/* galerkin.h: Galerkin radiosity, with or without hierarchical refinement, with or
 *	       without clusters, with Jacobi, Gauss-Seidel or Southwell iterations,
 *	       potential-dirven or not ...
 *
 * References:
 *
 */

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

#include "galerkin.h"
#include "galerkinP.h"
#include "scene.h"
#include "statistics.h"
#include "vertex.h"
#include "error.h"
#include "pools.h"
#include "basis.h"
#include "canvas.h"
#include "render.h"
#include "camera.h"
#include "topology.h"
#include "scene.h"
#include "cluster.h"
#include "monitor.h"
#include "scratch.h"
#include "options.h"
#include "writevrml.h"
#include "options.h"
#include "coefficients.h"
#include "mystrings.h"

GALERKIN_STATE gal;

/* installs cubature rules for triangles and quadrilaterals of the specified degree */
void SetCubatureRules(CUBARULE **trirule, CUBARULE **quadrule, CUBATURE_DEGREE degree)
{
  switch (degree) {
  case DEGREE_1:
    *trirule = &CRT1;
    *quadrule = &CRQ1;
    break;
  case DEGREE_2:
    *trirule = &CRT2;
    *quadrule = &CRQ2;
    break;
  case DEGREE_3:
    *trirule = &CRT3;
    *quadrule = &CRQ3;
    break;
  case DEGREE_4:
    *trirule = &CRT4;
    *quadrule = &CRQ4;
    break;
  case DEGREE_5:
    *trirule = &CRT5;
    *quadrule = &CRQ5;
    break;
  case DEGREE_6:
    *trirule = &CRT7;
    *quadrule = &CRQ6;
    break;
  case DEGREE_7:
    *trirule = &CRT7;
    *quadrule = &CRQ7;
    break;
  case DEGREE_8:
    *trirule = &CRT8;
    *quadrule = &CRQ8;
    break;
  case DEGREE_9:
    *trirule = &CRT9;
    *quadrule = &CRQ9;
    break;
  case DEGREE_3_PROD:
    *trirule = &CRT5;
    *quadrule = &CRQ3PG;
    break;
  case DEGREE_5_PROD:
    *trirule = &CRT7;
    *quadrule = &CRQ5PG;
    break;
  case DEGREE_7_PROD:
    *trirule = &CRT9;
    *quadrule = &CRQ7PG;
    break;
  default:
    Fatal(2, "SetCubatureRules", "Invalid degree %d", degree);
  }
}

static void GalerkinDefaults(void)
{
  gal.hierarchical = DEFAULT_GAL_HIERARCHICAL;
  gal.importance_driven = DEFAULT_GAL_IMPORTANCE_DRIVEN;
  gal.clustered = DEFAULT_GAL_CLUSTERED;
  gal.iteration_method = DEFAULT_GAL_ITERATION_METHOD;
  gal.lazy_linking = DEFAULT_GAL_LAZY_LINKING;
  gal.use_constant_radiance = DEFAULT_GAL_CONSTANT_RADIANCE;
  gal.use_ambient_radiance = DEFAULT_GAL_AMBIENT_RADIANCE;
  gal.shaftcullmode = DEFAULT_GAL_SHAFTCULLMODE;
  gal.rcv_degree = DEFAULT_GAL_RCV_CUBATURE_DEGREE;
  gal.src_degree = DEFAULT_GAL_SRC_CUBATURE_DEGREE;
  SetCubatureRules(&gal.rcv3rule, &gal.rcv4rule, gal.rcv_degree);
  SetCubatureRules(&gal.src3rule, &gal.src4rule, gal.src_degree);
  gal.clusRule = &CRV1;
  gal.rel_min_elem_area = DEFAULT_GAL_REL_MIN_ELEM_AREA;
  gal.rel_link_error_threshold = DEFAULT_GAL_REL_LINK_ERROR_THRESHOLD;
  gal.error_norm = DEFAULT_GAL_ERROR_NORM;
  gal.basis_type = DEFAULT_GAL_BASIS_TYPE;
  gal.exact_visibility = DEFAULT_GAL_EXACT_VISIBILITY;
  gal.multires_visibility = DEFAULT_GAL_MULTIRES_VISIBILITY;
  gal.rel_total_error_threshold = DEFAULT_GAL_REL_TOTAL_ERROR_THRESHOLD;
  gal.clustering_strategy = DEFAULT_GAL_CLUSTERING_STRATEGY;
  gal.scratch = (SGL_CONTEXT *)NULL;
  gal.scratch_fb_size = DEFAULT_GAL_SCRATCH_FB_SIZE;

  gal.iteration_nr = -1;	/* means "not initialized" */
}

static int true = TRUE, false = FALSE;

