#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <signal.h>
#include <time.h>
#include <math.h>

#include "brdf.h"
#include "mymath.h"
#include "patch.h"
#include "phong.h"
#include "pools.h"
#include "polygon.h"
#include "ray.h"
#include "render.h"
#include "scene.h"
#include "spherical.h"
#include "vertex.h"

#include "dest.h"
#include "destP.h"

#include "Common.h"
#include "EmissionStrategy.h"
#include "Impact.h"
#include "ImpactStore.h"
#include "PatchProbabilityAssigner.h"
#include "SamplingProcess.h"
#include "SamplingStrategy.h"
#include "ScratchPad.h"

#include "options.h"




#ifndef NOPOOLS
static POOL *dataPool = (POOL*)NULL;
#define NEWDATA()           (DENSITY_ESTIMATION_DATA_PTR)NewPoolCell(sizeof(DENSITY_ESTIMATION_DATA), 0, "density estimation data", &dataPool);
#define DISPOSEDATA(data)   Dispose((char*)data,&dataPool)
#else
#define NEWDATA()           (DENSITY_ESTIMATION_DATA_PTR)Alloc(sizeof(DENSITY_ESTIMATION_DATA))
#define DISPOSEDATA(data)   Free((char*)data,sizeof(DENSITY_ESTIMATION_DATA))
#endif




// Initialize the state of Density Estimation to default values...
static void DestDefaults(void)
{
  // Standard values
  deState.particleTracingMode = DEFAULT_DE_PARTICLE_TRACING_MODE;
  deState.totNP               = DEFAULT_DE_TOTNP;
  deState.kernelShape         = DEFAULT_DE_KERNEL_SHAPE;
  deState.simplifyMode        = DEFAULT_DE_SIMPLIFY_MODE;
  deState.c1                  = DEFAULT_DE_C1;
  deState.maxSubLevel         = DEFAULT_DE_MAX_SUB_LEVEL;
  deState.c2                  = DEFAULT_DE_C2;
  deState.maxSubLevelD        = DEFAULT_DE_MAX_SUB_LEVEL_D;
  deState.directionnal        = DEFAULT_DE_DIRECTIONNAL;
  deState.edgeCompensate      = DEFAULT_DE_EDGE_COMPENSATE;
  deState.largeKernel         = DEFAULT_DE_LARGE_KERNEL;
  deState.simplifyThreshold   = DEFAULT_DE_SIMPLIFY_THRESHOLD;

  // Initialise sampling processes
  if (deState.simpleBRDFSP.getSize()==0)
    deState.simpleBRDFSP.insert(new SamplingStep(1,DE_COSTHETA_EMISSION,HEMICUBE_NONE,DE_BRDF_SAMPLE,1.0,1.0),0);

  if (deState.adaptiveEquiWinSP.getSize()==0)
    {
      deState.adaptiveEquiWinSP.insert(new SamplingStep(1,DE_COSTHETA_EMISSION,HEMICUBE_NONE,DE_BRDF_SAMPLE,0.5,0.5),0);
      deState.adaptiveEquiWinSP.insert(new SamplingStep(10,DE_UNIFORM_EMISSION|DE_EQUI_WIN_EMISSION,HEMICUBE_PIXEL_COUNT,DE_BRDF_SAMPLE,0.5,0.5),1);
    }

  if (deState.expertSP.getSize()==0)
    {
      deState.expertSP.insert(new SamplingStep(1,DE_COSTHETA_EMISSION,HEMICUBE_NONE,DE_BRDF_SAMPLE,0.5,0.5),0);
      deState.expertSP.insert(new SamplingStep(5,DE_UNIFORM_EMISSION|DE_EQUI_WIN_EMISSION,HEMICUBE_PIXEL_COUNT,DE_BRDF_SAMPLE,0.5,0.5),1);  
    }
}

/* command line options */
static CMDLINEOPTDESC DestOptions[] = {
  {NULL	, 	0,	TYPELESS, 	NULL, 		DEFAULT_ACTION,
   NULL }  	/* sentinel */
};

static void ParseDestOptions(int *argc, char **argv)
{
  ParseOptions(DestOptions, argc, argv);
}

static void PrintDestOptions(FILE *fp)
{
  fprintf(fp, "\nDensity Estimation Options:\n");
  PrintOptions(fp, DestOptions);
}

