/****************************************************************************/
/*                                                                          */
/*  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_cell2.c
 *	      Author:  Hui Chen
 *		Date:  9/24/1992
 * 	 Description:  This file contains routines: s_seed_cell2().
 *		       This is device_independent.
 *
 *		       s_seed_cell2() is the routine which calculates the
 *		       surface area by "marching cubes" algorithm.
 *
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */


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

/****************************************/
/*					
/*	Procedure Name:  s_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 respectly
/*			 C_MeaVoxelList  *mea_surface_voxels[] - surface list
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine calculates the surface area 
/*			 by "marching cubes" algorithm.
 */


void s_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 area_cal2();       /* calculate the surface area inside the cell */

	extern int cube_table[256][4];

	C_MeaVoxelList	*vptr, *tempptr;
	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;

	float	surf_area;	/* variable contains calculated surface area */
	
	/* set up the increaments */

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

	surf_area = 0.;		/* initialization */

	/* now start calculating the surface area */

	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;

		    /* 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 surface area
		       inside the cell according the way specified for
		       that standard pattern
		    */
		    if ((cell_index != 0) && (cell_index != 255))
		    {
			/* 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 surface area */
			    surf_area += area_cal2(low, high, 
						   cube_table[cell_index][2],
						  rx, ry, rz);
		    }
		    tempptr = vptr;
		    vptr = vptr->next_voxel;
		    free(tempptr);
	        }
	    }
	}
	else {
	    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;

		    /* 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 surface area
		       inside the cell according the way specified for
		       that standard pattern
		    */
		    if ((cell_index != 0) && (cell_index != 255))
		    {
			/* 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 surface area */
			    switch(cube_table[cell_index][1]) {
				case 0: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  rx, ry, rz);
					break;
				case 1: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  ry, rx, rz);
					break;
				case 2: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  rx, rz, ry);
					break;
				case 3: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  rz, ry, rx);
					break;
				case 4: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  ry, rz, rx);
					break;
				case 5: surf_area += 
					 area_cal2(low, high, 
						   cube_table[cell_index][2],
						  rz, rx, ry);
					break;
				default: 
					printf("error in data file: axis_case is invalid.\n");
					break;
			    }
		    }
		    tempptr = vptr;
		    vptr = vptr->next_voxel;
		    free(tempptr);
	        }
	    }
	}
	/* since I didn't do "devided by 2. i.e. /2." in each triangle area
	   calculations in order to save time, I need to do it for the sum 
	*/
	surf_area = surf_area / 2.;
	printf("Surface area = %f in square of unit.\n", surf_area);
}


