This is computeContour.m in view mode; [Download] [Up]
/* computeContour.m
Just a collection of C functions, but made into a .m file to handle things
like "id"s and a little messaging. This is not an object implementation.
Receives grid data and contour informations and calls back CurveView object
for each contour it finds. A CurveView object passes its "id" so functions
here can do call backs.
92-11-30, Version 1.0, Izumi Ohzawa, izumi@pinoko.berkeley.edu.
Modified for ContourView. Added Color fills that work more-or-less OK for me.
------------------------------------------------------------------------------
Based on contour_plot.c from NeXTcontour1.4 by Thomas H. Pulliam,
pulliam@rft29.nas.nasa.gov
MS 202A-1 NASA Ames Research Center, Moffett Field, CA 94035.
This f2c converted file (contour_plot.c) is the only module taken from
NeXTcontour1.4. All other objects including ContourView[.h.m] class
has been written from scrach. Although ContourView is also in
NeXTcontour1.4, his and mine are different objects.
NeXTcontour_.__(.tar.Z) is available from FTP archive sites for NeXT apps.
I don't know how the original Fortran code looked like or where it came from,
other than that NeXTcontour1.4 is based on Pieter Bunings' PLOT3D package
for Computational Fluid Dynamics.
(int) nx,ny : index sizes of data
for (i = 0; i < nx*ny; i++) (float)x[i] : x coordinate
for (i = 0; i < nx*ny; i++) (float)y[i] : y coordinate
for (i = 0; i < nx*ny; i++) (float)f[i] : f[iy*nx+ix] value at (ix, iy)
*/
#import <objc/objc.h> /* for BOOL values YES, NO */
#import "ContourView.h"
#import "contour.h"
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
/* Common Block Declarations */
struct {
float xt[1000], yt[1000], zt[1000];
int ia[3000];
} pathbuf;
/* =================================================================================== */
/*
id vObj; Object that called this function (for call back).
int nx, ny; X, Y dimension of data.
float *x, *y, *f; arrays containing grid data [0..nx*ny -1]
int ncont; number of contour levels
CntrAttribute *cAttr; array containing attributes of contour levels [0..ncont-1]
*/
/* =================================================================================== */
int computeContour(id vObj, int nx, int ny, float *x, float *y, float *f,
int ncont, CntrAttribute *cAttr)
{
int x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset, ntmp_2,
ntmp_3, ntmp_4;
int ibeg, jbeg, idir;
float cont;
int i, j, iaend, icont;
int ij, np=0;
float wx=0.0, wy=0.0;
int ignext, lnstrt, iae, iia, iee, jee, imb, ima;
float fij;
int iss, jss;
float xyf, wxx=0.0, wyy=0.0;
if(nx == 1 || ny == 1)
return 0;
/* Don't ask me... This is some FORTRAN to C conversion shit having to do
* with addressing 0-offset C arrays as and 1-offset arrays.
* It can be cleaned up, but not now.
*/
x_dim1 = y_dim1 = f_dim1 = nx;
x_offset = nx + 1;
x -= x_offset;
y_offset = nx + 1;
y -= y_offset;
f_offset = nx + 1;
f -= f_offset;
/* Zero the marker array. */
for(iia = 0; iia < 3000; ++iia)
pathbuf.ia[iia] = 0;
/* lnstrt=1 (linestart) means we're starting a new line. */
lnstrt = 1;
ignext = 0;
/* ####### Loop through each contour level. ################################### */
for (icont = 0; icont < ncont; ++icont)
{
cont = cAttr[icont].level;
iss = 1;
iee = nx;
jss = 1;
jee = ny;
/* Flag points in IA where the the function increases through the contour
* level, not including the boundaries. This is so we have a list of at least
* one point on each contour line that doesn't intersect a boundary.
*/
L110:
iae = 0;
ntmp_2 = jee - 1;
for (j = jss + 1; j <= ntmp_2; ++j) {
imb = 0;
iaend = iae;
ntmp_3 = iee;
for (i = iss; i <= ntmp_3; ++i) {
if (f[i + j * f_dim1] <= cont) {
imb = 1;
} else if (imb == 1) {
++iae;
pathbuf.ia[iae - 1] = i * 1000 + j;
imb = 0;
/* Check if the IA array is full. If so, the subdividing algorithm goes like
* this: if we've marked at least one J row, drop back to the last completed
* J and call that the region. If we haven't even finished one J row, our
* region just extends to this I location.
*/
if (iae == 3000) {
if (j > jss + 1) {
iae = iaend;
jee = j;
} else {
/* Computing MIN */
ntmp_4 = j + 1;
jee = min(ntmp_4,jee);
iee = i;
}
goto L210;
}
}
/* L200: */
} /* for(i=iss; i <= ntmp_3; ++i) */
} /* for(j=jss+1; ... ) */
/* Search along the boundaries for contour line starts. IMA tells which
* boundary of the region we're on.
*/
L210:
ima = 1;
imb = 0;
ibeg = iss - 1;
jbeg = jss;
L1:
++ibeg;
if (ibeg == iee) {
ima = 2;
}
goto L5;
L2:
++jbeg;
if (jbeg == jee) {
ima = 3;
}
goto L5;
L3:
--ibeg;
if (ibeg == iss) {
ima = 4;
}
goto L5;
L4:
--jbeg;
if (jbeg == jss) {
ima = 5;
}
L5:
if (f[ibeg + jbeg * f_dim1] > cont) {
goto L7;
}
imb = 1;
L6:
switch ((int)ima) {
case 1: goto L1;
case 2: goto L2;
case 3: goto L3;
case 4: goto L4;
case 5: goto L91;
}
L7:
if (imb != 1) {
goto L6;
}
/* Got a start point. */
imb = 0;
i = ibeg; /* x index of starting point grid */
j = jbeg; /* y index of starting point grid */
fij = f[ibeg + jbeg * f_dim1]; /* z value of starting point */
/* Round the corner if necessary. */
switch ((int)ima) {
case 1: goto L21;
case 2: goto L11;
case 3: goto L12;
case 4: goto L13;
case 5: goto L51;
}
L11:
if (j != jss) {
goto L31;
}
goto L21;
L12:
if (i != iee) {
goto L41;
}
goto L31;
L13:
if (j != jee) {
goto L51;
}
goto L41;
/* This is how we start in the middle of the region, using IA. */
L10:
i = pathbuf.ia[iia - 1] / 1000;
j = pathbuf.ia[iia - 1] - i * 1000;
fij = f[i + j * f_dim1];
pathbuf.ia[iia - 1] = 0;
goto L21;
/* Look different directions to see which way the contour line went: */
/* 4 */
/* 1-|-3 */
/* 2 */
L20:
++j;
L21:
--i;
if (i < iss) {
goto L90;
}
idir = 1;
if (f[i + j * f_dim1] <= cont) {
goto L52;
}
fij = f[i + j * f_dim1];
goto L31;
L30:
--i;
L31:
--j;
if (j < jss) {
goto L90;
}
idir = 2;
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L41;
L40:
--j;
L41:
++i;
if (i > iee) {
goto L90;
}
idir = 3;
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L51;
L50:
++i;
L51:
++j;
idir = 4;
if (j > jee) {
goto L90;
}
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L21;
/* Wipe this point out of IA if it's in the list. */
L52:
if (iae == 0) {
goto L60;
}
ij = i * 1000 + j + 1000;
ntmp_3 = iae;
for (iia = 1; iia <= ntmp_3; ++iia) {
if (pathbuf.ia[iia - 1] == ij) {
pathbuf.ia[iia - 1] = 0;
goto L60;
}
/* L300: */
}
/* Do interpolation for X,Y coordinates. */
L60:
xyf = (cont - f[i + j * f_dim1]) / (fij - f[i + j * f_dim1]);
/* This tests for a contour point coinciding with a grid point. In this case
* the contour routine comes up with the same physical coordinate twice. If
* If we don't trap it, it can (in some cases significantly) increase the
* number of points in a contour line. Also, if this happens on the first
* point in a line, the second point could be misinterpreted as the end of a
* (circling) contour line.
*/
if(xyf == (float)0.0)
++ignext;
switch ((int)idir) {
case 1: /* east */
wxx = x[i + j * x_dim1] + xyf * (x[i + 1 + j * x_dim1] - x[i + j * x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + 1 + j * y_dim1] - y[i + j * y_dim1]);
break;
case 2: /* north */
wxx = x[i + j * x_dim1] + xyf * (x[i + (j + 1) * x_dim1] - x[i + j * x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + (j + 1) * y_dim1] - y[i + j * y_dim1]);
break;
case 3: /* west */
wxx = x[i + j * x_dim1] + xyf * (x[i - 1 + j * x_dim1] - x[i + j * x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i - 1 + j * y_dim1] - y[i + j * y_dim1]);
break;
case 4: /* south */
wxx = x[i + j * x_dim1] + xyf * (x[i + (j - 1) * x_dim1] - x[i + j * x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + (j - 1) * y_dim1] - y[i + j * y_dim1]);
break;
}
/* Figure out what to do with this point. */
if (lnstrt != 1) {
goto L66; /* if not the first point in a contour */
}
/* ###### This is the first point in a contour line. ###### */
np = 1; /* # of points in contour set to 1 */
pathbuf.xt[np - 1] = wxx;
pathbuf.yt[np - 1] = wyy;
wx = wxx; /* save starting point as (wx, wy) */
wy = wyy;
lnstrt = 0; /* clear first point flag, we got one */
goto L67;
/* Second point and after comes here */
/* Add a point to this line. Check for duplicate point first. */
L66:
if(ignext == 2) {
if(wxx == pathbuf.xt[np - 1] && wyy == pathbuf.yt[np - 1]) {
ignext = 0;
goto L67;
} else {
ignext = 1;
}
}
++np; /* increment # of points in contour */
pathbuf.xt[np - 1] = wxx;
pathbuf.yt[np - 1] = wyy;
/* See if the temporary array xt, yt are full. Sendf it out if so. */
if (np == 1000) {
[vObj accumContour:icont :np :pathbuf.xt :pathbuf.yt];
pathbuf.xt[0] = pathbuf.xt[np - 1]; /* last pt becomes first point to continue */
pathbuf.yt[0] = pathbuf.yt[np - 1];
np = 1;
}
/* Check to see if we're back to the initial point. */
if (wxx != wx) {
goto L67;
}
if (wyy == wy) {
goto L90;
}
/* Nope. Search for the next point on this line. */
L67:
switch ((int)idir) {
case 1: goto L50;
case 2: goto L20;
case 3: goto L30;
case 4: goto L40;
}
/* This is the end of a contour line. After this we'll start a new line.*/
L90:
lnstrt = 1; /* contour line start flag */
ignext = 0;
[vObj accumContour:icont :np :pathbuf.xt :pathbuf.yt];
/* If we're not done looking along the boundaries, go look there some more.*/
if (ima != 5) {
goto L6;
}
/* Otherwise, get the next start out of IA. */
L91:
if (iae != 0) {
ntmp_3 = iae;
for (iia = 1; iia <= ntmp_3; ++iia) {
if (pathbuf.ia[iia - 1] != 0) {
goto L10;
}
/* L400: */
}
}
/* And if there are no more of these, we're done with this region. If we've
* subdivided, update the region pointers and go back for more.
*/
if (iee == nx) {
if (jee != ny) {
jss = jee;
jee = ny;
goto L110;
}
} else {
iss = iee;
iee = nx;
goto L110;
}
/* L1000: */
/* ########### Loop back for the next contour level. ############################# */
} /* end for(icont=0; ...) */
return 0;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.