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

#include "scene.h"

#include "Common.h"
#include "Impact.h"
#include "SamplingProcess.h"


SamplingStep::SamplingStep(int times,int eflags,HemicubeType ht,int sflags,double pEnerg,double pPartic)
{
  nApp = times;
  emissionFlags = eflags;
  sampleFlags = sflags;
  hemiType = ht;
  es = EmissionStrategy::newEmissionStrategy(eflags,ht);
  ss = SamplingStrategy::newSamplingStrategy(sflags);
  partEnergy = pEnerg;
  partParticles = pPartic;
}

SamplingStep::SamplingStep()
{
  nApp = 1;
  emissionFlags = 0;
  sampleFlags = 0;
  hemiType = HEMICUBE_NONE;
  es = NULL;
  ss = NULL;
  partEnergy = 1.0;
  partParticles = 1.0;
}

SamplingStep::~SamplingStep()
{
  if (es != NULL)
    delete es;

  if (ss != NULL)
    delete ss;
}


char* SamplingStep::getDescription()
{
  char tmp[1024];

  
  sprintf(tmp,"Step (%i,%.2f,%.2f) - Emission : (%s%s%s%s%s)+(%s) - Reflection : (%s%s%s%s%s%s)",
	  nApp,partParticles,partEnergy,
	  ISSET(emissionFlags,DE_UNIFORM_EMISSION) ? "Uniform " : "",
	  ISSET(emissionFlags,DE_COSTHETA_EMISSION) ? "Cosinus " : "",
	  ISSET(emissionFlags,DE_AREA_EMISSION) ? "Area " : "",
	  ISSET(emissionFlags,DE_DIRECT_POTENTIAL_EMISSION) ? "Direct potential " : "",
	  ISSET(emissionFlags,DE_EQUI_WIN_EMISSION) ? "Equi win " : "",
	  (hemiType==HEMICUBE_NONE) ? "None" : ( (hemiType==HEMICUBE_VISIBILITY) ? "Visibility" : ( (hemiType==HEMICUBE_PIXEL_COUNT) ? "Pixel count" : "Form factor" )),
	  ISSET(sampleFlags,DE_UNIFORM_SAMPLE) ? "Uniform " : "",
	  ISSET(sampleFlags,DE_COSTHETA_SAMPLE) ? "Cosinus " : "",
	  ISSET(sampleFlags,DE_BRDF_SAMPLE) ? "Brdf " : "",
	  ISSET(sampleFlags,DE_AREA_SAMPLE) ? "Area " : "",
	  ISSET(sampleFlags,DE_DIRECT_POTENTIAL_SAMPLE) ? "DirectPotential " : "",
	  ISSET(sampleFlags,DE_EQUI_WIN_SAMPLE) ? "Equi win " : ""
	  );

  return strdup(tmp);
}



SamplingProcess::SamplingProcess() : Vector<SamplingStep>()
{}



SamplingProcess::~SamplingProcess()
{
  while (getSize() != 0)
    {
      delete remove(0);
    }
}


void SamplingProcess::doSampling()
{
  SamplingStep *ss;
  PATCHLIST    *pl;
  COLOR        emit;
  int          n,i;
  char         *str;

  // Find total power to distribute
  totAveragePower = 0.0;
  pl = Patches;
  ForAllPatches(cp,pl) {
    emit = EdfEmittance(cp->surface->material->edf,
			&(cp->midpoint), ALL_COMPONENTS);
    totAveragePower += COLORAVERAGE(emit) * cp->area;
  } EndForAll;


  // Get number of steps
  n = getSize();

  // Accumulate fraction of power and fraction of energy
  totEnergy = totParticles = 0.0;
  for (i=0;i<n;i++)
    {
      ss = get(i);
      totEnergy += ss->partEnergy;
      totParticles += ss->partParticles;
    }

  // Do all steps
  for (i=0;i<n;i++)
    {
      ss = get(i);
      str = ss->getDescription();
      printf("Step %i\n%s\n",i,str);
      delete str;
      doStep(ss);
    }
}



