From 215bc83d75ad383719501159a2253dc9a726de3c Mon Sep 17 00:00:00 2001 From: "Thomas G. Lockhart" Date: Tue, 29 Jul 1997 16:08:18 +0000 Subject: [PATCH] Remove #ifdef'd support for old i/o styles. Change box terminology from "length" to "width". Use length terminology in common with other geometric types (usually perimeter). Fix bugs in line arithmetic which resulted in bad intersection calculations. Deprecate temporary unstored slope fields. Check explicitly for intersections at endpoints to avoid rounding ugliness. Add center() routines for lseg, path, polygon. Add distance() routines for circle-polygon, polygon-polygon. Check explicitly for points and polygons contained within polygons using an axis-crossing algorithm. (Old code just checked bounding boxes). Add routine to convert circle-box. *whew* --- src/backend/utils/adt/geo_ops.c | 1148 ++++++++++++++++++++++--------- 1 file changed, 811 insertions(+), 337 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index 778ea5ef19..320270e7ab 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -7,11 +7,12 @@ * * * IDENTIFICATION - * $Header: /cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v 1.12 1997/06/03 14:01:22 thomas Exp $ + * $Header: /cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v 1.13 1997/07/29 16:08:18 thomas Exp $ * *------------------------------------------------------------------------- */ #include +#include #include #include /* for sprintf proto, etc. */ #include /* for strtod, etc. */ @@ -23,8 +24,12 @@ #include "utils/geo_decls.h" #include "utils/palloc.h" -#define OLD_FORMAT_IN 0 -#define OLD_FORMAT_OUT 0 +#ifndef PI +#define PI 3.1415926536 +#endif + +int point_inside( Point *p, int npts, Point plist[]); +int lseg_crossing( double x, double y, double px, double py); /* * Delimiters for input and output strings. @@ -46,17 +51,6 @@ static int digits8 = P_MAXDIG; -int geo_precision(int digits); - -int geo_precision(int digits) -{ - if (digits > P_MAXDIG) { - digits8 = P_MAXDIG; - } else if (digits > 0) { - digits8 = digits; - }; - return digits8; -} /* * Geometric data types are composed of points. @@ -83,6 +77,8 @@ int geo_precision(int digits) * and restore that order for text output - tgl 97/01/16 */ +int single_decode(char *str, float8 *x, char **ss); +int single_encode(float8 x, char *str); int pair_decode(char *str, float8 *x, float8 *y, char **s); int pair_encode(float8 x, float8 y, char *str); int pair_count(char *s, char delim); @@ -90,6 +86,32 @@ int path_decode(int opentype, int npts, char *str, int *isopen, char **ss, Point char *path_encode( bool closed, int npts, Point *pt); +int single_decode(char *str, float8 *x, char **s) +{ + char *cp; + + if (!PointerIsValid(str)) + return(FALSE); + + while (isspace( *str)) str++; + *x = strtod( str, &cp); +#ifdef GEODEBUG +fprintf( stderr, "single_decode- (%x) try decoding %s to %g\n", (cp-str), str, *x); +#endif + if (cp <= str) return(FALSE); + while (isspace( *cp)) cp++; + + if (s != NULL) *s = cp; + + return(TRUE); +} /* single_decode() */ + +int single_encode(float8 x, char *str) +{ + (void) sprintf(str, "%.*g", digits8, x); + return(TRUE); +} /* single_encode() */ + int pair_decode(char *str, float8 *x, float8 *y, char **s) { int has_delim; @@ -284,39 +306,17 @@ BOX *box_in(char *str) }; return(box); -} +} /* box_in() */ /* box_out - convert a box to external form. */ char *box_out(BOX *box) { -#if OLD_FORMAT_OUT - char *result = PALLOC(2*(P_MAXLEN+1)+2); - - char *cp; -#endif - if (!PointerIsValid(box)) return(NULL); -#if OLD_FORMAT_OUT - cp = result; - *cp++ = LDELIM; - if (! pair_encode( box->high.x, box->high.y, cp)) - elog (WARN, "Unable to format box", NULL); - cp += strlen(cp); - *cp++ = DELIM; - if (! pair_encode( box->low.x, box->low.y, cp)) - elog (WARN, "Unable to format box", NULL); - cp += strlen(cp); - *cp++ = RDELIM; - *cp = '\0'; - - return( result); -#else return( path_encode( -1, 2, (Point *) &(box->high))); -#endif -} +} /* box_out() */ /* box_construct - fill in a new box. @@ -498,23 +498,23 @@ double *box_area(BOX *box) { double *result = PALLOCTYPE(double); - *result = box_ln(box) * box_ht(box); + *result = box_wd(box) * box_ht(box); return(result); } -/* box_length - returns the length of the box +/* box_width - returns the width of the box * (horizontal magnitude). */ -double *box_length(BOX *box) +double *box_width(BOX *box) { double *result = PALLOCTYPE(double); *result = box->high.x - box->low.x; return(result); -} +} /* box_width() */ /* box_height - returns the height of the box @@ -565,14 +565,14 @@ Point *box_center(BOX *box) */ double box_ar(BOX *box) { - return( box_ln(box) * box_ht(box) ); + return( box_wd(box) * box_ht(box) ); } -/* box_ln - returns the length of the box +/* box_wd - returns the width (length) of the box * (horizontal magnitude). */ -double box_ln(BOX *box) +double box_wd(BOX *box) { return( box->high.x - box->low.x ); } @@ -670,8 +670,11 @@ line_construct_pm(Point *pt, double m) result->A = m; result->B = -1.0; result->C = pt->y - m * pt->x; + + result->m = m; + return(result); -} +} /* line_construct_pm() */ LINE * /* two points */ @@ -681,19 +684,40 @@ line_construct_pp(Point *pt1, Point *pt2) if (FPeq(pt1->x, pt2->x)) { /* vertical */ /* use "x = C" */ - result->m = 0.0; - result->A = -1.0; - result->B = 0.0; + result->A = -1; + result->B = 0; result->C = pt1->x; +#ifdef GEODEBUG +printf( "line_construct_pp- line is vertical\n"); +#endif + result->m = DBL_MAX; + + } else if (FPeq(pt1->y, pt2->y)) { /* horizontal */ + /* use "x = C" */ + result->A = 0; + result->B = -1; + result->C = pt1->y; +#ifdef GEODEBUG +printf( "line_construct_pp- line is horizontal\n"); +#endif + result->m = 0.0; + } else { /* use "mx - y + yinter = 0" */ - result->m = (pt1->y - pt2->y) / (pt1->x - pt2->x); - result->A = result->m; +#if FALSE + result->A = (pt1->y - pt2->y) / (pt1->x - pt2->x); +#endif + result->A = (pt2->y - pt1->y) / (pt2->x - pt1->x); result->B = -1.0; - result->C = pt1->y - result->m * pt1->x; - } + result->C = pt1->y - result->A * pt1->x; +#ifdef GEODEBUG +printf( "line_construct_pp- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n", + digits8, (pt2->x - pt1->x), digits8, (pt2->y - pt1->y)); +#endif + result->m = result->A; + }; return(result); -} +} /* line_construct_pp() */ /*---------------------------------------------------------- @@ -707,32 +731,53 @@ bool line_intersect(LINE *l1, LINE *l2) bool line_parallel(LINE *l1, LINE *l2) { +#if FALSE return( FPeq(l1->m, l2->m) ); -} +#endif + if (FPzero(l1->B)) { + return(FPzero(l2->B)); + }; + + return(FPeq(l2->A, l1->A*(l2->B / l1->B))); +} /* line_parallel() */ bool line_perp(LINE *l1, LINE *l2) { +#if FALSE if (l1->m) return( FPeq(l2->m / l1->m, -1.0) ); else if (l2->m) return( FPeq(l1->m / l2->m, -1.0) ); - return(1); /* both 0.0 */ -} +#endif + if (FPzero(l1->A)) { + return( FPzero(l2->B) ); + } else if (FPzero(l1->B)) { + return( FPzero(l2->A) ); + }; + + return( FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0) ); +} /* line_perp() */ bool line_vertical(LINE *line) { +#if FALSE return( FPeq(line->A, -1.0) && FPzero(line->B) ); -} +#endif + return( FPzero(line->B) ); +} /* line_vertical() */ bool line_horizontal(LINE *line) { +#if FALSE return( FPzero(line->m) ); -} +#endif + return( FPzero(line->A) ); +} /* line_horizontal() */ bool line_eq(LINE *l1, LINE *l2) { - double k; + double k; if (! FPzero(l2->A)) k = l1->A / l2->A; @@ -742,9 +787,10 @@ bool line_eq(LINE *l1, LINE *l2) k = l1->C / l2->C; else k = 1.0; + return( FPeq(l1->A, k * l2->A) && - FPeq(l1->B, k * l2->B) && - FPeq(l1->C, k * l2->C) ); + FPeq(l1->B, k * l2->B) && + FPeq(l1->C, k * l2->C) ); } @@ -772,14 +818,18 @@ line_distance(LINE *l1, LINE *l2) return(result); } -Point * /* point where l1, l2 intersect (if any) */ +/* line_interpt() + * Point where two lines l1, l2 intersect (if any) + */ +Point * line_interpt(LINE *l1, LINE *l2) { Point *result; - double x; + double x, y; if (line_parallel(l1, l2)) return(NULL); +#if FALSE if (line_vertical(l1)) result = point_construct(l2->m * l1->C + l2->C, l1->C); else if (line_vertical(l2)) @@ -788,8 +838,42 @@ line_interpt(LINE *l1, LINE *l2) x = (l1->C - l2->C) / (l2->A - l1->A); result = point_construct(x, l1->m * x + l1->C); } +#endif + + if (line_vertical(l1)) { +#if FALSE + x = l1->C; + y = -((l2->A * x + l2->C) / l2->B); +#endif + x = l1->C; + y = (l2->A * x + l2->C); + + } else if (line_vertical(l2)) { +#if FALSE + x = l2->C; + y = -((l1->A * x + l1->C) / l1->B); +#endif + x = l2->C; + y = (l1->A * x + l1->C); + + } else { +#if FALSE + x = (l2->B * l1->C - l1->B * l2->C) / (l2->A * l1->B - l1->A * l2->B); + y = -((l1->A * x + l1->C) / l1->B); +#endif + x = (l1->C - l2->C) / (l2->A - l1->A); + y = (l1->A * x + l1->C); + }; + result = point_construct(x, y); + +#ifdef GEODEBUG +printf( "line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n", + digits8, l1->A, digits8, l1->B, digits8, l1->C, digits8, l2->A, digits8, l2->B, digits8, l2->C); +printf( "line_interpt- lines intersect at (%.*g,%.*g)\n", digits8, x, digits8, y); +#endif return(result); -} +} /* line_interpt() */ + /*********************************************************************** ** @@ -821,12 +905,7 @@ PATH *path_in(char *str) char *s; int npts; int size; -#if OLD_FORMAT_IN - int oldstyle = FALSE; - double x, y; -#else int depth = 0; -#endif if (!PointerIsValid(str)) elog(WARN, "Bad (null) path external representation"); @@ -837,27 +916,11 @@ PATH *path_in(char *str) s = str; while (isspace( *s)) s++; -#if OLD_FORMAT_IN - /* identify old style format as having only one left delimiter in string... */ - oldstyle = ((*s == LDELIM) && (strrchr( s, LDELIM) == s)); - - /* old-style format? then first two fields are closed flag and point count... */ - if (oldstyle) { - s++; - if ((! pair_decode( s, &x, &y, &s)) || (*s++ != DELIM) - || ((x != 0) && (x != 1)) || (y <= 0)) - elog (WARN, "Bad path external representation '%s'",str); - isopen = (x == 0); - npts = y; - }; - -#else /* skip single leading paren */ if ((*s == LDELIM) && (strrchr( s, LDELIM) == s)) { s++; depth++; }; -#endif size = offsetof(PATH, p[0]) + (sizeof(path->p[0]) * npts); path = PALLOC(size); @@ -865,67 +928,23 @@ PATH *path_in(char *str) path->size = size; path->npts = npts; -#if OLD_FORMAT_IN - if (oldstyle) path->closed = (! isopen); - - if ((! path_decode(TRUE, npts, s, &isopen, &s, &(path->p[0]))) - || ! (oldstyle? (*s++ == RDELIM): (*s == '\0'))) - -#else if ((!path_decode(TRUE, npts, s, &isopen, &s, &(path->p[0]))) && (!((depth == 0) && (*s == '\0'))) && !((depth >= 1) && (*s == RDELIM))) -#endif elog (WARN, "Bad path external representation '%s'",str); -#if OLD_FORMAT_IN - if (oldstyle) { - while (isspace( *s)) s++; - if (*s != '\0') - elog (WARN, "Bad path external representation '%s'",str); - }; - - if (! oldstyle) path->closed = (! isopen); - -#else path->closed = (! isopen); -#endif return(path); -} +} /* path_in() */ char *path_out(PATH *path) { -#if OLD_FORMAT_OUT - int i; - char *result, *cp; -#endif - if (!PointerIsValid(path)) return NULL; -#if OLD_FORMAT_OUT - result = PALLOC(path->npts*(P_MAXLEN+3)+2); - - cp = result; - *cp++ = LDELIM; - if (! pair_encode( path->closed, path->npts, cp)) - elog (WARN, "Unable to format path", NULL); - cp += strlen(cp); - - for (i=0; inpts; i++) { - *cp++ = DELIM; - if (! pair_encode( path->p[i].x, path->p[i].y, cp)) - elog (WARN, "Unable to format path", NULL); - cp += strlen(cp); - }; - *cp++ = RDELIM; - *cp = '\0'; - return(result); -#else return( path_encode( path->closed, path->npts, (Point *) &(path->p[0]))); -#endif -} +} /* path_out() */ /*---------------------------------------------------------- @@ -1010,6 +1029,7 @@ path_close(PATH *path) return(result); } /* path_close() */ + PATH * path_open(PATH *path) { @@ -1122,29 +1142,32 @@ double *path_distance(PATH *p1, PATH *p2) double *path_length(PATH *path) { - double *result = PALLOCTYPE(double); - int ct, i; - - ct = path->npts - 1; - for (i = 0; i < ct; i++) + double *result; + int i; + + result = PALLOCTYPE(double); + + *result = 0; + for (i = 0; i < (path->npts - 1); i++) *result += point_dt(&path->p[i], &path->p[i+1]); - + return(result); -} +} /* path_length() */ double path_ln(PATH *path) { double result; - int ct, i; - - ct = path->npts - 1; - for (result = i = 0; i < ct; i++) + int i; + + result = 0; + for (i = 0; i < (path->npts - 1); i++) result += point_dt(&path->p[i], &path->p[i+1]); - + return(result); -} +} /* path_ln() */ + /*********************************************************************** ** ** Routines for 2D points. @@ -1166,10 +1189,8 @@ point_in(char *str) double x, y; char *s; - if (str == NULL) { + if (! PointerIsValid( str)) elog(WARN, "Bad (null) point external representation"); - return NULL; - } if (! pair_decode( str, &x, &y, &s) || (strlen(s) > 0)) elog (WARN, "Bad point external representation '%s'",str); @@ -1185,7 +1206,7 @@ point_in(char *str) char * point_out(Point *pt) { - if (!PointerIsValid(pt)) + if (! PointerIsValid(pt)) return(NULL); return( path_encode( -1, 1, pt)); @@ -1204,7 +1225,12 @@ Point *point_construct(double x, double y) Point *point_copy(Point *pt) { - Point *result = PALLOCTYPE(Point); + Point *result; + + if (! PointerIsValid( pt)) + return(NULL); + + result = PALLOCTYPE(Point); result->x = pt->x; result->y = pt->y; @@ -1336,7 +1362,7 @@ LSEG *lseg_in(char *str) lseg->m = point_sl(&lseg->p[0], &lseg->p[1]); return(lseg); -} +} /* lseg_in() */ char *lseg_out(LSEG *ls) @@ -1345,7 +1371,7 @@ char *lseg_out(LSEG *ls) return(NULL); return( path_encode( FALSE, 2, (Point *) &(ls->p[0]))); -} +} /* lseg_out() */ /* lseg_construct - @@ -1403,17 +1429,26 @@ bool lseg_intersect(LSEG *l1, LSEG *l2) bool lseg_parallel(LSEG *l1, LSEG *l2) { +#if FALSE return( FPeq(l1->m, l2->m) ); -} +#endif + return( FPeq( point_sl( &(l1->p[0]), &(l1->p[1])), + point_sl( &(l2->p[0]), &(l2->p[1]))) ); +} /* lseg_parallel() */ bool lseg_perp(LSEG *l1, LSEG *l2) { - if (! FPzero(l1->m)) - return( FPeq(l2->m / l1->m, -1.0) ); - else if (! FPzero(l2->m)) - return( FPeq(l1->m / l2->m, -1.0) ); + double m1, m2; + + m1 = point_sl( &(l1->p[0]), &(l1->p[1])); + m2 = point_sl( &(l2->p[0]), &(l2->p[1])); + + if (! FPzero(m1)) + return( FPeq(m2 / m1, -1.0) ); + else if (! FPzero(m2)) + return( FPeq(m1 / m2, -1.0) ); return(0); /* both 0.0 */ -} +} /* lseg_perp() */ bool lseg_vertical(LSEG *lseg) { @@ -1454,7 +1489,8 @@ double *lseg_distance(LSEG *l1, LSEG *l2) } /* distance between l1, l2 */ -double lseg_dt(LSEG *l1, LSEG *l2) +double +lseg_dt(LSEG *l1, LSEG *l2) { double *d, result; @@ -1467,15 +1503,37 @@ double lseg_dt(LSEG *l1, LSEG *l2) d = dist_ps(&l1->p[1], l2); result = Min(result, *d); PFREE(d); +#if FALSE +/* XXX Why are we checking distances from all endpoints to the other segment? + * One set of endpoints should be sufficient - tgl 97/07/03 + */ d = dist_ps(&l2->p[0], l1); result = Min(result, *d); PFREE(d); d = dist_ps(&l2->p[1], l1); result = Min(result, *d); PFREE(d); +#endif return(result); -} +} /* lseg_dt() */ + + +Point * +lseg_center(LSEG *lseg) +{ + Point *result; + + if (!PointerIsValid(lseg)) + return(NULL); + + result = PALLOCTYPE(Point); + + result->x = (lseg->p[0].x - lseg->p[1].x) / 2; + result->y = (lseg->p[0].y - lseg->p[1].y) / 2; + + return(result); +} /* lseg_center() */ /* lseg_interpt - @@ -1483,25 +1541,44 @@ double lseg_dt(LSEG *l1, LSEG *l2) * Find the intersection of the appropriate lines; if the * point is not on a given segment, there is no valid segment * intersection point at all. + * If there is an intersection, then check explicitly for matching + * endpoints since there may be rounding effects with annoying + * lsb residue. - tgl 1997-07-09 */ -Point *lseg_interpt(LSEG *l1, LSEG *l2) +Point * +lseg_interpt(LSEG *l1, LSEG *l2) { Point *result; LINE *tmp1, *tmp2; - + + if (!PointerIsValid(l1) || !PointerIsValid(l2)) + return(NULL); + tmp1 = line_construct_pp(&l1->p[0], &l1->p[1]); tmp2 = line_construct_pp(&l2->p[0], &l2->p[1]); result = line_interpt(tmp1, tmp2); - if (result) - if (! on_ps(result, l1)) { + if (PointerIsValid(result)) { + if (on_ps(result, l1)) { + if ((FPeq( l1->p[0].x, l2->p[0].x) && FPeq( l1->p[0].y, l2->p[0].y)) + || (FPeq( l1->p[0].x, l2->p[1].x) && FPeq( l1->p[0].y, l2->p[1].y))) { + result->x = l1->p[0].x; + result->y = l1->p[0].y; + + } else if ((FPeq( l1->p[1].x, l2->p[0].x) && FPeq( l1->p[1].y, l2->p[0].y)) + || (FPeq( l1->p[1].x, l2->p[1].x) && FPeq( l1->p[1].y, l2->p[1].y))) { + result->x = l1->p[1].x; + result->y = l1->p[1].y; + }; + } else { PFREE(result); result = NULL; - } - + }; + }; PFREE(tmp1); PFREE(tmp2); + return(result); -} +} /* lseg_interpt() */ /*********************************************************************** ** @@ -1536,27 +1613,49 @@ double *dist_ps(Point *pt, LSEG *lseg) LINE *ln; double *result, *tmpdist; Point *ip; - - /* construct a line that's perpendicular. See if the intersection of - the two lines is on the line segment. */ - if (lseg->p[1].x == lseg->p[0].x) + +/* + * Construct a line perpendicular to the input segment + * and through the input point + */ + if (lseg->p[1].x == lseg->p[0].x) { m = 0; - else if (lseg->p[1].y == lseg->p[0].y) /* slope is infinite */ + } else if (lseg->p[1].y == lseg->p[0].y) { /* slope is infinite */ m = (double)DBL_MAX; - else m = (-1) * (lseg->p[1].y - lseg->p[0].y) / - (lseg->p[1].x - lseg->p[0].x); + } else { +#if FALSE + m = (-1) * (lseg->p[1].y - lseg->p[0].y) / + (lseg->p[1].x - lseg->p[0].x); +#endif + m = ((lseg->p[0].y - lseg->p[1].y) / (lseg->p[1].x - lseg->p[0].x)); + }; ln = line_construct_pm(pt, m); - - if ((ip = interpt_sl(lseg, ln)) != NULL) + +#ifdef GEODEBUG +printf( "dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n", + ln->A, ln->B, ln->C, pt->x, pt->y, m); +#endif + +/* + * Calculate distance to the line segment + * or to the endpoints of the segment. + */ + + /* intersection is on the line segment? */ + if ((ip = interpt_sl(lseg, ln)) != NULL) { result = point_distance(pt, ip); - else /* intersection is not on line segment, so distance is min - of distance from point to an endpoint */ - { +#ifdef GEODEBUG +printf( "dist_ps- distance is %f to intersection point is (%f,%f)\n", + *result, ip->x, ip->y); +#endif + + /* otherwise, intersection is not on line segment */ + } else { result = point_distance(pt, &lseg->p[0]); tmpdist = point_distance(pt, &lseg->p[1]); if (*tmpdist < *result) *result = *tmpdist; PFREE (tmpdist); - } + }; if (ip != NULL) PFREE(ip); PFREE(ln); @@ -1567,7 +1666,7 @@ double *dist_ps(Point *pt, LSEG *lseg) /* ** Distance from a point to a path */ -double *dist_ppth(Point *pt, PATH *path) +double *dist_ppath(Point *pt, PATH *path) { double *result; double *tmp; @@ -1675,6 +1774,58 @@ double *dist_lb(LINE *line, BOX *box) } +double * +dist_cpoly(CIRCLE *circle, POLYGON *poly) +{ + double *result; + int i; + double *d; + LSEG seg; + + if (!PointerIsValid(circle) || !PointerIsValid(poly)) + elog (WARN, "Invalid (null) input for distance", NULL); + + if (point_inside( &(circle->center), poly->npts, poly->p)) { +#ifdef GEODEBUG +printf( "dist_cpoly- center inside of polygon\n"); +#endif + result = PALLOCTYPE(double); + + *result = 0; + return(result); + }; + + /* initialize distance with segment between first and last points */ + seg.p[0].x = poly->p[0].x; + seg.p[0].y = poly->p[0].y; + seg.p[1].x = poly->p[poly->npts-1].x; + seg.p[1].y = poly->p[poly->npts-1].y; + result = dist_ps( &(circle->center), &seg); +#ifdef GEODEBUG +printf( "dist_cpoly- segment 0/n distance is %f\n", *result); +#endif + + /* check distances for other segments */ + for (i = 0; (i < poly->npts - 1); i++) { + seg.p[0].x = poly->p[i].x; + seg.p[0].y = poly->p[i].y; + seg.p[1].x = poly->p[i+1].x; + seg.p[1].y = poly->p[i+1].y; + d = dist_ps( &(circle->center), &seg); +#ifdef GEODEBUG +printf( "dist_cpoly- segment %d distance is %f\n", (i+1), *d); +#endif + if (*d < *result) *result = *d; + PFREE(d); + }; + + *result -= circle->radius; + if (*result < 0) *result = 0; + + return(result); +} /* dist_cpoly() */ + + /*--------------------------------------------------------------------- * interpt_ * Intersection point of objects. @@ -1689,11 +1840,26 @@ Point *interpt_sl(LSEG *lseg, LINE *line) tmp = line_construct_pp(&lseg->p[0], &lseg->p[1]); p = line_interpt(tmp, line); - if (p) - if (! on_ps(p, lseg)) { +#ifdef GEODEBUG +printf( "interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n", + digits8, lseg->p[0].x, digits8, lseg->p[0].y, digits8, lseg->p[1].x, digits8, lseg->p[1].y); +printf( "interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n", + digits8, tmp->A, digits8, tmp->B, digits8, tmp->C); +#endif + if (PointerIsValid(p)) { +#ifdef GEODEBUG +printf( "interpt_sl- intersection point is (%.*g %.*g)\n", digits8, p->x, digits8, p->y); +#endif + if (on_ps(p, lseg)) { +#ifdef GEODEBUG +printf( "interpt_sl- intersection point is on segment\n"); +#endif + + } else { PFREE(p); p = NULL; - } + }; + }; PFREE(tmp); return(p); @@ -1716,21 +1882,32 @@ Point *close_pl(Point *pt, LINE *line) double invm; result = PALLOCTYPE(Point); +#if FALSE if (FPeq(line->A, -1.0) && FPzero(line->B)) { /* vertical */ +#endif + if (line_vertical(line)) { result->x = line->C; result->y = pt->y; return(result); + +#if FALSE } else if (FPzero(line->m)) { /* horizontal */ +#endif + } else if (line_horizontal(line)) { result->x = pt->x; result->y = line->C; return(result); } /* drop a perpendicular and find the intersection point */ +#if FALSE invm = -1.0 / line->m; +#endif + /* invert and flip the sign on the slope to get a perpendicular */ + invm = line->B / line->A; tmp = line_construct_pm(pt, invm); result = line_interpt(tmp, line); return(result); -} +} /* close_pl() */ /* close_ps - @@ -1758,26 +1935,36 @@ Point *close_ps(Point *pt, LSEG *lseg) result = point_copy(&lseg->p[yh]); if (result) return(result); +#if FALSE if (FPeq(lseg->p[0].x, lseg->p[1].x)) { /* vertical */ +#endif + if (lseg_vertical(lseg)) { result->x = lseg->p[0].x; result->y = pt->y; return(result); +#if FALSE } else if (FPzero(lseg->m)) { /* horizontal */ +#endif + } else if (lseg_horizontal(lseg)) { result->x = pt->x; result->y = lseg->p[0].y; return(result); } +#if FALSE invm = -1.0 / lseg->m; +#endif + invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1])); tmp = line_construct_pm(pt, invm); result = interpt_sl(lseg, tmp); return(result); -} +} /* close_ps() */ Point *close_pb(Point *pt, BOX *box) { /* think about this one for a while */ - + elog(WARN, "close_pb not implemented", NULL); + return(NULL); } @@ -1804,14 +1991,16 @@ Point *close_sl(LSEG *lseg, LINE *line) Point *close_sb(LSEG *lseg, BOX *box) { /* think about this one for a while */ - + elog(WARN, "close_sb not implemented", NULL); + return(NULL); } Point *close_lb(LINE *line, BOX *box) { /* think about this one for a while */ - + elog(WARN, "close_lb not implemented", NULL); + return(NULL); } @@ -1825,21 +2014,31 @@ Point *close_lb(LINE *line, BOX *box) */ bool on_pl(Point *pt, LINE *line) { + if (!PointerIsValid(pt) || !PointerIsValid(line)) + return(FALSE); + return( FPzero(line->A * pt->x + line->B * pt->y + line->C) ); } /* on_ps - * Determine colinearity by detecting a triangle inequality. + * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09 */ bool on_ps(Point *pt, LSEG *lseg) { + if (!PointerIsValid(pt) || !PointerIsValid(lseg)) + return(FALSE); + return( FPeq (point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]), point_dt(&lseg->p[0], &lseg->p[1])) ); } bool on_pb(Point *pt, BOX *box) { + if (!PointerIsValid(pt) || !PointerIsValid(box)) + return(FALSE); + return( pt->x <= box->high.x && pt->x >= box->low.x && pt->y <= box->high.y && pt->y >= box->low.y ); } @@ -1847,6 +2046,7 @@ bool on_pb(Point *pt, BOX *box) /* on_ppath - * Whether a point lies within (on) a polyline. * If open, we have to (groan) check each segment. + * (uses same algorithm as for point intersecting segment - tgl 1997-07-09) * If closed, we use the old O(n) ray method for point-in-polygon. * The ray is horizontal, from pt out to the right. * Each segment that crosses the ray counts as an @@ -1858,13 +2058,18 @@ bool on_pb(Point *pt, BOX *box) bool on_ppath(Point *pt, PATH *path) { +#if FALSE int above, next, /* is the seg above the ray? */ inter, /* # of times path crosses ray */ - hi, /* index inc of higher seg (0,1) */ - i, n; - double a, b, x, - yh, yl, xh, xl; - + hi; /* index inc of higher seg (0,1) */ + double x, yh, yl, xh, xl; +#endif + int i, n; + double a, b; + + if (!PointerIsValid(pt) || !PointerIsValid(path)) + return(FALSE); + if (! path->closed) { /*-- OPEN --*/ n = path->npts - 1; a = point_dt(pt, &path->p[0]); @@ -1878,6 +2083,8 @@ bool on_ppath(Point *pt, PATH *path) return(0); } + return(point_inside( pt, path->npts, path->p)); +#if FALSE inter = 0; /*-- CLOSED --*/ above = FPgt(path->p[0].y, pt->y) ? ABOVE : FPlt(path->p[0].y, pt->y) ? BELOW : UNDEF; @@ -1922,18 +2129,25 @@ bool on_ppath(Point *pt, PATH *path) } return( above == UNDEF || /* path is horizontal */ inter % 2); /* odd # of intersections */ -} +#endif +} /* on_ppath() */ bool on_sl(LSEG *lseg, LINE *line) { + if (!PointerIsValid(lseg) || !PointerIsValid(line)) + return(FALSE); + return( on_pl(&lseg->p[0], line) && on_pl(&lseg->p[1], line) ); -} +} /* on_sl() */ bool on_sb(LSEG *lseg, BOX *box) { + if (!PointerIsValid(lseg) || !PointerIsValid(box)) + return(FALSE); + return( on_pb(&lseg->p[0], box) && on_pb(&lseg->p[1], box) ); -} +} /* on_sb() */ /*--------------------------------------------------------------------- * inter_ @@ -1944,6 +2158,9 @@ bool inter_sl(LSEG *lseg, LINE *line) { Point *tmp; + if (!PointerIsValid(lseg) || !PointerIsValid(line)) + return(FALSE); + tmp = interpt_sl(lseg, line); if (tmp) { PFREE(tmp); @@ -2013,12 +2230,6 @@ POLYGON *poly_in(char *str) int size; int isopen; char *s; -#if OLD_FORMAT_IN - int oldstyle; - int oddcount; - int i; - double x1, x2; -#endif if (!PointerIsValid(str)) elog (WARN," Bad (null) polygon external representation"); @@ -2033,62 +2244,10 @@ POLYGON *poly_in(char *str) poly->size = size; poly->npts = npts; -#if OLD_FORMAT_IN - s = str; - while (isspace( *s)) s++; - /* identify old style format as having only one left delimiter in string... */ - oldstyle = ((*s == LDELIM) && (strrchr( s, LDELIM) == s)); - - if (oldstyle) { - s++; - while (isspace( *s)) s++; - - for (i=0; ip[i*2].x = x1; - poly->p[i*2+1].x = x2; - }; - oddcount = (npts % 2); - if (oddcount) { - if (! pair_decode( s, &x1, &x2, &s)) - elog (WARN, "Bad polygon external representation '%s'",str); - - if (*s == DELIM) s++; - poly->p[npts-1].x = x1; - poly->p[0].y = x2; - }; - for (i=0; ip[i*2+oddcount].y = x1; - poly->p[i*2+1+oddcount].y = x2; - }; - - if (*s == RDELIM) { - s++; - while (isspace( *s)) s++; - if (*s != '\0') - elog(WARN, "Bad polygon external representation '%s'", str); - - } else { - elog(WARN, "Bad polygon external representation '%s'", str); - }; - - } else { -#endif - if ((! path_decode(FALSE, npts, str, &isopen, &s, &(poly->p[0]))) - || (*s != '\0')) + if ((! path_decode(FALSE, npts, str, &isopen, &s, &(poly->p[0]))) + || (*s != '\0')) elog (WARN, "Bad polygon external representation '%s'",str); -#if OLD_FORMAT_IN - }; -#endif - make_bound_box(poly); return( poly); @@ -2101,33 +2260,11 @@ POLYGON *poly_in(char *str) *---------------------------------------------------------------*/ char *poly_out(POLYGON *poly) { -#if OLD_FORMAT_OUT - int i; - char *result, *cp; -#endif - if (!PointerIsValid(poly)) return NULL; -#if OLD_FORMAT_OUT - result = PALLOC(poly->npts*(P_MAXLEN+3)+2); - - cp = result; - *cp++ = LDELIM; - - for (i=0; inpts; i++) { - if (! pair_encode( poly->p[i].x, poly->p[i].y, cp)) - elog (WARN, "Unable to format polygon", NULL); - cp += strlen(cp); - *cp++ = DELIM; - }; - *(cp-1) = RDELIM; - *cp = '\0'; - return(result); -#else return( path_encode( TRUE, poly->npts, &(poly->p[0]))); -#endif -} +} /* poly_out() */ /*------------------------------------------------------- @@ -2173,20 +2310,29 @@ bool poly_overright(POLYGON *polya, POLYGON *polyb) /*------------------------------------------------------- * Is polygon A the same as polygon B? i.e. are all the * points the same? + * Check all points for matches in both forward and reverse + * direction since polygons are non-directional and are + * closed shapes. *-------------------------------------------------------*/ bool poly_same(POLYGON *polya, POLYGON *polyb) { - int i; + if (! PointerIsValid( polya) || ! PointerIsValid( polyb)) + return FALSE; + if (polya->npts != polyb->npts) return FALSE; + return(plist_same( polya->npts, polya->p, polyb->p)); + +#if FALSE for (i = 0; i < polya->npts; i++) { if ((polya->p[i].x != polyb->p[i].x) || (polya->p[i].y != polyb->p[i].y)) return FALSE; }; return TRUE; -} +#endif +} /* poly_same() */ /*----------------------------------------------------------------- * Determine if polygon A overlaps polygon B by determining if @@ -2197,23 +2343,114 @@ bool poly_overlap(POLYGON *polya, POLYGON *polyb) return box_overlap(&(polya->boundbox), &(polyb->boundbox)); } + /*----------------------------------------------------------------- * Determine if polygon A contains polygon B by determining if A's * bounding box contains B's bounding box. *-----------------------------------------------------------------*/ +#if FALSE bool poly_contain(POLYGON *polya, POLYGON *polyb) { return box_contain(&(polya->boundbox), &(polyb->boundbox)); } +#endif + +bool +poly_contain(POLYGON *polya, POLYGON *polyb) +{ + int i; + + if (!PointerIsValid(polya) || !PointerIsValid(polyb)) + return(FALSE); + + if (box_contain(&(polya->boundbox), &(polyb->boundbox))) { + for (i = 0; i < polyb->npts; i++) { + if (point_inside(&(polyb->p[i]), polya->npts, &(polya->p[0])) == 0) { +#if GEODEBUG +printf( "poly_contain- point (%f,%f) not in polygon\n", polyb->p[i].x, polyb->p[i].y); +#endif + return(FALSE); + }; + }; + for (i = 0; i < polya->npts; i++) { + if (point_inside(&(polya->p[i]), polyb->npts, &(polyb->p[0])) == 1) { +#if GEODEBUG +printf( "poly_contain- point (%f,%f) in polygon\n", polya->p[i].x, polya->p[i].y); +#endif + return(FALSE); + }; + }; + return(TRUE); + }; +#if GEODEBUG +printf( "poly_contain- bound box ((%f,%f),(%f,%f)) not inside ((%f,%f),(%f,%f))\n", + polyb->boundbox.low.x,polyb->boundbox.low.y,polyb->boundbox.high.x,polyb->boundbox.high.y, + polya->boundbox.low.x,polya->boundbox.low.y,polya->boundbox.high.x,polya->boundbox.high.y); +#endif + return(FALSE); +} /* poly_contain() */ + /*----------------------------------------------------------------- * Determine if polygon A is contained by polygon B by determining * if A's bounding box is contained by B's bounding box. *-----------------------------------------------------------------*/ +#if FALSE bool poly_contained(POLYGON *polya, POLYGON *polyb) { - return box_contained(&(polya->boundbox), &(polyb->boundbox)); + return(box_contained(&(polya->boundbox), &(polyb->boundbox))); } +#endif + +bool poly_contained(POLYGON *polya, POLYGON *polyb) +{ + return(poly_contain(polyb, polya)); +} /* poly_contained() */ + + +/* poly_contain_pt() + * Test to see if the point is inside the polygon. + * Code adapted from integer-based routines in + * Wn: A Server for the HTTP + * File: wn/image.c + * Version 1.15.1 + * Copyright (C) 1995 + * (code offered for use by J. Franks in Linux Journal letter.) + */ + +bool +poly_contain_pt( POLYGON *poly, Point *p) +{ + if (!PointerIsValid(poly) || !PointerIsValid(p)) + return(FALSE); + + return(point_inside(p, poly->npts, &(poly->p[0])) != 0); +} /* poly_contain_pt() */ + +bool +pt_contained_poly( Point *p, POLYGON *poly) +{ + if (!PointerIsValid(p) || !PointerIsValid(poly)) + return(FALSE); + + return(poly_contain_pt( poly, p)); +} /* pt_contained_poly() */ + + +double * +poly_distance( POLYGON *polya, POLYGON *polyb) +{ + double *result; + + if (!PointerIsValid(polya) || !PointerIsValid(polyb)) + return(NULL); + + result = PALLOCTYPE(double); + + *result = 0; + + return(result); +} /* poly_distance() */ /*********************************************************************** @@ -2402,8 +2639,6 @@ box_div(BOX *box, Point *p) ** ***********************************************************************/ -POLYGON *path_poly(PATH *path); - /* path_add() * Concatenate two paths (only if they are both open). */ @@ -2527,6 +2762,41 @@ path_div_pt(PATH *path, Point *point) } /* path_div_pt() */ +bool +path_contain_pt( PATH *path, Point *p) +{ + if (!PointerIsValid(path) || !PointerIsValid(p)) + return(FALSE); + + return( (path->closed? (point_inside(p, path->npts, &(path->p[0])) != 0): FALSE)); +} /* path_contain_pt() */ + +bool +pt_contained_path( Point *p, PATH *path) +{ + if (!PointerIsValid(p) || !PointerIsValid(path)) + return(FALSE); + + return( path_contain_pt( path, p)); +} /* pt_contained_path() */ + + +Point * +path_center(PATH *path) +{ + Point *result; + + if (!PointerIsValid(path)) + return(NULL); + + elog(WARN, "path_center not implemented", NULL); + + result = PALLOCTYPE(Point); + result = NULL; + + return(result); +} /* path_center() */ + POLYGON *path_poly(PATH *path) { POLYGON *poly; @@ -2597,7 +2867,7 @@ bool isoldpath(PATH *path) { if (!PointerIsValid(path) || (path->npts < 2)) - return(0); + return(FALSE); return(path->npts == (path->p[0].y+1)); } /* isoldpath() */ @@ -2610,7 +2880,7 @@ isoldpath(PATH *path) ***********************************************************************/ int4 -poly_npoints( POLYGON *poly) +poly_npoints(POLYGON *poly) { if (!PointerIsValid(poly)) return(0); @@ -2618,6 +2888,28 @@ poly_npoints( POLYGON *poly) return(poly->npts); } /* poly_npoints() */ + +Point * +poly_center(POLYGON *poly) +{ + Point *result; + CIRCLE *circle; + + if (!PointerIsValid(poly)) + return(NULL); + + if (PointerIsValid(circle = poly_circle(poly))) { + result = circle_center(circle); + PFREE(circle); + + } else { + result = NULL; + }; + + return(result); +} /* poly_center() */ + + BOX * poly_box(POLYGON *poly) { @@ -2631,6 +2923,7 @@ poly_box(POLYGON *poly) return(box); } /* poly_box() */ + /* box_poly() * Convert a box to a polygon. */ @@ -2664,6 +2957,7 @@ box_poly(BOX *box) return(poly); } /* box_poly() */ + PATH * poly_path(POLYGON *poly) { @@ -2773,53 +3067,6 @@ POLYGON } /* revertpoly() */ -/*------------------------------------------------------------------------- - * - * circle.c-- - * 2D geometric operations - * - * Copyright (c) 1994, Regents of the University of California - * - * - * IDENTIFICATION - * $Header: /cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v 1.12 1997/06/03 14:01:22 thomas Exp $ - * - *------------------------------------------------------------------------- - */ -#ifndef PI -#define PI 3.1415926536 -#endif - -int single_decode(char *str, float8 *x, char **ss); -int single_encode(float8 x, char *str); - -int single_decode(char *str, float8 *x, char **s) -{ - char *cp; - - if (!PointerIsValid(str)) - return(FALSE); - - while (isspace( *str)) str++; - *x = strtod( str, &cp); -#ifdef GEODEBUG -fprintf( stderr, "single_decode- (%x) try decoding %s to %g\n", (cp-str), str, *x); -#endif - if (cp <= str) return(FALSE); - while (isspace( *cp)) cp++; - - if (s != NULL) *s = cp; - - return(TRUE); -} - -int single_encode(float8 x, char *str) -{ - (void) sprintf(str, "%.*g", digits8, x); - return(TRUE); -} - - /*********************************************************************** ** ** Routines for circles. @@ -3191,6 +3438,30 @@ double *circle_distance(CIRCLE *circle1, CIRCLE *circle2) } /* circle_distance() */ +bool +circle_contain_pt(CIRCLE *circle, Point *point) +{ + bool within; + double *d; + + if (!PointerIsValid(circle) || !PointerIsValid(point)) + return(FALSE); + + d = point_distance(&(circle->center), point); + within = (*d <= circle->radius); + PFREE(d); + + return(within); +} /* circle_contain_pt() */ + + +bool +pt_contained_circle(Point *point, CIRCLE *circle) +{ + return(circle_contain_pt(circle,point)); +} /* circle_contain_pt() */ + + /* dist_pc - returns the distance between * a point and a circle. */ @@ -3199,6 +3470,7 @@ double *dist_pc(Point *point, CIRCLE *circle) double *result; result = PALLOCTYPE(double); + *result = (point_dt(point,&circle->center) - circle->radius); if (*result < 0) *result = 0; @@ -3261,6 +3533,50 @@ CIRCLE *circle(Point *center, float8 *radius) return(result); } + +BOX * +circle_box(CIRCLE *circle) +{ + BOX *box; + double delta; + + if (!PointerIsValid(circle)) + return(NULL); + + box = PALLOCTYPE(BOX); + + delta = circle->radius / sqrt(2.0e0); + + box->high.x = circle->center.x + delta; + box->low.x = circle->center.x - delta; + box->high.y = circle->center.y + delta; + box->low.y = circle->center.y - delta; + + return(box); +} /* circle_box() */ + +/* box_circle() + * Convert a box to a circle. + */ +CIRCLE * +box_circle(BOX *box) +{ + CIRCLE *circle; + + if (!PointerIsValid(box)) + return(NULL); + + circle = PALLOCTYPE(CIRCLE); + + circle->center.x = (box->high.x + box->low.x) / 2; + circle->center.y = (box->high.y + box->low.y) / 2; + + circle->radius = point_dt(&circle->center, &box->high); + + return(circle); +} /* box_circle() */ + + POLYGON *circle_poly(int npts, CIRCLE *circle) { POLYGON *poly; @@ -3271,7 +3587,7 @@ POLYGON *circle_poly(int npts, CIRCLE *circle) if (!PointerIsValid(circle)) return(NULL); - if (FPzero(circle->radius) || (npts <= 2)) + if (FPzero(circle->radius) || (npts < 2)) elog (WARN, "Unable to convert circle to polygon", NULL); size = offsetof(POLYGON, p[0]) + (sizeof(poly->p[0]) * npts); @@ -3305,7 +3621,7 @@ CIRCLE *poly_circle(POLYGON *poly) if (!PointerIsValid(poly)) return(NULL); - if (poly->npts <= 2) + if (poly->npts < 2) elog (WARN, "Unable to convert polygon to circle", NULL); circle = PALLOCTYPE(CIRCLE); @@ -3330,4 +3646,162 @@ CIRCLE *poly_circle(POLYGON *poly) elog (WARN, "Unable to convert polygon to circle", NULL); return(circle); -} +} /* poly_circle() */ + + +/*********************************************************************** + ** + ** Private routines for multiple types. + ** + ***********************************************************************/ + +#define HIT_IT INT_MAX + +int +point_inside( Point *p, int npts, Point plist[]) +{ + double x0, y0; + double px, py; + + int i; + double x, y; + int cross, crossnum; + +/* + * We calculate crossnum, which is twice the crossing number of a + * ray from the origin parallel to the positive X axis. + * A coordinate change is made to move the test point to the origin. + * Then the function lseg_crossing() is called to calculate the crossnum of + * one segment of the translated polygon with the ray which is the + * positive X-axis. + */ + + crossnum = 0; + i = 0; + if (npts <= 0) return 0; + + x0 = plist[0].x - p->x; + y0 = plist[0].y - p->y; + + px = x0; + py = y0; + for (i = 1; i < npts; i++) { + x = plist[i].x - p->x; + y = plist[i].y - p->y; + + if ( (cross = lseg_crossing( x, y, px, py)) == HIT_IT ) { + return 2; + } + crossnum += cross; + + px = x; + py = y; + } + if ( (cross = lseg_crossing( x0, y0, px, py)) == HIT_IT ) { + return 2; + } + crossnum += cross; + if ( crossnum != 0 ) { + return 1; + } + return 0; +} /* point_inside() */ + + +/* lseg_crossing() + * The function lseg_crossing() returns +2, or -2 if the segment from (x,y) + * to previous (x,y) crosses the positive X-axis positively or negatively. + * It returns +1 or -1 if one endpoint is on this ray, or 0 if both are. + * It returns 0 if the ray and the segment don't intersect. + * It returns HIT_IT if the segment contains (0,0) + */ + +int +lseg_crossing( double x, double y, double px, double py) +{ + double z; + int sgn; + + /* If (px,py) = (0,0) and not first call we have already sent HIT_IT */ + + if (FPzero( y)) { + if (FPzero( x)) { + return(HIT_IT); + + } else if (FPgt( x, 0)) { + if (FPzero( py)) return(FPgt( px, 0)? 0 : HIT_IT); + return(FPlt( py, 0)? 1 : -1); + + } else { /* x < 0 */ + if (FPzero( py)) return(FPlt( px, 0)? 0 : HIT_IT); + return(0); + }; + }; + + /* Now we know y != 0; set sgn to sign of y */ + sgn = (FPgt( y, 0)? 1 : -1); + if (FPzero( py)) return(FPlt( px, 0)? 0 : sgn); + + if (FPgt( (sgn * py), 0)) { /* y and py have same sign */ + return(0); + + } else { /* y and py have opposite signs */ + if (FPge( x, 0) && FPgt( px, 0)) return(2 * sgn); + if (FPlt( x, 0) && FPle( px, 0)) return(0); + + z = (x-px) * y - (y-py) * x; + if (FPzero( z)) return(HIT_IT); + return( FPgt( (sgn*z), 0)? 0 : 2 * sgn); + }; +} /* lseg_crossing() */ + + +bool +plist_same(int npts, Point p1[], Point p2[]) +{ + int i, ii, j; + + /* find match for first point */ + for (i = 0; i < npts; i++) { + if ((FPeq( p2[i].x, p1[0].x)) + && (FPeq( p2[i].y, p1[0].y))) { + + /* match found? then look forward through remaining points */ + for (ii = 1, j = i+1; ii < npts; ii++, j++) { + if (j >= npts) j = 0; + if ((!FPeq( p2[j].x, p1[ii].x)) + || (!FPeq( p2[j].y, p1[ii].y))) { +#ifdef GEODEBUG +printf( "plist_same- %d failed forward match with %d\n", j, ii); +#endif + break; + }; + }; +#ifdef GEODEBUG +printf( "plist_same- ii = %d/%d after forward match\n", ii, npts); +#endif + if (ii == npts) + return(TRUE); + + /* match not found forwards? then look backwards */ + for (ii = 1, j = i-1; ii < npts; ii++, j--) { + if (j < 0) j = (npts-1); + if ((!FPeq( p2[j].x, p1[ii].x)) + || (!FPeq( p2[j].y, p1[ii].y))) { +#ifdef GEODEBUG +printf( "plist_same- %d failed reverse match with %d\n", j, ii); +#endif + break; + }; + }; +#ifdef GEODEBUG +printf( "plist_same- ii = %d/%d after reverse match\n", ii, npts); +#endif + if (ii == npts) + return(TRUE); + }; + }; + + return(FALSE); +} /* plist_same() */ + -- 2.40.0