/* 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;
}




