1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS 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 * PostGIS 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 PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright (C) 2015-2020 Sandro Santilli <strk@kbt.io>
22 *
23 **********************************************************************/
24
25
26
27 #include "../postgis_config.h"
28
29 /*#define POSTGIS_DEBUG_LEVEL 1*/
30 #include "lwgeom_log.h"
31
32 #include "liblwgeom_internal.h"
33 #include "liblwgeom_topo_internal.h"
34 #include "lwgeom_geos.h"
35
36 #include <stdio.h>
37 #include <inttypes.h> /* for PRId64 */
38 #include <math.h>
39
40 #ifdef WIN32
41 # define LWTFMT_ELEMID "lld"
42 #else
43 # define LWTFMT_ELEMID PRId64
44 #endif
45
46 /*********************************************************************
47 *
48 * Backend iface
49 *
50 ********************************************************************/
51
lwt_CreateBackendIface(const LWT_BE_DATA * data)52 LWT_BE_IFACE* lwt_CreateBackendIface(const LWT_BE_DATA *data)
53 {
54 LWT_BE_IFACE *iface = lwalloc(sizeof(LWT_BE_IFACE));
55 iface->data = data;
56 iface->cb = NULL;
57 return iface;
58 }
59
lwt_BackendIfaceRegisterCallbacks(LWT_BE_IFACE * iface,const LWT_BE_CALLBACKS * cb)60 void lwt_BackendIfaceRegisterCallbacks(LWT_BE_IFACE *iface,
61 const LWT_BE_CALLBACKS* cb)
62 {
63 iface->cb = cb;
64 }
65
lwt_FreeBackendIface(LWT_BE_IFACE * iface)66 void lwt_FreeBackendIface(LWT_BE_IFACE* iface)
67 {
68 lwfree(iface);
69 }
70
71 /*********************************************************************
72 *
73 * Backend wrappers
74 *
75 ********************************************************************/
76
77 #define CHECKCB(be, method) do { \
78 if ( ! (be)->cb || ! (be)->cb->method ) \
79 lwerror("Callback " # method " not registered by backend"); \
80 } while (0)
81
82 #define CB0(be, method) \
83 CHECKCB(be, method);\
84 return (be)->cb->method((be)->data)
85
86 #define CB1(be, method, a1) \
87 CHECKCB(be, method);\
88 return (be)->cb->method((be)->data, a1)
89
90 #define CBT0(to, method) \
91 CHECKCB((to)->be_iface, method);\
92 return (to)->be_iface->cb->method((to)->be_topo)
93
94 #define CBT1(to, method, a1) \
95 CHECKCB((to)->be_iface, method);\
96 return (to)->be_iface->cb->method((to)->be_topo, a1)
97
98 #define CBT2(to, method, a1, a2) \
99 CHECKCB((to)->be_iface, method);\
100 return (to)->be_iface->cb->method((to)->be_topo, a1, a2)
101
102 #define CBT3(to, method, a1, a2, a3) \
103 CHECKCB((to)->be_iface, method);\
104 return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3)
105
106 #define CBT4(to, method, a1, a2, a3, a4) \
107 CHECKCB((to)->be_iface, method);\
108 return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4)
109
110 #define CBT5(to, method, a1, a2, a3, a4, a5) \
111 CHECKCB((to)->be_iface, method);\
112 return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4, a5)
113
114 #define CBT6(to, method, a1, a2, a3, a4, a5, a6) \
115 CHECKCB((to)->be_iface, method);\
116 return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4, a5, a6)
117
118 const char *
lwt_be_lastErrorMessage(const LWT_BE_IFACE * be)119 lwt_be_lastErrorMessage(const LWT_BE_IFACE* be)
120 {
121 CB0(be, lastErrorMessage);
122 }
123
124 LWT_BE_TOPOLOGY *
lwt_be_loadTopologyByName(LWT_BE_IFACE * be,const char * name)125 lwt_be_loadTopologyByName(LWT_BE_IFACE *be, const char *name)
126 {
127 CB1(be, loadTopologyByName, name);
128 }
129
130 static int
lwt_be_topoGetSRID(LWT_TOPOLOGY * topo)131 lwt_be_topoGetSRID(LWT_TOPOLOGY *topo)
132 {
133 CBT0(topo, topoGetSRID);
134 }
135
136 static double
lwt_be_topoGetPrecision(LWT_TOPOLOGY * topo)137 lwt_be_topoGetPrecision(LWT_TOPOLOGY *topo)
138 {
139 CBT0(topo, topoGetPrecision);
140 }
141
142 static int
lwt_be_topoHasZ(LWT_TOPOLOGY * topo)143 lwt_be_topoHasZ(LWT_TOPOLOGY *topo)
144 {
145 CBT0(topo, topoHasZ);
146 }
147
148 int
lwt_be_freeTopology(LWT_TOPOLOGY * topo)149 lwt_be_freeTopology(LWT_TOPOLOGY *topo)
150 {
151 CBT0(topo, freeTopology);
152 }
153
154 LWT_ISO_NODE *
lwt_be_getNodeById(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields)155 lwt_be_getNodeById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields)
156 {
157 CBT3(topo, getNodeById, ids, numelems, fields);
158 }
159
160 LWT_ISO_NODE *
lwt_be_getNodeWithinDistance2D(LWT_TOPOLOGY * topo,LWPOINT * pt,double dist,uint64_t * numelems,int fields,int64_t limit)161 lwt_be_getNodeWithinDistance2D(LWT_TOPOLOGY *topo,
162 LWPOINT *pt,
163 double dist,
164 uint64_t *numelems,
165 int fields,
166 int64_t limit)
167 {
168 CBT5(topo, getNodeWithinDistance2D, pt, dist, numelems, fields, limit);
169 }
170
171 static LWT_ISO_NODE *
lwt_be_getNodeWithinBox2D(const LWT_TOPOLOGY * topo,const GBOX * box,uint64_t * numelems,int fields,uint64_t limit)172 lwt_be_getNodeWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit)
173 {
174 CBT4(topo, getNodeWithinBox2D, box, numelems, fields, limit);
175 }
176
177 static LWT_ISO_EDGE *
lwt_be_getEdgeWithinBox2D(const LWT_TOPOLOGY * topo,const GBOX * box,uint64_t * numelems,int fields,uint64_t limit)178 lwt_be_getEdgeWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit)
179 {
180 CBT4(topo, getEdgeWithinBox2D, box, numelems, fields, limit);
181 }
182
183 static LWT_ISO_FACE *
lwt_be_getFaceWithinBox2D(const LWT_TOPOLOGY * topo,const GBOX * box,uint64_t * numelems,int fields,uint64_t limit)184 lwt_be_getFaceWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit)
185 {
186 CBT4(topo, getFaceWithinBox2D, box, numelems, fields, limit);
187 }
188
189 int
lwt_be_insertNodes(LWT_TOPOLOGY * topo,LWT_ISO_NODE * node,uint64_t numelems)190 lwt_be_insertNodes(LWT_TOPOLOGY *topo, LWT_ISO_NODE *node, uint64_t numelems)
191 {
192 CBT2(topo, insertNodes, node, numelems);
193 }
194
195 static int
lwt_be_insertFaces(LWT_TOPOLOGY * topo,LWT_ISO_FACE * face,uint64_t numelems)196 lwt_be_insertFaces(LWT_TOPOLOGY *topo, LWT_ISO_FACE *face, uint64_t numelems)
197 {
198 CBT2(topo, insertFaces, face, numelems);
199 }
200
201 static int
lwt_be_deleteFacesById(const LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t numelems)202 lwt_be_deleteFacesById(const LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t numelems)
203 {
204 CBT2(topo, deleteFacesById, ids, numelems);
205 }
206
207 static int
lwt_be_deleteNodesById(const LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t numelems)208 lwt_be_deleteNodesById(const LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t numelems)
209 {
210 CBT2(topo, deleteNodesById, ids, numelems);
211 }
212
213 LWT_ELEMID
lwt_be_getNextEdgeId(LWT_TOPOLOGY * topo)214 lwt_be_getNextEdgeId(LWT_TOPOLOGY* topo)
215 {
216 CBT0(topo, getNextEdgeId);
217 }
218
219 LWT_ISO_EDGE *
lwt_be_getEdgeById(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields)220 lwt_be_getEdgeById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields)
221 {
222 CBT3(topo, getEdgeById, ids, numelems, fields);
223 }
224
225 static LWT_ISO_FACE *
lwt_be_getFaceById(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields)226 lwt_be_getFaceById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields)
227 {
228 CBT3(topo, getFaceById, ids, numelems, fields);
229 }
230
231 static LWT_ISO_EDGE *
lwt_be_getEdgeByNode(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields)232 lwt_be_getEdgeByNode(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields)
233 {
234 CBT3(topo, getEdgeByNode, ids, numelems, fields);
235 }
236
237 static LWT_ISO_EDGE *
lwt_be_getEdgeByFace(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields,const GBOX * box)238 lwt_be_getEdgeByFace(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields, const GBOX *box)
239 {
240 CBT4(topo, getEdgeByFace, ids, numelems, fields, box);
241 }
242
243 static LWT_ISO_NODE *
lwt_be_getNodeByFace(LWT_TOPOLOGY * topo,const LWT_ELEMID * ids,uint64_t * numelems,int fields,const GBOX * box)244 lwt_be_getNodeByFace(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields, const GBOX *box)
245 {
246 CBT4(topo, getNodeByFace, ids, numelems, fields, box);
247 }
248
249 LWT_ISO_EDGE *
lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY * topo,LWPOINT * pt,double dist,uint64_t * numelems,int fields,int64_t limit)250 lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY *topo,
251 LWPOINT *pt,
252 double dist,
253 uint64_t *numelems,
254 int fields,
255 int64_t limit)
256 {
257 CBT5(topo, getEdgeWithinDistance2D, pt, dist, numelems, fields, limit);
258 }
259
260 int
lwt_be_insertEdges(LWT_TOPOLOGY * topo,LWT_ISO_EDGE * edge,uint64_t numelems)261 lwt_be_insertEdges(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edge, uint64_t numelems)
262 {
263 CBT2(topo, insertEdges, edge, numelems);
264 }
265
266 int
lwt_be_updateEdges(LWT_TOPOLOGY * topo,const LWT_ISO_EDGE * sel_edge,int sel_fields,const LWT_ISO_EDGE * upd_edge,int upd_fields,const LWT_ISO_EDGE * exc_edge,int exc_fields)267 lwt_be_updateEdges(LWT_TOPOLOGY* topo,
268 const LWT_ISO_EDGE* sel_edge, int sel_fields,
269 const LWT_ISO_EDGE* upd_edge, int upd_fields,
270 const LWT_ISO_EDGE* exc_edge, int exc_fields
271 )
272 {
273 CBT6(topo, updateEdges, sel_edge, sel_fields,
274 upd_edge, upd_fields,
275 exc_edge, exc_fields);
276 }
277
278 static int
lwt_be_updateNodes(LWT_TOPOLOGY * topo,const LWT_ISO_NODE * sel_node,int sel_fields,const LWT_ISO_NODE * upd_node,int upd_fields,const LWT_ISO_NODE * exc_node,int exc_fields)279 lwt_be_updateNodes(LWT_TOPOLOGY* topo,
280 const LWT_ISO_NODE* sel_node, int sel_fields,
281 const LWT_ISO_NODE* upd_node, int upd_fields,
282 const LWT_ISO_NODE* exc_node, int exc_fields
283 )
284 {
285 CBT6(topo, updateNodes, sel_node, sel_fields,
286 upd_node, upd_fields,
287 exc_node, exc_fields);
288 }
289
290 static uint64_t
lwt_be_updateFacesById(LWT_TOPOLOGY * topo,const LWT_ISO_FACE * faces,uint64_t numfaces)291 lwt_be_updateFacesById(LWT_TOPOLOGY* topo,
292 const LWT_ISO_FACE* faces, uint64_t numfaces
293 )
294 {
295 CBT2(topo, updateFacesById, faces, numfaces);
296 }
297
298 static int
lwt_be_updateEdgesById(LWT_TOPOLOGY * topo,const LWT_ISO_EDGE * edges,int numedges,int upd_fields)299 lwt_be_updateEdgesById(LWT_TOPOLOGY* topo,
300 const LWT_ISO_EDGE* edges, int numedges, int upd_fields
301 )
302 {
303 CBT3(topo, updateEdgesById, edges, numedges, upd_fields);
304 }
305
306 static int
lwt_be_updateNodesById(LWT_TOPOLOGY * topo,const LWT_ISO_NODE * nodes,int numnodes,int upd_fields)307 lwt_be_updateNodesById(LWT_TOPOLOGY* topo,
308 const LWT_ISO_NODE* nodes, int numnodes, int upd_fields
309 )
310 {
311 CBT3(topo, updateNodesById, nodes, numnodes, upd_fields);
312 }
313
314 int
lwt_be_deleteEdges(LWT_TOPOLOGY * topo,const LWT_ISO_EDGE * sel_edge,int sel_fields)315 lwt_be_deleteEdges(LWT_TOPOLOGY* topo,
316 const LWT_ISO_EDGE* sel_edge, int sel_fields
317 )
318 {
319 CBT2(topo, deleteEdges, sel_edge, sel_fields);
320 }
321
322 LWT_ELEMID
lwt_be_getFaceContainingPoint(LWT_TOPOLOGY * topo,LWPOINT * pt)323 lwt_be_getFaceContainingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt)
324 {
325 CBT1(topo, getFaceContainingPoint, pt);
326 }
327
328
329 int
lwt_be_updateTopoGeomEdgeSplit(LWT_TOPOLOGY * topo,LWT_ELEMID split_edge,LWT_ELEMID new_edge1,LWT_ELEMID new_edge2)330 lwt_be_updateTopoGeomEdgeSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_edge, LWT_ELEMID new_edge1, LWT_ELEMID new_edge2)
331 {
332 CBT3(topo, updateTopoGeomEdgeSplit, split_edge, new_edge1, new_edge2);
333 }
334
335 static int
lwt_be_updateTopoGeomFaceSplit(LWT_TOPOLOGY * topo,LWT_ELEMID split_face,LWT_ELEMID new_face1,LWT_ELEMID new_face2)336 lwt_be_updateTopoGeomFaceSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_face,
337 LWT_ELEMID new_face1, LWT_ELEMID new_face2)
338 {
339 CBT3(topo, updateTopoGeomFaceSplit, split_face, new_face1, new_face2);
340 }
341
342 static int
lwt_be_checkTopoGeomRemEdge(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id,LWT_ELEMID face_left,LWT_ELEMID face_right)343 lwt_be_checkTopoGeomRemEdge(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id,
344 LWT_ELEMID face_left, LWT_ELEMID face_right)
345 {
346 CBT3(topo, checkTopoGeomRemEdge, edge_id, face_left, face_right);
347 }
348
349 static int
lwt_be_checkTopoGeomRemNode(LWT_TOPOLOGY * topo,LWT_ELEMID node_id,LWT_ELEMID eid1,LWT_ELEMID eid2)350 lwt_be_checkTopoGeomRemNode(LWT_TOPOLOGY* topo, LWT_ELEMID node_id,
351 LWT_ELEMID eid1, LWT_ELEMID eid2)
352 {
353 CBT3(topo, checkTopoGeomRemNode, node_id, eid1, eid2);
354 }
355
356 static int
lwt_be_updateTopoGeomFaceHeal(LWT_TOPOLOGY * topo,LWT_ELEMID face1,LWT_ELEMID face2,LWT_ELEMID newface)357 lwt_be_updateTopoGeomFaceHeal(LWT_TOPOLOGY* topo,
358 LWT_ELEMID face1, LWT_ELEMID face2,
359 LWT_ELEMID newface)
360 {
361 CBT3(topo, updateTopoGeomFaceHeal, face1, face2, newface);
362 }
363
364 static int
lwt_be_updateTopoGeomEdgeHeal(LWT_TOPOLOGY * topo,LWT_ELEMID edge1,LWT_ELEMID edge2,LWT_ELEMID newedge)365 lwt_be_updateTopoGeomEdgeHeal(LWT_TOPOLOGY* topo,
366 LWT_ELEMID edge1, LWT_ELEMID edge2,
367 LWT_ELEMID newedge)
368 {
369 CBT3(topo, updateTopoGeomEdgeHeal, edge1, edge2, newedge);
370 }
371
372 static LWT_ELEMID *
lwt_be_getRingEdges(LWT_TOPOLOGY * topo,LWT_ELEMID edge,uint64_t * numedges,uint64_t limit)373 lwt_be_getRingEdges(LWT_TOPOLOGY *topo, LWT_ELEMID edge, uint64_t *numedges, uint64_t limit)
374 {
375 CBT3(topo, getRingEdges, edge, numedges, limit);
376 }
377
378
379 /* wrappers of backend wrappers... */
380
381 int
lwt_be_ExistsCoincidentNode(LWT_TOPOLOGY * topo,LWPOINT * pt)382 lwt_be_ExistsCoincidentNode(LWT_TOPOLOGY* topo, LWPOINT* pt)
383 {
384 uint64_t exists = 0;
385 lwt_be_getNodeWithinDistance2D(topo, pt, 0, &exists, 0, -1);
386 if (exists == UINT64_MAX)
387 {
388 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
389 return 0;
390 }
391 return exists;
392 }
393
394 int
lwt_be_ExistsEdgeIntersectingPoint(LWT_TOPOLOGY * topo,LWPOINT * pt)395 lwt_be_ExistsEdgeIntersectingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt)
396 {
397 uint64_t exists = 0;
398 lwt_be_getEdgeWithinDistance2D(topo, pt, 0, &exists, 0, -1);
399 if (exists == UINT64_MAX)
400 {
401 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
402 return 0;
403 }
404 return exists;
405 }
406
407 /************************************************************************
408 *
409 * Utility functions
410 *
411 ************************************************************************/
412
413 static LWGEOM *
_lwt_toposnap(LWGEOM * src,LWGEOM * tgt,double tol)414 _lwt_toposnap(LWGEOM *src, LWGEOM *tgt, double tol)
415 {
416 LWGEOM *tmp = src;
417 LWGEOM *tmp2;
418 int changed;
419 int iterations = 0;
420
421 int maxiterations = lwgeom_count_vertices(tgt);
422
423 /* GEOS snapping can be unstable */
424 /* See https://trac.osgeo.org/geos/ticket/760 */
425 do {
426 tmp2 = lwgeom_snap(tmp, tgt, tol);
427 ++iterations;
428 changed = ( lwgeom_count_vertices(tmp2) != lwgeom_count_vertices(tmp) );
429 LWDEBUGF(2, "After iteration %d, geometry changed ? %d (%d vs %d vertices)", iterations, changed, lwgeom_count_vertices(tmp2), lwgeom_count_vertices(tmp));
430 if ( tmp != src ) lwgeom_free(tmp);
431 tmp = tmp2;
432 } while ( changed && iterations <= maxiterations );
433
434 LWDEBUGF(1, "It took %d/%d iterations to properly snap",
435 iterations, maxiterations);
436
437 return tmp;
438 }
439
440 static void
_lwt_release_faces(LWT_ISO_FACE * faces,int num_faces)441 _lwt_release_faces(LWT_ISO_FACE *faces, int num_faces)
442 {
443 int i;
444 for ( i=0; i<num_faces; ++i ) {
445 if ( faces[i].mbr ) lwfree(faces[i].mbr);
446 }
447 lwfree(faces);
448 }
449
450 static void
_lwt_release_edges(LWT_ISO_EDGE * edges,int num_edges)451 _lwt_release_edges(LWT_ISO_EDGE *edges, int num_edges)
452 {
453 int i;
454 for ( i=0; i<num_edges; ++i ) {
455 if ( edges[i].geom ) lwline_free(edges[i].geom);
456 }
457 lwfree(edges);
458 }
459
460 static void
_lwt_release_nodes(LWT_ISO_NODE * nodes,int num_nodes)461 _lwt_release_nodes(LWT_ISO_NODE *nodes, int num_nodes)
462 {
463 int i;
464 for ( i=0; i<num_nodes; ++i ) {
465 if ( nodes[i].geom ) lwpoint_free(nodes[i].geom);
466 }
467 lwfree(nodes);
468 }
469
470 /************************************************************************
471 *
472 * API implementation
473 *
474 ************************************************************************/
475
476 LWT_TOPOLOGY *
lwt_LoadTopology(LWT_BE_IFACE * iface,const char * name)477 lwt_LoadTopology( LWT_BE_IFACE *iface, const char *name )
478 {
479 LWT_BE_TOPOLOGY* be_topo;
480 LWT_TOPOLOGY* topo;
481
482 be_topo = lwt_be_loadTopologyByName(iface, name);
483 if ( ! be_topo ) {
484 //lwerror("Could not load topology from backend: %s",
485 lwerror("%s", lwt_be_lastErrorMessage(iface));
486 return NULL;
487 }
488 topo = lwalloc(sizeof(LWT_TOPOLOGY));
489 topo->be_iface = iface;
490 topo->be_topo = be_topo;
491 topo->srid = lwt_be_topoGetSRID(topo);
492 topo->hasZ = lwt_be_topoHasZ(topo);
493 topo->precision = lwt_be_topoGetPrecision(topo);
494
495 return topo;
496 }
497
498 void
lwt_FreeTopology(LWT_TOPOLOGY * topo)499 lwt_FreeTopology( LWT_TOPOLOGY* topo )
500 {
501 if ( ! lwt_be_freeTopology(topo) ) {
502 lwnotice("Could not release backend topology memory: %s",
503 lwt_be_lastErrorMessage(topo->be_iface));
504 }
505 lwfree(topo);
506 }
507
508 /**
509 * @param checkFace if non zero will check the given face
510 * for really containing the point or determine the
511 * face when given face is -1. Use 0 to simply use
512 * the given face value, no matter what (effectively
513 * allowing adding a non-isolated point when used
514 * with face=-1).
515 */
516 static LWT_ELEMID
_lwt_AddIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID face,LWPOINT * pt,int skipISOChecks,int checkFace)517 _lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face,
518 LWPOINT* pt, int skipISOChecks, int checkFace )
519 {
520 LWT_ELEMID foundInFace = -1;
521
522 if ( lwpoint_is_empty(pt) )
523 {
524 lwerror("Cannot add empty point as isolated node");
525 return -1;
526 }
527
528
529 if ( ! skipISOChecks )
530 {
531 if ( lwt_be_ExistsCoincidentNode(topo, pt) ) /*x*/
532 {
533 lwerror("SQL/MM Spatial exception - coincident node");
534 return -1;
535 }
536 if ( lwt_be_ExistsEdgeIntersectingPoint(topo, pt) ) /*x*/
537 {
538 lwerror("SQL/MM Spatial exception - edge crosses node.");
539 return -1;
540 }
541 }
542
543 if ( checkFace && ( face == -1 || ! skipISOChecks ) )
544 {
545 foundInFace = lwt_be_getFaceContainingPoint(topo, pt); /*x*/
546 if ( foundInFace == -2 ) {
547 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
548 return -1;
549 }
550 if ( foundInFace == -1 ) foundInFace = 0;
551 }
552
553 if ( face == -1 ) {
554 face = foundInFace;
555 }
556 else if ( ! skipISOChecks && foundInFace != face ) {
557 #if 0
558 lwerror("SQL/MM Spatial exception - within face %d (not %d)",
559 foundInFace, face);
560 #else
561 lwerror("SQL/MM Spatial exception - not within face");
562 #endif
563 return -1;
564 }
565
566 LWT_ISO_NODE node;
567 node.node_id = -1;
568 node.containing_face = face;
569 node.geom = pt;
570 if ( ! lwt_be_insertNodes(topo, &node, 1) )
571 {
572 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
573 return -1;
574 }
575
576 return node.node_id;
577 }
578
579 LWT_ELEMID
lwt_AddIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID face,LWPOINT * pt,int skipISOChecks)580 lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face,
581 LWPOINT* pt, int skipISOChecks )
582 {
583 return _lwt_AddIsoNode( topo, face, pt, skipISOChecks, 1 );
584 }
585
586 /* Check that an edge does not cross an existing node or edge
587 *
588 * @param myself the id of an edge to skip, if any
589 * (for ChangeEdgeGeom). Can use 0 for none.
590 *
591 * Return -1 on cross or error, 0 if everything is fine.
592 * Note that before returning -1, lwerror is invoked...
593 */
594 static int
_lwt_CheckEdgeCrossing(LWT_TOPOLOGY * topo,LWT_ELEMID start_node,LWT_ELEMID end_node,const LWLINE * geom,LWT_ELEMID myself)595 _lwt_CheckEdgeCrossing( LWT_TOPOLOGY* topo,
596 LWT_ELEMID start_node, LWT_ELEMID end_node,
597 const LWLINE *geom, LWT_ELEMID myself )
598 {
599 uint64_t i, num_nodes, num_edges;
600 LWT_ISO_EDGE *edges;
601 LWT_ISO_NODE *nodes;
602 const GBOX *edgebox;
603 GEOSGeometry *edgegg;
604
605 initGEOS(lwnotice, lwgeom_geos_error);
606
607 edgegg = LWGEOM2GEOS(lwline_as_lwgeom(geom), 0);
608 if (!edgegg)
609 {
610 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
611 return -1;
612 }
613 edgebox = lwgeom_get_bbox( lwline_as_lwgeom(geom) );
614
615 /* loop over each node within the edge's gbox */
616 nodes = lwt_be_getNodeWithinBox2D( topo, edgebox, &num_nodes,
617 LWT_COL_NODE_ALL, 0 );
618 LWDEBUGF(1, "lwt_be_getNodeWithinBox2D returned %d nodes", num_nodes);
619 if (num_nodes == UINT64_MAX)
620 {
621 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
622 return -1;
623 }
624 for ( i=0; i<num_nodes; ++i )
625 {
626 const POINT2D *pt;
627 LWT_ISO_NODE* node = &(nodes[i]);
628 if ( node->node_id == start_node ) continue;
629 if ( node->node_id == end_node ) continue;
630 /* check if the edge contains this node (not on boundary) */
631 /* ST_RelateMatch(rec.relate, 'T********') */
632 pt = getPoint2d_cp(node->geom->point, 0);
633 int contains = ptarray_contains_point_partial(geom->points, pt, LW_FALSE, NULL) == LW_BOUNDARY;
634 if ( contains )
635 {
636 GEOSGeom_destroy(edgegg);
637 _lwt_release_nodes(nodes, num_nodes);
638 lwerror("SQL/MM Spatial exception - geometry crosses a node");
639 return -1;
640 }
641 }
642 if ( nodes ) _lwt_release_nodes(nodes, num_nodes);
643 /* may be NULL if num_nodes == 0 */
644
645 /* loop over each edge within the edge's gbox */
646 edges = lwt_be_getEdgeWithinBox2D( topo, edgebox, &num_edges, LWT_COL_EDGE_ALL, 0 );
647 LWDEBUGF(1, "lwt_be_getEdgeWithinBox2D returned %d edges", num_edges);
648 if (num_edges == UINT64_MAX)
649 {
650 GEOSGeom_destroy(edgegg);
651 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
652 return -1;
653 }
654 for ( i=0; i<num_edges; ++i )
655 {
656 LWT_ISO_EDGE* edge = &(edges[i]);
657 LWT_ELEMID edge_id = edge->edge_id;
658 GEOSGeometry *eegg;
659 char *relate;
660 int match;
661
662 if ( edge_id == myself ) continue;
663
664 if ( ! edge->geom ) {
665 _lwt_release_edges(edges, num_edges);
666 lwerror("Edge %d has NULL geometry!", edge_id);
667 return -1;
668 }
669
670 eegg = LWGEOM2GEOS( lwline_as_lwgeom(edge->geom), 0 );
671 if ( ! eegg ) {
672 GEOSGeom_destroy(edgegg);
673 _lwt_release_edges(edges, num_edges);
674 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
675 return -1;
676 }
677
678 LWDEBUGF(2, "Edge %d converted to GEOS", edge_id);
679
680 /* check if the edge crosses our edge (not boundary-boundary) */
681
682 relate = GEOSRelateBoundaryNodeRule(eegg, edgegg, 2);
683 if ( ! relate ) {
684 GEOSGeom_destroy(eegg);
685 GEOSGeom_destroy(edgegg);
686 _lwt_release_edges(edges, num_edges);
687 lwerror("GEOSRelateBoundaryNodeRule error: %s", lwgeom_geos_errmsg);
688 return -1;
689 }
690
691 LWDEBUGF(2, "Edge %d relate pattern is %s", edge_id, relate);
692
693 match = GEOSRelatePatternMatch(relate, "F********");
694 if ( match ) {
695 /* error or no interior intersection */
696 GEOSGeom_destroy(eegg);
697 GEOSFree(relate);
698 if ( match == 2 ) {
699 _lwt_release_edges(edges, num_edges);
700 GEOSGeom_destroy(edgegg);
701 lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg);
702 return -1;
703 }
704 else continue; /* no interior intersection */
705 }
706
707 match = GEOSRelatePatternMatch(relate, "1FFF*FFF2");
708 if ( match ) {
709 _lwt_release_edges(edges, num_edges);
710 GEOSGeom_destroy(edgegg);
711 GEOSGeom_destroy(eegg);
712 GEOSFree(relate);
713 if ( match == 2 ) {
714 lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg);
715 } else {
716 lwerror("SQL/MM Spatial exception - coincident edge %" LWTFMT_ELEMID,
717 edge_id);
718 }
719 return -1;
720 }
721
722 match = GEOSRelatePatternMatch(relate, "1********");
723 if ( match ) {
724 _lwt_release_edges(edges, num_edges);
725 GEOSGeom_destroy(edgegg);
726 GEOSGeom_destroy(eegg);
727 GEOSFree(relate);
728 if ( match == 2 ) {
729 lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg);
730 } else {
731 lwerror("Spatial exception - geometry intersects edge %"
732 LWTFMT_ELEMID, edge_id);
733 }
734 return -1;
735 }
736
737 match = GEOSRelatePatternMatch(relate, "T********");
738 if ( match ) {
739 _lwt_release_edges(edges, num_edges);
740 GEOSGeom_destroy(edgegg);
741 GEOSGeom_destroy(eegg);
742 GEOSFree(relate);
743 if ( match == 2 ) {
744 lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg);
745 } else {
746 lwerror("SQL/MM Spatial exception - geometry crosses edge %"
747 LWTFMT_ELEMID, edge_id);
748 }
749 return -1;
750 }
751
752 LWDEBUGF(2, "Edge %d analisys completed, it does no harm", edge_id);
753
754 GEOSFree(relate);
755 GEOSGeom_destroy(eegg);
756 }
757 if ( edges ) _lwt_release_edges(edges, num_edges);
758 /* would be NULL if num_edges was 0 */
759
760 GEOSGeom_destroy(edgegg);
761
762 return 0;
763 }
764
765
766 LWT_ELEMID
lwt_AddIsoEdge(LWT_TOPOLOGY * topo,LWT_ELEMID startNode,LWT_ELEMID endNode,const LWLINE * geom)767 lwt_AddIsoEdge( LWT_TOPOLOGY* topo, LWT_ELEMID startNode,
768 LWT_ELEMID endNode, const LWLINE* geom )
769 {
770 uint64_t num_nodes;
771 uint64_t i;
772 LWT_ISO_EDGE newedge;
773 LWT_ISO_NODE *endpoints;
774 LWT_ELEMID containing_face = -1;
775 LWT_ELEMID node_ids[2];
776 LWT_ISO_NODE updated_nodes[2];
777 int skipISOChecks = 0;
778 POINT2D p1, p2;
779
780 /* NOT IN THE SPECS:
781 * A closed edge is never isolated (as it forms a face)
782 */
783 if (startNode == endNode)
784 {
785 lwerror("Closed edges would not be isolated, try lwt_AddEdgeNewFaces");
786 return -1;
787 }
788
789 if ( ! skipISOChecks )
790 {
791 /* Acurve must be simple */
792 if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) )
793 {
794 lwerror("SQL/MM Spatial exception - curve not simple");
795 return -1;
796 }
797 }
798
799 /*
800 * Check for:
801 * existence of nodes
802 * nodes faces match
803 * Extract:
804 * nodes face id
805 * nodes geoms
806 */
807 num_nodes = 2;
808 node_ids[0] = startNode;
809 node_ids[1] = endNode;
810 endpoints = lwt_be_getNodeById( topo, node_ids, &num_nodes,
811 LWT_COL_NODE_ALL );
812 if (num_nodes == UINT64_MAX)
813 {
814 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
815 return -1;
816 }
817 else if ( num_nodes < 2 )
818 {
819 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
820 lwerror("SQL/MM Spatial exception - non-existent node");
821 return -1;
822 }
823 for ( i=0; i<num_nodes; ++i )
824 {
825 const LWT_ISO_NODE *n = &(endpoints[i]);
826 if ( n->containing_face == -1 )
827 {
828 _lwt_release_nodes(endpoints, num_nodes);
829 lwerror("SQL/MM Spatial exception - not isolated node");
830 return -1;
831 }
832 if ( containing_face == -1 ) containing_face = n->containing_face;
833 else if ( containing_face != n->containing_face )
834 {
835 _lwt_release_nodes(endpoints, num_nodes);
836 lwerror("SQL/MM Spatial exception - nodes in different faces");
837 return -1;
838 }
839
840 if ( ! skipISOChecks )
841 {
842 if ( n->node_id == startNode )
843 {
844 /* l) Check that start point of acurve match start node geoms. */
845 getPoint2d_p(geom->points, 0, &p1);
846 getPoint2d_p(n->geom->point, 0, &p2);
847 if ( ! p2d_same(&p1, &p2) )
848 {
849 _lwt_release_nodes(endpoints, num_nodes);
850 lwerror("SQL/MM Spatial exception - "
851 "start node not geometry start point.");
852 return -1;
853 }
854 }
855 else
856 {
857 /* m) Check that end point of acurve match end node geoms. */
858 getPoint2d_p(geom->points, geom->points->npoints-1, &p1);
859 getPoint2d_p(n->geom->point, 0, &p2);
860 if ( ! p2d_same(&p1, &p2) )
861 {
862 _lwt_release_nodes(endpoints, num_nodes);
863 lwerror("SQL/MM Spatial exception - "
864 "end node not geometry end point.");
865 return -1;
866 }
867 }
868 }
869 }
870
871 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
872
873 if ( ! skipISOChecks )
874 {
875 if ( _lwt_CheckEdgeCrossing( topo, startNode, endNode, geom, 0 ) )
876 {
877 /* would have called lwerror already, leaking :( */
878 return -1;
879 }
880 }
881
882 /*
883 * All checks passed, time to prepare the new edge
884 */
885
886 newedge.edge_id = lwt_be_getNextEdgeId( topo );
887 if ( newedge.edge_id == -1 ) {
888 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
889 return -1;
890 }
891
892 /* TODO: this should likely be an exception instead ! */
893 if ( containing_face == -1 ) containing_face = 0;
894
895 newedge.start_node = startNode;
896 newedge.end_node = endNode;
897 newedge.face_left = newedge.face_right = containing_face;
898 newedge.next_left = -newedge.edge_id;
899 newedge.next_right = newedge.edge_id;
900 newedge.geom = (LWLINE *)geom; /* const cast.. */
901
902 int ret = lwt_be_insertEdges(topo, &newedge, 1);
903 if ( ret == -1 ) {
904 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
905 return -1;
906 } else if ( ret == 0 ) {
907 lwerror("Insertion of split edge failed (no reason)");
908 return -1;
909 }
910
911 /*
912 * Update Node containing_face values
913 *
914 * the nodes anode and anothernode are no more isolated
915 * because now there is an edge connecting them
916 */
917 updated_nodes[0].node_id = startNode;
918 updated_nodes[0].containing_face = -1;
919 updated_nodes[1].node_id = endNode;
920 updated_nodes[1].containing_face = -1;
921 ret = lwt_be_updateNodesById(topo, updated_nodes, 2,
922 LWT_COL_NODE_CONTAINING_FACE);
923 if ( ret == -1 ) {
924 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
925 return -1;
926 }
927
928 return newedge.edge_id;
929 }
930
931 static LWCOLLECTION *
_lwt_EdgeSplit(LWT_TOPOLOGY * topo,LWT_ELEMID edge,LWPOINT * pt,int skipISOChecks,LWT_ISO_EDGE ** oldedge)932 _lwt_EdgeSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge, LWPOINT* pt, int skipISOChecks, LWT_ISO_EDGE** oldedge )
933 {
934 LWGEOM *split;
935 LWCOLLECTION *split_col;
936 uint64_t i;
937
938 /* Get edge */
939 i = 1;
940 LWDEBUG(1, "calling lwt_be_getEdgeById");
941 *oldedge = lwt_be_getEdgeById(topo, &edge, &i, LWT_COL_EDGE_ALL);
942 LWDEBUGF(1, "lwt_be_getEdgeById returned %p", *oldedge);
943 if ( ! *oldedge )
944 {
945 LWDEBUGF(1, "lwt_be_getEdgeById returned NULL and set i=%d", i);
946 if (i == UINT64_MAX)
947 {
948 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
949 return NULL;
950 }
951 else if ( i == 0 )
952 {
953 lwerror("SQL/MM Spatial exception - non-existent edge");
954 return NULL;
955 }
956 else
957 {
958 lwerror("Backend coding error: getEdgeById callback returned NULL "
959 "but numelements output parameter has value %d "
960 "(expected 0 or 1)", i);
961 return NULL;
962 }
963 }
964
965
966 /*
967 * - check if a coincident node already exists
968 */
969 if ( ! skipISOChecks )
970 {
971 LWDEBUG(1, "calling lwt_be_ExistsCoincidentNode");
972 if ( lwt_be_ExistsCoincidentNode(topo, pt) ) /*x*/
973 {
974 LWDEBUG(1, "lwt_be_ExistsCoincidentNode returned");
975 _lwt_release_edges(*oldedge, 1);
976 lwerror("SQL/MM Spatial exception - coincident node");
977 return NULL;
978 }
979 LWDEBUG(1, "lwt_be_ExistsCoincidentNode returned");
980 }
981
982 /* Split edge */
983 split = lwgeom_split((LWGEOM*)(*oldedge)->geom, (LWGEOM*)pt);
984 if ( ! split )
985 {
986 _lwt_release_edges(*oldedge, 1);
987 lwerror("could not split edge by point ?");
988 return NULL;
989 }
990 split_col = lwgeom_as_lwcollection(split);
991 if ( ! split_col ) {
992 _lwt_release_edges(*oldedge, 1);
993 lwgeom_free(split);
994 lwerror("lwgeom_as_lwcollection returned NULL");
995 return NULL;
996 }
997 if (split_col->ngeoms < 2) {
998 _lwt_release_edges(*oldedge, 1);
999 lwgeom_free(split);
1000 lwerror("SQL/MM Spatial exception - point not on edge");
1001 return NULL;
1002 }
1003
1004 #if 0
1005 {
1006 size_t sz;
1007 char *wkt = lwgeom_to_wkt((LWGEOM*)split_col, WKT_EXTENDED, 2, &sz);
1008 LWDEBUGF(1, "returning split col: %s", wkt);
1009 lwfree(wkt);
1010 }
1011 #endif
1012 return split_col;
1013 }
1014
1015 LWT_ELEMID
lwt_ModEdgeSplit(LWT_TOPOLOGY * topo,LWT_ELEMID edge,LWPOINT * pt,int skipISOChecks)1016 lwt_ModEdgeSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge,
1017 LWPOINT* pt, int skipISOChecks )
1018 {
1019 LWT_ISO_NODE node;
1020 LWT_ISO_EDGE* oldedge = NULL;
1021 LWCOLLECTION *split_col;
1022 const LWGEOM *oldedge_geom;
1023 const LWGEOM *newedge_geom;
1024 LWT_ISO_EDGE newedge1;
1025 LWT_ISO_EDGE seledge, updedge, excedge;
1026 int ret;
1027
1028 split_col = _lwt_EdgeSplit( topo, edge, pt, skipISOChecks, &oldedge );
1029 if ( ! split_col ) return -1; /* should have raised an exception */
1030 oldedge_geom = split_col->geoms[0];
1031 newedge_geom = split_col->geoms[1];
1032 /* Make sure the SRID is set on the subgeom */
1033 ((LWGEOM*)oldedge_geom)->srid = split_col->srid;
1034 ((LWGEOM*)newedge_geom)->srid = split_col->srid;
1035
1036 /* Add new node, getting new id back */
1037 node.node_id = -1;
1038 node.containing_face = -1; /* means not-isolated */
1039 node.geom = pt;
1040 if ( ! lwt_be_insertNodes(topo, &node, 1) )
1041 {
1042 _lwt_release_edges(oldedge, 1);
1043 lwcollection_free(split_col);
1044 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1045 return -1;
1046 }
1047 if (node.node_id == -1) {
1048 /* should have been set by backend */
1049 _lwt_release_edges(oldedge, 1);
1050 lwcollection_free(split_col);
1051 lwerror("Backend coding error: "
1052 "insertNodes callback did not return node_id");
1053 return -1;
1054 }
1055
1056 /* Insert the new edge */
1057 newedge1.edge_id = lwt_be_getNextEdgeId(topo);
1058 if ( newedge1.edge_id == -1 ) {
1059 _lwt_release_edges(oldedge, 1);
1060 lwcollection_free(split_col);
1061 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1062 return -1;
1063 }
1064 newedge1.start_node = node.node_id;
1065 newedge1.end_node = oldedge->end_node;
1066 newedge1.face_left = oldedge->face_left;
1067 newedge1.face_right = oldedge->face_right;
1068 newedge1.next_left = oldedge->next_left == -oldedge->edge_id ?
1069 -newedge1.edge_id : oldedge->next_left;
1070 newedge1.next_right = -oldedge->edge_id;
1071 newedge1.geom = lwgeom_as_lwline(newedge_geom);
1072 /* lwgeom_split of a line should only return lines ... */
1073 if ( ! newedge1.geom ) {
1074 _lwt_release_edges(oldedge, 1);
1075 lwcollection_free(split_col);
1076 lwerror("first geometry in lwgeom_split output is not a line");
1077 return -1;
1078 }
1079 ret = lwt_be_insertEdges(topo, &newedge1, 1);
1080 if ( ret == -1 ) {
1081 _lwt_release_edges(oldedge, 1);
1082 lwcollection_free(split_col);
1083 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1084 return -1;
1085 } else if ( ret == 0 ) {
1086 _lwt_release_edges(oldedge, 1);
1087 lwcollection_free(split_col);
1088 lwerror("Insertion of split edge failed (no reason)");
1089 return -1;
1090 }
1091
1092 /* Update the old edge */
1093 updedge.geom = lwgeom_as_lwline(oldedge_geom);
1094 /* lwgeom_split of a line should only return lines ... */
1095 if ( ! updedge.geom ) {
1096 _lwt_release_edges(oldedge, 1);
1097 lwcollection_free(split_col);
1098 lwerror("second geometry in lwgeom_split output is not a line");
1099 return -1;
1100 }
1101 updedge.next_left = newedge1.edge_id;
1102 updedge.end_node = node.node_id;
1103 ret = lwt_be_updateEdges(topo,
1104 oldedge, LWT_COL_EDGE_EDGE_ID,
1105 &updedge, LWT_COL_EDGE_GEOM|LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE,
1106 NULL, 0);
1107 if ( ret == -1 ) {
1108 _lwt_release_edges(oldedge, 1);
1109 lwcollection_free(split_col);
1110 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1111 return -1;
1112 } else if ( ret == 0 ) {
1113 _lwt_release_edges(oldedge, 1);
1114 lwcollection_free(split_col);
1115 lwerror("Edge being split (%d) disappeared during operations?", oldedge->edge_id);
1116 return -1;
1117 } else if ( ret > 1 ) {
1118 _lwt_release_edges(oldedge, 1);
1119 lwcollection_free(split_col);
1120 lwerror("More than a single edge found with id %d !", oldedge->edge_id);
1121 return -1;
1122 }
1123
1124 /* Update all next edge references to match new layout (ST_ModEdgeSplit) */
1125
1126 updedge.next_right = -newedge1.edge_id;
1127 excedge.edge_id = newedge1.edge_id;
1128 seledge.next_right = -oldedge->edge_id;
1129 seledge.start_node = oldedge->end_node;
1130 ret = lwt_be_updateEdges(topo,
1131 &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE,
1132 &updedge, LWT_COL_EDGE_NEXT_RIGHT,
1133 &excedge, LWT_COL_EDGE_EDGE_ID);
1134 if ( ret == -1 ) {
1135 _lwt_release_edges(oldedge, 1);
1136 lwcollection_free(split_col);
1137 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1138 return -1;
1139 }
1140
1141 updedge.next_left = -newedge1.edge_id;
1142 excedge.edge_id = newedge1.edge_id;
1143 seledge.next_left = -oldedge->edge_id;
1144 seledge.end_node = oldedge->end_node;
1145 ret = lwt_be_updateEdges(topo,
1146 &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE,
1147 &updedge, LWT_COL_EDGE_NEXT_LEFT,
1148 &excedge, LWT_COL_EDGE_EDGE_ID);
1149 if ( ret == -1 ) {
1150 _lwt_release_edges(oldedge, 1);
1151 lwcollection_free(split_col);
1152 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1153 return -1;
1154 }
1155
1156 /* Update TopoGeometries composition */
1157 ret = lwt_be_updateTopoGeomEdgeSplit(topo, oldedge->edge_id, newedge1.edge_id, -1);
1158 if ( ! ret ) {
1159 _lwt_release_edges(oldedge, 1);
1160 lwcollection_free(split_col);
1161 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1162 return -1;
1163 }
1164
1165 _lwt_release_edges(oldedge, 1);
1166 lwcollection_free(split_col);
1167
1168 /* return new node id */
1169 return node.node_id;
1170 }
1171
1172 LWT_ELEMID
lwt_NewEdgesSplit(LWT_TOPOLOGY * topo,LWT_ELEMID edge,LWPOINT * pt,int skipISOChecks)1173 lwt_NewEdgesSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge,
1174 LWPOINT* pt, int skipISOChecks )
1175 {
1176 LWT_ISO_NODE node;
1177 LWT_ISO_EDGE* oldedge = NULL;
1178 LWCOLLECTION *split_col;
1179 const LWGEOM *oldedge_geom;
1180 const LWGEOM *newedge_geom;
1181 LWT_ISO_EDGE newedges[2];
1182 LWT_ISO_EDGE seledge, updedge;
1183 int ret;
1184
1185 split_col = _lwt_EdgeSplit( topo, edge, pt, skipISOChecks, &oldedge );
1186 if ( ! split_col ) return -1; /* should have raised an exception */
1187 oldedge_geom = split_col->geoms[0];
1188 newedge_geom = split_col->geoms[1];
1189 /* Make sure the SRID is set on the subgeom */
1190 ((LWGEOM*)oldedge_geom)->srid = split_col->srid;
1191 ((LWGEOM*)newedge_geom)->srid = split_col->srid;
1192
1193 /* Add new node, getting new id back */
1194 node.node_id = -1;
1195 node.containing_face = -1; /* means not-isolated */
1196 node.geom = pt;
1197 if ( ! lwt_be_insertNodes(topo, &node, 1) )
1198 {
1199 _lwt_release_edges(oldedge, 1);
1200 lwcollection_free(split_col);
1201 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1202 return -1;
1203 }
1204 if (node.node_id == -1) {
1205 _lwt_release_edges(oldedge, 1);
1206 lwcollection_free(split_col);
1207 /* should have been set by backend */
1208 lwerror("Backend coding error: "
1209 "insertNodes callback did not return node_id");
1210 return -1;
1211 }
1212
1213 /* Delete the old edge */
1214 seledge.edge_id = edge;
1215 ret = lwt_be_deleteEdges(topo, &seledge, LWT_COL_EDGE_EDGE_ID);
1216 if ( ret == -1 ) {
1217 _lwt_release_edges(oldedge, 1);
1218 lwcollection_free(split_col);
1219 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1220 return -1;
1221 }
1222
1223 /* Get new edges identifiers */
1224 newedges[0].edge_id = lwt_be_getNextEdgeId(topo);
1225 if ( newedges[0].edge_id == -1 ) {
1226 _lwt_release_edges(oldedge, 1);
1227 lwcollection_free(split_col);
1228 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1229 return -1;
1230 }
1231 newedges[1].edge_id = lwt_be_getNextEdgeId(topo);
1232 if ( newedges[1].edge_id == -1 ) {
1233 _lwt_release_edges(oldedge, 1);
1234 lwcollection_free(split_col);
1235 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1236 return -1;
1237 }
1238
1239 /* Define the first new edge (to new node) */
1240 newedges[0].start_node = oldedge->start_node;
1241 newedges[0].end_node = node.node_id;
1242 newedges[0].face_left = oldedge->face_left;
1243 newedges[0].face_right = oldedge->face_right;
1244 newedges[0].next_left = newedges[1].edge_id;
1245 if ( oldedge->next_right == edge )
1246 newedges[0].next_right = newedges[0].edge_id;
1247 else if ( oldedge->next_right == -edge )
1248 newedges[0].next_right = -newedges[1].edge_id;
1249 else
1250 newedges[0].next_right = oldedge->next_right;
1251 newedges[0].geom = lwgeom_as_lwline(oldedge_geom);
1252 /* lwgeom_split of a line should only return lines ... */
1253 if ( ! newedges[0].geom ) {
1254 _lwt_release_edges(oldedge, 1);
1255 lwcollection_free(split_col);
1256 lwerror("first geometry in lwgeom_split output is not a line");
1257 return -1;
1258 }
1259
1260 /* Define the second new edge (from new node) */
1261 newedges[1].start_node = node.node_id;
1262 newedges[1].end_node = oldedge->end_node;
1263 newedges[1].face_left = oldedge->face_left;
1264 newedges[1].face_right = oldedge->face_right;
1265 newedges[1].next_right = -newedges[0].edge_id;
1266 if ( oldedge->next_left == -edge )
1267 newedges[1].next_left = -newedges[1].edge_id;
1268 else if ( oldedge->next_left == edge )
1269 newedges[1].next_left = newedges[0].edge_id;
1270 else
1271 newedges[1].next_left = oldedge->next_left;
1272 newedges[1].geom = lwgeom_as_lwline(newedge_geom);
1273 /* lwgeom_split of a line should only return lines ... */
1274 if ( ! newedges[1].geom ) {
1275 _lwt_release_edges(oldedge, 1);
1276 lwcollection_free(split_col);
1277 lwerror("second geometry in lwgeom_split output is not a line");
1278 return -1;
1279 }
1280
1281 /* Insert both new edges */
1282 ret = lwt_be_insertEdges(topo, newedges, 2);
1283 if ( ret == -1 ) {
1284 _lwt_release_edges(oldedge, 1);
1285 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1286 return -1;
1287 } else if ( ret == 0 ) {
1288 _lwt_release_edges(oldedge, 1);
1289 lwcollection_free(split_col);
1290 lwerror("Insertion of split edge failed (no reason)");
1291 return -1;
1292 }
1293
1294 /* Update all next edge references pointing to old edge id */
1295
1296 updedge.next_right = newedges[1].edge_id;
1297 seledge.next_right = edge;
1298 seledge.start_node = oldedge->start_node;
1299 ret = lwt_be_updateEdges(topo,
1300 &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE,
1301 &updedge, LWT_COL_EDGE_NEXT_RIGHT,
1302 NULL, 0);
1303 if ( ret == -1 ) {
1304 _lwt_release_edges(oldedge, 1);
1305 lwcollection_free(split_col);
1306 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1307 return -1;
1308 }
1309
1310 updedge.next_right = -newedges[0].edge_id;
1311 seledge.next_right = -edge;
1312 seledge.start_node = oldedge->end_node;
1313 ret = lwt_be_updateEdges(topo,
1314 &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE,
1315 &updedge, LWT_COL_EDGE_NEXT_RIGHT,
1316 NULL, 0);
1317 if ( ret == -1 ) {
1318 _lwt_release_edges(oldedge, 1);
1319 lwcollection_free(split_col);
1320 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1321 return -1;
1322 }
1323
1324 updedge.next_left = newedges[0].edge_id;
1325 seledge.next_left = edge;
1326 seledge.end_node = oldedge->start_node;
1327 ret = lwt_be_updateEdges(topo,
1328 &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE,
1329 &updedge, LWT_COL_EDGE_NEXT_LEFT,
1330 NULL, 0);
1331 if ( ret == -1 ) {
1332 _lwt_release_edges(oldedge, 1);
1333 lwcollection_free(split_col);
1334 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1335 return -1;
1336 }
1337
1338 updedge.next_left = -newedges[1].edge_id;
1339 seledge.next_left = -edge;
1340 seledge.end_node = oldedge->end_node;
1341 ret = lwt_be_updateEdges(topo,
1342 &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE,
1343 &updedge, LWT_COL_EDGE_NEXT_LEFT,
1344 NULL, 0);
1345 if ( ret == -1 ) {
1346 _lwt_release_edges(oldedge, 1);
1347 lwcollection_release(split_col);
1348 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1349 return -1;
1350 }
1351
1352 /* Update TopoGeometries composition */
1353 ret = lwt_be_updateTopoGeomEdgeSplit(topo, oldedge->edge_id, newedges[0].edge_id, newedges[1].edge_id);
1354 if ( ! ret ) {
1355 _lwt_release_edges(oldedge, 1);
1356 lwcollection_free(split_col);
1357 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1358 return -1;
1359 }
1360
1361 _lwt_release_edges(oldedge, 1);
1362 lwcollection_free(split_col);
1363
1364 /* return new node id */
1365 return node.node_id;
1366 }
1367
1368 /* Data structure used by AddEdgeX functions */
1369 typedef struct edgeend_t {
1370 /* Signed identifier of next clockwise edge (+outgoing,-incoming) */
1371 LWT_ELEMID nextCW;
1372 /* Identifier of face between myaz and next CW edge */
1373 LWT_ELEMID cwFace;
1374 /* Signed identifier of next counterclockwise edge (+outgoing,-incoming) */
1375 LWT_ELEMID nextCCW;
1376 /* Identifier of face between myaz and next CCW edge */
1377 LWT_ELEMID ccwFace;
1378 int was_isolated;
1379 double myaz; /* azimuth of edgeend geometry */
1380 } edgeend;
1381
1382 /*
1383 * Get first distinct vertex from endpoint
1384 * @param pa the pointarray to seek points in
1385 * @param ref the point we want to search a distinct one
1386 * @param from vertex index to start from (will really start from "from"+dir)
1387 * @param dir 1 to go forward
1388 * -1 to go backward
1389 * @return 0 if edge is collapsed (no distinct points)
1390 */
1391 static int
_lwt_FirstDistinctVertex2D(const POINTARRAY * pa,POINT2D * ref,int from,int dir,POINT2D * op)1392 _lwt_FirstDistinctVertex2D(const POINTARRAY* pa, POINT2D *ref, int from, int dir, POINT2D *op)
1393 {
1394 int i, toofar, inc;
1395 POINT2D fp;
1396
1397 if ( dir > 0 )
1398 {
1399 toofar = pa->npoints;
1400 inc = 1;
1401 }
1402 else
1403 {
1404 toofar = -1;
1405 inc = -1;
1406 }
1407
1408 LWDEBUGF(1, "first point is index %d", from);
1409 fp = *ref; /* getPoint2d_p(pa, from, &fp); */
1410 for ( i = from+inc; i != toofar; i += inc )
1411 {
1412 LWDEBUGF(1, "testing point %d", i);
1413 getPoint2d_p(pa, i, op); /* pick next point */
1414 if ( p2d_same(op, &fp) ) continue; /* equal to startpoint */
1415 /* this is a good one, neither same of start nor of end point */
1416 return 1; /* found */
1417 }
1418
1419 /* no distinct vertices found */
1420 return 0;
1421 }
1422
1423
1424 /*
1425 * Return non-zero on failure (lwerror is invoked in that case)
1426 * Possible failures:
1427 * -1 no two distinct vertices exist
1428 * -2 azimuth computation failed for first edge end
1429 */
1430 static int
_lwt_InitEdgeEndByLine(edgeend * fee,edgeend * lee,LWLINE * edge,POINT2D * fp,POINT2D * lp)1431 _lwt_InitEdgeEndByLine(edgeend *fee, edgeend *lee, LWLINE *edge,
1432 POINT2D *fp, POINT2D *lp)
1433 {
1434 POINTARRAY *pa = edge->points;
1435 POINT2D pt;
1436
1437 fee->nextCW = fee->nextCCW =
1438 lee->nextCW = lee->nextCCW = 0;
1439 fee->cwFace = fee->ccwFace =
1440 lee->cwFace = lee->ccwFace = -1;
1441
1442 /* Compute azimuth of first edge end */
1443 LWDEBUG(1, "computing azimuth of first edge end");
1444 if ( ! _lwt_FirstDistinctVertex2D(pa, fp, 0, 1, &pt) )
1445 {
1446 lwerror("Invalid edge (no two distinct vertices exist)");
1447 return -1;
1448 }
1449 if ( ! azimuth_pt_pt(fp, &pt, &(fee->myaz)) ) {
1450 lwerror("error computing azimuth of first edgeend [%.15g %.15g,%.15g %.15g]",
1451 fp->x, fp->y, pt.x, pt.y);
1452 return -2;
1453 }
1454 LWDEBUGF(1, "azimuth of first edge end [%.15g %.15g,%.15g %.15g] is %g",
1455 fp->x, fp->y, pt.x, pt.y, fee->myaz);
1456
1457 /* Compute azimuth of second edge end */
1458 LWDEBUG(1, "computing azimuth of second edge end");
1459 if ( ! _lwt_FirstDistinctVertex2D(pa, lp, pa->npoints-1, -1, &pt) )
1460 {
1461 lwerror("Invalid edge (no two distinct vertices exist)");
1462 return -1;
1463 }
1464 if ( ! azimuth_pt_pt(lp, &pt, &(lee->myaz)) ) {
1465 lwerror("error computing azimuth of last edgeend [%.15g %.15g,%.15g %.15g]",
1466 lp->x, lp->y, pt.x, pt.y);
1467 return -2;
1468 }
1469 LWDEBUGF(1, "azimuth of last edge end [%.15g %.15g,%.15g %.15g] is %g",
1470 lp->x, lp->y, pt.x, pt.y, lee->myaz);
1471
1472 return 0;
1473 }
1474
1475 /*
1476 * Find the first edges encountered going clockwise and counterclockwise
1477 * around a node, starting from the given azimuth, and take
1478 * note of the face on the both sides.
1479 *
1480 * @param topo the topology to act upon
1481 * @param node the identifier of the node to analyze
1482 * @param data input (myaz) / output (nextCW, nextCCW) parameter
1483 * @param other edgeend, if also incident to given node (closed edge).
1484 * @param myedge_id identifier of the edge id that data->myaz belongs to
1485 * @return number of incident edges found
1486 *
1487 */
1488 static int
_lwt_FindAdjacentEdges(LWT_TOPOLOGY * topo,LWT_ELEMID node,edgeend * data,edgeend * other,int myedge_id)1489 _lwt_FindAdjacentEdges( LWT_TOPOLOGY* topo, LWT_ELEMID node, edgeend *data,
1490 edgeend *other, int myedge_id )
1491 {
1492 LWT_ISO_EDGE *edges;
1493 uint64_t numedges = 1;
1494 uint64_t i;
1495 double minaz, maxaz;
1496 double az, azdif;
1497
1498 data->nextCW = data->nextCCW = 0;
1499 data->cwFace = data->ccwFace = -1;
1500
1501 if ( other ) {
1502 azdif = other->myaz - data->myaz;
1503 if ( azdif < 0 ) azdif += 2 * M_PI;
1504 minaz = maxaz = azdif;
1505 /* TODO: set nextCW/nextCCW/cwFace/ccwFace to other->something ? */
1506 LWDEBUGF(1, "Other edge end has cwFace=%d and ccwFace=%d",
1507 other->cwFace, other->ccwFace);
1508 } else {
1509 minaz = maxaz = -1;
1510 }
1511
1512 LWDEBUGF(1, "Looking for edges incident to node %" LWTFMT_ELEMID
1513 " and adjacent to azimuth %g", node, data->myaz);
1514
1515 /* Get incident edges */
1516 edges = lwt_be_getEdgeByNode( topo, &node, &numedges, LWT_COL_EDGE_ALL );
1517 if (numedges == UINT64_MAX)
1518 {
1519 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1520 return 0;
1521 }
1522
1523 LWDEBUGF(1, "getEdgeByNode returned %d edges, minaz=%g, maxaz=%g",
1524 numedges, minaz, maxaz);
1525
1526 /* For each incident edge-end (1 or 2): */
1527 for ( i = 0; i < numedges; ++i )
1528 {
1529 LWT_ISO_EDGE *edge;
1530 LWGEOM *g;
1531 LWGEOM *cleangeom;
1532 POINT2D p1, p2;
1533 POINTARRAY *pa;
1534
1535 edge = &(edges[i]);
1536
1537 if ( edge->edge_id == myedge_id ) continue;
1538
1539 g = lwline_as_lwgeom(edge->geom);
1540 /* NOTE: remove_repeated_points call could be replaced by
1541 * some other mean to pick two distinct points for endpoints */
1542 cleangeom = lwgeom_remove_repeated_points( g, 0 );
1543 pa = lwgeom_as_lwline(cleangeom)->points;
1544
1545 if ( pa->npoints < 2 ) {{
1546 LWT_ELEMID id = edge->edge_id;
1547 _lwt_release_edges(edges, numedges);
1548 lwgeom_free(cleangeom);
1549 lwerror("corrupted topology: edge %" LWTFMT_ELEMID
1550 " does not have two distinct points", id);
1551 return -1;
1552 }}
1553
1554 if ( edge->start_node == node ) {
1555 getPoint2d_p(pa, 0, &p1);
1556 if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, 0, 1, &p2) )
1557 {
1558 lwerror("Edge %d has no distinct vertices: [%.15g %.15g,%.15g %.15g]: ",
1559 edge->edge_id, p1.x, p1.y, p2.x, p2.y);
1560 return -1;
1561 }
1562 LWDEBUGF(1, "edge %" LWTFMT_ELEMID
1563 " starts on node %" LWTFMT_ELEMID
1564 ", edgeend is %g,%g-%g,%g",
1565 edge->edge_id, node, p1.x, p1.y, p2.x, p2.y);
1566 if ( ! azimuth_pt_pt(&p1, &p2, &az) ) {{
1567 LWT_ELEMID id = edge->edge_id;
1568 _lwt_release_edges(edges, numedges);
1569 lwgeom_free(cleangeom);
1570 lwerror("error computing azimuth of edge %d first edgeend [%.15g %.15g,%.15g %.15g]",
1571 id, p1.x, p1.y, p2.x, p2.y);
1572 return -1;
1573 }}
1574 azdif = az - data->myaz;
1575 LWDEBUGF(1, "azimuth of edge %" LWTFMT_ELEMID
1576 ": %g (diff: %g)", edge->edge_id, az, azdif);
1577
1578 if ( azdif < 0 ) azdif += 2 * M_PI;
1579 if ( minaz == -1 ) {
1580 minaz = maxaz = azdif;
1581 data->nextCW = data->nextCCW = edge->edge_id; /* outgoing */
1582 data->cwFace = edge->face_left;
1583 data->ccwFace = edge->face_right;
1584 LWDEBUGF(1, "new nextCW and nextCCW edge is %" LWTFMT_ELEMID
1585 ", outgoing, "
1586 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1587 " (face_right is new ccwFace, face_left is new cwFace)",
1588 edge->edge_id, edge->face_left,
1589 edge->face_right);
1590 } else {
1591 if ( azdif < minaz ) {
1592 data->nextCW = edge->edge_id; /* outgoing */
1593 data->cwFace = edge->face_left;
1594 LWDEBUGF(1, "new nextCW edge is %" LWTFMT_ELEMID
1595 ", outgoing, "
1596 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1597 " (previous had minaz=%g, face_left is new cwFace)",
1598 edge->edge_id, edge->face_left,
1599 edge->face_right, minaz);
1600 minaz = azdif;
1601 }
1602 else if ( azdif > maxaz ) {
1603 data->nextCCW = edge->edge_id; /* outgoing */
1604 data->ccwFace = edge->face_right;
1605 LWDEBUGF(1, "new nextCCW edge is %" LWTFMT_ELEMID
1606 ", outgoing, "
1607 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1608 " (previous had maxaz=%g, face_right is new ccwFace)",
1609 edge->edge_id, edge->face_left,
1610 edge->face_right, maxaz);
1611 maxaz = azdif;
1612 }
1613 }
1614 }
1615
1616 if ( edge->end_node == node ) {
1617 getPoint2d_p(pa, pa->npoints-1, &p1);
1618 if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, pa->npoints-1, -1, &p2) )
1619 {
1620 lwerror("Edge %d has no distinct vertices: [%.15g %.15g,%.15g %.15g]: ",
1621 edge->edge_id, p1.x, p1.y, p2.x, p2.y);
1622 return -1;
1623 }
1624 LWDEBUGF(1, "edge %" LWTFMT_ELEMID " ends on node %" LWTFMT_ELEMID
1625 ", edgeend is %g,%g-%g,%g",
1626 edge->edge_id, node, p1.x, p1.y, p2.x, p2.y);
1627 if ( ! azimuth_pt_pt(&p1, &p2, &az) ) {{
1628 LWT_ELEMID id = edge->edge_id;
1629 _lwt_release_edges(edges, numedges);
1630 lwgeom_free(cleangeom);
1631 lwerror("error computing azimuth of edge %d last edgeend [%.15g %.15g,%.15g %.15g]",
1632 id, p1.x, p1.y, p2.x, p2.y);
1633 return -1;
1634 }}
1635 azdif = az - data->myaz;
1636 LWDEBUGF(1, "azimuth of edge %" LWTFMT_ELEMID
1637 ": %g (diff: %g)", edge->edge_id, az, azdif);
1638 if ( azdif < 0 ) azdif += 2 * M_PI;
1639 if ( minaz == -1 ) {
1640 minaz = maxaz = azdif;
1641 data->nextCW = data->nextCCW = -edge->edge_id; /* incoming */
1642 data->cwFace = edge->face_right;
1643 data->ccwFace = edge->face_left;
1644 LWDEBUGF(1, "new nextCW and nextCCW edge is %" LWTFMT_ELEMID
1645 ", incoming, "
1646 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1647 " (face_right is new cwFace, face_left is new ccwFace)",
1648 edge->edge_id, edge->face_left,
1649 edge->face_right);
1650 } else {
1651 if ( azdif < minaz ) {
1652 data->nextCW = -edge->edge_id; /* incoming */
1653 data->cwFace = edge->face_right;
1654 LWDEBUGF(1, "new nextCW edge is %" LWTFMT_ELEMID
1655 ", incoming, "
1656 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1657 " (previous had minaz=%g, face_right is new cwFace)",
1658 edge->edge_id, edge->face_left,
1659 edge->face_right, minaz);
1660 minaz = azdif;
1661 }
1662 else if ( azdif > maxaz ) {
1663 data->nextCCW = -edge->edge_id; /* incoming */
1664 data->ccwFace = edge->face_left;
1665 LWDEBUGF(1, "new nextCCW edge is %" LWTFMT_ELEMID
1666 ", outgoing, from start point, "
1667 "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID
1668 " (previous had maxaz=%g, face_left is new ccwFace)",
1669 edge->edge_id, edge->face_left,
1670 edge->face_right, maxaz);
1671 maxaz = azdif;
1672 }
1673 }
1674 }
1675
1676 lwgeom_free(cleangeom);
1677 }
1678 if ( numedges ) _lwt_release_edges(edges, numedges);
1679
1680 LWDEBUGF(1, "edges adjacent to azimuth %g"
1681 " (incident to node %" LWTFMT_ELEMID ")"
1682 ": CW:%" LWTFMT_ELEMID "(%g) CCW:%" LWTFMT_ELEMID "(%g)",
1683 data->myaz, node, data->nextCW, minaz,
1684 data->nextCCW, maxaz);
1685
1686 if ( myedge_id < 1 && numedges && data->cwFace != data->ccwFace )
1687 {
1688 if ( data->cwFace != -1 && data->ccwFace != -1 ) {
1689 lwerror("Corrupted topology: adjacent edges %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID
1690 " bind different face (%" LWTFMT_ELEMID " and %" LWTFMT_ELEMID ")",
1691 data->nextCW, data->nextCCW,
1692 data->cwFace, data->ccwFace);
1693 return -1;
1694 }
1695 }
1696
1697 /* Return number of incident edges found */
1698 return numedges;
1699 }
1700
1701 /*
1702 * Get a point internal to the line and write it into the "ip"
1703 * parameter
1704 *
1705 * return 0 on failure (line is empty or collapsed), 1 otherwise
1706 */
1707 static int
_lwt_GetInteriorEdgePoint(const LWLINE * edge,POINT2D * ip)1708 _lwt_GetInteriorEdgePoint(const LWLINE* edge, POINT2D* ip)
1709 {
1710 uint32_t i;
1711 POINT2D fp, lp, tp;
1712 POINTARRAY *pa = edge->points;
1713
1714 if ( pa->npoints < 2 ) return 0; /* empty or structurally collapsed */
1715
1716 getPoint2d_p(pa, 0, &fp); /* save first point */
1717 getPoint2d_p(pa, pa->npoints-1, &lp); /* save last point */
1718 for (i=1; i<pa->npoints-1; ++i)
1719 {
1720 getPoint2d_p(pa, i, &tp); /* pick next point */
1721 if ( p2d_same(&tp, &fp) ) continue; /* equal to startpoint */
1722 if ( p2d_same(&tp, &lp) ) continue; /* equal to endpoint */
1723 /* this is a good one, neither same of start nor of end point */
1724 *ip = tp;
1725 return 1; /* found */
1726 }
1727
1728 /* no distinct vertex found */
1729
1730 /* interpolate if start point != end point */
1731
1732 if ( p2d_same(&fp, &lp) ) return 0; /* no distinct points in edge */
1733
1734 ip->x = fp.x + ( (lp.x - fp.x) * 0.5 );
1735 ip->y = fp.y + ( (lp.y - fp.y) * 0.5 );
1736
1737 return 1;
1738 }
1739
1740 static LWPOLY *
_lwt_MakeRingShell(LWT_TOPOLOGY * topo,LWT_ELEMID * signed_edge_ids,uint64_t num_signed_edge_ids)1741 _lwt_MakeRingShell(LWT_TOPOLOGY *topo, LWT_ELEMID *signed_edge_ids, uint64_t num_signed_edge_ids)
1742 {
1743 LWT_ELEMID *edge_ids;
1744 uint64_t numedges, i, j;
1745 LWT_ISO_EDGE *ring_edges;
1746
1747 /* Construct a polygon using edges of the ring */
1748 numedges = 0;
1749 edge_ids = lwalloc(sizeof(LWT_ELEMID)*num_signed_edge_ids);
1750 for (i=0; i<num_signed_edge_ids; ++i) {
1751 int absid = llabs(signed_edge_ids[i]);
1752 int found = 0;
1753 /* Do not add the same edge twice */
1754 for (j=0; j<numedges; ++j) {
1755 if ( edge_ids[j] == absid ) {
1756 found = 1;
1757 break;
1758 }
1759 }
1760 if ( ! found ) edge_ids[numedges++] = absid;
1761 }
1762 i = numedges;
1763 ring_edges = lwt_be_getEdgeById(topo, edge_ids, &i,
1764 LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM);
1765 lwfree( edge_ids );
1766 if (i == UINT64_MAX)
1767 {
1768 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1769 return NULL;
1770 }
1771 else if ( i != numedges )
1772 {
1773 lwfree( signed_edge_ids );
1774 _lwt_release_edges(ring_edges, i);
1775 lwerror("Unexpected error: %d edges found when expecting %d", i, numedges);
1776 return NULL;
1777 }
1778
1779 /* Should now build a polygon with those edges, in the order
1780 * given by GetRingEdges.
1781 */
1782 POINTARRAY *pa = NULL;
1783 for ( i=0; i<num_signed_edge_ids; ++i )
1784 {
1785 LWT_ELEMID eid = signed_edge_ids[i];
1786 LWDEBUGF(2, "Edge %d in ring is edge %" LWTFMT_ELEMID, i, eid);
1787 LWT_ISO_EDGE *edge = NULL;
1788 POINTARRAY *epa;
1789 for ( j=0; j<numedges; ++j )
1790 {
1791 if ( ring_edges[j].edge_id == llabs(eid) )
1792 {
1793 edge = &(ring_edges[j]);
1794 break;
1795 }
1796 }
1797 if ( edge == NULL )
1798 {
1799 _lwt_release_edges(ring_edges, numedges);
1800 lwerror("missing edge that was found in ring edges loop");
1801 return NULL;
1802 }
1803
1804 if ( pa == NULL )
1805 {
1806 pa = ptarray_clone_deep(edge->geom->points);
1807 if ( eid < 0 ) ptarray_reverse_in_place(pa);
1808 }
1809 else
1810 {
1811 if ( eid < 0 )
1812 {
1813 epa = ptarray_clone_deep(edge->geom->points);
1814 ptarray_reverse_in_place(epa);
1815 ptarray_append_ptarray(pa, epa, 0);
1816 ptarray_free(epa);
1817 }
1818 else
1819 {
1820 /* avoid a clone here */
1821 ptarray_append_ptarray(pa, edge->geom->points, 0);
1822 }
1823 }
1824 }
1825 _lwt_release_edges(ring_edges, numedges);
1826 POINTARRAY **points = lwalloc(sizeof(POINTARRAY*));
1827 points[0] = pa;
1828
1829 /* NOTE: the ring may very well have collapsed components,
1830 * which would make it topologically invalid
1831 */
1832 LWPOLY* shell = lwpoly_construct(0, 0, 1, points);
1833 return shell;
1834 }
1835
1836 /*
1837 * Add a split face by walking on the edge side.
1838 *
1839 * @param topo the topology to act upon
1840 * @param sedge edge id and walking side and direction
1841 * (forward,left:positive backward,right:negative)
1842 * @param face the face in which the edge identifier is known to be
1843 * @param mbr_only do not create a new face but update MBR of the current
1844 *
1845 * @return:
1846 * -1: if mbr_only was requested
1847 * 0: if the edge does not form a ring
1848 * -1: if it is impossible to create a face on the requested side
1849 * ( new face on the side is the universe )
1850 * -2: error
1851 * >0 : id of newly added face
1852 */
1853 static LWT_ELEMID
_lwt_AddFaceSplit(LWT_TOPOLOGY * topo,LWT_ELEMID sedge,LWT_ELEMID face,int mbr_only)1854 _lwt_AddFaceSplit( LWT_TOPOLOGY* topo,
1855 LWT_ELEMID sedge, LWT_ELEMID face,
1856 int mbr_only )
1857 {
1858 uint64_t numfaceedges, i, j;
1859 int newface_outside;
1860 uint64_t num_signed_edge_ids;
1861 LWT_ELEMID *signed_edge_ids;
1862 LWT_ISO_EDGE *edges;
1863 LWT_ISO_EDGE *forward_edges = NULL;
1864 int forward_edges_count = 0;
1865 LWT_ISO_EDGE *backward_edges = NULL;
1866 int backward_edges_count = 0;
1867
1868 signed_edge_ids = lwt_be_getRingEdges(topo, sedge, &num_signed_edge_ids, 0);
1869 if (!signed_edge_ids)
1870 {
1871 lwerror("Backend error (no ring edges for edge %" LWTFMT_ELEMID "): %s",
1872 sedge,
1873 lwt_be_lastErrorMessage(topo->be_iface));
1874 return -2;
1875 }
1876 LWDEBUGF(1, "getRingEdges returned %d edges", num_signed_edge_ids);
1877
1878 /* You can't get to the other side of an edge forming a ring */
1879 for (i=0; i<num_signed_edge_ids; ++i) {
1880 if ( signed_edge_ids[i] == -sedge ) {
1881 /* No split here */
1882 LWDEBUG(1, "not a ring");
1883 lwfree( signed_edge_ids );
1884 return 0;
1885 }
1886 }
1887
1888 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " split face %" LWTFMT_ELEMID " (mbr_only:%d)",
1889 sedge, face, mbr_only);
1890
1891 /* Construct a polygon using edges of the ring */
1892 LWPOLY *shell = _lwt_MakeRingShell(topo, signed_edge_ids,
1893 num_signed_edge_ids);
1894 if ( ! shell ) {
1895 lwfree( signed_edge_ids );
1896 /* ring_edges should be NULL */
1897 lwerror("Could not create ring shell: %s", lwt_be_lastErrorMessage(topo->be_iface));
1898 return -2;
1899 }
1900 const POINTARRAY *pa = shell->rings[0];
1901 if ( ! ptarray_is_closed(pa) )
1902 {
1903 lwpoly_free(shell);
1904 lwfree( signed_edge_ids );
1905 lwerror("Corrupted topology: ring of edge %" LWTFMT_ELEMID
1906 " is geometrically not-closed", sedge);
1907 return -2;
1908 }
1909
1910 int isccw = ptarray_isccw(pa);
1911 LWDEBUGF(1, "Ring of edge %" LWTFMT_ELEMID " is %sclockwise",
1912 sedge, isccw ? "counter" : "");
1913 const GBOX* shellbox = lwgeom_get_bbox(lwpoly_as_lwgeom(shell));
1914
1915 if ( face == 0 )
1916 {
1917 /* Edge split the universe face */
1918 if ( ! isccw )
1919 {
1920 lwpoly_free(shell);
1921 lwfree( signed_edge_ids );
1922 /* Face on the left side of this ring is the universe face.
1923 * Next call (for the other side) should create the split face
1924 */
1925 LWDEBUG(1, "The left face of this clockwise ring is the universe, "
1926 "won't create a new face there");
1927 return -1;
1928 }
1929 }
1930
1931 if ( mbr_only && face != 0 )
1932 {
1933 if ( isccw )
1934 {{
1935 LWT_ISO_FACE updface;
1936 updface.face_id = face;
1937 updface.mbr = (GBOX *)shellbox; /* const cast, we won't free it, later */
1938 int ret = lwt_be_updateFacesById( topo, &updface, 1 );
1939 if ( ret == -1 )
1940 {
1941 lwfree( signed_edge_ids );
1942 lwpoly_free(shell); /* NOTE: owns shellbox above */
1943 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1944 return -2;
1945 }
1946 if ( ret != 1 )
1947 {
1948 lwfree( signed_edge_ids );
1949 lwpoly_free(shell); /* NOTE: owns shellbox above */
1950 lwerror("Unexpected error: %d faces found when expecting 1", ret);
1951 return -2;
1952 }
1953 }}
1954 lwfree( signed_edge_ids );
1955 lwpoly_free(shell); /* NOTE: owns shellbox above */
1956 return -1; /* mbr only was requested */
1957 }
1958
1959 LWT_ISO_FACE *oldface = NULL;
1960 LWT_ISO_FACE newface;
1961 newface.face_id = -1;
1962 if ( face != 0 && ! isccw)
1963 {{
1964 /* Face created an hole in an outer face */
1965 uint64_t nfaces = 1;
1966 oldface = lwt_be_getFaceById(topo, &face, &nfaces, LWT_COL_FACE_ALL);
1967 if (nfaces == UINT64_MAX)
1968 {
1969 lwfree( signed_edge_ids );
1970 lwpoly_free(shell); /* NOTE: owns shellbox */
1971 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1972 return -2;
1973 }
1974 if ( nfaces != 1 )
1975 {
1976 lwfree( signed_edge_ids );
1977 lwpoly_free(shell); /* NOTE: owns shellbox */
1978 lwerror("Unexpected error: %d faces found when expecting 1", nfaces);
1979 return -2;
1980 }
1981 newface.mbr = oldface->mbr;
1982 }}
1983 else
1984 {
1985 newface.mbr = (GBOX *)shellbox; /* const cast, we won't free it, later */
1986 }
1987
1988 /* Insert the new face */
1989 int ret = lwt_be_insertFaces( topo, &newface, 1 );
1990 if ( ret == -1 )
1991 {
1992 lwfree( signed_edge_ids );
1993 lwpoly_free(shell); /* NOTE: owns shellbox */
1994 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
1995 return -2;
1996 }
1997 if ( ret != 1 )
1998 {
1999 lwfree( signed_edge_ids );
2000 lwpoly_free(shell); /* NOTE: owns shellbox */
2001 lwerror("Unexpected error: %d faces inserted when expecting 1", ret);
2002 return -2;
2003 }
2004 if ( oldface ) {
2005 newface.mbr = NULL; /* it is a reference to oldface mbr... */
2006 _lwt_release_faces(oldface, 1);
2007 }
2008
2009 /* Update side location of new face edges */
2010
2011 /* We want the new face to be on the left, if possible */
2012 if ( face != 0 && ! isccw ) { /* ring is clockwise in a real face */
2013 /* face shrinked, must update all non-contained edges and nodes */
2014 LWDEBUG(1, "New face is on the outside of the ring, updating rings in former shell");
2015 newface_outside = 1;
2016 /* newface is outside */
2017 } else {
2018 LWDEBUG(1, "New face is on the inside of the ring, updating forward edges in new ring");
2019 newface_outside = 0;
2020 /* newface is inside */
2021 }
2022
2023 /* Update edges bounding the old face */
2024 /* (1) fetch all edges where left_face or right_face is = oldface */
2025 int fields = LWT_COL_EDGE_EDGE_ID |
2026 LWT_COL_EDGE_FACE_LEFT |
2027 LWT_COL_EDGE_FACE_RIGHT |
2028 LWT_COL_EDGE_GEOM
2029 ;
2030 numfaceedges = 1;
2031 edges = lwt_be_getEdgeByFace( topo, &face, &numfaceedges, fields, newface.mbr );
2032 if (numfaceedges == UINT64_MAX)
2033 {
2034 lwfree(signed_edge_ids);
2035 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2036 return -2;
2037 }
2038 LWDEBUGF(1, "lwt_be_getEdgeByFace returned %d edges", numfaceedges);
2039
2040 if ( numfaceedges )
2041 {
2042 forward_edges = lwalloc(sizeof(LWT_ISO_EDGE)*numfaceedges);
2043 forward_edges_count = 0;
2044 backward_edges = lwalloc(sizeof(LWT_ISO_EDGE)*numfaceedges);
2045 backward_edges_count = 0;
2046
2047 /* (2) loop over the results and: */
2048 for ( i=0; i<numfaceedges; ++i )
2049 {
2050 LWT_ISO_EDGE *e = &(edges[i]);
2051 int found = 0;
2052 int contains;
2053 POINT2D ep;
2054
2055 /* (2.1) skip edges whose ID is in the list of boundary edges ? */
2056 for ( j=0; j<num_signed_edge_ids; ++j )
2057 {
2058 int seid = signed_edge_ids[j];
2059 if ( seid == e->edge_id )
2060 {
2061 /* IDEA: remove entry from signed_edge_ids ? */
2062 LWDEBUGF(1, "Edge %d is a forward edge of the new ring", e->edge_id);
2063 forward_edges[forward_edges_count].edge_id = e->edge_id;
2064 forward_edges[forward_edges_count++].face_left = newface.face_id;
2065 found++;
2066 if ( found == 2 ) break;
2067 }
2068 else if ( -seid == e->edge_id )
2069 {
2070 /* IDEA: remove entry from signed_edge_ids ? */
2071 LWDEBUGF(1, "Edge %d is a backward edge of the new ring", e->edge_id);
2072 backward_edges[backward_edges_count].edge_id = e->edge_id;
2073 backward_edges[backward_edges_count++].face_right = newface.face_id;
2074 found++;
2075 if ( found == 2 ) break;
2076 }
2077 }
2078 if ( found ) continue;
2079
2080 /* We need to check only a single point
2081 * (to avoid collapsed elements of the shell polygon
2082 * giving false positive).
2083 * The point but must not be an endpoint.
2084 */
2085 if ( ! _lwt_GetInteriorEdgePoint(e->geom, &ep) )
2086 {
2087 lwfree(signed_edge_ids);
2088 lwpoly_free(shell);
2089 lwfree(forward_edges); /* contents owned by edges */
2090 lwfree(backward_edges); /* contents owned by edges */
2091 _lwt_release_edges(edges, numfaceedges);
2092 lwerror("Could not find interior point for edge %d: %s",
2093 e->edge_id, lwgeom_geos_errmsg);
2094 return -2;
2095 }
2096
2097 /* IDEA: check that bounding box shortcut is taken, or use
2098 * shellbox to do it here */
2099 contains = ptarray_contains_point(pa, &ep) == LW_INSIDE;
2100 if ( contains == 2 )
2101 LWDEBUGF(1, "Edge %d %scontained in new ring", e->edge_id,
2102 (contains?"":"not "));
2103
2104 /* (2.2) skip edges (NOT, if newface_outside) contained in shell */
2105 if ( newface_outside )
2106 {
2107 if ( contains )
2108 {
2109 LWDEBUGF(1, "Edge %d contained in an hole of the new face",
2110 e->edge_id);
2111 continue;
2112 }
2113 }
2114 else
2115 {
2116 if ( ! contains )
2117 {
2118 LWDEBUGF(1, "Edge %d not contained in the face shell",
2119 e->edge_id);
2120 continue;
2121 }
2122 }
2123
2124 /* (2.3) push to forward_edges if left_face = oface */
2125 if ( e->face_left == face )
2126 {
2127 LWDEBUGF(1, "Edge %d has new face on the left side", e->edge_id);
2128 forward_edges[forward_edges_count].edge_id = e->edge_id;
2129 forward_edges[forward_edges_count++].face_left = newface.face_id;
2130 }
2131
2132 /* (2.4) push to backward_edges if right_face = oface */
2133 if ( e->face_right == face )
2134 {
2135 LWDEBUGF(1, "Edge %d has new face on the right side", e->edge_id);
2136 backward_edges[backward_edges_count].edge_id = e->edge_id;
2137 backward_edges[backward_edges_count++].face_right = newface.face_id;
2138 }
2139 }
2140
2141 /* Update forward edges */
2142 if ( forward_edges_count )
2143 {
2144 ret = lwt_be_updateEdgesById(topo, forward_edges,
2145 forward_edges_count,
2146 LWT_COL_EDGE_FACE_LEFT);
2147 if ( ret == -1 )
2148 {
2149 lwfree( signed_edge_ids );
2150 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2151 return -2;
2152 }
2153 if ( ret != forward_edges_count )
2154 {
2155 lwfree( signed_edge_ids );
2156 lwerror("Unexpected error: %d edges updated when expecting %d",
2157 ret, forward_edges_count);
2158 return -2;
2159 }
2160 }
2161
2162 /* Update backward edges */
2163 if ( backward_edges_count )
2164 {
2165 ret = lwt_be_updateEdgesById(topo, backward_edges,
2166 backward_edges_count,
2167 LWT_COL_EDGE_FACE_RIGHT);
2168 if ( ret == -1 )
2169 {
2170 lwfree( signed_edge_ids );
2171 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2172 return -2;
2173 }
2174 if ( ret != backward_edges_count )
2175 {
2176 lwfree( signed_edge_ids );
2177 lwerror("Unexpected error: %d edges updated when expecting %d",
2178 ret, backward_edges_count);
2179 return -2;
2180 }
2181 }
2182
2183 lwfree(forward_edges);
2184 lwfree(backward_edges);
2185
2186 }
2187
2188 _lwt_release_edges(edges, numfaceedges);
2189
2190 /* Update isolated nodes which are now in new face */
2191 uint64_t numisonodes = 1;
2192 fields = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM;
2193 LWT_ISO_NODE *nodes = lwt_be_getNodeByFace(topo, &face,
2194 &numisonodes, fields, newface.mbr);
2195 if (numisonodes == UINT64_MAX)
2196 {
2197 lwfree(signed_edge_ids);
2198 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2199 return -2;
2200 }
2201 if ( numisonodes ) {
2202 LWT_ISO_NODE *updated_nodes = lwalloc(sizeof(LWT_ISO_NODE)*numisonodes);
2203 int nodes_to_update = 0;
2204 for (i=0; i<numisonodes; ++i)
2205 {
2206 LWT_ISO_NODE *n = &(nodes[i]);
2207 const POINT2D *pt = getPoint2d_cp(n->geom->point, 0);
2208 int contains = ptarray_contains_point(pa, pt) == LW_INSIDE;
2209 LWDEBUGF(1, "Node %d is %scontained in new ring, newface is %s",
2210 n->node_id, contains ? "" : "not ",
2211 newface_outside ? "outside" : "inside" );
2212 if ( newface_outside )
2213 {
2214 if ( contains )
2215 {
2216 LWDEBUGF(1, "Node %d contained in an hole of the new face",
2217 n->node_id);
2218 continue;
2219 }
2220 }
2221 else
2222 {
2223 if ( ! contains )
2224 {
2225 LWDEBUGF(1, "Node %d not contained in the face shell",
2226 n->node_id);
2227 continue;
2228 }
2229 }
2230 updated_nodes[nodes_to_update].node_id = n->node_id;
2231 updated_nodes[nodes_to_update++].containing_face =
2232 newface.face_id;
2233 LWDEBUGF(1, "Node %d will be updated", n->node_id);
2234 }
2235 _lwt_release_nodes(nodes, numisonodes);
2236 if ( nodes_to_update )
2237 {
2238 int ret = lwt_be_updateNodesById(topo, updated_nodes,
2239 nodes_to_update,
2240 LWT_COL_NODE_CONTAINING_FACE);
2241 if ( ret == -1 ) {
2242 lwfree( signed_edge_ids );
2243 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2244 return -2;
2245 }
2246 }
2247 lwfree(updated_nodes);
2248 }
2249
2250 lwfree(signed_edge_ids);
2251 lwpoly_free(shell);
2252
2253 return newface.face_id;
2254 }
2255
2256 /**
2257 * @param modFace can be
2258 * 0 - have two new faces replace a splitted face
2259 * 1 - modify a splitted face, adding a new one
2260 * -1 - do not check at all for face splitting
2261 *
2262 */
2263 static LWT_ELEMID
_lwt_AddEdge(LWT_TOPOLOGY * topo,LWT_ELEMID start_node,LWT_ELEMID end_node,LWLINE * geom,int skipChecks,int modFace)2264 _lwt_AddEdge( LWT_TOPOLOGY* topo,
2265 LWT_ELEMID start_node, LWT_ELEMID end_node,
2266 LWLINE *geom, int skipChecks, int modFace )
2267 {
2268 LWT_ISO_EDGE newedge;
2269 LWGEOM *cleangeom;
2270 edgeend span; /* start point analisys */
2271 edgeend epan; /* end point analisys */
2272 POINT2D p1, pn, p2;
2273 POINTARRAY *pa;
2274 LWT_ELEMID node_ids[2];
2275 const LWPOINT *start_node_geom = NULL;
2276 const LWPOINT *end_node_geom = NULL;
2277 uint64_t num_nodes;
2278 LWT_ISO_NODE *endpoints;
2279 uint64_t i;
2280 int prev_left;
2281 int prev_right;
2282 LWT_ISO_EDGE seledge;
2283 LWT_ISO_EDGE updedge;
2284
2285 if ( ! skipChecks )
2286 {
2287 /* curve must be simple */
2288 if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) )
2289 {
2290 lwerror("SQL/MM Spatial exception - curve not simple");
2291 return -1;
2292 }
2293 }
2294
2295 newedge.start_node = start_node;
2296 newedge.end_node = end_node;
2297 newedge.geom = geom;
2298 newedge.face_left = -1;
2299 newedge.face_right = -1;
2300 /* TODO: should do the repeated points removal in 2D space */
2301 cleangeom = lwgeom_remove_repeated_points( lwline_as_lwgeom(geom), 0 );
2302
2303 pa = lwgeom_as_lwline(cleangeom)->points;
2304 if ( pa->npoints < 2 ) {
2305 lwgeom_free(cleangeom);
2306 lwerror("Invalid edge (no two distinct vertices exist)");
2307 return -1;
2308 }
2309
2310 /* Initialize endpoint info (some of that ) */
2311 span.cwFace = span.ccwFace =
2312 epan.cwFace = epan.ccwFace = -1;
2313
2314 /* Compute azimuth of first edge end on start node */
2315 getPoint2d_p(pa, 0, &p1);
2316 if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, 0, 1, &pn) )
2317 {
2318 lwgeom_free(cleangeom);
2319 lwerror("Invalid edge (no two distinct vertices exist)");
2320 return -1;
2321 }
2322 if ( ! azimuth_pt_pt(&p1, &pn, &span.myaz) ) {
2323 lwgeom_free(cleangeom);
2324 lwerror("error computing azimuth of first edgeend [%.15g %.15g,%.15g %.15g]",
2325 p1.x, p1.y, pn.x, pn.y);
2326 return -1;
2327 }
2328 LWDEBUGF(1, "edge's start node is %g,%g", p1.x, p1.y);
2329
2330 /* Compute azimuth of last edge end on end node */
2331 getPoint2d_p(pa, pa->npoints-1, &p2);
2332 if ( ! _lwt_FirstDistinctVertex2D(pa, &p2, pa->npoints-1, -1, &pn) )
2333 {
2334 lwgeom_free(cleangeom);
2335 /* This should never happen as we checked the edge while computing first edgend */
2336 lwerror("Invalid clean edge (no two distinct vertices exist) - should not happen");
2337 return -1;
2338 }
2339 lwgeom_free(cleangeom);
2340 if ( ! azimuth_pt_pt(&p2, &pn, &epan.myaz) ) {
2341 lwerror("error computing azimuth of last edgeend [%.15g %.15g,%.15g %.15g]",
2342 p2.x, p2.y, pn.x, pn.y);
2343 return -1;
2344 }
2345 LWDEBUGF(1, "edge's end node is %g,%g", p2.x, p2.y);
2346
2347 /*
2348 * Check endpoints existence, match with Curve geometry
2349 * and get face information (if any)
2350 */
2351
2352 if ( start_node != end_node ) {
2353 num_nodes = 2;
2354 node_ids[0] = start_node;
2355 node_ids[1] = end_node;
2356 } else {
2357 num_nodes = 1;
2358 node_ids[0] = start_node;
2359 }
2360
2361 endpoints = lwt_be_getNodeById( topo, node_ids, &num_nodes, LWT_COL_NODE_ALL );
2362 if (num_nodes == UINT64_MAX)
2363 {
2364 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2365 return -1;
2366 }
2367 for ( i=0; i<num_nodes; ++i )
2368 {
2369 LWT_ISO_NODE* node = &(endpoints[i]);
2370 if ( node->containing_face != -1 )
2371 {
2372 if ( newedge.face_left == -1 )
2373 {
2374 newedge.face_left = newedge.face_right = node->containing_face;
2375 }
2376 else if ( newedge.face_left != node->containing_face )
2377 {
2378 _lwt_release_nodes(endpoints, num_nodes);
2379 lwerror("SQL/MM Spatial exception - geometry crosses an edge"
2380 " (endnodes in faces %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID ")",
2381 newedge.face_left, node->containing_face);
2382 }
2383 }
2384
2385 LWDEBUGF(1, "Node %d, with geom %p (looking for %d and %d)",
2386 node->node_id, node->geom, start_node, end_node);
2387 if ( node->node_id == start_node ) {
2388 start_node_geom = node->geom;
2389 }
2390 if ( node->node_id == end_node ) {
2391 end_node_geom = node->geom;
2392 }
2393 }
2394
2395 if ( ! skipChecks )
2396 {
2397 if ( ! start_node_geom )
2398 {
2399 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
2400 lwerror("SQL/MM Spatial exception - non-existent node");
2401 return -1;
2402 }
2403 else
2404 {
2405 pa = start_node_geom->point;
2406 getPoint2d_p(pa, 0, &pn);
2407 if ( ! p2d_same(&pn, &p1) )
2408 {
2409 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
2410 lwerror("SQL/MM Spatial exception"
2411 " - start node not geometry start point."
2412 //" - start node not geometry start point (%g,%g != %g,%g).", pn.x, pn.y, p1.x, p1.y
2413 );
2414 return -1;
2415 }
2416 }
2417
2418 if ( ! end_node_geom )
2419 {
2420 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
2421 lwerror("SQL/MM Spatial exception - non-existent node");
2422 return -1;
2423 }
2424 else
2425 {
2426 pa = end_node_geom->point;
2427 getPoint2d_p(pa, 0, &pn);
2428 if ( ! p2d_same(&pn, &p2) )
2429 {
2430 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
2431 lwerror("SQL/MM Spatial exception"
2432 " - end node not geometry end point."
2433 //" - end node not geometry end point (%g,%g != %g,%g).", pn.x, pn.y, p2.x, p2.y
2434 );
2435 return -1;
2436 }
2437 }
2438
2439 if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes);
2440
2441 if ( _lwt_CheckEdgeCrossing( topo, start_node, end_node, geom, 0 ) )
2442 return -1;
2443
2444 } /* ! skipChecks */
2445
2446 /*
2447 * All checks passed, time to prepare the new edge
2448 */
2449
2450 newedge.edge_id = lwt_be_getNextEdgeId( topo );
2451 if ( newedge.edge_id == -1 ) {
2452 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2453 return -1;
2454 }
2455
2456 /* Find adjacent edges to each endpoint */
2457 int isclosed = start_node == end_node;
2458 int found;
2459 found = _lwt_FindAdjacentEdges( topo, start_node, &span,
2460 isclosed ? &epan : NULL, -1 );
2461 if ( found ) {
2462 span.was_isolated = 0;
2463 newedge.next_right = span.nextCW ? span.nextCW : -newedge.edge_id;
2464 prev_left = span.nextCCW ? -span.nextCCW : newedge.edge_id;
2465 LWDEBUGF(1, "New edge %d is connected on start node, "
2466 "next_right is %d, prev_left is %d",
2467 newedge.edge_id, newedge.next_right, prev_left);
2468 if ( newedge.face_right == -1 ) {
2469 newedge.face_right = span.cwFace;
2470 }
2471 if ( newedge.face_left == -1 ) {
2472 newedge.face_left = span.ccwFace;
2473 }
2474 } else {
2475 span.was_isolated = 1;
2476 newedge.next_right = isclosed ? -newedge.edge_id : newedge.edge_id;
2477 prev_left = isclosed ? newedge.edge_id : -newedge.edge_id;
2478 LWDEBUGF(1, "New edge %d is isolated on start node, "
2479 "next_right is %d, prev_left is %d",
2480 newedge.edge_id, newedge.next_right, prev_left);
2481 }
2482
2483 found = _lwt_FindAdjacentEdges( topo, end_node, &epan,
2484 isclosed ? &span : NULL, -1 );
2485 if ( found ) {
2486 epan.was_isolated = 0;
2487 newedge.next_left = epan.nextCW ? epan.nextCW : newedge.edge_id;
2488 prev_right = epan.nextCCW ? -epan.nextCCW : -newedge.edge_id;
2489 LWDEBUGF(1, "New edge %d is connected on end node, "
2490 "next_left is %d, prev_right is %d",
2491 newedge.edge_id, newedge.next_left, prev_right);
2492 if ( newedge.face_right == -1 ) {
2493 newedge.face_right = span.ccwFace;
2494 } else if ( modFace != -1 && newedge.face_right != epan.ccwFace ) {
2495 /* side-location conflict */
2496 lwerror("Side-location conflict: "
2497 "new edge starts in face"
2498 " %" LWTFMT_ELEMID " and ends in face"
2499 " %" LWTFMT_ELEMID,
2500 newedge.face_right, epan.ccwFace
2501 );
2502 return -1;
2503 }
2504 if ( newedge.face_left == -1 ) {
2505 newedge.face_left = span.cwFace;
2506 } else if ( modFace != -1 && newedge.face_left != epan.cwFace ) {
2507 /* side-location conflict */
2508 lwerror("Side-location conflict: "
2509 "new edge starts in face"
2510 " %" LWTFMT_ELEMID " and ends in face"
2511 " %" LWTFMT_ELEMID,
2512 newedge.face_left, epan.cwFace
2513 );
2514 return -1;
2515 }
2516 } else {
2517 epan.was_isolated = 1;
2518 newedge.next_left = isclosed ? newedge.edge_id : -newedge.edge_id;
2519 prev_right = isclosed ? -newedge.edge_id : newedge.edge_id;
2520 LWDEBUGF(1, "New edge %d is isolated on end node, "
2521 "next_left is %d, prev_right is %d",
2522 newedge.edge_id, newedge.next_left, prev_right);
2523 }
2524
2525 /*
2526 * If we don't have faces setup by now we must have encountered
2527 * a malformed topology (no containing_face on isolated nodes, no
2528 * left/right faces on adjacent edges or mismatching values)
2529 */
2530 if ( newedge.face_left != newedge.face_right )
2531 {
2532 lwerror("Left(%" LWTFMT_ELEMID ")/right(%" LWTFMT_ELEMID ")"
2533 "faces mismatch: invalid topology ?",
2534 newedge.face_left, newedge.face_right);
2535 return -1;
2536 }
2537 else if ( newedge.face_left == -1 && modFace > -1 )
2538 {
2539 lwerror("Could not derive edge face from linked primitives:"
2540 " invalid topology ?");
2541 return -1;
2542 }
2543
2544 /*
2545 * Insert the new edge, and update all linking
2546 */
2547
2548 int ret = lwt_be_insertEdges(topo, &newedge, 1);
2549 if ( ret == -1 ) {
2550 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2551 return -1;
2552 } else if ( ret == 0 ) {
2553 lwerror("Insertion of split edge failed (no reason)");
2554 return -1;
2555 }
2556
2557 int updfields;
2558
2559 /* Link prev_left to us
2560 * (if it's not us already) */
2561 if ( llabs(prev_left) != newedge.edge_id )
2562 {
2563 if ( prev_left > 0 )
2564 {
2565 /* its next_left_edge is us */
2566 updfields = LWT_COL_EDGE_NEXT_LEFT;
2567 updedge.next_left = newedge.edge_id;
2568 seledge.edge_id = prev_left;
2569 }
2570 else
2571 {
2572 /* its next_right_edge is us */
2573 updfields = LWT_COL_EDGE_NEXT_RIGHT;
2574 updedge.next_right = newedge.edge_id;
2575 seledge.edge_id = -prev_left;
2576 }
2577
2578 ret = lwt_be_updateEdges(topo,
2579 &seledge, LWT_COL_EDGE_EDGE_ID,
2580 &updedge, updfields,
2581 NULL, 0);
2582 if ( ret == -1 ) {
2583 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2584 return -1;
2585 }
2586 }
2587
2588 /* Link prev_right to us
2589 * (if it's not us already) */
2590 if ( llabs(prev_right) != newedge.edge_id )
2591 {
2592 if ( prev_right > 0 )
2593 {
2594 /* its next_left_edge is -us */
2595 updfields = LWT_COL_EDGE_NEXT_LEFT;
2596 updedge.next_left = -newedge.edge_id;
2597 seledge.edge_id = prev_right;
2598 }
2599 else
2600 {
2601 /* its next_right_edge is -us */
2602 updfields = LWT_COL_EDGE_NEXT_RIGHT;
2603 updedge.next_right = -newedge.edge_id;
2604 seledge.edge_id = -prev_right;
2605 }
2606
2607 ret = lwt_be_updateEdges(topo,
2608 &seledge, LWT_COL_EDGE_EDGE_ID,
2609 &updedge, updfields,
2610 NULL, 0);
2611 if ( ret == -1 ) {
2612 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2613 return -1;
2614 }
2615 }
2616
2617 /* NOT IN THE SPECS...
2618 * set containing_face = null for start_node and end_node
2619 * if they where isolated
2620 *
2621 */
2622 LWT_ISO_NODE updnode, selnode;
2623 updnode.containing_face = -1;
2624 if ( span.was_isolated )
2625 {
2626 selnode.node_id = start_node;
2627 ret = lwt_be_updateNodes(topo,
2628 &selnode, LWT_COL_NODE_NODE_ID,
2629 &updnode, LWT_COL_NODE_CONTAINING_FACE,
2630 NULL, 0);
2631 if ( ret == -1 ) {
2632 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2633 return -1;
2634 }
2635 }
2636 if ( epan.was_isolated )
2637 {
2638 selnode.node_id = end_node;
2639 ret = lwt_be_updateNodes(topo,
2640 &selnode, LWT_COL_NODE_NODE_ID,
2641 &updnode, LWT_COL_NODE_CONTAINING_FACE,
2642 NULL, 0);
2643 if ( ret == -1 ) {
2644 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2645 return -1;
2646 }
2647 }
2648
2649 /* Check face splitting, if required */
2650
2651 if ( modFace > -1 ) {
2652
2653 if ( ! isclosed && ( epan.was_isolated || span.was_isolated ) )
2654 {
2655 LWDEBUG(1, "New edge is dangling, so it cannot split any face");
2656 return newedge.edge_id; /* no split */
2657 }
2658
2659 int newface1 = -1;
2660
2661 /* IDEA: avoid building edge ring if input is closed, which means we
2662 * know in advance it splits a face */
2663
2664 if ( ! modFace )
2665 {
2666 newface1 = _lwt_AddFaceSplit( topo, -newedge.edge_id, newedge.face_left, 0 );
2667 if ( newface1 == 0 ) {
2668 LWDEBUG(1, "New edge does not split any face");
2669 return newedge.edge_id; /* no split */
2670 }
2671 }
2672
2673 int newface = _lwt_AddFaceSplit( topo, newedge.edge_id,
2674 newedge.face_left, 0 );
2675 if ( modFace )
2676 {
2677 if ( newface == 0 ) {
2678 LWDEBUG(1, "New edge does not split any face");
2679 return newedge.edge_id; /* no split */
2680 }
2681
2682 if ( newface < 0 )
2683 {
2684 /* face on the left is the universe face */
2685 /* must be forming a maximal ring in universal face */
2686 newface = _lwt_AddFaceSplit( topo, -newedge.edge_id,
2687 newedge.face_left, 0 );
2688 if ( newface < 0 ) return newedge.edge_id; /* no split */
2689 }
2690 else
2691 {
2692 _lwt_AddFaceSplit( topo, -newedge.edge_id, newedge.face_left, 1 );
2693 }
2694 }
2695
2696 /*
2697 * Update topogeometries, if needed
2698 */
2699 if ( newedge.face_left != 0 )
2700 {
2701 ret = lwt_be_updateTopoGeomFaceSplit(topo, newedge.face_left,
2702 newface, newface1);
2703 if ( ret == 0 ) {
2704 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2705 return -1;
2706 }
2707
2708 if ( ! modFace )
2709 {
2710 /* drop old face from the face table */
2711 ret = lwt_be_deleteFacesById(topo, &(newedge.face_left), 1);
2712 if ( ret == -1 ) {
2713 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2714 return -1;
2715 }
2716 }
2717 }
2718
2719 } // end of face split checking
2720
2721 return newedge.edge_id;
2722 }
2723
2724 LWT_ELEMID
lwt_AddEdgeModFace(LWT_TOPOLOGY * topo,LWT_ELEMID start_node,LWT_ELEMID end_node,LWLINE * geom,int skipChecks)2725 lwt_AddEdgeModFace( LWT_TOPOLOGY* topo,
2726 LWT_ELEMID start_node, LWT_ELEMID end_node,
2727 LWLINE *geom, int skipChecks )
2728 {
2729 return _lwt_AddEdge( topo, start_node, end_node, geom, skipChecks, 1 );
2730 }
2731
2732 LWT_ELEMID
lwt_AddEdgeNewFaces(LWT_TOPOLOGY * topo,LWT_ELEMID start_node,LWT_ELEMID end_node,LWLINE * geom,int skipChecks)2733 lwt_AddEdgeNewFaces( LWT_TOPOLOGY* topo,
2734 LWT_ELEMID start_node, LWT_ELEMID end_node,
2735 LWLINE *geom, int skipChecks )
2736 {
2737 return _lwt_AddEdge( topo, start_node, end_node, geom, skipChecks, 0 );
2738 }
2739
2740 static LWGEOM *
_lwt_FaceByEdges(LWT_TOPOLOGY * topo,LWT_ISO_EDGE * edges,int numfaceedges)2741 _lwt_FaceByEdges(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edges, int numfaceedges)
2742 {
2743 LWGEOM *outg;
2744 LWCOLLECTION *bounds;
2745 LWGEOM **geoms = lwalloc( sizeof(LWGEOM*) * numfaceedges );
2746 int i, validedges = 0;
2747
2748 for ( i=0; i<numfaceedges; ++i )
2749 {
2750 /* NOTE: skipping edges with same face on both sides, although
2751 * correct, results in a failure to build faces from
2752 * invalid topologies as expected by legacy tests.
2753 * TODO: update legacy tests expectances/unleash this skipping ?
2754 */
2755 /* if ( edges[i].face_left == edges[i].face_right ) continue; */
2756 geoms[validedges++] = lwline_as_lwgeom(edges[i].geom);
2757 }
2758 if ( ! validedges )
2759 {
2760 /* Face has no valid boundary edges, we'll return EMPTY, see
2761 * https://trac.osgeo.org/postgis/ticket/3221 */
2762 if ( numfaceedges ) lwfree(geoms);
2763 LWDEBUG(1, "_lwt_FaceByEdges returning empty polygon");
2764 return lwpoly_as_lwgeom(
2765 lwpoly_construct_empty(topo->srid, topo->hasZ, 0)
2766 );
2767 }
2768 bounds = lwcollection_construct(MULTILINETYPE,
2769 topo->srid,
2770 NULL, /* gbox */
2771 validedges,
2772 geoms);
2773 outg = lwgeom_buildarea( lwcollection_as_lwgeom(bounds) );
2774 lwcollection_release(bounds);
2775 lwfree(geoms);
2776 #if 0
2777 {
2778 size_t sz;
2779 char *wkt = lwgeom_to_wkt(outg, WKT_EXTENDED, 2, &sz);
2780 LWDEBUGF(1, "_lwt_FaceByEdges returning area: %s", wkt);
2781 lwfree(wkt);
2782 }
2783 #endif
2784 return outg;
2785 }
2786
2787 LWGEOM*
lwt_GetFaceGeometry(LWT_TOPOLOGY * topo,LWT_ELEMID faceid)2788 lwt_GetFaceGeometry(LWT_TOPOLOGY* topo, LWT_ELEMID faceid)
2789 {
2790 uint64_t numfaceedges;
2791 LWT_ISO_EDGE *edges;
2792 LWT_ISO_FACE *face;
2793 LWPOLY *out;
2794 LWGEOM *outg;
2795 uint64_t i, edgeid;
2796 int fields;
2797
2798 if (faceid == 0)
2799 {
2800 lwerror("SQL/MM Spatial exception - universal face has no geometry");
2801 return NULL;
2802 }
2803
2804 /* Construct the face geometry */
2805 numfaceedges = 1;
2806 fields = LWT_COL_EDGE_GEOM |
2807 LWT_COL_EDGE_EDGE_ID |
2808 LWT_COL_EDGE_FACE_LEFT |
2809 LWT_COL_EDGE_FACE_RIGHT
2810 ;
2811 edges = lwt_be_getEdgeByFace( topo, &faceid, &numfaceedges, fields, NULL );
2812 if (numfaceedges == UINT64_MAX)
2813 {
2814 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2815 return NULL;
2816 }
2817
2818 if ( numfaceedges == 0 )
2819 {
2820 i = 1;
2821 face = lwt_be_getFaceById(topo, &faceid, &i, LWT_COL_FACE_FACE_ID);
2822 if (i == UINT64_MAX)
2823 {
2824 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2825 return NULL;
2826 }
2827 if ( i == 0 ) {
2828 lwerror("SQL/MM Spatial exception - non-existent face.");
2829 return NULL;
2830 }
2831 lwfree( face );
2832 if ( i > 1 ) {
2833 lwerror("Corrupted topology: multiple face records have face_id=%"
2834 LWTFMT_ELEMID, faceid);
2835 return NULL;
2836 }
2837 /* Face has no boundary edges, we'll return EMPTY, see
2838 * https://trac.osgeo.org/postgis/ticket/3221 */
2839 lwnotice("Corrupted topology: face %"
2840 LWTFMT_ELEMID " has no associated edges.", faceid);
2841 out = lwpoly_construct_empty(topo->srid, topo->hasZ, 0);
2842 return lwpoly_as_lwgeom(out);
2843 }
2844 edgeid = edges[0].edge_id;
2845
2846 outg = _lwt_FaceByEdges( topo, edges, numfaceedges );
2847 _lwt_release_edges(edges, numfaceedges);
2848
2849 if ( ! outg )
2850 {
2851 /* Face did have edges but no polygon could be constructed
2852 * with that material, sounds like a corrupted topology..
2853 *
2854 * We'll return EMPTY, see
2855 * https://trac.osgeo.org/postgis/ticket/3221 */
2856 lwnotice("Corrupted topology: face %"
2857 LWTFMT_ELEMID " could not be constructed only from edges "
2858 "knowing about it (like edge %" LWTFMT_ELEMID ").",
2859 faceid, edgeid);
2860 out = lwpoly_construct_empty(topo->srid, topo->hasZ, 0);
2861 return lwpoly_as_lwgeom(out);
2862 }
2863
2864 return outg;
2865 }
2866
2867 /* Find which edge from the "edges" set defines the next
2868 * portion of the given "ring".
2869 *
2870 * The edge might be either forward or backward.
2871 *
2872 * @param ring The ring to find definition of.
2873 * It is assumed it does not contain duplicated vertices.
2874 * @param from offset of the ring point to start looking from
2875 * @param edges array of edges to search into
2876 * @param numedges number of edges in the edges array
2877 *
2878 * @return index of the edge defining the next ring portion or
2879 * -1 if no edge was found to be part of the ring
2880 */
2881 static int
_lwt_FindNextRingEdge(const POINTARRAY * ring,int from,const LWT_ISO_EDGE * edges,int numedges)2882 _lwt_FindNextRingEdge(const POINTARRAY *ring, int from,
2883 const LWT_ISO_EDGE *edges, int numedges)
2884 {
2885 int i;
2886 POINT2D p1;
2887
2888 /* Get starting ring point */
2889 getPoint2d_p(ring, from, &p1);
2890
2891 LWDEBUGF(1, "Ring's 'from' point (%d) is %g,%g", from, p1.x, p1.y);
2892
2893 /* find the edges defining the next portion of ring starting from
2894 * vertex "from" */
2895 for ( i=0; i<numedges; ++i )
2896 {
2897 const LWT_ISO_EDGE *isoe = &(edges[i]);
2898 LWLINE *edge = isoe->geom;
2899 POINTARRAY *epa = edge->points;
2900 POINT2D p2, pt;
2901 int match = 0;
2902 uint32_t j;
2903
2904 /* Skip if the edge is a dangling one */
2905 if ( isoe->face_left == isoe->face_right )
2906 {
2907 LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID
2908 " has same face (%" LWTFMT_ELEMID
2909 ") on both sides, skipping",
2910 isoe->edge_id, isoe->face_left);
2911 continue;
2912 }
2913
2914 if (epa->npoints < 2)
2915 {
2916 LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID
2917 " has only %"PRIu32" points",
2918 isoe->edge_id, epa->npoints);
2919 continue;
2920 }
2921
2922 #if 0
2923 size_t sz;
2924 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s",
2925 isoe->edge_id,
2926 lwgeom_to_wkt(lwline_as_lwgeom(edge), WKT_EXTENDED, 2, &sz));
2927 #endif
2928
2929 /* ptarray_remove_repeated_points ? */
2930
2931 getPoint2d_p(epa, 0, &p2);
2932 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'first' point is %g,%g",
2933 isoe->edge_id, p2.x, p2.y);
2934 LWDEBUGF(1, "Rings's 'from' point is still %g,%g", p1.x, p1.y);
2935 if ( p2d_same(&p1, &p2) )
2936 {
2937 LWDEBUG(1, "p2d_same(p1,p2) returned true");
2938 LWDEBUGF(1, "First point of edge %" LWTFMT_ELEMID
2939 " matches ring vertex %d", isoe->edge_id, from);
2940 /* first point matches, let's check next non-equal one */
2941 for ( j=1; j<epa->npoints; ++j )
2942 {
2943 getPoint2d_p(epa, j, &p2);
2944 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'next' point %d is %g,%g",
2945 isoe->edge_id, j, p2.x, p2.y);
2946 /* we won't check duplicated edge points */
2947 if ( p2d_same(&p1, &p2) ) continue;
2948 /* we assume there are no duplicated points in ring */
2949 getPoint2d_p(ring, from+1, &pt);
2950 LWDEBUGF(1, "Ring's point %d is %g,%g",
2951 from+1, pt.x, pt.y);
2952 match = p2d_same(&pt, &p2);
2953 break; /* we want to check a single non-equal next vertex */
2954 }
2955 #if POSTGIS_DEBUG_LEVEL > 0
2956 if ( match ) {
2957 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2958 " matches ring vertex %d", isoe->edge_id, from+1);
2959 } else {
2960 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2961 " does not match ring vertex %d", isoe->edge_id, from+1);
2962 }
2963 #endif
2964 }
2965
2966 if ( ! match )
2967 {
2968 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " did not match as forward",
2969 isoe->edge_id);
2970 getPoint2d_p(epa, epa->npoints-1, &p2);
2971 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'last' point is %g,%g",
2972 isoe->edge_id, p2.x, p2.y);
2973 if ( p2d_same(&p1, &p2) )
2974 {
2975 LWDEBUGF(1, "Last point of edge %" LWTFMT_ELEMID
2976 " matches ring vertex %d", isoe->edge_id, from);
2977 /* last point matches, let's check next non-equal one */
2978 for ( j=2; j<=epa->npoints; j++ )
2979 {
2980 getPoint2d_p(epa, epa->npoints - j, &p2);
2981 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'prev' point %d is %g,%g",
2982 isoe->edge_id, epa->npoints - j, p2.x, p2.y);
2983 /* we won't check duplicated edge points */
2984 if ( p2d_same(&p1, &p2) ) continue;
2985 /* we assume there are no duplicated points in ring */
2986 getPoint2d_p(ring, from+1, &pt);
2987 LWDEBUGF(1, "Ring's point %d is %g,%g",
2988 from+1, pt.x, pt.y);
2989 match = p2d_same(&pt, &p2);
2990 break; /* we want to check a single non-equal next vertex */
2991 }
2992 }
2993 #if POSTGIS_DEBUG_LEVEL > 0
2994 if ( match ) {
2995 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2996 " matches ring vertex %d", isoe->edge_id, from+1);
2997 } else {
2998 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2999 " does not match ring vertex %d", isoe->edge_id, from+1);
3000 }
3001 #endif
3002 }
3003
3004 if ( match ) return i;
3005
3006 }
3007
3008 return -1;
3009 }
3010
3011 /* Reverse values in array between "from" (inclusive)
3012 * and "to" (exclusive) indexes */
3013 static void
_lwt_ReverseElemidArray(LWT_ELEMID * ary,int from,int to)3014 _lwt_ReverseElemidArray(LWT_ELEMID *ary, int from, int to)
3015 {
3016 LWT_ELEMID t;
3017 while (from < to)
3018 {
3019 t = ary[from];
3020 ary[from++] = ary[to];
3021 ary[to--] = t;
3022 }
3023 }
3024
3025 /* Rotate values in array between "from" (inclusive)
3026 * and "to" (exclusive) indexes, so that "rotidx" is
3027 * the new value at "from" */
3028 static void
_lwt_RotateElemidArray(LWT_ELEMID * ary,int from,int to,int rotidx)3029 _lwt_RotateElemidArray(LWT_ELEMID *ary, int from, int to, int rotidx)
3030 {
3031 _lwt_ReverseElemidArray(ary, from, rotidx-1);
3032 _lwt_ReverseElemidArray(ary, rotidx, to-1);
3033 _lwt_ReverseElemidArray(ary, from, to-1);
3034 }
3035
3036
3037 int
lwt_GetFaceEdges(LWT_TOPOLOGY * topo,LWT_ELEMID face_id,LWT_ELEMID ** out)3038 lwt_GetFaceEdges(LWT_TOPOLOGY* topo, LWT_ELEMID face_id, LWT_ELEMID **out )
3039 {
3040 LWGEOM *face;
3041 LWPOLY *facepoly;
3042 LWT_ISO_EDGE *edges;
3043 uint64_t numfaceedges;
3044 int fields;
3045 uint32_t i;
3046 int nseid = 0; /* number of signed edge ids */
3047 int prevseid;
3048 LWT_ELEMID *seid; /* signed edge ids */
3049
3050 /* Get list of face edges */
3051 numfaceedges = 1;
3052 fields = LWT_COL_EDGE_EDGE_ID |
3053 LWT_COL_EDGE_GEOM |
3054 LWT_COL_EDGE_FACE_LEFT |
3055 LWT_COL_EDGE_FACE_RIGHT
3056 ;
3057 edges = lwt_be_getEdgeByFace( topo, &face_id, &numfaceedges, fields, NULL );
3058 if (numfaceedges == UINT64_MAX)
3059 {
3060 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3061 return -1;
3062 }
3063 if ( ! numfaceedges ) return 0; /* no edges in output */
3064
3065 /* order edges by occurrence in face */
3066
3067 face = _lwt_FaceByEdges(topo, edges, numfaceedges);
3068 if ( ! face )
3069 {
3070 /* _lwt_FaceByEdges should have already invoked lwerror in this case */
3071 _lwt_release_edges(edges, numfaceedges);
3072 return -1;
3073 }
3074
3075 if ( lwgeom_is_empty(face) )
3076 {
3077 /* no edges in output */
3078 _lwt_release_edges(edges, numfaceedges);
3079 lwgeom_free(face);
3080 return 0;
3081 }
3082
3083 /* force_lhr, if the face is not the universe */
3084 /* _lwt_FaceByEdges seems to guaranteed RHR */
3085 /* lwgeom_force_clockwise(face); */
3086 if ( face_id ) lwgeom_reverse_in_place(face);
3087
3088 #if 0
3089 {
3090 size_t sz;
3091 char *wkt = lwgeom_to_wkt(face, WKT_EXTENDED, 6, &sz);
3092 LWDEBUGF(1, "Geometry of face %" LWTFMT_ELEMID " is: %s",
3093 face_id, wkt);
3094 lwfree(wkt);
3095 }
3096 #endif
3097
3098 facepoly = lwgeom_as_lwpoly(face);
3099 if ( ! facepoly )
3100 {
3101 _lwt_release_edges(edges, numfaceedges);
3102 lwgeom_free(face);
3103 lwerror("Geometry of face %" LWTFMT_ELEMID " is not a polygon", face_id);
3104 return -1;
3105 }
3106
3107 nseid = prevseid = 0;
3108 seid = lwalloc( sizeof(LWT_ELEMID) * numfaceedges );
3109
3110 /* for each ring of the face polygon... */
3111 for ( i=0; i<facepoly->nrings; ++i )
3112 {
3113 const POINTARRAY *ring = facepoly->rings[i];
3114 int32_t j = 0;
3115 LWT_ISO_EDGE *nextedge;
3116 LWLINE *nextline;
3117
3118 LWDEBUGF(1, "Ring %d has %d points", i, ring->npoints);
3119
3120 while ( j < (int32_t) ring->npoints-1 )
3121 {
3122 LWDEBUGF(1, "Looking for edge covering ring %d from vertex %d",
3123 i, j);
3124
3125 int edgeno = _lwt_FindNextRingEdge(ring, j, edges, numfaceedges);
3126 if ( edgeno == -1 )
3127 {
3128 /* should never happen */
3129 _lwt_release_edges(edges, numfaceedges);
3130 lwgeom_free(face);
3131 lwfree(seid);
3132 lwerror("No edge (among %d) found to be defining geometry of face %"
3133 LWTFMT_ELEMID, numfaceedges, face_id);
3134 return -1;
3135 }
3136
3137 nextedge = &(edges[edgeno]);
3138 nextline = nextedge->geom;
3139
3140 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
3141 " covers ring %d from vertex %d to %d",
3142 nextedge->edge_id, i, j, j + nextline->points->npoints - 1);
3143
3144 #if 0
3145 {
3146 size_t sz;
3147 char *wkt = lwgeom_to_wkt(lwline_as_lwgeom(nextline), WKT_EXTENDED, 6, &sz);
3148 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s",
3149 nextedge->edge_id, wkt);
3150 lwfree(wkt);
3151 }
3152 #endif
3153
3154 j += nextline->points->npoints - 1;
3155
3156 /* Add next edge to the output array */
3157 seid[nseid++] = nextedge->face_left == face_id ?
3158 nextedge->edge_id :
3159 -nextedge->edge_id;
3160
3161 /* avoid checking again on next time turn */
3162 nextedge->face_left = nextedge->face_right = -1;
3163 }
3164
3165 /* now "scroll" the list of edges so that the one
3166 * with smaller absolute edge_id is first */
3167 /* Range is: [prevseid, nseid) -- [inclusive, exclusive) */
3168 if ( (nseid - prevseid) > 1 )
3169 {{
3170 LWT_ELEMID minid = 0;
3171 int minidx = 0;
3172 LWDEBUGF(1, "Looking for smallest id among the %d edges "
3173 "composing ring %d", (nseid-prevseid), i);
3174 for ( j=prevseid; j<nseid; ++j )
3175 {
3176 LWT_ELEMID id = llabs(seid[j]);
3177 LWDEBUGF(1, "Abs id of edge in pos %d is %" LWTFMT_ELEMID, j, id);
3178 if ( ! minid || id < minid )
3179 {
3180 minid = id;
3181 minidx = j;
3182 }
3183 }
3184 LWDEBUGF(1, "Smallest id is %" LWTFMT_ELEMID
3185 " at position %d", minid, minidx);
3186 if ( minidx != prevseid )
3187 {
3188 _lwt_RotateElemidArray(seid, prevseid, nseid, minidx);
3189 }
3190 }}
3191
3192 prevseid = nseid;
3193 }
3194
3195 lwgeom_free(face);
3196 _lwt_release_edges(edges, numfaceedges);
3197
3198 *out = seid;
3199 return nseid;
3200 }
3201
3202 int
lwt_ChangeEdgeGeom(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id,LWLINE * geom)3203 lwt_ChangeEdgeGeom(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, LWLINE *geom)
3204 {
3205 LWT_ISO_EDGE *oldedge;
3206 LWT_ISO_EDGE newedge;
3207 POINT2D p1, p2, pt;
3208 uint64_t i;
3209 int isclosed = 0;
3210
3211 /* curve must be simple */
3212 if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) )
3213 {
3214 lwerror("SQL/MM Spatial exception - curve not simple");
3215 return -1;
3216 }
3217
3218 i = 1;
3219 oldedge = lwt_be_getEdgeById(topo, &edge_id, &i, LWT_COL_EDGE_ALL);
3220 if ( ! oldedge )
3221 {
3222 LWDEBUGF(1, "lwt_ChangeEdgeGeom: "
3223 "lwt_be_getEdgeById returned NULL and set i=%d", i);
3224 if (i == UINT64_MAX)
3225 {
3226 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3227 return -1;
3228 }
3229 else if ( i == 0 )
3230 {
3231 lwerror("SQL/MM Spatial exception - non-existent edge %"
3232 LWTFMT_ELEMID, edge_id);
3233 return -1;
3234 }
3235 else
3236 {
3237 lwerror("Backend coding error: getEdgeById callback returned NULL "
3238 "but numelements output parameter has value %d "
3239 "(expected 0 or 1)", i);
3240 return -1;
3241 }
3242 }
3243
3244 LWDEBUGF(1, "lwt_ChangeEdgeGeom: "
3245 "old edge has %d points, new edge has %d points",
3246 oldedge->geom->points->npoints, geom->points->npoints);
3247
3248 /*
3249 * e) Check StartPoint consistency
3250 */
3251 getPoint2d_p(oldedge->geom->points, 0, &p1);
3252 getPoint2d_p(geom->points, 0, &pt);
3253 if ( ! p2d_same(&p1, &pt) )
3254 {
3255 _lwt_release_edges(oldedge, 1);
3256 lwerror("SQL/MM Spatial exception - "
3257 "start node not geometry start point.");
3258 return -1;
3259 }
3260
3261 /*
3262 * f) Check EndPoint consistency
3263 */
3264 if ( oldedge->geom->points->npoints < 2 )
3265 {
3266 _lwt_release_edges(oldedge, 1);
3267 lwerror("Corrupted topology: edge %" LWTFMT_ELEMID
3268 " has less than 2 vertices", oldedge->edge_id);
3269 return -1;
3270 }
3271 getPoint2d_p(oldedge->geom->points, oldedge->geom->points->npoints-1, &p2);
3272 if ( geom->points->npoints < 2 )
3273 {
3274 _lwt_release_edges(oldedge, 1);
3275 lwerror("Invalid edge: less than 2 vertices");
3276 return -1;
3277 }
3278 getPoint2d_p(geom->points, geom->points->npoints-1, &pt);
3279 if ( ! p2d_same(&pt, &p2) )
3280 {
3281 _lwt_release_edges(oldedge, 1);
3282 lwerror("SQL/MM Spatial exception - "
3283 "end node not geometry end point.");
3284 return -1;
3285 }
3286
3287 /* Not in the specs:
3288 * if the edge is closed, check we didn't change winding !
3289 * (should be part of isomorphism checking)
3290 */
3291 if ( oldedge->start_node == oldedge->end_node )
3292 {
3293 isclosed = 1;
3294 #if 1 /* TODO: this is actually bogus as a test */
3295 /* check for valid edge (distinct vertices must exist) */
3296 if ( ! _lwt_GetInteriorEdgePoint(geom, &pt) )
3297 {
3298 _lwt_release_edges(oldedge, 1);
3299 lwerror("Invalid edge (no two distinct vertices exist)");
3300 return -1;
3301 }
3302 #endif
3303
3304 if ( ptarray_isccw(oldedge->geom->points) !=
3305 ptarray_isccw(geom->points) )
3306 {
3307 _lwt_release_edges(oldedge, 1);
3308 lwerror("Edge twist at node POINT(%g %g)", p1.x, p1.y);
3309 return -1;
3310 }
3311 }
3312
3313 if ( _lwt_CheckEdgeCrossing(topo, oldedge->start_node,
3314 oldedge->end_node, geom, edge_id ) )
3315 {
3316 /* would have called lwerror already, leaking :( */
3317 _lwt_release_edges(oldedge, 1);
3318 return -1;
3319 }
3320
3321 LWDEBUG(1, "lwt_ChangeEdgeGeom: "
3322 "edge crossing check passed ");
3323
3324 /*
3325 * Not in the specs:
3326 * Check topological isomorphism
3327 */
3328
3329 /* Check that the "motion range" doesn't include any node */
3330 // 1. compute combined bbox of old and new edge
3331 GBOX mbox; /* motion box */
3332 lwgeom_add_bbox((LWGEOM*)oldedge->geom); /* just in case */
3333 lwgeom_add_bbox((LWGEOM*)geom); /* just in case */
3334 gbox_union(oldedge->geom->bbox, geom->bbox, &mbox);
3335 // 2. fetch all nodes in the combined box
3336 LWT_ISO_NODE *nodes;
3337 uint64_t numnodes;
3338 nodes = lwt_be_getNodeWithinBox2D(topo, &mbox, &numnodes,
3339 LWT_COL_NODE_ALL, 0);
3340 LWDEBUGF(1, "lwt_be_getNodeWithinBox2D returned %d nodes", numnodes);
3341 if (numnodes == UINT64_MAX)
3342 {
3343 _lwt_release_edges(oldedge, 1);
3344 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3345 return -1;
3346 }
3347 // 3. if any node beside endnodes are found:
3348 if ( numnodes > ( 1 + isclosed ? 0 : 1 ) )
3349 {{
3350 // 3.2. bail out if any node is in one and not the other
3351 for (i=0; i<numnodes; ++i)
3352 {
3353 LWT_ISO_NODE *n = &(nodes[i]);
3354 if ( n->node_id == oldedge->start_node ) continue;
3355 if ( n->node_id == oldedge->end_node ) continue;
3356 const POINT2D *pt = getPoint2d_cp(n->geom->point, 0);
3357 int ocont = ptarray_contains_point_partial(oldedge->geom->points, pt, isclosed, NULL) == LW_INSIDE;
3358 int ncont = ptarray_contains_point_partial(geom->points, pt, isclosed, NULL) == LW_INSIDE;
3359 if (ocont != ncont)
3360 {
3361 size_t sz;
3362 char *wkt = lwgeom_to_wkt(lwpoint_as_lwgeom(n->geom), WKT_ISO, 15, &sz);
3363 _lwt_release_nodes(nodes, numnodes);
3364 lwerror("Edge motion collision at %s", wkt);
3365 lwfree(wkt); /* would not necessarely reach this point */
3366 return -1;
3367 }
3368 }
3369 }}
3370 if ( numnodes ) _lwt_release_nodes(nodes, numnodes);
3371
3372 LWDEBUG(1, "nodes containment check passed");
3373
3374 /*
3375 * Check edge adjacency before
3376 * TODO: can be optimized to gather azimuths of all edge ends once
3377 */
3378
3379 edgeend span_pre, epan_pre;
3380 /* initialize span_pre.myaz and epan_pre.myaz with existing edge */
3381 int res = _lwt_InitEdgeEndByLine(&span_pre, &epan_pre, oldedge->geom, &p1, &p2);
3382 if (res)
3383 return -1; /* lwerror should have been raised */
3384 _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_pre,
3385 isclosed ? &epan_pre : NULL, edge_id );
3386 _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_pre,
3387 isclosed ? &span_pre : NULL, edge_id );
3388
3389 LWDEBUGF(1, "edges adjacent to old edge are %" LWTFMT_ELEMID
3390 " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID
3391 " and %" LWTFMT_ELEMID " (last point)",
3392 span_pre.nextCW, span_pre.nextCCW,
3393 epan_pre.nextCW, epan_pre.nextCCW);
3394
3395 /* update edge geometry */
3396 newedge.edge_id = edge_id;
3397 newedge.geom = geom;
3398 res = lwt_be_updateEdgesById(topo, &newedge, 1, LWT_COL_EDGE_GEOM);
3399 if (res == -1)
3400 {
3401 _lwt_release_edges(oldedge, 1);
3402 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3403 return -1;
3404 }
3405 if (!res)
3406 {
3407 _lwt_release_edges(oldedge, 1);
3408 lwerror("Unexpected error: %d edges updated when expecting 1", i);
3409 return -1;
3410 }
3411
3412 /*
3413 * Check edge adjacency after
3414 */
3415 edgeend span_post, epan_post;
3416 /* initialize epan_post.myaz and epan_post.myaz */
3417 res = _lwt_InitEdgeEndByLine(&span_post, &epan_post, geom, &p1, &p2);
3418 if (res)
3419 return -1; /* lwerror should have been raised */
3420 _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_post,
3421 isclosed ? &epan_post : NULL, edge_id );
3422 _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_post,
3423 isclosed ? &span_post : NULL, edge_id );
3424
3425 LWDEBUGF(1, "edges adjacent to new edge are %" LWTFMT_ELEMID
3426 " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID
3427 " and %" LWTFMT_ELEMID " (last point)",
3428 span_pre.nextCW, span_pre.nextCCW,
3429 epan_pre.nextCW, epan_pre.nextCCW);
3430
3431
3432 /* Bail out if next CW or CCW edge on start node changed */
3433 if ( span_pre.nextCW != span_post.nextCW ||
3434 span_pre.nextCCW != span_post.nextCCW )
3435 {{
3436 LWT_ELEMID nid = oldedge->start_node;
3437 _lwt_release_edges(oldedge, 1);
3438 lwerror("Edge changed disposition around start node %"
3439 LWTFMT_ELEMID, nid);
3440 return -1;
3441 }}
3442
3443 /* Bail out if next CW or CCW edge on end node changed */
3444 if ( epan_pre.nextCW != epan_post.nextCW ||
3445 epan_pre.nextCCW != epan_post.nextCCW )
3446 {{
3447 LWT_ELEMID nid = oldedge->end_node;
3448 _lwt_release_edges(oldedge, 1);
3449 lwerror("Edge changed disposition around end node %"
3450 LWTFMT_ELEMID, nid);
3451 return -1;
3452 }}
3453
3454 /*
3455 -- Update faces MBR of left and right faces
3456 -- TODO: think about ways to optimize this part, like see if
3457 -- the old edge geometry participated in the definition
3458 -- of the current MBR (for shrinking) or the new edge MBR
3459 -- would be larger than the old face MBR...
3460 --
3461 */
3462 const GBOX* oldbox = lwgeom_get_bbox(lwline_as_lwgeom(oldedge->geom));
3463 const GBOX* newbox = lwgeom_get_bbox(lwline_as_lwgeom(geom));
3464 if ( ! gbox_same(oldbox, newbox) )
3465 {
3466 uint64_t facestoupdate = 0;
3467 LWT_ISO_FACE faces[2];
3468 LWGEOM *nface1 = NULL;
3469 LWGEOM *nface2 = NULL;
3470 if ( oldedge->face_left > 0 )
3471 {
3472 nface1 = lwt_GetFaceGeometry(topo, oldedge->face_left);
3473 if ( ! nface1 )
3474 {
3475 lwerror("lwt_ChangeEdgeGeom could not construct face %"
3476 LWTFMT_ELEMID ", on the left of edge %" LWTFMT_ELEMID,
3477 oldedge->face_left, edge_id);
3478 return -1;
3479 }
3480 #if 0
3481 {
3482 size_t sz;
3483 char *wkt = lwgeom_to_wkt(nface1, WKT_EXTENDED, 2, &sz);
3484 LWDEBUGF(1, "new geometry of face left (%d): %s", (int)oldedge->face_left, wkt);
3485 lwfree(wkt);
3486 }
3487 #endif
3488 lwgeom_add_bbox(nface1);
3489 if ( ! nface1->bbox )
3490 {
3491 lwerror("Corrupted topology: face %d, left of edge %d, has no bbox",
3492 oldedge->face_left, edge_id);
3493 return -1;
3494 }
3495 faces[facestoupdate].face_id = oldedge->face_left;
3496 /* ownership left to nface */
3497 faces[facestoupdate++].mbr = nface1->bbox;
3498 }
3499 if ( oldedge->face_right > 0
3500 /* no need to update twice the same face.. */
3501 && oldedge->face_right != oldedge->face_left )
3502 {
3503 nface2 = lwt_GetFaceGeometry(topo, oldedge->face_right);
3504 if ( ! nface2 )
3505 {
3506 lwerror("lwt_ChangeEdgeGeom could not construct face %"
3507 LWTFMT_ELEMID ", on the right of edge %" LWTFMT_ELEMID,
3508 oldedge->face_right, edge_id);
3509 return -1;
3510 }
3511 #if 0
3512 {
3513 size_t sz;
3514 char *wkt = lwgeom_to_wkt(nface2, WKT_EXTENDED, 2, &sz);
3515 LWDEBUGF(1, "new geometry of face right (%d): %s", (int)oldedge->face_right, wkt);
3516 lwfree(wkt);
3517 }
3518 #endif
3519 lwgeom_add_bbox(nface2);
3520 if ( ! nface2->bbox )
3521 {
3522 lwerror("Corrupted topology: face %d, right of edge %d, has no bbox",
3523 oldedge->face_right, edge_id);
3524 return -1;
3525 }
3526 faces[facestoupdate].face_id = oldedge->face_right;
3527 faces[facestoupdate++].mbr = nface2->bbox; /* ownership left to nface */
3528 }
3529 LWDEBUGF(1, "%d faces to update", facestoupdate);
3530 if ( facestoupdate )
3531 {
3532 uint64_t updatedFaces = lwt_be_updateFacesById(topo, &(faces[0]), facestoupdate);
3533 if (updatedFaces != facestoupdate)
3534 {
3535 if (nface1)
3536 lwgeom_free(nface1);
3537 if (nface2)
3538 lwgeom_free(nface2);
3539 _lwt_release_edges(oldedge, 1);
3540 if (updatedFaces == UINT64_MAX)
3541 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3542 else
3543 lwerror("Unexpected error: %d faces found when expecting 1", i);
3544 return -1;
3545 }
3546 }
3547 if ( nface1 ) lwgeom_free(nface1);
3548 if ( nface2 ) lwgeom_free(nface2);
3549 }
3550 else
3551 {
3552 LWDEBUG(1, "BBOX of changed edge did not change");
3553 }
3554
3555 LWDEBUG(1, "all done, cleaning up edges");
3556
3557 _lwt_release_edges(oldedge, 1);
3558 return 0; /* success */
3559 }
3560
3561 /* Only return CONTAINING_FACE in the node object */
3562 static LWT_ISO_NODE *
_lwt_GetIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid)3563 _lwt_GetIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid)
3564 {
3565 LWT_ISO_NODE *node;
3566 uint64_t n = 1;
3567
3568 node = lwt_be_getNodeById( topo, &nid, &n, LWT_COL_NODE_CONTAINING_FACE );
3569 if (n == UINT64_MAX)
3570 {
3571 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3572 return 0;
3573 }
3574 if ( n < 1 ) {
3575 lwerror("SQL/MM Spatial exception - non-existent node");
3576 return 0;
3577 }
3578 if ( node->containing_face == -1 )
3579 {
3580 lwfree(node);
3581 lwerror("SQL/MM Spatial exception - not isolated node");
3582 return 0;
3583 }
3584
3585 return node;
3586 }
3587
3588 int
lwt_MoveIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid,LWPOINT * pt)3589 lwt_MoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid, LWPOINT *pt)
3590 {
3591 LWT_ISO_NODE *node;
3592 int ret;
3593
3594 node = _lwt_GetIsoNode( topo, nid );
3595 if ( ! node ) return -1;
3596
3597 if ( lwt_be_ExistsCoincidentNode(topo, pt) )
3598 {
3599 lwfree(node);
3600 lwerror("SQL/MM Spatial exception - coincident node");
3601 return -1;
3602 }
3603
3604 if ( lwt_be_ExistsEdgeIntersectingPoint(topo, pt) )
3605 {
3606 lwfree(node);
3607 lwerror("SQL/MM Spatial exception - edge crosses node.");
3608 return -1;
3609 }
3610
3611 /* TODO: check that the new point is in the same containing face !
3612 * See https://trac.osgeo.org/postgis/ticket/3232
3613 */
3614
3615 node->node_id = nid;
3616 node->geom = pt;
3617 ret = lwt_be_updateNodesById(topo, node, 1,
3618 LWT_COL_NODE_GEOM);
3619 if ( ret == -1 ) {
3620 lwfree(node);
3621 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3622 return -1;
3623 }
3624
3625 lwfree(node);
3626 return 0;
3627 }
3628
3629 int
lwt_RemoveIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid)3630 lwt_RemoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid)
3631 {
3632 LWT_ISO_NODE *node;
3633 int n = 1;
3634
3635 node = _lwt_GetIsoNode( topo, nid );
3636 if ( ! node ) return -1;
3637
3638 n = lwt_be_deleteNodesById( topo, &nid, n );
3639 if ( n == -1 )
3640 {
3641 lwfree(node);
3642 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3643 return -1;
3644 }
3645 if ( n != 1 )
3646 {
3647 lwfree(node);
3648 lwerror("Unexpected error: %d nodes deleted when expecting 1", n);
3649 return -1;
3650 }
3651
3652 /* TODO: notify to caller about node being removed ?
3653 * See https://trac.osgeo.org/postgis/ticket/3231
3654 */
3655
3656 lwfree(node);
3657 return 0; /* success */
3658 }
3659
3660 int
lwt_RemIsoEdge(LWT_TOPOLOGY * topo,LWT_ELEMID id)3661 lwt_RemIsoEdge(LWT_TOPOLOGY* topo, LWT_ELEMID id)
3662 {
3663 LWT_ISO_EDGE deledge;
3664 LWT_ISO_EDGE *edge;
3665 LWT_ELEMID nid[2];
3666 LWT_ISO_NODE upd_node[2];
3667 LWT_ELEMID containing_face;
3668 uint64_t n = 1;
3669 uint64_t i;
3670
3671 edge = lwt_be_getEdgeById( topo, &id, &n, LWT_COL_EDGE_START_NODE|
3672 LWT_COL_EDGE_END_NODE |
3673 LWT_COL_EDGE_FACE_LEFT |
3674 LWT_COL_EDGE_FACE_RIGHT );
3675 if ( ! edge )
3676 {
3677 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3678 return -1;
3679 }
3680 if ( ! n )
3681 {
3682 lwerror("SQL/MM Spatial exception - non-existent edge");
3683 return -1;
3684 }
3685 if ( n > 1 )
3686 {
3687 lwfree(edge);
3688 lwerror("Corrupted topology: more than a single edge have id %"
3689 LWTFMT_ELEMID, id);
3690 return -1;
3691 }
3692
3693 if ( edge[0].face_left != edge[0].face_right )
3694 {
3695 lwfree(edge);
3696 lwerror("SQL/MM Spatial exception - not isolated edge");
3697 return -1;
3698 }
3699 containing_face = edge[0].face_left;
3700
3701 nid[0] = edge[0].start_node;
3702 nid[1] = edge[0].end_node;
3703 lwfree(edge);
3704
3705 n = 2;
3706 edge = lwt_be_getEdgeByNode( topo, nid, &n, LWT_COL_EDGE_EDGE_ID );
3707 if ((n == UINT64_MAX) || (edge == NULL))
3708 {
3709 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3710 return -1;
3711 }
3712 for (i = 0; i < n; ++i)
3713 {
3714 if (edge[i].edge_id != id)
3715 {
3716 lwfree(edge);
3717 lwerror("SQL/MM Spatial exception - not isolated edge");
3718 return -1;
3719 }
3720 }
3721 lwfree(edge);
3722
3723 deledge.edge_id = id;
3724 n = lwt_be_deleteEdges( topo, &deledge, LWT_COL_EDGE_EDGE_ID );
3725 if (n == UINT64_MAX)
3726 {
3727 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3728 return -1;
3729 }
3730 if ( n != 1 )
3731 {
3732 lwerror("Unexpected error: %d edges deleted when expecting 1", n);
3733 return -1;
3734 }
3735
3736 upd_node[0].node_id = nid[0];
3737 upd_node[0].containing_face = containing_face;
3738 n = 1;
3739 if ( nid[1] != nid[0] ) {
3740 upd_node[1].node_id = nid[1];
3741 upd_node[1].containing_face = containing_face;
3742 ++n;
3743 }
3744 n = lwt_be_updateNodesById(topo, upd_node, n,
3745 LWT_COL_NODE_CONTAINING_FACE);
3746 if (n == UINT64_MAX)
3747 {
3748 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3749 return -1;
3750 }
3751
3752 /* TODO: notify to caller about edge being removed ?
3753 * See https://trac.osgeo.org/postgis/ticket/3248
3754 */
3755
3756 return 0; /* success */
3757 }
3758
3759 /* Used by _lwt_RemEdge to update edge face ref on healing
3760 *
3761 * @param of old face id (never 0 as you cannot remove face 0)
3762 * @param nf new face id
3763 * @return 0 on success, -1 on backend error
3764 */
3765 static int
_lwt_UpdateEdgeFaceRef(LWT_TOPOLOGY * topo,LWT_ELEMID of,LWT_ELEMID nf)3766 _lwt_UpdateEdgeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf)
3767 {
3768 LWT_ISO_EDGE sel_edge, upd_edge;
3769 int ret;
3770
3771 assert( of != 0 );
3772
3773 /* Update face_left for all edges still referencing old face */
3774 sel_edge.face_left = of;
3775 upd_edge.face_left = nf;
3776 ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_LEFT,
3777 &upd_edge, LWT_COL_EDGE_FACE_LEFT,
3778 NULL, 0);
3779 if ( ret == -1 ) return -1;
3780
3781 /* Update face_right for all edges still referencing old face */
3782 sel_edge.face_right = of;
3783 upd_edge.face_right = nf;
3784 ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_RIGHT,
3785 &upd_edge, LWT_COL_EDGE_FACE_RIGHT,
3786 NULL, 0);
3787 if ( ret == -1 ) return -1;
3788
3789 return 0;
3790 }
3791
3792 /* Used by _lwt_RemEdge to update node face ref on healing
3793 *
3794 * @param of old face id (never 0 as you cannot remove face 0)
3795 * @param nf new face id
3796 * @return 0 on success, -1 on backend error
3797 */
3798 static int
_lwt_UpdateNodeFaceRef(LWT_TOPOLOGY * topo,LWT_ELEMID of,LWT_ELEMID nf)3799 _lwt_UpdateNodeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf)
3800 {
3801 LWT_ISO_NODE sel, upd;
3802 int ret;
3803
3804 assert( of != 0 );
3805
3806 /* Update face_left for all edges still referencing old face */
3807 sel.containing_face = of;
3808 upd.containing_face = nf;
3809 ret = lwt_be_updateNodes(topo, &sel, LWT_COL_NODE_CONTAINING_FACE,
3810 &upd, LWT_COL_NODE_CONTAINING_FACE,
3811 NULL, 0);
3812 if ( ret == -1 ) return -1;
3813
3814 return 0;
3815 }
3816
3817 /* Used by lwt_RemEdgeModFace and lwt_RemEdgeNewFaces
3818 *
3819 * Returns -1 on error, identifier of the face that takes up the space
3820 * previously occupied by the removed edge if modFace is 1, identifier of
3821 * the created face (0 if none) if modFace is 0.
3822 */
3823 static LWT_ELEMID
_lwt_RemEdge(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id,int modFace)3824 _lwt_RemEdge( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, int modFace )
3825 {
3826 uint64_t i, nedges, nfaces, fields;
3827 LWT_ISO_EDGE *edge = NULL;
3828 LWT_ISO_EDGE *upd_edge = NULL;
3829 LWT_ISO_EDGE upd_edge_left[2];
3830 int nedge_left = 0;
3831 LWT_ISO_EDGE upd_edge_right[2];
3832 int nedge_right = 0;
3833 LWT_ISO_NODE upd_node[2];
3834 int nnode = 0;
3835 LWT_ISO_FACE *faces = NULL;
3836 LWT_ISO_FACE newface;
3837 LWT_ELEMID node_ids[2];
3838 LWT_ELEMID face_ids[2];
3839 int fnode_edges = 0; /* number of edges on the first node (excluded
3840 * the one being removed ) */
3841 int lnode_edges = 0; /* number of edges on the last node (excluded
3842 * the one being removed ) */
3843
3844 newface.face_id = 0;
3845
3846 i = 1;
3847 edge = lwt_be_getEdgeById(topo, &edge_id, &i, LWT_COL_EDGE_ALL);
3848 if (!edge)
3849 {
3850 LWDEBUGF(1, "lwt_be_getEdgeById returned NULL and set i=%d", i);
3851 if (i == UINT64_MAX)
3852 {
3853 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3854 return -1;
3855 }
3856 else if (i == 0)
3857 {
3858 lwerror("SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID, edge_id);
3859 return -1;
3860 }
3861 else
3862 {
3863 lwerror(
3864 "Backend coding error: getEdgeById callback returned NULL "
3865 "but numelements output parameter has value %d "
3866 "(expected 0 or 1)",
3867 i);
3868 return -1;
3869 }
3870 }
3871
3872 if ( ! lwt_be_checkTopoGeomRemEdge(topo, edge_id,
3873 edge->face_left, edge->face_right) )
3874 {
3875 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
3876 return -1;
3877 }
3878
3879 LWDEBUG(1, "Updating next_{right,left}_face of ring edges...");
3880
3881 /* Update edge linking */
3882
3883 nedges = 0;
3884 node_ids[nedges++] = edge->start_node;
3885 if ( edge->end_node != edge->start_node )
3886 {
3887 node_ids[nedges++] = edge->end_node;
3888 }
3889 fields = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_START_NODE |
3890 LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_NEXT_LEFT |
3891 LWT_COL_EDGE_NEXT_RIGHT;
3892 upd_edge = lwt_be_getEdgeByNode( topo, &(node_ids[0]), &nedges, fields );
3893 if (nedges == UINT64_MAX)
3894 {
3895 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3896 return -1;
3897 }
3898 nedge_left = nedge_right = 0;
3899 for ( i=0; i<nedges; ++i )
3900 {
3901 LWT_ISO_EDGE *e = &(upd_edge[i]);
3902 if ( e->edge_id == edge_id ) continue;
3903 if ( e->start_node == edge->start_node || e->end_node == edge->start_node )
3904 {
3905 ++fnode_edges;
3906 }
3907 if ( e->start_node == edge->end_node || e->end_node == edge->end_node )
3908 {
3909 ++lnode_edges;
3910 }
3911 if ( e->next_left == -edge_id )
3912 {
3913 upd_edge_left[nedge_left].edge_id = e->edge_id;
3914 upd_edge_left[nedge_left++].next_left =
3915 edge->next_left != edge_id ? edge->next_left : edge->next_right;
3916 }
3917 else if ( e->next_left == edge_id )
3918 {
3919 upd_edge_left[nedge_left].edge_id = e->edge_id;
3920 upd_edge_left[nedge_left++].next_left =
3921 edge->next_right != -edge_id ? edge->next_right : edge->next_left;
3922 }
3923
3924 if ( e->next_right == -edge_id )
3925 {
3926 upd_edge_right[nedge_right].edge_id = e->edge_id;
3927 upd_edge_right[nedge_right++].next_right =
3928 edge->next_left != edge_id ? edge->next_left : edge->next_right;
3929 }
3930 else if ( e->next_right == edge_id )
3931 {
3932 upd_edge_right[nedge_right].edge_id = e->edge_id;
3933 upd_edge_right[nedge_right++].next_right =
3934 edge->next_right != -edge_id ? edge->next_right : edge->next_left;
3935 }
3936 }
3937
3938 if ( nedge_left )
3939 {
3940 LWDEBUGF(1, "updating %d 'next_left' edges", nedge_left);
3941 /* update edges in upd_edge_left set next_left */
3942 int result = lwt_be_updateEdgesById(topo, &(upd_edge_left[0]), nedge_left, LWT_COL_EDGE_NEXT_LEFT);
3943 if (result == -1)
3944 {
3945 _lwt_release_edges(edge, 1);
3946 lwfree(upd_edge);
3947 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3948 return -1;
3949 }
3950 }
3951 if ( nedge_right )
3952 {
3953 LWDEBUGF(1, "updating %d 'next_right' edges", nedge_right);
3954 /* update edges in upd_edge_right set next_right */
3955 int result = lwt_be_updateEdgesById(topo, &(upd_edge_right[0]), nedge_right, LWT_COL_EDGE_NEXT_RIGHT);
3956 if (result == -1)
3957 {
3958 _lwt_release_edges(edge, 1);
3959 lwfree(upd_edge);
3960 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3961 return -1;
3962 }
3963 }
3964 LWDEBUGF(1, "releasing %d updateable edges in %p", nedges, upd_edge);
3965 lwfree(upd_edge);
3966
3967 /* Id of face that will take up all the space previously
3968 * taken by left and right faces of the edge */
3969 LWT_ELEMID floodface;
3970
3971 /* Find floodface, and update its mbr if != 0 */
3972 if ( edge->face_left == edge->face_right )
3973 {
3974 floodface = edge->face_right;
3975 }
3976 else
3977 {
3978 /* Two faces healed */
3979 if ( edge->face_left == 0 || edge->face_right == 0 )
3980 {
3981 floodface = 0;
3982 LWDEBUG(1, "floodface is universe");
3983 }
3984 else
3985 {
3986 /* we choose right face as the face that will remain
3987 * to be symmetric with ST_AddEdgeModFace */
3988 floodface = edge->face_right;
3989 LWDEBUGF(1, "floodface is %" LWTFMT_ELEMID, floodface);
3990 /* update mbr of floodface as union of mbr of both faces */
3991 face_ids[0] = edge->face_left;
3992 face_ids[1] = edge->face_right;
3993 nfaces = 2;
3994 fields = LWT_COL_FACE_ALL;
3995 faces = lwt_be_getFaceById(topo, face_ids, &nfaces, fields);
3996 if (nfaces == UINT64_MAX)
3997 {
3998 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3999 return -1;
4000 }
4001 GBOX *box1=NULL;
4002 GBOX *box2=NULL;
4003 for ( i=0; i<nfaces; ++i )
4004 {
4005 if ( faces[i].face_id == edge->face_left )
4006 {
4007 if ( ! box1 ) box1 = faces[i].mbr;
4008 else
4009 {
4010 i = edge->face_left;
4011 _lwt_release_edges(edge, 1);
4012 _lwt_release_faces(faces, nfaces);
4013 lwerror("corrupted topology: more than 1 face have face_id=%"
4014 LWTFMT_ELEMID, i);
4015 return -1;
4016 }
4017 }
4018 else if ( faces[i].face_id == edge->face_right )
4019 {
4020 if ( ! box2 ) box2 = faces[i].mbr;
4021 else
4022 {
4023 i = edge->face_right;
4024 _lwt_release_edges(edge, 1);
4025 _lwt_release_faces(faces, nfaces);
4026 lwerror("corrupted topology: more than 1 face have face_id=%"
4027 LWTFMT_ELEMID, i);
4028 return -1;
4029 }
4030 }
4031 else
4032 {
4033 i = faces[i].face_id;
4034 _lwt_release_edges(edge, 1);
4035 _lwt_release_faces(faces, nfaces);
4036 lwerror("Backend coding error: getFaceById returned face "
4037 "with non-requested id %" LWTFMT_ELEMID, i);
4038 return -1;
4039 }
4040 }
4041 if ( ! box1 ) {
4042 i = edge->face_left;
4043 _lwt_release_edges(edge, 1);
4044 _lwt_release_faces(faces, nfaces);
4045 lwerror("corrupted topology: no face have face_id=%"
4046 LWTFMT_ELEMID " (left face for edge %"
4047 LWTFMT_ELEMID ")", i, edge_id);
4048 return -1;
4049 }
4050 if ( ! box2 ) {
4051 i = edge->face_right;
4052 _lwt_release_edges(edge, 1);
4053 _lwt_release_faces(faces, nfaces);
4054 lwerror("corrupted topology: no face have face_id=%"
4055 LWTFMT_ELEMID " (right face for edge %"
4056 LWTFMT_ELEMID ")", i, edge_id);
4057 return -1;
4058 }
4059 gbox_merge(box2, box1); /* box1 is now the union of the two */
4060 newface.mbr = box1;
4061 if ( modFace )
4062 {
4063 newface.face_id = floodface;
4064 int result = lwt_be_updateFacesById(topo, &newface, 1);
4065 _lwt_release_faces(faces, 2);
4066 if (result == -1)
4067 {
4068 _lwt_release_edges(edge, 1);
4069 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4070 return -1;
4071 }
4072 if (result != 1)
4073 {
4074 _lwt_release_edges(edge, 1);
4075 lwerror("Unexpected error: %d faces updated when expecting 1", i);
4076 return -1;
4077 }
4078 }
4079 else
4080 {
4081 /* New face replaces the old two faces */
4082 newface.face_id = -1;
4083 int result = lwt_be_insertFaces(topo, &newface, 1);
4084 _lwt_release_faces(faces, 2);
4085 if (result == -1)
4086 {
4087 _lwt_release_edges(edge, 1);
4088 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4089 return -1;
4090 }
4091 if (result != 1)
4092 {
4093 _lwt_release_edges(edge, 1);
4094 lwerror("Unexpected error: %d faces inserted when expecting 1", result);
4095 return -1;
4096 }
4097 floodface = newface.face_id;
4098 }
4099 }
4100
4101 /* Update face references for edges and nodes still referencing
4102 * the removed face(s) */
4103
4104 if ( edge->face_left != floodface )
4105 {
4106 if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_left, floodface) )
4107 {
4108 _lwt_release_edges(edge, 1);
4109 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4110 return -1;
4111 }
4112 if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_left, floodface) )
4113 {
4114 _lwt_release_edges(edge, 1);
4115 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4116 return -1;
4117 }
4118 }
4119
4120 if ( edge->face_right != floodface )
4121 {
4122 if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_right, floodface) )
4123 {
4124 _lwt_release_edges(edge, 1);
4125 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4126 return -1;
4127 }
4128 if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_right, floodface) )
4129 {
4130 _lwt_release_edges(edge, 1);
4131 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4132 return -1;
4133 }
4134 }
4135
4136 /* Update topogeoms on heal */
4137 if ( ! lwt_be_updateTopoGeomFaceHeal(topo,
4138 edge->face_right, edge->face_left,
4139 floodface) )
4140 {
4141 _lwt_release_edges(edge, 1);
4142 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4143 return -1;
4144 }
4145 } /* two faces healed */
4146
4147 /* Delete the edge */
4148 int result = lwt_be_deleteEdges(topo, edge, LWT_COL_EDGE_EDGE_ID);
4149 if (result == -1)
4150 {
4151 _lwt_release_edges(edge, 1);
4152 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4153 return -1;
4154 }
4155
4156 /* If any of the edge nodes remained isolated, set
4157 * containing_face = floodface
4158 */
4159 if ( ! fnode_edges )
4160 {
4161 upd_node[nnode].node_id = edge->start_node;
4162 upd_node[nnode].containing_face = floodface;
4163 ++nnode;
4164 }
4165 if ( edge->end_node != edge->start_node && ! lnode_edges )
4166 {
4167 upd_node[nnode].node_id = edge->end_node;
4168 upd_node[nnode].containing_face = floodface;
4169 ++nnode;
4170 }
4171 if ( nnode )
4172 {
4173 int result = lwt_be_updateNodesById(topo, upd_node, nnode, LWT_COL_NODE_CONTAINING_FACE);
4174 if (result == -1)
4175 {
4176 _lwt_release_edges(edge, 1);
4177 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4178 return -1;
4179 }
4180 }
4181
4182 if ( edge->face_left != edge->face_right )
4183 /* or there'd be no face to remove */
4184 {
4185 LWT_ELEMID ids[2];
4186 int nids = 0;
4187 if ( edge->face_right != floodface )
4188 ids[nids++] = edge->face_right;
4189 if ( edge->face_left != floodface )
4190 ids[nids++] = edge->face_left;
4191 int result = lwt_be_deleteFacesById(topo, ids, nids);
4192 if (result == -1)
4193 {
4194 _lwt_release_edges(edge, 1);
4195 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4196 return -1;
4197 }
4198 }
4199
4200 _lwt_release_edges(edge, 1);
4201 return modFace ? floodface : newface.face_id;
4202 }
4203
4204 LWT_ELEMID
lwt_RemEdgeModFace(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id)4205 lwt_RemEdgeModFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id )
4206 {
4207 return _lwt_RemEdge( topo, edge_id, 1 );
4208 }
4209
4210 LWT_ELEMID
lwt_RemEdgeNewFace(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id)4211 lwt_RemEdgeNewFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id )
4212 {
4213 return _lwt_RemEdge( topo, edge_id, 0 );
4214 }
4215
4216 static LWT_ELEMID
_lwt_HealEdges(LWT_TOPOLOGY * topo,LWT_ELEMID eid1,LWT_ELEMID eid2,int modEdge)4217 _lwt_HealEdges( LWT_TOPOLOGY* topo, LWT_ELEMID eid1, LWT_ELEMID eid2,
4218 int modEdge )
4219 {
4220 LWT_ELEMID ids[2];
4221 LWT_ELEMID commonnode = -1;
4222 int caseno = 0;
4223 LWT_ISO_EDGE *node_edges;
4224 uint64_t num_node_edges;
4225 LWT_ISO_EDGE *edges;
4226 LWT_ISO_EDGE *e1 = NULL;
4227 LWT_ISO_EDGE *e2 = NULL;
4228 LWT_ISO_EDGE newedge, updedge, seledge;
4229 uint64_t nedges, i;
4230 int e1freenode;
4231 int e2sign, e2freenode;
4232 POINTARRAY *pa;
4233 char buf[256];
4234 char *ptr;
4235 size_t bufleft = 256;
4236
4237 ptr = buf;
4238
4239 /* NOT IN THE SPECS: see if the same edge is given twice.. */
4240 if ( eid1 == eid2 )
4241 {
4242 lwerror("Cannot heal edge %" LWTFMT_ELEMID
4243 " with itself, try with another", eid1);
4244 return -1;
4245 }
4246 ids[0] = eid1;
4247 ids[1] = eid2;
4248 nedges = 2;
4249 edges = lwt_be_getEdgeById(topo, ids, &nedges, LWT_COL_EDGE_ALL);
4250 if ((nedges == UINT64_MAX) || (edges == NULL))
4251 {
4252 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4253 return -1;
4254 }
4255 for ( i=0; i<nedges; ++i )
4256 {
4257 if ( edges[i].edge_id == eid1 ) {
4258 if ( e1 ) {
4259 _lwt_release_edges(edges, nedges);
4260 lwerror("Corrupted topology: multiple edges have id %"
4261 LWTFMT_ELEMID, eid1);
4262 return -1;
4263 }
4264 e1 = &(edges[i]);
4265 }
4266 else if ( edges[i].edge_id == eid2 ) {
4267 if ( e2 ) {
4268 _lwt_release_edges(edges, nedges);
4269 lwerror("Corrupted topology: multiple edges have id %"
4270 LWTFMT_ELEMID, eid2);
4271 return -1;
4272 }
4273 e2 = &(edges[i]);
4274 }
4275 }
4276 if ( ! e1 )
4277 {
4278 _lwt_release_edges(edges, nedges);
4279 lwerror(
4280 "SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID,
4281 eid1);
4282 return -1;
4283 }
4284 if ( ! e2 )
4285 {
4286 _lwt_release_edges(edges, nedges);
4287 lwerror(
4288 "SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID,
4289 eid2);
4290 return -1;
4291 }
4292
4293 /* NOT IN THE SPECS: See if any of the two edges are closed. */
4294 if ( e1->start_node == e1->end_node )
4295 {
4296 _lwt_release_edges(edges, nedges);
4297 lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %"
4298 LWTFMT_ELEMID, eid1, eid2);
4299 return -1;
4300 }
4301 if ( e2->start_node == e2->end_node )
4302 {
4303 _lwt_release_edges(edges, nedges);
4304 lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %"
4305 LWTFMT_ELEMID, eid2, eid1);
4306 return -1;
4307 }
4308
4309 /* Find common node */
4310
4311 if ( e1->end_node == e2->start_node )
4312 {
4313 commonnode = e1->end_node;
4314 caseno = 1;
4315 }
4316 else if ( e1->end_node == e2->end_node )
4317 {
4318 commonnode = e1->end_node;
4319 caseno = 2;
4320 }
4321 /* Check if any other edge is connected to the common node, if found */
4322 if ( commonnode != -1 )
4323 {
4324 num_node_edges = 1;
4325 node_edges = lwt_be_getEdgeByNode( topo, &commonnode,
4326 &num_node_edges, LWT_COL_EDGE_EDGE_ID );
4327 if (num_node_edges == UINT64_MAX)
4328 {
4329 _lwt_release_edges(edges, nedges);
4330 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4331 return -1;
4332 }
4333 for (i=0; i<num_node_edges; ++i)
4334 {
4335 int r;
4336 if ( node_edges[i].edge_id == eid1 ) continue;
4337 if ( node_edges[i].edge_id == eid2 ) continue;
4338 commonnode = -1;
4339 /* append to string, for error message */
4340 if ( bufleft > 0 ) {
4341 r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID,
4342 ( ptr==buf ? "" : "," ), node_edges[i].edge_id);
4343 if ( r >= (int) bufleft )
4344 {
4345 bufleft = 0;
4346 buf[252] = '.';
4347 buf[253] = '.';
4348 buf[254] = '.';
4349 buf[255] = '\0';
4350 }
4351 else
4352 {
4353 bufleft -= r;
4354 ptr += r;
4355 }
4356 }
4357 }
4358 lwfree(node_edges);
4359 }
4360
4361 if ( commonnode == -1 )
4362 {
4363 if ( e1->start_node == e2->start_node )
4364 {
4365 commonnode = e1->start_node;
4366 caseno = 3;
4367 }
4368 else if ( e1->start_node == e2->end_node )
4369 {
4370 commonnode = e1->start_node;
4371 caseno = 4;
4372 }
4373 /* Check if any other edge is connected to the common node, if found */
4374 if ( commonnode != -1 )
4375 {
4376 num_node_edges = 1;
4377 node_edges = lwt_be_getEdgeByNode( topo, &commonnode,
4378 &num_node_edges, LWT_COL_EDGE_EDGE_ID );
4379 if (num_node_edges == UINT64_MAX)
4380 {
4381 _lwt_release_edges(edges, nedges);
4382 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4383 return -1;
4384 }
4385 for (i=0; i<num_node_edges; ++i)
4386 {
4387 int r;
4388 if ( node_edges[i].edge_id == eid1 ) continue;
4389 if ( node_edges[i].edge_id == eid2 ) continue;
4390 commonnode = -1;
4391 /* append to string, for error message */
4392 if ( bufleft > 0 ) {
4393 r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID,
4394 ( ptr==buf ? "" : "," ), node_edges[i].edge_id);
4395 if ( r >= (int) bufleft )
4396 {
4397 bufleft = 0;
4398 buf[252] = '.';
4399 buf[253] = '.';
4400 buf[254] = '.';
4401 buf[255] = '\0';
4402 }
4403 else
4404 {
4405 bufleft -= r;
4406 ptr += r;
4407 }
4408 }
4409 }
4410 if ( num_node_edges ) lwfree(node_edges);
4411 }
4412 }
4413
4414 if ( commonnode == -1 )
4415 {
4416 _lwt_release_edges(edges, nedges);
4417 if ( ptr != buf )
4418 {
4419 lwerror("SQL/MM Spatial exception - other edges connected (%s)",
4420 buf);
4421 }
4422 else
4423 {
4424 lwerror("SQL/MM Spatial exception - non-connected edges");
4425 }
4426 return -1;
4427 }
4428
4429 if ( ! lwt_be_checkTopoGeomRemNode(topo, commonnode,
4430 eid1, eid2 ) )
4431 {
4432 _lwt_release_edges(edges, nedges);
4433 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4434 return -1;
4435 }
4436
4437 /* Construct the geometry of the new edge */
4438 switch (caseno)
4439 {
4440 case 1: /* e1.end = e2.start */
4441 pa = ptarray_clone_deep(e1->geom->points);
4442 //pa = ptarray_merge(pa, e2->geom->points);
4443 ptarray_append_ptarray(pa, e2->geom->points, 0);
4444 newedge.start_node = e1->start_node;
4445 newedge.end_node = e2->end_node;
4446 newedge.next_left = e2->next_left;
4447 newedge.next_right = e1->next_right;
4448 e1freenode = 1;
4449 e2freenode = -1;
4450 e2sign = 1;
4451 break;
4452 case 2: /* e1.end = e2.end */
4453 {
4454 POINTARRAY *pa2;
4455 pa2 = ptarray_clone_deep(e2->geom->points);
4456 ptarray_reverse_in_place(pa2);
4457 pa = ptarray_clone_deep(e1->geom->points);
4458 //pa = ptarray_merge(e1->geom->points, pa);
4459 ptarray_append_ptarray(pa, pa2, 0);
4460 ptarray_free(pa2);
4461 newedge.start_node = e1->start_node;
4462 newedge.end_node = e2->start_node;
4463 newedge.next_left = e2->next_right;
4464 newedge.next_right = e1->next_right;
4465 e1freenode = 1;
4466 e2freenode = 1;
4467 e2sign = -1;
4468 break;
4469 }
4470 case 3: /* e1.start = e2.start */
4471 pa = ptarray_clone_deep(e2->geom->points);
4472 ptarray_reverse_in_place(pa);
4473 //pa = ptarray_merge(pa, e1->geom->points);
4474 ptarray_append_ptarray(pa, e1->geom->points, 0);
4475 newedge.end_node = e1->end_node;
4476 newedge.start_node = e2->end_node;
4477 newedge.next_left = e1->next_left;
4478 newedge.next_right = e2->next_left;
4479 e1freenode = -1;
4480 e2freenode = -1;
4481 e2sign = -1;
4482 break;
4483 case 4: /* e1.start = e2.end */
4484 pa = ptarray_clone_deep(e2->geom->points);
4485 //pa = ptarray_merge(pa, e1->geom->points);
4486 ptarray_append_ptarray(pa, e1->geom->points, 0);
4487 newedge.end_node = e1->end_node;
4488 newedge.start_node = e2->start_node;
4489 newedge.next_left = e1->next_left;
4490 newedge.next_right = e2->next_right;
4491 e1freenode = -1;
4492 e2freenode = 1;
4493 e2sign = 1;
4494 break;
4495 default:
4496 pa = NULL;
4497 e1freenode = 0;
4498 e2freenode = 0;
4499 e2sign = 0;
4500 _lwt_release_edges(edges, nedges);
4501 lwerror("Coding error: caseno=%d should never happen", caseno);
4502 break;
4503 }
4504 newedge.geom = lwline_construct(topo->srid, NULL, pa);
4505
4506 if ( modEdge )
4507 {
4508 /* Update data of the first edge */
4509 newedge.edge_id = eid1;
4510 int result = lwt_be_updateEdgesById(topo,
4511 &newedge,
4512 1,
4513 LWT_COL_EDGE_NEXT_LEFT | LWT_COL_EDGE_NEXT_RIGHT | LWT_COL_EDGE_START_NODE |
4514 LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_GEOM);
4515 if (result == -1)
4516 {
4517 lwline_free(newedge.geom);
4518 _lwt_release_edges(edges, nedges);
4519 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4520 return -1;
4521 }
4522 else if (result != 1)
4523 {
4524 lwline_free(newedge.geom);
4525 if ( edges ) _lwt_release_edges(edges, nedges);
4526 lwerror("Unexpected error: %d edges updated when expecting 1", i);
4527 return -1;
4528 }
4529 }
4530 else
4531 {
4532 /* Add new edge */
4533 newedge.edge_id = -1;
4534 newedge.face_left = e1->face_left;
4535 newedge.face_right = e1->face_right;
4536 int result = lwt_be_insertEdges(topo, &newedge, 1);
4537 if (result == -1)
4538 {
4539 lwline_free(newedge.geom);
4540 _lwt_release_edges(edges, nedges);
4541 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4542 return -1;
4543 }
4544 else if (result == 0)
4545 {
4546 lwline_free(newedge.geom);
4547 _lwt_release_edges(edges, nedges);
4548 lwerror("Insertion of split edge failed (no reason)");
4549 return -1;
4550 }
4551 }
4552 lwline_free(newedge.geom);
4553
4554 /*
4555 -- Update next_left_edge/next_right_edge for
4556 -- any edge having them still pointing at the edge being removed
4557 -- (eid2 only when modEdge, or both otherwise)
4558 --
4559 -- NOTE:
4560 -- e#freenode is 1 when edge# end node was the common node
4561 -- and -1 otherwise. This gives the sign of possibly found references
4562 -- to its "free" (non connected to other edge) endnode.
4563 -- e2sign is -1 if edge1 direction is opposite to edge2 direction,
4564 -- or 1 otherwise.
4565 --
4566 */
4567
4568 /* update edges connected to e2's boundary from their end node */
4569 seledge.next_left = e2freenode * eid2;
4570 updedge.next_left = e2freenode * newedge.edge_id * e2sign;
4571 int result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0);
4572 if (result == -1)
4573 {
4574 _lwt_release_edges(edges, nedges);
4575 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4576 return -1;
4577 }
4578
4579 /* update edges connected to e2's boundary from their start node */
4580 seledge.next_right = e2freenode * eid2;
4581 updedge.next_right = e2freenode * newedge.edge_id * e2sign;
4582 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0);
4583 if (result == -1)
4584 {
4585 _lwt_release_edges(edges, nedges);
4586 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4587 return -1;
4588 }
4589
4590 if ( ! modEdge )
4591 {
4592 /* update edges connected to e1's boundary from their end node */
4593 seledge.next_left = e1freenode * eid1;
4594 updedge.next_left = e1freenode * newedge.edge_id;
4595 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0);
4596 if (result == -1)
4597 {
4598 _lwt_release_edges(edges, nedges);
4599 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4600 return -1;
4601 }
4602
4603 /* update edges connected to e1's boundary from their start node */
4604 seledge.next_right = e1freenode * eid1;
4605 updedge.next_right = e1freenode * newedge.edge_id;
4606 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0);
4607 if (result == -1)
4608 {
4609 _lwt_release_edges(edges, nedges);
4610 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4611 return -1;
4612 }
4613 }
4614
4615 /* delete the edges (only second on modEdge or both) */
4616 result = lwt_be_deleteEdges(topo, e2, LWT_COL_EDGE_EDGE_ID);
4617 if (result == -1)
4618 {
4619 _lwt_release_edges(edges, nedges);
4620 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4621 return -1;
4622 }
4623 if ( ! modEdge ) {
4624 i = lwt_be_deleteEdges(topo, e1, LWT_COL_EDGE_EDGE_ID);
4625 if (result == -1)
4626 {
4627 _lwt_release_edges(edges, nedges);
4628 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4629 return -1;
4630 }
4631 }
4632
4633 _lwt_release_edges(edges, nedges);
4634
4635 /* delete the common node */
4636 i = lwt_be_deleteNodesById( topo, &commonnode, 1 );
4637 if (result == -1)
4638 {
4639 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4640 return -1;
4641 }
4642
4643 /*
4644 --
4645 -- NOT IN THE SPECS:
4646 -- Drop composition rows involving second
4647 -- edge, as the first edge took its space,
4648 -- and all affected TopoGeom have been previously checked
4649 -- for being composed by both edges.
4650 */
4651 if ( ! lwt_be_updateTopoGeomEdgeHeal(topo,
4652 eid1, eid2, newedge.edge_id) )
4653 {
4654 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4655 return -1;
4656 }
4657
4658 return modEdge ? commonnode : newedge.edge_id;
4659 }
4660
4661 LWT_ELEMID
lwt_ModEdgeHeal(LWT_TOPOLOGY * topo,LWT_ELEMID e1,LWT_ELEMID e2)4662 lwt_ModEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 )
4663 {
4664 return _lwt_HealEdges( topo, e1, e2, 1 );
4665 }
4666
4667 LWT_ELEMID
lwt_NewEdgeHeal(LWT_TOPOLOGY * topo,LWT_ELEMID e1,LWT_ELEMID e2)4668 lwt_NewEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 )
4669 {
4670 return _lwt_HealEdges( topo, e1, e2, 0 );
4671 }
4672
4673 LWT_ELEMID
lwt_GetNodeByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4674 lwt_GetNodeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4675 {
4676 LWT_ISO_NODE *elem;
4677 uint64_t num;
4678 int flds = LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM; /* geom not needed */
4679 LWT_ELEMID id = 0;
4680 POINT2D qp; /* query point */
4681
4682 if ( ! getPoint2d_p(pt->point, 0, &qp) )
4683 {
4684 lwerror("Empty query point");
4685 return -1;
4686 }
4687 elem = lwt_be_getNodeWithinDistance2D(topo, pt, tol, &num, flds, 0);
4688 if (num == UINT64_MAX)
4689 {
4690 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4691 return -1;
4692 }
4693 else if ( num )
4694 {
4695 if ( num > 1 )
4696 {
4697 _lwt_release_nodes(elem, num);
4698 lwerror("Two or more nodes found");
4699 return -1;
4700 }
4701 id = elem[0].node_id;
4702 _lwt_release_nodes(elem, num);
4703 }
4704
4705 return id;
4706 }
4707
4708 LWT_ELEMID
lwt_GetEdgeByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4709 lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4710 {
4711 LWT_ISO_EDGE *elem;
4712 uint64_t num, i;
4713 int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; /* GEOM is not needed */
4714 LWT_ELEMID id = 0;
4715 LWGEOM *qp = lwpoint_as_lwgeom(pt); /* query point */
4716
4717 if ( lwgeom_is_empty(qp) )
4718 {
4719 lwerror("Empty query point");
4720 return -1;
4721 }
4722 elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol, &num, flds, 0);
4723 if (num == UINT64_MAX)
4724 {
4725 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4726 return -1;
4727 }
4728 for (i=0; i<num;++i)
4729 {
4730 LWT_ISO_EDGE *e = &(elem[i]);
4731 #if 0
4732 LWGEOM* geom;
4733 double dist;
4734
4735 if ( ! e->geom )
4736 {
4737 _lwt_release_edges(elem, num);
4738 lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
4739 " has null geometry", e->edge_id);
4740 continue;
4741 }
4742
4743 /* Should we check for intersection not being on an endpoint
4744 * as documented ? */
4745 geom = lwline_as_lwgeom(e->geom);
4746 dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
4747 if ( dist > tol ) continue;
4748 #endif
4749
4750 if ( id )
4751 {
4752 _lwt_release_edges(elem, num);
4753 lwerror("Two or more edges found");
4754 return -1;
4755 }
4756 else id = e->edge_id;
4757 }
4758
4759 if ( num ) _lwt_release_edges(elem, num);
4760
4761 return id;
4762 }
4763
4764 LWT_ELEMID
lwt_GetFaceByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4765 lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4766 {
4767 LWT_ELEMID id = 0;
4768 LWT_ISO_EDGE *elem;
4769 uint64_t num, i;
4770 int flds = LWT_COL_EDGE_EDGE_ID |
4771 LWT_COL_EDGE_GEOM |
4772 LWT_COL_EDGE_FACE_LEFT |
4773 LWT_COL_EDGE_FACE_RIGHT;
4774 LWGEOM *qp = lwpoint_as_lwgeom(pt);
4775
4776 id = lwt_be_getFaceContainingPoint(topo, pt);
4777 if ( id == -2 ) {
4778 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4779 return -1;
4780 }
4781
4782 if ( id > 0 )
4783 {
4784 return id;
4785 }
4786 id = 0; /* or it'll be -1 for not found */
4787
4788 LWDEBUG(1, "No face properly contains query point,"
4789 " looking for edges");
4790
4791 /* Not in a face, may be in universe or on edge, let's check
4792 * for distance */
4793 /* NOTE: we never pass a tolerance of 0 to avoid ever using
4794 * ST_Within, which doesn't include endpoints matches */
4795 elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol?tol:1e-5, &num, flds, 0);
4796 if (num == UINT64_MAX)
4797 {
4798 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4799 return -1;
4800 }
4801 for (i=0; i<num; ++i)
4802 {
4803 LWT_ISO_EDGE *e = &(elem[i]);
4804 LWT_ELEMID eface = 0;
4805 LWGEOM* geom;
4806 double dist;
4807
4808 if ( ! e->geom )
4809 {
4810 _lwt_release_edges(elem, num);
4811 lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
4812 " has null geometry", e->edge_id);
4813 continue;
4814 }
4815
4816 /* don't consider dangling edges */
4817 if ( e->face_left == e->face_right )
4818 {
4819 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
4820 " is dangling, won't consider it", e->edge_id);
4821 continue;
4822 }
4823
4824 geom = lwline_as_lwgeom(e->geom);
4825 dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
4826
4827 LWDEBUGF(1, "Distance from edge %" LWTFMT_ELEMID
4828 " is %g (tol=%g)", e->edge_id, dist, tol);
4829
4830 /* we won't consider edges too far */
4831 if ( dist > tol ) continue;
4832 if ( e->face_left == 0 ) {
4833 eface = e->face_right;
4834 }
4835 else if ( e->face_right == 0 ) {
4836 eface = e->face_left;
4837 }
4838 else {
4839 _lwt_release_edges(elem, num);
4840 lwerror("Two or more faces found");
4841 return -1;
4842 }
4843
4844 if ( id && id != eface )
4845 {
4846 _lwt_release_edges(elem, num);
4847 lwerror("Two or more faces found"
4848 #if 0 /* debugging */
4849 " (%" LWTFMT_ELEMID
4850 " and %" LWTFMT_ELEMID ")", id, eface
4851 #endif
4852 );
4853 return -1;
4854 }
4855 else id = eface;
4856 }
4857 if ( num ) _lwt_release_edges(elem, num);
4858
4859 return id;
4860 }
4861
4862 /* Return the smallest delta that can perturbate
4863 * the given value */
4864 static inline double
_lwt_minToleranceDouble(double d)4865 _lwt_minToleranceDouble( double d )
4866 {
4867 double ret = 3.6 * pow(10, - ( 15 - log10(d?d:1.0) ) );
4868 return ret;
4869 }
4870
4871 /* Return the smallest delta that can perturbate
4872 * the given point
4873 static inline double
4874 _lwt_minTolerancePoint2d( const POINT2D* p )
4875 {
4876 double max = FP_ABS(p->x);
4877 if ( max < FP_ABS(p->y) ) max = FP_ABS(p->y);
4878 return _lwt_minToleranceDouble(max);
4879 }
4880 */
4881
4882 /* Return the smallest delta that can perturbate
4883 * the maximum absolute value of a geometry ordinate
4884 */
4885 static double
_lwt_minTolerance(LWGEOM * g)4886 _lwt_minTolerance( LWGEOM *g )
4887 {
4888 const GBOX* gbox;
4889 double max;
4890 double ret;
4891
4892 gbox = lwgeom_get_bbox(g);
4893 if ( ! gbox ) return 0; /* empty */
4894 max = FP_ABS(gbox->xmin);
4895 if ( max < FP_ABS(gbox->xmax) ) max = FP_ABS(gbox->xmax);
4896 if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin);
4897 if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax);
4898
4899 ret = _lwt_minToleranceDouble(max);
4900
4901 return ret;
4902 }
4903
4904 #define _LWT_MINTOLERANCE( topo, geom ) ( \
4905 topo->precision ? topo->precision : _lwt_minTolerance(geom) )
4906
4907 typedef struct scored_pointer_t {
4908 void *ptr;
4909 double score;
4910 } scored_pointer;
4911
4912 static int
compare_scored_pointer(const void * si1,const void * si2)4913 compare_scored_pointer(const void *si1, const void *si2)
4914 {
4915 double a = ((scored_pointer *)si1)->score;
4916 double b = ((scored_pointer *)si2)->score;
4917 if ( a < b )
4918 return -1;
4919 else if ( a > b )
4920 return 1;
4921 else
4922 return 0;
4923 }
4924
4925 /*
4926 * @param findFace if non-zero the code will determine which face
4927 * contains the given point (unless it is known to be NOT
4928 * isolated)
4929 * @param moved if not-null will be set to 0 if the point was added
4930 * w/out any snapping or 1 otherwise.
4931 */
4932 static LWT_ELEMID
_lwt_AddPoint(LWT_TOPOLOGY * topo,LWPOINT * point,double tol,int findFace,int * moved)4933 _lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol, int
4934 findFace, int *moved)
4935 {
4936 uint64_t num, i;
4937 double mindist = FLT_MAX;
4938 LWT_ISO_NODE *nodes, *nodes2;
4939 LWT_ISO_EDGE *edges, *edges2;
4940 LWGEOM *pt = lwpoint_as_lwgeom(point);
4941 int flds;
4942 LWT_ELEMID id = 0;
4943 scored_pointer *sorted;
4944
4945 /* Get tolerance, if 0 was given */
4946 if (!tol)
4947 tol = _LWT_MINTOLERANCE(topo, pt);
4948
4949 LWDEBUGG(1, pt, "Adding point");
4950
4951 /*
4952 -- 1. Check if any existing node is closer than the given precision
4953 -- and if so pick the closest
4954 TODO: use WithinBox2D
4955 */
4956 flds = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM;
4957 nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
4958 if (num == UINT64_MAX)
4959 {
4960 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4961 return -1;
4962 }
4963 if ( num )
4964 {
4965 LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
4966 /* Order by distance if there are more than a single return */
4967 if ( num > 1 )
4968 {{
4969 sorted= lwalloc(sizeof(scored_pointer)*num);
4970 for (i=0; i<num; ++i)
4971 {
4972 sorted[i].ptr = nodes+i;
4973 sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
4974 LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
4975 ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
4976 }
4977 qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
4978 nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
4979 for (i=0; i<num; ++i)
4980 {
4981 nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
4982 }
4983 lwfree(sorted);
4984 lwfree(nodes);
4985 nodes = nodes2;
4986 }}
4987
4988 for ( i=0; i<num; ++i )
4989 {
4990 LWT_ISO_NODE *n = &(nodes[i]);
4991 LWGEOM *g = lwpoint_as_lwgeom(n->geom);
4992 double dist = lwgeom_mindistance2d(g, pt);
4993 /* TODO: move this check in the previous sort scan ... */
4994 /* must be closer than tolerated, unless distance is zero */
4995 if ( dist && dist >= tol ) continue;
4996 if ( ! id || dist < mindist )
4997 {
4998 id = n->node_id;
4999 mindist = dist;
5000 }
5001 }
5002 if ( id )
5003 {
5004 /* found an existing node */
5005 if ( nodes ) _lwt_release_nodes(nodes, num);
5006 if ( moved ) *moved = mindist == 0 ? 0 : 1;
5007 return id;
5008 }
5009 }
5010
5011 initGEOS(lwnotice, lwgeom_geos_error);
5012
5013 /*
5014 -- 2. Check if any existing edge falls within tolerance
5015 -- and if so split it by a point projected on it
5016 TODO: use WithinBox2D
5017 */
5018 flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
5019 edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
5020 if (num == UINT64_MAX)
5021 {
5022 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5023 return -1;
5024 }
5025 if ( num )
5026 {
5027 LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
5028
5029 /* Order by distance if there are more than a single return */
5030 if ( num > 1 )
5031 {{
5032 int j;
5033 sorted = lwalloc(sizeof(scored_pointer)*num);
5034 for (i=0; i<num; ++i)
5035 {
5036 sorted[i].ptr = edges+i;
5037 sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
5038 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
5039 ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
5040 }
5041 qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5042 edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
5043 for (j=0, i=0; i<num; ++i)
5044 {
5045 if ( sorted[i].score == sorted[0].score )
5046 {
5047 edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
5048 }
5049 else
5050 {
5051 lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
5052 }
5053 }
5054 num = j;
5055 lwfree(sorted);
5056 lwfree(edges);
5057 edges = edges2;
5058 }}
5059
5060 for (i=0; i<num; ++i)
5061 {
5062 /* The point is on or near an edge, split the edge */
5063 LWT_ISO_EDGE *e = &(edges[i]);
5064 LWGEOM *g = lwline_as_lwgeom(e->geom);
5065 LWGEOM *prj;
5066 int contains;
5067 LWT_ELEMID edge_id = e->edge_id;
5068
5069 LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
5070
5071 /* project point to line, split edge by point */
5072 prj = lwgeom_closest_point(g, pt);
5073 if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1;
5074 if ( lwgeom_has_z(pt) )
5075 {{
5076 /*
5077 -- This is a workaround for ClosestPoint lack of Z support:
5078 -- http://trac.osgeo.org/postgis/ticket/2033
5079 */
5080 LWGEOM *tmp;
5081 double z;
5082 POINT4D p4d;
5083 LWPOINT *prjpt;
5084 /* add Z to "prj" */
5085 tmp = lwgeom_force_3dz(prj, 0);
5086 prjpt = lwgeom_as_lwpoint(tmp);
5087 getPoint4d_p(point->point, 0, &p4d);
5088 z = p4d.z;
5089 getPoint4d_p(prjpt->point, 0, &p4d);
5090 p4d.z = z;
5091 ptarray_set_point4d(prjpt->point, 0, &p4d);
5092 lwgeom_free(prj);
5093 prj = tmp;
5094 }}
5095 const POINT2D *pt = getPoint2d_cp(lwgeom_as_lwpoint(prj)->point, 0);
5096 contains = ptarray_contains_point_partial(e->geom->points, pt, 0, NULL) == LW_BOUNDARY;
5097 if ( ! contains )
5098 {{
5099 double snaptol;
5100 LWGEOM *snapedge;
5101 LWLINE *snapline;
5102 POINT4D p1, p2;
5103
5104 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
5105 " does not contain projected point to it",
5106 edge_id);
5107
5108 /* In order to reduce the robustness issues, we'll pick
5109 * an edge that contains the projected point, if possible */
5110 if ( i+1 < num )
5111 {
5112 LWDEBUG(1, "But there's another to check");
5113 lwgeom_free(prj);
5114 continue;
5115 }
5116
5117 /*
5118 -- The tolerance must be big enough for snapping to happen
5119 -- and small enough to snap only to the projected point.
5120 -- Unfortunately ST_Distance returns 0 because it also uses
5121 -- a projected point internally, so we need another way.
5122 */
5123 snaptol = _lwt_minTolerance(prj);
5124 snapedge = _lwt_toposnap(g, prj, snaptol);
5125 snapline = lwgeom_as_lwline(snapedge);
5126
5127 LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
5128
5129 /* TODO: check if snapping did anything ? */
5130 #if POSTGIS_DEBUG_LEVEL > 0
5131 {
5132 size_t sz;
5133 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5134 char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
5135 LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
5136 lwfree(wkt1);
5137 lwfree(wkt2);
5138 }
5139 #endif
5140
5141
5142 /*
5143 -- Snapping currently snaps the first point below tolerance
5144 -- so may possibly move first point. See ticket #1631
5145 */
5146 getPoint4d_p(e->geom->points, 0, &p1);
5147 getPoint4d_p(snapline->points, 0, &p2);
5148 LWDEBUGF(1, "Edge first point is %g %g, "
5149 "snapline first point is %g %g",
5150 p1.x, p1.y, p2.x, p2.y);
5151 if ( p1.x != p2.x || p1.y != p2.y )
5152 {
5153 LWDEBUG(1, "Snapping moved first point, re-adding it");
5154 if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
5155 {
5156 lwgeom_free(prj);
5157 lwgeom_free(snapedge);
5158 _lwt_release_edges(edges, num);
5159 lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5160 return -1;
5161 }
5162 #if POSTGIS_DEBUG_LEVEL > 0
5163 {
5164 size_t sz;
5165 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5166 LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
5167 lwfree(wkt1);
5168 }
5169 #endif
5170 }
5171 #if POSTGIS_DEBUG_LEVEL > 0
5172 else {
5173 LWDEBUG(1, "Snapping did not move first point");
5174 }
5175 #endif
5176
5177 if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
5178 {
5179 /* TODO: should have invoked lwerror already, leaking memory */
5180 lwgeom_free(prj);
5181 lwgeom_free(snapedge);
5182 _lwt_release_edges(edges, num);
5183 lwerror("lwt_ChangeEdgeGeom failed");
5184 return -1;
5185 }
5186 lwgeom_free(snapedge);
5187 }}
5188 #if POSTGIS_DEBUG_LEVEL > 0
5189 else
5190 {{
5191 size_t sz;
5192 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5193 char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
5194 LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
5195 lwfree(wkt1);
5196 lwfree(wkt2);
5197 }}
5198 #endif
5199
5200 /* TODO: pass 1 as last argument (skipChecks) ? */
5201 id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
5202 if ( -1 == id )
5203 {
5204 /* TODO: should have invoked lwerror already, leaking memory */
5205 lwgeom_free(prj);
5206 _lwt_release_edges(edges, num);
5207 lwerror("lwt_ModEdgeSplit failed");
5208 return -1;
5209 }
5210
5211 lwgeom_free(prj);
5212
5213 /*
5214 * TODO: decimate the two new edges with the given tolerance ?
5215 *
5216 * the edge identifiers to decimate would be: edge_id and "id"
5217 * The problem here is that decimation of existing edges
5218 * may introduce intersections or topological inconsistencies,
5219 * for example:
5220 *
5221 * - A node may end up falling on the other side of the edge
5222 * - The decimated edge might intersect another existing edge
5223 *
5224 */
5225
5226 break; /* we only want to snap a single edge */
5227 }
5228 _lwt_release_edges(edges, num);
5229 }
5230 else
5231 {
5232 /* The point is isolated, add it as such */
5233 /* TODO: pass 1 as last argument (skipChecks) ? */
5234 id = _lwt_AddIsoNode(topo, -1, point, 0, findFace);
5235 if ( moved ) *moved = 0;
5236 if ( -1 == id )
5237 {
5238 /* should have invoked lwerror already, leaking memory */
5239 lwerror("lwt_AddIsoNode failed");
5240 return -1;
5241 }
5242 }
5243
5244 return id;
5245 }
5246
5247 LWT_ELEMID
lwt_AddPoint(LWT_TOPOLOGY * topo,LWPOINT * point,double tol)5248 lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol)
5249 {
5250 return _lwt_AddPoint(topo, point, tol, 1, NULL);
5251 }
5252
5253 /* Return identifier of an equal edge, 0 if none or -1 on error
5254 * (and lwerror gets called on error)
5255 */
5256 static LWT_ELEMID
_lwt_GetEqualEdge(LWT_TOPOLOGY * topo,LWLINE * edge)5257 _lwt_GetEqualEdge( LWT_TOPOLOGY *topo, LWLINE *edge )
5258 {
5259 LWT_ELEMID id;
5260 LWT_ISO_EDGE *edges;
5261 uint64_t num, i;
5262 const GBOX *qbox = lwgeom_get_bbox( lwline_as_lwgeom(edge) );
5263 GEOSGeometry *edgeg;
5264 const int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
5265
5266 edges = lwt_be_getEdgeWithinBox2D( topo, qbox, &num, flds, 0 );
5267 if (num == UINT64_MAX)
5268 {
5269 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5270 return -1;
5271 }
5272 if ( num )
5273 {
5274 initGEOS(lwnotice, lwgeom_geos_error);
5275
5276 edgeg = LWGEOM2GEOS( lwline_as_lwgeom(edge), 0 );
5277 if ( ! edgeg )
5278 {
5279 _lwt_release_edges(edges, num);
5280 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5281 return -1;
5282 }
5283 for (i=0; i<num; ++i)
5284 {
5285 LWT_ISO_EDGE *e = &(edges[i]);
5286 LWGEOM *g = lwline_as_lwgeom(e->geom);
5287 GEOSGeometry *gg;
5288 int equals;
5289 gg = LWGEOM2GEOS( g, 0 );
5290 if ( ! gg )
5291 {
5292 GEOSGeom_destroy(edgeg);
5293 _lwt_release_edges(edges, num);
5294 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5295 return -1;
5296 }
5297 equals = GEOSEquals(gg, edgeg);
5298 GEOSGeom_destroy(gg);
5299 if ( equals == 2 )
5300 {
5301 GEOSGeom_destroy(edgeg);
5302 _lwt_release_edges(edges, num);
5303 lwerror("GEOSEquals exception: %s", lwgeom_geos_errmsg);
5304 return -1;
5305 }
5306 if ( equals )
5307 {
5308 id = e->edge_id;
5309 GEOSGeom_destroy(edgeg);
5310 _lwt_release_edges(edges, num);
5311 return id;
5312 }
5313 }
5314 GEOSGeom_destroy(edgeg);
5315 _lwt_release_edges(edges, num);
5316 }
5317
5318 return 0;
5319 }
5320
5321 /*
5322 * Add a pre-noded pre-split line edge. Used by lwt_AddLine
5323 * Return edge id, 0 if none added (empty edge), -1 on error
5324 *
5325 * @param handleFaceSplit if non-zero the code will check
5326 * if the newly added edge would split a face and if so
5327 * would create new faces accordingly. Otherwise it will
5328 * set left_face and right_face to null (-1)
5329 */
5330 static LWT_ELEMID
_lwt_AddLineEdge(LWT_TOPOLOGY * topo,LWLINE * edge,double tol,int handleFaceSplit)5331 _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol,
5332 int handleFaceSplit )
5333 {
5334 LWCOLLECTION *col;
5335 LWPOINT *start_point, *end_point;
5336 LWGEOM *tmp = 0, *tmp2;
5337 LWT_ISO_NODE *node;
5338 LWT_ELEMID nid[2]; /* start_node, end_node */
5339 LWT_ELEMID id; /* edge id */
5340 POINT4D p4d;
5341 uint64_t nn, i;
5342 int moved=0, mm;
5343
5344 LWDEBUGG(1, lwline_as_lwgeom(edge), "_lwtAddLineEdge");
5345 LWDEBUGF(1, "_lwtAddLineEdge with tolerance %g", tol);
5346
5347 start_point = lwline_get_lwpoint(edge, 0);
5348 if ( ! start_point )
5349 {
5350 lwnotice("Empty component of noded line");
5351 return 0; /* must be empty */
5352 }
5353 nid[0] = _lwt_AddPoint( topo, start_point,
5354 _lwt_minTolerance(lwpoint_as_lwgeom(start_point)),
5355 handleFaceSplit, &mm );
5356 lwpoint_free(start_point); /* too late if lwt_AddPoint calls lwerror */
5357 if ( nid[0] == -1 ) return -1; /* lwerror should have been called */
5358 moved += mm;
5359
5360
5361 end_point = lwline_get_lwpoint(edge, edge->points->npoints-1);
5362 if ( ! end_point )
5363 {
5364 lwerror("could not get last point of line "
5365 "after successfully getting first point !?");
5366 return -1;
5367 }
5368 nid[1] = _lwt_AddPoint( topo, end_point,
5369 _lwt_minTolerance(lwpoint_as_lwgeom(end_point)),
5370 handleFaceSplit, &mm );
5371 moved += mm;
5372 lwpoint_free(end_point); /* too late if lwt_AddPoint calls lwerror */
5373 if ( nid[1] == -1 ) return -1; /* lwerror should have been called */
5374
5375 /*
5376 -- Added endpoints may have drifted due to tolerance, so
5377 -- we need to re-snap the edge to the new nodes before adding it
5378 */
5379 if ( moved )
5380 {
5381
5382 nn = nid[0] == nid[1] ? 1 : 2;
5383 node = lwt_be_getNodeById( topo, nid, &nn,
5384 LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM );
5385 if (nn == UINT64_MAX)
5386 {
5387 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5388 return -1;
5389 }
5390 start_point = NULL; end_point = NULL;
5391 for (i=0; i<nn; ++i)
5392 {
5393 if ( node[i].node_id == nid[0] ) start_point = node[i].geom;
5394 if ( node[i].node_id == nid[1] ) end_point = node[i].geom;
5395 }
5396 if ( ! start_point || ! end_point )
5397 {
5398 if ( nn ) _lwt_release_nodes(node, nn);
5399 lwerror("Could not find just-added nodes % " LWTFMT_ELEMID
5400 " and %" LWTFMT_ELEMID, nid[0], nid[1]);
5401 return -1;
5402 }
5403
5404 /* snap */
5405
5406 getPoint4d_p( start_point->point, 0, &p4d );
5407 lwline_setPoint4d(edge, 0, &p4d);
5408
5409 getPoint4d_p( end_point->point, 0, &p4d );
5410 lwline_setPoint4d(edge, edge->points->npoints-1, &p4d);
5411
5412 if ( nn ) _lwt_release_nodes(node, nn);
5413
5414 /* make valid, after snap (to handle collapses) */
5415 tmp = lwgeom_make_valid(lwline_as_lwgeom(edge));
5416
5417 col = lwgeom_as_lwcollection(tmp);
5418 if ( col )
5419 {{
5420
5421 LWCOLLECTION *colex = lwcollection_extract(col, LINETYPE);
5422
5423 /* Check if the so-snapped edge collapsed (see #1650) */
5424 if ( colex->ngeoms == 0 )
5425 {
5426 lwcollection_free(colex);
5427 lwgeom_free(tmp);
5428 LWDEBUG(1, "Made-valid snapped edge collapsed");
5429 return 0;
5430 }
5431
5432 tmp2 = lwgeom_clone_deep(colex->geoms[0]);
5433 lwgeom_free(tmp);
5434 tmp = tmp2;
5435 edge = lwgeom_as_lwline(tmp);
5436 lwcollection_free(colex);
5437 if ( ! edge )
5438 {
5439 /* should never happen */
5440 lwerror("lwcollection_extract(LINETYPE) returned a non-line?");
5441 return -1;
5442 }
5443 }}
5444 else
5445 {
5446 edge = lwgeom_as_lwline(tmp);
5447 if ( ! edge )
5448 {
5449 LWDEBUGF(1, "Made-valid snapped edge collapsed to %s",
5450 lwtype_name(lwgeom_get_type(tmp)));
5451 lwgeom_free(tmp);
5452 return 0;
5453 }
5454 }
5455 }
5456
5457 /* check if the so-snapped edge _now_ exists */
5458 id = _lwt_GetEqualEdge ( topo, edge );
5459 LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id);
5460 if ( id == -1 )
5461 {
5462 if ( tmp ) lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5463 return -1;
5464 }
5465 if ( id )
5466 {
5467 if ( tmp ) lwgeom_free(tmp); /* possibly takes "edge" down with it */
5468 return id;
5469 }
5470
5471 /* No previously existing edge was found, we'll add one */
5472
5473 /* Remove consecutive vertices below given tolerance
5474 * on edge addition */
5475 if ( tol )
5476 {{
5477 tmp2 = lwline_remove_repeated_points(edge, tol);
5478 LWDEBUGG(1, tmp2, "Repeated-point removed");
5479 edge = lwgeom_as_lwline(tmp2);
5480 if ( tmp ) lwgeom_free(tmp);
5481 tmp = tmp2;
5482
5483 /* check if the so-decimated edge collapsed to single-point */
5484 if ( nid[0] == nid[1] && edge->points->npoints == 2 )
5485 {
5486 lwgeom_free(tmp);
5487 LWDEBUG(1, "Repeated-point removed edge collapsed");
5488 return 0;
5489 }
5490
5491 /* check if the so-decimated edge _now_ exists */
5492 id = _lwt_GetEqualEdge ( topo, edge );
5493 LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id);
5494 if ( id == -1 )
5495 {
5496 lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5497 return -1;
5498 }
5499 if ( id )
5500 {
5501 lwgeom_free(tmp); /* takes "edge" down with it */
5502 return id;
5503 }
5504 }}
5505
5506
5507 /* TODO: skip checks ? */
5508 id = _lwt_AddEdge( topo, nid[0], nid[1], edge, 0, handleFaceSplit ? 1 : -1 );
5509 LWDEBUGF(1, "lwt_AddEdgeModFace returned %" LWTFMT_ELEMID, id);
5510 if ( id == -1 )
5511 {
5512 lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5513 return -1;
5514 }
5515 lwgeom_free(tmp); /* possibly takes "edge" down with it */
5516
5517 return id;
5518 }
5519
5520 /* Simulate split-loop as it was implemented in pl/pgsql version
5521 * of TopoGeo_addLinestring */
5522 static LWGEOM *
_lwt_split_by_nodes(const LWGEOM * g,const LWGEOM * nodes)5523 _lwt_split_by_nodes(const LWGEOM *g, const LWGEOM *nodes)
5524 {
5525 LWCOLLECTION *col = lwgeom_as_lwcollection(nodes);
5526 uint32_t i;
5527 LWGEOM *bg;
5528
5529 bg = lwgeom_clone_deep(g);
5530 if ( ! col->ngeoms ) return bg;
5531
5532 for (i=0; i<col->ngeoms; ++i)
5533 {
5534 LWGEOM *g2;
5535 g2 = lwgeom_split(bg, col->geoms[i]);
5536 lwgeom_free(bg);
5537 bg = g2;
5538 }
5539 bg->srid = nodes->srid;
5540
5541 return bg;
5542 }
5543
5544 static LWT_ELEMID*
_lwt_AddLine(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges,int handleFaceSplit)5545 _lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges,
5546 int handleFaceSplit)
5547 {
5548 LWGEOM *geomsbuf[1];
5549 LWGEOM **geoms;
5550 uint32_t ngeoms;
5551 LWGEOM *noded, *tmp;
5552 LWCOLLECTION *col;
5553 LWT_ELEMID *ids;
5554 LWT_ISO_EDGE *edges;
5555 LWT_ISO_NODE *nodes;
5556 uint64_t num, numedges = 0, numnodes = 0;
5557 uint64_t i;
5558 GBOX qbox;
5559
5560 *nedges = -1; /* error condition, by default */
5561
5562 /* Get tolerance, if 0 was given */
5563 if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)line );
5564 LWDEBUGF(1, "Working tolerance:%.15g", tol);
5565 LWDEBUGF(1, "Input line has srid=%d", line->srid);
5566
5567 /* Remove consecutive vertices below given tolerance upfront */
5568 if ( tol )
5569 {{
5570 LWLINE *clean = lwgeom_as_lwline(lwline_remove_repeated_points(line, tol));
5571 tmp = lwline_as_lwgeom(clean); /* NOTE: might collapse to non-simple */
5572 LWDEBUGG(1, tmp, "Repeated-point removed");
5573 }} else tmp=(LWGEOM*)line;
5574
5575 /* 1. Self-node */
5576 noded = lwgeom_node((LWGEOM*)tmp);
5577 if ( tmp != (LWGEOM*)line ) lwgeom_free(tmp);
5578 if ( ! noded ) return NULL; /* should have called lwerror already */
5579 LWDEBUGG(1, noded, "Noded");
5580
5581 qbox = *lwgeom_get_bbox( lwline_as_lwgeom(line) );
5582 LWDEBUGF(1, "Line BOX is %.15g %.15g, %.15g %.15g", qbox.xmin, qbox.ymin,
5583 qbox.xmax, qbox.ymax);
5584 gbox_expand(&qbox, tol);
5585 LWDEBUGF(1, "BOX expanded by %g is %.15g %.15g, %.15g %.15g",
5586 tol, qbox.xmin, qbox.ymin, qbox.xmax, qbox.ymax);
5587
5588 LWGEOM **nearby = 0;
5589 int nearbyindex = 0;
5590 int nearbycount = 0;
5591
5592 /* 2.0. Find edges falling within tol distance */
5593 edges = lwt_be_getEdgeWithinBox2D( topo, &qbox, &numedges, LWT_COL_EDGE_ALL, 0 );
5594 if (numedges == UINT64_MAX)
5595 {
5596 lwgeom_free(noded);
5597 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5598 return NULL;
5599 }
5600 LWDEBUGF(1, "Line has %d points, its bbox intersects %d edges bboxes",
5601 line->points->npoints, numedges);
5602 if ( numedges )
5603 {{
5604 /* collect those whose distance from us is < tol */
5605 nearbycount += numedges;
5606 nearby = lwalloc(numedges * sizeof(LWGEOM *));
5607 for (i=0; i<numedges; ++i)
5608 {
5609 LW_ON_INTERRUPT(return NULL);
5610 LWT_ISO_EDGE *e = &(edges[i]);
5611 LWGEOM *g = lwline_as_lwgeom(e->geom);
5612 LWDEBUGF(2, "Computing distance from edge %d having %d points", i, e->geom->points->npoints);
5613 double dist = lwgeom_mindistance2d(g, noded);
5614 /* must be closer than tolerated, unless distance is zero */
5615 if ( dist && dist >= tol ) continue;
5616 nearby[nearbyindex++] = g;
5617 }
5618 LWDEBUGF(1, "Found %d edges closer than tolerance (%g)", nearbyindex, tol);
5619 }}
5620 int nearbyedgecount = nearbyindex;
5621
5622 /* 2.1. Find isolated nodes falling within tol distance
5623 *
5624 * TODO: add backend-interface support for only getting isolated nodes
5625 */
5626 nodes = lwt_be_getNodeWithinBox2D( topo, &qbox, &numnodes, LWT_COL_NODE_ALL, 0 );
5627 if (numnodes == UINT64_MAX)
5628 {
5629 lwgeom_free(noded);
5630 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5631 return NULL;
5632 }
5633 LWDEBUGF(1, "Line bbox intersects %d nodes bboxes", numnodes);
5634 if ( numnodes )
5635 {{
5636 /* collect those whose distance from us is < tol */
5637 nearbycount = nearbyedgecount + numnodes;
5638 nearby = nearby ?
5639 lwrealloc(nearby, nearbycount * sizeof(LWGEOM *))
5640 :
5641 lwalloc(nearbycount * sizeof(LWGEOM *))
5642 ;
5643 int nn = 0;
5644 for (i=0; i<numnodes; ++i)
5645 {
5646 LWT_ISO_NODE *n = &(nodes[i]);
5647 if ( n->containing_face == -1 ) continue; /* skip not-isolated nodes */
5648 LWGEOM *g = lwpoint_as_lwgeom(n->geom);
5649 double dist = lwgeom_mindistance2d(g, noded);
5650 /* must be closer than tolerated, unless distance is zero */
5651 if ( dist && dist >= tol )
5652 {
5653 LWDEBUGF(1, "Node %d is %g units away, we tolerate only %g", n->node_id, dist, tol);
5654 continue;
5655 }
5656 nearby[nearbyindex++] = g;
5657 ++nn;
5658 }
5659 LWDEBUGF(1, "Found %d isolated nodes closer than tolerance (%g)", nn, tol);
5660 }}
5661 int nearbynodecount = nearbyindex - nearbyedgecount;
5662 nearbycount = nearbyindex;
5663
5664 LWDEBUGF(1, "Number of nearby elements is %d", nearbycount);
5665
5666 /* 2.2. Snap to nearby elements */
5667 if ( nearbycount )
5668 {{
5669 LWCOLLECTION *col;
5670 LWGEOM *elems;
5671
5672 col = lwcollection_construct(COLLECTIONTYPE, topo->srid,
5673 NULL, nearbycount, nearby);
5674 elems = lwcollection_as_lwgeom(col);
5675
5676 LWDEBUGG(1, elems, "Collected nearby elements");
5677
5678 tmp = _lwt_toposnap(noded, elems, tol);
5679 lwgeom_free(noded);
5680 noded = tmp;
5681 LWDEBUGG(1, noded, "Elements-snapped");
5682
5683 /* will not release the geoms array */
5684 lwcollection_release(col);
5685
5686 /*
5687 -- re-node to account for ST_Snap introduced self-intersections
5688 -- See http://trac.osgeo.org/postgis/ticket/1714
5689 -- TODO: consider running UnaryUnion once after all noding
5690 */
5691 tmp = lwgeom_unaryunion(noded);
5692 lwgeom_free(noded);
5693 noded = tmp;
5694 LWDEBUGG(1, noded, "Unary-unioned");
5695 }}
5696
5697 /* 2.3. Node with nearby edges */
5698 if ( nearbyedgecount )
5699 {{
5700 LWCOLLECTION *col;
5701 LWGEOM *iedges; /* just an alias for col */
5702 LWGEOM *diff, *xset;
5703
5704 LWDEBUGF(1, "Line intersects %d edges", nearbyedgecount);
5705
5706 col = lwcollection_construct(COLLECTIONTYPE, topo->srid,
5707 NULL, nearbyedgecount, nearby);
5708 iedges = lwcollection_as_lwgeom(col);
5709 LWDEBUGG(1, iedges, "Collected edges");
5710
5711 LWDEBUGF(1, "Diffing noded, with srid=%d "
5712 "and interesecting edges, with srid=%d",
5713 noded->srid, iedges->srid);
5714 diff = lwgeom_difference(noded, iedges);
5715 LWDEBUGG(1, diff, "Differenced");
5716
5717 LWDEBUGF(1, "Intersecting noded, with srid=%d "
5718 "and interesecting edges, with srid=%d",
5719 noded->srid, iedges->srid);
5720 xset = lwgeom_intersection(noded, iedges);
5721 LWDEBUGG(1, xset, "Intersected");
5722 lwgeom_free(noded);
5723
5724 /* We linemerge here because INTERSECTION, as of GEOS 3.8,
5725 * will result in shared segments being output as multiple
5726 * lines rather than a single line. Example:
5727
5728 INTERSECTION(
5729 'LINESTRING(0 0, 5 0, 8 0, 10 0,12 0)',
5730 'LINESTRING(5 0, 8 0, 10 0)'
5731 )
5732 ==
5733 MULTILINESTRING((5 0,8 0),(8 0,10 0))
5734
5735 * We will re-split in a subsequent step, by splitting
5736 * the final line with pre-existing nodes
5737 */
5738 LWDEBUG(1, "Linemerging intersection");
5739 tmp = lwgeom_linemerge(xset);
5740 LWDEBUGG(1, tmp, "Linemerged");
5741 lwgeom_free(xset);
5742 xset = tmp;
5743
5744 /*
5745 * Here we union the (linemerged) intersection with
5746 * the difference (new lines)
5747 */
5748 LWDEBUG(1, "Unioning difference and (linemerged) intersection");
5749 noded = lwgeom_union(diff, xset);
5750 LWDEBUGG(1, noded, "Diff-Xset Unioned");
5751 lwgeom_free(xset);
5752 lwgeom_free(diff);
5753
5754 /* will not release the geoms array */
5755 lwcollection_release(col);
5756 }}
5757
5758
5759 /* 2.4. Split by pre-existing nodes
5760 *
5761 * Pre-existing nodes are isolated nodes AND endpoints
5762 * of intersecting edges
5763 */
5764 if ( nearbyedgecount )
5765 {
5766 nearbycount += nearbyedgecount * 2; /* make space for endpoints */
5767 nearby = lwrealloc(nearby, nearbycount * sizeof(LWGEOM *));
5768 for (int i=0; i<nearbyedgecount; i++)
5769 {
5770 LWLINE *edge = lwgeom_as_lwline(nearby[i]);
5771 LWPOINT *startNode = lwline_get_lwpoint(edge, 0);
5772 LWPOINT *endNode = lwline_get_lwpoint(edge, edge->points->npoints-1);
5773 /* TODO: only add if within distance to noded AND if not duplicated */
5774 nearby[nearbyindex++] = lwpoint_as_lwgeom(startNode);
5775 nearbynodecount++;
5776 nearby[nearbyindex++] = lwpoint_as_lwgeom(endNode);
5777 nearbynodecount++;
5778 }
5779 }
5780 if ( nearbynodecount )
5781 {
5782 col = lwcollection_construct(MULTIPOINTTYPE, topo->srid,
5783 NULL, nearbynodecount,
5784 nearby + nearbyedgecount);
5785 LWGEOM *inodes = lwcollection_as_lwgeom(col);
5786 /* TODO: use lwgeom_split of lwgeom_union ... */
5787 tmp = _lwt_split_by_nodes(noded, inodes);
5788 lwgeom_free(noded);
5789 noded = tmp;
5790 LWDEBUGG(1, noded, "Node-split");
5791 /* will not release the geoms array */
5792 lwcollection_release(col);
5793 }
5794
5795
5796 LWDEBUG(1, "Freeing up nearby elements");
5797
5798 /* TODO: free up endpoints of nearbyedges */
5799 if ( nearby ) lwfree(nearby);
5800 if ( nodes ) _lwt_release_nodes(nodes, numnodes);
5801 if ( edges ) _lwt_release_edges(edges, numedges);
5802
5803 LWDEBUGG(1, noded, "Finally-noded");
5804
5805 /* 3. For each (now-noded) segment, insert an edge */
5806 col = lwgeom_as_lwcollection(noded);
5807 if ( col )
5808 {
5809 LWDEBUG(1, "Noded line was a collection");
5810 geoms = col->geoms;
5811 ngeoms = col->ngeoms;
5812 }
5813 else
5814 {
5815 LWDEBUG(1, "Noded line was a single geom");
5816 geomsbuf[0] = noded;
5817 geoms = geomsbuf;
5818 ngeoms = 1;
5819 }
5820
5821 LWDEBUGF(1, "Line was split into %d edges", ngeoms);
5822
5823 /* TODO: refactor to first add all nodes (re-snapping edges if
5824 * needed) and then check all edges for existing already
5825 * ( so to save a DB scan for each edge to be added )
5826 */
5827 ids = lwalloc(sizeof(LWT_ELEMID)*ngeoms);
5828 num = 0;
5829 for ( i=0; i<ngeoms; ++i )
5830 {
5831 LWT_ELEMID id;
5832 LWGEOM *g = geoms[i];
5833 g->srid = noded->srid;
5834
5835 #if POSTGIS_DEBUG_LEVEL > 0
5836 {
5837 size_t sz;
5838 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5839 LWDEBUGF(1, "Component %d of split line is: %s", i, wkt1);
5840 lwfree(wkt1);
5841 }
5842 #endif
5843
5844 id = _lwt_AddLineEdge( topo, lwgeom_as_lwline(g), tol, handleFaceSplit );
5845 LWDEBUGF(1, "_lwt_AddLineEdge returned %" LWTFMT_ELEMID, id);
5846 if ( id < 0 )
5847 {
5848 lwgeom_free(noded);
5849 lwfree(ids);
5850 return NULL;
5851 }
5852 if ( ! id )
5853 {
5854 LWDEBUGF(1, "Component %d of split line collapsed", i);
5855 continue;
5856 }
5857
5858 LWDEBUGF(1, "Component %d of split line is edge %" LWTFMT_ELEMID,
5859 i, id);
5860 ids[num++] = id; /* TODO: skip duplicates */
5861 }
5862
5863 LWDEBUGG(1, noded, "Noded before free");
5864 lwgeom_free(noded);
5865
5866 /* TODO: XXX remove duplicated ids if not done before */
5867
5868 *nedges = num;
5869 return ids;
5870 }
5871
5872 LWT_ELEMID*
lwt_AddLine(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges)5873 lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges)
5874 {
5875 return _lwt_AddLine(topo, line, tol, nedges, 1);
5876 }
5877
5878 LWT_ELEMID*
lwt_AddLineNoFace(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges)5879 lwt_AddLineNoFace(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges)
5880 {
5881 return _lwt_AddLine(topo, line, tol, nedges, 0);
5882 }
5883
5884 LWT_ELEMID*
lwt_AddPolygon(LWT_TOPOLOGY * topo,LWPOLY * poly,double tol,int * nfaces)5885 lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* poly, double tol, int* nfaces)
5886 {
5887 uint32_t i;
5888 *nfaces = -1; /* error condition, by default */
5889 int num;
5890 LWT_ISO_FACE *faces;
5891 uint64_t nfacesinbox;
5892 uint64_t j;
5893 LWT_ELEMID *ids = NULL;
5894 GBOX qbox;
5895 const GEOSPreparedGeometry *ppoly;
5896 GEOSGeometry *polyg;
5897
5898 /* Get tolerance, if 0 was given */
5899 if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)poly );
5900 LWDEBUGF(1, "Working tolerance:%.15g", tol);
5901
5902 /* Add each ring as an edge */
5903 for ( i=0; i<poly->nrings; ++i )
5904 {
5905 LWLINE *line;
5906 POINTARRAY *pa;
5907 LWT_ELEMID *eids;
5908 int nedges;
5909
5910 pa = ptarray_clone(poly->rings[i]);
5911 line = lwline_construct(topo->srid, NULL, pa);
5912 eids = lwt_AddLine( topo, line, tol, &nedges );
5913 if ( nedges < 0 ) {
5914 /* probably too late as lwt_AddLine invoked lwerror */
5915 lwline_free(line);
5916 lwerror("Error adding ring %d of polygon", i);
5917 return NULL;
5918 }
5919 lwline_free(line);
5920 lwfree(eids);
5921 }
5922
5923 /*
5924 -- Find faces covered by input polygon
5925 -- NOTE: potential snapping changed polygon edges
5926 */
5927 qbox = *lwgeom_get_bbox( lwpoly_as_lwgeom(poly) );
5928 gbox_expand(&qbox, tol);
5929 faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nfacesinbox,
5930 LWT_COL_FACE_ALL, 0 );
5931 if (nfacesinbox == UINT64_MAX)
5932 {
5933 lwfree(ids);
5934 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5935 return NULL;
5936 }
5937
5938 num = 0;
5939 if ( nfacesinbox )
5940 {
5941 polyg = LWGEOM2GEOS(lwpoly_as_lwgeom(poly), 0);
5942 if ( ! polyg )
5943 {
5944 _lwt_release_faces(faces, nfacesinbox);
5945 lwerror("Could not convert poly geometry to GEOS: %s", lwgeom_geos_errmsg);
5946 return NULL;
5947 }
5948 ppoly = GEOSPrepare(polyg);
5949 ids = lwalloc(sizeof(LWT_ELEMID)*nfacesinbox);
5950 for ( j=0; j<nfacesinbox; ++j )
5951 {
5952 LWT_ISO_FACE *f = &(faces[j]);
5953 LWGEOM *fg;
5954 GEOSGeometry *fgg, *sp;
5955 int covers;
5956
5957 /* check if a point on this face surface is covered by our polygon */
5958 fg = lwt_GetFaceGeometry( topo, f->face_id );
5959 if ( ! fg )
5960 {
5961 j = f->face_id; /* so we can destroy faces */
5962 GEOSPreparedGeom_destroy(ppoly);
5963 GEOSGeom_destroy(polyg);
5964 lwfree(ids);
5965 _lwt_release_faces(faces, nfacesinbox);
5966 lwerror("Could not get geometry of face %" LWTFMT_ELEMID, j);
5967 return NULL;
5968 }
5969 /* check if a point on this face's surface is covered by our polygon */
5970 fgg = LWGEOM2GEOS(fg, 0);
5971 lwgeom_free(fg);
5972 if ( ! fgg )
5973 {
5974 GEOSPreparedGeom_destroy(ppoly);
5975 GEOSGeom_destroy(polyg);
5976 _lwt_release_faces(faces, nfacesinbox);
5977 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5978 return NULL;
5979 }
5980 sp = GEOSPointOnSurface(fgg);
5981 GEOSGeom_destroy(fgg);
5982 if ( ! sp )
5983 {
5984 GEOSPreparedGeom_destroy(ppoly);
5985 GEOSGeom_destroy(polyg);
5986 _lwt_release_faces(faces, nfacesinbox);
5987 lwerror("Could not find point on face surface: %s", lwgeom_geos_errmsg);
5988 return NULL;
5989 }
5990 covers = GEOSPreparedCovers( ppoly, sp );
5991 GEOSGeom_destroy(sp);
5992 if (covers == 2)
5993 {
5994 GEOSPreparedGeom_destroy(ppoly);
5995 GEOSGeom_destroy(polyg);
5996 _lwt_release_faces(faces, nfacesinbox);
5997 lwerror("PreparedCovers error: %s", lwgeom_geos_errmsg);
5998 return NULL;
5999 }
6000 if ( ! covers )
6001 {
6002 continue; /* we're not composed by this face */
6003 }
6004
6005 /* TODO: avoid duplicates ? */
6006 ids[num++] = f->face_id;
6007 }
6008 GEOSPreparedGeom_destroy(ppoly);
6009 GEOSGeom_destroy(polyg);
6010 _lwt_release_faces(faces, nfacesinbox);
6011 }
6012
6013 /* possibly 0 if non face's surface point was found
6014 * to be covered by input polygon */
6015 *nfaces = num;
6016
6017 return ids;
6018 }
6019
6020 /*
6021 *---- polygonizer
6022 */
6023
6024 /* An array of pointers to EDGERING structures */
6025 typedef struct LWT_ISO_EDGE_TABLE_T {
6026 LWT_ISO_EDGE *edges;
6027 int size;
6028 } LWT_ISO_EDGE_TABLE;
6029
6030 static int
compare_iso_edges_by_id(const void * si1,const void * si2)6031 compare_iso_edges_by_id(const void *si1, const void *si2)
6032 {
6033 int a = ((LWT_ISO_EDGE *)si1)->edge_id;
6034 int b = ((LWT_ISO_EDGE *)si2)->edge_id;
6035 if ( a < b )
6036 return -1;
6037 else if ( a > b )
6038 return 1;
6039 else
6040 return 0;
6041 }
6042
6043 static LWT_ISO_EDGE *
_lwt_getIsoEdgeById(LWT_ISO_EDGE_TABLE * tab,LWT_ELEMID id)6044 _lwt_getIsoEdgeById(LWT_ISO_EDGE_TABLE *tab, LWT_ELEMID id)
6045 {
6046 LWT_ISO_EDGE key;
6047 key.edge_id = id;
6048
6049 void *match = bsearch( &key, tab->edges, tab->size,
6050 sizeof(LWT_ISO_EDGE),
6051 compare_iso_edges_by_id);
6052 return match;
6053 }
6054
6055 typedef struct LWT_EDGERING_ELEM_T {
6056 /* externally owned */
6057 LWT_ISO_EDGE *edge;
6058 /* 0 if false, 1 if true */
6059 int left;
6060 } LWT_EDGERING_ELEM;
6061
6062 /* A ring of edges */
6063 typedef struct LWT_EDGERING_T {
6064 /* Signed edge identifiers
6065 * positive ones are walked in their direction, negative ones
6066 * in the opposite direction */
6067 LWT_EDGERING_ELEM **elems;
6068 /* Number of edges in the ring */
6069 int size;
6070 int capacity;
6071 /* Bounding box of the ring */
6072 GBOX *env;
6073 /* Bounding box of the ring in GEOS format (for STRTree) */
6074 GEOSGeometry *genv;
6075 } LWT_EDGERING;
6076
6077 #define LWT_EDGERING_INIT(a) { \
6078 (a)->size = 0; \
6079 (a)->capacity = 1; \
6080 (a)->elems = lwalloc(sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \
6081 (a)->env = NULL; \
6082 (a)->genv = NULL; \
6083 }
6084
6085 #define LWT_EDGERING_PUSH(a, r) { \
6086 if ( (a)->size + 1 > (a)->capacity ) { \
6087 (a)->capacity *= 2; \
6088 (a)->elems = lwrealloc((a)->elems, sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \
6089 } \
6090 /* lwdebug(1, "adding elem %d (%p) of edgering %p", (a)->size, (r), (a)); */ \
6091 (a)->elems[(a)->size++] = (r); \
6092 }
6093
6094 #define LWT_EDGERING_CLEAN(a) { \
6095 int i; for (i=0; i<(a)->size; ++i) { \
6096 if ( (a)->elems[i] ) { \
6097 /* lwdebug(1, "freeing elem %d (%p) of edgering %p", i, (a)->elems[i], (a)); */ \
6098 lwfree((a)->elems[i]); \
6099 } \
6100 } \
6101 if ( (a)->elems ) { lwfree((a)->elems); (a)->elems = NULL; } \
6102 (a)->size = 0; \
6103 (a)->capacity = 0; \
6104 if ( (a)->env ) { lwfree((a)->env); (a)->env = NULL; } \
6105 if ( (a)->genv ) { GEOSGeom_destroy((a)->genv); (a)->genv = NULL; } \
6106 }
6107
6108 /* An array of pointers to EDGERING structures */
6109 typedef struct LWT_EDGERING_ARRAY_T {
6110 LWT_EDGERING **rings;
6111 int size;
6112 int capacity;
6113 GEOSSTRtree* tree;
6114 } LWT_EDGERING_ARRAY;
6115
6116 #define LWT_EDGERING_ARRAY_INIT(a) { \
6117 (a)->size = 0; \
6118 (a)->capacity = 1; \
6119 (a)->rings = lwalloc(sizeof(LWT_EDGERING *) * (a)->capacity); \
6120 (a)->tree = NULL; \
6121 }
6122
6123 /* WARNING: use of 'j' is intentional, not to clash with
6124 * 'i' used in LWT_EDGERING_CLEAN */
6125 #define LWT_EDGERING_ARRAY_CLEAN(a) { \
6126 int j; for (j=0; j<(a)->size; ++j) { \
6127 LWT_EDGERING_CLEAN((a)->rings[j]); \
6128 } \
6129 if ( (a)->capacity ) lwfree((a)->rings); \
6130 if ( (a)->tree ) { \
6131 GEOSSTRtree_destroy( (a)->tree ); \
6132 (a)->tree = NULL; \
6133 } \
6134 }
6135
6136 #define LWT_EDGERING_ARRAY_PUSH(a, r) { \
6137 if ( (a)->size + 1 > (a)->capacity ) { \
6138 (a)->capacity *= 2; \
6139 (a)->rings = lwrealloc((a)->rings, sizeof(LWT_EDGERING *) * (a)->capacity); \
6140 } \
6141 (a)->rings[(a)->size++] = (r); \
6142 }
6143
6144 typedef struct LWT_EDGERING_POINT_ITERATOR_T {
6145 LWT_EDGERING *ring;
6146 LWT_EDGERING_ELEM *curelem;
6147 int curelemidx;
6148 int curidx;
6149 } LWT_EDGERING_POINT_ITERATOR;
6150
6151 static int
_lwt_EdgeRingIterator_next(LWT_EDGERING_POINT_ITERATOR * it,POINT2D * pt)6152 _lwt_EdgeRingIterator_next(LWT_EDGERING_POINT_ITERATOR *it, POINT2D *pt)
6153 {
6154 LWT_EDGERING_ELEM *el = it->curelem;
6155 POINTARRAY *pa;
6156
6157 if ( ! el ) return 0; /* finished */
6158
6159 pa = el->edge->geom->points;
6160
6161 int tonext = 0;
6162 LWDEBUGF(3, "iterator fetching idx %d from pa of %d points", it->curidx, pa->npoints);
6163 getPoint2d_p(pa, it->curidx, pt);
6164 if ( el->left ) {
6165 it->curidx++;
6166 if ( it->curidx >= (int) pa->npoints ) tonext = 1;
6167 } else {
6168 it->curidx--;
6169 if ( it->curidx < 0 ) tonext = 1;
6170 }
6171
6172 if ( tonext )
6173 {
6174 LWDEBUG(3, "iterator moving to next element");
6175 it->curelemidx++;
6176 if ( it->curelemidx < it->ring->size )
6177 {
6178 el = it->curelem = it->ring->elems[it->curelemidx];
6179 it->curidx = el->left ? 0 : el->edge->geom->points->npoints - 1;
6180 }
6181 else
6182 {
6183 it->curelem = NULL;
6184 }
6185 }
6186
6187 return 1;
6188 }
6189
6190 /* Release return with lwfree */
6191 static LWT_EDGERING_POINT_ITERATOR *
_lwt_EdgeRingIterator_begin(LWT_EDGERING * er)6192 _lwt_EdgeRingIterator_begin(LWT_EDGERING *er)
6193 {
6194 LWT_EDGERING_POINT_ITERATOR *ret = lwalloc(sizeof(LWT_EDGERING_POINT_ITERATOR));
6195 ret->ring = er;
6196 if ( er->size ) ret->curelem = er->elems[0];
6197 else ret->curelem = NULL;
6198 ret->curelemidx = 0;
6199 ret->curidx = ret->curelem->left ? 0 : ret->curelem->edge->geom->points->npoints - 1;
6200 return ret;
6201 }
6202
6203 /* Identifier for a placeholder face that will be
6204 * used to mark hole rings */
6205 #define LWT_HOLES_FACE_PLACEHOLDER INT32_MIN
6206
6207 static int
_lwt_FetchNextUnvisitedEdge(LWT_TOPOLOGY * topo,LWT_ISO_EDGE_TABLE * etab,int from)6208 _lwt_FetchNextUnvisitedEdge(__attribute__((__unused__)) LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *etab, int from)
6209 {
6210 while (
6211 from < etab->size &&
6212 etab->edges[from].face_left != -1 &&
6213 etab->edges[from].face_right != -1
6214 ) from++;
6215 return from < etab->size ? from : -1;
6216 }
6217
6218 static LWT_ISO_EDGE *
_lwt_FetchAllEdges(LWT_TOPOLOGY * topo,int * numedges)6219 _lwt_FetchAllEdges(LWT_TOPOLOGY *topo, int *numedges)
6220 {
6221 LWT_ISO_EDGE *edge;
6222 int fields = LWT_COL_EDGE_ALL;
6223 uint64_t nelems = 1;
6224
6225 edge = lwt_be_getEdgeWithinBox2D( topo, NULL, &nelems, fields, 0);
6226 *numedges = nelems;
6227 if (nelems == UINT64_MAX)
6228 {
6229 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6230 return NULL;
6231 }
6232 return edge;
6233 }
6234
6235 /* Update the side face of given ring edges
6236 *
6237 * Edge identifiers are signed, those with negative identifier
6238 * need to be updated their right_face, those with positive
6239 * identifier need to be updated their left_face.
6240 *
6241 * @param face identifier of the face bound by the ring
6242 * @return 0 on success, -1 on error
6243 */
6244 static int
_lwt_UpdateEdgeRingSideFace(LWT_TOPOLOGY * topo,LWT_EDGERING * ring,LWT_ELEMID face)6245 _lwt_UpdateEdgeRingSideFace(LWT_TOPOLOGY *topo, LWT_EDGERING *ring,
6246 LWT_ELEMID face)
6247 {
6248 LWT_ISO_EDGE *forward_edges = NULL;
6249 int forward_edges_count = 0;
6250 LWT_ISO_EDGE *backward_edges = NULL;
6251 int backward_edges_count = 0;
6252 int i, ret;
6253
6254 /* Make a list of forward_edges and backward_edges */
6255
6256 forward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size);
6257 forward_edges_count = 0;
6258 backward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size);
6259 backward_edges_count = 0;
6260
6261 for ( i=0; i<ring->size; ++i )
6262 {
6263 LWT_EDGERING_ELEM *elem = ring->elems[i];
6264 LWT_ISO_EDGE *edge = elem->edge;
6265 LWT_ELEMID id = edge->edge_id;
6266 if ( elem->left )
6267 {
6268 LWDEBUGF(3, "Forward edge %d is %d", forward_edges_count, id);
6269 forward_edges[forward_edges_count].edge_id = id;
6270 forward_edges[forward_edges_count++].face_left = face;
6271 edge->face_left = face;
6272 }
6273 else
6274 {
6275 LWDEBUGF(3, "Backward edge %d is %d", forward_edges_count, id);
6276 backward_edges[backward_edges_count].edge_id = id;
6277 backward_edges[backward_edges_count++].face_right = face;
6278 edge->face_right = face;
6279 }
6280 }
6281
6282 /* Update forward edges */
6283 if ( forward_edges_count )
6284 {
6285 ret = lwt_be_updateEdgesById(topo, forward_edges,
6286 forward_edges_count,
6287 LWT_COL_EDGE_FACE_LEFT);
6288 if ( ret == -1 )
6289 {
6290 lwfree( forward_edges );
6291 lwfree( backward_edges );
6292 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6293 return -1;
6294 }
6295 if ( ret != forward_edges_count )
6296 {
6297 lwfree( forward_edges );
6298 lwfree( backward_edges );
6299 lwerror("Unexpected error: %d edges updated when expecting %d (forward)",
6300 ret, forward_edges_count);
6301 return -1;
6302 }
6303 }
6304
6305 /* Update backward edges */
6306 if ( backward_edges_count )
6307 {
6308 ret = lwt_be_updateEdgesById(topo, backward_edges,
6309 backward_edges_count,
6310 LWT_COL_EDGE_FACE_RIGHT);
6311 if ( ret == -1 )
6312 {
6313 lwfree( forward_edges );
6314 lwfree( backward_edges );
6315 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6316 return -1;
6317 }
6318 if ( ret != backward_edges_count )
6319 {
6320 lwfree( forward_edges );
6321 lwfree( backward_edges );
6322 lwerror("Unexpected error: %d edges updated when expecting %d (backward)",
6323 ret, backward_edges_count);
6324 return -1;
6325 }
6326 }
6327
6328 lwfree( forward_edges );
6329 lwfree( backward_edges );
6330
6331 return 0;
6332 }
6333
6334 /*
6335 * @param side 1 for left side, -1 for right side
6336 */
6337 static LWT_EDGERING *
_lwt_BuildEdgeRing(LWT_TOPOLOGY * topo,LWT_ISO_EDGE_TABLE * edges,LWT_ISO_EDGE * edge,int side)6338 _lwt_BuildEdgeRing(__attribute__((__unused__)) LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *edges,
6339 LWT_ISO_EDGE *edge, int side)
6340 {
6341 LWT_EDGERING *ring;
6342 LWT_EDGERING_ELEM *elem;
6343 LWT_ISO_EDGE *cur;
6344 int curside;
6345
6346 ring = lwalloc(sizeof(LWT_EDGERING));
6347 LWT_EDGERING_INIT(ring);
6348
6349 cur = edge;
6350 curside = side;
6351
6352 LWDEBUGF(2, "Building rings for edge %d (side %d)", cur->edge_id, side);
6353
6354 do {
6355 LWT_ELEMID next;
6356
6357 elem = lwalloc(sizeof(LWT_EDGERING_ELEM));
6358 elem->edge = cur;
6359 elem->left = ( curside == 1 );
6360
6361 /* Mark edge as "visited" */
6362 if ( elem->left ) cur->face_left = LWT_HOLES_FACE_PLACEHOLDER;
6363 else cur->face_right = LWT_HOLES_FACE_PLACEHOLDER;
6364
6365 LWT_EDGERING_PUSH(ring, elem);
6366 next = elem->left ? cur->next_left : cur->next_right;
6367
6368 LWDEBUGF(3, " next edge is %d", next);
6369
6370 if ( next > 0 ) curside = 1;
6371 else { curside = -1; next = -next; }
6372 cur = _lwt_getIsoEdgeById(edges, next);
6373 if ( ! cur )
6374 {
6375 lwerror("Could not find edge with id %d", next);
6376 break;
6377 }
6378 } while (cur != edge || curside != side);
6379
6380 LWDEBUGF(1, "Ring for edge %d has %d elems", edge->edge_id*side, ring->size);
6381
6382 return ring;
6383 }
6384
6385 static double
_lwt_EdgeRingSignedArea(LWT_EDGERING_POINT_ITERATOR * it)6386 _lwt_EdgeRingSignedArea(LWT_EDGERING_POINT_ITERATOR *it)
6387 {
6388 POINT2D P1;
6389 POINT2D P2;
6390 POINT2D P3;
6391 double sum = 0.0;
6392 double x0, x, y1, y2;
6393
6394 if ( ! _lwt_EdgeRingIterator_next(it, &P1) ) return 0.0;
6395 if ( ! _lwt_EdgeRingIterator_next(it, &P2) ) return 0.0;
6396
6397 LWDEBUG(2, "_lwt_EdgeRingSignedArea");
6398
6399 x0 = P1.x;
6400 while ( _lwt_EdgeRingIterator_next(it, &P3) )
6401 {
6402 x = P2.x - x0;
6403 y1 = P3.y;
6404 y2 = P1.y;
6405 sum += x * (y2-y1);
6406
6407 /* Move forwards! */
6408 P1 = P2;
6409 P2 = P3;
6410 }
6411
6412 return sum / 2.0;
6413 }
6414
6415
6416 /* Return 1 for true, 0 for false */
6417 static int
_lwt_EdgeRingIsCCW(LWT_EDGERING * ring)6418 _lwt_EdgeRingIsCCW(LWT_EDGERING *ring)
6419 {
6420 double sa;
6421
6422 LWDEBUGF(2, "_lwt_EdgeRingIsCCW, ring has %d elems", ring->size);
6423 LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring);
6424 sa = _lwt_EdgeRingSignedArea(it);
6425 LWDEBUGF(2, "_lwt_EdgeRingIsCCW, signed area is %g", sa);
6426 lwfree(it);
6427 if ( sa >= 0 ) return 0;
6428 else return 1;
6429 }
6430
6431 static int
_lwt_EdgeRingCrossingCount(const POINT2D * p,LWT_EDGERING_POINT_ITERATOR * it)6432 _lwt_EdgeRingCrossingCount(const POINT2D *p, LWT_EDGERING_POINT_ITERATOR *it)
6433 {
6434 int cn = 0; /* the crossing number counter */
6435 POINT2D v1, v2;
6436 #ifndef RELAX
6437 POINT2D v0;
6438 #endif
6439
6440 if ( ! _lwt_EdgeRingIterator_next(it, &v1) ) return cn;
6441 v0 = v1;
6442 while ( _lwt_EdgeRingIterator_next(it, &v2) )
6443 {
6444 double vt;
6445
6446 /* edge from vertex i to vertex i+1 */
6447 if
6448 (
6449 /* an upward crossing */
6450 ((v1.y <= p->y) && (v2.y > p->y))
6451 /* a downward crossing */
6452 || ((v1.y > p->y) && (v2.y <= p->y))
6453 )
6454 {
6455
6456 vt = (double)(p->y - v1.y) / (v2.y - v1.y);
6457
6458 /* P->x <intersect */
6459 if (p->x < v1.x + vt * (v2.x - v1.x))
6460 {
6461 /* a valid crossing of y=p->y right of p->x */
6462 ++cn;
6463 }
6464 }
6465 v1 = v2;
6466 }
6467
6468 LWDEBUGF(3, "_lwt_EdgeRingCrossingCount returning %d", cn);
6469
6470 #ifndef RELAX
6471 if ( memcmp(&v1, &v0, sizeof(POINT2D)) )
6472 {
6473 lwerror("_lwt_EdgeRingCrossingCount: V[n] != V[0] (%g %g != %g %g)",
6474 v1.x, v1.y, v0.x, v0.y);
6475 return -1;
6476 }
6477 #endif
6478
6479 return cn;
6480 }
6481
6482 /* Return 1 for true, 0 for false */
6483 static int
_lwt_EdgeRingContainsPoint(LWT_EDGERING * ring,POINT2D * p)6484 _lwt_EdgeRingContainsPoint(LWT_EDGERING *ring, POINT2D *p)
6485 {
6486 int cn = 0;
6487
6488 LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring);
6489 cn = _lwt_EdgeRingCrossingCount(p, it);
6490 lwfree(it);
6491 return (cn&1); /* 0 if even (out), and 1 if odd (in) */
6492 }
6493
6494 static GBOX *
_lwt_EdgeRingGetBbox(LWT_EDGERING * ring)6495 _lwt_EdgeRingGetBbox(LWT_EDGERING *ring)
6496 {
6497 int i;
6498
6499 if ( ! ring->env )
6500 {
6501 LWDEBUGF(2, "Computing GBOX for ring %p", ring);
6502 for (i=0; i<ring->size; ++i)
6503 {
6504 LWT_EDGERING_ELEM *elem = ring->elems[i];
6505 LWLINE *g = elem->edge->geom;
6506 const GBOX *newbox = lwgeom_get_bbox(lwline_as_lwgeom(g));
6507 if ( ! i ) ring->env = gbox_clone( newbox );
6508 else gbox_merge( newbox, ring->env );
6509 }
6510 }
6511
6512 return ring->env;
6513 }
6514
6515 static LWT_ELEMID
_lwt_EdgeRingGetFace(LWT_EDGERING * ring)6516 _lwt_EdgeRingGetFace(LWT_EDGERING *ring)
6517 {
6518 LWT_EDGERING_ELEM *el = ring->elems[0];
6519 return el->left ? el->edge->face_left : el->edge->face_right;
6520 }
6521
6522
6523 /*
6524 * Register a face on an edge side
6525 *
6526 * Create and register face to shell (CCW) walks,
6527 * register arbitrary negative face_id to CW rings.
6528 *
6529 * Push CCW rings to shells, CW rings to holes.
6530 *
6531 * The ownership of the "geom" and "ids" members of the
6532 * LWT_EDGERING pushed to the given LWT_EDGERING_ARRAYS
6533 * are transferred to caller.
6534 *
6535 * @param side 1 for left side, -1 for right side
6536 *
6537 * @param holes an array where holes will be pushed
6538 *
6539 * @param shells an array where shells will be pushed
6540 *
6541 * @param registered id of registered face. It will be a negative number
6542 * for holes or isolated edge strips (still registered in the face
6543 * table, but only temporary).
6544 *
6545 * @return 0 on success, -1 on error.
6546 *
6547 */
6548 static int
_lwt_RegisterFaceOnEdgeSide(LWT_TOPOLOGY * topo,LWT_ISO_EDGE * edge,int side,LWT_ISO_EDGE_TABLE * edges,LWT_EDGERING_ARRAY * holes,LWT_EDGERING_ARRAY * shells,LWT_ELEMID * registered)6549 _lwt_RegisterFaceOnEdgeSide(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edge,
6550 int side, LWT_ISO_EDGE_TABLE *edges,
6551 LWT_EDGERING_ARRAY *holes,
6552 LWT_EDGERING_ARRAY *shells,
6553 LWT_ELEMID *registered)
6554 {
6555 const LWT_BE_IFACE *iface = topo->be_iface;
6556 /* this is arbitrary, could be taken as parameter */
6557 static const int placeholder_faceid = LWT_HOLES_FACE_PLACEHOLDER;
6558 LWT_EDGERING *ring;
6559
6560 /* Get edge ring */
6561 ring = _lwt_BuildEdgeRing(topo, edges, edge, side);
6562
6563 LWDEBUG(2, "Ring built, calling EdgeRingIsCCW");
6564
6565 /* Compute winding (CW or CCW?) */
6566 int isccw = _lwt_EdgeRingIsCCW(ring);
6567
6568 if ( isccw )
6569 {
6570 /* Create new face */
6571 LWT_ISO_FACE newface;
6572
6573 LWDEBUGF(1, "Ring of edge %d is a shell (shell %d)", edge->edge_id * side, shells->size);
6574
6575 newface.mbr = _lwt_EdgeRingGetBbox(ring);
6576
6577 newface.face_id = -1;
6578 /* Insert the new face */
6579 int ret = lwt_be_insertFaces( topo, &newface, 1 );
6580 newface.mbr = NULL;
6581 if ( ret == -1 )
6582 {
6583 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6584 return -1;
6585 }
6586 if ( ret != 1 )
6587 {
6588 lwerror("Unexpected error: %d faces inserted when expecting 1", ret);
6589 return -1;
6590 }
6591 /* return new face_id */
6592 *registered = newface.face_id;
6593 LWT_EDGERING_ARRAY_PUSH(shells, ring);
6594
6595 /* update ring edges set new face_id on resp. side to *registered */
6596 ret = _lwt_UpdateEdgeRingSideFace(topo, ring, *registered);
6597 if ( ret )
6598 {
6599 lwerror("Errors updating edgering side face: %s",
6600 lwt_be_lastErrorMessage(iface));
6601 return -1;
6602 }
6603
6604 }
6605 else /* cw, so is an hole */
6606 {
6607 LWDEBUGF(1, "Ring of edge %d is a hole (hole %d)", edge->edge_id * side, holes->size);
6608 *registered = placeholder_faceid;
6609 LWT_EDGERING_ARRAY_PUSH(holes, ring);
6610 }
6611
6612 return 0;
6613 }
6614
6615 static void
_lwt_AccumulateCanditates(void * item,void * userdata)6616 _lwt_AccumulateCanditates(void* item, void* userdata)
6617 {
6618 LWT_EDGERING_ARRAY *candidates = userdata;
6619 LWT_EDGERING *sring = item;
6620 LWT_EDGERING_ARRAY_PUSH(candidates, sring);
6621 }
6622
6623 static LWT_ELEMID
_lwt_FindFaceContainingRing(LWT_TOPOLOGY * topo,LWT_EDGERING * ring,LWT_EDGERING_ARRAY * shells)6624 _lwt_FindFaceContainingRing(LWT_TOPOLOGY* topo, LWT_EDGERING *ring,
6625 LWT_EDGERING_ARRAY *shells)
6626 {
6627 LWT_ELEMID foundInFace = -1;
6628 int i;
6629 const GBOX *minenv = NULL;
6630 POINT2D pt;
6631 const GBOX *testbox;
6632 GEOSGeometry *ghole;
6633
6634 getPoint2d_p( ring->elems[0]->edge->geom->points, 0, &pt );
6635
6636 testbox = _lwt_EdgeRingGetBbox(ring);
6637
6638 /* Create a GEOS Point from a vertex of the hole ring */
6639 {
6640 LWPOINT *point = lwpoint_make2d(topo->srid, pt.x, pt.y);
6641 ghole = LWGEOM2GEOS( lwpoint_as_lwgeom(point), 1 );
6642 lwpoint_free(point);
6643 if ( ! ghole ) {
6644 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
6645 return -1;
6646 }
6647 }
6648
6649 /* Build STRtree of shell envelopes */
6650 if ( ! shells->tree )
6651 {
6652 static const int STRTREE_NODE_CAPACITY = 10;
6653 LWDEBUG(1, "Building STRtree");
6654 shells->tree = GEOSSTRtree_create(STRTREE_NODE_CAPACITY);
6655 if (shells->tree == NULL)
6656 {
6657 lwerror("Could not create GEOS STRTree: %s", lwgeom_geos_errmsg);
6658 return -1;
6659 }
6660 for (i=0; i<shells->size; ++i)
6661 {
6662 LWT_EDGERING *sring = shells->rings[i];
6663 const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring);
6664 LWDEBUGF(2, "GBOX of shell %p for edge %d is %g %g,%g %g",
6665 sring, sring->elems[0]->edge->edge_id, shellbox->xmin,
6666 shellbox->ymin, shellbox->xmax, shellbox->ymax);
6667 POINTARRAY *pa = ptarray_construct(0, 0, 2);
6668 POINT4D pt;
6669 LWLINE *diag;
6670 pt.x = shellbox->xmin;
6671 pt.y = shellbox->ymin;
6672 ptarray_set_point4d(pa, 0, &pt);
6673 pt.x = shellbox->xmax;
6674 pt.y = shellbox->ymax;
6675 ptarray_set_point4d(pa, 1, &pt);
6676 diag = lwline_construct(topo->srid, NULL, pa);
6677 /* Record just envelope in ggeom */
6678 /* making valid, probably not needed */
6679 sring->genv = LWGEOM2GEOS( lwline_as_lwgeom(diag), 1 );
6680 lwline_free(diag);
6681 GEOSSTRtree_insert(shells->tree, sring->genv, sring);
6682 }
6683 LWDEBUG(1, "STRtree build completed");
6684 }
6685
6686 LWT_EDGERING_ARRAY candidates;
6687 LWT_EDGERING_ARRAY_INIT(&candidates);
6688 GEOSSTRtree_query(shells->tree, ghole, &_lwt_AccumulateCanditates, &candidates);
6689 LWDEBUGF(1, "Found %d candidate shells containing first point of ring's originating edge %d",
6690 candidates.size, ring->elems[0]->edge->edge_id * ( ring->elems[0]->left ? 1 : -1 ) );
6691
6692 /* TODO: sort candidates by bounding box size */
6693
6694 for (i=0; i<candidates.size; ++i)
6695 {
6696 LWT_EDGERING *sring = candidates.rings[i];
6697 const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring);
6698 int contains = 0;
6699
6700 if ( sring->elems[0]->edge->edge_id == ring->elems[0]->edge->edge_id )
6701 {
6702 LWDEBUGF(1, "Shell %d is on other side of ring",
6703 _lwt_EdgeRingGetFace(sring));
6704 continue;
6705 }
6706
6707 /* The hole envelope cannot equal the shell envelope */
6708 if ( gbox_same(shellbox, testbox) )
6709 {
6710 LWDEBUGF(1, "Bbox of shell %d equals that of hole ring",
6711 _lwt_EdgeRingGetFace(sring));
6712 continue;
6713 }
6714
6715 /* Skip if ring box is not in shell box */
6716 if ( ! gbox_contains_2d(shellbox, testbox) )
6717 {
6718 LWDEBUGF(1, "Bbox of shell %d does not contain bbox of ring point",
6719 _lwt_EdgeRingGetFace(sring));
6720 continue;
6721 }
6722
6723 /* Skip test if a containing shell was already found
6724 * and this shell's bbox is not contained in the other */
6725 if ( minenv && ! gbox_contains_2d(minenv, shellbox) )
6726 {
6727 LWDEBUGF(2, "Bbox of shell %d (face %d) not contained by bbox "
6728 "of last shell found to contain the point",
6729 i, _lwt_EdgeRingGetFace(sring));
6730 continue;
6731 }
6732
6733 contains = _lwt_EdgeRingContainsPoint(sring, &pt);
6734 if ( contains )
6735 {
6736 /* Continue until all shells are tested, as we want to
6737 * use the one with the smallest bounding box */
6738 /* IDEA: sort shells by bbox size, stopping on first match */
6739 LWDEBUGF(1, "Shell %d contains hole of edge %d",
6740 _lwt_EdgeRingGetFace(sring),
6741 ring->elems[0]->edge->edge_id);
6742 minenv = shellbox;
6743 foundInFace = _lwt_EdgeRingGetFace(sring);
6744 }
6745 }
6746 if ( foundInFace == -1 ) foundInFace = 0;
6747
6748 candidates.size = 0; /* Avoid destroying the actual shell rings */
6749 LWT_EDGERING_ARRAY_CLEAN(&candidates);
6750
6751 GEOSGeom_destroy(ghole);
6752
6753 return foundInFace;
6754 }
6755
6756 /*
6757 * @return -1 on error (and report error),
6758 * 1 if faces beside the universal one exist
6759 * 0 otherwise
6760 */
6761 static int
_lwt_CheckFacesExist(LWT_TOPOLOGY * topo)6762 _lwt_CheckFacesExist(LWT_TOPOLOGY *topo)
6763 {
6764 LWT_ISO_FACE *faces;
6765 int fields = LWT_COL_FACE_FACE_ID;
6766 uint64_t nelems = 1;
6767 GBOX qbox;
6768
6769 qbox.xmin = qbox.ymin = -DBL_MAX;
6770 qbox.xmax = qbox.ymax = DBL_MAX;
6771 faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nelems, fields, 1);
6772 if (nelems == UINT64_MAX)
6773 {
6774 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6775 return -1;
6776 }
6777 if ( faces ) _lwt_release_faces(faces, nelems);
6778 return nelems;
6779 }
6780
6781 int
lwt_Polygonize(LWT_TOPOLOGY * topo)6782 lwt_Polygonize(LWT_TOPOLOGY* topo)
6783 {
6784 /*
6785 Fetch all edges
6786 Sort edges by edge_id
6787 Mark all edges' left and right face as -1
6788 Iteratively:
6789 Fetch next edge with left or right face == -1
6790 For each side with face == -1:
6791 Find ring on its side
6792 If ring is CCW:
6793 create a new face, assign to the ring edges' appropriate side
6794 If ring is CW (face needs to be same of external):
6795 assign a negative face_id the ring edges' appropriate side
6796 Now for each edge with a negative face_id on the side:
6797 Find containing face (mbr cache and all)
6798 Update with id of containing face
6799 */
6800
6801 const LWT_BE_IFACE *iface = topo->be_iface;
6802 LWT_ISO_EDGE *edge;
6803 int numfaces = -1;
6804 LWT_ISO_EDGE_TABLE edgetable;
6805 LWT_EDGERING_ARRAY holes, shells;
6806 int i;
6807 int err = 0;
6808
6809 LWT_EDGERING_ARRAY_INIT(&holes);
6810 LWT_EDGERING_ARRAY_INIT(&shells);
6811
6812 initGEOS(lwnotice, lwgeom_geos_error);
6813
6814 /*
6815 Check if Topology already contains some Face
6816 (ignoring the Universal Face)
6817 */
6818 numfaces = _lwt_CheckFacesExist(topo);
6819 if ( numfaces != 0 ) {
6820 if ( numfaces > 0 ) {
6821 /* Faces exist */
6822 lwerror("Polygonize: face table is not empty.");
6823 }
6824 /* Backend error, message should have been printed already */
6825 return -1;
6826 }
6827
6828
6829 edgetable.edges = _lwt_FetchAllEdges(topo, &(edgetable.size));
6830 if ( ! edgetable.edges ) {
6831 if (edgetable.size == 0) {
6832 /* not an error: no Edges */
6833 return 0;
6834 }
6835 /* error should have been printed already */
6836 return -1;
6837 }
6838
6839 /* Sort edges by ID (to allow btree searches) */
6840 qsort(edgetable.edges, edgetable.size, sizeof(LWT_ISO_EDGE), compare_iso_edges_by_id);
6841
6842 /* Mark all edges as unvisited */
6843 for (i=0; i<edgetable.size; ++i)
6844 edgetable.edges[i].face_left = edgetable.edges[i].face_right = -1;
6845
6846 i = 0;
6847 while (1)
6848 {
6849 i = _lwt_FetchNextUnvisitedEdge(topo, &edgetable, i);
6850 if ( i < 0 ) break; /* end of unvisited */
6851 edge = &(edgetable.edges[i]);
6852
6853 LWT_ELEMID newface = -1;
6854
6855 LWDEBUGF(1, "Next face-missing edge has id:%d, face_left:%d, face_right:%d",
6856 edge->edge_id, edge->face_left, edge->face_right);
6857 if ( edge->face_left == -1 )
6858 {
6859 err = _lwt_RegisterFaceOnEdgeSide(topo, edge, 1, &edgetable,
6860 &holes, &shells, &newface);
6861 if ( err ) break;
6862 LWDEBUGF(1, "New face on the left of edge %d is %d",
6863 edge->edge_id, newface);
6864 edge->face_left = newface;
6865 }
6866 if ( edge->face_right == -1 )
6867 {
6868 err = _lwt_RegisterFaceOnEdgeSide(topo, edge, -1, &edgetable,
6869 &holes, &shells, &newface);
6870 if ( err ) break;
6871 LWDEBUGF(1, "New face on the right of edge %d is %d",
6872 edge->edge_id, newface);
6873 edge->face_right = newface;
6874 }
6875 }
6876
6877 if ( err )
6878 {
6879 _lwt_release_edges(edgetable.edges, edgetable.size);
6880 LWT_EDGERING_ARRAY_CLEAN( &holes );
6881 LWT_EDGERING_ARRAY_CLEAN( &shells );
6882 lwerror("Errors fetching or registering face-missing edges: %s",
6883 lwt_be_lastErrorMessage(iface));
6884 return -1;
6885 }
6886
6887 LWDEBUGF(1, "Found %d holes and %d shells", holes.size, shells.size);
6888
6889 /* TODO: sort holes by pt.x, sort shells by bbox.xmin */
6890
6891 /* Assign shells to holes */
6892 for (i=0; i<holes.size; ++i)
6893 {
6894 LWT_ELEMID containing_face;
6895 LWT_EDGERING *ring = holes.rings[i];
6896
6897 containing_face = _lwt_FindFaceContainingRing(topo, ring, &shells);
6898 LWDEBUGF(1, "Ring %d contained by face %" LWTFMT_ELEMID, i, containing_face);
6899 if ( containing_face == -1 )
6900 {
6901 _lwt_release_edges(edgetable.edges, edgetable.size);
6902 LWT_EDGERING_ARRAY_CLEAN( &holes );
6903 LWT_EDGERING_ARRAY_CLEAN( &shells );
6904 lwerror("Errors finding face containing ring: %s",
6905 lwt_be_lastErrorMessage(iface));
6906 return -1;
6907 }
6908 int ret = _lwt_UpdateEdgeRingSideFace(topo, holes.rings[i], containing_face);
6909 if ( ret )
6910 {
6911 _lwt_release_edges(edgetable.edges, edgetable.size);
6912 LWT_EDGERING_ARRAY_CLEAN( &holes );
6913 LWT_EDGERING_ARRAY_CLEAN( &shells );
6914 lwerror("Errors updating edgering side face: %s",
6915 lwt_be_lastErrorMessage(iface));
6916 return -1;
6917 }
6918 }
6919
6920 LWDEBUG(1, "All holes assigned, cleaning up");
6921
6922 _lwt_release_edges(edgetable.edges, edgetable.size);
6923
6924 /* delete all shell and hole EDGERINGS */
6925 LWT_EDGERING_ARRAY_CLEAN( &holes );
6926 LWT_EDGERING_ARRAY_CLEAN( &shells );
6927
6928 return 0;
6929 }
6930