/****************************************************************************/
/*                                                                          */
/*  VolVis is a volume visualization system for investigating, manipulating */
/*  and rendering geometric and volumetric data.                            */
/*                                                                          */
/*  Copyright (C) 1993 by the Research Foundation of the State University   */
/*                            of New York                                   */
/*                                                                          */
/*  This program is free software; you can redistribute it and/or modify    */
/*  it under the terms of the GNU General Public License as published by    */
/*  the Free Software Foundation; either version 1, or (at your option)     */
/*  any later version.                                                      */
/*                                                                          */
/*  This program is distributed in the hope that it will be useful,         */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of          */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           */
/*  GNU General Public License for more details.                            */
/*                                                                          */
/*  You should have received a copy of the GNU General Public License       */
/*  along with this program; if not, write to the Free Software             */
/*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.               */
/*                                                                          */
/*  For information on VolVis, contact us at:                               */
/*                                                                          */
/*                volvis@cs.sunysb.edu                         (email)      */
/*                                                                          */
/*                Lisa Sobierajski & Ricardo Avila             (US Mail)    */
/*                Department of Computer Science                            */
/*                State University of New York at Stony Brook               */
/*                Stony Brook, New York  11794-4400                         */
/*                                                                          */
/****************************************************************************/




/*
 *			File: C_bounding_intersect.c
 *		      Author: Lisa Sobierajski
 *			Date: 
 *		 Description: Ray-Bounding Box Intersection Routine
 *	Modification History:
 *
 *		Who?		When?		Why?
 *	--------------------------------------------------------------------
 *
 */

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

# include "C_volvis.h"
# include "C_raytrace.h"

# define debug(x) 


C_IntersectionList *C_bounding_intersection_list( 
			ray_origin, ray_direction, world, max_distance )
