/*ScianPictures.c
  Eric Pepke
  June 25, 1990
  Code to draw 3-D pictures in scian.
*/

#include "Scian.h"
#include "ScianTypes.h"
#include "ScianColors.h"
#include "ScianPictures.h"
#include "ScianIDs.h"
#include "ScianErrors.h"
#include "ScianWindows.h"
#include "ScianDatasets.h"
#include "ScianArrays.h"
#include "ScianPreferences.h"
#include "ScianSpaces.h"
#include "ScianVisObjects.h"
#include "ScianPick.h"

static void AppendItemToPicture(PicPtr, PicItemPtr);
static PolyPtr ConvertPolyOntoPolys(PolysPtr polys, PolyPtr poly);

#define MAXCYLSIDES	512
#define CYLSIDES	16	/*Sides around cylinders, must be 2^n, 4 or more*/
#define MAXSUB		8	/*Maximum number of subdivisions*/
#ifndef SPHERESUB
#define SPHERESUB	2	/*Number of subdivisions of a sphere*/
#endif

#ifdef GRAPHICS
Linestyle dashedLine = 0xF0F0;
Linestyle dottedLine = 0xAAAA;
#endif

double cefs[MAXSUB];		/*Chord extension factor for each subdivision*/
double cefs_2[MAXSUB];		/*cefs / 2*/

ObjPtr picClass = 0;		/*Class for all pictures*/

Bool orderReversed = false;	/*True iff polygon order is reversed*/
Bool overrideColor = false;	/*True iff override color*/
int curColorMode;		/*The current color mode, CMODECMAP or CMODERGB*/
int curColorShading;		/*The current color shading*/
int curLightShading;		/*The current light shading*/ 
Bool curIsTransparent;		/*True iff transparent*/

long pictureTime = 0;		/*Picture drawing time*/

VertexPtr freedVertices = 0;	/*Freed vertices for the picking*/

#ifdef GRAPHICS
unsigned short greyPattern[16] = 
	{
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA
	};

unsigned short altGreyPattern[16] = 
	{
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555,
	    0xAAAA,
	    0x5555
	};

unsigned short funkyGreyPattern[16] = 
	{
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA
	};

unsigned short funkyAltGreyPattern[16] = 
	{
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555,
	    0xAAAA,
	    0xAAAA,
	    0x5555,
	    0x5555
	};

unsigned short oddPattern[16] = 
	{
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000
	};

unsigned short evenPattern[16] = 
	{
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF,
	    0x0000,
	    0xFFFF
	};

unsigned short tempPattern[16];

Matrix cardBasis=
    {
        {-0.5,  1.5, -1.5,  0.5},
        { 1.0, -2.5,  2.0, -0.5},
        {-0.5,  0.0,  0.5,  0.0},
        { 0.0,  1.0,  0.0,  0.0}
    };
#endif

typedef struct
    {
	Bool rWrite;
	Bool gWrite;
	Bool bWrite;
	Bool zWrite;
    } WriteMask;

#define NWRITEMASKS	10
static int curWriteMask = 0;

static WriteMask writeMasks[NWRITEMASKS];

#ifdef HAVE_PROTOTYPES
static void SetCurMask(void)
#else
static void SetCurMask()
#endif
/*Sets the mask according to the current mask*/
{
    if (rgbp)
    {
	RGBwritemask(writeMasks[curWriteMask] . rWrite ? 0xff : 0,
		     writeMasks[curWriteMask] . gWrite ? 0xff : 0,
		     writeMasks[curWriteMask] . bWrite ? 0xff : 0);
    }
    zwritemask(writeMasks[curWriteMask] . zWrite ? 0xffffffff : 0);
}

#ifdef HAVE_PROTOTYPES
void BeginMask(Bool rWrite, Bool gWrite, Bool bWrite, Bool zWrite)
#else
void BeginMask(rWrite, gWrite, bWrite, zWrite)
#endif
/*Sets up write masks to r, g, b, and z buffer.  Bracket with EndMask*/
{
    ++curWriteMask;

    writeMasks[curWriteMask] . rWrite = rWrite ? writeMasks[curWriteMask - 1] . rWrite : false;
    writeMasks[curWriteMask] . gWrite = gWrite ? writeMasks[curWriteMask - 1] . gWrite : false;
    writeMasks[curWriteMask] . bWrite = bWrite ? writeMasks[curWriteMask - 1] . bWrite : false;
    writeMasks[curWriteMask] . zWrite = zWrite ? writeMasks[curWriteMask - 1] . zWrite : false;

    SetCurMask();
}

#ifdef HAVE_PROTOTYPES
void EndMask(void)
#else
void EndMask()
#endif
/*End of a BeginMask block*/
{
    --curWriteMask;

    SetCurMask();
}

#ifdef HAVE_PROTOTYPES
void OverrideColor(Bool whether)
#else
void OverrideColor(whether)
Bool whether;
#endif
/*Overrides the colors in a picture*/
{
    overrideColor = whether;
}

static ObjPtr MakePictureDeformed(picture)
ObjPtr picture;
/*Makes a picture's PICDEFORMED variable.  Really doesn't have to do anything*/
{
    SetVar(picture, PICDEFORMED, ObjTrue);
}

#ifdef HAVE_PROTOTYPES
VertexPtr NewVertex(ObjPtr picture, int flags)
#else
VertexPtr NewVertex(picture, flags)
ObjPtr picture;
int flags;
#endif
/*Returns a new vertex, stuck on picture*/
{
    VertexPtr retVal;

    if (freedVertices)
    {
	retVal = freedVertices;
	freedVertices = freedVertices -> next;
    }
    else
    {
	retVal = newp(Vertex);
    }

    if (retVal)
    {
	retVal -> normal[0] = retVal -> normal[1] = 0.0;
	retVal -> normal[2] = 1.0;
	retVal -> position[0] = retVal -> position[1] = retVal -> position[2] = 0.0;
	retVal -> colorIndex = 0;

	retVal -> next = 0;
	retVal -> flags = flags;
	if (((PicPtr) picture) -> lastVertex)
	{
	    (((PicPtr) picture) -> lastVertex) -> next = retVal;
	    ((PicPtr) picture) -> lastVertex = retVal;
	}
	else
	{
	    ((PicPtr) picture) -> lastVertex = retVal;
	    ((PicPtr) picture) -> vertices = retVal;
	}
    }
    return retVal;
}

#ifdef HAVE_PROTOTYPES
VertexPtr InsertVertex(ObjPtr picture, VertexPtr after)
#else
VertexPtr InsertVertex(picture, after)
ObjPtr picture;
VertexPtr after;
#endif
/*Returns a new vertex, stuck on picture after after, with VF_SAMEPLACE*/
{
    VertexPtr retVal;

    if (freedVertices)
    {
	retVal = freedVertices;
	freedVertices = freedVertices -> next;
    }
    else
    {
	retVal = newp(Vertex);
    }

    if (retVal)
    {
	retVal -> normal[0] = retVal -> normal[1] = 0.0;
	retVal -> normal[2] = 1.0;
	retVal -> position[0] = retVal -> position[1] = retVal -> position[2] = 0.0;
	retVal -> colorIndex = 0;

	retVal -> next = after -> next;
	retVal -> flags = VF_SAMEPLACE;

	after -> next = retVal;

	if (((PicPtr) picture) -> lastVertex == after)
	{
	    ((PicPtr) picture) -> lastVertex = retVal;
	}
    }
    return retVal;
}

#ifdef HAVE_PROTOTYPES
void FreeVertex(VertexPtr vertex)
#else
void FreeVertex(vertex)
VertexPtr vertex;
#endif
/*Frees a vertex*/
{
    vertex -> next = freedVertices;
    freedVertices = vertex;
}

static ObjPtr RegisterPicField(pic, whichField)
ObjPtr pic;
int whichField;
/*Registers pic in a field*/
{
    curFields[whichField] . objectInfo = pic;
    return ObjTrue;
}

void InitPictures()
/*Initializes the pictures system*/
{
    int k;
    double theta;

    writeMasks[curWriteMask] . rWrite = true;
    writeMasks[curWriteMask] . gWrite = true;
    writeMasks[curWriteMask] . bWrite = true;
    writeMasks[curWriteMask] . zWrite = true;

#ifdef GRAPHICS
    deflinestyle(DASHEDLINE, dashedLine);
    deflinestyle(DOTTEDLINE, dottedLine);

    /*Define the patterns*/
    defpattern(GREYPAT, 16, greyPattern);
    defpattern(EVENPAT, 16, evenPattern);
    defpattern(ODDPAT, 16, oddPattern);
    defpattern(ALTGREYPAT, 16, altGreyPattern); 

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyGreyPattern[k] & evenPattern[k];
    defpattern(EVENGREYPAT, 16, tempPattern);

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyGreyPattern[k] & oddPattern[k];
    defpattern(ODDGREYPAT, 16, tempPattern);

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyAltGreyPattern[k] & evenPattern[k];
    defpattern(EVENALTGREYPAT, 16, tempPattern);

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyAltGreyPattern[k] & oddPattern[k];
    defpattern(ODDALTGREYPAT, 16, tempPattern);

    /*Define the cardinal spline basis*/
    defbasis(CARDBASIS, cardBasis);
#endif

    picClass = NewObject(NULLOBJ, 0);
    AddToReferenceList(picClass);
    SetMethod(picClass, REGISTERFIELD, RegisterPicField);

    DeclareIndirectDependency(picClass, PICDEFORMED, REPOBJ, PICDEFORMED);
    SetMethod(picClass, PICDEFORMED, MakePictureDeformed);

    /*Initialize cef*/
    theta = M_PI_4;
    for (k = 0; k < MAXSUB; ++k)
    {
	cefs[k] = 1.0 / rcos(theta);
	cefs_2[k] = cefs[k] * 0.5;
	theta *= 0.5;
    }
}

void KillPictures()
/*Kills the pictures system*/
{
    DeleteThing(picClass);
}

static ObjPtr CleanupPicture(pic)
PicPtr pic;
/*Cleans up pic*/
{
    while (pic -> vertices)
    {
	VertexPtr next;
	next = pic -> vertices -> next;
	Free(pic -> vertices);
	pic -> vertices = next;
    }
    pic -> lastVertex = 0;

    while (pic -> items)
    {
	PicItemPtr next;
	if (pic -> items -> type == POLYGONS)
	{
	    PolyPtr runner, nr;
	    runner = ((PolysPtr) pic -> items) -> polygons;
	    while (runner)
	    {
		nr = (PolyPtr) runner -> item . next;
		Free(runner);
		runner = nr;
	    }
	}
	next = pic -> items -> next;
	Free(pic -> items);
	pic -> items = next;
    }
    pic -> lastItem = 0;
    return ObjTrue;
}

ObjPtr NewPicture()
/*Returns a new empty picture or 0 if it could not make one.*/
{
    PicPtr retVal;
    retVal = (PicPtr) NewObject(picClass, sizeof(Picture) - sizeof(Thing));
    if (retVal)
    {
        SETOBJTYPE(retVal -> thing . flags, PICTURE);
	SetMethod((ObjPtr) retVal, CLEANUP, CleanupPicture);
	retVal -> vertices = 0;
	retVal -> lastVertex = 0;

	retVal -> items = 0;
	retVal -> lastItem = 0;
	return (ObjPtr) retVal;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolysPtr AppendPolysToPicture(ObjPtr picture)
#else
PolysPtr AppendPolysToPicture(picture)
ObjPtr picture;
#endif
/*Returns a new empty set of polygons 0 if it could not make one.*/
{
    PolysPtr item;
    item = (PolysPtr) Alloc(sizeof(Polygons));
    if (item)
    {
	item -> item . type = POLYGONS;
	item -> item . proximity = 0;
	item -> item . flags = 0;

	item -> picture = (PicPtr) picture;
	item -> polygons = 0;
	item -> lastItem = 0;
	item -> enclosed = false;

	AppendItemToPicture((PicPtr) picture, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendPolyToPolys(PolysPtr polys, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolyToPolys(polys, nVertices, vertices)
PolysPtr polys;
long nVertices;
Vertex vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to polygons polys.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYGON;
	item -> item . proximity = 0;
	item -> item . flags = 0;

	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex((ObjPtr) polys -> picture, 0);
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] = vertices[k] . position[d];
		item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
	    }
	    item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
	}
	if (polys -> polygons)
	{
	    /*It's not the first*/
	    polys -> lastItem -> item . next = (PicItemPtr) item;
	}
	else
	{
	    /*It is the first*/
	    polys -> polygons = item;
	}
	item -> item . next = 0;
	polys -> lastItem = item;
	return item;
    }
    else
    {
	OMErr();
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendSPolyToPolys(PolysPtr polys, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolyToPolys(polys, nVertices, vertices)
PolysPtr polys;
long nVertices;
VertexPtr vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of VertexPtr
  to polygons polys.  Returns true iff successful.  The vertices are 
  presumed to be shared.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYGON;
	item -> item . proximity = 0;
	item -> item . flags = 0;

	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = vertices[k];
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] = vertices[k] -> position[d];
		item -> vertices[k] -> normal[d] = vertices[k] -> normal[d];
	    }
	    item -> vertices[k] -> colorIndex = vertices[k] -> colorIndex;
	}
	if (polys -> polygons)
	{
	    /*It's not the first*/
	    polys -> lastItem -> item . next = (PicItemPtr) item;
	}
	else
	{
	    /*It is the first*/
	    polys -> polygons = item;
	}
	item -> item . next = 0;
	polys -> lastItem = item;
	return item;
    }
    else
    {
	OMErr();
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr TesselateSPolyToPicture(ObjPtr pic, long nVertices, VertexPtr vertices[])
#else
PolyPtr TesselateSPolyToPicture(pic, nVertices, vertices)
ObjPtr pic;
long nVertices;
VertexPtr vertices[];
#endif
/*Tesselates a polygon with nVertices where vertices is an array of VertexPtr
  to picture pic.  Returns the resulting item iff successful.  The vertices are 
  presumed to be shared.  It's a very simple-minded tesselation, not guaranteed
  to work for non convex polygons*/
{
    PolyPtr item;
    int k, d;
    register VertexPtr *v1, *v2;

    v1 = vertices;
    v2 = vertices + nVertices - 1;


    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYTRI;
	item -> item . proximity = 0;
	item -> item . flags = 0;

	item -> nVertices = nVertices;

	/*Fill the vertices*/
	k = 0;

	item -> vertices[k++] = *(v1++);
	while (v2 >= v1)
	{
	    item -> vertices[k++] = *(v1++);
	    if (v2 >= v1)
	    {
		item -> vertices[k++] = *(v2--);
	    }
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	OMErr();
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
static PolyPtr ConvertPolyOntoPolys(PolysPtr polys, PolyPtr poly)
#else
static PolyPtr ConvertPolyOntoPolys(polys, poly)
PolysPtr polys;
PolyPtr poly;
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to polygons polys.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (poly -> nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYGON;
	item -> item . proximity = poly -> item . proximity;
	item -> item . flags = 0;

	item -> nVertices = poly -> nVertices;
	for (k = 0; k < poly -> nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex((ObjPtr) polys -> picture, 0);
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] =
			poly -> vertices[k] -> position[d];
		item -> vertices[k] -> normal[d] =
			poly -> vertices[k] -> normal[d];
	    }
	    item -> vertices[k] -> colorIndex =
		poly -> vertices[k] -> colorIndex;
	}
	if (polys -> polygons)
	{
	    /*It's not the first*/
	    polys -> lastItem -> item . next = (PicItemPtr) item;
	}
	else
	{
	    /*It is the first*/
	    polys -> polygons = item;
	}
	item -> item . next = 0;
	polys -> lastItem = item;
	return item;
    }
    else
    {
	OMErr();
	return 0;
    }
}

static void AppendItemToPicture(pic, item)
PicPtr pic;
PicItemPtr item;
/*Appends item to the end of pic*/
{
    if (pic -> items)
    {
	/*It's not the first*/
	pic -> lastItem -> next = item;
    }
    else
    {
	/*It is the first*/
	pic -> items = item;
    }
    item -> next = 0;
    pic -> lastItem = item;
}

#ifdef HAVE_PROTOTYPES
RectMeshPtr AppendRectMeshToPicture(ObjPtr picture, long xDim, long yDim, Bool inCenter)
#else
RectMeshPtr AppendRectMeshToPicture(picture, xDim, yDim, inCenter)
ObjPtr picture;
long yDim, xDim;
Bool inCenter;
#endif
/*Appends a new rectangular mesh.  inCenter if nodes in center*/
{
    RectMeshPtr item;
    long nVertices;
    long k;
    Bool first;

    nVertices = (xDim * (2 * yDim - 1));
    item = (RectMeshPtr)
	   Alloc(sizeof(RectMesh) + (nVertices - 1) *
			sizeof(VertexPtr));

    if (inCenter)
    {
	/*Make the nodes with the center nodes canonical*/
	item -> vertices[0] = NewVertex(picture, 0);
	item -> vertices[1] = NewVertex(picture, VF_FIRSTCANON);
	for (k = 2; k < nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex(picture, (k & 1) ? VF_NEXTCANON : 0);
	}
    }
    else
    {
	/*Make as if vertex nodes are canonical*/
	item -> vertices[0] = NewVertex(picture, VF_FIRSTCANON);
	for (k = 1; k < nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex(picture, (k & 1) ? 0 : VF_NEXTCANON);
	}
    }

    item -> item . type = RECTMESH;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> inCenter = inCenter;
    item -> xDim = xDim;
    item -> yDim = yDim;

    AppendItemToPicture((PicPtr) picture, (PicItemPtr) item);
    return item;
}

