/*ScianVisContours.c
  June 14, 1991
  Eric Pepke
  Routines for contour visualization object

  The contour object takes a scalar field defined over a data form with
  topological dimension 2.  The spatial dimension of the dataset may be 
  either 2 or 3.  If the spatial dimension is 2, the Z is assumed to be
  zero.  The contour visualization object works for regular data forms,
  and curvilinear data forms.

*/

#include "Scian.h"
#include "ScianTypes.h"
#include "ScianArrays.h"
#include "ScianWindows.h"
#include "ScianTextBoxes.h"
#include "ScianObjWindows.h"
#include "ScianEvents.h"
#include "ScianIcons.h"
#include "ScianColors.h"
#include "ScianControls.h"
#include "ScianLists.h"
#include "ScianSpaces.h"
#include "ScianSliders.h"
#include "ScianIDs.h"
#include "ScianDatasets.h"
#include "ScianErrors.h"
#include "ScianVisObjects.h"
#include "ScianVisContours.h"
#include "ScianStyle.h"
#include "ScianPictures.h"
#include "ScianTitleBoxes.h"
#include "ScianButtons.h"
#include "ScianMethods.h"
#include "ScianDraw.h"
#include "ScianTemplates.h"
#include "ScianTemplateHelper.h"

ObjPtr contourLineClass;		/*Contour line object*/
ObjPtr contoursClass;

typedef struct xyzStruct
    {
	struct xyzStruct *next;		/*Pointer to next in list*/
	real x, y, z;			/*Coordinates*/
	real nx, ny, nz;		/*Normals*/
    } XYZ;

#ifdef PROTO
static Bool NextContour(ObjPtr, real *);
unsigned char FindQuadJoinings(int initState[2][2]);
void FollowContour(ObjPtr surface, unsigned char *joinings, 
			long i, long j, long iSize, long jSize,
			int startJoin, real contourValue, Bool reverse,
			int lineWidth);
XYZ *NewXYZNode(long i, long j, int dir, real contourValue);
#else
static Bool NextContour();
unsigned char FindQuadJoinings();
void FollowContour();
XYZ *NewXYZNode();
#endif



static ObjPtr ContoursInit(object)
ObjPtr object;
/*Initializes a contours object*/
{
    ObjPtr minMax;
    real bounds[6];
    real xySize, zSize;
    ObjPtr contourField, mainDataset, colorObj, deformObj;
    real min, max, *elements;
    double ddiff, majorWidth;
    int nTics;

    MakeVar(object, MAINDATASET);
    SetVar(object, COLOROBJ, colorObj = GetVar(object, MAINDATASET));
    if (colorObj)
    {
	SetVar(object, COLORS, ObjTrue);
    }
    SetVar(object, DEFORMOBJ, deformObj = GetVar(object, MAINDATASET));

    /*Heuristically determine whether to displace Z*/
    GetBounds(object, bounds);
    xySize = bounds[1] - bounds[0];
    xySize = MAX(xySize, bounds[3] - bounds[2]);
    xySize = MAX(xySize, bounds[5] - bounds[4]);
    if (GetVar(colorObj, DATA))
    {
	MakeVar(deformObj, MINMAX);
	minMax = GetVar(deformObj, MINMAX);
	zSize = ((real *) ELEMENTS(minMax))[1] - ((real *) ELEMENTS(minMax))[0];
	if (zSize < xySize * 0.5)
	{
	    SetVar(object, DEFORMSWITCH, NewInt(1));
	}
	else
	{
	    SetVar(object, DEFORMSWITCH, NewInt(0));
	}
    }

    /*Get the contour field*/
    contourField = GetObjectVar("ContoursInit", object, MAINDATASET);
    if (!contourField) return ObjFalse;
    while (mainDataset = GetVar(contourField, MAINDATASET))
    {
	contourField = mainDataset;
    }

    /*Find its min and max*/
    MakeVar(contourField, MINMAX);
    minMax = GetVar(contourField, MINMAX);
    if (minMax)
    {
	elements = ELEMENTS(minMax);
	min = elements[0];
	max = elements[1];
    }
    else
    {
	min = 0.0;
	max = 1.0;
    }

    ddiff = max - min;
    CalcGoodSteps(ddiff, 50, 10, &majorWidth, &nTics);

    SetVar(object, CONTOURCENTER, NewReal(0.0));
    SetVar(object, CONTOURSTEP, NewReal(majorWidth));
    SetVar(object, CONTOURMAX, NewReal(PLUSINF));
    SetVar(object, CONTOURDIV, NewInt(MAX(1, nTics)));
    SetVar(object, CONTOURMIN, NewReal(MINUSINF));

    return ObjTrue;
}

static ObjPtr SetContoursMainDataset(visObj, dataSet)
ObjPtr visObj, dataSet;
/*Sets the main data set of visObj to dataSet*/
{
    SetVar(visObj, MAINDATASET, dataSet);
    return ObjTrue;
}

