/**********************************************************************
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
*
* PostGIS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* PostGIS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PostGIS. If not, see .
*
**********************************************************************
*
* Copyright (C) 2011 Sandro Santilli
*
**********************************************************************/
#include "lwgeom_geos.h"
#include "liblwgeom_internal.h"
#include
#include
static int
lwgeom_ngeoms(const LWGEOM* n)
{
const LWCOLLECTION* c = lwgeom_as_lwcollection(n);
if ( c ) return c->ngeoms;
else return 1;
}
static const LWGEOM*
lwgeom_subgeom(const LWGEOM* g, int n)
{
const LWCOLLECTION* c = lwgeom_as_lwcollection(g);
if ( c ) return lwcollection_getsubgeom((LWCOLLECTION*)c, n);
else return g;
}
static void
lwgeom_collect_endpoints(const LWGEOM* lwg, LWMPOINT* col)
{
int i, n;
LWLINE* l;
switch (lwg->type)
{
case MULTILINETYPE:
for ( i = 0,
n = lwgeom_ngeoms(lwg);
i < n; ++i )
{
lwgeom_collect_endpoints(
lwgeom_subgeom(lwg, i),
col);
}
break;
case LINETYPE:
l = (LWLINE*)lwg;
col = lwmpoint_add_lwpoint(col,
lwline_get_lwpoint(l, 0));
col = lwmpoint_add_lwpoint(col,
lwline_get_lwpoint(l, l->points->npoints-1));
break;
default:
lwerror("lwgeom_collect_endpoints: invalid type %s",
lwtype_name(lwg->type));
break;
}
}
static LWMPOINT*
lwgeom_extract_endpoints(const LWGEOM* lwg)
{
LWMPOINT* col = lwmpoint_construct_empty(SRID_UNKNOWN,
FLAGS_GET_Z(lwg->flags),
FLAGS_GET_M(lwg->flags));
lwgeom_collect_endpoints(lwg, col);
return col;
}
/* Assumes initGEOS was called already */
/* May return LWPOINT or LWMPOINT */
static LWGEOM*
lwgeom_extract_unique_endpoints(const LWGEOM* lwg)
{
LWGEOM* ret;
GEOSGeometry *gepu;
LWMPOINT *epall = lwgeom_extract_endpoints(lwg);
GEOSGeometry *gepall = LWGEOM2GEOS((LWGEOM*)epall, 1);
lwmpoint_free(epall);
if ( ! gepall ) {
lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
return NULL;
}
/* UnaryUnion to remove duplicates */
/* TODO: do it all within pgis using indices */
gepu = GEOSUnaryUnion(gepall);
if ( ! gepu ) {
GEOSGeom_destroy(gepall);
lwerror("GEOSUnaryUnion: %s", lwgeom_geos_errmsg);
return NULL;
}
GEOSGeom_destroy(gepall);
ret = GEOS2LWGEOM(gepu, FLAGS_GET_Z(lwg->flags));
GEOSGeom_destroy(gepu);
if ( ! ret ) {
lwerror("Error during GEOS2LWGEOM");
return NULL;
}
return ret;
}
/* exported */
extern LWGEOM* lwgeom_node(const LWGEOM* lwgeom_in);
LWGEOM*
lwgeom_node(const LWGEOM* lwgeom_in)
{
GEOSGeometry *g1, *gn, *gm;
LWGEOM *ep, *lines;
LWCOLLECTION *col, *tc;
int pn, ln, np, nl;
if ( lwgeom_dimension(lwgeom_in) != 1 ) {
lwerror("Noding geometries of dimension != 1 is unsupported");
return NULL;
}
initGEOS(lwgeom_geos_error, lwgeom_geos_error);
g1 = LWGEOM2GEOS(lwgeom_in, 1);
if ( ! g1 ) {
lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
return NULL;
}
ep = lwgeom_extract_unique_endpoints(lwgeom_in);
if ( ! ep ) {
GEOSGeom_destroy(g1);
lwerror("Error extracting unique endpoints from input");
return NULL;
}
gn = GEOSNode(g1);
GEOSGeom_destroy(g1);
if ( ! gn ) {
lwgeom_free(ep);
lwerror("GEOSNode: %s", lwgeom_geos_errmsg);
return NULL;
}
gm = GEOSLineMerge(gn);
GEOSGeom_destroy(gn);
if ( ! gm ) {
lwgeom_free(ep);
lwerror("GEOSLineMerge: %s", lwgeom_geos_errmsg);
return NULL;
}
lines = GEOS2LWGEOM(gm, FLAGS_GET_Z(lwgeom_in->flags));
GEOSGeom_destroy(gm);
if ( ! lines ) {
lwgeom_free(ep);
lwerror("Error during GEOS2LWGEOM");
return NULL;
}
/*
* Reintroduce endpoints from input, using split-line-by-point.
* Note that by now we can be sure that each point splits at
* most _one_ segment as any point shared by multiple segments
* would already be a node. Also we can be sure that any of
* the segments endpoints won't split any other segment.
* We can use the above 2 assertions to early exit the loop.
*/
col = lwcollection_construct_empty(MULTILINETYPE, lwgeom_in->srid,
FLAGS_GET_Z(lwgeom_in->flags),
FLAGS_GET_M(lwgeom_in->flags));
np = lwgeom_ngeoms(ep);
for (pn=0; pn ln+1) {
tc->geoms[nl] = tc->geoms[nl-1];
--nl;
}
lwgeom_free(tc->geoms[ln]);
tc->geoms[ln] = col->geoms[0];
tc->geoms[ln+1] = col->geoms[1];
tc->ngeoms++;
} else {
lwgeom_free(lines);
/* transfer ownership rather than cloning */
lines = (LWGEOM*)lwcollection_clone_deep(col);
assert(col->ngeoms == 2);
lwgeom_free(col->geoms[0]);
lwgeom_free(col->geoms[1]);
}
/* reset the vector */
assert(col->ngeoms == 2);
col->ngeoms = 0;
break;
}
}
lwgeom_free(ep);
lwcollection_free(col);
lwgeom_set_srid(lines, lwgeom_in->srid);
return (LWGEOM*)lines;
}