#ifdef HAVE_PROTOTYPES
VertexPtr RectMeshVertex(RectMeshPtr rectMesh, long i, long j)
#else
VertexPtr RectMeshVertex(rectMesh, i, j)
RectMeshPtr rectMesh;
long i, j;
#endif
/*Returns vertex i, j of a rect mesh*/
{
    long offset;
    offset = RECTMESHVERTEX(rectMesh, i, j);
    return rectMesh -> vertices[offset]; 
}

#ifdef HAVE_PROTOTYPES
VertexPtr RectMeshCenter(RectMeshPtr rectMesh, long i, long j)
#else
VertexPtr RectMeshCenter(rectMesh, i, j)
RectMeshPtr rectMesh;
long i, j;
#endif
/*Returns center i, j of a rect mesh*/
{
    long offset;
    offset = RECTMESHCENTER(rectMesh, i, j);
    return rectMesh -> vertices[offset]; 
}

#define ISVECZERO(v) (((v)[0] == 0.0) && ((v)[1] == 0.0) && ((v)[2] == 0.0))

#ifdef HAVE_PROTOTYPES
void CalcRectNormals(RectMeshPtr rectMesh)
#else
void CalcRectNormals(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Calculates the normals in rectMesh and copies them to the vertices.
  Obviously, do this after the positions have been set.*/
{
    register long i, j, k, offset, offset2, offset3, iDim, jDim, nX;
    register VertexPtr *vertices;
    register long nVertices;
    float normal[3];		/*Normal at center vertex*/
    float position1[3];		/*Position of base vertex*/
    float vI[3], vJ[3];		/*Vectors for cross product*/
    float vI0[3], vI1[3], vJ0[3], vJ1[3];
    float vt[3];		/*Test vertex*/
    float cn[3];		/*Candidate normal*/

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;
    vertices = rectMesh -> vertices;
    nVertices = (iDim * (2 * jDim - 1));
    nX = jDim * 2 - 1;

#define nY 2

#ifdef OLDRECTNORMALS
    /*Make all the _vertex_ normals*/
    for (j = 0; j < jDim; ++j)
    {
	for (i = 0; i < iDim; ++i)
	{
	    vI[0] = vI[1] = vI[2] = 0.0;
	    vJ[0] = vJ[1] = vJ[2] = 0.0;

	    offset = RECTMESHVERTEX(rectMesh, i, j);

	    /*Get position here at vertex*/
	    position1[0] = vertices[offset] -> position[0];
	    position1[1] = vertices[offset] -> position[1];
	    position1[2] = vertices[offset] -> position[2];

	    /*Make i component*/
	    if (i >= iDim - 1)
	    {
		/*Get from i - 1*/
		offset2 = offset - nX;
		vI[0] += position1[0] - vertices[offset2] -> position[0];
		vI[1] += position1[1] - vertices[offset2] -> position[1];
		vI[2] += position1[2] - vertices[offset2] -> position[2];
	    }
	    else if (i <= 0)
	    {
		/*Get from i + 1*/
		offset2 = offset + nX;
		vI[0] += vertices[offset2] -> position[0] - position1[0];
		vI[1] += vertices[offset2] -> position[1] - position1[1];
		vI[2] += vertices[offset2] -> position[2] - position1[2];
	    }
	    else
	    {
		/*Get from i + 1 and i - 1*/
		offset2 = offset - nX;
		offset3 = offset + nX;
		vI[0] += vertices[offset3] -> position[0] - vertices[offset2] -> position[0];
		vI[1] += vertices[offset3] -> position[1] - vertices[offset2] -> position[1];
		vI[2] += vertices[offset3] -> position[2] - vertices[offset2] -> position[2];
	    }

	    /*Make j component*/
	    if (j >= jDim - 1)
	    {
		/*Get from j - 1*/
		offset2 = offset - nY;
		vJ[0] += position1[0] - vertices[offset2] -> position[0];
		vJ[1] += position1[1] - vertices[offset2] -> position[1];
		vJ[2] += position1[2] - vertices[offset2] -> position[2];
	    }
	    else if (j <= 0)
	    {
		/*Get from j + 1*/
		offset2 = offset + nY;
		vJ[0] += vertices[offset2] -> position[0] - position1[0];
		vJ[1] += vertices[offset2] -> position[1] - position1[1];
		vJ[2] += vertices[offset2] -> position[2] - position1[2];
	    }
	    else
	    {
		/*Get from j + 1 and j - 1*/
		offset2 = offset - nY;
		offset3 = offset + nY;
		vJ[0] += vertices[offset3] -> position[0] - vertices[offset2] -> position[0];
		vJ[1] += vertices[offset3] -> position[1] - vertices[offset2] -> position[1];
		vJ[2] += vertices[offset3] -> position[2] - vertices[offset2] -> position[2];
	    }

	    /*Cross i by j to get normal*/
	    CROSS(vI, vJ, normal);

	    if (normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0)
	    {
		normal[2] = 1.0;
	    }

	    /*Normalize it*/
	    NORMALIZE(normal);
	    vertices[offset] -> normal[0] += normal[0];
	    vertices[offset] -> normal[1] += normal[1];
	    vertices[offset] -> normal[2] += normal[2];
	}
    }

    /*Interpolate the centers with four corners*/
    for (i = 1; i < iDim - 1; ++i)
    {
	for (j = 1; j < jDim - 1; ++j)
	{
	    /*Lower left corner*/
	    offset = RECTMESHVERTEX(rectMesh, i - 1, j - 1);
	    normal[0] = vertices[offset] -> normal[0];
	    normal[1] = vertices[offset] -> normal[1];
	    normal[2] = vertices[offset] -> normal[2];

	    /*Upper left corner*/
	    offset = RECTMESHVERTEX(rectMesh, i - 1, j);
	    normal[0] += vertices[offset] -> normal[0];
	    normal[1] += vertices[offset] -> normal[1];
	    normal[2] += vertices[offset] -> normal[2];

	    /*Lower right corner*/
	    offset = RECTMESHVERTEX(rectMesh, i, j - 1);
	    normal[0] += vertices[offset] -> normal[0];
	    normal[1] += vertices[offset] -> normal[1];
	    normal[2] += vertices[offset] -> normal[2];

	    /*Upper right corner*/
	    offset = RECTMESHVERTEX(rectMesh, i, j);
	    normal[0] += vertices[offset] -> normal[0];
	    normal[1] += vertices[offset] -> normal[1];
	    normal[2] += vertices[offset] -> normal[2];

	    if (normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0)
	    {
		normal[2] = 1.0;
	    }

	    /*Normalize it*/
	    NORMALIZE(normal);

	    /*Stuff them*/
	    offset = RECTMESHCENTER(rectMesh, i, j);
	    vertices[offset] -> normal[0] += normal[0];
	    vertices[offset] -> normal[1] += normal[1];
	    vertices[offset] -> normal[2] += normal[2];
	}
    }

#else
    /*Zero all the vertices*/
    for (k = 0; k < nVertices; ++k)
    {
	vertices[k] -> normal[0] = 0.0;
	vertices[k] -> normal[1] = 0.0;
	vertices[k] -> normal[2] = 0.0;
    }

    /*Make all the center normals from edges*/
    for (j = 0; j < jDim - 1; ++j)
    {
	for (i = 0; i < iDim - 1; ++i)
	{
	    offset = RECTMESHVERTEX(rectMesh, i, j);

	    /*Get position here at vertex*/
	    position1[0] = vertices[offset] -> position[0];
	    position1[1] = vertices[offset] -> position[1];
	    position1[2] = vertices[offset] -> position[2];

	    /*Make i0 component*/
	    offset2 = offset + nX;
	    vI0[0] = vertices[offset2] -> position[0] - position1[0];
	    vI0[1] = vertices[offset2] -> position[1] - position1[1];
	    vI0[2] = vertices[offset2] -> position[2] - position1[2];

	    /*Make j0 component*/
	    offset2 = offset + nY;
	    vJ0[0] = vertices[offset2] -> position[0] - position1[0];
	    vJ0[1] = vertices[offset2] -> position[1] - position1[1];
	    vJ0[2] = vertices[offset2] -> position[2] - position1[2];

	    /*Get position at catty corner vertex*/
	    offset = offset + nX + nY;
	    position1[0] = vertices[offset] -> position[0];
	    position1[1] = vertices[offset] -> position[1];
	    position1[2] = vertices[offset] -> position[2];

	    /*Make i1 component*/
	    offset2 = offset - nX;
	    vI1[0] = position1[0] - vertices[offset2] -> position[0];
	    vI1[1] = position1[1] - vertices[offset2] -> position[1];
	    vI1[2] = position1[2] - vertices[offset2] -> position[2];

	    /*Make j1 component*/
	    offset2 = offset - nY;
	    vJ1[0] = position1[0] - vertices[offset2] -> position[0];
	    vJ1[1] = position1[1] - vertices[offset2] -> position[1];
	    vJ1[2] = position1[2] - vertices[offset2] -> position[2];


	    /*Now use the four candidates to get an approximate normal*/
	    normal[0] = normal[1] = normal[2] = 0.0;
	    CROSS(vI0, vJ0, cn);
	    normal[0] += cn[0];
	    normal[1] += cn[1];
	    normal[2] += cn[2];
	    CROSS(vI1, vJ0, cn);
	    normal[0] += cn[0];
	    normal[1] += cn[1];
	    normal[2] += cn[2];
	    CROSS(vI0, vJ1, cn);
	    normal[0] += cn[0];
	    normal[1] += cn[1];
	    normal[2] += cn[2];
	    CROSS(vI1, vJ1, cn);
	    normal[0] += cn[0];
	    normal[1] += cn[1];
	    normal[2] += cn[2];

	    if (normal[0] != 0.0 || normal[1] != 0.0 || normal[2] != 0.0)
	    {
		/*It isn't degenerate.  Normalize it.*/
		NORMALIZE(normal);

		/*Set it*/
		offset = RECTMESHCENTER(rectMesh, i, j);
		vertices[offset] -> normal[0] = normal[0];
		vertices[offset] -> normal[1] = normal[1];
		vertices[offset] -> normal[2] = normal[2];
	    }
	    /*If not set, a center will be considered degenerate.*/
	}
    }

    /*Make all the vertex normals from center*/
    for (j = 0; j < jDim - 1; ++j)
    {
	for (i = 0; i < iDim - 1; ++i)
	{
	    /*Get pointer for this vertex*/
	    offset = RECTMESHCENTER(rectMesh, i, j);

	    /*Get temporary normal*/
	    normal[0] = vertices[offset] -> normal[0];
	    normal[1] = vertices[offset] -> normal[1];
	    normal[2] = vertices[offset] -> normal[2];

	    /*Get pointer for lower left center*/
	    offset2 = RECTMESHVERTEX(rectMesh, i, j);

	    /*Distribute this normal to all four corners*/

	    vertices[offset2] -> normal[0] += normal[0];
	    vertices[offset2] -> normal[1] += normal[1];
	    vertices[offset2] -> normal[2] += normal[2];

	    offset2 += nX; 

	    vertices[offset2] -> normal[0] += normal[0];
	    vertices[offset2] -> normal[1] += normal[1];
	    vertices[offset2] -> normal[2] += normal[2];

	    offset2 += nY; 

	    vertices[offset2] -> normal[0] += normal[0];
	    vertices[offset2] -> normal[1] += normal[1];
	    vertices[offset2] -> normal[2] += normal[2];

	    offset2 -= nX; 
	    vertices[offset2] -> normal[0] += normal[0];
	    vertices[offset2] -> normal[1] += normal[1];
	    vertices[offset2] -> normal[2] += normal[2];
	}
    }

    /*Go through and fix degenerate normals*/
    for (k = 0; k < nVertices; ++k)
    {
	if (ISVECZERO(vertices[k] -> normal))
	{
	    vertices[k] -> normal[2] = 1.0;
	}
	else
	{
	    NORMALIZE(vertices[k] -> normal);
	}
    }
#endif

#undef nY

}

#ifdef HAVE_PROTOTYPES
void InterpRectCenters(RectMeshPtr rectMesh)
#else
void InterpRectCenters(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Interpolates rect centers from vertices*/
{
    register long i, j, offset, iDim, jDim, nX;
    register VertexPtr *vertices;
    float position[3];		/*Position of center vertex*/
    int colorIndex;		/*Color within color table*/

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;
    vertices = rectMesh -> vertices;
    nX = jDim * 2 - 1;

    for (i = 0; i < iDim - 1; ++i)
    {
	for (j = 0; j < jDim - 1; ++j)
	{
	    /*Lower left corner*/
	    offset = RECTMESHVERTEX(rectMesh, i, j);
	    position[0] = vertices[offset] -> position[0];
	    position[1] = vertices[offset] -> position[1];
	    position[2] = vertices[offset] -> position[2];
	    colorIndex = vertices[offset] -> colorIndex;

	    /*Upper left corner*/
	    offset = RECTMESHVERTEX(rectMesh, i, j + 1);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;

	    /*Lower right corner*/
	    offset = RECTMESHVERTEX(rectMesh, i + 1, j);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;

	    /*Upper right corner*/
	    offset = RECTMESHVERTEX(rectMesh, i + 1, j + 1);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;

	    /*Normalize results*/
	    position[0] *= 0.25;
	    position[1] *= 0.25;
	    position[2] *= 0.25;
	    colorIndex /= 4;

	    /*Stuff them*/
	    offset = RECTMESHCENTER(rectMesh, i, j);
	    vertices[offset] -> position[0] = position[0];
	    vertices[offset] -> position[1] = position[1];
	    vertices[offset] -> position[2] = position[2];
	    vertices[offset] -> colorIndex = colorIndex;
	}
    }
}

#ifdef HAVE_PROTOTYPES
void InterpRectVertices(RectMeshPtr rectMesh)
#else
void InterpRectVertices(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Interpolates and extrapolates rect vertices from centers*/
{
    register long i, j, k, offset, iDim, jDim, nX;
    register VertexPtr *vertices;
    float position[3];		/*Position of center vertex*/
    int colorIndex;		/*Color within color table*/
    int nColorIndices;

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;

    vertices = rectMesh -> vertices;
    nX = jDim * 2 - 1;

    /*Interpolate the centers with four corners*/
    for (i = 1; i < iDim - 1; ++i)
    {
	for (j = 1; j < jDim - 1; ++j)
	{
	    nColorIndices = 0;

	    /*Lower left corner*/
	    offset = RECTMESHCENTER(rectMesh, i - 1, j - 1);
	    position[0] = vertices[offset] -> position[0];
	    position[1] = vertices[offset] -> position[1];
	    position[2] = vertices[offset] -> position[2];
	    colorIndex = vertices[offset] -> colorIndex;
	    if (vertices[offset] -> colorIndex) ++nColorIndices;

	    /*Upper left corner*/
	    offset = RECTMESHCENTER(rectMesh, i - 1, j);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;
	    if (vertices[offset] -> colorIndex) ++nColorIndices;

	    /*Lower right corner*/
	    offset = RECTMESHCENTER(rectMesh, i, j - 1);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;
	    if (vertices[offset] -> colorIndex) ++nColorIndices;

	    /*Upper right corner*/
	    offset = RECTMESHCENTER(rectMesh, i, j);
	    position[0] += vertices[offset] -> position[0];
	    position[1] += vertices[offset] -> position[1];
	    position[2] += vertices[offset] -> position[2];
	    colorIndex += vertices[offset] -> colorIndex;
	    if (vertices[offset] -> colorIndex) ++nColorIndices;

	    /*Normalize results*/
	    position[0] *= 0.25;
	    position[1] *= 0.25;
	    position[2] *= 0.25;
#if 1
	    if (nColorIndices) colorIndex /= nColorIndices;
	    printf("%d\n", nColorIndices);
#else
	    colorIndex /= 4;
#endif

	    /*Stuff them*/
	    offset = RECTMESHVERTEX(rectMesh, i, j);
	    vertices[offset] -> position[0] = position[0];
	    vertices[offset] -> position[1] = position[1];
	    vertices[offset] -> position[2] = position[2];
	    vertices[offset] -> colorIndex = colorIndex;
	}
    }

    /*Extrapolate left and right sides*/
    for (j = 1; j < jDim - 1; ++j)
    {
	/*Left side*/
	offset = RECTMESHVERTEX(rectMesh, 1, j);
	for (k = 0; k < 3; ++k)
	{
	    position[k] = vertices[offset] -> position[k];
	}
	colorIndex = vertices[offset] -> colorIndex;
	if (iDim >= 3)
	{
	    /*Extend position only*/
	    offset = RECTMESHVERTEX(rectMesh, 2, j);
	    for (k = 0; k < 3; ++k)
	    {
		position[k] += position[k] - vertices[offset] -> position[k];
	    }
	}
	/*Stuff them back*/
	offset = RECTMESHVERTEX(rectMesh, 0, j);
	for (k = 0; k < 3; ++k)
	{
	    vertices[offset] -> position[k] = position[k];
	}
	vertices[offset] -> colorIndex = colorIndex;

	/*Right side*/
	offset = RECTMESHVERTEX(rectMesh, iDim - 2, j);
	for (k = 0; k < 3; ++k)
	{
	    position[k] = vertices[offset] -> position[k];
	}
	colorIndex = vertices[offset] -> colorIndex;
	if (iDim >= 3)
	{
	    /*Extend position only*/
	    offset = RECTMESHVERTEX(rectMesh, iDim - 3, j);
	    for (k = 0; k < 3; ++k)
	    {
		position[k] += position[k] - vertices[offset] -> position[k];
	    }
	}
	/*Stuff them back*/
	offset = RECTMESHVERTEX(rectMesh, iDim - 1, j);
	for (k = 0; k < 3; ++k)
	{
	    vertices[offset] -> position[k] = position[k];
	}
	vertices[offset] -> colorIndex = colorIndex;
    }	

    /*Extrapolate bottom and top sides*/
    for (i = 1; i < iDim - 1; ++i)
    {
	/*Left side*/
	offset = RECTMESHVERTEX(rectMesh, i, 1);
	for (k = 0; k < 3; ++k)
	{
	    position[k] = vertices[offset] -> position[k];
	}
	colorIndex = vertices[offset] -> colorIndex;
	if (jDim >= 3)
	{
	    /*Extend position only*/
	    offset = RECTMESHVERTEX(rectMesh, i, 2);
	    for (k = 0; k < 3; ++k)
	    {
		position[k] += position[k] - vertices[offset] -> position[k];
	    }
	}
	/*Stuff them back*/
	offset = RECTMESHVERTEX(rectMesh, i, 0);
	for (k = 0; k < 3; ++k)
	{
	    vertices[offset] -> position[k] = position[k];
	}
	vertices[offset] -> colorIndex = colorIndex;

	/*Right side*/
	offset = RECTMESHVERTEX(rectMesh, i, jDim - 2);
	for (k = 0; k < 3; ++k)
	{
	    position[k] = vertices[offset] -> position[k];
	}
	colorIndex = vertices[offset] -> colorIndex;
	if (jDim >= 3)
	{
	    /*Extend position only*/
	    offset = RECTMESHVERTEX(rectMesh, i, jDim - 3);
	    for (k = 0; k < 3; ++k)
	    {
		position[k] += position[k] - vertices[offset] -> position[k];
	    }
	}
	/*Stuff them back*/
	offset = RECTMESHVERTEX(rectMesh, i, jDim - 1);
	for (k = 0; k < 3; ++k)
	{
	    vertices[offset] -> position[k] = position[k];
	}
	vertices[offset] -> colorIndex = colorIndex;
    }

    /*Do lower left corner*/
    offset = RECTMESHCENTER(rectMesh, 0, 0);
    for (k = 0; k < 3; ++k)
    {
	position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, 1, 1);
    for (k = 0; k < 3; ++k)
    {
	position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, 0, 0);
    for (k = 0; k < 3; ++k)
    {
	vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do upper left corner*/
    offset = RECTMESHCENTER(rectMesh, 0, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
	position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, 1, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
	position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, 0, jDim - 1);
    for (k = 0; k < 3; ++k)
    {
	vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do lower right corner*/
    offset = RECTMESHCENTER(rectMesh, iDim - 2, 0);
    for (k = 0; k < 3; ++k)
    {
	position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, iDim - 2, 1);
    for (k = 0; k < 3; ++k)
    {
	position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, iDim - 1, 0);
    for (k = 0; k < 3; ++k)
    {
	vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do upper left corner*/
    offset = RECTMESHCENTER(rectMesh, iDim - 2, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
	position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, iDim - 2, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
	position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, iDim - 1, jDim - 1);
    for (k = 0; k < 3; ++k)
    {
	vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendPolyToPicture(ObjPtr pic, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolyToPicture(pic, nVertices, vertices)
ObjPtr pic;
long nVertices;
Vertex vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to picture pic.  Returns polygon iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYGON;
	item -> item . proximity = 0;
	item -> item . flags = 0;

	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex(pic, 0);
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] = vertices[k] . position[d];
		item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
	    }
	    item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
static PolyPtr ConvertPolyOntoPicture(ObjPtr pic, PolyPtr poly)
#else
static PolyPtr ConvertPolyOntoPicture(pic, poly)
ObjPtr pic;
PolyPtr poly;
#endif
/*Converts a polygon or polyline onto picture.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (poly -> nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = poly -> item . type;
	item -> lineWidth = poly -> lineWidth;
	item -> item . proximity = poly -> item . proximity;
	item -> item . flags = 0;

	item -> nVertices = poly -> nVertices;
	for (k = 0; k < poly -> nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex(pic, 0);
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] =
			poly -> vertices[k] -> position[d];
		item -> vertices[k] -> normal[d] =
			poly -> vertices[k] -> normal[d];
	    }
	    item -> vertices[k] -> colorIndex =
		poly -> vertices[k] -> colorIndex;
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendPolylineToPicture(ObjPtr pic, int lineWidth, int proximity, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolylineToPicture(pic, lineWidth, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
Vertex vertices[];
#endif
/*Append a polyline with nVertices where vertices is an array of Vertex
  to picture pic.  proximity is extra proximity for lines that are to
  go above something.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYLINE;
	item -> item . proximity = proximity;
	item -> item . flags = 0;

	item -> lineWidth = lineWidth;
	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = NewVertex(pic, 0);
	    for (d = 0; d < 3; ++d)
	    {
		item -> vertices[k] -> position[d] = vertices[k] . position[d];
		item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
	    }
	    item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendSPolylineToPicture(ObjPtr pic, int lineWidth, int proximity, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolylineToPicture(pic, lineWidth, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
VertexPtr vertices[];
#endif
/*Append a polyline with nVertices where vertices is an array of Vertex
  to picture pic.  proximity is extra proximity for lines that are to
  go above something.  Returns true iff successful.*/
{
    register PolyPtr item;
    register int k, d;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYLINE;
	item -> item . proximity = proximity;
	item -> item . flags = 0;

	item -> lineWidth = lineWidth;
	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = vertices[k];
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
PolyPtr AppendSPolypointToPicture(ObjPtr pic, int proximity, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolypointToPicture(pic, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
Vertex *vertices[];
#endif
/*Append a polypoint with nVertices where vertices is an array of VertexPtr
  to picture pic.  proximity is extra proximity for points that are to
  go above something.  Returns true iff successful.*/
{
    PolyPtr item;
    int k;

    item = (PolyPtr)
	   Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
	item -> item . type = POLYPOINT;
	item -> item . proximity = proximity;

	item -> lineWidth = 1;
	item -> nVertices = nVertices;
	for (k = 0; k < nVertices; ++k)
	{
	    item -> vertices[k] = vertices[k];
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
FrustumPtr AppendFrustumToPicture(ObjPtr pic, float end1[3], float rad1,
		float end2[3], float rad2, int colorIndex)
#else
FrustumPtr AppendFrustumToPicture(pic, end1, rad1, end2, rad2, colorIndex)
ObjPtr pic;
float end1[3];
float rad1;
float end2[3];
float rad2;
int colorIndex;
#endif
/*Appends a conical frustum with endpoints end1 and end2 and radii rad1 and
  rad2 to picture pic*/
{
    FrustumPtr item;
    int k;

    item = newp(Frustum);
    if (item)
    {
	item -> item . type = FRUSTUM;
	item -> rad1 = rad1;
	item -> rad2 = rad2;
	item -> colorIndex = colorIndex;
	for (k = 0; k < 3; ++k)
	{
	    item -> end1[k] = end1[k];
	    item -> end2[k] = end2[k];
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

#ifdef HAVE_PROTOTYPES
SpherePtr AppendSphereToPicture(ObjPtr pic, float center[3], float radius, int colorIndex)
#else
SpherePtr AppendSphereToPicture(pic, center, radius, colorIndex)
ObjPtr pic;
float center[3];
float radius;
int colorIndex;
#endif
/*Appends a sphere with center center and radius radius to pic*/
{
    SpherePtr item;
    int k;

    item = newp(Sphere);
    if (item)
    {
	item -> item . type = SPHERE;
	item -> radius = radius;
	item -> centerVertex = NewVertex(pic, 0);
	item -> centerVertex -> colorIndex = colorIndex;
	for (k = 0; k < 3; ++k)
	{
	    item -> centerVertex -> position[k] = center[k];
	}
	AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
	return item;
    }
    else
    {
	return 0;
    }
}

static float curCenter[3];	/*Current center*/
static float curRadius;		/*Current radius*/
static int curSub;		/*Current number of subdivisions*/
static int curNSides;		/*Current number of sides*/

#ifdef GRAPHICS
void DrawQuadrant(r1, r2, r3, level)
float r1[3];
float r2[3];
float r3[3];
int level;
/*Draws a quadrant of sphere at center curCenter with radius curRadius given 
  radial normal vectors r1, r2, and r3, using current color table, at recursion 
  level level.  Assumes that the appropriate color table is current*/
{
    if (level <= 0)
    {
	/*Maximum number of subdivisions has been reached*/
	float v1[3], v2[3], v3[3];
	v1[0] = curCenter[0] + r1[0] * curRadius;
	v1[1] = curCenter[1] + r1[1] * curRadius;
	v1[2] = curCenter[2] + r1[2] * curRadius;
	v2[0] = curCenter[0] + r2[0] * curRadius;
	v2[1] = curCenter[1] + r2[1] * curRadius;
	v2[2] = curCenter[2] + r2[2] * curRadius;
	v3[0] = curCenter[0] + r3[0] * curRadius;
	v3[1] = curCenter[1] + r3[1] * curRadius;
	v3[2] = curCenter[2] + r3[2] * curRadius;

        if (rgbp)
        {
            bgnpolygon();
	    N3F(r1);
	    v3f(v1);
	    N3F(r2);
	    v3f(v2);
	    N3F(r3);
	    v3f(v3);
	    endpolygon();
        }
	else
        {
	    bgnpolygon();
	    v3f(v1);
	    v3f(v2);
	    v3f(v3);
            endpolygon();
        }
    }
    else
    {
	/*Do another subdivision*/
	float r12[3], r23[3], r31[3];	/*Inner triangle subdivisions*/
	r12[0] = (r1[0] + r2[0]);
	r12[1] = (r1[1] + r2[1]);
	r12[2] = (r1[2] + r2[2]);
	r23[0] = (r2[0] + r3[0]);
	r23[1] = (r2[1] + r3[1]);
	r23[2] = (r2[2] + r3[2]);
	r31[0] = (r3[0] + r1[0]);
	r31[1] = (r3[1] + r1[1]);
	r31[2] = (r3[2] + r1[2]);
	NORMALIZE(r12);
	NORMALIZE(r23);
	NORMALIZE(r31);

	/*Draw the subdivisions*/
	DrawQuadrant(r1, r12, r31, level - 1);
	DrawQuadrant(r12, r2, r23, level - 1);
	DrawQuadrant(r23, r3, r31, level - 1);
	DrawQuadrant(r12, r23, r31, level - 1);
    }
}
#endif

static Bool formsAreEqual;

#ifdef HAVE_PROTOTYPES
void UnColorPicture(ObjPtr pic)
#else
void UnColorPicture(pic)
ObjPtr pic;
#endif
/*Uncolors picture pic*/
{
    PicItemPtr curItem;
    ObjPtr repObj;

    curItem = ((PicPtr) pic) -> items;

    while (curItem)
    {
	curItem = curItem -> next;
    }
}

#ifdef HAVE_PROTOTYPES
void ColorItemByObject(PicItemPtr curItem,ObjPtr colorObj, Bool interpolate)
#else
void ColorItemByObject(curItem, colorObj, interpolate)
PicItemPtr curItem;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors item curItem by colorObj*/
{
    long k;
    long nVertices;
    VertexPtr *vertices;

    switch(curItem -> type)
    {
	case POLYGON:
	case POLYLINE:
	case POLYPOINT:
	    nVertices = ((PolyPtr) curItem) -> nVertices;
	    vertices = ((PolyPtr) curItem) -> vertices;
	    break;
	case POLYGONS:
	    ColorItemsByObject((PicItemPtr) ((PolysPtr) curItem) -> polygons, 
		    colorObj, interpolate);
	    return;
	    break;
	case RECTMESH:
	    {
		int i, j;
		real sample;
		int index;
		long offset;
		real position[3];


		vertices = ((RectMeshPtr) curItem) -> vertices;
		if (((RectMeshPtr) curItem) -> inCenter)
		{
		    for (i = 0; i < ((RectMeshPtr) curItem) -> xDim - 1; ++i)
		    {
			for (j = 0; j < ((RectMeshPtr) curItem) -> yDim - 1; ++j)
			{
			    offset = RECTMESHCENTER((RectMeshPtr) curItem, i, j);
			    position[0] = vertices[offset] -> position[0];
			    position[1] = vertices[offset] -> position[1];
			    position[2] = vertices[offset] -> position[2];

			    sample = SampleSpatScalar(FIELD1, FIELD2,
				3, position, interpolate);

			    index = GetRealColorIndex(sample);
			    vertices[offset] -> colorIndex = index;
			}
		    }
		    InterpRectVertices((RectMeshPtr) curItem);
		}
		else
		{
		    for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
		    {
			for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
			{
			    offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);

			    position[0] = vertices[offset] -> position[0];
			    position[1] = vertices[offset] -> position[1];
			    position[2] = vertices[offset] -> position[2];

			    sample = SampleSpatScalar(FIELD1, FIELD2,
				3, position, interpolate);

			    index = GetRealColorIndex(sample);
			    vertices[offset] -> colorIndex = index;
			}
		    }
		    InterpRectCenters((RectMeshPtr) curItem);
		}
	    }
	    return;
	    break;
	default:
	    return;
    }

    if (formsAreEqual && curItem -> flags & PF_CANONVERTICES)
    {
	int nTraversalDims;		/*Number of dimensions to traverse*/
	int whichDim;			/*Dimension counter*/
	long *traversalDims = 0;	/*Dimensions of the dataset to traverse*/
	long *index = 0;		/*Counting index*/
	long nNodes;			/*# Total nodes in the dataset*/
	int nComponents;		/*# components of data form*/
	real sample;

	/*Get the information on traversing the dataset*/
	nTraversalDims = CountTraversalDims(FORMFIELD);
	if (nTraversalDims)
	{
	    traversalDims = (long *) Alloc(sizeof(long) * nTraversalDims);
	    index = (long *) Alloc(sizeof(long) * nTraversalDims);
	}
	else
	{
	    index = (long *) Alloc(sizeof(long));
	}
	GetTraversalDims(FORMFIELD, traversalDims);

	/*Calculate the number of nodes*/
	nNodes = 1;
	for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
	{
	    if (traversalDims[whichDim] > 0)
	    {
		nNodes *= traversalDims[whichDim];
	    }
	}

	/*Color as many vertices as we have nodes*/

	/*Zero the index*/
	for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
	{
	    index[whichDim] = 0;
	}

	/*Traverse all the points*/
	k = 0;
	do
	{
	    /*Sample the location at a point*/
	    sample = SelectFieldScalar(FIELD1, index);
	    vertices[k] -> colorIndex = GetRealColorIndex(sample);
	    ++k;

	    /*Advance to next point*/
	    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
	    {
		if (traversalDims[whichDim] > 0)
		{
		    if ((++index[whichDim]) >= traversalDims[whichDim])
		    {
			index[whichDim] = 0;
		    }
		    else
		    {
			break;
		    }
		}
	    }
	} while (whichDim < nTraversalDims); /*Break is based on advance*/

	/*Free up temporary storage*/
	SAFEFREE(traversalDims);
	SAFEFREE(index);
    }
    else
    {
	real position[3];
	real sample;
	int index;
	int j;
	/*Make it colored*/
	for (k = 0; k < nVertices; ++k)
	{
    
	    position[0] = vertices[k] -> position[0];
	    position[1] = vertices[k] -> position[1];
	    position[2] = vertices[k] -> position[2];
    
	    sample = SampleSpatScalar(FIELD1, FIELD2,
				3, position, interpolate);
	
	    index = GetRealColorIndex(sample);
	    vertices[k] -> colorIndex = index;
	}
    }
}

#ifdef HAVE_PROTOTYPES
void ColorItemWithIndex(PicItemPtr curItem, int index)
#else
void ColorItemWithIndex(curItem, index)
PicItemPtr curItem;
int index;
#endif
/*Colors item curItem with index*/
{
    long k;
	long nVertices;
	VertexPtr *vertices;

	switch(curItem -> type)
	{
	    case POLYGON:
	    case POLYLINE:
		nVertices = ((PolyPtr) curItem) -> nVertices;
		vertices = ((PolyPtr) curItem) -> vertices;
		break;
	    case POLYGONS:
		ColorItemsWithIndex((PicItemPtr) ((PolysPtr) curItem) -> polygons, 
			index);
		return;
		break;
	    case RECTMESH:
		{
		    int i, j;
		    real sample;
		    long offset;

		    vertices = ((RectMeshPtr) curItem) -> vertices;
		    if (((RectMeshPtr) curItem) -> inCenter)
		    {
			for (i = 0; i < ((RectMeshPtr) curItem) -> xDim - 1; ++i)
			{
			    for (j = 0; j < ((RectMeshPtr) curItem) -> yDim - 1; ++j)
			    {
				offset = RECTMESHCENTER((RectMeshPtr) curItem, i, j);

				vertices[offset] -> colorIndex = index;
			    }
			}
			InterpRectVertices((RectMeshPtr) curItem);
		    }
		    else
		    {
			for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
			{
			    for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
			    {
				offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);

				vertices[offset] -> colorIndex = index;
			    }
			}
			InterpRectCenters((RectMeshPtr) curItem);
		    }
		}
		return;
		break;
	    default:
		return;
	}

	/*Make it colored*/
	for (k = 0; k < nVertices; ++k)
	{
	    vertices[k] -> colorIndex = index;
	}
}

#ifdef HAVE_PROTOTYPES
void ColorItemsByObject(PicItemPtr curItem,ObjPtr colorObj, Bool interpolate)
#else
void ColorItemsByObject(curItem, colorObj, interpolate)
PicItemPtr curItem;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors items at curItem by colorObj*/
{
    ObjPtr repObj;
    ObjPtr colors;
    long k;

    MakeVar(colorObj, CPALETTE);
    colors = GetPaletteVar("ColorItemsByObject", colorObj, CPALETTE);
    if (!colors) return;
    SetPalette(colors);

    while (curItem)
    {
	ColorItemByObject(curItem, colorObj, interpolate);

	curItem = curItem -> next;
    }
}

#ifdef HAVE_PROTOTYPES
void ColorItemsWithIndex(PicItemPtr curItem, int index)
#else
void ColorItemsWithIndex(curItem, index)
PicItemPtr curItem;
int index;
#endif
/*Colors items starting at curItem with index*/
{
    ObjPtr repObj;
    ObjPtr colors;
    long k;

    while (curItem)
    {
	ColorItemWithIndex(curItem, index);

	curItem = curItem -> next;
    }
}

#ifdef HAVE_PROTOTYPES
void ColorPictureByObject(ObjPtr pic,ObjPtr colorObj, Bool interpolate)
#else
void ColorPictureByObject(pic, colorObj, interpolate)
ObjPtr pic;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors picture pic by colorObj.  if interpolate, interpolates on lookup*/
{
    register VertexPtr runner;
    register VertexPtr last = NULL;
    int nTraversalDims;		/*Number of dimensions to traverse*/
    int whichDim;			/*Dimension counter*/
    long *traversalDims = 0;	/*Dimensions of the dataset to traverse*/
    long *index = 0;		/*Counting index*/
    long nNodes;			/*# Total nodes in the dataset*/
    int nComponents;		/*# components of data form*/
    real sample;
    ObjPtr palette;

    ObjPtr mainDataset = NULLOBJ;

    SetCurField(FIELD1, colorObj);
    SetCurForm(FIELD2, colorObj);

    MakeVar(colorObj, CPALETTE);
    palette = GetPaletteVar("ColorPictureByObject", colorObj, CPALETTE);
    if (palette)
    {
	SetPalette(palette);
    }

    formsAreEqual = false;

    mainDataset = GetVar(pic, MAINDATASET);

    if (mainDataset)
    {
	SetCurForm(FORMFIELD, mainDataset);
	if (IdenticalFields(FORMFIELD, FIELD2))
	{
	    formsAreEqual = true;
	}
    }

    /*Make all vertices not colored*/
    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
	runner -> flags &= ~VF_COLORED;
	runner = runner -> next;
    }

    /*Create information for the canonical group*/
    if (formsAreEqual)
    {
	nTraversalDims = CountTraversalDims(FORMFIELD);
	if (nTraversalDims)
	{
	    traversalDims = (long *) Alloc(sizeof(long) * nTraversalDims);
	    index = (long *) Alloc(sizeof(long) * nTraversalDims);
	}
	else
	{
	    index = (long *) Alloc(sizeof(long));
	}
	GetTraversalDims(FORMFIELD, traversalDims);

	/*Calculate the number of nodes*/
	nNodes = 1;
	for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
	{
	    if (traversalDims[whichDim] > 0)
	    {
		nNodes *= traversalDims[whichDim];
	    }
	}
    }

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
	if (!(runner -> flags & VF_COLORED))
	{
	    runner -> flags |= VF_COLORED;

	    if (formsAreEqual && runner -> flags & VF_FIRSTCANON)
	    {
		/*Begin a canonical group*/
		
		/*Zero the index*/
		for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
		{
		    index[whichDim] = 0;
		}
	    }

	    if (formsAreEqual && runner -> flags & VF_NEXTCANON)
	    {
		/*Advance to next point*/
		for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
		{
		    if (traversalDims[whichDim] > 0)
		    {
			if ((++index[whichDim]) >= traversalDims[whichDim])
			{
			    index[whichDim] = 0;
			}
			else
			{
			    break;
			}
		    }
		}
		if (whichDim >= nTraversalDims)
		{
		    /*Wraparound*/
		    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
		    {
			index[whichDim] = 0;
		    }
		}	
	    }

	    if (last && (runner -> flags & VF_SAMEPLACE))
	    {
		runner -> colorIndex = last -> colorIndex;
	    }
	    else if (formsAreEqual &&
		((runner -> flags & VF_NEXTCANON) ||
		 (runner -> flags & VF_FIRSTCANON)))
	    {
		/*Sample from dataform*/
		sample = SelectFieldScalar(FIELD1, index);
		runner -> colorIndex = GetRealColorIndex(sample);
	    }
	    else
	    {
		/*Sample arbitrarily*/
		sample = SampleSpatScalar(FIELD1, FIELD2,
				3, runner -> position, interpolate);
	
		runner -> colorIndex = GetRealColorIndex(sample);
	    }
	}
	last = runner;
	runner = runner -> next;
    }

    if (formsAreEqual)
    {
	/*Free up temporary storage*/
	SAFEFREE(traversalDims);
	SAFEFREE(index);
    }
}

#ifdef HAVE_PROTOTYPES
void PickCanonicalPictureVertices(ObjPtr pic)
#else
void PickCanonicalPictureVertices(pic)
ObjPtr pic;
#endif
/*Picks the vertices in a picture*/
{
    VertexPtr runner;
 
    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
	PickVertex(runner);
	bgnpoint();
	v3f(runner -> position);
	endpoint();
	runner = runner -> next;
    }
}

#ifdef HAVE_PROTOTYPES
void DeformPictureByObject(ObjPtr pic)
#else
void DeformPictureByObject(pic)
ObjPtr pic;
#endif
/*Deforms picture pic.

  ASSUMES that SetupDeformation has been done on the visualization object
*/
{
    register VertexPtr runner;

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
	runner -> flags &= ~VF_DEFORMED;
	runner = runner -> next;
    }

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
	if (!(runner -> flags & VF_DEFORMED))
	{
	    runner -> flags |= VF_DEFORMED;
	    DEFORMVERTEX(runner);
	}
	runner = runner -> next;
    }

    CalcPictureNormals(pic);
}

#ifdef HAVE_PROTOTYPES
void CalcPolyNormals(PolyPtr poly)
#else
void CalcPolyNormals(poly)
PolyPtr poly;
#endif
/*Calculates the normals in a polygon*/
{
    float x[3], nml[3], v1[3], v2[3];
    register float mag;
    register int i, p, m, n;

    /*Set up z normal, just in case*/
    nml[0] = 0.0;
    nml[1] = 0.0;
    nml[2] = 1.0;

    /*Go around, starting at 1 for robustness*/
    for (i = 1; i <= poly -> nVertices; ++i)
    {
	p = i - 1;
	m = i % (poly -> nVertices);
	n = (i + 1) % (poly -> nVertices);
	
	v1[0] = poly -> vertices[p] -> position[0] -
		poly -> vertices[m] -> position[0];
	v1[1] = poly -> vertices[p] -> position[1] -
		poly -> vertices[m] -> position[1];
	v1[2] = poly -> vertices[p] -> position[2] -
		poly -> vertices[m] -> position[2];

	v2[0] = poly -> vertices[n] -> position[0] -
		poly -> vertices[m] -> position[0];
	v2[1] = poly -> vertices[n] -> position[1] -
		poly -> vertices[m] -> position[1];
	v2[2] = poly -> vertices[n] -> position[2] -
		poly -> vertices[m] -> position[2];

	CROSS(v2, v1, x);

	mag = MAG3(x);

	if (mag > 0.0)
	{
	    /*Normal stuff*/
	    mag = 1.0 / mag;
	    nml[0] = x[0] * mag;
	    nml[1] = x[1] * mag;
	    nml[2] = x[2] * mag;
	}
	poly -> vertices[m] -> normal[0] += nml[0];
	poly -> vertices[m] -> normal[1] += nml[1];
	poly -> vertices[m] -> normal[2] += nml[2];
    }
}

#ifdef HAVE_PROTOTYPES
void CalcPolytriNormals(PolyPtr poly)
#else
void CalcPolytriNormals(poly)
PolyPtr poly;
#endif
/*Calculates the normals in a triangular mesh strip*/
{
    float x[3], v1[3], v2[3];
    register float mag;
    register int i;

    i = 0;
    while (i < poly -> nVertices - 2)
    {
	/*Do a clockwise panel*/
	v1[0] = poly -> vertices[i + 1] -> position[0] -
		poly -> vertices[i] -> position[0];
	v1[1] = poly -> vertices[i + 1] -> position[1] -
		poly -> vertices[i] -> position[1];
	v1[2] = poly -> vertices[i + 1] -> position[2] -
		poly -> vertices[i] -> position[2];

	v2[0] = poly -> vertices[i + 2] -> position[0] -
		poly -> vertices[i] -> position[0];
	v2[1] = poly -> vertices[i + 2] -> position[1] -
		poly -> vertices[i] -> position[1];
	v2[2] = poly -> vertices[i + 2] -> position[2] -
		poly -> vertices[i] -> position[2];

	CROSS(v1, v2, x);

	mag = MAG3(x);
	if (mag > 0.0)
	{
	    /*Normal stuff*/
	    mag = 1.0 / mag;
	    poly -> vertices[i] -> normal[0] += x[0] * mag;
	    poly -> vertices[i] -> normal[1] += x[1] * mag;
	    poly -> vertices[i] -> normal[2] += x[2] * mag;

	    poly -> vertices[i + 1] -> normal[0] += x[0] * mag;
	    poly -> vertices[i + 1] -> normal[1] += x[1] * mag;
	    poly -> vertices[i + 1] -> normal[2] += x[2] * mag;

	    poly -> vertices[i + 2] -> normal[0] += x[0] * mag;
	    poly -> vertices[i + 2] -> normal[1] += x[1] * mag;
	    poly -> vertices[i + 2] -> normal[2] += x[2] * mag;
	}

	++i;

	if (i < poly -> nVertices - 2)
	{
	     /*Do a counterclockwise panel*/
	    v1[0] = poly -> vertices[i + 2] -> position[0] -
		    poly -> vertices[i] -> position[0];
	    v1[1] = poly -> vertices[i + 2] -> position[1] -
		    poly -> vertices[i] -> position[1];
	    v1[2] = poly -> vertices[i + 2] -> position[2] -
		    poly -> vertices[i] -> position[2];
    
	    v2[0] = poly -> vertices[i + 1] -> position[0] -
		    poly -> vertices[i] -> position[0];
	    v2[1] = poly -> vertices[i + 1] -> position[1] -
		    poly -> vertices[i] -> position[1];
	    v2[2] = poly -> vertices[i + 1] -> position[2] -
		    poly -> vertices[i] -> position[2];
    
	    CROSS(v1, v2, x);
	    mag = MAG3(x);
	    if (mag > 0.0)
	    {
		/*Normal stuff*/
		mag = 1.0 / mag;
		poly -> vertices[i] -> normal[0] += x[0] * mag;
		poly -> vertices[i] -> normal[1] += x[1] * mag;
		poly -> vertices[i] -> normal[2] += x[2] * mag;
    
		poly -> vertices[i + 1] -> normal[0] += x[0] * mag;
		poly -> vertices[i + 1] -> normal[1] += x[1] * mag;
		poly -> vertices[i + 1] -> normal[2] += x[2] * mag;
    
		poly -> vertices[i + 2] -> normal[0] += x[0] * mag;
		poly -> vertices[i + 2] -> normal[1] += x[1] * mag;
		poly -> vertices[i + 2] -> normal[2] += x[2] * mag;
	    }
    
	    ++i;
	}
    }
}

#ifdef HAVE_PROTOTYPES
void CalcPictureNormals(ObjPtr pic)
#else
void CalcPictureNormals(pic)
ObjPtr pic;
#endif
/*Calculates the normals within pic*/
{
    PicItemPtr curItem;
    VertexPtr vertices;

    /*Go through and zero the normals*/
    vertices = ((PicPtr) pic) -> vertices;
    while (vertices)
    {
	vertices -> normal[0] = 0.0;
	vertices -> normal[1] = 0.0;
	vertices -> normal[2] = 0.0;
	vertices = vertices -> next;
    }

    curItem = ((PicPtr) pic) -> items;
    while (curItem)
    {
	switch(curItem -> type)
	{
	    case POLYGON:
	  	CalcPolyNormals((PolyPtr) curItem);
		break;
	    case POLYTRI:
		CalcPolytriNormals((PolyPtr) curItem);
		break;
	    case POLYGONS:
		{
		    PolyPtr runner;
		    runner = ((PolysPtr) ((PicPtr) pic) -> items) -> polygons;
		    while (runner)
		    {
			CalcPolyNormals(runner);
			runner = (PolyPtr) runner -> item . next;
		    }
		}
		break;
	    case RECTMESH:
		CalcRectNormals((RectMeshPtr) curItem);
		break;
	}
	curItem = curItem -> next;
    }

    /*Go through and NORMALIZE the normals*/
    vertices = ((PicPtr) pic) -> vertices;
    while (vertices)
    {
	double f;
	f = MAG3(vertices -> normal);
	if (f > 0.0)
	{
	    f = 1.0 / f;
	    vertices -> normal[0] *= f;
	    vertices -> normal[1] *= f;
	    vertices -> normal[2] *= f;
	}
	else
	{
	    vertices -> normal[2] = 1.0;
	}
	vertices = vertices -> next;
    }

}

#ifdef GRAPHICS
/*Shading models*/
int shadeModels[2][4][3] =
    {
	{{    FLAT,    FLAT, GOURAUD},
	 {    FLAT,    FLAT, GOURAUD},
	 { GOURAUD,    FLAT, GOURAUD},
	 { GOURAUD, GOURAUD, GOURAUD}},
	{{    FLAT,    FLAT, GOURAUD},
	 {    FLAT,    FLAT, GOURAUD},
	 { GOURAUD, GOURAUD, GOURAUD},
	 { GOURAUD, GOURAUD, GOURAUD}}
    };

#define PrepareToDrawC						\
	color(curBeg + curNColors - 4);
	
#define PrepareToDrawRGB						\
	C3S(curColors[curNColors - 4]);

#define PrepareToDrawCDepthCue					\
	lshaderange(curBeg + 2, curBeg + curNColors - 4, curZMin, curZMax);

#define PrepareToDrawRGBDepthCue					\
	lRGBrange(0,		\
		  0,		\
		  0,		\
		  curColors[curNColors - 4][0],		\
		  curColors[curNColors - 4][1],		\
		  curColors[curNColors - 4][2], curZMin, curZMax);


float tN[3];

/*Traversal macros to make life a little easier*/

/*Polygon traversal macro.  Decomposes simply into triangles*/
#define TRAVERSEPOLY {							\
    register VertexPtr *v1, *v2;					\
    v1 = ((PolyPtr) curItem) -> vertices;				\
	    v2 = ((PolyPtr) curItem) -> vertices + 			\
	    (((PolyPtr) curItem) -> nVertices - 1);			\
    bgntmesh();								\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
									\
    while (v2 >= v1)							\
    {									\
	/*Emit v1 and increment*/					\
	V(v1);								\
	++v1;								\
	if (v2 >= v1)							\
	{								\
	    /*Emit v2 and decrement*/					\
	    V(v2);							\
	    --v2;							\
	}								\
    }									\
    endtmesh();}

/*Polytri traversal macro.  */
#define TRAVERSEPOLYTRI {						\
    register int k, n;							\
    register VertexPtr *v1;						\
    n = ((PolyPtr) curItem) -> nVertices;				\
    v1 = ((PolyPtr) curItem) -> vertices;				\
    bgntmesh();								\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
									\
    for (k = 0; k < n; ++k)						\
    {									\
	/*Emit v1 and increment*/					\
	V(v1); ++v1;							\
    }									\
    endtmesh();}

/*Planar traversal macro*/
#define TRAVERSEPLANARPOLY {						\
    register VertexPtr *v1;						\
    register int k;							\
    v1 = ((PolyPtr) curItem) -> vertices;				\
									\
    bgnpolygon();							\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)		\
    {									\
	/*Emit v1 and increment*/					\
	V(v1);								\
	++v1;								\
    }									\
    endpolygon();}

/*Polyline traversal macro*/
#define TRAVERSEPOLYLINE {						\
    register VertexPtr *v1;						\
    register int k;							\
    v1 = ((PolyPtr) curItem) -> vertices;				\
									\
    bgnline();								\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)		\
    {									\
	/*Emit v1 and increment*/					\
	V(v1);								\
	++v1;								\
    }									\
    endline();}

/*Polypoint traversal macro*/
#define TRAVERSEPOLYPOINT {						\
    register VertexPtr *v1;						\
    register int k;							\
    v1 = ((PolyPtr) curItem) -> vertices;				\
									\
    bgnpoint();								\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)		\
    {									\
	/*Emit v1 and increment*/					\
	if (0 == ((*v1) -> flags & VF_DONTDRAW)) {V(v1)}			\
	++v1;								\
    }									\
    endpoint();								\
}

