/*
 * ifs.c
 *
 * Philippe Bekaert - bekaert@cs.kuleuven.ac.be - 19 may 93
 */
#include <stdlib.h>
#include <math.h>
#include "geom.h"
#include "ifs.h"
#include "bounds.h"
#include "libcommon/sampling.h"
#include "libshade/viewing.h"
#include "libcommon/common.h"

static Methods *iIfsMethods = NULL, *iCIfsMethods = NULL;
static char ifsName[] = "ifs", CifsName[] = "condensed ifs";
unsigned long IfsTests, IfsHits, CIfsTests, CIfsHits;;

#define OPTIMAL_BV  /* for better bounding spheres */

/* return less than zero if first map has higher rang number than the second,
 *	  zero if two maps have equal rang numbers, 
 *        greater than zero if first map has lower rang number than the second */
static int rangcmp(map1, map2)
IfsMap *map1, *map2;
{
  	return map2->rang - map1->rang;
}

/* computes radius of smallest sphere with centre (x,y,z) which encloses
 * the n spheres with centre c[i] and radius r[i]: to be minimized */
static Float REval(n, c, r, x, y, z)
int n;
Vector *c; Float *r;
Float x, y, z;
{
	Vector C; 
	Float R, R1;
	int i;

	C.x = x; C.y = y; C.z = z;
	R = VecDist(&C, &c[0]) + r[0];
	for (i=1; i<n; i++)
		if ((R1 = VecDist(&C, &c[i]) + r[i]) > R) R = R1;
	return R;
}

static Float SmallestBounding(n, c, r, C)
int n;
Vector *c; Float *r;
Vector *C;
{
	Float R, R1, step;
	int i, dir, go_on;

/* initial guess for center point: gravitation centre of centre points */
	*C = c[0];
	for (i=1; i<n; i++)
		VecAdd(*C, c[i], C);
	VecScale(1./(float)n, *C, C);

/* calculate radius and minimize */
	R = REval(n, c, r, C->x, C->y, C->z);
	step = R/8.;
#ifdef OPTIMAL_BV
	do {
		dir = 0;
		do {
			go_on = FALSE;
			if (dir != -1 && (R1 = REval(n, c, r, C->x + step, C->y, C->z)) < R) 
			      	{ R = R1; C->x += step; dir = 1; go_on = TRUE;}
			if (dir !=  1 && (R1 = REval(n, c, r, C->x - step, C->y, C->z)) < R) 
		      		{ R = R1; C->x -= step; dir = -1; go_on = TRUE;}
		} while (go_on); 

		dir = 0;
		do {
			go_on = FALSE;
			if (dir != -1 && (R1 = REval(n, c, r, C->x, C->y + step, C->z)) < R) 
			      	{ R = R1; C->y += step; dir = 1; go_on = TRUE;}
			if (dir !=  1 && (R1 = REval(n, c, r, C->x, C->y - step, C->z)) < R) 
		      		{ R = R1; C->y -= step; dir = -1; go_on = TRUE;}
		} while (go_on); 

		dir = 0;
		do {
			go_on = FALSE;
			if (dir != -1 && (R1 = REval(n, c, r, C->x, C->y, C->z + step)) < R) 
			      	{ R = R1; C->z += step; dir = 1; go_on = TRUE;}
			if (dir !=  1 && (R1 = REval(n, c, r, C->x, C->y, C->z - step)) < R) 
		      		{ R = R1; C->z -= step; dir = -1; go_on = TRUE;}
		} while (go_on); 

		step /= 2.;
	} while (step > EPSILON);
#endif /* OPTIMAL_BV */

	return R;
}

static void IfsPartialBounding(nrmaps, maps, condensed, C, R, minsize, maxdepth)
int nrmaps, condensed, maxdepth;
IfsMap *maps;
Vector *C;
Float *R, minsize;
{
	int i, N, d;
	Float *r, R1, delta, maxcontractivity;
	Vector *c, C1;

/* largest Lipschitz constant of the maps */
	for (i = 0, maxcontractivity = 0.; i < nrmaps; i++)
		if (maps[i].contractivity > maxcontractivity) 
			maxcontractivity = maps[i].contractivity;

/* for IFS without condensation: make new sphere which contains the images
 * of this sphere and repeat as long the new sphere gets smaller. For IFS
 * with condensation set: make new sphere which contains the old and it's
 * images and repeat as long as it gets larger */
	N = condensed ? nrmaps+1 : nrmaps;
        c = (Vector *)Malloc(N * sizeof(Vector)); 
        r =  (Float *)Malloc(N * sizeof( Float)); 

        for (delta = R1 = *R, C1 = *C, d = 0;
	     delta > EPSILON && delta >= minsize * maxcontractivity && d <= maxdepth; 
	     delta *= maxcontractivity, d++) {
                if (!condensed)
                        if (R1 > *R) break; 
                *C = C1; *R = R1;

                for (i=0; i<nrmaps; i++) {
			c[i] = *C;
                        PointTransform(&c[i], &maps[i].matrix);
			r[i] = *R * maps[i].contractivity;
		}
                if (condensed) {
                        c[nrmaps] = *C;
			r[nrmaps] = *R;
		}

		R1 = SmallestBounding(N, c, r, &C1);
        }
	if (condensed) { *C = C1; *R = R1; }
        free(c); free(r);

        *R += EPSILON;          /* to be on the safe side */
}

