/* mrvisibility.c: Multi-Resolution Visibility ala Sillion & Drettakis, "Multi-Resolution
 * Control of Visibility Error", SIGGRAPH '95 p145 */

#include <math.h>
#include "mrvisibility.h"
#include "element.h"
#include "geom.h"
#include "sgl.h"
#include "error.h"
#include "pools.h"
#include "shadowcaching.h"

double GeomMultiResolutionVisibility(GEOM *geom, RAY *ray, float rcvdist, float srcsize, float minfeaturesize)
{
  VECTOR vtmp;
  float tmin, tmax, t, fsize, *bbx;
  ELEMENT *clus = (ELEMENT *)(geom->radiance_data);
  extern GEOM *excludedGeom1, *excludedGeom2;	/* in geom.c */
  HITREC hitstore;

  if (geom==excludedGeom1 || geom==excludedGeom2)
    return 1.;

  if (!geom->bounded) 
    Fatal(-1, "GeomMultiResolutionVisibility", "Don't know what to do with unbounded geoms");

  fsize = HUGE;
  tmin = rcvdist*EPSILON;
  tmax = rcvdist;
  bbx = geom->bounds;

  /* Check ray/bounding volume intersection and compute feature size of 
   * occluder. */
  VECTORADDSCALED(ray->pos, tmin, ray->dir, vtmp);
  if (OutOfBounds(&vtmp, bbx)) {
    if (!BoundsIntersectingSegment(ray, bbx, &tmin, &tmax)) 
      return 1.;  /* ray doesn't intersect the bounding box of the GEOM within 
		   * distance interval tmin ... tmax. */

    if (clus) {
      /* Compute feature size using equivalent blocker size of the occluder */
      t = (tmin + tmax) / 2.;	/* put the centre of the equivalent blocker halfway tmin and tmax */
      fsize = srcsize + rcvdist/t * (clus->bsize - srcsize);
    }
  }

  if (fsize < minfeaturesize) {
    double kappa, vol;
    vol = (bbx[MAX_X] - bbx[MIN_X] + EPSILON) * (bbx[MAX_Y] - bbx[MIN_Y] + EPSILON) * (bbx[MAX_Z] - bbx[MIN_Z] + EPSILON);
    kappa = clus->area / (4. * vol);
    return exp(-kappa * (tmax - tmin));
  } else {
    if (GeomIsAggregate(geom))
      return GeomListMultiResolutionVisibility(GeomPrimList(geom), ray, rcvdist, srcsize, minfeaturesize);
    else {
      HITREC *hit;
      if ((hit = PatchListIntersect(GeomPatchList(geom), ray, rcvdist*EPSILON, &rcvdist, HIT_FRONT|HIT_ANY, &hitstore))) {
	AddToShadowCache(hit->patch);
	return 0.;
      } else
	return 1.;
    }
  }
}

double GeomListMultiResolutionVisibility(GEOMLIST *occluderlist, RAY *ray, float rcvdist, float srcsize, float minfeaturesize)
{
  double vis = 1.;

  while (occluderlist) {
    double v;
    v = GeomMultiResolutionVisibility(occluderlist->geom, ray, rcvdist, srcsize, minfeaturesize);
    if (v < EPSILON)
      return 0.;
    else
      vis *= v;
    
    occluderlist = occluderlist->next;
  }

  return vis;
}

/* ************* Equivalent blcoker size determination. ************** */
#ifdef DEBUG
#include "render.h"
#endif /*DEBUG*/

/* sgl context for determining equivalent blocker sizes */
static SGL_CONTEXT *sgl;

/* needed for eroding and expanding */
static unsigned char *buf1, *buf2;

/* geoms will be rendered into a frame buffer of this size for determining
 * the equivalent blocker size. */
#define FBSIZE 30

/* Creates an sgl context needed for determining the equivalent blocker size
 * of some objects. */
