#include #include #include #include #include #include #include #include #include static char SccsId[] = "@(#)contour.c 1.12 05/20/05 USGS"; #ifndef DBL_MAX #define DBL_MAX 1.7976931348623157E+308 /* max decimal value of a */ /* "double" */ #endif /* DBL_MAX */ #ifndef FLT_MAX #define FLT_MAX 3.402823466E+38F /* max decimal value of a "float" */ #endif /* FLT_MAX */ #define NEG_HOLES 0 #define POS_HOLES 1 /* * This is the standard GMT grid header, pulled from gmt_grd.h */ struct GRD_HEADER { int nx; /* Number of columns */ int ny; /* Number of rows */ int node_offset; /* 0 for node grids, 1 for pixel grids */ double x_min; /* Minimum x coordinate */ double x_max; /* Maximum x coordinate */ double y_min; /* Minimum y coordinate */ double y_max; /* Maximum y coordinate */ double z_min; /* Minimum z value */ double z_max; /* Maximum z value */ double x_inc; /* x increment */ double y_inc; /* y increment */ double z_scale_factor; /* grd values must be multiplied by this */ double z_add_offset; /* After scaling, add this */ char x_units[80]; /* units in x-direction */ char y_units[80]; /* units in y-direction */ char z_units[80]; /* grid value units */ char title[80]; /* name of data set */ char command[320]; /* name of generating command */ char remark[160]; /* comments re this data set */ }; /* * In this contouring scheme, the grid is divided into a series of * boxes, one grid point on each corner (p0 - p3). A fifth point (p4) * is placed at the center of the box, and is the average of the * the other four points (this is an approximation, and could * probably be improved by fitting a surface). * * The box is then divided into four triangular regions (t0 - t3) * with eight total edges (e0 - e7), as illustrated below: p0 e1 p1 +------------------+ |\ /| | \ t1 / | | \ / | | \e4 e5/ | | \ / | | \ / | | \ / | | \ / | e0 | t0 \/ t2 | e2 | /\p4 | | / \ | | / \ | | / \ | | / \ | | /e7 e6\ | | / \ | | / t3 \ | |/ \| +------------------+ p3 e3 p2 * Contouring takes place by determining if a contour of a particular * interval passes through one of the four outside edges of the box * (e0 - e3), and then is traced through the box until it reaches * another outside edge. For example, if a contour passes through * e0, it must then intersect one of e4 or e7. If it intersects e4, * then it must intersect e1 or e5, and so on. We use linear * interpolation, so no edge will be penetrated by a contour of a * given level more than once. * * The four outside edges can each be given a unique global identifier * (see get_global_edge(), below). By contouring each box, then * connecting the the ends of the contours by means of the global * identifiers, complete contours may be traced through the grid. */ /* * These are the edges that make up each triangle (by index) */ int tri_egs[4][3] = { { 0, 4, 7 }, { 1, 4, 5 }, { 2, 5, 6 }, { 3, 6, 7 } }; /* * These are the points that make up each edge (by index) * Since we would ultimately like to create "right-handed" polygons * (i.e. when moving from point to point the interior of the polygon * is to the right), we want the points for the four outer edges to * be in a known order, so we choose them so that when standing on * the edge and facing inward, the first coordinate is on the right, * and the second is on the left. The point order of the interior * edges doesn't matter. */ struct { int p1; int p2; } eg_pts[] = { { 3, 0 }, { 0, 1 }, { 1, 2 }, { 2, 3 }, { 0, 4 }, { 1, 4 }, { 2, 4 }, { 3, 4 } }; /* * Each interior edge (4-7) belongs to two triangles. Given * one of those triangles, we want to know the other (so that * we may trace the contour from triangle to triangle). This * array is set up so that for a given edge index, and a given * triangle index, the dereferenced value is the other triangle's * index. Got that? */ int eg_tri[8][4] = { { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { 1, 0, -1, -1 }, { -1, 2, 1, -1 }, { -1, -1, 3, 2 }, { 3, -1, -1, 0 } }; typedef struct cpoint { double x; /* x coordinate the contour intersects the edge */ double y; /* y coordinate the contour intersects the edge */ int ti; /* triangle index -- temporary use only */ int ei; /* local edge index of the intersected edge */ int gei; /* global edge index of the intersected edge */ int order; /* indicating the segment's handedness */ struct cpoint *next; /* next point in the segment */ } CPoint; typedef struct cseg { double xmin; double xmax; double ymin; double ymax; int ci; int nesting; int neg_depth; int pos_depth; struct cseg *neg_holes; struct cseg *neg_next; struct cseg *pos_holes; struct cseg *pos_next; struct cseg *tmp_hole; struct cseg *next; CPoint *first; CPoint *last; } CSeg; typedef struct cinfo { double cv; /* value of the contour level */ CSeg *open; /* chain of open contours at this level */ CSeg *closed; /* chain of closed contours at this level */ } CInfo; CInfo *carray = NULL; CInfo contour0; int ncontours = 0; int avecontint = 0; double contour_delta; typedef struct { double x; double y; double z; } XYZPt; struct { double xmin; double xmax; double ymin; double ymax; } gedge; int nfixed = 0; int verbose = 0; int outstyle = 1; int debug = 0; float multiplier = 0.0; FILE *fout; SHPHandle hSHP; DBFHandle hDBF; int nShapeType; int dbfrecord = 0; int polygon_id = 1; int polygon = 2; int nx, ny; double dx, dy, xorig, yorig, zmin, zmax; float *grid; void contour_box(int ix, int iy); void get_p2(XYZPt pt[], double cv, CPoint *p1, CPoint **p2); void fill_cpoint(XYZPt pt[], int ei, int ti, double cv, int order, CPoint **p); void store_segment(CSeg *newseg, int ci); int get_global_edge(int ix, int iy, int ei); CPoint *get_new_cpoint(void); CPoint *dup_cpoint(CPoint *p); CSeg *get_new_cseg(void); void reverse_cseg(CSeg *seg); CSeg *dup_cseg(CSeg *s); CSeg *get_outer_boundry_polygon(void); void fixval(int i, int j, float *g); double epsilon(void); void write_contours(void); void write_polygons(void); void set_depressions(void); int is_contained(CSeg *out, CSeg *in, CPoint *cp); int is_depression(CSeg *out); void set_xyminmax(CSeg *s); void write_polystart(double cv); void write_polyend(void); void write_segstart(void); void write_segment(CSeg *s); void add_contour(double cval); void finish_contours(void); void finish_polygons(void); void connect_segments(CInfo *ca); int grid_edge(double x, double y); int get_polygon_id(void); void clean_up_polygons(void); void nest_polygons(void); int main(int ac, char **av) { int i, j, k; double cv, cval; FILE *fd; struct GRD_HEADER gh; int c; double cincr; double cfirst = 0.0; double clast = FLT_MAX; char *grdfile = NULL; char *contourfile = NULL; char *outfile = NULL; int errflg = 0; int help = 0; int cfile = 0; int cparam = 0; while( (c = getopt(ac, av, "c:C:df:g:hl:o:T:VZ:")) != EOF && !errflg ) { switch( c ) { case 'c': /* contour intervals from file */ if( cparam ) { errflg++; } else if(strlen(optarg) == 0 || (contourfile = strdup(optarg)) == 0) { fprintf(stderr, "Bad argument for flag \"-c\": \"%s\"\n", optarg); errflg++; } else { cfile++; } break; case 'C': /* contour intervals from parameters */ if( cfile ) { errflg++; } else if(strlen(optarg) == 0 || (cincr = atof(optarg)) <= 0) { fprintf(stderr, "Bad argument for flag \"-C\": \"%s\"\n", optarg); errflg++; } else { cparam++; } break; case 'd': /* debug output */ debug++; break; case 'f': /* first contour */ if( cfile ) { errflg++; } else if(strlen(optarg) == 0 || (cfirst = atof(optarg)) < 0) { fprintf(stderr, "Bad argument for flag \"-f\": \"%s\"\n", optarg); errflg++; } break; case 'g': /* grid file (input) */ if(strlen(optarg) == 0 || (grdfile = strdup(optarg)) == 0) { fprintf(stderr, "Bad argument for flag \"-g\": \"%s\"\n", optarg); errflg++; } break; case 'h': /* help -- print documentation */ help++; break; case 'l': /* maximum (last) contour */ if( cfile ) { errflg++; } else if(strlen(optarg) == 0 || (clast = atof(optarg)) <= 0) { fprintf(stderr, "Bad argument for flag \"-f\": \"%s\"\n", optarg); errflg++; } break; case 'o': /* output file */ if(strlen(optarg) == 0 || (outfile = strdup(optarg)) == 0) { fprintf(stderr, "Bad argument for flag \"-o\": \"%s\"\n", optarg); errflg++; } break; case 'T': /* output type: 1, 2, 3 */ if(strlen(optarg) == 0 || (outstyle = atoi(optarg)) <= 0 || outstyle > 4) { fprintf(stderr, "Bad argument for flag \"-T\": \"%s\"\n", optarg); errflg++; } break; case 'V': /* verbose */ verbose++; break; case 'Z': /* data value multiplier */ if(strlen(optarg) == 0) { fprintf(stderr, "Bad argument for flag \"-Z\": \"%s\"\n", optarg); errflg++; } else { multiplier = atof(optarg); } break; case '?': /* error */ default: errflg++; break; } } if( !help && !cfile && !cparam ) { fprintf(stderr, "\nERROR: must specify one of -C or -c.\n"); help++; } if( !help && !grdfile ) { fprintf(stderr, "\nERROR: must specify a grid file with -g.\n"); help++; } if( !help && outstyle == 4 && !outfile ) { fprintf(stderr, "\nERROR: " "must specify an output file with '-o' when using '-T 4'.\n"); help++; } if( help || errflg ) { fprintf(stderr, "\nUsage: contour -g -c | " "-C contour_interval\n" "\t[-f first_contour][-h][-l max_contour][-o ]\n" "\t[-T output_type][-V][-Z mult]\n"); if( errflg ) { fprintf(stderr, "\nUse the -h option for help\n"); exit(1); } } if( help ) { fprintf(stderr,"\nOptions:\n" " -g file\t A GMT grid file in grd format #1\n" " -c file\t A file listing the contour levels, one per line\n" " -C fval\t val specifies the contour interval if -c is not used\n" " -f fval\t Specifies the first contour to draw, which will then be\n" "\t\t incremented by the value given in the -C option\n" " -h\t\t Help; prints this message\n" " -l fval\t Specifies the upper limit to contour; note that this\n" "\t\t value itself is not guaranteed to be contoured, only that\n" "\t\t none greater than it will be\n" " -o file\t Specifies a file to which to write the output; the \n" "\t\t default is stdout\n" " -T ival\t Output type; four values are currently supported:\n" "\t 1\t Output most similar to that of GMT's grdcontour program;\n" "\t\t contour lines consist of a series of lon/lat pairs separated\n" "\t\t by lines beginning with '>' followed by the contour value\n" "\t 2\t Output closed polygons which may or may not overlap;\n" "\t\t contours that end at the edge of the grid will follow the\n" "\t\t grid border until they reach a matching contour (which is\n" "\t\t a likely place for contours to overlap if, indeed, they do)\n" "\t 3\t Output non-overlapping polygons, with possible holes;\n" "\t\t this format is used for creating 'clean' polygons for the\n" "\t\t generation of ArcView shape files\n" "\t 4\t Output ArcView shape files; it this option is selected\n" "\t\t an output file must be specified with the '-o' option, and\n" "\t\t the file name given will be the root name for the .shp, .shx\n" "\t\t and .dbf files.\n" ); exit(0); } /* * Using a fixed seed means we get the same random numbers * every time, but that gives us the same results every time * (which is good) and I don't think the systematic errors it * introduces are a problem. */ srand48(444); /* srand48(getpid()); */ /* * Read the grid header, set some useful variables */ if( (fd = fopen(grdfile, "rb")) == NULL ) { fprintf( stderr, "Unable to open %s\n", grdfile); exit(1); } if( fread(&gh, sizeof(gh), 1, fd) != 1 ) { fprintf( stderr, "Unable to read header from %s\n", grdfile); exit(1); } nx = gh.nx; ny = gh.ny; dx = gh.x_inc; dy = -1 * gh.y_inc; xorig = gh.x_min; yorig = gh.y_max; zmin = gh.z_min; zmax = gh.z_max; /* * These values are at the borders of the grid * Note: the strange calculation (i.e. (nx-2)*dx then later * adding dx again) is to duplicate the way this value will * be calculated later on, and since we compare this value to * another floating point value for equality, we don't want to * introduce even the smallest round-off difference. */ if( outstyle == 2 || outstyle == 3 || outstyle == 4 ) { gedge.xmin = xorig; gedge.xmax = (nx - 2) * dx + xorig; gedge.xmax += dx; gedge.ymax = yorig; gedge.ymin = (ny - 2) * dy + yorig; gedge.ymin += dy; } /* * Allocate the grid array and read the grid into it */ if( (grid = (float *)malloc(nx * ny * sizeof(float))) == NULL) { fprintf( stderr, "Can't allocate memory for grid\n"); exit(1); } if( fread(grid, sizeof(float), nx * ny, fd) != nx * ny ) { fprintf( stderr, "Error reading grid from %s\n", grdfile); exit(1); } fclose(fd); if( verbose ) fprintf(stderr, "Read grid file.\n"); /* * look for NaNs, set them to something sane */ for(j = 0; j < ny; j++) { for(i = 0; i < nx; i++) { if (isnanf(grid[i + j*nx]) ) fixval(i, j, grid); } } if( verbose ) fprintf(stderr, "Fixed %d NaNs\n", nfixed); /* * multiply by the value supplied with -Z */ if( multiplier != 0.0 ) { for(j = 0; j < ny; j++) { for(i = 0; i < nx; i++) { grid[i + j*nx] *= multiplier; } } zmin *= multiplier; zmax *= multiplier; } /* * get contour info */ if( cfile ) { if( (fd = fopen(contourfile, "r")) == NULL ) { fprintf( stderr, "Unable to open %s\n", contourfile); exit(1); } while( fscanf(fd, "%lf\n", &cval) != EOF ) { if( cval < zmin || cval > zmax ) continue; add_contour(cval); } if( ncontours == 0 ) { rewind(fd); if( fscanf(fd, "%lf\n", &cval) == EOF ) { fprintf(stderr, "No contours in contourfile %s\n", contourfile); exit(1); } add_contour(cval); } fclose(fd); if( verbose ) fprintf(stderr, "Read contour file.\n"); } else { for(cval = cfirst; cval <= clast; cval += cincr) { if( cval < zmin ) continue; if( cval > zmax ) break; add_contour(cval); } if( ncontours == 0 ) add_contour(cfirst); if( verbose ) fprintf(stderr, "Set up contour levels.\n"); } /* * Compute the average contour interval */ if( ncontours == 1 ) { avecontint = carray[0].cv; } else { for( i = 1; i < ncontours; i++ ) { avecontint += abs(carray[i].cv - carray[i - 1].cv); } avecontint /= (ncontours - 1); } if( avecontint == 0 ) avecontint = 1.0; contour_delta = avecontint * 1.e-4; /* * Open output file for writing, or use stdout... */ if( outfile ) { if( outstyle == 4 ) { if( (hSHP = SHPCreate(outfile, SHPT_POLYGON)) == NULL ) { fprintf( stderr, "Unable to create shapefile %s\n", outfile); exit(1); } SHPGetInfo(hSHP, NULL, &nShapeType, NULL, NULL); /* SHPClose(hSHP); */ if( (hDBF = DBFCreate(outfile)) == NULL ) { fprintf( stderr, "Unable to create .dbf file %s.dbf\n", outfile); exit(1); } if( DBFAddField(hDBF, "AREA", FTDouble, 13, 3) == -1 || DBFAddField(hDBF, "PERIMETER", FTDouble, 14, 3) == -1 || DBFAddField(hDBF, "PGAPOL_", FTDouble, 12, 0) == -1 || DBFAddField(hDBF, "PGAPOL_ID", FTDouble, 12, 0) == -1 || DBFAddField(hDBF, "GRID_CODE", FTDouble, 12, 0) == -1 || DBFAddField(hDBF, "VALUE", FTDouble, 14, 4) == -1 ) { fprintf(stderr, "DBFAddField failed...\n"); exit(1); } dbfrecord = DBFGetRecordCount(hDBF); /* DBFClose(hDBF); */ } else { if( (fout = fopen(outfile, "w")) == NULL ) { fprintf( stderr, "Unable to open %s\n", outfile); exit(1); } } } else { fout = stdout; } /* * We can't handle grid points that fall exactly on contours, * so we nudge them off, just a bit... */ for(k = 0; k < ncontours; k++) { cv = carray[k].cv; for(j = 0; j < ny; j++) { for(i = 0; i < nx; i++) { if( grid[i + j*nx] > cv - contour_delta && grid[i + j*nx] < cv + contour_delta ) { grid[i + j*nx] += epsilon(); /* * We want to recheck this point to make sure we changed * the value within the resolution of the machine's * floating point representation. */ i--; } } } } if( verbose ) fprintf(stderr, "Moved nodes off contour values.\n"); /* * Do the actual contouring */ for(j = 0; j < ny - 1; j++) { for(i = 0; i < nx - 1; i++) { contour_box(i, j); } } if( verbose ) fprintf(stderr, "Did contouring.\n"); /* * output the results */ if( outstyle == 1 ) { write_contours(); } else if( outstyle == 2 ) { finish_contours(); if( verbose ) fprintf(stderr, "Finished contours.\n"); write_contours(); if( verbose ) fprintf(stderr, "Wrote contours.\n"); } else { finish_polygons(); if( verbose ) fprintf(stderr, "Finished polygons.\n"); write_polygons(); if( verbose ) fprintf(stderr, "Wrote polygons.\n"); if(outstyle == 4) { SHPClose(hSHP); DBFClose(hDBF); } } if( verbose ) fprintf(stderr, "Done.\n"); free(grid); free(carray); return 0; } /* * ci => contour index * cv => contour value * ti => triangle index * ei => edge index */ void contour_box(int ix, int iy) { int lp1, lp2, ci, ei, order, tmp; int done[4]; double cv; CSeg *seg; XYZPt pt[5]; CPoint *p1, *p2, *last; double x = ix * dx + xorig; double y = iy * dy + yorig; pt[0].x = x; pt[1].x = x + dx; pt[2].x = x + dx; pt[3].x = x; pt[4].x = x + dx / 2.0; pt[0].y = y; pt[1].y = y; pt[2].y = y + dy; pt[3].y = y + dy; pt[4].y = y + dy / 2.0; pt[0].z = grid[iy * nx + ix]; pt[1].z = grid[iy * nx + ix + 1]; pt[2].z = grid[(iy + 1) * nx + ix + 1]; pt[3].z = grid[(iy + 1) * nx + ix]; pt[4].z = (pt[0].z + pt[1].z + pt[2].z + pt[3].z) / 4.0; /* * Loop over the contour values */ for(ci = 0; ci < ncontours; ci++) { cv = carray[ci].cv; if( pt[4].z > cv - contour_delta && pt[4].z < cv + contour_delta ) { pt[4].z += epsilon(); ci--; continue; } done[0] = done[1] = done[2] = done[3] = 0; /* * Loop over the four outside edges of this block */ for(ei = 0; ei < 4; ei++) { if( done[ei] ) continue; lp1 = eg_pts[ei].p1; lp2 = eg_pts[ei].p2; if( pt[lp1].z >= pt[lp2].z ) { /* normal order */ order = 1; } else { /* reverse order */ tmp = lp1; lp1 = lp2; lp2 = tmp; order = -1; } if( cv > pt[lp1].z || cv < pt[lp2].z ) continue; /* * Fill the first point with initial info; note that this * is an outside edge, so ei == ti */ fill_cpoint(pt, ei, ei, cv, order, &p1); last = p1; /* * global edges only matter on the anchor points of a segment */ p1->gei = get_global_edge(ix, iy, ei); /* *while p2 not on edge 0 -> 3, continue through the current box */ do { get_p2(pt, cv, last, &p2); /* add p2 to current segment */ last->next = p2; last = p2; } while (last->ei > 3); last->gei = get_global_edge(ix, iy, last->ei); /* * segment is now complete, store it */ done[last->ei] = 1; seg = get_new_cseg(); seg->first = p1; seg->last = last; seg->next = carray[ci].open; carray[ci].open = seg; store_segment(seg, ci); } } return; } void get_p2(XYZPt pt[], double cv, CPoint *p1, CPoint **p2) { int i, ei, lp1, lp2, tmp, order; for( i = 0; i < 3; i++ ) { if((ei = tri_egs[p1->ti][i]) == p1->ei) continue; lp1 = eg_pts[ei].p1; lp2 = eg_pts[ei].p2; if( pt[lp1].z >= pt[lp2].z ) { /* normal order */ order = 1; } else { /* reverse order */ tmp = lp1; lp1 = lp2; lp2 = tmp; order = -1; } if( cv > pt[lp1].z || cv < pt[lp2].z ) continue; fill_cpoint(pt, ei, p1->ti, cv, order, p2); break; } if( i >= 3 ) { fprintf( stderr, "get_p2: couldn't find another crossing point\n"); exit(1); } /* * He we set up the next iteration: if this edge is an interior * edge, we will need to do the next computation on the adjacent * triangle, rather than the current one. */ if( (*p2)->ei > 3 ) (*p2)->ti = eg_tri[(*p2)->ei][(*p2)->ti]; return; } void fill_cpoint(XYZPt pt[], int ei, int ti, double cv, int order, CPoint **p) { int lp1, lp2; double f; *p = get_new_cpoint(); (*p)->order = order; (*p)->ti = ti; (*p)->ei = ei; if(order > 0) { lp1 = eg_pts[ei].p1; lp2 = eg_pts[ei].p2; } else { lp1 = eg_pts[ei].p2; lp2 = eg_pts[ei].p1; } f = (pt[lp1].z - cv) / (pt[lp1].z - pt[lp2].z); if(ei == 0 || ei == 2) { /* vertical edge */ (*p)->x = pt[lp1].x; (*p)->y = pt[lp2].y * f + pt[lp1].y * (1.0 - f); } else if(ei == 1 || ei == 3) { /* horizontal edge */ (*p)->x = pt[lp2].x * f + pt[lp1].x * (1.0 - f); (*p)->y = pt[lp1].y; } else { /* diagonal edge */ (*p)->x = pt[lp2].x * f + pt[lp1].x * (1.0 - f); (*p)->y = pt[lp2].y * f + pt[lp1].y * (1.0 - f); } return; } /* * store_segment tries to add a segment to an existing segment, if * no current segment matches, it starts a new segment. If, by * adding the new segment, a segment becomes closed, the segment * is moved to the closed segment list. */ void store_segment(CSeg *newseg, int ci) { CSeg *seg, *prev; int didmerge = 1; if( newseg->first->order < 0 ) reverse_cseg(newseg); /* * Try to find a matching segment on the open segment list */ while( didmerge ) { for( prev = 0, seg = carray[ci].open; seg; prev = seg, seg = seg->next ) { if( seg == newseg ) continue; if( newseg->first->gei == seg->last->gei ) { /* * add the existing segment to the front of the new one * skip the first point in the new segment since it is a duplicate * WARNING: memory leak -- lose track of 1 CPoint */ seg->last->next = newseg->first->next; newseg->first = seg->first; break; } else if( newseg->last->gei == seg->first->gei ) { /* * add new segment to the front of the existing one * skip the first point in the old segment since it is a duplicate * WARNING: memory leak -- lose track of 1 CPoint */ newseg->last->next = seg->first->next; newseg->last = seg->last; break; } } if( seg ) { /* * We did a merge, so first we want to get rid of the old * segment, which is 'seg' */ if( prev ) prev->next = seg->next; else carray[ci].open = seg->next; free(seg); /* * Now see if the merged segment is closed */ if( newseg->first->gei == newseg->last->gei ) { /* * pop the segment out of the list and add it to the closed list */ carray[ci].open = newseg->next; newseg->next = carray[ci].closed; carray[ci].closed = newseg; /* * Yes, we did merge, but that segment is done, so there's * nothing left to do here */ didmerge = 0; } } else { didmerge = 0; } } return; } /* * egmap remaps the edge indices so that they work within the * get_global_edge function */ int egmap[] = { 0, 2, 1, 3 }; /* * f1eg and f2eg are functions that are used by get_global_edge * to compute global edge numbers */ int f1eg[] = { 0, 0, 0, 1 }; int f2eg[] = { 0, 0, 1, 1 }; /* * get_global_edge(ix, iy, ei) returns a unique numeric identifier for * each edge in the grid. For edges that overlap (i.e. the right edge * of a box overlaps the left edge of the box to its right), the * identifier will be the same. 'ix', 'iy' are the grid indices of the * upper left corner of the box, and ei is the index of the edge in * question. If you know a better way to do this, let me know. */ int get_global_edge(int ix, int iy, int ei) { int mi = egmap[ei]; return 2 * (mi + ix + iy * nx + (nx - 1) * f1eg[mi]) - 3 * f2eg[mi]; } #define NGRAB 1000 CPoint *get_new_cpoint(void) { static int cp = NGRAB; static CPoint *buf; if(cp == NGRAB) { buf = (CPoint *)malloc(NGRAB * sizeof(CPoint)); if (buf == NULL) { fprintf(stderr, "Out of memory in get_new_cpoint\n"); exit(1); } memset(buf, 0, NGRAB * sizeof(CPoint)); cp = 0; } return &buf[cp++]; } CPoint *dup_cpoint(CPoint *p) { CPoint *np = get_new_cpoint(); memcpy(np, p, sizeof(CPoint)); return np; } CSeg *get_new_cseg() { CSeg *s = (CSeg *)malloc(sizeof(CSeg)); if(s == NULL) { fprintf(stderr, "Out of memory in get_new_cseg\n"); exit(1); } s->nesting = 0; s->neg_depth = INT_MAX; s->pos_depth = INT_MAX; s->neg_holes = NULL; s->neg_next = NULL; s->pos_holes = NULL; s->pos_next = NULL; s->next = NULL; s->first = NULL; s->last = NULL; return s; } void reverse_cseg(CSeg *seg) { CPoint *s = seg->first; CPoint *tmp, *top; tmp = s->next; s->next = NULL; top = s; s = tmp; while( s->next ) { tmp = s->next; s->next = top; top = s; s = tmp; } s->next = top; seg->last = seg->first; seg->first = s; return; } CSeg *dup_cseg(CSeg *s) { CSeg *new = get_new_cseg(); CPoint *ncp = dup_cpoint(s->first); CPoint *cp = s->first->next; new->first = ncp; while( cp ) { ncp->next = dup_cpoint(cp); ncp = ncp->next; cp = cp->next; } new->last = ncp; new->nesting = s->nesting; new->xmin = s->xmin; new->xmax = s->xmax; new->ymin = s->ymin; new->ymax = s->ymax; return new; } CSeg *get_outer_boundry_polygon() { CSeg *s = get_new_cseg(); CPoint *p; p = get_new_cpoint(); s->first = p; p->x = gedge.xmin; p->y = gedge.ymax; p->next = get_new_cpoint(); p = p->next; p->x = gedge.xmax; p->y = gedge.ymax; p->next = get_new_cpoint(); p = p->next; p->x = gedge.xmax; p->y = gedge.ymin; p->next = get_new_cpoint(); p = p->next; p->x = gedge.xmin; p->y = gedge.ymin; p->next = get_new_cpoint(); p = p->next; p->x = gedge.xmin; p->y = gedge.ymax; s->last = p; return s; } void fixval(int i, int j, float *g) { int cnt = 0; int val = 0; float t; if( i == 0 ) { if( j == 0 ) { val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + (j + 1)*nx]) ? 0 : (cnt++, t); } else if( j == ny - 1 ) { val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); } else { val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); } } else if( i == nx - 1 ) { if( j == 0 ) { val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); } else if( j == ny - 1 ) { val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); } else { val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); } } else if( j == 0 ) { val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); } else if( j == ny - 1 ) { val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); } else { val += isnanf(t = g[i + (j - 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + (j + 1)*nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i + 1 + j *nx]) ? 0 : (cnt++, t); val += isnanf(t = g[i - 1 + j *nx]) ? 0 : (cnt++, t); } if( cnt > 0 ) g[i + j*nx] = val / cnt; else g[i + j*nx] = 0.0; nfixed++; return; } double epsilon(void) { return (drand48() - 0.5) * 20.0 * contour_delta; } void write_contours(void) { int i; double cv; CSeg *s; CPoint *cp; for( i = 0; i < ncontours; i++ ) { cv = carray[i].cv; for( s = carray[i].closed; s; s = s->next ) { fprintf(fout, "> %f\n", cv); for( cp = s->first; cp; cp = cp->next ) { fprintf(fout, "%f\t%f\n", cp->x, cp->y); } } for( s = carray[i].open; s; s = s->next ) { fprintf(fout, "> %f\n", cv); for( cp = s->first; cp; cp = cp->next ) { fprintf(fout, "%f\t%f\n", cp->x, cp->y); } } } return; } void clean_up_polygons(void) { int ci; int shorts = 0; int fixes = 0; CInfo *clist; CSeg *s, *sprev; CPoint *first; for( ci = -1; ci < ncontours; ci++ ) { if( ci == -1 ) clist = &contour0; else clist = &(carray[ci]); for( s = clist->closed; s; s = s->next ) { first = s->first; if( first->x != s->last->x || first->y != s->last->y ) { if(verbose) fprintf(stderr, "Points not identical in contour %d (%g %g)\n", ci, first->x - s->last->x, first->y - s->last->y); s->last->x = first->x; s->last->y = first->y; fixes++; } } /* * Get rid of short contours (i.e. < 4 points) */ for( s = clist->closed, sprev = NULL; s; ) { if( !s->first || !s->first->next || !s->first->next->next || !s->first->next->next->next ) { if( sprev ) { sprev->next = s->next; free(s); s = sprev->next; } else { clist->closed = s->next; free(s); s = clist->closed; } shorts++; } else { sprev = s; s = s->next; } } } if( verbose ) fprintf(stderr, "Fixed %d mismatches, discarded %d short polygons\n", fixes, shorts); return; } void write_polygons(void) { int i; CInfo *clist; CSeg *s, *t; for( i = -1; i < ncontours; i++ ) { if( i == -1 ) clist = &contour0; else clist = &(carray[i]); for( s = clist->closed; s; s = s->next ) { write_polystart(clist->cv); write_segment(s); for( t = s->neg_holes; t; t = t->neg_next ) { reverse_cseg(t); write_segstart(); write_segment(t); reverse_cseg(t); } for( t = s->pos_holes; t; t = t->pos_next ) { reverse_cseg(t); write_segstart(); write_segment(t); reverse_cseg(t); } write_polyend(); } } return; } double *xpts, *ypts; int nvertsmax, nverts; int panParts[5000]; int nparts; void write_polystart(double cv) { if( outstyle == 4 ) { nverts = 0; nvertsmax = 1000; panParts[0] = 0; nparts = 1; if( (xpts = (double *)malloc(nvertsmax * sizeof(double))) == NULL || (ypts = (double *)malloc(nvertsmax * sizeof(double))) == NULL ) { fprintf(stderr,"Out of memory in write_polystart\n"); exit(1); } DBFWriteDoubleAttribute(hDBF, dbfrecord, 0, 0); DBFWriteDoubleAttribute(hDBF, dbfrecord, 1, 0); DBFWriteDoubleAttribute(hDBF, dbfrecord, 2, polygon); DBFWriteDoubleAttribute(hDBF, dbfrecord, 3, polygon_id); if( multiplier != 0 ) { DBFWriteDoubleAttribute(hDBF, dbfrecord, 4, cv / multiplier); } else { DBFWriteDoubleAttribute(hDBF, dbfrecord, 4, cv); } DBFWriteDoubleAttribute(hDBF, dbfrecord, 5, cv); polygon++; polygon_id++; dbfrecord++; } else { fprintf(fout, "> %f\n", cv); } return; } void write_polyend(void) { SHPObject *psObject; if( outstyle == 4 ) { psObject = SHPCreateObject(nShapeType, -1, nparts, panParts, NULL, nverts, xpts, ypts, NULL, NULL); SHPWriteObject(hSHP, -1, psObject); SHPDestroyObject(psObject); free(xpts); free(ypts); } return; } void write_segstart(void) { if( outstyle == 4 ) { panParts[nparts++] = nverts; } else { fprintf(fout, "+\n"); } return; } void write_segment(CSeg *s) { CPoint *cp; if( outstyle == 4 ) { for( cp = s->first; cp; cp = cp->next ) { if( nverts == nvertsmax ) { nvertsmax *= 2; if((xpts = (double*)realloc(xpts, nvertsmax*sizeof(double))) == NULL ||(ypts = (double*)realloc(ypts, nvertsmax*sizeof(double))) == NULL) { fprintf(stderr, "Out of memory in write_segment\n"); exit(1); } } xpts[nverts] = cp->x; ypts[nverts] = cp->y; nverts++; } } else { for( cp = s->first; cp; cp = cp->next ) fprintf(fout, "%f\t%f\n", cp->x, cp->y); } return; } void set_xyminmax(CSeg *s) { CPoint *cp; while( s ) { cp = s->first; s->xmin = s->xmax = cp->x; s->ymin = s->ymax = cp->y; for( cp = cp->next; cp; cp = cp->next ) { if( cp->x < s->xmin ) s->xmin = cp->x; if( cp->x > s->xmax ) s->xmax = cp->x; if( cp->y < s->ymin ) s->ymin = cp->y; if( cp->y > s->ymax ) s->ymax = cp->y; } s = s->next; } return; } int is_contained(CSeg *out, CSeg *in, CPoint *cp) { CPoint *p1, *p2; int inside = 0; /* * A polygon can't be contained within another one if any point * is outside a box bounding the container */ if( cp->x < out->xmin || cp->x > out->xmax || cp->y < out->ymin || cp->y > out->ymax ) { return 0; } /* * Find a point on the polygon to be tested that does not lie * atop a point in the container -- a polygon can contain another * if they share a vertex, but not if they share an edge (this * is a rule for ArcView shapefiles). */ for( p1 = out->first, p2 = out->last; p1; p2 = p1, p1 = p1->next ) { if( p1->x == cp->x && p1->y == cp->y ) { cp = cp->next; if( !cp ) return 0; if( ( p2->x == cp->x && p2->y == cp->y ) || ( p1->next->x == cp->x && p1->next->y == cp->y )) { cp = 0; break; } } } if( !cp ) return 0; /* * Might have a new point, so check the box again (anything * to short-circuit this routine...) */ if( cp->x < out->xmin || cp->x > out->xmax || cp->y < out->ymin || cp->y > out->ymax ) { return 0; } /* * If any vertex of a polygon is inside the other, then the whole * polygon is contained in the other. * * Count the edges intersecting a ray shot from the point: an odd * number of intersections means the point is inside the polygon, * even means it is outside. I found this algorithm in some damn * book or another. */ for( p1 = out->first, p2 = out->last; p1; p2 = p1, p1 = p1->next ) { if( ( ( p1->y <= cp->y && cp->y < p2->y ) || ( p2->y <= cp->y && cp->y < p1->y ) ) && ( cp->x < (p2->x - p1->x) * (cp->y - p1->y) / (p2->y - p1->y) + p1->x ) ) inside = !inside; } return inside ? 1 : 0; } int is_depression(CSeg *s) { CPoint *p1, *p2, *pmin; double y1, y2, x1, x2; double x, xmin, yval; /* * Find the leftmost edge of the polygon (x == xmin) at some * fixed value of y. If, along that edge, the interior of the * polygon is to the left, the polygon is a depression, otherwise * it is a normal polygon. */ for( p1 = s->first, p2 = p1->next; p2; p1 = p2, p2 = p1->next ) { if( p1->y == p2->y ) continue; yval = (p1->y + p2->y) / 2.0; xmin = (p1->x + p2->x) / 2.0; pmin = p1; break; } if( !p2 ) { fprintf(stderr, "Error in is_depression: can't find a point\n"); return 0; } for( p1 = p2, p2 = p1->next; p2; p1 = p2, p2 = p1->next ) { if( p1->y == p2->y ) { continue; } else if( p1->y < p2->y ) { y1 = p1->y; x1 = p1->x; y2 = p2->y; x2 = p2->x; } else { y2 = p1->y; x2 = p1->x; y1 = p2->y; x1 = p2->x; } if( yval >= y2 || yval < y1 ) continue; x = x2 + (x1 - x2) * ((y2 - yval) / (y2 - y1)); if( x < xmin ) { pmin = p1; xmin = x; } } if( pmin->y < pmin->next->y ) return 0; else return 1; } void add_contour(double cval) { static int ncarray = 0; if( ncontours >= ncarray ) { if( ncarray == 0 ) ncarray = 50; ncarray *= 2; carray = (CInfo *)realloc(carray, ncarray * sizeof(CInfo)); if( carray == NULL ) { fprintf( stderr, "Can't reallocate memory for carray\n"); exit(1); } } carray[ncontours].cv = cval; carray[ncontours].open = NULL; carray[ncontours].closed = NULL; ncontours++; return; } void finish_contours(void) { int ci; for( ci = 0; ci < ncontours; ci++ ) { connect_segments( &(carray[ci]) ); } return; } void finish_polygons(void) { int ci, ciopen, ciclosed; CSeg *nseg, *s; /* * Find the lowest levels with open and with closed contours */ ciopen = -1; ciclosed = -1; for( ci = 0; ci < ncontours - 1; ci++ ) { if( ciopen == -1 && carray[ci].open != NULL ) ciopen = ci; if( ciclosed == -1 && carray[ci].closed != NULL ) ciclosed = ci; } if( ciopen == -1 ) { /* * If there are no open contours then the outside border of the map * won't be included. So we dummy up a set of points to create the * enclosing polygon. */ s = get_outer_boundry_polygon(); if( ciclosed <= 0 ) { s->next = contour0.closed; contour0.closed = s; } else { s->next = carray[ciclosed - 1].closed; carray[ciclosed - 1].closed = s; } } /* * Close the polygons -- this is similar to closing the contours * (i.e. finish_contours()) except that contours of the next larger * level can close off a polygon, so we duplicate and reverse them, * (which is, essentially, the same as going to the left of the * *first* point instead of the last) and put them on the open list * for the current level. */ /* * The lowest level segments go to contour0... */ for( s = carray[0].open; s; s = s->next ) { nseg = dup_cseg(s); reverse_cseg(nseg); nseg->next = contour0.open; contour0.open = nseg; } connect_segments(&contour0); /* * Need to set contour0.cv -- we know it is less than the * first contour, but we have no way of knowing by how * much, so we set it to the value of the first contour */ contour0.cv = carray[0].cv; for( ci = 0; ci < ncontours - 1; ci++ ) { for( s = carray[ci+1].open; s; s = s->next ) { nseg = dup_cseg(s); reverse_cseg(nseg); nseg->next = carray[ci].open; carray[ci].open = nseg; } connect_segments( &(carray[ci]) ); /* * Values inside a polygon are assumed to be halfway between * the inner and outer edges (contours) */ carray[ci].cv = (carray[ci].cv + carray[ci+1].cv) / 2.0; } /* * Simply close the highest level contours. * We leave the contour value alone, because, as with the lowest * contour, we know the enclosed area has a higher value, but we * have no way of knowing how much higher. */ connect_segments( &(carray[ncontours - 1]) ); clean_up_polygons(); set_xyminmax(contour0.closed); for( ci = 0; ci < ncontours; ci++ ) set_xyminmax(carray[ci].closed); set_depressions(); nest_polygons(); return; } void connect_segments(CInfo *ca) { int fei, lei; double mmpt, curval; CSeg *this, *fseg, *prev; CSeg *curseg, *curprv; CPoint *newpt; while( ca->open ) { fseg = ca->open; if( (lei = grid_edge(fseg->last->x, fseg->last->y)) < 0 ) { /* * This can only happen if an open contour does not reach the * edge of the grid, which is an unrecoverable error and indicates * a bug in the program... */ fprintf(stderr, "Error contouring val %f: gei1 %d " "gei2 %d\n", ca->cv, fseg->first->gei, fseg->last->gei); fprintf(stderr, "UNKNOWN edge x=%f y=%f xmin=%f xmax=%f ymin=%f " "ymax=%f\n", fseg->last->x, fseg->last->y, gedge.xmin, gedge.xmax, gedge.ymin, gedge.ymax); exit(1); } if( lei == 0 || lei == 2 ) mmpt = fseg->last->y; else mmpt = fseg->last->x; curval = DBL_MAX; curseg = NULL; curprv = NULL; /* * Inner loop over possible mates */ for( this = fseg, prev = NULL; this; prev = this, this = this->next ) { if( (fei = grid_edge(this->first->x, this->first->y)) < 0 ) { fprintf(stderr, "Error contouring val %f: gei1 %d " "gei2 %d\n", ca->cv, this->first->gei, this->last->gei); fprintf(stderr, "first edge %d, last edge %d\n", fei, lei); fprintf(stderr, "UNKNOWN edge x=%f y=%f xmin=%f xmax=%f ymin=%f " "ymax=%f\n", this->first->x, this->first->y, gedge.xmin, gedge.xmax, gedge.ymin, gedge.ymax); exit(1); } if( fei != lei ) continue; if( lei == 0 ) { /* want smallest y > y0 */ if( this->first->y < mmpt ) continue; if( this->first->y < curval || curseg == NULL ) { curval = this->first->y; curseg = this; curprv = prev; } } else if( lei == 1 ) { /* want smallest x > x0 */ if( this->first->x < mmpt ) continue; if( this->first->x < curval || curseg == NULL ) { curval = this->first->x; curseg = this; curprv = prev; } } else if( lei == 2 ) { /* want greatest y < y0 */ if( this->first->y > mmpt ) continue; if( this->first->y > curval || curseg == NULL ) { curval = this->first->y; curseg = this; curprv = prev; } } else { /* want greatest x < x0 */ if( this->first->x > mmpt ) continue; if( this->first->x > curval || curseg == NULL ) { curval = this->first->x; curseg = this; curprv = prev; } } } if( curseg ) { /* * found something? */ if( curseg == fseg ) { /* * connect last point to first point, then put on closed list */ newpt = dup_cpoint(fseg->first); fseg->last->next = newpt; fseg->last = newpt; newpt->next = NULL; ca->open = fseg->next; fseg->next = ca->closed; ca->closed = fseg; } else { /* * Connect the matching segment to the first segment, * remove it from the open list, then go back to the loop */ fseg->last->next = curseg->first; fseg->last = curseg->last; curprv->next = curseg->next; free(curseg); } } else { /* * No point to connect to on this edge; make a new point * located at the leftmost point on the edge and add it * to the segment; (the next loop iteration will view * this point as being on the next sequential edge -- see * the grid_edge() function to see how this works */ newpt = get_new_cpoint(); if( lei == 0 ) { newpt->x = gedge.xmin; newpt->y = gedge.ymax; } else if( lei == 1 ) { newpt->x = gedge.xmax; newpt->y = gedge.ymax; } else if( lei == 2 ) { newpt->x = gedge.xmax; newpt->y = gedge.ymin; } else { newpt->x = gedge.xmin; newpt->y = gedge.ymin; } fseg->last->next = newpt; fseg->last = newpt; newpt->next = NULL; } } return; } void set_depressions(void) { int ci; CSeg *s, *prev; for( s = contour0.closed; s; s = s->next ) { if( is_depression(s) ) { fprintf(stderr, "Found a depression in the lowest contour...\n"); exit(1); } } for( ci = 0; ci < ncontours; ci++ ) { for( s = carray[ci].closed, prev = NULL; s; ) { if( is_depression(s) ) { /* * Depressions in this level are simply polygons at the level below... * Reverse the polygon: */ reverse_cseg(s); /* * Remove the polygon from the current list */ if( prev ) prev->next = s->next; else carray[ci].closed = s->next; /* * Put the polygon on the next-lower list */ if( ci == 0 ) { s->next = contour0.closed; contour0.closed = s; } else { s->next = carray[ci-1].closed; carray[ci-1].closed = s; } if( prev ) s = prev->next; else s = carray[ci].closed; } else { prev = s; s = s->next; } } } return; } /* * Given the coordinates of a point on the edge of the grid, returns * the index of the edge upon which the point lies */ int grid_edge(double x, double y) { /* * Normally, you are on the x/y max/min of your coordinate, but * the leftmost point (looking inward from outside of the box) * should be considered on the following edge. */ if( x == gedge.xmin ) { if( y == gedge.ymax ) return 1; return 0; } else if( y == gedge.ymax ) { if( x == gedge.xmax ) return 2; return 1; } else if( x == gedge.xmax ) { if( y == gedge.ymin ) return 3; return 2; } else if( y == gedge.ymin ) { /* x == gedge.xmin has already been handled */ return 3; } else { return -1; } } void nest_polygons(void) { int i, j, best; CInfo *clist, *clist1, *clist2, *clist3; CInfo dummy; CSeg *s, *t, *potholes; if( verbose ) fprintf(stderr,"Starting to compute nesting..."); /* * Compute the self-nesting... */ for( i = -1; i < ncontours; i++ ) { if( i == -1 ) clist = &contour0; else clist = &(carray[i]); for( s = clist->closed; s; s = s->next ) { s->nesting = 0; for( t = clist->closed; t; t = t->next ) { if( t == s ) continue; if( is_contained(t, s, s->first) ) s->nesting++; } } } /* * Holes in any given polygon can be either polygons of one level * lower or one level higher. */ dummy.closed = NULL; for( i = -1; i < ncontours; i++ ) { if( i == -1 ) { clist1 = &dummy; clist2 = &contour0; } else if( i == 0 ) { clist1 = &contour0; clist2 = &(carray[i]); } else { clist1 = &(carray[i-1]); clist2 = &(carray[i]); } if( i == ncontours - 1 ) clist3 = &dummy; else clist3 = &(carray[i+1]); /* * Find the holes in the polygons. There can be negative holes (i.e. * depressions) and positive holes (i.e. the next-larger polygons). * This is a bit tricky: for each potential hole, we go through * the list of polygons we could be a hole in, looking for the * most deeply nested one inside of which we reside, then we * go through the holes, and place them on the (neg/pos)_hole list * of the polygon in which they form a hole (keeping only the * least-nested ones. */ for( j = NEG_HOLES; j <= POS_HOLES; j++ ) { if( j == NEG_HOLES ) { potholes = clist1->closed; } else { potholes = clist3->closed; } for( s = potholes; s; s = s->next ) { best = -1; s->tmp_hole = NULL; for( t = clist2->closed; t; t = t->next ) { /* * If I know I'm a hole inside a polygon of nesting level 2, * for instance, then there is no point in checking other * polygons at level 2 or greater... */ if( best >= t->nesting ) continue; if( is_contained(t, s, s->first) ) { best = t->nesting; s->tmp_hole = t; } } } for( s = potholes; s; s = s->next ) { if( s->tmp_hole != NULL ) { /* * Add the temp variable 't' for clarity... * (s->tmp_hole points to the polygon within which s is a hole) */ t = s->tmp_hole; /* * If s is a lower nesting level than the current holes on the * list (or the list is empty) start the list with s, and set * the nesting level; if s is at the same nesting level, push * it onto t's list of holes... */ if( j == NEG_HOLES ) { if( s->nesting < t->neg_depth ) { s->neg_next = NULL; t->neg_holes = s; t->neg_depth = s->nesting; } else if( s->nesting == t->neg_depth ) { s->neg_next = t->neg_holes; t->neg_holes = s; } } else { if( s->nesting < t->pos_depth ) { s->pos_next = NULL; t->pos_holes = s; t->pos_depth = s->nesting; } else if( s->nesting == t->pos_depth ) { s->pos_next = t->pos_holes; t->pos_holes = s; } } } } } } if( verbose ) fprintf(stderr,"done.\n"); return; }