// CPP program to disect a polyhedra in vrml format
// and output it in jvx format

#include <vector>
#include <list>
#include <map>
#include <valarray>
#include <iostream>
#include <strstream>
#include <string>
#include <cstring>
#include <cstdlib>

/*
#define PRINT_INT
#define PRINT_MERGE
#define PRINT_HIDE
*/

/**
 * couple of templates to treat vectors cyclically
 */

template<class T> void rotate_to_set(vector<T>&,const T,const T);
template<class T> vector<T> merge_at(vector<T>,vector<T>,const T,const T);
template<class T> void clean_spikes(vector<T> &);


template<class T> void rotate_to_set(vector<T>& vec,const T a,const T b)
{
	vector<T>::iterator i;
	vector<T>::iterator j;
	/* need to find i,j which are consecutive */
	i = vec.begin();
	while(i!=vec.end())
	{
		i = find(i,vec.end(),a);
		if(i == vec.end() ) break;
		if(i == vec.begin() )
		{
			j = i+1; if( *j == b) break;
			j = vec.end()-1; if( *j == b ) break;
		}
		else if(i == vec.end()-1)
		{
			j = i-1; if( *j == b) break;
			j = vec.begin(); if( *j == b ) break;
		}
		else
		{
			j = i+1; if( *j == b) break;
			j = i-1; if( *j == b) break;
		}
		++i;
	}
	if(i==vec.end() || j == vec.end()) 
	{
		cerr << "Rotate_to_set failed 198\n!";
		return;
	}
	if(i == vec.begin())
	{
		if(j == i + 1)
		{
			// the state we want do nothing
		}
		else if(j == vec.end()-1 )
		{
			rotate(vec.begin(),i+1,vec.end());
			reverse(vec.begin(),vec.end());
		}
		else
		{
			cerr << "Wierdness in cycle\n";
		}
	}
	else if(j == vec.begin())
	{
		if(i==j+1)
		{
			rotate(vec.begin(),i+1,vec.end());
			reverse(vec.begin(),vec.end());
		}
		else if(i==vec.end()-1)
		{		
			rotate(vec.begin(),i,vec.end());
		}
		else
		{
			cerr << "Wierdness in cycle\n";
		}
	}
	else if(i<j)
	{
			rotate(vec.begin(),i,vec.end());
	}
	else
	{
			rotate(vec.begin(),i+1,vec.end());
			reverse(vec.begin(),vec.end());
	}

}

template<class T> vector<T> merge_at(vector<T> vec1,vector<T> vec2,const T a,const T b)
{
	rotate_to_set(vec1,a,b);
	rotate_to_set(vec2,b,a);
//cout << "c1 (" << c1 << ")\n";
//cout << "c2 (" << c2 << ")\n";
	vector<T> c(0);
	vector<T>::iterator i;
	for(i=vec1.begin()+1;i!=vec1.end();++i) c.push_back(*i);
	for(i=vec2.begin()+1;i!=vec2.end();++i) c.push_back(*i);
	return c;
}

template<class T> void clean_spikes(vector<T> &vec)
{
	int size = vec.size();
	vector<T>::iterator i;
	bool more_to_do = true;

	while(more_to_do)
	{
	    more_to_do = false;
	    if(vec.size() <= 2)
	    {
		vec.clear();
		return;
	    }
	    for(i=vec.begin();i!=vec.end();++i)
	    {
		vector<T>::iterator previ;
		vector<T>::iterator nexti;
		previ = ( i==vec.begin() ? vec.end()-1 : i-1);
		nexti = (i==vec.end()-1 ? vec.begin() : i+1);
		if(*previ == *nexti)
		{
//cerr << "\nClean spikes (" << vec << ") -> (";
			if(i==vec.begin())
			{
				vec.erase(vec.end()-1);
				vec.erase(i);
			}
			else if(i==vec.end()-1)
			{
				vec.erase(vec.end()-1);
				vec.erase(vec.begin());
			}
			else
			{
				vec.erase(previ,nexti);
			}
//cerr << vec << ")\n";
			more_to_do = true;
			break;
		}
		else if(*i == *nexti)
		{
			vec.erase(i);
			more_to_do = true;
			break;
		}						
	   }
	}
/*
	if(vec.size() == 4)
	{
		if((vec[0] == vec[1] && vec[2] == vec[3])
		 ||(vec[0] == vec[2] && vec[1] == vec[3])
		 ||(vec[0] == vec[3] && vec[1] == vec[2]))
			vec.clear();
	}
*/
}
	
class Point
{
	double coords[3];
public:
	const static double limit=0.000001;
	Point(double xx=0.0,double yy=0.0,double zz=0.0)
	{
		coords[0]=xx; coords[1]=yy; coords[2]=zz;
	};
//	Point(const Point &p) 
	friend ostream& operator<<(ostream&,const Point&);
	friend Point operator+ (const Point&,const Point&);
	friend Point operator- (const Point&,const Point&);
	friend Point operator* (double,const Point&);
	friend Point operator^ (const Point&,const Point&); // vector product
	friend double dot(const Point&,const Point&); // dot product
	friend bool parallel(const Point&,const Point&); // whether two vectors are parallel
	friend double lensq(const Point&); // square of the length
	bool unit();
	void flip();
};
inline ostream& operator<<(ostream& os,const Point& p)
{
	os << p.coords[0] << " " << p.coords[1] << " " << p.coords[2];
}
inline Point operator+ (const Point& a,const Point& b)
{
	Point p(a.coords[0]+b.coords[0],a.coords[1]+b.coords[1],a.coords[2]+b.coords[2]);
	return p;
}
inline Point operator- (const Point& a,const Point& b)
{
	Point p(a.coords[0]-b.coords[0],a.coords[1]-b.coords[1],a.coords[2]-b.coords[2]);
	return p;
}
inline Point operator* (double d,const Point& a)
{
	Point p(d* a.coords[0],d*a.coords[1],d*a.coords[2]);
	return p;
}
inline double dot(const Point& a,const Point& b)
{
	return a.coords[0]*b.coords[0]+a.coords[1]*b.coords[1]+a.coords[2]*b.coords[2];
}
inline double lensq(const Point& a)
{
	return a.coords[0]*a.coords[0]+a.coords[1]*a.coords[1]+a.coords[2]*a.coords[2];
}
inline Point operator^ (const Point& a,const Point& b)
{
	Point p(a.coords[1]*b.coords[2]-a.coords[2]*b.coords[1],
		a.coords[2]*b.coords[0]-a.coords[0]*b.coords[2],
		a.coords[0]*b.coords[1]-a.coords[1]*b.coords[0]);
	return p;
}

inline bool Point::unit()
{
	double len = sqrt(coords[0]*coords[0] + coords[1] * coords[1] + coords[2] * coords[2]);
	if(len!=len || len < limit) 
	{
		coords[0] = coords[1] = coords[2] = 0.0;
		return false;
	}
	coords[0] /= len;
	coords[1] /= len;
	coords[2] /= len;
	return true;
}
inline void Point::flip()
{
	coords[0] = -coords[0];
	coords[1] = -coords[1];
	coords[2] = -coords[2];
}
inline bool parallel(const Point& a,const Point& b)
{
	return lensq(a^b) < Point::limit;
}

/**
 * A functional class 
 * returns true if lensq(p) < dist - Point::limit
 */

class DistPred
{
	double dist;
public:
	DistPred(double d) {dist = d;};
	bool operator() (const Point &p) { bool flag = lensq(p) < dist - Point::limit; 
//	cerr << "DIstPred " << flag << endl;
	return flag;}
};