/* Computes a sphere which is supposed to fit the attractor of the IFS quite tightly */
static void IfsBoundingSphere(ifs, R, C)
Ifs *ifs;
Float *R;
Vector *C;
{
        Float R1;
        int i, n, condensed;
	IfsMap *m;

	condensed = ifs->condensation ? TRUE : FALSE;
        if (!condensed) {
/* first guess: a sphere with centre in the origin which contains all it's
 * images: The radius must be equal or larger 
 *
 *      |B| / (1-s) 
 *
 * for each IFS map x -> Ax+B with contractivity s.
 * hierarchical IFS: use only the maps with highest rang number */
                C->x = C->y = C->z = 0.;
                for (i=0, *R=0.; i < ifs->nrmaps && ifs->maps[i].rang == ifs->maps->rang; i++) 
                        if ((R1 = VecNorm(ifs->maps[i].matrix.translate) /
                                (1. - ifs->maps[i].contractivity)) > *R) *R = R1;
        } else {
/*
 * construct a sphere around the condensation bounding box 
 */
                C->x = (ifs->condensation->bounds[HIGH][X] + ifs->condensation->bounds[LOW][X])/2.;
                C->y = (ifs->condensation->bounds[HIGH][Y] + ifs->condensation->bounds[LOW][Y])/2.;
                C->z = (ifs->condensation->bounds[HIGH][Z] + ifs->condensation->bounds[LOW][Z])/2.;
		*R = ifs->condensize/2.;
        }

/* process maps in groups of same rang */
	for (i = 0; i < ifs->nrmaps; i += n, condensed = TRUE) {
/* count the number of maps with same rang number as ifs->maps[i] */
		for (m = &ifs->maps[i], n = 0; i+n < ifs->nrmaps && m[n].rang == m->rang; n++) {}
/* construct a bounding sphere for the partial IFS defined by these maps of same rang */
		IfsPartialBounding(n, m, condensed, C, R, ifs->options.size, ifs->options.depth);
	}
}

/* global variables for depth first processing of the tree of images */
static IfsMap *maps;
static int nrmaps, maxdepth, condensation, oppenheimer_tree;
static Float minsize, isize;
static Float ibounds[2][3], cbounds[2][3], lbounds[2][3];
static Vector ipoint;

static IfsBB(map, depth, size, rang)
RSMatrix *map;
int depth;
Float size;
int rang;
{
	int i;

	if (depth >= maxdepth || (minsize > EPSILON && size <= minsize)) {
		Vector p;
		Float b[2][3];
		if (oppenheimer_tree) {	/* enlarge with bounds of image of leaves */
			BoundsCopy(lbounds, b);
			BoundsTransform(map, b);
			BoundsEnlarge(ibounds, b);	
		} else if (condensation) {	/* enlarge with bounds of condensation */
			BoundsCopy(cbounds, b);
			BoundsTransform(map, b);
			BoundsEnlarge(ibounds, b);	
		} else {		/* enlarge with sphere with centre at image of ipoint */
			p = ipoint;
			PointTransform(&p, map);
			BoundsEnlargeSphere(ibounds, &p, size/2);
		}
		return;
	}

	if (condensation) {		/* enlarge with bounds of image of condensation set */
		Float b[2][3];
		BoundsCopy(cbounds, b);
		BoundsTransform(map, b);
		BoundsEnlarge(ibounds, b);	
	}

	for (i=0; i<nrmaps; i++) {
		RSMatrix m;
		Float lambda[3], siz;

		if (maps[i].rang < rang) break;		/* for hierarchical IFS */
		MatrixMult(&maps[i].matrix, map, &m);	
	        if (maps[i].uniform)			
        		siz = size * maps[i].contractivity;
	        else {
        	        MatrixScaling(&m, lambda);
                	siz = isize * Lipschitz(lambda);
	        }
		IfsBB(&m, depth+1, siz, maps[i].rang);	
	}
}

static IfsBoundingBox(ifs, bounds)
Ifs *ifs;
Float bounds[2][3];
{
	RSMatrix unity_transform;

	MatrixInit(&unity_transform);
	maps = ifs->maps;
	nrmaps = ifs->nrmaps;
	condensation = (ifs->condensation != NULL);
	oppenheimer_tree = (ifs->leaves != NULL);
	maxdepth = ifs->options.depth;
	minsize = ifs->options.size;
	if (minsize < EPSILON)
		minsize = ifs->extent.r/20;
	isize = ifs->condensation ? ifs->condensize : ifs->size ;
	ipoint = ifs->extent.c;

	BoundsInit(ibounds);
	if (ifs->condensation) 
		BoundsCopy(ifs->condensation->bounds, cbounds);
	if (ifs->leaves)
		BoundsCopy(ifs->leaves->bounds, lbounds);

	IfsBB(&unity_transform, 0, isize, -99999999);
	BoundsCopy(ibounds, bounds);
}

