#include "STL_Distance.h"





/* Function Prototypes */

double SolveCircleCenter(int i);
void gregWireSphere(double r, double center[3]);
double norm(double v[3]);
bool BubbleSpheres(void);
int CanBubble(double r, GSphere* i, GSphere* j);//, int eiwl, int* Working_List);

void PrintWorkingList(int eiwl, int* Working_List);
void GenLoneSphere(GSphere* i, GSphere* new_sphere);
void GenMergeSphere(GSphere* i, GSphere* j, GSphere* new_sphere);

void PrintSphereData(GSphere x);
void PrintTestSphere(GSphere x);

void FindClosestSphere(GPoint TestPoint);


// math functions
void cross(double v1[3], double v2[3], double v3[3]);
void CalcRotationMatrix(double N[3], double R[3][3]);
void RotateVector(double R[3][3], double V[3]);
void CalcInverseRotation(double R[3][3]);
double abs_d(double x);
//int CheckIsDegenerate(double* v1,double* v2,double* v3);
void VectorProj(double* a, double* b, double*);
void VectBetweenPoints(double* a, double* b, double* c);
double ClosestTrianglePoint(int i,double* Point, double* closest);
void PointPlusVector(double* P, double* V, double* P2);
int f(double* V, double* C);
int GetMaxIndex(double* a);
void MakeEqual(double* a, double* b);
double CheckRemainingTriangles(double* Point, double* closest);
void RandomPointCollider();


void ConstructSphere(GSphere *sphere, GTriangle &triangle);

/**
 ** non-constant global variables:
 **/
static GLfloat spin = 0.0;
double Max[3];
double Min[3];
int num_triangles=0;
int num_sphere_levels=0;
GTriangle* Triangle_List;

GSphere* Sphere_List[MAX_SPHERE_LEVELS];
int sphere_level=0;
int num_spheres[MAX_SPHERE_LEVELS];

int numTestSpheres=0;

int numDegenerates=0;


int level_to_render = 0;

int getSphereRenderLevel()
{
	return level_to_render;
}

void changeSphereRenderLevel() {
	level_to_render ++;
	if( level_to_render >= num_sphere_levels)
		level_to_render = 0;
}

void MinMax(double x,double y, double z)
{
	if(Max[0]<x)
		Max[0]=x;
	if(Max[1]<y)
		Max[1]=y;
	if(Max[2]<z)
		Max[2]=z;

	if(Min[0]>x)
		Min[0]=x;
	if(Min[1]>y)
		Min[1]=y;
	if(Min[2]>z)
		Min[2]=z;
}

void cross(double v1[3], double v2[3], double v3[3])
{
	v3[0]=(v1[1]*v2[2])-(v1[2]*v2[1]);
	v3[1]=-((v1[0]*v2[2])-(v1[2]*v2[0]));
	v3[2]=(v1[0]*v2[1])-(v1[1]*v2[0]);
}

double dot(double v1[3], double v2[3])
{
	double val=0.0;
	val+=v1[0]*v2[0];
	val+=v1[1]*v2[1];
	val+=v1[2]*v2[2];
	return(val);
}

double norm(double v[3])
{
	return(sqrt(dot(v,v)));
}

void Normalize(double v[3])
{	
	double d=norm(v);
	if(d){
		for(int i=0;i<3;i++){
			v[i]/=d;
		}
	}
}


void LoadTriangles(CLEDSGeometry * STL);

void LoadSTL( char* filename , CLEDSGeometry * STL )
{

	FILE* file;

	if(file=fopen(filename,"r")){
		printf("File %s Opened OK!\n",filename);
	}
	else {
		printf("Failure Opening %s!\n",filename);
		exit(0);
	}
	
	char string1[15];
	
	//Read through and figure out how many triangles there are
	num_triangles=0;
	while(true){
		fscanf(file,"%s",string1);  //Just grab strings until you get something interesting
		if(!strncmp("endsolid",string1,8))
			break;
		if(!strncmp("loop",string1,4))
		{
			num_triangles++;
		}
	}
	
	num_spheres[0] = num_triangles;  //One sphere per triangle at first level
	Triangle_List  = new GTriangle[num_triangles];
	Sphere_List[0] = new GSphere[num_triangles];

	printf("STL file read.  Found %d triangles in file.\n",num_triangles);

	//Close the file and reopen it.
	fclose(file);

	LoadTriangles( STL );

}