/*Closed polyline traversal macro*/
#define TRAVERSECLOSEDPOLYLINE {						\
    register VertexPtr *v1;						\
    register int k;							\
    v1 = ((PolyPtr) curItem) -> vertices;				\
									\
    bgnclosedline();								\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)		\
    {									\
	/*Emit v1 and increment*/					\
	V(v1);								\
	++v1;								\
    }									\
    endclosedline();}

/*Closed polytri line traversal macro*/
#define TRAVERSECLOSEDPOLYTRI {						\
    register VertexPtr *v1;						\
    register int k;							\
    v1 = ((PolyPtr) curItem) -> vertices;				\
									\
    bgnclosedline();							\
    P(v1, ((PolyPtr) curItem) -> nVertices);				\
    V(v1);								\
    k = 1; ++v1;							\
    while (k < ((PolyPtr) curItem) -> nVertices)			\
    {									\
	V(v1);								\
	v1 += 2;							\
	k += 2;								\
    } --k; --v1;							\
    if (k < ((PolyPtr) curItem) -> nVertices) {k += 2; v1 += 2;}	\
    do									\
    {									\
	k -= 2;								\
	v1 -= 2;							\
	V(v1);								\
    } while (k > 0);							\
    endclosedline();}

/*Dots traversal macro*/
#define TRAVERSEDOTS	{						\
    VertexPtr vertices;							\
    vertices = ((PicPtr) pic) -> vertices;				\
    bgnpoint();								\
    P(&vertices, 1);							\
    while (vertices)							\
    {									\
	if (0 == (vertices -> flags & VF_DONTDRAW)) {V(&vertices)}	\
	vertices = vertices -> next;					\
    }									\
    endpoint();}

