/*************************************************************************/
/*                                                                       */
/*  Copyright (c) 1994 Stanford University                               */
/*                                                                       */
/*  All rights reserved.                                                 */
/*                                                                       */
/*  Permission is given to use, copy, and modify this software for any   */
/*  non-commercial purpose as long as this copyright notice is not       */
/*  removed.  All other uses, including redistribution in whole or in    */
/*  part, are forbidden without prior written permission.                */
/*                                                                       */
/*  This software is provided with absolutely no warranty and no         */
/*  support.                                                             */
/*                                                                       */
/*************************************************************************/

/*
 * NAME
 *	hutv.c
 *
 * DESCRIPTION
 *
 */


#include <stdio.h>
#include <math.h>
#include "rt.h"

/*
 *	eps_t is a small increment in t to be added to t_in to insure that the
 *	point is inside a cell not on the boundary.
 */

REAL	eps_t	= 1.0E-10;


/*
 * NAME
 *	prn_tv_stats - print out HU traversal statistics -- not implemented in
 *		this version
 *
 * SYNOPSIS
 *	VOID	prn_tv_stats()
 *
 * RETURNS
 *	Nothing.
 */

VOID	prn_tv_stats()
	{
	}


/*
 * NAME
 *	send_ray(r, v) - send ray to remote cluster
 *
 * SYNOPSIS
 *	INT	send_ray(r, v)
 *	RAY	*r;
 *	VOXEL	*v;
 *
 * DESCRIPTION
 *	Not implemented yet.
 *
 * RETURNS
 *	Nothing, yet.
 */

INT	send_ray(r, v)
RAY	*r;
VOXEL	*v;
	{
	}


/*
 * NAME
 *	lookup_hashtable -
 *
 * SYNOPSIS
 *	VOXEL	*lookup_hashtable(indx, g)
 *	INT	indx;
 *	GRID	*g;
 *
 * DESCRIPTION
 *	Hashtable is a table of all non-empty cells in the grid.
 *	hashtable[num_buckets] is indexed by index1D mod num_buckets,
 *	num_buckets and n, the # of cells per axis, should be relatively
 *	prime, grids at different levels may have different num_buckets.  This
 *	routine assumes that there exists a non-empty voxel in the table.
 *
 * RETURNS
 *
 */

VOXEL	*lookup_hashtable(indx, g)
INT	indx;
GRID	*g;
	{
	INT	i;
	VOXEL	*v;

	i = indx - ((indx / g->num_buckets) * g->num_buckets);
	v = g->hashtable[i];

	while (v->id != indx)
		{
		v = v->next;
		if (v == NULL)
			{
			fprintf(stderr, "hashtable error \n");
			exit(-1);
			}
		}

	return (v);
	}


/*
 * NAME
 *	lookup_emptycells -
 *
 * SYNOPSIS
 *	INT	lookup_emptycells(indx, g)
 *	INT	indx;
 *	GRID	*g;
 *
 * DESCRIPTION
 *	Emptycells is a packed array of bits, one bit per cell in the grid.  A
 *	1 indicates that the cell is empty and therefore there is no entry for
 *	that cell in the hashtable.
 *
 * RETURNS
 *
 */

INT	lookup_emptycells(indx, g)
INT	indx;
GRID	*g;
	{
	INT	i, w, r, num_bits;
	UINT	p, b;

	num_bits = sizeof(UINT)*8;

	w = indx / num_bits;
	r = indx - w * num_bits;
	p = g->emptycells[w];
	b = p & (MSB >> r);
	i = b > 0 ? EMPTY : NONEMPTY;

	return (i);
	}


/*
 * NAME
 *	pop_up_a_grid -
 *
 * SYNOPSIS
 *	VOID	pop_up_a_grid(r)
 *	RAY	*r;
 *
 * DESCRIPTION
 *	Pop_up_a_grid takes RAYINFO & a grid pointer off the top of the stack
 *	and discards it which puts the ray in the voxel containing the grid
 *	being left.
 *
 * RETURNS
 *	Nothing.
 */

