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

#include "Common.h"
#include "HemisphereScratchPad.h"


HemisphereScratchPad::HemisphereScratchPad(int rpos,int rdir)
{
  int                     i;

  DebugAssert(rdir>0,"Recursion level is a strictly positive number");
  
  maxr = rdir;
  maxnumsub = (1<<rdir);
  maxnumpar = maxnumsub + 1;
  maxnumspps = (maxnumsub*maxnumpar*2)+1;
  spps = new ScratchPad* [maxnumspps];
  
  for(i=0;i<maxnumspps;i++)
    spps[i] = new ScratchPad(rpos);
}


HemisphereScratchPad::~HemisphereScratchPad()
{
  int i;

  if (spps != NULL)
    {
      for (i=0;i<maxnumspps;i++)
	if (spps[i] != NULL)
	  delete spps[i];
      
      delete spps;
    }
}


void HemisphereScratchPad::configure(PATCH* p,const float grdsz,const int rdir)
{
  int i;
  
  curr = rdir;
  curnumsub = (1<<rdir);
  curnumpar = curnumsub + 1;
  curnumspps = (curnumsub*curnumpar*2)+1;

  chgsz = grdsz;

  for(i=0;i<curnumspps;i++)
    spps[i]->configure(p,grdsz);
}
  

void HemisphereScratchPad::setKernelSize(const float hp,const float hd)
{
  int   i;

  // Precompute some values
  h = (hd>M_PI_2) ? M_PI_2 : hd;
  h2 = h*h;
  h2inv = 1.0f / h2;

  for (i=0;i<curnumspps;i++)
    spps[i]->setKernelSize(hp);
}
  

void HemisphereScratchPad::addImpact(const VEC2D &pos,const VEC2D &dir,const COLOR &weight)
{
  float  tmp;
  float  dt,dp,deltat,deltap;
  int    i,d,it,ip,isr;
  VECTOR posd,posr,delta;
  COLOR  tw;


  // Find the vector of the impact
  to3d(dir.u,dir.v,posd);

  // Find the reflection, if any 
  tmp = M_PI_2 - dir.v;
  if (tmp<h)
    {
      isr = TRUE;
      tmp = M_PI - dir.v;
      to3d(dir.u,tmp,posr);
    }
  else
    isr = FALSE;

  // Compare it with the other ones
  i = 0;
  d = 1;
  dt = 0.0;
  deltat = M_PI_2 / ((float)(curnumsub));
  for (it=0;it<=curnumsub;it++)
    {
      dp = 0.0;
      deltap = M_2PI / ((float)d);
      for (ip=0;ip<d;ip++)
	{
	  to3d(dp,dt,delta);
	  
	  tmp = VECTORDOTPRODUCT(posd,delta);
	  if (tmp<0.0)
	    continue;
	  tmp = acos(tmp);

	  if (tmp<h)
	    {
	      tmp = (1.0-(tmp*tmp*h2inv));
	      tmp = M_2_PI*tmp*h2inv;
	      COLORSCALE(tmp,weight,tw);

	      if (isr)
		{
		  tmp = VECTORDOTPRODUCT(posr,delta);
		  if (tmp>0.0)
		    {
		      tmp = acos(tmp);
		      
		      if (tmp<h)
			{
			  tmp = (1.0-(tmp*tmp*h2inv));
			  tmp = M_2_PI*tmp*h2inv;
			  COLORADDSCALED(tw,tmp,weight,tw);
			}
		    }
		}
	      spps[i+ip]->addImpact(pos,tw);
	    }
	  dp += deltap;
	}
      
      i += d;
      d += (it==0) ? 3 : 4;
      dt += deltat;
    }

}


