/****************************************************************************/
/*                                                                          */
/*  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_function.c
 *		      Author: Lisa Sobierajski
 *			Date: 02/01/92
 *		 Description: Function Routines Of VolVis
 *	Modification History:
 *
 *		Who?		When?		Why?
 *	--------------------------------------------------------------------
 *
 */

#include <stdio.h>
#include <stdlib.h>	/* For malloc() and strtod() */
#include <fcntl.h>	/* For open()   */
#include <math.h>

#include "C_volvis.h"
#include "C_function.h"


int C_input_function( filename, volume )
char		*filename;
C_Volume	*volume;
{
	extern	C_World		world;

	C_FunNode	*C_create_function_tree();
	float		C_eval_function();

	FILE		*fp;
	char		string[200];
	char		*ptr;
	char		*new_ptr;
	float		fun_x_min, fun_x_max;
	float		fun_y_min, fun_y_max;
	float		fun_z_min, fun_z_max;
	float		map_to_0, map_to_255, segment_at;
	int		invert_function;
	register int	loop, xloop, yloop, zloop, multloop;
	register C_Byte	*dptr;
	register float	fval, xval, yval, zval;
	register float	x, y, z, tmp;
	int		ival;

	C_FunNode	*function_tree;
	
/*
	sprintf( string, "%s.fun", filename );
*/
	sprintf( string, "%s", filename);
	fp = fopen( string, "r" );

	if ( !fp )
	{
		C_error_message("Could not open function file!!!\n");
		fclose( fp );
		return C_ERROR;
	}

	fscanf( fp, "%s", string );
	printf("%s\n",string);
	fscanf( fp, "%s", string );
	printf("%s\n",string);
	volume->x_size_voxels = (int) strtol( string, &ptr, 10 ) + 2;
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n",string);
	fscanf( fp, "%s", string );
	printf("%s\n",string);
	volume->y_size_voxels = (int) strtol( string, &ptr, 10 ) + 2;
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n",string);
	fscanf( fp, "%s", string );
	printf("%s\n",string);
	volume->z_size_voxels = (int) strtol( string, &ptr, 10 ) + 2;
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	fscanf( fp, "%s", string );
	if ( ( strcmp( string, "METER" )) == 0 )
	{
		printf("METER\n");
		volume->unit_type = C_METER_UNIT;
	}
	else if ( ( strcmp( string, "MILLIMETER" )) == 0 )
	{
		printf("MILLIMETER\n");
		volume->unit_type = C_MILLIMETER_UNIT;
	}
	else if ( ( strcmp( string, "MICRON" )) == 0 )
	{
		printf("MICRON\n");
		volume->unit_type = C_MICRON_UNIT;
	}
	else if ( ( strcmp( string, "INCH" )) == 0 )
	{
		printf("INCH\n");
		volume->unit_type = C_INCH_UNIT;
	}
	else if ( ( strcmp( string, "FOOT" )) == 0 )
	{
		printf("FOOT\n");
		volume->unit_type = C_FOOT_UNIT;
	}
	else
	{
		C_error_message("Unknow Unit Type In Function File\n");
		fclose( fp );
		return C_ERROR;
	}
	
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	volume->x_size_units = (float) strtod( string, &ptr ) *
					(volume->x_size_voxels-1);
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	volume->y_size_units = (float) strtod( string, &ptr ) *
					(volume->y_size_voxels-1);
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	volume->z_size_units = (float) strtod( string, &ptr ) *
					(volume->z_size_voxels-1);
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_x_min = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_x_max = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_y_min = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_y_max = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_z_min = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
		
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fun_z_max = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 

	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	map_to_0 = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 

	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	map_to_255 = (float) strtod( string, &ptr );
	if ( ptr == string )
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 

	fscanf( fp, "%s", string );
	printf("%s\n", string );
	fscanf( fp, "%s", string );
	printf("%s\n", string );
	segment_at = (float) strtod( string, &ptr );
	if ( ptr == string )  
	{
		C_error_message("Invalid Function File Format!\n");
		fclose( fp );
		return C_ERROR;
	} 
	if ( segment_at > map_to_255 || segment_at < map_to_0 )
	{
		C_error_message(
		   "You must segment between the values mapped to 0 and 255" );
		fclose( fp );
		return C_ERROR;
	}

	fscanf( fp, "%s", string );
	fscanf( fp, "%s", string );
	if ( ( strcmp( string, "TRUE" )) == 0 )
	{
		invert_function = TRUE;
	}
	else if ( ( strcmp( string, "FALSE" )) == 0 )
	{
		invert_function = FALSE;
	}
	else
	{
		C_error_message("Invalid File Format\n");
		fclose( fp );
		return C_ERROR;
	}

	volume->icon_8bit = NULL;
	volume->icon_24bit = NULL;

	volume->data_type = C_SCALAR_DATA_8BIT;
	volume->data_modification = C_ORIGINAL_DATA;
	volume->data.scalar_data_8bit = C_New( C_ScalarData_8bit );
	
	volume->data.scalar_data_8bit->scalar_data = dptr = (C_Byte *)
		malloc( volume->x_size_voxels *
			volume->y_size_voxels *
			volume->z_size_voxels * sizeof(C_Byte) );

	function_tree = C_create_function_tree( fp );

	for (loop = 0; 
	     loop < volume->x_size_voxels * volume->y_size_voxels; loop++)
		*(dptr++) = 0;

	for ( zloop = 1; zloop < volume->z_size_voxels - 1; zloop++ )
	{
printf("zloop = %d\n",zloop);

	  for (loop = 0; loop < volume->x_size_voxels; loop++)
		*(dptr++) = 0;

	  z = fun_z_min + ((float)(zloop-1)/(float)(volume->z_size_voxels-2))*
			  (fun_z_max - fun_z_min);

	  for ( yloop = 1; yloop < volume->y_size_voxels - 1; yloop++ )
	  {
	    *(dptr++) = 0;
	    y = fun_y_min + ((float)(yloop-1)/(float)(volume->y_size_voxels-2))*
			    (fun_y_max - fun_y_min);

	    for ( xloop = 1; xloop < volume->x_size_voxels - 1; xloop++ )
	    {
		x = fun_x_min + ((float)(xloop-1)/
			        (float)(volume->x_size_voxels-2))*
				(fun_x_max - fun_x_min);


		fval = C_eval_function( function_tree, x, y, z );

		if ( invert_function )
			ival = (map_to_255 - fval) * 255.0 / 
				(map_to_255 - map_to_0);
		else
			ival = (fval - map_to_0) * 255.0 / 
				(map_to_255 - map_to_0);


		if ( ival < 0 )   ival = 0;
		if ( ival > 255 ) ival = 255;

		*(dptr++) = ival;
	    }
	    *(dptr++) = 0;
	  }
	  for (loop = 0; loop < volume->x_size_voxels; loop++)
		*(dptr++) = 0;
	}

	for (loop = 0; 
	     loop < volume->x_size_voxels * volume->y_size_voxels; loop++)
		*(dptr++) = 0;

	/* Initialize Any Additional Volume Entries */
        C_create_volume_matrices( volume );
        C_create_plane_equations( volume );
	
	volume->data_origin = C_FUN_VOXELIZED_DATA;

        volume->data.scalar_data_8bit->seg_type     = C_ISO_SURFACE;
        volume->data.scalar_data_8bit->seg.isovalue  = 
		( segment_at - map_to_0) * 255.0 / ( map_to_255 - map_to_0 );
	volume->data.scalar_data_8bit->seg.s_opacity = C_New( C_OpacityArray );

	for ( loop = 0; loop < 255; loop++ )
	    volume->data.scalar_data_8bit->seg.s_opacity->opacity[loop] = 
			loop / 255.0;
	
	fclose( fp );
	return C_OK;
}