/* compute a bunding box and bounding sphere for the IFS attractor */
static void IfsComputeExtent(ifs)
Ifs *ifs;
{
        Float R, sphavpa, boxavpa;
        Vector C, d;

/* compute bounds of condensation set, if any */
	if (ifs->condensation) {
                GeomComputeBounds(ifs->condensation);
		if (UNBOUNDED(ifs->condensation)) 
                	RLerror(RL_ABORT, "Can't compute bounds of IFS with unbounded condensation.\n");

                d.x = ifs->condensation->bounds[HIGH][X] - ifs->condensation->bounds[LOW][X];
                d.y = ifs->condensation->bounds[HIGH][Y] - ifs->condensation->bounds[LOW][Y];
                d.z = ifs->condensation->bounds[HIGH][Z] - ifs->condensation->bounds[LOW][Z];
		ifs->condensize = VecNorm(d);
        }

/* compute bounds of leaves, if any */
        if (ifs->leaves) {
                GeomComputeBounds(ifs->leaves);
                if (UNBOUNDED(ifs->leaves))
                        RLerror(RL_ABORT, "Oppenheimer tree with unbounded leaves.\n");

                d.x = ifs->leaves->bounds[HIGH][X] - ifs->leaves->bounds[LOW][X];
                d.y = ifs->leaves->bounds[HIGH][Y] - ifs->leaves->bounds[LOW][Y];
                d.z = ifs->leaves->bounds[HIGH][Z] - ifs->leaves->bounds[LOW][Z];
		ifs->leavesize = VecNorm(d);
        }

/* construct a sphere fitting tightly the attractor of the IFS */
        IfsBoundingSphere(ifs, &R, &C);

        ifs->extent.r = R;
        ifs->extent.rsq = R*R;
        ifs->extent.c = C;
        ifs->size = 2*R;

/* construct a box fitting tightly the attractor of the IFS */
	IfsBoundingBox(ifs, ifs->bounds);

/* compare average projected area of the box and the sphere. The avaerage projected
 * area for a convex volume is one fourth of the surface area of the volume [Glassner,
 * An Introduction to Ray Tracing, Academic Press, 1989 p212]. The volume
 * with the smallest average projected area will be used for building a bounding
 * volume hierarchy for the attractor of the IFS */
	sphavpa = 3.14 * ifs->extent.rsq;
	d.x = ifs->bounds[HIGH][X] - ifs->bounds[LOW][X];
	d.y = ifs->bounds[HIGH][Y] - ifs->bounds[LOW][Y];
	d.z = ifs->bounds[HIGH][Z] - ifs->bounds[LOW][Z];
	boxavpa = (d.x * d.y + d.y * d.z + d.z * d.x) / 2;
	
	if (sphavpa < boxavpa)	ifs->boxbounding = FALSE;
	else			ifs->boxbounding = TRUE;
	if (ifs->options.boxbounding) ifs->boxbounding = TRUE;
	if (ifs->options.spherebounding) ifs->boxbounding = FALSE;
}

/*
 * for debugging or listing the volumes that will be rendered
 */
static TransPrint(trans, out)
RSMatrix *trans;
FILE *out;
{
	int i;

	for (i=0; i<3; i++) 
		fprintf(out, "%lf %lf %lf   ", 
			(double)trans->matrix[i][0],
			(double)trans->matrix[i][1],
			(double)trans->matrix[i][2]  );
	fprintf(out, "%lf %lf %lf",
			(double)trans->translate.x,
			(double)trans->translate.y,
			(double)trans->translate.z  );
}

/* produce a list of volumes that will be rendered */
static FILE *listfp;	
static long volcount;

static ListIt(map, depth, size, rang)
RSMatrix *map;
int depth;
Float size;
int rang;
{
	int i;
	RSMatrix m;
	Float lambda[3];

	if (depth >= maxdepth || (minsize > EPSILON && size <= minsize)) {
		if (oppenheimer_tree)
			fprintf(listfp, "object leaf         ");
		else if (condensation)
			fprintf(listfp, "object condensation ");
		else
			fprintf(listfp, "object extent ");
		fprintf(listfp, "transform   ");
		TransPrint(map, listfp);
		fprintf(listfp, "\n");
		volcount++;
		return;
	}

	if (condensation) {
		fprintf(listfp, "object condensation ");
		fprintf(listfp, "transform   ");
		TransPrint(map, listfp);
		fprintf(listfp, "\n");
		volcount++;
	}

	for (i=0; i<nrmaps; i++) {
	        Float siz;
		if (maps[i].rang < rang) break;
		MatrixMult(&maps[i].matrix, map, &m);	
	        if (maps[i].uniform)			
        		siz = size * maps[i].contractivity;
	        else {
        	        MatrixScaling(&m, lambda);
                	siz = isize * Lipschitz(lambda);
	        }
		ListIt(&m, depth+1, siz, maps[i].rang);
	}
}