HemisphereScratchPad* HemisphereScratchPad::evaluateOutgoing()
{
  HemisphereScratchPad  *tp;
  int                   ia,ib;
  int                   da,db;
  int                   ita,itb;
  int                   ipa,ipb;
  float                 dta,dtb;
  float                 dpa,dpb;
  float                 deltata,deltatb;
  float                 deltapa,deltapb;
  float                 scale,fscale;
  VECTOR                dira,dirb,norm;
  COLOR                 col;


  // Initialise
  VECTORSET(norm,0.0,0.0,1.0);
  scale = M_2PI / ((float)curnumspps); 

  // Ajouter accum
  for (ib=0;ib<curnumspps;ib++)
    spps[ib]->addAccum();

  // Alloc another HemisphereScratchPad for temporary calculations.
  tp = new HemisphereScratchPad(spps[0]->curr,curr);
  tp->configure(spps[0]->patch,chgsz,curr);
  
  // Visit all directions
  ia = 0;
  da = 1;
  dta = 0.0;
  deltata = M_PI_2 / ((float)(curnumsub));
  for (ita=0;ita<=curnumsub;ita++)
    {
      fscale = cos(dta)*scale;
      dpa = 0.0;
      deltapa = M_2PI / ((float)da);
      for (ipa=0;ipa<da;ipa++)
	{
	  to3d(dpa,dta,dira);
	  VECTORSUBTRACT(origin,dira,dira);

	  // Contribute to all the other ones
	  ib = 0;
	  db = 1;
	  dtb = 0.0;
	  deltatb = M_PI_2 / ((float)(curnumsub));
	  for (itb=0;itb<=curnumsub;itb++)
	    {
	      dpb = 0.0;
	      deltapb = M_2PI / ((float)db);
	      for (ipb=0;ipb<db;ipb++)
		{
		  to3d(dpb,dtb,dirb);
	  
		  // dira and dirb are the two directions
		  col = BsdfEval(spps[0]->patch->surface->material->bsdf,
				 (&spps[0]->patch->midpoint),
				 NULL, NULL, &dira, &dirb, &norm,
				 DIFFUSE_COMPONENT|GLOSSY_COMPONENT);
		  COLORSCALE(fscale,col,col);

		  tp->spps[ib+ipb]->add(spps[ia+ipa],col);
      
		  dpb += deltapb;
		}
	      
	      ib += db;
	      db += (itb==0) ? 3 : 4;
	      dtb += deltatb;
	    }


	  dpa += deltapa;
	}
      
      ia += da;
      da += (ita==0) ? 3 : 4;
      dta += deltata;
    }

  return tp;
}



NonDiffuseSurface* HemisphereScratchPad::constructNonDiffuseCompressedSurface(const float t)
{
  HemisphereScratchPad *tp;
  NonDiffuseSurface    *nds;
  int                  i;

  // Reconstruct outgoing radiance
  tp = evaluateOutgoing();
 
  // Construct the non diffuse compressed surface
  nds = new NonDiffuseSurface(spps[0]->patch,curr);
  for (i=0;i<curnumspps;i++)
    nds->setDiffuseSurface(i,tp->spps[i]->constructDiffuseCompressedSurface(t));
    
  // Free temporary storage
  delete tp;

  return nds;
}



NonDiffuseSurface *HemisphereScratchPad::constructNonDiffuseDecimatedSurface(const float t)
{
  HemisphereScratchPad *tp;
  NonDiffuseSurface    *ndsp;
  int                           i;

  // Reconstruct outgoing radiance
  tp = evaluateOutgoing();

  // Construct the coefficient store array
  ndsp = new NonDiffuseSurface(spps[0]->patch,curr);
  for (i=0;i<curnumspps;i++)
    ndsp->setDiffuseSurface(i,tp->spps[i]->constructDiffuseDecimatedSurface(t));
    
  // Free temporary storage
  delete tp;

  return ndsp;
}





void HemisphereScratchPad::to3d(const float p,const float t,VECTOR& v)
{
  v.x = cos(p)*sin(t); 
  v.y = sin(p)*sin(t); 
  v.z = cos(t);        
}




int HemisphereScratchPad::numberOfPoints(int rdir)
{
  int nsub;

  DebugAssert(rdir>=0,"The recursion level is a positive number");

  nsub = (1<<rdir);
  return (nsub*(nsub+1)*2)+1;
}


