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