From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Wed, 25 Nov 2015 20:10:09 +0000 (+0000)
Subject: Fix to arc center calculation, from tiipponen, #3099
X-Git-Tag: 2.3.0beta1~367
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=b695c01409f607c73e6c39031476032a05e73f6a;p=postgis

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 03faf1101..f3510ea2e 100644
--- 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 37a711536..aa8fdfd17 100644
--- 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)
 {