/**
 *	Face: a ordered list of indicies to the points array
 */

class Face
{
public:
	vector<int> indices;
	Face(vector<int> in) {indices = in;}
	friend ostream& operator<<(ostream&,const Face&);
	friend void printFacePoints(ostream &,const Face&);
};
ostream& operator<<(ostream& os,const Face& f)
{
		for(int j=0;j < f.indices.size();++j)
			os << f.indices[j] << " ";
	return os;
}

/**
 * A colour value
 */

class ColourInfo
{
	int r,g,b;	// rgb colours 0..255
public:
	ColourInfo(int rr=-1,int gg=-1,int bb=-1) { r = rr; g=gg; b=bb; }
	ColourInfo(const ColourInfo &c) { r = c.r; g=c.g; b=c.b; }
	void setColour(int rr,int gg,int bb) { r = rr; g = gg; b = bb; }
	friend ostream& operator<<(ostream& os,const ColourInfo& c);
};
ostream& operator<<(ostream& os,const ColourInfo& c)
{
	os << c.r << " " << c.g << " " << c.b;
}

/**
 * A plane
 */

class Plane
{
public:
	Point normal;	// outward pointing normal
	double d;	// minimum distance from origin
	Plane() : normal() { d = 0.0; }	// uses default constructor for normal
	Plane(const Point &n,double dd) : normal(n) 
	{ 
		if( dd < 0 ) 
		{ d = -dd; normal.flip(); }
		else d = dd;
	}
	double dist(const Point &p) const { return (dot(p,normal) - d); } // the distance from the plane along normal
};

/**
 * A face with info on facial plane
 */

class NFace : public Face, public Plane 
{
	void calcPlane();
public:
	NFace(vector<int> in) : Plane() , Face(in) { calcPlane(); };
	NFace(vector<int> in,const Plane &p) : Plane(p) , Face(in) {};
	NFace(vector<int> in,const Point &p,double dd) : Plane(p,dd) , Face(in) {};
	NFace(vector<int> in,const NFace &f2) : Face(in), Plane(f2) {};
};

/**
 * A face with colour info and plane info
 */

class NCFace : public NFace, public ColourInfo
{
public:
	NCFace(vector<int> in) 
		: NFace(in) , ColourInfo() {};
	NCFace(vector<int> in,const Plane &p) 
		: NFace(in,p) , ColourInfo() {};
	NCFace(vector<int> in,const Point &p,double dd) 
		: NFace(in,p,dd) , ColourInfo() {};

	NCFace(vector<int> in,const ColourInfo &c)
		: NFace(in) , ColourInfo(c) {};
	NCFace(vector<int> in,const Plane &p,const ColourInfo &c)
		: NFace(in,p) , ColourInfo(c) {};
	NCFace(vector<int> in,const Point &p, double dd,const ColourInfo &c)
		: NFace(in,p,dd) , ColourInfo(c) {};

	NCFace(vector<int> in,const NCFace &f2)
		: NFace(in,f2) , ColourInfo(f2) {};

};

/**
 * The end class in the Face hieraracy
 * 	good determins whether we print the face at the end
 *	TODO: add list of faces which have been done
 */

class MyFace : public NCFace
{
public:
	bool good;	// whether we want to print in final ouput
	friend bool no_good(const MyFace &f);
	int originalIndex,Findex;
	vector<int> linkingFaces;

	void setFindex(int i) { Findex = i; };
	MyFace(vector<int> in) 
		: NCFace(in) { good = true; };
	MyFace(vector<int> in,const Plane &p) 
		: NCFace(in,p) { good = true; };
	MyFace(vector<int> in,const Point &p,double dd) 
		: NCFace(in,p,dd) { good = true; };

	MyFace(vector<int> in,const ColourInfo &c)
		: NCFace(in,c) {good = true; };
	MyFace(vector<int> in,const Plane &p,const ColourInfo &c)
		: NCFace(in,p,c) {good = true; };
	MyFace(vector<int> in,const Point &p, double dd,const ColourInfo &c)
		: NCFace(in,p,dd,c) {good = true; };

	MyFace(vector<int> in,const MyFace &f2)
		: NCFace(in,f2) {good = true; originalIndex = f2.originalIndex; };
};
bool no_good(const MyFace &f) { return !f.good; };


/**
 * EdgeKey used for constructing a hash table for edges
 */

class EdgeKey
{
public:
	int a,b;
	EdgeKey(int aa,int bb) { if(aa<bb) { a=aa; b=bb; } else {a=bb;b=aa;} };
	friend bool operator<(const EdgeKey& e1,const EdgeKey& e2);
};

inline bool operator<(const EdgeKey& e1,const EdgeKey& e2)
{
	return (e1.a < e2.a ) ||
		(e1.a == e2.a && e1.b < e2.b );
}

class Edge
{
public:
	bool good;
	vector<MyFace *> faces;
	Edge() {};
	void addFace(MyFace *f) { faces.push_back(f); }
};


/******************* Global variable *******************/

vector<Point> points;	// the list of points needs to be up here as used by NFace
list<MyFace> faces;	// the list of faces
list<MyFace> newfaces;	// new faces added
vector<MyFace> originalfaces;	// the list of faces
map<EdgeKey,Edge> edges;
string polyname;	// contains the name of the polyhedron

void NFace::calcPlane()
{
	if(indices.size()<3)
	{
		cerr << "Face with only two points found\n";
		d = 0.0;
	}
	Point &f1a = points[indices[0]];
	Point &f1b = points[indices[1]];
	Point &f1c = points[indices[2]];
	
	Point a1 = f1a-f1b;
	Point b1 = f1c-f1b;
	normal = a1^b1;	// normal for face1
	if( !normal.unit() )
	{
		cerr << "Normal for face is zero\n";
		d = 0.0;
		return;	
	}
	d = dot(normal,f1b);	// distance from origin for face1
	if(d<0) { d = -d; normal.flip(); }
}

/******************* input routines *******************/

void readPoints()
{
	Point buf;
	string str;
	char ch;

	getline(cin,str);
	getline(cin,str);
	getline(cin,str);
	// read until ) is reached
	while(cin>>ch && ch!=')');
	cin.get(ch);
	getline(cin,polyname);
	// now gobble until point[
	while(1)
	{

		getline(cin,str);
		if(str.find("point[") != string::npos) break;
	}

	while(1)
	{
		double x,y,z;
		if( (cin >> x ) )
		{
			cin >> y >> z >> ch;
//			cout << "read " << x << " " << y << " " << z << " " << ch << endl;
			points.push_back(Point(x,y,z));
		}
		else
		{
			cin.clear();
			cin>>ch;
			if(ch == '#' )
			{
				getline(cin,str);
			}
			else
			{
				break;
			}
		}					
	}
}

/**
 * readFaces(bool joinFlag)
 * if join flag is true then each face is like
 * 0,1,12,-1,1,6,12,-1,6,9,12,-1,9,2,12,-1,2,0,12,-1,
 */