static ObjPtr MakeContoursColored(visObject)
ObjPtr visObject;
/*Makes the contours colored*/
{
    SetVar(visObject, PICCOLORED, ObjTrue);
    if (GetPredicate(visObject, COLORS))
    {
	ObjPtr contourField, contourForm;
	ObjPtr colorField, colorForm;
	ObjPtr var;
	ObjPtr surface;
	ObjPtr palette;

	/*Get the contour field and its form*/
	contourField = GetObjectVar("MakeContoursColored", visObject, MAINDATASET);
	if (!contourField) return ObjFalse;
	contourForm = GetObjectVar("MakeContoursColored", contourField, DATAFORM);
	if (!contourForm) return ObjFalse;

	/*Get the color field and its form*/
	colorField = GetObjectVar("MakeContoursColored", visObject, COLOROBJ);
	if (!colorField) return ObjFalse;
	colorForm = GetObjectVar("MakeContoursColored", colorField, DATAFORM);
	if (!colorForm) return ObjFalse;

	/*Get the color palette*/
	MakeVar(colorField, CPALETTE);
	palette = GetPaletteVar("MakeContoursColored", colorField, CPALETTE);
	if (!palette)
	{
	    return ObjFalse;
	}
	SetPalette(palette);

	/*Get the surface to color*/
	surface = GetPictureVar("MakeContoursColored", visObject, SURFACE);
	if (!surface) return ObjFalse;

	/*Have to make it colored by the object.*/
	ColorPictureByObject(surface, colorField, GetPredicate(visObject, INTERPCOLORS));
    }
    return ObjTrue;
}

#ifdef PROTO
static Bool NextContour(ObjPtr contourList, real *contourValue)
#else
static Bool NextContour(contourList, contourValue)
ObjPtr contourList;
real *contourValue;
#endif
/*Stuffs the next contour ABOVE contourValue in contourList into contourValue
*/
{
    ObjPtr var;
    real test, closest;
    Bool closestSet = false;
    ThingListPtr runner;

    runner = LISTOF(contourList);
    while (runner)
    {
	var = GetRealVar("NextContour", runner -> thing, STARTVALUE);
	if (var)
	{
	    test = GetReal(var);
	    if (test > *contourValue)
	    {
		if (!closestSet || test < closest)
		{
		    /*Found the next one*/
		    closestSet = true;
		    closest = test;
		}
	    }
	}
	runner = runner -> next;
    }

    if (closestSet)
    {
	*contourValue = closest;
    }
    return closestSet;
}

/*Quadrilateral for contouring is numbered like this

  2       3
   +-----+
   |  2  |
j  |3   1|
   |  0  |
   +-----+
  0       1
      i
*/

unsigned char fixedJoinings[] =
	{
	    0x00,	/* 0000 */
	    0x30,	/* 0001 */
	    0x10,	/* 0010 */
	    0x70,	/* 0011 */
	    0x0B,	/* 0100 */
	    0x20,	/* 0101 */
	    0x1B,	/* 0110 */
	    0x60,	/* 0111 */
	    0x60,	/* 1000 */
	    0x36,	/* 1001 */
	    0x20,	/* 1010 */
	    0x0B,	/* 1011 */
	    0x70,	/* 1100 */
	    0x10,	/* 1101 */
	    0x30,	/* 1110 */
	    0x00	/* 1111 */
	};
unsigned char joinChars[] =
	{
	    ' ',	/* 0000 */
	    '\\',	/* 0001 */
	    '/',	/* 0010 */
	    '-',	/* 0011 */
	    '/',	/* 0100 */
	    '|',	/* 0101 */
	    '/',	/* 0110 */
	    '\\',	/* 0111 */
	    '\\',	/* 1000 */
	    '\\',	/* 1001 */
	    '|',	/* 1010 */
	    '/',	/* 1011 */
	    '-',	/* 1100 */
	    '/',	/* 1101 */
	    '\\',	/* 1110 */
	    ' '		/* 1111 */
	};

#ifdef PROTO
unsigned char FindQuadJoinings(int initState[2][2])
#else
unsigned char FindJoinings(initState)
int initState[2][2];
#endif
/*Finds the joinings for initial state initState*/
{
    int index;

    index = initState[0][0] +
	    (initState[1][0] << 1) +
	    (initState[0][1] << 2) +
	    (initState[1][1] << 3);

    return fixedJoinings[index];
}

#define SELJOIN(j, n) (((n) ? ((j) >> 4) : (j)) & 0xF)
#define WIPEJOIN(j, n) ((n) ? ((j) &= 0xF) : ((j) &= 0xF0))
#define SELEND(j, n) (((n) ? ((j) >> 2) : (j)) & 0x3)

