/****************************************************************************************
 *											*
 *	An Implementation of the Marching Cubes Algorithm with some variations		*
 *											*
 *			Program by Andreas Hadjiprocopis				*
 *			March/April 1995						*
 *											*
 ****************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <sys/fcntl.h>

#define	MAX_FILENAME		250
#define	MAX_CUBE_TRIANGLES	5
#define	TRUE			1
#define	FALSE			0
#define	ERROR			-1
#define	BIG_NUMBER		1000000
#define	WHITE			0
#define	BLACK			1
#define	DXF			0x2
#define	NFF			0x3
#define	IV			0x8
#define	ONE_BYTE		0x4
#define	TWO_BYTE		0x5
#define	CUBE_MOVES_XYZ		0x6
#define	CUBE_MOVES_YXZ		0x7

#define	DeduceVertexValue1(rv, ll, ul)	( (((rv)<=(ul))&&((rv)>=(ll))) ? (BLACK):(WHITE) )

typedef	struct	_point_structure {
	double		X, Y, Z;
} PointStructure;

typedef	struct	_cube_structure {
	PointStructure	Vertex[8], Edge[12],
			Triangle[MAX_CUBE_TRIANGLES][3];
	int		NumberOfTriangles,
			VertexValue[8];
} CubeStructure;

typedef	struct	_triangle_data_type {
	int		Vertex[3];
} TriangleType;

typedef struct	lookup_table {
	int		Vertex[8],
			NumberOfTriangles,
			Occupied;
	TriangleType	Triangle[MAX_CUBE_TRIANGLES];
} LookUpTableDataType;

LookUpTableDataType	LookUpTable[256], BasisLookUpTable[30];
int			ImageWidth, ImageHeight;
int			TotalEntries;
int			OutputFileFormat;
int			IsoRangeLow, IsoRangeHigh;
int			AveragingWindowWidth, AveragingWindowHeight, AveragingWindowThreshold;
int			CubeMovesOrder;
int			TotalTriangles;

/*
	      k = z
		|  / j = y
		| /
		|/
		------- i = x

		7		6
	BACK	x------e6-------x
		|		|
	    e7	|		|
	4      e11	5   e5	e10
FRONT	x--------e4-----x	|
	|	|	|	|
	|	x--------e2-----x
       e8	3	|	2
	|  e3	       e9  e1
	|		|
	x------e0-------x
	0		1
*/

int	DoMarchingCubes(int *current, int *next, int current_slice, double slice_separation, int cube_width, int cube_height);
int	AssignValuesToCube(int *current, int *next, int cube_width, int cube_height, double slice_separation, int i, int j, int k, CubeStructure *my_cube);
void	ConvertCharsToInts(unsigned char *char_data, int *int_data, int format, int size);
int	CalculateTriangles(CubeStructure *my_cube);
void	OutputTriangles(CubeStructure *my_cube);
void	OutputFileHeader(int num_slices, double slice_separation);
void	OutputFileTail(void);
int	DeduceVertexValue(int *data, int x_coord, int y_coord, int low_limit, int upper_limit);

void	CreateLookUpTable(void);
void	RotateAroundX_Axis(LookUpTableDataType *cube_in);
void	RotateAroundY_Axis(LookUpTableDataType *cube_in);
void	RotateAroundZ_Axis(LookUpTableDataType *cube_in);
void	CopyLookUpTables(LookUpTableDataType *table_target, LookUpTableDataType *table_source);
int	GiveCubeIndex(int *vertex_value);

