From 9665523267ac99f096073730d8c803fe89db47f3 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 11 Oct 2012 17:29:55 +0000 Subject: [PATCH] Move area core calculation to ptarray. git-svn-id: http://svn.osgeo.org/postgis/trunk@10404 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_ptarray.c | 30 +++++++++++++++++++ liblwgeom/liblwgeom_internal.h | 10 ++++++- liblwgeom/lwpoly.c | 55 ++++++++++------------------------ liblwgeom/measures.c | 1 + liblwgeom/ptarray.c | 44 ++++++++++++++++++++++++--- 5 files changed, 96 insertions(+), 44 deletions(-) diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c index 4ad95f622..80d3004d5 100644 --- 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 240580904..4b45f3c54 100644 --- 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 105fcc3d5..119bc61c3 100644 --- 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; inrings; 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; jnpoints-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 47aa5eb47..2bc52dcab 100644 --- 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 86c890403..d129d03a3 100644 --- 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) { -- 2.40.0