void LoadTriangles(CLEDSGeometry * STL)
{

	int curr_triangle = 0;

	CListIter<CLEDSliteFace *, UINT_8> CListIterIntFace;
	CLEDSliteFace *pFace;
	CLEDSliteFaceIter FaceIter;
	CLEDSliteVertex *p1stVtx, *pVtx;
	Point P,Q,R,S;
	
	//Initialize face list iterator
	CListIterIntFace.Init(STL->GetCListFaces());

	//Pull out the first face
	pFace = CListIterIntFace.PeekFirst();

	while ( pFace ) 
	{ // While the face is valid
		// Initialize face iterator
		FaceIter.Init(pFace);

		GTriangle *tri = &Triangle_List[curr_triangle];
		GSphere *sphere = &Sphere_List[0][curr_triangle];

		tri->Normal = pFace->GetNormal();
		
		// perhaps later it will come in use
		// convertPointToDouble( pFace->GetCentroid() , sphere->Orig);

		// Get first vertex in face
		pVtx = p1stVtx = (CLEDSVertex *)FaceIter.PeekFirstCVertex();
		tri->Point1 = pVtx->GetPoint();

		pVtx = (CLEDSVertex *)FaceIter.PeekNextCVertex();
		tri->Point2 = pVtx->GetPoint();
		
		pVtx = (CLEDSVertex *)FaceIter.PeekNextCVertex();
		tri->Point3 = pVtx->GetPoint();
		
		pVtx = (CLEDSVertex *)FaceIter.PeekNextCVertex();
		// assert that this is a triangle
		if( p1stVtx != pVtx ) {
			printf( "THE STL FILE IS NOT COMPOSED OF STRICTLY TRIANGLES");
			assert( false );
		}

		sphere->triangle = curr_triangle;
		sphere->Children[0] = sphere->Children[1] = NULL;
		sphere->Parent = NULL;

		ConstructSphere( sphere , *tri );

		curr_triangle++;
		pFace = CListIterIntFace.PeekNext();
	}

}



class WorkingSet {

public:
	GSphere* *items;
	int size;
	int maximum_size;
	int nextItem;

	WorkingSet(int max_size)
	{
		size = 0;
		maximum_size = max_size;
		items = new GSphere*[max_size];
		for(int i=0;i<max_size;i++)
			items[i] = NULL;
		nextItem = max_size;
	}

	~WorkingSet()
	{
		delete [] items;
	}

	/**
	 * starts a new iterator for the list - pointing to the first item in the list
	 */
	int startIterator()
	{
		int i = 0;
		while( i < size ) {
			if( items[i] != NULL)
				return i;
			else
				i++;
		}
		return i;
	}

	
	/**
	 * gets the next item in the working list and 
	 */
	GSphere* getNext(int *itr , int *prev)
	{
		if( *itr < 0  || *itr >= size )
			return NULL;
		GSphere* next = items[*itr];

		int i = *itr;
		while( i < size ) {
			if( items[i] != NULL) {
				*itr = i;
				return next;
			}
			else
				i++;
		}

		*itr = size;
		return NULL;
	}

	/**
	 * removes an item from the working set
	 */
	GSphere* removeFromWorkingSet( int index )
	{
		assert( size > index || size >= 0 );

		GSphere* i = items[index];
		items[index] = NULL;
		return i;
	}

	
	int getSize()
	{
		return size;
	}

	void addItem(GSphere* item)
	{
		assert( size < maximum_size && size >= 0 && item != NULL);
		
		items[size] = item;
		size++;
	}

};