void	main(int argc, char **argv)
{
	int		CurrentSlice, NextSlice, FileSize, i, j,
			Isovalue, InputFilesNumber, current_triangles[500],
			InputFileFormat, limit, CubeWidth, CubeHeight,
			*CurrentData, *NextData;
	double		SliceSeparation;
	char		**InputFilesList;
	unsigned char	*dummy;

	if( argc < 16 ){
		fprintf(stderr, "Usage: %s ImageWidth ImageHeight IsoRangeLow IsoRangeHigh SliceSeparation InputFileFormat[OneByte|TwoByte] OutputFileFormat[DXF|NFF|IV] CubeWidth CubeHeight AveragingWindowWidth AveragingWindowHeight AveragingWindowThreshold CubeMoveOrder[XYZ|YXZ] InputFilesList...\n", argv[0]);
		fprintf(stderr, "The input files list must contain at least two valid files.\n");
		fprintf(stderr, "The input file is binary and can have one byte or two bytes per pixel.\n");
		exit(-1);
	}
	if( (ImageWidth=atoi(argv[1])) <= 0 ){
		fprintf(stderr, "Image width must be a positive integer.\n");
		exit(-1);
	}
	if( (ImageHeight=atoi(argv[2])) <= 0 ){
		fprintf(stderr, "Image height must be a positive integer.\n");
		exit(-1);
	}
	IsoRangeLow = atoi(argv[3]);
	IsoRangeHigh= atoi(argv[4]);
	if( (IsoRangeLow>IsoRangeHigh) || (IsoRangeLow<0) || (IsoRangeHigh<0) ){
		fprintf(stderr, "IsoRangeLow/IsoRangeHigh must not be negative.\nNote order(Low<High).\n");
		exit(-1);
	}
	if( (SliceSeparation=atof(argv[5])) <= 0.0 ){
		fprintf(stderr, "The distance between adjacent slices (SliceSeparation) must be positive.\n");
		exit(-1);
	}

	if( toupper(argv[6][0]) == 'O' ) InputFileFormat = ONE_BYTE;
	else InputFileFormat = TWO_BYTE;
	if( toupper(argv[7][0]) == 'N' ) OutputFileFormat = NFF;
	else if( toupper(argv[7][0]) == 'D' ) OutputFileFormat = DXF;
	else OutputFileFormat = IV;

	if( InputFileFormat == TWO_BYTE )
		FileSize = ImageHeight * ImageWidth * 2;
	else
		FileSize = ImageHeight * ImageWidth;
	limit = (InputFileFormat == ONE_BYTE ? 256 : 65536);
	if( (IsoRangeLow>=limit) || (IsoRangeHigh>=limit) ){
		fprintf(stderr, "IsoRangeLow/IsoRangeHigh must not exceed image resolution.\n");
		exit(-1);
	}
	if( (CubeWidth=atoi(argv[8])) <= 0 ){
		fprintf(stderr, "Cube width must be positive.\n");
		exit(-1);
	}
	if( (CubeHeight=atoi(argv[9])) <= 0 ){
		fprintf(stderr, "Cube height must be positive.\n");
		exit(-1);
	}
	/* Negative or zero value on the averaging parameters indicates NO averaging */
	AveragingWindowWidth = atoi(argv[10]);
	AveragingWindowHeight= atoi(argv[11]);
	AveragingWindowThreshold = atoi(argv[12]);
	if( (AveragingWindowWidth<=0) ||
	    (AveragingWindowHeight<=0)||
	    (AveragingWindowThreshold<=0) )
			AveragingWindowWidth =
			AveragingWindowHeight=
			AveragingWindowThreshold = 0;
	if( AveragingWindowThreshold > ((AveragingWindowWidth*2+1)*(AveragingWindowHeight*2+1)) ){
		fprintf(stderr, "Averaging Window Threshold is too big.\n");
		exit(-1);
	}
	if( !strcmp("YXZ", argv[13]) )
		CubeMovesOrder = CUBE_MOVES_YXZ;
	else
		CubeMovesOrder = CUBE_MOVES_XYZ;

	InputFilesNumber = argc - 14; 
	if( (InputFilesList=(char **)calloc(InputFilesNumber+5, sizeof(char *))) == NULL ){
		fprintf(stderr, "Memory allocation error: failed for %d+5 items for InputFilesList as char **.\n", InputFilesNumber);
		exit(-1);
	}
	for(i=0;i<InputFilesNumber+5;i++)
		if( (InputFilesList[i]=(char *)calloc(MAX_FILENAME, sizeof(char))) == NULL ){
			fprintf(stderr, "Memory allocation error: failed for %d+5 items for InputFilesList[%d] as char *.\n", MAX_FILENAME, i);
			exit(-1);
		}
	for(i=14;i<argc;i++){
		strcpy(InputFilesList[i-14], argv[i]);
		if( (CurrentSlice=open(argv[i], O_RDONLY)) == -1 ){
			fprintf(stderr, "Could not open file '%s' for reading.\n", argv[i]);
			exit(-1);
		}
		close(CurrentSlice);
	}
	/* Check whether all input files are valid */

	if( (CurrentData=(int *)calloc(ImageHeight*ImageWidth+5, sizeof(int))) == NULL ){
		fprintf(stderr, "Memory allocation error: failed for %d+5 items for CurrentData as int *.\n", ImageHeight*ImageWidth);
		exit(-1);
	}
	if( (NextData=(int *)calloc(FileSize+5, sizeof(int))) == NULL ){
		fprintf(stderr, "Memory allocation error: failed for %d+5 items for NextData as int *.\n", FileSize);
		free(CurrentData);
		exit(-1);
	}
	if( (dummy=(unsigned char *)calloc(FileSize+5, sizeof(unsigned char))) == NULL ){
		fprintf(stderr, "Memory allocation error: failed for %d+5 items for dummy as unsigned char *.\n", FileSize);
		free(CurrentData);free(NextData);
		exit(-1);
	}

	if( (CurrentSlice=open(InputFilesList[0], O_RDONLY)) == -1 ){
		fprintf(stderr, "Could not open file '%s' for reading.\n", InputFilesList[0]);
		free(CurrentData);free(NextData);
		exit(-1);
	}
	fprintf(stderr, "Reading image file '%s'.\n", InputFilesList[0]);
	if( read(CurrentSlice, dummy, FileSize) != FileSize ){
		fprintf(stderr, "Error reading data from file '%s'.\n", InputFilesList[0]);
		free(CurrentData);free(NextData);free(dummy);
		close(CurrentSlice);
		exit(-1);
	}
	close(CurrentSlice);
	ConvertCharsToInts(dummy, CurrentData, InputFileFormat, FileSize);
	/* Create Look Up Table */
	CreateLookUpTable();

	OutputFileHeader(InputFilesNumber, SliceSeparation);

	for(i=1;i<InputFilesNumber;i++){
		if( (NextSlice=open(InputFilesList[i], O_RDONLY)) == -1 ){
			fprintf(stderr, "Could not open file '%s' for reading.\n", InputFilesList[i]);
			free(CurrentData);free(NextData);free(dummy);
			exit(-1);
		}
		fprintf(stderr, "Reading image file '%s'.\n", InputFilesList[i]);
		if( read(NextSlice, dummy, FileSize) != FileSize ){
			fprintf(stderr, "Error reading data from file '%s'.\n", InputFilesList[i]);
			free(CurrentData);free(NextData);free(dummy);
			close(NextSlice);
			exit(-1);
		}
		close(NextSlice);
		ConvertCharsToInts(dummy, NextData, InputFileFormat, FileSize);
		fprintf(stderr, "Processing image files '%s' and '%s' ", InputFilesList[i-1], InputFilesList[i]);
		fflush(stderr);
		current_triangles[i-1] = DoMarchingCubes(CurrentData, NextData, i-1, SliceSeparation, CubeWidth, CubeHeight);
		fprintf(stderr, ", %d triangles created.\n", current_triangles[i-1]);
		if( i != (InputFilesNumber-1) )
			for(j=0;j<FileSize;j++)
				CurrentData[j] = NextData[j];
  	}
	TotalTriangles = 0;
	for(i=0;i<InputFilesNumber-1;i++){
		fprintf(stderr, "Slice %d - Slice %d: # of triangles is %d\n", i+1, i+2, current_triangles[i]);
		TotalTriangles += current_triangles[i];
	}
	OutputFileTail();
	fprintf(stderr, "Total # of triangles is %d\nTotal # of slices is %d.\n\n", TotalTriangles, InputFilesNumber);
	if( AveragingWindowWidth <= 0 )
		fprintf(stderr, "Averaging Window Method was not used.\n");
	else
		fprintf(stderr, "Averaging Window Method was used,\nWindow Width:\t%d\nWindow Height:\t%d\nThreshold:\t%d out of %d.\n", AveragingWindowWidth*2+1, AveragingWindowHeight*2+1, AveragingWindowThreshold, (AveragingWindowHeight*2+1)*(AveragingWindowWidth*2+1));
}