static void IterationMethodOption(void *value)
{
  char *name = *(char **)value;

  if (strncasecmp(name, "jacobi", 2) == 0)
    gal.iteration_method = JACOBI;
  else if (strncasecmp(name, "gaussseidel", 2) == 0)
    gal.iteration_method = GAUSS_SEIDEL;
  else if (strncasecmp(name, "southwell", 2) == 0)
    gal.iteration_method = SOUTHWELL;
  else
    Error(NULL, "Invalid iteration method '%s'", name);
}

static void HierarchicalOption(void *value)
{
  int yesno = *(int *)value;
  gal.hierarchical = yesno;
}

static void LazyOption(void *value)
{
  int yesno = *(int *)value;
  gal.lazy_linking = yesno;
}

static void ClusteringOption(void *value)
{
  int yesno = *(int *)value;
  gal.clustered = yesno;
}

static void ImportanceOption(void *value)
{
  int yesno = *(int *)value;
  gal.importance_driven = yesno;
}

static void AmbientOption(void *value)
{
  int yesno = *(int *)value;
  gal.use_ambient_radiance = yesno;
}

static CMDLINEOPTDESC galerkinOptions[] = {
  {"-gr-iteration-method", 6,	Tstring,	NULL,	IterationMethodOption,
   "-gr-iteration-method <methodname>: Jacobi, GaussSeidel, Southwell"},
  {"-gr-hierarchical",	6,	TYPELESS,	(void *)&true, HierarchicalOption,
   "-gr-hierarchical    \t: do hierarchical refinement"},
  {"-gr-not-hierarchical", 10,	TYPELESS,	(void *)&false, HierarchicalOption,
   "-gr-not-hierarchical\t: don't do hierarchical refinement"},
  {"-gr-lazy-linking",	6,	TYPELESS,	(void *)&true,	LazyOption,
   "-gr-lazy-linking    \t: do lazy linking"},
  {"-gr-no-lazy-linking", 10,	TYPELESS,	(void *)&false,	LazyOption,
   "-gr-no-lazy-linking \t: don't do lazy linking"},
  {"-gr-clustering",	6,	TYPELESS,	(void *)&true,	ClusteringOption,
   "-gr-clustering      \t: do clustering"},
  {"-gr-no-clustering", 10,	TYPELESS,	(void *)&false,	ClusteringOption,
   "-gr-no-clustering   \t: don't do clustering"},
  {"-gr-importance",	6,	TYPELESS,	(void *)&true,	ImportanceOption,
   "-gr-importance      \t: do view-potential driven computations"},
  {"-gr-no-importance", 10,	TYPELESS,	(void *)&false,	ImportanceOption,
   "-gr-no-importance   \t: don't use view-potential"},
  {"-gr-ambient",	6,	TYPELESS,	(void *)&true,	AmbientOption,
   "-gr-ambient         \t: do visualisation with ambient term"},
  {"-gr-no-ambient", 10,	TYPELESS,	(void *)&false,	AmbientOption,
   "-gr-no-ambient      \t: do visualisation without ambient term"},
  {"-gr-link-error-threshold", 6, Tfloat,	&gal.rel_link_error_threshold, DEFAULT_ACTION,
   "-gr-link-error-threshold <float>: Relative link error threshold"},
  {"-gr-min-elem-area",	6,	Tfloat,		&gal.rel_min_elem_area, DEFAULT_ACTION,
   "-gr-min-elem-area <float> \t: Relative element area threshold"},
  {NULL,		0,	TYPELESS,	NULL,	DEFAULT_ACTION,
   NULL}
};

static void ParseGalerkinOptions(int *argc, char **argv)
{
  ParseOptions(galerkinOptions, argc, argv);
}

static void PrintGalerkinOptions(FILE *fp)
{
  fprintf(fp, "\nGalerkin radiosity options:\n");
  PrintOptions(fp, galerkinOptions);
}

/* for counting how much CPU time was used for the computations */
static void UpdateCpuSecs(void)
{
  clock_t t;

  t = clock();
  gal.cpu_secs += (float)(t - gal.lastclock)/(float)CLOCKS_PER_SEC;
  gal.lastclock = t;
}

/* for waking up now and then while the computations are going on */
static void wake_up(int sig)
{
  gal.wake_up = TRUE;
  signal(SIGALRM, wake_up);
  alarm(/*after*/ 1 /*second*/);

  UpdateCpuSecs();
}

/* radiance data for a PATCH is a surface element. */
static void *CreatePatchData(PATCH *patch)
{
  return patch->radiance_data = (void *)CreateToplevelElement(patch);
}

static void PrintPatchData(FILE *out, PATCH *patch)
{
  PrintElement(out, (ELEMENT *)patch->radiance_data);
}

