/*
 *
 * Description
 *	Pack an OFF object in memory
 *	``Packing'' means to rearrange the data structure so common vertices
 *	are shared.  Since in typical objects every vertex
 *	is shared by three or more polygons, this can often result in a
 *	substantial reduction in storage and increase in rendering
 *	efficiency (because shared vertices need only be transformed once).
 *
 * Input
 *	Obj		Pointer to object structure in which data is stored.
 * Output
 *	obj		Same as input pointer, but underlying object structure
 *			is modified.
 *
 * Diagnostics
 *	None.  (If unsuccessful, the object is returned unchanged.)
 *
 * Author
 *	Allen Akin
 *	Digital Equipment Corp.
 *	Workstation Systems Engineering
 *	Palo Alto, CA
 *
 * History
 *	20-Feb-87	Created
 *	11-Oct-89	Ability to pack normals was removed when conventions
 *			for dealing with normals were changed.
 *
 * Notes
 *	In order to do this properly, one should really treat all of the
 * 	attributes at a vertex simultaneously - a vertex is redundant only
 *	if all of the attributes at that vertex are redundant (vertex
 *	coordinate, color, vertex normal, etc.).  This routine will just
 *	pack the geometry property.  It is currently up to the application to
 *	regenerate per-vertex colors or normals.
 */

#include <stdio.h>
#include "off.h"



typedef struct {float x, y, z;} OFFxyz;
static qsort();
static turtlesort();
static int precedes();
static int follows();
static int equals();



OFFPackObj(obj)
OFFObjDesc *obj;
	{
	register OFFProperty *prop;

	for (prop = obj->FirstProp; prop; prop = prop->NextProp)
		{
		if (strcmp(prop->PropName, "geometry") == 0
		    && strcmp(prop->DataFormat, "fff") == 0)
			OFFPackProperty(prop);
		}
	}



OFFPackProperty(prop)
OFFProperty *prop;
	{
	register int newI;
	register int *vlist;
	register OFFxyz *newVertices;
	register int *ip;
	register int *indexMap;
	register int oldI;
	int nVertices;
	int nCounts;
	int nIndices;
	OFFxyz *oldVertices;
	short *oldCounts;
	short *oldIndices;
	char *newData;
	OFFxyz *xyzp;
	short *sp;
	int vertex;
	char *malloc();


	ip = (int *) prop->PropData;
	nVertices = *ip++;
	nCounts = *ip++;
	nIndices = *ip++;
	oldVertices = (OFFxyz *) ip;
	oldCounts = (short *) (oldVertices + nVertices);
	oldIndices = oldCounts + nCounts;


	newVertices = NULL;
	vlist = indexMap = NULL;
	if (!(newVertices = (OFFxyz *) malloc(nVertices * sizeof(OFFxyz)))
	 || !(vlist = (int *) malloc(nVertices * sizeof(int)))
	 || !(indexMap = (int *) malloc(nVertices * sizeof(int))) )
		{
		fprintf(stderr, "OFFPackProperty:  insufficient memory\n");
		if (newVertices)
			free((char *) newVertices);
		if (vlist)
			free((char *) vlist);
		if (indexMap)
			free((char *) indexMap);
		return;
		}


	/* Set up indices of original vertices: */

	for (newI = 0; newI < nVertices; ++newI)
		vlist[newI] = newI;


	/* Sort original vertices by swapping indices in vlist: */

	qsort(oldVertices, vlist, vlist + (nVertices - 1));


	/*
	 * The following code does two things:
	 *	(1) Copies vertices from oldVertices to newVertices,
	 *	    throwing away duplicates.
	 *	(2) Builds a mapping from old vertex indices to new
	 *	    vertex indices.
	 * Because vlist gives us ordered access to oldVertices, it's
	 * efficient to do these operations in parallel.
	 */

	newI = -1;
	for (vertex = 0, oldI = vlist[vertex]; vertex < nVertices; )
		{
		/*
		 * At this point, we have a new unique vertex.  Store it
		 * in the next position in the newVertices array, and
		 * update the indexMap to show how the old vertex index
		 * maps into the new one.
		 */

		newVertices[++newI] = oldVertices[oldI];
		indexMap[oldI] = newI;


		/*
		 * Now scan through succeeding oldVertices (continuing in
		 * sorted order).  If the old vertex equals the current new
		 * vertex, discard it and make the old vertex index map to
		 * the new vertex index.  If the old vertex doesn't equal
		 * the current new vertex, go back to the outer loop and
		 * store the new vertex.
		 */

		for (++vertex; vertex < nVertices; ++vertex)
			{
			oldI = vlist[vertex];
			if (!equals(oldVertices + oldI, newVertices + newI))
				break;
			indexMap[oldI] = newI;
			}
		}


	/*
	 * newI contains the index of the last-used element of newVertices.
	 * Use this information plus the old data sizes to compute the size
	 * of the new data block, and allocate it.
	 */
	
	++newI;
	if (!(newData = malloc(
	    3 * sizeof(int)
	    + newI * sizeof(OFFxyz)
	    + nCounts * sizeof(short)
	    + nIndices * sizeof(short)  )))
		{
		free((char *) newVertices);
		free((char *) vlist);
		free((char *) indexMap);
		fprintf(stderr, "OFFPackProperty:  insufficient memory\n");
		return;
		}


	/*
	 * Stuff the new values into the new data record, mapping old
	 * vertex indices into new ones as we go.
	 */

	ip = (int *) newData;
	*ip++ = newI;
	*ip++ = nCounts;
	*ip++ = nIndices;
	bcopy((char *) newVertices, (char *) ip, newI * sizeof(OFFxyz));
	xyzp = (OFFxyz *) ip;
	xyzp += newI;
	sp = (short *) xyzp;
	bcopy((char *) oldCounts, (char *) sp, nCounts * sizeof(int));
	sp += nCounts;
	while (--nIndices >= 0)
		*sp++ = indexMap[*oldIndices++ - 1] + 1;


	/*
	 * Throw away the old and auxiliary data structures, and store the
	 * new ones.
	 */

	free(prop->PropData);
	prop->PropData = newData;
	free((char *) newVertices);
	free((char *) vlist);
	free((char *) indexMap);
	}