static IfsList(ifs, filename)
Ifs *ifs;
char *filename;
{
	RSMatrix unity_transform;

/* initialise some global variables */
	MatrixInit(&unity_transform);
	maps = ifs->maps;
	nrmaps = ifs->nrmaps;
	condensation = (ifs->condensation != NULL);
	oppenheimer_tree = (ifs->leaves != NULL);
	maxdepth = ifs->options.depth;
	minsize = ifs->options.size;
	isize = ifs->condensation ? ifs->condensize : ifs->size ;
	if (maxdepth > 999999 && minsize < EPSILON) {
		minsize = isize/100;
		fprintf(stderr, "Warning: no stopcriterium for IFS list of volumes. Using 'minsize %lf'\n", (double)minsize);
		return;
	}

/* open output file */
	if ((listfp = fopen(filename, "w")) == NULL) {
		fprintf(stderr, "IfsList: Cannot open %s for writing\n", filename);
		perror(ifs->options.listfilename);
		return;
	}

	fprintf(listfp, "/* maxdepth=%d, initial_size=%lf, minsize=%lf */\n\n",
			 maxdepth, (double)isize, (double)minsize);

	if (!ifs->boxbounding) fprintf(listfp, "/* ");
	fprintf(listfp, "name extent box  %lf %lf %lf  %lf %lf %lf", 
		(double)ifs->bounds[LOW][X],
		(double)ifs->bounds[LOW][Y],
		(double)ifs->bounds[LOW][Z],
		(double)ifs->bounds[HIGH][X],
		(double)ifs->bounds[HIGH][Y],
		(double)ifs->bounds[HIGH][Z]  );
	if (!ifs->boxbounding) fprintf(listfp, "*/");
	fprintf(listfp, "\n");

	if (ifs->boxbounding) fprintf(listfp, "/* ");
	fprintf(listfp, "name extent sphere  %lf  %lf %lf %lf", 
		(double)ifs->extent.r,
		(double)ifs->extent.c.x, 
		(double)ifs->extent.c.y, 
		(double)ifs->extent.c.z        );
	if (ifs->boxbounding) fprintf(listfp, "*/");
	fprintf(listfp, "\n");

	if (condensation)
		fprintf(listfp, "name condensation /* fill in the condensation-set object */\n");
	if (oppenheimer_tree)
		fprintf(listfp, "name leaf /* fill in the leaf object */\n");

	volcount = 0;
	fprintf(listfp, "\ngrid 1 1 1 /* change to your personal taste */\n");
	ListIt(&unity_transform, 0, isize, -99999999);
	fprintf(listfp, "end /* %ld objects */\n", volcount);

	fclose(listfp);
}

/*
 * normal weighting functions 
 */
static Float HighPassNormalWeighting(Size, size)
Float Size, size;
{
        return (Size-size)/Size;
}

static Float LowPassNormalWeighting(Size, size)
Float Size, size;
{
        return size/Size;
}

static Float ConstantNormalWeighting(Size, size)
Float Size, size;
{
        return 1.;
}


/*
 * Create & return reference to an ifs.
 */