void RenderPatchDensityEstimation(PATCH *patch)
{
  Surface *sp;

  sp = (((DENSITY_ESTIMATION_DATA *)patch->radiance_data)->swptr);
  if (sp==NULL)
    {
      POINT pt[4];
      RGB   black;
      int   i;
      
      /*
      RGBSET(black,1.0,1.0,1.0);
      RenderSetColor(&black);
      for(i=0;i<patch->nrvertices;i++)
	{
	  VECTORCOPY(*(patch->vertex[i]->point),pt[i]);
	}
      RenderPatchOutline(patch);
      */
      
      RGBSET(black,0.0,0.0,0.0);
      RenderSetColor(&black);
      for(i=0;i<patch->nrvertices;i++)
	{
	  VECTORCOPY(*(patch->vertex[i]->point),pt[i]);
	}
      RenderPolygonFlat(patch->nrvertices,pt);
      
    }
  else
    sp->display();
}



// Initialize the Density Estimation algorithm...
void InitDensityEstimation(void)
{
  DENSITY_ESTIMATION_DATA  *data;
  PATCHLIST                *pl;

  // Initialize patch-associated data
  pl = Patches;
  ForAllPatches(patch,pl)
    {
      data = (DENSITY_ESTIMATION_DATA *)patch->radiance_data; // Get density estimation data 
      PatchCoordSys(patch,&(data->coordsys));                 // Coordinate system for sampling   
      data->swptr = NULL;                                     // Final result data
    }
  EndForAll;

  // Initialize timer
  deState.cpuParticleTracing = 0;
  deState.cpuReconstruction  = 0;
  deState.bytesResult        = 0;
  deState.lastClock          = clock();

  // Initialize to first phase
  deState.inWork             = FALSE;
  deState.curPhase           = DENSITY_ESTIMATION_PHASE_INITIALISATION;
  deState.isp                = NULL;
  deState.hsp                = NULL;
  deState.spp                = NULL;
}


void StartDensityEstimation()
{
  // Initialize storage method
  deState.isp = new ImpactStore(PatchListCount(Patches));

  // Check initialisation
  if (deState.isp == NULL)
    fprintf(stderr,"Density Estimation initialisation failed.\n"); 
  else
    {
      // Step to next phase
      deState.curPhase     = DENSITY_ESTIMATION_PHASE_PARTICLE_TRACING;
      deState.inWork       = TRUE;
    }
} 


void DoParticleTracing()
{
  // Execute the selected sampling process
  switch(deState.particleTracingMode)
    {
    case DE_PT_SIMPLE:
      deState.simpleBRDFSP.doSampling();
      break;
      
    case DE_PT_ADAPTIVE:
      deState.adaptiveEquiWinSP.doSampling();
      break;

    case DE_PT_EXPERT:
      deState.expertSP.doSampling();
      break;
    }

  // Tell the impact store we are finished with recording...
  deState.isp->flush();

  // Alloc hemisphere scratchpad if we are reconstruction directionnal information
  if (deState.directionnal==DE_MULTIDIRECTION_RECONSTRUCTION)
    deState.hsp = new HemisphereScratchPad(deState.maxSubLevel,deState.maxSubLevelD);
  
  // Alloc a scratchpad and assign a kernel shape. 
  deState.spp = new ScratchPad(deState.maxSubLevel);
  deState.spp->setKernel(Kernel::newKernel(deState.kernelShape));

  // Prepare density estimation / decimation cycle
  deState.curPhase     = DENSITY_ESTIMATION_PHASE_RECONSTRUCTION;
}


