/* gathering.c: Jocabi or Gauss-Seidel Galerkin radiosity */

#include "galerkinP.h"
#include "scene.h"
#include "vertex.h"
#include "camera.h"
#include "potential.h"
#include "statistics.h"
#include "error.h"
#include "coefficients.h"
#include "basis.h"

#include "interactionlist.h"

/* Lazy linking: delay creating the initial links for a patch until it has
 * radiance to distribute in the environment. Leads to faster feedback to the
 * user during the first Jacobi iterations. During the first iteration, only 
 * interactions with light sources are created. See Holschuch, EGRW '94 
 * (Darmstadt). */
static void PatchLazyCreateInteractions(PATCH *P)
{
  ELEMENT *el = TOPLEVEL_ELEMENT(P);

  if (!COLORNULL(el->radiance[0]) && !(el->flags & INTERACTIONS_CREATED)) {
    CreateInitialLinks(el, SOURCE);
    el->flags |= INTERACTIONS_CREATED;
  }
}

/* Converts the accumulated received radiance into exitant radiance, making the
 * hierarchical representation consistent and computes a new color for the patch. */
static void PatchUpdateRadiance(PATCH *P)
{
  PushPullRadiance(TOPLEVEL_ELEMENT(P));
  PatchRecomputeColor(P);
}

/* Creates initial interactions if necessary, recursively refines the interactions
 * of the toplevel element, gathers radiance over the resulting lowest level 
 * interactions and updates the radiance for the patch if doing 
 * Gauss-Seidel iterations. */
static void PatchGather(PATCH *P)
{
  ELEMENT *el = TOPLEVEL_ELEMENT(P);

  /* don't gather to patches without importance. This optimisation can not
   * be combined with lazy linking based on radiance. */
  if (gal.importance_driven && 
      el->potential.f < max_direct_potential * EPSILON)
    return;

  /* the form factors have been computed and stored with the source patch 
   * already before when doing non-importance-driven Jacobi iterations with lazy 
   * linking. */
  if (gal.iteration_method == GAUSS_SEIDEL || !gal.lazy_linking ||
      gal.importance_driven) {
    if (!(el->flags & INTERACTIONS_CREATED)) {
      CreateInitialLinks(el, RECEIVER);
      el->flags |= INTERACTIONS_CREATED;
    }
  }

  /* Refine the interactions and compute light transport at the leaves */
  RefineInteractions(el);

  /* Immediately convert received radiance into radiance, make the representation
   * consistent and recompute the color of the patch when doing Gauss-Seidel. 
   * The new radiance values are immediately used in subsequent steps of
   * the current iteration. */
  if (gal.iteration_method == GAUSS_SEIDEL)
    PatchUpdateRadiance(P);
}

/* Updates the potential of the element after a change of the camera, and as such
 * a potential change in directly received potential. */
void GatheringUpdateDirectPotential(ELEMENT *elem, float potential_increment)
{
  if (elem->regular_subelements) {
    int i;
    for (i=0; i<4; i++)
      GatheringUpdateDirectPotential(elem->regular_subelements[i], potential_increment);
  }
  elem->direct_potential.f += potential_increment;
  elem->potential.f += potential_increment;
}

/* Makes the representation of potential consistent after an iteration 
 * (potential is always propagated using Jacobi iterations). */
static float GatheringPushPullPotential(ELEMENT *elem, float down)
{
  float up;
  
  down += elem->received_potential.f / elem->area;
  elem->received_potential.f = 0.;

  up = 0.;

  if (!elem->regular_subelements && !elem->irregular_subelements) {
    up = down + elem->pog.patch->direct_potential;
  } 

  if (elem->regular_subelements) {
    int i;
    for (i=0; i<4; i++) 
      up += 0.25 * GatheringPushPullPotential(elem->regular_subelements[i], down);
  }

  if (elem->irregular_subelements) {
    ELEMENTLIST *subellist;
    for (subellist=elem->irregular_subelements; subellist; subellist=subellist->next) {
      ELEMENT *subel = subellist->element;
      if (!IsCluster(elem))
	down = 0.;	/* don't push to irregular surface subelements */
      up += subel->area / elem->area * GatheringPushPullPotential(subel, down);
    }
  }

  return (elem->potential.f = up);
}

static void PatchUpdatePotential(PATCH *patch)
{
  GatheringPushPullPotential(TOPLEVEL_ELEMENT(patch), 0.);
}

/* Recomputes the potential of the cluster and its subclusters based on the 
 * potential of the contained patches. */
static float ClusterUpdatePotential(ELEMENT *clus)
{
  if (IsCluster(clus)) {
    ELEMENTLIST *subcluslist;
    clus->potential.f = 0.;
    for (subcluslist=clus->irregular_subelements; subcluslist; subcluslist=subcluslist->next) {
      ELEMENT *subclus = subcluslist->element;
      clus->potential.f += subclus->area * ClusterUpdatePotential(subclus);
    }
    clus->potential.f /= clus->area;
  }
  return clus->potential.f;
}