int HemisphereScratchPad::indexDirection(int r,double phi,double theta)
{
  double convtheta;
  int    pi,t;


  DebugAssert(rdir>=0,"The recursion level is a positive number");

  // Be tolerant with numerical errors...
  if (phi<0.0)      phi = 0.0;
  if (phi>M_2PI)    phi = M_2PI;
  if (theta<0.0)    theta = 0.0;
  if (theta>M_PI_2) theta = M_PI_2;

  // Find sub-triangle edge length
  convtheta = ((double)(1<<r)) / M_PI_2;

  // Calculate various indexes
  t = (int) rint(theta*convtheta);                            // t index
  pi = (t==0) ? 0 : (int)rint(((double)(4*t))*phi*M_1_2PI);   // p index

  return pi + (2*t*(t-1));                                    // general index
}

void HemisphereScratchPad::addImpact(ScratchPad *spp,const int vr,const float vh,const VEC2D &pos,const VEC2D &dir,const COLOR &weight)
{
  float  tmp,angle,vh2inv;
  float  dt,dp,deltat,deltap;
  int    i,d,it,ip,isr;
  VECTOR indir,indirrefl,outdir,indirtp,delta;
  COLOR  tw,twtp;  
  int    vnumsub;

  // Find the outgoing direction of the impact.
  PatchPoint(spp->patch,pos.u,pos.v,&outdir);
  VECTORSUBTRACT(Camera.eyep,outdir,outdir);
  if (VECTORDOTPRODUCT(spp->patch->normal,outdir)<=0.0)
    return;
  VECTORNORMALIZE(outdir);

  // Find how far we discretise
  vnumsub = (1<<vr);
  vh2inv = 1.0f/(vh*vh);

  // Clear accumulation buffer
  COLORCLEAR(tw);

  // Find the base incoming direction of the impact
  to3d(dir.u,dir.v,indir);

  // Find the base incoming direction reflection of the impact, if any 
  tmp = M_PI_2 - dir.v;
  if (tmp<vh)
    {
      isr = TRUE;
      tmp = M_PI - dir.v;
      to3d(dir.u,tmp,indirrefl);
    }
  else
    isr = FALSE;

  // Compare it with the other ones
  i = 0;
  d = 1;
  dt = 0.0;
  deltat = M_PI_2 / ((float)(vnumsub));
  for (it=0;it<=vnumsub;it++)
    {
      dp = 0.0;
      deltap = M_2PI / ((float)d);
      for (ip=0;ip<d;ip++)
	{
	  to3d(dp,dt,delta);
	  
	  angle = VECTORDOTPRODUCT(indir,delta);
	  if (angle<0.0)
	    continue;
	  angle = acos(angle);

	  if (angle<vh)
	    {
	      angle = (1.0-(angle*angle*vh2inv));
	      tmp = M_2_PI*angle*vh2inv;

	      if (isr)
		{
		  angle = VECTORDOTPRODUCT(indir,delta);
		  if (angle>0.0)
		    {
		      angle = acos(angle);
		      
		      if (angle<vh)
			{
			  angle = (1.0-(angle*angle*vh2inv));
			  tmp += M_2_PI*angle*vh2inv;
			}
		    }
		}
	      
	      VECTORSUBTRACT(delta,origin,indirtp);
	      twtp = BsdfEval(spp->patch->surface->material->bsdf,
			      &(spp->patch->midpoint),
			      NULL, NULL,
			      &indirtp, &outdir, &(spp->patch->normal),
			      DIFFUSE_COMPONENT|GLOSSY_COMPONENT);
	      COLORPROD(twtp,weight,twtp);

	      tmp *= M_1_PI;
	      COLORADDSCALED(tw,tmp,twtp,tw);
	    }
	  dp += deltap;
	}
      
      i += d;
      d += (it==0) ? 3 : 4;
      dt += deltat;
    }

  spp->addImpact(pos,tw);
}