Ifs *
IfsCreate(maps, condensation, leaves, options)
IfsTransList *maps;
Geom *condensation, *leaves;
IfsOptions *options;
{
        Ifs *ifs;
        IfsTransList *map; 
        Float contractivity, maxcontractivity;
        Float lambda[3];
        int n;

/* first: count the number of maps */
        for (n=0, map=maps; map; n++, map=map->next) {}
        if (n<=0) {
                RLerror(RL_WARN, "I don't create IFS with no maps.\n");
                return (Ifs *)NULL;
        }

/* create Ifs structure */
        ifs = (Ifs *)share_malloc(sizeof(Ifs));
        ifs->condensation = condensation;       /* NULL = no cendensation */
        ifs->leaves = leaves;   /* if not NULL -> Oppenheimer tree */
        ifs->options = *options;

/* allocate space for keeping the maps */
        ifs->nrmaps = n;
        ifs->maps = (IfsMap *)share_malloc(n*sizeof(IfsMap));

/* the maps are in reverse order */
        ifs->animated = FALSE;
        ifs->uniform_maps = TRUE;
        maxcontractivity = 0.;
        for (n=ifs->nrmaps-1, map=maps; map; n--, map=map->next) {
                ifs->maps[n].trans = map->trans;
                ifs->maps[n].transtail = map->transtail;
                ifs->maps[n].rang = map->rang;

/* is the map animated ? */
                ifs->maps[n].animated = (map->trans->animated || map->trans->next!=(Trans *)NULL);
/* if so, the IFS is animated */
                if (ifs->maps[n].animated)
                        ifs->animated = TRUE;
/* if not, we can do some computations already */
                else {
/* "compose" the transformations - they are already composed during the 
 * parsing, so, just copy them. They are also garanteed to be non-singular. */
                        MatrixCopy(&map->trans->trans,  &ifs->maps[n].matrix);
                        MatrixCopy(&map->trans->itrans, &ifs->maps[n].imatrix);

/* compute the contractivity of the map */
                        MatrixScaling(&ifs->maps[n].matrix, lambda);
                        contractivity = Lipschitz(lambda);
                        if (contractivity >= 1.) 
                                RLerror(RL_WARN, "Ifs map %d is not contractive (contractivity = %lf)\n", n+1, contractivity);
                        if (contractivity > maxcontractivity)
                                maxcontractivity = contractivity;
                        ifs->maps[n].contractivity = contractivity;
                        ifs->maps[n].uniform = Uniform(lambda);
                        if (!ifs->maps[n].uniform) ifs->uniform_maps = FALSE;
                }
        }

/* some procedures may never end for non-contractive IFS's */
/*        if (maxcontractivity >= 1.0) {
                free(ifs->maps);
                free(ifs);
                return (Ifs *)NULL;
        }
 * just warn the user. we'll see ... */

/* for hierarchical IFS: sort maps to descending rang number (those maps 
 * that can follow any other map will be processed first) */
	qsort(ifs->maps, ifs->nrmaps, sizeof(IfsMap), rangcmp);

        if (!ifs->animated) ifs->contractivity = maxcontractivity;

/* the IFS is also animated if its condensation is -- no code to support this yet */
/*        ifs->animated |= ifs->condensation->animtrans; */

/* if the IFS is not animated, compute its extent here */
        if (!ifs->animated && !ifs->condensation)
                IfsComputeExtent(ifs);

/* normal weighting function for IFS without condensation */
        if (!ifs->condensation) {
                switch (ifs->options.normalweighting) {
                        case HIGHPASS_NORMAL_WEIGHTING:
                                ifs->normalweight = HighPassNormalWeighting; break;
                        case LOWPASS_NORMAL_WEIGHTING:
                                ifs->normalweight = LowPassNormalWeighting; break;
                        case CONSTANT_NORMAL_WEIGHTING:
                                ifs->normalweight = ConstantNormalWeighting; break;
                        default: {
                                RLerror(RL_WARN, "invalid normal weighting method.\n");
                                free(ifs->maps);
                                free(ifs);
                                return (Ifs *)NULL;
                        } 
                }
        } 

        ifs->frame = -1;        /* impossible value */

        return ifs;
}

/* resolve geometric associates -- untested */
static void IfsResolveAssoc(ifs)
Ifs *ifs;
{
        Trans tr;
        Float contractivity, maxcontractivity;
        Float lambda[3];
        int i;  

        maxcontractivity = 0.;
        ifs->uniform_maps = TRUE;
        for (i = 0; i < ifs->nrmaps; i++) {
                if (ifs->maps[i].animated) {
                        TransResolveAssoc(ifs->maps[i].trans);
                        TransComposeList(ifs->maps[i].trans, &tr);
			ifs->maps[i].matrix = tr.trans;
			ifs->maps[i].imatrix = tr.itrans;

/* compute the contractivity of the map */
                        MatrixScaling(&ifs->maps[i].matrix, lambda);
                        contractivity = Lipschitz(lambda);
                        if (contractivity >= 1.) 
                                RLerror(RL_WARN, "Ifs map %d is not contractive (contractivity = %lf)\n", i+1, contractivity);
                        ifs->maps[i].contractivity = contractivity;
                        ifs->maps[i].uniform = Uniform(lambda);
                } else 
                        contractivity = ifs->maps[i].contractivity;

                if (contractivity > maxcontractivity)
                        maxcontractivity = ifs->maps[i].contractivity;
                if (!ifs->maps[i].uniform) ifs->uniform_maps = FALSE;
        }

        ifs->contractivity = maxcontractivity;
/*        if (ifs->contractivity >= 1.)
                RLerror(RL_ABORT, "Can't go on with non-contractive IFS.\n");
*/
        ifs->frame = Sampling.framenum;
}

Methods *
IfsMethods()
{
	if (iIfsMethods == (Methods *)NULL) {
		iIfsMethods = MethodsCreate();
		iIfsMethods->create = (GeomCreateFunc *)IfsCreate;
		iIfsMethods->methods = IfsMethods;
		iIfsMethods->name = IfsName;
		iIfsMethods->intersect = IfsIntersect;
		iIfsMethods->normal = IfsNormal;
		iIfsMethods->uv = NULL;    /* IfsUV; */
		iIfsMethods->enter = IfsEnter;
		iIfsMethods->bounds = IfsBounds;
		iIfsMethods->stats = IfsStats;
		iIfsMethods->checkbounds = TRUE;
		iIfsMethods->closed = FALSE;
	}

	return iIfsMethods;
}

