1 /**********************************************************************
2  *
3  * rttopo - topology library
4  * http://git.osgeo.org/gitea/rttopo/librttopo
5  *
6  * rttopo is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * rttopo is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with rttopo.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
22  *
23  **********************************************************************/
24 
25 #include "rttopo_config.h"
26 #include "rtgeom_geos.h"
27 #include "librttopo_geom_internal.h"
28 
29 #include <string.h>
30 #include <assert.h>
31 
32 static int
rtgeom_ngeoms(const RTCTX * ctx,const RTGEOM * n)33 rtgeom_ngeoms(const RTCTX *ctx, const RTGEOM* n)
34 {
35   const RTCOLLECTION* c = rtgeom_as_rtcollection(ctx, n);
36   if ( c ) return c->ngeoms;
37   else return 1;
38 }
39 
40 static const RTGEOM*
rtgeom_subgeom(const RTCTX * ctx,const RTGEOM * g,int n)41 rtgeom_subgeom(const RTCTX *ctx, const RTGEOM* g, int n)
42 {
43   const RTCOLLECTION* c = rtgeom_as_rtcollection(ctx, g);
44   if ( c ) return rtcollection_getsubgeom(ctx, (RTCOLLECTION*)c, n);
45   else return g;
46 }
47 
48 
49 static void
rtgeom_collect_endpoints(const RTCTX * ctx,const RTGEOM * rtg,RTMPOINT * col)50 rtgeom_collect_endpoints(const RTCTX *ctx, const RTGEOM* rtg, RTMPOINT* col)
51 {
52   int i, n;
53   RTLINE* l;
54 
55   switch (rtg->type)
56   {
57     case RTMULTILINETYPE:
58       for ( i = 0,
59               n = rtgeom_ngeoms(ctx, rtg);
60             i < n; ++i )
61       {
62         rtgeom_collect_endpoints(ctx,
63           rtgeom_subgeom(ctx, rtg, i),
64           col);
65       }
66       break;
67     case RTLINETYPE:
68       l = (RTLINE*)rtg;
69       col = rtmpoint_add_rtpoint(ctx, col,
70         rtline_get_rtpoint(ctx, l, 0));
71       col = rtmpoint_add_rtpoint(ctx, col,
72         rtline_get_rtpoint(ctx, l, l->points->npoints-1));
73       break;
74     default:
75       rterror(ctx, "rtgeom_collect_endpoints: invalid type %s",
76         rttype_name(ctx, rtg->type));
77       break;
78   }
79 }
80 
81 static RTMPOINT*
rtgeom_extract_endpoints(const RTCTX * ctx,const RTGEOM * rtg)82 rtgeom_extract_endpoints(const RTCTX *ctx, const RTGEOM* rtg)
83 {
84   RTMPOINT* col = rtmpoint_construct_empty(ctx, SRID_UNKNOWN,
85                                 RTFLAGS_GET_Z(rtg->flags),
86                                 RTFLAGS_GET_M(rtg->flags));
87   rtgeom_collect_endpoints(ctx, rtg, col);
88 
89   return col;
90 }
91 
92 /* Assumes rtgeom_geos_ensure_init(ctx) was called */
93 /* May return RTPOINT or RTMPOINT */
94 static RTGEOM*
rtgeom_extract_unique_endpoints(const RTCTX * ctx,const RTGEOM * rtg)95 rtgeom_extract_unique_endpoints(const RTCTX *ctx, const RTGEOM* rtg)
96 {
97 #if RTGEOM_GEOS_VERSION < 33
98   rterror(ctx, "The GEOS version this postgis binary "
99           "was compiled against (%d) doesn't support "
100           "'GEOSUnaryUnion' function _r(ctx->gctx, 3.3.0+ required)",
101           RTGEOM_GEOS_VERSION);
102   return NULL;
103 #else /* RTGEOM_GEOS_VERSION >= 33 */
104   RTGEOM* ret;
105   GEOSGeometry *gepu;
106   RTMPOINT *epall = rtgeom_extract_endpoints(ctx, rtg);
107   GEOSGeometry *gepall = RTGEOM2GEOS(ctx, (RTGEOM*)epall, 1);
108   rtmpoint_free(ctx, epall);
109   if ( ! gepall ) {
110     rterror(ctx, "RTGEOM2GEOS: %s", rtgeom_get_last_geos_error(ctx));
111     return NULL;
112   }
113 
114   /* UnaryUnion to remove duplicates */
115   /* TODO: do it all within pgis using indices */
116   gepu = GEOSUnaryUnion_r(ctx->gctx, gepall);
117   if ( ! gepu ) {
118     GEOSGeom_destroy_r(ctx->gctx, gepall);
119     rterror(ctx, "GEOSUnaryUnion: %s", rtgeom_get_last_geos_error(ctx));
120     return NULL;
121   }
122   GEOSGeom_destroy_r(ctx->gctx, gepall);
123 
124   ret = GEOS2RTGEOM(ctx, gepu, RTFLAGS_GET_Z(rtg->flags));
125   GEOSGeom_destroy_r(ctx->gctx, gepu);
126   if ( ! ret ) {
127     rterror(ctx, "Error during GEOS2RTGEOM");
128     return NULL;
129   }
130 
131   return ret;
132 #endif /* RTGEOM_GEOS_VERSION >= 33 */
133 }
134 
135 /* exported */
136 extern RTGEOM* rtgeom_node(const RTCTX *ctx, const RTGEOM* rtgeom_in);
137 RTGEOM*
rtgeom_node(const RTCTX * ctx,const RTGEOM * rtgeom_in)138 rtgeom_node(const RTCTX *ctx, const RTGEOM* rtgeom_in)
139 {
140 #if RTGEOM_GEOS_VERSION < 33
141   rterror(ctx, "The GEOS version this postgis binary "
142           "was compiled against (%d) doesn't support "
143           "'GEOSUnaryUnion' function _r(ctx->gctx, 3.3.0+ required)",
144           RTGEOM_GEOS_VERSION);
145   return NULL;
146 #else /* RTGEOM_GEOS_VERSION >= 33 */
147   GEOSGeometry *g1, *gu, *gm;
148   RTGEOM *ep, *lines;
149   RTCOLLECTION *col, *tc;
150   int pn, ln, np, nl;
151 
152   if ( rtgeom_dimension(ctx, rtgeom_in) != 1 ) {
153     rterror(ctx, "Noding geometries of dimension != 1 is unsupported");
154     return NULL;
155   }
156 
157   rtgeom_geos_ensure_init(ctx);
158   g1 = RTGEOM2GEOS(ctx, rtgeom_in, 1);
159   if ( ! g1 ) {
160     rterror(ctx, "RTGEOM2GEOS: %s", rtgeom_get_last_geos_error(ctx));
161     return NULL;
162   }
163 
164   ep = rtgeom_extract_unique_endpoints(ctx, rtgeom_in);
165   if ( ! ep ) {
166     GEOSGeom_destroy_r(ctx->gctx, g1);
167     rterror(ctx, "Error extracting unique endpoints from input");
168     return NULL;
169   }
170 
171   /* Unary union input to fully node */
172   gu = GEOSUnaryUnion_r(ctx->gctx, g1);
173   GEOSGeom_destroy_r(ctx->gctx, g1);
174   if ( ! gu ) {
175     rtgeom_free(ctx, ep);
176     rterror(ctx, "GEOSUnaryUnion: %s", rtgeom_get_last_geos_error(ctx));
177     return NULL;
178   }
179 
180   /* Linemerge (in case of overlaps) */
181   gm = GEOSLineMerge_r(ctx->gctx, gu);
182   GEOSGeom_destroy_r(ctx->gctx, gu);
183   if ( ! gm ) {
184     rtgeom_free(ctx, ep);
185     rterror(ctx, "GEOSLineMerge: %s", rtgeom_get_last_geos_error(ctx));
186     return NULL;
187   }
188 
189   lines = GEOS2RTGEOM(ctx, gm, RTFLAGS_GET_Z(rtgeom_in->flags));
190   GEOSGeom_destroy_r(ctx->gctx, gm);
191   if ( ! lines ) {
192     rtgeom_free(ctx, ep);
193     rterror(ctx, "Error during GEOS2RTGEOM");
194     return NULL;
195   }
196 
197   /*
198    * Reintroduce endpoints from input, using split-line-by-point.
199    * Note that by now we can be sure that each point splits at
200    * most _one_ segment as any point shared by multiple segments
201    * would already be a node. Also we can be sure that any of
202    * the segments endpoints won't split any other segment.
203    * We can use the above 2 assertions to early exit the loop.
204    */
205 
206   col = rtcollection_construct_empty(ctx, RTMULTILINETYPE, rtgeom_in->srid,
207                                 RTFLAGS_GET_Z(rtgeom_in->flags),
208                                 RTFLAGS_GET_M(rtgeom_in->flags));
209 
210   np = rtgeom_ngeoms(ctx, ep);
211   for (pn=0; pn<np; ++pn) { /* for each point */
212 
213     const RTPOINT* p = (RTPOINT*)rtgeom_subgeom(ctx, ep, pn);
214 
215     nl = rtgeom_ngeoms(ctx, lines);
216     for (ln=0; ln<nl; ++ln) { /* for each line */
217 
218       const RTLINE* l = (RTLINE*)rtgeom_subgeom(ctx, lines, ln);
219 
220       int s = rtline_split_by_point_to(ctx, l, p, (RTMLINE*)col);
221 
222       if ( ! s ) continue; /* not on this line */
223 
224       if ( s == 1 ) {
225         /* found on this line, but not splitting it */
226         break;
227       }
228 
229       /* splits this line */
230 
231       /* replace this line with the two splits */
232       if ( rtgeom_is_collection(ctx, lines) ) {
233         tc = (RTCOLLECTION*)lines;
234         rtcollection_reserve(ctx, tc, nl + 1);
235         while (nl > ln+1) {
236           tc->geoms[nl] = tc->geoms[nl-1];
237           --nl;
238         }
239         rtgeom_free(ctx, tc->geoms[ln]);
240         tc->geoms[ln]   = col->geoms[0];
241         tc->geoms[ln+1] = col->geoms[1];
242         tc->ngeoms++;
243       } else {
244         rtgeom_free(ctx, lines);
245         /* transfer ownership rather than cloning */
246         lines = (RTGEOM*)rtcollection_clone_deep(ctx, col);
247         assert(col->ngeoms == 2);
248         rtgeom_free(ctx, col->geoms[0]);
249         rtgeom_free(ctx, col->geoms[1]);
250       }
251 
252       /* reset the vector */
253       assert(col->ngeoms == 2);
254       col->ngeoms = 0;
255 
256       break;
257     }
258 
259   }
260 
261   rtgeom_free(ctx, ep);
262   rtcollection_free(ctx, col);
263 
264   lines->srid = rtgeom_in->srid;
265   return (RTGEOM*)lines;
266 #endif /* RTGEOM_GEOS_VERSION >= 33 */
267 }
268 
269