/*Wire frame rectangle traversal macro*/
#define TRAVERSEWRECTMESH						\
    {									\
    register int i, j;							\
    int xDim, yDim;							\
    register VertexPtr *vertexPtr;					\
    int nX, nY;								\
									\
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;			\
    xDim = ((RectMeshPtr) curItem) -> xDim;				\
    yDim = ((RectMeshPtr) curItem) -> yDim;				\
    nX = (2 * yDim - 1);						\
    nY = 2;								\
									\
    /*Draw the fast moving axes*/					\
    for (i = 0; i < xDim; ++i)						\
    {									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices +		\
		RECTMESHVERTEX(((RectMeshPtr) curItem), i, 0);		\
	bgnline();							\
	P(vertexPtr, 1);						\
	for (j = 0; j < yDim; ++j)					\
	{								\
	    V(vertexPtr);						\
	    vertexPtr += nY;						\
	}								\
	endline();							\
    }									\
									\
    /*Draw the slow moving axes*/					\
    for (j = 0; j < yDim; ++j)						\
    {									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices +		\
			RECTMESHVERTEX(((RectMeshPtr) curItem), 0, j);	\
	bgnline();							\
	P(vertexPtr, 1);						\
	for (i = 0; i < xDim; ++i)					\
	{								\
	    V(vertexPtr);						\
	    vertexPtr += nX;						\
	}								\
	endline();							\
    }									\
    }

