/************************************************************************
 *									*
 *		    Copyright (c) 1988, David B. Wecker			*
 *			    All Rights Reserved				*
 *									*
 * This file is part of DBW_uRAY					*
 *									*
 * DBW_uRAY is distributed in the hope that it will be useful, but	*
 * WITHOUT ANY WARRANTY. No author or distributor accepts		*
 * responsibility to anyone for the consequences of using it or for	*
 * whether it serves any particular purpose or works at all, unless	*
 * he says so in writing. Refer to the DBW_uRAY General Public		*
 * License for full details.						*
 *									*
 * Everyone is granted permission to copy, modify and redistribute	*
 * DBW_uRAY, but only under the conditions described in the		*
 * DBW_uRAY General Public License. A copy of this license is		*
 * supposed to have been given to you along with DBW_uRAY so you	*
 * can know your rights and responsibilities. It should be in a file	*
 * named COPYING. Among other things, the copyright notice and this	*
 * notice must be preserved on all copies.				*
 ************************************************************************
 *									*
 * Authors:								*
 *	DBW - David B. Wecker						*
 *									*
 * Versions:								*
 *	V1.0 881023 DBW	- First released version			*
 *	V1.1 881110 DBW - Fixed scan coherence code			*
 *	V1.2 881125 DBW - Removed ALL scan coherence code (useless)	*
 *			  added "fat" extent boxes			*
 *									*
 ************************************************************************/

#include "uray.h"


/**************************************************************************/
/* This module builds a hierarchy of oct-tree extents to speed up tracing */
/**************************************************************************/


/* union an extent and a vector. (min of the mins and max of the maxs) */
void unionvector( re, v, nre )
EXTENT	*re,*nre;
VEC	v;
    {
    int	    i;

    for (i = 0; i < 3; i++) {
	nre->min[i] = MIN(re->min[i], v[i] );
	nre->max[i] = MAX(re->max[i], v[i] );
	}
    }

/* union 2 extents. (min of the mins and max of the maxs) */
void unionextent( e1, e2, ne )
EXTENT	*e1,*e2,*ne;
    {
    int	    i;

    for (i = 0; i < 3; i++) {
	ne->min[i] = MIN(e1->min[i], e2->min[i] );
	ne->max[i] = MAX(e1->max[i], e2->max[i] );
	}
    }

/* build all extents */
void getextent(np,e)
NODE	*np;
EXTENT	*e;
    {
    EXTENT	newe;
    EXTENT	*ext;
    SPHERE	*sph;
    RING	*flat;
    int		i;
    VEC		v1,v2;
    FLTDBL	rad;

    /* initialize the extent */
    for (i=0; i<3; i++) {
	e->min[i] =  MYHUGE;
	e->max[i] = -MYHUGE;
	}
  
    while (np) {
	switch (np->typ) {

	    /* if an extent, build all children, then return result */
	    case TYP_E:
	    ext = (EXTENT *)np;
	    getextent(ext->child,&newe);
	    for (i = 0; i < 3; i++) {
		ext->min[i] = newe.min[i];
		ext->max[i] = newe.max[i];
		}
	    break;

	    /* if a sphere, build the extent for this one object */
	    case TYP_S:
	    sph = (SPHERE *)np;
	    rad = sqrt(sph->rad);
	    for (i = 0; i < 3; i++) {
		newe.min[i] = sph->cen[i] - (rad + TOL);
		newe.max[i] = sph->cen[i] + (rad + TOL);
		}
	    break;

	    /* if a planar, build the extent for this one object */
	    case TYP_T:
	    case TYP_Q:
	    case TYP_R:
	    flat = (RING *)np;
	    for (i=0; i<3; i++) {
		newe.min[i] =  MYHUGE;
		newe.max[i] = -MYHUGE;
		}

	    /* do ring specific calculations */
	    if (np->typ == TYP_R) {
		rad = sqrt(flat->rad2);
		vcomb(rad,flat->v1,flat->p0,v1);
		vcomb(rad,flat->v2,v1,v2);
		unionvector(&newe,v2,&newe);
		vcomb(-rad,flat->v2,v1,v2);
		unionvector(&newe,v2,&newe);
		vcomb(-rad,flat->v1,flat->p0,v1);
		vcomb(rad,flat->v2,v1,v2);
		unionvector(&newe,v2,&newe);
		vcomb(-rad,flat->v2,v1,v2);
		unionvector(&newe,v2,&newe);
		}
	    else {
		vcomb(1.0,flat->p0,flat->v1,v1);
		vcomb(1.0,flat->p0,flat->v2,v2);
		unionvector(&newe,flat->p0,&newe);
		unionvector(&newe,v1,&newe);
		unionvector(&newe,v2,&newe);
		}

	    /* pad the extent (so no round off errors) */
	    for (i=0; i<3; i++) {
		newe.min[i] -= TOL;
		newe.max[i] += TOL;
		}
	    break;
	    }

	/* unionextent for next level up */
	unionextent(e,&newe,e);

	/* get next node */
	np = np->next;
	}
    }


