home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Nebula
/
nebula.bin
/
SourceCode
/
Classes
/
ContourView
/
ContourView.m
< prev
next >
Wrap
Text File
|
1992-12-03
|
26KB
|
858 lines
/* 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