#include <stdlib.h>
#include <math.h>
#include "lsix.h"

#include "ETriangle.H"
#include "SphTri.H"

#define SIGN(x) ((x) > 0.0 ? 1 : -1)

const int next[3] = {1,2,0} ;
const int prev[3] = {2,0,1} ;

#define INVALID(x) ((x)!=(x))


//********************************************************

ETriangle::ETriangle ()
{
  // void 
}

//********************************************************

ETriangle::ETriangle
( 
   const Vector & cVertex0 ,
   const Vector & cVertex1 ,
   const Vector & cVertex2 ,
   Scalar         iEpsilon 
) 
:  epsilon( iEpsilon )
{
   vertex[0]= cVertex0 ;
   vertex[1]= cVertex1 ;
   vertex[2]= cVertex2 ;

   UpdateNormals() ;
   UpdateArcsLengths() ;
   UpdateFormFactor() ;
   UpdateSolidAngle() ;
}

//*************************************************

const Vector & ETriangle::operator []
(  int i
)  const  
{
   return vertex[i] ;
}

//************************************************

void ETriangle::UpdateFromParent
(
   const ETriangle *  parent  ,
   const ETriangle *  intern  ,
   int                i
)
{
     epsilon = parent->epsilon ;

     int ip = prev[i] ,
         in = next[i] ;

     vertex[i]  = parent->vertex[i]  ;
     vertex[ip] = intern->vertex[ip] ;
     vertex[in] = intern->vertex[i]  ;

     normal[i]  = parent->normal[i]  ;
     normal[ip] = parent->normal[ip] ;
     normal[in] = intern->normal[ip] ;

     normal[in] = -1.0 * normal[in]  ;

     sideArcLen[i]  =  0.5 * parent->sideArcLen[i]  ;
     sideArcLen[ip] =  0.5 * parent->sideArcLen[ip] ; 
     sideArcLen[in] =  intern->sideArcLen[ip] ;

     inteArcLen[i]  = parent->inteArcLen[i] ; 
     inteArcLen[in] = M_PI - acos( normal[i]  | normal[in]  ) ;
     inteArcLen[ip] = M_PI - acos( normal[in] | normal[ip]  ) ;

     UpdateFormFactor() ;
     UpdateSolidAngle() ;  

 //   Draw() ;
}

//*************************************************

void ETriangle::UpdateNormals() 
{
   for( int i= 0 ; i<3 ; i++ )
      normal[i]= ( vertex[i] * vertex[next[i]] ).Normalized() ;
      
}

//*************************************************

void ETriangle::UpdateArcsLengths()
{
   for( int i= 0 ; i<3 ; i++ )
   {
       sideArcLen[i] = acos( vertex[i] | vertex[next[i]] ) ;

       if ( INVALID( sideArcLen[i])  )
       {  cout << endl << "cos == " << ( vertex[i] | vertex[next[i]] ) ;
       }

       inteArcLen[i] = M_PI - acos( normal[prev[i]] | normal[i] ) ; 
   }
}

//*************************************************

void ETriangle::UpdateFormFactor()
{
   
   formFactor =
     +  sideArcLen[0] * normal[0][Z]  
     +  sideArcLen[1] * normal[1][Z]
     +  sideArcLen[2] * normal[2][Z] ;

   formFactor = fabs( 0.5 * formFactor ) ;
}

//**************************************************

void ETriangle::UpdateSolidAngle()
{
   solidAngle =   inteArcLen[0] 
                + inteArcLen[1] 
                + inteArcLen[2] 
                - M_PI ;
}

//**************************************************

ETriangle * ETriangle::GetSubTriangles( ) const
{

  ETriangle * et = new ETriangle[4] ;

  Vector ivertex[3] ;

  ivertex[0]= ( vertex[0] + vertex[1] ).Normalized() ;
  ivertex[1]= ( vertex[1] + vertex[2] ).Normalized() ;
  ivertex[2]= ( vertex[2] + vertex[0] ).Normalized() ;

  ETriangle internal( ivertex[0], ivertex[1], ivertex[2], epsilon ) ;

  et[0].UpdateFromParent( this, &internal, 0 ) ;
  et[1].UpdateFromParent( this, &internal, 1 ) ;
  et[2].UpdateFromParent( this, &internal, 2 ) ;

  et[3] = internal ;

  return et ;
}

//**************************************************

Vector ETriangle::GetMaxArc
(
   int i0 
)  const
{
   int             i1 = next[i0] ;
   const Vector &  n  = normal[i0] ;
   Scalar          r  = sqrt( 1.0 - n[Z]*n[Z] ) ;
   

   Vector maxArc = Vector( (-n[X]*n[Z])/r, (-n[Y]*n[Z])/r, r ) ;

   Scalar z0 = CrossProdZ( vertex[i0], maxArc ) ,
          z1 = CrossProdZ( vertex[i1], maxArc ) ;

   if ( SIGN(z0) != SIGN(z1) ) 
      return maxArc ;
   else
   if ( vertex[i0][Z] > vertex[i1][Z] )
      return vertex[i0] ;
   else
      return vertex[i1] ;
          
}

//**************************************************

Vector ETriangle::GetMaxBoundary
( 
)  const
{
   Vector  maxArc[3] , 
           maxP      ;
   int     i         ;

   maxArc[0]= GetMaxArc( 0 ) ;
   maxArc[1]= GetMaxArc( 1 ) ;
   maxArc[2]= GetMaxArc( 2 ) ;

   maxP = maxArc[0] ;

   if ( maxArc[1][Z] > maxP[Z] ) maxP = maxArc[1] ; 
   if ( maxArc[2][Z] > maxP[Z] ) maxP = maxArc[2] ;

   return maxP ;

}

