\datethis % print date on listing

@* qhConvert.  Covert output of qhull to a form readable by MorseExtract.
Note it does not always work correctly.
In particular if vertices on the boundary span a simplex not in
the boundary it will not detect this and will cone the simplex to infinity.
So to make it work correctly you should enclose the region
in a tetrahedron.
Consequently, the last five vertices will not be in K, 
so on output they are preceded with an x.


@c
@<Header files to include@>@/
@<Data structures@>@/
@<Subroutine prototypes@>@/
@<Global variables@>@/
@<Subroutines@>@/
@<The main routine@>@/

@ @<The main routine@>=

//int mymain 
int main
(int argc, const char * argv[]) {
int i,j,k;
struct vertex *v;
struct edge *e;
struct triangle *f;
struct tetrahedron *p;

ne=nf=nbf=0;

if(argc<2)
{
	abort_message("Usage: xqhconvert vertex-file < qdelaunay output");
}

@<Set up vertices@>@;
@<Read in tetrahedra and determine edges and faces@>@;
@<Cone boundary to infinity@>@;
@<Output vertices@>@;
@<Output edges@>@;
@<Output faces@>@;
@<Output tetrahedra@>@;
return 0;
}

@ @<Global variables@>=
struct vertex *vertices;
struct tetrahedron *tetrahedra;
int nv,ne,nf,nt;  // number of vertices, edges, faces and tetrahedra
int nbv,nbf;  // number of boundary vertices, faces
FILE *df;
struct edge *laste,*firste;
struct triangle *lastf,*firstf;

@ @<Subroutine prototypes@>=

@ @<Header files to include@>=
#include <stdio.h>

@ Here are data structures

@<Data structures@>=

struct vertex
{
	struct edge *e;  // first on list of edges in star, to infinity if possible
	struct triangle *t;  // first on list of triangles in star
	int h; // value
	int flag;
};

struct edge
{
	int v[2];  // vertices always in order
	struct triangle *t; // a triangle in star, to infinity if possible
	struct edge *next[2];
	struct edge *nextid;  // edge with next id
	int id;
};

struct triangle
{
	int v[3];  // vertices in order
	struct edge *e[3];  // the edges
	int links[2];  // indices of the two cofaces
	struct triangle *next[3];
	struct triangle *nextid;
	int id;
};

struct tetrahedron
{
	int v[4]; // vertices in order
	struct triangle *f[4]; // faces
};


@ @<Set up vertices@>=
{
	int val;
	
	df = fopen(argv[1],"r");

	i = fscanf(df,"%d\n",&nv);
	if(i!=1) abort_message("bad number of vertices");
	nv++;  // make room for vertex at infinity
	v = vertices = (struct vertex*)malloc(nv*sizeof(struct vertex));
	for(i=0;i<nv;i++)
	{
		v->flag = 0;
		v->e = NULL;
		v->t = NULL;
		
		if(i<nv-1)
		{
			fscanf(df,"%d\n",&val);
			v->h = val; 
		}
		else v->h=0; 
		v++;
	}
	fclose(df);
	df = stdin;

	@<flag boundary vertices@>@;
}

@ @<flag boundary vertices@>=
{
	i = fscanf(df,"%d\n",&nbv);
	if(i!=1) abort_message("bad number of boundary vertices");
	if(nbv!=4) abort_message("there must exactly 4 boundary vertices");
	
	for(j=0;j<nbv;j++)
	{
		i = fscanf(df,"%d\n",&k);
		if(i!=1) abort_message("bad boundary vertex");
		if(k< 0 || k>=nv) abort_message("bad boundary vertex id");
		vertices[k].flag = 1;
	}
}

@ @<Read in tetrahedra and determine edges and faces@>=
{
	i = fscanf(df,"%d\n",&nt);
	if(i!=1) abort_message("bad number of tetrahedra");
	p = tetrahedra = (struct tetrahedron*)malloc(nt*sizeof(struct tetrahedron));
	
	for(j=0;j<nt;j++,p++)
	{
		@<Read in a tetrahedron and determine edges and faces@>@;
	}
}




@ @<Read in a tetrahedron and determine edges and faces@>=
{
	int ids[4];
	
	i = fscanf(df,"%d %d %d %d\n",ids,ids+1,ids+2,ids+3);
	if(i!=4) abort_message("bad tetrahedron vertex");
	
	@<sort tetra ids@>@;
	for(k=0;k<4;k++)
	{
		@<determine k-th face@>@;
	}
}


@ @<sort tetra ids@>=
{
	int ii,temp;
	
	for(k=0;k<3;k++)
	{
		for(ii=k+1;ii<4;ii++)
		{
			if(ids[ii]<ids[k])
			{
				temp = ids[k];
				ids[k] = ids[ii];
				ids[ii] = temp;
			}
		}
	}
	for(k=0;k<4;k++)
	{
		p->v[k]=ids[k];
	}
}

@ @<determine k-th face@>=
{
	int ii,jj;
	
	for(ii=0,jj=0;ii<4;ii++)
	{
		if(ii!=k)
		{
			ids[jj] = p->v[ii];
			jj++;
		}
	}
	
	p->f[k] = findaddface(ids,j);
}

@ @<Subroutines@>=

void abort_message(char *s)
{
	fprintf(stderr,"Error: %s",s);
	exit(1);
}


@ The following sees if the edge with vertex ids |vids| is already defined
and returns its address if it is.
But if it is not, it creates it.

@<Subroutines@>=
struct edge *findaddedge(int vids[2],struct triangle *t)
{
	struct edge *e;
	int i;
	