VOID	pop_up_a_grid(r)
RAY	*r;
	{
	RAYINFO *old_ri;

	old_ri = r->ri;
	r->ri  = old_ri->next;

	free_rayinfo(r);

	}


/*
 * NAME
 *	push_down_grid -
 *
 * SYNOPSIS
 *	VOID	push_down_grid(r, v)
 *	RAY	*r;
 *	VOXEL	*v;			/* Voxel containing the new grid.
 *
 * DESCRIPTION
 *	push_down_grid creates rayinfo for the new grid being entered and puts
 *	it on top of the stack.
 *
 * RETURNS
 *	Nothing.
 */

VOID	push_down_grid(r, v)
RAY	*r;
VOXEL	*v;
	{
	INT	n;		       /* # cells per axis in new grid.      */
	INT	small;
	INT	i_in, i_out, i, il, ih;
	REAL	wc[3];		       /* Entry point of grid in world coord.*/
	REAL	ti;
	REAL	del1, del2;
	REAL	t_in, t_out, tl, th;
	REAL	t[6], min[3];
	GRID	*new_g;
	RAYINFO *new_ri;
	RAYINFO *old_ri;

	old_ri = r->ri;
	new_g  = (GRID *)v->cell;

	new_ri = ma_rayinfo(r);
	r->ri  = new_ri;

	new_ri->next = old_ri;		/* Push new RAYINFO onto the stack.  */
	new_ri->grid = new_g;

	n = new_g->indx_inc[1];

	new_ri->delta[0] = old_ri->delta[0]/n;
	new_ri->delta[1] = old_ri->delta[1]/n;
	new_ri->delta[2] = old_ri->delta[2]/n;

	if (old_ri->t_in >= 0.0)
		{
		/* Ray origin outside the voxel. */

		new_ri->entry_plane = old_ri->entry_plane;
		new_ri->t_in = old_ri->t_in;
		ti = old_ri->t_in + eps_t;
		}
	else
		{
		/* Ray origin inside the voxel. */

		ti= 0.0;
		new_ri->entry_plane = -1;  /* no entry plane calculated */
		new_ri->t_in = old_ri->t_in;
		}

	wc[0] = ti*r->D[0] + r->P[0];
	wc[1] = ti*r->D[1] + r->P[1];
	wc[2] = ti*r->D[2] + r->P[2];


	new_ri->index3D[0] = (int)((wc[0] - new_g->min[0]) / new_g->cellsize[0]);

	if (new_ri->index3D[0] < 0)
		new_ri->index3D[0] = 0;

	if (new_ri->index3D[0] >= n)
		new_ri->index3D[0] = n - 1;


	new_ri->index3D[1] = (int)((wc[1] - new_g->min[1]) / new_g->cellsize[1]);

	if (new_ri->index3D[1] < 0)
		new_ri->index3D[1] = 0;

	if (new_ri->index3D[1] >= n)
		new_ri->index3D[1] = n - 1;


	new_ri->index3D[2] = (int)((wc[2] - new_g->min[2]) / new_g->cellsize[2]);

	if (new_ri->index3D[2] < 0)
		new_ri->index3D[2] = 0;

	if (new_ri->index3D[2] >= n)
		new_ri->index3D[2] = n - 1;


	new_ri->index1D = new_ri->index3D[0] + new_ri->index3D[1]*n + new_ri->index3D[2]*new_g->indx_inc[2];

	new_ri->indx_inc1D[0] = r->indx_inc3D[0];
	new_ri->indx_inc1D[1] = r->indx_inc3D[1]*n;
	new_ri->indx_inc1D[2] = r->indx_inc3D[2]*new_g->indx_inc[2];

	switch (new_ri->entry_plane)
		{
		case 3: /* max X */
			new_ri->entry_plane = 0;

		case 0: /* min X */
			new_ri->d[0] = new_ri->t_in + new_ri->delta[0];

			del1 = fmod((double)(old_ri->d[1] - new_ri->t_in),
			(double)new_ri->delta[1]);

			if (del1 <= eps_t)
				new_ri->d[1] = new_ri->t_in + new_ri->delta[1];
			else
				new_ri->d[1] = new_ri->t_in + del1;

			del2= fmod((double)(old_ri->d[2] - new_ri->t_in),
			(double)new_ri->delta[2]);

			if (del2 <= eps_t)
				new_ri->d[2] = new_ri->t_in + new_ri->delta[2];
			else
				new_ri->d[2] = new_ri->t_in + del2;

			small = new_ri->d[0] <= new_ri->d[1] ? 0 : 1;
			small = new_ri->d[small] <= new_ri->d[2] ? small : 2;

			new_ri->t_out = new_ri->d[small];
			new_ri->exit_plane = small;
			break;

		case 4: /* max Y */
			new_ri->entry_plane = 1;

		case 1: /* min Y */
			new_ri->d[1] = new_ri->t_in + new_ri->delta[1];

			del1 = fmod((double)(old_ri->d[0] - new_ri->t_in),
			(double)new_ri->delta[0]);

			if (del1 <= eps_t)
				new_ri->d[0] = new_ri->t_in + new_ri->delta[0];
			else
				new_ri->d[0] = new_ri->t_in + del1;

			del2 = fmod((double)(old_ri->d[2] - new_ri->t_in),
			(double)new_ri->delta[2]);

			if (del2 <= eps_t)
				new_ri->d[2] = new_ri->t_in + new_ri->delta[2];
			else
				new_ri->d[2] = new_ri->t_in + del2;

			small = new_ri->d[0] <= new_ri->d[1] ? 0 : 1;
			small = new_ri->d[small] <= new_ri->d[2] ? small : 2;

			new_ri->t_out = new_ri->d[small];
			new_ri->exit_plane = small;
			break;

		case 5: /* max Z */
			new_ri->entry_plane = 2;

		case 2: /* min Z */
			new_ri->d[2] = new_ri->t_in + new_ri->delta[2];

			del1 = fmod((double)(old_ri->d[0] - new_ri->t_in),
			(double)new_ri->delta[0]);

			if (del1 <= eps_t)
				new_ri->d[0] = new_ri->t_in + new_ri->delta[0];
			else
				new_ri->d[0] = new_ri->t_in + del1;

			del2 = fmod((double)(old_ri->d[1] - new_ri->t_in),
			(double)new_ri->delta[1]);

			if (del2 <= eps_t)
				new_ri->d[1] = new_ri->t_in + new_ri->delta[1];
			else
				new_ri->d[1] = new_ri->t_in + del2;

			small = new_ri->d[0] <= new_ri->d[1] ? 0 : 1;
			small = new_ri->d[small] <= new_ri->d[2] ? small : 2;

			new_ri->t_out = new_ri->d[small];
			new_ri->exit_plane = small;
			break;

		case -1: /* No entry plane calculated, origin inside voxel. */

			min[0] = new_g->min[0] + new_ri->index3D[0]*new_g->cellsize[0];
			min[1] = new_g->min[1] + new_ri->index3D[1]*new_g->cellsize[1];
			min[2] = new_g->min[2] + new_ri->index3D[2]*new_g->cellsize[2];

			if (r->D[0] == 0.0)
				{
				if (r->P[0] >= min[0] && r->P[0] <=  min[0] + new_g->cellsize[0])
					t[0] = -HUGE_REAL;
				else
					t[0] =	HUGE_REAL;
				}
			else
				t[0] = (min[0] - r->P[0]) / r->D[0];


			if (r->D[1] == 0.0)
				{
				if (r->P[1] >= min[1] && r->P[1] <=  min[1] + new_g->cellsize[1])
					t[1] = -HUGE_REAL;
				else
					t[1] =	HUGE_REAL;
				}
			else
				t[1] = (min[1] - r->P[1]) / r->D[1];


			if (r->D[2] == 0.0)
				{
				if (r->P[2] >= min[2] && r->P[2] <=  min[2] + new_g->cellsize[2])
					t[2] = -HUGE_REAL;
				else
					t[2] =	HUGE_REAL;
				}
			else
				t[2] = (min[2] - r->P[2]) / r->D[2];


			if (r->D[0] == 0.0)
				{
				if (r->P[0] >= min[0] && r->P[0] <=  min[0] + new_g->cellsize[0])
					t[3] = HUGE_REAL;
				else
					t[3] = HUGE_REAL;
				}
			else
				t[3] = (min[0] + new_g->cellsize[0] - r->P[0]) / r->D[0];


			if (r->D[1] == 0.0)
				{
				if (r->P[1] >= min[1] && r->P[1] <=  min[1] + new_g->cellsize[1])
					t[4] = HUGE_REAL;
				else
					t[4] = HUGE_REAL;
				}
			else
				t[4] = (min[1] + new_g->cellsize[1] - r->P[1]) / r->D[1];

			if (r->D[2] == 0.0)
				{
				if (r->P[2] >= min[2] && r->P[2] <=  min[2] + new_g->cellsize[2])
					t[5] = HUGE_REAL;
				else
					t[5] = HUGE_REAL;
				}
			else
				t[5] = (min[2] + new_g->cellsize[2] - r->P[2]) / r->D[2];

			t_in  = -HUGE_REAL;
			i_in  = -1;
			t_out =  HUGE_REAL;
			i_out = -1;

			for (i = 0; i < 3; i++)
				{
				if (t[i] < t[i + 3])
					{
					tl = t[i];
					il = i;
					th = t[i + 3];
					ih = i + 3;
					}
				else
					{
					tl = t[i + 3];
					il = i + 3;
					th = t[i];
					ih = i;
					}

				new_ri->d[i] = th;
				if (t_in < tl)			/* max min   */
					{
					t_in = tl;
					i_in = il;
					}

				if (t_out > th) 		/* min max   */
					{
					t_out = th;
					i_out = ih;
					}
				}

			if ((t_in > t_out) || (t_out < 0.0))
				{
				fprintf(stderr, "push_down_grid: Ray origin not in voxel \n");
				exit(-1);
				}

			if (i_in > 2)
				i_in -= 3;

			if (i_out > 2)
				i_out -= 3;

			new_ri->entry_plane = i_in;
			new_ri->t_in	    = t_in;
			new_ri->t_out	    = t_out;
			new_ri->exit_plane  = i_out;
			break;
		}
	}