void readFaces(int joinFlag)
{
	string str;
	vector<int> indices;
	char ch;
	int a;

	// first gobble until point[

	while(1)
	{

		getline(cin,str);
		if(str.find("coordIndex[") != string::npos) break;
	}
	if(joinFlag == 1)
	{
		while(1)
		{
			// 0,1,12,-1,1,6,12,-1,6,9,12,-1,9,2,12,-1,2,0,12,-1,
			getline(cin,str);
			istream *inptr;
			inptr = new istrstream(str.c_str());
			istream &input = *inptr;
			if( !(input >> a >> ch) ) break; // catch the end of the faces
			indices.push_back(a);
			while( input >> a >> ch)
			{
				if(a != -1) indices.push_back(a);
			}
//			delete input;
			vector<int> goodind;
			for(int i=0;i<indices.size();++i)
			{
				// first check that this index not already done
				int j;
				bool flag = true;
				for(j=0;j<i;++j)
				{
					if(indices[i] == indices[j]) flag = false;
				}
				// now check it only occurs twice
				int count = 0;
				for(j=0;j<indices.size();++j)
				{
					if(indices[i] == indices[j]) ++count;
				}
				if(flag && count <= 2)
					goodind.push_back(indices[i]);
			}
//cerr << "Original indices: ";
//for(int i=0;i<indices.size();++i) { cerr << indices[i] << " "; }
//cerr << "\nOriginal indices: ";
//for(int i=0;i<goodind.size();++i) { cerr << goodind[i] << " "; }
//cerr << endl;
			faces.push_back(MyFace(goodind));
			originalfaces.push_back(MyFace(goodind));
			MyFace &createdface = faces.back();
			createdface.originalIndex = originalfaces.size()-1;
			goodind.clear();
			indices.clear();
		}
	}
	else if(joinFlag == 2 || joinFlag == 3)
	{
		/* for star faces want to produce a non intersecting outline */

		while(1)
		{
			// 0,1,12,-1,1,6,12,-1,6,9,12,-1,9,2,12,-1,2,0,12,-1,
			getline(cin,str);
			istream *inptr;
			inptr = new istrstream(str.c_str());
			istream &input = *inptr;
//			istream input = istrstream(str.c_str());
			if( !(input >> a >> ch) ) break; // catch the end of the faces
			indices.push_back(a);
			while( input >> a >> ch)
			{
				if(a != -1) indices.push_back(a);
			}
//			delete input;
			vector<int> goodind;
			for(int i=0;i<indices.size();++i)
			{
				// first check that this index not already done
				int j;
				bool flag = true;
				for(j=0;j<i;++j)
				{
					if(indices[i] == indices[j]) flag = false;
				}
				// now check it only occurs twice
				int count = 0;
				for(j=0;j<indices.size();++j)
				{
					if(indices[i] == indices[j]) ++count;
				}
				if(flag && count <= 2)
					goodind.push_back(indices[i]);
			}
//cerr << "Original indices: ";
//for(int i=0;i<indices.size();++i) { cerr << indices[i] << " "; }
//cerr << "\nOriginal indices: ";
//for(int i=0;i<goodind.size();++i) { cerr << goodind[i] << " "; }
//cerr << endl;
			// Now have just the verticies

			// construct a coordinate system
			MyFace f(goodind);
			Point origin = f.d * f.normal;
			Point u = points[goodind[0]] - origin;
			Point v = f.normal ^ u;
			// set of intersection points
			vector<Point> intersectPoints;
			int size = goodind.size();
			for(int i=0;i<size;++i)
			{
				int nexti = i+1;
				if(nexti == size) nexti = 0;
				Point a = points[goodind[i]] - origin;
				Point b = points[goodind[nexti]] - origin;
				double ax = dot(a,u),ay = dot(a,v), bx = dot(b,u), by = dot(b,v);
				
				for(int j=i+1;j<size;++j)
				{
					int nextj = j+1;
					if(nextj == size) nextj = 0;
					Point c = points[goodind[j]] - origin;
					Point d = points[goodind[nextj]] - origin;
					double cx=dot(c,u),cy=dot(c,v),dx= dot(d,u),dy=dot(d,v);
					double matA = -bx+ax,matB = dx-cx,matC=-by+ay,matD=dy-cy;
					double det = matA*matD - matB*matC;
					if(det < Point::limit && det > -Point::limit) continue;
					double lambda = (matD*(ax-cx)-matB*(ay-cy))/det;
					double mu = (-matC*(ax-cx)+matA*(ay-cy))/det;
					if(lambda < 0.0+Point::limit || lambda > 1.0-Point::limit) continue;
					Point p = origin + (1-lambda)*a + lambda*b;
					intersectPoints.push_back(p);
				}
			}
			if(intersectPoints.size() == 0)
			{
				faces.push_back(MyFace(goodind));
				originalfaces.push_back(MyFace(goodind));
				MyFace &createdface = faces.back();
				createdface.originalIndex = originalfaces.size()-1;

			}
			else
			{
				vector<Point>::iterator k,l;
				vector<int> intersectInd(intersectPoints.size());
				vector<Point> junk;
				int kk,ll;

				double dist1 = lensq(intersectPoints[0]);
				double dist2 = lensq(intersectPoints[1]);
				double maxdist = (dist1>dist2 ? dist1 : dist2);
				double mindist = (dist1>dist2 ? dist2 : dist1);

				if( maxdist - mindist > Point::limit)
				{
					DistPred mydp(maxdist);	
					remove_copy_if(intersectPoints.begin(),intersectPoints.end(),back_inserter(junk),mydp);
					intersectPoints.clear();
					copy(junk.begin(),junk.end(),back_inserter(intersectPoints));
				}


				// add the intersectPoints to points
//cerr << "Interior pts " << endl;
				
				for(k=intersectPoints.begin(),kk=0;k!=intersectPoints.end();++k,++kk)
				{
					for(l=points.begin(),ll=0;l!=points.end();++l,++ll)
					{
						if(lensq(*k-*l) < Point::limit )
						{
							intersectInd[kk] = ll;
							break;
						}
					}
					if(l==points.end())
					{
						intersectInd[kk] = ll;
						points.push_back(intersectPoints[kk]);
					}
//cerr << intersectInd[kk] << ' ';
				}
//cerr << "\nFinal ";
				int lastvert = -1;
				int curextvert = 0;
				int curintvert = -1;
				vector<int> finalind;
				do
				{
//cerr << goodind[curextvert] << ' ';
					finalind.push_back(goodind[curextvert]);
					// find the closest interior vert to ext vert
					Point &a = points[goodind[curextvert]];
					double mindist = 10.0;
					int minind = -1;
//cerr << "curintvert " << curintvert << endl;
					for(int i =0;i<intersectPoints.size();++i)
					{
						if(i==curintvert
						 || intersectInd[i] == intersectInd[curintvert] ) 
							continue;
						if( lensq(a-intersectPoints[i]) < mindist)
						{
							mindist = lensq(a-intersectPoints[i]);
							minind = i;
						}
					}
					curintvert = minind;
//cerr << "minind " << minind << endl;
//cerr << intersectInd[minind] << ' ';
					finalind.push_back(intersectInd[minind]);

					// find the closest exterior to interior

					Point &b = intersectPoints[curintvert];
					mindist = 10.0;
					minind = -1;
					for(int i =0;i<goodind.size();++i)
					{
						if(i==curextvert) continue;
						if( lensq(b-points[goodind[i]]) < mindist)
						{
							mindist = lensq(b-points[goodind[i]]);
							minind = i;
						}
					}
					curextvert = minind;
				} while(curextvert != 0);
//cerr << endl;
				if(joinFlag == 2)
				{
					faces.push_back(MyFace(finalind));
					originalfaces.push_back(MyFace(finalind));
					MyFace &createdface = faces.back();
					createdface.originalIndex = originalfaces.size()-1;

				}
				else
				{
					int finalsize = finalind.size();
					vector<int> interiorverts;
					for(int i=1;i<finalsize;i+=2)
						interiorverts.push_back(finalind[i]);
					faces.push_back(MyFace(interiorverts));
					originalfaces.push_back(MyFace(interiorverts));
					MyFace &createdface = faces.back();
					createdface.originalIndex = originalfaces.size()-1;
					vector<int> starverts;
					starverts.push_back(finalind[0]);
					starverts.push_back(finalind[1]);
					starverts.push_back(finalind[finalsize-1]);
					faces.push_back(MyFace(starverts,createdface));
					for(int i=1;i<finalsize-2;i+=2)
					{
						starverts.clear();
						starverts.push_back(finalind[i]);
						starverts.push_back(finalind[i+1]);
						starverts.push_back(finalind[i+2]);
						faces.push_back(MyFace(starverts,createdface));
					}				

				}
				finalind.clear();
			}
			goodind.clear();
			indices.clear();
		}
	}
	else
	{
		while(cin >> a >> ch)
		{
			if(a != -1) indices.push_back(a);
			else
			{
				faces.push_back(MyFace(indices));
				originalfaces.push_back(MyFace(indices));
				indices.clear();
			}
		}
	}
	cin.clear();
}