class WSIterator {
	int next;
	int prev;
	WorkingSet *wset;
public:
	WSIterator( WorkingSet *w )
	{
		wset = w;
		next = 0;
		prev = -1;
	}

	~WSIterator()
	{
		
	}

	void start()
	{
		next = 0;
		prev = -1;
	}

	GSphere* getNext()
	{
		while( next < wset->size ) {
			if( wset->items[next] != NULL) {
				prev = next;
				GSphere *ret = wset->items[next];
				next++;
				return ret;
			}
			else
				next++;
		}
		return NULL;
	}

	/**
	 * return the index of the item just returned by this iterator
	 */
	int getPrev()
	{
		return prev;
	}


};

/**
 * Find the next largest level of spheres
 *
 * @return	whether or not to proceed with the tree
 */
bool BubbleSpheres(void)
{
	double r=0;  // bubble radius
	WorkingSet *Working_List = new WorkingSet( num_spheres[sphere_level] );  //elements in working list
	int i;
	// TODO: eliminate hardcoding
	// 1) Calculate the average sphere radius and multiply it by 2.5
		int current_sphere=0;
		for(i=0;i<num_spheres[sphere_level];i++)
			r += Sphere_List[sphere_level][i].r;

		r /= num_spheres[sphere_level];  //the average sphere radius
		
	// TODO: is this right?
		r *= 2.5; // the radius of the merge sphere should be at least the size of two and a half spheres right next to each other

	// 2) Setup spheres to test
		for(i=0;i<num_spheres[sphere_level];i++) 
			Working_List->addItem( &Sphere_List[sphere_level][i] ); //Setup the list of available spheres to bubble
	
		num_spheres[sphere_level+1] = 0;
		Sphere_List[sphere_level+1] = (GSphere*)malloc( sizeof(GSphere) * num_spheres[sphere_level] );
		
	
		int currentNewSphere = 0;
		int num_merged = 0;
	// 3) Go thru each pairing of spheres and determine whether or not they should be paired
		// these iterators are basically pointers to the next item in the list and when you call getNext they are advanced
		// thus you save them before every call to getNext
		GSphere* sphereI;
		GSphere* sphereJ;
		WSIterator *itrI = new WSIterator(Working_List);
		WSIterator *itrJ = new WSIterator(Working_List);
		itrI->start();
		while( (sphereI = itrI->getNext()) != NULL )
		{
			// loop thru all other entries and see if you can bubble them in
			itrJ->start();

			while( (sphereJ = itrJ->getNext()) != NULL )
			{
				if( sphereI != sphereJ && CanBubble(r,sphereI,sphereJ) )
				{
					GenMergeSphere(sphereI,sphereJ, &Sphere_List[sphere_level+1][currentNewSphere] );  //New sphere is a merging of these two
					currentNewSphere++;
					Working_List->removeFromWorkingSet( itrJ->getPrev() );
					Working_List->removeFromWorkingSet( itrI->getPrev() );
					num_merged++;
					break;
				}
			}
		}

	// 4) Loop thru the remaining spheres and generate lonely spheres
		itrI->start();
		while( (sphereI = itrI->getNext()) != NULL )
		{
			GenLoneSphere(sphereI,&Sphere_List[sphere_level+1][currentNewSphere] );
			currentNewSphere++;
		}

	// 5) Reallocate space for the Sphere_List
		num_spheres[sphere_level+1] = currentNewSphere;
		realloc( Sphere_List[sphere_level+1] , sizeof(GSphere) * num_spheres[sphere_level+1] );
		printf("total num new bubbles = %d\n",num_spheres[sphere_level+1]);
		printf("memory allocated for %d spheres at level %d of size %d\n",num_spheres[sphere_level+1],(sphere_level+1), sizeof(GSphere) * num_spheres[sphere_level+1]);
		printf("num merged = %d\n",num_merged);

		sphere_level++;
		num_sphere_levels++;

	// 6) Clean up
		delete Working_List;
		delete itrI;
		delete itrJ;

	return( num_merged > 0 );
}



