/*
 *	fl.c
 *
 *	fracland: a fractal landscape primitive for rayshade.
 *
 * 	Copyright (C) 1992, Lawrence K. Coffin.
 *	Parts Copyright (C) 1989, 1991, Craig E. Kolb.
 *	All rights reserved.
 * 
 *	Parts of this program are from "Numerical Recipes in C" by William H.
 *	Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.
 *
 *	This software may be freely copied, modified, and redistributed
 *	provided that this copyright notice is preserved on all copies.
 *
 *	You may not distribute this software, in whole or in part, as part of
 *	any commercial product without the express consent of the authors.
 * 
 *	There is no warranty or other guarantee of fitness of this software
 *	for any purpose.  It is provided solely "as is".
 */

#include "geom.h"
#include "hf.h"
#include "fl.h"
#include <math.h>
#if defined(AMIGA) && defined(__GNUC__)
#define pow mypow
extern double mypow(const double x, const double y);
#endif
#include <stdlib.h>

static Methods *iHfMethods = NULL;
static char hfName[] = "heighfield";

static void integrate_grid(), QueueTri();
static int DDA2D(), CheckCell();
static Float intHftri();
static float minalt(), maxalt();

typedef struct {
	int stepX, stepY;
	Float tDX, tDY;
	float minz, maxz;
	int outX, outY;
	Vector cp, pDX, pDY;
} Trav2D;

hfTri *CreateHfTriangle(), *GetQueuedTri();

extern unsigned long HFTests, HFHits;

float rv();

Hf *
FlCreate(seed,subdiv)
long int seed;
int subdiv;
{
	Hf *hf;
	FILE *fp;
	float val, *maxptr, *minptr;
	int i, j, n;
	int in,i2,x,y;
	float rc, temp;


	seednrand(seed);

	hf = (Hf *)Malloc(sizeof(Hf));
	/*
	 * Make the following an option someday.
	 */
	hf->BestSize = BESTSIZE;
	/*
	 * Store the inverse for faster computation.
	 */
	hf->iBestSize = 1. / (float)hf->BestSize;
	/*
	 * Get HF size.
	 */
	n = pow(2.0,(double)subdiv);
	hf->size = n+1;
		
	hf->data = (float **)share_malloc(hf->size * sizeof(float *));

	for(i = 0; i < hf->size; i++){
		hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
	}



        hf->data[0][0] = 0;
        hf->data[0][n] = 0;
        hf->data[n][0] = 0;
        hf->data[n][n] = 0;
        for (i=0; i<subdiv; i++){
                in = (double)n/pow(2.0,(double)i);
                i2 = in/2;
                rc = 1.0/pow(2.0,(double)i);
                for (x=0; x<n; x+=in){
                        for (y=0; y<n; y+=in){
				hf->data[x+i2][y   ] = (hf->data[x   ][y   ] + hf->data[x+in][y   ])/2 + rv()*rc;
				hf->data[x   ][y+i2] = (hf->data[x   ][y   ] + hf->data[x   ][y+in])/2 + rv()*rc;
				hf->data[x+in][y+i2] = (hf->data[x+in][y   ] + hf->data[x+in][y+in])/2 + rv()*rc;
				hf->data[x+i2][y+in] = (hf->data[x   ][y+in] + hf->data[x+in][y+in])/2 + rv()*rc;
				hf->data[x+i2][y+i2] = (hf->data[x][y] + hf->data[x+in][y] + hf->data[x][y+in] + hf->data[x+in][y+in])/4 + rv()*rc;
                        }
                }
        }
/*
 *  This flips the points around so the surfaces created here are the same
 *	as the surfaces created with the program "fracland".
 */
	for (x = 0; x < hf->size; x++){
		for (y = x; y < hf->size; y++){
			temp = hf->data[x][y];
			hf->data[x][y] = hf->data[y][x];
			hf->data[y][x] = temp;
		}
	}

	for (i = 0; i < hf->size; i++) {
		for (j = 0; j < hf->size; j++) {
			val = hf->data[i][j];
			if (val <= HF_UNSET) {
				hf->data[i][j] = HF_UNSET;
				/*
				 * Don't include the point in min/max
				 * calculations.
				 */
				continue;
			}
			if (val > hf->maxz)
				hf->maxz = val;
			if (val < hf->minz)
				hf->minz = val;
		}
	}

	/*
	 * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
	 */
	for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
				hf->levels++)
			;
	hf->levels++;
	hf->qsize = CACHESIZE;
	hf->q = (hfTri **)Calloc((unsigned)hf->qsize, sizeof(hfTri *));
	hf->qtail = 0;

	hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
	hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
	hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
	hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));

	hf->spacing[0] = hf->size -1;
	hf->lsize[0] = (int)hf->spacing[0];
	hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
	hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
	/*
	 * Compute initial bounding boxes
	 */
	for (i = 0; i < hf->lsize[0]; i++) {
		hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
		hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
		maxptr = hf->boundsmax[0][i];
		minptr = hf->boundsmin[0][i];
		for (j = 0; j < hf->lsize[0]; j++) {
			*maxptr++ = maxalt(i, j, hf->data) + EPSILON;
			*minptr++ = minalt(i, j, hf->data) - EPSILON;
		}
	}

	for (i = 1; i < hf->levels; i++) {
		hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
		hf->lsize[i] = (int)hf->spacing[i];
		if ((Float)hf->lsize[i] != hf->spacing[i])
			hf->lsize[i]++;
		hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
		hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
		for (j = 0; j < hf->lsize[i]; j++) {
			hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
							sizeof(float));
			hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
							sizeof(float));
		}
		integrate_grid(hf, i);
	}

	hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0.0;
	hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1.0;
	hf->boundbox[LOW][Z] = hf->minz;
	hf->boundbox[HIGH][Z] = hf->maxz;

	return hf;
}
/*
 * Generate a random value from -0.5 to 0.5
 */

