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


extern C_IntersectionList *C_bounding_intersection_list();

	
void C_rt_volume_scalar_opacity( 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(); 

	int			exiting_volume;
	int			exited_volume;
	C_Volume		*volume;
	C_Voxel_8bit		*dptr;
	C_StepType		prev_step_type;
	float			dot_product;	
	float			t, distance;
	C_FPosition		new_ray_origin;
	C_FPosition		point_one, point_two;
	C_FPosition		intersection_point;
	C_FPosition		junk_point;
	C_FVector		new_ray_direction;
	int			xstep, ystep, zstep,
				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;

	register float		A, B, C, D, E, F, G, H;
	float 			value1, value2;
	C_ScalarData_8bit	*scan_data;


	float			delta_t;
	float			scale_distance;
	float			remaining_opacity;
	float			accum_color_red;
	float			accum_color_green;
	float			accum_color_blue;
	float			curr_opacity;
        float                   diffx, diffy, diffz;
	float			delta_distance;



	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) );


	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 )
        {
                exited_volume = TRUE;
        }
        else
                exited_volume = 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 ( !exited_volume )
	{
		*closest_distance = iptr->distance;
			  
		shade_info->world_distance = iptr->distance;
		shade_info->volume_num = iptr->volume;
		
		C_transform_fposition( &(iptr->position), 
				       &(intersection_point),
				       &(volume->ltow_voxels) );

		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;

		point_one.x = iptr->position.x + new_ray_direction.x;
		point_one.y = iptr->position.y + new_ray_direction.y;
		point_one.z = iptr->position.z + new_ray_direction.z;

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

		scale_distance = sqrt( (double) (
			C_Square(point_two.x - intersection_point.x) +
			C_Square(point_two.y - intersection_point.y) +
			C_Square(point_two.z - intersection_point.z))); 
			

	  	exited_volume = FALSE;
	  	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 = value2 = 
		    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 );

		point_one.x = point_two.x = iptr->position.x;
		point_one.y = point_two.y = iptr->position.y;
		point_one.z = point_two.z = iptr->position.z;

		distance = 0.0;

		remaining_opacity = 1.0;

		accum_color_red = 0.0;
		accum_color_green = 0.0;
		accum_color_blue = 0.0;

	}

	if ( view->image_accuracy == C_HIGH_IMAGE_ACCURACY || 
	     view->image_accuracy == C_MEDIUM_IMAGE_ACCURACY  )
	{
	  if ( view->image_accuracy == C_HIGH_IMAGE_ACCURACY )
	  {
	    distance = 0.125;
	    delta_distance = 0.25;
	  }
	  else
	  {
	    distance = 0.5;
	    delta_distance = 1.0;
	  }

	  while ( voxel_x >= 0 && voxel_x < volume->x_size_voxels - 1 && 
	          voxel_y >= 0 && voxel_y < volume->y_size_voxels - 1 && 
	          voxel_z >= 0 && voxel_z < volume->z_size_voxels - 1    )
	  {
	    value2 = 
		    C_Trilin( 
		    	    ( point_two.x - (float) voxel_x ),
			    ( point_two.y - (float) voxel_y ),
			    ( point_two.z - (float) voxel_z ),
			      A, B, C, D, E, F, G, H );

	    if ( distance > 0.0 )
	    {
		curr_opacity = 
			scan_data->seg.s_opacity->opacity[(int)value2] *
				distance;

		remaining_opacity *= curr_opacity;

		if ( remaining_opacity < 0.0 )
			curr_opacity = 1.0;

		switch ( volume->color_type )
		{
		  case C_SINGLE_COLOR:

			accum_color_red += curr_opacity * 
				(float)(
				volume->color_data.single_color->red)
				/ 255.0;
			accum_color_green += curr_opacity * 
				(float)(
				volume->color_data.single_color->green)
				/ 255.0;
			accum_color_blue += curr_opacity * 
				(float)(
				volume->color_data.single_color->blue)
				/ 255.0;
	
			break;

		  case C_COLOR_ARRAY:

			accum_color_red += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)value2].red) /
				255.0;
			accum_color_green += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)value2].green) /
				255.0;
			accum_color_blue += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)value2].blue) /
				255.0;

			break;

		  case C_TEXTURE_COLOR:

			C_error_message("You can't have a texture!\n");
			return;
			break;
		}
	    }

	    point_one.x = point_two.x;
	    point_one.y = point_two.y;
	    point_one.z = point_two.z;

	    point_two.x += delta_distance * new_ray_direction.x;
	    point_two.y += delta_distance * new_ray_direction.y;
	    point_two.z += delta_distance * new_ray_direction.z;

	    distance = delta_distance;

	    voxel_x = C_Floor( point_two.x );
	    voxel_y = C_Floor( point_two.y );
	    voxel_z = C_Floor( point_two.z );

	    dptr = volume->data.scalar_data_8bit->scalar_data + 
		   voxel_z * zstep + 
		   voxel_y * ystep + 
	  	   voxel_x * xstep;

	    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 = value2;
	  }
	}
	else
	{
	  while ( exited_volume == FALSE )
	  {
	    value2 = 
		    C_Trilin( 
		    	    ( point_two.x - (float) voxel_x ),
			    ( point_two.y - (float) voxel_y ),
			    ( point_two.z - (float) voxel_z ),
			      A, B, C, D, E, F, G, H );

	    if ( distance > 0.0 )
	    {
		curr_opacity = 
			scan_data->seg.s_opacity->opacity[
				(int)((value1 + value2) / 2.0)] *
				distance;

		remaining_opacity *= curr_opacity;

		if ( remaining_opacity < 0.0 )
			curr_opacity = 1.0;

		switch ( volume->color_type )
		{
		  case C_SINGLE_COLOR:

			accum_color_red += curr_opacity * 
				(float)(
				volume->color_data.single_color->red)
				/ 255.0;
			accum_color_green += curr_opacity * 
				(float)(
				volume->color_data.single_color->green)
				/ 255.0;
			accum_color_blue += curr_opacity * 
				(float)(
				volume->color_data.single_color->blue)
				/ 255.0;
	
			break;

		  case C_COLOR_ARRAY:

			accum_color_red += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)((value1+value2)/2.0)].red) /
				255.0;
			accum_color_green += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)((value1+value2)/2.0)].green) /
				255.0;
			accum_color_blue += curr_opacity * 
				(float)(
				volume->color_data.color_array->color[
				(int)((value1+value2)/2.0)].blue) /
				255.0;

			break;

		  case C_TEXTURE_COLOR:

			C_error_message("You can't have a texture!\n");
			return;
			break;
		}
	    }

	    if ( exited_volume == FALSE && remaining_opacity < 0.0 )
	    {
		remaining_opacity = 0.0;
		exited_volume = TRUE;
	    }
	    else if ( exited_volume == FALSE )
	    {	
	      point_one.x = point_two.x;
	      point_one.y = point_two.y;
	      point_one.z = point_two.z;
  	      value1 = value2;

	      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)
		  {
		    exited_volume = TRUE;
		  }
		  else
		  {
		    delta_t = ((float)(voxel_x + (tstep_x < 0)) - 
				point_one.x) / new_ray_direction.x;
		    point_two.x = point_one.x + 
					delta_t * new_ray_direction.x;
		    point_two.y = point_one.y + 
					delta_t * new_ray_direction.y;
		    point_two.z = point_one.z + 
					delta_t * new_ray_direction.z;
		    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)
		  {
		    exited_volume = TRUE;
		  }
		  else
		  {
		    delta_t = ((float)(voxel_z + (tstep_z < 0)) - 
				point_one.z) / new_ray_direction.z;
		    point_two.x = point_one.x + 
					delta_t * new_ray_direction.x;
		    point_two.y = point_one.y + 
					delta_t * new_ray_direction.y;
		    point_two.z = point_one.z + 
					delta_t * new_ray_direction.z;
		    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)
		  {
		    exited_volume = TRUE;
		  }
		  else
		  {
		    delta_t = ((float)(voxel_y + (tstep_y < 0)) - 
				point_one.y) / new_ray_direction.y;
		    point_two.x = point_one.x + 
					delta_t * new_ray_direction.x;
		    point_two.y = point_one.y + 
					delta_t * new_ray_direction.y;
		    point_two.z = point_one.z + 
					delta_t * new_ray_direction.z;
		    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)
		  {
		    exited_volume = TRUE;
		  }
		  else
		  {
		    delta_t = ((float)(voxel_z + (tstep_z < 0)) - 
				point_one.z) / new_ray_direction.z;
		    point_two.x = point_one.x + 
					delta_t * new_ray_direction.x;
		    point_two.y = point_one.y + 
					delta_t * new_ray_direction.y;
		    point_two.z = point_one.z + 
					delta_t * new_ray_direction.z;
		    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);
		    }
		  }
	        }
	      }
	    }

	    distance = delta_t * scale_distance;
	  }
	}
		    
	shade_info->opacity = (1.0 - remaining_opacity);
	shade_info->opacity_red = accum_color_red;
	shade_info->opacity_green = accum_color_green;
	shade_info->opacity_blue = accum_color_blue;
		
	C_transform_fposition( &point_two, &junk_point, 
				&(volume->ltow_voxels) );

	shade_info->world_point.x = junk_point.x + 
		3.0*rt_info->ray_direction.x;
	shade_info->world_point.y = junk_point.y +
		3.0*rt_info->ray_direction.y;
	shade_info->world_point.z = junk_point.z +
		3.0*rt_info->ray_direction.z;
 
}