/* Does one step of the radiance computations, returns TRUE if the computations
 * have converged and FALSE if not. */
int DoGatheringIteration(void)
{
  if (gal.importance_driven) {
    if (gal.iteration_nr<=1 || Camera.changed) {
      UpdateDirectPotential();
      Camera.changed = FALSE;
    }
  }

  /* non importance-driven Jacobi iterations with lazy linking */
  if (gal.iteration_method != GAUSS_SEIDEL && gal.lazy_linking &&
      !gal.importance_driven)
    PatchListIterate(Patches, PatchLazyCreateInteractions);

  /* no visualisation with ambient term for gathering radiosity algorithms */
  COLORCLEAR(gal.ambient_radiance);

  /* one iteration = gather to all patches */
  PatchListIterate(Patches, PatchGather);

  /* update the radiosity after gathering to all patches with Jacobi, immediately
   * update with Gauss-Seidel so the new radiosity are already used for the
   * still-to-be-processed patches in the same iteration. */
  if (gal.iteration_method == JACOBI)
    PatchListIterate(Patches, PatchUpdateRadiance);

  if (gal.importance_driven)
    PatchListIterate(Patches, PatchUpdatePotential);

  return FALSE;	/* never done, until we have a better criterium. */
}

/* #ident <<< WMP */
#ifdef WMP_WEIGHTS
void ElementCopyDirectRadiance(ELEMENT *top)
{
   if (top->regular_subelements)
   {
      ForAllRegularSubelements(child, top)
	 {
	    ElementCopyDirectRadiance(child);
	 }
      EndForAll;

   }
   
   COPYCOEFFICIENTS(top->direct_radiance, top->radiance, top->basis_size);
}


/* ********************** Weight transport ************************* */
static void ClusterGatherWeight(ELEMENT *rcv, INTERACTION *link, COLOR *srcrad)
{
  COLOR *rcvrad = rcv->F;
  int alpha, a;
  double area_factor = 1.0;
  
  if(gal.clustering_strategy != ISOTROPIC)
  {
    Error("ClusterGatherWeight", "Wrong clustering strategy");
    return;
  }

  a = MIN(link->nrcv, rcv->basis_size);
  
  for (alpha=0; alpha<a; alpha++)
    COLORADDSCALED(rcvrad[alpha], area_factor*link->FK[alpha*link->nsrc], srcrad[0], rcvrad[alpha]);
}


void ComputeWeightOverLink(INTERACTION *link)
{
  int a, alpha;
  COLOR srcrad2;
  
  if((!IsLightSource(link->src)) || IsCluster(link->src)
     || gal.iteration_nr > 1)
  {
    /* No transport calculated */
    return;
  }
  
  if (gal.iteration_method != JACOBI)
  {
    Error("ComputeWeightOverLink", "Jacobi iteration required");
    return;
  }
  
  /* We need Le of the source !! */
  srcrad2 = SELFEMITTED_RADIANCE(link->src->pog.patch);
  COLORPROD(srcrad2, srcrad2, srcrad2);  /* square */

  /* number of basis function on receiver for this link */
  a = MIN(link->nrcv, link->rcv->basis_size);   

  for (alpha=0; alpha<a; alpha++) 
  {
    /* Multiply by scrrad^2 !! This is ok, since only a constant
       approximation of emitted light is made */
    COLORADDSCALED(link->rcv->F[alpha], link->FK[alpha], srcrad2, link->rcv->F[alpha]);
  }
}

void GatheringComputeWeights(ELEMENT *top)
{
  COLOR tmpCoeff[MAXBASISSIZE];
  
  CLEARCOEFFICIENTS(tmpCoeff, top->basis_size);

  ForAllInteractions(link, top->interactions)
    {
      ComputeWeightOverLink(link);
    }
  EndForAll;
  
  
  ForAllRegularSubelements(child, top)
    {
      Push(top, top->F, child, child->F);
      
      GatheringComputeWeights(child);
    }
  EndForAll;

  ForAllIrregularSubelements(child, top)
    {
      Push(top, top->F, child, child->F);
      
      GatheringComputeWeights(child);
    }
  EndForAll;

}
#endif
/* #ident >>> WMP */