/*
 * NAME
 *	step_grid -
 *
 * SYNOPSIS
 *	INT	step_grid(r)
 *	RAY	*r;
 *
 * DESCRIPTION
 *	Step to next voxel on grid and return index1D, updating RAYINFO in the
 *	process, return -1 if step carries off the current grid, index1D in
 *	RAYINFO is not valid if off the grid.  Assume that ray & RAYINFO are
 *	all initialized and that the current point along the ray is in the
 *	cell indentified by index1D & index3D[3].  The corresponding grid and
 *	RAYINFO are on the top of their stacks resp.
 *
 * RETURNS
 *	Nothing.
 */

INT	step_grid(r)
RAY	*r;
	{
	INT	n;			/* # cells per axis in grid.	     */
	INT	small;			/* Index of closest cell boundary.   */
	INT	*indx;
	RAY	*ra;
	GRID	*gr;
	RAYINFO *rinfo;

	ra    = r;
	rinfo = r->ri;
	gr    = rinfo->grid;
	indx  = gr->indx_inc;
	n     = indx[1];		/*  n = r->ri->grid->indx_inc[1];    */


	rinfo->t_in = rinfo->t_out;
	rinfo->index3D[rinfo->exit_plane] += r->indx_inc3D[rinfo->exit_plane];
	rinfo->entry_plane = rinfo->exit_plane;

	if (rinfo->index3D[rinfo->exit_plane] < 0 || rinfo->index3D[rinfo->exit_plane] >= n)
		{
		/* Out of range, off the grid. */

		return (-1);
		}
	else
		{
		/* Within current grid. */

		rinfo->d[rinfo->exit_plane] += rinfo->delta[rinfo->exit_plane];
		rinfo->index1D += rinfo->indx_inc1D[rinfo->exit_plane];

		/* Update exit info. */

		small = rinfo->d[0] <= rinfo->d[1] ? 0 : 1;
		small = rinfo->d[small] <= rinfo->d[2] ? small : 2;

		rinfo->t_out	  = rinfo->d[small];
		rinfo->exit_plane = small;

		return (rinfo->index1D);
		}
	}