void DoReconstructionStep()
{
  DENSITY_ESTIMATION_DATA     *data;
  int                         id;
  int                         n;
  int                         hasGlossy;
  int                         hasSpecular;
  int                         kernelSmallEnough;
  float                       chwsz;    /* Current window size for position  */
  float                       chgsz;    /* Current grid size for position */
  float                       chwszD;   /* Current window size for direction */
  Impact                      impact;
  ImpactStoreWindow           *isw;
  COLOR                       emit,ctmp,col;
  VEC2D                       cood,dir;
  VECTOR                      midpoint;
  BSDF                       *currentBSDF;

  // Initiate
  deState.curPatchList = Patches;
  deState.curPatch = PatchListNext(&deState.curPatchList);

  // Reconstruct and compress one patch info
  while (deState.curPatch != NULL)
    {
      // Easy access
      data = (DENSITY_ESTIMATION_DATA_PTR) (deState.curPatch->radiance_data);
      id = deState.curPatch->id-1;
      midpoint = deState.curPatch->midpoint;
      currentBSDF = deState.curPatch->surface->material->bsdf;

      // Get emission properties
      emit = EdfEmittance(deState.curPatch->surface->material->edf,
			  &midpoint, ALL_COMPONENTS);

      // If a patch received no impact, we don't care about this patch ...
      if ( (n = deState.isp->getNumberOfImpacts(id)) == 0)
	{
	  // ... except if it's a light source !
	  if (!COLORNULL(emit))
	    {
	      COLORSCALE(M_1_PI,emit,emit);
	      data->swptr = new ConstantSurface(deState.curPatch,emit);
	    }

	  deState.curPatch = PatchListNext(&deState.curPatchList);
	  continue;
	}

      // Choose kernel size
      chwsz  = sqrt( (deState.curPatch->area*deState.c1) / (float)n );
      chgsz  = sqrt( (deState.curPatch->area*deState.c2) / (float)n );
      chwszD = sqrt( (M_2PI*deState.c1) / (float) n );

      chwszD = M_PI_2/4.0;

      // Is the directionnal kernel small enough ?
      kernelSmallEnough = (chwszD<M_PI_2);

      // The patch has glossy component ?
      ctmp = BsdfReflectance(currentBSDF, NULL/*TODO*/, &midpoint, GLOSSY_COMPONENT);
      hasGlossy = ! COLORNULL(ctmp);

      // The patch has specular component ?
      ctmp = BsdfReflectance(currentBSDF, NULL/*TODO*/, &midpoint, SPECULAR_COMPONENT);
      hasSpecular = ! COLORNULL(ctmp);


      fprintf(stderr,
	      "Reconstructing patch %i receiving %i impacts "
	      "with a window size of h %.3f\n",
	      id, n, chwsz);

      if ((deState.directionnal == DE_MULTIDIRECTION_RECONSTRUCTION) &&
	  (kernelSmallEnough) &&
	  (hasGlossy))
	{
	  // Multi-direction reconstruction

	  // Configure ScratchPad.
	  deState.hsp->configure(deState.curPatch,chgsz,deState.maxSubLevelD);
	  deState.hsp->setKernelSize(chwsz,chwszD);

	  // Add the impacts.
	  isw = deState.isp->newImpactStoreWindow(id);
	  while (isw->readNext(impact)!=0)
	    {
	      // Unpack usefull information
	      impact.getCoordinates(cood);
	      impact.getDirection(dir);
	      impact.getWeight(col);

	      // Accumulate the impact
	      deState.hsp->addImpact(cood,dir,col);
	    }
	  delete isw;

	  if (deState.simplifyMode == DE_SIMPLIFY_DECIMATION)
	    data->swptr = deState.hsp->constructNonDiffuseDecimatedSurface(deState.simplifyThreshold);
	  else
	    data->swptr = deState.hsp->constructNonDiffuseCompressedSurface(deState.simplifyThreshold);
	}
      else
	if ((deState.directionnal==DE_UNIDIRECTION_RECONSTRUCTION)&&(kernelSmallEnough)&&(hasGlossy))
	  {
	    // Uni-direction reconstruction
	  
	    // Configure ScratchPad.
	    deState.spp->configure(deState.curPatch,chgsz);
	    deState.spp->setKernelSize(chwsz);
	    
	    // Add emission
	    COLORSCALE(M_1_PI,emit,emit);
	    deState.spp->setBaseRadiance(emit);
	    
	    // Add the impacts.
	    isw = deState.isp->newImpactStoreWindow(id);
	    while (isw->readNext(impact)!=0)
	      {
		// Unpack usefull information
		impact.getCoordinates(cood);
		impact.getDirection(dir);
		impact.getWeight(col);
	      
		// Calcul direction entree
		HemisphereScratchPad::addImpact(deState.spp,deState.maxSubLevelD,chwszD,cood,dir,col);
	      }
	    delete isw;

	    if (deState.simplifyMode == DE_SIMPLIFY_DECIMATION)
	      data->swptr = deState.spp->constructDiffuseDecimatedSurface(deState.simplifyThreshold);
	    else
	      data->swptr = deState.spp->constructDiffuseCompressedSurface(deState.simplifyThreshold);
	  }
	else
	  {
	    // No direction reconstruction
	  
	    // Configure ScratchPad.
	    deState.spp->configure(deState.curPatch,chgsz);
	    deState.spp->setKernelSize(chwsz);

	    // Add emission
	    COLORSCALE(M_1_PI,emit,emit);
	    deState.spp->setBaseRadiance(emit);
	    
	    // Find multiplicative color
	    // If has glossy, don't reconstruct specular 
	    // However, if has specular but non glossy, reconstruct all...
	    if (hasGlossy)
	      ctmp = 
		BsdfReflectance(currentBSDF, NULL/*&midpoint TODO*/, NULL/*TODO*/,
				DIFFUSE_COMPONENT|GLOSSY_COMPONENT);
	    else
	      {
		if (hasSpecular)
		  ctmp = BsdfReflectance(currentBSDF, NULL/*&midpoint,*/, NULL,/*TODO*/
					 DIFFUSE_COMPONENT);
		else
		  ctmp = BsdfReflectance(currentBSDF, NULL/*&midpoint,*/, NULL,/*TODO*/
					 DIFFUSE_COMPONENT);
	      }

	    // Irradiance -> Radiance
	    COLORSCALE(M_1_PI,ctmp,ctmp);

	    // Add the impacts.
	    isw = deState.isp->newImpactStoreWindow(id);
	    while (isw->readNext(impact)!=0)
	      {
		// Unpack usefull info from the impact
		impact.getCoordinates(cood);
		impact.getWeight(col);
		
		// Do the reflection
		COLORPROD(col,ctmp,col);

		// Accumulate the impact
		deState.spp->addImpact(cood,col);
	      }
	    delete isw;

	    // Simplify all
	    if (deState.simplifyMode == DE_SIMPLIFY_DECIMATION)
	      data->swptr = deState.spp->constructDiffuseDecimatedSurface(deState.simplifyThreshold);
	    else
	      data->swptr = deState.spp->constructDiffuseCompressedSurface(deState.simplifyThreshold);
	  }
         
      // Display the patch
      RenderPatchDensityEstimation(deState.curPatch);

      // Count bytes for stats
      if (data->swptr!=NULL)
	deState.bytesResult += data->swptr->getNumBytes();

      // Iterate
      deState.curPatch = PatchListNext(&deState.curPatchList);
    }

  // Nothing more to reconstruct, free data structures
  if (deState.directionnal==DE_MULTIDIRECTION_RECONSTRUCTION)
    {
      delete deState.hsp;
      deState.hsp = NULL;
    }

  // Also write a VRML file
  if (deState.simplifyMode == DE_SIMPLIFY_DECIMATION)
    {
      FILE *outFile;

      outFile = fopen("test.wrl","w");

      fprintf(outFile,"#VRML V2.0 utf8\n");
      deState.curPatchList = Patches;
      deState.curPatch = PatchListNext(&deState.curPatchList);
      while (deState.curPatch!=NULL)
	{
	  data = (DENSITY_ESTIMATION_DATA_PTR) (deState.curPatch->radiance_data);
	  if (data->swptr!=NULL)
	    data->swptr->writeVRML(outFile);
	  deState.curPatch = PatchListNext(&deState.curPatchList);
	}

      fclose(outFile);
    }

  deState.inWork = FALSE;
}



