/****************************************************************************/
/*                                                                          */
/*  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_anal.c
 *	      Author:  Hui Chen
 *		Date:  Feb. 9, 1992
 * 	 Description:  This file contains volume_cal1_anal(),
 * 		       It is device-independent.
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */


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

/****************************************/
/*					
/*	Procedure Name:  volume_cal1_anal
/*	  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.
/*			 You may notice that I used range [0, rx], [0, ry],
/*                       [0, rz] instead of [0, 1], [0, 1], [0, 1] in the
/*                       trilinear interpolation function because this way
/*                       this method will work for any case of cell unit size.
 */


float volume_cal1_anal(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 anal_vol();
	extern C_MeaInfo        mea_info;
	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 volume;

	int a, b, c, d, e, f, g, h;
	float c0, c1, c2, c3, d0, d1, d2, d3;
	float dx, dy, dz;

	a = *mea_vertices.v1;
	b = *mea_vertices.v2;
	c = *mea_vertices.v4;
	d = *mea_vertices.v3;
	e = *mea_vertices.v5;
	f = *mea_vertices.v6;
	g = *mea_vertices.v8;
	h = *mea_vertices.v7;

	/* calculate the coeficiences of the trilinear interpolation function */
	c0 = a * rx * ry;
	c1 = (b - a) * ry;
	c2 = (c - a) * rx;
	c3 = a - b - c + d;
	
	d0 = (e - a) * rx * ry;
	d1 = (a - b - e + f) * ry;
	d2 = (a - c - e + g) * rx;
	d3 = e - f - g + h - a + b + c - d;
	
	dx = rx /(mea_info.steps);
	dy = ry /(mea_info.steps);
	dz = rz /(mea_info.steps);

	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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,e1,0.0,e4,0.0,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];
		if (e2 > e4) {
		    if (e10 > e9) {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e2,0.0,e10);
		    }
		    else {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e2,0.0,e9);
		    }
		}
		else {
		    if (e10 > e9) {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e4,0.0,e10);
		    }
		    else {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e4,0.0,e9);
		    }
		}
		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];
		if (((e1+e3) > rx) && ((e4+e2) > ry)) {
		    if (e12 > e9) {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
		    			 0.0,rx,0.0,ry,0.0,e12);
		    }
		    else {
			volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
		    			 0.0,rx,0.0,ry,0.0,e9);
		    }
		}
		else {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
		    		 0.0,e1,0.0,e4,0.0,e9);
		    volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				  rx-e3,rx,ry-e2,ry,0.0,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];
		if (((e1+e7) > rx) && ((e4+e6) > ry) && ((e9+e12) > rz)) {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		}
		else {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
		    		 0.0,e1,0.0,e4,0.0,e9);
		    volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				  rx-e7,rx,ry-e6,ry,rz-e12,rz);
		}
		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];
		if (e6 > e8) {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				     0.0,rx,0.0,e6,0.0,rz);
		}
		else {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				     0.0,rx,0.0,e8,0.0,rz);
		}
		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];
		if ((((e6+e2) > ry) || ((e6+e4) > ry)) && 
		    (((e12+e10) >rz) || ((e12+e9) > rz))) {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		}
		else {
		    if (e2 > e4) {
		        if (e10 > e9) {
		    	    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e2,0.0,e10);
		        }
		        else {
		    	    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e2,0.0,e9);
		        }
		    }
		    else {
		        if (e10 > e9) {
			    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e4,0.0,e10);
		        }
		        else {
		    	    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 	 0.0,rx,0.0,e4,0.0,e9);
		        }
		    }
		    volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				  rx-e7,rx,ry-e6,ry,rz-e12,rz);
		}
		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];
		e3 = xcut[*mea_vertices.v4][*mea_vertices.v3];
		e4 = ycut[*mea_vertices.v4][*mea_vertices.v1];
		e11 = zcut[*mea_vertices.v4][*mea_vertices.v8];
		e6 = ycut[*mea_vertices.v7][*mea_vertices.v6];
		e7 = xcut[*mea_vertices.v7][*mea_vertices.v8];
		e12 = zcut[*mea_vertices.v7][*mea_vertices.v3];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,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];
		if (e6 > e8) {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				     0.0,rx,0.0,e6,0.0,rz);
		    if ((e4+e6) > ry) {
		        volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				      0.0,e3,ry-e6,ry,0.0,e11);
		    }
		    else {
		        volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				      0.0,e3,ry-e4,ry,0.0,e11);
		    }
		}
		else {
		    volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				     0.0,rx,0.0,e8,0.0,rz);
		    if ((e4+e8) > ry) {
		        volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				      0.0,e3,ry-e8,ry,0.0,e11);
		    }
		    else {
		        volume += anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				      0.0,e3,ry-e4,ry,0.0,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];
		e2 = ycut[*mea_vertices.v3][*mea_vertices.v2];
		e3 = xcut[*mea_vertices.v3][*mea_vertices.v4];
		e12 = zcut[*mea_vertices.v3][*mea_vertices.v7];
		e10 = zcut[*mea_vertices.v6][*mea_vertices.v2];
		e6 = ycut[*mea_vertices.v6][*mea_vertices.v7];
		e5 = xcut[*mea_vertices.v6][*mea_vertices.v5];
		e7 = xcut[*mea_vertices.v8][*mea_vertices.v7];
		e8 = ycut[*mea_vertices.v8][*mea_vertices.v5];
		e11 = zcut[*mea_vertices.v8][*mea_vertices.v4];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		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];
		volume = anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,
				 0.0,rx,0.0,ry,0.0,rz);
		break;
	default:
		printf("error in data file: pattern_num is invalid.\n");
		break;
	}
	return(volume);
}
