home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1995 August / NEBULA.mdf / SourceCode / Classes / ContourView / ContourView.m < prev    next >
Encoding:
Text File  |  1992-12-02  |  25.2 KB  |  858 lines
  1. /* ContourView.m
  2.    Contour Plot object with optional color fills.
  3.  
  4.      It is pretty easy to use.  If instantiated within IB, all you need to do
  5.      to get a plot is to call the following two method to get a default
  6.      plot.  If you understand the following method, you can use the
  7.      ContourView object easily.
  8.  
  9.      - setCartesianGridData:(float *)f :(float)xmin :(float)xmax
  10.                   :(float)ymin :(float)ymax
  11.                   ofSize:(int)nx :(int)ny
  12.                   withInterpolationTo:(int)n1x :(int)n1y;
  13.  
  14.     f[nx*ny] is a 1-d array containing 2-d grid data such that f[iy*nx+ix]
  15.     contains the value at (ix, iy).
  16.     Typicall, just 3 messages below will produce a contour plot with color
  17.     fills.
  18.  
  19.     [contourView setCartesianGridData: fdata :1.0 :5.0 :1.0 :5.0
  20.                 ofSize: 20 :20
  21.                 withInterpolationTo: 50 :50];
  22.     [contourView setFillEnable:YES];
  23.     [contourView display];
  24.  
  25.  
  26.  * ---------------------------------------------------------------------------
  27.  * Version History:
  28.  * V0.92 92-12-03 Izumi
  29.  *    It now correctly does color fills via sorting of contour drawing order.
  30.  * V0.91 92-05-24 Izumi
  31.  *    Does OK color fills, printing, copying PS to paste board.
  32.  * V0.90 92-05-18 Izumi Ohzawa, izumi@pinoko.berkeley.edu
  33.  *    The initial version.
  34. */
  35.  
  36. #import "ContourView.h"
  37. #import "contour.h"
  38. #import "splin2.h"
  39. #import <stdlib.h>        /* malloc */
  40. #import <math.h>
  41. #import <streams/streams.h>    /* NXOpenMemory() etc */
  42. #import <appkit/Window.h>
  43. #import <appkit/Pasteboard.h>
  44. #import <appkit/color.h>
  45. #import <appkit/graphics.h>    /* for definitions NX_DefaultDepth etc. */
  46. #import <dpsclient/dpsclient.h>
  47. #import <dpsclient/dpsNeXT.h>
  48. #import <dpsclient/wraps.h>    /* imports psops.h too */
  49.  
  50. #define N6            6    /* expand 3 points on 4 sides to close all contours */
  51. #define N3            3
  52. #define N2            2
  53. #define DEFAULTNCLEVELS        16    /* Default number of contour levels */
  54.  
  55.  
  56. static float pdash[] = {4.0, 4.0 };
  57. static float psolid[] = {0.0 };
  58.  
  59.  
  60. @implementation ContourView
  61.  
  62. - initFrame:(NXRect *)nf
  63. {
  64.     self = [super initFrame:nf];
  65.     ClearFlag = 0;
  66.     doFill = NO;
  67.     debug = NO;
  68.     frameON = YES;
  69.     frameLineWidth = 1.0;
  70.     minNumPoints = 6;
  71.     ndata = 0;
  72.     xd = NULL;
  73.     yd = NULL;
  74.     fd = NULL;
  75.     bipolar = 0;
  76.     basevalue = 0.0;
  77.     pXmin = bounds.origin.x;
  78.     pYmin = bounds.origin.y+0.25;
  79.     pXmax = bounds.size.width + bounds.origin.x -0.25;
  80.     pYmax = bounds.size.height + bounds.origin.y -0.25;
  81.  
  82.     nclevels = DEFAULTNCLEVELS;
  83.     cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
  84.     contourList = NULL;
  85.     numContours = 0;
  86.     [self setDefaultContourAttributes];
  87.  
  88.     positiveColor = NXConvertRGBToColor(0.0, 0.6, 0.0);    // dim green
  89.     negativeColor = NX_COLORRED;
  90.     backgroundColor = NX_COLORWHITE;
  91.     contourLineColor = NX_COLORBLACK;
  92.     frameColor = NX_COLORBLACK;
  93.     minorContourLineWidth = 1.0;
  94.     minorContourLineWidth = 2.0;
  95.     return self;
  96. }
  97.  
  98.  
  99. - free
  100. {
  101.     if(cA) free(cA);
  102.     if(xd) free(xd);
  103.     if(yd) free(yd);
  104.     if(fd) free(fd);
  105.     [super free];
  106.     return self;
  107. }
  108.  
  109. - (BOOL)acceptsFirstMouse { return YES; }
  110. - (BOOL)acceptsFirstResponder { return YES; }
  111.  
  112.  
  113. // This will change scale of plotting, but will not touch internally
  114. // stored grid data: xd[], yd[] arrays.  May be used for zooming?
  115. //
  116. - setScaleXY:(float)xmin :(float)xmax :(float)ymin :(float)ymax
  117. {
  118.    xdmin = xmin;
  119.    ydmin = ymin;
  120.    ppxu = (pXmax - pXmin)/(xmax - xmin);    // scaling factor points per unit of X
  121.    ppyu = (pYmax - pYmin)/(ymax - ymin);    // scaling factor points per unit of Y
  122.    return self;
  123. }
  124.  
  125.  
  126.  
  127. // Given array x[] with npts elements, returns min and max values in the array.
  128. //
  129. - findMinMax:(float *)x :(int)npts :(float *)vmin :(float *)vmax
  130. {
  131. int i;
  132. float fmin= MAXFLOAT;
  133. float fmax= MINFLOAT;
  134.     for(i=0; i<npts; i++)
  135.     {
  136.     if(x[i] < fmin) fmin = x[i];
  137.     if(x[i] > fmax) fmax = x[i];
  138.     }
  139.     *vmin = fmin;
  140.     *vmax = fmax;
  141.     return self;
  142. }
  143.  
  144. - setDebugEnable:(BOOL)de
  145. {
  146.     debug = de;
  147.     return self;
  148. }
  149.  
  150. - setFillEnable:(BOOL)fe
  151. {
  152.     doFill = fe;
  153.     return self;
  154. }
  155.  
  156.  
  157. // Set the number of contour levels.
  158. // It will be forced to an EVEN number if not already.
  159. // It will automatically set default contour level values.
  160. // If customized contour levels are needed,
  161. // Call - setContourLevelArray:(float *)caa :(int)nl instead.
  162. //
  163. - setNumberOfContourLevels:(int)nl
  164. {    
  165.     if(nl == nclevels) return self;    /* no change */
  166.     if(cA) free((void *)cA);
  167.     nclevels = nl/2*2;        // force it to be even
  168.     cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
  169.     [self setDefaultContourAttributes];
  170.     return self;
  171. }
  172.  
  173. // Set the range of contour level values used.
  174. // This can be used to override default autoscaling.
  175. // For real flexibility use -setContourLevelArray:: method.
  176. //
  177. - setMinMaxOfContourLevels:(float)min :(float)max
  178. {
  179.     fdmin = min;
  180.     fdmax = max;
  181.     [self setDefaultContourAttributes];
  182.     return self;
  183. }
  184.  
  185. - setContourLineColor:(NXColor)clc
  186. {
  187.     contourLineColor = clc;
  188.     [self setDefaultContourAttributes];
  189.     return self;
  190. }
  191.  
  192. - setFrameColor:(NXColor)fc
  193. {
  194.     frameColor = fc;
  195.     return self;
  196. }
  197.  
  198. - setBackgroundColor:(NXColor)bc
  199. {
  200.     backgroundColor = bc;
  201.     [self setDefaultContourAttributes];
  202.     return self;
  203. }
  204.  
  205. - setFillColors:(NXColor)pe :(NXColor)ne
  206. {
  207.     positiveColor = pe;
  208.     negativeColor = ne;
  209.     [self setDefaultContourAttributes];
  210.     return self;
  211. }
  212.  
  213.  
  214. // This method allows customization of contour attributes to use
  215. // overriding the default ones.
  216. //
  217. - setContourAttributeArray:(CntrAttribute *)caa :(int)nl;
  218. {
  219. // FIXME #### copy caa to cA.
  220.     return self;
  221. }
  222.  
  223. // Setup contour levels and other attributes to use based on the max and min of
  224. // values in the data.  This will be called automatically when grid
  225. // data are set.
  226. //
  227. - setDefaultContourAttributes
  228. {
  229. int i, ncl2;
  230. float absfmax, step;
  231. float rbase, gbase, bbase;    // RGB of background
  232. float rp, gp, bp;        // RGB of positive extreme
  233. float rn, gn, bn;        // RGB of negative extreme
  234. float rpstep, gpstep, bpstep;
  235. float rnstep, gnstep, bnstep;
  236. float r, g, b;
  237.  
  238.     NXConvertColorToRGB(positiveColor, &rp, &gp, &bp);
  239.     NXConvertColorToRGB(negativeColor, &rn, &gn, &bn);
  240.     NXConvertColorToRGB(backgroundColor, &rbase, &gbase, &bbase);
  241.  
  242.     if(fabs(fdmax) > fabs(fdmin))
  243.     absfmax = fabs(fdmax);
  244.     else
  245.     absfmax = fabs(fdmin);
  246.     ncl2 = nclevels/2;
  247.     if(fdmin<0.0 && fdmax>0.0)
  248.     {
  249.     /* bipolar data */
  250.     bipolar = 1;        /* instance var flag */
  251.     step = absfmax /(float)ncl2;
  252.         rpstep = (rp - rbase)/(float)ncl2;
  253.     gpstep = (gp - gbase)/(float)ncl2;
  254.     bpstep = (bp - bbase)/(float)ncl2;
  255.         rnstep = (rn - rbase)/(float)ncl2;
  256.     gnstep = (gn - gbase)/(float)ncl2;
  257.     bnstep = (bn - bbase)/(float)ncl2;
  258.     for(i=0; i<ncl2; i++)
  259.     {
  260.         /* positive */
  261.         cA[ncl2+i].dash = 0;
  262.         cA[ncl2+i].level = step*i + step/2.0;
  263.         r = rpstep*(i+1) + rbase;
  264.         g = gpstep*(i+1) + gbase;
  265.         b = bpstep*(i+1) + bbase;
  266.         cA[ncl2+i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
  267.         if(i) cA[ncl2+i].fillcolor_lo = cA[ncl2+i-1].fillcolor_hi;
  268.         else  cA[ncl2+i].fillcolor_lo = backgroundColor;
  269.         cA[ncl2+i].linecolor = contourLineColor;
  270.         /* negative */
  271.         cA[ncl2-1-i].dash = 1;
  272.         cA[ncl2-1-i].level = -(step*i + step/2.0);
  273.         r = rnstep*(i+1) + rbase;
  274.         g = gnstep*(i+1) + gbase;
  275.         b = bnstep*(i+1) + bbase;
  276.         cA[ncl2-1-i].fillcolor_lo = NXConvertRGBToColor(r, g, b);
  277.         if(i) cA[ncl2-1-i].fillcolor_hi = cA[ncl2-i].fillcolor_lo;
  278.         else  cA[ncl2-1-i].fillcolor_hi = backgroundColor;
  279.         cA[ncl2-1-i].linecolor = contourLineColor;
  280.     }
  281.     }
  282.     else
  283.     {
  284.     /* monopolar data */
  285.     bipolar = 0;        /* store flag in instance variable */
  286.     step = (fdmax - fdmin)/(float)(nclevels+1);
  287.         rpstep = (rp - rbase)/(float)nclevels;
  288.     gpstep = (gp - gbase)/(float)nclevels;
  289.     bpstep = (bp - bbase)/(float)nclevels;
  290.     for(i=0; i<nclevels; i++) {
  291.         cA[i].dash = 0;
  292.         cA[i].level = fdmin + step*(i+1);
  293.         r = rpstep*(i+1) + rbase;
  294.         g = gpstep*(i+1) + gbase;
  295.         b = bpstep*(i+1) + bbase;
  296.         cA[i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
  297.         if(i) cA[i].fillcolor_lo = cA[i-1].fillcolor_hi;
  298.         else  cA[i].fillcolor_lo = backgroundColor;
  299.         cA[i].linecolor = contourLineColor;
  300.     }
  301.     }
  302.     return self;
  303. }
  304.  
  305.  
  306.  
  307. // ==========================================================================
  308. // Set surface data for regular Cartesian grid with optional interpolation.
  309. // This is probably the most important method of this object.
  310. // Use this if you want to plot data on uniform XY grid.
  311. // Input data f[] of size (nx,ny) will be interpolated into (n1x, n1y) for
  312. //     smoother contour plot if desired.  If no interpolation is needed,
  313. //     make (n1x, n1y) same as (nx, ny).  This will turn off interpolation.
  314. // 
  315. //     (int)  nx, ny            : index sizes of input data
  316. //     (int) n1x, n1y        : size of interpolated grid
  317. //     for (i = 0; i < nx*ny; i++) (float)f[i]   : f value at (x,y)
  318. //     f[i] is a 1-D array.  f(x,y) is in f[x + nx*y]
  319. //
  320. //     (xmin,ymin) and (xmax,ymax) define the (X,Y) domain of the plot.
  321. //
  322. // We extend the final (interpolated) data 3 points on each border to
  323. // let the contour algorithm close all contours.  1-point wide region copy
  324. // the data of the original border, and 2 points margins outside that are
  325. // set to the base value.  This way, open contour lines that would have
  326. // terminated at a border will get closed outside the original domain.
  327. // This extra region for contour closure outside the original domain is
  328. // outside the clip path, thus will not show up in the final plot.
  329. // ---------------------------------------------------------------------------
  330. //
  331. - setCartesianGridData:(float *)f :(float)xmin :(float)xmax
  332.                   :(float)ymin :(float)ymax
  333.                   ofSize:(int)nnx :(int)nny
  334.                   withInterpolationTo:(int)n1x :(int)n1y
  335. {
  336. int i, ix, iy;
  337. float xstep, ystep;
  338. float *xt, *yt;    /* temporary input arrays for interpolation */
  339. float *fd1;
  340. float **f2, **fd2a;
  341.  
  342.     xdmin = xmin;
  343.     xdmax = xmax;
  344.     ydmin = ymin;
  345.     ydmax = ymax;
  346.     // do reallocation only if # of data points has changed.
  347.     nx = n1x+N6;
  348.     ny = n1y+N6;
  349.     nx1 = n1x;
  350.     ny1 = n1y;
  351.     if(ndata != ((n1x+N6)*(n1y+N6)))
  352.     {
  353.         ndata = (n1x+N6)*(n1y+N6);        /* new grid size */
  354.         if(xd) free(xd);
  355.     if(yd) free(yd);
  356.     if(fd) free(fd);
  357.         xd = (float *)malloc((size_t)ndata*sizeof(float));
  358.         yd = (float *)malloc((size_t)ndata*sizeof(float));
  359.         fd = (float *)malloc((size_t)ndata*sizeof(float));
  360.     }
  361.  
  362.     if(nnx == n1x && nny == n1y)        /* Same grid size indicates NO INTERPOLATION */
  363.     {
  364.     /* copy original domain to expanded domain while defining border areas */
  365.     for(i=0; i<ndata; i++) {
  366.         ix = i % nx;
  367.         iy = i / nx;
  368.         if(ix >= N3 && ix < (nx-N3) && iy >= N3 && iy < (ny-N3))
  369.         fd[i] = f[(iy-N3)*nnx+(ix-N3)];        /* original data */
  370.         else if(iy < N3 && ix >= N3 && ix <(nx-N3))
  371.         fd[i] = f[ix-N3];
  372.         else if(iy >= (ny-N3) && ix >= N3 && ix <(nx-N3))
  373.         fd[i] = f[(nny-1)*nnx+(ix-N3)];
  374.         else if(ix < N3 && iy >= N3 && iy <(ny-N3))
  375.         fd[i] = f[(iy-N3)*nnx];
  376.         else if(ix >=(nx-N3) && iy >= N3 && iy <(ny-N3))
  377.         fd[i] = f[(iy-N3)*nnx + nnx-1];
  378.         else if(ix < N3 && iy <N3)        /* lower left */
  379.         fd[i] = f[0];
  380.         else if(ix < N3 && iy >= (ny-N3))    /* uppper left */
  381.         fd[i] = f[(nny-1)*nnx];
  382.         else if(ix >= (nx-N3) && iy < N3)    /* lower right */
  383.         fd[i] = f[nnx-1];
  384.         else if(ix >= (nx-N3) && iy >= (ny-N3))
  385.         fd[i] = f[(nny-1)*nnx + nnx-1];
  386.  
  387.         /* wipe 2 pixel border to base level */
  388.         if(ix < N2 || ix >= (nx-N2) || iy <N2 || iy >= (ny-N2))
  389.         fd[i] = basevalue;        /* totally out */
  390.     }
  391.     }
  392.     else        /* Interpolate (nx, ny) to new size grid (n1x, n1y) */
  393.     {
  394.     // Allocate temporary input arrays for interpolation
  395.     xt = (float *)malloc((size_t)nnx*sizeof(float));
  396.     yt = (float *)malloc((size_t)nny*sizeof(float));
  397.     fd1 =  (float *)malloc((size_t)nx1*ny1*sizeof(float)); /* for interpolation */
  398.     f2 = fmatrix(1, nnx, 1, nny);
  399.     fd2a = fmatrix(1, nnx, 1, nny);
  400.     xstep = (xdmax-xdmin)/(float)(nnx-1);
  401.     ystep = (ydmax-ydmin)/(float)(nny-1);
  402.     for(ix=0; ix<nnx; ix++)
  403.         xt[ix] = xstep*ix + xdmin;
  404.     for(iy=0; iy<nny; iy++)
  405.         yt[iy] = ystep*iy + ydmin;
  406.     /* copy original data to 2-D array for interpolation */
  407.     for(ix=0; ix <nnx; ix++)
  408.         for(iy=0; iy <nny; iy++)
  409.         f2[ix+1][iy+1] = f[iy*nnx+ix];
  410.  
  411.     /* precompute second-derivative arrays for splin2()
  412.      * 2-nd derivative is returned in fd2a[][].
  413.      * (xt-1) and (yt-1) are for passing C's zero-offset array to
  414.      * splie2() which expects 1-offset arrays. */
  415.     splie2(xt-1, yt-1, f2, nnx, nny, fd2a);
  416.  
  417.     // Interpolate f(nnx,nny) to finer array fd1(n1x,n1y)
  418.         xstep = (xdmax-xdmin)/(float)(nx1-1);
  419.     ystep = (ydmax-ydmin)/(float)(ny1-1);
  420.     for(ix=0; ix <nx1; ix++)
  421.         for(iy=0; iy <ny1; iy++)
  422.         splin2(xt-1, yt-1, f2, fd2a, nnx, nny,
  423.             xstep*ix+xdmin, ystep*iy+ydmin, &fd1[iy*nx1+ix]);
  424.  
  425.     /* Copy interpolated array to expanded array while seting border values. */
  426.     for(i=0; i<ndata; i++) {
  427.         ix = i % nx;
  428.         iy = i / nx;
  429.         if(ix >= N3 && ix < (nx-N3) && iy >= N3 && iy < (ny-N3))
  430.         fd[i] = fd1[(iy-N3)*n1x+(ix-N3)];        /* copy interpolated data */
  431.         else if(iy < N3 && ix >= N3 && ix <(nx-N3))        /* bottom */
  432.         fd[i] = fd1[ix-N3];
  433.         else if(iy >= (ny-N3) && ix >= N3 && ix <(nx-N3))      /* top */
  434.         fd[i] = fd1[(n1y-1)*n1x+(ix-N3)];
  435.         else if(ix < N3 && iy >= N3 && iy <(ny-N3))        /* left */
  436.         fd[i] = fd1[(iy-N3)*n1x];
  437.         else if(ix >=(nx-N3) && iy >= N3 && iy <(ny-N3))    /* right */
  438.         fd[i] = fd1[(iy-N3)*n1x + n1x-1];
  439.         else if(ix < N3 && iy <N3)        /* lower left */
  440.         fd[i] = fd1[0];
  441.         else if(ix < N3 && iy >= (ny-N3))    /* uppper left */
  442.         fd[i] = fd1[(n1y-1)*n1x];
  443.         else if(ix >= (nx-N3) && iy < N3)    /* lower right */
  444.         fd[i] = fd1[n1x-1];
  445.         else if(ix >= (nx-N3) && iy >= (ny-N3))
  446.         fd[i] = fd1[(n1y-1)*n1x + n1x-1];
  447.  
  448.         /* wipe 2 pixel border to base level */
  449.         if(ix < N2 || ix >= (nx-N2) || iy <N2 || iy >= (ny-N2))
  450.         fd[i] = basevalue;        /* totally out */
  451.     }
  452.  
  453.         // Free temporary input array.
  454.     if(f2) free_fmatrix(f2, 1, nnx, 1, nny);
  455.     if(fd2a) free_fmatrix(fd2a, 1, nnx, 1, nny);
  456.     if(fd1) free((void *)fd1);
  457.         if(xt) free((void *)xt);
  458.         if(yt) free((void *)yt);
  459.     } /* END OF: if(nnx == n1x && nny == n1y) {} else {} */
  460.  
  461.     // Generate proper xd[], yd[] array for countour_().
  462.     xstep = (xdmax-xdmin)/(float)(nx1-1);
  463.     ystep = (ydmax-ydmin)/(float)(ny1-1);
  464.     for(iy=0; iy<ny; iy++)
  465.     for(ix=0; ix<nx; ix++)
  466.     {
  467.         xd[ix+iy*nx] = xstep*(ix-N3) + xdmin;
  468.         yd[ix+iy*nx] = ystep*(iy-N3) + ydmin;
  469.     }
  470.  
  471.     // Do appropriate scaling
  472.     [self setScaleXY:xdmin :xdmax :ydmin :ydmax];
  473.  
  474.     // Get min and max fd[]
  475.     [self findMinMax:fd :ndata :&fdmin :&fdmax];
  476.     [self setDefaultContourAttributes];
  477.  
  478.     return self;
  479. }
  480.  
  481.  
  482.  
  483. // Set XY grid points and values at each grid points.
  484. // Grid does not have to be Cartesian or regular.
  485. // Use this routine if you need flexibility.  No interpolation of
  486. // grid is provided.  If you want interpolation,
  487. // you have to do that yourself external to this object.
  488. //     (int)  nj, nk                       : index sizes of data
  489. //     for (i = 0; i < nj*nk; i++) (float)x[i]   : x coordinate
  490. //     for (i = 0; i < nj*nk; i++) (float)y[i]   : y coordinate
  491. //     for (i = 0; i < nj*nk; i++) (float)f[i]   : f value at (x[i], y[i])}
  492. //     f[i] is a 1-D array.  f(j,k) is in f[j + nj*k].  Same for x[], y[].
  493. //
  494. // NOTE: This methond does not support contour closure by expanding the domain.
  495. //
  496. - setGridAndValueData:(float *)x :(float *)y :(float *)f ofSize:(int)nj :(int)nk
  497. {
  498. int i;
  499.     // do reallocation only if # of data points has changed.
  500.     nx = nj;
  501.     ny = nk;
  502.     if(ndata != (nj*nk))
  503.     {
  504.         ndata = nj*nk;
  505.         if(xd) free(xd);
  506.     if(yd) free(yd);
  507.     if(fd) free(fd);
  508.         xd = (float *)malloc((size_t)ndata*sizeof(float));
  509.         yd = (float *)malloc((size_t)ndata*sizeof(float));
  510.         fd = (float *)malloc((size_t)ndata*sizeof(float));
  511.     }
  512.     for(i=0; i<ndata; i++)    // copy data into instance variables
  513.     {
  514.     xd[i] = x[i];
  515.     yd[i] = y[i];
  516.     fd[i] = f[i];
  517.     }
  518.  
  519.     // Do appropriate scaling
  520.     // Get min and max of xd[], yd[], and fd[]
  521.     [self findMinMax:xd :ndata :&xdmin :&xdmax];
  522.     [self findMinMax:yd :ndata :&ydmin :&ydmax];
  523.     [self findMinMax:fd :ndata :&fdmin :&fdmax];
  524.     [self setScaleXY:xdmin :xdmax :ydmin :ydmax];
  525.     [self setDefaultContourAttributes];
  526.  
  527.     return self;
  528. }
  529.  
  530.  
  531.  
  532. // Controls width and whether a rectangular frame is drawn.
  533. //
  534. - setFrameEnable:(BOOL)fflag lineWidth:(float)fw
  535. {
  536.  
  537.     frameON= fflag;    // Draw bounds frame
  538.     frameLineWidth=fw;
  539.     return self;
  540. }
  541.  
  542. - setContourLineWidthMinor:(float)clw andMajor:(float)clwm
  543. {
  544.     minorContourLineWidth = clw;
  545.     majorContourLineWidth = clwm;
  546.     return self;
  547. }
  548.  
  549.  
  550.  
  551.  
  552. // Contours with number of points less than nmp are not plotted.
  553. // This is used to eliminate tiny "speckles" in contour plot.
  554. //
  555. - setMinNumberOfPointsPerContour:(int)mnp
  556. {
  557.     minNumPoints = mnp;
  558.     return self;
  559. }
  560.  
  561. - clear:sender
  562. {
  563.     ClearFlag = 1;
  564.     [self display];
  565.     ClearFlag = 0;
  566.     return self;
  567. }
  568.  
  569. - copy:sender
  570. {
  571.     [self copyPScode:sender];
  572.     return self;
  573. }
  574.  
  575. // Copies the current view into the Pasteboard as PostScript.
  576. //
  577. - copyPScode:sender
  578. {
  579.     NXStream *psStream;
  580.     id        pb;
  581.     char     *data;
  582.     int       dataLen, maxDataLen;
  583.  
  584.   /* Open a stream on memory where we will collect the PostScript */
  585.     psStream = NXOpenMemory(NULL, 0, NX_WRITEONLY);
  586.     if (!psStream)
  587.       return self;
  588.  
  589.   /* Tell the Pasteboard we're going to copy PostScript */
  590.     pb = [Pasteboard new];
  591.     [pb declareTypes:&NXPostScriptPboardType num:1 owner:self];
  592.  
  593.   /* writes the PostScript for the whole plot as EPS into the stream */
  594.     [self copyPSCodeInside:&bounds to:psStream];
  595.  
  596.   /* get the buffered up PostScript out of the stream */
  597.     NXGetMemoryBuffer(psStream, &data, &dataLen, &maxDataLen);
  598.  
  599.   /* put the buffer in the Pasteboard, free the stream (and the buffer) */
  600.     [pb writeType:NXPostScriptPboardType data:data length:dataLen];
  601.     NXCloseMemory(psStream, NX_FREEBUFFER);
  602.     return self;
  603. }
  604.  
  605. // Call back method messaged from computeContour() function.
  606. - accumContour:(int)icont :(int)np :(float *)x :(float *)y
  607. {
  608. ContourPath *newcntr;    /* scratch pad pointer to new contourList items */
  609. int c_closed = 1;
  610.  
  611.     if(np < minNumPoints) return self;        /* too few points for contour */
  612.  
  613.     if (fabs(x[0] - x[np-1]) < 0.0001 && fabs(y[0] - y[np-1]) < 0.0001) {
  614.     c_closed = 1;
  615.     x[np-1] = x[0];  y[np-1] = y[0];    /* guarantee closure */
  616.     }
  617.     else
  618.     c_closed = 0;            /* contour not closed */
  619.  
  620.     /* ---- Save new contour in a linked list ----------------------------------- */
  621.     /* Allocate new item for contourList */
  622.     newcntr = (ContourPath *)malloc((size_t)sizeof(ContourPath));
  623.     newcntr->num_pts = np;        /* # of (x,y) points in contour */
  624.     newcntr->closed = c_closed;        /* contour is closed flag */
  625.     newcntr->levelindex = icont;    /* levelindex-th contour (indx to contour level) */
  626.     newcntr->level = cA[icont].level;    /* value of contour */
  627.     /* Allocate arrays for X and Y coordinates.
  628.     * 7 more points may be needed to close unclosed contours with a path
  629.     * that goes around the domain.
  630.     */
  631.     newcntr->x = (float *)malloc((size_t)(np*sizeof(float)));
  632.     newcntr->y = (float *)malloc((size_t)(np*sizeof(float)));
  633.     /* copy (x,y) coordinates of contour path */
  634.     memcpy((void *)newcntr->x, x, np * sizeof(float));
  635.     memcpy((void *)newcntr->y, y, np * sizeof(float));
  636.     /* insert new item into the list */
  637.     newcntr->next = contourList;    /* point to previous one as next */
  638.     contourList = newcntr;        /* for next time */
  639.     numContours++;            /* increment number of contours counter */
  640.     /* ---- New contour saved in a linked list ----------------------------------- */
  641.     return self;
  642. }
  643.  
  644.  
  645. // This will override the default -drawSelf::, which does no real drawing.
  646. // Should never be called directly.  It will be called by [self display];
  647. //
  648. - drawSelf:(const NXRect *)r:(int)count
  649. {
  650. int j;
  651. ContourPath *cntr;
  652. DPSContext curContext = DPSGetCurrentContext();
  653.  
  654.     if(NXDrawingStatus != NX_DRAWING)
  655.         DPSPrintf(curContext, "\n%% Start of ContourView PostScript drawing..\n");
  656.     NXSetColor(backgroundColor);
  657.     NXRectFill(&bounds);
  658.  
  659.     if(!fd || ClearFlag) return self;    // if no data set, do nothing.
  660.  
  661.     if(contourList) [self freeContourList];        /* to start out */
  662.     numContours = 0;
  663.  
  664.  
  665.     computeContour(self, nx, ny, xd, yd, fd, nclevels, cA);
  666.  
  667.     if(debug)
  668.     fprintf(stderr, "ContourView: # levels=%d,  # contours=%d\n", nclevels, numContours);
  669.  
  670.     if(doFill) {
  671.     /* Sort the order of contour drawing from outermost to innermost for filling */
  672.         SortedCntrPtr = (ContourPath **)(malloc((size_t)(sizeof(ContourPath *)*numContours)));
  673.         sort_contourList(contourList, xdmin, xdmax, ydmin, ydmax, numContours, SortedCntrPtr);
  674.     for(j=0; j<numContours; j++) {
  675.         cntr = SortedCntrPtr[j];
  676.         [self findInsideHighLow:cntr];        /* inside contour goin up or down */
  677.         [self plotContour:cntr :curContext];    /* This plots one contour curve */
  678.     }
  679.     }
  680.     else {
  681.     /* Do NOT do filling -- just go through the list */
  682.     cntr = contourList;
  683.     while(cntr) {
  684.         [self plotContour:cntr :curContext];    /* This plots one contour curve */
  685.         cntr = cntr->next;        /* point to the next one */
  686.     };
  687.     }
  688.  
  689.     if(frameON) {
  690.         if(NXDrawingStatus != NX_DRAWING)
  691.             DPSPrintf(curContext, "\n%% Frame rectangle\n");
  692.     // Frame rectangle
  693.     NXSetColor(frameColor);
  694.     PSsetdash(psolid, 0, 0.0);
  695.     PSsetlinewidth(frameLineWidth);
  696.     PSrectstroke(pXmin,pYmin,pXmax,pYmax);    /* frame rectangle */
  697.     }
  698.  
  699.     if(doFill) free(SortedCntrPtr);
  700.     [self freeContourList];
  701.  
  702.     if(NXDrawingStatus != NX_DRAWING)
  703.         DPSPrintf(curContext, "\n%% End of CurveView drawing.\n\n");
  704.  
  705.     return self;
  706. }
  707.  
  708.  
  709.  
  710. - plotContour:(ContourPath *)cntr :(DPSContext)cContext
  711. {
  712. int i, np;
  713. float xp, yp;
  714. static char *ocmsg[] = {"open", "closed",""};
  715.  
  716.     np = cntr->num_pts;        /* # points in this contour */
  717.     if(NXDrawingStatus != NX_DRAWING)
  718.         DPSPrintf(cContext, "\n%% ### Contour: level index=%d, value=%g, #pts=%d, [%s]\n",
  719.             cntr->levelindex, cntr->level, np, ocmsg[(cntr->closed)?1:0]);
  720.  
  721.     for(i=0; i<np; i++) {
  722.     xp = (cntr->x[i] - xdmin)*ppxu;
  723.     yp = (cntr->y[i] - ydmin)*ppyu;
  724.     if(i==0)
  725.         PSmoveto(xp, yp);
  726.     else
  727.         PSlineto(xp, yp);
  728.     }
  729.     if(cntr->closed)
  730.         PSclosepath();
  731.  
  732.     if(doFill) {
  733.     if(cntr->hi_inside)
  734.         NXSetColor(cA[cntr->levelindex].fillcolor_hi);
  735.     else
  736.         NXSetColor(cA[cntr->levelindex].fillcolor_lo);
  737.     PSgsave();
  738.     PSfill();
  739.     PSgrestore();
  740.     }
  741.  
  742.     if(cA[cntr->levelindex].dash)
  743.     PSsetdash(pdash, 2, 0.0);
  744.     /* fprintf(fpp, "[4 4] 0 setdash\n"); */ /* dashed curve */
  745.     else
  746.     PSsetdash(psolid, 0, 0.0);
  747.     /* fprintf(fpp, "[] 0 setdash\n");    */    /* solid curve */
  748.     NXSetColor(cA[cntr->levelindex].linecolor);
  749.     PSstroke();
  750.     return self;
  751. }
  752.  
  753. - findInsideHighLow:(ContourPath *)cntr
  754. {
  755. int i, i3, j, ii=0, ix=0, iy=0, npp, inside=0, gpfound=0, jx=0, jy=0;
  756. float xstep, ystep, xg=0.0, yg=0.0, fg= 0.0;
  757. static int ipx[] = { 0, -1,  1,  0, -1,  1,  0, -1,  1 };
  758. static int ipy[] = { 0,  0,  0, -1, -1, -1,  1,  1,  1 };
  759.  
  760.     xstep = (xdmax-xdmin)/(float)(nx1-1);
  761.     ystep = (ydmax-ydmin)/(float)(ny1-1);
  762.     npp = cntr->num_pts;
  763.  
  764.     // Typical setting for hi_inside flag
  765.     if(cntr->level > 0.0) cntr->hi_inside = 1;
  766.     else          cntr->hi_inside = 0;
  767.  
  768.     // Now try to determine it exactly.  Start with i=3 in an attempt to get
  769.     // a good point that is not on the edge.
  770.     for(i3=3; i3<(npp+3); i3++) {
  771.     i = i3 % npp;        /* i = [3, .. (npp-1), 0, 1, 2] */
  772.     if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
  773.         continue;        // this contour point is outside domain
  774.     ix  = (cntr->x[i] - xdmin)/xstep + 0.5;    // indices of closest grid point
  775.     iy  = (cntr->y[i] - ydmin)/ystep + 0.5;
  776.     // find indices (jx, jy) which is inside contour and inside domain
  777.     // by checking 3x3 points centered on (ix, iy)
  778.     for(j=0; j<9; j++) {
  779.         jx = ix + ipx[j];    
  780.         jy = iy + ipy[j];
  781.         xg = xstep*(float)jx + xdmin;
  782.         yg = ystep*(float)jy + ydmin;
  783.         if([self pointInDomain:xg :yg] && non_zero_winding(xg, yg, cntr->x, cntr->y, npp)==2)
  784.         {
  785.         inside = 1;
  786.         break;        // out of for(j..)
  787.         }
  788.     }
  789.     if(!inside) continue;    // try another point on contour
  790.  
  791.     ii = i;
  792.     gpfound = 1;
  793.     break;                // good point found.
  794.     }  /* end of for(i=0...) */
  795.  
  796.     // If inside point not found, use closed point that is outside
  797.     if(!inside) {
  798.     for(i=0; i<npp; i++) {
  799.         if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
  800.             continue;        // this contour point is outside domain
  801.         jx  = (cntr->x[i] - xdmin)/xstep + 0.5;    // indices of closest grid point
  802.         jy  = (cntr->y[i] - ydmin)/ystep + 0.5;
  803.         xg = xstep*(float)jx + xdmin;
  804.         yg = ystep*(float)jy + ydmin;
  805.         if( ![self pointInDomain:xg :yg])    // at least it must be inside domain
  806.         continue;
  807.  
  808.         ii = i;
  809.         gpfound = 1;
  810.         if(non_zero_winding(xg, yg, cntr->x, cntr->y, npp) == 2) inside = 1;
  811.         break;
  812.     }
  813.     }
  814.  
  815.     fg = fd[(jy+N3)*nx+jx+N3];        // grid value
  816.     if(inside) {
  817.     if(fg > cntr->level) cntr->hi_inside = 1;    // inside is high
  818.     else             cntr->hi_inside = 0;
  819.     } else {
  820.     if(fg > cntr->level) cntr->hi_inside = 0;    // outside is high
  821.     else             cntr->hi_inside = 1;
  822.     }
  823.  
  824.     if(debug && !gpfound)
  825.     fprintf(stderr,
  826.         "ContourView: Can't determine hi_inside: idx=%d, val=%g, #pt=%d, hi_ins=%d\n",
  827.             cntr->levelindex, cntr->level, npp, cntr->hi_inside);
  828.     return self;
  829. }
  830.  
  831. - (BOOL)pointInDomain:(float)xx :(float)yy
  832. {
  833.     if(xx >= xdmin && xx <= xdmax && yy >= ydmin && yy <= ydmax)
  834.     return YES;
  835.     else
  836.     return NO;
  837. }
  838.  
  839.  
  840. - freeContourList;
  841. {
  842. ContourPath *cntr;
  843.     while((cntr = contourList))
  844.     {
  845.         if(cntr->x) free(cntr->x);
  846.         if(cntr->y) free(cntr->y);
  847.         contourList = cntr->next;
  848.         free(cntr);
  849.     };
  850.         contourList = NULL;
  851.     return self;
  852. }
  853.  
  854.  
  855. @end
  856.  
  857.  
  858.