/*********************** output routines ****************/

void print_head()
{
	cout << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\" standalone=\"yes\"?>\n";
	cout << "<jvx-model>\n";
	cout << "	<meta generator=\"JavaView v.2.00.a13\"/>\n";
	cout << "	<meta date=\"Tue Jan 23 15:35:13 GMT+00:00 2001\"/>\n";
	cout << "	<version type=\"dump\">0.02</version>\n";
	cout << "	<title>" << polyname << "</title>\n";
	cout << " <authors>\n";
	cout << "  <author>\n";
	cout << "   <firstname>Richard</firstname>\n";
	cout << "   <lastname>Morris</lastname>\n";
	cout << "   <affiliation>\n";
	cout << "    <organization>Dept Statistics, University of Leeds</organization>\n";
	cout << "    <address>\n";
	cout << "     <line>Leeds</line>\n";
	cout << "     <line>LS2 9JT</line>\n";
	cout << "    </address>\n";
	cout << "   </affiliation>\n";
	cout << "   <email>webmaster@pfaf.org, rjm@amsta.leeds.ac.uk</email>\n";
	cout << "   <url>www.pfaf.org</url>\n";
	cout << "  </author>\n";
	cout << " </authors>\n";
	cout << "       <description>\n";
	cout << "		<abstract>A uniform polyhedra.</abstract>\n";
	cout << "	        <detail>\n";
	cout << "	                   This surface was evolved while working on the\n";
	cout << "	        </detail>\n";
	cout << "	        <keywords>\n";
	cout << "	            <keyword>uniform polyhedra</keyword>\n";
	cout << "		    <keyword>" << polyname << "</keyword>\n";
	cout << "	        </keywords>\n";
	cout << "	        <msc2000>\n";
	cout << "	                   <primary>52B10</primary>\n";
	cout << "	                   <secondary>52B99</secondary>\n";
	cout << "	                   <secondary>51M20</secondary>\n";
	cout << "	                </msc2000>\n";
	cout << "	                <software>Kalidio data for polyhedra triangulation by rjm</software>         - optional -\n";
	cout << "       </description>\n";
	cout << "	<geometries>\n";
	cout << "		<geometry name=\"" << polyname << "\">\n";
	cout << "			<pointSet point=\"hide\" dim=\"3\">\n";
}

void print_mid()
{
	cout << "			</pointSet>\n";
	cout << "			<faceSet face=\"show\" edge=\"show\" color=\"show\">\n";
}

void print_tail()
{
	cout << "			</faceSet>\n";
	cout << "		</geometry>\n";
	cout << "	</geometries>\n";
	cout << "</jvx-model>\n";
}

void print_points()
{
	typedef vector<Point>::const_iterator VI;

	cout << "				<points>\n";
	for(VI i=points.begin();i!=points.end();++i)
	{
		Point p = *i;
		cout << "<p>" << p << "</p>\n";
	}
	cout << "				</points>\n";
}

void print_faces()
{
	int ii,count=0;
	list<MyFace>::const_iterator i;

	cout << "				<faces>\n";
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{
		MyFace f = *i;
		if( f.good && f.indices.size() > 0 )
		{
			Face &ff = f;
			cout << "<f>" << ff << "</f>\n";
			++count;
		}
	}
	cout << "				</faces>\n";
	cerr << "Printed " << count << " faces\n";
}

void print_face_colours()
{
	int ii;
	list<MyFace>::const_iterator i;

	cout << "				<colors>\n";
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{
		MyFace f = *i;
		if( f.good )
		{
			ColourInfo &c = f;
			cout << "<c>" << c << "</c>\n";
		}
	}
	cout << "				</colors>\n";
}

void printFacePoints(ostream & os,const MyFace& f)
{
		os << "Face:\n";
		for(int j=0;j < f.indices.size();++j)
		{
			os << points[f.indices[j]] << endl;
		}
}

void print(const bool colourFlag)
{
	print_head();
	print_points();
	print_mid();
	print_faces();
	if(colourFlag)
		print_face_colours();
	print_tail();
}

/********************** intersection routines *******************/


/**
 *	bool intersectPlane(const Point &v,const double d,Face &f)
 *	returns true if the face iteresects the plane
 *	sideeffects: will add and delete faces from the list of faces if an intersection is found
 *		will also add points
 */

