/* globallines.c: generic global line tracing engine */

#include "globallines.h"

#include "scene.h"
#include "mymath.h"
#include "render.h"
#include "QMC/nied31.h"
#include "error.h"

/* World bounding sphere, available after InitGlobalLines() */
VECTOR glin_Center;
float glin_Radius;

/* Generate a bounding sphere */
static void GenerateBoundingSphere(void)
{
  float *bbx = WorldGrid->bounds;
  VECTOR diagonal;
  VECTORSET(glin_Center,
	    (bbx[MIN_X] +  bbx[MAX_X])/2,   
	    (bbx[MIN_Y] +  bbx[MAX_Y])/2,
	    (bbx[MIN_Z] +  bbx[MAX_Z])/2);
  VECTORSET(diagonal,
	    bbx[MAX_X] -  bbx[MIN_X],   
	    bbx[MAX_Y] -  bbx[MIN_Y],
	    bbx[MAX_Z] -  bbx[MIN_Z] );
  glin_Radius =  VECTORNORM(diagonal)/2;
}

static double *DefSample4D(int n)
{
  unsigned int *zeta = Nied31(n);
  static double xi[4];
  xi[0] = zeta[0] * RECIP;
  xi[1] = zeta[1] * RECIP;
  xi[2] = zeta[2] * RECIP;
  xi[3] = zeta[3] * RECIP;
  return xi;
}
static double* (*Sample4D)(int n) = DefSample4D;

/* Generates a line with two random points on bounding sphere.
 * Ray direction is not normalized: t values lay between 0 and 1
 * for points inside the scenes bounding sphere. */
static RAY GenerateGlobalLine(int n)
{
  RAY ray;
  VECTOR p1, p2;
  double alea;
  double direc1,direc2;
  double *xi;

  xi = Sample4D(n);

  alea = 1.0-2.0*xi[0];
  direc1 = 2.0*M_PI*xi[1];
  direc2 = acos(alea);

  VECTORSET(p1,cos(direc1)*sin(direc2),  
	       sin(direc1)*sin(direc2),
	       cos(direc2));
  VECTORSCALE(glin_Radius, p1,p1);
  VECTORADD(glin_Center,p1,p1);

  alea = 1.0-2.0*xi[2];
  direc1 = 2.0*M_PI*xi[3];
  direc2 = acos(alea);

  VECTORSET(p2,cos(direc1)*sin(direc2),  
	    sin(direc1)*sin(direc2),
	    cos(direc2));
  VECTORSCALE(glin_Radius, p2,p2);
  VECTORADD(glin_Center,p2,p2);

  ray.pos = p1; 
  VECTORSUBTRACT(p2, p1, ray.dir);
  /* ray direction is not normalized!!! */
  /*
  RenderSetColor(&Yellow);
  RenderLine(&p1, &p2);
  */
  return ray;
}

static HITLIST *TraceGlobalLine(RAY *ray)
{
  return AllGridIntersections((HITLIST *)NULL, WorldGrid, ray, -EPSILON, 1.+EPSILON, HIT_FRONT|HIT_BACK|HIT_PATCH|HIT_POINT);
  /* HIT_POINT is only needed to draw the points in ShowHit() */
}

/* "returns" TRUE if the hit record describes a frontal hit and "FALSE" if
 * it is a hit from the back. */
#define IsFrontalHit(hit)	(hit->flags & HIT_FRONT)

/* Returns -1 if hit1 comes before hit2, +1 if hit1 should come after
 * hit2 and 0 if both hits are coincident are are both frontal hits
 * or hits from the back (a situation that occurs e.g. for coplanar
 * polygons) */