Methods *
CIfsMethods()
{
	if (iCIfsMethods == (Methods *)NULL) {
		iCIfsMethods = MethodsCreate();
		iCIfsMethods->create = (GeomCreateFunc *)IfsCreate;
		iCIfsMethods->methods = CIfsMethods;
		iCIfsMethods->name = CIfsName;
		iCIfsMethods->intersect = CIfsIntersect;
		iCIfsMethods->bounds = IfsBounds;
		iCIfsMethods->stats = CIfsStats;
		iCIfsMethods->checkbounds = TRUE;
		iCIfsMethods->closed = TRUE;    /* depends on its condensation */
	}

	return iCIfsMethods;
}

/* Transformed Extends Lists */
static TEList *TENew()
{
        return (TEList *)Malloc(sizeof(TEList));
}

static void TEDispose(p)
TEList *p;
{
        free(p);
}

static void TEDisposeList(list)
TEList *list;
{
        TEList *p;

/* de-allocate nodes in reverse order of allocation */
        for (p=list; p->next; p=p->next) {}
        list=p;
        while (list) {
                p = list->prev;
                TEDispose(list);
                list = p;
        }
}

/* ray - sphere intersection */
static int IfsIntersectSphere(sph, ray, mindist, maxdist)
struct IfsExtent *sph;
Ray *ray;
Float mindist, *maxdist;
{
	Float xadj, yadj, zadj;
	Float b, d, t, s;

	xadj = sph->c.x - ray->pos.x;
	yadj = sph->c.y - ray->pos.y;
	zadj = sph->c.z - ray->pos.z;
        d = xadj * xadj + yadj * yadj + zadj * zadj - sph->rsq;

	b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
	t = b * b - d;
	if (t < 0.)
		return FALSE;
	t = (Float)sqrt((double)t);
	s = b - t;
	if (s > mindist) {
		if (s < *maxdist) {
			*maxdist = s;
			return TRUE;
		}
		return FALSE;
	}
	s = b + t;
	if (s > mindist && s < *maxdist) {
		*maxdist = s;
		return TRUE;
	}
	return FALSE;
}

static int IfsExtentIntersect(ifs, ray, mindist, maxdist)
Ifs *ifs;
Ray *ray;
Float mindist, *maxdist;
{
        Vector vtmp;

        VecAddScaled(ray->pos, mindist, ray->dir, &vtmp);
        if (ifs->boxbounding) {       /*  test against solid box */
                if (OutOfBounds(&vtmp, ifs->bounds))
                        return BoundsIntersect(ray, ifs->bounds, mindist, maxdist);
        } else {                       /* test against solid sphere */
                VecSub(vtmp, ifs->extent.c, &vtmp);
                if (vtmp.x*vtmp.x + vtmp.y*vtmp.y + vtmp.z*vtmp.z > ifs->extent.rsq)
                        return IfsIntersectSphere(&ifs->extent, ray, mindist, maxdist);
        } 
        *maxdist = mindist; 
        return TRUE;
}

/* for IFS without condensation */
static void IfsExtentNormal(ifs, ray, dist, normal)
Ifs *ifs;
Ray *ray;
Float dist;
Vector *normal;	
{
	Vector pos;

        pos.x = ray->pos.x + dist * ray->dir.x;
	pos.y = ray->pos.y + dist * ray->dir.y;
        pos.z = ray->pos.z + dist * ray->dir.z;

	if (ifs->boxbounding) {     /* this produces nicer normals, so ... */
		normal->x = pos.x - (ifs->bounds[HIGH][X] + ifs->bounds[LOW][X])/2.;
		normal->y = pos.y - (ifs->bounds[HIGH][Y] + ifs->bounds[LOW][Y])/2.;
		normal->z = pos.z - (ifs->bounds[HIGH][Z] + ifs->bounds[LOW][Z])/2.;
		VecNormalize(normal); 
	} else {
        	normal->x = (pos.x - ifs->extent.c.x) / ifs->extent.r;
	        normal->y = (pos.y - ifs->extent.c.y) / ifs->extent.r;
        	normal->z = (pos.z - ifs->extent.c.z) / ifs->extent.r;
        }
}

static int DoIfsIntersect();

/* an IFS without condensation set is a primitive object, so no hitlist here */
int 
IfsIntersect(ifs, ray, mindist, maxdist)
Ifs *ifs;
Ray *ray;
Float mindist, *maxdist;
{
        IfsTests++;
        if (DoIfsIntersect(ifs, ray, (HitList *)NULL, mindist, maxdist)) {
                IfsHits++;
                return TRUE;
        }
        return FALSE;
}

