Fix to arc center calculation, from tiipponen, #3099

git-svn-id: http://svn.osgeo.org/postgis/trunk@14425 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index 03faf110119450e6b2129c6a1c961730700e4d89..f3510ea2e9db0f72c9ace6b535fc02b67faa0ab0 100644 (file)
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -90,6 +90,29 @@ static void test_lw_segment_side(void)
 
 }
 
+static void test_lw_arc_center(void)
+{
+/* double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result); */
+       POINT2D c1;
+       double d1;
+       POINT2D p1, p2, p3;
+
+       p1.x = 2047538.600;
+       p1.y = 7268770.435;
+       p2.x = 2047538.598;
+       p2.y = 7268770.435;
+       p3.x = 2047538.596;
+       p3.y = 7268770.436;
+
+       d1 = lw_arc_center(&p1, &p2, &p3, &c1);
+
+       CU_ASSERT_DOUBLE_EQUAL(d1, 0.0046097720751, 0.0001);
+       CU_ASSERT_DOUBLE_EQUAL(c1.x, 2047538.599, 0.001);
+       CU_ASSERT_DOUBLE_EQUAL(c1.y, 7268770.4395, 0.001);
+       
+       // printf("lw_arc_center: (%12.12g, %12.12g) %12.12g\n", c1.x, c1.y, d1);
+}
+
 /*
 ** Test crossings side.
 */
@@ -983,4 +1006,5 @@ void algorithms_suite_setup(void)
        PG_ADD_TEST(suite,test_geohash_point_as_int);
        PG_ADD_TEST(suite,test_isclosed);
        PG_ADD_TEST(suite,test_lwgeom_simplify);
+       PG_ADD_TEST(suite,test_lw_arc_center);
 }

diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index 37a711536ab3243c21f50b2d6da29e5b4a0f6251..aa8fdfd171f03af5092bdf1998ae662e9634f109 100644 (file)
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -224,13 +224,13 @@ int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const P
 * point is coincident with either end point, they are taken as colinear.
 */
 double
-lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
+lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)        
 {
        POINT2D c;
        double cx, cy, cr;
-       double temp, bc, cd, det;
+       double dx21, dy21, dx31, dy31, h21, h31, d;
 
-    c.x = c.y = 0.0;
+       c.x = c.y = 0.0;
 
        LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
 
@@ -247,30 +247,35 @@ lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *
                return cr;
        }
 
-       temp = p2->x*p2->x + p2->y*p2->y;
-       bc = (p1->x*p1->x + p1->y*p1->y - temp) / 2.0;
-       cd = (temp - p3->x*p3->x - p3->y*p3->y) / 2.0;
-       det = (p1->x - p2->x)*(p2->y - p3->y)-(p2->x - p3->x)*(p1->y - p2->y);
+       /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
+       dx21 = p2->x - p1->x;
+       dy21 = p2->y - p1->y;
+       dx31 = p3->x - p1->x;
+       dy31 = p3->y - p1->y;
 
-       /* Check colinearity */
-       if (fabs(det) < EPSILON_SQLMM)
-               return -1.0;
+       h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
+       h31 = pow(dx31, 2.0) + pow(dy31, 2.0);
+
+       /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
+       d = 2 * (dx21 * dy31 - dx31 * dy21);
 
+       /* Check colinearity, |Cross product| = 0 */
+       if (fabs(d) < EPSILON_SQLMM)
+               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;
+       /* Calculate centroid coordinates and radius */
+       cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
+       cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
        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) + pow(cy - p1->y, 2));
+
        LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);
-       
+
        return cr;
 }
 
-
 int
 pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
 {