int DefCompareHits(HITREC *hit1, HITREC *hit2)
{
  int code=0;

  if (hit1->dist < hit2->dist-EPSILON)
    code = -1;
  else if (hit1->dist > hit2->dist+EPSILON)
    code = +1;
  else {
    /* coincident hitpoints */
    if (IsFrontalHit(hit1) && !IsFrontalHit(hit2))
      code = -1;	/* frontal hits come before hits from the back */
    else if (!IsFrontalHit(hit1) && IsFrontalHit(hit2))
      code = +1;
    else
      code = 0;		/* both hits are frontal or from back */
  }
  return code;
}
static int (*CompareHits)(HITREC*, HITREC*) = DefCompareHits;

/* Determines "spans" connecting a hit from the back with the next frontal
 * hit. Calls do_span for each found span. Requires that the hitlist has
 * been sorted using the default hitlist sorting routine. 
 * Returns number of spans processed. */
int ProcessSpans(HITLIST *hits, void (*do_span)(HITREC *start, HITREC *end))
{
  int nrspans = 0;
  HITREC *start = NULL;
  ForAllHits(hit, hits) {
    if (IsFrontalHit(hit)) {
      if (start) {
	do_span(start, hit); nrspans++;
	start = (HITREC *)NULL;
      }
      /* else (frontal hit without preceeding hit from the back) skip it */
    } else
      start = hit;
  } EndForAll;
  return nrspans;
}

/* example of a span processing routine */
void PrintSpan(HITREC *start, HITREC *end)
{
  fprintf(stderr, "start = "); PrintHit(stderr, start);
  fprintf(stderr, "end   = "); PrintHit(stderr, end);
}

/* default hit processing routine: calles ProcessSpans() with
 * PrintSpan as a span processing routine. */
void DefProcessHits(HITLIST *hits)
{
  ProcessSpans(hits, PrintSpan);
}
static void (*ProcessHits)(HITLIST *hits) = DefProcessHits;

static int inited = FALSE;

/* Initialises global line tracing.
 * - sample4d is a routine returning a 4D sample vector with given
 * index, used to generate a global line. This parameter can be 
 * NULL, in which case 4D Niederreiter sampling is used.
 * - compare_hits is a routine that compares to hit records. Use
 * 'DefCompareHits' in order to sort global line hit points in such
 * a way that ProcessSpans() can be used for processing the hits.
 * Also this parameter can be NULL, in which case the hits are
 * not sorted.
 * - process_hits is a pointer to a routine for processing the 
 * sorted global line hits. If sorted using DefCompareHits, process_hits
 * can call ProcessSpans with a routine for handling spans of mutually
 * visible points. */
void InitGlobalLines(double* (*sample4d)(int index),
		     int (*compare_hits)(HITREC *hit1, HITREC *hit2),
		     void (*process_hits)(HITLIST *hits))
{
  if (!World) {
    Error("InitGlobalLines", "No world!");
    inited = FALSE;
    return;
  }
  inited = TRUE;
  GenerateBoundingSphere();
  Sample4D = sample4d ? sample4d : DefSample4D;
  CompareHits = compare_hits;
  ProcessHits = process_hits;
}

/* Generates and traces a global line with given index number using
 * the two-points-on-a-bounding-sphere method, using sample4d()
 * for sampling the points on the sphere. If the global line 
 * intersects surfaces in the scene, the hit points are sorted using
 * compare_hits(). Eventually, process_hits() is called for processing
 * the sorted hit points. 
 * Returns FALSE is the generated global line did not hit the scene and
 * returns TRUE if it did. */
int DoGlobalLine(int index)
{
  RAY ray;
  HITLIST *hits;

  if (!inited) {
    Fatal(-1, "DoGlobalLine", "Call InitGlobalLines() first before using DoGlobalLine()");
  }

  ray = GenerateGlobalLine(index);
  hits = TraceGlobalLine(&ray);
  if (!hits)
    return FALSE;

  if (CompareHits) {
    hits = HitListSort(hits, CompareHits);
  }

  if (ProcessHits) {
    ProcessHits(hits);
  }

  DestroyHitlist(hits);

  return TRUE;
}