int	DoMarchingCubes(int *current, int *next, int current_slice, double slice_separation, int cube_width, int cube_height)
{
	int		i, j, index, triangles_number = 0;
	CubeStructure	my_cube;

	if( CubeMovesOrder == CUBE_MOVES_XYZ ){
		for(j=0;j<ImageHeight-cube_height;j+=cube_height)
			for(i=0;i<ImageWidth-cube_width;i+=cube_width){
				AssignValuesToCube(current, next, cube_width, cube_height, slice_separation, i, j, current_slice, &my_cube);
				index = GiveCubeIndex(my_cube.VertexValue);
				if( (index==0) || (index==255) ) continue;
				CalculateTriangles(&my_cube);
				OutputTriangles(&my_cube);
				triangles_number += my_cube.NumberOfTriangles;
			}
	} else {
		for(i=0;i<ImageWidth-cube_width;i+=cube_width)
			for(j=0;j<ImageHeight-cube_height;j+=cube_height){	
				AssignValuesToCube(current, next, cube_width, cube_height, slice_separation, i, j, current_slice, &my_cube);
				index = GiveCubeIndex(my_cube.VertexValue);
				if( (index==0) || (index==255) ) continue;
				CalculateTriangles(&my_cube);
				OutputTriangles(&my_cube);
				triangles_number += my_cube.NumberOfTriangles;
			}
	}

	return(triangles_number);
}

void	ConvertCharsToInts(unsigned char *char_data, int *int_data, int format, int size)
{
	int		i, j;

	if( format == ONE_BYTE ){
		for(i=0;i<size;i++)
			int_data[i] = (int )char_data[i];
		return;
	}

	j = 0;
	if( format == TWO_BYTE ){
		for(i=1;i<size-1;i+=2)
			int_data[j++] = 256 * ((int )char_data[i]) + (int )char_data[i+1];
		return;
	}
}