#ifdef PROTO
XYZ *NewXYZNode(long i, long j, int dir, real contourValue)
#else
XYZ *NewXYZNode(i, j, dir, contourValue)
long i, y;
int dir;
real contourValue;
#endif
/*Returns a new XYZ node for dir from i, j assuming contourValue.
  Expects dataset in FIELD1 and data form in FIELD2.*/
{
    XYZ *node;
    long iPlus, jPlus;
    real x1, y1, z1, f1, x2, y2, z2, f2, w1, w2;
    long curDim[2];
    int nComponents;
    real c[2][2][3];		/*Corner position*/
    real cf[2][2];		/*Corner function*/
    real v1[3], v2[3];		/*Vectors*/
    real x[3];			/*Cross product*/

    /*Make the node*/
    node = (XYZ *) Alloc(sizeof(XYZ));
    node -> next = NULL;

    /*Get the components in the data field*/
    nComponents = GetNComponents(FIELD2);

    /*Fill the corners*/
    for (iPlus = 0; iPlus < 2; ++iPlus)
    {
	for (jPlus = 0; jPlus < 2; ++jPlus)
	{
	    curDim[0] = i + iPlus;
	    curDim[1] = j + jPlus;
	    c[iPlus][jPlus][0] = SelectFieldComponent(FIELD2, 0, curDim);
	    c[iPlus][jPlus][1] = SelectFieldComponent(FIELD2, 1, curDim);
	    c[iPlus][jPlus][2] = nComponents >= 3 ?						
			SelectFieldComponent(FIELD2, 2, curDim) : 0.0;
	    cf[iPlus][jPlus] = SelectFieldScalar(FIELD1, curDim);	
	}
    }

    /*First endpoint*/							
    switch(dir)								
    {									
	case 0:								
	    /*Do 0, 0*/							
	    iPlus = 0;						
	    jPlus = 0;						
	    break;							
	case 1:								
	    /*Do 1, 0*/							
	    iPlus = 1;					
	    jPlus = 0;						
	    break;							
	case 2:								
	    /*Do 1, 1*/							
	    iPlus = 1;					
	    jPlus = 1;					
	    break;							
	case 3:								
	    /*Do 0, 1*/							
	    iPlus = 0;						
	    jPlus = 1;					
	    break;							
    }									
    x1 = c[iPlus][jPlus][0];			
    y1 = c[iPlus][jPlus][1];			
    z1 = c[iPlus][jPlus][2];			
    f1 = cf[iPlus][jPlus];			
									
    /*Second endpoint*/							
    switch(dir)								
    {									
	case 0:								
	    /*Do 1, 0*/							
	    iPlus = 1;					
	    jPlus = 0;						
	    break;							
	case 1:								
	    /*Do 1, 1*/							
	    iPlus = 1;					
	    jPlus = 1;					
	    break;							
	case 2:								
	    /*Do 0, 1*/							
	    iPlus = 0;						
	    jPlus = 1;					
	    break;							
	case 3:								
	    /*Do 0, 0*/							
	    iPlus = 0;						
	    jPlus = 0;						
	    break;							
    }									
    x2 = c[iPlus][jPlus][0];			
    y2 = c[iPlus][jPlus][1];			
    z2 = c[iPlus][jPlus][2];			
    f2 = cf[iPlus][jPlus];			
																		
    /*Now interpolate*/
    if ((f1 <= contourValue) && (contourValue <= f2))
    {
	real rDiff = 1.0 / (f2 - f1);
	w1 = (f2 - contourValue) * rDiff;
	w2 = (contourValue - f1) * rDiff;
    }
    else if ((f2 <= contourValue) && (contourValue <= f1))
    {
	real rDiff = 1.0 / (f1 - f2);
	w2 = (f1 - contourValue) * rDiff;
	w1 = (contourValue - f2) * rDiff;
    }
    else
    {
	w1 = w2 = 0.0;
	fprintf(stderr, "Contour error at (%d %d)\n", i, j);
    }

    node -> x = x1 * w1 + x2 * w2;
    node -> y = y1 * w1 + y2 * w2;
    node -> z = z1 * w1 + z2 * w2;

    /*Now find the normals*/
    v1[0] = c[0][1][0] - c[0][0][0];
    v1[1] = c[0][1][1] - c[0][0][1];
    v1[2] = c[0][1][2] - c[0][0][2];

    v2[0] = c[1][0][0] - c[0][0][0];
    v2[1] = c[1][0][1] - c[0][0][1];
    v2[2] = c[1][0][2] - c[0][0][2];

    CROSS(v1, v2, x);
    NORMALIZE(x);
    node -> nx = x[0];
    node -> ny = x[1];
    node -> nz = x[2];

    return node;
}


#ifdef PROTO
void FollowContour(ObjPtr surface, unsigned char *joinings, 
			long i, long j, long iSize, long jSize,
			int startJoin,
			real contourValue, Bool reverse, int lineWidth)
#else
void FollowContour(surface, joinings, i, j, iSize, jSize, startJoin, 
	contourValue, reverse, lineWidth)
