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


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

/****************************************/
/*					
/*	Procedure Name:  anal_area
/*	  Return Value:  float
/*     Input Variables:  
/*			 float c0,c1,c2,c3,d0,d1,d2,d3 
/*				- coeficiences of the trilinear interpolation function
/*			 float dx,dy,dz - size of a tiny area
/*			 float x0,x1,y0,y1,z0,z1 - boundary of the subarea
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine calculates the surface area inside
/*			 a subarea.
 */

float anal_area(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,x0,x1,y0,y1,z0,z1)
float c0,c1,c2,c3;
float d0,d1,d2,d3;
float dx,dy,dz;
float x0,x1,y0,y1,z0,z1;
{
	extern C_MeaInfo        mea_info;
	int s_count;	/* surface count */
	float surf;	/* surface as a floating point number */	
	float x, y, z;	/* x, y, z coordinate of tiny area */
	int loopx, loopy, loopz; /* loop variables */
	int stepx, stepy, stepz; /* maximum steps of looping */
	int low, high;	/* segmentation values */
	int value;	/* trilinear interpolated intensity value */
	int flag;
	int x_flag;

	low = mea_info.low;
	high = mea_info.high;
	stepx = (int) ((x1-x0)/dx);
	stepy = (int) ((y1-y0)/dy);
	stepz = (int) ((z1-z0)/dz);

	s_count = 0;
	for (loopz=0; loopz<=stepz; loopz++) {
	  z = z0 + loopz * dz;
	  for (loopy=0; loopy<=stepy; loopy++) {
	     y = y0 + loopy * dy;
	     value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,x0,y,z);	
	     if (C_In_Range(value, low, high))
	        x_flag = 1;
	     else
	        x_flag = 0;
	     for (loopx=1; loopx<=stepx; loopx++) {
	        x = x0 + loopx * dx;
		value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,x,y,z);	
		if (C_In_Range(value, low, high)) {
		    	if (x_flag == 0) {
			    flag = 1;
			}
			else {
		            flag = 0;
		        }
			value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,
						   x,y,z-dz);	
			if (!(C_In_Range(value, low, high)))
		            flag += 1;
			value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,
						   x,y,z+dz);	
			if (!(C_In_Range(value, low, high)))
		            flag += 1;
			value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,
						   x,y-dy,z);	
			if (!(C_In_Range(value, low, high)))
		            flag += 1;
			value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,
						   x,y+dy,z);	
			if (!(C_In_Range(value, low, high)))
		            flag += 1;
			value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,
						   x+dx,y,z);	
			if (!(C_In_Range(value, low, high)))
		            flag += 1;
			s_count += flag;
			x_flag = 1;
		}
		else {
			x_flag = 0;
		}
	     }
	   }
	}
	surf = ((float) s_count);
	return(surf);
}

/****************************************/
/*					
/*	Procedure Name:  anal_vol
/*	  Return Value:  float
/*     Input Variables:  
/*			 float c0,c1,c2,c3,d0,d1,d2,d3 
/*				- coeficiences of the trilinear interpolation function
/*			 float dx,dy,dz - size of a tiny area
/*			 float x0,x1,y0,y1,z0,z1 - boundary of the subarea
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine calculates the volume inside
/*			 a subarea.
 */

float anal_vol(c0,c1,c2,c3,d0,d1,d2,d3,dx,dy,dz,x0,x1,y0,y1,z0,z1) 
float c0,c1,c2,c3;
float d0,d1,d2,d3;
float dx,dy,dz;
float x0,x1,y0,y1,z0,z1;
{
	extern C_MeaInfo        mea_info;
	int v_count;	/* volume count */
	float vol;	
	float x, y, z;	/* x, y, z coordinate of tiny area */
	int loopx, loopy, loopz; /* loop variables */
	int stepx, stepy, stepz; /* maximum steps of looping */
	int low, high;	/* segmentation values */
	int value;	/* trilinear interpolated intensity value */

	low = mea_info.low;
	high = mea_info.high;
	stepx = (int) ((x1-x0)/dx);
	stepy = (int) ((y1-y0)/dy);
	stepz = (int) ((z1-z0)/dz);

	v_count = 0;
	if (mea_info.complementary == 0) {
	  /* count the number of tiny areas inside the object for this case */
	  for (loopz=0; loopz<=stepz; loopz++) {
	    z = z0 + loopz * dz;
	    for (loopy=0; loopy<=stepy; loopy++) {
	      y = y0 + loopy * dy;
	      for (loopx=0; loopx<=stepx; loopx++) {
	        x = x0 + loopx * dx;
		value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,x,y,z);	
		if (C_In_Range(value, low, high)) {
		    v_count += 1;
		}
	      } 
	    }
	  }
	}
	else {
	  /* count the number of tiny areas outside the object for this case */
	  for (loopz=0; loopz<=stepz; loopz++) {
	    z = z0 + loopz * dz;
	    for (loopy=0; loopy<=stepy; loopy++) {
	      y = y0 + loopy * dy;
	      for (loopx=0; loopx<=stepx; loopx++) {
	        x = x0 + loopx * dx;
		value = C_Mea_Trilin(c0,c1,c2,c3,d0,d1,d2,d3,x,y,z);	
		if (!(C_In_Range(value, low, high))) {
		    v_count += 1;
		}
	      } 
	    }
	  }
	}
	/* since the final volume will be devided by 6. due to the other 
	   part of computation, we do mutiply 6. here.
	 */
	vol = ((float) v_count) * dx * dy * dz * 6.;
	return(vol);
}