void BlockerInit(void)
{
  sgl = sglOpen(FBSIZE, FBSIZE);
  sglDepthTesting(TRUE);

  buf1 = (unsigned char *)Alloc(FBSIZE * FBSIZE);
  buf2 = (unsigned char *)Alloc(FBSIZE * FBSIZE);
}

/* Destroys the sgl context created by BlockerInit(). */
void BlockerTerminate(void)
{
  Free((char *)buf2, FBSIZE * FBSIZE);
  Free((char *)buf1, FBSIZE * FBSIZE);

  sglClose(sgl);
}

static void BlockerRenderPatch(PATCH *patch)
{
  VECTOR v[4]; int i;

  /* no backface culling */

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

  sglSetColor((SGL_PIXEL)patch);
  sglPolygon(patch->nrvertices, v);
}

static void BlockerRenderGeomRecursive(GEOM *geom)
{
  if (GeomIsAggregate(geom)) {
    GeomListIterate(GeomPrimList(geom), BlockerRenderGeomRecursive);
  } else {
    PatchListIterate(GeomPatchList(geom), BlockerRenderPatch);
  }
}

/* Renders the GEOM as seen form the given direction into the sgl frame buffer. */
static float *BlockerRenderGeom(GEOM *geom, VECTOR viewdir)
{
  POINT eye, centre;
  VECTOR up = {0., 0., 1.};
  static BOUNDINGBOX bbx;
  float *gbounds;
  TRANSFORM lookat;
  SGL_CONTEXT *prev_sgl_context;

  /* set up an orthographics projection, centre is the centre of the boundingbox of
   * the geom. eye is the centre - viewdir. */
  gbounds = GeomBounds(geom);
  VECTORSET(centre, (gbounds[MIN_X] + gbounds[MAX_X])/2., 
	            (gbounds[MIN_Y] + gbounds[MAX_Y])/2., 
	            (gbounds[MIN_Z] + gbounds[MAX_Z])/2.);
  VECTORSUBTRACT(centre, viewdir, eye);
  if (fabs(VECTORDOTPRODUCT(up, viewdir)) > 1.-EPSILON)
    VECTORSET(up, 0., 1., 0.);
  lookat = LookAt(eye, centre, up);

  BoundsTransform(gbounds, &lookat, bbx);

  prev_sgl_context = sglMakeCurrent(sgl);
  sglLoadMatrix(Ortho(bbx[MIN_X], bbx[MAX_X], bbx[MIN_Y], bbx[MAX_Y], -bbx[MAX_Z], -bbx[MIN_Z]));
  sglMultMatrix(lookat);

  sglClear((SGL_PIXEL)NULL, SGL_ZMAX);
  BlockerRenderGeomRecursive(geom);

  sglMakeCurrent(prev_sgl_context);
  return bbx;
}

/* Copies the pixels in the viewport of the sgl frame buffer used for blocker size
 * determination contiguously into 'buf'. Copies a zero for each background pixel 
 * and a 1 for each non background pixel. */
static void BlockerCopyFrameBuffer(unsigned char *buf)
{
  int i, j; 
  SGL_PIXEL *pix;
  unsigned char *p;

  for (p=buf, j=0; j<sgl->vp_height; j++) {
    pix = sgl->fbuf + j * sgl->width;
    for (i=0; i<sgl->vp_width; i++, pix++) {
      *p++ = (*pix == (SGL_PIXEL)0) ? 0 : 1;
    }
  }
}

#ifdef DEBUG
static void ShowBuffer(unsigned char *buf, int width, int height)
{
  int i, j;
  RGB rgb[width];
  unsigned char *p = buf;

  for (j=0; j<height; j++) {
    for (i=0; i<width; i++) {
      if (*p++) 
	rgb[i] = Black;
      else
	rgb[i] = Yellow;
    }
    RenderPixels(0, height-j-1, width, rgb);
  }
}
#endif /*DEBUG*/

/* returns FALSE if buf1 is all background. Result of erosion
 * operation is stored in buf2. */
