This is ContourView.m in view mode; [Download] [Up]
/* 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
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.