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 963dcfb2aaeba34d96d5d2da3fa505a94a8b42c2..fee2bdf43239e3db13882322fd2236ddc9aa714f 100644 (file)
--- 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 fee34f62832ab4d869b74e2aab4da3f95d8cdc0e..85d03057a5b3057d81c05b2fc87b6d0fdbe81137 100644 (file)
--- 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 a705f15c98598ce8ee31332ba46807b26b979774..8cd3d9bd71ed5541d444317e06e06388e7db8d7c 100644 (file)
--- 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);
+