Skip to content

Commit

Permalink
Bounding box calculation for straight CircString fails (postgis#668)
Browse files Browse the repository at this point in the history 
git-svn-id: http://svn.osgeo.org/postgis/trunk@7416 b70326c6-7e19-0410-871a-916f4a2858ee
  • Loading branch information
@pramsey
pramsey committed Jun 17, 2011
1 parent 96d94fb commit c8d72b7
Show file tree
Hide file tree
Showing 6 changed files with 73 additions and 59 deletions.
35 changes: 22 additions & 13 deletions liblwgeom/g_box.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,22 +347,31 @@ static int lwcircle_calculate_gbox_cartesian(POINT4D p1, POINT4D p2, POINT4D p3,
/* angles from center once a1 is rotated to zero */
double r2, r3;
double xe = 0.0, ye = 0.0;
POINT4D *center;
POINT4D center;
int i;

LWDEBUG(2, "lwcircle_calculate_gbox called.");

radius = lwcircle_center(&p1, &p2, &p3, &center);
if (radius < 0.0) return LW_FAILURE;
if (radius < 0.0) /* p1,p2,p3 are colinear */
{
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);
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;
Expand Down Expand Up @@ -424,26 +433,26 @@ static int lwcircle_calculate_gbox_cartesian(POINT4D p1, POINT4D p2, POINT4D p3,
/* right extent */
case 0:
angle = 0.0;
xe = center->x + radius;
ye = center->y;
xe = center.x + radius;
ye = center.y;
break;
/* top extent */
case 1:
angle = M_PI_2;
xe = center->x;
ye = center->y + radius;
xe = center.x;
ye = center.y + radius;
break;
/* left extent */
case 2:
angle = M_PI;
xe = center->x - radius;
ye = center->y;
xe = center.x - radius;
ye = center.y;
break;
/* bottom extent */
case 3:
angle = -1 * M_PI_2;
xe = center->x;
ye = center->y - radius;
xe = center.x;
ye = center.y - radius;
break;
/* first point */
case 4:
Expand Down
3 changes: 1 addition & 2 deletions liblwgeom/liblwgeom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2186,10 +2186,9 @@ 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(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result);
LWGEOM *lwgeom_segmentize(LWGEOM *geom, uint32 perQuad);
LWGEOM *lwgeom_desegmentize(LWGEOM *geom);
double lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result);

#endif /* !defined _LIBLWGEOM_H */

54 changes: 30 additions & 24 deletions liblwgeom/lwcircstring.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,31 +290,40 @@ lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)
/* angles from center once a1 is rotated to zero */
double r2, r3;
double xe = 0.0, ye = 0.0;
POINT4D *center;
POINT4D center;
int i;
BOX3D *box;

LWDEBUG(2, "lwcircle_compute_box3d called.");

center = lwalloc(sizeof(POINT4D));
box = lwalloc(sizeof(BOX3D));
radius = lwcircle_center(p1, p2, p3, &center);
if (radius < 0.0) return NULL;
if (radius < 0.0)
{
box->xmin = FP_MIN(p1->x, p3->x);
box->ymin = FP_MIN(p1->y, p3->y);
box->zmin = FP_MIN(p1->z, p3->z);
box->xmax = FP_MAX(p1->x, p3->x);
box->ymax = FP_MAX(p1->y, p3->y);
box->zmax = FP_MAX(p1->z, p3->z);
return box;
}

/*
top = center->y + radius;
left = center->x - radius;
top = center.y + radius;
left = center.x - radius;
LWDEBUGF(3, "lwcircle_compute_box3d: center (%.16f, %.16f)", center->x, center->y);
LWDEBUGF(3, "lwcircle_compute_box3d: center (%.16f, %.16f)", center.x, center.y);
*/

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);
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;
Expand Down Expand Up @@ -376,26 +385,26 @@ lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)
/* right extent */
case 0:
angle = 0.0;
xe = center->x + radius;
ye = center->y;
xe = center.x + radius;
ye = center.y;
break;
/* top extent */
case 1:
angle = M_PI_2;
xe = center->x;
ye = center->y + radius;
xe = center.x;
ye = center.y + radius;
break;
/* left extent */
case 2:
angle = M_PI;
xe = center->x - radius;
ye = center->y;
xe = center.x - radius;
ye = center.y;
break;
/* bottom extent */
case 3:
angle = -1 * M_PI_2;
xe = center->x;
ye = center->y - radius;
xe = center.x;
ye = center.y - radius;
break;
/* first point */
case 4:
Expand Down Expand Up @@ -448,26 +457,23 @@ lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)
LWDEBUGF(3, "lwcircle_compute_box3d: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2);

/*
x1 = center->x + x1 * radius;
x2 = center->x + x2 * radius;
y1 = center->y + y1 * radius;
y2 = center->y + y2 * radius;
x1 = center.x + x1 * radius;
x2 = center.x + x2 * radius;
y1 = center.y + y1 * radius;
y2 = center.y + y2 * radius;
*/
z1 = (FP_LT(p1->z, p2->z)) ? p1->z : p2->z;
z1 = (FP_LT(z1, p3->z)) ? z1 : p3->z;
z2 = (FP_GT(p1->z, p2->z)) ? p1->z : p2->z;
z2 = (FP_GT(z2, p3->z)) ? z2 : p3->z;

box = lwalloc(sizeof(BOX3D));
box->xmin = x1;
box->xmax = x2;
box->ymin = y1;
box->ymax = y2;
box->zmin = z1;
box->zmax = z2;

lwfree(center);

return box;
}

