/****************************************************************************/
/*                                                                          */
/*  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_volume_cal1.c
 *	      Author:  Hui Chen
 *		Date:  Dec. 9, 1991
 * 	 Description:  This file contains volume_cal1(),
 *		       which calculates the volume 
 *		       inside a cell.
 * 		       It is device-independent.
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */


#include <math.h>
#include "C_volvis.h"
#include "C_measure.h"

/****************************************/
/*					
/*	Procedure Name:  volume_cal1
/*	  Return Value:  float
/*     Input Variables:  int	pattern_num	- pattern index
/*                       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
/*                       float rx, ry, rz       - units/voxels ratios in
/*                                                x, y, z axis, respectively
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
*/

float volume_cal1(pattern_num, xcut, ycut, zcut, rx, ry, rz)
int pattern_num;
float xcut[256][256];
float ycut[256][256];
float zcut[256][256];
float rx, ry, rz;
{
	extern float volume_cal1_anal();
	float volume_cal1_cell();

	extern C_MeaInfo        mea_info;
        float   volume;

        if (mea_info.mea_method == C_CELL_METHOD) {
	   volume = volume_cal1_cell(pattern_num, xcut, ycut, zcut, rx, ry, rz);
	}
	else {
	   volume = volume_cal1_anal(pattern_num, xcut, ycut, zcut, rx, ry, rz);
	}
	return(volume);
}

/****************************************/
/*					
/*	Procedure Name:  volume_cal1_cell
/*	  Return Value:  float
/*     Input Variables:  int	pattern_num	- pattern index
/*                       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
/*                       float rx, ry, rz       - units/voxels ratios in
/*                                                x, y, z axis, respectively
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  C_Mea_Vertices           mea_vertices
/* 				 /*These are pointers to the eight vertices
/*				   of the cell in one of the standard patterns */
/*     	   Description:  This routine calculates the volume inside a cell.
/*			 For each standard pattern, the volume is defined
/*			 by one or more tetrahedrons. Applying
/* 			 appropriate formula to compute the tetrahedral volume. 
 */