/*Wire frame monocolor rectangle traversal macro*/
#define TRAVERSEWMONOCOLORRECTMESH						\
    {									\
    register int i, j;							\
    int xDim, yDim;							\
    register VertexPtr *vertexPtr;					\
    int nX, nY;								\
									\
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;			\
    xDim = ((RectMeshPtr) curItem) -> xDim;				\
    yDim = ((RectMeshPtr) curItem) -> yDim;				\
    nX = (2 * yDim - 1);						\
    nY = 2;								\
									\
    /*Draw the fast moving axes*/					\
    for (i = 0; i < xDim; ++i)						\
    {									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices +		\
		RECTMESHVERTEX(((RectMeshPtr) curItem), i, 0);		\
	for (j = 0; j < yDim - 1; ++j)					\
	{								\
	    bgnline();							\
	    P1(vertexPtr, 1);						\
	    V(vertexPtr);						\
	    vertexPtr += nY;						\
	    V(vertexPtr);						\
	    endline();							\
	}								\
    }									\
									\
    /*Draw the slow moving axes*/					\
    for (j = 0; j < yDim; ++j)						\
    {									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices +		\
			RECTMESHVERTEX(((RectMeshPtr) curItem), 0, j);	\
	for (i = 0; i < xDim - 1; ++i)					\
	{								\
	    bgnline();							\
	    P2(vertexPtr, 1);						\
	    V(vertexPtr);						\
	    vertexPtr += nX;						\
	    V(vertexPtr);						\
	    endline();							\
	}								\
    }									\
    }


/*Rectangular mesh traversal macro*/
#define TRAVERSERECTMESH						\
    {									\
	register int i, j;						\
	int xDim, yDim;							\
	register VertexPtr *vertexPtr;					\
	int nX, nY;							\
									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices;		\
	xDim = ((RectMeshPtr) curItem) -> xDim;				\
	yDim = ((RectMeshPtr) curItem) -> yDim;				\
	nY = 2;								\
	nX = (yDim * 2 - 1);						\
	for (i = 0; i < xDim - 1; ++i)					\
	{								\
	    bgntmesh();						\
	    P(vertexPtr, 1);					\
	    V(vertexPtr);					\
	    vertexPtr += nX;					\
	    V(vertexPtr);					\
	    for (j = 1; j < yDim; ++j)				\
	    {							\
		vertexPtr += nY - nX;				\
		V(vertexPtr);					\
		vertexPtr += nX;				\
		V(vertexPtr);					\
	    }							\
	    vertexPtr -= nX - 1;				\
	    endtmesh();						\
	}								\
    }

/*Monocolor Rectangular mesh traversal macro*/
#define TRAVERSEMCRECTMESH						\
    {									\
	register int i, j;						\
	int xDim, yDim;							\
	register VertexPtr *vertexPtr;					\
	int nX, nY;							\
									\
	vertexPtr = ((RectMeshPtr) curItem) -> vertices;		\
	xDim = ((RectMeshPtr) curItem) -> xDim;				\
	yDim = ((RectMeshPtr) curItem) -> yDim;				\
	nY = 2;								\
	nX = (yDim * 2 - 1);						\
	for (i = 0; i < xDim - 1; ++i)					\
	{								\
	    for (j = 0; j < yDim - 1; ++j)				\
	    {								\
		bgntmesh();						\
		P(vertexPtr + 1);						\
		V(vertexPtr);						\
		vertexPtr += nX;					\
		V(vertexPtr);						\
		vertexPtr += nY - nX;					\
		V(vertexPtr);						\
		vertexPtr += nX;					\
		V(vertexPtr);						\
		endtmesh();						\
		vertexPtr -= nX;					\
	    }								\
	    vertexPtr += nY - 1;  /*-1 for missing center*/		\
	}								\
    }
#endif