/**
 * Returns whether or not two spheres can be bubbled.
 * dS = runion^3 / ( ri^3 + rj^3 )
 * @param	i,j		the two spheres needing to be merged
 * @param	r		the 
 * @return	(ds <= r)
 */
int CanBubble(double r, GSphere* i, GSphere* j)//, int eiwl, int* Working_List)
{
	// TODO: you may want to change this back
	/*double runion=0,ds=0;
	double temp[3],temp2,temp3;
	int k=0;
	for(k=0;k<3;k++){
		temp[k]= j->Orig[k] - i->Orig[k];
	}
	temp2 = i->r;
	temp3 = j->r;
	runion = ( temp2 + temp3 + norm(temp) ) / 2.0;
	ds = (runion*runion*runion)/((temp2*temp2*temp2)+(temp3*temp3*temp3));
	if(ds<=r) 
		return(TRUE);
	else{
		return(FALSE);
	}*/

	double runion=0;
	double temp[3],temp2,temp3;
	int k=0;
	temp[0]= j->Orig[0] - i->Orig[0];
	temp[1]= j->Orig[1] - i->Orig[1];
	temp[2]= j->Orig[2] - i->Orig[2];
	
	temp2 = i->r;
	temp3 = j->r;
	runion = ( temp2 + temp3 + norm(temp) ) / 2.0;

	return( runion <= r );
}


/**
 * Generate one sphere from two mergesd spheres
 * will set Sphere_List[sphere_level+1][current_sphere] = the newly merged sphere
 * @param	i,j					the two spheres to be merged
 * @param	current_sphere		the sphere that will contain the merged spheres
 */
void GenMergeSphere(GSphere* i, GSphere* j, GSphere* new_sphere)
{
//	printf("\t\tIn GenMergeSphere, sphere_level=%d\n",sphere_level);
	double runion=0,ds=0;
	double temp[3],temp2,temp3;
	int k=0;
	
	// Calculate the radius of the sphere
		for(k=0;k<3;k++)
			temp[k] = j->Orig[k] - i->Orig[k];

		temp2 = i->r;
		temp3 = j->r;
		runion = (temp2+temp3+norm(temp))/2.0;
		new_sphere->r = runion;

	// Calculate the origin of the new sphere
		Normalize(temp);  //give me a unit vector in the correct direction
		for(k=0;k<3;k++){
			new_sphere->Orig[k] = i->Orig[k] + ((runion-temp2)*temp[k]);
		}
	
	//Update Parent/Child Data
		new_sphere->Children[0] = i;
		new_sphere->Children[1] = j;
		new_sphere->Parent = NULL;
		new_sphere->triangle = -1;

		i->Parent = new_sphere;
		j->Parent = new_sphere;

}


void PrintTestSphere(GSphere x)
{
	printf("Test Sphere:");
	PrintSphereData(x);
}

void PrintSphereData(GSphere x)
{
	printf("Sphere radius %f at:\n",x.r);
	for(int k=0;k<3;k++){
		printf("%f\t",x.Orig[k]);
	}
	printf("\n");
}

/**
 * Generate a new sphere and place it in the location new_sphere
 * the new sphere will be exactly the same as sphere i
 */
void GenLoneSphere(GSphere* i, GSphere* new_sphere)
{
	new_sphere->r = i->r;
	for(int k=0;k<3;k++)
	{
		new_sphere->Orig[k] = i->Orig[k];
	}
	new_sphere->Children[0] = i;
	new_sphere->Children[1] = i;
	new_sphere->Parent = NULL;
	new_sphere->triangle = -1;

	i->Parent = new_sphere;
	i->Parent = new_sphere;
}


/*This Creates the entire sphere tree in one shot */
void CreateSphereTree(void)
{
	printf("Creating Sphere Tree\n");
	int i=0;
	bool no_more_merged = false;
	while(i<MAX_SPHERE_LEVELS-1)
	{
		no_more_merged = !(BubbleSpheres());
		printf("Sphere Level %d of max= %d Finished!\n", i+1, MAX_SPHERE_LEVELS-1 );

		if(num_spheres[sphere_level] <= 1 || no_more_merged)
			break;

		i++;
	}
	printf("Sphere Tree Created!\n");
}



