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;
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_FACE_LEFT |
2808 LWT_COL_EDGE_FACE_RIGHT
2809 ;
2810 edges = lwt_be_getEdgeByFace( topo, &faceid, &numfaceedges, fields, NULL );
2811 if (numfaceedges == UINT64_MAX)
2812 {
2813 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2814 return NULL;
2815 }
2816
2817 if ( numfaceedges == 0 )
2818 {
2819 i = 1;
2820 face = lwt_be_getFaceById(topo, &faceid, &i, LWT_COL_FACE_FACE_ID);
2821 if (i == UINT64_MAX)
2822 {
2823 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
2824 return NULL;
2825 }
2826 if ( i == 0 ) {
2827 lwerror("SQL/MM Spatial exception - non-existent face.");
2828 return NULL;
2829 }
2830 lwfree( face );
2831 if ( i > 1 ) {
2832 lwerror("Corrupted topology: multiple face records have face_id=%"
2833 LWTFMT_ELEMID, faceid);
2834 return NULL;
2835 }
2836 /* Face has no boundary edges, we'll return EMPTY, see
2837 * https://trac.osgeo.org/postgis/ticket/3221 */
2838 out = lwpoly_construct_empty(topo->srid, topo->hasZ, 0);
2839 return lwpoly_as_lwgeom(out);
2840 }
2841
2842 outg = _lwt_FaceByEdges( topo, edges, numfaceedges );
2843 _lwt_release_edges(edges, numfaceedges);
2844
2845 return outg;
2846 }
2847
2848 /* Find which edge from the "edges" set defines the next
2849 * portion of the given "ring".
2850 *
2851 * The edge might be either forward or backward.
2852 *
2853 * @param ring The ring to find definition of.
2854 * It is assumed it does not contain duplicated vertices.
2855 * @param from offset of the ring point to start looking from
2856 * @param edges array of edges to search into
2857 * @param numedges number of edges in the edges array
2858 *
2859 * @return index of the edge defining the next ring portion or
2860 * -1 if no edge was found to be part of the ring
2861 */
2862 static int
_lwt_FindNextRingEdge(const POINTARRAY * ring,int from,const LWT_ISO_EDGE * edges,int numedges)2863 _lwt_FindNextRingEdge(const POINTARRAY *ring, int from,
2864 const LWT_ISO_EDGE *edges, int numedges)
2865 {
2866 int i;
2867 POINT2D p1;
2868
2869 /* Get starting ring point */
2870 getPoint2d_p(ring, from, &p1);
2871
2872 LWDEBUGF(1, "Ring's 'from' point (%d) is %g,%g", from, p1.x, p1.y);
2873
2874 /* find the edges defining the next portion of ring starting from
2875 * vertex "from" */
2876 for ( i=0; i<numedges; ++i )
2877 {
2878 const LWT_ISO_EDGE *isoe = &(edges[i]);
2879 LWLINE *edge = isoe->geom;
2880 POINTARRAY *epa = edge->points;
2881 POINT2D p2, pt;
2882 int match = 0;
2883 uint32_t j;
2884
2885 /* Skip if the edge is a dangling one */
2886 if ( isoe->face_left == isoe->face_right )
2887 {
2888 LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID
2889 " has same face (%" LWTFMT_ELEMID
2890 ") on both sides, skipping",
2891 isoe->edge_id, isoe->face_left);
2892 continue;
2893 }
2894
2895 if (epa->npoints < 2)
2896 {
2897 LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID
2898 " has only %"PRIu32" points",
2899 isoe->edge_id, epa->npoints);
2900 continue;
2901 }
2902
2903 #if 0
2904 size_t sz;
2905 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s",
2906 isoe->edge_id,
2907 lwgeom_to_wkt(lwline_as_lwgeom(edge), WKT_EXTENDED, 2, &sz));
2908 #endif
2909
2910 /* ptarray_remove_repeated_points ? */
2911
2912 getPoint2d_p(epa, 0, &p2);
2913 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'first' point is %g,%g",
2914 isoe->edge_id, p2.x, p2.y);
2915 LWDEBUGF(1, "Rings's 'from' point is still %g,%g", p1.x, p1.y);
2916 if ( p2d_same(&p1, &p2) )
2917 {
2918 LWDEBUG(1, "p2d_same(p1,p2) returned true");
2919 LWDEBUGF(1, "First point of edge %" LWTFMT_ELEMID
2920 " matches ring vertex %d", isoe->edge_id, from);
2921 /* first point matches, let's check next non-equal one */
2922 for ( j=1; j<epa->npoints; ++j )
2923 {
2924 getPoint2d_p(epa, j, &p2);
2925 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'next' point %d is %g,%g",
2926 isoe->edge_id, j, p2.x, p2.y);
2927 /* we won't check duplicated edge points */
2928 if ( p2d_same(&p1, &p2) ) continue;
2929 /* we assume there are no duplicated points in ring */
2930 getPoint2d_p(ring, from+1, &pt);
2931 LWDEBUGF(1, "Ring's point %d is %g,%g",
2932 from+1, pt.x, pt.y);
2933 match = p2d_same(&pt, &p2);
2934 break; /* we want to check a single non-equal next vertex */
2935 }
2936 #if POSTGIS_DEBUG_LEVEL > 0
2937 if ( match ) {
2938 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2939 " matches ring vertex %d", isoe->edge_id, from+1);
2940 } else {
2941 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2942 " does not match ring vertex %d", isoe->edge_id, from+1);
2943 }
2944 #endif
2945 }
2946
2947 if ( ! match )
2948 {
2949 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " did not match as forward",
2950 isoe->edge_id);
2951 getPoint2d_p(epa, epa->npoints-1, &p2);
2952 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'last' point is %g,%g",
2953 isoe->edge_id, p2.x, p2.y);
2954 if ( p2d_same(&p1, &p2) )
2955 {
2956 LWDEBUGF(1, "Last point of edge %" LWTFMT_ELEMID
2957 " matches ring vertex %d", isoe->edge_id, from);
2958 /* last point matches, let's check next non-equal one */
2959 for ( j=2; j<=epa->npoints; j++ )
2960 {
2961 getPoint2d_p(epa, epa->npoints - j, &p2);
2962 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'prev' point %d is %g,%g",
2963 isoe->edge_id, epa->npoints - j, p2.x, p2.y);
2964 /* we won't check duplicated edge points */
2965 if ( p2d_same(&p1, &p2) ) continue;
2966 /* we assume there are no duplicated points in ring */
2967 getPoint2d_p(ring, from+1, &pt);
2968 LWDEBUGF(1, "Ring's point %d is %g,%g",
2969 from+1, pt.x, pt.y);
2970 match = p2d_same(&pt, &p2);
2971 break; /* we want to check a single non-equal next vertex */
2972 }
2973 }
2974 #if POSTGIS_DEBUG_LEVEL > 0
2975 if ( match ) {
2976 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2977 " matches ring vertex %d", isoe->edge_id, from+1);
2978 } else {
2979 LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID
2980 " does not match ring vertex %d", isoe->edge_id, from+1);
2981 }
2982 #endif
2983 }
2984
2985 if ( match ) return i;
2986
2987 }
2988
2989 return -1;
2990 }
2991
2992 /* Reverse values in array between "from" (inclusive)
2993 * and "to" (exclusive) indexes */
2994 static void
_lwt_ReverseElemidArray(LWT_ELEMID * ary,int from,int to)2995 _lwt_ReverseElemidArray(LWT_ELEMID *ary, int from, int to)
2996 {
2997 LWT_ELEMID t;
2998 while (from < to)
2999 {
3000 t = ary[from];
3001 ary[from++] = ary[to];
3002 ary[to--] = t;
3003 }
3004 }
3005
3006 /* Rotate values in array between "from" (inclusive)
3007 * and "to" (exclusive) indexes, so that "rotidx" is
3008 * the new value at "from" */
3009 static void
_lwt_RotateElemidArray(LWT_ELEMID * ary,int from,int to,int rotidx)3010 _lwt_RotateElemidArray(LWT_ELEMID *ary, int from, int to, int rotidx)
3011 {
3012 _lwt_ReverseElemidArray(ary, from, rotidx-1);
3013 _lwt_ReverseElemidArray(ary, rotidx, to-1);
3014 _lwt_ReverseElemidArray(ary, from, to-1);
3015 }
3016
3017
3018 int
lwt_GetFaceEdges(LWT_TOPOLOGY * topo,LWT_ELEMID face_id,LWT_ELEMID ** out)3019 lwt_GetFaceEdges(LWT_TOPOLOGY* topo, LWT_ELEMID face_id, LWT_ELEMID **out )
3020 {
3021 LWGEOM *face;
3022 LWPOLY *facepoly;
3023 LWT_ISO_EDGE *edges;
3024 uint64_t numfaceedges;
3025 int fields;
3026 uint32_t i;
3027 int nseid = 0; /* number of signed edge ids */
3028 int prevseid;
3029 LWT_ELEMID *seid; /* signed edge ids */
3030
3031 /* Get list of face edges */
3032 numfaceedges = 1;
3033 fields = LWT_COL_EDGE_EDGE_ID |
3034 LWT_COL_EDGE_GEOM |
3035 LWT_COL_EDGE_FACE_LEFT |
3036 LWT_COL_EDGE_FACE_RIGHT
3037 ;
3038 edges = lwt_be_getEdgeByFace( topo, &face_id, &numfaceedges, fields, NULL );
3039 if (numfaceedges == UINT64_MAX)
3040 {
3041 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3042 return -1;
3043 }
3044 if ( ! numfaceedges ) return 0; /* no edges in output */
3045
3046 /* order edges by occurrence in face */
3047
3048 face = _lwt_FaceByEdges(topo, edges, numfaceedges);
3049 if ( ! face )
3050 {
3051 /* _lwt_FaceByEdges should have already invoked lwerror in this case */
3052 _lwt_release_edges(edges, numfaceedges);
3053 return -1;
3054 }
3055
3056 if ( lwgeom_is_empty(face) )
3057 {
3058 /* no edges in output */
3059 _lwt_release_edges(edges, numfaceedges);
3060 lwgeom_free(face);
3061 return 0;
3062 }
3063
3064 /* force_lhr, if the face is not the universe */
3065 /* _lwt_FaceByEdges seems to guaranteed RHR */
3066 /* lwgeom_force_clockwise(face); */
3067 if ( face_id ) lwgeom_reverse_in_place(face);
3068
3069 #if 0
3070 {
3071 size_t sz;
3072 char *wkt = lwgeom_to_wkt(face, WKT_EXTENDED, 6, &sz);
3073 LWDEBUGF(1, "Geometry of face %" LWTFMT_ELEMID " is: %s",
3074 face_id, wkt);
3075 lwfree(wkt);
3076 }
3077 #endif
3078
3079 facepoly = lwgeom_as_lwpoly(face);
3080 if ( ! facepoly )
3081 {
3082 _lwt_release_edges(edges, numfaceedges);
3083 lwgeom_free(face);
3084 lwerror("Geometry of face %" LWTFMT_ELEMID " is not a polygon", face_id);
3085 return -1;
3086 }
3087
3088 nseid = prevseid = 0;
3089 seid = lwalloc( sizeof(LWT_ELEMID) * numfaceedges );
3090
3091 /* for each ring of the face polygon... */
3092 for ( i=0; i<facepoly->nrings; ++i )
3093 {
3094 const POINTARRAY *ring = facepoly->rings[i];
3095 int32_t j = 0;
3096 LWT_ISO_EDGE *nextedge;
3097 LWLINE *nextline;
3098
3099 LWDEBUGF(1, "Ring %d has %d points", i, ring->npoints);
3100
3101 while ( j < (int32_t) ring->npoints-1 )
3102 {
3103 LWDEBUGF(1, "Looking for edge covering ring %d from vertex %d",
3104 i, j);
3105
3106 int edgeno = _lwt_FindNextRingEdge(ring, j, edges, numfaceedges);
3107 if ( edgeno == -1 )
3108 {
3109 /* should never happen */
3110 _lwt_release_edges(edges, numfaceedges);
3111 lwgeom_free(face);
3112 lwfree(seid);
3113 lwerror("No edge (among %d) found to be defining geometry of face %"
3114 LWTFMT_ELEMID, numfaceedges, face_id);
3115 return -1;
3116 }
3117
3118 nextedge = &(edges[edgeno]);
3119 nextline = nextedge->geom;
3120
3121 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
3122 " covers ring %d from vertex %d to %d",
3123 nextedge->edge_id, i, j, j + nextline->points->npoints - 1);
3124
3125 #if 0
3126 {
3127 size_t sz;
3128 char *wkt = lwgeom_to_wkt(lwline_as_lwgeom(nextline), WKT_EXTENDED, 6, &sz);
3129 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s",
3130 nextedge->edge_id, wkt);
3131 lwfree(wkt);
3132 }
3133 #endif
3134
3135 j += nextline->points->npoints - 1;
3136
3137 /* Add next edge to the output array */
3138 seid[nseid++] = nextedge->face_left == face_id ?
3139 nextedge->edge_id :
3140 -nextedge->edge_id;
3141
3142 /* avoid checking again on next time turn */
3143 nextedge->face_left = nextedge->face_right = -1;
3144 }
3145
3146 /* now "scroll" the list of edges so that the one
3147 * with smaller absolute edge_id is first */
3148 /* Range is: [prevseid, nseid) -- [inclusive, exclusive) */
3149 if ( (nseid - prevseid) > 1 )
3150 {{
3151 LWT_ELEMID minid = 0;
3152 int minidx = 0;
3153 LWDEBUGF(1, "Looking for smallest id among the %d edges "
3154 "composing ring %d", (nseid-prevseid), i);
3155 for ( j=prevseid; j<nseid; ++j )
3156 {
3157 LWT_ELEMID id = llabs(seid[j]);
3158 LWDEBUGF(1, "Abs id of edge in pos %d is %" LWTFMT_ELEMID, j, id);
3159 if ( ! minid || id < minid )
3160 {
3161 minid = id;
3162 minidx = j;
3163 }
3164 }
3165 LWDEBUGF(1, "Smallest id is %" LWTFMT_ELEMID
3166 " at position %d", minid, minidx);
3167 if ( minidx != prevseid )
3168 {
3169 _lwt_RotateElemidArray(seid, prevseid, nseid, minidx);
3170 }
3171 }}
3172
3173 prevseid = nseid;
3174 }
3175
3176 lwgeom_free(face);
3177 _lwt_release_edges(edges, numfaceedges);
3178
3179 *out = seid;
3180 return nseid;
3181 }
3182
3183 int
lwt_ChangeEdgeGeom(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id,LWLINE * geom)3184 lwt_ChangeEdgeGeom(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, LWLINE *geom)
3185 {
3186 LWT_ISO_EDGE *oldedge;
3187 LWT_ISO_EDGE newedge;
3188 POINT2D p1, p2, pt;
3189 uint64_t i;
3190 int isclosed = 0;
3191
3192 /* curve must be simple */
3193 if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) )
3194 {
3195 lwerror("SQL/MM Spatial exception - curve not simple");
3196 return -1;
3197 }
3198
3199 i = 1;
3200 oldedge = lwt_be_getEdgeById(topo, &edge_id, &i, LWT_COL_EDGE_ALL);
3201 if ( ! oldedge )
3202 {
3203 LWDEBUGF(1, "lwt_ChangeEdgeGeom: "
3204 "lwt_be_getEdgeById returned NULL and set i=%d", i);
3205 if (i == UINT64_MAX)
3206 {
3207 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3208 return -1;
3209 }
3210 else if ( i == 0 )
3211 {
3212 lwerror("SQL/MM Spatial exception - non-existent edge %"
3213 LWTFMT_ELEMID, edge_id);
3214 return -1;
3215 }
3216 else
3217 {
3218 lwerror("Backend coding error: getEdgeById callback returned NULL "
3219 "but numelements output parameter has value %d "
3220 "(expected 0 or 1)", i);
3221 return -1;
3222 }
3223 }
3224
3225 LWDEBUGF(1, "lwt_ChangeEdgeGeom: "
3226 "old edge has %d points, new edge has %d points",
3227 oldedge->geom->points->npoints, geom->points->npoints);
3228
3229 /*
3230 * e) Check StartPoint consistency
3231 */
3232 getPoint2d_p(oldedge->geom->points, 0, &p1);
3233 getPoint2d_p(geom->points, 0, &pt);
3234 if ( ! p2d_same(&p1, &pt) )
3235 {
3236 _lwt_release_edges(oldedge, 1);
3237 lwerror("SQL/MM Spatial exception - "
3238 "start node not geometry start point.");
3239 return -1;
3240 }
3241
3242 /*
3243 * f) Check EndPoint consistency
3244 */
3245 if ( oldedge->geom->points->npoints < 2 )
3246 {
3247 _lwt_release_edges(oldedge, 1);
3248 lwerror("Corrupted topology: edge %" LWTFMT_ELEMID
3249 " has less than 2 vertices", oldedge->edge_id);
3250 return -1;
3251 }
3252 getPoint2d_p(oldedge->geom->points, oldedge->geom->points->npoints-1, &p2);
3253 if ( geom->points->npoints < 2 )
3254 {
3255 _lwt_release_edges(oldedge, 1);
3256 lwerror("Invalid edge: less than 2 vertices");
3257 return -1;
3258 }
3259 getPoint2d_p(geom->points, geom->points->npoints-1, &pt);
3260 if ( ! p2d_same(&pt, &p2) )
3261 {
3262 _lwt_release_edges(oldedge, 1);
3263 lwerror("SQL/MM Spatial exception - "
3264 "end node not geometry end point.");
3265 return -1;
3266 }
3267
3268 /* Not in the specs:
3269 * if the edge is closed, check we didn't change winding !
3270 * (should be part of isomorphism checking)
3271 */
3272 if ( oldedge->start_node == oldedge->end_node )
3273 {
3274 isclosed = 1;
3275 #if 1 /* TODO: this is actually bogus as a test */
3276 /* check for valid edge (distinct vertices must exist) */
3277 if ( ! _lwt_GetInteriorEdgePoint(geom, &pt) )
3278 {
3279 _lwt_release_edges(oldedge, 1);
3280 lwerror("Invalid edge (no two distinct vertices exist)");
3281 return -1;
3282 }
3283 #endif
3284
3285 if ( ptarray_isccw(oldedge->geom->points) !=
3286 ptarray_isccw(geom->points) )
3287 {
3288 _lwt_release_edges(oldedge, 1);
3289 lwerror("Edge twist at node POINT(%g %g)", p1.x, p1.y);
3290 return -1;
3291 }
3292 }
3293
3294 if ( _lwt_CheckEdgeCrossing(topo, oldedge->start_node,
3295 oldedge->end_node, geom, edge_id ) )
3296 {
3297 /* would have called lwerror already, leaking :( */
3298 _lwt_release_edges(oldedge, 1);
3299 return -1;
3300 }
3301
3302 LWDEBUG(1, "lwt_ChangeEdgeGeom: "
3303 "edge crossing check passed ");
3304
3305 /*
3306 * Not in the specs:
3307 * Check topological isomorphism
3308 */
3309
3310 /* Check that the "motion range" doesn't include any node */
3311 // 1. compute combined bbox of old and new edge
3312 GBOX mbox; /* motion box */
3313 lwgeom_add_bbox((LWGEOM*)oldedge->geom); /* just in case */
3314 lwgeom_add_bbox((LWGEOM*)geom); /* just in case */
3315 gbox_union(oldedge->geom->bbox, geom->bbox, &mbox);
3316 // 2. fetch all nodes in the combined box
3317 LWT_ISO_NODE *nodes;
3318 uint64_t numnodes;
3319 nodes = lwt_be_getNodeWithinBox2D(topo, &mbox, &numnodes,
3320 LWT_COL_NODE_ALL, 0);
3321 LWDEBUGF(1, "lwt_be_getNodeWithinBox2D returned %d nodes", numnodes);
3322 if (numnodes == UINT64_MAX)
3323 {
3324 _lwt_release_edges(oldedge, 1);
3325 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3326 return -1;
3327 }
3328 // 3. if any node beside endnodes are found:
3329 if ( numnodes > ( 1 + isclosed ? 0 : 1 ) )
3330 {{
3331 // 3.2. bail out if any node is in one and not the other
3332 for (i=0; i<numnodes; ++i)
3333 {
3334 LWT_ISO_NODE *n = &(nodes[i]);
3335 if ( n->node_id == oldedge->start_node ) continue;
3336 if ( n->node_id == oldedge->end_node ) continue;
3337 const POINT2D *pt = getPoint2d_cp(n->geom->point, 0);
3338 int ocont = ptarray_contains_point_partial(oldedge->geom->points, pt, isclosed, NULL) == LW_INSIDE;
3339 int ncont = ptarray_contains_point_partial(geom->points, pt, isclosed, NULL) == LW_INSIDE;
3340 if (ocont != ncont)
3341 {
3342 size_t sz;
3343 char *wkt = lwgeom_to_wkt(lwpoint_as_lwgeom(n->geom), WKT_ISO, 15, &sz);
3344 _lwt_release_nodes(nodes, numnodes);
3345 lwerror("Edge motion collision at %s", wkt);
3346 lwfree(wkt); /* would not necessarely reach this point */
3347 return -1;
3348 }
3349 }
3350 }}
3351 if ( numnodes ) _lwt_release_nodes(nodes, numnodes);
3352
3353 LWDEBUG(1, "nodes containment check passed");
3354
3355 /*
3356 * Check edge adjacency before
3357 * TODO: can be optimized to gather azimuths of all edge ends once
3358 */
3359
3360 edgeend span_pre, epan_pre;
3361 /* initialize span_pre.myaz and epan_pre.myaz with existing edge */
3362 int res = _lwt_InitEdgeEndByLine(&span_pre, &epan_pre, oldedge->geom, &p1, &p2);
3363 if (res)
3364 return -1; /* lwerror should have been raised */
3365 _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_pre,
3366 isclosed ? &epan_pre : NULL, edge_id );
3367 _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_pre,
3368 isclosed ? &span_pre : NULL, edge_id );
3369
3370 LWDEBUGF(1, "edges adjacent to old edge are %" LWTFMT_ELEMID
3371 " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID
3372 " and %" LWTFMT_ELEMID " (last point)",
3373 span_pre.nextCW, span_pre.nextCCW,
3374 epan_pre.nextCW, epan_pre.nextCCW);
3375
3376 /* update edge geometry */
3377 newedge.edge_id = edge_id;
3378 newedge.geom = geom;
3379 res = lwt_be_updateEdgesById(topo, &newedge, 1, LWT_COL_EDGE_GEOM);
3380 if (res == -1)
3381 {
3382 _lwt_release_edges(oldedge, 1);
3383 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3384 return -1;
3385 }
3386 if (!res)
3387 {
3388 _lwt_release_edges(oldedge, 1);
3389 lwerror("Unexpected error: %d edges updated when expecting 1", i);
3390 return -1;
3391 }
3392
3393 /*
3394 * Check edge adjacency after
3395 */
3396 edgeend span_post, epan_post;
3397 /* initialize epan_post.myaz and epan_post.myaz */
3398 res = _lwt_InitEdgeEndByLine(&span_post, &epan_post, geom, &p1, &p2);
3399 if (res)
3400 return -1; /* lwerror should have been raised */
3401 _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_post,
3402 isclosed ? &epan_post : NULL, edge_id );
3403 _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_post,
3404 isclosed ? &span_post : NULL, edge_id );
3405
3406 LWDEBUGF(1, "edges adjacent to new edge are %" LWTFMT_ELEMID
3407 " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID
3408 " and %" LWTFMT_ELEMID " (last point)",
3409 span_pre.nextCW, span_pre.nextCCW,
3410 epan_pre.nextCW, epan_pre.nextCCW);
3411
3412
3413 /* Bail out if next CW or CCW edge on start node changed */
3414 if ( span_pre.nextCW != span_post.nextCW ||
3415 span_pre.nextCCW != span_post.nextCCW )
3416 {{
3417 LWT_ELEMID nid = oldedge->start_node;
3418 _lwt_release_edges(oldedge, 1);
3419 lwerror("Edge changed disposition around start node %"
3420 LWTFMT_ELEMID, nid);
3421 return -1;
3422 }}
3423
3424 /* Bail out if next CW or CCW edge on end node changed */
3425 if ( epan_pre.nextCW != epan_post.nextCW ||
3426 epan_pre.nextCCW != epan_post.nextCCW )
3427 {{
3428 LWT_ELEMID nid = oldedge->end_node;
3429 _lwt_release_edges(oldedge, 1);
3430 lwerror("Edge changed disposition around end node %"
3431 LWTFMT_ELEMID, nid);
3432 return -1;
3433 }}
3434
3435 /*
3436 -- Update faces MBR of left and right faces
3437 -- TODO: think about ways to optimize this part, like see if
3438 -- the old edge geometry participated in the definition
3439 -- of the current MBR (for shrinking) or the new edge MBR
3440 -- would be larger than the old face MBR...
3441 --
3442 */
3443 const GBOX* oldbox = lwgeom_get_bbox(lwline_as_lwgeom(oldedge->geom));
3444 const GBOX* newbox = lwgeom_get_bbox(lwline_as_lwgeom(geom));
3445 if ( ! gbox_same(oldbox, newbox) )
3446 {
3447 uint64_t facestoupdate = 0;
3448 LWT_ISO_FACE faces[2];
3449 LWGEOM *nface1 = NULL;
3450 LWGEOM *nface2 = NULL;
3451 if ( oldedge->face_left > 0 )
3452 {
3453 nface1 = lwt_GetFaceGeometry(topo, oldedge->face_left);
3454 if ( ! nface1 )
3455 {
3456 lwerror("lwt_ChangeEdgeGeom could not construct face %"
3457 LWTFMT_ELEMID ", on the left of edge %" LWTFMT_ELEMID,
3458 oldedge->face_left, edge_id);
3459 return -1;
3460 }
3461 #if 0
3462 {
3463 size_t sz;
3464 char *wkt = lwgeom_to_wkt(nface1, WKT_EXTENDED, 2, &sz);
3465 LWDEBUGF(1, "new geometry of face left (%d): %s", (int)oldedge->face_left, wkt);
3466 lwfree(wkt);
3467 }
3468 #endif
3469 lwgeom_add_bbox(nface1);
3470 if ( ! nface1->bbox )
3471 {
3472 lwerror("Corrupted topology: face %d, left of edge %d, has no bbox",
3473 oldedge->face_left, edge_id);
3474 return -1;
3475 }
3476 faces[facestoupdate].face_id = oldedge->face_left;
3477 /* ownership left to nface */
3478 faces[facestoupdate++].mbr = nface1->bbox;
3479 }
3480 if ( oldedge->face_right > 0
3481 /* no need to update twice the same face.. */
3482 && oldedge->face_right != oldedge->face_left )
3483 {
3484 nface2 = lwt_GetFaceGeometry(topo, oldedge->face_right);
3485 if ( ! nface2 )
3486 {
3487 lwerror("lwt_ChangeEdgeGeom could not construct face %"
3488 LWTFMT_ELEMID ", on the right of edge %" LWTFMT_ELEMID,
3489 oldedge->face_right, edge_id);
3490 return -1;
3491 }
3492 #if 0
3493 {
3494 size_t sz;
3495 char *wkt = lwgeom_to_wkt(nface2, WKT_EXTENDED, 2, &sz);
3496 LWDEBUGF(1, "new geometry of face right (%d): %s", (int)oldedge->face_right, wkt);
3497 lwfree(wkt);
3498 }
3499 #endif
3500 lwgeom_add_bbox(nface2);
3501 if ( ! nface2->bbox )
3502 {
3503 lwerror("Corrupted topology: face %d, right of edge %d, has no bbox",
3504 oldedge->face_right, edge_id);
3505 return -1;
3506 }
3507 faces[facestoupdate].face_id = oldedge->face_right;
3508 faces[facestoupdate++].mbr = nface2->bbox; /* ownership left to nface */
3509 }
3510 LWDEBUGF(1, "%d faces to update", facestoupdate);
3511 if ( facestoupdate )
3512 {
3513 uint64_t updatedFaces = lwt_be_updateFacesById(topo, &(faces[0]), facestoupdate);
3514 if (updatedFaces != facestoupdate)
3515 {
3516 if (nface1)
3517 lwgeom_free(nface1);
3518 if (nface2)
3519 lwgeom_free(nface2);
3520 _lwt_release_edges(oldedge, 1);
3521 if (updatedFaces == UINT64_MAX)
3522 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3523 else
3524 lwerror("Unexpected error: %d faces found when expecting 1", i);
3525 return -1;
3526 }
3527 }
3528 if ( nface1 ) lwgeom_free(nface1);
3529 if ( nface2 ) lwgeom_free(nface2);
3530 }
3531 else
3532 {
3533 lwnotice("BBOX of changed edge did not change");
3534 }
3535
3536 LWDEBUG(1, "all done, cleaning up edges");
3537
3538 _lwt_release_edges(oldedge, 1);
3539 return 0; /* success */
3540 }
3541
3542 /* Only return CONTAINING_FACE in the node object */
3543 static LWT_ISO_NODE *
_lwt_GetIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid)3544 _lwt_GetIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid)
3545 {
3546 LWT_ISO_NODE *node;
3547 uint64_t n = 1;
3548
3549 node = lwt_be_getNodeById( topo, &nid, &n, LWT_COL_NODE_CONTAINING_FACE );
3550 if (n == UINT64_MAX)
3551 {
3552 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3553 return 0;
3554 }
3555 if ( n < 1 ) {
3556 lwerror("SQL/MM Spatial exception - non-existent node");
3557 return 0;
3558 }
3559 if ( node->containing_face == -1 )
3560 {
3561 lwfree(node);
3562 lwerror("SQL/MM Spatial exception - not isolated node");
3563 return 0;
3564 }
3565
3566 return node;
3567 }
3568
3569 int
lwt_MoveIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid,LWPOINT * pt)3570 lwt_MoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid, LWPOINT *pt)
3571 {
3572 LWT_ISO_NODE *node;
3573 int ret;
3574
3575 node = _lwt_GetIsoNode( topo, nid );
3576 if ( ! node ) return -1;
3577
3578 if ( lwt_be_ExistsCoincidentNode(topo, pt) )
3579 {
3580 lwfree(node);
3581 lwerror("SQL/MM Spatial exception - coincident node");
3582 return -1;
3583 }
3584
3585 if ( lwt_be_ExistsEdgeIntersectingPoint(topo, pt) )
3586 {
3587 lwfree(node);
3588 lwerror("SQL/MM Spatial exception - edge crosses node.");
3589 return -1;
3590 }
3591
3592 /* TODO: check that the new point is in the same containing face !
3593 * See https://trac.osgeo.org/postgis/ticket/3232
3594 */
3595
3596 node->node_id = nid;
3597 node->geom = pt;
3598 ret = lwt_be_updateNodesById(topo, node, 1,
3599 LWT_COL_NODE_GEOM);
3600 if ( ret == -1 ) {
3601 lwfree(node);
3602 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3603 return -1;
3604 }
3605
3606 lwfree(node);
3607 return 0;
3608 }
3609
3610 int
lwt_RemoveIsoNode(LWT_TOPOLOGY * topo,LWT_ELEMID nid)3611 lwt_RemoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid)
3612 {
3613 LWT_ISO_NODE *node;
3614 int n = 1;
3615
3616 node = _lwt_GetIsoNode( topo, nid );
3617 if ( ! node ) return -1;
3618
3619 n = lwt_be_deleteNodesById( topo, &nid, n );
3620 if ( n == -1 )
3621 {
3622 lwfree(node);
3623 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3624 return -1;
3625 }
3626 if ( n != 1 )
3627 {
3628 lwfree(node);
3629 lwerror("Unexpected error: %d nodes deleted when expecting 1", n);
3630 return -1;
3631 }
3632
3633 /* TODO: notify to caller about node being removed ?
3634 * See https://trac.osgeo.org/postgis/ticket/3231
3635 */
3636
3637 lwfree(node);
3638 return 0; /* success */
3639 }
3640
3641 int
lwt_RemIsoEdge(LWT_TOPOLOGY * topo,LWT_ELEMID id)3642 lwt_RemIsoEdge(LWT_TOPOLOGY* topo, LWT_ELEMID id)
3643 {
3644 LWT_ISO_EDGE deledge;
3645 LWT_ISO_EDGE *edge;
3646 LWT_ELEMID nid[2];
3647 LWT_ISO_NODE upd_node[2];
3648 LWT_ELEMID containing_face;
3649 uint64_t n = 1;
3650 uint64_t i;
3651
3652 edge = lwt_be_getEdgeById( topo, &id, &n, LWT_COL_EDGE_START_NODE|
3653 LWT_COL_EDGE_END_NODE |
3654 LWT_COL_EDGE_FACE_LEFT |
3655 LWT_COL_EDGE_FACE_RIGHT );
3656 if ( ! edge )
3657 {
3658 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3659 return -1;
3660 }
3661 if ( ! n )
3662 {
3663 lwerror("SQL/MM Spatial exception - non-existent edge");
3664 return -1;
3665 }
3666 if ( n > 1 )
3667 {
3668 lwfree(edge);
3669 lwerror("Corrupted topology: more than a single edge have id %"
3670 LWTFMT_ELEMID, id);
3671 return -1;
3672 }
3673
3674 if ( edge[0].face_left != edge[0].face_right )
3675 {
3676 lwfree(edge);
3677 lwerror("SQL/MM Spatial exception - not isolated edge");
3678 return -1;
3679 }
3680 containing_face = edge[0].face_left;
3681
3682 nid[0] = edge[0].start_node;
3683 nid[1] = edge[0].end_node;
3684 lwfree(edge);
3685
3686 n = 2;
3687 edge = lwt_be_getEdgeByNode( topo, nid, &n, LWT_COL_EDGE_EDGE_ID );
3688 if ((n == UINT64_MAX) || (edge == NULL))
3689 {
3690 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3691 return -1;
3692 }
3693 for (i = 0; i < n; ++i)
3694 {
3695 if (edge[i].edge_id != id)
3696 {
3697 lwfree(edge);
3698 lwerror("SQL/MM Spatial exception - not isolated edge");
3699 return -1;
3700 }
3701 }
3702 lwfree(edge);
3703
3704 deledge.edge_id = id;
3705 n = lwt_be_deleteEdges( topo, &deledge, LWT_COL_EDGE_EDGE_ID );
3706 if (n == UINT64_MAX)
3707 {
3708 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3709 return -1;
3710 }
3711 if ( n != 1 )
3712 {
3713 lwerror("Unexpected error: %d edges deleted when expecting 1", n);
3714 return -1;
3715 }
3716
3717 upd_node[0].node_id = nid[0];
3718 upd_node[0].containing_face = containing_face;
3719 n = 1;
3720 if ( nid[1] != nid[0] ) {
3721 upd_node[1].node_id = nid[1];
3722 upd_node[1].containing_face = containing_face;
3723 ++n;
3724 }
3725 n = lwt_be_updateNodesById(topo, upd_node, n,
3726 LWT_COL_NODE_CONTAINING_FACE);
3727 if (n == UINT64_MAX)
3728 {
3729 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3730 return -1;
3731 }
3732
3733 /* TODO: notify to caller about edge being removed ?
3734 * See https://trac.osgeo.org/postgis/ticket/3248
3735 */
3736
3737 return 0; /* success */
3738 }
3739
3740 /* Used by _lwt_RemEdge to update edge face ref on healing
3741 *
3742 * @param of old face id (never 0 as you cannot remove face 0)
3743 * @param nf new face id
3744 * @return 0 on success, -1 on backend error
3745 */
3746 static int
_lwt_UpdateEdgeFaceRef(LWT_TOPOLOGY * topo,LWT_ELEMID of,LWT_ELEMID nf)3747 _lwt_UpdateEdgeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf)
3748 {
3749 LWT_ISO_EDGE sel_edge, upd_edge;
3750 int ret;
3751
3752 assert( of != 0 );
3753
3754 /* Update face_left for all edges still referencing old face */
3755 sel_edge.face_left = of;
3756 upd_edge.face_left = nf;
3757 ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_LEFT,
3758 &upd_edge, LWT_COL_EDGE_FACE_LEFT,
3759 NULL, 0);
3760 if ( ret == -1 ) return -1;
3761
3762 /* Update face_right for all edges still referencing old face */
3763 sel_edge.face_right = of;
3764 upd_edge.face_right = nf;
3765 ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_RIGHT,
3766 &upd_edge, LWT_COL_EDGE_FACE_RIGHT,
3767 NULL, 0);
3768 if ( ret == -1 ) return -1;
3769
3770 return 0;
3771 }
3772
3773 /* Used by _lwt_RemEdge to update node face ref on healing
3774 *
3775 * @param of old face id (never 0 as you cannot remove face 0)
3776 * @param nf new face id
3777 * @return 0 on success, -1 on backend error
3778 */
3779 static int
_lwt_UpdateNodeFaceRef(LWT_TOPOLOGY * topo,LWT_ELEMID of,LWT_ELEMID nf)3780 _lwt_UpdateNodeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf)
3781 {
3782 LWT_ISO_NODE sel, upd;
3783 int ret;
3784
3785 assert( of != 0 );
3786
3787 /* Update face_left for all edges still referencing old face */
3788 sel.containing_face = of;
3789 upd.containing_face = nf;
3790 ret = lwt_be_updateNodes(topo, &sel, LWT_COL_NODE_CONTAINING_FACE,
3791 &upd, LWT_COL_NODE_CONTAINING_FACE,
3792 NULL, 0);
3793 if ( ret == -1 ) return -1;
3794
3795 return 0;
3796 }
3797
3798 /* Used by lwt_RemEdgeModFace and lwt_RemEdgeNewFaces
3799 *
3800 * Returns -1 on error, identifier of the face that takes up the space
3801 * previously occupied by the removed edge if modFace is 1, identifier of
3802 * the created face (0 if none) if modFace is 0.
3803 */
3804 static LWT_ELEMID
_lwt_RemEdge(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id,int modFace)3805 _lwt_RemEdge( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, int modFace )
3806 {
3807 uint64_t i, nedges, nfaces, fields;
3808 LWT_ISO_EDGE *edge = NULL;
3809 LWT_ISO_EDGE *upd_edge = NULL;
3810 LWT_ISO_EDGE upd_edge_left[2];
3811 int nedge_left = 0;
3812 LWT_ISO_EDGE upd_edge_right[2];
3813 int nedge_right = 0;
3814 LWT_ISO_NODE upd_node[2];
3815 int nnode = 0;
3816 LWT_ISO_FACE *faces = NULL;
3817 LWT_ISO_FACE newface;
3818 LWT_ELEMID node_ids[2];
3819 LWT_ELEMID face_ids[2];
3820 int fnode_edges = 0; /* number of edges on the first node (excluded
3821 * the one being removed ) */
3822 int lnode_edges = 0; /* number of edges on the last node (excluded
3823 * the one being removed ) */
3824
3825 newface.face_id = 0;
3826
3827 i = 1;
3828 edge = lwt_be_getEdgeById(topo, &edge_id, &i, LWT_COL_EDGE_ALL);
3829 if (!edge)
3830 {
3831 LWDEBUGF(1, "lwt_be_getEdgeById returned NULL and set i=%d", i);
3832 if (i == UINT64_MAX)
3833 {
3834 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3835 return -1;
3836 }
3837 else if (i == 0)
3838 {
3839 lwerror("SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID, edge_id);
3840 return -1;
3841 }
3842 else
3843 {
3844 lwerror(
3845 "Backend coding error: getEdgeById callback returned NULL "
3846 "but numelements output parameter has value %d "
3847 "(expected 0 or 1)",
3848 i);
3849 return -1;
3850 }
3851 }
3852
3853 if ( ! lwt_be_checkTopoGeomRemEdge(topo, edge_id,
3854 edge->face_left, edge->face_right) )
3855 {
3856 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
3857 return -1;
3858 }
3859
3860 LWDEBUG(1, "Updating next_{right,left}_face of ring edges...");
3861
3862 /* Update edge linking */
3863
3864 nedges = 0;
3865 node_ids[nedges++] = edge->start_node;
3866 if ( edge->end_node != edge->start_node )
3867 {
3868 node_ids[nedges++] = edge->end_node;
3869 }
3870 fields = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_START_NODE |
3871 LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_NEXT_LEFT |
3872 LWT_COL_EDGE_NEXT_RIGHT;
3873 upd_edge = lwt_be_getEdgeByNode( topo, &(node_ids[0]), &nedges, fields );
3874 if (nedges == UINT64_MAX)
3875 {
3876 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3877 return -1;
3878 }
3879 nedge_left = nedge_right = 0;
3880 for ( i=0; i<nedges; ++i )
3881 {
3882 LWT_ISO_EDGE *e = &(upd_edge[i]);
3883 if ( e->edge_id == edge_id ) continue;
3884 if ( e->start_node == edge->start_node || e->end_node == edge->start_node )
3885 {
3886 ++fnode_edges;
3887 }
3888 if ( e->start_node == edge->end_node || e->end_node == edge->end_node )
3889 {
3890 ++lnode_edges;
3891 }
3892 if ( e->next_left == -edge_id )
3893 {
3894 upd_edge_left[nedge_left].edge_id = e->edge_id;
3895 upd_edge_left[nedge_left++].next_left =
3896 edge->next_left != edge_id ? edge->next_left : edge->next_right;
3897 }
3898 else if ( e->next_left == edge_id )
3899 {
3900 upd_edge_left[nedge_left].edge_id = e->edge_id;
3901 upd_edge_left[nedge_left++].next_left =
3902 edge->next_right != -edge_id ? edge->next_right : edge->next_left;
3903 }
3904
3905 if ( e->next_right == -edge_id )
3906 {
3907 upd_edge_right[nedge_right].edge_id = e->edge_id;
3908 upd_edge_right[nedge_right++].next_right =
3909 edge->next_left != edge_id ? edge->next_left : edge->next_right;
3910 }
3911 else if ( e->next_right == edge_id )
3912 {
3913 upd_edge_right[nedge_right].edge_id = e->edge_id;
3914 upd_edge_right[nedge_right++].next_right =
3915 edge->next_right != -edge_id ? edge->next_right : edge->next_left;
3916 }
3917 }
3918
3919 if ( nedge_left )
3920 {
3921 LWDEBUGF(1, "updating %d 'next_left' edges", nedge_left);
3922 /* update edges in upd_edge_left set next_left */
3923 int result = lwt_be_updateEdgesById(topo, &(upd_edge_left[0]), nedge_left, LWT_COL_EDGE_NEXT_LEFT);
3924 if (result == -1)
3925 {
3926 _lwt_release_edges(edge, 1);
3927 lwfree(upd_edge);
3928 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3929 return -1;
3930 }
3931 }
3932 if ( nedge_right )
3933 {
3934 LWDEBUGF(1, "updating %d 'next_right' edges", nedge_right);
3935 /* update edges in upd_edge_right set next_right */
3936 int result = lwt_be_updateEdgesById(topo, &(upd_edge_right[0]), nedge_right, LWT_COL_EDGE_NEXT_RIGHT);
3937 if (result == -1)
3938 {
3939 _lwt_release_edges(edge, 1);
3940 lwfree(upd_edge);
3941 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3942 return -1;
3943 }
3944 }
3945 LWDEBUGF(1, "releasing %d updateable edges in %p", nedges, upd_edge);
3946 lwfree(upd_edge);
3947
3948 /* Id of face that will take up all the space previously
3949 * taken by left and right faces of the edge */
3950 LWT_ELEMID floodface;
3951
3952 /* Find floodface, and update its mbr if != 0 */
3953 if ( edge->face_left == edge->face_right )
3954 {
3955 floodface = edge->face_right;
3956 }
3957 else
3958 {
3959 /* Two faces healed */
3960 if ( edge->face_left == 0 || edge->face_right == 0 )
3961 {
3962 floodface = 0;
3963 LWDEBUG(1, "floodface is universe");
3964 }
3965 else
3966 {
3967 /* we choose right face as the face that will remain
3968 * to be symmetric with ST_AddEdgeModFace */
3969 floodface = edge->face_right;
3970 LWDEBUGF(1, "floodface is %" LWTFMT_ELEMID, floodface);
3971 /* update mbr of floodface as union of mbr of both faces */
3972 face_ids[0] = edge->face_left;
3973 face_ids[1] = edge->face_right;
3974 nfaces = 2;
3975 fields = LWT_COL_FACE_ALL;
3976 faces = lwt_be_getFaceById(topo, face_ids, &nfaces, fields);
3977 if (nfaces == UINT64_MAX)
3978 {
3979 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
3980 return -1;
3981 }
3982 GBOX *box1=NULL;
3983 GBOX *box2=NULL;
3984 for ( i=0; i<nfaces; ++i )
3985 {
3986 if ( faces[i].face_id == edge->face_left )
3987 {
3988 if ( ! box1 ) box1 = faces[i].mbr;
3989 else
3990 {
3991 i = edge->face_left;
3992 _lwt_release_edges(edge, 1);
3993 _lwt_release_faces(faces, nfaces);
3994 lwerror("corrupted topology: more than 1 face have face_id=%"
3995 LWTFMT_ELEMID, i);
3996 return -1;
3997 }
3998 }
3999 else if ( faces[i].face_id == edge->face_right )
4000 {
4001 if ( ! box2 ) box2 = faces[i].mbr;
4002 else
4003 {
4004 i = edge->face_right;
4005 _lwt_release_edges(edge, 1);
4006 _lwt_release_faces(faces, nfaces);
4007 lwerror("corrupted topology: more than 1 face have face_id=%"
4008 LWTFMT_ELEMID, i);
4009 return -1;
4010 }
4011 }
4012 else
4013 {
4014 i = faces[i].face_id;
4015 _lwt_release_edges(edge, 1);
4016 _lwt_release_faces(faces, nfaces);
4017 lwerror("Backend coding error: getFaceById returned face "
4018 "with non-requested id %" LWTFMT_ELEMID, i);
4019 return -1;
4020 }
4021 }
4022 if ( ! box1 ) {
4023 i = edge->face_left;
4024 _lwt_release_edges(edge, 1);
4025 _lwt_release_faces(faces, nfaces);
4026 lwerror("corrupted topology: no face have face_id=%"
4027 LWTFMT_ELEMID " (left face for edge %"
4028 LWTFMT_ELEMID ")", i, edge_id);
4029 return -1;
4030 }
4031 if ( ! box2 ) {
4032 i = edge->face_right;
4033 _lwt_release_edges(edge, 1);
4034 _lwt_release_faces(faces, nfaces);
4035 lwerror("corrupted topology: no face have face_id=%"
4036 LWTFMT_ELEMID " (right face for edge %"
4037 LWTFMT_ELEMID ")", i, edge_id);
4038 return -1;
4039 }
4040 gbox_merge(box2, box1); /* box1 is now the union of the two */
4041 newface.mbr = box1;
4042 if ( modFace )
4043 {
4044 newface.face_id = floodface;
4045 int result = lwt_be_updateFacesById(topo, &newface, 1);
4046 _lwt_release_faces(faces, 2);
4047 if (result == -1)
4048 {
4049 _lwt_release_edges(edge, 1);
4050 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4051 return -1;
4052 }
4053 if (result != 1)
4054 {
4055 _lwt_release_edges(edge, 1);
4056 lwerror("Unexpected error: %d faces updated when expecting 1", i);
4057 return -1;
4058 }
4059 }
4060 else
4061 {
4062 /* New face replaces the old two faces */
4063 newface.face_id = -1;
4064 int result = lwt_be_insertFaces(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 inserted when expecting 1", result);
4076 return -1;
4077 }
4078 floodface = newface.face_id;
4079 }
4080 }
4081
4082 /* Update face references for edges and nodes still referencing
4083 * the removed face(s) */
4084
4085 if ( edge->face_left != floodface )
4086 {
4087 if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_left, floodface) )
4088 {
4089 _lwt_release_edges(edge, 1);
4090 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4091 return -1;
4092 }
4093 if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_left, floodface) )
4094 {
4095 _lwt_release_edges(edge, 1);
4096 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4097 return -1;
4098 }
4099 }
4100
4101 if ( edge->face_right != floodface )
4102 {
4103 if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_right, floodface) )
4104 {
4105 _lwt_release_edges(edge, 1);
4106 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4107 return -1;
4108 }
4109 if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_right, floodface) )
4110 {
4111 _lwt_release_edges(edge, 1);
4112 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4113 return -1;
4114 }
4115 }
4116
4117 /* Update topogeoms on heal */
4118 if ( ! lwt_be_updateTopoGeomFaceHeal(topo,
4119 edge->face_right, edge->face_left,
4120 floodface) )
4121 {
4122 _lwt_release_edges(edge, 1);
4123 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4124 return -1;
4125 }
4126 } /* two faces healed */
4127
4128 /* Delete the edge */
4129 int result = lwt_be_deleteEdges(topo, edge, LWT_COL_EDGE_EDGE_ID);
4130 if (result == -1)
4131 {
4132 _lwt_release_edges(edge, 1);
4133 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4134 return -1;
4135 }
4136
4137 /* If any of the edge nodes remained isolated, set
4138 * containing_face = floodface
4139 */
4140 if ( ! fnode_edges )
4141 {
4142 upd_node[nnode].node_id = edge->start_node;
4143 upd_node[nnode].containing_face = floodface;
4144 ++nnode;
4145 }
4146 if ( edge->end_node != edge->start_node && ! lnode_edges )
4147 {
4148 upd_node[nnode].node_id = edge->end_node;
4149 upd_node[nnode].containing_face = floodface;
4150 ++nnode;
4151 }
4152 if ( nnode )
4153 {
4154 int result = lwt_be_updateNodesById(topo, upd_node, nnode, LWT_COL_NODE_CONTAINING_FACE);
4155 if (result == -1)
4156 {
4157 _lwt_release_edges(edge, 1);
4158 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4159 return -1;
4160 }
4161 }
4162
4163 if ( edge->face_left != edge->face_right )
4164 /* or there'd be no face to remove */
4165 {
4166 LWT_ELEMID ids[2];
4167 int nids = 0;
4168 if ( edge->face_right != floodface )
4169 ids[nids++] = edge->face_right;
4170 if ( edge->face_left != floodface )
4171 ids[nids++] = edge->face_left;
4172 int result = lwt_be_deleteFacesById(topo, ids, nids);
4173 if (result == -1)
4174 {
4175 _lwt_release_edges(edge, 1);
4176 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4177 return -1;
4178 }
4179 }
4180
4181 _lwt_release_edges(edge, 1);
4182 return modFace ? floodface : newface.face_id;
4183 }
4184
4185 LWT_ELEMID
lwt_RemEdgeModFace(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id)4186 lwt_RemEdgeModFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id )
4187 {
4188 return _lwt_RemEdge( topo, edge_id, 1 );
4189 }
4190
4191 LWT_ELEMID
lwt_RemEdgeNewFace(LWT_TOPOLOGY * topo,LWT_ELEMID edge_id)4192 lwt_RemEdgeNewFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id )
4193 {
4194 return _lwt_RemEdge( topo, edge_id, 0 );
4195 }
4196
4197 static LWT_ELEMID
_lwt_HealEdges(LWT_TOPOLOGY * topo,LWT_ELEMID eid1,LWT_ELEMID eid2,int modEdge)4198 _lwt_HealEdges( LWT_TOPOLOGY* topo, LWT_ELEMID eid1, LWT_ELEMID eid2,
4199 int modEdge )
4200 {
4201 LWT_ELEMID ids[2];
4202 LWT_ELEMID commonnode = -1;
4203 int caseno = 0;
4204 LWT_ISO_EDGE *node_edges;
4205 uint64_t num_node_edges;
4206 LWT_ISO_EDGE *edges;
4207 LWT_ISO_EDGE *e1 = NULL;
4208 LWT_ISO_EDGE *e2 = NULL;
4209 LWT_ISO_EDGE newedge, updedge, seledge;
4210 uint64_t nedges, i;
4211 int e1freenode;
4212 int e2sign, e2freenode;
4213 POINTARRAY *pa;
4214 char buf[256];
4215 char *ptr;
4216 size_t bufleft = 256;
4217
4218 ptr = buf;
4219
4220 /* NOT IN THE SPECS: see if the same edge is given twice.. */
4221 if ( eid1 == eid2 )
4222 {
4223 lwerror("Cannot heal edge %" LWTFMT_ELEMID
4224 " with itself, try with another", eid1);
4225 return -1;
4226 }
4227 ids[0] = eid1;
4228 ids[1] = eid2;
4229 nedges = 2;
4230 edges = lwt_be_getEdgeById(topo, ids, &nedges, LWT_COL_EDGE_ALL);
4231 if ((nedges == UINT64_MAX) || (edges == NULL))
4232 {
4233 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4234 return -1;
4235 }
4236 for ( i=0; i<nedges; ++i )
4237 {
4238 if ( edges[i].edge_id == eid1 ) {
4239 if ( e1 ) {
4240 _lwt_release_edges(edges, nedges);
4241 lwerror("Corrupted topology: multiple edges have id %"
4242 LWTFMT_ELEMID, eid1);
4243 return -1;
4244 }
4245 e1 = &(edges[i]);
4246 }
4247 else if ( edges[i].edge_id == eid2 ) {
4248 if ( e2 ) {
4249 _lwt_release_edges(edges, nedges);
4250 lwerror("Corrupted topology: multiple edges have id %"
4251 LWTFMT_ELEMID, eid2);
4252 return -1;
4253 }
4254 e2 = &(edges[i]);
4255 }
4256 }
4257 if ( ! e1 )
4258 {
4259 _lwt_release_edges(edges, nedges);
4260 lwerror(
4261 "SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID,
4262 eid1);
4263 return -1;
4264 }
4265 if ( ! e2 )
4266 {
4267 _lwt_release_edges(edges, nedges);
4268 lwerror(
4269 "SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID,
4270 eid2);
4271 return -1;
4272 }
4273
4274 /* NOT IN THE SPECS: See if any of the two edges are closed. */
4275 if ( e1->start_node == e1->end_node )
4276 {
4277 _lwt_release_edges(edges, nedges);
4278 lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %"
4279 LWTFMT_ELEMID, eid1, eid2);
4280 return -1;
4281 }
4282 if ( e2->start_node == e2->end_node )
4283 {
4284 _lwt_release_edges(edges, nedges);
4285 lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %"
4286 LWTFMT_ELEMID, eid2, eid1);
4287 return -1;
4288 }
4289
4290 /* Find common node */
4291
4292 if ( e1->end_node == e2->start_node )
4293 {
4294 commonnode = e1->end_node;
4295 caseno = 1;
4296 }
4297 else if ( e1->end_node == e2->end_node )
4298 {
4299 commonnode = e1->end_node;
4300 caseno = 2;
4301 }
4302 /* Check if any other edge is connected to the common node, if found */
4303 if ( commonnode != -1 )
4304 {
4305 num_node_edges = 1;
4306 node_edges = lwt_be_getEdgeByNode( topo, &commonnode,
4307 &num_node_edges, LWT_COL_EDGE_EDGE_ID );
4308 if (num_node_edges == UINT64_MAX)
4309 {
4310 _lwt_release_edges(edges, nedges);
4311 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4312 return -1;
4313 }
4314 for (i=0; i<num_node_edges; ++i)
4315 {
4316 int r;
4317 if ( node_edges[i].edge_id == eid1 ) continue;
4318 if ( node_edges[i].edge_id == eid2 ) continue;
4319 commonnode = -1;
4320 /* append to string, for error message */
4321 if ( bufleft > 0 ) {
4322 r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID,
4323 ( ptr==buf ? "" : "," ), node_edges[i].edge_id);
4324 if ( r >= (int) bufleft )
4325 {
4326 bufleft = 0;
4327 buf[252] = '.';
4328 buf[253] = '.';
4329 buf[254] = '.';
4330 buf[255] = '\0';
4331 }
4332 else
4333 {
4334 bufleft -= r;
4335 ptr += r;
4336 }
4337 }
4338 }
4339 lwfree(node_edges);
4340 }
4341
4342 if ( commonnode == -1 )
4343 {
4344 if ( e1->start_node == e2->start_node )
4345 {
4346 commonnode = e1->start_node;
4347 caseno = 3;
4348 }
4349 else if ( e1->start_node == e2->end_node )
4350 {
4351 commonnode = e1->start_node;
4352 caseno = 4;
4353 }
4354 /* Check if any other edge is connected to the common node, if found */
4355 if ( commonnode != -1 )
4356 {
4357 num_node_edges = 1;
4358 node_edges = lwt_be_getEdgeByNode( topo, &commonnode,
4359 &num_node_edges, LWT_COL_EDGE_EDGE_ID );
4360 if (num_node_edges == UINT64_MAX)
4361 {
4362 _lwt_release_edges(edges, nedges);
4363 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4364 return -1;
4365 }
4366 for (i=0; i<num_node_edges; ++i)
4367 {
4368 int r;
4369 if ( node_edges[i].edge_id == eid1 ) continue;
4370 if ( node_edges[i].edge_id == eid2 ) continue;
4371 commonnode = -1;
4372 /* append to string, for error message */
4373 if ( bufleft > 0 ) {
4374 r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID,
4375 ( ptr==buf ? "" : "," ), node_edges[i].edge_id);
4376 if ( r >= (int) bufleft )
4377 {
4378 bufleft = 0;
4379 buf[252] = '.';
4380 buf[253] = '.';
4381 buf[254] = '.';
4382 buf[255] = '\0';
4383 }
4384 else
4385 {
4386 bufleft -= r;
4387 ptr += r;
4388 }
4389 }
4390 }
4391 if ( num_node_edges ) lwfree(node_edges);
4392 }
4393 }
4394
4395 if ( commonnode == -1 )
4396 {
4397 _lwt_release_edges(edges, nedges);
4398 if ( ptr != buf )
4399 {
4400 lwerror("SQL/MM Spatial exception - other edges connected (%s)",
4401 buf);
4402 }
4403 else
4404 {
4405 lwerror("SQL/MM Spatial exception - non-connected edges");
4406 }
4407 return -1;
4408 }
4409
4410 if ( ! lwt_be_checkTopoGeomRemNode(topo, commonnode,
4411 eid1, eid2 ) )
4412 {
4413 _lwt_release_edges(edges, nedges);
4414 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4415 return -1;
4416 }
4417
4418 /* Construct the geometry of the new edge */
4419 switch (caseno)
4420 {
4421 case 1: /* e1.end = e2.start */
4422 pa = ptarray_clone_deep(e1->geom->points);
4423 //pa = ptarray_merge(pa, e2->geom->points);
4424 ptarray_append_ptarray(pa, e2->geom->points, 0);
4425 newedge.start_node = e1->start_node;
4426 newedge.end_node = e2->end_node;
4427 newedge.next_left = e2->next_left;
4428 newedge.next_right = e1->next_right;
4429 e1freenode = 1;
4430 e2freenode = -1;
4431 e2sign = 1;
4432 break;
4433 case 2: /* e1.end = e2.end */
4434 {
4435 POINTARRAY *pa2;
4436 pa2 = ptarray_clone_deep(e2->geom->points);
4437 ptarray_reverse_in_place(pa2);
4438 pa = ptarray_clone_deep(e1->geom->points);
4439 //pa = ptarray_merge(e1->geom->points, pa);
4440 ptarray_append_ptarray(pa, pa2, 0);
4441 ptarray_free(pa2);
4442 newedge.start_node = e1->start_node;
4443 newedge.end_node = e2->start_node;
4444 newedge.next_left = e2->next_right;
4445 newedge.next_right = e1->next_right;
4446 e1freenode = 1;
4447 e2freenode = 1;
4448 e2sign = -1;
4449 break;
4450 }
4451 case 3: /* e1.start = e2.start */
4452 pa = ptarray_clone_deep(e2->geom->points);
4453 ptarray_reverse_in_place(pa);
4454 //pa = ptarray_merge(pa, e1->geom->points);
4455 ptarray_append_ptarray(pa, e1->geom->points, 0);
4456 newedge.end_node = e1->end_node;
4457 newedge.start_node = e2->end_node;
4458 newedge.next_left = e1->next_left;
4459 newedge.next_right = e2->next_left;
4460 e1freenode = -1;
4461 e2freenode = -1;
4462 e2sign = -1;
4463 break;
4464 case 4: /* e1.start = e2.end */
4465 pa = ptarray_clone_deep(e2->geom->points);
4466 //pa = ptarray_merge(pa, e1->geom->points);
4467 ptarray_append_ptarray(pa, e1->geom->points, 0);
4468 newedge.end_node = e1->end_node;
4469 newedge.start_node = e2->start_node;
4470 newedge.next_left = e1->next_left;
4471 newedge.next_right = e2->next_right;
4472 e1freenode = -1;
4473 e2freenode = 1;
4474 e2sign = 1;
4475 break;
4476 default:
4477 pa = NULL;
4478 e1freenode = 0;
4479 e2freenode = 0;
4480 e2sign = 0;
4481 _lwt_release_edges(edges, nedges);
4482 lwerror("Coding error: caseno=%d should never happen", caseno);
4483 break;
4484 }
4485 newedge.geom = lwline_construct(topo->srid, NULL, pa);
4486
4487 if ( modEdge )
4488 {
4489 /* Update data of the first edge */
4490 newedge.edge_id = eid1;
4491 int result = lwt_be_updateEdgesById(topo,
4492 &newedge,
4493 1,
4494 LWT_COL_EDGE_NEXT_LEFT | LWT_COL_EDGE_NEXT_RIGHT | LWT_COL_EDGE_START_NODE |
4495 LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_GEOM);
4496 if (result == -1)
4497 {
4498 lwline_free(newedge.geom);
4499 _lwt_release_edges(edges, nedges);
4500 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4501 return -1;
4502 }
4503 else if (result != 1)
4504 {
4505 lwline_free(newedge.geom);
4506 if ( edges ) _lwt_release_edges(edges, nedges);
4507 lwerror("Unexpected error: %d edges updated when expecting 1", i);
4508 return -1;
4509 }
4510 }
4511 else
4512 {
4513 /* Add new edge */
4514 newedge.edge_id = -1;
4515 newedge.face_left = e1->face_left;
4516 newedge.face_right = e1->face_right;
4517 int result = lwt_be_insertEdges(topo, &newedge, 1);
4518 if (result == -1)
4519 {
4520 lwline_free(newedge.geom);
4521 _lwt_release_edges(edges, nedges);
4522 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4523 return -1;
4524 }
4525 else if (result == 0)
4526 {
4527 lwline_free(newedge.geom);
4528 _lwt_release_edges(edges, nedges);
4529 lwerror("Insertion of split edge failed (no reason)");
4530 return -1;
4531 }
4532 }
4533 lwline_free(newedge.geom);
4534
4535 /*
4536 -- Update next_left_edge/next_right_edge for
4537 -- any edge having them still pointing at the edge being removed
4538 -- (eid2 only when modEdge, or both otherwise)
4539 --
4540 -- NOTE:
4541 -- e#freenode is 1 when edge# end node was the common node
4542 -- and -1 otherwise. This gives the sign of possibly found references
4543 -- to its "free" (non connected to other edge) endnode.
4544 -- e2sign is -1 if edge1 direction is opposite to edge2 direction,
4545 -- or 1 otherwise.
4546 --
4547 */
4548
4549 /* update edges connected to e2's boundary from their end node */
4550 seledge.next_left = e2freenode * eid2;
4551 updedge.next_left = e2freenode * newedge.edge_id * e2sign;
4552 int result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0);
4553 if (result == -1)
4554 {
4555 _lwt_release_edges(edges, nedges);
4556 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4557 return -1;
4558 }
4559
4560 /* update edges connected to e2's boundary from their start node */
4561 seledge.next_right = e2freenode * eid2;
4562 updedge.next_right = e2freenode * newedge.edge_id * e2sign;
4563 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0);
4564 if (result == -1)
4565 {
4566 _lwt_release_edges(edges, nedges);
4567 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4568 return -1;
4569 }
4570
4571 if ( ! modEdge )
4572 {
4573 /* update edges connected to e1's boundary from their end node */
4574 seledge.next_left = e1freenode * eid1;
4575 updedge.next_left = e1freenode * newedge.edge_id;
4576 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0);
4577 if (result == -1)
4578 {
4579 _lwt_release_edges(edges, nedges);
4580 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4581 return -1;
4582 }
4583
4584 /* update edges connected to e1's boundary from their start node */
4585 seledge.next_right = e1freenode * eid1;
4586 updedge.next_right = e1freenode * newedge.edge_id;
4587 result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0);
4588 if (result == -1)
4589 {
4590 _lwt_release_edges(edges, nedges);
4591 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4592 return -1;
4593 }
4594 }
4595
4596 /* delete the edges (only second on modEdge or both) */
4597 result = lwt_be_deleteEdges(topo, e2, LWT_COL_EDGE_EDGE_ID);
4598 if (result == -1)
4599 {
4600 _lwt_release_edges(edges, nedges);
4601 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4602 return -1;
4603 }
4604 if ( ! modEdge ) {
4605 i = lwt_be_deleteEdges(topo, e1, LWT_COL_EDGE_EDGE_ID);
4606 if (result == -1)
4607 {
4608 _lwt_release_edges(edges, nedges);
4609 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4610 return -1;
4611 }
4612 }
4613
4614 _lwt_release_edges(edges, nedges);
4615
4616 /* delete the common node */
4617 i = lwt_be_deleteNodesById( topo, &commonnode, 1 );
4618 if (result == -1)
4619 {
4620 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4621 return -1;
4622 }
4623
4624 /*
4625 --
4626 -- NOT IN THE SPECS:
4627 -- Drop composition rows involving second
4628 -- edge, as the first edge took its space,
4629 -- and all affected TopoGeom have been previously checked
4630 -- for being composed by both edges.
4631 */
4632 if ( ! lwt_be_updateTopoGeomEdgeHeal(topo,
4633 eid1, eid2, newedge.edge_id) )
4634 {
4635 lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface));
4636 return -1;
4637 }
4638
4639 return modEdge ? commonnode : newedge.edge_id;
4640 }
4641
4642 LWT_ELEMID
lwt_ModEdgeHeal(LWT_TOPOLOGY * topo,LWT_ELEMID e1,LWT_ELEMID e2)4643 lwt_ModEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 )
4644 {
4645 return _lwt_HealEdges( topo, e1, e2, 1 );
4646 }
4647
4648 LWT_ELEMID
lwt_NewEdgeHeal(LWT_TOPOLOGY * topo,LWT_ELEMID e1,LWT_ELEMID e2)4649 lwt_NewEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 )
4650 {
4651 return _lwt_HealEdges( topo, e1, e2, 0 );
4652 }
4653
4654 LWT_ELEMID
lwt_GetNodeByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4655 lwt_GetNodeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4656 {
4657 LWT_ISO_NODE *elem;
4658 uint64_t num;
4659 int flds = LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM; /* geom not needed */
4660 LWT_ELEMID id = 0;
4661 POINT2D qp; /* query point */
4662
4663 if ( ! getPoint2d_p(pt->point, 0, &qp) )
4664 {
4665 lwerror("Empty query point");
4666 return -1;
4667 }
4668 elem = lwt_be_getNodeWithinDistance2D(topo, pt, tol, &num, flds, 0);
4669 if (num == UINT64_MAX)
4670 {
4671 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4672 return -1;
4673 }
4674 else if ( num )
4675 {
4676 if ( num > 1 )
4677 {
4678 _lwt_release_nodes(elem, num);
4679 lwerror("Two or more nodes found");
4680 return -1;
4681 }
4682 id = elem[0].node_id;
4683 _lwt_release_nodes(elem, num);
4684 }
4685
4686 return id;
4687 }
4688
4689 LWT_ELEMID
lwt_GetEdgeByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4690 lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4691 {
4692 LWT_ISO_EDGE *elem;
4693 uint64_t num, i;
4694 int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; /* GEOM is not needed */
4695 LWT_ELEMID id = 0;
4696 LWGEOM *qp = lwpoint_as_lwgeom(pt); /* query point */
4697
4698 if ( lwgeom_is_empty(qp) )
4699 {
4700 lwerror("Empty query point");
4701 return -1;
4702 }
4703 elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol, &num, flds, 0);
4704 if (num == UINT64_MAX)
4705 {
4706 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4707 return -1;
4708 }
4709 for (i=0; i<num;++i)
4710 {
4711 LWT_ISO_EDGE *e = &(elem[i]);
4712 #if 0
4713 LWGEOM* geom;
4714 double dist;
4715
4716 if ( ! e->geom )
4717 {
4718 _lwt_release_edges(elem, num);
4719 lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
4720 " has null geometry", e->edge_id);
4721 continue;
4722 }
4723
4724 /* Should we check for intersection not being on an endpoint
4725 * as documented ? */
4726 geom = lwline_as_lwgeom(e->geom);
4727 dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
4728 if ( dist > tol ) continue;
4729 #endif
4730
4731 if ( id )
4732 {
4733 _lwt_release_edges(elem, num);
4734 lwerror("Two or more edges found");
4735 return -1;
4736 }
4737 else id = e->edge_id;
4738 }
4739
4740 if ( num ) _lwt_release_edges(elem, num);
4741
4742 return id;
4743 }
4744
4745 LWT_ELEMID
lwt_GetFaceByPoint(LWT_TOPOLOGY * topo,LWPOINT * pt,double tol)4746 lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
4747 {
4748 LWT_ELEMID id = 0;
4749 LWT_ISO_EDGE *elem;
4750 uint64_t num, i;
4751 int flds = LWT_COL_EDGE_EDGE_ID |
4752 LWT_COL_EDGE_GEOM |
4753 LWT_COL_EDGE_FACE_LEFT |
4754 LWT_COL_EDGE_FACE_RIGHT;
4755 LWGEOM *qp = lwpoint_as_lwgeom(pt);
4756
4757 id = lwt_be_getFaceContainingPoint(topo, pt);
4758 if ( id == -2 ) {
4759 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4760 return -1;
4761 }
4762
4763 if ( id > 0 )
4764 {
4765 return id;
4766 }
4767 id = 0; /* or it'll be -1 for not found */
4768
4769 LWDEBUG(1, "No face properly contains query point,"
4770 " looking for edges");
4771
4772 /* Not in a face, may be in universe or on edge, let's check
4773 * for distance */
4774 /* NOTE: we never pass a tolerance of 0 to avoid ever using
4775 * ST_Within, which doesn't include endpoints matches */
4776 elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol?tol:1e-5, &num, flds, 0);
4777 if (num == UINT64_MAX)
4778 {
4779 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4780 return -1;
4781 }
4782 for (i=0; i<num; ++i)
4783 {
4784 LWT_ISO_EDGE *e = &(elem[i]);
4785 LWT_ELEMID eface = 0;
4786 LWGEOM* geom;
4787 double dist;
4788
4789 if ( ! e->geom )
4790 {
4791 _lwt_release_edges(elem, num);
4792 lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID
4793 " has null geometry", e->edge_id);
4794 continue;
4795 }
4796
4797 /* don't consider dangling edges */
4798 if ( e->face_left == e->face_right )
4799 {
4800 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
4801 " is dangling, won't consider it", e->edge_id);
4802 continue;
4803 }
4804
4805 geom = lwline_as_lwgeom(e->geom);
4806 dist = lwgeom_mindistance2d_tolerance(geom, qp, tol);
4807
4808 LWDEBUGF(1, "Distance from edge %" LWTFMT_ELEMID
4809 " is %g (tol=%g)", e->edge_id, dist, tol);
4810
4811 /* we won't consider edges too far */
4812 if ( dist > tol ) continue;
4813 if ( e->face_left == 0 ) {
4814 eface = e->face_right;
4815 }
4816 else if ( e->face_right == 0 ) {
4817 eface = e->face_left;
4818 }
4819 else {
4820 _lwt_release_edges(elem, num);
4821 lwerror("Two or more faces found");
4822 return -1;
4823 }
4824
4825 if ( id && id != eface )
4826 {
4827 _lwt_release_edges(elem, num);
4828 lwerror("Two or more faces found"
4829 #if 0 /* debugging */
4830 " (%" LWTFMT_ELEMID
4831 " and %" LWTFMT_ELEMID ")", id, eface
4832 #endif
4833 );
4834 return -1;
4835 }
4836 else id = eface;
4837 }
4838 if ( num ) _lwt_release_edges(elem, num);
4839
4840 return id;
4841 }
4842
4843 /* Return the smallest delta that can perturbate
4844 * the given value */
4845 static inline double
_lwt_minToleranceDouble(double d)4846 _lwt_minToleranceDouble( double d )
4847 {
4848 double ret = 3.6 * pow(10, - ( 15 - log10(d?d:1.0) ) );
4849 return ret;
4850 }
4851
4852 /* Return the smallest delta that can perturbate
4853 * the given point
4854 static inline double
4855 _lwt_minTolerancePoint2d( const POINT2D* p )
4856 {
4857 double max = FP_ABS(p->x);
4858 if ( max < FP_ABS(p->y) ) max = FP_ABS(p->y);
4859 return _lwt_minToleranceDouble(max);
4860 }
4861 */
4862
4863 /* Return the smallest delta that can perturbate
4864 * the maximum absolute value of a geometry ordinate
4865 */
4866 static double
_lwt_minTolerance(LWGEOM * g)4867 _lwt_minTolerance( LWGEOM *g )
4868 {
4869 const GBOX* gbox;
4870 double max;
4871 double ret;
4872
4873 gbox = lwgeom_get_bbox(g);
4874 if ( ! gbox ) return 0; /* empty */
4875 max = FP_ABS(gbox->xmin);
4876 if ( max < FP_ABS(gbox->xmax) ) max = FP_ABS(gbox->xmax);
4877 if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin);
4878 if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax);
4879
4880 ret = _lwt_minToleranceDouble(max);
4881
4882 return ret;
4883 }
4884
4885 #define _LWT_MINTOLERANCE( topo, geom ) ( \
4886 topo->precision ? topo->precision : _lwt_minTolerance(geom) )
4887
4888 typedef struct scored_pointer_t {
4889 void *ptr;
4890 double score;
4891 } scored_pointer;
4892
4893 static int
compare_scored_pointer(const void * si1,const void * si2)4894 compare_scored_pointer(const void *si1, const void *si2)
4895 {
4896 double a = ((scored_pointer *)si1)->score;
4897 double b = ((scored_pointer *)si2)->score;
4898 if ( a < b )
4899 return -1;
4900 else if ( a > b )
4901 return 1;
4902 else
4903 return 0;
4904 }
4905
4906 /*
4907 * @param findFace if non-zero the code will determine which face
4908 * contains the given point (unless it is known to be NOT
4909 * isolated)
4910 * @param moved if not-null will be set to 0 if the point was added
4911 * w/out any snapping or 1 otherwise.
4912 */
4913 static LWT_ELEMID
_lwt_AddPoint(LWT_TOPOLOGY * topo,LWPOINT * point,double tol,int findFace,int * moved)4914 _lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol, int
4915 findFace, int *moved)
4916 {
4917 uint64_t num, i;
4918 double mindist = FLT_MAX;
4919 LWT_ISO_NODE *nodes, *nodes2;
4920 LWT_ISO_EDGE *edges, *edges2;
4921 LWGEOM *pt = lwpoint_as_lwgeom(point);
4922 int flds;
4923 LWT_ELEMID id = 0;
4924 scored_pointer *sorted;
4925
4926 /* Get tolerance, if 0 was given */
4927 if (!tol)
4928 tol = _LWT_MINTOLERANCE(topo, pt);
4929
4930 LWDEBUGG(1, pt, "Adding point");
4931
4932 /*
4933 -- 1. Check if any existing node is closer than the given precision
4934 -- and if so pick the closest
4935 TODO: use WithinBox2D
4936 */
4937 flds = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM;
4938 nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
4939 if (num == UINT64_MAX)
4940 {
4941 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4942 return -1;
4943 }
4944 if ( num )
4945 {
4946 LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
4947 /* Order by distance if there are more than a single return */
4948 if ( num > 1 )
4949 {{
4950 sorted= lwalloc(sizeof(scored_pointer)*num);
4951 for (i=0; i<num; ++i)
4952 {
4953 sorted[i].ptr = nodes+i;
4954 sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
4955 LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
4956 ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
4957 }
4958 qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
4959 nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
4960 for (i=0; i<num; ++i)
4961 {
4962 nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
4963 }
4964 lwfree(sorted);
4965 lwfree(nodes);
4966 nodes = nodes2;
4967 }}
4968
4969 for ( i=0; i<num; ++i )
4970 {
4971 LWT_ISO_NODE *n = &(nodes[i]);
4972 LWGEOM *g = lwpoint_as_lwgeom(n->geom);
4973 double dist = lwgeom_mindistance2d(g, pt);
4974 /* TODO: move this check in the previous sort scan ... */
4975 /* must be closer than tolerated, unless distance is zero */
4976 if ( dist && dist >= tol ) continue;
4977 if ( ! id || dist < mindist )
4978 {
4979 id = n->node_id;
4980 mindist = dist;
4981 }
4982 }
4983 if ( id )
4984 {
4985 /* found an existing node */
4986 if ( nodes ) _lwt_release_nodes(nodes, num);
4987 if ( moved ) *moved = mindist == 0 ? 0 : 1;
4988 return id;
4989 }
4990 }
4991
4992 initGEOS(lwnotice, lwgeom_geos_error);
4993
4994 /*
4995 -- 2. Check if any existing edge falls within tolerance
4996 -- and if so split it by a point projected on it
4997 TODO: use WithinBox2D
4998 */
4999 flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
5000 edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
5001 if (num == UINT64_MAX)
5002 {
5003 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5004 return -1;
5005 }
5006 if ( num )
5007 {
5008 LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
5009
5010 /* Order by distance if there are more than a single return */
5011 if ( num > 1 )
5012 {{
5013 int j;
5014 sorted = lwalloc(sizeof(scored_pointer)*num);
5015 for (i=0; i<num; ++i)
5016 {
5017 sorted[i].ptr = edges+i;
5018 sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
5019 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
5020 ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
5021 }
5022 qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5023 edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
5024 for (j=0, i=0; i<num; ++i)
5025 {
5026 if ( sorted[i].score == sorted[0].score )
5027 {
5028 edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
5029 }
5030 else
5031 {
5032 lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
5033 }
5034 }
5035 num = j;
5036 lwfree(sorted);
5037 lwfree(edges);
5038 edges = edges2;
5039 }}
5040
5041 for (i=0; i<num; ++i)
5042 {
5043 /* The point is on or near an edge, split the edge */
5044 LWT_ISO_EDGE *e = &(edges[i]);
5045 LWGEOM *g = lwline_as_lwgeom(e->geom);
5046 LWGEOM *prj;
5047 int contains;
5048 LWT_ELEMID edge_id = e->edge_id;
5049
5050 LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
5051
5052 /* project point to line, split edge by point */
5053 prj = lwgeom_closest_point(g, pt);
5054 if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1;
5055 if ( lwgeom_has_z(pt) )
5056 {{
5057 /*
5058 -- This is a workaround for ClosestPoint lack of Z support:
5059 -- http://trac.osgeo.org/postgis/ticket/2033
5060 */
5061 LWGEOM *tmp;
5062 double z;
5063 POINT4D p4d;
5064 LWPOINT *prjpt;
5065 /* add Z to "prj" */
5066 tmp = lwgeom_force_3dz(prj);
5067 prjpt = lwgeom_as_lwpoint(tmp);
5068 getPoint4d_p(point->point, 0, &p4d);
5069 z = p4d.z;
5070 getPoint4d_p(prjpt->point, 0, &p4d);
5071 p4d.z = z;
5072 ptarray_set_point4d(prjpt->point, 0, &p4d);
5073 lwgeom_free(prj);
5074 prj = tmp;
5075 }}
5076 const POINT2D *pt = getPoint2d_cp(lwgeom_as_lwpoint(prj)->point, 0);
5077 contains = ptarray_contains_point_partial(e->geom->points, pt, 0, NULL) == LW_BOUNDARY;
5078 if ( ! contains )
5079 {{
5080 double snaptol;
5081 LWGEOM *snapedge;
5082 LWLINE *snapline;
5083 POINT4D p1, p2;
5084
5085 LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
5086 " does not contain projected point to it",
5087 edge_id);
5088
5089 /* In order to reduce the robustness issues, we'll pick
5090 * an edge that contains the projected point, if possible */
5091 if ( i+1 < num )
5092 {
5093 LWDEBUG(1, "But there's another to check");
5094 lwgeom_free(prj);
5095 continue;
5096 }
5097
5098 /*
5099 -- The tolerance must be big enough for snapping to happen
5100 -- and small enough to snap only to the projected point.
5101 -- Unfortunately ST_Distance returns 0 because it also uses
5102 -- a projected point internally, so we need another way.
5103 */
5104 snaptol = _lwt_minTolerance(prj);
5105 snapedge = _lwt_toposnap(g, prj, snaptol);
5106 snapline = lwgeom_as_lwline(snapedge);
5107
5108 LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
5109
5110 /* TODO: check if snapping did anything ? */
5111 #if POSTGIS_DEBUG_LEVEL > 0
5112 {
5113 size_t sz;
5114 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5115 char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
5116 LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
5117 lwfree(wkt1);
5118 lwfree(wkt2);
5119 }
5120 #endif
5121
5122
5123 /*
5124 -- Snapping currently snaps the first point below tolerance
5125 -- so may possibly move first point. See ticket #1631
5126 */
5127 getPoint4d_p(e->geom->points, 0, &p1);
5128 getPoint4d_p(snapline->points, 0, &p2);
5129 LWDEBUGF(1, "Edge first point is %g %g, "
5130 "snapline first point is %g %g",
5131 p1.x, p1.y, p2.x, p2.y);
5132 if ( p1.x != p2.x || p1.y != p2.y )
5133 {
5134 LWDEBUG(1, "Snapping moved first point, re-adding it");
5135 if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
5136 {
5137 lwgeom_free(prj);
5138 lwgeom_free(snapedge);
5139 _lwt_release_edges(edges, num);
5140 lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5141 return -1;
5142 }
5143 #if POSTGIS_DEBUG_LEVEL > 0
5144 {
5145 size_t sz;
5146 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5147 LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
5148 lwfree(wkt1);
5149 }
5150 #endif
5151 }
5152 #if POSTGIS_DEBUG_LEVEL > 0
5153 else {
5154 LWDEBUG(1, "Snapping did not move first point");
5155 }
5156 #endif
5157
5158 if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
5159 {
5160 /* TODO: should have invoked lwerror already, leaking memory */
5161 lwgeom_free(prj);
5162 lwgeom_free(snapedge);
5163 _lwt_release_edges(edges, num);
5164 lwerror("lwt_ChangeEdgeGeom failed");
5165 return -1;
5166 }
5167 lwgeom_free(snapedge);
5168 }}
5169 #if POSTGIS_DEBUG_LEVEL > 0
5170 else
5171 {{
5172 size_t sz;
5173 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5174 char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
5175 LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
5176 lwfree(wkt1);
5177 lwfree(wkt2);
5178 }}
5179 #endif
5180
5181 /* TODO: pass 1 as last argument (skipChecks) ? */
5182 id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
5183 if ( -1 == id )
5184 {
5185 /* TODO: should have invoked lwerror already, leaking memory */
5186 lwgeom_free(prj);
5187 _lwt_release_edges(edges, num);
5188 lwerror("lwt_ModEdgeSplit failed");
5189 return -1;
5190 }
5191
5192 lwgeom_free(prj);
5193
5194 /*
5195 * TODO: decimate the two new edges with the given tolerance ?
5196 *
5197 * the edge identifiers to decimate would be: edge_id and "id"
5198 * The problem here is that decimation of existing edges
5199 * may introduce intersections or topological inconsistencies,
5200 * for example:
5201 *
5202 * - A node may end up falling on the other side of the edge
5203 * - The decimated edge might intersect another existing edge
5204 *
5205 */
5206
5207 break; /* we only want to snap a single edge */
5208 }
5209 _lwt_release_edges(edges, num);
5210 }
5211 else
5212 {
5213 /* The point is isolated, add it as such */
5214 /* TODO: pass 1 as last argument (skipChecks) ? */
5215 id = _lwt_AddIsoNode(topo, -1, point, 0, findFace);
5216 if ( moved ) *moved = 0;
5217 if ( -1 == id )
5218 {
5219 /* should have invoked lwerror already, leaking memory */
5220 lwerror("lwt_AddIsoNode failed");
5221 return -1;
5222 }
5223 }
5224
5225 return id;
5226 }
5227
5228 LWT_ELEMID
lwt_AddPoint(LWT_TOPOLOGY * topo,LWPOINT * point,double tol)5229 lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol)
5230 {
5231 return _lwt_AddPoint(topo, point, tol, 1, NULL);
5232 }
5233
5234 /* Return identifier of an equal edge, 0 if none or -1 on error
5235 * (and lwerror gets called on error)
5236 */
5237 static LWT_ELEMID
_lwt_GetEqualEdge(LWT_TOPOLOGY * topo,LWLINE * edge)5238 _lwt_GetEqualEdge( LWT_TOPOLOGY *topo, LWLINE *edge )
5239 {
5240 LWT_ELEMID id;
5241 LWT_ISO_EDGE *edges;
5242 uint64_t num, i;
5243 const GBOX *qbox = lwgeom_get_bbox( lwline_as_lwgeom(edge) );
5244 GEOSGeometry *edgeg;
5245 const int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
5246
5247 edges = lwt_be_getEdgeWithinBox2D( topo, qbox, &num, flds, 0 );
5248 if (num == UINT64_MAX)
5249 {
5250 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5251 return -1;
5252 }
5253 if ( num )
5254 {
5255 initGEOS(lwnotice, lwgeom_geos_error);
5256
5257 edgeg = LWGEOM2GEOS( lwline_as_lwgeom(edge), 0 );
5258 if ( ! edgeg )
5259 {
5260 _lwt_release_edges(edges, num);
5261 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5262 return -1;
5263 }
5264 for (i=0; i<num; ++i)
5265 {
5266 LWT_ISO_EDGE *e = &(edges[i]);
5267 LWGEOM *g = lwline_as_lwgeom(e->geom);
5268 GEOSGeometry *gg;
5269 int equals;
5270 gg = LWGEOM2GEOS( g, 0 );
5271 if ( ! gg )
5272 {
5273 GEOSGeom_destroy(edgeg);
5274 _lwt_release_edges(edges, num);
5275 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5276 return -1;
5277 }
5278 equals = GEOSEquals(gg, edgeg);
5279 GEOSGeom_destroy(gg);
5280 if ( equals == 2 )
5281 {
5282 GEOSGeom_destroy(edgeg);
5283 _lwt_release_edges(edges, num);
5284 lwerror("GEOSEquals exception: %s", lwgeom_geos_errmsg);
5285 return -1;
5286 }
5287 if ( equals )
5288 {
5289 id = e->edge_id;
5290 GEOSGeom_destroy(edgeg);
5291 _lwt_release_edges(edges, num);
5292 return id;
5293 }
5294 }
5295 GEOSGeom_destroy(edgeg);
5296 _lwt_release_edges(edges, num);
5297 }
5298
5299 return 0;
5300 }
5301
5302 /*
5303 * Add a pre-noded pre-split line edge. Used by lwt_AddLine
5304 * Return edge id, 0 if none added (empty edge), -1 on error
5305 *
5306 * @param handleFaceSplit if non-zero the code will check
5307 * if the newly added edge would split a face and if so
5308 * would create new faces accordingly. Otherwise it will
5309 * set left_face and right_face to null (-1)
5310 */
5311 static LWT_ELEMID
_lwt_AddLineEdge(LWT_TOPOLOGY * topo,LWLINE * edge,double tol,int handleFaceSplit)5312 _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol,
5313 int handleFaceSplit )
5314 {
5315 LWCOLLECTION *col;
5316 LWPOINT *start_point, *end_point;
5317 LWGEOM *tmp = 0, *tmp2;
5318 LWT_ISO_NODE *node;
5319 LWT_ELEMID nid[2]; /* start_node, end_node */
5320 LWT_ELEMID id; /* edge id */
5321 POINT4D p4d;
5322 uint64_t nn, i;
5323 int moved=0, mm;
5324
5325 LWDEBUGG(1, lwline_as_lwgeom(edge), "_lwtAddLineEdge");
5326 LWDEBUGF(1, "_lwtAddLineEdge with tolerance %g", tol);
5327
5328 start_point = lwline_get_lwpoint(edge, 0);
5329 if ( ! start_point )
5330 {
5331 lwnotice("Empty component of noded line");
5332 return 0; /* must be empty */
5333 }
5334 nid[0] = _lwt_AddPoint( topo, start_point,
5335 _lwt_minTolerance(lwpoint_as_lwgeom(start_point)),
5336 handleFaceSplit, &mm );
5337 lwpoint_free(start_point); /* too late if lwt_AddPoint calls lwerror */
5338 if ( nid[0] == -1 ) return -1; /* lwerror should have been called */
5339 moved += mm;
5340
5341
5342 end_point = lwline_get_lwpoint(edge, edge->points->npoints-1);
5343 if ( ! end_point )
5344 {
5345 lwerror("could not get last point of line "
5346 "after successfully getting first point !?");
5347 return -1;
5348 }
5349 nid[1] = _lwt_AddPoint( topo, end_point,
5350 _lwt_minTolerance(lwpoint_as_lwgeom(end_point)),
5351 handleFaceSplit, &mm );
5352 moved += mm;
5353 lwpoint_free(end_point); /* too late if lwt_AddPoint calls lwerror */
5354 if ( nid[1] == -1 ) return -1; /* lwerror should have been called */
5355
5356 /*
5357 -- Added endpoints may have drifted due to tolerance, so
5358 -- we need to re-snap the edge to the new nodes before adding it
5359 */
5360 if ( moved )
5361 {
5362
5363 nn = nid[0] == nid[1] ? 1 : 2;
5364 node = lwt_be_getNodeById( topo, nid, &nn,
5365 LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM );
5366 if (nn == UINT64_MAX)
5367 {
5368 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5369 return -1;
5370 }
5371 start_point = NULL; end_point = NULL;
5372 for (i=0; i<nn; ++i)
5373 {
5374 if ( node[i].node_id == nid[0] ) start_point = node[i].geom;
5375 if ( node[i].node_id == nid[1] ) end_point = node[i].geom;
5376 }
5377 if ( ! start_point || ! end_point )
5378 {
5379 if ( nn ) _lwt_release_nodes(node, nn);
5380 lwerror("Could not find just-added nodes % " LWTFMT_ELEMID
5381 " and %" LWTFMT_ELEMID, nid[0], nid[1]);
5382 return -1;
5383 }
5384
5385 /* snap */
5386
5387 getPoint4d_p( start_point->point, 0, &p4d );
5388 lwline_setPoint4d(edge, 0, &p4d);
5389
5390 getPoint4d_p( end_point->point, 0, &p4d );
5391 lwline_setPoint4d(edge, edge->points->npoints-1, &p4d);
5392
5393 if ( nn ) _lwt_release_nodes(node, nn);
5394
5395 /* make valid, after snap (to handle collapses) */
5396 tmp = lwgeom_make_valid(lwline_as_lwgeom(edge));
5397
5398 col = lwgeom_as_lwcollection(tmp);
5399 if ( col )
5400 {{
5401
5402 col = lwcollection_extract(col, LINETYPE);
5403
5404 /* Check if the so-snapped edge collapsed (see #1650) */
5405 if ( col->ngeoms == 0 )
5406 {
5407 lwcollection_free(col);
5408 lwgeom_free(tmp);
5409 LWDEBUG(1, "Made-valid snapped edge collapsed");
5410 return 0;
5411 }
5412
5413 tmp2 = lwgeom_clone_deep( col->geoms[0] );
5414 lwgeom_free(tmp);
5415 tmp = tmp2;
5416 edge = lwgeom_as_lwline(tmp);
5417 lwcollection_free(col);
5418 if ( ! edge )
5419 {
5420 /* should never happen */
5421 lwerror("lwcollection_extract(LINETYPE) returned a non-line?");
5422 return -1;
5423 }
5424 }}
5425 else
5426 {
5427 edge = lwgeom_as_lwline(tmp);
5428 if ( ! edge )
5429 {
5430 LWDEBUGF(1, "Made-valid snapped edge collapsed to %s",
5431 lwtype_name(lwgeom_get_type(tmp)));
5432 lwgeom_free(tmp);
5433 return 0;
5434 }
5435 }
5436 }
5437
5438 /* check if the so-snapped edge _now_ exists */
5439 id = _lwt_GetEqualEdge ( topo, edge );
5440 LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id);
5441 if ( id == -1 )
5442 {
5443 if ( tmp ) lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5444 return -1;
5445 }
5446 if ( id )
5447 {
5448 if ( tmp ) lwgeom_free(tmp); /* possibly takes "edge" down with it */
5449 return id;
5450 }
5451
5452 /* No previously existing edge was found, we'll add one */
5453
5454 /* Remove consecutive vertices below given tolerance
5455 * on edge addition */
5456 if ( tol )
5457 {{
5458 tmp2 = lwline_remove_repeated_points(edge, tol);
5459 LWDEBUGG(1, tmp2, "Repeated-point removed");
5460 edge = lwgeom_as_lwline(tmp2);
5461 if ( tmp ) lwgeom_free(tmp);
5462 tmp = tmp2;
5463
5464 /* check if the so-decimated edge collapsed to single-point */
5465 if ( nid[0] == nid[1] && edge->points->npoints == 2 )
5466 {
5467 lwgeom_free(tmp);
5468 LWDEBUG(1, "Repeated-point removed edge collapsed");
5469 return 0;
5470 }
5471
5472 /* check if the so-decimated edge _now_ exists */
5473 id = _lwt_GetEqualEdge ( topo, edge );
5474 LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id);
5475 if ( id == -1 )
5476 {
5477 lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5478 return -1;
5479 }
5480 if ( id )
5481 {
5482 lwgeom_free(tmp); /* takes "edge" down with it */
5483 return id;
5484 }
5485 }}
5486
5487
5488 /* TODO: skip checks ? */
5489 id = _lwt_AddEdge( topo, nid[0], nid[1], edge, 0, handleFaceSplit ? 1 : -1 );
5490 LWDEBUGF(1, "lwt_AddEdgeModFace returned %" LWTFMT_ELEMID, id);
5491 if ( id == -1 )
5492 {
5493 lwgeom_free(tmp); /* probably too late, due to internal lwerror */
5494 return -1;
5495 }
5496 lwgeom_free(tmp); /* possibly takes "edge" down with it */
5497
5498 return id;
5499 }
5500
5501 /* Simulate split-loop as it was implemented in pl/pgsql version
5502 * of TopoGeo_addLinestring */
5503 static LWGEOM *
_lwt_split_by_nodes(const LWGEOM * g,const LWGEOM * nodes)5504 _lwt_split_by_nodes(const LWGEOM *g, const LWGEOM *nodes)
5505 {
5506 LWCOLLECTION *col = lwgeom_as_lwcollection(nodes);
5507 uint32_t i;
5508 LWGEOM *bg;
5509
5510 bg = lwgeom_clone_deep(g);
5511 if ( ! col->ngeoms ) return bg;
5512
5513 for (i=0; i<col->ngeoms; ++i)
5514 {
5515 LWGEOM *g2;
5516 g2 = lwgeom_split(bg, col->geoms[i]);
5517 lwgeom_free(bg);
5518 bg = g2;
5519 }
5520 bg->srid = nodes->srid;
5521
5522 return bg;
5523 }
5524
5525 static LWT_ELEMID*
_lwt_AddLine(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges,int handleFaceSplit)5526 _lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges,
5527 int handleFaceSplit)
5528 {
5529 LWGEOM *geomsbuf[1];
5530 LWGEOM **geoms;
5531 uint32_t ngeoms;
5532 LWGEOM *noded, *tmp;
5533 LWCOLLECTION *col;
5534 LWT_ELEMID *ids;
5535 LWT_ISO_EDGE *edges;
5536 LWT_ISO_NODE *nodes;
5537 uint64_t num, numedges = 0, numnodes = 0;
5538 uint64_t i;
5539 GBOX qbox;
5540
5541 *nedges = -1; /* error condition, by default */
5542
5543 /* Get tolerance, if 0 was given */
5544 if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)line );
5545 LWDEBUGF(1, "Working tolerance:%.15g", tol);
5546 LWDEBUGF(1, "Input line has srid=%d", line->srid);
5547
5548 /* Remove consecutive vertices below given tolerance upfront */
5549 if ( tol )
5550 {{
5551 LWLINE *clean = lwgeom_as_lwline(lwline_remove_repeated_points(line, tol));
5552 tmp = lwline_as_lwgeom(clean); /* NOTE: might collapse to non-simple */
5553 LWDEBUGG(1, tmp, "Repeated-point removed");
5554 }} else tmp=(LWGEOM*)line;
5555
5556 /* 1. Self-node */
5557 noded = lwgeom_node((LWGEOM*)tmp);
5558 if ( tmp != (LWGEOM*)line ) lwgeom_free(tmp);
5559 if ( ! noded ) return NULL; /* should have called lwerror already */
5560 LWDEBUGG(1, noded, "Noded");
5561
5562 qbox = *lwgeom_get_bbox( lwline_as_lwgeom(line) );
5563 LWDEBUGF(1, "Line BOX is %.15g %.15g, %.15g %.15g", qbox.xmin, qbox.ymin,
5564 qbox.xmax, qbox.ymax);
5565 gbox_expand(&qbox, tol);
5566 LWDEBUGF(1, "BOX expanded by %g is %.15g %.15g, %.15g %.15g",
5567 tol, qbox.xmin, qbox.ymin, qbox.xmax, qbox.ymax);
5568
5569 LWGEOM **nearby = 0;
5570 int nearbyindex = 0;
5571 int nearbycount = 0;
5572
5573 /* 2.0. Find edges falling within tol distance */
5574 edges = lwt_be_getEdgeWithinBox2D( topo, &qbox, &numedges, LWT_COL_EDGE_ALL, 0 );
5575 if (numedges == UINT64_MAX)
5576 {
5577 lwgeom_free(noded);
5578 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5579 return NULL;
5580 }
5581 LWDEBUGF(1, "Line has %d points, its bbox intersects %d edges bboxes",
5582 line->points->npoints, numedges);
5583 if ( numedges )
5584 {{
5585 /* collect those whose distance from us is < tol */
5586 nearbycount += numedges;
5587 nearby = lwalloc(numedges * sizeof(LWGEOM *));
5588 for (i=0; i<numedges; ++i)
5589 {
5590 LW_ON_INTERRUPT(return NULL);
5591 LWT_ISO_EDGE *e = &(edges[i]);
5592 LWGEOM *g = lwline_as_lwgeom(e->geom);
5593 LWDEBUGF(2, "Computing distance from edge %d having %d points", i, e->geom->points->npoints);
5594 double dist = lwgeom_mindistance2d(g, noded);
5595 /* must be closer than tolerated, unless distance is zero */
5596 if ( dist && dist >= tol ) continue;
5597 nearby[nearbyindex++] = g;
5598 }
5599 LWDEBUGF(1, "Found %d edges closer than tolerance (%g)", nearbyindex, tol);
5600 }}
5601 int nearbyedgecount = nearbyindex;
5602
5603 /* 2.1. Find nodes falling within tol distance *
5604 * TODO: check if we should be only considering _isolated_ nodes! */
5605 nodes = lwt_be_getNodeWithinBox2D( topo, &qbox, &numnodes, LWT_COL_NODE_ALL, 0 );
5606 if (numnodes == UINT64_MAX)
5607 {
5608 lwgeom_free(noded);
5609 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5610 return NULL;
5611 }
5612 LWDEBUGF(1, "Line bbox intersects %d nodes bboxes", numnodes);
5613 if ( numnodes )
5614 {{
5615 /* collect those whose distance from us is < tol */
5616 nearbycount = nearbyedgecount + numnodes;
5617 nearby = nearby ?
5618 lwrealloc(nearby, nearbycount * sizeof(LWGEOM *))
5619 :
5620 lwalloc(nearbycount * sizeof(LWGEOM *))
5621 ;
5622 int nn = 0;
5623 for (i=0; i<numnodes; ++i)
5624 {
5625 LWT_ISO_NODE *n = &(nodes[i]);
5626 LWGEOM *g = lwpoint_as_lwgeom(n->geom);
5627 double dist = lwgeom_mindistance2d(g, noded);
5628 /* must be closer than tolerated, unless distance is zero */
5629 if ( dist && dist >= tol )
5630 {
5631 LWDEBUGF(1, "Node %d is %g units away, we tolerate only %g", n->node_id, dist, tol);
5632 continue;
5633 }
5634 nearby[nearbyindex++] = g;
5635 ++nn;
5636 }
5637 LWDEBUGF(1, "Found %d nodes closer than tolerance (%g)", nn, tol);
5638 }}
5639 int nearbynodecount = nearbyindex - nearbyedgecount;
5640 nearbycount = nearbyindex;
5641
5642 LWDEBUGF(1, "Number of nearby elements is %d", nearbycount);
5643
5644 /* 2.2. Snap to nearby elements */
5645 if ( nearbycount )
5646 {{
5647 LWCOLLECTION *col;
5648 LWGEOM *elems;
5649
5650 col = lwcollection_construct(COLLECTIONTYPE, topo->srid,
5651 NULL, nearbycount, nearby);
5652 elems = lwcollection_as_lwgeom(col);
5653
5654 LWDEBUGG(1, elems, "Collected nearby elements");
5655
5656 tmp = _lwt_toposnap(noded, elems, tol);
5657 lwgeom_free(noded);
5658 noded = tmp;
5659 LWDEBUGG(1, noded, "Elements-snapped");
5660
5661 /* will not release the geoms array */
5662 lwcollection_release(col);
5663
5664 /*
5665 -- re-node to account for ST_Snap introduced self-intersections
5666 -- See http://trac.osgeo.org/postgis/ticket/1714
5667 -- TODO: consider running UnaryUnion once after all noding
5668 */
5669 tmp = lwgeom_unaryunion(noded);
5670 lwgeom_free(noded);
5671 noded = tmp;
5672 LWDEBUGG(1, noded, "Unary-unioned");
5673 }}
5674
5675 /* 2.3. Node with nearby edges */
5676 if ( nearbyedgecount )
5677 {{
5678 LWCOLLECTION *col;
5679 LWGEOM *iedges; /* just an alias for col */
5680 LWGEOM *diff, *xset;
5681
5682 LWDEBUGF(1, "Line intersects %d edges", nearbyedgecount);
5683
5684 col = lwcollection_construct(COLLECTIONTYPE, topo->srid,
5685 NULL, nearbyedgecount, nearby);
5686 iedges = lwcollection_as_lwgeom(col);
5687 LWDEBUGG(1, iedges, "Collected edges");
5688
5689 LWDEBUGF(1, "Diffing noded, with srid=%d "
5690 "and interesecting edges, with srid=%d",
5691 noded->srid, iedges->srid);
5692 diff = lwgeom_difference(noded, iedges);
5693 LWDEBUGG(1, diff, "Differenced");
5694
5695 LWDEBUGF(1, "Intersecting noded, with srid=%d "
5696 "and interesecting edges, with srid=%d",
5697 noded->srid, iedges->srid);
5698 xset = lwgeom_intersection(noded, iedges);
5699 LWDEBUGG(1, xset, "Intersected");
5700 lwgeom_free(noded);
5701
5702 /* We linemerge here because INTERSECTION, as of GEOS 3.8,
5703 * will result in shared segments being output as multiple
5704 * lines rather than a single line. Example:
5705
5706 INTERSECTION(
5707 'LINESTRING(0 0, 5 0, 8 0, 10 0,12 0)',
5708 'LINESTRING(5 0, 8 0, 10 0)'
5709 )
5710 ==
5711 MULTILINESTRING((5 0,8 0),(8 0,10 0))
5712
5713 * We will re-split in a subsequent step, by splitting
5714 * the final line with pre-existing nodes
5715 */
5716 LWDEBUG(1, "Linemerging intersection");
5717 tmp = lwgeom_linemerge(xset);
5718 LWDEBUGG(1, tmp, "Linemerged");
5719 lwgeom_free(xset);
5720 xset = tmp;
5721
5722 /*
5723 * Here we union the (linemerged) intersection with
5724 * the difference (new lines)
5725 */
5726 LWDEBUG(1, "Unioning difference and (linemerged) intersection");
5727 noded = lwgeom_union(diff, xset);
5728 LWDEBUGG(1, noded, "Diff-Xset Unioned");
5729 lwgeom_free(xset);
5730 lwgeom_free(diff);
5731
5732 /* will not release the geoms array */
5733 lwcollection_release(col);
5734 }}
5735
5736
5737 /* 2.4. Split by pre-existing nodes
5738 *
5739 * Pre-existing nodes are isolated nodes AND endpoints
5740 * of intersecting edges
5741 */
5742 if ( nearbyedgecount )
5743 {
5744 nearbycount += nearbyedgecount * 2; /* make space for endpoints */
5745 nearby = lwrealloc(nearby, nearbycount * sizeof(LWGEOM *));
5746 for (int i=0; i<nearbyedgecount; i++)
5747 {
5748 LWLINE *edge = lwgeom_as_lwline(nearby[i]);
5749 LWPOINT *startNode = lwline_get_lwpoint(edge, 0);
5750 LWPOINT *endNode = lwline_get_lwpoint(edge, edge->points->npoints-1);
5751 /* TODO: only add if within distance to noded AND if not duplicated */
5752 nearby[nearbyindex++] = lwpoint_as_lwgeom(startNode);
5753 nearbynodecount++;
5754 nearby[nearbyindex++] = lwpoint_as_lwgeom(endNode);
5755 nearbynodecount++;
5756 }
5757 }
5758 if ( nearbynodecount )
5759 {
5760 col = lwcollection_construct(MULTIPOINTTYPE, topo->srid,
5761 NULL, nearbynodecount,
5762 nearby + nearbyedgecount);
5763 LWGEOM *inodes = lwcollection_as_lwgeom(col);
5764 /* TODO: use lwgeom_split of lwgeom_union ... */
5765 tmp = _lwt_split_by_nodes(noded, inodes);
5766 lwgeom_free(noded);
5767 noded = tmp;
5768 LWDEBUGG(1, noded, "Node-split");
5769 /* will not release the geoms array */
5770 lwcollection_release(col);
5771 }
5772
5773
5774 LWDEBUG(1, "Freeing up nearby elements");
5775
5776 /* TODO: free up endpoints of nearbyedges */
5777 if ( nearby ) lwfree(nearby);
5778 if ( nodes ) _lwt_release_nodes(nodes, numnodes);
5779 if ( edges ) _lwt_release_edges(edges, numedges);
5780
5781 LWDEBUGG(1, noded, "Finally-noded");
5782
5783 /* 3. For each (now-noded) segment, insert an edge */
5784 col = lwgeom_as_lwcollection(noded);
5785 if ( col )
5786 {
5787 LWDEBUG(1, "Noded line was a collection");
5788 geoms = col->geoms;
5789 ngeoms = col->ngeoms;
5790 }
5791 else
5792 {
5793 LWDEBUG(1, "Noded line was a single geom");
5794 geomsbuf[0] = noded;
5795 geoms = geomsbuf;
5796 ngeoms = 1;
5797 }
5798
5799 LWDEBUGF(1, "Line was split into %d edges", ngeoms);
5800
5801 /* TODO: refactor to first add all nodes (re-snapping edges if
5802 * needed) and then check all edges for existing already
5803 * ( so to save a DB scan for each edge to be added )
5804 */
5805 ids = lwalloc(sizeof(LWT_ELEMID)*ngeoms);
5806 num = 0;
5807 for ( i=0; i<ngeoms; ++i )
5808 {
5809 LWT_ELEMID id;
5810 LWGEOM *g = geoms[i];
5811 g->srid = noded->srid;
5812
5813 #if POSTGIS_DEBUG_LEVEL > 0
5814 {
5815 size_t sz;
5816 char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5817 LWDEBUGF(1, "Component %d of split line is: %s", i, wkt1);
5818 lwfree(wkt1);
5819 }
5820 #endif
5821
5822 id = _lwt_AddLineEdge( topo, lwgeom_as_lwline(g), tol, handleFaceSplit );
5823 LWDEBUGF(1, "_lwt_AddLineEdge returned %" LWTFMT_ELEMID, id);
5824 if ( id < 0 )
5825 {
5826 lwgeom_free(noded);
5827 lwfree(ids);
5828 return NULL;
5829 }
5830 if ( ! id )
5831 {
5832 LWDEBUGF(1, "Component %d of split line collapsed", i);
5833 continue;
5834 }
5835
5836 LWDEBUGF(1, "Component %d of split line is edge %" LWTFMT_ELEMID,
5837 i, id);
5838 ids[num++] = id; /* TODO: skip duplicates */
5839 }
5840
5841 LWDEBUGG(1, noded, "Noded before free");
5842 lwgeom_free(noded);
5843
5844 /* TODO: XXX remove duplicated ids if not done before */
5845
5846 *nedges = num;
5847 return ids;
5848 }
5849
5850 LWT_ELEMID*
lwt_AddLine(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges)5851 lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges)
5852 {
5853 return _lwt_AddLine(topo, line, tol, nedges, 1);
5854 }
5855
5856 LWT_ELEMID*
lwt_AddLineNoFace(LWT_TOPOLOGY * topo,LWLINE * line,double tol,int * nedges)5857 lwt_AddLineNoFace(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges)
5858 {
5859 return _lwt_AddLine(topo, line, tol, nedges, 0);
5860 }
5861
5862 LWT_ELEMID*
lwt_AddPolygon(LWT_TOPOLOGY * topo,LWPOLY * poly,double tol,int * nfaces)5863 lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* poly, double tol, int* nfaces)
5864 {
5865 uint32_t i;
5866 *nfaces = -1; /* error condition, by default */
5867 int num;
5868 LWT_ISO_FACE *faces;
5869 uint64_t nfacesinbox;
5870 uint64_t j;
5871 LWT_ELEMID *ids = NULL;
5872 GBOX qbox;
5873 const GEOSPreparedGeometry *ppoly;
5874 GEOSGeometry *polyg;
5875
5876 /* Get tolerance, if 0 was given */
5877 if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)poly );
5878 LWDEBUGF(1, "Working tolerance:%.15g", tol);
5879
5880 /* Add each ring as an edge */
5881 for ( i=0; i<poly->nrings; ++i )
5882 {
5883 LWLINE *line;
5884 POINTARRAY *pa;
5885 LWT_ELEMID *eids;
5886 int nedges;
5887
5888 pa = ptarray_clone(poly->rings[i]);
5889 line = lwline_construct(topo->srid, NULL, pa);
5890 eids = lwt_AddLine( topo, line, tol, &nedges );
5891 if ( nedges < 0 ) {
5892 /* probably too late as lwt_AddLine invoked lwerror */
5893 lwline_free(line);
5894 lwerror("Error adding ring %d of polygon", i);
5895 return NULL;
5896 }
5897 lwline_free(line);
5898 lwfree(eids);
5899 }
5900
5901 /*
5902 -- Find faces covered by input polygon
5903 -- NOTE: potential snapping changed polygon edges
5904 */
5905 qbox = *lwgeom_get_bbox( lwpoly_as_lwgeom(poly) );
5906 gbox_expand(&qbox, tol);
5907 faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nfacesinbox,
5908 LWT_COL_FACE_ALL, 0 );
5909 if (nfacesinbox == UINT64_MAX)
5910 {
5911 lwfree(ids);
5912 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5913 return NULL;
5914 }
5915
5916 num = 0;
5917 if ( nfacesinbox )
5918 {
5919 polyg = LWGEOM2GEOS(lwpoly_as_lwgeom(poly), 0);
5920 if ( ! polyg )
5921 {
5922 _lwt_release_faces(faces, nfacesinbox);
5923 lwerror("Could not convert poly geometry to GEOS: %s", lwgeom_geos_errmsg);
5924 return NULL;
5925 }
5926 ppoly = GEOSPrepare(polyg);
5927 ids = lwalloc(sizeof(LWT_ELEMID)*nfacesinbox);
5928 for ( j=0; j<nfacesinbox; ++j )
5929 {
5930 LWT_ISO_FACE *f = &(faces[j]);
5931 LWGEOM *fg;
5932 GEOSGeometry *fgg, *sp;
5933 int covers;
5934
5935 /* check if a point on this face surface is covered by our polygon */
5936 fg = lwt_GetFaceGeometry( topo, f->face_id );
5937 if ( ! fg )
5938 {
5939 j = f->face_id; /* so we can destroy faces */
5940 GEOSPreparedGeom_destroy(ppoly);
5941 GEOSGeom_destroy(polyg);
5942 lwfree(ids);
5943 _lwt_release_faces(faces, nfacesinbox);
5944 lwerror("Could not get geometry of face %" LWTFMT_ELEMID, j);
5945 return NULL;
5946 }
5947 /* check if a point on this face's surface is covered by our polygon */
5948 fgg = LWGEOM2GEOS(fg, 0);
5949 lwgeom_free(fg);
5950 if ( ! fgg )
5951 {
5952 GEOSPreparedGeom_destroy(ppoly);
5953 GEOSGeom_destroy(polyg);
5954 _lwt_release_faces(faces, nfacesinbox);
5955 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5956 return NULL;
5957 }
5958 sp = GEOSPointOnSurface(fgg);
5959 GEOSGeom_destroy(fgg);
5960 if ( ! sp )
5961 {
5962 GEOSPreparedGeom_destroy(ppoly);
5963 GEOSGeom_destroy(polyg);
5964 _lwt_release_faces(faces, nfacesinbox);
5965 lwerror("Could not find point on face surface: %s", lwgeom_geos_errmsg);
5966 return NULL;
5967 }
5968 covers = GEOSPreparedCovers( ppoly, sp );
5969 GEOSGeom_destroy(sp);
5970 if (covers == 2)
5971 {
5972 GEOSPreparedGeom_destroy(ppoly);
5973 GEOSGeom_destroy(polyg);
5974 _lwt_release_faces(faces, nfacesinbox);
5975 lwerror("PreparedCovers error: %s", lwgeom_geos_errmsg);
5976 return NULL;
5977 }
5978 if ( ! covers )
5979 {
5980 continue; /* we're not composed by this face */
5981 }
5982
5983 /* TODO: avoid duplicates ? */
5984 ids[num++] = f->face_id;
5985 }
5986 GEOSPreparedGeom_destroy(ppoly);
5987 GEOSGeom_destroy(polyg);
5988 _lwt_release_faces(faces, nfacesinbox);
5989 }
5990
5991 /* possibly 0 if non face's surface point was found
5992 * to be covered by input polygon */
5993 *nfaces = num;
5994
5995 return ids;
5996 }
5997
5998 /*
5999 *---- polygonizer
6000 */
6001
6002 /* An array of pointers to EDGERING structures */
6003 typedef struct LWT_ISO_EDGE_TABLE_T {
6004 LWT_ISO_EDGE *edges;
6005 int size;
6006 } LWT_ISO_EDGE_TABLE;
6007
6008 static int
compare_iso_edges_by_id(const void * si1,const void * si2)6009 compare_iso_edges_by_id(const void *si1, const void *si2)
6010 {
6011 int a = ((LWT_ISO_EDGE *)si1)->edge_id;
6012 int b = ((LWT_ISO_EDGE *)si2)->edge_id;
6013 if ( a < b )
6014 return -1;
6015 else if ( a > b )
6016 return 1;
6017 else
6018 return 0;
6019 }
6020
6021 static LWT_ISO_EDGE *
_lwt_getIsoEdgeById(LWT_ISO_EDGE_TABLE * tab,LWT_ELEMID id)6022 _lwt_getIsoEdgeById(LWT_ISO_EDGE_TABLE *tab, LWT_ELEMID id)
6023 {
6024 LWT_ISO_EDGE key;
6025 key.edge_id = id;
6026
6027 void *match = bsearch( &key, tab->edges, tab->size,
6028 sizeof(LWT_ISO_EDGE),
6029 compare_iso_edges_by_id);
6030 return match;
6031 }
6032
6033 typedef struct LWT_EDGERING_ELEM_T {
6034 /* externally owned */
6035 LWT_ISO_EDGE *edge;
6036 /* 0 if false, 1 if true */
6037 int left;
6038 } LWT_EDGERING_ELEM;
6039
6040 /* A ring of edges */
6041 typedef struct LWT_EDGERING_T {
6042 /* Signed edge identifiers
6043 * positive ones are walked in their direction, negative ones
6044 * in the opposite direction */
6045 LWT_EDGERING_ELEM **elems;
6046 /* Number of edges in the ring */
6047 int size;
6048 int capacity;
6049 /* Bounding box of the ring */
6050 GBOX *env;
6051 /* Bounding box of the ring in GEOS format (for STRTree) */
6052 GEOSGeometry *genv;
6053 } LWT_EDGERING;
6054
6055 #define LWT_EDGERING_INIT(a) { \
6056 (a)->size = 0; \
6057 (a)->capacity = 1; \
6058 (a)->elems = lwalloc(sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \
6059 (a)->env = NULL; \
6060 (a)->genv = NULL; \
6061 }
6062
6063 #define LWT_EDGERING_PUSH(a, r) { \
6064 if ( (a)->size + 1 > (a)->capacity ) { \
6065 (a)->capacity *= 2; \
6066 (a)->elems = lwrealloc((a)->elems, sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \
6067 } \
6068 /* lwdebug(1, "adding elem %d (%p) of edgering %p", (a)->size, (r), (a)); */ \
6069 (a)->elems[(a)->size++] = (r); \
6070 }
6071
6072 #define LWT_EDGERING_CLEAN(a) { \
6073 int i; for (i=0; i<(a)->size; ++i) { \
6074 if ( (a)->elems[i] ) { \
6075 /* lwdebug(1, "freeing elem %d (%p) of edgering %p", i, (a)->elems[i], (a)); */ \
6076 lwfree((a)->elems[i]); \
6077 } \
6078 } \
6079 if ( (a)->elems ) { lwfree((a)->elems); (a)->elems = NULL; } \
6080 (a)->size = 0; \
6081 (a)->capacity = 0; \
6082 if ( (a)->env ) { lwfree((a)->env); (a)->env = NULL; } \
6083 if ( (a)->genv ) { GEOSGeom_destroy((a)->genv); (a)->genv = NULL; } \
6084 }
6085
6086 /* An array of pointers to EDGERING structures */
6087 typedef struct LWT_EDGERING_ARRAY_T {
6088 LWT_EDGERING **rings;
6089 int size;
6090 int capacity;
6091 GEOSSTRtree* tree;
6092 } LWT_EDGERING_ARRAY;
6093
6094 #define LWT_EDGERING_ARRAY_INIT(a) { \
6095 (a)->size = 0; \
6096 (a)->capacity = 1; \
6097 (a)->rings = lwalloc(sizeof(LWT_EDGERING *) * (a)->capacity); \
6098 (a)->tree = NULL; \
6099 }
6100
6101 /* WARNING: use of 'j' is intentional, not to clash with
6102 * 'i' used in LWT_EDGERING_CLEAN */
6103 #define LWT_EDGERING_ARRAY_CLEAN(a) { \
6104 int j; for (j=0; j<(a)->size; ++j) { \
6105 LWT_EDGERING_CLEAN((a)->rings[j]); \
6106 } \
6107 if ( (a)->capacity ) lwfree((a)->rings); \
6108 if ( (a)->tree ) { \
6109 GEOSSTRtree_destroy( (a)->tree ); \
6110 (a)->tree = NULL; \
6111 } \
6112 }
6113
6114 #define LWT_EDGERING_ARRAY_PUSH(a, r) { \
6115 if ( (a)->size + 1 > (a)->capacity ) { \
6116 (a)->capacity *= 2; \
6117 (a)->rings = lwrealloc((a)->rings, sizeof(LWT_EDGERING *) * (a)->capacity); \
6118 } \
6119 (a)->rings[(a)->size++] = (r); \
6120 }
6121
6122 typedef struct LWT_EDGERING_POINT_ITERATOR_T {
6123 LWT_EDGERING *ring;
6124 LWT_EDGERING_ELEM *curelem;
6125 int curelemidx;
6126 int curidx;
6127 } LWT_EDGERING_POINT_ITERATOR;
6128
6129 static int
_lwt_EdgeRingIterator_next(LWT_EDGERING_POINT_ITERATOR * it,POINT2D * pt)6130 _lwt_EdgeRingIterator_next(LWT_EDGERING_POINT_ITERATOR *it, POINT2D *pt)
6131 {
6132 LWT_EDGERING_ELEM *el = it->curelem;
6133 POINTARRAY *pa;
6134
6135 if ( ! el ) return 0; /* finished */
6136
6137 pa = el->edge->geom->points;
6138
6139 int tonext = 0;
6140 LWDEBUGF(3, "iterator fetching idx %d from pa of %d points", it->curidx, pa->npoints);
6141 getPoint2d_p(pa, it->curidx, pt);
6142 if ( el->left ) {
6143 it->curidx++;
6144 if ( it->curidx >= (int) pa->npoints ) tonext = 1;
6145 } else {
6146 it->curidx--;
6147 if ( it->curidx < 0 ) tonext = 1;
6148 }
6149
6150 if ( tonext )
6151 {
6152 LWDEBUG(3, "iterator moving to next element");
6153 it->curelemidx++;
6154 if ( it->curelemidx < it->ring->size )
6155 {
6156 el = it->curelem = it->ring->elems[it->curelemidx];
6157 it->curidx = el->left ? 0 : el->edge->geom->points->npoints - 1;
6158 }
6159 else
6160 {
6161 it->curelem = NULL;
6162 }
6163 }
6164
6165 return 1;
6166 }
6167
6168 /* Release return with lwfree */
6169 static LWT_EDGERING_POINT_ITERATOR *
_lwt_EdgeRingIterator_begin(LWT_EDGERING * er)6170 _lwt_EdgeRingIterator_begin(LWT_EDGERING *er)
6171 {
6172 LWT_EDGERING_POINT_ITERATOR *ret = lwalloc(sizeof(LWT_EDGERING_POINT_ITERATOR));
6173 ret->ring = er;
6174 if ( er->size ) ret->curelem = er->elems[0];
6175 else ret->curelem = NULL;
6176 ret->curelemidx = 0;
6177 ret->curidx = ret->curelem->left ? 0 : ret->curelem->edge->geom->points->npoints - 1;
6178 return ret;
6179 }
6180
6181 /* Identifier for a placeholder face that will be
6182 * used to mark hole rings */
6183 #define LWT_HOLES_FACE_PLACEHOLDER INT32_MIN
6184
6185 static int
_lwt_FetchNextUnvisitedEdge(LWT_TOPOLOGY * topo,LWT_ISO_EDGE_TABLE * etab,int from)6186 _lwt_FetchNextUnvisitedEdge(__attribute__((__unused__)) LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *etab, int from)
6187 {
6188 while (
6189 from < etab->size &&
6190 etab->edges[from].face_left != -1 &&
6191 etab->edges[from].face_right != -1
6192 ) from++;
6193 return from < etab->size ? from : -1;
6194 }
6195
6196 static LWT_ISO_EDGE *
_lwt_FetchAllEdges(LWT_TOPOLOGY * topo,int * numedges)6197 _lwt_FetchAllEdges(LWT_TOPOLOGY *topo, int *numedges)
6198 {
6199 LWT_ISO_EDGE *edge;
6200 int fields = LWT_COL_EDGE_ALL;
6201 uint64_t nelems = 1;
6202
6203 edge = lwt_be_getEdgeWithinBox2D( topo, NULL, &nelems, fields, 0);
6204 *numedges = nelems;
6205 if (nelems == UINT64_MAX)
6206 {
6207 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6208 return NULL;
6209 }
6210 return edge;
6211 }
6212
6213 /* Update the side face of given ring edges
6214 *
6215 * Edge identifiers are signed, those with negative identifier
6216 * need to be updated their right_face, those with positive
6217 * identifier need to be updated their left_face.
6218 *
6219 * @param face identifier of the face bound by the ring
6220 * @return 0 on success, -1 on error
6221 */
6222 static int
_lwt_UpdateEdgeRingSideFace(LWT_TOPOLOGY * topo,LWT_EDGERING * ring,LWT_ELEMID face)6223 _lwt_UpdateEdgeRingSideFace(LWT_TOPOLOGY *topo, LWT_EDGERING *ring,
6224 LWT_ELEMID face)
6225 {
6226 LWT_ISO_EDGE *forward_edges = NULL;
6227 int forward_edges_count = 0;
6228 LWT_ISO_EDGE *backward_edges = NULL;
6229 int backward_edges_count = 0;
6230 int i, ret;
6231
6232 /* Make a list of forward_edges and backward_edges */
6233
6234 forward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size);
6235 forward_edges_count = 0;
6236 backward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size);
6237 backward_edges_count = 0;
6238
6239 for ( i=0; i<ring->size; ++i )
6240 {
6241 LWT_EDGERING_ELEM *elem = ring->elems[i];
6242 LWT_ISO_EDGE *edge = elem->edge;
6243 LWT_ELEMID id = edge->edge_id;
6244 if ( elem->left )
6245 {
6246 LWDEBUGF(3, "Forward edge %d is %d", forward_edges_count, id);
6247 forward_edges[forward_edges_count].edge_id = id;
6248 forward_edges[forward_edges_count++].face_left = face;
6249 edge->face_left = face;
6250 }
6251 else
6252 {
6253 LWDEBUGF(3, "Backward edge %d is %d", forward_edges_count, id);
6254 backward_edges[backward_edges_count].edge_id = id;
6255 backward_edges[backward_edges_count++].face_right = face;
6256 edge->face_right = face;
6257 }
6258 }
6259
6260 /* Update forward edges */
6261 if ( forward_edges_count )
6262 {
6263 ret = lwt_be_updateEdgesById(topo, forward_edges,
6264 forward_edges_count,
6265 LWT_COL_EDGE_FACE_LEFT);
6266 if ( ret == -1 )
6267 {
6268 lwfree( forward_edges );
6269 lwfree( backward_edges );
6270 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6271 return -1;
6272 }
6273 if ( ret != forward_edges_count )
6274 {
6275 lwfree( forward_edges );
6276 lwfree( backward_edges );
6277 lwerror("Unexpected error: %d edges updated when expecting %d (forward)",
6278 ret, forward_edges_count);
6279 return -1;
6280 }
6281 }
6282
6283 /* Update backward edges */
6284 if ( backward_edges_count )
6285 {
6286 ret = lwt_be_updateEdgesById(topo, backward_edges,
6287 backward_edges_count,
6288 LWT_COL_EDGE_FACE_RIGHT);
6289 if ( ret == -1 )
6290 {
6291 lwfree( forward_edges );
6292 lwfree( backward_edges );
6293 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6294 return -1;
6295 }
6296 if ( ret != backward_edges_count )
6297 {
6298 lwfree( forward_edges );
6299 lwfree( backward_edges );
6300 lwerror("Unexpected error: %d edges updated when expecting %d (backward)",
6301 ret, backward_edges_count);
6302 return -1;
6303 }
6304 }
6305
6306 lwfree( forward_edges );
6307 lwfree( backward_edges );
6308
6309 return 0;
6310 }
6311
6312 /*
6313 * @param side 1 for left side, -1 for right side
6314 */
6315 static LWT_EDGERING *
_lwt_BuildEdgeRing(LWT_TOPOLOGY * topo,LWT_ISO_EDGE_TABLE * edges,LWT_ISO_EDGE * edge,int side)6316 _lwt_BuildEdgeRing(__attribute__((__unused__)) LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *edges,
6317 LWT_ISO_EDGE *edge, int side)
6318 {
6319 LWT_EDGERING *ring;
6320 LWT_EDGERING_ELEM *elem;
6321 LWT_ISO_EDGE *cur;
6322 int curside;
6323
6324 ring = lwalloc(sizeof(LWT_EDGERING));
6325 LWT_EDGERING_INIT(ring);
6326
6327 cur = edge;
6328 curside = side;
6329
6330 LWDEBUGF(2, "Building rings for edge %d (side %d)", cur->edge_id, side);
6331
6332 do {
6333 LWT_ELEMID next;
6334
6335 elem = lwalloc(sizeof(LWT_EDGERING_ELEM));
6336 elem->edge = cur;
6337 elem->left = ( curside == 1 );
6338
6339 /* Mark edge as "visited" */
6340 if ( elem->left ) cur->face_left = LWT_HOLES_FACE_PLACEHOLDER;
6341 else cur->face_right = LWT_HOLES_FACE_PLACEHOLDER;
6342
6343 LWT_EDGERING_PUSH(ring, elem);
6344 next = elem->left ? cur->next_left : cur->next_right;
6345
6346 LWDEBUGF(3, " next edge is %d", next);
6347
6348 if ( next > 0 ) curside = 1;
6349 else { curside = -1; next = -next; }
6350 cur = _lwt_getIsoEdgeById(edges, next);
6351 if ( ! cur )
6352 {
6353 lwerror("Could not find edge with id %d", next);
6354 break;
6355 }
6356 } while (cur != edge || curside != side);
6357
6358 LWDEBUGF(1, "Ring for edge %d has %d elems", edge->edge_id*side, ring->size);
6359
6360 return ring;
6361 }
6362
6363 static double
_lwt_EdgeRingSignedArea(LWT_EDGERING_POINT_ITERATOR * it)6364 _lwt_EdgeRingSignedArea(LWT_EDGERING_POINT_ITERATOR *it)
6365 {
6366 POINT2D P1;
6367 POINT2D P2;
6368 POINT2D P3;
6369 double sum = 0.0;
6370 double x0, x, y1, y2;
6371
6372 if ( ! _lwt_EdgeRingIterator_next(it, &P1) ) return 0.0;
6373 if ( ! _lwt_EdgeRingIterator_next(it, &P2) ) return 0.0;
6374
6375 LWDEBUG(2, "_lwt_EdgeRingSignedArea");
6376
6377 x0 = P1.x;
6378 while ( _lwt_EdgeRingIterator_next(it, &P3) )
6379 {
6380 x = P2.x - x0;
6381 y1 = P3.y;
6382 y2 = P1.y;
6383 sum += x * (y2-y1);
6384
6385 /* Move forwards! */
6386 P1 = P2;
6387 P2 = P3;
6388 }
6389
6390 return sum / 2.0;
6391 }
6392
6393
6394 /* Return 1 for true, 0 for false */
6395 static int
_lwt_EdgeRingIsCCW(LWT_EDGERING * ring)6396 _lwt_EdgeRingIsCCW(LWT_EDGERING *ring)
6397 {
6398 double sa;
6399
6400 LWDEBUGF(2, "_lwt_EdgeRingIsCCW, ring has %d elems", ring->size);
6401 LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring);
6402 sa = _lwt_EdgeRingSignedArea(it);
6403 LWDEBUGF(2, "_lwt_EdgeRingIsCCW, signed area is %g", sa);
6404 lwfree(it);
6405 if ( sa >= 0 ) return 0;
6406 else return 1;
6407 }
6408
6409 static int
_lwt_EdgeRingCrossingCount(const POINT2D * p,LWT_EDGERING_POINT_ITERATOR * it)6410 _lwt_EdgeRingCrossingCount(const POINT2D *p, LWT_EDGERING_POINT_ITERATOR *it)
6411 {
6412 int cn = 0; /* the crossing number counter */
6413 POINT2D v1, v2;
6414 #ifndef RELAX
6415 POINT2D v0;
6416 #endif
6417
6418 if ( ! _lwt_EdgeRingIterator_next(it, &v1) ) return cn;
6419 v0 = v1;
6420 while ( _lwt_EdgeRingIterator_next(it, &v2) )
6421 {
6422 double vt;
6423
6424 /* edge from vertex i to vertex i+1 */
6425 if
6426 (
6427 /* an upward crossing */
6428 ((v1.y <= p->y) && (v2.y > p->y))
6429 /* a downward crossing */
6430 || ((v1.y > p->y) && (v2.y <= p->y))
6431 )
6432 {
6433
6434 vt = (double)(p->y - v1.y) / (v2.y - v1.y);
6435
6436 /* P->x <intersect */
6437 if (p->x < v1.x + vt * (v2.x - v1.x))
6438 {
6439 /* a valid crossing of y=p->y right of p->x */
6440 ++cn;
6441 }
6442 }
6443 v1 = v2;
6444 }
6445
6446 LWDEBUGF(3, "_lwt_EdgeRingCrossingCount returning %d", cn);
6447
6448 #ifndef RELAX
6449 if ( memcmp(&v1, &v0, sizeof(POINT2D)) )
6450 {
6451 lwerror("_lwt_EdgeRingCrossingCount: V[n] != V[0] (%g %g != %g %g)",
6452 v1.x, v1.y, v0.x, v0.y);
6453 return -1;
6454 }
6455 #endif
6456
6457 return cn;
6458 }
6459
6460 /* Return 1 for true, 0 for false */
6461 static int
_lwt_EdgeRingContainsPoint(LWT_EDGERING * ring,POINT2D * p)6462 _lwt_EdgeRingContainsPoint(LWT_EDGERING *ring, POINT2D *p)
6463 {
6464 int cn = 0;
6465
6466 LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring);
6467 cn = _lwt_EdgeRingCrossingCount(p, it);
6468 lwfree(it);
6469 return (cn&1); /* 0 if even (out), and 1 if odd (in) */
6470 }
6471
6472 static GBOX *
_lwt_EdgeRingGetBbox(LWT_EDGERING * ring)6473 _lwt_EdgeRingGetBbox(LWT_EDGERING *ring)
6474 {
6475 int i;
6476
6477 if ( ! ring->env )
6478 {
6479 LWDEBUGF(2, "Computing GBOX for ring %p", ring);
6480 for (i=0; i<ring->size; ++i)
6481 {
6482 LWT_EDGERING_ELEM *elem = ring->elems[i];
6483 LWLINE *g = elem->edge->geom;
6484 const GBOX *newbox = lwgeom_get_bbox(lwline_as_lwgeom(g));
6485 if ( ! i ) ring->env = gbox_clone( newbox );
6486 else gbox_merge( newbox, ring->env );
6487 }
6488 }
6489
6490 return ring->env;
6491 }
6492
6493 static LWT_ELEMID
_lwt_EdgeRingGetFace(LWT_EDGERING * ring)6494 _lwt_EdgeRingGetFace(LWT_EDGERING *ring)
6495 {
6496 LWT_EDGERING_ELEM *el = ring->elems[0];
6497 return el->left ? el->edge->face_left : el->edge->face_right;
6498 }
6499
6500
6501 /*
6502 * Register a face on an edge side
6503 *
6504 * Create and register face to shell (CCW) walks,
6505 * register arbitrary negative face_id to CW rings.
6506 *
6507 * Push CCW rings to shells, CW rings to holes.
6508 *
6509 * The ownership of the "geom" and "ids" members of the
6510 * LWT_EDGERING pushed to the given LWT_EDGERING_ARRAYS
6511 * are transferred to caller.
6512 *
6513 * @param side 1 for left side, -1 for right side
6514 *
6515 * @param holes an array where holes will be pushed
6516 *
6517 * @param shells an array where shells will be pushed
6518 *
6519 * @param registered id of registered face. It will be a negative number
6520 * for holes or isolated edge strips (still registered in the face
6521 * table, but only temporary).
6522 *
6523 * @return 0 on success, -1 on error.
6524 *
6525 */
6526 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)6527 _lwt_RegisterFaceOnEdgeSide(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edge,
6528 int side, LWT_ISO_EDGE_TABLE *edges,
6529 LWT_EDGERING_ARRAY *holes,
6530 LWT_EDGERING_ARRAY *shells,
6531 LWT_ELEMID *registered)
6532 {
6533 const LWT_BE_IFACE *iface = topo->be_iface;
6534 /* this is arbitrary, could be taken as parameter */
6535 static const int placeholder_faceid = LWT_HOLES_FACE_PLACEHOLDER;
6536 LWT_EDGERING *ring;
6537
6538 /* Get edge ring */
6539 ring = _lwt_BuildEdgeRing(topo, edges, edge, side);
6540
6541 LWDEBUG(2, "Ring built, calling EdgeRingIsCCW");
6542
6543 /* Compute winding (CW or CCW?) */
6544 int isccw = _lwt_EdgeRingIsCCW(ring);
6545
6546 if ( isccw )
6547 {
6548 /* Create new face */
6549 LWT_ISO_FACE newface;
6550
6551 LWDEBUGF(1, "Ring of edge %d is a shell (shell %d)", edge->edge_id * side, shells->size);
6552
6553 newface.mbr = _lwt_EdgeRingGetBbox(ring);
6554
6555 newface.face_id = -1;
6556 /* Insert the new face */
6557 int ret = lwt_be_insertFaces( topo, &newface, 1 );
6558 newface.mbr = NULL;
6559 if ( ret == -1 )
6560 {
6561 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6562 return -1;
6563 }
6564 if ( ret != 1 )
6565 {
6566 lwerror("Unexpected error: %d faces inserted when expecting 1", ret);
6567 return -1;
6568 }
6569 /* return new face_id */
6570 *registered = newface.face_id;
6571 LWT_EDGERING_ARRAY_PUSH(shells, ring);
6572
6573 /* update ring edges set new face_id on resp. side to *registered */
6574 ret = _lwt_UpdateEdgeRingSideFace(topo, ring, *registered);
6575 if ( ret )
6576 {
6577 lwerror("Errors updating edgering side face: %s",
6578 lwt_be_lastErrorMessage(iface));
6579 return -1;
6580 }
6581
6582 }
6583 else /* cw, so is an hole */
6584 {
6585 LWDEBUGF(1, "Ring of edge %d is a hole (hole %d)", edge->edge_id * side, holes->size);
6586 *registered = placeholder_faceid;
6587 LWT_EDGERING_ARRAY_PUSH(holes, ring);
6588 }
6589
6590 return 0;
6591 }
6592
6593 static void
_lwt_AccumulateCanditates(void * item,void * userdata)6594 _lwt_AccumulateCanditates(void* item, void* userdata)
6595 {
6596 LWT_EDGERING_ARRAY *candidates = userdata;
6597 LWT_EDGERING *sring = item;
6598 LWT_EDGERING_ARRAY_PUSH(candidates, sring);
6599 }
6600
6601 static LWT_ELEMID
_lwt_FindFaceContainingRing(LWT_TOPOLOGY * topo,LWT_EDGERING * ring,LWT_EDGERING_ARRAY * shells)6602 _lwt_FindFaceContainingRing(LWT_TOPOLOGY* topo, LWT_EDGERING *ring,
6603 LWT_EDGERING_ARRAY *shells)
6604 {
6605 LWT_ELEMID foundInFace = -1;
6606 int i;
6607 const GBOX *minenv = NULL;
6608 POINT2D pt;
6609 const GBOX *testbox;
6610 GEOSGeometry *ghole;
6611
6612 getPoint2d_p( ring->elems[0]->edge->geom->points, 0, &pt );
6613
6614 testbox = _lwt_EdgeRingGetBbox(ring);
6615
6616 /* Create a GEOS Point from a vertex of the hole ring */
6617 {
6618 LWPOINT *point = lwpoint_make2d(topo->srid, pt.x, pt.y);
6619 ghole = LWGEOM2GEOS( lwpoint_as_lwgeom(point), 1 );
6620 lwpoint_free(point);
6621 if ( ! ghole ) {
6622 lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
6623 return -1;
6624 }
6625 }
6626
6627 /* Build STRtree of shell envelopes */
6628 if ( ! shells->tree )
6629 {
6630 static const int STRTREE_NODE_CAPACITY = 10;
6631 LWDEBUG(1, "Building STRtree");
6632 shells->tree = GEOSSTRtree_create(STRTREE_NODE_CAPACITY);
6633 if (shells->tree == NULL)
6634 {
6635 lwerror("Could not create GEOS STRTree: %s", lwgeom_geos_errmsg);
6636 return -1;
6637 }
6638 for (i=0; i<shells->size; ++i)
6639 {
6640 LWT_EDGERING *sring = shells->rings[i];
6641 const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring);
6642 LWDEBUGF(2, "GBOX of shell %p for edge %d is %g %g,%g %g",
6643 sring, sring->elems[0]->edge->edge_id, shellbox->xmin,
6644 shellbox->ymin, shellbox->xmax, shellbox->ymax);
6645 POINTARRAY *pa = ptarray_construct(0, 0, 2);
6646 POINT4D pt;
6647 LWLINE *diag;
6648 pt.x = shellbox->xmin;
6649 pt.y = shellbox->ymin;
6650 ptarray_set_point4d(pa, 0, &pt);
6651 pt.x = shellbox->xmax;
6652 pt.y = shellbox->ymax;
6653 ptarray_set_point4d(pa, 1, &pt);
6654 diag = lwline_construct(topo->srid, NULL, pa);
6655 /* Record just envelope in ggeom */
6656 /* making valid, probably not needed */
6657 sring->genv = LWGEOM2GEOS( lwline_as_lwgeom(diag), 1 );
6658 lwline_free(diag);
6659 GEOSSTRtree_insert(shells->tree, sring->genv, sring);
6660 }
6661 LWDEBUG(1, "STRtree build completed");
6662 }
6663
6664 LWT_EDGERING_ARRAY candidates;
6665 LWT_EDGERING_ARRAY_INIT(&candidates);
6666 GEOSSTRtree_query(shells->tree, ghole, &_lwt_AccumulateCanditates, &candidates);
6667 LWDEBUGF(1, "Found %d candidate shells containing first point of ring's originating edge %d",
6668 candidates.size, ring->elems[0]->edge->edge_id * ( ring->elems[0]->left ? 1 : -1 ) );
6669
6670 /* TODO: sort candidates by bounding box size */
6671
6672 for (i=0; i<candidates.size; ++i)
6673 {
6674 LWT_EDGERING *sring = candidates.rings[i];
6675 const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring);
6676 int contains = 0;
6677
6678 if ( sring->elems[0]->edge->edge_id == ring->elems[0]->edge->edge_id )
6679 {
6680 LWDEBUGF(1, "Shell %d is on other side of ring",
6681 _lwt_EdgeRingGetFace(sring));
6682 continue;
6683 }
6684
6685 /* The hole envelope cannot equal the shell envelope */
6686 if ( gbox_same(shellbox, testbox) )
6687 {
6688 LWDEBUGF(1, "Bbox of shell %d equals that of hole ring",
6689 _lwt_EdgeRingGetFace(sring));
6690 continue;
6691 }
6692
6693 /* Skip if ring box is not in shell box */
6694 if ( ! gbox_contains_2d(shellbox, testbox) )
6695 {
6696 LWDEBUGF(1, "Bbox of shell %d does not contain bbox of ring point",
6697 _lwt_EdgeRingGetFace(sring));
6698 continue;
6699 }
6700
6701 /* Skip test if a containing shell was already found
6702 * and this shell's bbox is not contained in the other */
6703 if ( minenv && ! gbox_contains_2d(minenv, shellbox) )
6704 {
6705 LWDEBUGF(2, "Bbox of shell %d (face %d) not contained by bbox "
6706 "of last shell found to contain the point",
6707 i, _lwt_EdgeRingGetFace(sring));
6708 continue;
6709 }
6710
6711 contains = _lwt_EdgeRingContainsPoint(sring, &pt);
6712 if ( contains )
6713 {
6714 /* Continue until all shells are tested, as we want to
6715 * use the one with the smallest bounding box */
6716 /* IDEA: sort shells by bbox size, stopping on first match */
6717 LWDEBUGF(1, "Shell %d contains hole of edge %d",
6718 _lwt_EdgeRingGetFace(sring),
6719 ring->elems[0]->edge->edge_id);
6720 minenv = shellbox;
6721 foundInFace = _lwt_EdgeRingGetFace(sring);
6722 }
6723 }
6724 if ( foundInFace == -1 ) foundInFace = 0;
6725
6726 candidates.size = 0; /* Avoid destroying the actual shell rings */
6727 LWT_EDGERING_ARRAY_CLEAN(&candidates);
6728
6729 GEOSGeom_destroy(ghole);
6730
6731 return foundInFace;
6732 }
6733
6734 /*
6735 * @return -1 on error (and report error),
6736 * 1 if faces beside the universal one exist
6737 * 0 otherwise
6738 */
6739 static int
_lwt_CheckFacesExist(LWT_TOPOLOGY * topo)6740 _lwt_CheckFacesExist(LWT_TOPOLOGY *topo)
6741 {
6742 LWT_ISO_FACE *faces;
6743 int fields = LWT_COL_FACE_FACE_ID;
6744 uint64_t nelems = 1;
6745 GBOX qbox;
6746
6747 qbox.xmin = qbox.ymin = -DBL_MAX;
6748 qbox.xmax = qbox.ymax = DBL_MAX;
6749 faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nelems, fields, 1);
6750 if (nelems == UINT64_MAX)
6751 {
6752 lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
6753 return -1;
6754 }
6755 if ( faces ) _lwt_release_faces(faces, nelems);
6756 return nelems;
6757 }
6758
6759 int
lwt_Polygonize(LWT_TOPOLOGY * topo)6760 lwt_Polygonize(LWT_TOPOLOGY* topo)
6761 {
6762 /*
6763 Fetch all edges
6764 Sort edges by edge_id
6765 Mark all edges' left and right face as -1
6766 Iteratively:
6767 Fetch next edge with left or right face == -1
6768 For each side with face == -1:
6769 Find ring on its side
6770 If ring is CCW:
6771 create a new face, assign to the ring edges' appropriate side
6772 If ring is CW (face needs to be same of external):
6773 assign a negative face_id the ring edges' appropriate side
6774 Now for each edge with a negative face_id on the side:
6775 Find containing face (mbr cache and all)
6776 Update with id of containing face
6777 */
6778
6779 const LWT_BE_IFACE *iface = topo->be_iface;
6780 LWT_ISO_EDGE *edge;
6781 int numfaces = -1;
6782 LWT_ISO_EDGE_TABLE edgetable;
6783 LWT_EDGERING_ARRAY holes, shells;
6784 int i;
6785 int err = 0;
6786
6787 LWT_EDGERING_ARRAY_INIT(&holes);
6788 LWT_EDGERING_ARRAY_INIT(&shells);
6789
6790 initGEOS(lwnotice, lwgeom_geos_error);
6791
6792 /*
6793 Check if Topology already contains some Face
6794 (ignoring the Universal Face)
6795 */
6796 numfaces = _lwt_CheckFacesExist(topo);
6797 if ( numfaces != 0 ) {
6798 if ( numfaces > 0 ) {
6799 /* Faces exist */
6800 lwerror("Polygonize: face table is not empty.");
6801 }
6802 /* Backend error, message should have been printed already */
6803 return -1;
6804 }
6805
6806
6807 edgetable.edges = _lwt_FetchAllEdges(topo, &(edgetable.size));
6808 if ( ! edgetable.edges ) {
6809 if (edgetable.size == 0) {
6810 /* not an error: no Edges */
6811 return 0;
6812 }
6813 /* error should have been printed already */
6814 return -1;
6815 }
6816
6817 /* Sort edges by ID (to allow btree searches) */
6818 qsort(edgetable.edges, edgetable.size, sizeof(LWT_ISO_EDGE), compare_iso_edges_by_id);
6819
6820 /* Mark all edges as unvisited */
6821 for (i=0; i<edgetable.size; ++i)
6822 edgetable.edges[i].face_left = edgetable.edges[i].face_right = -1;
6823
6824 i = 0;
6825 while (1)
6826 {
6827 i = _lwt_FetchNextUnvisitedEdge(topo, &edgetable, i);
6828 if ( i < 0 ) break; /* end of unvisited */
6829 edge = &(edgetable.edges[i]);
6830
6831 LWT_ELEMID newface = -1;
6832
6833 LWDEBUGF(1, "Next face-missing edge has id:%d, face_left:%d, face_right:%d",
6834 edge->edge_id, edge->face_left, edge->face_right);
6835 if ( edge->face_left == -1 )
6836 {
6837 err = _lwt_RegisterFaceOnEdgeSide(topo, edge, 1, &edgetable,
6838 &holes, &shells, &newface);
6839 if ( err ) break;
6840 LWDEBUGF(1, "New face on the left of edge %d is %d",
6841 edge->edge_id, newface);
6842 edge->face_left = newface;
6843 }
6844 if ( edge->face_right == -1 )
6845 {
6846 err = _lwt_RegisterFaceOnEdgeSide(topo, edge, -1, &edgetable,
6847 &holes, &shells, &newface);
6848 if ( err ) break;
6849 LWDEBUGF(1, "New face on the right of edge %d is %d",
6850 edge->edge_id, newface);
6851 edge->face_right = newface;
6852 }
6853 }
6854
6855 if ( err )
6856 {
6857 _lwt_release_edges(edgetable.edges, edgetable.size);
6858 LWT_EDGERING_ARRAY_CLEAN( &holes );
6859 LWT_EDGERING_ARRAY_CLEAN( &shells );
6860 lwerror("Errors fetching or registering face-missing edges: %s",
6861 lwt_be_lastErrorMessage(iface));
6862 return -1;
6863 }
6864
6865 LWDEBUGF(1, "Found %d holes and %d shells", holes.size, shells.size);
6866
6867 /* TODO: sort holes by pt.x, sort shells by bbox.xmin */
6868
6869 /* Assign shells to holes */
6870 for (i=0; i<holes.size; ++i)
6871 {
6872 LWT_ELEMID containing_face;
6873 LWT_EDGERING *ring = holes.rings[i];
6874
6875 containing_face = _lwt_FindFaceContainingRing(topo, ring, &shells);
6876 LWDEBUGF(1, "Ring %d contained by face %" LWTFMT_ELEMID, i, containing_face);
6877 if ( containing_face == -1 )
6878 {
6879 _lwt_release_edges(edgetable.edges, edgetable.size);
6880 LWT_EDGERING_ARRAY_CLEAN( &holes );
6881 LWT_EDGERING_ARRAY_CLEAN( &shells );
6882 lwerror("Errors finding face containing ring: %s",
6883 lwt_be_lastErrorMessage(iface));
6884 return -1;
6885 }
6886 int ret = _lwt_UpdateEdgeRingSideFace(topo, holes.rings[i], containing_face);
6887 if ( ret )
6888 {
6889 _lwt_release_edges(edgetable.edges, edgetable.size);
6890 LWT_EDGERING_ARRAY_CLEAN( &holes );
6891 LWT_EDGERING_ARRAY_CLEAN( &shells );
6892 lwerror("Errors updating edgering side face: %s",
6893 lwt_be_lastErrorMessage(iface));
6894 return -1;
6895 }
6896 }
6897
6898 LWDEBUG(1, "All holes assigned, cleaning up");
6899
6900 _lwt_release_edges(edgetable.edges, edgetable.size);
6901
6902 /* delete all shell and hole EDGERINGS */
6903 LWT_EDGERING_ARRAY_CLEAN( &holes );
6904 LWT_EDGERING_ARRAY_CLEAN( &shells );
6905
6906 return 0;
6907 }
6908