/****************************************************************************/
/*                                                                          */
/*  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_parc_subdivide.c
 *		      Author: Rick Avila
 *			Date: 08/19/92
 *		 Description: PARC Subdivision Algorithm
 *	Modification History:
 *
 *		Who?		When?		Why?
 *	--------------------------------------------------------------------
 *
 */

#include <stdio.h>
#include <math.h>

#include "C_volvis.h"
#include "C_navigator.h"

# define C_ffloor( x )	(float)((int)(x))
# define C_fceil( x )	((x == (int)(x))?(x):(float)((int)(x+1)))

extern	C_World		world;
extern	C_View		view;

C_parc_subdivide( volume, level )
C_Volume	*volume;
int		level;
{
	int		i;
	C_ParcInfo	*p_info;
	C_ParcSubNode	*list_ptr;
	C_ParcSubNode	*temp_ptr;

	/* Check If Volume Has Already Been PARC Subdivided */
	if( volume->p_info )
	{
		/* Set List Pointer To Head Of List */
		list_ptr = volume->p_info->sub_list;

		/* Release All Nodes Previously Allocated */
		while( list_ptr )
		{
			temp_ptr = list_ptr;
			list_ptr = list_ptr->next;
			free( temp_ptr );
		}

		/* Free the Parc Information Structure */
		free( volume->p_info );

		volume->p_info = NULL;
	}

	p_info = (C_ParcInfo *)malloc( sizeof(C_ParcInfo) );
	if( !p_info )
	{
		C_error_message("Not Enough Memory To PARC Subdivide");
		volume->p_info = NULL;
		return;
	}

	p_info->sub_thresh_low = 0.0;
	p_info->sub_level = level;
	p_info->sub_count = 0;
	p_info->sub_list = NULL;

	volume->p_info = p_info;

	C_message("Starting PARC Subdivision");

	parc_subdivide( volume );

	C_message("PARC Subdivision Complete");

}