void TerminateDensityEstimation(void)
{
  if (deState.hsp != NULL)
    {
      delete deState.hsp;
      deState.hsp = NULL;
    }

  if (deState.spp != NULL)
    {
      delete deState.spp;
      deState.spp = NULL;
    }

  if (deState.isp != NULL)
    {
      delete deState.isp;
      deState.isp = NULL;
    }
}



// Do one iteration step of the Density Estimation algorithm...
int DoDensityEstimationStep(void)
{
  void                        (*prev_alrm_handler)(int signr); 
  unsigned                    prev_alarm_left;


  // Install a timer that will wake us up after one second for checking for 
  // user events.
  prev_alrm_handler = signal(SIGALRM,WakeUpDensityEstimation);
  prev_alarm_left = alarm(1);
  deState.wakeUp = FALSE;
  

  // Different phases, different works...
  switch(deState.curPhase)
    {
    case DENSITY_ESTIMATION_PHASE_INITIALISATION:
      // Initialize patch data
      StartDensityEstimation();
      break;

    case DENSITY_ESTIMATION_PHASE_PARTICLE_TRACING:
      // Do particle tracing
      DoParticleTracing();
      break;

    case DENSITY_ESTIMATION_PHASE_RECONSTRUCTION:
      // Do reconstruction and decimation/compression
      DoReconstructionStep();
      break;

    default:
      DebugAssert(0,"Incoherent state");
      break;
    }

  /* Reinstall the previous alarm handler. */
  signal(SIGALRM, prev_alrm_handler);
  alarm(prev_alarm_left);

  return !(deState.inWork); 
}