static void DestroyPatchData(PATCH *patch)
{
  DestroyToplevelElement((ELEMENT *)patch->radiance_data);
  patch->radiance_data = (void *)NULL;
}

void PatchRecomputeColor(PATCH *patch)
{
  COLOR rho = REFLECTIVITY(patch);
  COLOR rad_vis;

  /* compute the patches color based on its radiance + ambient radiance
   * if desired. */
  if (gal.use_ambient_radiance) {
    COLORPROD(rho, gal.ambient_radiance, rad_vis);
    COLORADD(rad_vis, RADIANCE(patch), rad_vis);
    RadianceToRGB(rad_vis, &patch->color);
  } else {
    RadianceToRGB(RADIANCE(patch), &patch->color);
  }
  PatchComputeVertexColors(patch);
}

static void PatchInit(PATCH *patch)
{
  COLOR rho = REFLECTIVITY(patch), Ed = SELFEMITTED_RADIANCE(patch);

  if (gal.use_constant_radiance) {
    /* see Neumann et al, "The Constant Radiosity Step", Eurographics Rendering Workshop
     * '95, Dublin, Ireland, June 1995, p 336-344. */
    COLORPROD(rho, gal.constant_radiance, RADIANCE(patch));
    COLORADD(RADIANCE(patch), Ed, RADIANCE(patch));
    if (gal.iteration_method == SOUTHWELL)
      COLORSUBTRACT(RADIANCE(patch), gal.constant_radiance, UNSHOT_RADIANCE(patch));
  } else {
    RADIANCE(patch) = Ed;
    if (gal.iteration_method == SOUTHWELL)
      UNSHOT_RADIANCE(patch) = RADIANCE(patch);
  }

  if (gal.importance_driven) {
    switch (gal.iteration_method) {
    case GAUSS_SEIDEL:
    case JACOBI:
      POTENTIAL(patch).f = patch->direct_potential;
      break;
    case SOUTHWELL:
      POTENTIAL(patch).f = UNSHOT_POTENTIAL(patch).f = patch->direct_potential;
      break;
    default:
      Fatal(-1, "PatchInit", "Invalid iteration method");
    }
  }

  PatchRecomputeColor(patch);
}

void InitGalerkin(void)
{
  gal.iteration_nr = gal.step_nr = 0;
  gal.cpu_secs = 0.;

  InitBasis();

  gal.constant_radiance = estimated_average_radiance;
  if (gal.use_constant_radiance) {
    COLORCLEAR(gal.ambient_radiance);
  } else
    gal.ambient_radiance = estimated_average_radiance;

  PatchListIterate(Patches, PatchInit);

  gal.top_geom = ClusteredWorldGeom;
  gal.top_cluster = GalerkinCreateClusterHierarchy(gal.top_geom);

  /* create a scratch software renderer for various operations on clusters */
  ScratchInit();

  /* global variables used for formfactor computation optimisation */
  gal.fflastrcv = gal.fflastsrc = (ELEMENT *)NULL;

  /* global variables for scratch rendering */
  gal.lastclusid = -1; VECTORSET(gal.lasteye, HUGE, HUGE, HUGE);
}

static int DoGalerkinOneStep(void)
{
  int done = FALSE;
  void (*prev_alrm_handler)(int signr); 
  unsigned prev_alarm_left;

  if (gal.iteration_nr < 0) {
    Error("DoGalerkinOneStep", "method not initialized");
    return TRUE;	/* done, don't continue! */
  }

  /* install a timer that will wake us up ofter one second for checking for 
   * user events */
  prev_alrm_handler = signal(SIGALRM, wake_up);
  prev_alarm_left = alarm(/*after*/ 1 /*second*/);
  gal.wake_up = FALSE;

  gal.iteration_nr++;
  gal.lastclock = clock();

  /* and now the real work */
  switch (gal.iteration_method) {
  case JACOBI:
  case GAUSS_SEIDEL:
    if (gal.clustered)
      done = DoClusteredGatheringIteration();
    else
      done = DoGatheringIteration();
    break;
  case SOUTHWELL:
    done = DoShootingStep();
    break;
  default:
    Fatal(2, "DoGalerkinOneStep", "Invalid iteration method %d\n", gal.iteration_method);
  }

  UpdateCpuSecs();

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

  return done;
}

static void GalerkinUpdatePatchPotential(PATCH *patch, float potential_increment)
{
  if (gal.iteration_method == SOUTHWELL)
    ShootingUpdateDirectPotential(TOPLEVEL_ELEMENT(patch), potential_increment);
  else if (gal.iteration_method == JACOBI ||
	   gal.iteration_method == GAUSS_SEIDEL)
    GatheringUpdateDirectPotential(TOPLEVEL_ELEMENT(patch), potential_increment);
}