/* an IFS with condensation is like an aggregate */
int 
CIfsIntersect(ifs, ray, hitlist, mindist, maxdist)
Ifs *ifs;
Ray *ray;
HitList *hitlist;
Float mindist, *maxdist;
{
        CIfsTests++;
        if (DoIfsIntersect(ifs, ray, hitlist, mindist, maxdist)) {
                CIfsHits++;
                return TRUE;
        }
        return FALSE;
}

/*
 * Ray/ifs intersection test.
 */
static int DoIfsIntersect(ifs, ray, hitlist, mindist, maxdist)
Ifs *ifs;
Ray *ray;
HitList *hitlist;
Float mindist, *maxdist;
{
        TEList *list, *closest, *p, *q, *new;
        Float nmindist, nmaxdist, distfac, pixelsize;
        Ray newray;
        Vector normal;
        int n, intersection_found, construct_children;
        Float lambda[3];
        Geom *obj;
	HitNode *np;
	IfsMap *map;

/*        if (ifs->animated && !equal(ifs->time, ray->time)) {
                IfsResolveAssoc(ifs);
                IfsComputeExtent(ifs);
                ifs->time = ray->time;
        }
*/
/* does the ray intersect the IFS's extent ? */
        distfac = 1.;
        newray = *ray;
        nmaxdist = *maxdist * distfac;
        if (!IfsExtentIntersect(ifs, &newray, EPSILON, &nmaxdist))
                return FALSE;

/* put the extent on the list */
        new = TENew();
        MatrixInit(&new->matrix);
        MatrixInit(&new->imatrix);
        new->ray = newray;
        new->distfac = distfac;
        new->dist = nmaxdist / distfac;
        new->size = ifs->size;
        new->depth = 0;
        new->rang = -99999999;
        if (!ifs->condensation) {
                IfsExtentNormal(ifs, &newray, nmaxdist, &normal);
                NormalTransform(&normal, &new->imatrix);
                VecScale((*ifs->normalweight)(ifs->size, new->size), normal, &new->normal);
        }
        new->next = new->prev = (TEList *)NULL;
        list = new;
        intersection_found = FALSE;

/* while list not empty */
        while (list) {
/* find closest entry in list */
                closest = list;
                for (p = list->next; p; p = p->next) 
                        if (p->dist < closest->dist) closest = p;

                construct_children = TRUE;
                pixelsize = PixelSize(ray, closest->dist);

/* does the (inverse transformed) ray intersect the condensation set (if any) ? */
                if (ifs->condensation) {
/* Oppenheimer trees: if TE small enough, than test against leave, not against condensation */
                        obj = ifs->condensation;
			if (ifs->condensize * closest->size / ifs->size < ifs->options.size ||
			    closest->depth >= ifs->options.depth) {
				if (ifs->leaves) obj = ifs->leaves;
				construct_children = FALSE;
			}

/* instantiate the object */
                        nmindist = mindist * closest->distfac;
                        nmaxdist = *maxdist * closest->distfac;
                        if (intersect(obj, &closest->ray, hitlist, nmindist, &nmaxdist)) {
                                *maxdist = nmaxdist / closest->distfac;
                                intersection_found = TRUE;

if (*maxdist <= closest->dist) fprintf(stderr, "impossible situation: closest = %p, closest->dist = %lf\n", closest, closest->dist);
/*				if (closest->depth > 0) 
				  	MatrixInvert(&closest->matrix, &closest->imatrix); 
*/
				np = &hitlist->data[hitlist->nodes-1];
				if (np->dotrans) {
					MatrixMult(&closest->matrix, &np->trans.trans, &np->trans.trans);
					MatrixMult(&np->trans.itrans, &closest->imatrix, &np->trans.itrans);
				} else {
		                        MatrixCopy(&closest->matrix, &np->trans.trans);
                		        MatrixCopy(&closest->imatrix, &np->trans.itrans);
				}
				np->dotrans = TRUE;

/* unless we don't need the *closest* intersection (shadow-rays!!) 
 * we cannot stop here because the generated hierarchy of extents is
 * *not* well-behaved: it is possible that we find a closer intersection
 * later.  */
                                if (ray->type == SHADOW_RAY) {
                                        TEDisposeList(list);
                                        return TRUE;
                                }

/* But we can remove all the transformed extents which are behind this one. Note
 * that at least one TE (closest) stays in the list. */
                                 for (p = list; p;) {
                                        if (p->dist > *maxdist + EPSILON) {
                                                if (p->next) p->next->prev = p->prev;
                                                if (p->prev) p->prev->next = p->next;
                                                else list = p->next;      
                                                q = p->next;
                                                TEDispose(p);
                                                p = q;
                                        } else
                                                p = p->next;
                                }
                        }
                }


/* or is the (transformed) extent smaller than one pixel ? */
                if (closest->size < pixelsize || closest->size < ifs->options.size || closest->depth >= ifs->options.depth) {
                        if (!ifs->condensation &&
                            closest->dist >= mindist && closest->dist <= *maxdist) {
                                intersection_found = TRUE;
                                *maxdist = closest->dist;
/* remember hierarchically computed normal and intersection position */
                                ifs->normal = closest->normal;
                                VecNormalize(&ifs->normal);
                                VecAddScaled(ray->pos, *maxdist, ray->dir, &ifs->intersectpos);
/* and we can stop here */
                                TEDisposeList(list);
                                return TRUE;
                        }

/* anyway, we don't construct the children of the transformed extent */
                        construct_children = FALSE;
                }


/* construct the (transformed) extent's "children" and add any that intersect to the list */
               if (construct_children) for (n=0, map=ifs->maps; n<ifs->nrmaps; n++, map++) {
/* for hierarchical IFS: use only maps of equal or higher rang. The maps are sorted to
 * descending rang number */
                        if (map->rang < closest->rang) break;
                        distfac = closest->distfac;
                        newray = closest->ray;
                        distfac *= RayTransform(&newray, &map->imatrix);
                        nmaxdist = *maxdist * distfac;

                        if (IfsExtentIntersect(ifs, &newray, EPSILON*distfac, &nmaxdist)) {
/* add the "transformed extent" to the list */
                                new = TENew();
                                new->ray = newray;
                                new->distfac = distfac;
                                new->dist = nmaxdist / distfac;
                                new->depth = closest->depth+1;
                                new->rang = map->rang;
/* pre-multiplication !! */
				if (ifs->condensation || !ifs->uniform_maps)
	                                MatrixMult(&map->matrix, &closest->matrix, &new->matrix);
                                MatrixMult(&closest->imatrix, &map->imatrix, &new->imatrix);

/* compute size of transformed extent */
                                if (map->uniform)
                                        new->size = closest->size * map->contractivity;
                                else {
                                        MatrixScaling(&new->matrix, lambda);
                                        new->size = ifs->size * Lipschitz(lambda);
                                }

/* if IFS without condensation, compute hierarchically composed normal */
                                if (!ifs->condensation) {
                                        IfsExtentNormal(ifs, &newray, nmaxdist, &normal);
                                        NormalTransform(&normal, &new->imatrix);
                                        VecAddScaled(closest->normal,
                                                     (*ifs->normalweight)(ifs->size, new->size),
                                                     normal,
                                                     &new->normal);
                                }
/* add to list */

                                new->prev = list->prev; list->prev = new;
                                new->next = list; list = new;
                        }
                }

/* remove closest from list */
                if (closest->next) closest->next->prev = closest->prev;
                if (closest->prev) closest->prev->next = closest->next;
                else list = closest->next;      /* closest was the first item */
                TEDispose(closest);
        }

        return intersection_found;
}