float C_eval_function( root, x, y, z )
C_FunNode		*root;
float			x, y, z;
{
	float		return_val;
	float		left_val, right_val;
	int		loop;

	if ( !root )
	{
		printf("Oh No - NULL root\n");
		return 0.0;
	}

	switch( root->node_type )
	{
	  case C_CONST:
		return root->node_value;
		break;

	  case C_VAR_X:
		return x;
		break;

	  case C_VAR_Y:
		return y;
		break;

	  case C_VAR_Z:
		return z;
		break;

	  case C_POWER:
		left_val = C_eval_function( root->left_child, x, y, z );	
		right_val = C_eval_function( root->right_child, x, y, z );	
		
		return_val = 1.0;

		for (loop = 0; loop < fabs((double)right_val); loop++ )
			return_val *= left_val;

		if ( right_val < 0.0 )
			return_val = 1.0 / return_val;
		return return_val;
		break;

	  case C_ADD:
		left_val = C_eval_function( root->left_child, x, y, z );	
		right_val = C_eval_function( root->right_child, x, y, z );	
		return ( left_val + right_val );
		break;

	  case C_SUB:
		left_val = C_eval_function( root->left_child, x, y, z );	
		right_val = C_eval_function( root->right_child, x, y, z );	
		return ( left_val - right_val );
		break;

	  case C_MUL:
		left_val = C_eval_function( root->left_child, x, y, z );	
		right_val = C_eval_function( root->right_child, x, y, z );	
		return ( left_val * right_val );
		break;

	  case C_DIV:
		left_val = C_eval_function( root->left_child, x, y, z );	
		right_val = C_eval_function( root->right_child, x, y, z );	
		if ( right_val != 0.0 )
			return ( left_val / right_val );
		else
			return C_BIG_FLOAT;
		break;

	  case C_SIN:
		left_val = C_eval_function( root->left_child, x, y, z );	
		return	(float) sin( (double) left_val );
		break;

	  case C_COS:
		left_val = C_eval_function( root->left_child, x, y, z );	
		return	(float) cos( (double) left_val );
		break;

	  case C_TAN:
		left_val = C_eval_function( root->left_child, x, y, z );	
		return	(float) tan( (double) left_val );
		break;

	  default:
		printf("Unknown node type in function tree\n");
		return 0.0;
	}
}