static void TerminateGalerkin(void)
{
  ScratchTerminate();
  GalerkinDestroyClusterHierarchy(gal.top_cluster);
}

static COLOR GetRadiance(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *leaf;
  COLOR rad;

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  leaf = RegularLeafElementAtPoint(TOPLEVEL_ELEMENT(patch), &u, &v);

  rad = RadianceAtPoint(leaf, leaf->radiance, u, v);

  if (gal.use_ambient_radiance) {
    /* add ambient radiance */
    COLOR rho = REFLECTIVITY(patch);
    COLOR ambirad;
    COLORPROD(rho, gal.ambient_radiance, ambirad);
    COLORADD(rad, ambirad, rad);
  }

  return rad;
}

/* #ident <<< WMP */
#ifdef WMP_WEIGHTS
COLOR GAL_GetAccuracy(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *leaf;
  COLOR gradient, rad[4], grad[6], col;
  double avgRad;

  COLORSETMONOCHROME(col, 1.0);

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  leaf = RegularLeafElementAtPoint(TOPLEVEL_ELEMENT(patch), &u, &v);

  /* Take the radiance at the corner points */

  if(patch->nrvertices == 3)
  {
    rad[0] = RadianceAtPoint(leaf, leaf->direct_radiance, 0, 0);
    rad[1] = RadianceAtPoint(leaf, leaf->direct_radiance, 1, 0);
    rad[2] = RadianceAtPoint(leaf, leaf->direct_radiance, 0, 1);

    COLORSUBTRACT(rad[1], rad[0], grad[0]);
    COLORSUBTRACT(rad[2], rad[0], grad[1]);
    COLORSUBTRACT(rad[2], rad[1], grad[2]);

    COLORABS(grad[0], grad[0]);
    COLORABS(grad[1], grad[1]);
    COLORABS(grad[2], grad[2]);

    /* Keep largest gradient for different spectral components */
    COLORMAX(grad[0], grad[1], gradient);
    COLORMAX(grad[2], gradient, gradient);
  }
  else
  {
    /* 4 vertices */
    rad[0] = RadianceAtPoint(leaf, leaf->direct_radiance, 0, 0);
    rad[1] = RadianceAtPoint(leaf, leaf->direct_radiance, 1, 0);
    rad[2] = RadianceAtPoint(leaf, leaf->direct_radiance, 0, 1);
    rad[3] = RadianceAtPoint(leaf, leaf->direct_radiance, 1, 1);

    COLORSUBTRACT(rad[1], rad[0], grad[0]);
    COLORSUBTRACT(rad[2], rad[0], grad[1]);
    COLORSUBTRACT(rad[3], rad[0], grad[2]);

    COLORSUBTRACT(rad[2], rad[1], grad[3]);
    COLORSUBTRACT(rad[3], rad[1], grad[4]);
    COLORSUBTRACT(rad[3], rad[2], grad[5]);

    COLORABS(grad[0], grad[0]);
    COLORABS(grad[1], grad[1]);
    COLORABS(grad[2], grad[2]);
    COLORABS(grad[3], grad[3]);
    COLORABS(grad[4], grad[4]);
    COLORABS(grad[5], grad[5]);

    /* Keep largest gradient for different spectral components */
    COLORMAX(grad[0], grad[1], grad[0]);
    COLORMAX(grad[2], grad[3], grad[2]);
    COLORMAX(grad[4], grad[5], grad[4]);
    COLORMAX(grad[0], grad[2], gradient);
    COLORMAX(grad[4], gradient, gradient);
  }

  /* We want relative gradients */

  avgRad = COLORAVERAGE(leaf->direct_radiance[0]);

  if(avgRad < 0.1)
  {
    return col;
  }

  COLORSCALEINVERSE(avgRad, gradient, gradient);

/*  printf("Col %g\n", COLORMAXCOMPONENT(gradient)); */

  COLORSCALE(1.0, gradient, gradient);

  COLORSUBTRACT(col, gradient, gradient);
  COLORCLIPPOSITIVE(gradient, gradient);

  return gradient;
}

COLOR GAL_GetSeWeight(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *leaf;
  COLOR weight, rho2;

  rho2 = REFLECTIVITY(patch);
  COLORPROD(rho2, rho2,rho2);

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  leaf = RegularLeafElementAtPoint(TOPLEVEL_ELEMENT(patch), &u, &v);

  weight = RadianceAtPoint(leaf, leaf->F, u, v);

  COLORPROD(weight, rho2, weight);
  
  return weight;
}

COLOR GAL_GetIndirectRadiance(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *leaf;
  COLOR rad, drad;

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  leaf = RegularLeafElementAtPoint(TOPLEVEL_ELEMENT(patch), &u, &v);

  rad = RadianceAtPoint(leaf, leaf->radiance, u, v);
  drad = RadianceAtPoint(leaf, leaf->direct_radiance, u, v);

  COLORSUBTRACT(rad, drad, rad);

  return rad;
}


