/*
* $Id: gcntr.c,v 1.1.1.1 2005/09/18 22:04:32 dhmunro Exp $
* General purpose contour tracer for quadrilateral meshes.
* Handles single level contours, or region between a pair of levels.
*/
/* Copyright (c) 2005, The Regents of the University of California.
* All rights reserved.
* This file is part of yorick (http://yorick.sourceforge.net).
* Read the accompanying LICENSE file for details.
*/
#include "gist.h"
/* note: The static functions in this file that actually do the work
* are independent of gist.h, except for the data type GpReal.
* However, they do make use of the region and triangulation arrays
* explained in gist.h in the definition of the GaQuadMesh struct. */
/* What is a contour?
*
* Given a quadrilateral mesh (x,y), and values of a z at the points
* of that mesh, we seek a set of polylines connecting points at a
* particular value of z. Each point on such a contour curve lies
* on an edge of the mesh, at a point linearly interpolated to the
* contour level z0 between the given values of z at the endpoints
* of the edge.
*
* Identifying these points is easy. Figuring out how to connect them
* into a curve -- or possibly a set of disjoint curves -- is difficult.
* Each disjoint curve may be either a closed circuit, or it may begin
* and end on a mesh boundary.
*
* One of the problems with a quadrilateral mesh is that when the z
* values at one pair of diagonally opposite points lie below z0, and
* the values at the other diagonal pair of the same zone lie above z0,
* all four edges of the zone are cut, and there is an ambiguity in
* how we should connect the points. I call this a saddle zone.
* The problem is that two disjoint curves cut through a saddle zone
* (I reject the alternative of connecting the opposite points to make
* a single self-intersecting curve, since those make ugly contour plots
* -- I've tried it). The real problem with saddle zones is that you
* need to communicate the connectivity decision you make back to the
* calling routine, since for the next contour level, we need to tell
* the contour tracer to make the same decision as on the previous
* level. The input/output triangulation array is the solution to this
* nasty problem.
*
* Another complicating factor is that there may be logical holes in
* the mesh -- zones which do not exist. We want our contours to stop
* if they hit the edge of such a zone, just as if they'd hit the edge
* of the whole mesh. The input region array addresses this issue.
*
* Yet another complication: We may want a list of closed polygons which
* outline the region between two contour levels z0 and z1. These may
* include sections of the mesh boundary (including edges of logical
* holes defined by the region array), in addition to sections of the
* contour curves at one or both levels. This introduces a huge
* topological problem -- if one of the closed contours (possibly
* including an interior logical hole in the mesh, but not any part of
* the boundary of the whole mesh) encloses a region which is not
* between z0 and z1, that curve must be connected by a slit (or "branch
* cut") to the enclosing curve, so that the list of disjoint polygons
* we return is each simply connected.
*
* Okay, one final stunning difficulty: For the two level case, no
* individual polygon should have more than a few thousand sides, since
* huge filled polygons place an inordinate load on rendering software,
* which needs an amount of scratch space proportional to the number
* of sides it needs to fill. So in the two level case, we want to
* chunk the mesh into rectangular pieces of no more than, say, 30x30
* zones, which keeps each returned polygon to less than a few thousand
* sides (the worst case is very very bad -- you can easily write down
* a function and two level values which produce a polygon that cuts
* every edge of the mesh twice).
*/
/*
* Here is the numbering scheme for points, edges, and zones in
* the mesh -- note that each ij corresponds to one point, one zone,
* one i-edge (i=constant edge) and one j-edge (j=constant edge):
*
* (ij-1)-------(ij)-------(ij)
* | |
* | |
* | |
* (ij-1) (ij) (ij)
* | |
* | |
* | |
* (ij-iX-1)----(ij-iX)----(ij-iX)
*
* At each point, the function value is either 0, 1, or 2, depending
* on whether it is below z0, between z0 and z1, or above z1.
* Each zone either exists (1) or not (0).
* From these three bits of data, all of the curve connectivity follows.
*
* The tracing algorithm is naturally edge-based: Either you are at a
* point where a level cuts an edge, ready to step across a zone to
* another edge, or you are drawing the edge itself, if it happens to
* be a boundary with at least one section between z0 and z1.
*
* In either case, the edge is a directed edge -- either the zone
* you are advancing into is to its left or right, or you are actually
* drawing it. I always trace curves keeping the region between z0 and
* z1 to the left of the curve. If I'm tracing a boundary, I'm always
* moving CCW (counter clockwise) around the zone that exists. And if
* I'm about to cross a zone, I'll make the direction of the edge I'm
* sitting on be such that the zone I'm crossing is to its left.
*
* I start tracing each curve near its lower left corner (mesh oriented
* as above), which is the first point I encounter scanning through the
* mesh in order. When I figure the 012 z values and zonal existence,
* I also mark the potential starting points: Each edge may harbor a
* potential starting point corresponding to either direction, so there
* are four start possibilities at each ij point. Only the following
* possibilities need to be marked as potential starting edges:
*
* +-+-+-+
* | | | |
* A-0-C-+ One or both levels cut E and have z=1 above them, and
* | EZ| | 0A is cut and either 0C is cut or CD is cut.
* +-B-D-+ Or, one or both levels cut E and E is a boundary edge.
* | | | | (and Z exists)
* +-+-+-+
*
* +-+-+-+
* | | | |
* +-A-0-C One or both levels cut E and have z=1 below them, and
* | |ZE | 0A is cut and either 0C is cut or CD is cut.
* +-+-B-D Or, one or both levels cut E and E is a boundary edge.
* | | | | (and Z exists)
* +-+-+-+
*
* +-+-+-+
* | | | |
* +-+-+-+ E is a boundary edge, Z exists, at some point on E
* | |Z| | lies between the levels.
* +-+E+-+
* | | | |
* +-+-+-+
*
* +-+-+-+
* | | | |
* +-+E+-+ E is a boundary edge, Z exists, at some point on E
* | |Z| | lies between the levels.
* +-+-+-+
* | | | |
* +-+-+-+
*
* During the first tracing pass, the start mark is erased whenever
* any non-starting edge is encountered, reducing the number of points
* that need to be considered for the second pass. The first pass
* makes the basic connectivity decisions. It figures out how many
* disjoint curves there will be, and identifies slits for the two level
* case or open contours for the single level case, and removes all but
* the actual start markers. A second tracing pass can perform the
* actual final trace.
*/
/* ------------------------------------------------------------------------ */
/* the data about edges, zones, and points -- boundary or not, exists
* or not, z value 0, 1, or 2 -- is kept in a mesh sized data array */
typedef short Cdata;
/* here is the minimum structure required to tell where we are in the
* mesh sized data array */
typedef struct Csite Csite;
struct Csite {
long edge; /* ij of current edge */
long left; /* +-1 or +-imax as the zone is to right, left, below,
* or above the edge */
long imax; /* imax for the mesh */
long jmax; /* jmax for the mesh */
long n; /* number of points marked on this curve so far */
long count; /* count of start markers visited */
GpReal zlevel[2]; /* contour levels, zlevel[1]<=zlevel[0]
* signals single level case */
short *triangle; /* triangulation array for the mesh */
int *reg; /* region array for the mesh */
long edge0, left0; /* starting site on this curve for closure */
int level0; /* starting level for closure */
long edge00; /* site needing START_ROW mark */
/* making the actual marks requires a bunch of other stuff */
const GpReal *x, *y, *z; /* mesh coordinates and function values */
GpReal *xcp, *ycp; /* output contour points */
};
/* the Cdata array consists of the following bits:
* Z_VALUE (2 bits) 0, 1, or 2 function value at point
* ZONE_EX 1 zone exists, 0 zone doesn't exist
* I_BNDY this i-edge (i=constant edge) is a mesh boundary
* J_BNDY this j-edge (i=constant edge) is a mesh boundary
* I0_START this i-edge is a start point into zone to left
* I1_START this i-edge is a start point into zone to right
* J0_START this j-edge is a start point into zone below
* J1_START this j-edge is a start point into zone above
* START_ROW next start point is in current row (accelerates 2nd pass)
* SLIT_UP marks this i-edge as the beginning of a slit upstroke
* SLIT_DN marks this i-edge as the beginning of a slit downstroke
* OPEN_END marks an i-edge start point whose other endpoint is
* on a boundary for the single level case
* ALL_DONE marks final start point
*/
#define Z_VALUE 0x0003
#define ZONE_EX 0x0004
#define I_BNDY 0x0008
#define J_BNDY 0x0010
#define I0_START 0x0020
#define I1_START 0x0040
#define J0_START 0x0080
#define J1_START 0x0100
#define START_ROW 0x0200
#define SLIT_UP 0x0400
#define SLIT_DN 0x0800
#define OPEN_END 0x1000
#define ALL_DONE 0x2000
/* some helpful macros to find points relative to a given directed
* edge -- points are designated 0, 1, 2, 3 CCW around zone with 0 and
* 1 the endpoints of the current edge */
#define FORWARD(left,ix) ((left)>0?((left)>1?1:-(ix)):((left)<-1?-1:(ix)))
#define POINT0(edge,fwd) ((edge)-((fwd)>0?fwd:0))
#define POINT1(edge,fwd) ((edge)+((fwd)<0?fwd:0))
#define IS_JEDGE(edge,left) ((left)>0?((left)>1?1:0):((left)<-1?1:0))
#define ANY_START (I0_START|I1_START|J0_START|J1_START)
#define START_MARK(left) \
((left)>0?((left)>1?J1_START:I1_START):((left)<-1?J0_START:I0_START))
/* ------------------------------------------------------------------------ */
/* these actually mark points */
static int zone_crosser(Csite *site, Cdata *data, int level, int pass2);
static int edge_walker(Csite *site, Cdata *data, int pass2);
static int slit_cutter(Csite *site, Cdata *data, int up, int pass2);
/* this calls the first three to trace the next disjoint curve
* -- return value is number of points on this curve, or
* 0 if there are no more curves this pass
* -(number of points) on first pass if:
* this is two level case, and the curve closed on a hole
* this is single level case, curve is open, and will start from
* a different point on the second pass
* -- in both cases, this curve will be combined with another
* on the second pass */
static long curve_tracer(Csite *site, Cdata *data, int pass2);
/* this initializes the data array for curve_tracer */
static void data_init(Csite *site, Cdata *data, int region, long nchunk);
/* ------------------------------------------------------------------------ */
/* zone_crosser assumes you are sitting at a cut edge about to cross
* the current zone. It always marks the initial point, crosses at
* least one zone, and marks the final point. On non-boundary i-edges,
* it is responsible for removing start markers on the first pass. */
static int zone_crosser(Csite *site, Cdata *data, int level, int pass2)
{
long edge= site->edge;
long left= site->left;
long n= site->n;
long fwd= FORWARD(left,site->imax);
long p0, p1;
int jedge= IS_JEDGE(edge,left);
long edge0= site->edge0;
long left0= site->left0;
int level0= site->level0==level;
int two_levels= site->zlevel[1]>site->zlevel[0];
short *triangle= site->triangle;
const GpReal *x= pass2? site->x : 0;
const GpReal *y= pass2? site->y : 0;
const GpReal *z= pass2? site->z : 0;
GpReal zlevel= pass2? site->zlevel[level] : 0.0;
GpReal *xcp= pass2? site->xcp : 0;
GpReal *ycp= pass2? site->ycp : 0;
int z0, z1, z2, z3;
int keep_left= 0; /* flag to try to minimize curvature in saddles */
int done= 0;
if (level) level= 2;
for (;;) {
/* set edge endpoints */
p0= POINT0(edge,fwd);
p1= POINT1(edge,fwd);
/* always mark cut on current edge */
if (pass2) {
/* second pass actually computes and stores the point */
GpReal zcp= (zlevel-z[p0])/(z[p1]-z[p0]);
xcp[n]= zcp*(x[p1]-x[p0]) + x[p0];
ycp[n]= zcp*(y[p1]-y[p0]) + y[p0];
}
if (!done && !jedge) {
if (n) {
/* if this is not the first point on the curve, and we're
* not done, and this is an i-edge, check several things */
if (!two_levels && !pass2 && (data[edge]&OPEN_END)) {
/* reached an OPEN_END mark, skip the n++ */
done= 4; /* same return value 4 used below */
break;
}
/* check for curve closure -- if not, erase any start mark */
if (edge==edge0 && left==left0) {
/* may signal closure on a downstroke */
if (level0) done= (!pass2 && two_levels && left<0)? 5 : 3;
} else if (!pass2) {
Cdata start= data[edge]&(fwd>0?I0_START:I1_START);
if (start) { data[edge]&=~start; site->count--; }
if (!two_levels) {
start= data[edge]&(fwd>0?I1_START:I0_START);
if (start) { data[edge]&=~start; site->count--; }
}
}
}
}
n++;
if (done) break;
/* cross current zone to another cut edge */
z0= (data[p0]&Z_VALUE) != level; /* 1 if fill toward p0 */
z1= !z0; /* know level cuts edge */
z2= (data[p1+left]&Z_VALUE) != level;
z3= (data[p0+left]&Z_VALUE) != level;
if (z0==z2) {
if (z1==z3) {
/* this is a saddle zone, need triangle to decide
* -- set triangle if not already decided for this zone */
long zone= edge + (left>0? left : 0);
if (triangle) {
if (!triangle[zone]) {
if (keep_left) triangle[zone]= jedge? -1 : 1;
else triangle[zone]= jedge? 1 : -1;
}
if (triangle[zone]>0? !jedge : jedge) goto bkwd;
} else {
if (keep_left) goto bkwd;
}
}
/* bend forward (right along curve) */
keep_left= 1;
jedge= !jedge;
edge= p1 + (left>0? left : 0);
{ long tmp=fwd; fwd=-left; left=tmp; }
} else if (z1==z3) {
bkwd:
/* bend backward (left along curve) */
keep_left= 0;
jedge= !jedge;
edge= p0 + (left>0? left : 0);
{ long tmp=fwd; fwd=left; left=-tmp; }
} else {
/* straight across to opposite edge */
edge+= left;
}
/* after crossing zone, edge/left/fwd is oriented CCW relative to
* the next zone, assuming we will step there */
/* now that we've taken a step, check for the downstroke
* of a slit on the second pass (upstroke checked above)
* -- taking step first avoids a race condition */
if (pass2 && two_levels && !jedge) {
if (left>0) {
if (data[edge]&SLIT_UP) done= 6;
} else {
if (data[edge]&SLIT_DN) done= 5;
}
}
if (!done) {
/* finally, check if we are on a boundary */
if (data[edge] & (jedge?J_BNDY:I_BNDY)) {
done= two_levels? 2 : 4;
/* flip back into the zone that exists */
left= -left;
fwd= -fwd;
if (!pass2 && (edge!=edge0||left!=left0)) {
Cdata start= data[edge]&START_MARK(left);
if (start) { data[edge]&=~start; site->count--; }
}
}
}
}
site->edge= edge;
site->n= n;
site->left= left;
return done>4? slit_cutter(site, data, done-5, pass2) : done;
}
/* edge_walker assumes that the current edge is being drawn CCW
* around the current zone. Since only boundary edges are drawn
* and we always walk around with the filled region to the left,
* no edge is ever drawn CW. We attempt to advance to the next
* edge on this boundary, but if current second endpoint is not
* between the two contour levels, we exit back to zone_crosser.
* Note that we may wind up marking no points.
* -- edge_walker is never called for single level case */
static int edge_walker(Csite *site, Cdata *data, int pass2)
{
long edge= site->edge;
long left= site->left;
long n= site->n;
long fwd= FORWARD(left,site->imax);
long p0= POINT0(edge,fwd);
long p1= POINT1(edge,fwd);
int jedge= IS_JEDGE(edge,left);
long edge0= site->edge0;
long left0= site->left0;
int level0= site->level0==2;
int marked;
const GpReal *x= pass2? site->x : 0;
const GpReal *y= pass2? site->y : 0;
GpReal *xcp= pass2? site->xcp : 0;
GpReal *ycp= pass2? site->ycp : 0;
int z0, z1, heads_up= 0;
for (;;) {
/* mark endpoint 0 only if value is 1 there, and this is a
* two level task */
z0= data[p0]&Z_VALUE;
z1= data[p1]&Z_VALUE;
marked= 0;
if (z0==1) {
/* mark current boundary point */
if (pass2) {
xcp[n]= x[p0];
ycp[n]= y[p0];
}
marked= 1;
} else if (!n) {
/* if this is the first point is not between the levels
* must do the job of the zone_crosser and mark the first cut here,
* so that it will be marked again by zone_crosser as it closes */
if (pass2) {
GpReal zcp= site->zlevel[(z0!=0)];
zcp= (zcp-site->z[p0])/(site->z[p1]-site->z[p0]);
xcp[n]= zcp*(x[p1]-x[p0]) + x[p0];
ycp[n]= zcp*(y[p1]-y[p0]) + y[p0];
}
marked= 1;
}
if (n) {
/* check for closure */
if (level0 && edge==edge0 && left==left0) {
site->edge= edge;
site->left= left;
site->n= n+marked;
/* if the curve is closing on a hole, need to make a downslit */
if (fwd<0 && !(data[edge]&(jedge?J_BNDY:I_BNDY)))
return slit_cutter(site, data, 0, pass2);
return 3;
} else if (pass2) {
if (heads_up || (fwd<0 && (data[edge]&SLIT_DN))) {
site->edge= edge;
site->left= left;
site->n= n+marked;
return slit_cutter(site, data, heads_up, pass2);
}
} else {
/* if this is not first point, clear start mark for this edge */
Cdata start= data[edge]&START_MARK(left);
if (start) { data[edge]&=~start; site->count--; }
}
}
if (marked) n++;
/* if next endpoint not between levels, need to exit to zone_crosser */
if (z1!=1) {
site->edge= edge;
site->left= left;
site->n= n;
return (z1!=0); /* return level closest to p1 */
}
/* step to p1 and find next edge
* -- turn left if possible, else straight, else right
* -- check for upward slit beginning at same time */
edge= p1 + (left>0? left : 0);
if (pass2 && jedge && fwd>0 && (data[edge]&SLIT_UP)) {
jedge= !jedge;
heads_up= 1;
} else if (data[edge]&(jedge?I_BNDY:J_BNDY)) {
long tmp=fwd; fwd=left; left=-tmp;
jedge= !jedge;
} else {
edge= p1 + (fwd>0? fwd : 0);
if (pass2 && !jedge && fwd>0 && (data[edge]&SLIT_UP)) {
heads_up= 1;
} else if (!(data[edge]&(jedge?J_BNDY:I_BNDY))) {
edge= p1 - (left<0? left : 0);
jedge= !jedge;
{ long tmp=fwd; fwd=-left; left=tmp; }
}
}
p0= p1;
p1= POINT1(edge,fwd);
}
}
/* -- slit_cutter is never called for single level case */
static int slit_cutter(Csite *site, Cdata *data, int up, int pass2)
{
long imax= site->imax;
long n= site->n;
const GpReal *x= pass2? site->x : 0;
const GpReal *y= pass2? site->y : 0;
GpReal *xcp= pass2? site->xcp : 0;
GpReal *ycp= pass2? site->ycp : 0;
if (up) {
/* upward stroke of slit proceeds up left side of slit until
* it hits a boundary or a point not between the contour levels
* -- this never happens on the first pass */
long p1= site->edge;
int z1;
for (;;) {
z1= data[p1]&Z_VALUE;
if (z1 != 1) {
site->edge= p1;
site->left= -1;
site->n= n;
return (z1!=0);
} else if (data[p1]&J_BNDY) {
/* this is very unusual case of closing on a mesh hole */
site->edge= p1;
site->left= -imax;
site->n= n;
return 2;
}
xcp[n]= x[p1];
ycp[n]= y[p1];
n++;
p1+= imax;
}
} else {
/* downward stroke proceeds down right side of slit until it
* hits a boundary or point not between the contour levels */
long p0= site->edge;
int z0;
/* at beginning of first pass, mark first i-edge with SLIT_DN */
data[p0]|= SLIT_DN;
p0-= imax;
for (;;) {
z0= data[p0]&Z_VALUE;
if (!pass2) {
if (z0!=1 || (data[p0]&I_BNDY) || (data[p0+1]&J_BNDY)) {
/* at end of first pass, mark final i-edge with SLIT_UP */
data[p0+imax]|= SLIT_UP;
/* one extra count for splicing at outer curve */
site->n= n+1;
return 4; /* return same special value as for OPEN_END */
}
} else {
if (z0 != 1) {
site->edge= p0+imax;
site->left= 1;
site->n= n;
return (z0!=0);
} else if (data[p0+1]&J_BNDY) {
site->edge= p0+1;
site->left= imax;
site->n= n;
return 2;
} else if (data[p0]&I_BNDY) {
site->edge= p0;
site->left= 1;
site->n= n;
return 2;
}
}
if (pass2) {
xcp[n]= x[p0];
ycp[n]= y[p0];
n++;
} else {
/* on first pass need to count for upstroke as well */
n+= 2;
}
p0-= imax;
}
}
}
/* ------------------------------------------------------------------------ */
/* curve_tracer finds the next starting point, then traces the curve,
* returning the number of points on this curve
* -- in a two level trace, the return value is negative on the
* first pass if the curve closed on a hole
* -- in a single level trace, the return value is negative on the
* first pass if the curve is an incomplete open curve
* -- a return value of 0 indicates no more curves */
static long curve_tracer(Csite *site, Cdata *data, int pass2)
{
long imax= site->imax;
long edge0= site->edge0;
long left0= site->left0;
long edge00= site->edge00;
int two_levels= site->zlevel[1]>site->zlevel[0];
int level, level0, mark_row;
long n;
/* it is possible for a single i-edge to serve as two actual start
* points, one to the right and one to the left
* -- for the two level case, this happens on the first pass for
* a doubly cut edge, or on a chunking boundary
* -- for single level case, this is impossible, but a similar
* situation involving open curves is handled below
* a second two start possibility is when the edge0 zone does not
* exist and both the i-edge and j-edge boundaries are cut
* yet another possibility is three start points at a junction
* of chunk cuts
* -- sigh, several other rare possibilities,
* allow for general case, just go in order i1, i0, j1, j0 */
int two_starts;
if (left0==1) two_starts= data[edge0]&(I0_START|J1_START|J0_START);
else if (left0==-1) two_starts= data[edge0]&(J1_START|J0_START);
else if (left0==imax) two_starts= data[edge0]&J0_START;
else two_starts= 0;
if (pass2 || edge0==0) {
/* zip up to row marked on first pass (or by data_init if edge0==0)
* -- but not for double start case */
if (!two_starts) {
/* final start point marked by ALL_DONE marker */
int first= (edge0==0 && !pass2);
long e0= edge0;
if (data[edge0]&ALL_DONE) return 0;
while (!(data[edge0]&START_ROW)) edge0+= imax;
if (e0==edge0) edge0++; /* two starts handled specially */
if (first)
/* if this is the very first start point, we want to remove
* the START_ROW marker placed by data_init */
data[edge0 - edge0%imax]&= ~START_ROW;
}
} else {
/* first pass ends when all potential start points visited */
if (site->count<=0) {
/* place ALL_DONE marker for second pass */
data[edge00]|= ALL_DONE;
/* reset initial site for second pass */
site->edge0= site->edge00= site->left0= 0;
return 0;
}
if (!two_starts) edge0++;
}
if (two_starts) {
/* trace second curve with this start immediately */
if (left0==1 && (data[edge0]&I0_START)) {
left0= -1;
level= (data[edge0]&I_BNDY)? 2 : 0;
} else if ((left0==1 || left0==-1) && (data[edge0]&J1_START)) {
left0= imax;
level= 2;
} else {
left0= -imax;
level= 2;
}
} else {
/* usual case is to scan for next start marker
* -- on second pass, this is at most one row of mesh, but first
* pass hits nearly every point of the mesh, since it can't
* know in advance which potential start marks removed */
while (!(data[edge0]&ANY_START)) edge0++;
if (data[edge0]&I1_START) left0= 1;
else if (data[edge0]&I0_START) left0= -1;
else if (data[edge0]&J1_START) left0= imax;
else /*data[edge0]&J0_START*/ left0= -imax;
if (data[edge0]&(I1_START|I0_START))
level= (data[edge0]&I_BNDY)? 2 : 0;
else
level= 2;
}
/* this start marker will not be unmarked, but it has been visited */
if (!pass2) site->count--;
/* if this curve starts on a non-boundary i-edge, we need to
* determine the level */
if (!level && two_levels)
level= left0>0?
((data[edge0-imax]&Z_VALUE)!=0) : ((data[edge0]&Z_VALUE)!=0);
/* initialize site for this curve */
site->edge= site->edge0= edge0;
site->left= site->left0= left0;
site->level0= level0= level; /* for open curve detection only */
/* single level case just uses zone_crosser */
if (!two_levels) level= 0;
/* to generate the curve, alternate between zone_crosser and
* edge_walker until closure or first call to edge_walker in
* single level case */
site->n= 0;
for (;;) {
if (level<2) level= zone_crosser(site, data, level, pass2);
else if (level<3) level= edge_walker(site, data, pass2);
else break;
}
n= site->n;
/* single level case may have ended at a boundary rather than closing
* -- need to recognize this case here in order to place the
* OPEN_END mark for zone_crosser, remove this start marker,
* and be sure not to make a START_ROW mark for this case
* two level case may close with slit_cutter, in which case start
* must also be removed and no START_ROW mark made
* -- change sign of return n to inform caller */
if (!pass2 && level>3 && (two_levels || level0==0)) {
if (!two_levels) data[edge0]|= OPEN_END;
data[edge0]&= ~(left0>0? I1_START : I0_START);
mark_row= 0; /* do not mark START_ROW */
n= -n;
} else {
if (two_levels) mark_row= !two_starts;
else mark_row= 1;
}
/* on first pass, must apply START_ROW mark in column above previous
* start marker
* -- but skip if we just did second of two start case */
if (!pass2 && mark_row) {
data[edge0 - (edge0-edge00)%imax]|= START_ROW;
site->edge00= edge0;
}
return n;
}
/* ------------------------------------------------------------------------ */
static void data_init(Csite *site, Cdata *data, int region, long nchunk)
{
long imax= site->imax;
long jmax= site->jmax;
long ijmax= imax*jmax;
const GpReal *z= site->z;
GpReal zlev0= site->zlevel[0];
GpReal zlev1= site->zlevel[1];
int two_levels= zlev1>zlev0;
int *reg= site->reg;
long count= 0;
int started= 0;
int ibndy, jbndy, i_was_chunk;
long icsize= imax-1;
long jcsize= jmax-1;
long ichunk, jchunk, irem, jrem, i, j, ij;
if (nchunk && two_levels) {
/* figure out chunk sizes
* -- input nchunk is square root of maximum allowed zones per chunk
* -- start points for single level case are wrong, so don't try it */
long inum= (nchunk*nchunk)/(jmax-1);
long jnum= (nchunk*nchunk)/(imax-1);
if (inum < nchunk) inum= nchunk;
if (jnum < nchunk) jnum= nchunk;
/* ijnum= actual number of chunks,
* ijrem= number of those chunks needing one more zone (ijcsize+1) */
inum= (imax-2)/inum + 1;
icsize= (imax-1)/inum;
irem= (imax-1)%inum;
jnum= (jmax-2)/jnum + 1;
jcsize= (jmax-1)/jnum;
jrem= (jmax-1)%jnum;
/* convert ijrem into value of i or j at which to begin adding an
* extra zone */
irem= (inum-irem)*icsize;
jrem= (inum-jrem)*jcsize;
} else {
irem= imax;
jrem= jmax;
}
/* do everything in a single pass through the data array to
* minimize cache faulting (z, reg, and data are potentially
* very large arrays)
* access to the z and reg arrays is strictly sequential,
* but we need two rows (+-imax) of the data array at a time */
if (z[0]>zlev0) data[0]= (two_levels && z[0]>zlev1)? 2 : 1;
else data[0]= 0;
jchunk= 0;
for (j=ij=0 ; j<jmax ; j++) {
ichunk= i_was_chunk= 0;
for (i=0 ; i<imax ; i++,ij++) {
/* transfer zonal existence from reg to data array
* -- get these for next row so we can figure existence of
* points and j-edges for this row */
data[ij+imax+1]= 0;
if (reg) {
if (region?(reg[ij+imax+1]==region):(reg[ij+imax+1]!=0))
data[ij+imax+1]= ZONE_EX;
} else {
if (i<imax-1 && j<jmax-1) data[ij+imax+1]= ZONE_EX;
}
/* translate z values to 0, 1, 2 flags */
if (ij<imax) data[ij+1]= 0;
if (ij<ijmax-1 && z[ij+1]>zlev0)
data[ij+1]|= (two_levels && z[ij+1]>zlev1)? 2 : 1;
/* apply edge boundary marks */
ibndy= i==ichunk || (data[ij]&ZONE_EX)!=(data[ij+1]&ZONE_EX);
jbndy= j==jchunk || (data[ij]&ZONE_EX)!=(data[ij+imax]&ZONE_EX);
if (ibndy) data[ij]|= I_BNDY;
if (jbndy) data[ij]|= J_BNDY;
/* apply i-edge start marks
* -- i-edges are only marked when actually cut
* -- no mark is necessary if one of the j-edges which share
* the lower endpoint is also cut
* -- no I0 mark necessary unless filled region below some cut,
* no I1 mark necessary unless filled region above some cut */
if (j) {
int v0= (data[ij]&Z_VALUE);
int vb= (data[ij-imax]&Z_VALUE);
if (v0!=vb) { /* i-edge is cut */
if (ibndy) {
if (data[ij]&ZONE_EX) {
data[ij]|= I0_START;
count++;
}
if (data[ij+1]&ZONE_EX) {
data[ij]|= I1_START;
count++;
}
} else {
int va= (data[ij-1]&Z_VALUE);
int vc= (data[ij+1]&Z_VALUE);
int vd= (data[ij-imax+1]&Z_VALUE);
if (v0!=1 && va!=v0 && (vc!=v0 || vd!=v0) &&
(data[ij]&ZONE_EX)) {
data[ij]|= I0_START;
count++;
}
if (vb!=1 && va==vb && (vc==vb || vd==vb) &&
(data[ij+1]&ZONE_EX)) {
data[ij]|= I1_START;
count++;
}
}
}
}
/* apply j-edge start marks
* -- j-edges are only marked when they are boundaries
* -- all cut boundary edges marked
* -- for two level case, a few uncut edges must be marked
*/
if (i && jbndy) {
int v0= (data[ij]&Z_VALUE);
int vb= (data[ij-1]&Z_VALUE);
if (v0!=vb) {
if (data[ij]&ZONE_EX) {
data[ij]|= J0_START;
count++;
}
if (data[ij+imax]&ZONE_EX) {
data[ij]|= J1_START;
count++;
}
} else if (two_levels && v0==1) {
if (data[ij+imax]&ZONE_EX) {
if (i_was_chunk || !(data[ij+imax-1]&ZONE_EX)) {
/* lower left is a drawn part of boundary */
data[ij]|= J1_START;
count++;
}
} else if (data[ij]&ZONE_EX) {
if (data[ij+imax-1]&ZONE_EX) {
/* weird case of open hole at lower left */
data[ij]|= J0_START;
count++;
}
}
}
}
i_was_chunk= (i==ichunk);
if (i_was_chunk) ichunk+= icsize + (ichunk>=irem);
}
if (j==jchunk) jchunk+= jcsize + (jchunk>=jrem);
/* place first START_ROW marker */
if (count && !started) {
data[ij-imax]|= START_ROW;
started= 1;
}
}
/* place immediate stop mark if nothing found */
if (!count) data[0]|= ALL_DONE;
/* initialize site */
site->edge0= site->edge00= site->edge= 0;
site->left0= site->left= 0;
site->n= 0;
site->count= count;
}
/* ------------------------------------------------------------------------ */
/* here are the interface routines for Gist */
extern int GaGetScratchS(long n);
extern short *gasScratch;
static Csite gc_site;
static long gc_common(GaQuadMesh *mesh, int region, const GpReal *zz,
long nchunk, long *nparts);
long GcInit1(GaQuadMesh *mesh, int region, const GpReal *zz,
GpReal lev, long *nparts)
{
gc_site.zlevel[0]= gc_site.zlevel[1]= lev;
return gc_common(mesh, region, zz, 0L, nparts);
}
long GcInit2(GaQuadMesh *mesh, int region, const GpReal *zz,
GpReal levs[2], long nchunk, long *nparts)
{
gc_site.zlevel[0]= levs[0];
gc_site.zlevel[1]= levs[1];
*nparts= 0;
if (gc_site.zlevel[0]==gc_site.zlevel[1]) return 0;
if (gc_site.zlevel[0] > gc_site.zlevel[1]) {
GpReal tmp= gc_site.zlevel[0];
gc_site.zlevel[0]= gc_site.zlevel[1];
gc_site.zlevel[1]= tmp;
}
return gc_common(mesh, region, zz, nchunk, nparts);
}
static long gc_common(GaQuadMesh *mesh, int region, const GpReal *zz,
long nchunk, long *nparts)
{
long ntotal= 0;
long n;
gc_site.imax= mesh->iMax;
gc_site.jmax= mesh->jMax;
gc_site.x= mesh->x;
gc_site.y= mesh->y;
gc_site.reg= mesh->reg;
gc_site.triangle= mesh->triangle;
gc_site.z= zz;
gc_site.xcp= gc_site.ycp= 0;
gc_site.n= gc_site.count= 0;
*nparts= 0;
/* get scratch space for data array */
if (GaGetScratchS(gc_site.imax*(gc_site.jmax+1)+1)) return 0;
/* initialize the data array */
data_init(&gc_site, gasScratch, region, nchunk);
/* make first pass to compute required sizes for GcTrace second pass */
for (;;) {
n= curve_tracer(&gc_site, gasScratch, 0);
if (!n) break;
if (n>0) {
(*nparts)++;
ntotal+= n;
} else {
ntotal-= n;
}
}
return ntotal;
}
long GcTrace(long *n, GpReal *px, GpReal *py)
{
long np, ntotal= 0;
/* make second pass to fill outputs */
for (;;) {
gc_site.xcp= px;
gc_site.ycp= py;
np= curve_tracer(&gc_site, gasScratch, 1);
if (!np) break;
if (np>0) {
*(n++)= np;
px+= np;
py+= np;
ntotal+= np;
} else {
/* serious bug */
ntotal= -1;
break;
}
}
/* data array is likely big -- go ahead and free it */
GaFreeScratch();
return ntotal;
}
/* ------------------------------------------------------------------------ */
syntax highlighted by Code2HTML, v. 0.9.1