void SamplingProcess::doStep(SamplingStep *ss)
{
  PATCHLIST                   *listPatch;
  COLOR                       emit;
  double                      fractionForPatch;
  double                      fractionEnergy;
  double                      numToSendForStep;
  double                      numToSend;
  double                      times;
  int                         i;

  // Convert to double
  times = (double) ss->nApp;

  // Calculate the number of particles to distribute at each sub-step.
  numToSendForStep = ( (double)deState.totNP * ss->partParticles ) / (totParticles*times);

  // Calculate the fraction of energy to distribute at each sub-step.
  fractionEnergy = ss->partEnergy / (totEnergy*times);

  // Do each sub-step
  for (i=0;i<ss->nApp;i++)
    {
      // Print substep info
      printf("Substep %i\n",i);

      // Iterate the patch list to find light sources
      listPatch = Patches; 
      ForAllPatches (curPatch,listPatch)
	{
	  emit = EdfEmittance(curPatch->surface->material->edf,
			      &(curPatch->midpoint), ALL_COMPONENTS);
	  if (!COLORNULL(emit))
	    {
	      // Calculate the number of particles to send from curPatch
	      fractionForPatch = COLORAVERAGE(emit) * curPatch->area  / totAveragePower;
	      numToSend = floor( fractionForPatch * numToSendForStep );

	      // Print some information
	      printf("Sending from patch id %i ...\n" 
		     "Fraction of energy to distribute: %.3f\n"
		     "Number of particles to send: %.0f\n",
		     curPatch->id,fractionEnergy*fractionForPatch,numToSend);

	      // Update emission strategy
	      ss->es->Update(curPatch);
	 
	      // Update sampling strategy of this step
	      ss->ss->Update();
	      
	      // Process the current patch
	      doPatch(ss,curPatch,numToSend,fractionEnergy);
	    }
	} EndForAll;
    }
}