/*
 * NAME
 *	next_voxel -
 *
 * SYNOPSIS
 *	INT	next_voxel(r)
 *	RAY	*r;
 *
 * DESCRIPTION
 *	Next_voxel() may move up or down in the heirarchy and make appropriate
 *	adjustments to the stack.  The routine returns the index1D of the next
 *	voxel.	An index == -1 indicates the ray exited the space.  Assume
 *	that ray & RAYINFO are all initialized and that the current point
 *	along the ray is in the cell indentified by index1D & index3D[3].  The
 *	corresponding grid and RAYINFO are on the top of their stacks resp.
 *
 * RETURNS
 *
 */

INT	next_voxel(r)
RAY	*r;
	{
	INT	indx;
	GRID	*gr;
	VOXEL	*v;
	RAYINFO *rinfo;

	while ((indx = step_grid(r)) == -1)
		{
		/* Step carried off grid. */

		rinfo = r->ri;
		gr = rinfo->grid;

		if (gr->subdiv_level != 0)
			{
			pop_up_a_grid(r);
			indx = r->ri->index1D;
			}
		else
			{
			/* Top level & off grid. */
			/* Exited world space.	 */

			return (-1);
			}
		}

	return (indx);
	}


/*
 * NAME
 *	next_nonempty_voxel -
 *
 * SYNOPSIS
 *	VOXEL	*next_nonempty_voxel(r)
 *	RAY	*r;
 *
 * DESCRIPTION
 *	Next_nonempty_voxel() may move up or down in the heirarchy and make
 *	appropriate adjustments to the stack.  The routine returns a pointer
 *	to the next nonempty voxel.  A zero pointer indicates the ray exited
 *	the space.  Assume that ray & RAYINFO are all initialized and that the
 *	current point along the ray is in the cell indentified by index1D &
 *	index3D[3].  The corresponding grid and RAYINFO are on the top of
 *	their stacks resp.
 *
 * RETURNS
 *	A pointer to the next nonempty voxel.
 */

