/****************************************************************************/
/*                                                                          */
/*  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_v_voxel.c
 *	      Author:  Hui Chen
 *		Date:  Dec. 7, 1991
 * 	 Description:  This file contains routine: v_voxel_method(),
 *		       which calculates the volume by counting 
 *		       the number of voxels in the object.	
 *		       It is device_independent.
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */

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

/****************************************/
/*					
/*	Procedure Name:  v_voxel_method
/*	  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 volume by counting 
/*		         the number of voxels in the object.	
 */


float v_voxel_method(data,xsize,ysize,zsize,low,high,rx,ry,rz)
C_Voxel_8bit	*data;
int xsize,ysize,zsize,low,high;
float rx, ry, rz;
{
	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 */
	float	vol_count;	/* count of voxels on the volume  */
	int	inside_count;	/* count of voxels inside the surface */
	int	onsurf_count;	/* count of voxels on the surface  */

	/* set up the increaments */

	zstep = xsize*ysize;
	ystep = xsize;
	xstep = 1;

	/* initialize the count to 0 */

	vol_count = 0;
	inside_count  = 0;
	onsurf_count  = 0;

	/* now start counting the voxels in the object */

	for (z=1; z<(zsize-1); z++) {
	    for (y=1; y<(ysize-1); y++)
	    {
		/* initialize dptr for this row */
		dptr = data + (z*zstep) + (y*ystep) + xstep;
		for (x=1; x<(xsize-1); x++)
		{
		    if (C_In_Range( *dptr, low, high)) {
		       if ( C_On_Surface(dptr,low,high,
                                         xstep,ystep,zstep)) {
			  /* count the surface voxel as half */
			  vol_count += 0.5;
			  onsurf_count ++;
		       }
		       else {
		          vol_count += 1.0;
			  inside_count ++;
		       }
		    }
		    dptr += xstep;
		}
	    }
	}
	return(vol_count);
}

/****************************************/
/*					
/*	Procedure Name:  v_seed_voxel
/*	  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
/*			 C_MeaVoxelList  *mea_surface_voxels[] - surface list
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none 
/*     	   Description:  This routine calculates the volume by counting 
/*		         the number of voxels in the object.	
/*			 See my project report for detail.
 */

float v_seed_voxel(mea_surface_voxels, data,xsize,ysize,zsize,low,high,rx,ry,rz)
C_MeaVoxelList  *mea_surface_voxels[C_MEA_HASH_MAX];
C_Voxel_8bit  *data;
int      xsize,ysize,zsize;
int      low, high;
float    rx, ry, rz;
{
	C_MeaVoxelList  *vptr, *vptr0, *vptr1;
	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	*/
	float	vol_count;	/* count of voxels on the volume  */
	int	inside_count;	/* count of voxels inside the surface */
	int	onsurf_count;	/* count of voxels on the surface  */
	int	index;
	int	hash_index;
	int 	flag;
	int     i;

	/* set up the increaments */

	zstep = xsize*ysize;
	ystep = xsize;
	xstep = 1;

	/* initialize the count to 0 */

	vol_count = 0;
	inside_count  = 0;
	onsurf_count  = 0;

	/* now start counting the voxels in the object */

	for (i=0; i<C_MEA_HASH_MAX; i++) {
            vptr = mea_surface_voxels[i];
            while (vptr != NULL) {
		dptr = data + vptr->voxel.index;
                vptr0 = vptr;
                vptr = vptr->next_voxel;
		free(vptr0);
		if (C_LEFT_Face(dptr, low, high, xstep)) {
		    vol_count += 0.5;
		    onsurf_count ++;
		    dptr ++;
		    while (C_In_Range( *(dptr), low, high)) {	
		       if ( C_On_Surface(dptr,low,high,
                                         xstep,ystep,zstep)) {
			  vol_count += 0.5;
			  onsurf_count ++;
			  index = dptr - data;
			  hash_index = index % C_MEA_HASH_MAX; 
            		  vptr0 = mea_surface_voxels[hash_index];
            		  vptr1 = mea_surface_voxels[hash_index];
			  flag = 0;
			  while ((vptr1 != NULL) && (flag == 0)){
				if (vptr1->voxel.index == index) {
				    flag = 1;
				    if (vptr0 == vptr1) {
				        mea_surface_voxels[hash_index] = 
					mea_surface_voxels[hash_index]->next_voxel;
				    }
				    else {
				        vptr0->next_voxel = vptr1->next_voxel;
				    }
				    free(vptr1);
				}
				else { 
				    vptr0 = vptr1;
				    vptr1 = vptr1->next_voxel;
				}
			  }
		       }
		       else {
		          vol_count += 1.0;
			  inside_count ++;
		       }
		    dptr ++;
		    }
		}
		else if (C_RIGHT_Face(dptr, low, high, xstep)) {
		    vol_count += 0.5;
		    onsurf_count ++;
		    dptr --;
		    while (C_In_Range( *(dptr), low, high)) {	
		       if ( C_On_Surface(dptr,low,high,
                                         xstep,ystep,zstep)) {
			  vol_count += 0.5;
			  onsurf_count ++;
			  index = dptr - data;
			  hash_index = index % C_MEA_HASH_MAX; 
            		  vptr0 = mea_surface_voxels[hash_index];
            		  vptr1 = mea_surface_voxels[hash_index];
			  flag = 0;
			  while ((vptr1 != NULL) && (flag == 0)) {
				if (vptr1->voxel.index == index) {
				    flag = 1;
				    if (vptr0 == vptr1) {
				        mea_surface_voxels[hash_index] = 
					mea_surface_voxels[hash_index]->next_voxel;
				    }
				    else {
				        vptr0->next_voxel = vptr1->next_voxel;
				    }
				    free(vptr1);
				}
				else { 
				    vptr0 = vptr1;
				    vptr1 = vptr1->next_voxel;
				}
			  }
		       }
		       else {
		          vol_count += 1.0;
			  inside_count ++;
		       }
		    dptr --;
		    }
		}
	    }
	}
	return(vol_count);
}
