From 83b6ddaa76847c45eebe0229035dac276a0cc373 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 23 Jun 2011 13:57:16 +0000 Subject: [PATCH] PostGIS Box2D (and && operator) gives wrong result for ST_CircularString type (#578) git-svn-id: http://svn.osgeo.org/postgis/trunk@7452 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/g_box.c | 245 +++++++++++---------------------------- liblwgeom/liblwgeom.h | 2 +- liblwgeom/lwsegmentize.c | 13 ++- 3 files changed, 77 insertions(+), 183 deletions(-) diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index 6e87af17b..85292c246 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -10,6 +10,7 @@ **********************************************************************/ #include "liblwgeom_internal.h" +#include "lwalgorithm.h" #include #include @@ -336,202 +337,92 @@ size_t gbox_serialized_size(uchar flags) } +static inline int signum(double n) +{ + if( n < 0 ) return -1; + if( n > 0 ) return 1; + return 0; +} /* ******************************************************************************** ** Compute cartesian bounding GBOX boxes from LWGEOM. */ -static int lwcircle_calculate_gbox_cartesian(POINT4D p1, POINT4D p2, POINT4D p3, GBOX *gbox) +static int lwcircle_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox) { - double x1, x2, y1, y2, z1, z2, m1, m2; - double angle, radius, sweep; - /* angles from center */ - double a1, a2, a3; - /* angles from center once a1 is rotated to zero */ - double r2, r3; - double xe = 0.0, ye = 0.0; + POINT2D xmin, ymin, xmax, ymax; POINT4D center; - int i; + int p2_side; + double radius; LWDEBUG(2, "lwcircle_calculate_gbox called."); - radius = lwcircle_center(&p1, &p2, &p3, ¢er); - if (radius < 0.0) /* p1,p2,p3 are colinear */ + radius = lwcircle_center(p1, p2, p3, ¢er); + + /* Negative radius signals straight line, p1/p2/p3 are colinear */ + if (radius < 0.0) { - gbox->xmin = FP_MIN(p1.x, p3.x); - gbox->ymin = FP_MIN(p1.y, p3.y); - gbox->zmin = FP_MIN(p1.z, p3.z); - gbox->xmax = FP_MAX(p1.x, p3.x); - gbox->ymax = FP_MAX(p1.y, p3.y); - gbox->zmax = FP_MAX(p1.z, p3.z); + gbox->xmin = FP_MIN(p1->x, p3->x); + gbox->ymin = FP_MIN(p1->y, p3->y); + gbox->zmin = FP_MIN(p1->z, p3->z); + gbox->xmax = FP_MAX(p1->x, p3->x); + gbox->ymax = FP_MAX(p1->y, p3->y); + gbox->zmax = FP_MAX(p1->z, p3->z); return LW_SUCCESS; } - - x1 = MAXFLOAT; - x2 = -1 * MAXFLOAT; - y1 = MAXFLOAT; - y2 = -1 * MAXFLOAT; - - a1 = atan2(p1.y - center.y, p1.x - center.x); - a2 = atan2(p2.y - center.y, p2.x - center.x); - a3 = atan2(p3.y - center.y, p3.x - center.x); - - /* Rotate a2 and a3 such that a1 = 0 */ - r2 = a2 - a1; - r3 = a3 - a1; - - /* - * There are six cases here I'm interested in - * Clockwise: - * 1. a1-a2 < 180 == r2 < 0 && (r3 > 0 || r3 < r2) - * 2. a1-a2 > 180 == r2 > 0 && (r3 > 0 && r3 < r2) - * 3. a1-a2 > 180 == r2 > 0 && (r3 > r2 || r3 < 0) - * Counter-clockwise: - * 4. a1-a2 < 180 == r2 > 0 && (r3 < 0 || r3 > r2) - * 5. a1-a2 > 180 == r2 < 0 && (r3 < 0 && r3 > r2) - * 6. a1-a2 > 180 == r2 < 0 && (r3 < r2 || r3 > 0) - * 3 and 6 are invalid cases where a3 is the midpoint. - * BBOX is fundamental, so these cannot error out and will instead - * calculate the sweep using a3 as the middle point. - */ - - /* clockwise 1 */ - if (FP_LT(r2, 0) && (FP_GT(r3, 0) || FP_LT(r3, r2))) - { - sweep = (FP_GT(r3, 0)) ? (r3 - 2 * M_PI) : r3; - } - /* clockwise 2 */ - else if (FP_GT(r2, 0) && FP_GT(r3, 0) && FP_LT(r3, r2)) - { - sweep = (FP_GT(r3, 0)) ? (r3 - 2 * M_PI) : r3; - } - /* counter-clockwise 4 */ - else if (FP_GT(r2, 0) && (FP_LT(r3, 0) || FP_GT(r3, r2))) - { - sweep = (FP_LT(r3, 0)) ? (r3 + 2 * M_PI) : r3; - } - /* counter-clockwise 5 */ - else if (FP_LT(r2, 0) && FP_LT(r3, 0) && FP_GT(r3, r2)) - { - sweep = (FP_LT(r3, 0)) ? (r3 + 2 * M_PI) : r3; - } - /* clockwise invalid 3 */ - else if (FP_GT(r2, 0) && (FP_GT(r3, r2) || FP_LT(r3, 0))) - { - sweep = (FP_GT(r2, 0)) ? (r2 - 2 * M_PI) : r2; - } - /* clockwise invalid 6 */ - else + + /* Matched start/end points imply circle */ + if ( p1->x == p3->x && p1->y == p3->y ) { - sweep = (FP_LT(r2, 0)) ? (r2 + 2 * M_PI) : r2; + gbox->xmin = center.x - radius; + gbox->ymin = center.y - radius; + gbox->zmin = FP_MIN(p1->z,p2->z); + gbox->mmin = FP_MIN(p1->m,p2->m); + gbox->xmax = center.x + radius; + gbox->ymax = center.y + radius; + gbox->zmax = FP_MAX(p1->z,p2->z); + gbox->mmax = FP_MAX(p1->m,p2->m); + return LW_SUCCESS; } - LWDEBUGF(3, "a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep); - - angle = 0.0; - for (i=0; i < 6; i++) - { - switch (i) - { - /* right extent */ - case 0: - angle = 0.0; - xe = center.x + radius; - ye = center.y; - break; - /* top extent */ - case 1: - angle = M_PI_2; - xe = center.x; - ye = center.y + radius; - break; - /* left extent */ - case 2: - angle = M_PI; - xe = center.x - radius; - ye = center.y; - break; - /* bottom extent */ - case 3: - angle = -1 * M_PI_2; - xe = center.x; - ye = center.y - radius; - break; - /* first point */ - case 4: - angle = a1; - xe = p1.x; - ye = p1.y; - break; - /* last point */ - case 5: - angle = a3; - xe = p3.x; - ye = p3.y; - break; - } - /* determine if the extents are outside the arc */ - if (i < 4) - { - if (FP_GT(sweep, 0.0)) - { - if (FP_LT(a3, a1)) - { - if (FP_GT(angle, (a3 + 2 * M_PI)) || FP_LT(angle, a1)) continue; - } - else - { - if (FP_GT(angle, a3) || FP_LT(angle, a1)) continue; - } - } - else - { - if (FP_GT(a3, a1)) - { - if (FP_LT(angle, (a3 - 2 * M_PI)) || FP_GT(angle, a1)) continue; - } - else - { - if (FP_LT(angle, a3) || FP_GT(angle, a1)) continue; - } - } - } - - LWDEBUGF(3, "lwcircle_calculate_gbox: potential extreame %d (%.16f, %.16f)", i, xe, ye); - - x1 = (FP_LT(x1, xe)) ? x1 : xe; - y1 = (FP_LT(y1, ye)) ? y1 : ye; - x2 = (FP_GT(x2, xe)) ? x2 : xe; - y2 = (FP_GT(y2, ye)) ? y2 : ye; - } - - LWDEBUGF(3, "lwcircle_calculate_gbox: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2); + /* First approximation, bounds of start/end points */ + gbox->xmin = FP_MIN(p1->x, p3->x); + gbox->ymin = FP_MIN(p1->y, p3->y); + gbox->zmin = FP_MIN(p1->z, p3->z); + gbox->mmin = FP_MIN(p1->m, p3->m); + gbox->xmax = FP_MAX(p1->x, p3->x); + gbox->ymax = FP_MAX(p1->y, p3->y); + gbox->zmax = FP_MAX(p1->z, p3->z); + gbox->mmax = FP_MAX(p1->m, p3->m); + + /* Create points for the possible extrema */ + xmin.x = center.x - radius; + xmin.y = center.y; + ymin.x = center.x; + ymin.y = center.y - radius; + xmax.x = center.x + radius; + xmax.y = center.y; + ymax.x = center.x; + ymax.y = center.y + radius; + + + /* Divide the circle into two parts, one on each side of a line + joining p1 and p3. The circle extrema on the same side of that line + as p2 is on, are also the extrema of the bbox. */ + + p2_side = signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2)); - z1 = FP_MIN(p1.z, p2.z); - z1 = FP_MIN(z1, p3.z); - z2 = FP_MAX(p1.z, p2.z); - z2 = FP_MAX(z2, p3.z); + if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &xmin)) ) + gbox->xmin = xmin.x; - m1 = FP_MIN(p1.m, p2.m); - m1 = FP_MIN(m1, p3.m); - m2 = FP_MAX(p1.m, p2.m); - m2 = FP_MAX(m2, p3.m); + if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &ymin)) ) + gbox->ymin = ymin.y; - gbox->xmin = x1; - gbox->xmax = x2; - gbox->ymin = y1; - gbox->ymax = y2; + if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &xmax)) ) + gbox->xmax = xmax.x; - if ( FLAGS_GET_Z(gbox->flags) ) - { - gbox->zmin = z1; - gbox->zmax = z2; - } - if ( FLAGS_GET_M(gbox->flags) ) - { - gbox->mmin = m1; - gbox->mmax = m2; - } + if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &ymax)) ) + gbox->ymax = ymax.y; return LW_SUCCESS; } @@ -602,7 +493,7 @@ static int lwcircstring_calculate_gbox_cartesian(LWCIRCSTRING *curve, GBOX *gbox getPoint4d_p(curve->points, i-1, &p2); getPoint4d_p(curve->points, i, &p3); - if (lwcircle_calculate_gbox_cartesian(p1, p2, p3, &tmp) == LW_FAILURE) + if (lwcircle_calculate_gbox_cartesian(&p1, &p2, &p3, &tmp) == LW_FAILURE) continue; gbox_merge(&tmp, gbox); diff --git a/liblwgeom/liblwgeom.h b/liblwgeom/liblwgeom.h index ae01abba4..cf53d3a2f 100644 --- a/liblwgeom/liblwgeom.h +++ b/liblwgeom/liblwgeom.h @@ -2186,7 +2186,7 @@ void error_if_srid_mismatch(int srid1, int srid2); ******************************************************************************/ int lwgeom_has_arc(const LWGEOM *geom); -double lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result); +double lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result); LWGEOM *lwgeom_segmentize(LWGEOM *geom, uint32 perQuad); LWGEOM *lwgeom_desegmentize(LWGEOM *geom); diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c index 6e2a41990..2affdf3a3 100644 --- a/liblwgeom/lwsegmentize.c +++ b/liblwgeom/lwsegmentize.c @@ -85,7 +85,7 @@ lwgeom_has_arc(const LWGEOM *geom) * point is coincident with either end point, they are taken as colinear. */ double -lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result) +lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result) { POINT4D c; double cx, cy, cr; @@ -93,18 +93,18 @@ lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result) c.x = c.y = c.z = c.m = 0.0; - LWDEBUGF(2, "lwcircle_center called (%.16f, %.16f), (%.16f, %.16f), (%.16f, %.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y); + LWDEBUGF(2, "lwcircle_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y); /* Closed circle */ - if (fabs(p1->x - p3->x) < EPSILON_SQLMM - && fabs(p1->y - p3->y) < EPSILON_SQLMM) + if (fabs(p1->x - p3->x) < EPSILON_SQLMM && + fabs(p1->y - p3->y) < EPSILON_SQLMM) { cx = p1->x + (p2->x - p1->x) / 2.0; cy = p1->y + (p2->y - p1->y) / 2.0; c.x = cx; c.y = cy; *result = c; - cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y)); + cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0)); return cr; } @@ -125,6 +125,9 @@ lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result) c.y = cy; *result = c; cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y)); + + LWDEBUGF(2, "lwcircle_center center is (%.16f,%.16f)", result->x, result->y); + return cr; } -- 2.50.1