bool intersectPlane(MyFace &f,const MyFace &f2)
{
	// find first point not on line

	int j;	// primary loop variable
	int size = f.indices.size();
	int signs[size];
	double dists[size];
	bool plusFlag=false,minusFlag=false;
	for(j=0;j < size;++j)
	{
		double dist0;
		dist0 = f2.dist(points[f.indices[j]]);
		dists[j] = dist0;
		if( dist0 > Point::limit ) { signs[j] = 1; plusFlag = true; }
		else if( dist0 < - Point::limit ) { signs[j] = -1; minusFlag = true; }
		else signs[j] = 0;
	}
	if(!plusFlag || !minusFlag) return false;

	// check second face actuall crosses first face 
/*
	plusFlag=false,minusFlag=false;
	for(j=0;j < f2.indices.size();++j)
	{
		double dist0;
		dist0 = f.dist(points[f2.indices[j]]);
		if( dist0 > Point::limit ) { plusFlag = true; if(minusFlag) break; }
		else if( dist0 < - Point::limit ) { minusFlag = true; if(plusFlag) break; }
	}
	if(!plusFlag || !minusFlag) return false;
*/
	// We know there is some intersection

	f.good = false;

#ifdef PRINT_INT
cerr << "Intersection with plane face " << f.indices << " points ";
for(j=0;j < size;++j) cerr << dists[j] << " " ;
cerr << endl;
cerr << "Plane (" << f.normal << ") = " << f.d << " (" << f2.normal << ") = " << f2.d << endl;
//	printFacePoints(cerr,f);


//cerr << "Initial ";
//for(j=0;j<size;++j) cerr << signs[j] << ":" << f2.dist(points[f.indices[j]]) << " ";
//cerr << endl;
#endif

	// go through the list and add any new verticies

	int numAdditionalVertices = 0;
	int nextj;
	for(j=0;j<size;++j)
	{
		if(j==size-1) nextj = 0;
		else	      nextj = j+1;
		if((signs[j] == -1 && signs[nextj] == 1 ) 
		 ||(signs[j] == 1 && signs[nextj] == -1 ))
			++numAdditionalVertices;
	}
#ifdef PRINT_INT
//cerr << "additional " << numAdditionalVertices << endl;
#endif
	int size2 = size+numAdditionalVertices;
	int signs2[size2+1];
	int indicies2[size2+1];
	int i = 0;
	for(j=0;j<size;++j)
	{
		if(j==size-1) nextj = 0;
		else	      nextj = j+1;
		signs2[i] = signs[j];
		indicies2[i] = f.indices[j];
		++i;
		if((signs[j] == -1 && signs[nextj] == 1 ) 
		 ||(signs[j] == 1 && signs[nextj] == -1 ))
		{
			double dist1 = dists[j];
			double dist2 = dists[nextj];
			Point p = (dist1/(dist1-dist2)) * points[f.indices[nextj]]
				+ (dist2/(dist2-dist1)) * points[f.indices[j]];
			vector<Point>::const_iterator pi;
			int k;
			for(pi=points.begin(),k=0;pi!=points.end();++pi,++k)
			{
				Point q = p - *pi;
				if( lensq(q) < Point::limit) break;
			}
			if(pi!=points.end())
			{
				indicies2[i] = k;
#ifdef PRINT_INT
cerr << "existing point " << *pi << ' ' << f.dist(points[indicies2[i]]) << ' ' << f2.dist(points[indicies2[i]]) << ' ' << indicies2[i] << endl;
#endif
			}
			else
			{
				indicies2[i] = points.size();
				points.push_back(p);
#ifdef PRINT_INT
cerr << "new point " << p << ' ' << f.dist(points[indicies2[i]]) << ' ' << f2.dist(points[indicies2[i]]) << ' ' << indicies2[i] << endl;
#endif
			}
			signs2[i] = 0;
			++i;
		}	
	}
//cerr << "Before cycle ";
//for(j=0;j<size2;++j)
//{ cerr << signs2[j] << ":" << indicies2[j] << ":" << f2.dist(points[indicies2[j]]) << " "; }
//cerr << endl;
	// cycle indicies so that it starts with a zero
	while(signs2[0]!=0)
	{
		int firstsign = signs2[0];
		int firstindex = indicies2[0];
		for(j=0;j<size2-1;++j)
		{
			signs2[j] = signs2[j+1];
			indicies2[j] = indicies2[j+1];
		}
		signs2[size2-1] = firstsign;
		indicies2[size2-1] = firstindex;
	}	
	signs2[size2] = 0;
	indicies2[size2] = indicies2[0];
//cerr << "After cycle  ";
//for(j=0;j<size2+1;++j) cerr << signs2[j] << " ";
//cerr << endl;
	
	int num_zero = 0;
	for(j=0;j<size2+1;++j) if(signs2[j]==0) num_zero++;

	// The sequence now starts and ends with a zero
	// use a state machine to construct faces

#ifdef NOT_DEF
	int state = 0;
	vector<int> va;
	for(j=0;j<size2+1;++j)
	{
		switch(state)
		{
		case 0:
			if(signs2[j] != 0)
			{
				// start a new facet
				va.clear();
				va.push_back(indicies2[j-1]);
				va.push_back(indicies2[j]);
			}
			break;					
		case -1: case 1:
			if(signs2[j] == 0)
			{
				// end facet
				va.push_back(indicies2[j]);
				newfaces.push_back(MyFace(va,f));
#ifdef PRINT_INT
				cerr << "creating new face: " << va << endl;
//				for(int i=0;i<va.size();++i)
//				{
//					double dist1 = f2.dist(points[va[i]]);
//					cerr << va[i] << " " << dist1 << " ";
//				}
//				cerr << endl;
#endif
			}
			else
				va.push_back(indicies2[j]);
			break;
		}
		state = signs2[j];
	}
#endif

	int state = 0;
	vector<int> plusind,minusind;
	for(j=0;j<size2+1;++j)
	{
		switch(state) // the state of the previous vert
		{
		case 0:
			if(signs2[j] == 1)
			{
				// start a new facet
				
				plusind.push_back(indicies2[j-1]);
				plusind.push_back(indicies2[j]);
			}
			else if(signs2[j] == -1)
			{
				// start a new facet
				
				minusind.push_back(indicies2[j-1]);
				minusind.push_back(indicies2[j]);
			}
			break;					
		case -1: case 1:
			if(signs2[j] == 0)
			{
				// end facet
				if(state == 1)
					plusind.push_back(indicies2[j]);
				else
					minusind.push_back(indicies2[j]);

//				newfaces.push_back(MyFace(va,f));
#ifdef PRINT_INT
//				cerr << "creating new face: " << va << endl;
//				for(int i=0;i<va.size();++i)
//				{
//					double dist1 = f2.dist(points[va[i]]);
//					cerr << va[i] << " " << dist1 << " ";
//				}
//				cerr << endl;
#endif
			}
			else if(signs2[j] == 1)
				plusind.push_back(indicies2[j]);
			else
				minusind.push_back(indicies2[j]);
			break;
		}
		state = signs2[j];
	}
	clean_spikes(plusind);
	clean_spikes(minusind);
	if(plusind.size() > 2)
	{
		newfaces.push_back(MyFace(plusind,f));
//cerr << "creating new plus face: " << plusind << endl;
	}
	if(minusind.size() > 2)
	{
		newfaces.push_back(MyFace(minusind,f));
//cerr << "creating new minus face: " << minusind << endl;
	}
}

void intersectAllFaces()
{
	int ii,jj;
	list<MyFace>::iterator i;
	vector<MyFace>::const_iterator j;
cerr << "Number of original faces " << originalfaces.size() << endl;
	for(j=originalfaces.begin(),jj=0;j!=originalfaces.end();++j,++jj)
	{
		int size = faces.size();
cerr << "Intersecting with face " << jj << " current number of faces " << size << endl;
		for(i=faces.begin(),ii=0;ii<size;++i,++ii)
		{
			if(jj== i->originalIndex) continue;
//			if(i->good )
			{
//				cerr << "Intersecting " << ii << " " << jj << endl;
				intersectPlane(*i,*j);
			}
		}
//cerr << "Size before splice " << faces.size() << endl;
		faces.splice(faces.end(),newfaces);
//cerr << "Size after splice " << faces.size() << endl;
		faces.remove_if(no_good);
//cerr << "Size after remove " << faces.size() << endl;
	}
cerr << "Final Number of faces " << faces.size() << endl;
}

/*********************************  hiding algorithms **************************/

/**
 *  hideInteriorFaces
 *  hides those faces which do not have 'num' verticies on the sphere
 *  very rudementary routine
 */

void hideInteriorFaces(const int num)
{
	list<MyFace>::iterator i;
	for(i=faces.begin();i!=faces.end();++i)
	{
		MyFace &f = *i;
		int j,count=0;	// primary loop variable
		int size = f.indices.size();
		for(j=0;j < size;++j)
		{
			Point &p = points[f.indices[j]];
			if( lensq(p) > 1 - Point::limit ) ++count;
		}
		if(count<num) f.good = false;
	}
}

/**
 * This version checks that at least two of the points on the facet where
 * lie on one of the original edges
 */

void hideInteriorFaces2(const int num)
{
	list<MyFace>::iterator i;
	for(i=faces.begin();i!=faces.end();++i)
	{
		MyFace &f = *i;
		const MyFace &of = originalfaces[f.originalIndex];

		int j,count=0;	// primary loop variable
		int size = f.indices.size();
		for(j=0;j < size;++j)
		{
			const Point &p = points[f.indices[j]];
			int osize = of.indices.size();
			for(int k=0; k < osize; ++k)
			{
				int nextk = k +1;
				if(nextk == osize) nextk = 0;
				const Point &a = points[of.indices[k]];
				const Point &b = points[of.indices[nextk]];
				if( parallel(b-a,p-a) ) ++count;
			}
		}
		if(count<2) f.good = false;
	}
}

