From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 28 Sep 2012 18:23:52 +0000 (+0000)
Subject: Measurement support for arcs (#2018)
X-Git-Tag: 2.1.0beta2~614
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=a59a8fb625387b6b31b7413b8e730b2782d75556;p=postgis

Measurement support for arcs (#2018)


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

diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c
index 963dcfb2a..fee2bdf43 100644
--- a/liblwgeom/cunit/cu_measures.c
+++ b/liblwgeom/cunit/cu_measures.c
@@ -542,6 +542,48 @@ test_lw_dist2d_arc_arc(void)
 	CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
 }
 
+static void
+test_lw_arc_length(void)
+{
+/* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
+
+	POINT2D A1, A2, A3;
+	double d;
+	
+	/* Unit semicircle at 0,0 */
+	A1.x = -1; A1.y = 0;
+	A2.x = 0 ; A2.y = 1;
+	A3.x = 1 ; A3.y = 0;
+
+	/* Arc above the unit semicircle */
+	d = lw_arc_length(&A1, &A2, &A3);
+	CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
+	d = lw_arc_length(&A3, &A2, &A1);
+	CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
+
+	/* Unit semicircle at 0,0 */
+	A1.x = 0; A1.y = 1;
+	A2.x = 1; A2.y = 0;
+	A3.x = 0; A3.y = -1;
+
+	/* Arc to right of the unit semicircle */
+	d = lw_arc_length(&A1, &A2, &A3);
+	CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
+	d = lw_arc_length(&A3, &A2, &A1);
+	CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
+
+	/* Unit 3/4 circle at 0,0 */
+	A1.x = -1; A1.y = 0;
+	A2.x = 1; A2.y = 0;
+	A3.x = 0; A3.y = -1;
+
+	/* Arc to right of the unit semicircle */
+	d = lw_arc_length(&A1, &A2, &A3);
+	CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);
+	d = lw_arc_length(&A3, &A2, &A1);
+	CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);	
+}
+
 /*
 ** Used by test harness to register the tests in this file.
 */
@@ -555,6 +597,7 @@ CU_TestInfo measures_tests[] =
 	PG_TEST(test_lw_dist2d_pt_arc),
 	PG_TEST(test_lw_dist2d_seg_arc),
 	PG_TEST(test_lw_dist2d_arc_arc),
+	PG_TEST(test_lw_arc_length),
 	CU_TEST_INFO_NULL
 };
 CU_SuiteInfo measures_suite = {"PostGIS Measures Suite",  NULL,  NULL, measures_tests};
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index fee34f628..85d03057a 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -860,6 +860,69 @@ lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
 		return LW_FALSE;
 }
 
+/**
+* Returns the length of a circular arc segment
+*/
+double
+lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
+{
+	POINT2D C;
+	double radius_A, circumference_A;
+	int a2_side, clockwise;
+	double a1, a2, a3;
+	double angle;
+	
+	if ( lw_arc_is_pt(A1, A2, A3) )
+		return 0.0;
+	
+	radius_A = lwcircle_center(A1, A2, A3, &C);
+
+	/* Co-linear! Return linear distance! */
+	if ( radius_A < 0 ) 
+	{
+		return sqrt((A1->x-A3->x)*(A1->x-A3->x) + (A1->y-A3->y)*(A1->y-A3->y));
+	}
+	
+	/* Closed circle! Return the circumference! */
+	circumference_A = M_PI * 2 * radius_A;
+	if ( p2d_same(A1, A3) )
+		return circumference_A;
+	
+	/* Determine the orientation of the arc */
+	a2_side = signum(lw_segment_side(A1, A3, A2));
+
+	/* The side of the A1/A3 line that A2 falls on dictates the sweep  
+	   direction from A1 to A3. */
+	if ( a2_side == -1 ) 
+		clockwise = LW_TRUE;
+	else 
+		clockwise = LW_FALSE;
+		
+	/* Angles of each point that defines the arc section */
+	a1 = atan2(A1->y - C.y, A1->x - C.x);
+	a2 = atan2(A2->y - C.y, A2->x - C.x);
+	a3 = atan2(A3->y - C.y, A3->x - C.x);
+
+	/* What's the sweep from A1 to A3? */
+	if ( clockwise )
+	{
+		if ( a1 > a3 )
+			angle = a1 - a3;
+		else
+			angle = 2*M_PI + a1 - a3; 
+	}
+	else
+	{
+		if ( a3 > a1 )
+			angle = a3 - a1;
+		else
+			angle = 2*M_PI + a3 - a1; 			
+	}
+
+	/* Length as proportion of circumference */
+	return circumference_A* (angle / (2*M_PI));
+}
+
 /**
 * Calculate the shortest distance between an arc and an edge.
 * Line/circle approach from http://stackoverflow.com/questions/1073336/circle-line-collision-detection 
diff --git a/liblwgeom/measures.h b/liblwgeom/measures.h
index a705f15c9..8cd3d9bd7 100644
--- a/liblwgeom/measures.h
+++ b/liblwgeom/measures.h
@@ -76,3 +76,8 @@ int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, c
 int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D* B3, DISTPTS *dl);
 void lw_dist2d_distpts_init(DISTPTS *dl, int mode);
 
+/*
+* Length primitives
+*/
+double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
+