/* Update the CPU time. */
static void UpdateCpuSecs(void)
{
  float   f;
  clock_t t;

  t = clock();
  f = (float)(t - deState.lastClock)/(float)CLOCKS_PER_SEC;

  switch(deState.curPhase)
    {
    case DENSITY_ESTIMATION_PHASE_INITIALISATION:
      break;

    case DENSITY_ESTIMATION_PHASE_PARTICLE_TRACING:
      deState.cpuParticleTracing += f;
      break;

    case DENSITY_ESTIMATION_PHASE_RECONSTRUCTION:
      deState.cpuReconstruction += f;
      break;

    default:
      DebugAssert(0,"Unknown state.");
      break;
    }

  deState.lastClock = t;
}


/* Sets up a timer to wake us up every second. */
void WakeUpDensityEstimation(int)
{
  deState.wakeUp = TRUE;
  signal(SIGALRM,WakeUpDensityEstimation);
  alarm(1);
  UpdateCpuSecs();
}



static void *CreatePatchData(PATCH *patch)
{
  DENSITY_ESTIMATION_DATA_PTR data = NEWDATA();

  patch->radiance_data = data;
  return data;
}


static void PrintPatchData(FILE*,PATCH*)
{ }


static void DestroyPatchData(PATCH *patch)
{
  if (patch->radiance_data)
    {
      Surface *sp;

      sp = ((DENSITY_ESTIMATION_DATA*)patch->radiance_data)->swptr;
      if (sp!=NULL)
	delete sp;

      DISPOSEDATA(patch->radiance_data);
    }
  patch->radiance_data = (void*)NULL;
}




static COLOR GetRadianceDensityEstimation(PATCH *,double,double,VECTOR)
{ return zerorad; }


static char* GetStatsDensityEstimation(void)
{
  static char stats[1000];

  sprintf(stats,
	  "Density Estimation Statistics:\n\n" 
	  "Time spend in particle tracing : %.1f\n" 
	  "Time spend in reconstruction : %.1f\n"
          "Bytes needed to store the result : %i\n",
	  deState.cpuParticleTracing,deState.cpuReconstruction,deState.bytesResult);
  return stats;
}




static void RenderSceneDensityEstimation(void)
{
  if (!deState.inWork)
    { PatchListIterate(Patches,RenderPatchDensityEstimation); }
}



static void *CreateGeomData(GEOM *)
{ return NULL; }

static void PrintGeomData(FILE *,GEOM *)
{ }

static void DestroyGeomData(GEOM *)
{ }



RADIANCEMETHOD DensityEstimation = {
  "DensityEstimation", 4,
  "Density Estimation",
  "destButton",
  DestDefaults,
  CreateDestControlPanel,
  ParseDestOptions,
  PrintDestOptions,
  InitDensityEstimation,
  DoDensityEstimationStep,
  TerminateDensityEstimation,
  GetRadianceDensityEstimation,
  CreatePatchData,
  PrintPatchData,
  DestroyPatchData,
  (void (*)(PATCH *, float))NULL,
  CreateGeomData,
  PrintGeomData,
  DestroyGeomData,
  ShowDestControlPanel,
  HideDestControlPanel,
  GetStatsDensityEstimation,
  RenderSceneDensityEstimation
};