int	AssignValuesToCube(int *current, int *next, int cube_width, int cube_height, double slice_separation, int i, int j, int k, CubeStructure *my_cube)
{
	float cw, ch;

	/* If no averaging is desired */
	if( AveragingWindowWidth <= 0 ){
		my_cube->VertexValue[0] = DeduceVertexValue1(current[j*ImageWidth+i], IsoRangeLow, IsoRangeHigh); 
		my_cube->VertexValue[1] = DeduceVertexValue1(current[j*ImageWidth+i+cube_width], IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[2] = DeduceVertexValue1(current[(j+cube_height)*ImageWidth+i+cube_width], IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[3] = DeduceVertexValue1(current[(j+cube_height)*ImageWidth+i], IsoRangeLow, IsoRangeHigh); 

		my_cube->VertexValue[4] = DeduceVertexValue1(next[j*ImageWidth+i], IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[5] = DeduceVertexValue1(next[j*ImageWidth+i+cube_width], IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[6] = DeduceVertexValue1(next[(j+cube_height)*ImageWidth+i+cube_width], IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[7] = DeduceVertexValue1(next[(j+cube_height)*ImageWidth+i], IsoRangeLow, IsoRangeHigh );
	} else {
		my_cube->VertexValue[0] = DeduceVertexValue(current, i, j, IsoRangeLow, IsoRangeHigh); 
		my_cube->VertexValue[1] = DeduceVertexValue(current, i+cube_width, j, IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[2] = DeduceVertexValue(current, i+cube_width, j+cube_height, IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[3] = DeduceVertexValue(current, i, j+cube_height, IsoRangeLow, IsoRangeHigh); 

		my_cube->VertexValue[4] = DeduceVertexValue(next, i, j, IsoRangeLow, IsoRangeHigh); 
		my_cube->VertexValue[5] = DeduceVertexValue(next, i+cube_width, j, IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[6] = DeduceVertexValue(next, i+cube_width, j+cube_height, IsoRangeLow, IsoRangeHigh);
		my_cube->VertexValue[7] = DeduceVertexValue(next, i, j+cube_height, IsoRangeLow, IsoRangeHigh); 
	}

	cw = (float )cube_width; ch = (float )cube_height;
	my_cube->Edge[0].X = i+cw/2;my_cube->Edge[0].Y = j;my_cube->Edge[0].Z = k*slice_separation;
	my_cube->Edge[1].X = i+cw;my_cube->Edge[1].Y = j+ch/2;my_cube->Edge[1].Z = k*slice_separation;
	my_cube->Edge[2].X = i+cw/2;my_cube->Edge[2].Y = j+ch;my_cube->Edge[2].Z = k*slice_separation;
	my_cube->Edge[3].X = i;my_cube->Edge[3].Y = j+ch/2;my_cube->Edge[3].Z = k*slice_separation;

	my_cube->Edge[4].X = i+cw/2;my_cube->Edge[4].Y = j;my_cube->Edge[4].Z = (k+1)*slice_separation;
	my_cube->Edge[5].X = i+cw;my_cube->Edge[5].Y = j+ch/2;my_cube->Edge[5].Z = (k+1)*slice_separation;
	my_cube->Edge[6].X = i+cw/2;my_cube->Edge[6].Y = j+ch;my_cube->Edge[6].Z = (k+1)*slice_separation;
	my_cube->Edge[7].X = i;my_cube->Edge[7].Y = j+ch/2;my_cube->Edge[7].Z = (k+1)*slice_separation;

	my_cube->Edge[8].X = i;my_cube->Edge[8].Y = j;my_cube->Edge[8].Z = (k+0.5)*slice_separation;
	my_cube->Edge[9].X = i+cw;my_cube->Edge[9].Y = j;my_cube->Edge[9].Z = (k+0.5)*slice_separation;
	my_cube->Edge[10].X = i+cw;my_cube->Edge[10].Y = j+ch;my_cube->Edge[10].Z = (k+0.5)*slice_separation;
	my_cube->Edge[11].X = i;my_cube->Edge[11].Y = j+ch;my_cube->Edge[11].Z = (k+0.5)*slice_separation;

	my_cube->NumberOfTriangles = 0;

	return(TRUE);
}

int	CalculateTriangles(CubeStructure *my_cube)
{
	int		index, i, j;
	int		vertex_value[8];

	index = GiveCubeIndex(my_cube->VertexValue);

	my_cube->NumberOfTriangles = LookUpTable[index].NumberOfTriangles;

	for(i=0;i<LookUpTable[index].NumberOfTriangles;i++)
		for(j=0;j<3;j++){
			my_cube->Triangle[i][j].X = my_cube->Edge[LookUpTable[index].Triangle[i].Vertex[j]].X;
			my_cube->Triangle[i][j].Y = my_cube->Edge[LookUpTable[index].Triangle[i].Vertex[j]].Y;
			my_cube->Triangle[i][j].Z = my_cube->Edge[LookUpTable[index].Triangle[i].Vertex[j]].Z;
		}
	return(TRUE);
}


void	OutputTriangles(CubeStructure *my_cube)
{
	int		i, j;

	if( OutputFileFormat == NFF ){
		for(i=0;i<my_cube->NumberOfTriangles;i++){
			printf("f Red 1 0 0 0 0\n");
			printf("p 3\n");
			for(j=0;j<3;j++)
				printf("%2.1lf %2.1lf %2.1lf\n", my_cube->Triangle[i][j].X, my_cube->Triangle[i][j].Y, my_cube->Triangle[i][j].Z);
		}
	}

	if( OutputFileFormat == DXF ){
		for(i=0;i<my_cube->NumberOfTriangles;i++){
			printf("  0\n3DFACE\n  8\n0\n  5\n1\n");
			for(j=0;j<3;j++)
				printf(" %d\n%2.1lf\n %d\n%2.1lf\n %d\n%2.1lf\n",
					 j+10,
					 my_cube->Triangle[i][j].X,
					 j+20,
					 my_cube->Triangle[i][j].Y,
					 j+30,
					 my_cube->Triangle[i][j].Z );
			printf(" 13\n%2.1lf\n 23\n%2.1lf\n 33\n%2.1lf\n",
					 my_cube->Triangle[i][2].X,
					 my_cube->Triangle[i][2].Y,
					 my_cube->Triangle[i][2].Z);
		}
	}

	if( OutputFileFormat == IV ){
		for(i=0;i<my_cube->NumberOfTriangles;i++){
			printf("\t\t\t\t\t%2.1lf %2.1lf %2.1lf ,\n\t\t\t\t\t%2.1lf %2.1lf %2.1lf ,\n\t\t\t\t\t%2.1lf %2.1lf %2.1lf ,\n",
								my_cube->Triangle[i][0].X,
								my_cube->Triangle[i][0].Y,
								my_cube->Triangle[i][0].Z,
								my_cube->Triangle[i][1].X,
								my_cube->Triangle[i][1].Y,
								my_cube->Triangle[i][1].Z,
								my_cube->Triangle[i][2].X,
								my_cube->Triangle[i][2].Y,
								my_cube->Triangle[i][2].Z );
		}
	}
}

void	OutputFileHeader(int num_slices, double slice_separation)
{

	if( OutputFileFormat == NFF ){
		printf("v\nfrom %lf %lf %lf\nat %lf %lf 0.0\nup 0 0 1\nangle 70\nhither 1\nresolution %d %d\n",
			ImageWidth*1.35, ImageHeight*0.5, num_slices*slice_separation*7.0,
			ImageWidth*0.75, ImageHeight*0.5,
			ImageWidth, ImageHeight);
		printf("b White\nl 0 %lf %lf\nl %lf 0 %lf\nl %lf %lf %lf\nl %lf %lf %lf\n",
			  ImageHeight*1.4, num_slices*slice_separation*1.4,
			  ImageWidth*1.4, num_slices*1.4,
			  ImageWidth*1.4, ImageHeight*1.4, -num_slices*1.4,
			  ImageWidth*1.4, ImageHeight*1.4,  num_slices*1.4);
	}

	if( OutputFileFormat == DXF ){
		printf("  0\nSECTION\n  2\nENTITIES\n");
	}

	if( OutputFileFormat == IV ){
		printf("#Inventor V1.0 ascii\nSeparator {\n\tRotationXYZ {\n\t\taxis\tX\n\t\tangle\t-1.5708\n\t}\n\tSeparator {\n\t\tSeparator {\n\t\t\tCoordinate3 {\n\t\t\t\tpoint[\n");
	}
}

void	OutputFileTail(void)
{
	int	i, index;

	if( OutputFileFormat == NFF ){
		return;
	}

	if( OutputFileFormat == DXF ){
		printf("  0\nENDSEC\n  0\nEOF\n");
		fflush(stdout);
		return;
	}

	if( OutputFileFormat == IV ){
		printf(" ]\n\t\t\t}\n");
		printf("\t\t\tShapeHints {\n\t\t\t\thints\t(ORDERED | CONVEX)\n\t\t\t}\n");
		printf("\t\t\tNormal {\n\t\t\t\tvector [-0.707107 -0.707107 0]\n\t\t\t}\n");
		printf("\t\t\tNormalBinding {\n\t\t\t\tvalue PER_VERTEX_INDEXED\n\t\t\t}\n");
		printf("\t\t\tIndexedFaceSet {\n\t\t\t\tcoordIndex [\n");
		index = 0;
		for(i=0;i<TotalTriangles/2-1;i++){
			printf("\t\t\t\t\t%d , %d , %d , -1 ,\n\t\t\t\t\t%d , %d , %d , -1 ,\n",
				index, index+1, index+2, index+3, index+4, index+5 );
			index += 6;
		}
		printf("\t\t\t\t\t%d , %d , %d , -1,\n\t\t\t\t\t%d , %d , %d , -1 ]\n",
                                index, index+1, index+2, index+3, index+4, index+5 );
		printf("\t\t\t}\n\t\t}\n\t}\n}\n");
	}
}


int	DeduceVertexValue(int *data, int x_coord, int y_coord, int low_limit, int upper_limit)
{
	int		sum, i, j;

	/* if ther edge is white there is no need for averaging */
	if( DeduceVertexValue1(data[x_coord+y_coord*ImageWidth], low_limit, upper_limit) == WHITE )
		return( WHITE );

	/* Only need to average if the vertex value is  BLACK */
	/* If we are at the edges */
	if( (x_coord<AveragingWindowWidth) || (x_coord>=(ImageWidth-AveragingWindowWidth)) ||
	    (y_coord<AveragingWindowHeight) || (y_coord>=(ImageHeight-AveragingWindowHeight)) )
		return( DeduceVertexValue1(data[x_coord+y_coord*ImageWidth], low_limit, upper_limit) );

	sum = 0;
	for(j=y_coord-AveragingWindowHeight;j<=y_coord+AveragingWindowHeight;j++)
		for(i=x_coord-AveragingWindowWidth;i<=x_coord+AveragingWindowWidth;i++)
			sum += DeduceVertexValue1(data[i+j*ImageWidth], low_limit, upper_limit);

	if( sum >= AveragingWindowThreshold ) return( BLACK );

	return( WHITE );
}


void	CreateLookUpTable(void)
{
	int			i, j, k, basic;
	LookUpTableDataType	temporary_cube;

	/* Fill up the LookUpTable with the 15 cube arrangments as described by
	   the SIGGRAPH 87 paper(Computer Graphics(ACM), Volume 21, Number 4, July 1987.
	   The N0-BLACK/NO-WHITE case is case 0 */

	TotalEntries = 0;
	/* Put all white first */
	for(basic=0;basic<15;basic++)
		for(i=0;i<8;i++)
			BasisLookUpTable[basic].Vertex[i] = WHITE;

	for(i=0;i<256;i++)
		LookUpTable[i].Occupied = FALSE;

basic = 0;
	BasisLookUpTable[basic].NumberOfTriangles = 0;

basic = 1;
	BasisLookUpTable[basic].Vertex[0] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 1;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 8;

basic = 2;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[1] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 2;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 1;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 3;

basic = 3;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[5] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 2;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 8;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 4;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 5;
	
basic = 4;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[6] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 2;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 8;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 10;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 5;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 6;

basic = 5;
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 9;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 3;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 11;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 3;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 9;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 10;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 11;
	
basic = 6;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[6] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 1;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 3;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 10;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 5;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 6;

basic = 7;
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[4] =
	BasisLookUpTable[basic].Vertex[6] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 7;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 6;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 5;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 10;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 1;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 9;

basic = 8;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 2;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 11;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 10;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 10;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 8;

basic = 9;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] =
	BasisLookUpTable[basic].Vertex[7] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 7;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 6;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 8;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 0;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 6;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 6;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 0;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 10;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 1;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 10;

basic = 10;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[4] =
	BasisLookUpTable[basic].Vertex[6] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 7;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 4;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 3;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 0;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 4;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 1;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 2;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 5;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 5;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 6;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 2;

basic = 11;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] =
	BasisLookUpTable[basic].Vertex[6] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 8;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 11;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 11;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 6;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 5;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 5;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 1;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 0;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 5;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 0;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 11;

basic = 12;
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] =
	BasisLookUpTable[basic].Vertex[4] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 9;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 11;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 3;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 11;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 9;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 10;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 7;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 4;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 8;

basic = 13;
	BasisLookUpTable[basic].Vertex[0] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[5] =
	BasisLookUpTable[basic].Vertex[7] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 8;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 1;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 2;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 10;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 5;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 9;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 4;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 7;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 6;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 11;

basic = 14;
	BasisLookUpTable[basic].Vertex[1] =
	BasisLookUpTable[basic].Vertex[2] =
	BasisLookUpTable[basic].Vertex[3] =
	BasisLookUpTable[basic].Vertex[7] = BLACK;
	BasisLookUpTable[basic].NumberOfTriangles = 4;
	BasisLookUpTable[basic].Triangle[0].Vertex[0] = 0;
	BasisLookUpTable[basic].Triangle[0].Vertex[1] = 3;
	BasisLookUpTable[basic].Triangle[0].Vertex[2] = 7;
	BasisLookUpTable[basic].Triangle[1].Vertex[0] = 7;
	BasisLookUpTable[basic].Triangle[1].Vertex[1] = 10;
	BasisLookUpTable[basic].Triangle[1].Vertex[2] = 0;
	BasisLookUpTable[basic].Triangle[2].Vertex[0] = 10;
	BasisLookUpTable[basic].Triangle[2].Vertex[1] = 1;
	BasisLookUpTable[basic].Triangle[2].Vertex[2] = 0;
	BasisLookUpTable[basic].Triangle[3].Vertex[0] = 7;
	BasisLookUpTable[basic].Triangle[3].Vertex[1] = 6;
	BasisLookUpTable[basic].Triangle[3].Vertex[2] = 10;


	/* Now find 15 more entries by complementing the 0 to 14 entries */
	/* Triangles stay as they are */
	for(basic=0;basic<15;basic++){
		CopyLookUpTables(&(BasisLookUpTable[basic+15]), &(BasisLookUpTable[basic]));
		for(i=0;i<8;i++)
			if( BasisLookUpTable[basic].Vertex[i] == BLACK )
				BasisLookUpTable[basic+15].Vertex[i] = WHITE;
			else
				BasisLookUpTable[basic+15].Vertex[i] = BLACK;
	}

	TotalEntries = 30;
	/* initialise occupied entries of look up table */
	for(i=0;i<256;i++)
		LookUpTable[i].Occupied = FALSE;

	/* Add all these entries to the normal look up table, the index is given by the
           the expression Sum( 2^(Black_Vertex_Number) ) */
	for(i=0;i<TotalEntries;i++)
		CopyLookUpTables( &(LookUpTable[GiveCubeIndex(BasisLookUpTable[i].Vertex)]), &(BasisLookUpTable[i]));

	/* And go up to 256 by rotating the first 30 entries. */
	/* Except, of course, from the first and 15th(NO-BLACK/NO-WHITE) cases */

	for(basic=1;basic<TotalEntries;basic++){
		if( basic == 15 ) continue;
		if( TotalEntries >= 256 ) break;
		CopyLookUpTables(&temporary_cube, &(BasisLookUpTable[basic]));
		for(i=0;i<4;i++){
			for(j=0;j<4;j++){
				RotateAroundZ_Axis(&temporary_cube);
				if( !LookUpTable[GiveCubeIndex(temporary_cube.Vertex)].Occupied ){
					CopyLookUpTables( &(LookUpTable[GiveCubeIndex(temporary_cube.Vertex)]), &temporary_cube );
					TotalEntries++;
				}
			}
			RotateAroundX_Axis(&temporary_cube);
		}
		RotateAroundY_Axis(&temporary_cube);
		for(j=0;j<4;j++){
			RotateAroundZ_Axis(&temporary_cube);
			if( !LookUpTable[GiveCubeIndex(temporary_cube.Vertex)].Occupied ){
				CopyLookUpTables( &(LookUpTable[GiveCubeIndex(temporary_cube.Vertex)]), &temporary_cube );
				TotalEntries++;
			}
		}
		RotateAroundY_Axis(&temporary_cube);
		RotateAroundY_Axis(&temporary_cube);
		for(j=0;j<4;j++){
			RotateAroundZ_Axis(&temporary_cube);
			if( !LookUpTable[GiveCubeIndex(temporary_cube.Vertex)].Occupied ){
				CopyLookUpTables( &(LookUpTable[GiveCubeIndex(temporary_cube.Vertex)]), &temporary_cube );
				TotalEntries++;
			}
		}			
	}

}

void	RotateAroundX_Axis(LookUpTableDataType *cube_in)
{
	int			i, j;
	LookUpTableDataType	temp_cube;

	temp_cube.Vertex[4] = cube_in->Vertex[0];
	temp_cube.Vertex[5] = cube_in->Vertex[1];
	temp_cube.Vertex[1] = cube_in->Vertex[2];
	temp_cube.Vertex[0] = cube_in->Vertex[3];
	temp_cube.Vertex[7] = cube_in->Vertex[4];
	temp_cube.Vertex[6] = cube_in->Vertex[5];
	temp_cube.Vertex[2] = cube_in->Vertex[6];
	temp_cube.Vertex[3] = cube_in->Vertex[7];

	temp_cube.NumberOfTriangles = cube_in->NumberOfTriangles;

	for(i=0;i<cube_in->NumberOfTriangles;i++)
		for(j=0;j<3;j++)
			switch( cube_in->Triangle[i].Vertex[j] ){
				case 0: temp_cube.Triangle[i].Vertex[j] = 4;break;
				case 1: temp_cube.Triangle[i].Vertex[j] = 9;break;
				case 2: temp_cube.Triangle[i].Vertex[j] = 0;break;
				case 3: temp_cube.Triangle[i].Vertex[j] = 8;break;
				case 4: temp_cube.Triangle[i].Vertex[j] = 6;break;
				case 5: temp_cube.Triangle[i].Vertex[j] = 10;break;
				case 6: temp_cube.Triangle[i].Vertex[j] = 2;break;
				case 7: temp_cube.Triangle[i].Vertex[j] = 11;break;
				case 8: temp_cube.Triangle[i].Vertex[j] = 7;break;
				case 9: temp_cube.Triangle[i].Vertex[j] = 5;break;
				case 10: temp_cube.Triangle[i].Vertex[j] = 1;break;
				case 11: temp_cube.Triangle[i].Vertex[j] = 3;break;
				default: printf("Nonsense X.\n");break;
			}
	
	CopyLookUpTables(cube_in, &temp_cube);
}


void	RotateAroundY_Axis(LookUpTableDataType *cube_in)
{
	int			i, j;
	LookUpTableDataType	temp_cube;

	temp_cube.Vertex[4] = cube_in->Vertex[0];
	temp_cube.Vertex[0] = cube_in->Vertex[1];
	temp_cube.Vertex[3] = cube_in->Vertex[2];
	temp_cube.Vertex[7] = cube_in->Vertex[3];
	temp_cube.Vertex[5] = cube_in->Vertex[4];
	temp_cube.Vertex[1] = cube_in->Vertex[5];
	temp_cube.Vertex[2] = cube_in->Vertex[6];
	temp_cube.Vertex[6] = cube_in->Vertex[7];

	temp_cube.NumberOfTriangles = cube_in->NumberOfTriangles;

	for(i=0;i<cube_in->NumberOfTriangles;i++)
		for(j=0;j<3;j++)
			switch( cube_in->Triangle[i].Vertex[j] ){
				case 0: temp_cube.Triangle[i].Vertex[j] = 8;break;
				case 1: temp_cube.Triangle[i].Vertex[j] = 3;break;
				case 2: temp_cube.Triangle[i].Vertex[j] = 11;break;
				case 3: temp_cube.Triangle[i].Vertex[j] = 7;break;
				case 4: temp_cube.Triangle[i].Vertex[j] = 9;break;
				case 5: temp_cube.Triangle[i].Vertex[j] = 1;break;
				case 6: temp_cube.Triangle[i].Vertex[j] = 10;break;
				case 7: temp_cube.Triangle[i].Vertex[j] = 5;break;
				case 8: temp_cube.Triangle[i].Vertex[j] = 4;break;
				case 9: temp_cube.Triangle[i].Vertex[j] = 0;break;
				case 10: temp_cube.Triangle[i].Vertex[j] = 2;break;
				case 11: temp_cube.Triangle[i].Vertex[j] = 6;break;
				default: printf("Nonsense Y.\n");break;
			}

	CopyLookUpTables(cube_in, &temp_cube);
}

void	RotateAroundZ_Axis(LookUpTableDataType *cube_in)
{
	int			i, j;
	LookUpTableDataType	temp_cube;

	temp_cube.Vertex[3] = cube_in->Vertex[0];
	temp_cube.Vertex[0] = cube_in->Vertex[1];
	temp_cube.Vertex[1] = cube_in->Vertex[2];
	temp_cube.Vertex[2] = cube_in->Vertex[3];
	temp_cube.Vertex[7] = cube_in->Vertex[4];
	temp_cube.Vertex[4] = cube_in->Vertex[5];
	temp_cube.Vertex[5] = cube_in->Vertex[6];
	temp_cube.Vertex[6] = cube_in->Vertex[7];

	temp_cube.NumberOfTriangles = cube_in->NumberOfTriangles;

	for(i=0;i<cube_in->NumberOfTriangles;i++)
		for(j=0;j<3;j++)
			switch( cube_in->Triangle[i].Vertex[j] ){
				case 0: temp_cube.Triangle[i].Vertex[j] = 3;break;
				case 1: temp_cube.Triangle[i].Vertex[j] = 0;break;
				case 2: temp_cube.Triangle[i].Vertex[j] = 1;break;
				case 3: temp_cube.Triangle[i].Vertex[j] = 2;break;
				case 4: temp_cube.Triangle[i].Vertex[j] = 7;break;
				case 5: temp_cube.Triangle[i].Vertex[j] = 4;break;
				case 6: temp_cube.Triangle[i].Vertex[j] = 5;break;
				case 7: temp_cube.Triangle[i].Vertex[j] = 6;break;
				case 8: temp_cube.Triangle[i].Vertex[j] = 11;break;
				case 9: temp_cube.Triangle[i].Vertex[j] = 8;break;
				case 10: temp_cube.Triangle[i].Vertex[j] = 9;break;
				case 11: temp_cube.Triangle[i].Vertex[j] = 10;break;
				default: printf("Nonsense Z.\n");break;
			}

	CopyLookUpTables(cube_in, &temp_cube);
}


void	CopyLookUpTables(LookUpTableDataType *table_target, LookUpTableDataType *table_source)
{
	int		i, j;

	for(i=0;i<8;i++)
		table_target->Vertex[i] = table_source->Vertex[i];

	table_target->NumberOfTriangles = table_source->NumberOfTriangles;
	for(i=0;i<table_source->NumberOfTriangles;i++)
		for(j=0;j<3;j++)
			table_target->Triangle[i].Vertex[j] = table_source->Triangle[i].Vertex[j];

	table_target->Occupied = TRUE;
}

int	GiveCubeIndex(int *vertex_value)
{
	int		i, pow = 1, sum = 0;

	for(i=0;i<8;i++){
		if( vertex_value[i] == BLACK )
			sum += pow;
		pow *= 2;
	}
	return(sum);
}