/****************************************************************************/
/*                                                                          */
/*  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_cell1.c
 *	      Author:  Hui Chen
 *		Date:  Feb. 7, 1992
 * 	 Description:  This file contains routines: s_cell1(),
 *						    compute_linear_cut().
 *		       They are device_independent.
 *
 *		       s_cell1() is the routine which calculates the
 *		       surface area by "marching cubes" algorithm.
 *
 *		       compute_linear_cut() is the routine which calculates
 *		       the intersection of surface with each edge of cell
 *		       by linear interpolation.  
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */

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

/****************************************/
/*					
/*	Procedure Name:  s_cell1
/*	  Return Value:  float
/*     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
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine calculates the surface area 
/*			 by "marching cubes" algorithm.
 */

float s_cell1(data,xsize,ysize,zsize,low,high,rx,ry,rz)
C_Voxel_8bit  *data;
int	 xsize,ysize,zsize;
int      low, high;
float	 rx, ry, rz;
{
	extern void rotate_cube();	/* rotate a cell to match a standard pattern */
	extern float area_cal1();       /* calculate the surface area inside a cell  */
	void compute_linear_cut();

	extern C_MeaInfo        mea_info;
	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	x, y, z;	/* integer position loop variables */
	int	cell_index;	/* the index value of the cell     */

	float	xcut[256][256]; /* tables which contains the linear */
	float	ycut[256][256]; /* interpolated intersections(cuts) */
	float	zcut[256][256]; /* for edges along x, y, z-axis     */  
				/* direction, respectively	    */

	float	surf_area;	/* variable contains calculated surface area */
	

	/* pre-compute the intersections of surface with edges
	   and save results in global table cut[][] for later usages */

	compute_linear_cut(low,high, rx, ry, rz, xcut, ycut, zcut);

	/* 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 (z=0; z<(zsize-1); z++)
	        for (y=0; y<(ysize-1); y++)
	        {		
		/* initialize dptr for this row */

		    dptr = data + (z*zstep) + (y*ystep);

		    for (x=0; x<(xsize-1); x++)
		    {
		    /* 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_cal1(cube_table[cell_index][2],
						  xcut, xcut, xcut,
						  rx, ry, rz);
		        }
		        dptr += xstep;
		    }
	        }
	}
	else {
	    for (z=0; z<(zsize-1); z++)
	        for (y=0; y<(ysize-1); y++)
	        {		
		/* initialize dptr for this row */

		    dptr = data + (z*zstep) + (y*ystep);

		    for (x=0; x<(xsize-1); x++)
		    {
		    /* 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_cal1(cube_table[cell_index][2],
						  xcut, ycut, zcut,
						  rx, ry, rz);
					break;
				case 1: surf_area += 
					 area_cal1(cube_table[cell_index][2],
						  ycut, xcut, zcut,
						  ry, rx, rz);
					break;
				case 2: surf_area += 
					 area_cal1(cube_table[cell_index][2],
						  xcut, zcut, ycut,
						  rx, rz, ry);
					break;
				case 3: surf_area += 
					 area_cal1(cube_table[cell_index][2],
						  zcut, ycut, xcut,
						  rz, ry, rx);
					break;
				case 4: surf_area += 
					 area_cal1(cube_table[cell_index][2],
						  ycut, zcut, xcut,
						  ry, rz, rx);
					break;
				case 5: surf_area += 
					 area_cal1(cube_table[cell_index][2],
						  zcut, xcut, ycut,
						  rz, rx, ry);
					break;
				default: 
					printf("error in data file: axis_case is invalid.\n");
					break;
			    }
		        }
		        dptr += xstep;
		    }
	        }
	}
	if (mea_info.mea_method == C_CELL_METHOD) {
	/* 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.;
	}
	else {
	/* convert the count of the tiny areas on the surface to a reasonable
	   value of surface area in square units.
	 */
		surf_area = surf_area * (rx*ry + rx*rz + ry*rz) /
		      (3.0 * (float)(mea_info.steps) * (float)(mea_info.steps));
	}
	return(surf_area);
}


/****************************************/
/*					
/*	Procedure Name:  compute_linear_cut
/*	  Return Value:  none
/*     Input Variables:  int	low	- low threshold value used for
/*					  segmentation
/*			 int	high	- high threshold value used for
/*					  segmentation
/*			 float 	rx	- units/voxels ratio in x-axis
/*			 float 	ry	- units/voxels ratio in y-axis
/*			 float 	rz	- units/voxels ratio in z-axis
/*    Output Variables:  float  xcut[256][256]	- table which contains the
/*                                                linear interpolated
/*                                                intersections(cuts) for edges
/*                                                along x-axis direction
/*                       float  ycut[256][256]  - table which contains the
/*                                                linear interpolated
/*                                                intersections(cuts) for edges
/*                                                along y-axis direction
/*                       float  zcut[256][256]  - table which contains the
/*                                                linear interpolated
/*                                                intersections(cuts) for edges
/*                                                along z-axis direction
/*    Update Variables:  none
/*    Global Variables:  none 	
/*     	   Description:  This routine pre-calculates the intersections of 
/*			 surface(defined by two threshold values) with 
/*			 each edge of cells by linear interpolation, 
/*			 and saves the results in the table cut[][] for
/*			 later usages.  
 */

void compute_linear_cut(low,high, rx, ry, rz, xcut, ycut, zcut)
int  low, high;
float rx, ry, rz;
float  xcut[256][256];
float  ycut[256][256];
float  zcut[256][256];
{
	int i, j;	/* loop control variables */
	    
	if (rx == ry) {
	    if (rx == rz) {
		if (rx == 1.0) { /* rx = ry = rz =1. */
		    /* calculate the elements in the table for edges with 
	   	       one vertex has intensity inside the range and the 
	   	       other has intensity below the low threshold value */
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   xcut[i][j] = 1. - xcut[j][i];
	    	    }

		    /* calculate the elements in the table for edges with 
		       one vertex has intensity inside the range and the 
	   	       other has intensity upper the high threshold value */
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   xcut[j][i] = 1. - xcut[i][j];
	    	    }
		}
		else { /* rx = ry = rz */
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = (((float)(low - j))/((float)(i - j)))*rx;
			   xcut[i][j] = rx - xcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = (((float)(high - i))/((float)(j - i)))*rx;
			   xcut[j][i] = rx - xcut[i][j];
	    	    }
		}
	    } /* end of ( rx = ry = rz ) */
	    else {
		if (rx == 1.0) { /* rx = ry = 1. */
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   xcut[i][j] = 1. - xcut[j][i];
			   ycut[j][i] = xcut[j][i];
			   ycut[i][j] = xcut[i][j];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   xcut[j][i] = 1. - xcut[i][j];
			   ycut[i][j] = xcut[i][j];
			   ycut[j][i] = xcut[j][i];
	    	    }
		}
		else { /* rx = ry */
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = (((float)(low - j))/((float)(i - j)))*rx;
			   xcut[i][j] = rx - xcut[j][i];
			   ycut[j][i] = xcut[j][i];
			   ycut[i][j] = xcut[i][j];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = (((float)(high - i))/((float)(j - i)))*rx;
			   xcut[j][i] = rx - xcut[i][j];
			   ycut[i][j] = xcut[i][j];
			   ycut[j][i] = xcut[j][i];
	    	    }
		}
		if (rz == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   zcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   zcut[i][j] = 1. - zcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   zcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   zcut[j][i] = 1. - zcut[i][j];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   zcut[j][i] = (((float)(low - j))/((float)(i - j)))*rz;
			   zcut[i][j] = rz - zcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   zcut[i][j] = (((float)(high - i))/((float)(j - i)))*rz;
			   zcut[j][i] = rz - zcut[i][j];
	    	    }
		}
	    } /* end of ( rx = ry, but rz <> rx ) */
	} /* end of ( rx = ry ) */
	else {
	    if (rx == rz) { /* rx = rz, but rx <> ry */
		if (rx == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   xcut[i][j] = 1. - xcut[j][i];
			   zcut[j][i] = xcut[j][i];
			   zcut[i][j] = xcut[i][j];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   xcut[j][i] = 1. - xcut[i][j];
			   zcut[i][j] = xcut[i][j];
			   zcut[j][i] = xcut[j][i];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = (((float)(low - j))/((float)(i - j)))*rx;
			   xcut[i][j] = rx - xcut[j][i];
			   zcut[j][i] = xcut[j][i];
			   zcut[i][j] = xcut[i][j];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = (((float)(high - i))/((float)(j - i)))*rx;
			   xcut[j][i] = rx - xcut[i][j];
			   zcut[i][j] = xcut[i][j];
			   zcut[j][i] = xcut[j][i];
	    	    }
		}
		if (ry == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   ycut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   ycut[i][j] = 1. - ycut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   ycut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   ycut[j][i] = 1. - ycut[i][j];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   ycut[j][i] = (((float)(low - j))/((float)(i - j)))*ry;
			   ycut[i][j] = ry - ycut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   ycut[i][j] = (((float)(high - i))/((float)(j - i)))*ry;
			   ycut[j][i] = ry - ycut[i][j];
	    	    }
		}
	    } /* end of ( rx = rz, but rx <> ry ) */
	    else {
		if (rx == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   xcut[i][j] = 1. - xcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   xcut[j][i] = 1. - xcut[i][j];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   xcut[j][i] = (((float)(low - j))/((float)(i - j)))*rx;
			   xcut[i][j] = rx - xcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   xcut[i][j] = (((float)(high - i))/((float)(j - i)))*rx;
			   xcut[j][i] = rx - xcut[i][j];
	    	    }
		}
		if (ry == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   ycut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   ycut[i][j] = 1. - ycut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   ycut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   ycut[j][i] = 1. - ycut[i][j];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   ycut[j][i] = (((float)(low - j))/((float)(i - j)))*ry;
			   ycut[i][j] = ry - ycut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   ycut[i][j] = (((float)(high - i))/((float)(j - i)))*ry;
			   ycut[j][i] = ry - ycut[i][j];
	    	    }
		}
		if (rz == 1.0) {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   zcut[j][i] = ((float)(low - j)) / ((float)(i - j));
			   zcut[i][j] = 1. - zcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   zcut[i][j] = ((float)(high - i)) / ((float)(j - i));
			   zcut[j][i] = 1. - zcut[i][j];
	    	    }
		}
		else {
 		    for (i=low; i<=high; i++) 
	    		for (j=0; j<low; j++) {
			   zcut[j][i] = (((float)(low - j))/((float)(i - j)))*rz;
			   zcut[i][j] = rz - zcut[j][i];
	    	    }
		    for (i=low; i<=high; i++) 
	    		for (j=high+1; j<=255; j++) {
			   zcut[i][j] = (((float)(high - i))/((float)(j - i)))*rz;
			   zcut[j][i] = rz - zcut[i][j];
	    	    }
		}
	    } /* end of ( rx <> ry <> rz ) */
	} /* end of ( rx <> ry ) */
}
