home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1995 August / NEBULA.bin / SourceCode / Classes / ContourView / computeContour.m < prev    next >
Encoding:
Text File  |  1992-12-02  |  10.6 KB  |  466 lines
  1. /* computeContour.m
  2.  
  3.   Just a collection of C functions, but made into a .m file to handle things
  4.   like "id"s and a little messaging.  This is not an object implementation.
  5.  
  6.   Receives grid data and contour informations and calls back CurveView object
  7.   for each contour it finds.  A CurveView object passes its "id" so functions
  8.   here can do call backs.
  9.  
  10.   92-11-30, Version 1.0, Izumi Ohzawa, izumi@pinoko.berkeley.edu.
  11.   Modified for ContourView.  Added Color fills that work more-or-less OK for me.
  12.  
  13.   ------------------------------------------------------------------------------
  14.   Based on contour_plot.c from NeXTcontour1.4 by Thomas H. Pulliam,
  15.   pulliam@rft29.nas.nasa.gov
  16.   MS 202A-1 NASA Ames Research Center, Moffett Field, CA 94035.
  17.  
  18.   This f2c converted file (contour_plot.c) is the only module taken from
  19.   NeXTcontour1.4.  All other objects including ContourView[.h.m] class
  20.   has been written from scrach.  Although ContourView is also in
  21.   NeXTcontour1.4, his and mine are different objects.
  22.  
  23.   NeXTcontour_.__(.tar.Z) is available from FTP archive sites for NeXT apps.
  24.   I don't know how the original Fortran code looked like or where it came from,
  25.   other than that NeXTcontour1.4 is based on Pieter Bunings' PLOT3D package
  26.   for Computational Fluid Dynamics.
  27.  
  28.      (int)  nx,ny                       : index sizes of data
  29.      for (i = 0; i < nx*ny; i++) (float)x[i]   : x coordinate
  30.      for (i = 0; i < nx*ny; i++) (float)y[i]   : y coordinate
  31.      for (i = 0; i < nx*ny; i++) (float)f[i]   : f[iy*nx+ix] value at (ix, iy)
  32. */
  33.  
  34. #import <objc/objc.h>        /* for BOOL values YES, NO */
  35.  
  36. #import "ContourView.h"
  37. #import "contour.h"
  38.  
  39. #define min(a,b) ((a) <= (b) ? (a) : (b))
  40. #define max(a,b) ((a) >= (b) ? (a) : (b))
  41.  
  42.  
  43.  
  44.  
  45. /* Common Block Declarations */
  46.  
  47. struct {
  48.     float xt[1000], yt[1000], zt[1000];
  49.     int ia[3000];
  50. } pathbuf;
  51.  
  52.  
  53. /* =================================================================================== */
  54. /*
  55. id  vObj;        Object that called this function (for call back).
  56. int nx, ny;           X, Y dimension of data.
  57. float *x, *y, *f;     arrays containing grid data [0..nx*ny -1]
  58. int ncont;        number of contour levels
  59. CntrAttribute *cAttr;   array containing attributes of contour levels [0..ncont-1]
  60. */
  61. /* =================================================================================== */
  62.  
  63. int computeContour(id vObj, int nx, int ny, float *x, float *y, float *f,
  64.         int ncont, CntrAttribute *cAttr)
  65. {
  66. int x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset, ntmp_2, 
  67.     ntmp_3, ntmp_4;
  68.  
  69. int ibeg, jbeg, idir;
  70. float cont;
  71. int i, j, iaend, icont;
  72. int ij, np=0;
  73. float wx=0.0, wy=0.0;
  74. int ignext, lnstrt, iae, iia, iee, jee, imb, ima;
  75. float fij;
  76. int iss, jss;
  77. float xyf, wxx=0.0, wyy=0.0;
  78.  
  79.     if(nx == 1 || ny == 1)
  80.     return 0;
  81.  
  82.     /* Don't ask me...  This is some FORTRAN to C conversion shit having to do
  83.      * with addressing 0-offset C arrays as and 1-offset arrays.
  84.      * It can be cleaned up, but not now.
  85.      */
  86.     x_dim1 = y_dim1 = f_dim1 = nx;
  87.     x_offset = nx + 1;
  88.     x -= x_offset;
  89.     y_offset = nx + 1;
  90.     y -= y_offset;
  91.     f_offset = nx + 1;
  92.     f -= f_offset;
  93.  
  94.     /*   Zero the marker array. */
  95.     for(iia = 0; iia < 3000; ++iia)
  96.     pathbuf.ia[iia] = 0;
  97.  
  98.     /* lnstrt=1 (linestart) means we're starting a new line. */
  99.     lnstrt = 1;
  100.     ignext = 0;
  101.  
  102.     /* #######  Loop through each contour level. ################################### */
  103.     for (icont = 0; icont < ncont; ++icont)
  104.     {
  105.     cont = cAttr[icont].level;
  106.     iss = 1;
  107.     iee = nx;
  108.     jss = 1;
  109.     jee = ny;
  110.  
  111. /*  Flag points in IA where the the function increases through the contour
  112.  *  level, not including the boundaries.  This is so we have a list of at least
  113.  *   one point on each contour line that doesn't intersect a boundary.
  114. */
  115.  
  116. L110:
  117.     iae = 0;
  118.     ntmp_2 = jee - 1;
  119.     for (j = jss + 1; j <= ntmp_2; ++j) {
  120.         imb = 0;
  121.         iaend = iae;
  122.         ntmp_3 = iee;
  123.         for (i = iss; i <= ntmp_3; ++i) {
  124.         if (f[i + j * f_dim1] <= cont) {
  125.             imb = 1;
  126.         } else if (imb == 1) {
  127.             ++iae;
  128.             pathbuf.ia[iae - 1] = i * 1000 + j;
  129.             imb = 0;
  130.  
  131. /*  Check if the IA array is full.  If so, the subdividing algorithm goes like
  132.  *  this:  if we've marked at least one J row, drop back to the last completed
  133.  *  J and call that the region.  If we haven't even finished one J row, our
  134.  *   region just extends to this I location.
  135. */
  136.  
  137.             if (iae == 3000) {
  138.             if (j > jss + 1) {
  139.                 iae = iaend;
  140.                 jee = j;
  141.             } else {
  142. /* Computing MIN */
  143.                 ntmp_4 = j + 1;
  144.                 jee = min(ntmp_4,jee);
  145.                 iee = i;
  146.             }
  147.             goto L210;
  148.             }
  149.  
  150.         }
  151. /* L200: */
  152.         }  /* for(i=iss; i <= ntmp_3; ++i) */
  153.     }   /* for(j=jss+1; ... ) */
  154.  
  155. /*   Search along the boundaries for contour line starts.  IMA tells which
  156.  *   boundary of the region we're on.
  157. */
  158.  
  159. L210:
  160.     ima = 1;
  161.     imb = 0;
  162.     ibeg = iss - 1;
  163.     jbeg = jss;
  164.  
  165. L1:
  166.     ++ibeg;
  167.     if (ibeg == iee) {
  168.         ima = 2;
  169.     }
  170.     goto L5;
  171. L2:
  172.     ++jbeg;
  173.     if (jbeg == jee) {
  174.         ima = 3;
  175.     }
  176.     goto L5;
  177. L3:
  178.     --ibeg;
  179.     if (ibeg == iss) {
  180.         ima = 4;
  181.     }
  182.     goto L5;
  183. L4:
  184.     --jbeg;
  185.     if (jbeg == jss) {
  186.         ima = 5;
  187.     }
  188. L5:
  189.     if (f[ibeg + jbeg * f_dim1] > cont) {
  190.         goto L7;
  191.     }
  192.     imb = 1;
  193. L6:
  194.     switch ((int)ima) {
  195.         case 1:  goto L1;
  196.         case 2:  goto L2;
  197.         case 3:  goto L3;
  198.         case 4:  goto L4;
  199.         case 5:  goto L91;
  200.     }
  201. L7:
  202.     if (imb != 1) {
  203.         goto L6;
  204.     }
  205.  
  206. /*   Got a start point. */
  207.  
  208.     imb = 0;
  209.     i = ibeg;            /* x index of starting point grid */
  210.     j = jbeg;            /* y index of starting point grid */
  211.     fij = f[ibeg + jbeg * f_dim1];    /* z value of starting point */
  212.  
  213. /*   Round the corner if necessary. */
  214.  
  215.     switch ((int)ima) {
  216.         case 1:  goto L21;
  217.         case 2:  goto L11;
  218.         case 3:  goto L12;
  219.         case 4:  goto L13;
  220.         case 5:  goto L51;
  221.     }
  222. L11:
  223.     if (j != jss) {
  224.         goto L31;
  225.     }
  226.     goto L21;
  227. L12:
  228.     if (i != iee) {
  229.         goto L41;
  230.     }
  231.     goto L31;
  232. L13:
  233.     if (j != jee) {
  234.         goto L51;
  235.     }
  236.     goto L41;
  237.  
  238. /*   This is how we start in the middle of the region, using IA. */
  239.  
  240. L10:
  241.     i = pathbuf.ia[iia - 1] / 1000;
  242.     j = pathbuf.ia[iia - 1] - i * 1000;
  243.     fij = f[i + j * f_dim1];
  244.     pathbuf.ia[iia - 1] = 0;
  245.     goto L21;
  246. /*  Look different directions to see which way the contour line went: */
  247. /*   4   */
  248. /* 1-|-3 */
  249. /*   2   */
  250.  
  251. L20:
  252.     ++j;
  253. L21:
  254.     --i;
  255.     if (i < iss) {
  256.         goto L90;
  257.     }
  258.     idir = 1;
  259.     if (f[i + j * f_dim1] <= cont) {
  260.         goto L52;
  261.     }
  262.     fij = f[i + j * f_dim1];
  263.     goto L31;
  264. L30:
  265.     --i;
  266. L31:
  267.     --j;
  268.     if (j < jss) {
  269.         goto L90;
  270.     }
  271.     idir = 2;
  272.     if (f[i + j * f_dim1] <= cont) {
  273.         goto L60;
  274.     }
  275.     fij = f[i + j * f_dim1];
  276.     goto L41;
  277. L40:
  278.     --j;
  279. L41:
  280.     ++i;
  281.     if (i > iee) {
  282.         goto L90;
  283.     }
  284.     idir = 3;
  285.     if (f[i + j * f_dim1] <= cont) {
  286.         goto L60;
  287.     }
  288.     fij = f[i + j * f_dim1];
  289.     goto L51;
  290. L50:
  291.     ++i;
  292. L51:
  293.     ++j;
  294.     idir = 4;
  295.     if (j > jee) {
  296.         goto L90;
  297.     }
  298.     if (f[i + j * f_dim1] <= cont) {
  299.         goto L60;
  300.     }
  301.     fij = f[i + j * f_dim1];
  302.     goto L21;
  303.  
  304. /*   Wipe this point out of IA if it's in the list. */
  305.  
  306. L52:
  307.     if (iae == 0) {
  308.         goto L60;
  309.     }
  310.     ij = i * 1000 + j + 1000;
  311.     ntmp_3 = iae;
  312.     for (iia = 1; iia <= ntmp_3; ++iia) {
  313.         if (pathbuf.ia[iia - 1] == ij) {
  314.         pathbuf.ia[iia - 1] = 0;
  315.         goto L60;
  316.         }
  317. /* L300: */
  318.     }
  319.  
  320. /*   Do interpolation for X,Y coordinates. */
  321.  
  322. L60:
  323.     xyf = (cont - f[i + j * f_dim1]) / (fij - f[i + j * f_dim1]);
  324.  
  325. /*  This tests for a contour point coinciding with a grid point.  In this case
  326.  *  the contour routine comes up with the same physical coordinate twice.  If
  327.  *  If we don't trap it, it can (in some cases significantly) increase the
  328.  *  number of points in a contour line.  Also, if this happens on the first
  329.  *  point in a line, the second point could be misinterpreted as the end of a
  330.  *   (circling) contour line.
  331. */
  332.     if(xyf == (float)0.0)
  333.         ++ignext;
  334.  
  335.     switch ((int)idir) {
  336.         case 1:    /* east */
  337.         wxx = x[i + j * x_dim1] + xyf * (x[i + 1 + j * x_dim1] - x[i + j * x_dim1]);
  338.         wyy = y[i + j * y_dim1] + xyf * (y[i + 1 + j * y_dim1] - y[i + j * y_dim1]);
  339.         break;
  340.         case 2:    /* north */
  341.         wxx = x[i + j * x_dim1] + xyf * (x[i + (j + 1) * x_dim1] - x[i + j * x_dim1]);
  342.         wyy = y[i + j * y_dim1] + xyf * (y[i + (j + 1) * y_dim1] - y[i + j * y_dim1]);
  343.         break;
  344.         case 3:    /* west */
  345.         wxx = x[i + j * x_dim1] + xyf * (x[i - 1 + j * x_dim1] - x[i + j * x_dim1]);
  346.         wyy = y[i + j * y_dim1] + xyf * (y[i - 1 + j * y_dim1] - y[i + j * y_dim1]);
  347.         break;
  348.         case 4:    /* south */
  349.         wxx = x[i + j * x_dim1] + xyf * (x[i + (j - 1) * x_dim1] - x[i + j * x_dim1]);
  350.         wyy = y[i + j * y_dim1] + xyf * (y[i + (j - 1) * y_dim1] - y[i + j * y_dim1]);
  351.         break;
  352.     }
  353.  
  354.     /*   Figure out what to do with this point. */
  355.     if (lnstrt != 1) {
  356.         goto L66;        /* if not the first point in a contour */
  357.     }
  358.  
  359.         /* ######  This is the first point in a contour line. ###### */
  360.     np = 1;                /* # of points in contour set to 1 */
  361.     pathbuf.xt[np - 1] = wxx;
  362.     pathbuf.yt[np - 1] = wyy;
  363.     wx = wxx;            /* save starting point as (wx, wy) */
  364.     wy = wyy;
  365.     lnstrt = 0;            /* clear first point flag, we got one */
  366.     goto L67;
  367.  
  368. /*   Second point and after comes here */
  369. /*   Add a point to this line.  Check for duplicate point first. */
  370.  
  371. L66:
  372.     if(ignext == 2) {
  373.         if(wxx == pathbuf.xt[np - 1] && wyy == pathbuf.yt[np - 1]) {
  374.         ignext = 0;
  375.         goto L67;
  376.         } else {
  377.         ignext = 1;
  378.         }
  379.     }
  380.  
  381.     ++np;                /* increment # of points in contour */
  382.     pathbuf.xt[np - 1] = wxx;
  383.     pathbuf.yt[np - 1] = wyy;
  384.  
  385. /*   See if the temporary array xt, yt are full.  Sendf it out if so. */
  386.  
  387.     if (np == 1000) {
  388.         [vObj accumContour:icont :np :pathbuf.xt :pathbuf.yt];
  389.         pathbuf.xt[0] = pathbuf.xt[np - 1];    /* last pt becomes first point to continue */
  390.         pathbuf.yt[0] = pathbuf.yt[np - 1];
  391.         np = 1;
  392.     }
  393.  
  394. /*   Check to see if we're back to the initial point. */
  395.  
  396.     if (wxx != wx) {
  397.         goto L67;
  398.     }
  399.     if (wyy == wy) {
  400.         goto L90;
  401.     }
  402.  
  403. /*   Nope.  Search for the next point on this line. */
  404.  
  405. L67:
  406.     switch ((int)idir) {
  407.         case 1:  goto L50;
  408.         case 2:  goto L20;
  409.         case 3:  goto L30;
  410.         case 4:  goto L40;
  411.     }
  412.  
  413. /*  This is the end of a contour line.  After this we'll start a new line.*/
  414.  
  415. L90:
  416.     lnstrt = 1;        /* contour line start flag */
  417.     ignext = 0;
  418.     [vObj accumContour:icont :np :pathbuf.xt :pathbuf.yt];
  419.  
  420.  
  421. /*  If we're not done looking along the boundaries, go look there some more.*/
  422.  
  423.     if (ima != 5) {
  424.         goto L6;
  425.     }
  426.  
  427. /*   Otherwise, get the next start out of IA. */
  428.  
  429. L91:
  430.     if (iae != 0) {
  431.         ntmp_3 = iae;
  432.         for (iia = 1; iia <= ntmp_3; ++iia) {
  433.         if (pathbuf.ia[iia - 1] != 0) {
  434.             goto L10;
  435.         }
  436. /* L400: */
  437.         }
  438.     }
  439.  
  440. /*  And if there are no more of these, we're done with this region.  If we've
  441.  *   subdivided, update the region pointers and go back for more.
  442. */
  443.  
  444.     if (iee == nx) {
  445.         if (jee != ny) {
  446.         jss = jee;
  447.         jee = ny;
  448.         goto L110;
  449.         }
  450.     } else {
  451.         iss = iee;
  452.         iee = nx;
  453.         goto L110;
  454.     }
  455.  
  456. /* L1000: */
  457. /*  ########### Loop back for the next contour level. ############################# */
  458.     }   /* end for(icont=0; ...) */
  459.  
  460.     return 0;
  461. }
  462.  
  463.  
  464.  
  465.  
  466.