#define VecEqual(a, b)  ( (equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z)) )

/*
 * Compute normal to ifs at position of last intersection found.
 */
int
IfsNormal(ifs, pos, nrm, gnrm)
Ifs *ifs;
Vector *pos, *nrm, *gnrm;
{
/* this situation does not happen in rayshade.4.0.6 */
        if (!VecEqual(*pos, ifs->intersectpos))
                RLerror(RL_ABORT, "Can't compute ifs normal in other point than last intersection.\n");

/* normal is the hierarchically computed normal */
        if (ifs->normal.x == 0. && ifs->normal.y == 0. && ifs->normal.z == 0.)
                ifs->normal.z = 1.;     /* may happen when the whole IFS is smaller than one pixel */
        *nrm = *gnrm = ifs->normal;
        return FALSE;
}


/*
 * Determine if ray enters (TRUE) or leaves (FALSE) ifs at position of
 * last intersection.
 */
int
IfsEnter(ifs, ray, mind, hitd)
Ifs *ifs;
Ray *ray;
Float mind, hitd;
{
        return TRUE;
}

void
IfsBounds(ifs, bounds)
Ifs *ifs;
Float bounds[2][3];
{
        if (ifs->animated) 
                IfsResolveAssoc(ifs);

        if (ifs->animated || ifs->condensation)
                IfsComputeExtent(ifs);

        BoundsCopy(ifs->bounds, bounds);

/* produce a list of volumes that will be rendered */
	if (ifs->options.listfilename) 
		IfsList(ifs, ifs->options.listfilename);
}

char *
IfsName()
{
	return ifsName;
}

char *
CIfsName()
{
	return CifsName;
}

void
IfsStats(tests, hits)
unsigned long *tests, *hits;
{
	*tests = IfsTests;
	*hits = IfsHits;
}

void
CIfsStats(tests, hits)
unsigned long *tests, *hits;
{
	*tests = CIfsTests;
	*hits = CIfsHits;
}

void
IfsMethodRegister(meth)
UserMethodType meth;
{
	if (iIfsMethods)
		iIfsMethods->user = meth;
}