#1898: Nathan Wagner's patch that adds a flag 2 to allow ST_DelaunayTriangles to dump out a TIN. Just commit and see if winnie has same issue with shp2pgsql-gui checks

git-svn-id: http://svn.osgeo.org/postgis/trunk@11361 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Regina Obe 2013-05-06 06:48:20 +00:00
parent 4d9a4acc9c
commit 40b455d466
5 changed files with 99 additions and 9 deletions

View file

@ -898,12 +898,13 @@ Return a <ulink
url="http://en.wikipedia.org/wiki/Delaunay_triangulation">Delaunay
triangulation</ulink> around the vertices of the input geometry.
Output is a COLLECTION of polygons (for flags=0) or a MULTILINESTRING
(for flags=1). The tolerance, if any, is used to snap input vertices
(for flags=1) or TIN (for flags=2). The tolerance, if any, is used to snap input vertices
togheter.
</para>
<para>Availability: 2.1.0 - requires GEOS &gt;= 3.4.0.</para>
<para>&Z_support;</para>
<para>&T_support;</para>
</refsection>
<refsection>

View file

@ -17,6 +17,8 @@
#include <stdlib.h>
LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
#undef LWGEOM_PROFILE_BUILDAREA
#define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
@ -1298,9 +1300,10 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl
#endif /* POSTGIS_GEOS_VERSION < 32 */
}
LWGEOM*
lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edgeOnly)
{
/*
* output = 1 for edges, 2 for TIN, 0 for polygons
*/
LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
#if POSTGIS_GEOS_VERSION < 34
lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
return NULL;
@ -1308,6 +1311,11 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
GEOSGeometry *g1, *g3;
LWGEOM *lwgeom_result;
if (output < 0 || output > 2) {
lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
return NULL;
}
initGEOS(lwnotice, lwgeom_geos_error);
g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in);
@ -1317,7 +1325,8 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
return NULL;
}
g3 = GEOSDelaunayTriangulation(g1, tolerance, edgeOnly);
/* if output != 1 we want polys */
g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
/* Don't need input geometry anymore */
GEOSGeom_destroy(g1);
@ -1331,12 +1340,21 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
/* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
if (output == 2) {
lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
} else {
lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
}
GEOSGeom_destroy(g3);
if (lwgeom_result == NULL)
{
lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
if (lwgeom_result == NULL) {
if (output != 2) {
lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
} else {
lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
}
return NULL;
}

View file

@ -13,9 +13,12 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "lwgeom_geos.h"
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
@ -23,6 +26,68 @@ LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
return (LWTIN*)lwcollection_add_lwgeom((LWCOLLECTION*)mobj, (LWGEOM*)obj);
}
LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
int type = GEOSGeomTypeId(geom);
int hasZ;
int SRID = GEOSGetSRID(geom);
/* GEOS's 0 is equivalent to our unknown as for SRID values */
if ( SRID == 0 ) SRID = SRID_UNKNOWN;
if ( want3d ) {
hasZ = GEOSHasZ(geom);
if ( ! hasZ ) {
LWDEBUG(3, "Geometry has no Z, won't provide one");
want3d = 0;
}
}
switch (type) {
LWTRIANGLE **geoms;
uint32_t i, ngeoms;
case GEOS_GEOMETRYCOLLECTION:
LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
ngeoms = GEOSGetNumGeometries(geom);
geoms = NULL;
if ( ngeoms ) {
geoms = lwalloc(ngeoms * sizeof *geoms);
if (!geoms) {
lwerror("lwtin_from_geos: can't allocate geoms");
return NULL;
}
for (i=0; i<ngeoms; i++) {
const GEOSGeometry *poly, *ring;
const GEOSCoordSequence *cs;
POINTARRAY *pa;
poly = GEOSGetGeometryN(geom, i);
ring = GEOSGetExteriorRing(poly);
cs = GEOSGeom_getCoordSeq(ring);
pa = ptarray_from_GEOSCoordSeq(cs, want3d);
geoms[i] = lwtriangle_construct(SRID, NULL, pa);
}
}
return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
case GEOS_POLYGON:
case GEOS_MULTIPOINT:
case GEOS_MULTILINESTRING:
case GEOS_MULTIPOLYGON:
case GEOS_LINESTRING:
case GEOS_LINEARRING:
case GEOS_POINT:
lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
break;
default:
lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
return NULL;
}
/* shouldn't get here */
return NULL;
}
void lwtin_free(LWTIN *tin)
{

View file

@ -3,3 +3,6 @@ SELECT 1, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry)
SELECT 2, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry));
SELECT 3, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2));
SELECT 4, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2, 1));
SELECT 5, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry, 0.0, 2));
SELECT 6, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 0.0, 2));
SELECT 7, ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2.0, 2));

View file

@ -2,3 +2,6 @@
2|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,8 9,5 5)),POLYGON((5 5,8 9,7 9,5 5)))
3|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,7 9,5 5)))
4|MULTILINESTRING((5 5,7 9),(5 5,6 0),(6 0,7 9))
5|TIN(((5 5,6 0,7 9,5 5)))
6|TIN(((5 5,6 0,8 9,5 5)),((5 5,8 9,7 9,5 5)))
7|TIN(((5 5,6 0,7 9,5 5)))