float volume_cal1_cell(pattern_num, xcut, ycut, zcut, rx, ry, rz)
int pattern_num;
float xcut[256][256];
float ycut[256][256];
float zcut[256][256];
float rx, ry, rz;
{
	extern C_Mea_Vertices		mea_vertices;

	float e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12;
	/* the intersection values of 12 edges of a cell */

	float t_volume;

	float sq_a, sq_b, sq_c, sq_r, sq_p, sq_q;
	
	switch(pattern_num) {
	case 1:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		t_volume = e1 * e4 * e9;
		break;
	case 2:
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		e10 = zcut[*mea_vertices.v2][*mea_vertices.v6];
		t_volume = (e4*e9 + e2*e10 + e2*e9)*rx;
		break;
	case 3:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		e2 = ycut[*mea_vertices.v3][*mea_vertices.v2];
		e3 = xcut[*mea_vertices.v3][*mea_vertices.v4];
		e12 = zcut[*mea_vertices.v3][*mea_vertices.v7];
		t_volume = (e1*e4*e9 + e2*e3*e12);
		break;
	case 4:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		e6 = ycut[*mea_vertices.v7][*mea_vertices.v6];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		t_volume = (e1*e4*e9 + e6*e7*e12);
		break;
	case 5:
		e1 = xcut[*mea_vertices.v2][*mea_vertices.v1];
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e8 = ycut[*mea_vertices.v5][*mea_vertices.v8];
		e9 = zcut[*mea_vertices.v5][*mea_vertices.v1];
		t_volume = ((e2+e6+e8)*rx*rz + e8*e9*rx + e2*e9*rx 
			   + e1*e2*(rz-e9));
		break;
	case 6:
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		e10 = zcut[*mea_vertices.v2][*mea_vertices.v6];
		e6 = ycut[*mea_vertices.v7][*mea_vertices.v6];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		t_volume = (e6*e7*e12 + (e4*e9 + e2*e10 + e2*e9)*rx);
		break;
	case 7:
		e1 = xcut[*mea_vertices.v2][*mea_vertices.v1];
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e10 = zcut[*mea_vertices.v2][*mea_vertices.v6];
		/* triangle #1 */
		e3 = xcut[*mea_vertices.v4][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v4][*mea_vertices.v1];
		e11 = zcut[*mea_vertices.v4][*mea_vertices.v8];
		/* triangle #2 */
		e6 = ycut[*mea_vertices.v7][*mea_vertices.v6];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		/* triangle #3 */
		t_volume = (e1*e2*e10 + e3*e4*e11 + e6*e7*e12);
		break;
	case 8:
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e8 = ycut[*mea_vertices.v5][*mea_vertices.v8];
		t_volume = (e4+e2+e6+e6+e8+e4)*rx*rz;
		break;
	case 9:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e10 = zcut[*mea_vertices.v6][*mea_vertices.v2];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e7 = xcut[*mea_vertices.v8][*mea_vertices.v7];
		e11 = zcut[*mea_vertices.v8][*mea_vertices.v4];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		t_volume = (e7*ry*(rz+e11) + e7*e4*(rz-e11) + e1*(ry+e4)*rz +
			    (rx-e7)*e6*rz + (rx-e1)*e6*e10 + rx*ry*rz);
		break;
	case 10:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
 		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		e3 = xcut[*mea_vertices.v4][*mea_vertices.v3];
		e11 = zcut[*mea_vertices.v4][*mea_vertices.v8];
		e10 = zcut[*mea_vertices.v6][*mea_vertices.v2];
		e5 = xcut[*mea_vertices.v6][*mea_vertices.v5];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		t_volume = (e3*e11 + e1*e9 + e3*e9 
			  + e5*e10 + e7*e12 + e5*e12)*ry;
		break;
	case 11:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e10 = zcut[*mea_vertices.v6][*mea_vertices.v2];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		e8 = ycut[*mea_vertices.v5][*mea_vertices.v8];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
/*
		sq_a = C_sq( e8 );
		sq_b = C_sq( e1 ) + C_sq( e8 ) + C_sq( rz );
		sq_c = C_sq( rx ) + C_sq( ry-e8 ) + C_sq( e12 );
		sq_r = C_sq( e1 ) + C_sq( rz );
		sq_p = C_sq( rx-e1 ) + C_sq( ry ) + C_sq( rz-e12 );
		sq_q = C_sq( rx ) + C_sq( ry ) + C_sq( e12 );
*/
		sq_a = ( e8*e8 );
		sq_b = ( e1*e1 ) + ( e8*e8 ) + ( rz*rz );
		sq_c = ( rx*rx ) + ( ry-e8 )*(ry-e8) + ( e12*e12 );
		sq_r = ( e1*e1 ) + ( rz*rz );
		sq_p = ( rx-e1 )*(rx-e1) + ( ry*ry ) + ( rz-e12 )*(rz-e12);
		sq_q = ( rx*rx ) + ( ry*ry ) + ( e12*e12 );

          	t_volume = ((float)(sqrt((double)(sq_a*sq_b*sq_q + sq_c*sq_r*sq_p))))/2.; 
/*
          	t_volume = (float)(sqrt(sq_a*sq_b*sq_q + sq_c*sq_r*sq_p)/2.); 
*/
		
		t_volume = t_volume +
			   ((rx+e7)*ry*e12 + (rx-e7)*e8*e12 + (rx-e1)*ry*e10 +
			    e1*(e4+e8)*rz + rx*ry*rz);
		break;
	case 12:
		e1 = xcut[*mea_vertices.v2][*mea_vertices.v1];
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e8 = ycut[*mea_vertices.v5][*mea_vertices.v8];
		e9 = zcut[*mea_vertices.v5][*mea_vertices.v1];
		e3 = xcut[*mea_vertices.v4][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v4][*mea_vertices.v1];
		e11 = zcut[*mea_vertices.v4][*mea_vertices.v8];
		t_volume = ((e2+e6+e8)*rx*rz + e8*e9*rx + e2*e9*rx 
			   + e1*e2*(rz-e9) + e3*e4*e11);
		break;
	case 13:
		e1 = xcut[*mea_vertices.v1][*mea_vertices.v2];
		e4 = ycut[*mea_vertices.v1][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v1][*mea_vertices.v5];
		/* triangle #1 */
		e2 = ycut[*mea_vertices.v3][*mea_vertices.v2];
		e3 = xcut[*mea_vertices.v3][*mea_vertices.v4];
		e12 = zcut[*mea_vertices.v3][*mea_vertices.v7];
		/* triangle #2 */
		e10 = zcut[*mea_vertices.v6][*mea_vertices.v2];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e5 = xcut[*mea_vertices.v6][*mea_vertices.v5];
		/* triangle #3 */
		e7 = xcut[*mea_vertices.v8][*mea_vertices.v7];
		e8 = ycut[*mea_vertices.v8][*mea_vertices.v5];
		e11 = zcut[*mea_vertices.v8][*mea_vertices.v4];
		/* triangle #4 */
		t_volume = (e1*e4*e9 + e2*e3*e12 + e5*e6*e10 + e7*e8*e11);
		break;
	case 14:
		e1 = xcut[*mea_vertices.v2][*mea_vertices.v1];
		e2 = ycut[*mea_vertices.v2][*mea_vertices.v3];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e7 = xcut[*mea_vertices.v8][*mea_vertices.v7];
		e11 = zcut[*mea_vertices.v8][*mea_vertices.v4];
		e9 = zcut[*mea_vertices.v5][*mea_vertices.v1];
/*
		sq_a = C_sq( ry ) + C_sq( e11 );
		sq_b = C_sq( rx-e1 ) + C_sq( ry ) + C_sq( rz-e11 );
		sq_c = C_sq( rx ) + C_sq( ry-e6 ) + C_sq( e11 );
		sq_r = C_sq( rx-e1 ) + C_sq( rz );
		sq_p = C_sq( e1 ) + C_sq( e6 ) + C_sq( rz );
		sq_q = C_sq( rx ) + C_sq( e6 );
*/
		sq_a = ( ry*ry ) + ( e11*e11 );
		sq_b = ( rx-e1 )*(rx-e1) + ( ry*ry ) + ( rz-e11 )*(rz-e11);
		sq_c = ( rx*rx ) + ( ry-e6 )*(ry-e6) + ( e11*e11 );
		sq_r = ( rx-e1 )*(rx-e1) + ( rz*rz );
		sq_p = ( e1*e1 ) + ( e6*e6 ) + ( rz*rz );
		sq_q = ( rx*rx ) + ( e6*e6 );

          	t_volume = ((float)(sqrt((double)(sq_a*sq_b*sq_q + sq_c*sq_r*sq_p))))/2.; 
/*
          	t_volume = (float)(sqrt(sq_a*sq_b*sq_q + sq_c*sq_r*sq_p)/2.); 
*/
		
		t_volume = t_volume + 
			   (e7*(ry-e6)*e11 + (rx-e1)*ry*e9 + rx*ry*e11 +
			    rx*e6*rz + e1*(e6+e2)*rz);
		break;
	default:
		printf("error in data file: pattern_num is invalid.\n");
		break;
	}
	return(t_volume);
}