	for(e = vertices[vids[0]].e; e!=NULL; e = e->next[1-i])
	{
		for(i=0;i<2;i++)
		{
			if(e->v[i]!=vids[i]) break;
		}
		if(i==2) 
		{
			e->t = t;  // so boundary edges point to infinity
			return e;
		}
	}
	@<make edge@>@;
	return e;
}


@ @<make edge@>=
{
	e = (struct edge *)malloc(sizeof(struct edge));
	for(i=0;i<2;i++)
	{
		e->v[i] = vids[i];
		e->next[i] = vertices[vids[i]].e;
		vertices[vids[i]].e = e;
	}
	e->t=t;
	e->id = ne++;
	if(e->id)
	{
		laste->nextid = e;
	}
	else firste = e;
	laste = e;
	e->nextid = NULL;
}

@ The following sees if the face with vertex ids |vids| is already defined
and returns its address if it is.
But if it is not, it creates it.

@<Subroutines@>=
struct triangle *findaddface(int vids[3],int tid)
{
	struct triangle *t;
	int i;
	
	for(t = vertices[vids[0]].t; t!=NULL; t = t->next[i])
	{
		for(i=0;i<3;i++)
		{
			if(t->v[i]!=vids[i]) break;
		}
		if(i==3) 
		{
			t->links[1]=tid;
			return t;
		}
		for(i=0;i<3;i++)
		{
			if(t->v[i]==vids[0]) break;
		}
	}
	@<make triangle and determine edges@>@;
	return t;
}

@ @<make triangle and determine edges@>=
{
	int ids[2],j,k;
	
	t = (struct triangle *)malloc(sizeof(struct triangle));
	for(i=0;i<3;i++)
	{
		t->v[i] = vids[i];
		t->next[i] = vertices[vids[i]].t;
		vertices[vids[i]].t = t;
		for(j=0,k=0;j<3;j++)
		{
			if(j!=i)
			{
				ids[k]=vids[j];
				k++;
			}
		}
		t->e[i] = findaddedge(ids,t);
	}
	t->links[0]=tid;
	t->links[1]=-1;
	t->id = nf++;
	if(t->id)
	{
		lastf->nextid = t;
	}
	else firstf = t;
	lastf = t;
	t->nextid = NULL;
}


@ @<Cone boundary to infinity@>=
{
	int bids[3];
	
	bids[2]=nv-1;
	
	for(i=0;i<nv-1;i++)
	{
		if(vertices[i].flag)
		{
			for(f = vertices[i].t; f!=NULL; f = f->next[j])
			{
				if(f->v[0]==i)
				{
					for(j=1;j<3;j++)
					{
						if(vertices[f->v[j]].flag==0) break;
					}
					if(j==3) 
					{
						f->links[1] = nt+nbf;
						bids[0]=f->v[0];
						bids[1]=f->v[1];
						findaddface(bids,nt+nbf);
						bids[1]=f->v[2];
						findaddface(bids,nt+nbf);
						bids[0]=f->v[1];
						findaddface(bids,nt+nbf);
						nbf++;
					}
					j=0;
				}
				else if (f->v[1]==i) j=1;
				else if (f->v[2]==i) j=2;
				else abort_message("error 1");
			}
		}
	}
}











@ @<Output vertices@>=
{
	printf("%d\n",nv);
	v = vertices;
	for(i=0;i<nv;i++) 
	{
		if(i>=nv-5) printf("x");
		printf("%d\n",v->h);
		v++;
	}
}

@ @<Output edges@>=
{
	printf("%d\n",ne);
	v = vertices;
	for(i=0;i<nv-1;i++) 
	{
		for(e = v->e; e!=NULL; e = e->next[j])
		{
			if(e->v[0]==i)
			{
				//printf("%d %d %d\n",e->id,e->v[0],e->v[1]);
				j=0;
			}
			else j=1;
		}
		v++;
	}
	for(e=firste,i=0;e!=NULL;e=e->nextid,i++)
	{
		printf("%d %d\n",e->v[0],e->v[1]);
		if(e->id!=i) abort_message("edges out of order");
	}
}

@ @<Output faces@>=
{
	printf("%d\n",nf);
	v = vertices;
	for(i=0;i<nv-2;i++) 
	{
		for(f = v->t; f!=NULL; f = f->next[j])
		{
			if(f->v[0]==i)
			{
				//printf("%d %d %d %d\n",f->id,f->e[0]->id,f->e[1]->id,f->e[2]->id);
				j=0;
			}
			else if (f->v[1]==i) j=1;
			else if (f->v[2]==i) j=2;
			else abort_message("error 1");
		}
		v++;
	}
	for(f=firstf,i=0;f!=NULL;f=f->nextid,i++)
	{
		printf("%d %d %d\n",f->e[0]->id,f->e[1]->id,f->e[2]->id);
		if(f->id!=i) abort_message("faces out of order");
	}

}

@ @<Output tetrahedra@>=
{
	struct tetrahedron *t;
	
	printf("%d\n",nt+nbf);
	for(i=0;i<nt;i++) 
	{
		t = tetrahedra+i;
		printf("%d %d %d %d\n",t->f[0]->id,t->f[1]->id,t->f[2]->id,t->f[3]->id);
	}
	v = vertices;
	for(i=0;i<nv-2;i++) 
	{
		if(v->flag)
		{
			for(f = v->t; f!=NULL; f = f->next[j])
			{
				if(f->v[0]==i)
				{
					if(vertices[f->v[1]].flag && vertices[f->v[2]].flag)
						printf("%d %d %d %d\n",f->e[0]->t->id,f->e[1]->t->id,f->e[2]->t->id,f->id);
					j=0;
				}
				else if (f->v[1]==i) j=1;
				else if (f->v[2]==i) j=2;
				else abort_message("error 3");
			}
		}
		v++;
	}
}