/**
 * This version checks that at one point lies at a vertex
 * and one lies on one of the original edges
 * if num is 5 or 6 then break on first face found
 */

void hideInteriorFaces3(const int num)
{
	list<MyFace>::iterator i;

	for(i=faces.begin();i!=faces.end();++i)
	{
		MyFace &f = *i;
		const MyFace &of = originalfaces[f.originalIndex];

		int j,count1=0,count2=0;	// primary loop variable
		int size = f.indices.size();
		for(j=0;j < size;++j)		// for each point
		{
			int pc1=0,pc2=0;
			const Point &p = points[f.indices[j]];
			int osize = of.indices.size();

			for(int k=0; k < osize; ++k)
			{
				int nextk = k +1;
				if(nextk == osize) nextk = 0;
				const Point &a = points[of.indices[k]];
				const Point &b = points[of.indices[nextk]];
				if( lensq(p-a) < Point::limit && lensq(p) > 1 - Point::limit ) ++pc1;
				if( parallel(b-a,p-a) ) 
				{
					++pc2;
//Point ba = b-a; Point pa = p-a;  ba.unit(); pa.unit(); Point n = ba^pa;
//cerr << "Parallel " << a << " b " << b << " p " << p <<endl;
//cerr << ba << " " << pa << " " << n << endl;
				}

			}
			if(pc1) ++count1;
			if(pc2) ++count2;
		}
//		cerr << "counts " << count1 << ' ' << count2 << endl;
		if(count1 < 1 || count2 <2) f.good = false;
	}

	if(num == 5 || num == 6 || num == 7 || num == 8) 
	{
		bool first = true;
		for(i=faces.begin();i!=faces.end();++i)
		{
			MyFace &f = *i;
			if(f.good)
			{
				if(first)
				{
					first = false;
					if(num == 6 || num == 8) f.normal.flip();
				}
				else
					f.good = false;
			}
		}
	}
}

/**
 *	given a face and an edge find the adjacent face which should be shown
 */

bool showAdjacentFace(MyFace &f,const Edge &e,const EdgeKey &ek,int num)
{
#ifdef PRINT_HIDE	
cerr << "Doing edge: " << ek.a << ' ' << ek.b << endl;
cerr << "First face: " << f.indices << endl;
#endif
	if(e.faces.size() < 2) return false;

	valarray<double> angles(e.faces.size());
	vector<int>::const_iterator i;
	Point firstpoint;
	const Point &a = points[ek.a];
	const Point &b = points[ek.b];
	const Point &n = f.normal;
	Point v = b-a;	// vector along edge
	v.unit();
	Point u = n^v;	// vector in face pirpendicular to edge
	u.unit();

	for(i=f.indices.begin();i!=f.indices.end();++i)
	{
#ifdef PRINT_HIDE	
//cerr << "Index " << *i << ' ' << points[*i] << endl;
#endif
		if(*i != ek.a && *i != ek.b)
		{
			firstpoint = points[*i];
			break;
		}
	}
	if(i==f.indices.end())
	{
		cerr << "Whoopse: apparently only two verticies\n";
		cerr << f.indices << endl;
		exit(1);
	}					

	Point v2 = firstpoint-a;
	if(dot(u,v2) < 0) u.flip();	// make sure u points in right direction

//cerr << "v " << v << "\nn " << n << "\nu " << u << "\nv2 " << v2 << endl;				
//Point goodnorm =  (points[0]-points[20])^(points[26]-points[20]);
//goodnorm.unit();
//cerr << "good norm " << goodnorm << endl;


	/* now have coorinate frame n,u 
		find the face with the smalest angle
	*/

	int ll,min_l,original_index;
	float min_angle = 4*M_PI;
	Point thispoint;
#ifdef PRINT_HIDE	
cerr << "faces on edges " << &e << ' ' << e.faces.size() << endl;
#endif
	vector<MyFace *>::const_iterator l;
	for(l=e.faces.begin(),ll=0;l!=e.faces.end();++l,++ll)
	{
		MyFace *fp = *l;
		if(fp == &f)
		{
			angles[ll] = 4*M_PI;
			original_index = ll;
			continue;
		}
#ifdef PRINT_HIDE	
cerr << "Adjacent face no " << ll << ' ' << fp << endl;
cerr << fp->indices << endl;
#endif
		for(i=fp->indices.begin();i!=fp->indices.end();++i)
		{
//cerr << "Index " << *i << ' ' << points[*i] << endl;
			if(*i != ek.a && *i != ek.b)
			{
				thispoint = points[*i];
				break;
			}
		}
		if(i==fp->indices.end())
		{
			cerr << "Whoopse2: apparently only two verticies\n";
			cerr << fp->indices << endl;
			exit(1);
		}					

		Point v3 = thispoint-a;
		angles[ll] = atan2(dot(v3,n),dot(v3,u));
		if(angles[ll] < 0.0 + Point::limit ) angles[ll] += 2*M_PI;
//cerr << "v3 " << v3 << " dot " << dot(v3,n) << ' ' << dot(v3,u) << " Angle " << angles[ll] << endl;
		if(angles[ll] < min_angle)
		{
			min_angle = angles[ll];
			min_l = ll;
		}
	}
	MyFace *fp = e.faces[min_l];
	if(fp->good) return false;
	fp->good = true;
	NFace &nf = *fp;
#ifdef PRINT_HIDE	
cerr << "Added face " << nf << endl;
#endif
	// maintain outward pointing normals

	if(num == 5 || num == 6 || num == 7 || num == 8)
	{
		if(min_angle > M_PI_2 - Point::limit &&
			min_angle < M_PI_2 + Point::limit )
		{
			if(dot(u,fp->normal) < 0 )
				fp->normal.flip();
		}
		else if(min_angle > M_PI + M_PI_2 - Point::limit &&
			min_angle < M_PI + M_PI_2 + Point::limit )
		{
			if(dot(u,fp->normal) > 0 )
			fp->normal.flip();
		}
		else if(min_angle < M_PI_2 || min_angle > M_PI + M_PI_2 )
		{
			if(dot(n,fp->normal) > 0 )
			fp->normal.flip();
		}
		else
		{
			if(dot(n,fp->normal) < 0 )
			fp->normal.flip();
		}
	}
//cerr << "Normal " << fp->normal << endl;
	return true;
}

/**
 *	Loop through list of edges
 *	if only one of the faces is good then find the other adjacent face
 *	which makes the smallest angle an mark that face as good
 *	repeat until this does not occur (just do it once to check)
 */