//**************************************************

Vector ETriangle::GetMax
( 
)  const
{
   
   int s[3] ;

   s[0] = SIGN( normal[0][Z] ) ;
   s[1] = SIGN( normal[1][Z] ) ;
   s[2] = SIGN( normal[2][Z] ) ;

   if ( s[0] == s[1] && s[1] == s[2] )
      return Vector( 0.0, 0.0, 1.0 ) ;
   else
      return GetMaxBoundary() ;
}

//*************************************************

void ETriangle::GetSamplesSA
(  
   unsigned int n  ,
   Sample       samples[]
)  const
{
   SphericalTriangle st( vertex[0], vertex[1], vertex[2] ) ;

   unsigned rn = (unsigned) floor( sqrt((double)n) ) ;
   
   unsigned fullStratified = ((rn*rn == n) && (n>1))  ;

 
   Scalar invSa = 1.0 / solidAngle ;

   for( int i= 0 ; i < n ; i++ )
   {
      Scalar u,v ;

      if ( fullStratified )
      {
          unsigned iu = i % rn ,
                   iv = i / rn ;

          u = ((double)iu + drand48()) / (double) rn ;
          v = ((double)iv + drand48()) / (double) rn ;
      }
      else
      {
         u = ((double)i + drand48()) / (double) n ;
         v = drand48() ;
      }

      samples[i].position      = st.Sample( u,v ) ;
      samples[i].probability   = invSa / samples[i].position[Z] ;
   }

   return ;
}

//************************************************

void ETriangle::GetSamples
(  
   unsigned int n  ,
   Sample       samples[]
)  const
{

   Vector max = GetMax() ;

   Scalar k = ( ( GetSolidAngle()* max[Z] )  / GetFormFactor() ) - 1.0 ;

   if ( k < epsilon )
      GetSamplesSA( n , samples ) ;
   else
      GetSamplesAdap( n , samples ) ;
}

//************************************************

void ETriangle::GetSamplesAdap
(  
   unsigned int n  ,
   Sample       samples[]
)  const
{

   ETriangle *   et     = GetSubTriangles() ;
   Scalar        ff[4]                      ,
                 sumff                      ;
   unsigned int  i                          ,
                 ns[4]                      ;

   for( i= 0 ; i < 4 ; i++ )
   {   ff[i]= et[i].GetFormFactor() ;
       sumff += ff[i] ;
   }

   GetNSamples( n, 4, ff, ns ) ;
   
   unsigned int base = 0 ;

   for( i= 0 ; i<4 ; i++ )
   {
      if ( ns[i] > 0 )
      {   
          et[i].GetSamples( ns[i], samples + base ) ;

          Scalar fi = ff[i] / sumff ;

          for ( unsigned j = base ; j < base + ns[i] ; j++ )
              samples[j].probability *= fi ;
          
          base += ns[i] ;
      }
   }

   delete[] et ;
}

//********************************************************

void GetNSamples
(
   unsigned    nSamples  ,
   unsigned    nEntries  ,
   double      prob[]    ,
   unsigned    ns[]
)
{
   double     sumProb = 0.0    ,
              sumRem  = 0.0    ;
   double *   rem     = new Scalar[nEntries] ;

   unsigned   i                ,
              remN   = nSamples       ;


   for( i= 0 ; i < nEntries ; i++ )
      sumProb += prob[i] ;
   

   for( i= 0 ; i < nEntries ; i++ )     
   {
      double exact   = prob[i] * (Scalar)nSamples / sumProb ,
             rounded = floor( exact )    ;

      ns[i]  = (unsigned int) rounded ;
      rem[i] = exact - rounded ;

      remN   -= ns[i] ;
      sumRem += rem[i] ;
   }

   while ( remN > 0 )
   {
      Scalar       r  = drand48() * sumRem ;
      unsigned int i  = 0 ;
 
      while ( r > rem[i] )
      {
          r -= rem[i] ;
          i++ ;
      }

      if ( i >= nEntries )
      {
         cerr << endl << "*** something's wrong ..... *** "
              << endl << " i        = " << i
              << endl << " sumRem   = " << sumRem
              << endl << " nSamples = " << nSamples
              << endl << flush ; 

         for( int j= 0 ; j < nEntries ; j++ )              
         {   
             cerr << endl << " rem[" << j << "] = " << rem[j] << " ns = " << ns[j] << " prob = " << prob[j] ;  
         }
         exit(1) ;
      }
         
      ns[i]++ ;
      remN-- ;
      sumRem-= rem[i] ;    
      rem[i]= 0.0 ;
   }

   delete [] rem ;
}


//******************************************************************

ostream & operator << 
( 
   ostream &         os ,  
   const ETriangle & et 
)
{
    os << endl 
       << "-------------------------------" << endl
       << "i    Arc Len.  Int.Angle NorLen" << endl
       << "-------------------------------" << endl ;


    for( int i = 0 ; i < 3 ; i++ )
    {
       os << i << "  " << et.sideArcLen[i] 
               << "  " << et.inteArcLen[i] 
               << "  " << et.normal[i] 
               << endl ;
    }

    SphericalTriangle st( et.vertex[0], et.vertex[1], et.vertex[2] ) ;


    os  << " -------------------------------" << endl
        << " Solid Angle = " << et.GetSolidAngle() << endl 
        << " S.A. (Arvo) = " << st.SolidAngle() << endl
        << " Form Factor = " << et.GetFormFactor() << endl
        << " Max.Z       = " << et.GetMax() << endl
        << " -----------------------------------" << endl
        << flush ; 

    return os ;
}