C_FPosition		*ray_origin;
C_FVector		*ray_direction;
C_World			*world;
float			max_distance;
{
	extern C_transform_fposition(); 

	int			vloop, loop;
	float			x_min_limit, y_min_limit, z_min_limit;
	float			x_max_limit, y_max_limit, z_max_limit;
	float			x_plane,
				y_plane,
				z_plane;
	C_Volume		*volume;
	C_IntersectionList	*intersections,
				*tempptr,
				*iptr;
	float			t, distance;
	C_FPosition		local_point;
	C_FPosition		world_point;
	C_FPosition		new_ray_origin;
	C_FVector		new_ray_direction;
	int			entered_volume;
	int			entered_from_inside;
	C_FaceType		which_x_face,
				which_y_face,
				which_z_face;
	C_CutGeom		*box_cut;

/***************************************************************************/
/****									****/
/****  Create the intersections list...					****/
/****									****/
/****  In this section of the code, an  intersection calculation is     ****/
/****  performed between the ray and each of the volumes (which includes****/ 
/****  both volume and light sources).  The following information is    ****/
/****  stored for each intersection:					****/
/****									****/
/****		position - an (x,y,z) location (in the local coordinates****/	
/****                      of the intersected volume, in voxels)        ****/
/****			   indicating the position of the intersection. ****/
/****									****/
/****		distance - the distance at which the intersection       ****/
/****			   occurred (in units, relative to the world    ****/
/****			   coordinate system).				****/
/****									****/
/****		volume   - the index of the intersected volume.		****/
/****	  								****/
/****		next_intersection - link to the next structure.		****/
/****									****/
/****	Each intersection calculation is done by first, testing if the  ****/
/****   ray.origin is inside of the volume being checked.  If not, then ****/
/****   an intersection calculation is performed for each of three      ****/
/****   possible faces of the bounding cube.  (These faces are          ****/
/****   determined according to the ray.direction - if the normal of    ****/
/****   a face is less than 90 degrees away from the ray.direction,     ****/
/****   then the ray would intersect this plane on the way out, not the ****/
/****   way in.)  If an intersection is found, then it is added to the  ****/
/****   intersections list, in sorted order (closest intersection at the****/
/****   head of the list.)						****/
/****									****/
/***************************************************************************/

	/****    initialize intersection and iptr to NULL    ****/

	intersections = iptr = NULL;


	/****    loop through every volume    ****/

	for (vloop = 0; vloop < world->num_volumes; vloop++)
	{
	  /****    for simplicity, let volume = world->volume[vloop]    ****/

	  volume = world->volume[vloop];

	  if ( volume->visible )
	  {

	    /****	   we have not yet entered this volume  ****/

	    entered_volume = FALSE;
	    entered_from_inside = FALSE;

	    switch (volume->data_type)
	    {
	      case C_SCALAR_DATA_8BIT:
			x_min_limit = 0.0;
			y_min_limit = 0.0;
			z_min_limit = 0.0;
			x_max_limit = volume->x_size_voxels - 1; 
			y_max_limit = volume->y_size_voxels - 1;
			z_max_limit = volume->z_size_voxels - 1;
			break;

	      case C_GEOMETRIC_DATA:
			x_min_limit = 0.0;
			y_min_limit = 0.0;
			z_min_limit = 0.0;
			x_max_limit = volume->x_size_voxels - 1; 
			y_max_limit = volume->y_size_voxels - 1;
			z_max_limit = volume->z_size_voxels - 1;
			break;
	      default:
			C_error_message("Internal Error: unknown data type!\n");
			return NULL; 
			break;
	    }

	    box_cut = NULL;

	    /****    modify limits due to simple box cuts    ****/
	    for ( loop = 0; loop < volume->num_cuts; loop++ )
	      if ( volume->local_cut[loop]->cut_type == C_SIMPLE_BOX_CUT )
	      {
		x_min_limit = volume->local_cut[loop]->cut.box_cut[0];
		y_min_limit = volume->local_cut[loop]->cut.box_cut[1];
		z_min_limit = volume->local_cut[loop]->cut.box_cut[2];
		x_max_limit = volume->local_cut[loop]->cut.box_cut[3];
		y_max_limit = volume->local_cut[loop]->cut.box_cut[4];
		z_max_limit = volume->local_cut[loop]->cut.box_cut[5];
		
		box_cut = volume->local_cut[loop];
		break;
	      }
	    

	    /****    transform origin, to check if it is in this volume   ****/

	    C_transform_fvector( ray_direction, ray_origin, 
				 &new_ray_direction, &new_ray_origin,
	 		         &(volume->wtol_voxels) );

	    if (  new_ray_origin.x >= x_min_limit &&
	      	  new_ray_origin.x < x_max_limit &&
		  new_ray_origin.y >= y_min_limit &&
		  new_ray_origin.y < y_max_limit &&
		  new_ray_origin.z >= z_min_limit &&
		  new_ray_origin.z < z_max_limit ) 
	    {

	      /****    the origin is in this volume, add this     ****/
	      /****    intersection to the intersections list     ****/

       	      tempptr = ( C_IntersectionList * ) 
	      		malloc( sizeof ( C_IntersectionList ) );

	      tempptr->position.x  = new_ray_origin.x;
	      tempptr->position.y  = new_ray_origin.y;
	      tempptr->position.z  = new_ray_origin.z;
	      tempptr->distance    = 0.0;
	      tempptr->volume	   = vloop;
	      tempptr->start_step  = C_UNDEFINED_STEP;
	      tempptr->start_cut   = NULL;
	      tempptr->end_cut     = NULL;

	      entered_volume = TRUE;
	      entered_from_inside = TRUE;

	      if ( intersections == NULL )
	      {
	          tempptr->next_intersection = NULL;
	          intersections  = tempptr;
	      }
	      else 
	      {
	          tempptr->next_intersection = intersections;
	          intersections = tempptr;
	      }
	    }
	    else 
	    {

	      /****    The origin was not in the volume, now we must    ****/
	      /****    intersect this ray with the bounding cube of     ****/
	      /****    this volume.				        ****/
  
	      /****    Check the dot product between the C_MINUS_X_FACE ****/
	      /****    of the bounding cube, and the ray_direction. If  ****/
	      /****    if is < 0 (the angle between them is more than   ****/
	      /****    90 degrees) then intersect with this plane.  If  ****/
	      /****    the dot product is > 0 (the angle between them   ****/
	      /****    is less than 90 degrees) then intersect this ray ****/
	      /****    with the C_PLUS_X_FACE.  If the dot product is 0,****/
	      /****    then don't intersect with either (the rays are   ****/
	      /****    parallel.)				        ****/
 
  
	      if ( new_ray_direction.x < 0 )
	      {

		/****   We want to intersect with C_PLUS_X_FACE   ****/

		x_plane = x_max_limit;

		/****   We need to negate the dot product, since  ****/
		/****   it was calculated for the C_MINUS_X_FACE. ****/

		which_x_face = C_MINUS_X_FACE;
	      }
	      else
	      {
	        x_plane = x_min_limit;
		which_x_face = C_PLUS_X_FACE;
	      }


	      if ( new_ray_direction.y < 0 )
	      {
                /****   We want to intersect with C_PLUS_Y_FACE   ****/

		y_plane = y_max_limit;

                /****   We need to negate the dot product, since  ****/
                /****   it was calculated for the C_MINUS_Y_FACE. ****/

		which_y_face = C_MINUS_Y_FACE;
	      }
	      else
	      {
	        y_plane = y_min_limit;
		which_y_face = C_PLUS_Y_FACE;
	      }


	      if ( new_ray_direction.z < 0 )
	      {
                /****   We want to intersect with C_PLUS_Z_FACE   ****/

	 	z_plane = z_max_limit;

                /****   We need to negate the dot product, since  ****/
                /****   it was calculated for the C_MINUS_Z_FACE. ****/

		which_z_face = C_MINUS_Z_FACE;
	      }
	      else
	      {
	        z_plane = z_min_limit; 
	 	which_z_face = C_PLUS_Z_FACE;
	      }
		
	      if ( new_ray_direction.x != 0 ) 
	      {
/****
printf("checking x intersection\n");
****/
		/****     Find the t at which this ray intersects    ****/
		/****     this plane.				     ****/

		t = (x_plane - new_ray_origin.x) / new_ray_direction.x;

		/****    Check if this happens at a positive t value ****/
		/****    which means that this happens forward along ****/
		/****    the ray as opposed to backwards             ****/

		if ( t >= 0.0 )
		{
/****
printf("plane, neworigin.x origin.x= %f %f %f\n", x_plane, new_ray_origin.x,
		ray_origin->x );
printf("pos distance\n");
****/
		    /****   find which point this t corresponds to   ****/

		    local_point.x = new_ray_direction.x * t + new_ray_origin.x;
		    local_point.y = new_ray_direction.y * t + new_ray_origin.y;
		    local_point.z = new_ray_direction.z * t + new_ray_origin.z;
/**
printf("t, local point: %f   %f %f %f\n", t, local_point.x,
	local_point.y, local_point.z );
*/


		    /****   check where on the plane this ray     ****/
		    /****   intersects - it must intersect on the ****/
		    /****   face of the cube.			  ****/
		    if ( local_point.y >= y_min_limit && 
		         local_point.y < y_max_limit &&
		         local_point.z >= z_min_limit && 
		         local_point.z < z_max_limit )
		    {
/**
printf("intersection x\n");
**?
			/****  It intersects on the face, now    ****/
			/****  add it to the intersections list  ****/
			/****  (if it is less than max_distance) ****/
			/****  in sorted order, with the lowest  ****/
			/****  distance at the head of the list. ****/


		        /****   transform this point to world coordinates ****/

		        C_transform_fposition( &local_point, &world_point, 
			    		       &(volume->ltow_voxels) );

			distance = (float) sqrt( (double) 
			  ( C_Square( ray_origin->x - world_point.x ) +
			    C_Square( ray_origin->y - world_point.y ) +
			    C_Square( ray_origin->z - world_point.z ) )); 

			if ( distance < max_distance )
			{

			  tempptr = ( C_IntersectionList * ) 
			      malloc( sizeof ( C_IntersectionList ) );

			  tempptr->position.x  = local_point.x;
			  tempptr->position.y  = local_point.y;
			  tempptr->position.z  = local_point.z;
			  tempptr->distance    = distance;
	      		  tempptr->volume      = vloop;
			  tempptr->start_step  = C_X_STEP;
			  tempptr->start_cut   = NULL;
	      		  tempptr->end_cut     = NULL;

			  entered_volume = TRUE;

			  if ( intersections == NULL )
			  {
			    tempptr->next_intersection = NULL;
			    intersections  = tempptr;
			  }
			  else if ( distance <= intersections->distance )
			  {
			    tempptr->next_intersection = intersections;
			    intersections = tempptr;
			  }
			  else 
			  {
			    iptr = intersections;
			    
			    while ( iptr->next_intersection != NULL &&
				    iptr->next_intersection->distance < 
				    distance )
			    {
				iptr = iptr->next_intersection;
			    }


			    tempptr->next_intersection = 
					iptr->next_intersection;
			    iptr->next_intersection = tempptr;
			  }
			}
		    }
		}

	      }
	    
	    
	      if (!entered_volume)
	      {

	        /****    Check the dot product between the C_MINUS_Y_FACE ****/
	        /****    of the bounding cube, and the ray_direction. If  ****/
	        /****    if is < 0 (the angle between them is more than   ****/
	        /****    90 degrees) then intersect with this plane.  If  ****/
	        /****    the dot product is > 0 (the angle between them   ****/
	        /****    is more than 90 degrees) then intersect this ray ****/
	        /****    with the C_PLUS_Y_FACE.  If the dot product is 0,****/
	        /****    then don't intersect with either (the rays are   ****/
	        /****    parallel.)					  ****/

	        if ( new_ray_direction.y != 0 ) 
	        {
/****
printf("checking y intersection\n");
****/
                  /****     Find the t at which this ray intersects    ****/
                  /****     this plane.                                ****/

		  t = (y_plane - new_ray_origin.y) / new_ray_direction.y;

                  /****    Check if this happens at a positive t value ****/
                  /****    which means that this happens forward along ****/
                  /****    the ray as opposed to backwards             ****/
  
		  if ( t >= 0.0 )
		  {
/****
printf("pos distance\n");
****/
                    /****   find which point this t corresponds to   ****/

		    local_point.x = new_ray_direction.x * t + new_ray_origin.x;
		    local_point.y = new_ray_direction.y * t + new_ray_origin.y;
		    local_point.z = new_ray_direction.z * t + new_ray_origin.z;



                    /****   check where on the plane this ray     ****/
                    /****   intersects - it must intersect on the ****/
                    /****   face of the cube.                     ****/

		    if ( local_point.x >= x_min_limit && 
		         local_point.x < x_max_limit &&
		         local_point.z >= z_min_limit && 
		         local_point.z < z_max_limit )
		    {
/***
printf("intersection y\n");
***/
			/****  It intersects on the face, now    ****/
			/****  add it to the intersections list  ****/
			/****  (if it is less than max_distance) ****/
			/****  in sorted order, with the lowest  ****/
			/****  distance at the head of the list. ****/

                        /****   transform this point to world coordinates ****/

		        C_transform_fposition( &local_point, &world_point, 
			    		       &(volume->ltow_voxels) );

			distance = (float) sqrt( (double) 
			  ( C_Square( ray_origin->x - world_point.x ) +
			    C_Square( ray_origin->y - world_point.y ) +
			    C_Square( ray_origin->z - world_point.z ) )); 

			if ( distance < max_distance )
			{
			  tempptr = ( C_IntersectionList * ) 
			      malloc( sizeof ( C_IntersectionList ) );

			  tempptr->position.x  = local_point.x;
			  tempptr->position.y  = local_point.y;
			  tempptr->position.z  = local_point.z;
			  tempptr->distance    = distance;
	      		  tempptr->volume      = vloop;
			  tempptr->start_step  = C_Y_STEP;
			  tempptr->start_cut   = box_cut;
	      		  tempptr->end_cut     = NULL;

			  entered_volume = TRUE;

			  if ( intersections == NULL )
			  {
			    tempptr->next_intersection = NULL;
			    intersections  = tempptr;
			  }
			  else if ( distance <= intersections->distance )
			  {
			    tempptr->next_intersection = intersections;
			    intersections = tempptr;
			  }
			  else 
			  {
			    iptr = intersections;
			    
			    while ( iptr->next_intersection != NULL &&
				    iptr->next_intersection->distance < 
				    distance )
			    {
				iptr = iptr->next_intersection;
			    }


			    tempptr->next_intersection = 
					iptr->next_intersection;
			    iptr->next_intersection = tempptr;
			  }
			}
		    }
		  }
	        }
	      
	        if (!entered_volume)
	        { 
	          /****  Check the dot product between the C_MINUS_Z_FACE ****/
	          /****  of the bounding cube, and the ray_direction. If  ****/
	          /****  if is < 0 (the angle between them is more than   ****/
	          /****  90 degrees) then intersect with this plane.  If  ****/
	          /****  the dot product is > 0 (the angle between them   ****/
	          /****  is more than 90 degrees) then intersect this ray ****/
	          /****  with the C_PLUS_Z_FACE.  If the dot product is 0,****/
	          /****  then don't intersect with either (the rays are   ****/
	          /****  parallel.)				          ****/

	          if ( new_ray_direction.z != 0 ) 
	          {
/****
printf("checking z intersection\n");
****/
                    /****     Find the t at which this ray intersects    ****/
                    /****     this plane.                                ****/

		    t = (z_plane - new_ray_origin.z) / new_ray_direction.z;

                    /****    Check if this happens at a positive t value ****/
                    /****    which means that this happens forward along ****/
                    /****    the ray as opposed to backwards             ****/

		    if ( t >= 0.0 )
		    {
/****
printf("pos distance\n");
****/
                      /****   find which point this t corresponds to   ****/
  
		      local_point.x = new_ray_direction.x*t + new_ray_origin.x;
		      local_point.y = new_ray_direction.y*t + new_ray_origin.y;
		      local_point.z = new_ray_direction.z*t + new_ray_origin.z;
  


                      /****   check where on the plane this ray     ****/
                      /****   intersects - it must intersect on the ****/
                      /****   face of the cube.                     ****/

		      if ( local_point.y >= y_min_limit && 
		           local_point.y < y_max_limit &&
		           local_point.x >= x_min_limit && 
		           local_point.x < x_max_limit )
		      {
/***
printf("intersection z\n");
***/
			/****  It intersects on the face, now    ****/
			/****  add it to the intersections list  ****/
			/****  (if it is less than max_distance) ****/
			/****  in sorted order, with the lowest  ****/
			/****  distance at the head of the list. ****/

                        /****   transform this point to local coordinates ****/
  
		        C_transform_fposition( &local_point, &world_point, 
			    		       &(volume->ltow_voxels) );

			distance = (float) sqrt( (double) 
			  ( C_Square( ray_origin->x - world_point.x ) +
			    C_Square( ray_origin->y - world_point.y ) +
			    C_Square( ray_origin->z - world_point.z ) )); 

			if ( distance < max_distance )
			{

			  tempptr = ( C_IntersectionList * ) 
			      malloc( sizeof ( C_IntersectionList ) );

			  tempptr->position.x  = local_point.x;
			  tempptr->position.y  = local_point.y;
			  tempptr->position.z  = local_point.z;
			  tempptr->distance    = distance;
	      		  tempptr->volume      = vloop;
			  tempptr->start_step  = C_Z_STEP;
			  tempptr->start_cut   = box_cut;
	      		  tempptr->end_cut     = NULL;

			  entered_volume = TRUE;

			  if ( intersections == NULL )
		  	  {
			    tempptr->next_intersection = NULL;
			    intersections  = tempptr;
			  }
			  else if ( distance <= intersections->distance )
			  {
			    tempptr->next_intersection = intersections;
			    intersections = tempptr;
			  }
			  else 
			  {
			    iptr = intersections;
			    
			    while ( iptr->next_intersection != NULL &&
				    iptr->next_intersection->distance < 
				    distance )
			    {
				iptr = iptr->next_intersection;
			    }


			    tempptr->next_intersection = 
					iptr->next_intersection;
			    iptr->next_intersection = tempptr;
			  }
			}
		      }
		    }
	          }
	        }
	      }
	    }
	      
		

	    /****   If we entered the volume, then we must exit it   ****/

	    if ( entered_volume ) 
	    {
	      /****   Loosen the limits so that we are guaranteed that  ****/
	      /****   every ENTRY will have a corresponding EXIT        ****/

	      x_min_limit -= 0.01;
	      y_min_limit -= 0.01;
	      z_min_limit -= 0.01;
	      x_max_limit += 0.01;
	      y_max_limit += 0.01;
	      z_max_limit += 0.01;

	      if ( entered_from_inside )
	      {

	        x_plane = x_min_limit;
  
	        if ( new_ray_direction.x > 0 )
	        {

		  /****   We want to intersect with C_PLUS_X_FACE   ****/

		  x_plane = x_max_limit;

		  /****   We need to negate the dot product, since  ****/
		  /****   it was calculated for the C_MINUS_X_FACE. ****/

	        }
	      } 
	      else
	      {
		if ( which_x_face == C_MINUS_X_FACE )
			x_plane = x_min_limit;
		else
			x_plane = x_max_limit;

	      }

	      if ( new_ray_direction.x != 0 ) 
	      {
		  /****     Find the t at which this ray intersects    ****/
		  /****     this plane.				       ****/

		  t = (x_plane - new_ray_origin.x) / new_ray_direction.x;

		  /****    Check if this happens at a positive t value ****/
		  /****    which means that this happens forward along ****/
		  /****    the ray as opposed to backwards             ****/

		  if ( t >= 0.0 )
		  {

/****
printf("pos distance\n");
****/
		    /****   find which point this t corresponds to   ****/

		    local_point.x = new_ray_direction.x * t + new_ray_origin.x;
		    local_point.y = new_ray_direction.y * t + new_ray_origin.y;
		    local_point.z = new_ray_direction.z * t + new_ray_origin.z;

		    /****   check where on the plane this ray     ****/
		    /****   intersects - it must intersect on the ****/
		    /****   face of the cube.			  ****/

		    if ( local_point.y >= y_min_limit && 
		         local_point.y < y_max_limit &&
		         local_point.z >= z_min_limit && 
		         local_point.z < z_max_limit )
		    {

/****
printf("intersection\n");
****/
			/****  It intersects on the face, now    ****/
			/****  add it to the intersections list  ****/
			/****  (if it is less than max_distance) ****/
			/****  in sorted order, with the lowest  ****/
			/****  distance at the head of the list. ****/

		        /****   transform this point to local coordinates ****/

		        C_transform_fposition( &local_point, &world_point, 
			    		       &(volume->ltow_voxels) );

			distance = (float) sqrt( (double) 
			  ( C_Square( ray_origin->x - world_point.x ) +
			    C_Square( ray_origin->y - world_point.y ) +
			    C_Square( ray_origin->z - world_point.z ) )); 


			if (distance > max_distance )
			{
				distance             = max_distance;
			        tempptr->start_cut   = NULL;
				
			}
			else
				tempptr->end_cut   = box_cut;

			tempptr->max_distance    = distance;
				

			entered_volume = FALSE;

		    }
		  }
	      }
	

	      if ( entered_volume && entered_from_inside )
	      {

	          y_plane = y_min_limit; 

	          if ( new_ray_direction.y > 0 )
	          {
                    /****   We want to intersect with C_PLUS_Y_FACE   ****/

		    y_plane = y_max_limit;


                    /****   We need to negate the dot product, since  ****/
                    /****   it was calculated for the C_MINUS_Y_FACE. ****/

	          }
	      }
	      else if ( entered_volume )
	      {
		  if ( which_y_face == C_MINUS_Y_FACE )
			  y_plane = y_min_limit;
		  else
			  y_plane = y_max_limit;

	      }

	      if ( entered_volume && new_ray_direction.y != 0 ) 
	      {
                  /****     Find the t at which this ray intersects    ****/
                  /****     this plane.                                ****/

		  t = (y_plane - new_ray_origin.y) / new_ray_direction.y;

                  /****    Check if this happens at a positive t value ****/
                  /****    which means that this happens forward along ****/
                  /****    the ray as opposed to backwards             ****/

		  if ( t >= 0.0 )
		  {
/****
printf("pos distance\n");
****/
                    /****   find which point this t corresponds to   ****/

		    local_point.x = new_ray_direction.x * t + new_ray_origin.x;
		    local_point.y = new_ray_direction.y * t + new_ray_origin.y;
		    local_point.z = new_ray_direction.z * t + new_ray_origin.z;



                    /****   check where on the plane this ray     ****/
                    /****   intersects - it must intersect on the ****/
                    /****   face of the cube.                     ****/

		    if ( local_point.x >= x_min_limit && 
		         local_point.x < x_max_limit &&
		         local_point.z >= z_min_limit && 
		         local_point.z < z_max_limit )
		    {
/****
printf("intersection\n");
****/
		      /****  It intersects on the face, now    ****/
		      /****  add it to the intersections list  ****/
		      /****  (if it is less than max_distance) ****/
		      /****  in sorted order, with the lowest  ****/
		      /****  distance at the head of the list. ****/

                      /****   transform this point to local coordinates ****/

		      C_transform_fposition( &local_point, &world_point, 
			    		     &(volume->ltow_voxels) );

		      distance = (float) sqrt( (double) 
		        ( C_Square( ray_origin->x - world_point.x ) +
		          C_Square( ray_origin->y - world_point.y ) +
		          C_Square( ray_origin->z - world_point.z ) )); 


		      if (distance > max_distance )
		      {
				distance         = max_distance;
				tempptr->end_cut = NULL;
				
		      }
		      else
				tempptr->end_cut = box_cut;

		      tempptr->max_distance    = distance;

		      entered_volume = FALSE;

		    }
		  }
	      }
		
		
		
	      
	      if ( entered_volume && entered_from_inside )
	      {
	        z_plane = z_min_limit; 

	        if ( new_ray_direction.z > 0 )
	        {
                      /****   We want to intersect with C_PLUS_Z_FACE   ****/
  
		      z_plane = z_max_limit;

                      /****   We need to negate the dot product, since  ****/
                      /****   it was calculated for the C_MINUS_Z_FACE. ****/

	        }
	      }
	      else if ( entered_volume )
	      {
		  if ( which_z_face == C_MINUS_Z_FACE )
			  z_plane = z_min_limit;
		  else
			  z_plane = z_max_limit;

	      }

	      if ( entered_volume && new_ray_direction.z != 0 ) 
	      {
                /****     Find the t at which this ray intersects    ****/
                /****     this plane.                                ****/

		t = (z_plane - new_ray_origin.z) / new_ray_direction.z;

                /****    Check if this happens at a positive t value ****/
                /****    which means that this happens forward along ****/
                /****    the ray as opposed to backwards             ****/

		if ( t >= 0.0 )
		{
/****
printf("pos distance\n");
****/
                  /****   find which point this t corresponds to   ****/
  
		  local_point.x = new_ray_direction.x * t + new_ray_origin.x;
		  local_point.y = new_ray_direction.y * t + new_ray_origin.y;
		  local_point.z = new_ray_direction.z * t + new_ray_origin.z;



                  /****   check where on the plane this ray     ****/
                  /****   intersects - it must intersect on the ****/
                  /****   face of the cube.                     ****/

		  if ( local_point.y >= y_min_limit && 
		       local_point.y < y_max_limit &&
		       local_point.x >= x_min_limit && 
		       local_point.x < x_max_limit )
		  {
/****
printf("intersection\n");
****/
		      /****  It intersects on the face, now    ****/
		      /****  add it to the intersections list  ****/
		      /****  (if it is less than max_distance) ****/
		      /****  in sorted order, with the lowest  ****/
		      /****  distance at the head of the list. ****/
  
                      /****   transform this point to local coordinates ****/

		      C_transform_fposition( &local_point, &world_point, 
		   		       &(volume->ltow_voxels) );

		      distance = (float) sqrt( (double) 
		        ( C_Square( ray_origin->x - world_point.x ) +
		          C_Square( ray_origin->y - world_point.y ) +
		          C_Square( ray_origin->z - world_point.z ) )); 
  
  
		      if (distance > max_distance )
		      {
				distance         = max_distance;
				tempptr->end_cut = NULL;
		      }
		      else
				tempptr->end_cut = box_cut;
  
		      tempptr->max_distance    = distance;

		      entered_volume = FALSE;

		  }
	        }
	      }
  

/***
	      if ( tempptr )
	      {
		printf("position: %f %f %f\n", tempptr->position.x,
			  tempptr->position.y, tempptr->position.z ); 
		printf("direction: %f %f %f\n", new_ray_direction.x,
			new_ray_direction.y, new_ray_direction.z );
		if ( entered_volume )
			printf("Didn't come back out!!!\n");
		else
			printf("distance = %f\n", tempptr->distance );
	      }
***/



	      if ( entered_volume )
	      {
		  printf("*************************************OH NO!!!\n");
		  printf("*************************************OH NO!!!\n");
		  printf("*************************************OH NO!!!\n");
		  tempptr->max_distance = C_BIG_FLOAT;
	      }
	    }
  
	  }

	}

	return intersections;
}
  
  