void hideInteriorFaces4(const int num)
{
	map<EdgeKey,Edge>::const_iterator k;
	bool more_to_do = true;

	while(more_to_do)
	{
		int num_faces_added = 0;
		more_to_do = false;

		for(k=edges.begin();k!=edges.end();++k)
		{
			int count=0;	// how many good faces are adjacent
			MyFace *firstface;
			vector<MyFace *> fvec = k->second.faces;
			vector<MyFace *>::const_iterator l;
			for(l=fvec.begin();l!=fvec.end();++l)
			{
				MyFace *fp = *l;
				if(fp->good)
				{
					firstface = fp;
					++count;
				}
			}
	
			if(count == 1)
			{
				more_to_do = true;
				++num_faces_added;
				valarray<double> angles(fvec.size());
				vector<int>::const_iterator i;
				Point firstpoint;
				const Point &a = points[k->first.a];
				const Point &b = points[k->first.b];
				const Point &n = firstface->normal;
				Point v = b-a;
				v.unit();
				Point u = n^v;
				u.unit();
#ifdef PRINT_HIDE	
cerr << "Doing edge: " << k->first.a << ' ' << k->first.b << ' ' << a << ',' << b << endl;
				// find a point on firstface which is not on edge
Face ff = *firstface;
cerr << "First face: " << ff << endl;
#endif
				for(i=firstface->indices.begin();i!=firstface->indices.end();++i)
				{
#ifdef PRINT_HIDE	
//cerr << "Index " << *i << ' ' << points[*i] << endl;
#endif
					if(*i != k->first.a && *i != k->first.b)
					{
						firstpoint = points[*i];
						break;
					}
				}
				if(i==firstface->indices.end())
				{
					cerr << "Whoopse: apparently only two verticies\n";
					cerr << firstface->indices << endl;
					exit(1);
				}					
				Point v2 = firstpoint-a;
				if(dot(u,v2) < 0) u.flip();
//cerr << "v " << v << "\nn " << n << "\nu " << u << "\nv2 " << v2 << endl;				
//Point goodnorm =  (points[0]-points[20])^(points[26]-points[20]);
//goodnorm.unit();
//cerr << "good norm " << goodnorm << endl;
				/* now have coorinate frame n,u */

				int ll,min_l;
				float min_angle = 4*M_PI;
				Point thispoint;
				for(l=fvec.begin(),ll=0;l!=fvec.end();++l,++ll)
				{
					MyFace *fp = *l;
#ifdef PRINT_HIDE	
cerr << "Adjacent face no " << ll;
cerr << fp->indices << endl;
#endif
					for(i=fp->indices.begin();i!=fp->indices.end();++i)
					{
cerr << "Index " << *i << ' ' << points[*i] << endl;
						if(*i != k->first.a && *i != k->first.b)
						{
							thispoint = points[*i];
							break;
						}
					}
					if(i==fp->indices.end())
					{
						cerr << "Whoopse2: apparently only two verticies\n";
						cerr << firstface->indices << endl;
						exit(1);
					}					

					Point v3 = thispoint-a;
					angles[ll] = atan2(dot(v3,n),dot(v3,u));
					if(angles[ll] < 0.0 + Point::limit ) angles[ll] += 2*M_PI;
//cerr << "v3 " << v3 << " dot " << dot(v3,n) << ' ' << dot(v3,u) << " Angle " << angles[ll] << endl;
					if(angles[ll] < min_angle)
					{
						min_angle = angles[ll];
						min_l = ll;
					}
				}
				MyFace *fp = fvec[min_l];
				fp->good = true;
				NFace &nf = *fp;
#ifdef PRINT_HIDE	
cerr << "Added face " << nf << endl;
#endif
				// maintain outward pointing normals

				if(num == 5 || num == 6)
				{
					if(min_angle > M_PI_2 - Point::limit &&
						min_angle < M_PI_2 + Point::limit )
					{
						if(dot(u,fp->normal) < 0 )
							fp->normal.flip();
					}
					else if(min_angle > M_PI + M_PI_2 - Point::limit &&
						min_angle < M_PI + M_PI_2 + Point::limit )
					{
						if(dot(u,fp->normal) > 0 )
							fp->normal.flip();
					}
					else if(min_angle < M_PI_2 || min_angle > M_PI + M_PI_2 )
					{
						if(dot(n,fp->normal) > 0 )
							fp->normal.flip();
					}
					else
					{
						if(dot(n,fp->normal) < 0 )
							fp->normal.flip();
					}

				}
//cerr << "Normal " << fp->normal << endl;
			}
			 // end if statement
			
		} // end loop through edges
		cerr << "Added " << num_faces_added << " faces adjacent to edges\n";
//		break;
	} // end while
}

/**
 * 	loop through faces
 *		adding adjacents where necessary
 */

void hideInteriorFaces7(const int num)
{
	list<MyFace>::iterator k;
	bool more_to_do = true;


	// clear the linking faces

	for(k=faces.begin();k!=faces.end();++k)
	{
	    	MyFace &fp = *k;
		int i;
		int size = fp.indices.size();
		for(i=0;i<size;++i)
			fp.linkingFaces.push_back(-1);
	}

	while(more_to_do)
	{
		int num_faces_added = 0;
		more_to_do = false;

		for(k=faces.begin();k!=faces.end();++k)
		{
		    	MyFace &fp = *k;
			if( fp.good )
			{
				int size = fp.indices.size();
				for(int i=0;i<size;++i)
				{
					if(fp.linkingFaces[i] >= 0) continue;

					int nexti = i+1; if(nexti==size) nexti = 0;
					EdgeKey ek(fp.indices[i],fp.indices[nexti]);
					Edge e = edges[ek];
					if(showAdjacentFace(fp,e,ek,num))
					{
						more_to_do = true;
						++num_faces_added;
					}
				}
			}
		}
		cerr << "Added " << num_faces_added << " faces adjacent to edges\n";
//		break;
	} // end while
}

/**************************** merging routines *****************************/

void mergeFaces(MyFace &f1,MyFace &f2,const int a,const int b)
{
#ifdef PRINT_MERGE
cerr << "Merge faces on " << a << ' ' << b << ' ';
cerr << "(" << f1.indices << ") ";
cerr << " (" << f2.indices << ")\n";
#endif
	vector<int> newind = merge_at(f1.indices,f2.indices,a,b);
#ifdef PRINT_MERGE
cerr << "before spike -> (" << newind << ")\n";
#endif
	clean_spikes(newind);
	if(newind.size() <= 2 ) return;

	faces.push_back(MyFace(newind,f1));
	MyFace &newface = faces.back();
	f1.good = false;
	f2.good = false;
	newface.good = true;
#ifdef PRINT_MERGE
cerr << "Mergeinds -> (" << newface.indices << ")\n";
#endif

	/* now need to adjust each of the edges to point to this new face */

	int i,size=newind.size();
	for(i=0;i<size;++i)
	{
		int nexti = i+1;
		if(nexti == size) nexti = 0;

		Edge &e = edges[EdgeKey(newind[i],newind[nexti])];
		vector<MyFace *>::iterator j;
		int jj;
		for(j=e.faces.begin(),jj=0;j!=e.faces.end();++j,++jj)
		{
			MyFace *fp1 = &f1;
			MyFace *fp2 = &f2;
			MyFace *fp3 = &newface;
			if(*j == fp1 || *j == fp2)
			{
//cerr << "Swapping face on " << newind[i] << ' ' << newind[nexti] << " (" << (*j)->indices << ") with (" << newface.indices << ")\n";
				// e.faces[jj] = fp3;
				*j = fp3;
			}
		}
	}
}

/**
 * Link faces
 * links together adjacent faces which share a common edge
 */