ObjPtr surface;
unsigned char *joinings; 
long i; 
long j; 
long iSize; 
long jSize;
int startJoin;
real contourValue;
Bool reverse;
int lineWidth;
#endif
/*Follows a contour within surface using joinings starting at startJoin at i, j.
startJoin is 0 for LSN and 1 for MSN;
iSize and jSize gives the size of the array.  Dataset is assumed to be in
FIELD1 and dataform in FIELD2.  This routine runs contours both ways.
Side effects: removes runs from joinings.  If reverse is true, reverses the
sense of the normal*/
{
    unsigned int joining, a, b, iCur, jCur, aCur, bCur, match, curJoin;
    long count, k;
    XYZ *first, *last, *runner;
    long curDim[2];
    VertexPtr vertices;
    Bool firstTime;

    /*Set joining to the joining to start*/
    joining = SELJOIN(joinings[i + j * iSize], startJoin);

    /*Set a and b to the two ends of the joining*/
    a = SELEND(joining, 1);
    b = SELEND(joining, 0);

    /*Return if null*/
    if (a == b) return;

    /*Want to go from a to b.  Reverse if neccessary.*/
    if ((i == 0 && b == 3) ||		/*Off left side*/
	(i == iSize - 2 && b == 1) ||	/*Off right side*/
	(j == 0 && b == 0) ||		/*Off bottom side*/
	(j == jSize - 2 && b == 2))	/*Off top side*/ 
    {
	int t;
	t = a;
	a = b;
	b = t;
    }

    /*Make the first XYZ in the contour line*/
    first = last = NewXYZNode(i, j, a, contourValue);
    last -> next = NULL;

    iCur = i;
    jCur = j;
    aCur = a;
    bCur = b;
    curJoin = startJoin;
    /*Now run the contour, from aCur to bCur within iCur, jCur*/
    for(;;)
    {
	/*Emit this joining onto the list*/
	last -> next = NewXYZNode(iCur, jCur, bCur, contourValue);
	last = last -> next;

	/*Delete this joining*/
	WIPEJOIN(joinings[iCur + jCur * iSize], curJoin);

	/*Find a next place to go*/
	switch(bCur)
	{
	    case 0:
		--jCur;
		match = 2;
		break;
	    case 1:
		++iCur;
		match = 3;
		break;
	    case 2:
		++jCur;
		match = 0;
		break;
	    case 3:
		--iCur;
		match = 1;
		break;
	}

	if (iCur < 0 || jCur < 0 || iCur >= iSize - 1 || jCur >= jSize - 1)
	{
	    /*Fallen off the edge*/
	    break;
	}

	/*Get the new joining and see if it matches up.  Try 0 first.*/
	joining = SELJOIN(joinings[iCur + jCur * iSize], 0);
	if (joining)
	{
	    if (SELEND(joining, 1) == match)
	    {
		aCur = SELEND(joining, 1);
		bCur = SELEND(joining, 0);
		curJoin = 0;
		continue;
	    }
	    if (SELEND(joining, 0) == match)
	    {
		aCur = SELEND(joining, 0);
		bCur = SELEND(joining, 1);
		curJoin = 0;
		continue;
	    }
	}

	/*OK, so now try one*/
	joining = SELJOIN(joinings[iCur + jCur * iSize], 1);
	if (joining)
	{
	    if (SELEND(joining, 1) == match)
	    {
		aCur = SELEND(joining, 1);
		bCur = SELEND(joining, 0);
		curJoin = 1;
		continue;
	    }
	    if (SELEND(joining, 0) == match)
	    {
		aCur = SELEND(joining, 0);
		bCur = SELEND(joining, 1);
		curJoin = 1;
		continue;
	    }
	}

	break;
    }

    /*Follow the contour backward from the beginning, just in case*/
    iCur = i;
    jCur = j;
    bCur = a;
    curJoin = startJoin;
    firstTime = true;
    /*Now run the contour, from aCur to bCur within iCur, jCur*/
    for(;;)
    {
	if (firstTime)
	{
	    firstTime = false;
	}
	else
	{
	    XYZ *temp;
	    /*Emit this joining onto the list*/
	    temp = NewXYZNode(iCur, jCur, bCur, contourValue);
	    temp -> next = first;
	    first = temp;
	}

	/*Delete this joining*/
	WIPEJOIN(joinings[iCur + jCur * iSize], curJoin);

	/*Find a next place to go*/
	switch(bCur)
	{
	    case 0:
		--jCur;
		match = 2;
		break;
	    case 1:
		++iCur;
		match = 3;
		break;
	    case 2:
		++jCur;
		match = 0;
		break;
	    case 3:
		--iCur;
		match = 1;
		break;
	}

	if (iCur < 0 || jCur < 0 || iCur >= iSize - 1 || jCur >= jSize - 1)
	{
	    /*Fallen off the edge*/
	    break;
	}

	/*Get the new joining and see if it matches up.  Try 0 first.*/
	joining = SELJOIN(joinings[iCur + jCur * iSize], 0);
	if (joining)
	{
	    if (SELEND(joining, 1) == match)
	    {
		aCur = SELEND(joining, 1);
		bCur = SELEND(joining, 0);
		curJoin = 0;
		continue;
	    }
	    if (SELEND(joining, 0) == match)
	    {
		aCur = SELEND(joining, 0);
		bCur = SELEND(joining, 1);
		curJoin = 0;
		continue;
	    }
	}

	/*OK, so now try one*/
	joining = SELJOIN(joinings[iCur + jCur * iSize], 1);
	if (joining)
	{
	    if (SELEND(joining, 1) == match)
	    {
		aCur = SELEND(joining, 1);
		bCur = SELEND(joining, 0);
		curJoin = 1;
		continue;
	    }
	    if (SELEND(joining, 0) == match)
	    {
		aCur = SELEND(joining, 0);
		bCur = SELEND(joining, 1);
		curJoin = 1;
		continue;
	    }
	}
	break;
    }

    /*Count the elements in the list*/
    runner = first;
    count = 0;
    while (runner)
    {
	++count;
	runner = runner -> next;
    }

    /*Make the vertices*/
    vertices = (VertexPtr) Alloc(count * sizeof(Vertex));
    if (vertices)
    {
	/*Run through the list, putting it on the vertices*/
	runner = first;
	for (k = 0; k < count; ++k)
	{
	    vertices[k] . position[0] = runner -> x;
	    vertices[k] . position[1] = runner -> y;
	    vertices[k] . position[2] = runner -> z;
	    if (reverse)
	    {
		vertices[k] . normal[0] = -runner -> nx;
		vertices[k] . normal[1] = -runner -> ny;
		vertices[k] . normal[2] = -runner -> nz;
	    }
	    else
	    {
		vertices[k] . normal[0] = runner -> nx;
		vertices[k] . normal[1] = runner -> ny;
		vertices[k] . normal[2] = runner -> nz;
	    }
	    vertices[k] . colorIndex = 0;
	    runner = runner -> next;
	}

	/*Make a contour line*/
	AppendPolylineToPicture(surface, lineWidth, 2, count, vertices);

	/*Free the vertices*/
	Free(vertices);
    }

    /*Eliminate the XYZ list*/
    while (first)
    {
	last = first -> next;
	Free(first);
	first = last;
    }
}

