From: Paul Ramsey Date: Wed, 22 Aug 2012 21:04:34 +0000 (+0000) Subject: ST_Segmentize(geography, maxseglength) (#1962) X-Git-Tag: 2.1.0beta2~689 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e508329bcbad84320bbb328e29d33d481ae6ab94;p=postgis ST_Segmentize(geography, maxseglength) (#1962) git-svn-id: http://svn.osgeo.org/postgis/trunk@10195 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index 74a70cd04..fa2df93ce 100644 --- a/NEWS +++ b/NEWS @@ -8,6 +8,7 @@ PostGIS 2.1.0 * New Features * + - ST_Segmentize(geography) (Paul Ramsey / OpenGeo) - ST_DelaunayTriangles (Sandro Santilli / Vizzuality) - ST_NearestValue, ST_Neighborhood (Bborie Park / UC Davis) - ST_PixelAsPoint, ST_PixelAsPoints (Bborie Park / UC Davis) diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 239dbb7c1..4a012d62d 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -1137,6 +1137,26 @@ static void test_vector_rotate(void) CU_ASSERT_DOUBLE_EQUAL(n.x, 1.0, 0.00000001); } +static void test_lwgeom_segmentize_sphere(void) +{ + LWGEOM *lwg1, *lwg2; + LWLINE *lwl; + double max = 100000.0 / WGS84_RADIUS; + //char *wkt; + + /* Simple case */ + lwg1 = lwgeom_from_wkt("LINESTRING(0 20, 5 20)", LW_PARSER_CHECK_NONE); + lwg2 = lwgeom_segmentize_sphere(lwg1, max); + lwl = (LWLINE*)lwg2; + //wkt = lwgeom_to_ewkt(lwg2); + CU_ASSERT_EQUAL(lwl->points->npoints, 7); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + //lwfree(wkt); + + return; +} + /* ** Used by test harness to register the tests in this file. */ @@ -1161,6 +1181,7 @@ CU_TestInfo geodetic_tests[] = PG_TEST(test_gbox_utils), PG_TEST(test_vector_angle), PG_TEST(test_vector_rotate), + PG_TEST(test_lwgeom_segmentize_sphere), CU_TEST_INFO_NULL }; CU_SuiteInfo geodetic_suite = {"Geodetic Suite", NULL, NULL, geodetic_tests}; diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 41e8de1d9..ee83055b4 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1064,6 +1064,11 @@ double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT * return d_nearest; } +/** +* Calculate the distance between two edges. +* IMPORTANT: this test does not check for edge intersection!!! (distance == 0) +* You have to check for intersection before calling this function. +*/ double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2) { double d; @@ -1573,6 +1578,172 @@ void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside) } +/** +* Create a new point array with no segment longer than the input segment length (expressed in radians!) +* @param pa_in - input point array pointer +* @param max_seg_length - maximum output segment length in radians +*/ +static POINTARRAY* +ptarray_segmentize_sphere(const POINTARRAY *pa_in, double max_seg_length) +{ + POINTARRAY *pa_out; + int hasz = ptarray_has_z(pa_in); + int hasm = ptarray_has_m(pa_in); + int pa_in_offset = 0; /* input point offset */ + POINT4D p1, p2, p; + POINT3D q1, q2, q, qn; + GEOGRAPHIC_POINT g1, g2, g; + double d; + + /* Just crap out on crazy input */ + if ( ! pa_in ) + lwerror("ptarray_segmentize_sphere: null input pointarray"); + if ( max_seg_length <= 0.0 ) + lwerror("ptarray_segmentize_sphere: maximum segment length must be positive"); + + /* Empty starting array */ + pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints); + + /* Add first point */ + getPoint4d_p(pa_in, pa_in_offset, &p1); + ptarray_append_point(pa_out, &p1, LW_FALSE); + geographic_point_init(p1.x, p1.y, &g1); + pa_in_offset++; + + while ( pa_in_offset < pa_in->npoints ) + { + getPoint4d_p(pa_in, pa_in_offset, &p2); + geographic_point_init(p2.x, p2.y, &g2); + + /* Skip duplicate points (except in case of 2-point lines!) */ + if ( (pa_in->npoints > 2) && p4d_same(&p1, &p2) ) + continue; + + /* How long is this edge? */ + d = sphere_distance(&g1, &g2); + + /* We need to segmentize this edge */ + if ( d > max_seg_length ) + { + int nsegs = 1 + d / max_seg_length; + int i; + double dx, dy, dz, dzz = 0, dmm = 0; + + geog2cart(&g1, &q1); + geog2cart(&g2, &q2); + + dx = (q2.x - q1.x) / nsegs; + dy = (q2.y - q1.y) / nsegs; + dz = (q2.z - q1.z) / nsegs; + + /* The independent Z/M values on the ptarray */ + if ( hasz ) dzz = (p2.z - p1.z) / nsegs; + if ( hasm ) dmm = (p2.m - p1.m) / nsegs; + + q = q1; + p = p1; + + for ( i = 0; i < nsegs - 1; i++ ) + { + /* Move one increment forwards */ + q.x += dx; q.y += dy; q.z += dz; + qn = q; + normalize(&qn); + + /* Back to spherical coordinates */ + cart2geog(&qn, &g); + /* Back to lon/lat */ + p.x = rad2deg(g.lon); + p.y = rad2deg(g.lat); + if ( hasz ) + p.z += dzz; + if ( hasm ) + p.m += dmm; + ptarray_append_point(pa_out, &p, LW_FALSE); + } + + ptarray_append_point(pa_out, &p2, LW_FALSE); + } + /* This edge is already short enough */ + else + { + ptarray_append_point(pa_out, &p2, (pa_in->npoints==2)?LW_TRUE:LW_FALSE); + } + + /* Move one offset forward */ + p1 = p2; + g1 = g2; + pa_in_offset++; + } + + return pa_out; +} + +/** +* Create a new, densified geometry where no segment is longer than max_seg_length. +* Input geometry is not altered, output geometry must be freed by caller. +* @param lwg_in = input geometry +* @param max_seg_length = maximum segment length in radians +*/ +LWGEOM* +lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length) +{ + POINTARRAY *pa_out; + LWLINE *lwline; + LWPOLY *lwpoly_in, *lwpoly_out; + LWCOLLECTION *lwcol_in, *lwcol_out; + int i; + + /* Reflect NULL */ + if ( ! lwg_in ) + return NULL; + + /* Clone empty */ + if ( lwgeom_is_empty(lwg_in) ) + return lwgeom_clone(lwg_in); + + switch (lwg_in->type) + { + case MULTIPOINTTYPE: + case POINTTYPE: + return lwgeom_clone_deep(lwg_in); + break; + case LINETYPE: + lwline = lwgeom_as_lwline(lwg_in); + pa_out = ptarray_segmentize_sphere(lwline->points, max_seg_length); + return lwline_as_lwgeom(lwline_construct(lwg_in->srid, NULL, pa_out)); + break; + case POLYGONTYPE: + lwpoly_in = lwgeom_as_lwpoly(lwg_in); + lwpoly_out = lwpoly_construct_empty(lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in)); + for ( i = 0; i < lwpoly_in->nrings; i++ ) + { + pa_out = ptarray_segmentize_sphere(lwpoly_in->rings[i], max_seg_length); + lwpoly_add_ring(lwpoly_out, pa_out); + } + return lwpoly_as_lwgeom(lwpoly_out); + break; + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + lwcol_in = lwgeom_as_lwcollection(lwg_in); + lwcol_out = lwcollection_construct_empty(lwg_in->type, lwg_in->srid, lwgeom_has_z(lwg_in), lwgeom_has_m(lwg_in)); + for ( i = 0; i < lwcol_in->ngeoms; i++ ) + { + lwcollection_add_lwgeom(lwcol_out, lwgeom_segmentize_sphere(lwcol_in->geoms[i], max_seg_length)); + } + return lwcollection_as_lwgeom(lwcol_out); + break; + default: + lwerror("lwgeom_segmentize_sphere: unsupported input geometry type: %d - %s", + lwg_in->type, lwtype_name(lwg_in->type)); + break; + } + + lwerror("lwgeom_segmentize_sphere got to the end of the function, should not happen"); + return NULL; +} + /** * Returns the area of the ring (ring must be closed) in square radians (surface of @@ -2337,7 +2508,7 @@ int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2) if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) && (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) ) { - lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON and POINT types are currently supported"); + lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON covers POINT tests are currently supported"); return LW_FALSE; } diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index de6de292a..2efdcf53b 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -683,6 +683,13 @@ CREATE OR REPLACE FUNCTION ST_CoveredBy(text, text) $$ SELECT ST_CoveredBy($1::geometry, $2::geometry); $$ LANGUAGE 'sql' IMMUTABLE ; +-- Availability: 2.1.0 +CREATE OR REPLACE FUNCTION ST_Segmentize(geog geography, max_segment_length float8) + RETURNS geography + AS 'MODULE_PATHNAME','geography_segmentize' + LANGUAGE 'c' IMMUTABLE STRICT + COST 100; + -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Intersects(geography, geography) RETURNS boolean diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 76ac3bdd7..f3ec0c106 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -41,6 +41,7 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS); Datum geography_perimeter(PG_FUNCTION_ARGS); Datum geography_project(PG_FUNCTION_ARGS); Datum geography_azimuth(PG_FUNCTION_ARGS); +Datum geography_segmentize(PG_FUNCTION_ARGS); /* ** geography_distance_uncached(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid) @@ -955,3 +956,45 @@ Datum geography_azimuth(PG_FUNCTION_ARGS) } + +/* +** geography_segmentize(GSERIALIZED *g1, double max_seg_length) +** returns densified geometry with no segment longer than max +*/ +PG_FUNCTION_INFO_V1(geography_segmentize); +Datum geography_segmentize(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + double max_seg_length; + uint32_t type1; + + /* Get our geometry object loaded into memory. */ + g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + type1 = gserialized_get_type(g1); + + /* Convert max_seg_length from metric units to radians */ + max_seg_length = PG_GETARG_FLOAT8(1) / WGS84_RADIUS; + + /* We can't densify points or points, reflect them back */ + if ( type1 == POINTTYPE || type1 == MULTIPOINTTYPE || gserialized_is_empty(g1) ) + PG_RETURN_POINTER(g1); + + /* Deserialize */ + lwgeom1 = lwgeom_from_gserialized(g1); + + /* Calculate the densified geometry */ + lwgeom2 = lwgeom_segmentize_sphere(lwgeom1, max_seg_length); + g2 = geography_serialize(lwgeom2); + + /* Clean up */ + lwgeom_free(lwgeom1); + lwgeom_free(lwgeom2); + PG_FREE_IF_COPY(g1, 0); + + PG_RETURN_POINTER(g2); +} + +