VOXEL	*next_nonempty_voxel(r)
RAY	*r;
	{
	INT	indx;
	VOXEL	*v;
	GRID	*gr;
	RAYINFO *rinfo;

	indx = next_voxel(r);
	if (indx < 0)
		return (NULL);

	rinfo = r->ri;
	gr    = rinfo->grid;

	while (lookup_emptycells(indx, gr) == EMPTY)
		{
		indx = next_voxel(r);

		if (indx < 0) {
			return (NULL);
			}

		rinfo = r->ri;
		gr = rinfo->grid;
		}

	/* Found a nonempty cell. */

	v = lookup_hashtable(indx, gr);

	return (v);
	}


/*
 * NAME
 *	next_nonempty_leaf -
 *
 * SYNOPSIS
 *	VOXEL	*next_nonempty_leaf(r , step, status)
 *	RAY	*r;
 *	INT	step;
 *	INT	*status;
 *
 *	Step = 1 implies step & move up/down in the heirarchy to a leaf node.
 *	Step = 0 implies move through the heirarchy to a nonempty leaf node,
 *	stepping to the next voxel if necessary.
 *
 * DESCRIPTION
 *	Next_nonempty_leaf steps to the next nonempty leaf node moving up or
 *	down in the heirarchy as required and makes appropriate adjustments to
 *	the stack.  If the routine is called with step = 0, then the ray must
 *	be in a nonempty voxel to begin with and will move to the first
 *	nonempty leaf without stepping if possible.
 *
 *	The routine returns a pointer to the next nonempty voxel.  A zero
 *	pointer indicates that a nonempty voxel was not found and the status
 *	word can be consulted to determine why.  If the voxel is in a remote
 *	cluster, the ray is sent to that cluster and NULL is returned.
 *
 *	Assume that ray & RAYINFO are all initialized and that the current
 *	point along the ray is in the cell indentified by index1D &
 *	index3D[3].  The corresponding RAYINFO are on the top of their stacks
 *	resp.
 *
 * RETURNS
 *	Nothing.
 */

