/****************************************************************************/
/*                                                                          */
/*  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_tracking.c
 *	      Author:  Hui Chen
 *		Date:  9/23/1992
 * 	 Description:  This file contains routines implementing 
 *		       surface tracking algorithm.
 *		       This is device_independent.
 *
 *
 * Modification History:
 * 
 *	who?		when?		why?
 * -----------------------------------------------
 *
 */


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

/****************************************/
/*					
/*	Procedure Name:  mea_valid_seed
/*	  Return Value:  int
/*     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
/*			 int		seedx	- x coordinate of seed point
/*			 int		seedy	- y coordinate of seed point
/*			 int		seedz	- z coordinate of seed point
/*    Output Variables:  none
/*    Update Variables:  none
/*    Global Variables:  none
/*     	   Description:  This routine checks if the seed point is on the 
/*			 surface of the object. 
 */

int mea_valid_seed(data,xsize,ysize,zsize,low,high,seedx, seedy, seedz)
C_Voxel_8bit  *data;
int	 xsize,ysize,zsize;
int      low, high;
int	 seedx, seedy, seedz;
{
	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	*/

	/* set up the increaments */

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

	dptr = data + (seedz*zstep) + (seedy*ystep) + (seedx * xstep);

	if (!(C_In_Range(*dptr, low, high))) {
	    printf("Seed is NOT in range!!!\n");
	    return(0);
	}

/*
	if (!(C_In_Range(*(dptr-xstep), low, high))) {
	    printf("Seed is on C_MEA_LEFT.\n");
	}
	if (!(C_In_Range(*(dptr+xstep), low, high))) {
	    printf("Seed is on C_MEA_RIGHT.\n");
	}
	if (!(C_In_Range(*(dptr-ystep), low, high))) {
	    printf("Seed is on C_MEA_DOWN.\n");
	}
	if (!(C_In_Range(*(dptr+ystep), low, high))) {
	    printf("Seed is on C_MEA_UP.\n");
	}
	if (!(C_In_Range(*(dptr-zstep), low, high))) {
	    printf("Seed is on C_MEA_FRONT.\n");
	}
	if (!(C_In_Range(*(dptr+zstep), low, high))) {
	    printf("Seed is on C_MEA_BACK.\n");
	}
*/
	if (!(C_On_Surface(dptr, low, high, xstep, ystep, zstep))) {
	    printf("Seed is NOT on surface!!!\n");
	    return(0);
	}
}

/****************************************/
/*					
/*	Procedure Name:  mea_tracking
/*	  Return Value:  none
/*     Input Variables:  int		method  - measurement method
/*			 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
/*			 int		seedx	- x coordinate of seed point
/*			 int		seedy	- y coordinate of seed point
/*			 int		seedz	- z coordinate of seed point
/*    Output Variables:  none
/*    Update Variables:  C_MeaVoxelList  *mea_surface_voxels[]
/*				 - output of the algorithm, a surface list
/*    Global Variables:  none
/*     	   Description:  This routine implements the surface tracking 
/*			 algorithm. See my report for the detail.
 */