Expand Down
36 changes: 16 additions & 20 deletions liblwgeom/lwsegmentize.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,12 +85,14 @@ 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(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *result)
{
POINT4D *c;
POINT4D c;
double cx, cy, cr;
double temp, bc, cd, det;

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);

/* Closed circle */
Expand All @@ -99,9 +101,8 @@ lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result)
{
cx = p1->x + (p2->x - p1->x) / 2.0;
cy = p1->y + (p2->y - p1->y) / 2.0;
c = lwalloc(sizeof(POINT2D));
c->x = cx;
c->y = cy;
c.x = cx;
c.y = cy;
*result = c;
cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
return cr;
Expand All @@ -114,17 +115,14 @@ lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result)

/* Check colinearity */
if (fabs(det) < EPSILON_SQLMM)
{
*result = NULL;
return -1.0;
}


det = 1.0 / det;
cx = (bc*(p2->y - p3->y)-cd*(p1->y - p2->y))*det;
cy = ((p1->x - p2->x)*cd-(p2->x - p3->x)*bc)*det;
c = lwalloc(sizeof(POINT4D));
c->x = cx;
c->y = cy;
c.x = cx;
c.y = cy;
*result = c;
cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
return cr;
Expand Down Expand Up @@ -156,7 +154,7 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad)
uint32 ptcount;
uchar *pt;

POINT4D *center;
POINT4D center;
double radius = 0.0,
sweep = 0.0,
angle = 0.0,
Expand All @@ -174,11 +172,11 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad)
return NULL;
}

LWDEBUGF(3, "lwcircle_segmentize, (%.16f, %.16f) radius=%.16f", center->x, center->y, radius);
LWDEBUGF(3, "lwcircle_segmentize, (%.16f, %.16f) radius=%.16f", center.x, center.y, radius);

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);
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);

LWDEBUGF(3, "a1 = %.16f, a2 = %.16f, a3 = %.16f", a1, a2, a3);

Expand Down Expand Up @@ -228,8 +226,8 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad)
angle += increment;
if (increment > 0.0 && angle > M_PI) angle -= 2*M_PI;
if (increment < 0.0 && angle < -1*M_PI) angle -= 2*M_PI;
pbuf.x = center->x + radius*cos(angle);
pbuf.y = center->y + radius*sin(angle);
pbuf.x = center.x + radius*cos(angle);
pbuf.y = center.y + radius*sin(angle);
if ((sweep > 0 && angle < a2) || (sweep < 0 && angle > a2))
{
if ((sweep > 0 && a1 < a2) || (sweep < 0 && a1 > a2))
Expand Down Expand Up @@ -278,8 +276,6 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad)
pt = getPoint_internal(result, ptcount - 1);
memcpy(pt, (uchar *)p3, ptsize);

lwfree(center);

return result;
}

Expand Down
3 changes: 3 additions & 0 deletions regress/tickets.sql
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,5 +384,8 @@ select '#938', 'POLYGON EMPTY'::geometry::box2d;

DROP TABLE foo;

-- #668 --
select '#668',box2d('CIRCULARSTRING(10 2,12 2,14 2)'::geometry) as b;

-- Clean up
DELETE FROM spatial_ref_sys;
1 change: 1 addition & 0 deletions regress/tickets_expected
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,3 +111,4 @@ ERROR: First argument must be a LINESTRING
#884|1|f
#884|2|t
#938|
#668|BOX(10 2,14 2)

0 comments on commit c8d72b7

Please sign in to comment.