/****************************************************************************/
/*                                                                          */
/*  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_area_cal2.c
 *	      Author:  Hui Chen
 *		Date:  Feb. 9, 1991
 * 	 Description:  This file contains area_cal2(),
 *		       which calculates the surface area
 *		       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:  area_cal2
/*	  Return Value:  float
/*     Input Variables:  int	pattern_num	- pattern index
/*			 float rx, ry, rz	- units/voxels ratios in 
/*						  x, y, z axis, respectively
/*			 low, high
/*    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 surface area inside a cell.
/*			 For each standard pattern, the surface area is defined
/*			 by one or more triangls. For each triangle, two edge 
/*			 vectors are specified using left hand coordinates. 
/*			 By applying appropriate formula to compute the 
/*			 triangle area. 
/*			 See detail in the project report, the numbering systems/*			 are the same.
 */


float area_cal2(low, high, pattern_num, rx, ry, rz)
int low, high;
int pattern_num;
float rx, ry, rz;
{
	float compute_cut();

	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 tri_area;
	
	switch(pattern_num) {
	case 1:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		tri_area = triangle((-e1),e4,0.,(-e1),0.,e9);
		break;
	case 2:
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		e10 = compute_cut(low, high, rz, *mea_vertices.v2, *mea_vertices.v6);
		/* triangle #1 */
		tri_area = triangle(0.,(-e4),e9,rx,(e2-e4),0.);
		/* triangle #2 */
		tri_area = tri_area + triangle(0.,e2,(-e10),(-rx),0.,(e9-e10));
		break;
	case 3:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		/* triangle #1 */
		tri_area = triangle((-e1),e4,0.,(-e1),0.,e9);
		e2 = compute_cut(low, high, ry, *mea_vertices.v3, *mea_vertices.v2);
		e3 = compute_cut(low, high, rx, *mea_vertices.v3, *mea_vertices.v4);
		e12 = compute_cut(low, high, rz, *mea_vertices.v3, *mea_vertices.v7);
		/* triangle #2 */
		tri_area = tri_area + triangle(e3,(-e2),0.,e3,0.,e12);
		break;
	case 4:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		/* triangle #1 */
		tri_area = triangle((-e1),e4,0.,(-e1),0.,e9);
		e6 = compute_cut(low, high, ry, *mea_vertices.v7, *mea_vertices.v6);
		e7 = compute_cut(low, high, rx, *mea_vertices.v7, *mea_vertices.v8);
		e12 = compute_cut(low, high, rz, *mea_vertices.v7, *mea_vertices.v3);
		/* triangle #2 */
		tri_area = tri_area + triangle(e7,(-e6),0.,e7,0.,(-e12));
		break;
	case 5:
		e1 = compute_cut(low, high, rx, *mea_vertices.v2, *mea_vertices.v1);
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e8 = compute_cut(low, high, ry, *mea_vertices.v5, *mea_vertices.v8);
		e9 = compute_cut(low, high, rz, *mea_vertices.v5, *mea_vertices.v1);
		/* triangle #1 */
		tri_area = triangle(e1,e2,0.,(e1-rx),0.,(rz-e9));
		/* triangle #2 */
		tri_area = tri_area + triangle(0.,e8,e9,rx,e2,(e9-rz));
		/* triangle #3 */
		tri_area = tri_area + triangle(0.,(e2-e6),(-rz),(-rx),(e8-e6),0.);
		break;
	case 6:
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		e10 = compute_cut(low, high, rz, *mea_vertices.v2, *mea_vertices.v6);
		/* triangle #1 */
		tri_area = triangle(0.,(-e4),e9,rx,(e2-e4),0.);
		/* triangle #2 */
		tri_area = tri_area + triangle(0.,e2,(-e10),(-rx),0.,(e9-e10));
		e6 = compute_cut(low, high, ry, *mea_vertices.v7, *mea_vertices.v6);
		e7 = compute_cut(low, high, rx, *mea_vertices.v7, *mea_vertices.v8);
		e12 = compute_cut(low, high, rz, *mea_vertices.v7, *mea_vertices.v3);
		/* triangle #3 */
		tri_area = tri_area + triangle(e7,(-e6),0.,e7,0.,(-e12));
		break;
	case 7:
		e1 = compute_cut(low, high, rx, *mea_vertices.v2, *mea_vertices.v1);
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e10 = compute_cut(low, high, rz, *mea_vertices.v2, *mea_vertices.v6);
		/* triangle #1 */
		tri_area = triangle(e1,e2,0.,e1,0.,e10);
		e3 = compute_cut(low, high, rx, *mea_vertices.v4, *mea_vertices.v3);
		e4 = compute_cut(low, high, ry, *mea_vertices.v4, *mea_vertices.v1);
		e11 = compute_cut(low, high, rz, *mea_vertices.v4, *mea_vertices.v8);
		/* triangle #2 */
		tri_area = tri_area + triangle(e3,e4,0.,0.,e4,e11);
		e6 = compute_cut(low, high, ry, *mea_vertices.v7, *mea_vertices.v6);
		e7 = compute_cut(low, high, rx, *mea_vertices.v7, *mea_vertices.v8);
		e12 = compute_cut(low, high, rz, *mea_vertices.v7, *mea_vertices.v3);
		/* triangle #3 */
		tri_area = tri_area + triangle(e7,(-e6),0.,e7,0.,(-e12));
		break;
	case 8:
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e8 = compute_cut(low, high, ry, *mea_vertices.v5, *mea_vertices.v8);
		/* triangle #1 */
		tri_area = triangle(0.,(e6-e2),rz,(-rx),(e4-e2),0.);
		/* triangle #2 */
		tri_area = tri_area + triangle(0.,(e4-e8),(-rz),rx,(e6-e8),0.);
		break;
	case 9:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e10 = compute_cut(low, high, rz, *mea_vertices.v6, *mea_vertices.v2);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e7 = compute_cut(low, high, rx, *mea_vertices.v8, *mea_vertices.v7);
		e11 = compute_cut(low, high, rz, *mea_vertices.v8, *mea_vertices.v4);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		/* triangle #1 */
		tri_area = triangle((e1-rx),0.,(e10-rz),0.,e6,e10);
		/* triangle #2 */
		tri_area = tri_area + triangle((e7-rx),(ry-e6),0.,(e1-rx),(-e6),(-rz));
		/* triangle #3 */
		tri_area = tri_area + triangle((-e1),e4,0.,(e7-e1),ry,rz);
		/* triangle #4 */
		tri_area = tri_area + triangle(0.,(ry-e4),(rz-e11),e7,(ry-e4),rz);
		break;
	case 10:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
 		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		e3 = compute_cut(low, high, rx, *mea_vertices.v4, *mea_vertices.v3);
		e11 = compute_cut(low, high, rz, *mea_vertices.v4, *mea_vertices.v8);
		e10 = compute_cut(low, high, rz, *mea_vertices.v6, *mea_vertices.v2);
		e5 = compute_cut(low, high, rx, *mea_vertices.v6, *mea_vertices.v5);
		e12 = compute_cut(low, high, rz, *mea_vertices.v7, *mea_vertices.v3);
		e7 = compute_cut(low, high, rx, *mea_vertices.v7, *mea_vertices.v8);
		/* triangle #1 */
		tri_area = triangle((-e1),0.,e9,(e3-e1),ry,0.);
		/* triangle #2 */
		tri_area = tri_area + triangle(e3,0.,(-e11),0.,(-ry),(e9-e11));
		/* triangle #3 */
		tri_area = tri_area + triangle((-e5),0.,e10,0.,ry,(e12-e10));
		/* triangle #4 */
		tri_area = tri_area + triangle(e7,0.,(-e12),(e5-e7),(-ry),0.);
		break;
	case 11:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e10 = compute_cut(low, high, rz, *mea_vertices.v6, *mea_vertices.v2);
		e12 = compute_cut(low, high, rz, *mea_vertices.v7, *mea_vertices.v3);
		e7 = compute_cut(low, high, rx, *mea_vertices.v7, *mea_vertices.v8);
		e8 = compute_cut(low, high, ry, *mea_vertices.v5, *mea_vertices.v8);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		/* triangle #1 */
		tri_area = triangle(e1,(-e4),0.,0.,(e8-e4),rz);
		/* triangle #2 */
		tri_area = tri_area + triangle((e1-rx),0.,(e10-rz),0.,ry,(e12-e10));
		/* triangle #3 */
		tri_area = tri_area + triangle((e7-rx),(e8-ry),0.,e7,0.,(-e12));
		/* triangle #4 */
		tri_area = tri_area + triangle(e1,(-e8),(-rz),rx,(ry-e8),(-e12));
		break;
	case 12:
		e1 = compute_cut(low, high, rx, *mea_vertices.v2, *mea_vertices.v1);
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e8 = compute_cut(low, high, ry, *mea_vertices.v5, *mea_vertices.v8);
		e9 = compute_cut(low, high, rz, *mea_vertices.v5, *mea_vertices.v1);
		/* triangle #1 */
		tri_area = triangle(e1,e2,0.,(e1-rx),0.,(rz-e9));
		/* triangle #2 */
		tri_area = tri_area + triangle(0.,e8,e9,rx,e2,(e9-rz));
		/* triangle #3 */
		tri_area = tri_area + triangle(0.,(e2-e6),(-rz),(-rx),(e8-e6),0.);
		e3 = compute_cut(low, high, rx, *mea_vertices.v4, *mea_vertices.v3);
		e4 = compute_cut(low, high, ry, *mea_vertices.v4, *mea_vertices.v1);
		e11 = compute_cut(low, high, rz, *mea_vertices.v4, *mea_vertices.v8);
		/* triangle #4 */
		tri_area = tri_area + triangle(e3,e4,0.,0.,e4,e11);
		break;
	case 13:
		e1 = compute_cut(low, high, rx, *mea_vertices.v1, *mea_vertices.v2);
		e4 = compute_cut(low, high, ry, *mea_vertices.v1, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v1, *mea_vertices.v5);
		/* triangle #1 */
		tri_area = triangle((-e1),e4,0.,(-e1),0.,e9);
		e2 = compute_cut(low, high, ry, *mea_vertices.v3, *mea_vertices.v2);
		e3 = compute_cut(low, high, rx, *mea_vertices.v3, *mea_vertices.v4);
		e12 = compute_cut(low, high, rz, *mea_vertices.v3, *mea_vertices.v7);
		/* triangle #2 */
		tri_area = tri_area + triangle(e3,(-e2),0.,e3,0.,e12);
		e10 = compute_cut(low, high, rz, *mea_vertices.v6, *mea_vertices.v2);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e5 = compute_cut(low, high, rx, *mea_vertices.v6, *mea_vertices.v5);
		e7 = compute_cut(low, high, rx, *mea_vertices.v8, *mea_vertices.v7);
		e8 = compute_cut(low, high, ry, *mea_vertices.v8, *mea_vertices.v5);
		e11 = compute_cut(low, high, rz, *mea_vertices.v8, *mea_vertices.v4);
		/* triangle #3 */
		tri_area = tri_area + triangle(e5,e6,0.,e5,0.,(-e10));
		/* triangle #4 */
		tri_area = tri_area + triangle(e7,e8,0.,0.,e8,(-e11));
		break;
	case 14:
		e1 = compute_cut(low, high, rx, *mea_vertices.v2, *mea_vertices.v1);
		e2 = compute_cut(low, high, ry, *mea_vertices.v2, *mea_vertices.v3);
		e6 = compute_cut(low, high, ry, *mea_vertices.v6, *mea_vertices.v7);
		e7 = compute_cut(low, high, rx, *mea_vertices.v8, *mea_vertices.v7);
		e11 = compute_cut(low, high, rz, *mea_vertices.v8, *mea_vertices.v4);
		e9 = compute_cut(low, high, rz, *mea_vertices.v5, *mea_vertices.v1);
		/* triangle #1 */
		tri_area = triangle((-e1),(-e2),0.,0.,(e6-e2),rz);
		/* triangle #2 */
		tri_area = tri_area + triangle((rx-e1),0.,(e9-rz),0.,ry,(e11-e9));
		/* triangle #3 */
		tri_area = tri_area + triangle((rx-e7),(e6-ry),0.,(-e7),0.,(-e11));
		/* triangle #4 */
		tri_area = tri_area + triangle((-e1),(-e6),(-rz),(-rx),(ry-e6),(-e11));
		break;
	default:
		printf("error in data file: pattern_num is invalid.\n");
		break;
	}
	return(tri_area);
}


float compute_cut(low, high, r, a, b)
int low, high;
float r;
int a, b;
{
	float cut;

	if (b<low)                 
           cut = (((float)(a-low))/((float)(a-b)))*r; 
        else        
           cut = (((float)(high-a))/((float)(b-a)))*r; 

	return (cut);
}