/* a=b */
void MakeEqual(double* a, double* b)
{
	for(int i=0;i<3;i++){
		a[i]=b[i];
	}
}

int GetMaxIndex(double* a)
{
	double x=a[0];
	int indmax=0;
	for(int i=1;i<3;i++){
		if(abs_d(a[i])>x){
			x=a[i];
			indmax=i;
		}
	}
	return(indmax);
}

/*projection of a onto b */
void VectorProj(double* a, double* b, double* C)
{
	double bb=dot(b,b);
	double ab=dot(a,b);
	ab/=bb;  //to save a scrap of memory and/or processing
	for(int i=0;i<3;i++){
		C[i]=ab*b[i];
	}
}

/* c=b-a */
void VectBetweenPoints(double* a, double* b, double* c)
{
	for(int i=0;i<3;i++){
		c[i]=b[i]-a[i];
	}
}

/*P2=P+V*/
void PointPlusVector(double* P, double* V, double* P2)
{
	for(int i=0;i<3;i++){
		P2[i]=P[i]+V[i];
	}
}

int f(double* V, double* C)
{
	double A[3];
	VectorProj(V,C,A);
	double s=dot(A,C);
	return(sign(s));
}


void printChildren(int i, int* kids){
	printf("%d's Children are: ",i);
	for(int j=0;j<2;j++)
	{
		printf("%d ",kids[j]);
	}
	printf("\n");
}




std::list<GSphere*> *testList;
void initTestList()
{
	testList = new std::list<GSphere*>;
}

void clearTestList()
{
	delete testList;
}

double d_minimum;

/**
 * d = |p-c| + r
 */
inline void update_dmin( double d )
{
	if( d_minimum > d )
		d_minimum = d;
}

double calcDistanceToCenter( GSphere *sphere , GPoint *point)
{

	double temp[3];
	temp[0] = sphere->Orig[0] - point->X[0];
	temp[1] = sphere->Orig[1] - point->X[1];
	temp[2] = sphere->Orig[2] - point->X[2];
	return norm(temp);
}


/**
 * Find closest sphere and returns the new dmin if it is closer
 */
bool addSphereToList( GSphere *sphere, GPoint *TestPoint)
{
	double d;
	assert( sphere != NULL  );
	
	d = calcDistanceToCenter( sphere , TestPoint );
	if( (d - sphere->r) <= d_minimum )
	{
		testList->push_back( sphere );
		update_dmin(d + sphere->r);
		return true;
	}
	else
		return false;
}


/**
 * Breadth first search of the tree. After being called the 
 * @param	top_level	the highest level of the tree - the root
 * @param	TestPoint	the point to test
 */
void FindClosestSpheres( int top_level , GPoint *TestPoint )
{
	bool complete,c1,c2;
	complete = false;
	GSphere *current;
	int i,size;
	
	d_minimum = calcDistanceToCenter( &Sphere_List[top_level][0] , TestPoint) + Sphere_List[top_level][0].r;
	
	for(i=0;i<num_spheres[top_level];i++)
		addSphereToList( &Sphere_List[top_level][i] , TestPoint );
	
	while( !complete )
	{
		complete = true;
		size = testList->size();
		for(i=0;i<size;i++)
		{
			assert( !testList->empty() );
			current = testList->front();
			testList->pop_front();
			
			if( current->Children[0] == NULL )
			{
				assert( current->Children[1] == NULL );
				assert( current->triangle != -1 );
				addSphereToList( current , TestPoint );
			} 
			else if( current->Children[0]  == current->Children[1] )
			{
				if( addSphereToList( current->Children[0] , TestPoint) )
					complete = false;
			}
			else
			{
				c1 = addSphereToList( current->Children[0] , TestPoint);
				c2 = addSphereToList( current->Children[1] , TestPoint);
				if( c1 || c2 )
					complete = false;
			}
		}
	}
}