VOXEL	*next_nonempty_leaf(r , step, status)
RAY	*r;
INT	step;
INT	*status;
	{
	INT	indx;
	GRID	*ng;
	VOXEL	*v;
	RAYINFO *rinfo;

	if (step == STEP)
		v = next_nonempty_voxel(r);
	else
		{
		/* Only used by init_ray when step == 0. */

		rinfo = r->ri;
		v = lookup_hashtable(rinfo->index1D, rinfo->grid);
		}

	if (v == NULL)
		{
		*status =  EXITED_WORLD;
		return (v);
		}

	/* Nonempty voxel either a GRID or a nonempty leaf. */

	while (v->celltype == REMOTE_GRID || v->celltype == GSM_GRID)
		{
		if (v->celltype == REMOTE_GRID)
			{
			send_ray(r, v);
			*status = SENT_TO_REMOTE;
			return (NULL);
			}

		push_down_grid(r, v);

		rinfo = r->ri;
		indx  = rinfo->index1D;

		if (lookup_emptycells(indx, rinfo->grid) != EMPTY)
			{
			v = lookup_hashtable(indx, rinfo->grid);
			if (v->celltype != REMOTE_GRID && v->celltype != GSM_GRID)
				{
				/* Nonempty leaf. */

				if (v->celltype == REMOTE_VOXEL)
					{
					send_ray(r, v);
					*status = SENT_TO_REMOTE;
					return (NULL);
					}

				return (v);
				}

			/* Nonempty grid so do another while loop. */
			}
		else
			{
			v = next_nonempty_leaf(r, STEP, status);

			return (v);
			}
		}

	return (v);
	}


/*
 * NAME
 *	init_ray -
 *
 * SYNOPSIS
 *	VOXEL	*init_ray(r, g)
 *	RAY	*r;
 *	GRID	*g;			// World grid.
 *
 * DESCRIPTION
 *	Add rayinfo to a ray & intersect with world grid and find the initial
 *	nonempty leaf voxel.
 *
 * RETURNS
 *	If the ray hits a nonempty leaf voxel a pointer to that voxel is
 *	returned, otherwise NULL is returned.
 */

