/****************************************************************************/
/*                                                                          */
/*  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_mea_seed_v_cell2.c
 *	      Author:  Hui Chen
 *		Date:  Dec. 7, 1991
 * 	 Description:  This file contains routines: v_seed_cell2(),
 *		       It is device_independent.
 *
 *		       v_seed_cell2() is the routine which calculates the
 *		       volume by "marching cubes" algorithm.
 *
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */

#include "C_volvis.h"
#include "C_measure.h"

/****************************************/
/*					
/*	Procedure Name:  v_seed_cell2
/*	  Return Value:  none
/*     Input Variables:  C_Voxel_8bit	*data	- pointer to raw data array
/*			 int		xsize	- x resolution of data array
/*			 int		ysize	- y resolution of data array
/*			 int		zsize	- z resolution of data array
/*			 int		low	- low threshold value used for
/*						  segmentation
/*			 int		high	- high threshold value used for
/*						  segmentation
/*                       float rx, ry, rz       - units/voxels ratios in
/*                                                x, y, z axis, respectively
/*			 C_MeaVoxelList  *mea_surface_voxels[] - surface list
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine calculates the volume 
/*			 by "marching cubes" algorithm.
 */


void v_seed_cell2(mea_surface_voxels,data,xsize,ysize,zsize,low,high,rx,ry,rz)
C_MeaVoxelList  *mea_surface_voxels[C_MEA_HASH_MAX];
C_Voxel_8bit  *data;
int	 xsize,ysize,zsize;
int      low, high;
float	 rx, ry, rz;
{
	extern void rotate_cube();	/* rotate the cell(cube) */
	extern float volume_cal1();     /* calculate the volume inside the cell */

	extern int cube_table[256][4];

	C_MeaVoxelList  *vptr, *vptr0, *vptr1;
	C_Voxel_8bit *dptr;		/* pointer to raw data		*/
	int	xstep;		/* increment for a step in x	*/
	int	ystep;		/* increment for a step in y	*/
	int	zstep;		/* increment for a step in z	*/
	int	cell_index;	/* the index value of the cell     */
	int	i;		/* loop control variable */

        int     index;
        int     hash_index;
        int     flag;

	float   vol_one_cell; 	/* volume of a single full ocupied cell */

	float	volume;		/* variable contains calculated volume */
	

	vol_one_cell = rx * ry * rz * 6.0; /* since we will do /6.0 to the 
					      sum result at the end */

	/* set up the increaments */

	zstep = xsize*ysize;
	ystep = xsize;
	xstep = 1;

	volume = 0.;		/* initialization */

	/* now start calculating the volume */

	if ((rx == ry) && (rx == rz)) {
	    for (i=0; i<C_MEA_HASH_MAX; i++) {
                vptr = mea_surface_voxels[i];
                while (vptr != NULL)
	        {		
		/* initialize dptr for this voxel */

		    dptr = data + vptr->voxel.index;
                    vptr0 = vptr;
                    vptr = vptr->next_voxel;
                    free(vptr0);

		    flag = 0;
                    if (C_In_Range( *(dptr), low, high)) {
                        if ((C_LEFT_Face(dptr, low, high, xstep)) ||
                            (C_RIGHT_Face(dptr, low, high, xstep)))
                            flag = 1;
                    }
                    else {
                        flag = 1;
                    }

                    if (flag == 1) {
		    /* calculate the index value of the cell depending on
		       the intensities of the eight vertices
		    */
		        index_cal((dptr),low,high,xstep,ystep,zstep,cell_index);

		    /* in case of surface cut through the cell,
		       rotate the cell to match one of the standard 
		       patterns and then calculate the volume 
		       inside the cell according the way specified for
		       that standard pattern
		    */
		        if (cell_index == 255) 
		            volume += vol_one_cell;
		        else if (cell_index != 0)
		        {
			/* rotate the cell */
		            rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			/* calculate the volume */
			    if (cube_table[cell_index][3] == 0)
			       volume += volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz);
			    else
			       volume += (vol_one_cell - 
					 	volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz));
		        }
		    }
		    if (C_In_Range( *(dptr), low, high)) {
			if (C_LEFT_Face(dptr, low, high, xstep)) {
			    dptr ++;
			    while (C_In_Range( *(dptr), low, high)) {
		        	index_cal((dptr),low,high,xstep,ystep,zstep,
						cell_index);

		        	if (cell_index == 255) 
			    	    volume += vol_one_cell;
				else if (cell_index != 0)
		        	{
			    	    rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			    	    if (cube_table[cell_index][3] == 0)
			       	        volume += volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz);
			            else
			                volume += (vol_one_cell - 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz));
		        	}

				if ( C_On_Surface(dptr,low,high,
						  xstep,ystep,zstep)) {
				    index = dptr - data;
                          	    hash_index = index % C_MEA_HASH_MAX;
                          	    vptr0 = mea_surface_voxels[hash_index];
                          	    vptr1 = mea_surface_voxels[hash_index];
                          	    flag = 0;
                          	    while ((vptr1 != NULL) && (flag == 0)){
                               	        if (vptr1->voxel.index == index) {
                                            flag = 1;
                                            if (vptr0 == vptr1) {
                                        	mea_surface_voxels[hash_index] =
                                        	mea_surface_voxels[hash_index]->next_voxel;
					    }
					    else {
                                        	vptr0->next_voxel = 
							vptr1->next_voxel;
                               		    }
                                    	    free(vptr1);
                                	}
					else {
					     vptr0 = vptr1;
                                   	     vptr1 = vptr1->next_voxel;
                                        }

				    }
				}
				dptr ++;
			    }
			}
			else if (C_RIGHT_Face(dptr, low, high, xstep)) {
			    dptr --;
			    while (C_In_Range( *(dptr), low, high)) {
		        	index_cal((dptr),low,high,xstep,ystep,zstep,
						cell_index);

		        	if (cell_index == 255) 
			    	    volume += vol_one_cell;
				else if (cell_index != 0)
		        	{
			    	    rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			    	    if (cube_table[cell_index][3] == 0)
			       	        volume += volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz);
			            else
			                volume += (vol_one_cell - 
					  	volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz));
		        	}

				if ( C_On_Surface(dptr,low,high,
						  xstep,ystep,zstep)) {
				    index = dptr - data;
                          	    hash_index = index % C_MEA_HASH_MAX;
                          	    vptr0 = mea_surface_voxels[hash_index];
                          	    vptr1 = mea_surface_voxels[hash_index];
                          	    flag = 0;
                          	    while ((vptr1 != NULL) && (flag == 0)){
                               	        if (vptr1->voxel.index == index) {
                                            flag = 1;
                                            if (vptr0 == vptr1) {
                                        	mea_surface_voxels[hash_index] =
                                        	mea_surface_voxels[hash_index]->next_voxel;
					    }
					    else {
                                        	vptr0->next_voxel = 
							vptr1->next_voxel;
                               		    }
                                    	    free(vptr1);
                                	}
					else {
					     vptr0 = vptr1;
                                   	     vptr1 = vptr1->next_voxel;
                                        }

				    }
				}
				dptr --;
			    }
			}
		    }
	        }
	    }
	}
	else {
	    for (i=0; i<C_MEA_HASH_MAX; i++) {
                vptr = mea_surface_voxels[i];
                while (vptr != NULL)
	        {		
		/* initialize dptr for this row */

		    dptr = data + vptr->voxel.index;
                    vptr0 = vptr;
                    vptr = vptr->next_voxel;
                    free(vptr0);

		    flag = 0;
                    if (C_In_Range( *(dptr), low, high)) {
                        if ((C_LEFT_Face(dptr, low, high, xstep)) ||
                            (C_RIGHT_Face(dptr, low, high, xstep)))
                            flag = 1;
                    }
                    else {
                        flag = 1;
                    }

                    if (flag == 1) {
		    /* calculate the index value of the cell depending on
		       the intensities of the eight vertices
		    */
		        index_cal((dptr),low,high,xstep,ystep,zstep,cell_index);

		    /* in case of surface cut through the cell,
		       rotate the cell to match one of the standard 
		       patterns and then calculate the volume 
		       inside the cell according the way specified for
		       that standard pattern
		    */
		        if (cell_index == 255) 
		            volume += vol_one_cell;
		        else if (cell_index != 0)
		        {
			/* rotate the cell */
		            rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			/* calculate the volume */
			    if (cube_table[cell_index][3] == 0) {
			        switch(cube_table[cell_index][1]) {
				case 0: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  rx, ry, rz);
					break;
				case 1: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  ry, rx, rz);
					break;
				case 2: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  rx, rz, ry);
					break;
				case 3: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  rz, ry, rx);
					break;
				case 4: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  ry, rz, rx);
					break;
				case 5: volume += volume_cal2(low, high, 
						  cube_table[cell_index][2],
						  rz, rx, ry);
					break;
				default: 
					printf("error in data file: axis_case is invalid.\n");
					break;
			        }
			    }
			    else {
			        switch(cube_table[cell_index][1]) {
				case 0: volume += (vol_one_cell - 
						   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   rx, ry, rz));
					break;
				case 1: volume += (vol_one_cell -
					 	   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   ry, rx, rz));
					break;
				case 2: volume += (vol_one_cell -
					 	   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   rx, rz, ry));
					break;
				case 3: volume += (vol_one_cell -
					 	   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   rz, ry, rx));
					break;
				case 4: volume += (vol_one_cell -
					  	   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   ry, rz, rx));
					break;
				case 5: volume += (vol_one_cell -
					 	   volume_cal2(low, high, 
						   cube_table[cell_index][2],
						   rz, rx, ry));
					break;
				default: 
					printf("error in data file: axis_case is invalid.\n");
					break;
			        }
			    }
		        }
		    }
		    if (C_In_Range( *(dptr), low, high)) {
			if (C_LEFT_Face(dptr, low, high, xstep)) {
			    dptr ++;
			    while (C_In_Range( *(dptr), low, high)) {
		        	index_cal((dptr),low,high,xstep,ystep,zstep,
						cell_index);

		        	if (cell_index == 255) 
			    	    volume += vol_one_cell;
				else if (cell_index != 0)
		        	{
			    	    rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			    	    if (cube_table[cell_index][3] == 0) {
			    		switch(cube_table[cell_index][1]) {
					case 0: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz);
						break;
					case 1: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rx, rz);
						break;
					case 2: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, rz, ry);
						break;
					case 3: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, ry, rx);
						break;
					case 4: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rz, rx);
						break;
					case 5: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, rx, ry);
						break;
					default: 
						printf("error in data file: axis_case is invalid.\n");
						break;
			    		}
				    }
				    else {
			    		switch(cube_table[cell_index][1]) {
					case 0: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz));
						break;
					case 1: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rx, rz));
						break;
					case 2: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, rz, ry));
						break;
					case 3: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, ry, rx));
						break;
					case 4: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rz, rx));
						break;
					case 5: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, rx, ry));
						break;
					default: 
						printf("error in data file: axis_case is invalid.\n");
						break;
			    		}
				    }
		        	}

				if ( C_On_Surface(dptr,low,high,
						  xstep,ystep,zstep)) {
				    index = dptr - data;
                          	    hash_index = index % C_MEA_HASH_MAX;
                          	    vptr0 = mea_surface_voxels[hash_index];
                          	    vptr1 = mea_surface_voxels[hash_index];
                          	    flag = 0;
                          	    while ((vptr1 != NULL) && (flag == 0)){
                               	        if (vptr1->voxel.index == index) {
                                            flag = 1;
                                            if (vptr0 == vptr1) {
                                        	mea_surface_voxels[hash_index] =
                                        	mea_surface_voxels[hash_index]->next_voxel;
					    }
					    else {
                                        	vptr0->next_voxel = 
							vptr1->next_voxel;
                               		    }
                                    	    free(vptr1);
                                	}
					else {
					     vptr0 = vptr1;
                                   	     vptr1 = vptr1->next_voxel;
                                        }

				    }
				}
				dptr ++;
			    }
			}
			else if (C_RIGHT_Face(dptr, low, high, xstep)) {
			    dptr --;
			    while (C_In_Range( *(dptr), low, high)) {
		        	index_cal((dptr),low,high,xstep,ystep,zstep,
						cell_index);

		        	if (cell_index == 255) 
			    	    volume += vol_one_cell;
				else if (cell_index != 0)
		        	{
			    	    rotate_cube(cube_table[cell_index][0],
				    (dptr),
				    (dptr+xstep),
				    (dptr+ystep+xstep), 
				    (dptr+ystep),
				    (dptr+zstep),
				    (dptr+zstep+xstep),
				    (dptr+zstep+ystep+xstep),
				    (dptr+zstep+ystep));

			    	    if (cube_table[cell_index][3] == 0) {
			    		switch(cube_table[cell_index][1]) {
					case 0: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz);
						break;
					case 1: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rx, rz);
						break;
					case 2: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, rz, ry);
						break;
					case 3: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, ry, rx);
						break;
					case 4: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rz, rx);
						break;
					case 5: volume += 
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, rx, ry);
						break;
					default: 
						printf("error in data file: axis_case is invalid.\n");
						break;
			    		}
				    }
				    else {
			    		switch(cube_table[cell_index][1]) {
					case 0: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, ry, rz));
						break;
					case 1: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rx, rz));
						break;
					case 2: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rx, rz, ry));
						break;
					case 3: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, ry, rx));
						break;
					case 4: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  ry, rz, rx));
						break;
					case 5: volume += (vol_one_cell -
						volume_cal2(low, high, 
						cube_table[cell_index][2],
						  rz, rx, ry));
						break;
					default: 
						printf("error in data file: axis_case is invalid.\n");
						break;
			    		}
				    }
		        	}

				if ( C_On_Surface(dptr,low,high,
						  xstep,ystep,zstep)) {
				    index = dptr - data;
                          	    hash_index = index % C_MEA_HASH_MAX;
                          	    vptr0 = mea_surface_voxels[hash_index];
                          	    vptr1 = mea_surface_voxels[hash_index];
                          	    flag = 0;
                          	    while ((vptr1 != NULL) && (flag == 0)){
                               	        if (vptr1->voxel.index == index) {
                                            flag = 1;
                                            if (vptr0 == vptr1) {
                                        	mea_surface_voxels[hash_index] =
                                        	mea_surface_voxels[hash_index]->next_voxel;
					    }
					    else {
                                        	vptr0->next_voxel = 
							vptr1->next_voxel;
                               		    }
                                    	    free(vptr1);
                                	}
					else {
					     vptr0 = vptr1;
                                   	     vptr1 = vptr1->next_voxel;
                                        }

				    }
				}
				dptr --;
			    }
			}
		    }
                }
	    }
	}
	/* since I didn't do "devided by 6. i.e. /6." in each tetrahedron 
	   volume calculations in order to save time, I need to do it for 
	   the sum here.
        */
	volume = volume / 6.0;
	printf("Volume = %f in cube of unit.\n", volume);
}