void linkFaces()
{
	map<EdgeKey,Edge>::iterator k;
	bool more_to_do = true;

		/* work out the good edges */

		for(k=edges.begin();k!=edges.end();++k)
		{

			vector<MyFace *> fvec = k->second.faces;
			vector<MyFace *>::const_iterator l,m;
			int count = 0;
			for(l=fvec.begin();l!=fvec.end();++l)
			{
				MyFace *fp1 = *l;
				if(fp1->good) ++count;
			}
			if(count==2) k->second.good = true;
			else	     k->second.good = false;
		}

	while(more_to_do)
	{
		int num_faces_merged = 0;
		more_to_do = false;


		for(k=edges.begin();k!=edges.end();++k)
		{
			if( ! k->second.good ) continue;
			int count=0;	// how many good faces are adjacent
			vector<MyFace *> fvec = k->second.faces;
			vector<MyFace *>::const_iterator l,m;
			for(l=fvec.begin();l!=fvec.end();++l)
			{
				MyFace *fp1 = *l;

				if(fp1->good)
				    for(m=l+1;m!=fvec.end();++m)
				    {
					MyFace *fp2 = *m;
					if(fp1 == fp2)
					{
						k->second.good = false;
						continue;
					}
					if(fp2->good && fp1->originalIndex == fp2->originalIndex)
					{
						mergeFaces(*fp1,*fp2,k->first.a,k->first.b);
						more_to_do = true;
						k->second.good = false;
						++num_faces_merged;
					}
				   }
			}
		}
		cerr << "Merged " << num_faces_merged << endl;
	}
}


void buildEdges()
{
	list<MyFace>::iterator i;
	int ii;
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{

		MyFace &f = *i;
		f.setFindex(ii);
		int j,nextj;	// primary loop variable
		int size = f.indices.size();
		for(j=0;j < size;++j)
		{
			nextj = j+1;
			if(nextj == size) nextj = 0;
			EdgeKey ek(f.indices[j],f.indices[nextj]);
			edges[ek].addFace(&f);
		}
	}
	map<EdgeKey,Edge>::const_iterator k;
	for(k=edges.begin();k!=edges.end();++k)
	{
//		cerr << k->first.a << ',' << k->first.b << ':';
		vector<MyFace *> fvec = k->second.faces;
		vector<MyFace *>::const_iterator l;
		for(l=fvec.begin();l!=fvec.end();++l)
		{
			MyFace *fp = *l;
//			cerr << ' ' << fp->Findex;
		}
//		cerr << endl;
	}
}

void colourFaces()
{
	list<MyFace>::iterator i;
	Point v1(1,1,1);
	Point vx(1,0,0);
	Point vy(0,1,0);
	Point vz(0,0,1);
	
	for(i=faces.begin();i!=faces.end();++i)
	{
		MyFace &f = *i;
		Point n = f.normal;
		if(dot(n,v1) < 0 ) n.flip();
		int r = int( 128 + dot(n,vx) * 127 );
		int g = int( 128 + dot(n,vy) * 127 );
		int b = int( 128 + dot(n,vz) * 127 );
//cerr << "Colour " << r << " " << g << " " << b << " normal " << n << endl;
		f.setColour(r,g,b);
	}
}

/*
 * colour the first num faces one colour, etc
 */

void colourFaces2(int num)
{
	list<MyFace>::iterator i;
	int ii,r,g,b;
	Point v1(1,1,1);
	Point vx(1,0,0);
	Point vy(0,1,0);
	Point vz(0,0,1);
	
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{
		MyFace &f = *i;
		Point n = f.normal;
		if(dot(n,v1) < 0 ) n.flip();
		if(ii==num) ii = 0;
		if(ii==0)
		{
			r = int( 128 + dot(n,vx) * 127 );
			g = int( 128 + dot(n,vy) * 127 );
			b = int( 128 + dot(n,vz) * 127 );
cerr << "Colour " << r << " " << g << " " << b << " normal " << n << endl;
		}
		f.setColour(r,g,b);
	}
}

/**
 * colour the faces by the type of polyhedra
 *	tri = red
 *	square = green
 *	pent = blue
 *	hex = yellow
 *	oct = purple
 *	dec = cyan
 *	other = white
 */

void colourFaces3(int num)
{
	list<MyFace>::iterator i;
	int ii,r,g,b;
	Point v1(1,1,1);
	Point vx(1,0,0);
	Point vy(0,1,0);
	Point vz(0,0,1);
	
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{
		MyFace &f = *i;
		switch(f.indices.size())
		{
		case 3: r = 255; g = 0; b = 0; break;
		case 4: r = 0; g = 255; b = 0; break;
		case 5: r = 0; g = 0; b = 255; break;
		case 6: r = 0; g = 255; b = 255; break;
		case 8: r = 255; g = 0; b = 255; break;
		case 10: r = 255; g = 255; b = 0; break;
		default: r = 255; g = 255; b = 255; break;
		}
		f.setColour(r,g,b);
	}
}

void setOriginalIndicies()
{
	list<MyFace>::iterator i;
	int ii;
	
	for(i=faces.begin(),ii=0;i!=faces.end();++i,++ii)
	{
		i->originalIndex = ii;
	}
}

void print_help()
{
	cout << "convert opts < stdin > stdout\n";
	cout << "Where the input is that produced by kaleido's -w option\n";
	cout << "and the output is in JVX format\n";
	cout << "\t-i\tIntersect all the faces with each other\n";
	cout << "\t-j\tJoin star faces together\n";
	cout << "\t-c\tColour the faces (uses direction of normal)\n";
	cout << "\t-h num\tRemove interior facets\n";
	cout << "\t-h 1\tRemoves those not adjacent to a vertex\n";
	cout << "\t-h 2\tRemoves those not adjacent to two verticies\n";
	cout << "\t-h 3\tRemoves those where two of verticies do not line on an edge of the polyhedron\n";
	cout << "\t-m file.WRL\tMatches the vertices to those in file\n";
	cout << "\t\tUseful for making coumpounds\n";
	cout << "\t-?\tPrints this screen\n";
	cout << "\t--help\tPrints this screen\n";
	exit(0);
}

main(int argc,char **argv)
{
	bool intersectFlag = false,colourFlag=false,linkFlag=false;
	int hideFlag = -1,colourNum = 0,joinFlag = 0;

	for(int i=1;i<argc;++i)
	{
		if(!strcmp(argv[i],"-?")) print_help();
		if(!strcmp(argv[i],"--help")) print_help();
		if(!strcmp(argv[i],"-i")) intersectFlag = true;
		if(!strcmp(argv[i],"-j")) joinFlag = atoi(argv[++i]);
		if(!strcmp(argv[i],"-c")) colourFlag = true;
		if(!strcmp(argv[i],"-C")) colourNum = atoi(argv[++i]);
		if(!strcmp(argv[i],"-h")) hideFlag = atoi(argv[++i]);
		if(!strcmp(argv[i],"-l")) linkFlag = true;
	}
	readPoints();
	readFaces(joinFlag);
//	setOriginalIndicies();
	if(colourFlag)
		colourFaces();
	if(colourNum > 0)
		colourFaces2(colourNum);
	if(colourNum == -1)
		colourFaces3(colourNum);
	if(intersectFlag)
		intersectAllFaces();


	switch(hideFlag)
	{
	case 1: case 2:
		hideInteriorFaces(hideFlag);
		break;
	case 3:
		hideInteriorFaces2(hideFlag);
		break;
	case 4:
		buildEdges();
		hideInteriorFaces3(hideFlag);
		hideInteriorFaces4(hideFlag);
		break;
	case 5: case 6:
		buildEdges();
		hideInteriorFaces3(hideFlag);
		hideInteriorFaces4(hideFlag);
		break;
	case 7: case 8:
		buildEdges();
		hideInteriorFaces3(hideFlag);
		hideInteriorFaces7(hideFlag);
		break;
	}
	if(linkFlag)
	{
		if(hideFlag < 4 || hideFlag > 8 )
			buildEdges();
		linkFaces();
	}
	print(colourFlag || (colourNum));
}