/* ContourView.m
   Contour Plot object with optional color fills.
 
     It is pretty easy to use.  If instantiated within IB, all you need to do
     to get a plot is to call the following two method to get a default
     plot.  If you understand the following method, you can use the
     ContourView object easily.

     - setCartesianGridData:(float *)f :(float)xmin :(float)xmax
				  :(float)ymin :(float)ymax
				  ofSize:(int)nx :(int)ny
				  withInterpolationTo:(int)n1x :(int)n1y;
 
    f[nx*ny] is a 1-d array containing 2-d grid data such that f[iy*nx+ix]
    contains the value at (ix, iy).
    Typicall, just 3 messages below will produce a contour plot with color
    fills.

    [contourView setCartesianGridData: fdata :1.0 :5.0 :1.0 :5.0
				ofSize: 20 :20
				withInterpolationTo: 50 :50];
    [contourView setFillEnable:YES];
    [contourView display];
 
 
 * ---------------------------------------------------------------------------
 * Version History:
 * V0.92 92-12-03 Izumi
 *    It now correctly does color fills via sorting of contour drawing order.
 * V0.91 92-05-24 Izumi
 *    Does OK color fills, printing, copying PS to paste board.
 * V0.90 92-05-18 Izumi Ohzawa, izumi@pinoko.berkeley.edu
 *    The initial version.
*/

#import "ContourView.h"
#import "contour.h"
#import "splin2.h"
#import <stdlib.h>		/* malloc */
#import <math.h>
#import <streams/streams.h>	/* NXOpenMemory() etc */
#import <appkit/Window.h>
#import <appkit/Pasteboard.h>
#import <appkit/color.h>
#import <appkit/graphics.h>	/* for definitions NX_DefaultDepth etc. */
#import <dpsclient/dpsclient.h>
#import <dpsclient/dpsNeXT.h>
#import <dpsclient/wraps.h>	/* imports psops.h too */

#define N6			6	/* expand 3 points on 4 sides to close all contours */
#define N3			3
#define N2			2
#define DEFAULTNCLEVELS		16	/* Default number of contour levels */


static float pdash[] = {4.0, 4.0 };
static float psolid[] = {0.0 };


@implementation ContourView

- initFrame:(NXRect *)nf
{
    self = [super initFrame:nf];
    ClearFlag = 0;
    doFill = NO;
    debug = NO;
    frameON = YES;
    frameLineWidth = 1.0;
    minNumPoints = 6;
    ndata = 0;
    xd = NULL;
    yd = NULL;
    fd = NULL;
    bipolar = 0;
    basevalue = 0.0;
    pXmin = bounds.origin.x;
    pYmin = bounds.origin.y+0.25;
    pXmax = bounds.size.width + bounds.origin.x -0.25;
    pYmax = bounds.size.height + bounds.origin.y -0.25;

    nclevels = DEFAULTNCLEVELS;
    cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
    contourList = NULL;
    numContours = 0;
    [self setDefaultContourAttributes];

    positiveColor = NXConvertRGBToColor(0.0, 0.6, 0.0);    // dim green
    negativeColor = NX_COLORRED;
    backgroundColor = NX_COLORWHITE;
    contourLineColor = NX_COLORBLACK;
    frameColor = NX_COLORBLACK;
    minorContourLineWidth = 1.0;
    minorContourLineWidth = 2.0;
    return self;
}


- free
{
    if(cA) free(cA);
    if(xd) free(xd);
    if(yd) free(yd);
    if(fd) free(fd);
    [super free];
    return self;
}

- (BOOL)acceptsFirstMouse { return YES; }
- (BOOL)acceptsFirstResponder { return YES; }


// This will change scale of plotting, but will not touch internally
// stored grid data: xd[], yd[] arrays.  May be used for zooming?
//
- setScaleXY:(float)xmin :(float)xmax :(float)ymin :(float)ymax
{
   xdmin = xmin;
   ydmin = ymin;
   ppxu = (pXmax - pXmin)/(xmax - xmin);	// scaling factor points per unit of X
   ppyu = (pYmax - pYmin)/(ymax - ymin);	// scaling factor points per unit of Y
   return self;
}



