Move area core calculation to ptarray.

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

diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 4ad95f6227d657a0a4a8aab0541cf8e1dbe72f11..80d3004d5de028706144004ef83750286f04afbc 100644 (file)
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -313,6 +313,34 @@ static void test_ptarray_isccw(void)
        lwpoly_free(poly);
 }
 
+static void test_ptarray_signed_area() 
+{
+       LWLINE *line;
+       double area;
+
+       /* parallelogram */
+       line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,1 1, 2 1, 1 0, 0 0)"));
+       area = ptarray_signed_area(line->points);
+       CU_ASSERT_DOUBLE_EQUAL(area, 1.0, 0.0000001);
+       lwline_free(line);
+
+       /* square */
+       line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,0 2, 2 2, 2 0, 0 0)"));
+       area = ptarray_signed_area(line->points);
+       CU_ASSERT_DOUBLE_EQUAL(area, 4.0, 0.0000001);
+       lwline_free(line);
+
+       /* square backwares*/
+       line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,2 0, 2 2, 0 2, 0 0)"));
+       area = ptarray_signed_area(line->points);
+       //printf("%g\n",area);
+       CU_ASSERT_DOUBLE_EQUAL(area, -4.0, 0.0000001);
+       lwline_free(line);
+       
+}
+
+
+
 static void test_ptarray_desegmentize() 
 {
        LWGEOM *in, *out;
@@ -580,6 +608,7 @@ static void test_ptarray_contains_point_arc()
        lwline_free(lwline);
 }
 
+
 /*
 ** Used by the test harness to register the tests in this file.
 */
@@ -589,6 +618,7 @@ CU_TestInfo ptarray_tests[] =
        PG_TEST(test_ptarray_append_ptarray),
        PG_TEST(test_ptarray_locate_point),
        PG_TEST(test_ptarray_isccw),
+       PG_TEST(test_ptarray_signed_area),
        PG_TEST(test_ptarray_desegmentize),
        PG_TEST(test_ptarray_insert_point),
        PG_TEST(test_ptarray_contains_point),

diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index 2405809043ca7fdae81218efe41983b408a7f145..4b45f3c540805ab3511a195a1a1e509aa70156a3 100644 (file)
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -121,11 +121,18 @@
 #define SIZE_GET(varsize) (((varsize) >> 2) & 0x3FFFFFFF)
 #define SIZE_SET(varsize, size) (((varsize) & 0x00000003)|(((size) & 0x3FFFFFFF) << 2 ))
 
-/*
+/**
 * Tolerance used to determine equality.
 */
 #define EPSILON_SQLMM 1e-8
 
+/**
+* Constants for point-in-polygon return values
+*/
+#define LW_INSIDE 1
+#define LW_BOUNDARY 0
+#define LW_OUTSIDE -1
+
 /*
 * Internal prototypes
 */
@@ -299,6 +306,7 @@ void ptarray_affine(POINTARRAY *pa, const AFFINE *affine);
 char ptarray_isccw(const POINTARRAY *pa);
 int ptarray_has_z(const POINTARRAY *pa);
 int ptarray_has_m(const POINTARRAY *pa);
