/****************************************************************************/
/*                                                                          */
/*  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                         */
/*                                                                          */
/****************************************************************************/



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

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

# define debug(x)


void C_rt_volume_iso_surface( world, view, rt_info, iptr, 
			      closest_distance, shade_info  )
C_World				*world;
C_View				*view;
C_RayTraceInfo			*rt_info;
C_IntersectionList		*iptr;
float				*closest_distance;
C_ShadeInfo			*shade_info;
{
	extern C_transform_fposition(); 
	extern C_transform_fvector(); 
	extern void C_shade_point();
	extern void C_trilin_line_intersection();
	

	int			loop;
	C_Volume		*volume;
	C_Voxel_8bit		*dptr;
	float			isovalue;
	C_StepType		prev_step_type;
	float			dot_product;	
	float			t, distance;
	C_FPosition		new_ray_origin;
	C_FPosition		trilin_origin;
	C_FPosition		point_one, point_two;
	C_FPosition		midpoint;
	C_FPosition		intersection_point;
	C_FPosition		junk_point;
	C_FVector		new_ray_direction;
	float			current_intersection;
	float			num_steps;
	int			xstep, ystep, zstep,
				found_intersection,
				voxel_x, voxel_y, voxel_z;
	float			tmax_x, tmax_y, tmax_z,
				tdelta_x, tdelta_y, tdelta_z;
	int			tstep_x, tstep_y, tstep_z;
	float			diffx, diffy, diffz;

	register float		A, B, C, D, E, F, G, H;

	C_LineIntersectInfo 	intersect_solution;

	float
				x_weight1, y_weight1, z_weight1,
				x_weight2, y_weight2, z_weight2,
				mid_x_weight, mid_y_weight, mid_z_weight,
				midpointvalue,
				value1, value2;
	C_ScalarData_8bit	*scan_data;

	int			start_in_range;
	int			fudge_steps;


	prev_step_type = iptr->start_step;

	volume = world->volume[iptr->volume];

	scan_data = volume->data.scalar_data_8bit;

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

	isovalue = scan_data->seg.isovalue;
	
	xstep = 1;
	ystep = volume->x_size_voxels;
	zstep = volume->x_size_voxels * volume->y_size_voxels;
		
        voxel_x = C_Floor( iptr->position.x );
        voxel_y = C_Floor( iptr->position.y );
        voxel_z = C_Floor( iptr->position.z );


	if ( voxel_x >= volume->x_size_voxels - 1 ) 
	{
		diffx = ( volume->x_size_voxels - 1.01 ) - iptr->position.x;

		if ( new_ray_direction.x != 0.0 )
			t = diffx / new_ray_direction.x;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
	if ( voxel_y >= volume->y_size_voxels - 1 ) 
	{
		diffy = ( volume->y_size_voxels - 1.01 ) - iptr->position.y;

		if ( new_ray_direction.y != 0.0 )
			t = diffy / new_ray_direction.y;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
	if ( voxel_z >= volume->z_size_voxels - 1 ) 
	{
		diffz = ( volume->z_size_voxels - 1.01 ) - iptr->position.z;

		if ( new_ray_direction.z != 0.0 )
			t = diffz / new_ray_direction.z;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
	if ( voxel_x < 0 ) 
	{
		diffx = -iptr->position.x + 0.01;

		if ( new_ray_direction.x != 0.0 )
			t = diffx / new_ray_direction.x;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
	if ( voxel_y < 0 ) 
	{
		diffy = -iptr->position.y + 0.01;

		if ( new_ray_direction.y != 0.0 )
			t = diffy / new_ray_direction.y;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
	if ( voxel_z < 0 ) 
	{
		diffz = -iptr->position.z + 0.01;

		if ( new_ray_direction.z != 0.0 )
			t = diffz / new_ray_direction.z;
		else
			t = -1.0;

		if ( t > 0.0 )
		{
			iptr->position.x += new_ray_direction.x * t;
			iptr->position.y += new_ray_direction.y * t;
			iptr->position.z += new_ray_direction.z * t;

        		voxel_x = C_Floor( iptr->position.x );
        		voxel_y = C_Floor( iptr->position.y );
        		voxel_z = C_Floor( iptr->position.z );
		}
	}
		
		
		
	if (   voxel_x >= volume->x_size_voxels - 1 ||
	       voxel_y >= volume->y_size_voxels - 1 ||	
	       voxel_z >= volume->z_size_voxels - 1 ||
	       voxel_x < 0 || voxel_y < 0 || voxel_z < 0 )
	{
		found_intersection = TRUE;
	}
	else
	  	found_intersection = FALSE;

	tstep_x = C_Sign( new_ray_direction.x );
	tstep_y = C_Sign( new_ray_direction.y );
	tstep_z = C_Sign( new_ray_direction.z );


	if (new_ray_direction.x != 0.0)
	{
	  	tmax_x = fabs((float)((voxel_x+(tstep_x==1)) - 
		     iptr->position.x) / new_ray_direction.x);
	  	tdelta_x = fabs(1.0 / new_ray_direction.x);
  	}
	else
  	{
	    	tmax_x = C_BIG_FLOAT;
	    	tdelta_x = C_BIG_FLOAT;
	}

	if (new_ray_direction.y != 0.0)
	{
	  	tmax_y = fabs((float)((voxel_y+(tstep_y==1)) - 
		     iptr->position.y) / new_ray_direction.y);
	  	tdelta_y = fabs(1.0 / new_ray_direction.y);
	}
	else
	{
	  	tmax_y = C_BIG_FLOAT;
	  	tdelta_y = C_BIG_FLOAT;
	}
	
	if (new_ray_direction.z != 0.0)
	{
	  	tmax_z = fabs((float)((voxel_z+(tstep_z==1)) - 
		   iptr->position.z) / new_ray_direction.z);
	  	tdelta_z = fabs(1.0 / new_ray_direction.z);
	}
	else
	{
	  	tmax_z = C_BIG_FLOAT;
	  	tdelta_z = C_BIG_FLOAT;
	}
	

	if (!found_intersection)
	{
	  	dptr = scan_data->scalar_data + 
	       		voxel_x * xstep + 
	       		voxel_y * ystep + 
	       		voxel_z * zstep;

	  	A = *(dptr);
	  	B = *(dptr + xstep);
	  	C = *(dptr + ystep);
	  	D = *(dptr + xstep + ystep);
	  	E = *(dptr + zstep);
	  	F = *(dptr + xstep + zstep);
	  	G = *(dptr + ystep + zstep);
	  	H = *(dptr + xstep + ystep + zstep); 

		value1 = C_Trilin( 
			( iptr->position.x - (float) voxel_x ),
			( iptr->position.y - (float) voxel_y ),
			( iptr->position.z - (float) voxel_z ),
			A, B, C, D, E, F, G, H );

		if ( value1 > isovalue )
			start_in_range = TRUE;
		else
			start_in_range = FALSE;
	}

	if ( start_in_range && iptr->start_cut != NULL )
	{
printf("in a place I shouldn't be\n");
	  shade_info->world_distance = iptr->distance;
	  shade_info->volume_num = iptr->volume;
	  C_transform_fposition( &(iptr->position), &intersection_point,
				 &(volume->ltow_voxels) );
	  shade_info->world_point.x = intersection_point.x;
	  shade_info->world_point.y = intersection_point.y;
	  shade_info->world_point.z = intersection_point.z;
	  shade_info->local_point.x = iptr->position.x;
	  shade_info->local_point.y = iptr->position.y;
	  shade_info->local_point.z = iptr->position.z;
	  shade_info->world_normal.x = C_BIG_FLOAT;
	  shade_info->exiting_volume = FALSE;
	  shade_info->start_on_cut = TRUE;
		  
	  found_intersection = TRUE;
	}


	while ( found_intersection == FALSE )
	{
	    if ( (!start_in_range &&
		 ( A >= isovalue ||
	           B >= isovalue ||
	           C >= isovalue ||
	           D >= isovalue ||
	           E >= isovalue ||
	           F >= isovalue ||
	           G >= isovalue ||
	           H >= isovalue )) ||
 		 (start_in_range &&
		 ( A < isovalue ||
	           B < isovalue ||
	           C < isovalue ||
	           D < isovalue ||
	           E < isovalue ||
	           F < isovalue ||
	           G < isovalue ||
	           H < isovalue )))
	    {
	     found_intersection = TRUE;

	     switch ( view->image_accuracy )
	     {
	      case C_LOW_IMAGE_ACCURACY:
	      case C_MEDIUM_IMAGE_ACCURACY:
	        switch ( prev_step_type )
	        {
	          case C_X_STEP:
	  	
		    num_steps = ((float)(voxel_x + (tstep_x < 0)) 
				  - iptr->position.x) / 
				  new_ray_direction.x;
  
		    point_one.x = voxel_x + (tstep_x < 0);
		    point_one.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		    point_one.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
  
	      	    x_weight1 = point_one.x - voxel_x;
	      	    y_weight1 = point_one.y - voxel_y;
	      	    z_weight1 = point_one.z - voxel_z;
  
		    if (tstep_x < 0)
		    {
			  value1 = C_Bilin(y_weight1, z_weight1,
					   B, D, F, H );
		    }
		    else
		    {
			  value1 = C_Bilin(y_weight1, z_weight1,
					 A, C, E, G );
		    }
					     
  
		    break;
  
	          case C_Y_STEP:
  	
	  	    num_steps = ((float)(voxel_y + (tstep_y < 0)) 
				  - iptr->position.y) / 
				  new_ray_direction.y;
  
		    point_one.x = iptr->position.x + 
			  num_steps * new_ray_direction.x;
		    point_one.y = voxel_y + (tstep_y < 0);
		    point_one.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
  
	      	    x_weight1 = point_one.x - voxel_x;
	      	    y_weight1 = point_one.y - voxel_y;
	      	    z_weight1 = point_one.z - voxel_z;
  
		    if (tstep_y < 0)
		    {
			  value1 = C_Bilin(x_weight1, z_weight1,
					   C, D, G, H );
		    }
		    else
		    {
			  value1 = C_Bilin(x_weight1, z_weight1,
					   A, B, E, F );
		    }
					     
	  	
		    break;
  
	          case C_Z_STEP:
  
		    num_steps = ((float)(voxel_z + (tstep_z < 0)) 
				  - iptr->position.z) / 
				  new_ray_direction.z;
  
		    point_one.x = iptr->position.x +
			  num_steps * new_ray_direction.x;
		    point_one.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		    point_one.z = voxel_z + (tstep_z < 0);
	  	
	      	    x_weight1 = point_one.x - voxel_x;
	      	    y_weight1 = point_one.y - voxel_y;
	      	    z_weight1 = point_one.z - voxel_z;
  
		    if (tstep_z < 0)
		    {
			  value1 = C_Bilin(x_weight1, y_weight1,
					   E, F, G, H );
		    }
		    else
		    {
			  value1 = C_Bilin(x_weight1, y_weight1,
					   A, B, C, D );
		    }
					     
		    break;
  	
	          case C_UNDEFINED_STEP:
  
		    if ( new_ray_direction.x != 0.0 )
		    {
  
		      num_steps = ((float)(voxel_x + (tstep_x < 0)) 
				  - iptr->position.x) / 
				  new_ray_direction.x;
  
		      point_one.x = voxel_x + (tstep_x < 0);
		      point_one.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		      point_one.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
		    }
  
		    if ( !(C_In_Range(point_one.x,voxel_x,voxel_x+1) 
		        && C_In_Range(point_one.y,voxel_y,voxel_y+1)
		        && C_In_Range(point_one.z,voxel_z,voxel_z+1)) ||
		        (new_ray_direction.x == 0.0) )
		    {
		      if (new_ray_direction.y != 0.0 )
		      {
	  	
		        num_steps = ((float)(voxel_y + (tstep_y < 0)) 
				  - iptr->position.y) / 
				  new_ray_direction.y;
  
		        point_one.x = iptr->position.x + 
		  	  num_steps * new_ray_direction.x;
		        point_one.y = voxel_y + (tstep_y < 0);
		        point_one.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
		      }
  
		      if ( !(C_In_Range(point_one.x,voxel_x,voxel_x+1) 
                         && C_In_Range(point_one.y,voxel_y,voxel_y+1)
                         && C_In_Range(point_one.z,voxel_z,voxel_z+1)) ||
			  (new_ray_direction.y == 0.0) )
                      {
  
		        num_steps = ((float)(voxel_z + (tstep_z < 0)) 
				  - iptr->position.z) / 
				  new_ray_direction.z;
  
		        point_one.x = iptr->position.x +
			  num_steps * new_ray_direction.x;
		        point_one.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		        point_one.z = voxel_z + (tstep_z < 0);
  
	      	        x_weight1 = point_one.x - voxel_x;
	      	        y_weight1 = point_one.y - voxel_y;
	      	        z_weight1 = point_one.z - voxel_z;
  
		        if (tstep_z < 0)
		        {
			  value1 = C_Bilin(x_weight1, y_weight1,
					   E, F, G, H );
		        }
		        else
		        {
			  value1 = C_Bilin(x_weight1, y_weight1,
					   A, B, C, D );
		        }
					     
	  	
		      }
		      else
		      {
	      	        x_weight1 = point_one.x - voxel_x;
	      	        y_weight1 = point_one.y - voxel_y;
	      	        z_weight1 = point_one.z - voxel_z;
  
		        if (tstep_y < 0)
		        {
			  value1 = C_Bilin(x_weight1, z_weight1,
					   C, D, G, H );
		        }
		        else
		        {
			  value1 = C_Bilin(x_weight1, z_weight1,
					   A, B, E, F );
		        }
					     
		      }
		    }
		    else
		    {
	      	      x_weight1 = point_one.x - voxel_x;
	      	      y_weight1 = point_one.y - voxel_y;
	      	      z_weight1 = point_one.z - voxel_z;
  
		      if (tstep_x < 0)
		      {
			  value1 = C_Bilin(y_weight1, z_weight1,
					   B, D, F, H );
		      }
		      else
		      {
			  value1 = C_Bilin(y_weight1, z_weight1,
					   A, C, E, G );
		      }
					     
		    }
		    break;
	  	
	        }
  
	        if ( tmax_x <= tmax_y && 
		     tmax_x <= tmax_z )
	        {
		    num_steps = ((float)(voxel_x + (tstep_x > 0)) 
				  - iptr->position.x) / 
				  new_ray_direction.x;
  
		    point_two.x = voxel_x + (tstep_x > 0);
		    point_two.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		    point_two.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
  
	      	    x_weight2 = point_two.x - voxel_x;
	      	    y_weight2 = point_two.y - voxel_y;
	      	    z_weight2 = point_two.z - voxel_z;
  
		    if (tstep_x > 0)
		    {
			  value2 = C_Bilin(y_weight2, z_weight2,
					   B, D, F, H );
		    }
		    else
		    {
			  value2 = C_Bilin(y_weight2, z_weight2,
					   A, C, E, G );
		    }
					     
	  	
	        }
	        else if ( tmax_y <= tmax_x &&
		          tmax_y <= tmax_z )
	        {
		    num_steps = ((float)(voxel_y + (tstep_y > 0)) 
				  - iptr->position.y) / 
				  new_ray_direction.y;
  
		    point_two.x = iptr->position.x + 
			  num_steps * new_ray_direction.x;
		    point_two.y = voxel_y + (tstep_y > 0);
		    point_two.z = iptr->position.z +
			  num_steps * new_ray_direction.z;
  
	      	    x_weight2 = point_two.x - voxel_x;
	      	    y_weight2 = point_two.y - voxel_y;
	      	    z_weight2 = point_two.z - voxel_z;
  
		    if (tstep_y > 0)
		    {
			  value2 = C_Bilin(x_weight2, z_weight2,
					   C, D, G, H );
		    }
		    else
		    {
			  value2 = C_Bilin(x_weight2, z_weight2,
					   A, B, E, F );
		    }
					     
	  	
	        }
	        else
	        {
		    num_steps = ((float)(voxel_z + (tstep_z > 0)) 
				  - iptr->position.z) / 
				  new_ray_direction.z;
  
		    point_two.x = iptr->position.x +
			  num_steps * new_ray_direction.x;
		    point_two.y = iptr->position.y + 
			  num_steps * new_ray_direction.y;
		    point_two.z = voxel_z + (tstep_z > 0);
	  	
	      	    x_weight2 = point_two.x - voxel_x;
	      	    y_weight2 = point_two.y - voxel_y;
	      	    z_weight2 = point_two.z - voxel_z;
  
		    if (tstep_z > 0)
		    {
			  value2 = C_Bilin(x_weight2, y_weight2,
					   E, F, G, H );
		    }
		    else
		    {
			  value2 = C_Bilin(x_weight2, y_weight2,
					   A, B, C, D );
		    }
					     
	        }
  
  
	        if ( start_in_range )
	        {
		  t      = value1;
		  value1 = value2;
		  value2 = t;
  
		  t           = point_one.x;
		  point_one.x = point_two.x;
		  point_two.x = t;
  
		  t           = point_one.y;
		  point_one.y = point_two.y;
		  point_two.y = t;
  
		  t           = point_one.z;
		  point_one.z = point_two.z;
		  point_two.z = t;
	        }
		
  
	        if (  value1 < isovalue && 
		     value2 >= isovalue )  
	        {
	          if (( point_one.x != point_two.x ||
		        point_one.y != point_two.y ||
		        point_one.z != point_two.z ) )
	          {
if ( view->image_accuracy == C_LOW_IMAGE_ACCURACY ||
     view->image_accuracy == C_MEDIUM_IMAGE_ACCURACY )
{
 if ( value2 != value1)
 {
  t = (float)( isovalue - value1 ) / (float)(value2 - value1); 
  point_two.x = point_one.x + t * new_ray_direction.x;
  point_two.y = point_one.y + t * new_ray_direction.y;
  point_two.z = point_one.z + t * new_ray_direction.z;
 }
  
}
else
{
		    for (loop = 0; loop < C_BINSEARCH_LIMIT; loop++)
		    {
		      midpoint.x = ( point_one.x + point_two.x )/2.0;
		      midpoint.y = ( point_one.y + point_two.y )/2.0;
		      midpoint.z = ( point_one.z + point_two.z )/2.0;
		
		      mid_x_weight = ( x_weight1 + x_weight2)/2.0;
		      mid_y_weight = ( y_weight1 + y_weight2)/2.0;
		      mid_z_weight = ( z_weight1 + z_weight2)/2.0;

		      midpointvalue = C_Trilin( 
			mid_x_weight, mid_y_weight, mid_z_weight,
			A, B, C, D, E, F, G, H );

		      if ( midpointvalue >= isovalue )
		      {
			point_two.x = midpoint.x;
			point_two.y = midpoint.y;
			point_two.z = midpoint.z;
			value2 = midpointvalue;
			x_weight2 = mid_x_weight;
			y_weight2 = mid_y_weight;
			z_weight2 = mid_z_weight;
		      }
		      else
		      {
			point_one.x = midpoint.x;
			point_one.y = midpoint.y;
			point_one.z = midpoint.z;
			value1 = midpointvalue;
			x_weight1 = mid_x_weight;
			y_weight1 = mid_y_weight;
			z_weight1 = mid_z_weight;
		      }
		    }
}
	          }
		}
		else
			found_intersection = FALSE;

		
		break;

	       case C_HIGH_IMAGE_ACCURACY:

		trilin_origin.x = new_ray_origin.x - voxel_x;
		trilin_origin.y = new_ray_origin.y - voxel_y;
		trilin_origin.z = new_ray_origin.z - voxel_z;

		C_trilin_line_intersection( &trilin_origin, 
			&new_ray_direction, (double)A, (double)B, (double)C, 
			(double)D, (double)E, (double)F, (double)G, (double)H, 
			(double)isovalue, &intersect_solution );

/***
printf("\n\nlooking for intersection!!!\n");
printf("intersection was found (search) at: %f %f %f\n", point_two.x,
	point_two.y, point_two.z );
printf("   origin: %f %f %f\n", 
	trilin_origin.x, trilin_origin.y, trilin_origin.z);
printf("direction: %f %f %f\n",
	new_ray_direction.x, new_ray_direction.y, new_ray_direction.z );
printf("A-H: %f %f %f %f %f %f %f %f\n", A, B, C, D, E, F, G, H );
printf("threshold: %d\n", isovalue );
****/

		if ( intersect_solution.num_intersections > 0 )
		{
		    for (loop=0;  loop<intersect_solution.num_intersections;
			 loop++)
		    {
			point_two.x = 
				intersect_solution.local_position[loop].x +
				voxel_x; 
			point_two.y = 
				intersect_solution.local_position[loop].y +
				voxel_y;
			point_two.z = 
				intersect_solution.local_position[loop].z +
				voxel_z;

/****
printf("intersection occurred at: %f %f %f\n",
	point_two.x, point_two.y, point_two.z );
printf("voxel location: %d %d %d\n", voxel_x, voxel_y, voxel_z );
printf("Real Solution: %f %f %f\n", 
	intersect_solution.local_position[loop].x,
	intersect_solution.local_position[loop].y,
	intersect_solution.local_position[loop].z );
****/

			if ((C_In_Range(point_two.x, voxel_x, voxel_x + 1)) &&
			    (C_In_Range(point_two.y, voxel_y, voxel_y + 1)) &&
			    (C_In_Range(point_two.z, voxel_z, voxel_z + 1)))
				break;
/***
			if ((C_In_Range(point_two.x, 0, volume->x_size_voxels)) &&
			    (C_In_Range(point_two.y, 0,volume->y_size_voxels )) &&
			    (C_In_Range(point_two.z, 0,volume->z_size_voxels )))
				break;
****/
		    } 


		    if ( loop >= intersect_solution.num_intersections )
			found_intersection = FALSE;
/*
else
printf("/n/n********************intersection occurred at: %f %f %f\n",
	point_two.x, point_two.y, point_two.z );
*/

		}
		else
			found_intersection = FALSE;

		break;

	      }
	      
	      if ( found_intersection )
	      {

	        C_transform_fposition( &point_two, 
			&intersection_point, &(volume->ltow_voxels));

		diffx = intersection_point.x - rt_info->ray_origin.x;
		diffy = intersection_point.y - rt_info->ray_origin.y;
		diffz = intersection_point.z - rt_info->ray_origin.z;

	        current_intersection = (float) sqrt ( 
		    C_Square( diffx ) +
		    C_Square( diffy ) +
		    C_Square( diffz ) ); 

		if ((rt_info->ray_type != C_PRIMARY_RAY &&
		     rt_info->ray_type != C_PARC_RAY) && 
		     (prev_step_type == C_UNDEFINED_STEP ||
		      current_intersection < 
		     (volume->x_size_units/(volume->x_size_voxels-1)))) 
		{
			current_intersection = C_BIG_FLOAT;
			if ( start_in_range )
				start_in_range = FALSE;
			else
				start_in_range = TRUE;

			found_intersection = FALSE;
		}
		else
		{
			found_intersection = TRUE;
		}

	        if ( found_intersection &&
		     current_intersection < *closest_distance &&
		     current_intersection > iptr->distance &&
		     current_intersection < iptr->max_distance )
	        {

	
		  *closest_distance = 
			current_intersection;
		  
		  shade_info->world_distance = current_intersection;
		  shade_info->volume_num = iptr->volume;
		  shade_info->world_point.x = intersection_point.x;
		  shade_info->world_point.y = intersection_point.y;
		  shade_info->world_point.z = intersection_point.z;
		  shade_info->local_point.x = point_two.x;
		  shade_info->local_point.y = point_two.y;
		  shade_info->local_point.z = point_two.z;
		  shade_info->world_normal.x = C_BIG_FLOAT;
		  if ( start_in_range )
		  	shade_info->exiting_volume = TRUE;
		  else
		  	shade_info->exiting_volume = FALSE;
	  	  shade_info->start_on_cut = FALSE;
		  


	        }
		else if ( found_intersection &&
			  current_intersection <= iptr->distance )
		{
			found_intersection = FALSE;
			if ( start_in_range )
                                start_in_range = FALSE;
                        else
                                start_in_range = TRUE;	
		}
			 
	      }	
	    
	    }
	    if ( found_intersection == FALSE )
	    {	
	      if (tmax_x < tmax_y)
	      {
	        if (tmax_x < tmax_z)
	        {
		  voxel_x += tstep_x;
		  if (voxel_x < 0 || 
		      voxel_x >= volume->x_size_voxels-1)
		  {
		    found_intersection = TRUE;
		  }
		  else
		  {
		    tmax_x += tdelta_x;
		    dptr += tstep_x * xstep;
		    prev_step_type = C_X_STEP;
		    if (tstep_x > 0)
		    {
			A = B;
			C = D;
		  	E = F;
			G = H;
			B = *(dptr + xstep);
			D = *(dptr + xstep + ystep);
			F = *(dptr + xstep + zstep);
			H = *(dptr + xstep + ystep + zstep);
		    }
		    else
		    {
			B = A;
			D = C;
			F = E;
			H = G;
			A = *(dptr);
			C = *(dptr + ystep);
			E = *(dptr + zstep);
			G = *(dptr + ystep + zstep);
		    } 
		  }
	        }
	        else
	        {
		  voxel_z += tstep_z;
		  if (voxel_z < 0 || 
		      voxel_z >= volume->z_size_voxels-1)
		  {
		    found_intersection = TRUE;
		  }
		  else
		  {
		    tmax_z += tdelta_z;
		    dptr += tstep_z * zstep;
		    prev_step_type = C_Z_STEP;
		    if (tstep_z > 0)
		    {
			A = E;
			B = F;
			C = G;
			D = H;
			E = *(dptr + zstep);
			F = *(dptr + xstep + zstep);
			G = *(dptr + ystep + zstep);
			H = *(dptr + xstep + ystep + zstep);
		    }
		    else
		    {
			E = A;
			F = B;
			G = C;
			H = D;
			A = *(dptr);
			B = *(dptr + xstep);
			C = *(dptr + ystep);
			D = *(dptr + xstep + ystep);
		    }
		  }
	        }
	      }
	      else
	      {
	        if (tmax_y < tmax_z)
	        {
		  voxel_y += tstep_y;
		  if (voxel_y < 0 || 
		      voxel_y >= volume->y_size_voxels-1)
		  {
		    found_intersection = TRUE;
		  }
		  else
		  {
		    tmax_y += tdelta_y;
		    dptr += tstep_y * ystep;
		    prev_step_type = C_Y_STEP;
		    if (tstep_y > 0)
		    {
			A = C;
			B = D;
			E = G;
			F = H;
			C = *(dptr + ystep);
			D = *(dptr + xstep + ystep);
			G = *(dptr + ystep + zstep);
			H = *(dptr + xstep + ystep + zstep);
		    }
		    else
		    {
			C = A;
			D = B;
			G = E;
			H = F;
			A = *(dptr);
			B = *(dptr + xstep);
			E = *(dptr + zstep);
			F = *(dptr + xstep + zstep);
		    }
		  }
	        }
	        else
	        {
		  voxel_z += tstep_z;
		  if (voxel_z < 0 || 
		      voxel_z >= volume->z_size_voxels-1)
		  {
		    found_intersection = TRUE;
		  }
		  else
		  {
		    tmax_z += tdelta_z;
		    dptr += tstep_z * zstep;
		    prev_step_type = C_Z_STEP;
		    if (tstep_z > 0)
		    {
			A = E;
			B = F;
			C = G;
			D = H;
			E = *(dptr + zstep);
			F = *(dptr + xstep + zstep);
			G = *(dptr + ystep + zstep);
			H = *(dptr + xstep + ystep + zstep);
		    }
		    else
		    {
			E = A;
			F = B;
			G = C;
			H = D;
			A = *(dptr);
			B = *(dptr + xstep);
			C = *(dptr + ystep);
			D = *(dptr + xstep + ystep);
		    }
		  }
	        }
	      }
	    }
	}
}