/* this is the actual routine that builds the extent hierarchy */
/* The tree can have upto SRTMAX nodes in each bucket and can be upto */
/* SRTLEV levels deep. These have been tuned for optimal results */

#define SRTMAX	 8	/* maximum allowed in one bucket */
#define SRTLEV	16	/* sort recursion level */

NODE	*sortext(ext,lev)
int	lev;
EXTENT	*ext;
    {
    int		i,j,x,y,z;
    EXTENT	*tmp,*tmp2,*curr,*next;
    FLTDBL	cen,ave[3],min[3],max[3];
    struct {
	int	count;
	EXTENT	*bkt;
	} cube[2][2][2], *cp = &cube[0][0][0];

    /* init the 8 extents (to build) */
    for (i=0; i< 8; i++) {
	cp[i].bkt   = NULL;
	cp[i].count = 0;
	}
    for (i=0; i<3; i++) {
	min[i] =  MYHUGE;
	max[i] = -MYHUGE;
	}

    /* compute the extent range */
    curr = ext;
    while (curr) {
	for (i=0; i < 3; i++) {
	    if (curr->min[i]<min[i]) min[i] = curr->min[i];
	    if (curr->max[i]>max[i]) max[i] = curr->max[i];
	    }
	curr = (EXTENT *)curr->next;
	}

    /* now compute centers of buckets */
    for (i=0; i<3; i++) ave[i] = 0.5 * (min[i] + max[i]);

    /* now place eveybody in buckets */
    curr = ext;
    while (curr) {
	next = (EXTENT *)curr->next;
	for (i=0; i<3; i++) {
	    cen = 0.5 * (curr->max[i] + curr->min[i]);
	    if (cen >= ave[i]) j = 1;
	    else	       j = 0;
	    switch (i) {
		case 0:	x = j; break;
		case 1: y = j; break;
		case 2: z = j; break;
		}
	    }
	curr->next	    = (NODE *)cube[x][y][z].bkt;
	cube[x][y][z].bkt   = curr;
	cube[x][y][z].count++;
	curr		    = next;
	}

    /* put the buckets back together */
    curr = NULL;
    for (x=0; x<2; x++)
	for (y=0; y<2; y++)
	    for (z=0; z<2; z++) {
		if (!cube[x][y][z].count) continue;

		/* do we need to split up this bucket */
		if (cube[x][y][z].count > SRTMAX && lev < SRTLEV) {
		    cube[x][y][z].bkt = tmp = 
			(EXTENT *)sortext(cube[x][y][z].bkt,lev+1);
		    cube[x][y][z].count = 0;
		    while (tmp) {
			cube[x][y][z].count++;
			tmp = (EXTENT *)tmp->next;
			}
		    }
		if (cube[x][y][z].count != 1) {
		    tmp		= (EXTENT *)doalloc(TYP_E);
		    tmp->child	= (NODE *)cube[x][y][z].bkt;
		    for (i=0; i<3; i++) {
			tmp->min[i] =  MYHUGE;
			tmp->max[i] = -MYHUGE;
			}
		    tmp2 = (EXTENT *)tmp->child;
		    while (tmp2) {
			unionextent(tmp,tmp2,tmp);
			tmp2 = (EXTENT *)tmp2->next;
			}
		    cube[x][y][z].bkt = tmp;
		    }
		cube[x][y][z].bkt->next	= (NODE *)curr;
		curr			= cube[x][y][z].bkt;
		}

    return((NODE *)curr);
    }

/* top level extent building routine */
void doextents() {
    VEC	    scr_v;
    EXTENT  e;

    /* Get all the leafs of the extent hierarchy */
    printf("Creating object extents\n");
    getextent(nodes,&e);

    /* now sort objects (building extent tree) */
    printf("Creating extent tree\n");
    nodes = sortext(nodes,0);

    printf("Using %d extents for %d objects\n\n",extcnt,objcnt);
    }