/* wat als ge clustering aan of afzet tijdens de berekeningen ??? */
int DoClusteredGatheringIteration(void)
{
  static double user_error_threshold, current_error_threshold;

  if (gal.importance_driven) {
    if (gal.iteration_nr<=1 || Camera.changed) {
      UpdateDirectPotential();
      ForAllPatches(P, Patches) {
	ELEMENT *top = TOPLEVEL_ELEMENT(P);
	float potential_increment = P->direct_potential - top->direct_potential.f;
	GatheringUpdateDirectPotential(top, potential_increment);
      } EndForAll;
      ClusterUpdatePotential(gal.top_cluster);
      Camera.changed = FALSE;
    }
  }

  printf("Gal iteration %i\n", gal.iteration_nr);

  /* initial linking stage is replaced by the creation of a self-link between
   * the whole scene and itself */
  if (gal.iteration_nr <= 1) {
    CreateInitialLinkWithTopCluster(gal.top_cluster, RECEIVER);
    current_error_threshold = 0.01;
  } else
    /* sliding error threshold */
    current_error_threshold /= 1.4;

  user_error_threshold = gal.rel_link_error_threshold;
  /*
  if (current_error_threshold < user_error_threshold)
    current_error_threshold = user_error_threshold;
  gal.rel_link_error_threshold = current_error_threshold; 
  */
  /* refines and computes light transport over the refined links */
  RefineInteractions(gal.top_cluster);

  gal.rel_link_error_threshold = user_error_threshold;

  /* push received radiance down the hierarchy to the leaf elements, where
   * it is multiplied with the reflectivity and the selfemitted radiance added,
   * and finally pulls back up for a consistent multiresolution representation
   * of radiance over all levels. */
  PushPullRadiance(gal.top_cluster);

  /* #ident <<< WMP */
#ifdef WMP_WEIGHTS
  if(gal.iteration_nr == 1)
  {
     /* First iteration done, now we can copy the direct radiance.
	Direct radiance is only copied for leaf elements.
	Later on, direct radiance is only pushed down the hierarchy
	to finer elements
     */

     ForAllPatches(p, Patches) {
	   ElementCopyDirectRadiance(TOPLEVEL_ELEMENT(p));
     } EndForAll;

     /* Now compute the weights in the hierarchy */

     CLEARCOEFFICIENTS(gal.top_cluster->F, gal.top_cluster->basis_size);
     GatheringComputeWeights(gal.top_cluster);
  }

#endif     
  /* #ident >>> WMP */

  if (gal.importance_driven)
    GatheringPushPullPotential(gal.top_cluster, 0.);

  /* no visualisation with ambient term for gathering radiosity algorithms */
  COLORCLEAR(gal.ambient_radiance);

  /* update the display colors of the patches */
  PatchListIterate(Patches, PatchRecomputeColor);

  return FALSE; /* never done */
}

/* updates the radiance on a surface element due to a change of the material 
 * properties */
static void ElementUpdateMaterial(ELEMENT *top,
				  COLOR oldrho, COLOR newrho, COLOR olde, COLOR newe)
{
  if (top->regular_subelements) {
    CLEARCOEFFICIENTS(top->radiance, top->basis_size);
    ForAllRegularSubelements(child, top) {
      COLOR Btmp[MAXBASISSIZE];
      ElementUpdateMaterial(child, oldrho, newrho, olde, newe);
      Pull(top, Btmp, child, child->radiance);
      ADDCOEFFICIENTS(top->radiance, Btmp, top->basis_size);
    } EndForAll;
  } else {
    int i;

    COLORSUBTRACT(top->radiance[0], olde, top->radiance[0]);
    for (i=0; i<top->basis_size; i++) {
      COLORDIV(top->radiance[i], oldrho, top->radiance[i]);
      COLORPROD(top->radiance[i], newrho, top->radiance[i]);
    }
    COLORADD(top->radiance[0], newe, top->radiance[0]);
  }
}

static void ClusterUpdateMaterial(ELEMENT *clus)
{
  if (clus->irregular_subelements) {
    CLEARCOEFFICIENTS(clus->radiance, clus->basis_size);
    ForAllIrregularSubelements(child, clus) {
      COLOR Btmp[MAXBASISSIZE];
      ClusterUpdateMaterial(child);
      Pull(clus, Btmp, child, child->radiance);
      ADDCOEFFICIENTS(clus->radiance, Btmp, clus->basis_size);
    } EndForAll;
  }
}

void GatheringUpdateMaterial(MATERIAL *oldmaterial, MATERIAL *newmaterial)
{
  /*TODO*/
  Warning("GatheringUpdateMaterial", "Out of order");
#ifdef NEVER
  POINT p = {0.,0.,0.};
  COLOR newrho = BsdfDiffuseReflectance(newmaterial->bsdf, &p),
        oldrho = BsdfDiffuseReflectance(oldmaterial->bsdf, &p),
        newe = EdfDiffuseRadiance(newmaterial->edf, &p),
        olde = EdfDiffuseRadiance(oldmaterial->edf, &p);
  {
    static int wgiv=0;
    if (!wgiv) {
      Warning("GatheringUpdateMaterial", "will only work for non-textured materials");
      wgiv = 1;
    }
  }
  ForAllPatches(p, Patches) {
    if (p->surface->material == oldmaterial)
      ElementUpdateMaterial(TOPLEVEL_ELEMENT(p), oldrho, newrho, olde, newe);
  } EndForAll;
  ClusterUpdateMaterial(gal.top_cluster);
#endif
}