/**
 * Determine the list of possible closest triangles
 * This updates the Test_Sphere_List to be a list of all the leafs
 * to check.
 */
octreeVtx * FindClosestDistance(GPoint *gTestPoint)
{
	initTestList();

	Point TestPoint;
	convertGPointToPoint( *gTestPoint , &TestPoint );
	
	double closest[3];  //The closest point
	
	// 1) find all the relevant spheres
	FindClosestSpheres( num_sphere_levels-1 , gTestPoint);

	ASSERT(testList->size() != 0);	// (Steve) if this is the case..something went wrong.
	
	// 2) Go thru list and find the closest point
	GSphere* tempS = testList->front();
	testList->pop_front();
	
	double dmin,tempD,tempP[3];
	
	// (Steve) Init'd to NULL
	octreeVtx *minvertex = NULL;
	octreeVtx *tempVertex = NULL;
	
	
	Vector foo; // (Steve)
	minvertex = minDistTri( Triangle_List[tempS->triangle].Point1 , Triangle_List[tempS->triangle].Point2 , Triangle_List[tempS->triangle].Point3 , TestPoint, foo);
	
	dmin = minvertex->dist;
	std::list<GSphere*>::iterator itr;
	
	for(itr = testList->begin(); itr != testList->end(); ++itr) 
	{
		tempS = *itr;
		assert( tempS->triangle != -1 );
		
		
		tempVertex = minDistTri( Triangle_List[tempS->triangle].Point1 , Triangle_List[tempS->triangle].Point2 , Triangle_List[tempS->triangle].Point3 , TestPoint, foo);
		
		tempD = tempVertex->dist;
		
		// (Steve) Added null check
		if(minvertex == NULL || (tempVertex != NULL && abs_d(dmin) > abs_d(tempD))) {
			MakeEqual(closest,tempP);
			dmin = tempD;
			// TODO: managing memory leaks
			// cause he allocated space for the vertex in the minDistTri function
			freeOctreeVtx( minvertex );
			minvertex = tempVertex;
		}
	}
	
	clearTestList();
	
	return(minvertex);
}


void CalcRotationMatrix(double N[3], double R[3][3])
{
		//http://astronomy.swin.edu.au/~pbourke/geometry/rotate/
		/*  U = (a,b,c)
		 d = sqrt(b^2 + c^2) 
		cos(t) = c/d, sin(t)=b/d   (shows the simple proof of that)
		Rx =
		[    1,    0,    0]
		[    0,  c/d, -b/d]
		[    0,  b/d,  c/d]
		Ry =
		[  d,  0, -a]
		[  0,  1,  0]
		[  a,  0,  d]
		R =
		[      d, -a*b/d, -a*c/d]
		[      0,    c/d,   -b/d]
		[      a,      b,      c]
		*/
	double a,b,c,d;
	a=N[0];
	b=N[1];
	c=N[2];
	d=sqrt(b*b+c*c);
	double temp[3][3]={{ d, -a*b/d, -a*c/d},{0,    c/d,   -b/d},{a,      b,      c}};
	for(int i=0;i<3;i++){
		for(int j=0;j<3;j++){
			R[i][j]=temp[i][j];
		}
	}
}



void CalcInverseRotation(double R[3][3])
{
	//Its just the transpose
	double temp[3][3];
	int i,j;
	for(i=0;i<3;i++){
		for(j=0;j<3;j++){
			temp[j][i]=R[i][j];
	}}
	for(i=0;i<3;i++){
		for(j=0;j<3;j++){
			R[i][j]=temp[i][j];
	}}
}
void RotateVector(double R[3][3], double V[3])
{
	int i,j;
	double temp[3];
	for(i=0;i<3;i++){
		temp[i]=0;
		for(j=0;j<3;j++){
			temp[i]+=R[i][j]*V[j];
		}
	}
	for(i=0;i<3;i++) V[i]=temp[i];
}