static ObjPtr MakeContoursSurface(visObject)
ObjPtr visObject;
/*Makes the surface in a contours object.  Also colors it.*/
{
    ObjPtr dataset;		/*The dataset the vis object represents*/
    long datasetFlags;		/*Flags of the dataset*/
    ObjPtr var;			/*Random variable*/
    ObjPtr picture;		/*The picture to be made*/

    dataset = GetObjectVar("MakeContoursSurface", visObject, MAINDATASET);
    if (!dataset)
    {
	return ObjFalse;
    }

    datasetFlags = GetDatasetInfo(dataset);

    if (0 == datasetFlags & DS_HASFORM)
    {
	ReportError("MakeContoursSurface", "No data form");
	return ObjFalse;
    }


    /*Make the new picture*/
    picture = NewPicture();
    if (!picture) return ObjFalse;

    if (datasetFlags & DS_UNSTRUCTURED)
    {
	return ObjFalse;
    }
    else
    {
	/*It's a structured dataset.*/
	int topDim, nComponents;
	long curDim[2], iDims[2];
	ObjPtr dims;
	RectMeshPtr rectMesh;
	real temp[3];
	int iDim, jDim;
	Bool reverse;

	/*It must have dimension 2*/
	topDim = GetTopDim(dataset);
	if (topDim != 2)
	{
	    ReportError("MakeContoursSurface", "Topological dimension must be 2.");
	    return ObjFalse;
	}

	/*Get the actual topological dimensions*/
	dims = GetDatasetFormDims(dataset);
	if (!dims || !IsRealArray(dims) || RANK(dims) != 1 || DIMS(dims)[0] != 2)
	{
	    ReportError("MakeContoursSurface", "No topological dimensions");
	    return ObjFalse;
	}

	iDims[0] = ((real *) ELEMENTS(dims))[0];
	iDims[1] = ((real *) ELEMENTS(dims))[1];

	/*Register the dataset and its dataform*/
	if (!SetCurField(FIELD1, dataset))
	{
	    return ObjFalse;
	}
	if (!SetCurForm(FIELD2, dataset))
	{
	    return ObjFalse;
	}
	
	reverse = GetPredicate(visObject, REVERSESENSE) ^
		  GetPredicate(dataset, ISLEFTHANDED);
	if (reverse)
	{
	    iDim = 1;
	    jDim = 0;
	}
	else
	{
	    iDim = 0;
	    jDim = 1;
	}

	var = GetVar(visObject, CONTOURCENTER);
	if (var)
	{
	    /*Create all the contour lines*/
	    ObjPtr minMax;
	    real *elements;
	    ThingListPtr runner;
	    double min, max, mid;
	    double curContour;
	    unsigned char *joinings;
	    long i, j, k, iPlus, jPlus;
	    int initState[2][2];
	    double test;
	    double center;		/*Center of contour*/
	    double step;		/*Plus and minus step*/
	    int nDivisions;		/*Number of divisions*/
	    double contourMin, contourMax;/*Min and max of contours*/
	    long contourIndex;

	    center = GetReal(var);

	    /*Get the step*/
	    var = GetRealVar("MakeContoursSurface", visObject, CONTOURSTEP);
	    if (!var)
	    {
		return ObjFalse;
	    }
	    step = GetReal(var);
	    var = GetIntVar("MakeContoursSurface", visObject, CONTOURDIV);
	    if (!var)
	    {
		return ObjFalse;
	    }
	    nDivisions = GetInt(var);
	    if (nDivisions > 0)
	    {
	        step /= nDivisions;
	    }
	    else
	    {
		nDivisions = 1;
	    }

	    /*Now the min and max*/
	    var = GetRealVar("MakeContoursSurface", visObject, CONTOURMAX);
	    if (!var)
	    {
		return ObjFalse;
	    }
	    contourMax = GetReal(var);
	    var = GetRealVar("MakeContoursSurface", visObject, CONTOURMIN);
	    if (!var)
	    {
		return ObjFalse;
	    }
	    contourMin = GetReal(var);

	    /*Get the min and max of the dataset*/
	    MakeVar(dataset, MINMAX);
	    minMax = GetVar(dataset, MINMAX);
	    if (minMax)
	    {
		elements = ELEMENTS(minMax);
		min = elements[0];
	 	max = elements[1];
	    }
	    else
	    {
		min = 0.0;
		max = 1.0;
	    }
	    mid = (min + max) * 0.5;

	    /*Make the joinings array*/
	    joinings = (unsigned char *) Alloc(iDims[0] * iDims[1]);
	    for (k = 0; k < iDims[0] * iDims[1]; ++k)
	    {
		joinings[k] = 0;
	    }

	    /*Clip min and max with contour min and max*/
	    min = MAX(min, contourMin);
	    max = MIN(max, contourMax);

	    /*Find the index of the first contour*/
	    contourIndex = (min - center) / ABS(step);
	    curContour = center + contourIndex * ABS(step);

	    /*Go through the contours, finding joinings*/	    
	    while (curContour <= max)
	    {
		/*Generate joinings for this contour*/
		for (j = 0; j < iDims[jDim] - 1; ++j)
		{
		    for (i = 0; i < iDims[iDim] - 1; ++i)
		    {
			/*Set up initState*/
			for (jPlus = 0; jPlus < 2; ++jPlus)
			{
			    for (iPlus = 0; iPlus < 2; ++iPlus)
			    {
				curDim[iDim] = i + iPlus;
				curDim[jDim] = j + jPlus;
				test = SelectFieldScalar(FIELD1, curDim);
				if (test == missingData)
				{
				    goto cellMissing;
				}
				initState[(iDim == 0) ? iPlus : jPlus]
					 [(jDim == 1) ? jPlus : iPlus] =
				    (curContour > mid) ?
					(test <= curContour ? 1 : 0) :
					(test <= curContour ? 0 : 1);
			    }
			}
			curDim[iDim] = i;
			curDim[jDim] = j;
			joinings[curDim[1] * iDims[0] + curDim[0]] =
				FindQuadJoinings(initState);
cellMissing:;
		    }
		}

		/*Follow all the contours*/
		for (j = 0; j < iDims[jDim] - 1; ++j)
		{
		    for (i = 0; i < iDims[iDim] - 1; ++i)
		    {
			int joining;
			curDim[iDim] = i;
			curDim[jDim] = j;
			/*Must have a joining involving the bottom edge*/
			joining = joinings[curDim[1] * iDims[0] + curDim[0]];
			if (joining & 0xF0)
			{
			    FollowContour(picture, joinings, curDim[0], curDim[1], iDims[0], iDims[1], 1, curContour, reverse ? false : true, 0);
			}
			joining = joinings[curDim[1] * iDims[0] + curDim[0]];
			if (joining & 0x0F)
			{
			    FollowContour(picture, joinings, curDim[0], curDim[1], iDims[0], iDims[1], 0, curContour, reverse ? false : true, 0);
			}
		    }  
		}

		/*Next contour*/
		++contourIndex;
		curContour = center + contourIndex * ABS(step);
	    }

	    Free(joinings);
	}
    }

    SetVar(visObject, SURFACE, picture);
    SetVar(picture, REPOBJ, visObject);
    return ObjTrue;
}