// Given array x[] with npts elements, returns min and max values in the array.
//
- findMinMax:(float *)x :(int)npts :(float *)vmin :(float *)vmax
{
int i;
float fmin= MAXFLOAT;
float fmax= MINFLOAT;
    for(i=0; i<npts; i++)
    {
	if(x[i] < fmin) fmin = x[i];
	if(x[i] > fmax) fmax = x[i];
    }
    *vmin = fmin;
    *vmax = fmax;
    return self;
}

- setDebugEnable:(BOOL)de
{
    debug = de;
    return self;
}

- setFillEnable:(BOOL)fe
{
    doFill = fe;
    return self;
}


// Set the number of contour levels.
// It will be forced to an EVEN number if not already.
// It will automatically set default contour level values.
// If customized contour levels are needed,
// Call - setContourLevelArray:(float *)caa :(int)nl instead.
//
- setNumberOfContourLevels:(int)nl
{    
    if(nl == nclevels) return self;	/* no change */
    if(cA) free((void *)cA);
    nclevels = nl/2*2;		// force it to be even
    cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
    [self setDefaultContourAttributes];
    return self;
}

// Set the range of contour level values used.
// This can be used to override default autoscaling.
// For real flexibility use -setContourLevelArray:: method.
//
- setMinMaxOfContourLevels:(float)min :(float)max
{
    fdmin = min;
    fdmax = max;
    [self setDefaultContourAttributes];
    return self;
}

- setContourLineColor:(NXColor)clc
{
    contourLineColor = clc;
    [self setDefaultContourAttributes];
    return self;
}

- setFrameColor:(NXColor)fc
{
    frameColor = fc;
    return self;
}

- setBackgroundColor:(NXColor)bc
{
    backgroundColor = bc;
    [self setDefaultContourAttributes];
    return self;
}

- setFillColors:(NXColor)pe :(NXColor)ne
{
    positiveColor = pe;
    negativeColor = ne;
    [self setDefaultContourAttributes];
    return self;
}


// This method allows customization of contour attributes to use
// overriding the default ones.
//
- setContourAttributeArray:(CntrAttribute *)caa :(int)nl;
{
// FIXME #### copy caa to cA.
    return self;
}