void mea_tracking(method, mea_surface_voxels, data,xsize,ysize,zsize,low,high,seedx, seedy, seedz)
int	 method; 	
C_MeaVoxelList	*mea_surface_voxels[C_MEA_HASH_MAX];
C_Voxel_8bit  *data;
int	 xsize,ysize,zsize;
int      low, high;
int	 seedx, seedy, seedz;
{
	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	hash_index;
	int	direction;	/* direction of current face	*/
	int	direct[2];	/* direction of faces current face can reach */
	int	index;		/* index of current face    */
	int	save_index;
	int	ind[2];		/* index of faces current face can reach */ 
	int	flag;
	int	i;
	int	j;

	C_MeaFaceList	*mea_marked_faces[C_MEA_HASH_MAX]; /* list of marked faces */
	C_MeaFaceList	*mea_queue; /* queue that contains faces to be processed */
	C_MeaFaceList	*tempptr, *ptr, *ptr0;
	C_MeaVoxelList	*vtempptr;

	/* initialize pointers */
	mea_queue = ptr = ptr0 = NULL;

	for (i=0; i<C_MEA_HASH_MAX; i++) {
	    mea_surface_voxels[i] = NULL;
	    mea_marked_faces[i] = NULL;
	}
	
	/* set up the increaments */

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

	index = (seedz*zstep) + (seedy*ystep) + (seedx * xstep);
	dptr = data + index;

	if (!(C_In_Range(*(dptr-xstep), low, high))) {
	    direction = C_MEA_LEFT;
	}
	else if (!(C_In_Range(*(dptr+xstep), low, high))) {
	    direction = C_MEA_RIGHT;
	}
	else if (!(C_In_Range(*(dptr-ystep), low, high))) {
	    direction = C_MEA_DOWN;
	}
	else if (!(C_In_Range(*(dptr+ystep), low, high))) {
	    direction = C_MEA_UP;
	}
	else if (!(C_In_Range(*(dptr-zstep), low, high))) {
	    direction = C_MEA_FRONT;
	}
	else if (!(C_In_Range(*(dptr+zstep), low, high))) {
	    direction = C_MEA_BACK;
	}
	else {
	    C_error_message("Error. Seed point is not on surface. \n");
	}

	/* put first face into queue */
	tempptr = ( C_MeaFaceList * )
                    malloc( sizeof ( C_MeaFaceList ) );
	if (tempptr == NULL) {
	    printf("Could not malloc C_MeaFaceList!!!\n");
	    exit(C_ERROR);
	}
	tempptr->face.index = index;
	tempptr->face.direction = direction;
	tempptr->next_face = NULL;
	mea_queue = tempptr;
 
	/* put two copies of first face into marked faces list */
	hash_index = index % C_MEA_HASH_MAX;

	tempptr = ( C_MeaFaceList * )
                    malloc( sizeof ( C_MeaFaceList ) );
	if (tempptr == NULL) {
	    printf("Could not malloc C_MeaFaceList!!!\n");
	    exit(C_ERROR);
	}
	tempptr->face.index = index;
	tempptr->face.direction = direction;
	tempptr->next_face = NULL;
	mea_marked_faces[hash_index] = tempptr;

	tempptr = ( C_MeaFaceList * )
                    malloc( sizeof ( C_MeaFaceList ) );
	if (tempptr == NULL) {
	    printf("Could not malloc C_MeaFaceList!!!\n");
	    exit(C_ERROR);
	}
	tempptr->face.index = index;
	tempptr->face.direction = direction;
	tempptr->next_face = mea_marked_faces[hash_index];
	mea_marked_faces[hash_index] = tempptr;

	/* Surface tracking starts */
	while (mea_queue != NULL) {
	    ptr = mea_queue;
	    ptr0 = mea_queue;
	    /* remove a face f from queue */
	    while (ptr->next_face != NULL) {
		ptr0 = ptr;
		ptr = ptr->next_face;
	    }
	    if (ptr == ptr0) {
		mea_queue = NULL;
	    }
	    else {
	        ptr0->next_face = NULL;
	    }
	    index = ptr->face.index;
	    direction = ptr->face.direction;
	    free(ptr);	

	    dptr = data + index;

	    /* find the two faces f1 and f2 that f's out-edge can reach */
	    switch (direction) {
		case C_MEA_RIGHT:
				if (C_In_Range(*(dptr+xstep+zstep), 
					 low, high)) {
	    				direct[0] = C_MEA_FRONT;
					ind[0] = index + xstep + zstep; 
				}
				else if (C_In_Range(*(dptr+zstep), low, high)) {
	    				direct[0] = C_MEA_RIGHT;
					ind[0] = index + zstep; 
				}
				else {
	    				direct[0] = C_MEA_BACK;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr+ystep+xstep), 
					 low, high)) {
	    				direct[1] = C_MEA_DOWN;
					ind[1] = index + xstep + ystep; 
				}
				else if (C_In_Range(*(dptr+ystep), low, high)) {
	    				direct[1] = C_MEA_RIGHT;
					ind[1] = index + ystep; 
				}
				else {
	    				direct[1] = C_MEA_UP;
					ind[1] = index; 
				}
				break;
		case C_MEA_LEFT:
				if (C_In_Range(*(dptr-xstep-zstep), 
					 low, high)) {
	    				direct[0] = C_MEA_BACK;
					ind[0] = index - xstep - zstep; 
				}
				else if (C_In_Range(*(dptr-zstep), low, high)) {
	    				direct[0] = C_MEA_LEFT;
					ind[0] = index - zstep; 
				}
				else {
	    				direct[0] = C_MEA_FRONT;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr-ystep-xstep), 
					 low, high)) {
	    				direct[1] = C_MEA_UP;
					ind[1] = index - xstep - ystep; 
				}
				else if (C_In_Range(*(dptr-ystep), low, high)) {
	    				direct[1] = C_MEA_LEFT;
					ind[1] = index - ystep; 
				}
				else {
	    				direct[1] = C_MEA_DOWN;
					ind[1] = index; 
				}
				break;
		case C_MEA_UP:
				if (C_In_Range(*(dptr-xstep+ystep), 
					 low, high)) {
	    				direct[0] = C_MEA_RIGHT;
					ind[0] = index - xstep + ystep; 
				}
				else if (C_In_Range(*(dptr-xstep), low, high)) {
	    				direct[0] = C_MEA_UP;
					ind[0] = index - xstep; 
                                }
				else {
	    				direct[0] = C_MEA_LEFT;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr+ystep+zstep), 
					 low, high)) {
	    				direct[1] = C_MEA_FRONT;
					ind[1] = index + ystep + zstep; 
				}
				else if (C_In_Range(*(dptr+zstep), low, high)) {
	    				direct[1] = C_MEA_UP;
					ind[1] = index + zstep; 
				}
				else {
	    				direct[1] = C_MEA_BACK;
					ind[1] = index; 
				}
				
				break;
		case C_MEA_DOWN:
				if (C_In_Range(*(dptr+xstep-ystep), 
					 low, high)) {
	    				direct[0] = C_MEA_LEFT;
					ind[0] = index + xstep - ystep; 
				}
				else if (C_In_Range(*(dptr+xstep), low, high)) {
	    				direct[0] = C_MEA_DOWN;
					ind[0] = index + xstep; 
                                }
				else {
	    				direct[0] = C_MEA_RIGHT;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr-ystep-zstep), 
					 low, high)) {
	    				direct[1] = C_MEA_BACK;
					ind[1] = index - ystep - zstep; 
				}
				else if (C_In_Range(*(dptr-zstep), low, high)) {
	    				direct[1] = C_MEA_DOWN;
					ind[1] = index - zstep; 
				}
				else {
	    				direct[1] = C_MEA_FRONT;
					ind[1] = index; 
				}
				break;
		case C_MEA_BACK:
				if (C_In_Range(*(dptr-ystep+zstep), 
					 low, high)) {
	    				direct[0] = C_MEA_UP;
					ind[0] = index - ystep + zstep; 
				}
				else if (C_In_Range(*(dptr-ystep), low, high)) {
	    				direct[0] = C_MEA_BACK;
					ind[0] = index - ystep; 
				}
				else {
	    				direct[0] = C_MEA_DOWN;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr-xstep+zstep), 
					 low, high)) {
	    				direct[1] = C_MEA_RIGHT;
					ind[1] = index - xstep + zstep; 
				}
				else if (C_In_Range(*(dptr-xstep), low, high)) {
	    				direct[1] = C_MEA_BACK;
					ind[1] = index - xstep; 
                                }
				else {
	    				direct[1] = C_MEA_LEFT;
					ind[1] = index; 
				}
				break;
		case C_MEA_FRONT:
				if (C_In_Range(*(dptr+ystep-zstep), 
					 low, high)) {
	    				direct[0] = C_MEA_DOWN;
					ind[0] = index + ystep - zstep; 
				}
				else if (C_In_Range(*(dptr+ystep), low, high)) {
	    				direct[0] = C_MEA_FRONT;
					ind[0] = index + ystep; 
				}
				else {
	    				direct[0] = C_MEA_UP;
					ind[0] = index; 
				}

				if (C_In_Range(*(dptr+xstep-zstep), 
					 low, high)) {
	    				direct[1] = C_MEA_LEFT;
					ind[1] = index + xstep - zstep; 
				}
				else if (C_In_Range(*(dptr+xstep), low, high)) {
	    				direct[1] = C_MEA_FRONT;
					ind[1] = index + xstep; 
                                }
				else {
	    				direct[1] = C_MEA_RIGHT;
					ind[1] = index; 
				}
				break;
		default:
				printf("Error: invalid face direction!!! \n");
				break;
	    }

	    /* output the voxel f belongs to */
	    save_index = index;
	    if (method == C_VOXEL_METHOD) {
	        hash_index = index % C_MEA_HASH_MAX;
	        if (mea_surface_voxels[hash_index] == NULL) {
	            vtempptr = ( C_MeaVoxelList * )
                       	         malloc( sizeof ( C_MeaVoxelList ) );
		    if (vtempptr == NULL) {
	    		printf("Could not malloc C_MeaVoxelList!!!\n");
	    		exit(C_ERROR);
		    }
	            vtempptr->voxel.index = index;
	            vtempptr->voxel.value = *dptr;
	            vtempptr->next_voxel = NULL;
		    mea_surface_voxels[hash_index] = vtempptr;
	        }
	        else {
	            vtempptr = ( C_MeaVoxelList * )
                   	         malloc( sizeof ( C_MeaVoxelList ) );
		    if (vtempptr == NULL) {
	    		printf("Could not malloc C_MeaVoxelList!!!\n");
	    		exit(C_ERROR);
		    }
	            vtempptr->voxel.index = index;
	            vtempptr->voxel.value = *dptr;
	            vtempptr->next_voxel = mea_surface_voxels[hash_index]; 
		    mea_surface_voxels[hash_index] = vtempptr;
	        }
	    }
	    else {
		for (j=0; j<8; j++) {
                    switch (j) {
                        case 0: index = save_index;
	    			dptr = data + index;
                                break;
                        case 1: index = save_index - xstep;
	    			dptr = data + index;
                                break;
                        case 2: index = save_index - xstep - zstep;
	    			dptr = data + index;
                                break;
                        case 3: index = save_index - zstep;
	    			dptr = data + index;
                                break;
                        case 4: index = save_index - ystep;
	    			dptr = data + index;
                                break;
                        case 5: index = save_index - xstep - ystep;
	    			dptr = data + index;
                                break;
                        case 6: index = save_index - xstep - ystep - zstep;
	    			dptr = data + index;
                                break;
			case 7: index = save_index - ystep - zstep;
	    			dptr = data + index;
                                break;
                    }
	            hash_index = index % C_MEA_HASH_MAX;
	            if (mea_surface_voxels[hash_index] == NULL) {
	                vtempptr = ( C_MeaVoxelList * )
                       	             malloc( sizeof ( C_MeaVoxelList ) );
		        if (vtempptr == NULL) {
	    		    printf("Could not malloc C_MeaVoxelList!!!\n");
	    		    exit(C_ERROR);
		        }
	                vtempptr->voxel.index = index;
	                vtempptr->voxel.value = *dptr;
	                vtempptr->next_voxel = NULL;
		        mea_surface_voxels[hash_index] = vtempptr;
	            }
	            else {
	                vtempptr = ( C_MeaVoxelList * )
                   	             malloc( sizeof ( C_MeaVoxelList ) );
		        if (vtempptr == NULL) {
	    		    printf("Could not malloc C_MeaVoxelList!!!\n");
	    		    exit(C_ERROR);
		        }
	                vtempptr->voxel.index = index;
	                vtempptr->voxel.value = *dptr;
	                vtempptr->next_voxel = mea_surface_voxels[hash_index]; 
		        mea_surface_voxels[hash_index] = vtempptr;
		    }
	        }
	    }

	    /* if f1 or f2 is in the marked faces list, delete it. otherwise
		put it in the list */
	    for (i=0; i<2; i++) {
		hash_index = ind[i] % C_MEA_HASH_MAX;
		ptr = ptr0 = mea_marked_faces[hash_index];
		flag = 0;
		while ((ptr != NULL) && (flag == 0)) {
		    if ((ptr->face.index == ind[i]) && 
			(ptr->face.direction == direct[i])) {
			flag = 1;
			if (ptr == ptr0) {
			    mea_marked_faces[hash_index] = 
					mea_marked_faces[hash_index]->next_face;
			    free(ptr);
			}
			else {
			    ptr0->next_face = ptr->next_face;
			    free(ptr);
			}
		    }
		    else {
			ptr0 = ptr;
			ptr = ptr->next_face;
		    }
		}

		if (flag == 0) {
		    if (mea_marked_faces[hash_index] == NULL) {
		        tempptr = ( C_MeaFaceList * )
                           	    malloc( sizeof ( C_MeaFaceList ) );
			if (tempptr == NULL) {
	    		    printf("Could not malloc C_MeaFaceList!!!\n");
	    		    exit(C_ERROR);
			}
		        tempptr->face.index = ind[i];
		        tempptr->face.direction = direct[i];
		        tempptr->next_face = NULL;
		        mea_marked_faces[hash_index] = tempptr;
		    }
		    else {
		        tempptr = ( C_MeaFaceList * )
                           	    malloc( sizeof ( C_MeaFaceList ) );
			if (tempptr == NULL) {
	    		    printf("Could not malloc C_MeaFaceList!!!\n");
	    		    exit(C_ERROR);
			}
		        tempptr->face.index = ind[i];
		        tempptr->face.direction = direct[i];
		        tempptr->next_face = mea_marked_faces[hash_index];
		        mea_marked_faces[hash_index] = tempptr;
		    }
		    if (mea_queue == NULL) {
		        tempptr = ( C_MeaFaceList * )
                   	            malloc( sizeof ( C_MeaFaceList ) );
			if (tempptr == NULL) {
	    		    printf("Could not malloc C_MeaFaceList!!!\n");
	    		    exit(C_ERROR);
			}
		        tempptr->face.index = ind[i];
		        tempptr->face.direction = direct[i];
		        tempptr->next_face = NULL;
		        mea_queue = tempptr;
		    }
		    else {
		        tempptr = ( C_MeaFaceList * )
                   	            malloc( sizeof ( C_MeaFaceList ) );
			if (tempptr == NULL) {
	    		    printf("Could not malloc C_MeaFaceList!!!\n");
	    		    exit(C_ERROR);
			}
		        tempptr->face.index = ind[i];
		        tempptr->face.direction = direct[i];
		        tempptr->next_face = mea_queue;
		        mea_queue = tempptr;
		    }
		}
	    }
	}
}