float rv(){
        float y;
        y = (float) nrand() - 0.5;
        return (y);
}

Methods *
FlMethods()
{
	if (iHfMethods == (Methods *)NULL) {
		iHfMethods = MethodsCreate();
		iHfMethods->create = (GeomCreateFunc *)FlCreate;
		iHfMethods->methods = FlMethods;
		iHfMethods->name = HfName;
		iHfMethods->intersect = HfIntersect;
		iHfMethods->normal = HfNormal;
		iHfMethods->uv = HfUV;
		iHfMethods->bounds = HfBounds;
		iHfMethods->stats = HfStats;
		iHfMethods->checkbounds = TRUE;
		iHfMethods->closed = FALSE;
	}
	return iHfMethods;
}


/*
 * Build min/max altitude value arrays for the given grid level.
 */
static void
integrate_grid(hf, level)
Hf *hf;
int level;
{
	int i, j, k, l, ii, ji;
	float max_alt, min_alt;
	float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
	int insize, fromsize, fact;

	maxinto = hf->boundsmax[level];
	mininto = hf->boundsmin[level];
	insize = hf->lsize[level];
	frommax = hf->boundsmax[level-1];
	frommin = hf->boundsmin[level-1];
	fact = hf->BestSize;
	fromsize = hf->lsize[level-1];

	ii = 0;

	for (i = 0; i < insize; i++) {
		ji = 0;
		for (j = 0; j < insize; j++) {
			max_alt = HF_UNSET;
			min_alt = -HF_UNSET;
			for (k = 0; k <= fact; k++) {
				if (ii+k >= fromsize)
					continue;
				maxptr = &frommax[ii+k][ji];
				minptr = &frommin[ii+k][ji];
				for (l = 0; l <= fact; l++,maxptr++,minptr++) {
					if (ji+l >= fromsize)
						continue;
					if (*maxptr > max_alt)
						max_alt = *maxptr;
					if (*minptr < min_alt)
						min_alt = *minptr;
				}
			}
			maxinto[i][j] = max_alt + EPSILON;
			mininto[i][j] = min_alt - EPSILON;
			ji += fact;
		}
		ii += fact;
	}
}
/*
 * Return maximum height of cell indexed by y,x.  This could be done
 * as a macro, but many C compliers will choke on it.
 */
static float
minalt(y,x,data)
int x, y;
float **data;
{
	float  min_alt;

	min_alt = min(data[y][x], data[y+1][x]);
	min_alt = min(min_alt, data[y][x+1]);
	min_alt = min(min_alt, data[y+1][x+1]);
	return min_alt;
}

/*
 * Return maximum cell height, as above.
 */
static float
maxalt(y,x,data)
int x, y;
float **data;
{
	float  max_alt;

	max_alt = max(data[y][x], data[y+1][x]);
	max_alt = max(max_alt, data[y][x+1]);
	max_alt = max(max_alt, data[y+1][x+1]);
	return max_alt;
}