#define PRECEDES(a,b) precedes(&(a), &(b))
#define FOLLOWS(a,b) follows(&(a), &(b))
#define QTHRESH 10	/* threshold for linear sort versus quicksort */


static turtlesort(v, first, last)
OFFxyz *v;
register int *first;
register int *last;
	{
	int temp;
	OFFxyz min_key;
	OFFxyz cur_key;
	register int *current;

	for (; first < last; ++first)
		{
		min_key = v[*first];
		for (current = first + 1; current <= last; ++current)
			{
			cur_key = v[*current];
			if (PRECEDES(cur_key, min_key))
				{
				min_key = cur_key;
				temp = *first;
				*first = *current;
				*current = temp;
				}
			}
		}
	}


static qsort(v, first, last)
OFFxyz *v;
register int *first;
register int *last;
	{
	OFFxyz pivot_key;
	int *pivot;
	register int *high;
	register int *low;
	int temp;

reenter:
	if (first + QTHRESH >= last)
		{
		turtlesort(v, first, last);
		return;
		}
	

	/* Find median of first, middle, and last keys: */

	{
	int *mid;
	OFFxyz f, m, l;

	mid = first + ((last - first) >> 1);

	f = v[*first];
	m = v[*mid];
	l = v[*last];

	if (PRECEDES(f, l))
		{
		if (PRECEDES(f, m))
			pivot = (PRECEDES(m, l))? mid: last;
		else
			pivot = first;
		}
	else
		{
		if (FOLLOWS(f, m))
			pivot = (PRECEDES(m, l))? last: mid;
		else
			pivot = first;
		}
	}


	/* Swap the pivot element into sentinel position at the beginning: */

	pivot_key = v[*pivot];

	temp = *first;
	*first = *pivot;
	*pivot = temp;


	/* Partition the array around the pivot: */

	low = first;
	high = last;
	while (low < high)
		{
		while (FOLLOWS(v[*high], pivot_key))
			--high;
		while (low < high && !(FOLLOWS(v[*low], pivot_key)))
			++low;
		if (low < high)
			{
			temp = *low;
			*low = *high;
			*high = temp;
			}
		}


	/* Swap the pivot into its proper place in the middle: */

	temp = *first;
	*first = *high;
	*high = temp;


	/*
	 * Recurse to handle subarrays.  Use direct recursion on shortest
	 * subarray, then unwind tail recursion to handle second subarray.
	 * (This minimizes worst-case stack requirements.)  (Note:  the
	 * casts below eliminate two divisions.)
	 */

	if ((char *)high - (char *)first < (char *)last - (char *)high)
		{
		qsort(v, first, high - 1);
		first = high + 1;
		goto reenter;
		}
	else
		{
		qsort(v, high + 1, last);
		last = high - 1;
		goto reenter;
		}
	}



static int precedes(a, b)
register OFFxyz *a;
register OFFxyz *b;
	{
	if (a->x < b->x)
		return 1;
	if (a->x > b->x)
		return 0;
	if (a->y < b->y)
		return 1;
	if (a->y > b->y)
		return 0;
	return a->z < b->z;
	}



static int follows(a, b)
register OFFxyz *a;
register OFFxyz *b;
	{
	if (a->x > b->x)
		return 1;
	if (a->x < b->x)
		return 0;
	if (a->y > b->y)
		return 1;
	if (a->y < b->y)
		return 0;
	return a->z > b->z;
	}



static int equals(a, b)
register OFFxyz *a;
register OFFxyz *b;
	{
	return a->x == b->x && a->y == b->y && a->z == b->z;
	}