/****************************************/
/*					
/*	Procedure Name:  mea_clean_list
/*	  Return Value:  none
/*     Input Variables:  C_Voxel_8bit	*data	- pointer to raw data array
/*    Output Variables:  none
/*    Update Variables:  C_MeaVoxelList  *mea_surface_voxels[]
/*				 - output of the algorithm, a surface list
/*    Global Variables:  none
/*     	   Description:  This routine deletes the repeated voxels in surface  
/*			 list. 
*/

void mea_clean_list(mea_surface_voxels, data)
C_MeaVoxelList	*mea_surface_voxels[C_MEA_HASH_MAX];
C_Voxel_8bit  *data;
{
	C_MeaVoxelList	*vptr, *vptr0, *tempptr;
	C_Voxel_8bit		*dptr;
	int		i;

	for (i=0; i<C_MEA_HASH_MAX; i++) {
	    vptr = mea_surface_voxels[i];
	    while (vptr != NULL) {
		dptr = data + vptr->voxel.index;
		if (*dptr == vptr->voxel.value) {
		    *dptr = 255;
		}
		else {
		    vptr->voxel.value = 255;
		}
		vptr = vptr->next_voxel;
	    }
	}

	for (i=0; i<C_MEA_HASH_MAX; i++) {
	    vptr = mea_surface_voxels[i];
	    vptr0 = mea_surface_voxels[i];
	    while (vptr != NULL) {
		if (vptr->voxel.value == 255) {
		    if (vptr == vptr0) {
			tempptr = vptr;
			mea_surface_voxels[i] = 
					mea_surface_voxels[i]->next_voxel;
			vptr0 = vptr0->next_voxel;
			vptr = vptr->next_voxel;
			free(tempptr);
		    }
		    else {
			tempptr = vptr;
			vptr0->next_voxel = vptr->next_voxel;
			vptr = vptr->next_voxel;
			free(tempptr);
		    }
		}
		else {
		    dptr = data + vptr->voxel.index;
		    *dptr = vptr->voxel.value;
	            vptr0 = vptr;
	            vptr = vptr->next_voxel;
		}
	    }
	}
}