VOXEL	*init_ray(r, g)
RAY	*r;
GRID	*g;
	{
	INT	status;
	INT	indx, grid_id;
	INT	i_in, i_out, i, il, ih;
	REAL	t_in, t_out, tl, th;
	REAL	t[6];
	VOXEL	*v;
	GRID	*gr;
	RAYINFO *ri;

	reset_rayinfo(r);

	if (r->D[0] == 0.0)
		{
		if (r->P[0] >= g->min[0] && r->P[0] <=	g->min[0] + g->cellsize[0])
			t[0] = -HUGE_REAL;
		else
			t[0] =	HUGE_REAL;
		}
	else
		t[0] = (g->min[0] - r->P[0]) / r->D[0];


	if (r->D[1] == 0.0)
		{
		if (r->P[1] >= g->min[1] && r->P[1] <=	g->min[1] + g->cellsize[1])
			t[1] = -HUGE_REAL;
		else
			t[1] =	HUGE_REAL;
		}
	else
		t[1] = (g->min[1] - r->P[1]) / r->D[1];


	if (r->D[2] == 0.0)
		{
		if (r->P[2] >= g->min[2] && r->P[2] <=	g->min[2] + g->cellsize[2])
			t[2] = -HUGE_REAL;
		else
			t[2] =	HUGE_REAL;
		}
	else
		t[2] = (g->min[2] - r->P[2]) / r->D[2];


	if (r->D[0] == 0.0)
		{
		if (r->P[0] >= g->min[0] && r->P[0] <=	g->min[0] + g->cellsize[0])
			t[3] = HUGE_REAL;
		else
			t[3] = HUGE_REAL;
		}
	else
		t[3] = (g->min[0] + g->cellsize[0] - r->P[0]) / r->D[0];


	if (r->D[1] == 0.0)
		{
		if (r->P[1] >= g->min[1] && r->P[1] <=	g->min[1] + g->cellsize[1])
			t[4] = HUGE_REAL;
		else
			t[4] = HUGE_REAL;
		}
	else
		t[4] = (g->min[1] + g->cellsize[1] - r->P[1]) / r->D[1];


	if (r->D[2] == 0.0)
		{
		if (r->P[2] >= g->min[2] && r->P[2] <=	g->min[2] + g->cellsize[2])
			t[5] = HUGE_REAL;
		else
			t[5] = HUGE_REAL;
		}
	else
		t[5] = (g->min[2] + g->cellsize[2] - r->P[2]) / r->D[2];


	t_in  = -HUGE_REAL;
	i_in  = -1;
	t_out =  HUGE_REAL;
	i_out = -1;

	for (i = 0 ; i < 3; i++)
		{
		if (t[i] < t[i + 3])
			{
			tl = t[i];
			il = i;
			th = t[i + 3];
			ih = i+3;
			}
		else
			{
			tl = t[i + 3];
			il = i + 3;
			th = t[i];
			ih = i;
			}

		if (t_in < tl)					/* max min   */
			{
			t_in = tl;
			i_in = il;
			}

		if (t_out > th) 				/* min max   */
			{
			t_out = th;
			i_out = ih;
			}
		}

	if (t_in >= t_out || t_out < 0.0)
		{
		return (NULL);			/* Ray missed world cube.    */
		}

	ri	 = ma_rayinfo(r);
	r->ri	 = ri;
	ri->grid = g;


	/* Store exit t[]s. */

	if (t[0] >= t[3])
		ri->d[0] = t[0];
	else
		ri->d[0] = t[3];


	if (t[1] >= t[4])
		ri->d[1] = t[1];
	else
		ri->d[1] = t[4];


	if (t[2] >= t[5])
		ri->d[2] = t[2];
	else
		ri->d[2] = t[5];


	if (i_in > 2)
		i_in  -= 3;

	if (i_out > 2)
		i_out -= 3;


	ri->entry_plane = i_in;
	ri->t_in	= t_in;
	ri->t_out	= t_out;
	ri->exit_plane	= i_out;

	ri->delta[0] =	r->D[0] == 0.0 ? HUGE_REAL : g->cellsize[0] / ABS(r->D[0]);
	ri->delta[1] =	r->D[1] == 0.0 ? HUGE_REAL : g->cellsize[1] / ABS(r->D[1]);
	ri->delta[2] =	r->D[2] == 0.0 ? HUGE_REAL : g->cellsize[2] / ABS(r->D[2]);

	ri->index3D[0] = 0;
	ri->index3D[1] = 0;
	ri->index3D[2] = 0;

	r->indx_inc3D[0] = r->D[0] >= 0.0 ? 1 : -1;
	r->indx_inc3D[1] = r->D[1] >= 0.0 ? 1 : -1;
	r->indx_inc3D[2] = r->D[2] >= 0.0 ? 1 : -1;

	ri->index1D = 0;

	ri->indx_inc1D[0] = r->D[0] >= 0.0 ? 1 : -1;
	ri->indx_inc1D[1] = r->D[1] >= 0.0 ? 1 : -1;
	ri->indx_inc1D[2] = r->D[2] >= 0.0 ? 1 : -1;
	ri->next = NULL;

	/*
	 *	At this point the ray is in the top grid in the hierarchy it
	 *	must now descend to a leaf cell.
	 */

	if ((v = next_nonempty_leaf(r, NO_STEP, &status)) != NULL)
		{
		ri	= r->ri;
		indx	= ri->index1D;
		gr	= ri->grid;
		grid_id = gr->id;
		}
	else
		{
		return (NULL);
		}

	return (v);
	}