COLOR GAL_GetDirectRadiance(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *leaf;
  COLOR rad, emitted;
  EDF *edf;

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  leaf = RegularLeafElementAtPoint(TOPLEVEL_ELEMENT(patch), &u, &v);

  rad = RadianceAtPoint(leaf, leaf->direct_radiance, u, v);

  /* Subtract the emitted radiance */
  
  edf = leaf->pog.patch->surface->material->edf;

  if(edf)
  {
     emitted = EdfEmittance(edf, NULL, DIFFUSE_COMPONENT);
     COLORSCALEINVERSE(M_PI, emitted, emitted);
     COLORSUBTRACT(rad, emitted, rad);
  }  

  return rad;
}

/* returns direct radiance: assumes gatherin without clustering */
static COLOR OldGetDirectRadiance(PATCH *patch, double u, double v, VECTOR dir)
{
  ELEMENT *elem;
  COLOR rad;
  COLOR rho = REFLECTIVITY(patch);

  COLORCLEAR(rad);

  if (!(gal.iteration_method == JACOBI || gal.iteration_method == GAUSS_SEIDEL) || gal.clustered) {
    static int wgiv = 0;
    if (!wgiv) {
      Error("GetDirectRadiance", "only implemented for gathering without clustering");
      wgiv = 1;
      return rad;
    }
  }

  if (patch->jacobian)
    BilinearToUniform(patch, &u, &v);

  elem = TOPLEVEL_ELEMENT(patch);
  while (elem) {
    int nrlightsrc = 0;
    COLOR rcvdirectrad[MAXBASISSIZE];
    CLEARCOEFFICIENTS(rcvdirectrad, elem->basis_size);

    ForAllInteractions(link, elem->interactions) {
      ELEMENT *src = link->src;	 /* assume gathering without clusering */
      if (!IsCluster(src) && src->pog.patch->surface->material->edf) {
	/* src is light source: gather direct illumination */
	COLOR Ed = SELFEMITTED_RADIANCE(src->pog.patch);
	if (link->nsrc==1 && link->nrcv==1) {
	  COLORADDSCALED(rcvdirectrad[0], link->K.f, Ed, rcvdirectrad[0]);
	} else {
	  int a, alpha;
	  a = MIN(link->nrcv, link->rcv->basis_size);
	  for (alpha=0; alpha<a; alpha++) {
	    COLORADDSCALED(rcvdirectrad[alpha], link->K.p[alpha*link->nsrc], Ed, rcvdirectrad[alpha]);
	  }
	}
	nrlightsrc++;
      }
    } EndForAll;

    if (nrlightsrc > 0) {
      COLOR directrad = RadianceAtPoint(elem, rcvdirectrad, u, v);
      COLORADDSCALED(rad, 1./elem->area, directrad, rad);
    }

    if (elem->regular_subelements) {
      elem = RegularSubelementAtPoint(elem, &u, &v);
    } else
      elem = (ELEMENT *)NULL;	/* break out of the loop */
  }

  COLORPROD(rho, rad, rad);
  return rad;
}
#endif /* WMP_WEIGHTS */
/* #ident >>> WMP */

static char *GetGalerkinStats(void)
{
  static char stats[2000];
  char *p;
  int n;

  p = stats;
  sprintf(p, "Galerkin Radiosity Statistics:\n\n%n", &n); p += n;
  sprintf(p, "Iteration: %d\n\n%n", gal.iteration_nr, &n); p += n;
  sprintf(p, "Nr. elements: %d\n%n", GetNumberOfElements(), &n); p += n;
  sprintf(p, "clusters: %d\n%n", GetNumberOfClusters(), &n); p += n;
  sprintf(p, "surface elements: %d\n\n%n", GetNumberOfSurfaceElements(), &n); p += n;
  sprintf(p, "Nr. interactions: %d\n%n", GetNumberOfInteractions(), &n); p += n;
  sprintf(p, "cluster to cluster: %d\n%n", GetNumberOfClusterToClusterInteractions(), &n); p += n;  
  sprintf(p, "cluster to surface: %d\n%n", GetNumberOfClusterToSurfaceInteractions(), &n); p += n;  
  sprintf(p, "surface to cluster: %d\n%n", GetNumberOfSurfaceToClusterInteractions(), &n); p += n;  
  sprintf(p, "surface to surface: %d\n\n%n", GetNumberOfSurfaceToSurfaceInteractions(), &n); p += n;  
  sprintf(p, "CPU time: %g secs.\n%n", gal.cpu_secs, &n); p += n;
  sprintf(p, "Memory usage: %ld KBytes.\n\n%n", GetMemoryUsage()/1024L, &n); p += n;
  sprintf(p, "Minimum element area: %g m^2\n%n", total_area * (double)gal.rel_min_elem_area, &n); p += n;
  sprintf(p, "Link error threshold: %g %s\n\n%n", 
	  (double)(gal.error_norm == RADIANCE_ERROR ?
	   M_PI * (gal.rel_link_error_threshold * ColorLuminance(max_selfemitted_radiance)) :
	   gal.rel_link_error_threshold * ColorLuminance(max_selfemitted_power)),
	  (gal.error_norm == RADIANCE_ERROR ? "lux" : "lumen"), 
	  &n); p += n;

  return stats;
}