parc_subdivide( volume )
C_Volume	*volume;
{
	int		level;		/* Level Of Subdivision */
	int		divisions;	/* Divisions On A Dimension */
	float		x_div_size;	/* Size Of X Division */
	float		y_div_size;	/* Size Of Y Division */
	float		z_div_size;	/* Size Of Z Division */
	int		x, y, z;	/* Subvolume Indeces */
	int		xx, yy, zz;	/* Floating Point Position */
	int		xvoxels;	/* Data Voxel Reolution In X */
	int		yvoxels;	/* Data Voxel Reolution In Y */
	int		zvoxels;	/* Data Voxel Reolution In Z */
	int		vox_plane;	/* # Voxels In XY Plane Of Data */
	int		sub_plane;	/* # Subdiv In A Plane Of Data */
	int		status;		/* Subvolume Status */
	int		sub_count;	/* Number Of Surface Subvolumes */
	C_ParcInfo	*parc_info;	/* Ptr To the PARC Structure */
	unsigned int	index_count;	/* Index To A Sobvolume */
	C_FPosition	start;		/* Starting Coordinates Of Subvolume */
	C_FPosition	end;		/* Ending Coordinates Of Subvolume */
	C_FPosition	half;		/* Halfway Position */
	unsigned char	*sub_data;	/* Ptr To Start Of Subvolume Data */
	unsigned char	*sub_data_ptr;	/* Ptr To Subvolume Data */
	C_Voxel_8bit	*scalar_data;	/* Ptr To Volume Scalar Data */
	C_Voxel_8bit	value;		/* Voxel Value */
	float		threshold;	/* Low Voxel Threshold Value */
	C_ParcSubNode	*sub_list;	/* Ptr To List Of Subvolume Data */
	C_ParcSubNode	*sub_node;	/* Ptr To List Node */
	unsigned char	faces;
	char		msg[C_MAX_STRING];

	/* Get The Volume Pointer & Subdivision Level */
	parc_info = volume->p_info;
	level  = parc_info->sub_level;

	/* Error Check The Level Number */
	if( (level < 0) || (level > 8) )
	{
		C_error_message("Unreasonable PARC Subdivision Level");
		return;
	}

	/* Get The Threshold Value */
	switch( volume->data_type )
	{
	    case C_GEOMETRIC_DATA:
	      C_error_message("Can't subdivide GEOMETRIC object!");
	      return;
	      break; /* Unnecessary */

	    case C_SCALAR_DATA_8BIT:

	      /* Get Pointer To Scan Data Array */
	      scalar_data = volume->data.scalar_data_8bit->scalar_data;

	      /* Check For Proper Segmentation Type!!!! */
	      switch( volume->data.scalar_data_8bit->seg_type )
	      {
		case C_ISO_SURFACE:
		   threshold = volume->data.scalar_data_8bit->seg.isovalue;
		   break;
		case C_SCALAR_OPACITY:
		   threshold = 1;
		   break;
		default:
		   C_error_message("Internal Error: PARC SUBDIVISION");
	      }
	}

	/* Create A Single Subvolume If Subdivision Level == 1 */
	if( level == 0 )
	{
		faces = 0;
		C_Set_Neg_X_Face( faces );
		C_Set_Pos_X_Face( faces );
		C_Set_Neg_Y_Face( faces );
		C_Set_Pos_Y_Face( faces );
		C_Set_Neg_Z_Face( faces );
		C_Set_Pos_Z_Face( faces );

		/* Allocate A Node For Subvolume List */
		sub_node = (C_ParcSubNode *)malloc( sizeof(C_ParcSubNode) );
		if( !sub_node )
		{
		   C_error_message("Unable To Allocate Subvolume List Node");
		   return;
		}

		sub_node->index = 0;
		sub_node->subnode_faces = faces;
		sub_node->normal.x = 0.0;
		sub_node->normal.y = 0.0;
		sub_node->normal.z = 0.0;
		sub_node->next = NULL;

		/* Set The Low Threshold In PARC Structure */
		parc_info->sub_thresh_low = threshold;

		/* Set Subvolume Count */
		parc_info->sub_count = 1;

		/* Set Subvolume List Pointer */
		parc_info->sub_list = sub_node;

		sprintf( msg, "SubVolume Count: 1" );
		C_message( msg );

		return;
	}

	divisions = (int)pow( (2.0), (double)level );

	/* Allocate Memory To Store Subdivision Data */
	sub_data = (unsigned char *)malloc( divisions*divisions*divisions );
	if( !sub_data )
	{
		C_error_message("Unable To Allocate PARC Subdivision Data");
		return;
	}

	switch( volume->data_type )
	{
	    case C_SCALAR_DATA_8BIT:

	      xvoxels = volume->x_size_voxels;
	      yvoxels = volume->y_size_voxels;
	      zvoxels = volume->z_size_voxels;

	      /* Number Of Voxels In An XY Plane Of Volume Data */
	      vox_plane = xvoxels * yvoxels;

	      /* Number Of Subdivisions In A Plane Of Subvolumes */
	      sub_plane = divisions * divisions;

	      /* Set Pointer Into Subvolume Data Array */
	      sub_data_ptr = sub_data;

	      x_div_size = (float)(xvoxels-1)/(float)(divisions);
	      y_div_size = (float)(yvoxels-1)/(float)(divisions);
	      z_div_size = (float)(zvoxels-1)/(float)(divisions);

	      /* Loop Through Every Subvolume And Check For Object */
	      /* Surface Within Subvolume			   */
	      for( z=0; z < divisions; z++ )
	      for( y=0; y < divisions; y++ )
	      for( x=0; x < divisions; x++ )
	      {
		/* Compute Boundary Of Subvolume */
		start.x = (float)(x) * x_div_size;
		end.x = start.x + x_div_size;

		start.y = (float)(y) * y_div_size;
		end.y = start.y + y_div_size;

		start.z = (float)(z) * z_div_size;
		end.z = start.z + z_div_size;

		status = C_NOTHING_FOUND;

		for( zz = C_ffloor(start.z); zz <= C_fceil(end.z); zz++ )
		for( yy = C_ffloor(start.y); yy <= C_fceil(end.y); yy++ )
		for( xx = C_ffloor(start.x); xx <= C_fceil(end.x); xx++ )
		{
		   value = *(scalar_data + xx + yy*xvoxels + zz*vox_plane);

		   /* Check For High Value Voxel */
		   if( value >= threshold )
		   {
			if( status == C_NOTHING_FOUND )
			{
				status = C_HIGH_FOUND;
			}
			else if( status == C_LOW_FOUND )
			{
				status = C_SURFACE_FOUND;
				break;
			}
		   }
		   else /* Check For Low Value Voxel */
		   {
			if( status == C_NOTHING_FOUND )
			{
				status = C_LOW_FOUND;
			}
			else if( status == C_HIGH_FOUND )
			{
				status = C_SURFACE_FOUND;
				break;
			}
		   }
		}

		*(sub_data_ptr++) = status;
	      }

	      /* No Subvolumes In List Yet! */
	      sub_list = NULL;

	      sub_data_ptr = sub_data;

	      /* No Surface Subvolumes Counted Yet */
	      sub_count = 0;
	      index_count = 0;

	      /* Loop Through Every Subvolume And Check For Surfaces */
	      /* Present At Subvolume Extent */
	      for( z=0; z < divisions; z++ )
	      for( y=0; y < divisions; y++ )
	      for( x=0; x < divisions; x++ )
	      {

		faces = 0;

		if( *(sub_data_ptr) == C_HIGH_FOUND ||
		    *(sub_data_ptr) == C_SURFACE_FOUND )
		{
		   if ( x == 0 )
			C_Set_Neg_X_Face( faces );
		   
		   if ( y == 0 )
			C_Set_Neg_Y_Face( faces );
		   
		   if ( z == 0 )
			C_Set_Neg_Z_Face( faces );
		   
		   if ( x == (divisions-1) )
			C_Set_Pos_X_Face( faces );
		   
		   if ( y == (divisions-1) )
			C_Set_Pos_Y_Face( faces );
		   
		   if ( z == (divisions-1) )
			C_Set_Pos_Z_Face( faces );
		   
		   if (  *(sub_data_ptr-1) == C_LOW_FOUND )
			C_Set_Neg_X_Face( faces );

		   if (  *(sub_data_ptr+1) == C_LOW_FOUND )
			C_Set_Pos_X_Face( faces );

		   if (  *(sub_data_ptr-divisions) == C_LOW_FOUND )
			C_Set_Neg_Y_Face( faces );

		   if (  *(sub_data_ptr+divisions) == C_LOW_FOUND )
			C_Set_Pos_Y_Face( faces );

		   if (  *(sub_data_ptr-sub_plane) == C_LOW_FOUND )
			C_Set_Neg_Z_Face( faces );

		   if (  *(sub_data_ptr+sub_plane) == C_LOW_FOUND )
			C_Set_Pos_Z_Face( faces );
		}

		if( *(sub_data_ptr) == C_SURFACE_FOUND  || faces )
		{
		    /* Allocate A Node For Subvolume List */
		    sub_node = (C_ParcSubNode *)malloc( sizeof(C_ParcSubNode) );
		    if( !sub_node )
		    {
		      C_error_message("Unable To Allocate Subvolume List Node");
		      return;
		    }

		    sub_node->index = index_count;
		    sub_node->subnode_faces = faces;
		    sub_node->normal.x = 0.0;
		    sub_node->normal.y = 0.0;
		    sub_node->normal.z = 0.0;
		    sub_node->next = sub_list;
		    sub_list = sub_node;

		    /* Increment The Subvolume Count */
		    sub_count++;
		}

		/* Go To Next Subvolume */
		sub_data_ptr++;

		/* Increment Subvolume Index Count */
		index_count++;
	      }

	      sprintf( msg, "SubVolume Count: %d\n", sub_count );
	      C_message( msg );

	      /* Set The Low Threshold In PARC Structure */
	      parc_info->sub_thresh_low = threshold;

	      /* Set Subvolume Count */
	      parc_info->sub_count = sub_count;

	      /* Set Subvolume List Pointer */
	      parc_info->sub_list = sub_list;

	      /* Free Subvolume Data Array */
	      free( sub_data );
	}
	
}

