[postgis-tickets] r14425 - Fix to arc center calculation, from tiipponen, #3099

[Paul Ramsey]

Author: pramsey
Date: 2015-11-25 12:10:09 -0800 (Wed, 25 Nov 2015)
New Revision: 14425

Modified:
   trunk/liblwgeom/cunit/cu_algorithm.c
   trunk/liblwgeom/lwalgorithm.c
Log:
Fix to arc center calculation, from tiipponen, #3099


Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c	2015-11-25 18:33:54 UTC (rev 14424)
+++ trunk/liblwgeom/cunit/cu_algorithm.c	2015-11-25 20:10:09 UTC (rev 14425)
@@ -90,6 +90,29 @@
 
 }
 
+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 @@
 	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);
 }

Modified: trunk/liblwgeom/lwalgorithm.c
===================================================================
--- trunk/liblwgeom/lwalgorithm.c	2015-11-25 18:33:54 UTC (rev 14424)
+++ trunk/liblwgeom/lwalgorithm.c	2015-11-25 20:10:09 UTC (rev 14425)
@@ -224,13 +224,13 @@
 * 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 @@
 		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)
+	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)
 {