static void RenderElementHierarchy(ELEMENT *elem)
{
  if (!elem->regular_subelements)
    RenderElement(elem);
  else
    ITERATE_REGULAR_SUBELEMENTS(elem, RenderElementHierarchy);
}

static void GalerkinRenderPatch(PATCH *patch)
{
  RenderElementHierarchy(TOPLEVEL_ELEMENT(patch));
}

void GalerkinRender(void)
{
  if (renderopts.frustum_culling)
    RenderWorldOctree(GalerkinRenderPatch);
  else
    PatchListIterate(Patches, GalerkinRenderPatch);
}

#include "vertex.h"

static void GalerkinRenderPatchPotential(PATCH *patch)
{
  POINT v[4];
  int i;
  RGB color;

  /*  RenderElementHierarchy(TOPLEVEL_ELEMENT(patch)); */
  float gray = TOPLEVEL_ELEMENT(patch)->potential.f;
  if (gray > 1.) gray = 1.;
  if (gray < 0.) gray = 0.;

  RGBSET(color, gray, gray, gray);
  RenderSetColor(&color);

  for (i=0; i<patch->nrvertices; i++)
    v[i] = *patch->vertex[i]->point;

  RenderPolygonFlat(patch->nrvertices, v);

  RenderSetColor(&renderopts.outline_color);
  for (i=0; i<patch->nrvertices; i++) {
    RenderLine(&v[i], &v[(i+1)%patch->nrvertices]);
  }
}

void GalerkinRenderPotential(void)
{
  PatchListIterate(Patches, GalerkinRenderPatchPotential); 
}

static void BuildIdToPatchTable(PATCH *patch, PATCH **id2patch)
{
  id2patch[patch->id] = patch;
}

void GalerkinRaycast(FILE *fp)
{
  PATCH **id2patch;
  unsigned long *ids, *id;
  long i, j, x, y, nrpixels, maxpatchid, lostpixels;
  RGB *pix;
  RAY ray;
  float v, h, xsample, ysample, dist;
  HITREC *hit;

  CanvasPushMode(CANVASMODE_RENDER);

  /* get the patch IDs for each pixel. */
  ids = RenderIds(&x, &y);
  if (!ids) return;
  nrpixels = x*y; lostpixels=0;

  /* build a table to convert a patch ID to the corresponding PATCH * */
  maxpatchid = PatchGetNextID()-1;  
  id2patch = (PATCH **)Alloc((int)(maxpatchid+1) * sizeof(PATCH *));
  for (i=0; i<=maxpatchid; i++) id2patch[i] = (PATCH *)NULL;
  PatchListIterate1A(Patches, BuildIdToPatchTable, (void *)id2patch);

  pix = (RGB *)Alloc((int)x * sizeof(RGB));

  /* h and v are the horizontal resp. vertical distance between two 
   * neighbooring pixels on the screen. */
  h = 2. * tan(Camera.hfov * M_PI/180.) / (float)x;
  v = 2. * tan(Camera.vfov * M_PI/180.) / (float)y;

  /* the ray origin is the same for all eye-rays */
  ray.pos = Camera.eyep; 

  if (fp) 
    fprintf(fp, "P6\n%ld %ld\n255\n", x, y);

  for (j=y-1, ysample=-v*(float)(y-1)/2.; j>=0; j--, ysample+=v) {
    id = ids + j*x;
    for (i=0, xsample=-h*(float)(x-1)/2.; i<x; i++, id++, xsample+=h) {
      unsigned long the_id = (*id)&0xffffff;
      PATCH *the_patch = the_id <= maxpatchid ? id2patch[the_id] : (PATCH *)NULL;
      HITREC myhit;

      if (the_id == 0) {	/* background visible through pixel */
	pix[i] = Black;
	continue;
      }

      /* construct ray direction to center of pixel */
      VECTORCOMB3(Camera.Z, xsample, Camera.X, ysample, Camera.Y, ray.dir);
      VECTORNORMALIZE(ray.dir);

      pix[i] = Black;
      dist = HUGE;
      hit = (HITREC *)NULL;
      if (the_patch)		/* most probably the nearest patch */
      {
	hit = PatchIntersect(the_patch, &ray, 0., &dist, HIT_FRONT, &myhit);
      }

      if (!hit) {
	dist = HUGE;
	hit = GeomListDiscretisationIntersect(ClusteredWorld, &ray, 0., &dist, HIT_FRONT, &myhit);
      }

      if (hit) {
	double u, v;
	VECTOR pt, dir;
	COLOR rad;

	VECTORADDSCALED(ray.pos, dist, ray.dir, pt);
	PatchUV(hit->patch, &pt, &u, &v);
	VECTORSCALE(-1, ray.dir, dir);

	rad = GetRadiance(hit->patch, u, v, dir);

	/* #ident <<< WMP */
	/* rad = GAL_GetAccuracy(hit->patch, u, v, dir); */
	/* rad = GAL_GetSeWeight(hit->patch, u, v, dir); */
	/* #ident >>> WMP */

	RadianceToRGB(rad, &pix[i]);
	/*
	if (hit->patch != the_patch) {
	  pix[i] = Red;
	  lostpixels++;
	}
	*/
      } else {
	pix[i] = Black;
	/*	
	if (the_patch) {
	  pix[i] = Red;
	  lostpixels++;
	}
	*/
      }
    }

    if (fp) {
      for (i=0; i<x; i++) {
	fprintf(fp, "%c%c%c", 
		(char)(pix[i].r * 255.), 
		(char)(pix[i].g * 255.), 
		(char)(pix[i].b * 255.));
      }
    }

    RenderPixels(0, j, x, 1, pix);
  }

  if (lostpixels > 0)
    Warning(NULL, "%d lost pixels", lostpixels);

  Free((char *)pix, (int)x * sizeof(RGB));
  Free((char *)id2patch, (int)(maxpatchid+1) * sizeof(PATCH *));
  Free((char *)ids, (int)nrpixels * sizeof(unsigned long));

  CanvasPullMode();
}