static int Erode(unsigned char *buf1, unsigned char *buf2, int width, int height)
{
  int i, j, allbackground=TRUE, allforeground=TRUE;
  unsigned char *left, *right, *above, *below;

  left = buf1-1;
  right = buf1+1;
  below = buf1-width;
  above = buf1+width;

  for (j=0; j<height; j++) {
    for (i=0; i<width; i++, buf1++, buf2++, left++, right++, below++, above++) {
      if (*buf1) 
	allbackground=FALSE;
      else
	allforeground=FALSE;

      *buf2 = ((*buf1) &&
	       (j==0 || *below) &&
	       (j==height-1 || *above) &&
	       (i==0 || *left) &&
	       (i==width-1 || *right));
    }
  }

  return !(allbackground);
}

/* returns FALSE if buf1 is all foreground. Result of expansion
 * operation is stored in buf2. */
static int Expand(unsigned char *buf1, unsigned char *buf2, int width, int height)
{
  int i, j, allbackground=TRUE, allforeground=TRUE;
  unsigned char *left, *right, *above, *below;

  left = buf1-1;
  right = buf1+1;
  below = buf1-width;
  above = buf1+width;

  for (j=0; j<height; j++) {
    for (i=0; i<width; i++, buf1++, buf2++, left++, right++, below++, above++) {
      if (*buf1) 
	allbackground=FALSE;
      else
	allforeground=FALSE;

      *buf2 = ((*buf1) ||
	       (j>0 && *below) ||
	       (j<height-1 && *above) ||
	       (i>0 && *left) ||
	       (i<width-1 && *right));
    }
  }

  return !(allforeground);
}

/* Determines the equivalent blocker size of the GEOM when viewed in the given
 * direction. */
double GeomBlcokerSizeInDirection(GEOM *geom, VECTOR viewdir)
{
  int width, height, nrerosions, nrexpansions;
  double pixelsize;
  float *bbx;
  
  /* render the patches in the geom into the scratch frame buffer. Use an 
   * orthographic projection along the direction 'viewdir' */
  bbx = BlockerRenderGeom(geom, viewdir);

  width = sgl->vp_width;
  height = sgl->vp_height;

  /* erode until nothing remains, count how many erosions were necessary */
  nrerosions = 0;
  BlockerCopyFrameBuffer(buf1);
  while (Erode(buf1, buf2, width, height) && nrerosions < MAX(width, height)/2) {
    /* swap the buffers */
    unsigned char *tmp;
    tmp = buf1; buf1 = buf2; buf2 = tmp;
    nrerosions++;
  }

  /* expand until no background remains, count how many expansions were necessary */
  nrexpansions = 0;
  BlockerCopyFrameBuffer(buf1);
  while (Expand(buf1, buf2, width, height) && nrexpansions < MAX(width, height)/2) {
    /* swap the buffers */
    unsigned char *tmp;
    tmp = buf1; buf1 = buf2; buf2 = tmp;
    nrexpansions++;
  }

  pixelsize = MAX((bbx[MAX_X] - bbx[MIN_X])/width, (bbx[MAX_Y] - bbx[MIN_Y])/height);
  return (double)(2*MAX(nrerosions, nrexpansions)) * pixelsize;
}

/* Determines a single equivalent blocker size for the GEOM which can be used for
 * any direction (approximately). */
double GeomBlcokerSize(GEOM *geom)
{
  double blocker_size, s;
  VECTOR viewdir;
  
  VECTORSET(viewdir, 1., 0., 0.);
  s = GeomBlcokerSizeInDirection(geom, viewdir);
  blocker_size = s;

  VECTORSET(viewdir, 0., 1., 0.);
  s = GeomBlcokerSizeInDirection(geom, viewdir);
  blocker_size = MAX(blocker_size, s);

  VECTORSET(viewdir, 0., 0., 1.);
  s = GeomBlcokerSizeInDirection(geom, viewdir);
  blocker_size = MAX(blocker_size, s);

  return blocker_size;
}