/****************************************************************************/
/*                                                                          */
/*  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_normal( volume, shade_info ) 
C_Volume		*volume;
C_ShadeInfo		*shade_info;

{
	extern 			C_transform_fposition(); 
	extern 			C_transform_fvector(); 

	extern C_View		view;

	int			exiting_volume;
	C_Voxel_8bit		*dptr;
	float			xratio, yratio, zratio;
	C_FPosition		new_point;
	C_FPosition		point_one;
	C_FPosition		junk_point;
	C_FVector		normal;
	int			xstep, ystep, zstep;

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

	C_GeoData		*table_entry;
	int			voxel_x, voxel_y, voxel_z;


	switch ( volume->data_type )
	{
	    case C_SCALAR_DATA_8BIT:

	        scan_data = volume->data.scalar_data_8bit;

		exiting_volume = shade_info->exiting_volume;

		if ( shade_info->world_normal.x == C_BIG_FLOAT )
		{

		  xstep = 1;
		  ystep = volume->x_size_voxels;
		  zstep = volume->x_size_voxels * volume->y_size_voxels;

		  if ( view.image_accuracy == C_LOW_IMAGE_ACCURACY )
		  {
		      voxel_x = C_Round( shade_info->local_point.x );
		      voxel_y = C_Round( shade_info->local_point.y );
		      voxel_z = C_Round( shade_info->local_point.z );

		      if ( voxel_x < 1 )
			voxel_x = 1;
		      else if ( voxel_x > volume->x_size_voxels - 1 )
			voxel_x = volume->x_size_voxels - 1;

		      if ( voxel_y < 1 )
			voxel_y = 1;
		      else if ( voxel_y > volume->y_size_voxels - 1 )
			voxel_y = volume->y_size_voxels - 1;

		      if ( voxel_z < 1 )
			voxel_z = 1;
		      else if ( voxel_z > volume->z_size_voxels - 1 )
			voxel_z = volume->z_size_voxels - 1;

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

		      normal.x = (int)*(dptr - xstep) - (int)*(dptr + xstep);
		      normal.y = (int)*(dptr - ystep) - (int)*(dptr + ystep);
		      normal.z = (int)*(dptr - zstep) - (int)*(dptr + zstep);

		      if (normal.x == 0.0 && normal.y == 0.0 && normal.z == 0.0)
		      {
printf("Ack!!! - zero normal!!!\n");
		        normal.x = shade_info->local_point.x -
		                (volume->x_size_voxels/2.0);
		        normal.y = shade_info->local_point.y -
		                (volume->y_size_voxels/2.0);
		        normal.z = shade_info->local_point.z -
		                (volume->z_size_voxels/2.0);
		      }




		      junk_point.x = shade_info->local_point.x * 
				   volume->x_size_units /
				   (float)(volume->x_size_voxels-1);

		      junk_point.y = shade_info->local_point.y * 
				   volume->y_size_units /
				   (float)(volume->y_size_voxels-1);

		      junk_point.z = shade_info->local_point.z * 
				   volume->z_size_units /
				   (float)(volume->z_size_voxels-1);


                      C_transform_fvector( &normal, &junk_point,
			&(shade_info->world_normal), &new_point, 
		 	&(volume->ltow_units) );


		  }
		  else
		  {

		      normal_step_size = (1.0*volume->x_size_units) / 
				(3.0*volume->x_size_voxels); 

		      xratio = ((float) (volume->x_size_voxels-1) /
			     volume->x_size_units);

		      yratio = ((float) (volume->y_size_voxels-1) /
			     volume->y_size_units);

		      zratio = ((float) (volume->z_size_voxels-1) /
			     volume->z_size_units);

		      xnormal_step_size = 
			    normal_step_size * xratio;

		      ynormal_step_size = 
			    normal_step_size * yratio; 

		      znormal_step_size = 
			    normal_step_size * zratio; 


                      if ( shade_info->local_point.x < xnormal_step_size )
		      {
                        normal_step_size = shade_info->local_point.x / xratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }

                      if ( shade_info->local_point.y < ynormal_step_size )
		      {
                        normal_step_size = shade_info->local_point.y / yratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }

                      if (shade_info->local_point.z < znormal_step_size )
		      {
                        normal_step_size = shade_info->local_point.z / zratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }

                      if ( ((float)(volume->x_size_voxels-1) - 
				shade_info->local_point.x) <
		          xnormal_step_size )
		      {
                        normal_step_size = ((float)(volume->x_size_voxels - 
			    1.0001) - shade_info->local_point.x) / xratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }

                      if ( ((float)(volume->y_size_voxels-1) - 
				shade_info->local_point.y) <
		          ynormal_step_size )
		      {
                        normal_step_size = ((float)(volume->y_size_voxels - 
			    1.0001) - shade_info->local_point.y) / yratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }

                      if ( ((float)(volume->z_size_voxels-1) - 
				shade_info->local_point.z) <
		          znormal_step_size )
		      {
                        normal_step_size = ((float)(volume->z_size_voxels - 
			    1.0001) - shade_info->local_point.z) / zratio;

			xnormal_step_size = 
			    normal_step_size * xratio;

			ynormal_step_size = 
			    normal_step_size * yratio;

			znormal_step_size = 
			    normal_step_size * zratio;
		      }


                      point_one.x = shade_info->local_point.x + 
					xnormal_step_size;
                      point_one.y = shade_info->local_point.y;
                      point_one.z = shade_info->local_point.z;

                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value1 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      point_one.x = shade_info->local_point.x - 
					xnormal_step_size;
                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value2 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      normal.x = value2 - value1;


                      point_one.x = shade_info->local_point.x;
                      point_one.y = shade_info->local_point.y + 
					ynormal_step_size;

                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value1 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      point_one.y = shade_info->local_point.y - 
					ynormal_step_size;

                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value2 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      normal.y = value2 - value1;


                      point_one.y = shade_info->local_point.y;
                      point_one.z = shade_info->local_point.z + 
					znormal_step_size;

                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value1 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      point_one.z = shade_info->local_point.z - 
					znormal_step_size;

                      dptr = (scan_data->scalar_data) +
                             ((int) point_one.x) * xstep +
                             ((int) point_one.y) * ystep +
                             ((int) point_one.z) * zstep;

                      value2 = C_Trilin( point_one.x - (int)(point_one.x),
                                       point_one.y - (int)(point_one.y),
                                       point_one.z - (int)(point_one.z),
                                       *(dptr),
                                       *(dptr + xstep),
                                       *(dptr + ystep),
                                       *(dptr + xstep + ystep),
                                       *(dptr + zstep),
                                       *(dptr + xstep + zstep),
                                       *(dptr + ystep + zstep),
                                       *(dptr + xstep + ystep + zstep) );

                      normal.z = value2 - value1;

		      if (normal.x == 0.0 && normal.y == 0.0 && normal.z == 0.0)
		      {
printf("Ack!!! - zero normal!!!\n");
		        normal.x = shade_info->local_point.x -
		                (volume->x_size_voxels/2.0);
		        normal.y = shade_info->local_point.y -
		                (volume->y_size_voxels/2.0);
		        normal.z = shade_info->local_point.z -
		                (volume->z_size_voxels/2.0);
		      }




		      junk_point.x = shade_info->local_point.x * 
				   volume->x_size_units /
				   (float)(volume->x_size_voxels-1);

		      junk_point.y = shade_info->local_point.y * 
				   volume->y_size_units /
				   (float)(volume->y_size_voxels-1);

		      junk_point.z = shade_info->local_point.z * 
				   volume->z_size_units /
				   (float)(volume->z_size_voxels-1);


                      C_transform_fvector( &normal, &junk_point,
			&(shade_info->world_normal), &new_point, 
		 	&(volume->ltow_units) );

		        
		    }
		}

	        if ( shade_info->exiting_volume )
	        {
	    	    shade_info->world_normal.x = 
		    	-shade_info->world_normal.x;
		    shade_info->world_normal.y = 
			-shade_info->world_normal.y;
		    shade_info->world_normal.z = 
			-shade_info->world_normal.z;
	        }

	    break;

	    case C_GEOMETRIC_DATA:


		    table_entry = volume->data.geo_data;

		    switch ( table_entry->geo_type )
		    {
		      case C_SPHERE:
		      case C_HOLLOW_SPHERE:
		      case C_CYLINDER:
		      case C_HOLLOW_CYLINDER:
		      case C_POLYGON:
		      case C_MULTIPOLY:


			if ( (shade_info->world_normal.x * 
					shade_info->ray_direction.x + 
			      shade_info->world_normal.y * 
					shade_info->ray_direction.y + 
			      shade_info->world_normal.z * 
					shade_info->ray_direction.z) >0)
			{
				shade_info->world_normal.x = 
					-shade_info->world_normal.x;
				shade_info->world_normal.y = 
					-shade_info->world_normal.y;
				shade_info->world_normal.z = 
					-shade_info->world_normal.z;
				exiting_volume = TRUE;
			}
			else
				exiting_volume = FALSE;

			break;

		      default:

			C_error_message("Unkown Geo Object Type!!!\n");
			
			break;
		    }

		    break;

	        default:

		    C_error_message("Error - no such data type\n");

		    break;

	}

}