static void GalerkinUpdateMaterial(MATERIAL *oldmaterial, MATERIAL *newmaterial)
{
  if (gal.iteration_method == SOUTHWELL)
    ShootingUpdateMaterial(oldmaterial, newmaterial);
  else if (gal.iteration_method == JACOBI ||
	   gal.iteration_method == GAUSS_SEIDEL)
    GatheringUpdateMaterial(oldmaterial, newmaterial);
  else
    fprintf(stderr, "GalerkinUpdateMaterial: not yet implemented.\n");

  RenderScene();
}

/* *********************************************************** */
/* VRML output */
static FILE *vrmlfp;
static int nwrit, vid;

static void WriteVertexCoord(POINT *p)
{
  if (nwrit>0)
    fprintf(vrmlfp, ", ");
  nwrit++;
  if (nwrit%4 == 0)
    fprintf(vrmlfp, "\n\t  ");
  fprintf(vrmlfp, "%g %g %g", p->x, p->y, p->z);
  vid++;
}

static void WriteVertexCoords(ELEMENT *elem)
{
  POINT v[8];
  int i, nverts = ElementVertices(elem, v);
  for (i=0; i<nverts; i++) {
    WriteVertexCoord(&v[i]);
  }
}

static void WriteCoords(void)
{
  nwrit = vid = 0;
  fprintf(vrmlfp, "\tcoord Coordinate {\n\t  point [ ");
  ForAllLeafElements(gal.top_cluster, WriteVertexCoords);
  fprintf(vrmlfp, " ] ");
  fprintf(vrmlfp, "\n\t}\n");
}

static void WriteVertexColor(RGB *color)
{
  /* not yet written */
  if (nwrit>0)
    fprintf(vrmlfp, ", ");
  nwrit++;
  if (nwrit%4 == 0)
    fprintf(vrmlfp, "\n\t  ");
  fprintf(vrmlfp, "%.3g %.3g %.3g", color->r, color->g, color->b);
  vid++;
}

static void ElementWriteVertexColors(ELEMENT *element)
{
  COLOR vertrad[4];
  int i;

  if (element->pog.patch->nrvertices == 3) {
    vertrad[0] = RadianceAtPoint(element, element->radiance, 0., 0.);
    vertrad[1] = RadianceAtPoint(element, element->radiance, 1., 0.);
    vertrad[2] = RadianceAtPoint(element, element->radiance, 0., 1.);
  } else {
    vertrad[0] = RadianceAtPoint(element, element->radiance, 0., 0.);
    vertrad[1] = RadianceAtPoint(element, element->radiance, 1., 0.);
    vertrad[2] = RadianceAtPoint(element, element->radiance, 1., 1.);
    vertrad[3] = RadianceAtPoint(element, element->radiance, 0., 1.);
  }

  if (gal.use_ambient_radiance) {
    COLOR rho = REFLECTIVITY(element->pog.patch), ambient;
    
    COLORPROD(rho, gal.ambient_radiance, ambient);
    for (i=0; i<element->pog.patch->nrvertices; i++) {
      COLORADD(vertrad[i], ambient, vertrad[i]);
    }
  }

  for (i=0; i<element->pog.patch->nrvertices; i++) {
    RGB col;
    RadianceToRGB(vertrad[i], &col);
    WriteVertexColor(&col);
  }
}