static void DrawCPolygon(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWCPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYLINE
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolytri(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYTRI
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWCPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYTRI
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolyline(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polyline*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYLINE
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolypoint(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polyline*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYPOINT
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYTRI
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWRGBPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYLINE
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWRGBPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYTRI
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolyline(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYLINE
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolypoint(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYPOINT
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}
static void DrawCPPolygon(curItem)
PicItemPtr curItem;
/*Draws a color map shaded planarpolygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPPolygon(curItem)
PicItemPtr curItem;
/*Draws an RGB shaded planar polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}


static void DrawWCRectMesh(curItem)
PicItemPtr curItem;
/*Draws a color map wire frame rect mesh*/
{
#ifdef GRAPHICS
	if (curLightShading == DEPTHCUELIGHT)
	{
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PCDN
#define V VCDN
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P1 PCDM
#define P2 PCDMSKIP
#define V VCDM
			TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PCDS
#define V VCDS
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		}
	}
	else
	{
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PCNN
#define V VCNN
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P1 PCNM
#define P2 PCNMSKIP
#define V VCNM
			TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PCNS
#define V VCNS
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		}
	}
#endif
}

static void DrawWRGBRectMesh(curItem)
PicItemPtr curItem;
/*Draws an RGB wire frame rectangular mesh*/
{
#ifdef GRAPHICS
	if (curLightShading == DEPTHCUELIGHT)
	{
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PRDN
#define V VRDN
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P1 PRDM
#define P2 PRDMSKIP
#define V VRDM
			TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PRDS
#define V VRDS
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		}
	}
	else
	{
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PRNN
#define V VRNN
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P1 PRNM
#define P2 PRNMSKIP
#define V VRNM
			TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PRNS
#define V VRNS
			TRAVERSEWRECTMESH;
#undef P
#undef V
			break;
		}
	}
#endif
}

static void DrawCRectMesh(curItem)
PicItemPtr curItem;
/*Draws a rectangular mesh in color map mode*/
{
#ifdef GRAPHICS
	switch (curLightShading)
	{
	    case NOLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PCNN
#define V VCNN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPCNM
#define V VCNM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PCNS
#define V VCNS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	    case MONOLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PCMN
#define V VCMN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPCNM
#define V VCMM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PCMS
#define V VCMS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	    case SMOOTHLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PCSN
#define V VCSN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPCNM
#define V VCSM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PCSS
#define V VCSS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	}
#endif
}

int yaga = 0;

static void DrawRGBRectMesh(curItem)
PicItemPtr curItem;
/*Draws a rectangular mesh in RGB mode*/
{
#ifdef GRAPHICS
	switch (curLightShading)
	{
	    case NOLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PRNN
#define V VRNN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPRNM
#define V VRNM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PRNS
#define V VRNS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	    case MONOLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PRMN
#define V VRMN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPRNM
#define V VRMM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PRMS
#define V VRMS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	    case SMOOTHLIGHT:
		switch (curColorShading)
		{
		    case NOCOLOR:
#define P PRSN
#define V VRSN
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		    case MONOCOLOR:
#define P MCRPRNM
#define V VRSM
			TRAVERSEMCRECTMESH;
#undef P
#undef V
			break;
		    case SMOOTHCOLOR:
#define P PRSS
#define V VRSS
			TRAVERSERECTMESH;
#undef P
#undef V
			break;
		}
		break;
	}
#endif
}

void DrawCSphere(curItem)
PicItemPtr curItem;
/*Draws a sphere in cmap mode*/
{
#ifdef GRAPHICS
			{
			    Sphere curSphere;
			    float vx[3], vy[3], vz[3]; /*Vertices stuck out on those axes*/
			    
			    curSphere = *((SpherePtr) curItem);
	
	        	    if (!curIsTransparent)
	        	    {
	       			backface(TRUE);
	        	    }
			    /*Draw the eight quadrants of the sphere*/
			    /* +x +y +z */
			    vx[0] = 1.0;
			    vx[1] = 0.0;
			    vx[2] = 0.0;
			    vy[0] = 0.0;
			    vy[1] = 1.0;
			    vy[2] = 0.0;
			    vz[0] = 0.0;
			    vz[1] = 0.0;
			    vz[2] = 1.0;
	
			    curRadius = curSphere. radius;
			    curCenter[0] = curSphere . centerVertex -> position[0];
			    curCenter[1] = curSphere . centerVertex -> position[1];
			    curCenter[2] = curSphere . centerVertex -> position[2];
			    DrawQuadrant(vx, vy, vz, SPHERESUB);
	
			    /* +y -x +z*/
			    vx[0] = -1.0;
			    DrawQuadrant(vy, vx, vz, SPHERESUB);
	
			    /* -x -y -z*/
			    vy[1] = -1.0;
			    DrawQuadrant(vx, vy, vz, SPHERESUB);
	
			    /* -y +x +z*/
			    vx[0] = 1.0;
			    DrawQuadrant(vy, vx, vz, SPHERESUB);
	
			    /* +x -y -z*/
			    vz[2] = -1.0;
			    DrawQuadrant(vx, vy, vz, SPHERESUB);
	
			    /* +y +x -z*/
			    vy[1] = 1.0;
			    DrawQuadrant(vy, vx, vz, SPHERESUB);
	
			    /* -x +y -z*/
			    vx[0] = -1.0;
			    DrawQuadrant(vx, vy, vz, SPHERESUB);
	
			    /* -y, -x, -z*/
			    vy[1] = -1.0;
			    DrawQuadrant(vy, vx, vz, SPHERESUB);
	
	        	    if (!curIsTransparent)
	        	    {
	       			backface(FALSE);
	        	    }
			}
#endif
}

void DrawSFrustum(curItem)
PicItemPtr curItem;
/*Draws a skeletal frustum*/
{
#ifdef GRAPHICS
    bgnline();
    v3f(((FrustumPtr) curItem) -> end1);
    v3f(((FrustumPtr) curItem) -> end2);
    endline();
#endif
}

void DrawCFrustum(curItem)
PicItemPtr curItem;
/*Draws a color map frustum*/
{
#ifdef GRAPHICS
    Frustum curFrustum;
    float a[3];		/*Axial unit vector*/
    float r[CYLSIDES][3];
    double length;	/*Length of cylinder*/
    double nf;		/*Normalization factor*/
			/*Radial units around cylinder*/
    int end;		/*Counter for which end*/
    int b;		/*Counter around radius*/
    int d;		/*Delta for subdivision*/
    int k;		/*Random counter*/

    if (!curIsTransparent)
    {
	backface(TRUE);
    }

    /*Copy the frustum into curFrustum for speed*/
    curFrustum = *((FrustumPtr) curItem);

    /*First make an axial unit vector*/
    a[0] = curFrustum . end2[0] -
	   curFrustum . end1[0];
    a[1] = curFrustum . end2[1] -
	   curFrustum . end1[1];
    a[2] = curFrustum . end2[2] -
	   curFrustum . end1[2];
    length = sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
    nf = 1.0 / length;
    a[0] *= nf; a[1] *= nf; a[2] *= nf;

    /*See if it's nearly colinear with i*/
    if (ABS(a[0]) > 0.8)
    {
	/*It is, so cross by j to get first r*/
	r[0][0] = -a[2];
	r[0][1] = 0.0;
	r[0][2] = a[0];
    }
    else
    {
	/*It isn't, so cross by i to get first r*/
	r[0][0] = 0.0;
	r[0][1] = a[2];
	r[0][2] = -a[1];
    }
    
    /*Normalize r*/
    length = sqrt(r[0][0]*r[0][0] + r[0][1]*r[0][1] + r[0][2]*r[0][2]);
    nf = 1.0 / length;
    r[0][0] *= nf; r[0][1] *= nf, r[0][2] *= nf;

    /*Cross a with first radial unit to get orthogonal unit*/
    CROSS(a, r[0], r[CYLSIDES / 4]);

    /*Fill out point radii 3 and 4 by reversing 1 and 2*/
    r[CYLSIDES / 2][0] = -r[0][0];
    r[CYLSIDES / 2][1] = -r[0][1];
    r[CYLSIDES / 2][2] = -r[0][2];

    r[CYLSIDES / 4 * 3][0] = -r[CYLSIDES / 4][0];
    r[CYLSIDES / 4 * 3][1] = -r[CYLSIDES / 4][1];
    r[CYLSIDES / 4 * 3][2] = -r[CYLSIDES / 4][2];

    /*Subdivide the sides of the cylinder*/
    k = 0;
    for (d = CYLSIDES / 4; d > 1; d /= 2)
    {
	register double cef_2;	/*Chord extension factor/2*/
	
	cef_2 = cefs_2[k];
	++k;
	
	for (b = 0; b < CYLSIDES; b += d)
	{
	    register int m, e;	/*Midpoint and endpoint*/

	    /*Interpolate center of chord*/
	    e = (b + d) & (CYLSIDES - 1);	/*Clip circular*/
	    m = b + d / 2;
	    r[m][0] = (r[b][0] + r[e][0]) * cef_2;
	    r[m][1] = (r[b][1] + r[e][1]) * cef_2;
	    r[m][2] = (r[b][2] + r[e][2]) * cef_2;
	}
    }

    /*Draw the cylinder*/
    
#ifndef PSYCHEDELIC
    if (rgbp == false)
    {
	float v[3];
	
	bgntmesh();
	SetRealColor(r[CYLSIDES - 1][0]);
	v[0] = curFrustum . rad2 * r[CYLSIDES - 1][0] + 
	       curFrustum . end2[0];
	v[1] = curFrustum . rad2 * r[CYLSIDES - 1][1] + 
	       curFrustum . end2[1];
	v[2] = curFrustum . rad2 * r[CYLSIDES - 1][2] + 
	       curFrustum . end2[2];
	v3f(v);
	v[0] = curFrustum . rad1 * r[CYLSIDES - 1][0] + 
	       curFrustum . end1[0];
	v[1] = curFrustum . rad1 * r[CYLSIDES - 1][1] + 
	       curFrustum . end1[1];
	v[2] = curFrustum . rad1 * r[CYLSIDES - 1][2] + 
	       curFrustum . end1[2];
	v3f(v);

	for (k = 0; k < CYLSIDES; ++k)
	{
	    SetRealColor(r[k][0]);
	    v[0] = curFrustum . rad2 * r[k][0] + 
		   curFrustum . end2[0];
	    v[1] = curFrustum . rad2 * r[k][1] + 
		   curFrustum . end2[1];
	    v[2] = curFrustum . rad2 * r[k][2] + 
		   curFrustum . end2[2];
	    v3f(v);
	    v[0] = curFrustum . rad1 * r[k][0] + 
		   curFrustum . end1[0];
	    v[1] = curFrustum . rad1 * r[k][1] + 
		   curFrustum . end1[1];
	    v[2] = curFrustum . rad1 * r[k][2] + 
		   curFrustum . end1[2];
	    v3f(v);
	}
	endtmesh();

	/*Do top cap*/
	SetRealColor(a[0]);
	bgnpolygon();
	for (k = 0; k < CYLSIDES; ++k)
	{
	    v[0] = curFrustum . rad2 * r[k][0] + 
		   curFrustum . end2[0];
	    v[1] = curFrustum . rad2 * r[k][1] + 
		   curFrustum . end2[1];
	    v[2] = curFrustum . rad2 * r[k][2] + 
		   curFrustum . end2[2];
	    v3f(v);
	}
	endpolygon();

	/*Do bottom cap*/
	SetRealColor(-a[0]);
	bgnpolygon();
	for (k = CYLSIDES - 1; k >= 0; --k)
	{
	    v[0] = curFrustum . rad1 * r[k][0] + 
		   curFrustum . end1[0];
	    v[1] = curFrustum . rad1 * r[k][1] + 
		   curFrustum . end1[1];
	    v[2] = curFrustum . rad1 * r[k][2] + 
		   curFrustum . end1[2];
	    v3f(v);
	}
	endpolygon();
    }
    else
#endif
    {
	float v[3];

	bgntmesh();
	N3F(r[CYLSIDES - 1]);
	v[0] = curFrustum . rad2 * r[CYLSIDES - 1][0] + 
	       curFrustum . end2[0];
	v[1] = curFrustum . rad2 * r[CYLSIDES - 1][1] + 
	       curFrustum . end2[1];
	v[2] = curFrustum . rad2 * r[CYLSIDES - 1][2] + 
	       curFrustum . end2[2];
	v3f(v);
	v[0] = curFrustum . rad1 * r[CYLSIDES - 1][0] + 
	       curFrustum . end1[0];
	v[1] = curFrustum . rad1 * r[CYLSIDES - 1][1] + 
	       curFrustum . end1[1];
	v[2] = curFrustum . rad1 * r[CYLSIDES - 1][2] + 
	       curFrustum . end1[2];
	v3f(v);

	for (k = 0; k < CYLSIDES; ++k)
	{
	    N3F(r[k]);
	    v[0] = curFrustum . rad2 * r[k][0] + 
		   curFrustum . end2[0];
	    v[1] = curFrustum . rad2 * r[k][1] + 
		   curFrustum . end2[1];
	    v[2] = curFrustum . rad2 * r[k][2] + 
		   curFrustum . end2[2];
	    v3f(v);
	    v[0] = curFrustum . rad1 * r[k][0] + 
		   curFrustum . end1[0];
	    v[1] = curFrustum . rad1 * r[k][1] + 
		   curFrustum . end1[1];
	    v[2] = curFrustum . rad1 * r[k][2] + 
		   curFrustum . end1[2];
	    v3f(v);
	}
	endtmesh();

	/*Do top cap*/
	N3F(a);
	bgnpolygon();
	for (k = 0; k < CYLSIDES; ++k)
	{
	    v[0] = curFrustum . rad2 * r[k][0] + 
		   curFrustum . end2[0];
	    v[1] = curFrustum . rad2 * r[k][1] + 
		   curFrustum . end2[1];
	    v[2] = curFrustum . rad2 * r[k][2] + 
		   curFrustum . end2[2];
	    v3f(v);
	}
	endpolygon();

	/*Do bottom cap*/
	a[0] = -a[0];
	a[1] = -a[1];
	a[2] = -a[2];
	N3F(a);
	bgnpolygon();
	for (k = CYLSIDES - 1; k >= 0; --k)
	{
	    v[0] = curFrustum . rad1 * r[k][0] + 
		   curFrustum . end1[0];
	    v[1] = curFrustum . rad1 * r[k][1] + 
		   curFrustum . end1[1];
	    v[2] = curFrustum . rad1 * r[k][2] + 
		   curFrustum . end1[2];
	    v3f(v);
	}
	endpolygon();
    }
    if (!curIsTransparent)
    {
	backface(FALSE);
    }
#endif
}

void DrawCPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons in cmode*/
{
#ifdef GRAPHICS
    PolyPtr curItem;

    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
	DrawCPolygon(curItem);
	curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

void DrawRGBPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons in RGB mode*/
{
#ifdef GRAPHICS
    PolyPtr curItem;

    if (!curIsTransparent && ((PolysPtr) item) -> enclosed && !orderReversed)	
    {
	backface(TRUE);
    }
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
	DrawRGBPolygon(curItem);
	curItem = (PolyPtr) curItem -> item . next;
    }
    if (!curIsTransparent && ((PolysPtr) item) -> enclosed && !orderReversed)	
    {
	backface(FALSE);
    }
#endif
}

void DrawWCPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons as wire frames*/
{
#ifdef GRAPHICS
    PolyPtr curItem;
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
	DrawWCPolygon(curItem);
	curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

void DrawWRGBPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons as wire frames*/
{
#ifdef GRAPHICS
    PolyPtr curItem;
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
	DrawWRGBPolygon(curItem);
	curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

#ifdef GRAPHICS
typedef void (* DrawFunc)();

/*Drawing Table*/

DrawFunc drawFuncs[NPICITEMTYPES][DQ_FULL + 1][2] =
    {
    /*    Wireframe C     Wireframe RGB       Full C          Full RGB           */
	{{DrawWCPolygon,  DrawWRGBPolygon }, {DrawCPolygon,   DrawRGBPolygon }},
	{{DrawSFrustum,   DrawSFrustum    }, {DrawCFrustum,   DrawCFrustum   }},
	{{DrawCSphere,    DrawCSphere     }, {DrawCSphere,    DrawCSphere    }},
	{{DrawWCRectMesh, DrawWRGBRectMesh}, {DrawCRectMesh,  DrawRGBRectMesh}},
	{{DrawWCPolygons, DrawWRGBPolygons}, {DrawCPolygons,  DrawRGBPolygons}},
	{{DrawCPolyline,  DrawRGBPolyline }, {DrawCPolyline,  DrawRGBPolyline}},
	{{DrawCPolypoint, DrawRGBPolypoint}, {DrawCPolypoint, DrawRGBPolypoint}},
	{{DrawWCPolytri,  DrawWRGBPolytri }, {DrawCPolytri,   DrawRGBPolytri }}
    };

#endif

#ifdef HAVE_PROTOTYPES
void PrepareToDraw(Bool isTransparent, int lightShading, int colorShading, ObjPtr colors)
#else
void PrepareToDraw(isTransparent, lightShading, colorShading, colors)
Bool isTransparent;
int lightShading;
int colorShading;
ObjPtr colors;
#endif
/*Prepares to do some drawing*/
{
    /*Store light and color shading*/
    curLightShading = lightShading;
    curColorShading = colorShading;

    if ((drawingQuality != DQ_FULL && curLightShading != DEPTHCUELIGHT)
	|| overrideColor)
    {
	curLightShading = NOLIGHT;
    }

    /*And the transparency*/
    curIsTransparent = isTransparent;
    SetPalette(colors);

    /*Set the color mode*/
    if (rgbp == false)
    {
	curColorMode = CMODECMAP;
    }
    else
    {
	curColorMode = CMODERGB;
    }

#ifdef GRAPHICS
    /*Set shading model*/
    shademodel(shadeModels[curColorMode][curLightShading][curColorShading]);

    if ((curLightShading == DEPTHCUELIGHT) && hasDepthCue)
    {
	/*Set depth cueing*/
	depthcue(TRUE);
    }
    else if (overrideColor == false)
#endif
    {
	/*Set the default color for skeleton*/
	if (curSpace)
	{
	    ObjPtr panel, color, var;
	    real rgb[3];
    
	    panel = GetVar(curSpace, BACKPANEL);
	    if (panel)
	    {
		var = GetVar(panel, BACKGROUND);
		if (var && IsRealArray(var) && RANK(var) == 1 && DIMS(var)[0] == 3)
		{
		    Array2CArray(rgb, var);
		}
		else if (var && IsInt(var))
		{
		    rgb[0] = uiColors[GetInt(var)][0] / 255.0;
		    rgb[1] = uiColors[GetInt(var)][1] / 255.0;
		    rgb[2] = uiColors[GetInt(var)][2] / 255.0;
		}
		else
		{
		    rgb[0] = rgb[1] = rgb[2] = 0.0;
		}
	    }
	    if (rgb[0] + rgb[1] + rgb[2] > 1.5)
	    {
		SetUIColor(UIBLACK);
	    }
	    else
	    {
		SetUIColor(UIWHITE);
	    }
	}
    }
}

void DoneDrawing()
/*Done drawing a picture*/
{
#ifdef GRAPHICS
    if (hasDepthCue)
    {
	depthcue(FALSE);
    }
#endif
}

void ZeroPictureTime()
/*Zeroes the picture time*/
{
    pictureTime = 0;
}

long GetPictureTime()
/*Gets the picture time*/
{
    return pictureTime;
}

#ifdef HAVE_PROTOTYPES
void DrawPicture(ObjPtr pic, Bool isTransparent, int lightShading, int colorShading, ObjPtr colors)
#else
void DrawPicture(pic, isTransparent, lightShading, colorShading, colors)
ObjPtr pic;
Bool isTransparent;
int lightShading;
int colorShading;
ObjPtr colors;
#endif
/*Draws picture pic.  isTransparent if transparent*/
{
#ifdef GRAPHICS
    PicItemPtr curItem;
    ObjPtr repObj;
    DrawFunc drawFunc;
    int qualIndex;
    long beginTime;
    struct tms buffer;

    beginTime = times(&buffer);

    /*Calculate the quality index*/
    qualIndex = drawingQuality;

    /*Make the picture deformed*/
    MakeVar(pic, PICDEFORMED);

    /*Prepare to draw*/
    PrepareToDraw(isTransparent, lightShading, colorShading, colors);

    switch (curColorMode)
    {
	case CMODECMAP:
		if (lightShading == DEPTHCUELIGHT)
		{
		    PrepareToDrawCDepthCue;
		}
		else
		{
		    PrepareToDrawC;
		}
		break;
	case CMODERGB:
		if (lightShading == DEPTHCUELIGHT)
		{
		    PrepareToDrawRGBDepthCue;
		}
		else
		{
		    PrepareToDrawRGB;
		}
		break;
    }
    if (qualIndex == DQ_DOTS)
    {
	switch (curColorMode)
	{
	    case CMODECMAP:
#define TRAVERSE TRAVERSEDOTS
#include "ScianStdCDraw.h"
#undef TRAVERSE
		break;
	    case CMODERGB:
#define TRAVERSE TRAVERSEDOTS
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
		break;
	}
    }
    else
    {
	curItem = ((PicPtr) pic) -> items;
    
	switch (curColorMode)
	{
	    case CMODECMAP:
		while (curItem)
		{
		    int proximity;
		    int k;
		    if (curItem -> type < 0 || curItem -> type >= NPICITEMTYPES)
		    {
			char errStr[256];
			sprintf(errStr, "Item type %d not valid", curItem -> type);
			ReportError("DrawPicture", errStr);
			continue;
		    }
		    proximity = curItem -> proximity;
		    if (proximity > 5 || proximity < -5)
		    {
			ReportError("DrawPicture", "Bad proximity");
			continue;
		    }
		    if (proximity > 0)
		    {
			for (k = 0; k < proximity; ++k)
			{
			    NudgeCloser();
			}
		    }
		    else if (proximity < 0)
		    {
			for (k = 0; k > proximity; --k)
			{
			    NudgeFarther();
			}
		    }
		    drawFunc = drawFuncs[curItem -> type][qualIndex][0];
		    if (drawFunc)
		    {
			(*drawFunc)(curItem);
		    }
		    if (proximity > 0)
		    {
			for (k = 0; k < proximity; ++k)
			{
			    NudgeFarther();
			}
		    }
		    else if (proximity < 0)
		    {
			for (k = 0; k > proximity; --k)
			{
			    NudgeCloser();
			}
		    }
		    curItem = curItem -> next;
		}
		break;
	    case CMODERGB:
		if (overrideColor)
		{
		    lmcolor(LMC_NULL);
		}
		else if (lightShading == NOLIGHT)
		{
		    lmcolor(LMC_COLOR);
		}
		else
		{
		    lmcolor(LMC_DIFFUSE);
		}
    
		while (curItem)
		{
		    int proximity;
		    int k;
		    if (curItem -> type < 0 || curItem -> type >= NPICITEMTYPES)
		    {
			char errStr[256];
			sprintf(errStr, "Item type %d not valid", curItem -> type);
			ReportError("DrawPicture", errStr);
			continue;
		    }
		    proximity = curItem -> proximity;
		    if (proximity > 5 || proximity < -5)
		    {
			ReportError("DrawPicture", "Bad proximity");
			continue;
		    }
		    if (proximity > 0)
		    {
			for (k = 0; k < proximity; ++k)
			{
			    NudgeCloser();
			}
		    }
		    else if (proximity < 0)
		    {
			for (k = 0; k > proximity; --k)
			{
			    NudgeFarther();
			}
		    }
		    drawFunc = drawFuncs[curItem -> type][qualIndex][1];
		    (*drawFunc)(curItem);
		    if (proximity > 0)
		    {
			for (k = 0; k < proximity; ++k)
			{
			    NudgeFarther();
			}
		    }
		    else if (proximity < 0)
		    {
			for (k = 0; k > proximity; --k)
			{
			    NudgeCloser();
			}
		    }
		    curItem = curItem -> next;
		}
		break;
	}
    }

    DoneDrawing();

    pictureTime += (times(&buffer) - beginTime);
#endif
}

#ifdef HAVE_PROTOTYPES
void SavePictureGeometry(char *picName, ObjPtr pic, ObjPtr colors)
#else
void SavePictureGeometry(picName, pic, colors)
char *picName;
ObjPtr pic;
ObjPtr colors;
#endif
/*Saves the geometry of picture pic into a file derived from picName*/
{
    PicItemPtr curItem;
    ObjPtr repObj;
    FILE *outFile;
    char fileName[256];
    int curColorIndex;
    float curNormal[3];
    int k;

    printf("Saving picture for %s\n", picName);
    strcpy(fileName, picName);
    strcat(fileName, ".c");

    outFile = fopen(fileName, "w");
    if (!outFile)
    {
	perror(fileName);
	return;
    }

    /*Make the picture deformed*/
    MakeVar(pic, PICDEFORMED);

    /*Prepare to draw*/
    SetPalette(colors);

    curItem = ((PicPtr) pic) -> items;

    curColorIndex = -1;
    curNormal[0] = curNormal[1] = curNormal[2] = 0.0;

    fprintf(outFile, "    o = genobj();\n");
    fprintf(outFile, "    makeobj(o);\n");

    while (curItem)
    {
	switch(curItem -> type)
	{
	    case POLYGON:
		{
		    PolyPtr poly;
		    VertexPtr vertex;

		    poly = (PolyPtr) curItem;

		    fprintf(outFile, "    bgnpolygon();\n");
		    for (k = 0; k < poly -> nVertices; ++k)
		    {
			vertex = poly -> vertices[k];

			if (vertex -> colorIndex != curColorIndex)
			{
			    curColorIndex = vertex -> colorIndex;

			    fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
			    fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
			    fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
			    fprintf(outFile, "    c3s(c);\n");
			}

			if (vertex -> normal[0] != curNormal[0] ||
			    vertex -> normal[1] != curNormal[1] ||
			    vertex -> normal[2] != curNormal[2])
			{
			    curNormal[0] = vertex -> normal[0];
			    curNormal[1] = vertex -> normal[1];
			    curNormal[2] = vertex -> normal[2];

			    fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
			    fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
			    fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
			    fprintf(outFile, "    n3f(n);\n");
			}

			fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			fprintf(outFile, "    v3f(v);\n");
		    }
		    fprintf(outFile, "    endpolygon();\n\n");
		}
		break;
	    case RECTMESH:
		{
		    register long i, j, offset, iDim, jDim, nX;
		    register VertexPtr *vertices;
		    register VertexPtr vertex;
		    float position[3];		/*Position of center vertex*/
		    int colorIndex;		/*Color within color table*/
		    RectMeshPtr rectMesh;

		    rectMesh = (RectMeshPtr) curItem;

		    iDim = rectMesh -> xDim;
		    jDim = rectMesh -> yDim;
		    vertices = rectMesh -> vertices;
		    nX = jDim * 2 - 1;
		
		    for (i = 0; i < iDim - 1; ++i)
		    {
			for (j = 0; j < jDim - 1; ++j)
			{
			    fprintf(outFile, "    bgnpolygon();\n");

			    /*Lower left corner*/
			    offset = RECTMESHVERTEX(rectMesh, i, j);
			    vertex = vertices[offset];
    
			    if (vertex -> colorIndex != curColorIndex)
			    {
				curColorIndex = vertex -> colorIndex;
    
				fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
				fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
				fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
				fprintf(outFile, "    c3s(c);\n");
			    }
    
			    if (vertex -> normal[0] != curNormal[0] ||
				vertex -> normal[1] != curNormal[1] ||
				vertex -> normal[2] != curNormal[2])
			    {
				curNormal[0] = vertex -> normal[0];
				curNormal[1] = vertex -> normal[1];
				curNormal[2] = vertex -> normal[2];
    
				fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
				fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
				fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
				fprintf(outFile, "    n3f(n);\n");
			    }
    
			    fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			    fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			    fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			    fprintf(outFile, "    v3f(v);\n");

			    /*Lower right corner*/
			    offset = RECTMESHVERTEX(rectMesh, i + 1, j);
			    vertex = vertices[offset];
    
			    if (vertex -> colorIndex != curColorIndex)
			    {
				curColorIndex = vertex -> colorIndex;
    
				fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
				fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
				fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
				fprintf(outFile, "    c3s(c);\n");
			    }
    
			    if (vertex -> normal[0] != curNormal[0] ||
				vertex -> normal[1] != curNormal[1] ||
				vertex -> normal[2] != curNormal[2])
			    {
				curNormal[0] = vertex -> normal[0];
				curNormal[1] = vertex -> normal[1];
				curNormal[2] = vertex -> normal[2];
    
				fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
				fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
				fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
				fprintf(outFile, "    n3f(n);\n");
			    }
    
			    fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			    fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			    fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			    fprintf(outFile, "    v3f(v);\n");
		
			    /*Upper right corner*/
			    offset = RECTMESHVERTEX(rectMesh, i + 1, j + 1);
			    vertex = vertices[offset];
    
			    if (vertex -> colorIndex != curColorIndex)
			    {
				curColorIndex = vertex -> colorIndex;
    
				fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
				fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
				fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
				fprintf(outFile, "    c3s(c);\n");
			    }
    
			    if (vertex -> normal[0] != curNormal[0] ||
				vertex -> normal[1] != curNormal[1] ||
				vertex -> normal[2] != curNormal[2])
			    {
				curNormal[0] = vertex -> normal[0];
				curNormal[1] = vertex -> normal[1];
				curNormal[2] = vertex -> normal[2];
    
				fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
				fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
				fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
				fprintf(outFile, "    n3f(n);\n");
			    }
    
			    fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			    fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			    fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			    fprintf(outFile, "    v3f(v);\n");

			    /*Upper left corner*/
			    offset = RECTMESHVERTEX(rectMesh, i, j + 1);
			    vertex = vertices[offset];
    
			    if (vertex -> colorIndex != curColorIndex)
			    {
				curColorIndex = vertex -> colorIndex;
    
				fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
				fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
				fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
				fprintf(outFile, "    c3s(c);\n");
			    }
    
			    if (vertex -> normal[0] != curNormal[0] ||
				vertex -> normal[1] != curNormal[1] ||
				vertex -> normal[2] != curNormal[2])
			    {
				curNormal[0] = vertex -> normal[0];
				curNormal[1] = vertex -> normal[1];
				curNormal[2] = vertex -> normal[2];
    
				fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
				fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
				fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
				fprintf(outFile, "    n3f(n);\n");
			    }
    
			    fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			    fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			    fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			    fprintf(outFile, "    v3f(v);\n");
		
			    fprintf(outFile, "    endpolygon();\n\n");
			}
		    }
		}
		break;
	    case POLYGONS:
		{
		    PolyPtr poly;
		    VertexPtr vertex;

		    for (poly = ((PolysPtr) curItem) -> polygons;
			 poly;
			 poly = (PolyPtr) (poly -> item . next))
		    {
			fprintf(outFile, "    bgnpolygon();\n");
			for (k = 0; k < poly -> nVertices; ++k)
			{
			    vertex = poly -> vertices[k];
    
			    if (vertex -> colorIndex != curColorIndex)
			    {
				curColorIndex = vertex -> colorIndex;
    
				fprintf(outFile, "    c[0] = %d;\n", curColors[curColorIndex][0]);
				fprintf(outFile, "    c[1] = %d;\n", curColors[curColorIndex][1]);
				fprintf(outFile, "    c[2] = %d;\n", curColors[curColorIndex][2]);
				fprintf(outFile, "    c3s(c);\n");
			    }
    
			    if (vertex -> normal[0] != curNormal[0] ||
				vertex -> normal[1] != curNormal[1] ||
				vertex -> normal[2] != curNormal[2])
			    {
				curNormal[0] = vertex -> normal[0];
				curNormal[1] = vertex -> normal[1];
				curNormal[2] = vertex -> normal[2];
    
				fprintf(outFile, "    n[0] = %#g;\n", curNormal[0]);
				fprintf(outFile, "    n[1] = %#g;\n", curNormal[1]);
				fprintf(outFile, "    n[2] = %#g;\n", curNormal[2]);
				fprintf(outFile, "    n3f(n);\n");
			    }
    
			    fprintf(outFile, "    v[0] = %#g;\n", vertex -> position[0]);
			    fprintf(outFile, "    v[1] = %#g;\n", vertex -> position[1]);
			    fprintf(outFile, "    v[2] = %#g;\n", vertex -> position[2]);
			    fprintf(outFile, "    v3f(v);\n");
			}
			fprintf(outFile, "    endpolygon();\n\n");

			if (poly == ((PolysPtr) curItem) -> lastItem)
			{
			    break;
			}
		    }
		}
	}
	curItem = curItem -> next;
    }
    fprintf(outFile, "    closeobj();\n");
    fclose(outFile);
}
static PolysPtr curPolys;
static ObjPtr curPicture;
static int curColorIndex;

void QuadrantPolys(r1, r2, r3, level)
float r1[3];
float r2[3];
float r3[3];
int level;
/*Emits a quadrant of sphere at center curCenter with radius curRadius given 
  radial normal vectors r1, r2, and r3, using current color table, at recursion 
  level level into curPolys*/
{
    if (level >= curSub)
    {
	/*Maximum number of subdivisions has been reached*/
	Vertex *v[3];
	v[0] = NewVertex(curPicture, VF_SAMEPLACE);
	v[1] = NewVertex(curPicture, VF_SAMEPLACE);
	v[2] = NewVertex(curPicture, VF_SAMEPLACE);
	v[0] -> position[0] = curCenter[0] + r1[0] * curRadius;
	v[0] -> position[1] = curCenter[1] + r1[1] * curRadius;
	v[0] -> position[2] = curCenter[2] + r1[2] * curRadius;
	v[1] -> position[0] = curCenter[0] + r2[0] * curRadius;
	v[1] -> position[1] = curCenter[1] + r2[1] * curRadius;
	v[1] -> position[2] = curCenter[2] + r2[2] * curRadius;
	v[2] -> position[0] = curCenter[0] + r3[0] * curRadius;
	v[2] -> position[1] = curCenter[1] + r3[1] * curRadius;
	v[2] -> position[2] = curCenter[2] + r3[2] * curRadius;

	v[0] -> normal[0] = r1[0];
	v[0] -> normal[1] = r1[1];
	v[0] -> normal[2] = r1[2];
	v[1] -> normal[0] = r2[0];
	v[1] -> normal[1] = r2[1];
	v[1] -> normal[2] = r2[2];
	v[2] -> normal[0] = r3[0];
	v[2] -> normal[1] = r3[1];
	v[2] -> normal[2] = r3[2];

	v[0] -> colorIndex = curColorIndex;
	v[1] -> colorIndex = curColorIndex;
	v[2] -> colorIndex = curColorIndex;

	AppendSPolyToPolys(curPolys, 3, v);
    }
    else
    {
	/*Do another subdivision*/
	float r12[3], r23[3], r31[3];	/*Inner triangle subdivisions*/
	r12[0] = (r1[0] + r2[0]);
	r12[1] = (r1[1] + r2[1]);
	r12[2] = (r1[2] + r2[2]);
	r23[0] = (r2[0] + r3[0]);
	r23[1] = (r2[1] + r3[1]);
	r23[2] = (r2[2] + r3[2]);
	r31[0] = (r3[0] + r1[0]);
	r31[1] = (r3[1] + r1[1]);
	r31[2] = (r3[2] + r1[2]);
	NORMALIZE(r12);
	NORMALIZE(r23);
	NORMALIZE(r31);

	/*Draw the subdivisions*/
	QuadrantPolys(r1, r12, r31, level + 1);
	QuadrantPolys(r12, r2, r23, level + 1);
	QuadrantPolys(r23, r3, r31, level + 1);
	QuadrantPolys(r12, r23, r31, level + 1);
    }
}

#ifdef HAVE_PROTOTYPES
PolysPtr ConvertSphereOntoPicture(ObjPtr picture, VertexPtr center, real radius, int flags, int nSubdivisions)
#else
PolysPtr ConvertSphereOntoPicture(picture, center, radius, flags, nSubdivisions)
ObjPtr picture;
VertexPtr center;
real radius;
int flags;
int nSubdivisions;
#endif
/*Returns a polys object of sphere*/
{
    float vx[3], vy[3], vz[3]; /*Vertices stuck out on those axes*/
    VertexPtr center2;

    curSub = nSubdivisions;	    
    curPolys = AppendPolysToPicture(picture);
    curPicture = picture;
    curPolys -> enclosed = true;

    /*Make the center point*/

    center2 = NewVertex(picture, flags);
    center2 -> colorIndex = center -> colorIndex;
    center2 -> position[0] = center -> position[0];
    center2 -> position[1] = center -> position[1];
    center2 -> position[2] = center -> position[2];


    /*Make the eight quadrants of the sphere*/
    /* +x +y +z */
    vx[0] = 1.0;
    vx[1] = 0.0;
    vx[2] = 0.0;
    vy[0] = 0.0;
    vy[1] = 1.0;
    vy[2] = 0.0;
    vz[0] = 0.0;
    vz[1] = 0.0;
    vz[2] = 1.0;
	
    curRadius = radius;
    curCenter[0] = center -> position[0];
    curCenter[1] = center -> position[1];
    curCenter[2] = center -> position[2];
    curColorIndex = center -> colorIndex;
    QuadrantPolys(vx, vy, vz, 0);
	
    /* +y -x +z*/
    vx[0] = -1.0;
    QuadrantPolys(vy, vx, vz, 0);
	
    /* -x -y -z*/
    vy[1] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);
	
    /* -y +x +z*/
    vx[0] = 1.0;
    QuadrantPolys(vy, vx, vz, 0);
	
    /* +x -y -z*/
    vz[2] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);
	
    /* +y +x -z*/
    vy[1] = 1.0;
    QuadrantPolys(vy, vx, vz, 0);
	
    /* -x +y -z*/
    vx[0] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);
	
    /* -y, -x, -z*/
    vy[1] = -1.0;
    QuadrantPolys(vy, vx, vz, 0);

    return curPolys;
}

#ifdef HAVE_PROTOTYPES
void ConvertFrustumOntoPicture(ObjPtr picture, float end1[3], float rad1,
		float end2[3], float rad2, int nSubdivisions, Bool capEnds)
#else
void ConvertFrustumOntoPicture(picture, end1, rad1, end2, rad2, nSubdivisions, capEnds)
float end1[3], end2[3];
float rad1, rad2;
ObjPtr picture;
int nSubdivisions;
Bool capEnds;
#endif
/*Converts a frustum onto a picture*/
{
    float a[3];		/*Axial unit vector*/
    float r[MAXCYLSIDES][3];
    float s[3], n[3], f[3];
    Vertex v[MAXCYLSIDES + 1];
    double length;	/*Length of cylinder*/
    double nf;		/*Normalization factor*/
			/*Radial units around cylinder*/
    int end;		/*Counter for which end*/
    int b;		/*Counter around radius*/
    int d;		/*Delta for subdivision*/
    int k;		/*Random counter*/
    RectMeshPtr rectMesh;
    VertexPtr vertex;

    curNSides = 4;
    while (nSubdivisions--) curNSides *= 2;

    if (end2[0] == end1[0] &&
	end2[1] == end1[1] &&
	end2[2] == end1[2])
    {
	/*can't do it*/
	return;
    }

    /*First make an axial unit vector*/
    a[0] = end2[0] -
	   end1[0];
    a[1] = end2[1] -
	   end1[1];
    a[2] = end2[2] -
	   end1[2];
    length = sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
    nf = 1.0 / length;
    a[0] *= nf; a[1] *= nf; a[2] *= nf;

    /*See if it's nearly colinear with i*/
    if (ABS(a[0]) > 0.8)
    {
	/*It is, so cross by j to get first r*/
	r[0][0] = -a[2];
	r[0][1] = 0.0;
	r[0][2] = a[0];
    }
    else
    {
	/*It isn't, so cross by i to get first r*/
	r[0][0] = 0.0;
	r[0][1] = a[2];
	r[0][2] = -a[1];
    }
    
    /*Normalize r*/
    length = sqrt(r[0][0]*r[0][0] + r[0][1]*r[0][1] + r[0][2]*r[0][2]);
    nf = 1.0 / length;
    r[0][0] *= nf; r[0][1] *= nf, r[0][2] *= nf;

    /*Cross a with first radial unit to get orthogonal unit*/
    CROSS(a, r[0], r[curNSides / 4]);

    /*Fill out point radii 3 and 4 by reversing 1 and 2*/
    r[curNSides / 2][0] = -r[0][0];
    r[curNSides / 2][1] = -r[0][1];
    r[curNSides / 2][2] = -r[0][2];

    r[curNSides / 4 * 3][0] = -r[curNSides / 4][0];
    r[curNSides / 4 * 3][1] = -r[curNSides / 4][1];
    r[curNSides / 4 * 3][2] = -r[curNSides / 4][2];

    /*Subdivide the sides of the cylinder*/
    k = 0;
    for (d = curNSides / 4; d > 1; d /= 2)
    {
	register double cef_2;	/*Chord extension factor/2*/
	
	cef_2 = cefs_2[k];
	++k;
	
	for (b = 0; b < curNSides; b += d)
	{
	    register int m, e;	/*Midpoint and endpoint*/

	    /*Interpolate center of chord*/
	    e = (b + d) & (curNSides - 1);	/*Clip circular*/
	    m = b + d / 2;
	    r[m][0] = (r[b][0] + r[e][0]) * cef_2;
	    r[m][1] = (r[b][1] + r[e][1]) * cef_2;
	    r[m][2] = (r[b][2] + r[e][2]) * cef_2;
	}
    }

    /*Convert the cylinder*/
    rectMesh = AppendRectMeshToPicture(picture, 2, curNSides + 1, false);

    f[0] = rad2 * r[curNSides - 1][0] + end2[0] - 
	   (rad1 * r[curNSides - 1][0] + end1[0]);
    f[1] = rad2 * r[curNSides - 1][1] + end2[1] - 
	   (rad1 * r[curNSides - 1][1] + end1[1]);
    f[2] = rad2 * r[curNSides - 1][2] + end2[2] - 
	   (rad1 * r[curNSides - 1][2] + end1[2]);
    CROSS(f, r[curNSides - 1], s);
    CROSS(s, f, n);
    NORMALIZE(n);

    vertex = RectMeshVertex(rectMesh, 0, 0);
    vertex -> position[0] = rad2 * r[curNSides - 1][0] + end2[0];
    vertex -> position[1] = rad2 * r[curNSides - 1][1] + end2[1];
    vertex -> position[2] = rad2 * r[curNSides - 1][2] + end2[2];

#if 0
    vertex -> normal[0] = r[curNSides - 1][0];
    vertex -> normal[1] = r[curNSides - 1][1];
    vertex -> normal[2] = r[curNSides - 1][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif

    vertex = RectMeshVertex(rectMesh, 1, 0);
    vertex -> position[0] = rad1 * r[curNSides - 1][0] + end1[0];
    vertex -> position[1] = rad1 * r[curNSides - 1][1] + end1[1];
    vertex -> position[2] = rad1 * r[curNSides - 1][2] + end1[2];

#if 0
    vertex -> normal[0] = r[curNSides - 1][0];
    vertex -> normal[1] = r[curNSides - 1][1];
    vertex -> normal[2] = r[curNSides - 1][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif

    for (k = 0; k < curNSides; ++k)
    {
	f[0] = rad2 * r[k][0] + end2[0] - 
	       (rad1 * r[k][0] + end1[0]);
	f[1] = rad2 * r[k][1] + end2[1] - 
	       (rad1 * r[k][1] + end1[1]);
	f[2] = rad2 * r[k][2] + end2[2] - 
	       (rad1 * r[k][2] + end1[2]);
	CROSS(f, r[k], s);
	CROSS(s, f, n);
	NORMALIZE(n);

	vertex = RectMeshVertex(rectMesh, 0, k + 1);	
	vertex -> position[0] = rad2 * r[k][0] + end2[0];
	vertex -> position[1] = rad2 * r[k][1] + end2[1];
	vertex -> position[2] = rad2 * r[k][2] + end2[2];

#if 0
	vertex -> normal[0] = r[k][0];
	vertex -> normal[1] = r[k][1];
	vertex -> normal[2] = r[k][2];
#else
	vertex -> normal[0] = n[0];
	vertex -> normal[1] = n[1];
	vertex -> normal[2] = n[2];
#endif

	vertex = RectMeshVertex(rectMesh, 1, k + 1);	
	vertex -> position[0] = rad1 * r[k][0] + end1[0];
	vertex -> position[1] = rad1 * r[k][1] + end1[1];
	vertex -> position[2] = rad1 * r[k][2] + end1[2];

#if 0
	vertex -> normal[0] = r[k][0];
	vertex -> normal[1] = r[k][1];
	vertex -> normal[2] = r[k][2];
#else
	vertex -> normal[0] = n[0];
	vertex -> normal[1] = n[1];
	vertex -> normal[2] = n[2];
#endif
    }

    InterpRectCenters(rectMesh);

    if (capEnds)
    {
	/*Do top cap*/
	for (k = 0; k < curNSides; ++k)
	{
	    v[k] . position[0] = rad2 * r[k][0] + end2[0];
	    v[k] . position[1] = rad2 * r[k][1] + end2[1];
	    v[k] . position[2] = rad2 * r[k][2] + end2[2];
	    v[k] . normal[0] = a[0];
	    v[k] . normal[1] = a[1];
	    v[k] . normal[2] = a[2];
	}
	AppendPolyToPicture(picture, curNSides, v);
    
	/*Do bottom cap*/
	for (k = 0; k < curNSides; ++k)
	{
	    v[k] . position[0] = rad1 * r[k][0] + end1[0];
	    v[k] . position[1] = rad1 * r[k][1] + end1[1];
	    v[k] . position[2] = rad1 * r[k][2] + end1[2];
	    v[k] . normal[0] = -a[0];
	    v[k] . normal[1] = -a[1];
	    v[k] . normal[2] = -a[2];
	}
	AppendPolyToPicture(picture, curNSides, v);
    }
}

#ifdef HAVE_PROTOTYPES
PicItemPtr ColorItemsByItems(PicItemPtr destItem, ObjPtr owner, PicItemPtr curItem)
#else
PicItemPtr ColorItemsByItems(dest, owner, curItem)
ObjPtr dest, pic;
ObjPtr owner;
PicItemPtr curItem;
#endif
/*Colors a picture with another picture, which must be the source picture*/
{
    ObjPtr var;
    Bool capEnds;

    capEnds = GetPredicate(owner, CAPENDSP);

    while (curItem)
    {
	if (!destItem)
	{
	    ReportError("ColorItemsByItems", "Pictures don't match");
	    break;
	}
	switch (curItem -> type)
	{
	    case POLYGON:
	    case POLYLINE:
		ColorItemWithIndex(destItem, 
			((PolyPtr) curItem) -> vertices[0] -> colorIndex);
		break;
	    case FRUSTUM:

		ColorItemWithIndex(destItem,
			((FrustumPtr) curItem) -> colorIndex);
		if (capEnds)
		{
		    /*Cap the ends*/
		    destItem = destItem -> next;
		    if (!destItem) break;
		    ColorItemWithIndex(destItem,
			((FrustumPtr) curItem) -> colorIndex);
		    destItem = destItem -> next;
		    if (!destItem) break;
		    ColorItemWithIndex(destItem,
			((FrustumPtr) curItem) -> colorIndex);
		}
		break;
	    case SPHERE:

		ColorItemWithIndex(destItem,
			((SpherePtr) curItem) -> centerVertex -> colorIndex);
		break;
	    case POLYGONS:

		{
		    PolyPtr polygons;
		    polygons = ((PolysPtr) curItem) -> polygons;
		    if (polygons)
		    {
			ColorItemWithIndex((PicItemPtr) polygons,
				polygons -> vertices[0] -> colorIndex);
		    }
		}
		break;
	    case RECTMESH:

		ColorItemWithIndex(destItem,
			((RectMeshPtr) curItem) -> vertices[0] -> colorIndex);
		break;

	}
	curItem = curItem -> next;
	destItem = destItem -> next;
    }
    return destItem;
}

#ifdef HAVE_PROTOTYPES
void ColorPictureByPicture(ObjPtr dest, ObjPtr owner, ObjPtr pic)
#else
void ColorPictureByPicture(dest, owner, pic)
ObjPtr dest, pic;
ObjPtr owner;
#endif
/*Colors a picture with another picture, which must be the source picture*/
{
    PicItemPtr curItem, destItem;

    curItem = ((PicPtr) pic) -> items;
    destItem = ((PicPtr) dest) -> items;
    ColorItemsByItems(destItem, owner, curItem);
}

#ifdef HAVE_PROTOTYPES
ObjPtr ConvertOntoPicture(ObjPtr retVal, ObjPtr pic, ObjPtr owner)
#else
ObjPtr ConvertOntoPicture(retVal, pic, owner)
ObjPtr retVal;
ObjPtr pic;
ObjPtr owner;
#endif
/*Returns a copy of pic converted into internal form on the end of retVal.  
  Owner can give info on conversion*/
{
    PicItemPtr curItem;
    ObjPtr var;
    int sphereSub;
    int frustumSub;
    Bool capEnds;

    var = GetIntVar("ConvertOntoPicture", owner, SPHERESUBDIV);
    if (var)
    {
	sphereSub = GetInt(var);
    }
    else
    {
	sphereSub = 2;
    }

    var = GetIntVar("ConvertOntoPicture", owner, FRUSTUMSUBDIV);
    if (var)
    {
	frustumSub = GetInt(var);
    }
    else
    {
	frustumSub = 2;
    }

    capEnds = GetPredicate(owner, CAPENDSP);

    curItem = ((PicPtr) pic) -> items;
    while(curItem)
    {
	switch (curItem -> type)
	{
	    case POLYGON:
		ConvertPolyOntoPicture(retVal, (PolyPtr) curItem);
		break;
	    case POLYLINE:
		ConvertPolyOntoPicture(retVal, (PolyPtr) curItem);
		break;
	    case FRUSTUM:
		ConvertFrustumOntoPicture(retVal, 
			((FrustumPtr) curItem) -> end1,
			((FrustumPtr) curItem) -> rad1,
			((FrustumPtr) curItem) -> end2,
			((FrustumPtr) curItem) -> rad2,
			frustumSub, capEnds);
		break;
	    case SPHERE:
		ConvertSphereOntoPicture(retVal,
			((SpherePtr) curItem) -> centerVertex,
			((SpherePtr) curItem) -> radius,
			0,
			sphereSub);
		break;
	    case POLYGONS:
		{
		    PolyPtr polygons;
		    PolysPtr polySet;
		    polySet = AppendPolysToPicture(retVal);
		    polygons = ((PolysPtr) curItem) -> polygons;
		    while (polygons)
		    {
			ConvertPolyOntoPolys(polySet, polygons);
			polygons = (PolyPtr) (polygons -> item . next);
		    }
		}
		break;
	    case RECTMESH:
		{
		    RectMeshPtr rectMesh;
		    long k;
		    rectMesh = AppendRectMeshToPicture(retVal,
			((RectMeshPtr) curItem) -> xDim, 
			((RectMeshPtr) curItem) -> yDim,
			((RectMeshPtr) curItem) -> inCenter);
		    for (k = 0; k < (((RectMeshPtr) curItem) -> xDim * (2 * ((RectMeshPtr) curItem) -> yDim - 1)); ++k)
		    {
			((RectMeshPtr) rectMesh) -> vertices[k] -> normal[0] =
			((RectMeshPtr) curItem) -> vertices[k] -> normal[0];
			((RectMeshPtr) rectMesh) -> vertices[k] -> normal[1] =
			((RectMeshPtr) curItem) -> vertices[k] -> normal[1];
			((RectMeshPtr) rectMesh) -> vertices[k] -> normal[2] =
			((RectMeshPtr) curItem) -> vertices[k] -> normal[2];
			((RectMeshPtr) rectMesh) -> vertices[k] -> position[0] =
			((RectMeshPtr) curItem) -> vertices[k] -> position[0];
			((RectMeshPtr) rectMesh) -> vertices[k] -> position[1] =
			((RectMeshPtr) curItem) -> vertices[k] -> position[1];
			((RectMeshPtr) rectMesh) -> vertices[k] -> position[2] =
			((RectMeshPtr) curItem) -> vertices[k] -> position[2];
			((RectMeshPtr) rectMesh) -> vertices[k] -> colorIndex =
			((RectMeshPtr) curItem) -> vertices[k] -> colorIndex;
		    }
		}
		break;
	}
	curItem = curItem -> next;
    }
    return retVal;
}

#ifdef HAVE_PROTOTYPES
ObjPtr ConvertPicture(ObjPtr pic, ObjPtr owner)
#else
ObjPtr ConvertPicture(pic, owner)
ObjPtr pic;
ObjPtr owner;
#endif
/*Returns a copy of pic converted into internal form.  Owner can give info on conversion*/
{
    ObjPtr retVal;
    retVal = NewPicture();
    ConvertOntoPicture(retVal, pic, owner);
    return retVal;
}

void GetPictureBounds(pic, bounds)
ObjPtr pic;
real bounds[6];
/*Returns the bounds of pic in bounds*/
{
    PicItemPtr curItem;
    bounds[0] = bounds[2] = bounds[4] = PLUSINF;
    bounds[1] = bounds[3] = bounds[5] = MINUSINF;

    curItem = ((PicPtr) pic) -> items;

    while (curItem)
    {
	switch(curItem -> type)
	{
	    case SPHERE:
		if (((SpherePtr) curItem) -> centerVertex -> position[0] -
		    ((SpherePtr) curItem) -> radius < bounds[0])
		{
		    bounds[0] = ((SpherePtr) curItem) -> centerVertex -> position[0] -
				((SpherePtr) curItem) -> radius;
		}
		if (((SpherePtr) curItem) -> centerVertex -> position[0] +
		    ((SpherePtr) curItem) -> radius > bounds[1])
		{
		    bounds[1] = ((SpherePtr) curItem) -> centerVertex -> position[0] +
				((SpherePtr) curItem) -> radius;
		}
		if (((SpherePtr) curItem) -> centerVertex -> position[1] -
		    ((SpherePtr) curItem) -> radius < bounds[2])
		{
		    bounds[2] = ((SpherePtr) curItem) -> centerVertex -> position[1] -
				((SpherePtr) curItem) -> radius;
		}
		if (((SpherePtr) curItem) -> centerVertex -> position[1] +
		    ((SpherePtr) curItem) -> radius > bounds[3])
		{
		    bounds[3] = ((SpherePtr) curItem) -> centerVertex -> position[1] +
				((SpherePtr) curItem) -> radius;
		}
		if (((SpherePtr) curItem) -> centerVertex -> position[2] -
		    ((SpherePtr) curItem) -> radius < bounds[4])
		{
		    bounds[4] = ((SpherePtr) curItem) -> centerVertex -> position[2] -
				((SpherePtr) curItem) -> radius;
		}
		if (((SpherePtr) curItem) -> centerVertex -> position[2] +
		    ((SpherePtr) curItem) -> radius > bounds[5])
		{
		    bounds[5] = ((SpherePtr) curItem) -> centerVertex -> position[2] +
				((SpherePtr) curItem) -> radius;
		}
		break;
	    case FRUSTUM:
		if (((FrustumPtr) curItem) -> end1[0] - 
		    ((FrustumPtr) curItem) -> rad1 < bounds[0])
		{
		    bounds[0] = ((FrustumPtr) curItem) -> end1[0] - 
		    		((FrustumPtr) curItem) -> rad1;
		}
		if (((FrustumPtr) curItem) -> end1[0] + 
		    ((FrustumPtr) curItem) -> rad1 > bounds[1])
		{
		    bounds[1] = ((FrustumPtr) curItem) -> end1[0] + 
		    		((FrustumPtr) curItem) -> rad1;
		}
		if (((FrustumPtr) curItem) -> end1[1] - 
		    ((FrustumPtr) curItem) -> rad1 < bounds[2])
		{
		    bounds[2] = ((FrustumPtr) curItem) -> end1[1] - 
		    		((FrustumPtr) curItem) -> rad1;
		}
		if (((FrustumPtr) curItem) -> end1[1] + 
		    ((FrustumPtr) curItem) -> rad1 > bounds[3])
		{
		    bounds[3] = ((FrustumPtr) curItem) -> end1[1] + 
		    		((FrustumPtr) curItem) -> rad1;
		}
		if (((FrustumPtr) curItem) -> end1[2] - 
		    ((FrustumPtr) curItem) -> rad1 < bounds[4])
		{
		    bounds[4] = ((FrustumPtr) curItem) -> end1[2] - 
		    		((FrustumPtr) curItem) -> rad1;
		}
		if (((FrustumPtr) curItem) -> end1[2] + 
		    ((FrustumPtr) curItem) -> rad1 > bounds[5])
		{
		    bounds[5] = ((FrustumPtr) curItem) -> end1[2] + 
		    		((FrustumPtr) curItem) -> rad1;
		}

		if (((FrustumPtr) curItem) -> end2[0] - 
		    ((FrustumPtr) curItem) -> rad2 < bounds[0])
		{
		    bounds[0] = ((FrustumPtr) curItem) -> end2[0] - 
		    		((FrustumPtr) curItem) -> rad2;
		}
		if (((FrustumPtr) curItem) -> end2[0] + 
		    ((FrustumPtr) curItem) -> rad2 > bounds[1])
		{
		    bounds[1] = ((FrustumPtr) curItem) -> end2[0] + 
		    		((FrustumPtr) curItem) -> rad2;
		}
		if (((FrustumPtr) curItem) -> end2[1] - 
		    ((FrustumPtr) curItem) -> rad2 < bounds[2])
		{
		    bounds[2] = ((FrustumPtr) curItem) -> end2[1] - 
		    		((FrustumPtr) curItem) -> rad2;
		}
		if (((FrustumPtr) curItem) -> end2[1] + 
		    ((FrustumPtr) curItem) -> rad2 > bounds[3])
		{
		    bounds[3] = ((FrustumPtr) curItem) -> end2[1] + 
		    		((FrustumPtr) curItem) -> rad2;
		}
		if (((FrustumPtr) curItem) -> end2[2] - 
		    ((FrustumPtr) curItem) -> rad2 < bounds[4])
		{
		    bounds[4] = ((FrustumPtr) curItem) -> end2[2] - 
		    		((FrustumPtr) curItem) -> rad2;
		}
		if (((FrustumPtr) curItem) -> end2[2] + 
		    ((FrustumPtr) curItem) -> rad2 > bounds[5])
		{
		    bounds[5] = ((FrustumPtr) curItem) -> end2[2] + 
		    		((FrustumPtr) curItem) -> rad2;
		}
		break;
	    case POLYGON:
	    case POLYLINE:
		{
		    int k;
		    VertexPtr *curVertex;

		    curVertex = ((PolyPtr) curItem) -> vertices;
		    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)
		    {
			if ((*curVertex) -> position[0] < bounds[0])
			{
			    bounds[0] = (*curVertex) -> position[0];
			}
			if ((*curVertex) -> position[0] > bounds[1])
			{
			    bounds[1] = (*curVertex) -> position[0];
			}
			if ((*curVertex) -> position[1] < bounds[2])
			{
			    bounds[2] = (*curVertex) -> position[1];
			}
			if ((*curVertex) -> position[1] > bounds[3])
			{
			    bounds[3] = (*curVertex) -> position[1];
			}
			if ((*curVertex) -> position[2] < bounds[4])
			{
			    bounds[4] = (*curVertex) -> position[2];
			}
			if ((*curVertex) -> position[2] > bounds[5])
			{
			    bounds[5] = (*curVertex) -> position[2];
			}
			++curVertex;
		    }
		}
		break;
	    case POLYGONS:
		{
		    PolyPtr polygons;
		    polygons = ((PolysPtr) curItem) -> polygons;
		    while (polygons)
		    {
			int k;
			VertexPtr *curVertex;

			curVertex = polygons -> vertices;
			for (k = 0; k < polygons -> nVertices; ++k)
			{
			    if ((*curVertex) -> position[0] < bounds[0])
			    {
				bounds[0] = (*curVertex) -> position[0];
			    }
			    if ((*curVertex) -> position[0] > bounds[1])
			    {
				bounds[1] = (*curVertex) -> position[0];
			    }
			    if ((*curVertex) -> position[1] < bounds[2])
			    {
				bounds[2] = (*curVertex) -> position[1];
			    }
			    if ((*curVertex) -> position[1] > bounds[3])
			    {
				bounds[3] = (*curVertex) -> position[1];
			    }
			    if ((*curVertex) -> position[2] < bounds[4])
			    {
				bounds[4] = (*curVertex) -> position[2];
			    }
			    if ((*curVertex) -> position[2] > bounds[5])
			    {
				bounds[5] = (*curVertex) -> position[2];
			    }
			    ++curVertex;
			}
			polygons = (PolyPtr) (polygons -> item . next);
		    }
		}
		break;
	    case RECTMESH:
		{
		    register long i, j, offset;
		    VertexPtr *vertices;

		    vertices = ((RectMeshPtr) curItem) -> vertices;

		    for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
		    {
			for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
			{
			    offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);
			    if (vertices[offset] -> position[0] < bounds[0])
			    {
				bounds[0] = vertices[offset] -> position[0];
			    }
			    if (vertices[offset] -> position[0] > bounds[1])
			    {
				bounds[1] = vertices[offset] -> position[0];
			    }
			    if (vertices[offset] -> position[1] < bounds[2])
			    {
				bounds[2] = vertices[offset] -> position[1];
			    }
			    if (vertices[offset] -> position[1] > bounds[3])
			    {
				bounds[3] = vertices[offset] -> position[1];
			    }
			    if (vertices[offset] -> position[2] < bounds[4])
			    {
				bounds[4] = vertices[offset] -> position[2];
			    }
			    if (vertices[offset] -> position[2] > bounds[5])
			    {
				bounds[5] = vertices[offset] -> position[2];
			    }
			}
		    }
		}
		break;
	}
	curItem = curItem -> next;
    }
}

#ifdef HAVE_PROTOTYPES
void OpenGraphObj(ObjPtr object)
#else
ObjPtr OpenGraphObj(object)
ObjPtr object;
#endif
/*Opens an existing graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    makeobj(GetInt(object));
#endif
#endif
}

#ifdef HAVE_PROTOTYPES
void CloseGraphObj(ObjPtr object)
#else
void CloseGraphObj(object)
ObjPtr object;
#endif
/*Closes a graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    closeobj();
#endif
#endif
}

#ifdef HAVE_PROTOTYPES
void DrawGraphObj(ObjPtr object)
#else
void DrawGraphObj(object)
ObjPtr object;
#endif
/*Draw a graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    callobj(GetInt(object));
#endif
#endif
}

#ifdef GRAPHICS
#ifdef GRAPHOBJ
static ObjPtr CleanupGraphObj(object)
ObjPtr object;
/*Cleans up a graphical object*/
{
    delobj(GetInt(object));
    return ObjTrue;
}
#endif
#endif

#ifdef HAVE_PROTOTYPES
ObjPtr NewOpenGraphObj(void)
#else
ObjPtr NewOpenGraphObj()
#endif
/*Returns a pointer to a new, opened graphical object.
*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    Object obj;
    ObjPtr retVal;

    obj = genobj();

    retVal = NewInt(obj);
    SetMethod(retVal, CLEANUP, CleanupGraphObj);
    OpenGraphObj(retVal);
    return retVal;
#endif
#endif
}