static ObjPtr AddContourControls(contour, panelContents)
ObjPtr contour, panelContents;
/*Adds controls appropriate to a contour object to panelContents*/
{
    ObjPtr titleBox, button, radio, var, corral, icon, name, contourField, mainDataset;
    ObjPtr textBox, defaultIcon;
    int width, left, top, bottom, right, mid;
    ObjPtr control;
    ObjPtr minMax;
    char number[50];
    real step, origin, min, max;
    int divisions;

    width = CWINWIDTH - 2 * CORRALBORDER - CWINCORRALWIDTH;

    /*Get the contour field*/
    contourField = GetObjectVar("AddContourControls", contour, MAINDATASET);
    if (!contourField) return ObjFalse;
    while (mainDataset = GetVar(contourField, MAINDATASET))
    {
	contourField = mainDataset;
    }

    /*Put in the contour corral at the top left*/
    left = MAJORBORDER;
    top = MAJORBORDER;
    corral = NewIconCorral(NULLOBJ,
			   left, left + ONECORRALWIDTH,
			   CWINHEIGHT - MAJORBORDER - ONECORRALHEIGHT,
			   CWINHEIGHT - MAJORBORDER, 0);
    SetVar(corral, SINGLECORRAL, ObjTrue);
    SetVar(corral, TOPDOWN, ObjTrue);
    SetVar(corral, NAME, NewString("Contour Field"));
    SetVar(corral, HELPSTRING,
	NewString("This corral shows the dataset that is being used to make \
the contour display.  The locations of contours and color cells are calculated \
using this field.  The color of color cells is calculated using the Color Field, \
available in the Color control set."));
    PrefixList(panelContents, corral);
    SetVar(corral, PARENT, panelContents);
    SetVar(corral, REPOBJ, contour);
    SetMethod(corral, DROPINCONTENTS, DropInMainDatasetCorral);

    /*Create the contour source text box*/
    textBox = NewTextBox(left, left + ONECORRALWIDTH, 
			 CWINHEIGHT - MAJORBORDER - ONECORRALHEIGHT - TEXTBOXSEP - TEXTBOXHEIGHT,
			 CWINHEIGHT - MAJORBORDER - ONECORRALHEIGHT - TEXTBOXSEP,
			 0, "Contour Field Text", "Contour Field");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);
    SetTextAlign(textBox, CENTERALIGN);

    top = CWINHEIGHT - 2 * MAJORBORDER - ONECORRALHEIGHT - TEXTBOXSEP - TEXTBOXHEIGHT;

    /*Put in an icon that represents the field*/
    name = GetVar(contourField, NAME);
    defaultIcon = GetVar(contourField, DEFAULTICON);
    if (defaultIcon)
    {
	icon = NewObject(defaultIcon, 0);
	SetVar(icon, NAME, name);
    }
    else
    {
	icon = NewIcon(0, 0, ICONQUESTION, GetString(name));
    }
    SetVar(icon, ICONLOC, NULLOBJ);
    SetVar(icon, REPOBJ, contourField);
    DropIconInCorral(corral, icon);

    var = GetRealVar("AddContourControls", contour, CONTOURSTEP);
    if (var)
    {
	step = GetReal(var);
    }
    else
    {
	SetVar(contour, CONTOURSTEP, NewReal(1.0));
	step = 1.0;
    }

    var = GetRealVar("AddContourControls", contour, CONTOURCENTER);
    if (var)
    {
	origin = GetReal(var);
    }
    else
    {
	SetVar(contour, CONTOURCENTER, NewReal(0.0));
	origin = 0.0;
    }

    var = GetRealVar("AddContourControls", contour, CONTOURMIN);
    if (var)
    {
	min = GetReal(var);
    }
    else
    {
	SetVar(contour, CONTOURMIN, NewReal(minusInf));
	min = minusInf;
    }

    var = GetRealVar("AddContourControls", contour, CONTOURMAX);
    if (var)
    {
	max = GetReal(var);
    }
    else
    {
	SetVar(contour, CONTOURMAX, NewReal(plusInf));
	max = plusInf;
    }

    var = GetIntVar("AddContourControls", contour, CONTOURDIV);
    if (var)
    {
	divisions = GetInt(var);
    }
    else
    {
	SetVar(contour, CONTOURDIV, NewInt(1));
	divisions = 1;
    }

    /*Make the origin*/
    textBox = TemplateTextBox(VisContoursTemplate, "Origin Text", 0, "Origin:");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);

    PrintNumber(number, origin); 
    textBox = TemplateTextBox(VisContoursTemplate, "Origin",
			 EDITABLE + WITH_PIT + ONE_LINE, number);
    SetVar(textBox, PARENT, panelContents);
    PrefixList(panelContents, textBox);
    SetTextAlign(textBox, RIGHTALIGN);
    AssocTextRealControlWithVar(textBox, contour, CONTOURCENTER, minusInf, plusInf, TR_NE_BOTTOM | TR_NE_TOP);
    SetVar(textBox, HELPSTRING,
	NewString("This number specifies the origin of the contours.  \
Contours are calculated starting at the origin, increasing and decreasing by the steps \
specified in the major and minor contours parameters, until the minimum or maximum is reached.\
For most ordinary contours that cover the entire visualization range, this does not \
need to be changed."));

    /*Make the minimum*/
    textBox = TemplateTextBox(VisContoursTemplate, "Minimum Text", 
				0, "Minimum:");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);

    PrintNumber(number, min); 
    textBox = TemplateTextBox(VisContoursTemplate, "Minimum",
			 EDITABLE + WITH_PIT + ONE_LINE, number);
    SetVar(textBox, PARENT, panelContents);
    PrefixList(panelContents, textBox);
    SetTextAlign(textBox, RIGHTALIGN);
    AssocTextRealControlWithVar(textBox, contour, CONTOURMIN, minusInf, plusInf, TR_NE_BOTTOM);
    SetVar(textBox, HELPSTRING,
	NewString("This number specifies the minimum contour value.  \
Contours are calculated starting at the origin, increasing and decreasing by the steps \
specified in the major and minor contours parameters, until the minimum or maximum is reached.  \
For most ordinary contours that cover the entire visualization range, this does not \
need to be changed.\n"));

    /*Make the maximum*/
    textBox = TemplateTextBox(VisContoursTemplate, "Maximum Text", 0, "Maximum:");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);

    PrintNumber(number, max); 
    textBox = TemplateTextBox(VisContoursTemplate, "Maximum",
			EDITABLE + WITH_PIT + ONE_LINE, number);
    SetVar(textBox, PARENT, panelContents);
    PrefixList(panelContents, textBox);
    SetTextAlign(textBox, RIGHTALIGN);
    AssocTextRealControlWithVar(textBox, contour, CONTOURMAX, minusInf, plusInf, TR_NE_BOTTOM);
    SetVar(textBox, HELPSTRING,
	NewString("This number specifies the maximum contour value.  \
Contours are calculated starting at the origin, increasing and decreasing by the steps \
specified in the major and minor contours parameters, until the minimum or maximum is reached.  \
For most ordinary contours that cover the entire visualization range, this does not \
need to be changed.\n"));

    /*Make the major contour controls*/
    textBox = TemplateTextBox(VisContoursTemplate, "Step Text", 
			 0, "Major Step:");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);

    PrintNumber(number, step); 
    textBox = TemplateTextBox(VisContoursTemplate, "Step", 
			 EDITABLE + WITH_PIT + ONE_LINE, number);
    SetVar(textBox, PARENT, panelContents);
    PrefixList(panelContents, textBox);
    SetTextAlign(textBox, RIGHTALIGN);
    AssocTextRealControlWithVar(textBox, contour, CONTOURSTEP, 0.0, plusInf, TR_NE_BOTTOM | TR_NE_TOP);
    SetVar(textBox, HELPSTRING,
	NewString("This number specifies the step size of the major contours.  \
It must be a real number greater than 0.  \
Contours are calculated starting at the origin, increasing and decreasing by the steps \
specified in the major and minor contours parameters, until the minimum or maximum is reached."));

    /*Make the minor contour controls*/
    textBox = TemplateTextBox(VisContoursTemplate, "Divisions Text", 
			 0, "Minor Divisions:");
    PrefixList(panelContents, textBox);
    SetVar(textBox, PARENT, panelContents);

    PrintNumber(number, (real) divisions); 
    textBox = TemplateTextBox(VisContoursTemplate, "Divisions", 
			 EDITABLE + WITH_PIT + ONE_LINE, number);
    SetVar(textBox, PARENT, panelContents);
    PrefixList(panelContents, textBox);
    SetTextAlign(textBox, RIGHTALIGN);
    AssocTextIntControlWithVar(textBox, contour, CONTOURDIV, 1.0, plusInf, TR_INT_ONLY | TR_NE_TOP);
    SetVar(textBox, HELPSTRING,
	NewString("This text box contains a number giving the number of divisions between \
major contours.  \
This must be an integer greater than zero.  For example, a value of 1 gives no minor divisions \
between major contours.  A value of 10 gives one major contour every 10 divisions.  \
Contours are calculated starting at the origin, increasing and decreasing by the steps \
specified in the major and minor contours parameters, until the minimum or maximum is reached."));

    return ObjTrue;
}