// used to randomly allocate a color to a sphere
int randomNum = 1;
/**
 * Render on sphere
 */
void renderSphere(GSphere* sphere)
{
	if( sphere == NULL )
		return;

	GLfloat red = 0;
	GLfloat blue = 0;
	GLfloat green = 0;

	// Determine Color for Sphere
	if (randomNum % 3 == 1) {
		red = 1.0;
	}
	else if (randomNum % 3 == 2) {
		green = 1.0;
	}
	else {
		blue = 1.0;
	}
	randomNum++;

	glPushMatrix();

	// Set Current Color
	glColor3f (red, green, blue);

	// Translate to sphere Position
	glTranslatef ((GLfloat) sphere->Orig[0], (GLfloat) sphere->Orig[1], (GLfloat) sphere->Orig[2]);

	// Draw sphere - the last two numbers are the number of stacks and slices
	glutWireSphere((GLdouble) sphere->r , 20 , 20);

	glPopMatrix();
}

/**
 * Render all the spheres in the tree.
 */
void renderSphereTree( )
{
	randomNum = 0;
	int numspheres = num_spheres[level_to_render];
	for(int j=0;j<numspheres;j++)
	{

		renderSphere( &Sphere_List[level_to_render][j] );
	}

}





// Different centroid code

void ConstructSphere(GSphere *sphere, GTriangle &triangle)
{

	double radVector[3];

	// calculate centroid of the triangle
	sphere->Orig[0] = ( triangle.Point1[0] + triangle.Point2[0] + triangle.Point3[0] ) / 3.0;
	sphere->Orig[1] = ( triangle.Point1[1] + triangle.Point2[1] + triangle.Point3[1] ) / 3.0;
	sphere->Orig[2] = ( triangle.Point1[2] + triangle.Point2[2] + triangle.Point3[2] ) / 3.0;

	
	double r1,r2,r3;
	radVector[0] = triangle.Point1[0] - sphere->Orig[0];
	radVector[1] = triangle.Point1[1] - sphere->Orig[1];
	radVector[2] = triangle.Point1[2] - sphere->Orig[2];

	r1 = norm(radVector);
	
	radVector[0] = triangle.Point2[0] - sphere->Orig[0];
	radVector[1] = triangle.Point2[1] - sphere->Orig[1];
	radVector[2] = triangle.Point2[2] - sphere->Orig[2];
	
	r2 = norm(radVector);
	
	radVector[0] = triangle.Point3[0] - sphere->Orig[0];
	radVector[1] = triangle.Point3[1] - sphere->Orig[1];
	radVector[2] = triangle.Point3[2] - sphere->Orig[2];
	
	r3 = norm(radVector);

	// get the maximum possible radius
	if( r1 > r2 && r1 > r3) 
		sphere->r = r1;
	else if( r2 > r1 && r2 > r3) 
		sphere->r = r2;
	else if( r3 > r2 && r3 > r1) 
		sphere->r = r3;
	
}

//********************************************************************

void triangle_centroid_3d ( double centroid[] ,const double x[], const double y[], const double z[] )

//********************************************************************
//
//  Purpose:
//
//    POLYGON_CENTROID_3D computes the centroid of a polygon in 3D.
//
//  Method:
//
//    The centroid is the area-weighted sum of the centroids of
//    disjoint triangles that make up the polygon.
// 
//  Modified:
//
//    20 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Adrian Bowyer and John Woodwark,
//    A Programmer's Geometry,
//    Butterworths, 1983.
//
//  Parameters:
//
//
//    Input, double X[3], Y[3], Z[3], the coordinates of the vertices.
//
//    Output, double POLYGON_CENTROID_3D[3], the coordinates of the centroid.
//
{

	centroid[0] = ( x[0] + x[1] + x[2] ) / 3.0E+00;
	centroid[1] = ( y[0] + y[1] + y[2] ) / 3.0E+00;
	centroid[2] = ( z[0] + z[1] + z[2] ) / 3.0E+00;
}