+double ptarray_signed_area(const POINTARRAY *pa);
 
 /*
 * Clone support

diff --git a/liblwgeom/lwpoly.c b/liblwgeom/lwpoly.c
index 105fcc3d5fd15fde40b6fffd9de103cc2babee7a..119bc61c375ae642e037c9c977a1695630af9ee2 100644 (file)
--- a/liblwgeom/lwpoly.c
+++ b/liblwgeom/lwpoly.c
@@ -384,54 +384,31 @@ LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist)
 }
 
 /**
- * Find the area of the outer ring - sum (area of inner rings).
- * Could use a more numerically stable calculator...
- */
+* Find the area of the outer ring - sum (area of inner rings).
+*/
 double
 lwpoly_area(const LWPOLY *poly)
 {
-       double poly_area=0.0;
+       double poly_area = 0.0;
        int i;
-       POINT2D pp;
-       POINT2D cp;
-       POINT2D np;
-        double x0 = cp.x;
-
-       LWDEBUGF(2, "in lwpoly_area (%d rings)", poly->nrings);
+       
+       if ( ! poly ) 
+               lwerror("lwpoly_area called with null polygon pointer!");
 
-       for (i=0; i<poly->nrings; i++)
+       for ( i=0; i < poly->nrings; i++ )
        {
-               int j;
                POINTARRAY *ring = poly->rings[i];
                double ringarea = 0.0;
 
-               LWDEBUGF(4, " rings %d has %d points", i, ring->npoints);
-
-               if ( ! ring->npoints ) continue; /* empty ring */
-
-               getPoint2d_p(ring, 0, &cp);
-               getPoint2d_p(ring, 1, &np);
-                x0 = cp.x;
-                np.x -= x0;
-                for (j=0; j<ring->npoints-1; j++)
-               {
-                        pp.y = cp.y;
-                        cp.x = np.x;
-                        cp.y = np.y;
-                       getPoint2d_p(ring, j+1, &np);
-                        np.x -= x0;
-                        ringarea += cp.x * (np.y - pp.y);
-               }
-
-               ringarea  /= 2.0;
-
-               LWDEBUGF(4, " ring 1 has area %lf",ringarea);
-
-               ringarea  = fabs(ringarea);
-               if (i != 0)     /*outer */
-                       ringarea  = -1.0*ringarea ; /* its a hole */
-
-               poly_area += ringarea;
+               /* Empty or messed-up ring. */
+               if ( ring->npoints < 3 ) 
+                       continue; 
+               
+               ringarea = fabs(ptarray_signed_area(ring));
+               if ( i == 0 ) /* Outer ring, positive area! */
+                       poly_area += ringarea; 
+               else /* Inner ring, negative area! */
+                       poly_area -= ringarea; 
        }
 
        return poly_area;

diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index 47aa5eb47f9e2c48306b3efd44663423f6bd9573..2bc52dcabd87df38b9757b79c0c5d27e8e091100 100644 (file)
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -977,6 +977,7 @@ lw_dist2d_ptarrayarc_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DIST
        }
        return LW_TRUE;
 }
+
 /**
 * 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/ptarray.c b/liblwgeom/ptarray.c
index 86c8904032f8debd5983bdbb294766aad15b82bb..d129d03a344dc0b61e07e66a7a071828cb5b97d3 100644 (file)
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -795,14 +795,14 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt)
        if ( (pa->npoints % 2) == 0 )
        {
                lwerror("ptarray_contains_point_arc called with even number of points");
-               return -1;
+               return LW_OUTSIDE;
        }
 
        /* Check for not an arc ring (always have >= 3 points) */
        if ( pa->npoints < 3 )
        {
                lwerror("ptarray_contains_point_arc called too-short pointarray");
-               return -1;
+               return LW_OUTSIDE;
        }
 
        /* Check for unclosed case */
@@ -811,7 +811,7 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt)
        if ( ! p2d_same(seg1, seg3) )
        {
                lwerror("ptarray_contains_point_arc called on unclosed ring");
-               return -1;
+               return LW_OUTSIDE;
        } 
        /* OK, it's closed. Is it just one circle? */
        else if ( pa->npoints == 3 )
@@ -832,7 +832,7 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt)
                else if ( d < radius ) 
                        return 1; /* Inside circle */
                else 
-                       return -1; /* Outside circle */
+                       return LW_OUTSIDE; /* Outside circle */
        } 
        else if ( p2d_same(seg1, pt) )
        {
@@ -921,6 +921,42 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt)
        return 1;
 }
 
+/**
+* Returns the area in cartesian units. Area is negative if ring is oriented CCW, 
+* positive if it is oriented CW and zero if the ring is degenerate or flat.
+* http://en.wikipedia.org/wiki/Shoelace_formula
+*/
+double
+ptarray_signed_area(const POINTARRAY *pa)
+{
+       const POINT2D *P1;
+       const POINT2D *P2;
+       const POINT2D *P3;
+       double sum = 0.0;
+       double x0, x, y1, y2;
+       int i;
+       
+       if (! pa || pa->npoints < 3 )
+               return 0.0;
+               
+       P1 = getPoint2d_cp(pa, 0);
+       P2 = getPoint2d_cp(pa, 1);
+       x0 = P1->x;
+       for ( i = 1; i < pa->npoints - 1; i++ )
+       {
+               P3 = getPoint2d_cp(pa, i+1);
+               x = P2->x - x0;
+               y1 = P3->y;
+               y2 = P1->y;
+               sum += x * (y2-y1);
+               
+               /* Move forwards! */
+               P1 = P2;
+               P2 = P3;
+       }
+       return sum / 2.0;       
+}
+
 POINTARRAY*
 ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)
 {