ObjPtr MakeContourLines(contour)
ObjPtr contour;
/*Makes some contour lines*/
{
    ObjPtr line, lines;
    SetVar(contour, CONTOURLINES, lines = NewList());

    return ObjTrue;
}

void InitContours()
/*Initializes the contour objects*/
{
    ObjPtr icon, color;

    /*Class for a contour line*/
    contourLineClass = NewObject(NULLOBJ, 0);
    color = NewRealArray(1, 3L);
    ((real *) ELEMENTS(color))[0] = 1.0;
    ((real *) ELEMENTS(color))[1] = 1.0;
    ((real *) ELEMENTS(color))[2] = 1.0;
    SetVar(contourLineClass, BASECOLOR, color);
    AddToReferenceList(contourLineClass);

    /*Class for a contour object*/
    contoursClass = NewObject(visDeformed, 0);
    AddToReferenceList(contoursClass);
    SetVar(contoursClass, NAME, NewString("Contours"));
    SetMethod(contoursClass, INITIALIZE, ContoursInit);
    SetVar(contoursClass, SHINVAL, NewReal(80.0));
    SetVar(contoursClass, SPECVAL, NewReal(0.2));
    SetVar(contoursClass, DEFAULTICON, icon = NewObject(visIcon, 0));
    SetMethod(contoursClass, CONTOURLINES, MakeContourLines);
    SetVar(contoursClass, DRAWSURFACE, ObjFalse);
    SetVar(contoursClass, DRAWWIREFRAME, ObjTrue);
    SetVar(icon, WHICHICON, NewInt(ICONCONTOURS));
    SetVar(icon, NAME, NewString("Contours"));
    SetVar(icon, HELPSTRING,
	NewString("This icon represents a contour visualization object.  \
The contour object shows contours and shaded color surfaces of 2-dimensional \
scalar fields defined over structured or nonstructured grids."));
    DeclareIndirectDependency(contoursClass, SURFACE, MAINDATASET, CHANGED);
    DeclareDependency(contoursClass, SURFACE, COLORCELLS);
    SetMethod(contoursClass, SURFACE, MakeContoursSurface);
    DeclareDependency(contoursClass, SURFACE, REVERSESENSE);
    DeclareDependency(contoursClass, SURFACE, CONTOURCENTER);
    DeclareDependency(contoursClass, SURFACE, CONTOURSTEP);
    DeclareDependency(contoursClass, SURFACE, CONTOURMAX);
    DeclareDependency(contoursClass, SURFACE, CONTOURDIV);
    DeclareDependency(contoursClass, SURFACE, CONTOURMIN);

    SetMethod(contoursClass, PICCOLORED, MakeContoursColored);
    SetMethod(contoursClass, SETMAINDATASET, SetContoursMainDataset);

    SetMethod(contoursClass, ADDCONTROLS, AddContourControls);
    icon = NewIcon(0, 0, ICONCONTOURS, "Contours");
    SetVar(contoursClass, CONTROLICON, icon);

    DefineVisMapping(DS_HASFORM | DS_HASFIELD, 2, 3, 1, contoursClass);
    DefineVisMapping(DS_HASFORM | DS_HASFIELD, 2, 2, 1, contoursClass);
    DefineVisMapping(DS_HASFORM | DS_HASFIELD | DS_UNSTRUCTURED, 2, 3, 1, contoursClass);
    DefineVisMapping(DS_HASFORM | DS_HASFIELD | DS_UNSTRUCTURED, 2, 2, 1, contoursClass);
}

void KillContours()
/*Kills the contours*/
{
    DeleteThing(contoursClass);
    DeleteThing(contourLineClass);
}