static void WriteVertexColors(void)
{
  vid = nwrit = 0;
  fprintf(vrmlfp, "\tcolor Color {\n\t  color [ ");
  ForAllLeafElements(gal.top_cluster, ElementWriteVertexColors);
  fprintf(vrmlfp, " ] ");
  fprintf(vrmlfp, "\n\t}\n");
}

static void WriteColors(void)
{
  if (!renderopts.smooth_shading)
    Warning(NULL, "I assume you want a smooth shaded model ...");
  fprintf(vrmlfp, "\tcolorPerVertex %s\n", "TRUE");
  WriteVertexColors();
}

static void WriteCoordIndex(int index)
{
  nwrit++;
  if (nwrit%20 == 0) fprintf(vrmlfp, "\n\t  ");
  fprintf(vrmlfp, "%d ", index);
}

static void ElementWriteCoordIndices(ELEMENT *elem)
{
  int i;
  for (i=0; i<elem->pog.patch->nrvertices; i++)
    WriteCoordIndex(vid++);
  WriteCoordIndex(-1);
}

static void WriteCoordIndices(void)
{
  vid = nwrit = 0;
  fprintf(vrmlfp, "\tcoordIndex [ ");
  ForAllLeafElements(gal.top_cluster, ElementWriteCoordIndices);
  fprintf(vrmlfp, " ]\n"); 
}

static void GalerkinWriteVRML(FILE *fp)
{
  WriteVRMLHeader(fp);

  vrmlfp = fp;
  WriteCoords();
  WriteColors();
  WriteCoordIndices();

  WriteVRMLTrailer(fp);
}

/* **************************************************************** */
/* MGF mesh output: vertices are not shared!! */
static FILE *mgfp;
static int mgfirstface;

static void IteratePrimitiveGeoms(GEOMLIST *geomlist, void (*func)(GEOM *))
{
  ForAllGeoms(geom, geomlist) {
    GEOMLIST *primlist = GeomPrimList(geom);
    if (!primlist)
      func(geom);
    else {
      IteratePrimitiveGeoms(primlist, func);
    }
  } EndForAll;
}

static void MgfWriteVertex(POINT *p)
{
  fprintf(mgfp, "v v%d =\n\tp %g %g %g\n", 
	  vid, p->x, p->y, p->z);
  vid++;
}

static void MgfWriteLeafElement(ELEMENT *elem)
{
  POINT v[8];
  int i, firstvid=vid;

  int nverts = ElementVertices(elem, v);
  for (i=0; i<nverts; i++) {
    MgfWriteVertex(&v[i]);
  }

  fprintf(mgfp, "f ");
  for (i=0; i<nverts; i++) {
    fprintf(mgfp, "v%d ", firstvid+i);
  }
  fprintf(mgfp, "\n");
}

static void MgfWritePatch(PATCH *patch)
{
  if (mgfirstface)
    fprintf(mgfp, "m %s\n", patch->surface->material->name);
  mgfirstface=FALSE;
  ForAllLeafElements(TOPLEVEL_ELEMENT(patch), MgfWriteLeafElement);
}

static void MgfWriteGeom(GEOM *geom)
{
  mgfirstface=TRUE;
  PatchListIterate(GeomPatchList(geom), MgfWritePatch);
}

void GalerkinWriteMGF(FILE *fp)
{
  mgfp = fp;
  vid = 1;
  IteratePrimitiveGeoms(World, MgfWriteGeom);
}

/* **************************************************************** */
RADIANCEMETHOD GalerkinRadiosity = {
  "Galerkin", 3,
  "Galerkin Radiosity",
  "galerkinButton",
  GalerkinDefaults, 
  ParseGalerkinOptions, 
  PrintGalerkinOptions,
  InitGalerkin,
  DoGalerkinOneStep,
  TerminateGalerkin,
  GetRadiance,
  CreatePatchData,
  PrintPatchData,
  DestroyPatchData,
  CreateGalerkinControlPanel,
  (void (*)(void *))NULL,
  ShowGalerkinControlPanel,
  HideGalerkinControlPanel,
  GetGalerkinStats,
  GalerkinRender,
  (void (*)(void))NULL,
  GalerkinUpdateMaterial,
  GalerkinWriteVRML
};