// Setup contour levels and other attributes to use based on the max and min of
// values in the data.  This will be called automatically when grid
// data are set.
//
- setDefaultContourAttributes
{
int i, ncl2;
float absfmax, step;
float rbase, gbase, bbase;	// RGB of background
float rp, gp, bp;		// RGB of positive extreme
float rn, gn, bn;		// RGB of negative extreme
float rpstep, gpstep, bpstep;
float rnstep, gnstep, bnstep;
float r, g, b;

    NXConvertColorToRGB(positiveColor, &rp, &gp, &bp);
    NXConvertColorToRGB(negativeColor, &rn, &gn, &bn);
    NXConvertColorToRGB(backgroundColor, &rbase, &gbase, &bbase);

    if(fabs(fdmax) > fabs(fdmin))
	absfmax = fabs(fdmax);
    else
	absfmax = fabs(fdmin);
    ncl2 = nclevels/2;
    if(fdmin<0.0 && fdmax>0.0)
    {
	/* bipolar data */
	bipolar = 1;		/* instance var flag */
	step = absfmax /(float)ncl2;
        rpstep = (rp - rbase)/(float)ncl2;
	gpstep = (gp - gbase)/(float)ncl2;
	bpstep = (bp - bbase)/(float)ncl2;
        rnstep = (rn - rbase)/(float)ncl2;
	gnstep = (gn - gbase)/(float)ncl2;
	bnstep = (bn - bbase)/(float)ncl2;
	for(i=0; i<ncl2; i++)
	{
	    /* positive */
	    cA[ncl2+i].dash = 0;
	    cA[ncl2+i].level = step*i + step/2.0;
	    r = rpstep*(i+1) + rbase;
	    g = gpstep*(i+1) + gbase;
	    b = bpstep*(i+1) + bbase;
	    cA[ncl2+i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
	    if(i) cA[ncl2+i].fillcolor_lo = cA[ncl2+i-1].fillcolor_hi;
	    else  cA[ncl2+i].fillcolor_lo = backgroundColor;
	    cA[ncl2+i].linecolor = contourLineColor;
	    /* negative */
	    cA[ncl2-1-i].dash = 1;
	    cA[ncl2-1-i].level = -(step*i + step/2.0);
	    r = rnstep*(i+1) + rbase;
	    g = gnstep*(i+1) + gbase;
	    b = bnstep*(i+1) + bbase;
	    cA[ncl2-1-i].fillcolor_lo = NXConvertRGBToColor(r, g, b);
	    if(i) cA[ncl2-1-i].fillcolor_hi = cA[ncl2-i].fillcolor_lo;
	    else  cA[ncl2-1-i].fillcolor_hi = backgroundColor;
	    cA[ncl2-1-i].linecolor = contourLineColor;
	}
    }
    else
    {
	/* monopolar data */
	bipolar = 0;		/* store flag in instance variable */
	step = (fdmax - fdmin)/(float)(nclevels+1);
        rpstep = (rp - rbase)/(float)nclevels;
	gpstep = (gp - gbase)/(float)nclevels;
	bpstep = (bp - bbase)/(float)nclevels;
	for(i=0; i<nclevels; i++) {
	    cA[i].dash = 0;
	    cA[i].level = fdmin + step*(i+1);
	    r = rpstep*(i+1) + rbase;
	    g = gpstep*(i+1) + gbase;
	    b = bpstep*(i+1) + bbase;
	    cA[i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
	    if(i) cA[i].fillcolor_lo = cA[i-1].fillcolor_hi;
	    else  cA[i].fillcolor_lo = backgroundColor;
	    cA[i].linecolor = contourLineColor;
	}
    }
    return self;
}



// ==========================================================================
// Set surface data for regular Cartesian grid with optional interpolation.
// This is probably the most important method of this object.
// Use this if you want to plot data on uniform XY grid.
// Input data f[] of size (nx,ny) will be interpolated into (n1x, n1y) for
//     smoother contour plot if desired.  If no interpolation is needed,
//     make (n1x, n1y) same as (nx, ny).  This will turn off interpolation.
// 
//     (int)  nx, ny            : index sizes of input data
//     (int) n1x, n1y		: size of interpolated grid
//     for (i = 0; i < nx*ny; i++) (float)f[i]   : f value at (x,y)
//     f[i] is a 1-D array.  f(x,y) is in f[x + nx*y]
//
//     (xmin,ymin) and (xmax,ymax) define the (X,Y) domain of the plot.
//
// We extend the final (interpolated) data 3 points on each border to
// let the contour algorithm close all contours.  1-point wide region copy
// the data of the original border, and 2 points margins outside that are
// set to the base value.  This way, open contour lines that would have
// terminated at a border will get closed outside the original domain.
// This extra region for contour closure outside the original domain is
// outside the clip path, thus will not show up in the final plot.
// ---------------------------------------------------------------------------
//
- setCartesianGridData:(float *)f :(float)xmin :(float)xmax
				  :(float)ymin :(float)ymax
				  ofSize:(int)nnx :(int)nny
				  withInterpolationTo:(int)n1x :(int)n1y
{
int i, ix, iy;
float xstep, ystep;
float *xt, *yt;	/* temporary input arrays for interpolation */
float *fd1;
float **f2, **fd2a;

    xdmin = xmin;
    xdmax = xmax;
    ydmin = ymin;
    ydmax = ymax;
    // do reallocation only if # of data points has changed.
    nx = n1x+N6;
    ny = n1y+N6;
    nx1 = n1x;
    ny1 = n1y;
    if(ndata != ((n1x+N6)*(n1y+N6)))
    {
        ndata = (n1x+N6)*(n1y+N6);		/* new grid size */
    	if(xd) free(xd);
	if(yd) free(yd);
	if(fd) free(fd);
    	xd = (float *)malloc((size_t)ndata*sizeof(float));
    	yd = (float *)malloc((size_t)ndata*sizeof(float));
    	fd = (float *)malloc((size_t)ndata*sizeof(float));
    }

    if(nnx == n1x && nny == n1y)		/* Same grid size indicates NO INTERPOLATION */
    {
	/* copy original domain to expanded domain while defining border areas */
	for(i=0; i<ndata; i++) {
	    ix = i % nx;
	    iy = i / nx;
	    if(ix >= N3 && ix < (nx-N3) && iy >= N3 && iy < (ny-N3))
		fd[i] = f[(iy-N3)*nnx+(ix-N3)];		/* original data */
	    else if(iy < N3 && ix >= N3 && ix <(nx-N3))
		fd[i] = f[ix-N3];
	    else if(iy >= (ny-N3) && ix >= N3 && ix <(nx-N3))
		fd[i] = f[(nny-1)*nnx+(ix-N3)];
	    else if(ix < N3 && iy >= N3 && iy <(ny-N3))
		fd[i] = f[(iy-N3)*nnx];
	    else if(ix >=(nx-N3) && iy >= N3 && iy <(ny-N3))
		fd[i] = f[(iy-N3)*nnx + nnx-1];
	    else if(ix < N3 && iy <N3)		/* lower left */
		fd[i] = f[0];
	    else if(ix < N3 && iy >= (ny-N3))	/* uppper left */
		fd[i] = f[(nny-1)*nnx];
	    else if(ix >= (nx-N3) && iy < N3)	/* lower right */
		fd[i] = f[nnx-1];
	    else if(ix >= (nx-N3) && iy >= (ny-N3))
		fd[i] = f[(nny-1)*nnx + nnx-1];

	    /* wipe 2 pixel border to base level */
	    if(ix < N2 || ix >= (nx-N2) || iy <N2 || iy >= (ny-N2))
		fd[i] = basevalue;		/* totally out */
	}
    }
    else		/* Interpolate (nx, ny) to new size grid (n1x, n1y) */
    {
	// Allocate temporary input arrays for interpolation
	xt = (float *)malloc((size_t)nnx*sizeof(float));
	yt = (float *)malloc((size_t)nny*sizeof(float));
	fd1 =  (float *)malloc((size_t)nx1*ny1*sizeof(float)); /* for interpolation */
	f2 = fmatrix(1, nnx, 1, nny);
	fd2a = fmatrix(1, nnx, 1, nny);
	xstep = (xdmax-xdmin)/(float)(nnx-1);
	ystep = (ydmax-ydmin)/(float)(nny-1);
	for(ix=0; ix<nnx; ix++)
	    xt[ix] = xstep*ix + xdmin;
	for(iy=0; iy<nny; iy++)
	    yt[iy] = ystep*iy + ydmin;
	/* copy original data to 2-D array for interpolation */
	for(ix=0; ix <nnx; ix++)
	    for(iy=0; iy <nny; iy++)
		f2[ix+1][iy+1] = f[iy*nnx+ix];

	/* precompute second-derivative arrays for splin2()
	 * 2-nd derivative is returned in fd2a[][].
	 * (xt-1) and (yt-1) are for passing C's zero-offset array to
	 * splie2() which expects 1-offset arrays. */
	splie2(xt-1, yt-1, f2, nnx, nny, fd2a);

	// Interpolate f(nnx,nny) to finer array fd1(n1x,n1y)
    	xstep = (xdmax-xdmin)/(float)(nx1-1);
	ystep = (ydmax-ydmin)/(float)(ny1-1);
	for(ix=0; ix <nx1; ix++)
	    for(iy=0; iy <ny1; iy++)
		splin2(xt-1, yt-1, f2, fd2a, nnx, nny,
			xstep*ix+xdmin, ystep*iy+ydmin, &fd1[iy*nx1+ix]);

	/* Copy interpolated array to expanded array while seting border values. */
	for(i=0; i<ndata; i++) {
	    ix = i % nx;
	    iy = i / nx;
	    if(ix >= N3 && ix < (nx-N3) && iy >= N3 && iy < (ny-N3))
		fd[i] = fd1[(iy-N3)*n1x+(ix-N3)];		/* copy interpolated data */
	    else if(iy < N3 && ix >= N3 && ix <(nx-N3))    	/* bottom */
		fd[i] = fd1[ix-N3];
	    else if(iy >= (ny-N3) && ix >= N3 && ix <(nx-N3))  	/* top */
		fd[i] = fd1[(n1y-1)*n1x+(ix-N3)];
	    else if(ix < N3 && iy >= N3 && iy <(ny-N3))		/* left */
		fd[i] = fd1[(iy-N3)*n1x];
	    else if(ix >=(nx-N3) && iy >= N3 && iy <(ny-N3))	/* right */
		fd[i] = fd1[(iy-N3)*n1x + n1x-1];
	    else if(ix < N3 && iy <N3)		/* lower left */
		fd[i] = fd1[0];
	    else if(ix < N3 && iy >= (ny-N3))	/* uppper left */
		fd[i] = fd1[(n1y-1)*n1x];
	    else if(ix >= (nx-N3) && iy < N3)	/* lower right */
		fd[i] = fd1[n1x-1];
	    else if(ix >= (nx-N3) && iy >= (ny-N3))
		fd[i] = fd1[(n1y-1)*n1x + n1x-1];

	    /* wipe 2 pixel border to base level */
	    if(ix < N2 || ix >= (nx-N2) || iy <N2 || iy >= (ny-N2))
		fd[i] = basevalue;		/* totally out */
	}

        // Free temporary input array.
	if(f2) free_fmatrix(f2, 1, nnx, 1, nny);
	if(fd2a) free_fmatrix(fd2a, 1, nnx, 1, nny);
	if(fd1) free((void *)fd1);
        if(xt) free((void *)xt);
        if(yt) free((void *)yt);
    } /* END OF: if(nnx == n1x && nny == n1y) {} else {} */

    // Generate proper xd[], yd[] array for countour_().
    xstep = (xdmax-xdmin)/(float)(nx1-1);
    ystep = (ydmax-ydmin)/(float)(ny1-1);
    for(iy=0; iy<ny; iy++)
	for(ix=0; ix<nx; ix++)
	{
	    xd[ix+iy*nx] = xstep*(ix-N3) + xdmin;
	    yd[ix+iy*nx] = ystep*(iy-N3) + ydmin;
	}

    // Do appropriate scaling
    [self setScaleXY:xdmin :xdmax :ydmin :ydmax];

    // Get min and max fd[]
    [self findMinMax:fd :ndata :&fdmin :&fdmax];
    [self setDefaultContourAttributes];

    return self;
}



// Set XY grid points and values at each grid points.
// Grid does not have to be Cartesian or regular.
// Use this routine if you need flexibility.  No interpolation of
// grid is provided.  If you want interpolation,
// you have to do that yourself external to this object.
//     (int)  nj, nk                       : index sizes of data
//     for (i = 0; i < nj*nk; i++) (float)x[i]   : x coordinate
//     for (i = 0; i < nj*nk; i++) (float)y[i]   : y coordinate
//     for (i = 0; i < nj*nk; i++) (float)f[i]   : f value at (x[i], y[i])}
//     f[i] is a 1-D array.  f(j,k) is in f[j + nj*k].  Same for x[], y[].
//
// NOTE: This methond does not support contour closure by expanding the domain.
//
- setGridAndValueData:(float *)x :(float *)y :(float *)f ofSize:(int)nj :(int)nk
{
int i;
    // do reallocation only if # of data points has changed.
    nx = nj;
    ny = nk;
    if(ndata != (nj*nk))
    {
        ndata = nj*nk;
    	if(xd) free(xd);
	if(yd) free(yd);
	if(fd) free(fd);
    	xd = (float *)malloc((size_t)ndata*sizeof(float));
    	yd = (float *)malloc((size_t)ndata*sizeof(float));
    	fd = (float *)malloc((size_t)ndata*sizeof(float));
    }
    for(i=0; i<ndata; i++)	// copy data into instance variables
    {
	xd[i] = x[i];
	yd[i] = y[i];
	fd[i] = f[i];
    }

    // Do appropriate scaling
    // Get min and max of xd[], yd[], and fd[]
    [self findMinMax:xd :ndata :&xdmin :&xdmax];
    [self findMinMax:yd :ndata :&ydmin :&ydmax];
    [self findMinMax:fd :ndata :&fdmin :&fdmax];
    [self setScaleXY:xdmin :xdmax :ydmin :ydmax];
    [self setDefaultContourAttributes];

    return self;
}



// Controls width and whether a rectangular frame is drawn.
//
- setFrameEnable:(BOOL)fflag lineWidth:(float)fw
{

    frameON= fflag;	// Draw bounds frame
    frameLineWidth=fw;
    return self;
}

- setContourLineWidthMinor:(float)clw andMajor:(float)clwm
{
    minorContourLineWidth = clw;
    majorContourLineWidth = clwm;
    return self;
}




// Contours with number of points less than nmp are not plotted.
// This is used to eliminate tiny "speckles" in contour plot.
//
- setMinNumberOfPointsPerContour:(int)mnp
{
    minNumPoints = mnp;
    return self;
}

- clear:sender
{
    ClearFlag = 1;
    [self display];
    ClearFlag = 0;
    return self;
}

- copy:sender
{
    [self copyPScode:sender];
    return self;
}

// Copies the current view into the Pasteboard as PostScript.
//
- copyPScode:sender
{
    NXStream *psStream;
    id        pb;
    char     *data;
    int       dataLen, maxDataLen;

  /* Open a stream on memory where we will collect the PostScript */
    psStream = NXOpenMemory(NULL, 0, NX_WRITEONLY);
    if (!psStream)
      return self;

  /* Tell the Pasteboard we're going to copy PostScript */
    pb = [Pasteboard new];
    [pb declareTypes:&NXPostScriptPboardType num:1 owner:self];

  /* writes the PostScript for the whole plot as EPS into the stream */
    [self copyPSCodeInside:&bounds to:psStream];

  /* get the buffered up PostScript out of the stream */
    NXGetMemoryBuffer(psStream, &data, &dataLen, &maxDataLen);

  /* put the buffer in the Pasteboard, free the stream (and the buffer) */
    [pb writeType:NXPostScriptPboardType data:data length:dataLen];
    NXCloseMemory(psStream, NX_FREEBUFFER);
    return self;
}

// Call back method messaged from computeContour() function.
- accumContour:(int)icont :(int)np :(float *)x :(float *)y
{
ContourPath *newcntr;	/* scratch pad pointer to new contourList items */
int c_closed = 1;

    if(np < minNumPoints) return self;		/* too few points for contour */

    if (fabs(x[0] - x[np-1]) < 0.0001 && fabs(y[0] - y[np-1]) < 0.0001) {
	c_closed = 1;
	x[np-1] = x[0];  y[np-1] = y[0];	/* guarantee closure */
    }
    else
	c_closed = 0;			/* contour not closed */

    /* ---- Save new contour in a linked list ----------------------------------- */
    /* Allocate new item for contourList */
    newcntr = (ContourPath *)malloc((size_t)sizeof(ContourPath));
    newcntr->num_pts = np;		/* # of (x,y) points in contour */
    newcntr->closed = c_closed;		/* contour is closed flag */
    newcntr->levelindex = icont;	/* levelindex-th contour (indx to contour level) */
    newcntr->level = cA[icont].level;	/* value of contour */
    /* Allocate arrays for X and Y coordinates.
	* 7 more points may be needed to close unclosed contours with a path
	* that goes around the domain.
	*/
    newcntr->x = (float *)malloc((size_t)(np*sizeof(float)));
    newcntr->y = (float *)malloc((size_t)(np*sizeof(float)));
    /* copy (x,y) coordinates of contour path */
    memcpy((void *)newcntr->x, x, np * sizeof(float));
    memcpy((void *)newcntr->y, y, np * sizeof(float));
    /* insert new item into the list */
    newcntr->next = contourList;	/* point to previous one as next */
    contourList = newcntr;		/* for next time */
    numContours++;			/* increment number of contours counter */
    /* ---- New contour saved in a linked list ----------------------------------- */
    return self;
}


// This will override the default -drawSelf::, which does no real drawing.
// Should never be called directly.  It will be called by [self display];
//
- drawSelf:(const NXRect *)r:(int)count
{
int j;
ContourPath *cntr;
DPSContext curContext = DPSGetCurrentContext();

    if(NXDrawingStatus != NX_DRAWING)
        DPSPrintf(curContext, "\n%% Start of ContourView PostScript drawing..\n");
    NXSetColor(backgroundColor);
    NXRectFill(&bounds);

    if(!fd || ClearFlag) return self;	// if no data set, do nothing.

    if(contourList) [self freeContourList];		/* to start out */
    numContours = 0;


    computeContour(self, nx, ny, xd, yd, fd, nclevels, cA);

    if(debug)
	fprintf(stderr, "ContourView: # levels=%d,  # contours=%d\n", nclevels, numContours);

    if(doFill) {
	/* Sort the order of contour drawing from outermost to innermost for filling */
        SortedCntrPtr = (ContourPath **)(malloc((size_t)(sizeof(ContourPath *)*numContours)));
        sort_contourList(contourList, xdmin, xdmax, ydmin, ydmax, numContours, SortedCntrPtr);
	for(j=0; j<numContours; j++) {
	    cntr = SortedCntrPtr[j];
	    [self findInsideHighLow:cntr];		/* inside contour goin up or down */
	    [self plotContour:cntr :curContext];	/* This plots one contour curve */
	}
    }
    else {
	/* Do NOT do filling -- just go through the list */
	cntr = contourList;
	while(cntr) {
	    [self plotContour:cntr :curContext];	/* This plots one contour curve */
	    cntr = cntr->next;		/* point to the next one */
	};
    }

    if(frameON) {
        if(NXDrawingStatus != NX_DRAWING)
            DPSPrintf(curContext, "\n%% Frame rectangle\n");
	// Frame rectangle
	NXSetColor(frameColor);
	PSsetdash(psolid, 0, 0.0);
	PSsetlinewidth(frameLineWidth);
	PSrectstroke(pXmin,pYmin,pXmax,pYmax);	/* frame rectangle */
    }

    if(doFill) free(SortedCntrPtr);
    [self freeContourList];

    if(NXDrawingStatus != NX_DRAWING)
        DPSPrintf(curContext, "\n%% End of CurveView drawing.\n\n");

    return self;
}



- plotContour:(ContourPath *)cntr :(DPSContext)cContext
{
int i, np;
float xp, yp;
static char *ocmsg[] = {"open", "closed",""};

    np = cntr->num_pts;		/* # points in this contour */
    if(NXDrawingStatus != NX_DRAWING)
        DPSPrintf(cContext, "\n%% ### Contour: level index=%d, value=%g, #pts=%d, [%s]\n",
			cntr->levelindex, cntr->level, np, ocmsg[(cntr->closed)?1:0]);

    for(i=0; i<np; i++) {
	xp = (cntr->x[i] - xdmin)*ppxu;
	yp = (cntr->y[i] - ydmin)*ppyu;
	if(i==0)
	    PSmoveto(xp, yp);
	else
	    PSlineto(xp, yp);
    }
    if(cntr->closed)
	    PSclosepath();

    if(doFill) {
	if(cntr->hi_inside)
	    NXSetColor(cA[cntr->levelindex].fillcolor_hi);
	else
	    NXSetColor(cA[cntr->levelindex].fillcolor_lo);
	PSgsave();
	PSfill();
	PSgrestore();
    }

    if(cA[cntr->levelindex].dash)
	PSsetdash(pdash, 2, 0.0);
	/* fprintf(fpp, "[4 4] 0 setdash\n"); */ /* dashed curve */
    else
	PSsetdash(psolid, 0, 0.0);
	/* fprintf(fpp, "[] 0 setdash\n");	*/	/* solid curve */
    NXSetColor(cA[cntr->levelindex].linecolor);
    PSstroke();
    return self;
}

- findInsideHighLow:(ContourPath *)cntr
{
int i, i3, j, ii=0, ix=0, iy=0, npp, inside=0, gpfound=0, jx=0, jy=0;
float xstep, ystep, xg=0.0, yg=0.0, fg= 0.0;
static int ipx[] = { 0, -1,  1,  0, -1,  1,  0, -1,  1 };
static int ipy[] = { 0,  0,  0, -1, -1, -1,  1,  1,  1 };

    xstep = (xdmax-xdmin)/(float)(nx1-1);
    ystep = (ydmax-ydmin)/(float)(ny1-1);
    npp = cntr->num_pts;

    // Typical setting for hi_inside flag
    if(cntr->level > 0.0) cntr->hi_inside = 1;
    else		  cntr->hi_inside = 0;

    // Now try to determine it exactly.  Start with i=3 in an attempt to get
    // a good point that is not on the edge.
    for(i3=3; i3<(npp+3); i3++) {
	i = i3 % npp;		/* i = [3, .. (npp-1), 0, 1, 2] */
	if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
		continue;		// this contour point is outside domain
	ix  = (cntr->x[i] - xdmin)/xstep + 0.5;	// indices of closest grid point
	iy  = (cntr->y[i] - ydmin)/ystep + 0.5;
	// find indices (jx, jy) which is inside contour and inside domain
	// by checking 3x3 points centered on (ix, iy)
	for(j=0; j<9; j++) {
	    jx = ix + ipx[j];	
	    jy = iy + ipy[j];
	    xg = xstep*(float)jx + xdmin;
	    yg = ystep*(float)jy + ydmin;
	    if([self pointInDomain:xg :yg] && non_zero_winding(xg, yg, cntr->x, cntr->y, npp)==2)
	    {
		inside = 1;
		break;		// out of for(j..)
	    }
	}
	if(!inside) continue;	// try another point on contour

	ii = i;
	gpfound = 1;
	break;				// good point found.
    }  /* end of for(i=0...) */

    // If inside point not found, use closed point that is outside
    if(!inside) {
	for(i=0; i<npp; i++) {
	    if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
		    continue;		// this contour point is outside domain
	    jx  = (cntr->x[i] - xdmin)/xstep + 0.5;	// indices of closest grid point
	    jy  = (cntr->y[i] - ydmin)/ystep + 0.5;
	    xg = xstep*(float)jx + xdmin;
	    yg = ystep*(float)jy + ydmin;
	    if( ![self pointInDomain:xg :yg])	// at least it must be inside domain
		continue;

	    ii = i;
	    gpfound = 1;
	    if(non_zero_winding(xg, yg, cntr->x, cntr->y, npp) == 2) inside = 1;
	    break;
	}
    }

    fg = fd[(jy+N3)*nx+jx+N3];		// grid value
    if(inside) {
	if(fg > cntr->level) cntr->hi_inside = 1;	// inside is high
	else		     cntr->hi_inside = 0;
    } else {
	if(fg > cntr->level) cntr->hi_inside = 0;	// outside is high
	else		     cntr->hi_inside = 1;
    }

    if(debug && !gpfound)
	fprintf(stderr,
		"ContourView: Can't determine hi_inside: idx=%d, val=%g, #pt=%d, hi_ins=%d\n",
			cntr->levelindex, cntr->level, npp, cntr->hi_inside);
    return self;
}

- (BOOL)pointInDomain:(float)xx :(float)yy
{
    if(xx >= xdmin && xx <= xdmax && yy >= ydmin && yy <= ydmax)
	return YES;
    else
	return NO;
}


- freeContourList;
{
ContourPath *cntr;
	while((cntr = contourList))
	{
	    if(cntr->x) free(cntr->x);
	    if(cntr->y) free(cntr->y);
	    contourList = cntr->next;
	    free(cntr);
	};
        contourList = NULL;
	return self;
}


@end


