diff --git a/Makefile.am b/Makefile.am index 3e99af8..47d1910 100644 --- a/Makefile.am +++ b/Makefile.am @@ -25,6 +25,7 @@ EXTRA_DIST = makefile.vc CMakeLists.txt autogen.sh \ tests/shp_test.cc \ tests/shpdxf_test.cc \ tests/shpgeo_test.cc \ + tests/shputils_test.cc \ tests/test1.sh tests/test2.sh tests/test3.sh \ tests/expect1.out tests/expect2.out tests/expect3.out \ tests/shape_eg_data/3dpoints.dbf \ diff --git a/contrib/shpgeo.c b/contrib/shpgeo.c index a350d13..ee33762 100644 --- a/contrib/shpgeo.c +++ b/contrib/shpgeo.c @@ -1485,12 +1485,12 @@ SHPObject *SHPUnCompound(const SHPObject *psCShape, int *ringNumber) *ringNumber = -1; else *ringNumber = ring; - /* I am strictly assuming that all R- parts of a complex object - * directly follow their R+, so when we hit a new R+ its a - * new part of a compound object - * a SHPClean may be needed to enforce this as it is not part - * of ESRI's definition of a SHPfile - */ + /* I am strictly assuming that all R- parts of a complex object + * directly follow their R+, so when we hit a new R+ its a + * new part of a compound object + * a SHPClean may be needed to enforce this as it is not part + * of ESRI's definition of a SHPfile + */ #ifdef DEBUG2 printf("(SHPUnCompound) asked for ring %d, lastring is %d \n", lRing, ring); diff --git a/shputils.c b/shputils.c index 830a9ce..378372b 100644 --- a/shputils.c +++ b/shputils.c @@ -32,6 +32,7 @@ #include "shapefil.h" +#include #include #include #include @@ -515,6 +516,182 @@ void check_theme_bnd() puts("WARNING: Theme is outside the clip area."); /** SKIP THEME **/ } +/* -------------------------------------------------------------------- */ +/* Compute where a line segment (x0,y0)-(x1,y1) first crosses */ +/* the axis-aligned clip rectangle. Returns true if an */ +/* intersection was found, storing it in (*xi, *yi). */ +/* -------------------------------------------------------------------- */ +bool compute_clip_intersection(double x0, double y0, double x1, double y1, + double clipxmin, double clipymin, + double clipxmax, double clipymax, double *xi, + double *yi) +{ + const double dx = x1 - x0; + const double dy = y1 - y0; + double tmin = 2.0; /* > 1 means no valid intersection yet */ + double best_xi = 0, best_yi = 0; + + /* Left edge (x = clipxmin) */ + if (dx != 0.0) + { + const double t = (clipxmin - x0) / dx; + if (t > 0.0 && t < 1.0) + { + const double y = y0 + t * dy; + if (y >= clipymin && y <= clipymax && t < tmin) + { + tmin = t; + best_xi = clipxmin; + best_yi = y; + } + } + } + + /* Right edge (x = clipxmax) */ + if (dx != 0.0) + { + const double t = (clipxmax - x0) / dx; + if (t > 0.0 && t < 1.0) + { + const double y = y0 + t * dy; + if (y >= clipymin && y <= clipymax && t < tmin) + { + tmin = t; + best_xi = clipxmax; + best_yi = y; + } + } + } + + /* Bottom edge (y = clipymin) */ + if (dy != 0.0) + { + const double t = (clipymin - y0) / dy; + if (t > 0.0 && t < 1.0) + { + const double x = x0 + t * dx; + if (x >= clipxmin && x <= clipxmax && t < tmin) + { + tmin = t; + best_xi = x; + best_yi = clipymin; + } + } + } + + /* Top edge (y = clipymax) */ + if (dy != 0.0) + { + const double t = (clipymax - y0) / dy; + if (t > 0.0 && t < 1.0) + { + const double x = x0 + t * dx; + if (x >= clipxmin && x <= clipxmax && t < tmin) + { + tmin = t; + best_xi = x; + best_yi = clipymax; + } + } + } + + if (tmin <= 1.0) + { + *xi = best_xi; + *yi = best_yi; + return true; + } + return false; +} + +/* -------------------------------------------------------------------- */ +/* Sutherland-Hodgman polygon clipping against a single edge. */ +/* edge_code: 0=left, 1=right, 2=bottom, 3=top */ +/* Returns the new vertex count. */ +/* -------------------------------------------------------------------- */ +static int sh_clip_edge(const double *inX, const double *inY, const double *inZ, + const double *inM, int nIn, double *outX, double *outY, + double *outZ, double *outM, int edge_code, + double edge_val) +{ + int nOut = 0; + for (int i = 0; i < nIn; i++) + { + const int j = (i + 1) % nIn; + double si, sj; /* signed distance to edge */ + + switch (edge_code) + { + case 0: + si = inX[i] - edge_val; + sj = inX[j] - edge_val; + break; /* left: inside when x >= val */ + case 1: + si = edge_val - inX[i]; + sj = edge_val - inX[j]; + break; /* right: inside when x <= val */ + case 2: + si = inY[i] - edge_val; + sj = inY[j] - edge_val; + break; /* bottom: inside when y >= val */ + case 3: + si = edge_val - inY[i]; + sj = edge_val - inY[j]; + break; /* top: inside when y <= val */ + default: + return 0; + } + + const bool i_in = (si >= 0); + const bool j_in = (sj >= 0); + + if (i_in && j_in) + { + /* Both inside -> output j */ + outX[nOut] = inX[j]; + outY[nOut] = inY[j]; + if (outZ) + outZ[nOut] = inZ[j]; + if (outM) + outM[nOut] = inM[j]; + nOut++; + } + else if (i_in && !j_in) + { + /* Leaving -> output intersection */ + const double t = si / (si - sj); + outX[nOut] = inX[i] + t * (inX[j] - inX[i]); + outY[nOut] = inY[i] + t * (inY[j] - inY[i]); + if (outZ) + outZ[nOut] = inZ[i] + t * (inZ[j] - inZ[i]); + if (outM) + outM[nOut] = inM[i] + t * (inM[j] - inM[i]); + nOut++; + } + else if (!i_in && j_in) + { + /* Entering -> output intersection, then j */ + const double t = si / (si - sj); + outX[nOut] = inX[i] + t * (inX[j] - inX[i]); + outY[nOut] = inY[i] + t * (inY[j] - inY[i]); + if (outZ) + outZ[nOut] = inZ[i] + t * (inZ[j] - inZ[i]); + if (outM) + outM[nOut] = inM[i] + t * (inM[j] - inM[i]); + nOut++; + outX[nOut] = inX[j]; + outY[nOut] = inY[j]; + if (outZ) + outZ[nOut] = inZ[j]; + if (outM) + outM[nOut] = inM[j]; + nOut++; + } + /* else both outside -> nothing */ + } + return nOut; +} + int clip_boundary() { /*** FIRST check the boundary of the feature ***/ @@ -588,55 +765,311 @@ int clip_boundary() if (icut) { /** CUT **/ - /*** Check each vertex in the feature with the Boundary and "CUT" ***/ - /*** THIS CODE WAS NOT COMPLETED! READ NOTE AT THE BOTTOM ***/ + const bool is_polygon = (psCShape->nSHPType == SHPT_POLYGON || + psCShape->nSHPType == SHPT_POLYGONZ || + psCShape->nSHPType == SHPT_POLYGONM); + + if (is_polygon && !ierase) + { + const int nParts = psCShape->nParts > 0 ? psCShape->nParts : 1; + const int totalOrig = psCShape->nVertices; + /* Worst case: each ring gains a few vertices from clipping */ + const int maxTotal = totalOrig * 3 + nParts * 10; + + double *outX = (double *)malloc(maxTotal * sizeof(double)); + double *outY = (double *)malloc(maxTotal * sizeof(double)); + double *outZ = psCShape->padfZ + ? (double *)malloc(maxTotal * sizeof(double)) + : NULL; + double *outM = psCShape->padfM + ? (double *)malloc(maxTotal * sizeof(double)) + : NULL; + int *outParts = (int *)malloc(nParts * sizeof(int)); + int outNParts = 0; + int outTotal = 0; + + for (int p = 0; p < nParts; p++) + { + /* Determine ring range */ + const int ringStart = + psCShape->panPartStart ? psCShape->panPartStart[p] : 0; + const int ringEnd = (p + 1 < nParts && psCShape->panPartStart) + ? psCShape->panPartStart[p + 1] + : totalOrig; + int nRing = ringEnd - ringStart; + + /* Strip closing vertex if present */ + if (nRing >= 4 && + psCShape->padfX[ringStart] == + psCShape->padfX[ringStart + nRing - 1] && + psCShape->padfY[ringStart] == + psCShape->padfY[ringStart + nRing - 1]) + nRing--; + + if (nRing < 3) + continue; + + /* Allocate SH ping-pong buffers for this ring */ + const int maxRing = (nRing + 4) * 2 + 2; + double *bufX[2], *bufY[2], *bufZ[2], *bufM[2]; + for (int b = 0; b < 2; b++) + { + bufX[b] = (double *)malloc(maxRing * sizeof(double)); + bufY[b] = (double *)malloc(maxRing * sizeof(double)); + bufZ[b] = psCShape->padfZ + ? (double *)malloc(maxRing * sizeof(double)) + : NULL; + bufM[b] = psCShape->padfM + ? (double *)malloc(maxRing * sizeof(double)) + : NULL; + } + + memcpy(bufX[0], psCShape->padfX + ringStart, + nRing * sizeof(double)); + memcpy(bufY[0], psCShape->padfY + ringStart, + nRing * sizeof(double)); + if (psCShape->padfZ) + memcpy(bufZ[0], psCShape->padfZ + ringStart, + nRing * sizeof(double)); + if (psCShape->padfM) + memcpy(bufM[0], psCShape->padfM + ringStart, + nRing * sizeof(double)); + + int nCur = nRing; + int src = 0; + + /* Clip against 4 edges: left, right, bottom, top */ + const int edges[] = {0, 1, 2, 3}; + const double vals[] = {cxmin, cxmax, cymin, cymax}; + for (int e = 0; e < 4 && nCur > 0; e++) + { + const int dst = 1 - src; + nCur = + sh_clip_edge(bufX[src], bufY[src], bufZ[src], bufM[src], + nCur, bufX[dst], bufY[dst], bufZ[dst], + bufM[dst], edges[e], vals[e]); + src = dst; + } + + if (nCur >= 3) + { + /* Record part start */ + outParts[outNParts++] = outTotal; + + /* Copy clipped ring + closing vertex */ + memcpy(outX + outTotal, bufX[src], nCur * sizeof(double)); + memcpy(outY + outTotal, bufY[src], nCur * sizeof(double)); + if (outZ) + memcpy(outZ + outTotal, bufZ[src], + nCur * sizeof(double)); + if (outM) + memcpy(outM + outTotal, bufM[src], + nCur * sizeof(double)); + outTotal += nCur; + + /* Close the ring */ + outX[outTotal] = outX[outTotal - nCur]; + outY[outTotal] = outY[outTotal - nCur]; + if (outZ) + outZ[outTotal] = outZ[outTotal - nCur]; + if (outM) + outM[outTotal] = outM[outTotal - nCur]; + outTotal++; + } + + for (int b = 0; b < 2; b++) + { + free(bufX[b]); + free(bufY[b]); + free(bufZ[b]); + free(bufM[b]); + } + } + + if (outNParts > 0) + { + free(psCShape->padfX); + free(psCShape->padfY); + psCShape->padfX = outX; + psCShape->padfY = outY; + if (psCShape->padfZ) + { + free(psCShape->padfZ); + psCShape->padfZ = outZ; + } + if (psCShape->padfM) + { + free(psCShape->padfM); + psCShape->padfM = outM; + } + free(psCShape->panPartStart); + psCShape->panPartStart = outParts; + psCShape->nParts = outNParts; + psCShape->nVertices = outTotal; + return (1); /** WRITE RECORD **/ + } + + free(outX); + free(outY); + free(outZ); + free(outM); + free(outParts); + return (0); /** SKIP RECORD **/ + } + + const int nOrig = psCShape->nVertices; + const int maxVerts = nOrig * 3 + 2; + double *newX = (double *)malloc(maxVerts * sizeof(double)); + double *newY = (double *)malloc(maxVerts * sizeof(double)); + double *newZ = psCShape->padfZ + ? (double *)malloc(maxVerts * sizeof(double)) + : NULL; + double *newM = psCShape->padfM + ? (double *)malloc(maxVerts * sizeof(double)) + : NULL; int i2 = 0; - bool prev_outside = false; - for (int j2 = 0; j2 < psCShape->nVertices; j2++) + + for (int j2 = 0; j2 < nOrig; j2++) { - bool inside = + bool cur_inside = psCShape->padfX[j2] >= cxmin && psCShape->padfX[j2] <= cxmax && psCShape->padfY[j2] >= cymin && psCShape->padfY[j2] <= cymax; - if (ierase) - inside = !inside; - if (inside) + cur_inside = !cur_inside; + + if (j2 > 0) { - if (i2 != j2) + bool prev_inside = psCShape->padfX[j2 - 1] >= cxmin && + psCShape->padfX[j2 - 1] <= cxmax && + psCShape->padfY[j2 - 1] >= cymin && + psCShape->padfY[j2 - 1] <= cymax; + if (ierase) + prev_inside = !prev_inside; + + if (cur_inside != prev_inside) { - if (prev_outside) + double xi, yi; + if (compute_clip_intersection( + psCShape->padfX[j2 - 1], psCShape->padfY[j2 - 1], + psCShape->padfX[j2], psCShape->padfY[j2], cxmin, + cymin, cxmax, cymax, &xi, &yi)) { - /*** AddIntersection(i2); ***/ /*** Add intersection ***/ - prev_outside = false; + newX[i2] = xi; + newY[i2] = yi; + if (newZ || newM) + { + const double dx = + psCShape->padfX[j2] - psCShape->padfX[j2 - 1]; + const double dy = + psCShape->padfY[j2] - psCShape->padfY[j2 - 1]; + const double t = + (fabs(dx) > fabs(dy)) + ? (xi - psCShape->padfX[j2 - 1]) / dx + : (yi - psCShape->padfY[j2 - 1]) / dy; + if (newZ) + newZ[i2] = psCShape->padfZ[j2 - 1] + + t * (psCShape->padfZ[j2] - + psCShape->padfZ[j2 - 1]); + if (newM) + newM[i2] = psCShape->padfM[j2 - 1] + + t * (psCShape->padfM[j2] - + psCShape->padfM[j2 - 1]); + } + i2++; } - psCShape->padfX[i2] = - psCShape->padfX[j2]; /** move vertex **/ - psCShape->padfY[i2] = psCShape->padfY[j2]; } - i2++; - } - else - { - if ((!prev_outside) && (j2 > 0)) + else if (!cur_inside && !prev_inside) { - /*** AddIntersection(i2); ***/ /*** Add intersection (Watch out for j2==i2-1) ***/ - /*** Also a polygon may overlap twice and will split into a several parts ***/ - prev_outside = true; + /* Both outside: check if segment crosses through box */ + double xi1, yi1, xi2, yi2; + if (compute_clip_intersection( + psCShape->padfX[j2 - 1], psCShape->padfY[j2 - 1], + psCShape->padfX[j2], psCShape->padfY[j2], cxmin, + cymin, cxmax, cymax, &xi1, &yi1) && + compute_clip_intersection( + psCShape->padfX[j2], psCShape->padfY[j2], + psCShape->padfX[j2 - 1], psCShape->padfY[j2 - 1], + cxmin, cymin, cxmax, cymax, &xi2, &yi2)) + { + newX[i2] = xi1; + newY[i2] = yi1; + if (newZ || newM) + { + const double dx = + psCShape->padfX[j2] - psCShape->padfX[j2 - 1]; + const double dy = + psCShape->padfY[j2] - psCShape->padfY[j2 - 1]; + const double t = + (fabs(dx) > fabs(dy)) + ? (xi1 - psCShape->padfX[j2 - 1]) / dx + : (yi1 - psCShape->padfY[j2 - 1]) / dy; + if (newZ) + newZ[i2] = psCShape->padfZ[j2 - 1] + + t * (psCShape->padfZ[j2] - + psCShape->padfZ[j2 - 1]); + if (newM) + newM[i2] = psCShape->padfM[j2 - 1] + + t * (psCShape->padfM[j2] - + psCShape->padfM[j2 - 1]); + } + i2++; + newX[i2] = xi2; + newY[i2] = yi2; + if (newZ || newM) + { + const double dx = + psCShape->padfX[j2] - psCShape->padfX[j2 - 1]; + const double dy = + psCShape->padfY[j2] - psCShape->padfY[j2 - 1]; + const double t = + (fabs(dx) > fabs(dy)) + ? (xi2 - psCShape->padfX[j2 - 1]) / dx + : (yi2 - psCShape->padfY[j2 - 1]) / dy; + if (newZ) + newZ[i2] = psCShape->padfZ[j2 - 1] + + t * (psCShape->padfZ[j2] - + psCShape->padfZ[j2 - 1]); + if (newM) + newM[i2] = psCShape->padfM[j2 - 1] + + t * (psCShape->padfM[j2] - + psCShape->padfM[j2 - 1]); + } + i2++; + } } } - } - printf("Vertices:%d OUT:%d Number of Parts:%d\n", - psCShape->nVertices, i2, psCShape->nParts); + if (cur_inside) + { + newX[i2] = psCShape->padfX[j2]; + newY[i2] = psCShape->padfY[j2]; + if (newZ) + newZ[i2] = psCShape->padfZ[j2]; + if (newM) + newM[i2] = psCShape->padfM[j2]; + i2++; + } + } + free(psCShape->padfX); + free(psCShape->padfY); + psCShape->padfX = newX; + psCShape->padfY = newY; + if (psCShape->padfZ) + { + free(psCShape->padfZ); + psCShape->padfZ = newZ; + } + if (psCShape->padfM) + { + free(psCShape->padfM); + psCShape->padfM = newM; + } psCShape->nVertices = i2; if (i2 < 2) return (0); /** SKIP RECORD **/ - /*** (WE ARE NOT CREATING INTERSECTIONS and some lines could be reduced to one point) **/ - // if (i2 == 0) return(0); /** SKIP RECORD **/ - // else return (1); /** WRITE RECORD **/ } /** End CUT **/ @@ -669,6 +1102,88 @@ double findunit(char *unit) return (unitfactor); } +/* -------------------------------------------------------------------- */ +/* Parse comma-separated integer selection values from a string. */ +/* Returns the number of values parsed. */ +/* -------------------------------------------------------------------- */ +long int parse_select_values(const char *input, long int *values, + int max_values) +{ + long int count = 0; + const char *cp = input; + long int val = atol(cp); + while (val > 0 && count < max_values) + { + values[count] = val; + while (*cp >= '0' && *cp <= '9') + cp++; + while (*cp > '\0' && (*cp < '0' || *cp > '9')) + cp++; + val = atol(cp); + count++; + } + return count; +} + +/* -------------------------------------------------------------------- */ +/* Apply factor and shift to the coordinates of a shape object. */ +/* -------------------------------------------------------------------- */ +void transform_coordinates(SHPObject *psShape, double dfFactor, double dfXShift, + double dfYShift) +{ + for (int j = 0; j < psShape->nVertices; j++) + { + psShape->padfX[j] = psShape->padfX[j] * dfFactor + dfXShift; + psShape->padfY[j] = psShape->padfY[j] * dfFactor + dfYShift; + } +} + +/* -------------------------------------------------------------------- */ +/* Copy one DBF record from hDBFin(iRecord) to hDBFout(jRecord) */ +/* using the field mapping in fieldmap[]. */ +/* Returns true on success, false if any write fails. */ +/* -------------------------------------------------------------------- */ +bool copy_dbf_record(DBFHandle hDBFin, int iRecord, DBFHandle hDBFout, + int jRecord, const int *fieldmap, int nFields) +{ + int w, d; + for (int i = 0; i < nFields; i++) + { + if (fieldmap[i] < 0) + continue; + switch (DBFGetFieldInfo(hDBFin, i, NULL, &w, &d)) + { + case FTString: + case FTLogical: + case FTDate: + { + const char *val = DBFReadStringAttribute(hDBFin, iRecord, i); + if (val == NULL) + return false; + if (!DBFWriteStringAttribute(hDBFout, jRecord, fieldmap[i], + val)) + return false; + break; + } + case FTInteger: + if (!DBFWriteIntegerAttribute( + hDBFout, jRecord, fieldmap[i], + DBFReadIntegerAttribute(hDBFin, iRecord, i))) + return false; + break; + case FTDouble: + if (!DBFWriteDoubleAttribute( + hDBFout, jRecord, fieldmap[i], + DBFReadDoubleAttribute(hDBFin, iRecord, i))) + return false; + break; + case FTInvalid: + break; + } + } + return true; +} + /* -------------------------------------------------------------------- */ /* Display a usage message. */ /* -------------------------------------------------------------------- */ @@ -722,75 +1237,18 @@ void error() " { }"); puts( " { }"); - puts(" Note: CUT is not complete and does not create intersections."); - puts(" For more information read programmer comment."); - - /**** Clip functions for Polygon and Cut is not supported - There are several web pages that describe methods of doing this function. - It seem easy to implement until you start writing code. I don't have the - time to add these functions but a did leave a simple cut routine in the - program that can be called by using CUT instead of TOUCH in the - CLIP or ERASE functions. It does not add the intersection of the line and - the clip box, so polygons could look incomplete and lines will come up short. - - Information about clipping lines with a box: - http://www.csclub.uwaterloo.ca/u/mpslager/articles/sutherland/wr.html - Information about finding the intersection of two lines: - http://www.whisqu.se/per/docs/math28.htm - - THE CODE LOOKS LIKE THIS: - ******************************************************** - void Intersect_Lines(float x0,float y0,float x1,float y1, - float x2,float y2,float x3,float y3, - float *xi,float *yi) - { -// this function computes the intersection of the sent lines -// and returns the intersection point, note that the function assumes -// the lines intersect. the function can handle vertical as well -// as horizontal lines. note the function isn't very clever, it simply -// applies the math, but we don't need speed since this is a -// pre-processing step -// The Intersect_lines program came from (http://www.whisqu.se/per/docs/math28.htm) - -float a1,b1,c1, // constants of linear equations -a2,b2,c2, -det_inv, // the inverse of the determinant of the coefficientmatrix -m1,m2; // the slopes of each line - -// compute slopes, note the cludge for infinity, however, this will -// be close enough -if ((x1-x0)!=0) -m1 = (y1-y0)/(x1-x0); -else -m1 = (float)1e+10; // close enough to infinity - - -if ((x3-x2)!=0) -m2 = (y3-y2)/(x3-x2); -else -m2 = (float)1e+10; // close enough to infinity - -// compute constants -a1 = m1; -a2 = m2; -b1 = -1; -b2 = -1; -c1 = (y0-m1*x0); -c2 = (y2-m2*x2); -// compute the inverse of the determinate -det_inv = 1/(a1*b2 - a2*b1); -// use Kramers rule to compute xi and yi -*xi=((b1*c2 - b2*c1)*det_inv); -*yi=((a2*c1 - a1*c2)*det_inv); -} // end Intersect_Lines - **********************************************************/ + puts(" Note: CUT clips lines, polylines, and polygons to the"); + puts(" boundary box (Sutherland-Hodgman algorithm)."); exit(1); } -int main(int argc, char **argv) +/* -------------------------------------------------------------------- */ +/* Parse command-line arguments into global state. */ +/* Calls error()/exit(1) on invalid arguments. */ +/* -------------------------------------------------------------------- */ +void parse_arguments(int argc, char **argv) { - // Check command line usage. if (argc < 2) error(); snprintf(infile, sizeof(infile), "%s", argv[1]); @@ -798,13 +1256,9 @@ int main(int argc, char **argv) { snprintf(outfile, sizeof(outfile), "%s", argv[2]); if (strncasecmp2(outfile, "LIST", 0) == 0) - { ilist = true; - } if (strncasecmp2(outfile, "ALL", 0) == 0) - { iall = true; - } } if (ilist || iall || argc == 2) { @@ -813,7 +1267,6 @@ int main(int argc, char **argv) outfile[0] = '\0'; } - // Look for other functions on the command line. (SELECT, UNIT) for (int i = 3; i < argc; i++) { if ((strncasecmp2(argv[i], "SEL", 3) == 0) || @@ -828,23 +1281,9 @@ int main(int argc, char **argv) i++; if (i >= argc) error(); - selcount = 0; - snprintf(temp, sizeof(temp), "%s", argv[i]); - cpt = temp; - tj = atoi(cpt); - ti = 0; - while (tj > 0) - { - selectvalues[selcount] = tj; - while (*cpt >= '0' && *cpt <= '9') - cpt++; - while (*cpt > '\0' && (*cpt < '0' || *cpt > '9')) - cpt++; - tj = atoi(cpt); - selcount++; - } + selcount = parse_select_values(argv[i], selectvalues, 150); iselect = true; - } /*** End SEL & UNSEL ***/ + } else if ((strncasecmp2(argv[i], "CLIP", 4) == 0) || (strncasecmp2(argv[i], "ERASE", 5) == 0)) { @@ -877,7 +1316,7 @@ int main(int argc, char **argv) cymin, cxmax, cymax); } else - { /*** xmin,ymin,xmax,ymax ***/ + { sscanf(argv[i], "%lf", &cymin); i++; if (i >= argc) @@ -902,7 +1341,7 @@ int main(int argc, char **argv) else error(); iclip = true; - } /*** End CLIP & ERASE ***/ + } else if (strncasecmp2(argv[i], "FACTOR", 0) == 0) { i++; @@ -934,7 +1373,7 @@ int main(int argc, char **argv) } printf("Output file coordinate values will be factored by %lg\n", factor); - } /*** End FACTOR ***/ + } else if (strncasecmp2(argv[i], "SHIFT", 5) == 0) { i++; @@ -947,34 +1386,86 @@ int main(int argc, char **argv) sscanf(argv[i], "%lf", &yshift); iunit = true; printf("X Shift: %lg Y Shift: %lg\n", xshift, yshift); - } /*** End SHIFT ***/ + } else { printf("ERROR: Unknown function %s\n", argv[i]); error(); } } +} + +/* -------------------------------------------------------------------- */ +/* Process all records: select, clip, copy DBF fields, */ +/* transform coordinates, and write shapes. */ +/* -------------------------------------------------------------------- */ +void process_records(void) +{ + int jRecord = DBFGetRecordCount(hDBFappend); + const int nFields = DBFGetFieldCount(hDBF); + + for (int iRecord = 0; iRecord < nEntities; iRecord++) + { + if (iselect && selectrec(iRecord) == 0) + { + psCShape = NULL; + continue; + } + + psCShape = SHPReadObject(hSHP, iRecord); + if (psCShape == NULL) + { + fprintf(stderr, "ERROR: Unable to read shape %d\n", iRecord); + continue; + } + + if (iclip && clip_boundary() == 0) + { + SHPDestroyObject(psCShape); + psCShape = NULL; + continue; + } + + if (!copy_dbf_record(hDBF, iRecord, hDBFappend, jRecord, pt, nFields)) + { + fprintf(stderr, "Warning: Failed to copy DBF record %d\n", iRecord); + } + jRecord++; + + if (iunit) + transform_coordinates(psCShape, factor, xshift, yshift); + + SHPComputeExtents(psCShape); + SHPWriteObject(hSHPappend, -1, psCShape); + + SHPDestroyObject(psCShape); + psCShape = NULL; + } +} - // If there is no data in this file let the user know. - openfiles(); /* Open the infile and the outfile for shape and dbf. */ +#ifndef SHPUTILS_NO_MAIN +int main(int argc, char **argv) +{ + parse_arguments(argc, argv); + + openfiles(); if (DBFGetFieldCount(hDBF) == 0) { puts("There are no fields in this table!"); exit(1); } - // Print out the file bounds. + /* Print out the file bounds. */ { const int iRecord = DBFGetRecordCount(hDBF); SHPGetInfo(hSHP, NULL, NULL, adfBoundsMin, adfBoundsMax); - printf( "Input Bounds: (%lg,%lg) - (%lg,%lg) Entities: %d DBF: %d\n", adfBoundsMin[0], adfBoundsMin[1], adfBoundsMax[0], adfBoundsMax[1], nEntities, iRecord); if (strcmp(outfile, "") == 0) - { /* Describe the shapefile; No other functions */ + { ti = DBFGetFieldCount(hDBF); showitems(); exit(0); @@ -995,116 +1486,25 @@ int main(int argc, char **argv) adfBoundsMin[0], adfBoundsMin[1], adfBoundsMax[0], adfBoundsMax[1], nEntitiesAppend, jRecord); } - /* -------------------------------------------------------------------- */ - /* Find matching fields in the append file or add new items. */ - /* -------------------------------------------------------------------- */ + mergefields(); - /* -------------------------------------------------------------------- */ - /* Find selection field if needed. */ - /* -------------------------------------------------------------------- */ + if (iselect) findselect(); - /* -------------------------------------------------------------------- */ - /* Read all the records */ - /* -------------------------------------------------------------------- */ - int jRecord = DBFGetRecordCount(hDBFappend); - for (int iRecord = 0; iRecord < nEntities; - iRecord++) /** DBFGetRecordCount(hDBF) **/ - { - /* -------------------------------------------------------------------- */ - /* SELECT for values if needed. (Can the record be skipped.) */ - /* -------------------------------------------------------------------- */ - if (iselect) - if (selectrec(iRecord) == 0) - goto SKIP_RECORD; /** SKIP RECORD **/ - - /* -------------------------------------------------------------------- */ - /* Read a Shape record */ - /* -------------------------------------------------------------------- */ - psCShape = SHPReadObject(hSHP, iRecord); - - /* -------------------------------------------------------------------- */ - /* Clip coordinates of shapes if needed. */ - /* -------------------------------------------------------------------- */ - if (iclip) - if (clip_boundary() == 0) - goto SKIP_RECORD; /** SKIP RECORD **/ - - /* -------------------------------------------------------------------- */ - /* Read a DBF record and copy each field. */ - /* -------------------------------------------------------------------- */ - for (int i = 0; i < DBFGetFieldCount(hDBF); i++) - { - /* -------------------------------------------------------------------- */ - /* Store the record according to the type and formatting */ - /* information implicit in the DBF field description. */ - /* -------------------------------------------------------------------- */ - if (pt[i] > -1) /* if the current field exists in output file */ - { - switch (DBFGetFieldInfo(hDBF, i, NULL, &iWidth, &iDecimals)) - { - case FTString: - case FTLogical: - case FTDate: - DBFWriteStringAttribute( - hDBFappend, jRecord, pt[i], - (DBFReadStringAttribute(hDBF, iRecord, i))); - break; - - case FTInteger: - DBFWriteIntegerAttribute( - hDBFappend, jRecord, pt[i], - (DBFReadIntegerAttribute(hDBF, iRecord, i))); - break; - - case FTDouble: - DBFWriteDoubleAttribute( - hDBFappend, jRecord, pt[i], - (DBFReadDoubleAttribute(hDBF, iRecord, i))); - break; - - case FTInvalid: - break; - } - } - } - jRecord++; - /* -------------------------------------------------------------------- */ - /* Change FACTOR and SHIFT coordinates of shapes if needed. */ - /* -------------------------------------------------------------------- */ - if (iunit) - { - for (int j = 0; j < psCShape->nVertices; j++) - { - psCShape->padfX[j] = psCShape->padfX[j] * factor + xshift; - psCShape->padfY[j] = psCShape->padfY[j] * factor + yshift; - } - } - - /* -------------------------------------------------------------------- */ - /* Write the Shape record after recomputing current extents. */ - /* -------------------------------------------------------------------- */ - SHPComputeExtents(psCShape); - SHPWriteObject(hSHPappend, -1, psCShape); + process_records(); - SKIP_RECORD: - SHPDestroyObject(psCShape); - psCShape = NULL; - // j=0; + /* Print out the # of Entities and the file bounds. */ + { + const int jRecord = DBFGetRecordCount(hDBFappend); + SHPGetInfo(hSHPappend, &nEntitiesAppend, &nShapeTypeAppend, + adfBoundsMin, adfBoundsMax); + printf( + "Output Bounds: (%lg,%lg) - (%lg,%lg) Entities: %d DBF: %d\n\n", + adfBoundsMin[0], adfBoundsMin[1], adfBoundsMax[0], adfBoundsMax[1], + nEntitiesAppend, jRecord); } - /* -------------------------------------------------------------------- */ - /* Print out the # of Entities and the file bounds. */ - /* -------------------------------------------------------------------- */ - jRecord = DBFGetRecordCount(hDBFappend); - SHPGetInfo(hSHPappend, &nEntitiesAppend, &nShapeTypeAppend, adfBoundsMin, - adfBoundsMax); - - printf("Output Bounds: (%lg,%lg) - (%lg,%lg) Entities: %d DBF: %d\n\n", - adfBoundsMin[0], adfBoundsMin[1], adfBoundsMax[0], adfBoundsMax[1], - nEntitiesAppend, jRecord); - SHPClose(hSHP); SHPClose(hSHPappend); DBFClose(hDBF); @@ -1123,3 +1523,4 @@ int main(int argc, char **argv) return 0; } +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c706ede..ef9a10a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -108,6 +108,29 @@ if (BUILD_SHARED_LIBS AND ${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.21" AND (WIN ) endif() +add_executable( + shputils_test + ${PROJECT_SOURCE_DIR}/shputils_test.cc + ${CMAKE_SOURCE_DIR}/shputils.c +) +target_include_directories(shputils_test PRIVATE ${CMAKE_SOURCE_DIR}) +target_compile_definitions(shputils_test PRIVATE SHPUTILS_NO_MAIN) +target_link_libraries(shputils_test PRIVATE ${PACKAGE} gtest benchmark) +add_test( + NAME shputils_test + COMMAND shputils_test + WORKING_DIRECTORY "${PROJECT_SOURCE_DIR}" +) +target_compile_features(shputils_test PUBLIC cxx_std_17) +set_target_properties(shputils_test PROPERTIES FOLDER tests CXX_EXTENSIONS OFF) +if (BUILD_SHARED_LIBS AND ${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.21" AND (WIN32 OR CYGWIN)) + add_custom_command( + TARGET shputils_test POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_if_different $ $ + COMMAND_EXPAND_LISTS + ) +endif() + configure_file( ${CMAKE_SOURCE_DIR}/cmake/shapelib.gta.runsettings.in ${CMAKE_BINARY_DIR}/shapelib.gta.runsettings diff --git a/tests/shputils_test.cc b/tests/shputils_test.cc new file mode 100644 index 0000000..fe7e3c2 --- /dev/null +++ b/tests/shputils_test.cc @@ -0,0 +1,1674 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +extern "C" +{ +#include "shapefil.h" + + long int parse_select_values(const char *input, long int *values, + int max_values); + void transform_coordinates(SHPObject *psShape, double dfFactor, + double dfXShift, double dfYShift); + bool copy_dbf_record(DBFHandle hDBFin, int iRecord, DBFHandle hDBFout, + int jRecord, const int *fieldmap, int nFields); + bool compute_clip_intersection(double x0, double y0, double x1, double y1, + double clipxmin, double clipymin, + double clipxmax, double clipymax, double *xi, + double *yi); + void process_records(void); + + /* Global state needed for clip_boundary / check_theme_bnd tests */ + extern SHPObject *psCShape; + extern SHPHandle hSHP, hSHPappend; + extern DBFHandle hDBF, hDBFappend; + extern int *pt; + extern int nEntities, nShapeType; + extern double cxmin, cymin, cxmax, cymax; + extern bool iclip, ierase, itouch, iinside, icut, iselect, iunit; + extern double adfBoundsMin[4], adfBoundsMax[4]; + extern double factor, xshift, yshift; + + int clip_boundary(); + void check_theme_bnd(); +} + +namespace fs = std::filesystem; + +namespace +{ + +// --------------------------------------------------------------------------- +// clip_boundary tests +// --------------------------------------------------------------------------- + +class ClipBoundaryTest : public ::testing::Test +{ + protected: + void SetUp() override + { + // Set default clip box + cxmin = 0.0; + cymin = 0.0; + cxmax = 100.0; + cymax = 100.0; + + // Reset flags + ierase = false; + itouch = false; + iinside = false; + icut = false; + } + + void TearDown() override + { + if (psCShape) + { + SHPDestroyObject(psCShape); + psCShape = nullptr; + } + } + + void MakePoint(double x, double y) + { + psCShape = SHPCreateSimpleObject(SHPT_POINT, 1, &x, &y, nullptr); + } + + void MakeLine(double x1, double y1, double x2, double y2) + { + double xs[] = {x1, x2}; + double ys[] = {y1, y2}; + psCShape = SHPCreateSimpleObject(SHPT_ARC, 2, xs, ys, nullptr); + } + + void MakePolyline(const std::vector> &pts) + { + std::vector xs, ys; + for (const auto &[x, y] : pts) + { + xs.push_back(x); + ys.push_back(y); + } + psCShape = SHPCreateSimpleObject(SHPT_ARC, static_cast(xs.size()), + xs.data(), ys.data(), nullptr); + } + + void MakePolygon(const std::vector> &pts) + { + std::vector xs, ys; + for (const auto &[x, y] : pts) + { + xs.push_back(x); + ys.push_back(y); + } + psCShape = + SHPCreateSimpleObject(SHPT_POLYGON, static_cast(xs.size()), + xs.data(), ys.data(), nullptr); + } + + void + MakePolygonWithHole(const std::vector> &outer, + const std::vector> &hole) + { + std::vector xs, ys; + for (const auto &[x, y] : outer) + { + xs.push_back(x); + ys.push_back(y); + } + for (const auto &[x, y] : hole) + { + xs.push_back(x); + ys.push_back(y); + } + int parts[] = {0, static_cast(outer.size())}; + psCShape = SHPCreateObject(SHPT_POLYGON, -1, 2, parts, nullptr, + static_cast(xs.size()), xs.data(), + ys.data(), nullptr, nullptr); + } +}; + +TEST_F(ClipBoundaryTest, FeatureTotallyOutsideClip) +{ + MakePoint(200.0, 200.0); + EXPECT_EQ(0, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, FeatureTotallyOutsideErase) +{ + ierase = true; + MakePoint(200.0, 200.0); + EXPECT_EQ(1, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, FeatureTotallyInsideClip) +{ + MakePoint(50.0, 50.0); + EXPECT_EQ(1, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, FeatureTotallyInsideErase) +{ + ierase = true; + MakePoint(50.0, 50.0); + EXPECT_EQ(0, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, FeatureCrossesClipInsideMode) +{ + // Feature straddles the clip boundary + iinside = true; + MakeLine(-10.0, 50.0, 50.0, 50.0); + // Inside mode: feature not fully inside -> skip + EXPECT_EQ(0, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, FeatureCrossesEraseInsideMode) +{ + ierase = true; + iinside = true; + MakeLine(-10.0, 50.0, 50.0, 50.0); + // Erase + inside: feature not fully inside clip box -> write + EXPECT_EQ(1, clip_boundary()); +} + +// --------------------------------------------------------------------------- +// compute_clip_intersection tests +// +// Clip box used throughout: (0,0) to (100,100) +// +// 100 +----------+ +// | | +// | clip box | +// | | +// 0 +----------+ +// 0 100 +// --------------------------------------------------------------------------- + +TEST(ComputeClipIntersectionTest, CrossesLeftEdge) +{ + // +----------+ + // B <----X---- A | A=(50,50) B=(-50,50) + // (-50) |0 (50) | X = intersection at (0,50) + // +----------+ + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(50, 50, -50, 50, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(0.0, xi); + EXPECT_DOUBLE_EQ(50.0, yi); +} + +TEST(ComputeClipIntersectionTest, CrossesRightEdge) +{ + // +----------+ + // | A ----X----> B A=(50,50) B=(150,50) + // | (50) |100 (150) X = intersection at (100,50) + // +----------+ + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(50, 50, 150, 50, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(100.0, xi); + EXPECT_DOUBLE_EQ(50.0, yi); +} + +TEST(ComputeClipIntersectionTest, CrossesBottomEdge) +{ + // +----------+ + // | A | A=(50,50) + // | | | + // +----X-----+ X = intersection at (50,0) + // | + // B B=(50,-50) + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(50, 50, 50, -50, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(50.0, xi); + EXPECT_DOUBLE_EQ(0.0, yi); +} + +TEST(ComputeClipIntersectionTest, CrossesTopEdge) +{ + // B B=(50,150) + // | + // +----X-----+ X = intersection at (50,100) + // | | | + // | A | A=(50,50) + // +----------+ + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(50, 50, 50, 150, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(50.0, xi); + EXPECT_DOUBLE_EQ(100.0, yi); +} + +TEST(ComputeClipIntersectionTest, OutsideToInsideCrossesLeft) +{ + double xi, yi; + // From outside (-50,50) to inside (50,50) crosses left edge at (0,50) + EXPECT_TRUE( + compute_clip_intersection(-50, 50, 50, 50, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(0.0, xi); + EXPECT_DOUBLE_EQ(50.0, yi); +} + +TEST(ComputeClipIntersectionTest, DiagonalCross) +{ + // B B=(150,150) + // / + // +-------X X = intersection at corner (100,100) + // | / | + // | A | A=(50,50) + // +-------+ + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(50, 50, 150, 150, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(100.0, xi); + EXPECT_DOUBLE_EQ(100.0, yi); +} + +TEST(ComputeClipIntersectionTest, FullyInsideNoIntersection) +{ + double xi = -1, yi = -1; + // Both endpoints inside -> no intersection + EXPECT_FALSE( + compute_clip_intersection(20, 20, 80, 80, 0, 0, 100, 100, &xi, &yi)); +} + +TEST(ComputeClipIntersectionTest, FullyOutsideNoIntersection) +{ + double xi = -1, yi = -1; + // Both endpoints outside on same side -> no intersection + EXPECT_FALSE(compute_clip_intersection(150, 150, 200, 200, 0, 0, 100, 100, + &xi, &yi)); +} + +TEST(ComputeClipIntersectionTest, TouchesCorner) +{ + // A A=(-50,150) + // \ + // X---------+ X = intersection at corner (0,100) + // | \ | + // | B | B=(50,50) (inside) + // | | + // +---------+ + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(-50, 150, 50, 50, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(0.0, xi); + EXPECT_DOUBLE_EQ(100.0, yi); +} + +TEST(ComputeClipIntersectionTest, GrazesCornerDiagonal) +{ + // B B=(50,150) + // / + // / + // X------+ X = intersection at corner (0,100) + // /| | + // / | | + // A +------+ A=(-100,0) + double xi, yi; + EXPECT_TRUE( + compute_clip_intersection(-100, 0, 50, 150, 0, 0, 100, 100, &xi, &yi)); + EXPECT_DOUBLE_EQ(0.0, xi); + EXPECT_DOUBLE_EQ(100.0, yi); +} + +// --------------------------------------------------------------------------- +// CUT mode tests +// +// Clip box: (0,0) to (100,100) +// CUT trims geometry to the box, computing intersection vertices. +// --------------------------------------------------------------------------- + +TEST_F(ClipBoundaryTest, CutLineExitingRight) +{ + // +----------+ + // | A=======X-------B A=(50,50) B=(150,50) + // | (50) 100 (150) result: A..X + // +----------+ + icut = true; + MakeLine(50.0, 50.0, 150.0, 50.0); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nVertices); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(100.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[1]); +} + +TEST_F(ClipBoundaryTest, CutLineEnteringFromLeft) +{ + // +----------+ + // A----X=======B | A=(-50,50) B=(50,50) + // (-50) 0 (50) | result: X..B + // +----------+ + icut = true; + MakeLine(-50.0, 50.0, 50.0, 50.0); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nVertices); + EXPECT_DOUBLE_EQ(0.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[1]); +} + +TEST_F(ClipBoundaryTest, CutLineCrossingThrough) +{ + // +----------+ + // A-----X1========X2-----B A=(-50,50) B=(150,50) + // (-50) 0 100 (150) both outside, segment crosses through + // +----------+ result: X1..X2 + icut = true; + MakeLine(-50.0, 50.0, 150.0, 50.0); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nVertices); + EXPECT_DOUBLE_EQ(0.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(100.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[1]); +} + +TEST_F(ClipBoundaryTest, CutLineFullyInside) +{ + // +----------+ + // | A====B | A=(20,20) B=(80,80) + // | | entirely inside, no clipping needed + // +----------+ result: A..B unchanged + icut = true; + MakeLine(20.0, 20.0, 80.0, 80.0); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nVertices); + EXPECT_DOUBLE_EQ(20.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(20.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(80.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(80.0, psCShape->padfY[1]); +} + +TEST_F(ClipBoundaryTest, CutLineFullyOutside) +{ + // +----------+ + // | | A=(150,150) B=(200,200) + // | | A----B entirely outside clip box + // +----------+ result: skip record + icut = true; + MakeLine(150.0, 150.0, 200.0, 200.0); + EXPECT_EQ(0, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, CutMultiVertexPolyline) +{ + // +----------+ + // A-----X1===B====X2-----C A=(-10,50) B=(50,50) C=(150,50) + // (-10) 0 (50) 100 (150) + // +----------+ result: X1..B..X2 (3 vertices) + icut = true; + MakePolyline({{-10, 50}, {50, 50}, {150, 50}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(3, psCShape->nVertices); + EXPECT_DOUBLE_EQ(0.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[1]); + EXPECT_DOUBLE_EQ(100.0, psCShape->padfX[2]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[2]); +} + +TEST_F(ClipBoundaryTest, CutEraseLineInside) +{ + // +----------+ + // | A~~~~~~~X=======B A=(50,50) B=(150,50) + // | (erased) 100 (150) erase mode: keep OUTSIDE part + // +----------+ result: X..B + icut = true; + ierase = true; + MakeLine(50.0, 50.0, 150.0, 50.0); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nVertices); + EXPECT_DOUBLE_EQ(100.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[0]); + EXPECT_DOUBLE_EQ(150.0, psCShape->padfX[1]); + EXPECT_DOUBLE_EQ(50.0, psCShape->padfY[1]); +} + +TEST_F(ClipBoundaryTest, CutReducesToOneVertex) +{ + icut = true; + // Single vertex inside, but the line mostly outside -> only 1 vertex + // -> should return 0 (skip) + MakePolyline({{-100, 50}, {1, 50}, {-100, 60}}); + int result = clip_boundary(); + if (psCShape->nVertices < 2) + EXPECT_EQ(0, result); +} + +// --------------------------------------------------------------------------- +// CUT mode polygon tests +// +// Clip box: (0,0) to (100,100) +// Polygon CUT uses Sutherland-Hodgman to clip the polygon to the box. +// All input polygons are closed (first vertex == last vertex). +// --------------------------------------------------------------------------- + +static double RingArea(const double *x, const double *y, int n) +{ + double sum = 0; + for (int i = 0; i < n - 1; i++) + sum += x[i] * y[i + 1] - x[i + 1] * y[i]; + return sum / 2.0; /* signed: >0 CCW, <0 CW */ +} + +static double PolygonArea(const SHPObject *obj) +{ + return fabs(RingArea(obj->padfX, obj->padfY, obj->nVertices)); +} + +static double TotalPolygonArea(const SHPObject *obj) +{ + double total = 0; + for (int p = 0; p < obj->nParts; p++) + { + const int start = obj->panPartStart[p]; + const int end = + (p + 1 < obj->nParts) ? obj->panPartStart[p + 1] : obj->nVertices; + total += + fabs(RingArea(obj->padfX + start, obj->padfY + start, end - start)); + } + return total; +} + +static double NetPolygonArea(const SHPObject *obj) +{ + double total = 0; + for (int p = 0; p < obj->nParts; p++) + { + const int start = obj->panPartStart[p]; + const int end = + (p + 1 < obj->nParts) ? obj->panPartStart[p + 1] : obj->nVertices; + total += RingArea(obj->padfX + start, obj->padfY + start, end - start); + } + return fabs(total); +} + +TEST_F(ClipBoundaryTest, CutPolygonFullyInside) +{ + // 100 +----------+ + // | +------+ | (20,20)→(80,20)→(80,80)→(20,80) + // | | poly | | fully inside → unchanged + // | +------+ | + // 0 +----------+ + // 0 100 + icut = true; + MakePolygon({{20, 20}, {80, 20}, {80, 80}, {20, 80}, {20, 20}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(20.0, psCShape->padfX[0]); + EXPECT_DOUBLE_EQ(20.0, psCShape->padfY[0]); +} + +TEST_F(ClipBoundaryTest, CutPolygonFullyOutside) +{ + // 100 +----------+ poly at (200,200)-(300,300) + // | | entirely outside → skip + // | | + // | | + // 0 +----------+ + // 0 100 + icut = true; + MakePolygon({{200, 200}, {300, 200}, {300, 300}, {200, 300}, {200, 200}}); + EXPECT_EQ(0, clip_boundary()); +} + +TEST_F(ClipBoundaryTest, CutPolygonClipsToRightEdge) +{ + // 100 +----------+ + // 75 | +-----X-----+ (50,25)→(150,25)→(150,75)→(50,75) + // | |clpd | out | X = clip points at x=100 + // 25 | +-----X-----+ result = 50×50 rectangle, area = 2500 + // 0 +----------+ + // 0 50 100 150 + icut = true; + MakePolygon({{50, 25}, {150, 25}, {150, 75}, {50, 75}, {50, 25}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + EXPECT_NEAR(2500.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutTriangleClipsToRightEdge) +{ + // 100 +----------+ + // 80 | *. | (50,20)→(150,50)→(50,80) + // | | `. | apex at (150,50) is outside clip box + // 50 | | `X---* X = clip points at x=100 (y=35 and y=65) + // | | .' | result = trapezoid, area = 2250 + // 20 | *' | + // 0 +----------+ + // 0 50 100 150 + icut = true; + MakePolygon({{50, 20}, {150, 50}, {50, 80}, {50, 20}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + EXPECT_NEAR(2250.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutPolygonCrossesCorner) +{ + // 150 +----------+ + // | | (50,50)→(150,50)→(150,150)→(50,150) + // 100 +----X----* | polygon overlaps top-right of clip box + // | |clpd| | * = clip corner (100,100) + // 50 | +----X-----+ X = clip points (100,50) and (50,100) + // | | result = 50×50 square, area = 2500 + // 0 +----------+ + // 0 50 100 150 + icut = true; + MakePolygon({{50, 50}, {150, 50}, {150, 150}, {50, 150}, {50, 50}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + bool hasCorner = false; + for (int i = 0; i < psCShape->nVertices; i++) + if (psCShape->padfX[i] == 100.0 && psCShape->padfY[i] == 100.0) + hasCorner = true; + EXPECT_TRUE(hasCorner) << "Clipped polygon must include corner (100,100)"; + EXPECT_NEAR(2500.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutPolygonCrossesTwoEdges) +{ + // 100 +----------+ + // 75 +--X-------X--+ (-50,25)→(150,25)→(150,75)→(-50,75) + // 50 | |clipped| | clips at x=0 (left) and x=100 (right) + // 25 +--X-------X--+ area = 100×50 = 5000 + // 0 +----------+ + // -50 0 100 150 + icut = true; + MakePolygon({{-50, 25}, {150, 25}, {150, 75}, {-50, 75}, {-50, 25}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + EXPECT_NEAR(5000.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutPolygonEnclosesBox) +{ + // 150 .--+----------+--. (-50,-50)→(150,-50)→(150,150)→(-50,150) + // 100 | +----------+ | encloses entire clip box + // 50 | | result | | result = clip box + // 0 | +----------+ | area = 100×100 = 10000 + // -50 '--+----------+--' + // -50 0 100 150 + icut = true; + MakePolygon({{-50, -50}, {150, -50}, {150, 150}, {-50, 150}, {-50, -50}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_GE(psCShape->nVertices, 5); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + EXPECT_NEAR(10000.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutDiamondCrossesFourEdges) +{ + // 120 * (50,-20)→(120,50)→(50,120)→(-20,50) + // 100 +----/--+--\---+ crosses all 4 edges + // | / | \ | result: octagon, 9 verts (8 + closing) + // 50 * clip box * area = 10000 - 4×(30²/2) = 8200 + // | \ | / | + // 0 +----\--+--/---+ + // -20 * + // -20 0 50 100 120 + icut = true; + MakePolygon({{50, -20}, {120, 50}, {50, 120}, {-20, 50}, {50, -20}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(9, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + for (int i = 0; i < psCShape->nVertices; i++) + { + EXPECT_GE(psCShape->padfX[i], cxmin - 1e-10); + EXPECT_LE(psCShape->padfX[i], cxmax + 1e-10); + EXPECT_GE(psCShape->padfY[i], cymin - 1e-10); + EXPECT_LE(psCShape->padfY[i], cymax + 1e-10); + } + EXPECT_NEAR(8200.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutConcaveNotchCrossesRightEdge) +{ + // 100 +----------+ + // 80 | +-----X--+ (50,20)→(120,20)→(120,40)→(70,40)→ + // 60 | | +--X--+ (70,60)→(120,60)→(120,80)→(50,80) + // 50 | | |ntch | notch exits/re-enters at x=100 + // 40 | | +--X--+ C-shape, 9 verts, area = 2400 + // 20 | +-----X--+ + // 0 +----------+ + // 0 50 70 100 120 + icut = true; + MakePolygon({{50, 20}, + {120, 20}, + {120, 40}, + {70, 40}, + {70, 60}, + {120, 60}, + {120, 80}, + {50, 80}, + {50, 20}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(9, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + for (int i = 0; i < psCShape->nVertices; i++) + { + EXPECT_GE(psCShape->padfX[i], cxmin - 1e-10); + EXPECT_LE(psCShape->padfX[i], cxmax + 1e-10); + } + EXPECT_NEAR(2400.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutChevronCrossesBottomEdge) +{ + // 100 +----------+ + // 80 | A--------B | A=(20,80) B=(80,80) + // 50 | \ / | chevron dips below y=0 + // | \ / | + // 0 +----X1-X2---+ X1≈(44,0) X2≈(56,0) + // -20 * C=(50,-20) outside + // 0 20 44 56 80 100 trapezoid, area = 2880 + icut = true; + MakePolygon({{20, 80}, {50, -20}, {80, 80}, {20, 80}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(5, psCShape->nVertices); + EXPECT_DOUBLE_EQ(psCShape->padfX[0], + psCShape->padfX[psCShape->nVertices - 1]); + EXPECT_DOUBLE_EQ(psCShape->padfY[0], + psCShape->padfY[psCShape->nVertices - 1]); + for (int i = 0; i < psCShape->nVertices; i++) + { + EXPECT_GE(psCShape->padfY[i], cymin - 1e-10); + EXPECT_LE(psCShape->padfY[i], cymax + 1e-10); + } + EXPECT_NEAR(2880.0, PolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutPolygonWithHoleFullyInside) +{ + // 100 +----------+ + // 90 | +------+ | outer: (10,10)→(90,10)→(90,90)→(10,90) + // 70 | | +--+ | | hole: (30,30)→(30,70)→(70,70)→(70,30) CW + // 50 | | | | | | both fully inside clip box + // 30 | | +--+ | | net area = 80×80 − 40×40 = 4800 + // 10 | +------+ | + // 0 +----------+ + // 0 10 30 70 90 100 + icut = true; + MakePolygonWithHole({{10, 10}, {90, 10}, {90, 90}, {10, 90}, {10, 10}}, + {{30, 30}, {30, 70}, {70, 70}, {70, 30}, {30, 30}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nParts); + EXPECT_EQ(10, psCShape->nVertices); // 5 outer + 5 hole + EXPECT_NEAR(4800.0, NetPolygonArea(psCShape), 1.0); +} + +TEST_F(ClipBoundaryTest, CutPolygonWithHoleCrossesTop) +{ + // 130 +------+ + // 110 | +--+ | outer: (10,10)→(90,10)→(90,130)→(10,130) + // 100 +-X-+--+-X-+ hole: (30,70)→(30,110)→(70,110)→(70,70) CW + // | | | | | | X = both rings clip at y=100 + // 70 | | +--+ | | outer clipped → 80×90 = 7200 + // | | | | hole clipped → 40×30 = 1200 + // 10 | +------+ | net area = 7200 − 1200 = 6000 + // 0 +----------+ + // 0 10 30 70 90 100 + icut = true; + MakePolygonWithHole({{10, 10}, {90, 10}, {90, 130}, {10, 130}, {10, 10}}, + {{30, 70}, {30, 110}, {70, 110}, {70, 70}, {30, 70}}); + EXPECT_EQ(1, clip_boundary()); + EXPECT_EQ(2, psCShape->nParts); + EXPECT_EQ(10, psCShape->nVertices); // 5 outer + 5 hole + // Each ring must be closed + for (int p = 0; p < psCShape->nParts; p++) + { + const int start = psCShape->panPartStart[p]; + const int end = (p + 1 < psCShape->nParts) + ? psCShape->panPartStart[p + 1] + : psCShape->nVertices; + EXPECT_DOUBLE_EQ(psCShape->padfX[start], psCShape->padfX[end - 1]) + << "ring " << p; + EXPECT_DOUBLE_EQ(psCShape->padfY[start], psCShape->padfY[end - 1]) + << "ring " << p; + } + EXPECT_NEAR(6000.0, NetPolygonArea(psCShape), 1.0); +} + +// --------------------------------------------------------------------------- +// check_theme_bnd tests +// --------------------------------------------------------------------------- + +class CheckThemeBndTest : public ::testing::Test +{ + protected: + void SetUp() override + { + cxmin = 0.0; + cymin = 0.0; + cxmax = 100.0; + cymax = 100.0; + ierase = false; + iclip = true; + nEntities = 10; + } +}; + +TEST_F(CheckThemeBndTest, ThemeTotallyInsideClip) +{ + adfBoundsMin[0] = 10.0; + adfBoundsMin[1] = 10.0; + adfBoundsMax[0] = 90.0; + adfBoundsMax[1] = 90.0; + check_theme_bnd(); + // Theme inside clip -> clip not needed + EXPECT_FALSE(iclip); + EXPECT_EQ(10, nEntities); +} + +TEST_F(CheckThemeBndTest, ThemeTotallyInsideErase) +{ + ierase = true; + adfBoundsMin[0] = 10.0; + adfBoundsMin[1] = 10.0; + adfBoundsMax[0] = 90.0; + adfBoundsMax[1] = 90.0; + check_theme_bnd(); + // Theme inside erase -> skip all + EXPECT_EQ(0, nEntities); +} + +TEST_F(CheckThemeBndTest, ThemeTotallyOutsideClip) +{ + adfBoundsMin[0] = 200.0; + adfBoundsMin[1] = 200.0; + adfBoundsMax[0] = 300.0; + adfBoundsMax[1] = 300.0; + check_theme_bnd(); + // Theme outside clip -> skip all + EXPECT_EQ(0, nEntities); +} + +TEST_F(CheckThemeBndTest, ThemeTotallyOutsideErase) +{ + ierase = true; + adfBoundsMin[0] = 200.0; + adfBoundsMin[1] = 200.0; + adfBoundsMax[0] = 300.0; + adfBoundsMax[1] = 300.0; + check_theme_bnd(); + // Theme outside erase -> clip not needed, write theme + EXPECT_FALSE(iclip); +} + +// --------------------------------------------------------------------------- +// parse_select_values tests +// --------------------------------------------------------------------------- + +TEST(ParseSelectValuesTest, SingleValue) +{ + long int values[150] = {}; + const long int count = parse_select_values("42", values, 150); + EXPECT_EQ(1, count); + EXPECT_EQ(42, values[0]); +} + +TEST(ParseSelectValuesTest, MultipleCommaSeparated) +{ + long int values[150] = {}; + const long int count = parse_select_values("3,5,9,13,17,27", values, 150); + EXPECT_EQ(6, count); + EXPECT_EQ(3, values[0]); + EXPECT_EQ(5, values[1]); + EXPECT_EQ(9, values[2]); + EXPECT_EQ(13, values[3]); + EXPECT_EQ(17, values[4]); + EXPECT_EQ(27, values[5]); +} + +TEST(ParseSelectValuesTest, SpaceSeparated) +{ + long int values[150] = {}; + const long int count = parse_select_values("10 20 30", values, 150); + EXPECT_EQ(3, count); + EXPECT_EQ(10, values[0]); + EXPECT_EQ(20, values[1]); + EXPECT_EQ(30, values[2]); +} + +TEST(ParseSelectValuesTest, EmptyString) +{ + long int values[150] = {}; + const long int count = parse_select_values("", values, 150); + EXPECT_EQ(0, count); +} + +TEST(ParseSelectValuesTest, NoDigits) +{ + long int values[150] = {}; + const long int count = parse_select_values("abc", values, 150); + EXPECT_EQ(0, count); +} + +TEST(ParseSelectValuesTest, MaxValuesLimit) +{ + long int values[3] = {}; + const long int count = parse_select_values("1,2,3,4,5", values, 3); + EXPECT_EQ(3, count); + EXPECT_EQ(1, values[0]); + EXPECT_EQ(2, values[1]); + EXPECT_EQ(3, values[2]); +} + +TEST(ParseSelectValuesTest, MixedDelimiters) +{ + long int values[150] = {}; + const long int count = parse_select_values("100, 200; 300", values, 150); + EXPECT_EQ(3, count); + EXPECT_EQ(100, values[0]); + EXPECT_EQ(200, values[1]); + EXPECT_EQ(300, values[2]); +} + +// --------------------------------------------------------------------------- +// transform_coordinates tests +// --------------------------------------------------------------------------- + +class TransformCoordinatesTest : public ::testing::Test +{ + protected: + void TearDown() override + { + if (shape_) + { + SHPDestroyObject(shape_); + shape_ = nullptr; + } + } + + SHPObject *shape_ = nullptr; +}; + +TEST_F(TransformCoordinatesTest, IdentityTransform) +{ + double x[] = {1.0, 2.0, 3.0}; + double y[] = {4.0, 5.0, 6.0}; + shape_ = SHPCreateSimpleObject(SHPT_ARC, 3, x, y, nullptr); + transform_coordinates(shape_, 1.0, 0.0, 0.0); + EXPECT_DOUBLE_EQ(1.0, shape_->padfX[0]); + EXPECT_DOUBLE_EQ(2.0, shape_->padfX[1]); + EXPECT_DOUBLE_EQ(3.0, shape_->padfX[2]); + EXPECT_DOUBLE_EQ(4.0, shape_->padfY[0]); + EXPECT_DOUBLE_EQ(5.0, shape_->padfY[1]); + EXPECT_DOUBLE_EQ(6.0, shape_->padfY[2]); +} + +TEST_F(TransformCoordinatesTest, ScaleOnly) +{ + double x[] = {10.0, 20.0}; + double y[] = {30.0, 40.0}; + shape_ = SHPCreateSimpleObject(SHPT_ARC, 2, x, y, nullptr); + transform_coordinates(shape_, 2.0, 0.0, 0.0); + EXPECT_DOUBLE_EQ(20.0, shape_->padfX[0]); + EXPECT_DOUBLE_EQ(40.0, shape_->padfX[1]); + EXPECT_DOUBLE_EQ(60.0, shape_->padfY[0]); + EXPECT_DOUBLE_EQ(80.0, shape_->padfY[1]); +} + +TEST_F(TransformCoordinatesTest, ShiftOnly) +{ + double x[] = {1.0}; + double y[] = {2.0}; + shape_ = SHPCreateSimpleObject(SHPT_POINT, 1, x, y, nullptr); + transform_coordinates(shape_, 1.0, 100.0, 200.0); + EXPECT_DOUBLE_EQ(101.0, shape_->padfX[0]); + EXPECT_DOUBLE_EQ(202.0, shape_->padfY[0]); +} + +TEST_F(TransformCoordinatesTest, ScaleAndShift) +{ + double x[] = {5.0, 10.0}; + double y[] = {15.0, 20.0}; + shape_ = SHPCreateSimpleObject(SHPT_ARC, 2, x, y, nullptr); + // factor = 3.0, xshift = 1.0, yshift = 2.0 + // new_x = x * 3.0 + 1.0, new_y = y * 3.0 + 2.0 + transform_coordinates(shape_, 3.0, 1.0, 2.0); + EXPECT_DOUBLE_EQ(16.0, shape_->padfX[0]); + EXPECT_DOUBLE_EQ(31.0, shape_->padfX[1]); + EXPECT_DOUBLE_EQ(47.0, shape_->padfY[0]); + EXPECT_DOUBLE_EQ(62.0, shape_->padfY[1]); +} + +TEST_F(TransformCoordinatesTest, FeetToMetersConversion) +{ + // 1200 (feet unit) / 3937 (meter unit) ≈ 0.3048 + const double feet_to_meters = 1200.0 / 3937.0; + double x[] = {100.0}; + double y[] = {200.0}; + shape_ = SHPCreateSimpleObject(SHPT_POINT, 1, x, y, nullptr); + transform_coordinates(shape_, feet_to_meters, 0.0, 0.0); + EXPECT_NEAR(30.48, shape_->padfX[0], 0.001); + EXPECT_NEAR(60.96, shape_->padfY[0], 0.001); +} + +// --------------------------------------------------------------------------- +// copy_dbf_record tests +// --------------------------------------------------------------------------- + +static std::string GenerateUniqueName(std::string_view stem) +{ + const auto now = std::chrono::system_clock::now(); + const auto timestamp = + std::chrono::duration_cast( + now.time_since_epoch()) + .count(); + std::ostringstream oss; + oss << stem << '_' << timestamp; + return oss.str(); +} + +class CopyDbfRecordTest : public ::testing::Test +{ + protected: + void SetUp() override + { + dir_ = fs::current_path() / GenerateUniqueName("copydbf"); + fs::create_directories(dir_); + } + + void TearDown() override + { + if (hIn_) + DBFClose(hIn_); + if (hOut_) + DBFClose(hOut_); + fs::remove_all(dir_); + } + + fs::path dir_; + DBFHandle hIn_ = nullptr; + DBFHandle hOut_ = nullptr; +}; + +TEST_F(CopyDbfRecordTest, CopiesIntegerField) +{ + const auto srcPath = (dir_ / "src.dbf").string(); + hIn_ = DBFCreate(srcPath.c_str()); + ASSERT_NE(nullptr, hIn_); + EXPECT_NE(-1, DBFAddField(hIn_, "NUM", FTInteger, 10, 0)); + EXPECT_TRUE(DBFWriteIntegerAttribute(hIn_, 0, 0, 42)); + DBFClose(hIn_); + hIn_ = DBFOpen(srcPath.c_str(), "rb"); + ASSERT_NE(nullptr, hIn_); + + const auto dstPath = (dir_ / "dst.dbf").string(); + hOut_ = DBFCreate(dstPath.c_str()); + ASSERT_NE(nullptr, hOut_); + EXPECT_NE(-1, DBFAddField(hOut_, "NUM", FTInteger, 10, 0)); + + int fieldmap[] = {0}; + EXPECT_TRUE(copy_dbf_record(hIn_, 0, hOut_, 0, fieldmap, 1)); + + EXPECT_EQ(42, DBFReadIntegerAttribute(hOut_, 0, 0)); +} + +TEST_F(CopyDbfRecordTest, CopiesStringField) +{ + const auto srcPath = (dir_ / "src.dbf").string(); + hIn_ = DBFCreate(srcPath.c_str()); + ASSERT_NE(nullptr, hIn_); + EXPECT_NE(-1, DBFAddField(hIn_, "NAME", FTString, 20, 0)); + EXPECT_TRUE(DBFWriteStringAttribute(hIn_, 0, 0, "hello")); + DBFClose(hIn_); + hIn_ = DBFOpen(srcPath.c_str(), "rb"); + ASSERT_NE(nullptr, hIn_); + + const auto dstPath = (dir_ / "dst.dbf").string(); + hOut_ = DBFCreate(dstPath.c_str()); + ASSERT_NE(nullptr, hOut_); + EXPECT_NE(-1, DBFAddField(hOut_, "NAME", FTString, 20, 0)); + + int fieldmap[] = {0}; + EXPECT_TRUE(copy_dbf_record(hIn_, 0, hOut_, 0, fieldmap, 1)); + + const char *raw = DBFReadStringAttribute(hOut_, 0, 0); + ASSERT_NE(nullptr, raw); + const std::string result = raw; + const auto trimmed = result.substr(0, result.find_last_not_of(' ') + 1); + EXPECT_EQ("hello", trimmed); +} + +TEST_F(CopyDbfRecordTest, CopiesDoubleField) +{ + const auto srcPath = (dir_ / "src.dbf").string(); + hIn_ = DBFCreate(srcPath.c_str()); + ASSERT_NE(nullptr, hIn_); + EXPECT_NE(-1, DBFAddField(hIn_, "VAL", FTDouble, 18, 6)); + EXPECT_TRUE(DBFWriteDoubleAttribute(hIn_, 0, 0, 3.14)); + DBFClose(hIn_); + hIn_ = DBFOpen(srcPath.c_str(), "rb"); + ASSERT_NE(nullptr, hIn_); + + const auto dstPath = (dir_ / "dst.dbf").string(); + hOut_ = DBFCreate(dstPath.c_str()); + ASSERT_NE(nullptr, hOut_); + EXPECT_NE(-1, DBFAddField(hOut_, "VAL", FTDouble, 18, 6)); + + int fieldmap[] = {0}; + EXPECT_TRUE(copy_dbf_record(hIn_, 0, hOut_, 0, fieldmap, 1)); + + EXPECT_DOUBLE_EQ(3.14, DBFReadDoubleAttribute(hOut_, 0, 0)); +} + +TEST_F(CopyDbfRecordTest, SkipsUnmappedField) +{ + const auto srcPath = (dir_ / "src.dbf").string(); + hIn_ = DBFCreate(srcPath.c_str()); + ASSERT_NE(nullptr, hIn_); + EXPECT_NE(-1, DBFAddField(hIn_, "A", FTInteger, 10, 0)); + EXPECT_NE(-1, DBFAddField(hIn_, "B", FTInteger, 10, 0)); + EXPECT_TRUE(DBFWriteIntegerAttribute(hIn_, 0, 0, 10)); + EXPECT_TRUE(DBFWriteIntegerAttribute(hIn_, 0, 1, 20)); + DBFClose(hIn_); + hIn_ = DBFOpen(srcPath.c_str(), "rb"); + ASSERT_NE(nullptr, hIn_); + + const auto dstPath = (dir_ / "dst.dbf").string(); + hOut_ = DBFCreate(dstPath.c_str()); + ASSERT_NE(nullptr, hOut_); + EXPECT_NE(-1, DBFAddField(hOut_, "B", FTInteger, 10, 0)); + + int fieldmap[] = {-1, 0}; + EXPECT_TRUE(copy_dbf_record(hIn_, 0, hOut_, 0, fieldmap, 2)); + + EXPECT_EQ(20, DBFReadIntegerAttribute(hOut_, 0, 0)); +} + +TEST_F(CopyDbfRecordTest, RemapsFieldOrder) +{ + const auto srcPath = (dir_ / "src.dbf").string(); + hIn_ = DBFCreate(srcPath.c_str()); + ASSERT_NE(nullptr, hIn_); + EXPECT_NE(-1, DBFAddField(hIn_, "X", FTInteger, 10, 0)); + EXPECT_NE(-1, DBFAddField(hIn_, "Y", FTInteger, 10, 0)); + EXPECT_TRUE(DBFWriteIntegerAttribute(hIn_, 0, 0, 100)); + EXPECT_TRUE(DBFWriteIntegerAttribute(hIn_, 0, 1, 200)); + DBFClose(hIn_); + hIn_ = DBFOpen(srcPath.c_str(), "rb"); + ASSERT_NE(nullptr, hIn_); + + const auto dstPath = (dir_ / "dst.dbf").string(); + hOut_ = DBFCreate(dstPath.c_str()); + ASSERT_NE(nullptr, hOut_); + EXPECT_NE(-1, DBFAddField(hOut_, "Y", FTInteger, 10, 0)); + EXPECT_NE(-1, DBFAddField(hOut_, "X", FTInteger, 10, 0)); + + int fieldmap[] = {1, 0}; + EXPECT_TRUE(copy_dbf_record(hIn_, 0, hOut_, 0, fieldmap, 2)); + + EXPECT_EQ(200, DBFReadIntegerAttribute(hOut_, 0, 0)); // Y + EXPECT_EQ(100, DBFReadIntegerAttribute(hOut_, 0, 1)); // X +} + +// --------------------------------------------------------------------------- +// process_records integration tests +// --------------------------------------------------------------------------- + +class ProcessRecordsTest : public ::testing::Test +{ + protected: + void SetUp() override + { + dir_ = fs::current_path() / GenerateUniqueName("procrec"); + fs::create_directories(dir_); + + // Reset global state + iclip = false; + ierase = false; + itouch = false; + iinside = false; + icut = false; + iselect = false; + iunit = false; + factor = 0; + xshift = 0; + yshift = 0; + psCShape = nullptr; + pt = nullptr; + } + + void TearDown() override + { + if (hSHP) + { + SHPClose(hSHP); + hSHP = nullptr; + } + if (hSHPappend) + { + SHPClose(hSHPappend); + hSHPappend = nullptr; + } + if (hDBF) + { + DBFClose(hDBF); + hDBF = nullptr; + } + if (hDBFappend) + { + DBFClose(hDBFappend); + hDBFappend = nullptr; + } + free(pt); + pt = nullptr; + fs::remove_all(dir_); + } + + // Create a simple point shapefile with nPts points and an integer field + void CreateTestData(int nPts) + { + const auto srcShp = (dir_ / "src.shp").string(); + hSHP = SHPCreate(srcShp.c_str(), SHPT_POINT); + ASSERT_NE(nullptr, hSHP); + for (int i = 0; i < nPts; i++) + { + double x = static_cast(i * 10); + double y = static_cast(i * 20); + auto *obj = SHPCreateSimpleObject(SHPT_POINT, 1, &x, &y, nullptr); + SHPWriteObject(hSHP, -1, obj); + SHPDestroyObject(obj); + } + SHPClose(hSHP); + + const auto srcDbf = (dir_ / "src.dbf").string(); + hDBF = DBFCreate(srcDbf.c_str()); + ASSERT_NE(nullptr, hDBF); + EXPECT_NE(-1, DBFAddField(hDBF, "VAL", FTInteger, 10, 0)); + for (int i = 0; i < nPts; i++) + EXPECT_TRUE(DBFWriteIntegerAttribute(hDBF, i, 0, (i + 1) * 100)); + DBFClose(hDBF); + + // Reopen as read-only + hSHP = SHPOpen(srcShp.c_str(), "rb"); + ASSERT_NE(nullptr, hSHP); + SHPGetInfo(hSHP, &nEntities, &nShapeType, nullptr, nullptr); + + hDBF = DBFOpen(srcDbf.c_str(), "rb"); + ASSERT_NE(nullptr, hDBF); + + // Create output files + const auto dstShp = (dir_ / "dst.shp").string(); + hSHPappend = SHPCreate(dstShp.c_str(), SHPT_POINT); + ASSERT_NE(nullptr, hSHPappend); + + const auto dstDbf = (dir_ / "dst.dbf").string(); + hDBFappend = DBFCreate(dstDbf.c_str()); + ASSERT_NE(nullptr, hDBFappend); + EXPECT_NE(-1, DBFAddField(hDBFappend, "VAL", FTInteger, 10, 0)); + + // Identity field map + pt = static_cast(malloc(sizeof(int))); + pt[0] = 0; + } + + fs::path dir_; +}; + +TEST_F(ProcessRecordsTest, CopiesAllRecords) +{ + CreateTestData(3); + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_EQ(3, outEntities); + EXPECT_EQ(3, DBFGetRecordCount(hDBFappend)); + EXPECT_EQ(100, DBFReadIntegerAttribute(hDBFappend, 0, 0)); + EXPECT_EQ(200, DBFReadIntegerAttribute(hDBFappend, 1, 0)); + EXPECT_EQ(300, DBFReadIntegerAttribute(hDBFappend, 2, 0)); +} + +TEST_F(ProcessRecordsTest, AppliesTransform) +{ + CreateTestData(1); + iunit = true; + factor = 2.0; + xshift = 5.0; + yshift = 10.0; + + process_records(); + + // Original point: (0, 0), transformed: (0*2+5, 0*2+10) = (5, 10) + auto *obj = SHPReadObject(hSHPappend, 0); + ASSERT_NE(nullptr, obj); + EXPECT_DOUBLE_EQ(5.0, obj->padfX[0]); + EXPECT_DOUBLE_EQ(10.0, obj->padfY[0]); + SHPDestroyObject(obj); +} + +TEST_F(ProcessRecordsTest, ClipFiltersOutside) +{ + CreateTestData(3); + // Points at (0,0), (10,20), (20,40) + // Clip box includes only (10,20) and (20,40) + iclip = true; + cxmin = 5.0; + cymin = 5.0; + cxmax = 25.0; + cymax = 45.0; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_EQ(2, outEntities); + EXPECT_EQ(2, DBFGetRecordCount(hDBFappend)); +} + +TEST_F(ProcessRecordsTest, EmptyInput) +{ + CreateTestData(0); + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_EQ(0, outEntities); + EXPECT_EQ(0, DBFGetRecordCount(hDBFappend)); +} + +// --------------------------------------------------------------------------- +// High-level CUT tests using real shapefiles from shape_eg_data/ +// +// pline.shp: Arc 460 shapes bounds (1296367..1302699, 228199..237185) +// brklinz.shp: ArcZ 122 shapes bounds (6294338..6296322, 1978444..1979694) +// csah.shp: Arc 124 shapes bounds (257105..335498, 5176099..5226768) +// polygon.shp: Polygon 474 shapes bounds (471127..489292, 4751545..4765611) +// --------------------------------------------------------------------------- + +static const auto kTestData = fs::path{"shape_eg_data"}; + +class CutRealDataTest : public ::testing::Test +{ + protected: + void SetUp() override + { + dir_ = fs::current_path() / GenerateUniqueName("cutreal"); + fs::create_directories(dir_); + + iclip = true; + ierase = false; + itouch = false; + iinside = false; + icut = true; + iselect = false; + iunit = false; + factor = 0; + xshift = 0; + yshift = 0; + psCShape = nullptr; + pt = nullptr; + } + + void TearDown() override + { + if (hSHP) + { + SHPClose(hSHP); + hSHP = nullptr; + } + if (hSHPappend) + { + SHPClose(hSHPappend); + hSHPappend = nullptr; + } + if (hDBF) + { + DBFClose(hDBF); + hDBF = nullptr; + } + if (hDBFappend) + { + DBFClose(hDBFappend); + hDBFappend = nullptr; + } + free(pt); + pt = nullptr; + fs::remove_all(dir_); + } + + // Open a real shapefile as input, create empty output, build field map + void OpenRealData(const fs::path &basename) + { + const auto shpPath = (kTestData / basename).string(); + auto dbfBase = basename; + dbfBase.replace_extension(".dbf"); + const auto dbfPath = (kTestData / dbfBase).string(); + + hSHP = SHPOpen(shpPath.c_str(), "rb"); + ASSERT_NE(nullptr, hSHP) << "Cannot open " << shpPath; + SHPGetInfo(hSHP, &nEntities, &nShapeType, nullptr, nullptr); + + hDBF = DBFOpen(dbfPath.c_str(), "rb"); + ASSERT_NE(nullptr, hDBF) << "Cannot open " << dbfPath; + + const int nFields = DBFGetFieldCount(hDBF); + + // Create output files + const auto dstShp = (dir_ / "out.shp").string(); + hSHPappend = SHPCreate(dstShp.c_str(), nShapeType); + ASSERT_NE(nullptr, hSHPappend); + + const auto dstDbf = (dir_ / "out.dbf").string(); + hDBFappend = DBFCreate(dstDbf.c_str()); + ASSERT_NE(nullptr, hDBFappend); + + // Copy field definitions and build identity field map + pt = static_cast(malloc(nFields * sizeof(int))); + for (int i = 0; i < nFields; i++) + { + char name[12]; + int width, decimals; + const auto type = DBFGetFieldInfo(hDBF, i, name, &width, &decimals); + EXPECT_NE(-1, DBFAddField(hDBFappend, name, type, width, decimals)); + pt[i] = i; + } + } + + fs::path dir_; +}; + +TEST_F(CutRealDataTest, CutPlineClipsToSubset) +{ + // pline.shp bounds: + // X: 1296367 .. 1302699 Y: 228199 .. 237185 + // + // 237185 +========================+ + // | pline.shp | + // | +----------+ | + // | | clip box | | + // | +----------+ | + // 228199 +========================+ + // 1296367 1302699 + // + // Clip box covers roughly the center third: + OpenRealData("pline.shp"); + const int totalIn = nEntities; + EXPECT_EQ(460, totalIn); + + cxmin = 1298000; + cymin = 231000; + cxmax = 1301000; + cymax = 235000; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, 0); // some shapes survived + EXPECT_LT(outEntities, totalIn); // but not all + + // Every output shape must have valid vertices inside or on the clip box + for (int i = 0; i < outEntities; i++) + { + auto *obj = SHPReadObject(hSHPappend, i); + ASSERT_NE(nullptr, obj); + EXPECT_GE(obj->nVertices, 2) << "CUT polyline must have >= 2 vertices"; + for (int v = 0; v < obj->nVertices; v++) + { + EXPECT_GE(obj->padfX[v], cxmin - 1e-6) + << "shape " << i << " vertex " << v << " X below clip min"; + EXPECT_LE(obj->padfX[v], cxmax + 1e-6) + << "shape " << i << " vertex " << v << " X above clip max"; + EXPECT_GE(obj->padfY[v], cymin - 1e-6) + << "shape " << i << " vertex " << v << " Y below clip min"; + EXPECT_LE(obj->padfY[v], cymax + 1e-6) + << "shape " << i << " vertex " << v << " Y above clip max"; + } + SHPDestroyObject(obj); + } + + EXPECT_EQ(outEntities, DBFGetRecordCount(hDBFappend)); +} + +TEST_F(CutRealDataTest, CutBrklinzPreservesZ) +{ + // brklinz.shp: ArcZ with 122 shapes, has Z coordinates + // Clip a sub-region and verify Z values are preserved (non-zero) + OpenRealData("brklinz.shp"); + const int totalIn = nEntities; + EXPECT_EQ(122, totalIn); + + cxmin = 6294500; + cymin = 1978500; + cxmax = 6296000; + cymax = 1979500; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, 0); + + // Check that at least one output shape has non-zero Z + bool foundZ = false; + for (int i = 0; i < outEntities && !foundZ; i++) + { + auto *obj = SHPReadObject(hSHPappend, i); + ASSERT_NE(nullptr, obj); + if (obj->padfZ != nullptr) + { + for (int v = 0; v < obj->nVertices; v++) + { + if (obj->padfZ[v] != 0.0) + { + foundZ = true; + break; + } + } + } + SHPDestroyObject(obj); + } + EXPECT_TRUE(foundZ) << "Expected non-zero Z in ArcZ output"; +} + +TEST_F(CutRealDataTest, CutCsahWithTightBox) +{ + // csah.shp bounds: + // X: 257105 .. 335498 Y: 5176099 .. 5226768 + // + // Use a very tight clip box in the center -> few shapes survive + OpenRealData("csah.shp"); + EXPECT_EQ(124, nEntities); + + cxmin = 290000; + cymin = 5195000; + cxmax = 300000; + cymax = 5205000; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, 0); + EXPECT_LT(outEntities, 124); + EXPECT_EQ(outEntities, DBFGetRecordCount(hDBFappend)); +} + +TEST_F(CutRealDataTest, CutTinyBoxYieldsNothing) +{ + // Clip box completely outside the data extent -> zero output + OpenRealData("pline.shp"); + + cxmin = 0; + cymin = 0; + cxmax = 1; + cymax = 1; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_EQ(0, outEntities); + EXPECT_EQ(0, DBFGetRecordCount(hDBFappend)); +} + +TEST_F(CutRealDataTest, CutHugeBoxPreservesAll) +{ + // Clip box encloses entire data extent -> no clipping needed + // (some null/deleted shapes may still be skipped by SHPReadObject) + OpenRealData("csah.shp"); + const int totalIn = nEntities; + + cxmin = 200000; + cymin = 5100000; + cxmax = 400000; + cymax = 5300000; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, totalIn / 2); // most shapes survive + EXPECT_EQ(outEntities, DBFGetRecordCount(hDBFappend)); +} + +TEST_F(CutRealDataTest, CutPolygonClipsToSubset) +{ + // polygon.shp bounds: + // X: 471127 .. 489292 Y: 4751545 .. 4765611 + // + // Clip a sub-region and verify polygon output + OpenRealData("polygon.shp"); + const int totalIn = nEntities; + EXPECT_EQ(474, totalIn); + + cxmin = 478000; + cymin = 4756000; + cxmax = 484000; + cymax = 4762000; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, 0); + EXPECT_LT(outEntities, totalIn); + + for (int i = 0; i < outEntities; i++) + { + auto *obj = SHPReadObject(hSHPappend, i); + ASSERT_NE(nullptr, obj); + EXPECT_GE(obj->nVertices, 4) << "CUT polygon must have >= 4 vertices"; + // Check each ring is closed + for (int p = 0; p < obj->nParts; p++) + { + const int start = obj->panPartStart[p]; + const int end = (p + 1 < obj->nParts) ? obj->panPartStart[p + 1] + : obj->nVertices; + EXPECT_DOUBLE_EQ(obj->padfX[start], obj->padfX[end - 1]) + << "shape " << i << " ring " << p << " not closed (X)"; + EXPECT_DOUBLE_EQ(obj->padfY[start], obj->padfY[end - 1]) + << "shape " << i << " ring " << p << " not closed (Y)"; + } + // All vertices inside or on clip box + for (int v = 0; v < obj->nVertices; v++) + { + EXPECT_GE(obj->padfX[v], cxmin - 1e-6) + << "shape " << i << " vertex " << v; + EXPECT_LE(obj->padfX[v], cxmax + 1e-6) + << "shape " << i << " vertex " << v; + EXPECT_GE(obj->padfY[v], cymin - 1e-6) + << "shape " << i << " vertex " << v; + EXPECT_LE(obj->padfY[v], cymax + 1e-6) + << "shape " << i << " vertex " << v; + } + SHPDestroyObject(obj); + } + + EXPECT_EQ(outEntities, DBFGetRecordCount(hDBFappend)); + + // Total clipped area must be positive. Overlapping shapes can + // exceed the clip box area, so only check a reasonable upper bound. + const double clipBoxArea = (cxmax - cxmin) * (cymax - cymin); + double totalArea = 0; + for (int i = 0; i < outEntities; i++) + { + auto *obj = SHPReadObject(hSHPappend, i); + ASSERT_NE(nullptr, obj); + totalArea += TotalPolygonArea(obj); + SHPDestroyObject(obj); + } + EXPECT_GT(totalArea, 0.0); + EXPECT_LT(totalArea, clipBoxArea * 10.0); +} + +TEST_F(CutRealDataTest, CutPolygonHugeBoxPreservesAll) +{ + // Clip box encloses entire polygon.shp extent + OpenRealData("polygon.shp"); + const int totalIn = nEntities; + + cxmin = 470000; + cymin = 4750000; + cxmax = 490000; + cymax = 4770000; + + process_records(); + + int outEntities = 0; + SHPGetInfo(hSHPappend, &outEntities, nullptr, nullptr, nullptr); + EXPECT_GT(outEntities, totalIn / 2); + EXPECT_EQ(outEntities, DBFGetRecordCount(hDBFappend)); + + // Total area must be positive (shapes preserved) + double totalArea = 0; + for (int i = 0; i < outEntities; i++) + { + auto *obj = SHPReadObject(hSHPappend, i); + ASSERT_NE(nullptr, obj); + totalArea += TotalPolygonArea(obj); + SHPDestroyObject(obj); + } + EXPECT_GT(totalArea, 0.0); +} + +} // namespace + +int main(int argc, char **argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +}