From 9baacd722f96e29b79615168bfcef9200a546269 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 28 Mar 2013 22:11:09 +0000 Subject: [PATCH] #2018, Distance calculation support for arc features (circstring, compoundcurve, curvepolygon) git-svn-id: http://svn.osgeo.org/postgis/trunk@11219 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 2 + liblwgeom/cunit/cu_measures.c | 113 +++++- liblwgeom/cunit/cu_ptarray.c | 77 ++-- liblwgeom/liblwgeom.h.in | 8 +- liblwgeom/liblwgeom_internal.h | 6 +- liblwgeom/lwcompound.c | 89 ++++- liblwgeom/lwcurvepoly.c | 20 + liblwgeom/measures.c | 641 +++++++++++++++++++++++++-------- liblwgeom/measures.h | 16 +- liblwgeom/ptarray.c | 63 +++- regress/tickets_expected | 2 +- 11 files changed, 821 insertions(+), 216 deletions(-) diff --git a/NEWS b/NEWS index 37b7e24a..7fdf9d35 100644 --- a/NEWS +++ b/NEWS @@ -59,6 +59,8 @@ PostGIS 2.1.0 - Casts to/from PostgreSQL geotypes (point/path/polygon). - #310, ST_DumpPoints converted to a C function (Nathan Wagner) - #2011, ST_DumpValues to output raster as array (Bborie Park / UC Davis) + - #2018, ST_Distance support for CircularString, CurvePolygon, MultiCurve, + MultiSurface, CompoundCurve - #2030, n-raster (and n-band) ST_MapAlgebra (Bborie Park / UC Davis) - Added geomval array variant of ST_SetValues() to set many pixel values of a band using a set of geometries and corresponding values in one call diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c index b8bdf6e1..1689cf2c 100644 --- a/liblwgeom/cunit/cu_measures.c +++ b/liblwgeom/cunit/cu_measures.c @@ -28,82 +28,167 @@ static LWGEOM* lwgeom_from_text(const char *str) return r.geom; } -static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res) +#define DIST2DTEST(str1, str2, res) do_test_mindistance2d_tolerance(str1, str2, res, __LINE__) + +static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line) { LWGEOM *lw1; LWGEOM *lw2; double distance; + char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n"; + char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n %s\n\n"; lw1 = lwgeom_from_wkt(in1, LW_PARSER_CHECK_NONE); lw2 = lwgeom_from_wkt(in2, LW_PARSER_CHECK_NONE); + if ( ! lw1 ) + { + printf(msg2, in1); + exit(1); + } + if ( ! lw2 ) + { + printf(msg2, in2); + exit(1); + } + distance = lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0); - CU_ASSERT_EQUAL(distance, expected_res); lwgeom_free(lw1); lwgeom_free(lw2); + if ( fabs(distance - expected_res) > 0.00001 ) + { + printf(msg1, distance, expected_res, line); + CU_FAIL(); + } + else + { + CU_PASS(); + } + } + static void test_mindistance2d_tolerance(void) { /* ** Simple case. */ - do_test_mindistance2d_tolerance("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5); + DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5); /* ** Point vs Geometry Collection. */ - do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0); + DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0); /* ** Point vs Geometry Collection Collection. */ - do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0); + DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0); /* ** Point vs Geometry Collection Collection Collection. */ - do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0); + DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0); /* ** Point vs Geometry Collection Collection Collection Multipoint. */ - do_test_mindistance2d_tolerance("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0); + DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0); /* ** Geometry Collection vs Geometry Collection */ - do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0); + DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0); /* ** Geometry Collection Collection vs Geometry Collection Collection */ - do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0); + DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0); /* ** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint */ - do_test_mindistance2d_tolerance("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0); + DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0); /* ** Linestring vs its start point */ - do_test_mindistance2d_tolerance("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0); + DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0); /* ** Linestring vs its end point */ - do_test_mindistance2d_tolerance("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0); + DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0); /* ** Linestring vs its start point (tricky number, see #1459) */ - do_test_mindistance2d_tolerance("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0); + DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0); /* ** Linestring vs its end point (tricky number, see #1459) */ - do_test_mindistance2d_tolerance("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0); + DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0); + + /* + * Circular string and point + */ + DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1); + DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1); + + /* + * Circular string and Circular string + */ + DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1); + + /* + * CurvePolygon and Point + */ + static char *cs1 = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))"; + DIST2DTEST(cs1, "POINT(3 14)", 1); + DIST2DTEST(cs1, "POINT(3 8)", 0); + DIST2DTEST(cs1, "POINT(6 5)", 1); + DIST2DTEST(cs1, "POINT(6 4)", 0); + + /* + * CurvePolygon and Linestring + */ + DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484); + DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0); + DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0); + DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786); + DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913); + + /* + * CurvePolygon and Polygon + */ + DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415); + DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0); + DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0); + + /* + * CurvePolygon and CurvePolygon + */ + DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666); + DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0); + + /* + * MultiSurface and CurvePolygon + */ + static char *cs2 = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))"; + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0); + DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555); + + /* + * MultiCurve and Linestring + */ + DIST2DTEST("LINESTRING(0.5 1,0.5 3)", "MULTICURVE(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3),(0 0, 0 5))", 0.5); } diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c index d56446f3..3bd74497 100644 --- a/liblwgeom/cunit/cu_ptarray.c +++ b/liblwgeom/cunit/cu_ptarray.c @@ -409,7 +409,7 @@ static void test_ptarray_desegmentize() static void test_ptarray_contains_point() { -/* int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) */ +/* int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt, int *winding_number) */ LWLINE *lwline; POINTARRAY *pa; @@ -496,81 +496,102 @@ static void test_ptarray_contains_point() lwline_free(lwline); } - -static void test_ptarray_contains_point_arc() +static void test_ptarrayarc_contains_point() { - /* int ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) */ + /* int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt) */ LWLINE *lwline; POINTARRAY *pa; POINT2D pt; int rv; - /* Collection of semi-circles surrounding unit square */ + /*** Collection of semi-circles surrounding unit square ***/ lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)")); pa = lwline->points; /* Point in middle of square */ pt.x = 0; pt.y = 0; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point in left lobe */ pt.x = -1.1; pt.y = 0.1; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point on boundary of left lobe */ pt.x = -1; pt.y = 0; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point on boundary vertex */ pt.x = -1; pt.y = 1; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_BOUNDARY); /* Point outside */ pt.x = -1.5; pt.y = 1.5; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_OUTSIDE); - /* Two-edge ring made up of semi-circles (really, a circle) */ + /*** Two-edge ring made up of semi-circles (really, a circle) ***/ lwline_free(lwline); - lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0, 0 -1, -1 0)")); pa = lwline->points; /* Point outside */ pt.x = -1.5; pt.y = 1.5; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); + CU_ASSERT_EQUAL(rv, LW_OUTSIDE); + + /* Point more outside */ + pt.x = 2.5; + pt.y = 1.5; + rv = ptarrayarc_contains_point(pa, &pt); + CU_ASSERT_EQUAL(rv, LW_OUTSIDE); + + /* Point more outside */ + pt.x = 2.5; + pt.y = 2.5; + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_OUTSIDE); /* Point inside at middle */ pt.x = 0; pt.y = 0; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point inside offset from middle */ pt.x = 0.01; pt.y = 0.01; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point on edge vertex */ pt.x = 0; pt.y = 1; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_BOUNDARY); - /* One-edge ring, closed circle */ + /*** Two-edge ring, closed ***/ + lwline_free(lwline); + lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1 6, 6 1, 9 7, 6 10, 1 6)")); + pa = lwline->points; + + /* Point to left of ring */ + pt.x = 20; + pt.y = 4; + rv = ptarrayarc_contains_point(pa, &pt); + CU_ASSERT_EQUAL(rv, LW_OUTSIDE); + + /*** One-edge ring, closed circle ***/ lwline_free(lwline); lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, -1 0)")); pa = lwline->points; @@ -578,37 +599,37 @@ static void test_ptarray_contains_point_arc() /* Point inside */ pt.x = 0; pt.y = 0; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_INSIDE); /* Point outside */ pt.x = 0; pt.y = 2; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_OUTSIDE); /* Point on boundary */ pt.x = 0; pt.y = 1; - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); CU_ASSERT_EQUAL(rv, LW_BOUNDARY); - /* Overshort ring */ + /*** Overshort ring ***/ lwline_free(lwline); lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0)")); pa = lwline->points; cu_error_msg_reset(); - rv = ptarray_contains_point_arc(pa, &pt); + rv = ptarrayarc_contains_point(pa, &pt); //printf("%s\n", cu_error_msg); - CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called with even number of points", cu_error_msg); + CU_ASSERT_STRING_EQUAL("ptarrayarc_contains_point called with even number of points", cu_error_msg); - /* Unclosed ring */ + /*** Unclosed ring ***/ lwline_free(lwline); - lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 1 0)")); + lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 2 0)")); pa = lwline->points; cu_error_msg_reset(); - rv = ptarray_contains_point_arc(pa, &pt); - CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called on unclosed ring", cu_error_msg); + rv = ptarrayarc_contains_point(pa, &pt); + CU_ASSERT_STRING_EQUAL("ptarrayarc_contains_point called on unclosed ring", cu_error_msg); lwline_free(lwline); } @@ -627,7 +648,7 @@ CU_TestInfo ptarray_tests[] = PG_TEST(test_ptarray_desegmentize), PG_TEST(test_ptarray_insert_point), PG_TEST(test_ptarray_contains_point), - PG_TEST(test_ptarray_contains_point_arc), + PG_TEST(test_ptarrayarc_contains_point), CU_TEST_INFO_NULL }; CU_SuiteInfo ptarray_suite = {"ptarray", NULL, NULL, ptarray_tests }; diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index 84c698f1..32065ba4 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -919,9 +919,11 @@ extern int lwcurvepoly_add_ring(LWCURVEPOLY *poly, LWGEOM *ring); */ extern int lwcompound_add_lwgeom(LWCOMPOUND *comp, LWGEOM *geom); - - - +/** +* Construct an equivalent curve polygon from a polygon. Curve polygons +* can have linear rings as their rings, so this works fine (in theory?) +*/ +extern LWCURVEPOLY* lwcurvepoly_construct_from_lwpoly(LWPOLY *lwpoly); /****************************************************************** diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index 99f48594..5cff2944 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -381,7 +381,11 @@ double lw_seg_length(const POINT2D *A1, const POINT2D *A2); double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3); int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring); int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt); -int ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt); +int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt); +int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number); +int ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number); +int lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt); +int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt); /** * Split a line by a point and push components to the provided multiline. diff --git a/liblwgeom/lwcompound.c b/liblwgeom/lwcompound.c index d4ec2414..fc01b13c 100644 --- a/liblwgeom/lwcompound.c +++ b/liblwgeom/lwcompound.c @@ -24,19 +24,31 @@ lwcompound_is_closed(const LWCOMPOUND *compound) size_t size; int npoints=0; - if (!FLAGS_GET_Z(compound->flags)) + if ( lwgeom_has_z((LWGEOM*)compound) ) + { + size = sizeof(POINT3D); + } + else + { size = sizeof(POINT2D); - else size = sizeof(POINT3D); + } - if (compound->geoms[compound->ngeoms - 1]->type == CIRCSTRINGTYPE) + if ( compound->geoms[compound->ngeoms - 1]->type == CIRCSTRINGTYPE ) + { npoints = ((LWCIRCSTRING *)compound->geoms[compound->ngeoms - 1])->points->npoints; + } else if (compound->geoms[compound->ngeoms - 1]->type == LINETYPE) + { npoints = ((LWLINE *)compound->geoms[compound->ngeoms - 1])->points->npoints; + } if ( memcmp(getPoint_internal( (POINTARRAY *)compound->geoms[0]->data, 0), getPoint_internal( (POINTARRAY *)compound->geoms[compound->ngeoms - 1]->data, npoints - 1), - size) ) return LW_FALSE; + size) ) + { + return LW_FALSE; + } return LW_TRUE; } @@ -106,3 +118,72 @@ lwcompound_construct_empty(int srid, char hasz, char hasm) return ret; } +int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt) +{ + switch( geom->type ) + { + case LINETYPE: + return ptarray_contains_point(((LWLINE*)geom)->points, pt); + case CIRCSTRINGTYPE: + return ptarrayarc_contains_point(((LWCIRCSTRING*)geom)->points, pt); + case COMPOUNDTYPE: + return lwcompound_contains_point((LWCOMPOUND*)geom, pt); + } + lwerror("lwgeom_contains_point failed"); + return LW_FAILURE; +} + +int +lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt) +{ + int i; + LWLINE *lwline; + LWCIRCSTRING *lwcirc; + int wn = 0; + int winding_number = 0; + int result; + + for ( i = 0; i < comp->ngeoms; i++ ) + { + LWGEOM *lwgeom = comp->geoms[i]; + if ( lwgeom->type == LINETYPE ) + { + lwline = lwgeom_as_lwline(lwgeom); + if ( comp->ngeoms == 1 ) + { + return ptarray_contains_point(lwline->points, pt); + } + else + { + /* Don't check closure while doing p-i-p test */ + result = ptarray_contains_point_partial(lwline->points, pt, LW_FALSE, &winding_number); + } + } + else if ( lwgeom->type == CIRCSTRINGTYPE ) + { + lwcirc = lwgeom_as_lwcircstring(lwgeom); + if ( comp->ngeoms == 1 ) + { + return ptarrayarc_contains_point(lwcirc->points, pt); + } + else + { + /* Don't check closure while doing p-i-p test */ + result = ptarrayarc_contains_point_partial(lwcirc->points, pt, LW_FALSE, &winding_number); + } + } + + /* Propogate boundary condition */ + if ( result == LW_BOUNDARY ) + return LW_BOUNDARY; + + wn += winding_number; + } + + /* Outside */ + if (wn == 0) + return LW_OUTSIDE; + + /* Inside */ + return LW_INSIDE; +} diff --git a/liblwgeom/lwcurvepoly.c b/liblwgeom/lwcurvepoly.c index 0235c754..0168e240 100644 --- a/liblwgeom/lwcurvepoly.c +++ b/liblwgeom/lwcurvepoly.c @@ -36,6 +36,26 @@ lwcurvepoly_construct_empty(int srid, char hasz, char hasm) return ret; } +LWCURVEPOLY * +lwcurvepoly_construct_from_lwpoly(LWPOLY *lwpoly) +{ + LWCURVEPOLY *ret; + int i; + ret = lwalloc(sizeof(LWCURVEPOLY)); + ret->type = CURVEPOLYTYPE; + ret->flags = lwpoly->flags; + ret->srid = lwpoly->srid; + ret->nrings = lwpoly->nrings; + ret->maxrings = lwpoly->nrings; /* Allocate room for sub-members, just in case. */ + ret->rings = lwalloc(ret->maxrings * sizeof(LWGEOM*)); + ret->bbox = lwpoly->bbox; + for ( i = 0; i < ret->nrings; i++ ) + { + ret->rings[i] = lwline_as_lwgeom(lwline_construct(ret->srid, NULL, ptarray_clone_deep(lwpoly->rings[i]))); + } + return ret; +} + int lwcurvepoly_add_ring(LWCURVEPOLY *poly, LWGEOM *ring) { int i; diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index 30e53e42..00021fa8 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -209,6 +209,27 @@ lw_dist2d_comp(LWGEOM *lw1, LWGEOM *lw2, DISTPTS *dl) return lw_dist2d_recursive(lw1, lw2, dl); } +static int +lw_dist2d_is_collection(const LWGEOM *g) +{ + + switch (g->type) + { + case MULTIPOINTTYPE: + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + case MULTICURVETYPE: + case MULTISURFACETYPE: + case COMPOUNDTYPE: + return LW_TRUE; + break; + + default: + return LW_FALSE; + } +} + /** This is a recursive function delivering every possible combinatin of subgeometries */ @@ -224,13 +245,13 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) LWDEBUGF(2, "lw_dist2d_comp is called with type1=%d, type2=%d", lwg1->type, lwg2->type); - if (lwgeom_is_collection(lwg1)) + if (lw_dist2d_is_collection(lwg1)) { LWDEBUG(3, "First geometry is collection"); c1 = lwgeom_as_lwcollection(lwg1); n1 = c1->ngeoms; } - if (lwgeom_is_collection(lwg2)) + if (lw_dist2d_is_collection(lwg2)) { LWDEBUG(3, "Second geometry is collection"); c2 = lwgeom_as_lwcollection(lwg2); @@ -240,7 +261,7 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) for ( i = 0; i < n1; i++ ) { - if (lwgeom_is_collection(lwg1)) + if (lw_dist2d_is_collection(lwg1)) { g1 = c1->geoms[i]; } @@ -251,7 +272,7 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) if (lwgeom_is_empty(g1)) return LW_TRUE; - if (lwgeom_is_collection(g1)) + if (lw_dist2d_is_collection(g1)) { LWDEBUG(3, "Found collection inside first geometry collection, recursing"); if (!lw_dist2d_recursive(g1, lwg2, dl)) return LW_FALSE; @@ -259,7 +280,7 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) } for ( j = 0; j < n2; j++ ) { - if (lwgeom_is_collection(lwg2)) + if (lw_dist2d_is_collection(lwg2)) { g2 = c2->geoms[j]; } @@ -267,7 +288,7 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) { g2 = (LWGEOM*)lwg2; } - if (lwgeom_is_collection(g2)) + if (lw_dist2d_is_collection(g2)) { LWDEBUG(3, "Found collection inside second geometry collection, recursing"); if (!lw_dist2d_recursive(g1, g2, dl)) return LW_FALSE; @@ -286,14 +307,17 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/ if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE; - if ((dl->mode==DIST_MAX)||(g1->type==POINTTYPE) ||(g2->type==POINTTYPE)||lw_dist2d_check_overlap(g1, g2)) + if ( (dl->mode != DIST_MAX) && + (! lw_dist2d_check_overlap(g1, g2)) && + (g1->type == LINETYPE || g1->type == POLYGONTYPE) && + (g2->type == LINETYPE || g2->type == POLYGONTYPE) ) { - if (!lw_dist2d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE; - if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/ + if (!lw_dist2d_distribute_fast(g1, g2, dl)) return LW_FALSE; } else { - if (!lw_dist2d_distribute_fast(g1, g2, dl)) return LW_FALSE; + if (!lw_dist2d_distribute_bruteforce2(g1, g2, dl)) return LW_FALSE; + if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/ } } } @@ -301,9 +325,130 @@ int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl) } +int +lw_dist2d_distribute_bruteforce2(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl) +{ + + int t1 = lwg1->type; + int t2 = lwg2->type; + + switch ( t1 ) + { + case POINTTYPE: + { + dl->twisted = 1; + switch ( t2 ) + { + case POINTTYPE: + return lw_dist2d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl); + case LINETYPE: + return lw_dist2d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl); + case POLYGONTYPE: + return lw_dist2d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2, dl); + case CIRCSTRINGTYPE: + return lw_dist2d_point_circstring((LWPOINT *)lwg1, (LWCIRCSTRING *)lwg2, dl); + case CURVEPOLYTYPE: + return lw_dist2d_point_curvepoly((LWPOINT *)lwg1, (LWCURVEPOLY *)lwg2, dl); + default: + lwerror("Unsupported geometry type: %s", lwtype_name(t2)); + } + } + case LINETYPE: + { + dl->twisted = 1; + switch ( t2 ) + { + case POINTTYPE: + dl->twisted=(-1); + return lw_dist2d_point_line((LWPOINT *)lwg2, (LWLINE *)lwg1, dl); + case LINETYPE: + return lw_dist2d_line_line((LWLINE *)lwg1, (LWLINE *)lwg2, dl); + case POLYGONTYPE: + return lw_dist2d_line_poly((LWLINE *)lwg1, (LWPOLY *)lwg2, dl); + case CIRCSTRINGTYPE: + return lw_dist2d_line_circstring((LWLINE *)lwg1, (LWCIRCSTRING *)lwg2, dl); + case CURVEPOLYTYPE: + return lw_dist2d_line_curvepoly((LWLINE *)lwg1, (LWCURVEPOLY *)lwg2, dl); + default: + lwerror("Unsupported geometry type: %s", lwtype_name(t2)); + } + } + case CIRCSTRINGTYPE: + { + dl->twisted = 1; + switch ( t2 ) + { + case POINTTYPE: + dl->twisted = -1; + return lw_dist2d_point_circstring((LWPOINT *)lwg2, (LWCIRCSTRING *)lwg1, dl); + case LINETYPE: + dl->twisted = -1; + return lw_dist2d_line_circstring((LWLINE *)lwg2, (LWCIRCSTRING *)lwg1, dl); + case POLYGONTYPE: + return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg1, (LWPOLY *)lwg2, dl); + case CIRCSTRINGTYPE: + return lw_dist2d_circstring_circstring((LWCIRCSTRING *)lwg1, (LWCIRCSTRING *)lwg2, dl); + case CURVEPOLYTYPE: + return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg1, (LWCURVEPOLY *)lwg2, dl); + default: + lwerror("Unsupported geometry type: %s", lwtype_name(t2)); + } + } + case POLYGONTYPE: + { + dl->twisted = -1; + switch ( t2 ) + { + case POINTTYPE: + return lw_dist2d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1, dl); + case LINETYPE: + return lw_dist2d_line_poly((LWLINE *)lwg2, (LWPOLY *)lwg1, dl); + case CIRCSTRINGTYPE: + return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg2, (LWPOLY *)lwg1, dl); + case POLYGONTYPE: + dl->twisted = 1; + return lw_dist2d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2, dl); + case CURVEPOLYTYPE: + dl->twisted = 1; + return lw_dist2d_poly_curvepoly((LWPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl); + default: + lwerror("Unsupported geometry type: %s", lwtype_name(t2)); + } + } + case CURVEPOLYTYPE: + { + dl->twisted = (-1); + switch ( t2 ) + { + case POINTTYPE: + return lw_dist2d_point_curvepoly((LWPOINT *)lwg2, (LWCURVEPOLY *)lwg1, dl); + case LINETYPE: + return lw_dist2d_line_curvepoly((LWLINE *)lwg2, (LWCURVEPOLY *)lwg1, dl); + case POLYGONTYPE: + return lw_dist2d_poly_curvepoly((LWPOLY *)lwg2, (LWCURVEPOLY *)lwg1, dl); + case CIRCSTRINGTYPE: + return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg2, (LWCURVEPOLY *)lwg1, dl); + case CURVEPOLYTYPE: + dl->twisted = 1; + return lw_dist2d_curvepoly_curvepoly((LWCURVEPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl); + default: + lwerror("Unsupported geometry type: %s", lwtype_name(t2)); + } + } + default: + { + lwerror("Unsupported geometry type: %s", lwtype_name(t1)); + } + } + + /*You shouldn't being able to get here*/ + lwerror("unspecified error in function lw_dist2d_distribute_bruteforce"); + return LW_FALSE; +} -/** + +/** This function distributes the "old-type" brut-force tasks depending on type */ int @@ -494,12 +639,20 @@ point to line calculation int lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl) { - POINT2D p; - POINTARRAY *pa = line->points; + const POINT2D *p; LWDEBUG(2, "lw_dist2d_point_line is called"); - getPoint2d_p(point->point, 0, &p); - return lw_dist2d_pt_ptarray(&p, pa, dl); + p = getPoint2d_cp(point->point, 0); + return lw_dist2d_pt_ptarray(p, line->points, dl); } + +int +lw_dist2d_point_circstring(LWPOINT *point, LWCIRCSTRING *circ, DISTPTS *dl) +{ + const POINT2D *p; + p = getPoint2d_cp(point->point, 0); + return lw_dist2d_pt_ptarrayarc(p, circ->points, dl); +} + /** * 1. see if pt in outer boundary. if no, then treat the outer ring like a line * 2. if in the boundary, test to see if its in a hole. @@ -508,23 +661,23 @@ lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl) int lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl) { - POINT2D p; + const POINT2D *p; int i; LWDEBUG(2, "lw_dist2d_point_poly called"); - getPoint2d_p(point->point, 0, &p); + p = getPoint2d_cp(point->point, 0); if (dl->mode == DIST_MAX) { LWDEBUG(3, "looking for maxdistance"); - return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl); + return lw_dist2d_pt_ptarray(p, poly->rings[0], dl); } /* Return distance to outer ring if not inside it */ - if ( ptarray_contains_point(poly->rings[0], &p) == LW_OUTSIDE ) + if ( ptarray_contains_point(poly->rings[0], p) == LW_OUTSIDE ) { LWDEBUG(3, "first point not inside outer-ring"); - return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl); + return lw_dist2d_pt_ptarray(p, poly->rings[0], dl); } /* @@ -533,27 +686,70 @@ lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl) * see if its inside. If not, distance==0. * Otherwise, distance = pt to ring distance */ - for (i=1; inrings; i++) + for ( i = 1; i < poly->nrings; i++) { /* Inside a hole. Distance = pt -> ring */ - if ( ptarray_contains_point(poly->rings[i], &p) != LW_OUTSIDE ) + if ( ptarray_contains_point(poly->rings[i], p) != LW_OUTSIDE ) { LWDEBUG(3, " inside an hole"); - return lw_dist2d_pt_ptarray(&p, poly->rings[i], dl); + return lw_dist2d_pt_ptarray(p, poly->rings[i], dl); + } + } + + LWDEBUG(3, " inside the polygon"); + if (dl->mode == DIST_MIN) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = p->x; + dl->p1.y = dl->p2.y = p->y; + } + return LW_TRUE; /* Is inside the polygon */ +} + +int +lw_dist2d_point_curvepoly(LWPOINT *point, LWCURVEPOLY *poly, DISTPTS *dl) +{ + const POINT2D *p; + int i; + + p = getPoint2d_cp(point->point, 0); + + if (dl->mode == DIST_MAX) + lwerror("lw_dist2d_point_curvepoly cannot calculate max distance"); + + /* Return distance to outer ring if not inside it */ + if ( lwgeom_contains_point(poly->rings[0], p) == LW_OUTSIDE ) + { + return lw_dist2d_recursive((LWGEOM*)point, poly->rings[0], dl); + } + + /* + * Inside the outer ring. + * Scan though each of the inner rings looking to + * see if its inside. If not, distance==0. + * Otherwise, distance = pt to ring distance + */ + for ( i = 1; i < poly->nrings; i++) + { + /* Inside a hole. Distance = pt -> ring */ + if ( lwgeom_contains_point(poly->rings[i], p) != LW_OUTSIDE ) + { + LWDEBUG(3, " inside a hole"); + return lw_dist2d_recursive((LWGEOM*)point, poly->rings[i], dl); } } LWDEBUG(3, " inside the polygon"); if (dl->mode == DIST_MIN) { - dl->distance=0.0; - dl->p1.x=p.x; - dl->p1.y=p.y; - dl->p2.x=p.x; - dl->p2.y=p.y; + dl->distance = 0.0; + dl->p1.x = dl->p2.x = p->x; + dl->p1.y = dl->p2.y = p->y; } + return LW_TRUE; /* Is inside the polygon */ } + /** line to line calculation @@ -566,19 +762,127 @@ lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl) LWDEBUG(2, "lw_dist2d_line_line is called"); return lw_dist2d_ptarray_ptarray(pa1, pa2, dl); } -/** -line to polygon calculation -*/ int -lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl) +lw_dist2d_line_circstring(LWLINE *line1, LWCIRCSTRING *line2, DISTPTS *dl) { - LWDEBUG(2, "lw_dist2d_line_poly is called"); - return lw_dist2d_ptarray_poly(line->points, poly, dl); + return lw_dist2d_ptarray_ptarrayarc(line1->points, line2->points, dl); } /** + * line to polygon calculation + * Brute force. + * Test line-ring distance against each ring. + * If there's an intersection (distance==0) then return 0 (crosses boundary). + * Otherwise, test to see if any point is inside outer rings of polygon, + * but not in inner rings. + * If so, return 0 (line inside polygon), + * otherwise return min distance to a ring (could be outside + * polygon or inside a hole) + */ +int +lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl) +{ + const POINT2D *pt; + int i; + + LWDEBUGF(2, "lw_dist2d_line_poly called (%d rings)", poly->nrings); + + pt = getPoint2d_cp(line->points, 0); + if ( ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE ) + { + return lw_dist2d_ptarray_ptarray(line->points, poly->rings[0], dl); + } + + for (i=1; inrings; i++) + { + if (!lw_dist2d_ptarray_ptarray(line->points, poly->rings[i], dl)) return LW_FALSE; + + LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", + i, dl->distance, dl->tolerance); + /* just a check if the answer is already given */ + if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; + } + + /* + * No intersection, have to check if a point is + * inside polygon + */ + pt = getPoint2d_cp(line->points, 0); + + /* + * Outside outer ring, so min distance to a ring + * is the actual min distance + + if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) + { + return ; + } */ + + /* + * Its in the outer ring. + * Have to check if its inside a hole + */ + for (i=1; inrings; i++) + { + if ( ptarray_contains_point(poly->rings[i], pt) != LW_OUTSIDE ) + { + /* + * Its inside a hole, then the actual + * distance is the min ring distance + */ + return LW_TRUE; + } + } + if (dl->mode == DIST_MIN) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; + } + return LW_TRUE; /* Not in hole, so inside polygon */ +} + +int +lw_dist2d_line_curvepoly(LWLINE *line, LWCURVEPOLY *poly, DISTPTS *dl) +{ + const POINT2D *pt = getPoint2d_cp(line->points, 0); + int i; + + if ( lwgeom_contains_point(poly->rings[0], pt) == LW_OUTSIDE ) + { + return lw_dist2d_recursive((LWGEOM*)line, poly->rings[0], dl); + } + + for ( i = 1; i < poly->nrings; i++ ) + { + if ( ! lw_dist2d_recursive((LWGEOM*)line, poly->rings[i], dl) ) + return LW_FALSE; + + if ( dl->distance<=dl->tolerance && dl->mode == DIST_MIN ) + return LW_TRUE; + } + + for ( i=1; i < poly->nrings; i++ ) + { + if ( lwgeom_contains_point(poly->rings[i],pt) != LW_OUTSIDE ) + { + /* Its inside a hole, then the actual */ + return LW_TRUE; + } + } + + if (dl->mode == DIST_MIN) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; + } + return LW_TRUE; /* Not in hole, so inside polygon */ +} + +/** Function handling polygon to polygon calculation 1 if we are looking for maxdistance, just check the outer rings. 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings @@ -590,7 +894,7 @@ int lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) { - POINT2D pt; + const POINT2D *pt; int i; LWDEBUG(2, "lw_dist2d_poly_poly called"); @@ -598,65 +902,61 @@ lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) /*1 if we are looking for maxdistance, just check the outer rings.*/ if (dl->mode == DIST_MAX) { - return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); + return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl); } /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings here it would be possible to handle the information about wich one is inside wich one and only search for the smaller ones in the bigger ones holes.*/ - getPoint2d_p(poly1->rings[0], 0, &pt); - if ( ptarray_contains_point(poly2->rings[0], &pt) == LW_OUTSIDE ) + pt = getPoint2d_cp(poly1->rings[0], 0); + if ( ptarray_contains_point(poly2->rings[0], pt) == LW_OUTSIDE ) { - getPoint2d_p(poly2->rings[0], 0, &pt); - if ( ptarray_contains_point(poly1->rings[0], &pt) == LW_OUTSIDE ) + pt = getPoint2d_cp(poly2->rings[0], 0); + if ( ptarray_contains_point(poly1->rings[0], pt) == LW_OUTSIDE ) { - return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); + return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl); } } /*3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1*/ - getPoint2d_p(poly2->rings[0], 0, &pt); + pt = getPoint2d_cp(poly2->rings[0], 0); for (i=1; inrings; i++) { /* Inside a hole */ - if ( ptarray_contains_point(poly1->rings[i], &pt) != LW_OUTSIDE ) + if ( ptarray_contains_point(poly1->rings[i], pt) != LW_OUTSIDE ) { - return lw_dist2d_ptarray_ptarray(poly1->rings[i], poly2->rings[0],dl); + return lw_dist2d_ptarray_ptarray(poly1->rings[i], poly2->rings[0], dl); } } /*4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2*/ - getPoint2d_p(poly1->rings[0], 0, &pt); + pt = getPoint2d_cp(poly1->rings[0], 0); for (i=1; inrings; i++) { /* Inside a hole */ - if ( ptarray_contains_point(poly2->rings[i], &pt) != LW_OUTSIDE ) + if ( ptarray_contains_point(poly2->rings[i], pt) != LW_OUTSIDE ) { - return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[i],dl); + return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[i], dl); } } /*5 If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside.*/ - getPoint2d_p(poly1->rings[0], 0, &pt); - if ( ptarray_contains_point(poly2->rings[0], &pt) != LW_OUTSIDE ) - { - dl->distance=0.0; - dl->p1.x=pt.x; - dl->p1.y=pt.y; - dl->p2.x=pt.x; - dl->p2.y=pt.y; + pt = getPoint2d_cp(poly1->rings[0], 0); + if ( ptarray_contains_point(poly2->rings[0], pt) != LW_OUTSIDE ) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; return LW_TRUE; } - getPoint2d_p(poly2->rings[0], 0, &pt); - if ( ptarray_contains_point(poly1->rings[0], &pt) != LW_OUTSIDE ) + pt = getPoint2d_cp(poly2->rings[0], 0); + if ( ptarray_contains_point(poly1->rings[0], pt) != LW_OUTSIDE ) { - dl->distance=0.0; - dl->p1.x=pt.x; - dl->p1.y=pt.y; - dl->p2.x=pt.x; - dl->p2.y=pt.y; + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; return LW_TRUE; } @@ -665,28 +965,155 @@ lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) return LW_FALSE; } +int +lw_dist2d_poly_curvepoly(LWPOLY *poly1, LWCURVEPOLY *curvepoly2, DISTPTS *dl) +{ + LWCURVEPOLY *curvepoly1 = lwcurvepoly_construct_from_lwpoly(poly1); + int rv = lw_dist2d_curvepoly_curvepoly(curvepoly1, curvepoly2, dl); + lwfree(curvepoly1); + return rv; +} + +int +lw_dist2d_circstring_poly(LWCIRCSTRING *circ, LWPOLY *poly, DISTPTS *dl) +{ + LWCURVEPOLY *curvepoly = lwcurvepoly_construct_from_lwpoly(poly); + int rv = lw_dist2d_line_curvepoly((LWLINE*)circ, curvepoly, dl); + lwfree(curvepoly); + return rv; +} + + +int +lw_dist2d_circstring_curvepoly(LWCIRCSTRING *circ, LWCURVEPOLY *poly, DISTPTS *dl) +{ + return lw_dist2d_line_curvepoly((LWLINE*)circ, poly, dl); +} + +int +lw_dist2d_circstring_circstring(LWCIRCSTRING *line1, LWCIRCSTRING *line2, DISTPTS *dl) +{ + return lw_dist2d_ptarrayarc_ptarrayarc(line1->points, line2->points, dl); +} + +static const POINT2D * +lw_curvering_getfirstpoint2d_cp(LWGEOM *geom) +{ + switch( geom->type ) + { + case LINETYPE: + return getPoint2d_cp(((LWLINE*)geom)->points, 0); + case CIRCSTRINGTYPE: + return getPoint2d_cp(((LWCIRCSTRING*)geom)->points, 0); + case COMPOUNDTYPE: + { + LWCOMPOUND *comp = (LWCOMPOUND*)geom; + LWLINE *line = (LWLINE*)(comp->geoms[0]); + return getPoint2d_cp(line->points, 0); + } + default: + lwerror("lw_curvering_getfirstpoint2d_cp: unknown type"); + } + return NULL; +} + +int +lw_dist2d_curvepoly_curvepoly(LWCURVEPOLY *poly1, LWCURVEPOLY *poly2, DISTPTS *dl) +{ + const POINT2D *pt; + int i; + + LWDEBUG(2, "lw_dist2d_curvepoly_curvepoly called"); + + /*1 if we are looking for maxdistance, just check the outer rings.*/ + if (dl->mode == DIST_MAX) + { + return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl); + } + + + /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings + here it would be possible to handle the information about wich one is inside wich one and only search for the smaller ones in the bigger ones holes.*/ + pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]); + if ( lwgeom_contains_point(poly2->rings[0], pt) == LW_OUTSIDE ) + { + pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]); + if ( lwgeom_contains_point(poly1->rings[0], pt) == LW_OUTSIDE ) + { + return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl); + } + } + + /*3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1*/ + pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]); + for (i = 1; i < poly1->nrings; i++) + { + /* Inside a hole */ + if ( lwgeom_contains_point(poly1->rings[i], pt) != LW_OUTSIDE ) + { + return lw_dist2d_recursive(poly1->rings[i], poly2->rings[0], dl); + } + } + + /*4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2*/ + pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]); + for (i = 1; i < poly2->nrings; i++) + { + /* Inside a hole */ + if ( lwgeom_contains_point(poly2->rings[i], pt) != LW_OUTSIDE ) + { + return lw_dist2d_recursive(poly1->rings[0], poly2->rings[i], dl); + } + } + + + /*5 If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside.*/ + pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]); + if ( lwgeom_contains_point(poly2->rings[0], pt) != LW_OUTSIDE ) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; + return LW_TRUE; + } + + pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]); + if ( lwgeom_contains_point(poly1->rings[0], pt) != LW_OUTSIDE ) + { + dl->distance = 0.0; + dl->p1.x = dl->p2.x = pt->x; + dl->p1.y = dl->p2.y = pt->y; + return LW_TRUE; + } + + lwerror("Unspecified error in function lw_dist2d_curvepoly_curvepoly"); + return LW_FALSE; +} + + + /** * search all the segments of pointarray to see which one is closest to p1 * Returns minimum distance between point and pointarray */ int -lw_dist2d_pt_ptarray(POINT2D *p, POINTARRAY *pa,DISTPTS *dl) +lw_dist2d_pt_ptarray(const POINT2D *p, POINTARRAY *pa,DISTPTS *dl) { int t; - POINT2D start, end; + const POINT2D *start, *end; int twist = dl->twisted; LWDEBUG(2, "lw_dist2d_pt_ptarray is called"); - getPoint2d_p(pa, 0, &start); + start = getPoint2d_cp(pa, 0); - if ( !lw_dist2d_pt_pt(p, &start, dl) ) return LW_FALSE; + if ( !lw_dist2d_pt_pt(p, start, dl) ) return LW_FALSE; for (t=1; tnpoints; t++) { dl->twisted=twist; - getPoint2d_p(pa, t, &end); - if (!lw_dist2d_pt_seg(p, &start, &end,dl)) return LW_FALSE; + end = getPoint2d_cp(pa, t); + if (!lw_dist2d_pt_seg(p, start, end, dl)) return LW_FALSE; if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/ start = end; @@ -746,82 +1173,6 @@ lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) } -/** - * Brute force. - * Test line-ring distance against each ring. - * If there's an intersection (distance==0) then return 0 (crosses boundary). - * Otherwise, test to see if any point is inside outer rings of polygon, - * but not in inner rings. - * If so, return 0 (line inside polygon), - * otherwise return min distance to a ring (could be outside - * polygon or inside a hole) - */ - -int -lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl) -{ - POINT2D pt; - int i; - - LWDEBUGF(2, "lw_dist2d_ptarray_poly called (%d rings)", poly->nrings); - - getPoint2d_p(pa, 0, &pt); - if ( ptarray_contains_point(poly->rings[0], &pt) == LW_OUTSIDE ) - { - return lw_dist2d_ptarray_ptarray(pa,poly->rings[0],dl); - } - - for (i=1; inrings; i++) - { - if (!lw_dist2d_ptarray_ptarray(pa, poly->rings[i], dl)) return LW_FALSE; - - LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", - i, dl->distance, dl->tolerance); - /* just a check if the answer is already given */ - if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; - } - - /* - * No intersection, have to check if a point is - * inside polygon - */ - getPoint2d_p(pa, 0, &pt); - - /* - * Outside outer ring, so min distance to a ring - * is the actual min distance - - if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) - { - return ; - } */ - - /* - * Its in the outer ring. - * Have to check if its inside a hole - */ - for (i=1; inrings; i++) - { - if ( ptarray_contains_point(poly->rings[i], &pt) != LW_OUTSIDE ) - { - /* - * Its inside a hole, then the actual - * distance is the min ring distance - */ - return LW_TRUE; - } - } - if (dl->mode == DIST_MIN) - { - dl->distance=0.0; - dl->p1.x=pt.x; - dl->p1.y=pt.y; - dl->p2.x=pt.x; - dl->p2.y=pt.y; - } - return LW_TRUE; /* Not in hole, so inside polygon */ -} - /** @@ -955,7 +1306,7 @@ lw_dist2d_ptarrayarc_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DIST else { A1 = getPoint2d_cp(pa, 0); - for ( t=1; t < pa->npoints; t++ ) /* For each segment in pa */ + for ( t=1; t < pa->npoints; t += 2 ) /* For each segment in pa */ { A2 = getPoint2d_cp(pa, t); A3 = getPoint2d_cp(pa, t+1); diff --git a/liblwgeom/measures.h b/liblwgeom/measures.h index 0a3958dd..a61f11e9 100644 --- a/liblwgeom/measures.h +++ b/liblwgeom/measures.h @@ -40,6 +40,7 @@ typedef struct */ int lw_dist2d_comp(LWGEOM *lw1, LWGEOM *lw2, DISTPTS *dl); int lw_dist2d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); +int lw_dist2d_distribute_bruteforce2(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl); int lw_dist2d_check_overlap(LWGEOM *lwg1,LWGEOM *lwg2); int lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); @@ -47,7 +48,7 @@ int lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl); /* * Brute force functions */ -int lw_dist2d_pt_ptarray(POINT2D *p, POINTARRAY *pa, DISTPTS *dl); +int lw_dist2d_pt_ptarray(const POINT2D *p, POINTARRAY *pa, DISTPTS *dl); int lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl); int lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl); int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl); @@ -55,10 +56,19 @@ int lw_dist2d_ptarrayarc_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, int lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl); int lw_dist2d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl); int lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl); -int lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); +int lw_dist2d_point_circstring(LWPOINT *point, LWCIRCSTRING *circ, DISTPTS *dl); int lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl); -int lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); +int lw_dist2d_point_curvepoly(LWPOINT *point, LWCURVEPOLY *poly, DISTPTS *dl); +int lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); +int lw_dist2d_line_circstring(LWLINE *line1, LWCIRCSTRING *line2, DISTPTS *dl); int lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl); +int lw_dist2d_line_curvepoly(LWLINE *line, LWCURVEPOLY *poly, DISTPTS *dl); +int lw_dist2d_circstring_circstring(LWCIRCSTRING *line1, LWCIRCSTRING *line2, DISTPTS *dl); +int lw_dist2d_circstring_poly(LWCIRCSTRING *circ, LWPOLY *poly, DISTPTS *dl); +int lw_dist2d_circstring_curvepoly(LWCIRCSTRING *circ, LWCURVEPOLY *poly, DISTPTS *dl); +int lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); +int lw_dist2d_poly_curvepoly(LWPOLY *poly1, LWCURVEPOLY *curvepoly2, DISTPTS *dl); +int lw_dist2d_curvepoly_curvepoly(LWCURVEPOLY *poly1, LWCURVEPOLY *poly2, DISTPTS *dl); /* * New faster distance calculations diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c index 4877edc3..6baad62d 100644 --- a/liblwgeom/ptarray.c +++ b/liblwgeom/ptarray.c @@ -695,6 +695,12 @@ ptarray_is_closed_z(const POINTARRAY *in) */ int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) +{ + return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL); +} + +int +ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number) { int wn = 0; int i; @@ -705,7 +711,7 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) seg1 = getPoint2d_cp(pa, 0); seg2 = getPoint2d_cp(pa, pa->npoints-1); - if ( ! p2d_same(seg1, seg2) ) + if ( check_closed && ! p2d_same(seg1, seg2) ) lwerror("ptarray_contains_point called on unclosed ring"); for ( i=1; i < pa->npoints; i++ ) @@ -763,6 +769,10 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) seg1 = seg2; } + /* Sent out the winding number for calls that are building on this as a primitive */ + if ( winding_number ) + *winding_number = wn; + /* Outside */ if (wn == 0) { @@ -781,8 +791,15 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) * Return 1 if the point is inside the POINTARRAY, -1 if it is outside, * and 0 if it is on the boundary. */ + +int +ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt) +{ + return ptarrayarc_contains_point_partial(pa, pt, LW_TRUE /* Check closed*/, NULL); +} + int -ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) +ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number) { int wn = 0; int i, side; @@ -794,27 +811,27 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) /* Check for not an arc ring (always have odd # of points) */ if ( (pa->npoints % 2) == 0 ) { - lwerror("ptarray_contains_point_arc called with even number of points"); + lwerror("ptarrayarc_contains_point called with even number of points"); 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"); + lwerror("ptarrayarc_contains_point called too-short pointarray"); return LW_OUTSIDE; } /* Check for unclosed case */ seg1 = getPoint2d_cp(pa, 0); seg3 = getPoint2d_cp(pa, pa->npoints-1); - if ( ! p2d_same(seg1, seg3) ) + if ( check_closed && ! p2d_same(seg1, seg3) ) { - lwerror("ptarray_contains_point_arc called on unclosed ring"); + lwerror("ptarrayarc_contains_point called on unclosed ring"); return LW_OUTSIDE; } /* OK, it's closed. Is it just one circle? */ - else if ( pa->npoints == 3 ) + else if ( p2d_same(seg1, seg3) && pa->npoints == 3 ) { double radius, d; POINT2D c; @@ -822,21 +839,21 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) /* Wait, it's just a point, so it can't contain anything */ if ( lw_arc_is_pt(seg1, seg2, seg3) ) - return -1; + return LW_OUTSIDE; /* See if the point is within the circle radius */ radius = lw_arc_center(seg1, seg2, seg3, &c); d = distance2d_pt_pt(pt, &c); if ( FP_EQUALS(d, radius) ) - return 0; /* Boundary of circle */ + return LW_BOUNDARY; /* Boundary of circle */ else if ( d < radius ) - return 1; /* Inside circle */ + return LW_INSIDE; /* Inside circle */ else return LW_OUTSIDE; /* Outside circle */ } - else if ( p2d_same(seg1, pt) ) + else if ( p2d_same(seg1, pt) || p2d_same(seg3, pt) ) { - return 0; /* Boundary case */ + return LW_BOUNDARY; /* Boundary case */ } /* Start on the ring */ @@ -848,7 +865,7 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) /* Catch an easy boundary case */ if( p2d_same(seg3, pt) ) - return 0; + return LW_BOUNDARY; /* Skip arcs that have no size */ if ( lw_arc_is_pt(seg1, seg2, seg3) ) @@ -864,13 +881,21 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) seg1 = seg3; continue; } + + /* Outside of horizontal range, and not between end points we also skip */ + if ( (pt->x > gbox.xmax || pt->x < gbox.xmin) && + (pt->y > FP_MAX(seg1->y, seg3->y) || pt->y < FP_MIN(seg1->y, seg3->y)) ) + { + seg1 = seg3; + continue; + } side = lw_arc_side(seg1, seg2, seg3, pt); /* On the boundary */ if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) ) { - return 0; + return LW_BOUNDARY; } /* Going "up"! Point to left of arc. */ @@ -894,7 +919,7 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) /* On the boundary! */ if ( d == radius ) - return 0; + return LW_BOUNDARY; /* Within the arc! */ if ( d < radius ) @@ -911,14 +936,18 @@ ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) seg1 = seg3; } + /* Sent out the winding number for calls that are building on this as a primitive */ + if ( winding_number ) + *winding_number = wn; + /* Outside */ if (wn == 0) { - return -1; + return LW_OUTSIDE; } /* Inside */ - return 1; + return LW_INSIDE; } /** diff --git a/regress/tickets_expected b/regress/tickets_expected index 50a3acb7..91bda278 100644 --- a/regress/tickets_expected +++ b/regress/tickets_expected @@ -15,7 +15,7 @@ ERROR: lwgeom_longitude_shift: unsupported geom type: CircularString #73|GEOMETRYCOLLECTION(CIRCULARSTRING(1 1,2 3,4 5,6 7,5 6)) #80|MULTILINESTRING((0 0,1 1)) #83|MULTICURVE(CIRCULARSTRING(220268 150415,220227 150505,220227 150406)) -ERROR: Unsupported geometry type: CircularString +#85|0 #112|GEOMETRYCOLLECTION(POINT(-10 50)) NOTICE: ST_Locate_Between_Measures and ST_Locate_Along_Measure are deprecated. Use ST_LocateAlong and ST_LocateBetween. ERROR: Geometry argument does not have an 'M' ordinate