void SamplingProcess::doPatch(SamplingStep* ss,PATCH *patch,double numtosend,double fractionEnergy)
{
  DENSITY_ESTIMATION_DATA     *thisData;
  PATCH                       *thisPatch;
  BSDF                        *thisBsdf;
  RAY                         rin,rout;
  HITREC                      hin,hout;
  COLOR                       weight,ctmp;
  VECTOR                      vtmp;
  double                      pdf,weightCoef;
  double                      phi,theta,cosTheta;

  VEC2D                       v2dtp;
  DVEC2D		      dv2dtp;

  // find constant factor
  weightCoef = (fractionEnergy * patch->area) / numtosend;

  // Send the particles from the patch
  while (numtosend>0.0)
    {
      // Choose a ray.
      pdf = ss->es->SendRay(patch,rin,hin);

      // Choose failed
      if (pdf == 0.0)
	continue;

      // Check to see if we are on the good side
      cosTheta = VECTORDOTPRODUCT(rin.dir,patch->normal);
      if (cosTheta<=0.0)
	continue;

      // Compute initial weight.
      weight = EdfEval(patch->surface->material->edf,
		       &(patch->midpoint),
		       &rin.dir, &(patch->normal), ALL_COMPONENTS);
      pdf = (cosTheta * weightCoef) / pdf;
      COLORSCALE(pdf,weight,weight);
  
      while (hin.patch != NULL)
	{
	  if (hin.flags & HIT_BACK)
	    {
	      // Humm... No time to properly implement BTDFs...
	      // SO THIS IS POTENTIALLY 'WRONG' !!!
	      // Try to ignore and continue if potentially transparent material
	      // if (hin.patch->surface->material->btdf != NULL)
	      //
	      // Another try (jp)
	      if (! COLORNULL(
                    BsdfTransmittance(hin.patch->surface->material->bsdf,
				      /*&(hin.patch->midpoint)*/ NULL, NULL /*TODO*/, ALL_COMPONENTS))
		  )
		{
		  double tp;
		  
		  rin.pos = hin.point;
		  rin.dir = SampleHemisphereCosTheta(&( ((DENSITY_ESTIMATION_DATA *) (hin.patch->radiance_data))->coordsys),drand48(),drand48(),&tp);
		  RayCastingFunction(hin.patch,&rin,hin);
		  continue;
		}
	      else
		break;
	    }

	  // Access data more easily.
	  thisPatch = hin.patch;
	  thisBsdf  = thisPatch->surface->material->bsdf;
	  thisData  = (DENSITY_ESTIMATION_DATA *) (thisPatch->radiance_data);

	  // Hitting the back of a surface or a patch that does not have a brdf...
	  if (thisBsdf == NULL)
	    {
	      // Humm... No time to properly implement BTDFs...  SO
	      // THIS IS POTENTIALLY 'WRONG' !  Try to ignore and
	      // continue if potentially transparent material.	      
	      // if (thisPatch->surface->material->btdf != NULL)
	      //
	      // Another try (jp)
	      if (! COLORNULL(
                    BsdfTransmittance(hin.patch->surface->material->bsdf,
				      /*&(hin.patch->midpoint),*/ 
				      NULL, NULL /*TODO*/, ALL_COMPONENTS))
		  )
		{
		  VECTOR opn;
		  COORDSYS cs;
		  double tp;
		  
		  VECTORSUBTRACT(origin,thisPatch->normal,opn);
		  VectorCoordSys(&opn,&cs);
		  
		  rin.pos = hin.point;
		  rin.dir = SampleHemisphereCosTheta(&cs,drand48(),drand48(),&tp);
		  RayCastingFunction(hin.patch,&rin,hin);
		  continue;
		}
	      else
		break;
	    }

	  // Store impact position.
	  PatchUV(thisPatch,&(hin.point),&dv2dtp.u, &dv2dtp.v);
	  v2dtp.u = dv2dtp.u; v2dtp.v = dv2dtp.v; 
	  impact.setCoordinates(v2dtp);

	  
	  // Store impact direction
	  VECTORSUBTRACT(origin,rin.dir,vtmp);
	  VectorToSphericalCoord(&vtmp,&(thisData->coordsys),&phi,&theta);
	  v2dtp.u = phi; v2dtp.v = theta;
	  impact.setDirection(v2dtp);	      

	  // Store impact weight
	  impact.setWeight(weight);
	  
	  // Record in impact store.
	  deState.isp->recordImpact((thisPatch->id)-1,impact);
	  
	  // Probalistically terminate the particle using russian roulette 
	  if (COLORAVERAGE(weight)<0.000001)
	    {
	      // Do we terminate the particle ?
	      pdf = drand48();
	      if (pdf<0.5)
		break;
	      else
		COLORSCALEINVERSE(1.0-0.5,weight,weight);
	    }

	  // Use particle bouncer to bounce the particle.
	  if ( (pdf = ss->ss->BounceRay(rin,hin,rout,hout)) == 0.0 )
	    break;

	  // Evaluate the weight of the particle.

	  // 1 --> Evaluate cosTheta / pdf .
	  cosTheta = VECTORDOTPRODUCT(rout.dir,thisPatch->normal);
	  if (cosTheta < 0.0) break;
	  COLORSCALE(cosTheta / pdf,weight,weight);

	  // 2 --> Evaluate the brdf.
	  ctmp = BsdfEval(thisBsdf, &(hin.patch->midpoint),
			  NULL, NULL,
			  &rin.dir,&rout.dir,&(hin.normal), ALL_COMPONENTS);
	  COLORPROD(ctmp,weight,weight);
	
	  // Get ready for another bounce
	  hin = hout;
	  rin = rout;
	}
      
      numtosend -= 1.0;
    }
}