[postgis-commits] svn - r3079 - in branches/1.3: lwgeom regress
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Wed Oct 8 02:05:29 PDT 2008
Author: mleslie
Date: 2008-10-08 02:05:24 -0700 (Wed, 08 Oct 2008)
New Revision: 3079
Modified:
branches/1.3/lwgeom/lwcurve.c
branches/1.3/regress/sql-mm-circularstring.sql
branches/1.3/regress/sql-mm-circularstring_expected.in
Log:
Rebuilt the box3d generation for circular strings to account for a special large-arc case. Fix for issue 58, includes regression test for the case.
Modified: branches/1.3/lwgeom/lwcurve.c
===================================================================
--- branches/1.3/lwgeom/lwcurve.c 2008-10-08 08:59:29 UTC (rev 3078)
+++ branches/1.3/lwgeom/lwcurve.c 2008-10-08 09:05:24 UTC (rev 3079)
@@ -34,7 +34,7 @@
{
LWCURVE *result;
- /*
+ /*
* The first arc requires three points. Each additional
* arc requires two more points. Thus the minimum point count
* is three, and the count must be odd.
@@ -272,10 +272,10 @@
{
double x1, x2, y1, y2, z1, z2;
double angle, radius, sweep;
- /*
- double top, left;
- */
+ // 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;
POINT4D *center;
int i;
@@ -285,18 +285,16 @@
lwnotice("lwcircle_compute_box3d called.");
#endif
+ /*
+ * Determine where the center of the circle lives, and the radius.
+ */
center = lwalloc(sizeof(POINT4D));
radius = lwcircle_center(p1, p2, p3, ¢er);
if(radius < 0.0) return NULL;
- /*
- top = center->y + radius;
- left = center->x - radius;
-
#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcircle_compute_box3d: top=%.16f, left=%.16f", top, left);
+ lwnotice("lwcircle_compute_box3d: center (%.16f, %.16f)", center->x, center->y);
#endif
- */
x1 = MAXFLOAT;
x2 = -1 * MAXFLOAT;
@@ -307,30 +305,55 @@
a2 = atan2(p2->y - center->y, p2->x - center->x);
a3 = atan2(p3->y - center->y, p3->x - center->x);
- /* Determine sweep angle */
- if(a1 > a2 && a2 > a3)
+ // 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(r2 < 0 && (r3 > 0 || r3 < r2))
{
- sweep = a3 - a1;
+ sweep = (r3 > 0) ? (r3 - 2 * M_PI) : r3;
}
- /* Counter-clockwise */
- else if(a1 < a2 && a2 < a3)
+ // clockwise 2
+ else if(r2 > 0 && r3 > 0 && r3 < r2)
{
- sweep = a3 - a1;
+ sweep = (r3 > 0) ? (r3 - 2 * M_PI) : r3;
}
- /* Clockwise, wrap */
- else if((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3))
+ // counter-clockwise 4
+ else if(r2 > 0 && (r3 < 0 || r3 > r2))
{
- sweep = a3 - a1 + 2*M_PI;
+ sweep = (r3 < 0) ? (r3 + 2 * M_PI) : r3;
}
- /* Counter-clockwise, wrap */
- else if((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3))
+ // counter-clockwisk 5
+ else if(r2 < 0 && r3 < 0 && r3 > r2)
{
- sweep = a3 - a1 - 2*M_PI;
- }
- else
+ sweep = (r3 < 0) ? (r3 + 2 * M_PI) : r3;
+ }
+ // clockwise invalid 3
+ else if(r2 > 0 && (r3 > r2 || r3 < 0))
{
- sweep = 0.0;
+ sweep = (r2 > 0) ? (r2 - 2 * M_PI) : r2;
}
+ // clockwise invalid 6
+ else
+ {
+ sweep = (r2 < 0) ? (r2 + 2 * M_PI) : r2;
+ }
#ifdef PGIS_DEBUG
lwnotice("a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep);
@@ -340,41 +363,68 @@
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(sweep > 0.0 && (angle > a3 || angle < a1)) continue;
- if(sweep < 0.0 && (angle < a3 || angle > a1)) continue;
+ if(sweep > 0.0)
+ {
+ if(a3 < a1)
+ {
+ if(angle > (a3 + 2 * M_PI) || angle < a1) continue;
+ }
+ else
+ {
+ if(angle > a3 || angle < a1) continue;
+ }
+ }
+ else
+ {
+ if(a3 > a1)
+ {
+ if(angle < (a3 - 2 * M_PI) || angle > a1) continue;
+ }
+ else
+ {
+ if(angle < a3 || angle > a1) continue;
+ }
+ }
}
#ifdef PGIS_DEBUG
Modified: branches/1.3/regress/sql-mm-circularstring.sql
===================================================================
--- branches/1.3/regress/sql-mm-circularstring.sql 2008-10-08 08:59:29 UTC (rev 3078)
+++ branches/1.3/regress/sql-mm-circularstring.sql 2008-10-08 09:05:24 UTC (rev 3079)
@@ -384,3 +384,4 @@
SELECT DropGeometryColumn('public', 'circularstring', 'the_geom_3dm');
SELECT DropGeometryColumn('public', 'circularstring', 'the_geom_2d');
DROP TABLE public.circularstring;
+SELECT ST_asText(ST_box2d('CIRCULARSTRING(220268.439465645 150415.359530563,220227.333322076 150505.561285879,220227.353105332 150406.434743975)'::geometry));
Modified: branches/1.3/regress/sql-mm-circularstring_expected.in
===================================================================
--- branches/1.3/regress/sql-mm-circularstring_expected.in 2008-10-08 08:59:29 UTC (rev 3078)
+++ branches/1.3/regress/sql-mm-circularstring_expected.in 2008-10-08 09:05:24 UTC (rev 3079)
@@ -203,3 +203,4 @@
public.circularstring.the_geom_3dz effectively removed.
public.circularstring.the_geom_3dm effectively removed.
public.circularstring.the_geom_2d effectively removed.
+POLYGON((220187.375 150406.421875,220187.375 150506.71875,220288.828125 150506.71875,220288.828125 150406.421875,220187.375 150406.421875))
More information about the postgis-commits
mailing list