C_FunNode *C_create_function_tree( fp )
FILE		*fp;
{
	char		string[20];
	char		*ptr;
	C_FunNode	*node;
	
	fscanf( fp, "%s", string );

	node = C_New( C_FunNode );
	
	if (strcmp( string, "CONST" ) == 0 )
	{
printf("constant\n");
		node->node_type = C_CONST;
		fscanf( fp, "%s", string );
		node->node_value = (float) strtod( string, &ptr ); 
		if ( ptr == string )
		{
			C_error_message("Invalid Function File\n");
			return NULL;
		}
printf("value = %f\n", node->node_value);
		node->left_child = NULL;
		node->right_child = NULL;
	}
	else if (strcmp( string, "X" ) == 0 )
	{
printf("X\n");
		node->node_type = C_VAR_X;
		node->left_child = NULL;
		node->right_child = NULL;
	}	
	else if (strcmp( string, "Y" ) == 0 )
	{
printf("Y\n");
		node->node_type = C_VAR_Y;
		node->left_child = NULL;
		node->right_child = NULL;
	}	
	else if (strcmp( string, "Z" ) == 0 )
	{
printf("Z\n");
		node->node_type = C_VAR_Z;
		node->left_child = NULL;
		node->right_child = NULL;
	}	
	else if (strcmp( string, "POWER" ) == 0 )
	{
printf("Power\n");
		node->node_type = C_POWER;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = C_create_function_tree( fp );
		if ( !node->right_child )
			return NULL;
	}	
	else if (strcmp( string, "ADD" ) == 0 )
	{
printf("Add\n");
		node->node_type = C_ADD;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = C_create_function_tree( fp );
		if ( !node->right_child )
			return NULL;
	}	
	else if (strcmp( string, "SUB" ) == 0 )
	{
printf("Sub\n");
		node->node_type = C_SUB;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = C_create_function_tree( fp );
		if ( !node->right_child )
			return NULL;
	}	
	else if (strcmp( string, "MUL" ) == 0 )
	{
printf("Mul\n");
		node->node_type = C_MUL;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = C_create_function_tree( fp );
		if ( !node->right_child )
			return NULL;
	}	
	else if (strcmp( string, "DIV" ) == 0 )
	{
printf("Div\n");
		node->node_type = C_DIV;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = C_create_function_tree( fp );
		if ( !node->right_child )
			return NULL;
	}	
	else if (strcmp( string, "SIN" ) == 0 )
	{
printf("Sin\n");
		node->node_type = C_SIN;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = NULL;
	}	
	else if (strcmp( string, "COS" ) == 0 )
	{
printf("Cos\n");
		node->node_type = C_COS;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = NULL;
	}	
	else if (strcmp( string, "TAN" ) == 0 )
	{
printf("Tan\n");
		node->node_type = C_TAN;
		node->left_child = C_create_function_tree( fp );
		if ( !node->left_child )
			return NULL;
		node->right_child = NULL;
	}	

	return node;
}
