1 /*
2
3 gg_relations_ext.c -- Gaia spatial relations [advanced]
4
5 version 5.0, 2020 August 1
6
7 Author: Sandro Furieri a.furieri@lqt.it
8
9 ------------------------------------------------------------------------------
10
11 Version: MPL 1.1/GPL 2.0/LGPL 2.1
12
13 The contents of this file are subject to the Mozilla Public License Version
14 1.1 (the "License"); you may not use this file except in compliance with
15 the License. You may obtain a copy of the License at
16 http://www.mozilla.org/MPL/
17
18 Software distributed under the License is distributed on an "AS IS" basis,
19 WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
20 for the specific language governing rights and limitations under the
21 License.
22
23 The Original Code is the SpatiaLite library
24
25 The Initial Developer of the Original Code is Alessandro Furieri
26
27 Portions created by the Initial Developer are Copyright (C) 2008-2021
28 the Initial Developer. All Rights Reserved.
29
30 Contributor(s):
31
32 Alternatively, the contents of this file may be used under the terms of
33 either the GNU General Public License Version 2 or later (the "GPL"), or
34 the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
35 in which case the provisions of the GPL or the LGPL are applicable instead
36 of those above. If you wish to allow use of your version of this file only
37 under the terms of either the GPL or the LGPL, and not to allow others to
38 use your version of this file under the terms of the MPL, indicate your
39 decision by deleting the provisions above and replace them with the notice
40 and other provisions required by the GPL or the LGPL. If you do not delete
41 the provisions above, a recipient may use your version of this file under
42 the terms of any one of the MPL, the GPL or the LGPL.
43
44 */
45
46 #include <sys/types.h>
47 #include <stdlib.h>
48 #include <stdio.h>
49 #include <string.h>
50 #include <float.h>
51
52 #if defined(_WIN32) && !defined(__MINGW32__)
53 #include "config-msvc.h"
54 #else
55 #include "config.h"
56 #endif
57
58 #ifndef OMIT_GEOS /* including GEOS */
59 #ifdef GEOS_REENTRANT
60 #ifdef GEOS_ONLY_REENTRANT
61 #define GEOS_USE_ONLY_R_API /* only fully thread-safe GEOS API */
62 #endif
63 #endif
64 #include <geos_c.h>
65 #endif
66
67 #include <spatialite_private.h>
68 #include <spatialite/sqlite.h>
69
70 #include <spatialite/gaiageo.h>
71
72 #ifndef OMIT_GEOS /* including GEOS */
73
74 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaOffsetCurve(gaiaGeomCollPtr geom,double radius,int points,int left_right)75 gaiaOffsetCurve (gaiaGeomCollPtr geom, double radius, int points,
76 int left_right)
77 {
78 /*
79 // builds a geometry that is the OffsetCurve of GEOM
80 // (which is expected to be of the LINESTRING type)
81 //
82 */
83 gaiaGeomCollPtr geo = NULL;
84 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
85 GEOSGeometry *g1;
86 GEOSGeometry *g2;
87 gaiaPointPtr pt;
88 gaiaLinestringPtr ln;
89 gaiaPolygonPtr pg;
90 int pts = 0;
91 int lns = 0;
92 int pgs = 0;
93 int closed = 0;
94 gaiaResetGeosMsg ();
95 if (!geom)
96 return NULL;
97
98 if (left_right < 0)
99 left_right = 0; /* silencing stupid compiler warnings */
100
101 /* checking the input geometry for validity */
102 pt = geom->FirstPoint;
103 while (pt)
104 {
105 /* counting how many POINTs are there */
106 pts++;
107 pt = pt->Next;
108 }
109 ln = geom->FirstLinestring;
110 while (ln)
111 {
112 /* counting how many LINESTRINGs are there */
113 lns++;
114 if (gaiaIsClosed (ln))
115 closed++;
116 ln = ln->Next;
117 }
118 pg = geom->FirstPolygon;
119 while (pg)
120 {
121 /* counting how many POLYGON are there */
122 pgs++;
123 pg = pg->Next;
124 }
125 if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
126 return NULL;
127
128 /* all right: this one simply is a LINESTRING */
129 geom->DeclaredType = GAIA_LINESTRING;
130
131 g1 = gaiaToGeos (geom);
132 g2 = GEOSOffsetCurve (g1, radius, points, GEOSBUF_JOIN_ROUND, 5.0);
133 GEOSGeom_destroy (g1);
134 if (!g2)
135 return NULL;
136 if (geom->DimensionModel == GAIA_XY_Z)
137 geo = gaiaFromGeos_XYZ (g2);
138 else if (geom->DimensionModel == GAIA_XY_M)
139 geo = gaiaFromGeos_XYM (g2);
140 else if (geom->DimensionModel == GAIA_XY_Z_M)
141 geo = gaiaFromGeos_XYZM (g2);
142 else
143 geo = gaiaFromGeos_XY (g2);
144 GEOSGeom_destroy (g2);
145 if (geo == NULL)
146 return NULL;
147 geo->Srid = geom->Srid;
148 #else
149 if (geom == NULL || radius == 0.0 || points == 0 || left_right == 0)
150 geom = NULL; /* silencing stupid compiler warnings */
151 #endif
152 return geo;
153 }
154
155 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaOffsetCurve_r(const void * p_cache,gaiaGeomCollPtr geom,double radius,int points,int left_right)156 gaiaOffsetCurve_r (const void *p_cache, gaiaGeomCollPtr geom, double radius,
157 int points, int left_right)
158 {
159 /*
160 // builds a geometry that is the OffsetCurve of GEOM
161 // (which is expected to be of the LINESTRING type)
162 //
163 */
164 gaiaGeomCollPtr geo;
165 GEOSGeometry *g1;
166 GEOSGeometry *g2;
167 gaiaPointPtr pt;
168 gaiaLinestringPtr ln;
169 gaiaPolygonPtr pg;
170 int pts = 0;
171 int lns = 0;
172 int pgs = 0;
173 int closed = 0;
174 struct splite_internal_cache *cache =
175 (struct splite_internal_cache *) p_cache;
176 GEOSContextHandle_t handle = NULL;
177 if (cache == NULL)
178 return NULL;
179 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
180 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
181 return NULL;
182 handle = cache->GEOS_handle;
183 if (handle == NULL)
184 return NULL;
185 gaiaResetGeosMsg_r (cache);
186 if (!geom)
187 return NULL;
188
189 if (left_right < 0)
190 left_right = 0; /* silencing stupid compiler warnings */
191
192 /* checking the input geometry for validity */
193 pt = geom->FirstPoint;
194 while (pt)
195 {
196 /* counting how many POINTs are there */
197 pts++;
198 pt = pt->Next;
199 }
200 ln = geom->FirstLinestring;
201 while (ln)
202 {
203 /* counting how many LINESTRINGs are there */
204 lns++;
205 if (gaiaIsClosed (ln))
206 closed++;
207 ln = ln->Next;
208 }
209 pg = geom->FirstPolygon;
210 while (pg)
211 {
212 /* counting how many POLYGON are there */
213 pgs++;
214 pg = pg->Next;
215 }
216 if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
217 return NULL;
218
219 /* all right: this one simply is a LINESTRING */
220 geom->DeclaredType = GAIA_LINESTRING;
221
222 g1 = gaiaToGeos_r (cache, geom);
223 g2 = GEOSOffsetCurve_r (handle, g1, radius, points,
224 GEOSBUF_JOIN_ROUND, 5.0);
225 GEOSGeom_destroy_r (handle, g1);
226 if (!g2)
227 return NULL;
228 if (geom->DimensionModel == GAIA_XY_Z)
229 geo = gaiaFromGeos_XYZ_r (cache, g2);
230 else if (geom->DimensionModel == GAIA_XY_M)
231 geo = gaiaFromGeos_XYM_r (cache, g2);
232 else if (geom->DimensionModel == GAIA_XY_Z_M)
233 geo = gaiaFromGeos_XYZM_r (cache, g2);
234 else
235 geo = gaiaFromGeos_XY_r (cache, g2);
236 GEOSGeom_destroy_r (handle, g2);
237 if (geo == NULL)
238 return NULL;
239 geo->Srid = geom->Srid;
240 return geo;
241 }
242
243 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSingleSidedBuffer(gaiaGeomCollPtr geom,double radius,int points,int left_right)244 gaiaSingleSidedBuffer (gaiaGeomCollPtr geom, double radius, int points,
245 int left_right)
246 {
247 /*
248 // builds a geometry that is the SingleSided BUFFER of GEOM
249 // (which is expected to be of the LINESTRING type)
250 //
251 */
252 gaiaGeomCollPtr geo = NULL;
253 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
254 GEOSGeometry *g1;
255 GEOSGeometry *g2;
256 GEOSBufferParams *params = NULL;
257 gaiaPointPtr pt;
258 gaiaLinestringPtr ln;
259 gaiaPolygonPtr pg;
260 int pts = 0;
261 int lns = 0;
262 int pgs = 0;
263 int closed = 0;
264 gaiaResetGeosMsg ();
265 if (!geom)
266 return NULL;
267
268 /* checking the input geometry for validity */
269 pt = geom->FirstPoint;
270 while (pt)
271 {
272 /* counting how many POINTs are there */
273 pts++;
274 pt = pt->Next;
275 }
276 ln = geom->FirstLinestring;
277 while (ln)
278 {
279 /* counting how many LINESTRINGs are there */
280 lns++;
281 if (gaiaIsClosed (ln))
282 closed++;
283 ln = ln->Next;
284 }
285 pg = geom->FirstPolygon;
286 while (pg)
287 {
288 /* counting how many POLYGON are there */
289 pgs++;
290 pg = pg->Next;
291 }
292 if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
293 return NULL;
294
295 /* all right: this one simply is a LINESTRING */
296 geom->DeclaredType = GAIA_LINESTRING;
297
298 g1 = gaiaToGeos (geom);
299 /* setting up Buffer params */
300 params = GEOSBufferParams_create ();
301 GEOSBufferParams_setEndCapStyle (params, GEOSBUF_CAP_ROUND);
302 GEOSBufferParams_setJoinStyle (params, GEOSBUF_JOIN_ROUND);
303 GEOSBufferParams_setMitreLimit (params, 5.0);
304 GEOSBufferParams_setQuadrantSegments (params, points);
305 GEOSBufferParams_setSingleSided (params, 1);
306
307 /* creating the SingleSided Buffer */
308 if (left_right == 0)
309 {
310 /* right-sided requires NEGATIVE radius */
311 radius *= -1.0;
312 }
313 g2 = GEOSBufferWithParams (g1, params, radius);
314 GEOSGeom_destroy (g1);
315 GEOSBufferParams_destroy (params);
316 if (!g2)
317 return NULL;
318 if (geom->DimensionModel == GAIA_XY_Z)
319 geo = gaiaFromGeos_XYZ (g2);
320 else if (geom->DimensionModel == GAIA_XY_M)
321 geo = gaiaFromGeos_XYM (g2);
322 else if (geom->DimensionModel == GAIA_XY_Z_M)
323 geo = gaiaFromGeos_XYZM (g2);
324 else
325 geo = gaiaFromGeos_XY (g2);
326 GEOSGeom_destroy (g2);
327 if (geo == NULL)
328 return NULL;
329 geo->Srid = geom->Srid;
330 #else
331 if (geom == NULL || radius == 0.0 || points == 0 || left_right == 0)
332 geom = NULL; /* silencing stupid compiler warnings */
333 #endif
334 return geo;
335 }
336
337 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSingleSidedBuffer_r(const void * p_cache,gaiaGeomCollPtr geom,double radius,int points,int left_right)338 gaiaSingleSidedBuffer_r (const void *p_cache, gaiaGeomCollPtr geom,
339 double radius, int points, int left_right)
340 {
341 /*
342 // builds a geometry that is the SingleSided BUFFER of GEOM
343 // (which is expected to be of the LINESTRING type)
344 //
345 */
346 gaiaGeomCollPtr geo;
347 GEOSGeometry *g1;
348 GEOSGeometry *g2;
349 GEOSBufferParams *params = NULL;
350 int quadsegs = 30;
351 gaiaPointPtr pt;
352 gaiaLinestringPtr ln;
353 gaiaPolygonPtr pg;
354 int pts = 0;
355 int lns = 0;
356 int pgs = 0;
357 int closed = 0;
358 struct splite_internal_cache *cache =
359 (struct splite_internal_cache *) p_cache;
360 GEOSContextHandle_t handle = NULL;
361 if (cache == NULL)
362 return NULL;
363 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
364 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
365 return NULL;
366 handle = cache->GEOS_handle;
367 if (handle == NULL)
368 return NULL;
369 gaiaResetGeosMsg_r (cache);
370 if (!geom)
371 return NULL;
372
373 /* checking the input geometry for validity */
374 pt = geom->FirstPoint;
375 while (pt)
376 {
377 /* counting how many POINTs are there */
378 pts++;
379 pt = pt->Next;
380 }
381 ln = geom->FirstLinestring;
382 while (ln)
383 {
384 /* counting how many LINESTRINGs are there */
385 lns++;
386 if (gaiaIsClosed (ln))
387 closed++;
388 ln = ln->Next;
389 }
390 pg = geom->FirstPolygon;
391 while (pg)
392 {
393 /* counting how many POLYGON are there */
394 pgs++;
395 pg = pg->Next;
396 }
397 if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
398 return NULL;
399
400 /* all right: this one simply is a LINESTRING */
401 geom->DeclaredType = GAIA_LINESTRING;
402
403 g1 = gaiaToGeos_r (cache, geom);
404 /* setting up Buffer params */
405 params = GEOSBufferParams_create_r (handle);
406 GEOSBufferParams_setEndCapStyle_r (handle, params,
407 cache->buffer_end_cap_style);
408 GEOSBufferParams_setJoinStyle_r (handle, params, cache->buffer_join_style);
409 GEOSBufferParams_setMitreLimit_r (handle, params,
410 cache->buffer_mitre_limit);
411 if (points > 0)
412 quadsegs = points;
413 else if (cache->buffer_quadrant_segments > 0)
414 quadsegs = cache->buffer_quadrant_segments;
415 GEOSBufferParams_setQuadrantSegments_r (handle, params, quadsegs);
416 GEOSBufferParams_setSingleSided_r (handle, params, 1);
417
418 /* creating the SingleSided Buffer */
419 if (left_right == 0)
420 {
421 /* right-sided requires NEGATIVE radius */
422 radius *= -1.0;
423 }
424 g2 = GEOSBufferWithParams_r (handle, g1, params, radius);
425 GEOSGeom_destroy_r (handle, g1);
426 GEOSBufferParams_destroy_r (handle, params);
427 if (!g2)
428 return NULL;
429 if (geom->DimensionModel == GAIA_XY_Z)
430 geo = gaiaFromGeos_XYZ_r (cache, g2);
431 else if (geom->DimensionModel == GAIA_XY_M)
432 geo = gaiaFromGeos_XYM_r (cache, g2);
433 else if (geom->DimensionModel == GAIA_XY_Z_M)
434 geo = gaiaFromGeos_XYZM_r (cache, g2);
435 else
436 geo = gaiaFromGeos_XY_r (cache, g2);
437 GEOSGeom_destroy_r (handle, g2);
438 if (geo == NULL)
439 return NULL;
440 geo->Srid = geom->Srid;
441 return geo;
442 }
443
444 GAIAGEO_DECLARE int
gaiaHausdorffDistance(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double * xdist)445 gaiaHausdorffDistance (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
446 double *xdist)
447 {
448 /*
449 / computes the (discrete) Hausdorff distance intercurring
450 / between GEOM-1 and GEOM-2
451 */
452 int ret = 0;
453 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
454 double dist;
455 GEOSGeometry *g1;
456 GEOSGeometry *g2;
457 gaiaResetGeosMsg ();
458 if (!geom1 || !geom2)
459 return 0;
460 g1 = gaiaToGeos (geom1);
461 g2 = gaiaToGeos (geom2);
462 ret = GEOSHausdorffDistance (g1, g2, &dist);
463 GEOSGeom_destroy (g1);
464 GEOSGeom_destroy (g2);
465 if (ret)
466 *xdist = dist;
467 #else
468 if (geom1 == NULL || geom2 == NULL || xdist == NULL)
469 geom1 = NULL; /* silencing stupid compiler warnings */
470 #endif
471 return ret;
472 }
473
474 GAIAGEO_DECLARE int
gaiaHausdorffDistance_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double * xdist)475 gaiaHausdorffDistance_r (const void *p_cache, gaiaGeomCollPtr geom1,
476 gaiaGeomCollPtr geom2, double *xdist)
477 {
478 /*
479 / computes the (discrete) Hausdorff distance intercurring
480 / between GEOM-1 and GEOM-2
481 */
482 double dist;
483 int ret;
484 GEOSGeometry *g1;
485 GEOSGeometry *g2;
486 struct splite_internal_cache *cache =
487 (struct splite_internal_cache *) p_cache;
488 GEOSContextHandle_t handle = NULL;
489 if (cache == NULL)
490 return 0;
491 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
492 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
493 return 0;
494 handle = cache->GEOS_handle;
495 if (handle == NULL)
496 return 0;
497 gaiaResetGeosMsg_r (cache);
498 if (!geom1 || !geom2)
499 return 0;
500 g1 = gaiaToGeos_r (cache, geom1);
501 g2 = gaiaToGeos_r (cache, geom2);
502 ret = GEOSHausdorffDistance_r (handle, g1, g2, &dist);
503 GEOSGeom_destroy_r (handle, g1);
504 GEOSGeom_destroy_r (handle, g2);
505 if (ret)
506 *xdist = dist;
507 return ret;
508 }
509
510 #ifdef GEOS_370 /* only if GEOS_370 support is available */
511
512 GAIAGEO_DECLARE int
gaiaHausdorffDistanceDensify(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double densify_fract,double * xdist)513 gaiaHausdorffDistanceDensify (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
514 double densify_fract, double *xdist)
515 {
516 /*
517 / computes the (discrete) Hausdorff distance intercurring
518 / between GEOM-1 and GEOM-2
519 */
520 int ret = 0;
521 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
522 double dist;
523 GEOSGeometry *g1;
524 GEOSGeometry *g2;
525 gaiaResetGeosMsg ();
526 if (!geom1 || !geom2)
527 return 0;
528 g1 = gaiaToGeos (geom1);
529 g2 = gaiaToGeos (geom2);
530 ret = GEOSHausdorffDistanceDensify (g1, g2, densify_fract, &dist);
531 GEOSGeom_destroy (g1);
532 GEOSGeom_destroy (g2);
533 if (ret)
534 *xdist = dist;
535 #else
536 if (geom1 == NULL || geom2 == NULL || xdist == NULL)
537 geom1 = NULL; /* silencing stupid compiler warnings */
538 #endif
539 return ret;
540 }
541
542 GAIAGEO_DECLARE int
gaiaHausdorffDistanceDensify_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double densify_fract,double * xdist)543 gaiaHausdorffDistanceDensify_r (const void *p_cache, gaiaGeomCollPtr geom1,
544 gaiaGeomCollPtr geom2, double densify_fract,
545 double *xdist)
546 {
547 /*
548 / computes the (discrete) Hausdorff distance intercurring
549 / between GEOM-1 and GEOM-2
550 */
551 double dist;
552 int ret;
553 GEOSGeometry *g1;
554 GEOSGeometry *g2;
555 struct splite_internal_cache *cache =
556 (struct splite_internal_cache *) p_cache;
557 GEOSContextHandle_t handle = NULL;
558 if (cache == NULL)
559 return 0;
560 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
561 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
562 return 0;
563 handle = cache->GEOS_handle;
564 if (handle == NULL)
565 return 0;
566 gaiaResetGeosMsg_r (cache);
567 if (!geom1 || !geom2)
568 return 0;
569 g1 = gaiaToGeos_r (cache, geom1);
570 g2 = gaiaToGeos_r (cache, geom2);
571 ret = GEOSHausdorffDistanceDensify_r (handle, g1, g2, densify_fract, &dist);
572 GEOSGeom_destroy_r (handle, g1);
573 GEOSGeom_destroy_r (handle, g2);
574 if (ret)
575 *xdist = dist;
576 return ret;
577 }
578
579 GAIAGEO_DECLARE int
gaiaFrechetDistance(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double * xdist)580 gaiaFrechetDistance (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
581 double *xdist)
582 {
583 /*
584 / computes the (discrete) Frechet distance intercurring
585 / between GEOM-1 and GEOM-2
586 */
587 int ret = 0;
588 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
589 double dist;
590 GEOSGeometry *g1;
591 GEOSGeometry *g2;
592 gaiaResetGeosMsg ();
593 if (!geom1 || !geom2)
594 return 0;
595 g1 = gaiaToGeos (geom1);
596 g2 = gaiaToGeos (geom2);
597 ret = GEOSFrechetDistance (g1, g2, &dist);
598 GEOSGeom_destroy (g1);
599 GEOSGeom_destroy (g2);
600 if (ret)
601 *xdist = dist;
602 #else
603 if (geom1 == NULL || geom2 == NULL || xdist == NULL)
604 geom1 = NULL; /* silencing stupid compiler warnings */
605 #endif
606 return ret;
607 }
608
609 GAIAGEO_DECLARE int
gaiaFrechetDistance_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double * xdist)610 gaiaFrechetDistance_r (const void *p_cache, gaiaGeomCollPtr geom1,
611 gaiaGeomCollPtr geom2, double *xdist)
612 {
613 /*
614 / computes the (discrete) Frechet distance intercurring
615 / between GEOM-1 and GEOM-2
616 */
617 double dist;
618 int ret;
619 GEOSGeometry *g1;
620 GEOSGeometry *g2;
621 struct splite_internal_cache *cache =
622 (struct splite_internal_cache *) p_cache;
623 GEOSContextHandle_t handle = NULL;
624 if (cache == NULL)
625 return 0;
626 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
627 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
628 return 0;
629 handle = cache->GEOS_handle;
630 if (handle == NULL)
631 return 0;
632 gaiaResetGeosMsg_r (cache);
633 if (!geom1 || !geom2)
634 return 0;
635 g1 = gaiaToGeos_r (cache, geom1);
636 g2 = gaiaToGeos_r (cache, geom2);
637 ret = GEOSFrechetDistance_r (handle, g1, g2, &dist);
638 GEOSGeom_destroy_r (handle, g1);
639 GEOSGeom_destroy_r (handle, g2);
640 if (ret)
641 *xdist = dist;
642 return ret;
643 }
644
645 GAIAGEO_DECLARE int
gaiaFrechetDistanceDensify(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double densify_fract,double * xdist)646 gaiaFrechetDistanceDensify (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
647 double densify_fract, double *xdist)
648 {
649 /*
650 / computes the (discrete) Frechet distance intercurring
651 / between GEOM-1 and GEOM-2
652 */
653 int ret = 0;
654 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
655 double dist;
656 GEOSGeometry *g1;
657 GEOSGeometry *g2;
658 gaiaResetGeosMsg ();
659 if (!geom1 || !geom2)
660 return 0;
661 g1 = gaiaToGeos (geom1);
662 g2 = gaiaToGeos (geom2);
663 ret = GEOSFrechetDistanceDensify (g1, g2, densify_fract, &dist);
664 GEOSGeom_destroy (g1);
665 GEOSGeom_destroy (g2);
666 if (ret)
667 *xdist = dist;
668 #else
669 if (geom1 == NULL || geom2 == NULL || xdist == NULL)
670 geom1 = NULL; /* silencing stupid compiler warnings */
671 #endif
672 return ret;
673 }
674
675 GAIAGEO_DECLARE int
gaiaFrechetDistanceDensify_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double densify_fract,double * xdist)676 gaiaFrechetDistanceDensify_r (const void *p_cache, gaiaGeomCollPtr geom1,
677 gaiaGeomCollPtr geom2, double densify_fract,
678 double *xdist)
679 {
680 /*
681 / computes the (discrete) Frechet distance intercurring
682 / between GEOM-1 and GEOM-2
683 */
684 double dist;
685 int ret;
686 GEOSGeometry *g1;
687 GEOSGeometry *g2;
688 struct splite_internal_cache *cache =
689 (struct splite_internal_cache *) p_cache;
690 GEOSContextHandle_t handle = NULL;
691 if (cache == NULL)
692 return 0;
693 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
694 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
695 return 0;
696 handle = cache->GEOS_handle;
697 if (handle == NULL)
698 return 0;
699 gaiaResetGeosMsg_r (cache);
700 if (!geom1 || !geom2)
701 return 0;
702 g1 = gaiaToGeos_r (cache, geom1);
703 g2 = gaiaToGeos_r (cache, geom2);
704 ret = GEOSFrechetDistanceDensify_r (handle, g1, g2, densify_fract, &dist);
705 GEOSGeom_destroy_r (handle, g1);
706 GEOSGeom_destroy_r (handle, g2);
707 if (ret)
708 *xdist = dist;
709 return ret;
710 }
711
712 #endif /* end GEOS_370 conditional */
713
714 static gaiaGeomCollPtr
geom_as_lines(gaiaGeomCollPtr geom)715 geom_as_lines (gaiaGeomCollPtr geom)
716 {
717 /* transforms a Geometry into a LINESTRING/MULTILINESTRING (if possible) */
718 gaiaGeomCollPtr result;
719 gaiaLinestringPtr ln;
720 gaiaLinestringPtr new_ln;
721 gaiaPolygonPtr pg;
722 gaiaRingPtr rng;
723 int iv;
724 int ib;
725 double x;
726 double y;
727 double z;
728 double m;
729
730 if (!geom)
731 return NULL;
732 if (geom->FirstPoint != NULL)
733 {
734 /* invalid: GEOM contains at least one POINT */
735 return NULL;
736 }
737
738 switch (geom->DimensionModel)
739 {
740 case GAIA_XY_Z_M:
741 result = gaiaAllocGeomCollXYZM ();
742 break;
743 case GAIA_XY_Z:
744 result = gaiaAllocGeomCollXYZ ();
745 break;
746 case GAIA_XY_M:
747 result = gaiaAllocGeomCollXYM ();
748 break;
749 default:
750 result = gaiaAllocGeomColl ();
751 break;
752 };
753 result->Srid = geom->Srid;
754 ln = geom->FirstLinestring;
755 while (ln)
756 {
757 /* copying any Linestring */
758 new_ln = gaiaAddLinestringToGeomColl (result, ln->Points);
759 for (iv = 0; iv < ln->Points; iv++)
760 {
761 if (ln->DimensionModel == GAIA_XY_Z)
762 {
763 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
764 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
765 }
766 else if (ln->DimensionModel == GAIA_XY_M)
767 {
768 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
769 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
770 }
771 else if (ln->DimensionModel == GAIA_XY_Z_M)
772 {
773 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
774 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
775 }
776 else
777 {
778 gaiaGetPoint (ln->Coords, iv, &x, &y);
779 gaiaSetPoint (new_ln->Coords, iv, x, y);
780 }
781 }
782 ln = ln->Next;
783 }
784 pg = geom->FirstPolygon;
785 while (pg)
786 {
787 /* copying any Polygon Ring (as Linestring) */
788 rng = pg->Exterior;
789 new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
790 for (iv = 0; iv < rng->Points; iv++)
791 {
792 /* exterior Ring */
793 if (rng->DimensionModel == GAIA_XY_Z)
794 {
795 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
796 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
797 }
798 else if (rng->DimensionModel == GAIA_XY_M)
799 {
800 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
801 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
802 }
803 else if (rng->DimensionModel == GAIA_XY_Z_M)
804 {
805 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
806 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
807 }
808 else
809 {
810 gaiaGetPoint (rng->Coords, iv, &x, &y);
811 gaiaSetPoint (new_ln->Coords, iv, x, y);
812 }
813 }
814 for (ib = 0; ib < pg->NumInteriors; ib++)
815 {
816 rng = pg->Interiors + ib;
817 new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
818 for (iv = 0; iv < rng->Points; iv++)
819 {
820 /* any interior Ring */
821 if (rng->DimensionModel == GAIA_XY_Z)
822 {
823 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
824 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
825 }
826 else if (rng->DimensionModel == GAIA_XY_M)
827 {
828 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
829 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
830 }
831 else if (rng->DimensionModel == GAIA_XY_Z_M)
832 {
833 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
834 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
835 }
836 else
837 {
838 gaiaGetPoint (rng->Coords, iv, &x, &y);
839 gaiaSetPoint (new_ln->Coords, iv, x, y);
840 }
841 }
842 }
843 pg = pg->Next;
844 }
845 return result;
846 }
847
848 static void
add_shared_linestring(gaiaGeomCollPtr geom,gaiaDynamicLinePtr dyn)849 add_shared_linestring (gaiaGeomCollPtr geom, gaiaDynamicLinePtr dyn)
850 {
851 /* adding a LINESTRING from Dynamic Line */
852 int count = 0;
853 gaiaLinestringPtr ln;
854 gaiaPointPtr pt;
855 int iv;
856
857 if (!geom)
858 return;
859 if (!dyn)
860 return;
861 pt = dyn->First;
862 while (pt)
863 {
864 /* counting how many Points are there */
865 count++;
866 pt = pt->Next;
867 }
868 if (count == 0)
869 return;
870 ln = gaiaAddLinestringToGeomColl (geom, count);
871 iv = 0;
872 pt = dyn->First;
873 while (pt)
874 {
875 /* copying points into the LINESTRING */
876 if (ln->DimensionModel == GAIA_XY_Z)
877 {
878 gaiaSetPointXYZ (ln->Coords, iv, pt->X, pt->Y, pt->Z);
879 }
880 else if (ln->DimensionModel == GAIA_XY_M)
881 {
882 gaiaSetPointXYM (ln->Coords, iv, pt->X, pt->Y, pt->M);
883 }
884 else if (ln->DimensionModel == GAIA_XY_Z_M)
885 {
886 gaiaSetPointXYZM (ln->Coords, iv, pt->X, pt->Y, pt->Z, pt->M);
887 }
888 else
889 {
890 gaiaSetPoint (ln->Coords, iv, pt->X, pt->Y);
891 }
892 iv++;
893 pt = pt->Next;
894 }
895 }
896
897 static void
append_shared_path(gaiaDynamicLinePtr dyn,gaiaLinestringPtr ln,int order)898 append_shared_path (gaiaDynamicLinePtr dyn, gaiaLinestringPtr ln, int order)
899 {
900 /* appends a Shared Path item to Dynamic Line */
901 int iv;
902 double x;
903 double y;
904 double z;
905 double m;
906 if (order)
907 {
908 /* reversed order */
909 for (iv = ln->Points - 1; iv >= 0; iv--)
910 {
911 if (ln->DimensionModel == GAIA_XY_Z)
912 {
913 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
914 if (x == dyn->Last->X && y == dyn->Last->Y
915 && z == dyn->Last->Z)
916 ;
917 else
918 gaiaAppendPointZToDynamicLine (dyn, x, y, z);
919 }
920 else if (ln->DimensionModel == GAIA_XY_M)
921 {
922 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
923 if (x == dyn->Last->X && y == dyn->Last->Y
924 && m == dyn->Last->M)
925 ;
926 else
927 gaiaAppendPointMToDynamicLine (dyn, x, y, m);
928 }
929 else if (ln->DimensionModel == GAIA_XY_Z_M)
930 {
931 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
932 if (x == dyn->Last->X && y == dyn->Last->Y
933 && z == dyn->Last->Z && m == dyn->Last->M)
934 ;
935 else
936 gaiaAppendPointZMToDynamicLine (dyn, x, y, z, m);
937 }
938 else
939 {
940 gaiaGetPoint (ln->Coords, iv, &x, &y);
941 if (x == dyn->Last->X && y == dyn->Last->Y)
942 ;
943 else
944 gaiaAppendPointToDynamicLine (dyn, x, y);
945 }
946 }
947 }
948 else
949 {
950 /* conformant order */
951 for (iv = 0; iv < ln->Points; iv++)
952 {
953 if (ln->DimensionModel == GAIA_XY_Z)
954 {
955 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
956 if (x == dyn->Last->X && y == dyn->Last->Y
957 && z == dyn->Last->Z)
958 ;
959 else
960 gaiaAppendPointZToDynamicLine (dyn, x, y, z);
961 }
962 else if (ln->DimensionModel == GAIA_XY_M)
963 {
964 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
965 if (x == dyn->Last->X && y == dyn->Last->Y
966 && m == dyn->Last->M)
967 ;
968 else
969 gaiaAppendPointMToDynamicLine (dyn, x, y, m);
970 }
971 else if (ln->DimensionModel == GAIA_XY_Z_M)
972 {
973 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
974 if (x == dyn->Last->X && y == dyn->Last->Y
975 && z == dyn->Last->Z && m == dyn->Last->M)
976 ;
977 else
978 gaiaAppendPointZMToDynamicLine (dyn, x, y, z, m);
979 }
980 else
981 {
982 gaiaGetPoint (ln->Coords, iv, &x, &y);
983 if (x == dyn->Last->X && y == dyn->Last->Y)
984 ;
985 else
986 gaiaAppendPointToDynamicLine (dyn, x, y);
987 }
988 }
989 }
990 }
991
992 static void
prepend_shared_path(gaiaDynamicLinePtr dyn,gaiaLinestringPtr ln,int order)993 prepend_shared_path (gaiaDynamicLinePtr dyn, gaiaLinestringPtr ln, int order)
994 {
995 /* prepends a Shared Path item to Dynamic Line */
996 int iv;
997 double x;
998 double y;
999 double z;
1000 double m;
1001 if (order)
1002 {
1003 /* reversed order */
1004 for (iv = 0; iv < ln->Points; iv++)
1005 {
1006 if (ln->DimensionModel == GAIA_XY_Z)
1007 {
1008 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
1009 if (x == dyn->First->X && y == dyn->First->Y
1010 && z == dyn->First->Z)
1011 ;
1012 else
1013 gaiaPrependPointZToDynamicLine (dyn, x, y, z);
1014 }
1015 else if (ln->DimensionModel == GAIA_XY_M)
1016 {
1017 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
1018 if (x == dyn->First->X && y == dyn->First->Y
1019 && m == dyn->First->M)
1020 ;
1021 else
1022 gaiaPrependPointMToDynamicLine (dyn, x, y, m);
1023 }
1024 else if (ln->DimensionModel == GAIA_XY_Z_M)
1025 {
1026 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
1027 if (x == dyn->First->X && y == dyn->First->Y
1028 && z == dyn->First->Z && m == dyn->First->M)
1029 ;
1030 else
1031 gaiaPrependPointZMToDynamicLine (dyn, x, y, z, m);
1032 }
1033 else
1034 {
1035 gaiaGetPoint (ln->Coords, iv, &x, &y);
1036 if (x == dyn->First->X && y == dyn->First->Y)
1037 ;
1038 else
1039 gaiaPrependPointToDynamicLine (dyn, x, y);
1040 }
1041 }
1042 }
1043 else
1044 {
1045 /* conformant order */
1046 for (iv = ln->Points - 1; iv >= 0; iv--)
1047 {
1048 if (ln->DimensionModel == GAIA_XY_Z)
1049 {
1050 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
1051 if (x == dyn->First->X && y == dyn->First->Y
1052 && z == dyn->First->Z)
1053 ;
1054 else
1055 gaiaPrependPointZToDynamicLine (dyn, x, y, z);
1056 }
1057 else if (ln->DimensionModel == GAIA_XY_M)
1058 {
1059 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
1060 if (x == dyn->First->X && y == dyn->First->Y
1061 && m == dyn->First->M)
1062 ;
1063 else
1064 gaiaPrependPointMToDynamicLine (dyn, x, y, m);
1065 }
1066 else if (ln->DimensionModel == GAIA_XY_Z_M)
1067 {
1068 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
1069 if (x == dyn->First->X && y == dyn->First->Y
1070 && z == dyn->First->Z && m == dyn->First->M)
1071 ;
1072 else
1073 gaiaPrependPointZMToDynamicLine (dyn, x, y, z, m);
1074 }
1075 else
1076 {
1077 gaiaGetPoint (ln->Coords, iv, &x, &y);
1078 if (x == dyn->First->X && y == dyn->First->Y)
1079 ;
1080 else
1081 gaiaPrependPointToDynamicLine (dyn, x, y);
1082 }
1083 }
1084 }
1085 }
1086
1087 static gaiaGeomCollPtr
arrange_shared_paths(gaiaGeomCollPtr geom)1088 arrange_shared_paths (gaiaGeomCollPtr geom)
1089 {
1090 /* final aggregation step for shared paths */
1091 gaiaLinestringPtr ln;
1092 gaiaLinestringPtr *ln_array;
1093 gaiaGeomCollPtr result;
1094 gaiaDynamicLinePtr dyn;
1095 int count;
1096 int i;
1097 int i2;
1098 int iv;
1099 double x;
1100 double y;
1101 double z;
1102 double m;
1103 int ok;
1104 int ok2;
1105
1106 if (!geom)
1107 return NULL;
1108 count = 0;
1109 ln = geom->FirstLinestring;
1110 while (ln)
1111 {
1112 /* counting how many Linestrings are there */
1113 count++;
1114 ln = ln->Next;
1115 }
1116 if (count == 0)
1117 return NULL;
1118
1119 ln_array = malloc (sizeof (gaiaLinestringPtr) * count);
1120 i = 0;
1121 ln = geom->FirstLinestring;
1122 while (ln)
1123 {
1124 /* populating the Linestring references array */
1125 ln_array[i++] = ln;
1126 ln = ln->Next;
1127 }
1128
1129 /* allocating a new Geometry [MULTILINESTRING] */
1130 switch (geom->DimensionModel)
1131 {
1132 case GAIA_XY_Z_M:
1133 result = gaiaAllocGeomCollXYZM ();
1134 break;
1135 case GAIA_XY_Z:
1136 result = gaiaAllocGeomCollXYZ ();
1137 break;
1138 case GAIA_XY_M:
1139 result = gaiaAllocGeomCollXYM ();
1140 break;
1141 default:
1142 result = gaiaAllocGeomColl ();
1143 break;
1144 };
1145 result->Srid = geom->Srid;
1146 result->DeclaredType = GAIA_MULTILINESTRING;
1147
1148 ok = 1;
1149 while (ok)
1150 {
1151 /* looping until we have processed any input item */
1152 ok = 0;
1153 for (i = 0; i < count; i++)
1154 {
1155 if (ln_array[i] != NULL)
1156 {
1157 /* starting a new LINESTRING */
1158 dyn = gaiaAllocDynamicLine ();
1159 ln = ln_array[i];
1160 ln_array[i] = NULL;
1161 for (iv = 0; iv < ln->Points; iv++)
1162 {
1163 /* inserting the 'seed' path */
1164 if (ln->DimensionModel == GAIA_XY_Z)
1165 {
1166 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
1167 gaiaAppendPointZToDynamicLine (dyn, x, y, z);
1168 }
1169 else if (ln->DimensionModel == GAIA_XY_M)
1170 {
1171 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
1172 gaiaAppendPointMToDynamicLine (dyn, x, y, m);
1173 }
1174 else if (ln->DimensionModel == GAIA_XY_Z_M)
1175 {
1176 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z,
1177 &m);
1178 gaiaAppendPointZMToDynamicLine (dyn, x, y, z,
1179 m);
1180 }
1181 else
1182 {
1183 gaiaGetPoint (ln->Coords, iv, &x, &y);
1184 gaiaAppendPointToDynamicLine (dyn, x, y);
1185 }
1186 }
1187 ok2 = 1;
1188 while (ok2)
1189 {
1190 /* looping until we have checked any other item */
1191 ok2 = 0;
1192 for (i2 = 0; i2 < count; i2++)
1193 {
1194 /* expanding the 'seed' path */
1195 if (ln_array[i2] == NULL)
1196 continue;
1197 ln = ln_array[i2];
1198 /* checking the first vertex */
1199 iv = 0;
1200 if (ln->DimensionModel == GAIA_XY_Z)
1201 {
1202 gaiaGetPointXYZ (ln->Coords, iv, &x, &y,
1203 &z);
1204 }
1205 else if (ln->DimensionModel == GAIA_XY_M)
1206 {
1207 gaiaGetPointXYM (ln->Coords, iv, &x, &y,
1208 &m);
1209 }
1210 else if (ln->DimensionModel == GAIA_XY_Z_M)
1211 {
1212 gaiaGetPointXYZM (ln->Coords, iv, &x,
1213 &y, &z, &m);
1214 }
1215 else
1216 {
1217 gaiaGetPoint (ln->Coords, iv, &x, &y);
1218 }
1219 if (x == dyn->Last->X && y == dyn->Last->Y)
1220 {
1221 /* appending this item to the 'seed' (conformant order) */
1222 append_shared_path (dyn, ln, 0);
1223 ln_array[i2] = NULL;
1224 ok2 = 1;
1225 continue;
1226 }
1227 if (x == dyn->First->X && y == dyn->First->Y)
1228 {
1229 /* prepending this item to the 'seed' (reversed order) */
1230 prepend_shared_path (dyn, ln, 1);
1231 ln_array[i2] = NULL;
1232 ok2 = 1;
1233 continue;
1234 }
1235 /* checking the last vertex */
1236 iv = ln->Points - 1;
1237 if (ln->DimensionModel == GAIA_XY_Z)
1238 {
1239 gaiaGetPointXYZ (ln->Coords, iv, &x, &y,
1240 &z);
1241 }
1242 else if (ln->DimensionModel == GAIA_XY_M)
1243 {
1244 gaiaGetPointXYM (ln->Coords, iv, &x, &y,
1245 &m);
1246 }
1247 else if (ln->DimensionModel == GAIA_XY_Z_M)
1248 {
1249 gaiaGetPointXYZM (ln->Coords, iv, &x,
1250 &y, &z, &m);
1251 }
1252 else
1253 {
1254 gaiaGetPoint (ln->Coords, iv, &x, &y);
1255 }
1256 if (x == dyn->Last->X && y == dyn->Last->Y)
1257 {
1258 /* appending this item to the 'seed' (reversed order) */
1259 append_shared_path (dyn, ln, 1);
1260 ln_array[i2] = NULL;
1261 ok2 = 1;
1262 continue;
1263 }
1264 if (x == dyn->First->X && y == dyn->First->Y)
1265 {
1266 /* prepending this item to the 'seed' (conformant order) */
1267 prepend_shared_path (dyn, ln, 0);
1268 ln_array[i2] = NULL;
1269 ok2 = 1;
1270 continue;
1271 }
1272 }
1273 }
1274 add_shared_linestring (result, dyn);
1275 gaiaFreeDynamicLine (dyn);
1276 ok = 1;
1277 break;
1278 }
1279 }
1280 }
1281 free (ln_array);
1282 return result;
1283 }
1284
1285 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSharedPaths(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)1286 gaiaSharedPaths (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1287 {
1288 /*
1289 // builds a geometry containing Shared Paths commons to GEOM1 & GEOM2
1290 // (which are expected to be of the LINESTRING/MULTILINESTRING type)
1291 //
1292 */
1293 gaiaGeomCollPtr result = NULL;
1294 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1295 gaiaGeomCollPtr geo;
1296 gaiaGeomCollPtr line1;
1297 gaiaGeomCollPtr line2;
1298 GEOSGeometry *g1;
1299 GEOSGeometry *g2;
1300 GEOSGeometry *g3;
1301 gaiaResetGeosMsg ();
1302 if (!geom1)
1303 return NULL;
1304 if (!geom2)
1305 return NULL;
1306 /* transforming input geoms as Lines */
1307 line1 = geom_as_lines (geom1);
1308 line2 = geom_as_lines (geom2);
1309 if (line1 == NULL || line2 == NULL)
1310 {
1311 if (line1)
1312 gaiaFreeGeomColl (line1);
1313 if (line2)
1314 gaiaFreeGeomColl (line2);
1315 return NULL;
1316 }
1317
1318 g1 = gaiaToGeos (line1);
1319 g2 = gaiaToGeos (line2);
1320 gaiaFreeGeomColl (line1);
1321 gaiaFreeGeomColl (line2);
1322 g3 = GEOSSharedPaths (g1, g2);
1323 GEOSGeom_destroy (g1);
1324 GEOSGeom_destroy (g2);
1325 if (!g3)
1326 return NULL;
1327 if (geom1->DimensionModel == GAIA_XY_Z)
1328 geo = gaiaFromGeos_XYZ (g3);
1329 else if (geom1->DimensionModel == GAIA_XY_M)
1330 geo = gaiaFromGeos_XYM (g3);
1331 else if (geom1->DimensionModel == GAIA_XY_Z_M)
1332 geo = gaiaFromGeos_XYZM (g3);
1333 else
1334 geo = gaiaFromGeos_XY (g3);
1335 GEOSGeom_destroy (g3);
1336 if (geo == NULL)
1337 return NULL;
1338 geo->Srid = geom1->Srid;
1339 result = arrange_shared_paths (geo);
1340 gaiaFreeGeomColl (geo);
1341 #else
1342 if (geom1 == NULL || geom2 == NULL)
1343 geom1 = NULL; /* silencing stupid compiler warnings */
1344 #endif
1345 return result;
1346 }
1347
1348 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSharedPaths_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)1349 gaiaSharedPaths_r (const void *p_cache, gaiaGeomCollPtr geom1,
1350 gaiaGeomCollPtr geom2)
1351 {
1352 /*
1353 // builds a geometry containing Shared Paths commons to GEOM1 & GEOM2
1354 // (which are expected to be of the LINESTRING/MULTILINESTRING type)
1355 //
1356 */
1357 gaiaGeomCollPtr geo;
1358 gaiaGeomCollPtr result;
1359 gaiaGeomCollPtr line1;
1360 gaiaGeomCollPtr line2;
1361 GEOSGeometry *g1;
1362 GEOSGeometry *g2;
1363 GEOSGeometry *g3;
1364 struct splite_internal_cache *cache =
1365 (struct splite_internal_cache *) p_cache;
1366 GEOSContextHandle_t handle = NULL;
1367 if (cache == NULL)
1368 return NULL;
1369 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1370 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1371 return NULL;
1372 handle = cache->GEOS_handle;
1373 if (handle == NULL)
1374 return NULL;
1375 gaiaResetGeosMsg_r (cache);
1376 if (!geom1)
1377 return NULL;
1378 if (!geom2)
1379 return NULL;
1380 /* transforming input geoms as Lines */
1381 line1 = geom_as_lines (geom1);
1382 line2 = geom_as_lines (geom2);
1383 if (line1 == NULL || line2 == NULL)
1384 {
1385 if (line1)
1386 gaiaFreeGeomColl (line1);
1387 if (line2)
1388 gaiaFreeGeomColl (line2);
1389 return NULL;
1390 }
1391
1392 g1 = gaiaToGeos_r (cache, line1);
1393 g2 = gaiaToGeos_r (cache, line2);
1394 gaiaFreeGeomColl (line1);
1395 gaiaFreeGeomColl (line2);
1396 g3 = GEOSSharedPaths_r (handle, g1, g2);
1397 GEOSGeom_destroy_r (handle, g1);
1398 GEOSGeom_destroy_r (handle, g2);
1399 if (!g3)
1400 return NULL;
1401 if (geom1->DimensionModel == GAIA_XY_Z)
1402 geo = gaiaFromGeos_XYZ_r (cache, g3);
1403 else if (geom1->DimensionModel == GAIA_XY_M)
1404 geo = gaiaFromGeos_XYM_r (cache, g3);
1405 else if (geom1->DimensionModel == GAIA_XY_Z_M)
1406 geo = gaiaFromGeos_XYZM_r (cache, g3);
1407 else
1408 geo = gaiaFromGeos_XY_r (cache, g3);
1409 GEOSGeom_destroy_r (handle, g3);
1410 if (geo == NULL)
1411 return NULL;
1412 geo->Srid = geom1->Srid;
1413 result = arrange_shared_paths (geo);
1414 gaiaFreeGeomColl (geo);
1415 return result;
1416 }
1417
1418 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineInterpolatePoint(gaiaGeomCollPtr geom,double fraction)1419 gaiaLineInterpolatePoint (gaiaGeomCollPtr geom, double fraction)
1420 {
1421 /*
1422 * attempts to intepolate a point on line at dist "fraction"
1423 *
1424 * the fraction is expressed into the range from 0.0 to 1.0
1425 */
1426 gaiaGeomCollPtr result = NULL;
1427 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1428 int pts = 0;
1429 int lns = 0;
1430 int pgs = 0;
1431 gaiaPointPtr pt;
1432 gaiaLinestringPtr ln;
1433 gaiaPolygonPtr pg;
1434 GEOSGeometry *g;
1435 GEOSGeometry *g_pt;
1436 double length;
1437 double projection;
1438 gaiaResetGeosMsg ();
1439 if (!geom)
1440 return NULL;
1441
1442 /* checking if a single Linestring has been passed */
1443 pt = geom->FirstPoint;
1444 while (pt)
1445 {
1446 pts++;
1447 pt = pt->Next;
1448 }
1449 ln = geom->FirstLinestring;
1450 while (ln)
1451 {
1452 lns++;
1453 ln = ln->Next;
1454 }
1455 pg = geom->FirstPolygon;
1456 while (pg)
1457 {
1458 pgs++;
1459 pg = pg->Next;
1460 }
1461 if (pts == 0 && lns == 1 && pgs == 0)
1462 ;
1463 else
1464 return NULL;
1465
1466 g = gaiaToGeos (geom);
1467 if (GEOSLength (g, &length))
1468 {
1469 /* transforming fraction to length */
1470 if (fraction < 0.0)
1471 fraction = 0.0;
1472 if (fraction > 1.0)
1473 fraction = 1.0;
1474 projection = length * fraction;
1475 }
1476 else
1477 {
1478 GEOSGeom_destroy (g);
1479 return NULL;
1480 }
1481 g_pt = GEOSInterpolate (g, projection);
1482 GEOSGeom_destroy (g);
1483 if (!g_pt)
1484 return NULL;
1485 if (geom->DimensionModel == GAIA_XY_Z)
1486 result = gaiaFromGeos_XYZ (g_pt);
1487 else if (geom->DimensionModel == GAIA_XY_M)
1488 result = gaiaFromGeos_XYM (g_pt);
1489 else if (geom->DimensionModel == GAIA_XY_Z_M)
1490 result = gaiaFromGeos_XYZM (g_pt);
1491 else
1492 result = gaiaFromGeos_XY (g_pt);
1493 GEOSGeom_destroy (g_pt);
1494 if (result == NULL)
1495 return NULL;
1496 result->Srid = geom->Srid;
1497 #else
1498 if (geom == NULL || fraction == 0.0)
1499 geom = NULL; /* silencing stupid compiler warnings */
1500 #endif
1501 return result;
1502 }
1503
1504 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineInterpolatePoint_r(const void * p_cache,gaiaGeomCollPtr geom,double fraction)1505 gaiaLineInterpolatePoint_r (const void *p_cache, gaiaGeomCollPtr geom,
1506 double fraction)
1507 {
1508 /*
1509 * attempts to intepolate a point on line at dist "fraction"
1510 *
1511 * the fraction is expressed into the range from 0.0 to 1.0
1512 */
1513 int pts = 0;
1514 int lns = 0;
1515 int pgs = 0;
1516 gaiaGeomCollPtr result;
1517 gaiaPointPtr pt;
1518 gaiaLinestringPtr ln;
1519 gaiaPolygonPtr pg;
1520 GEOSGeometry *g;
1521 GEOSGeometry *g_pt;
1522 double length;
1523 double projection;
1524 struct splite_internal_cache *cache =
1525 (struct splite_internal_cache *) p_cache;
1526 GEOSContextHandle_t handle = NULL;
1527 if (cache == NULL)
1528 return NULL;
1529 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1530 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1531 return NULL;
1532 handle = cache->GEOS_handle;
1533 if (handle == NULL)
1534 return NULL;
1535 gaiaResetGeosMsg_r (cache);
1536 if (!geom)
1537 return NULL;
1538
1539 /* checking if a single Linestring has been passed */
1540 pt = geom->FirstPoint;
1541 while (pt)
1542 {
1543 pts++;
1544 pt = pt->Next;
1545 }
1546 ln = geom->FirstLinestring;
1547 while (ln)
1548 {
1549 lns++;
1550 ln = ln->Next;
1551 }
1552 pg = geom->FirstPolygon;
1553 while (pg)
1554 {
1555 pgs++;
1556 pg = pg->Next;
1557 }
1558 if (pts == 0 && lns == 1 && pgs == 0)
1559 ;
1560 else
1561 return NULL;
1562
1563 g = gaiaToGeos_r (cache, geom);
1564 if (GEOSLength_r (handle, g, &length))
1565 {
1566 /* transforming fraction to length */
1567 if (fraction < 0.0)
1568 fraction = 0.0;
1569 if (fraction > 1.0)
1570 fraction = 1.0;
1571 projection = length * fraction;
1572 }
1573 else
1574 {
1575 GEOSGeom_destroy_r (handle, g);
1576 return NULL;
1577 }
1578 g_pt = GEOSInterpolate_r (handle, g, projection);
1579 GEOSGeom_destroy_r (handle, g);
1580 if (!g_pt)
1581 return NULL;
1582 if (geom->DimensionModel == GAIA_XY_Z)
1583 result = gaiaFromGeos_XYZ_r (cache, g_pt);
1584 else if (geom->DimensionModel == GAIA_XY_M)
1585 result = gaiaFromGeos_XYM_r (cache, g_pt);
1586 else if (geom->DimensionModel == GAIA_XY_Z_M)
1587 result = gaiaFromGeos_XYZM_r (cache, g_pt);
1588 else
1589 result = gaiaFromGeos_XY_r (cache, g_pt);
1590 GEOSGeom_destroy_r (handle, g_pt);
1591 if (result == NULL)
1592 return NULL;
1593 result->Srid = geom->Srid;
1594 return result;
1595 }
1596
1597 static gaiaGeomCollPtr
gaiaLineInterpolateEquidistantPointsCommon(struct splite_internal_cache * cache,gaiaGeomCollPtr geom,double distance)1598 gaiaLineInterpolateEquidistantPointsCommon (struct splite_internal_cache *cache,
1599 gaiaGeomCollPtr geom,
1600 double distance)
1601 {
1602 /*
1603 * attempts to intepolate a set of points on line at regular distances
1604 */
1605 int pts = 0;
1606 int lns = 0;
1607 int pgs = 0;
1608 gaiaGeomCollPtr result;
1609 gaiaGeomCollPtr xpt;
1610 gaiaPointPtr pt;
1611 gaiaLinestringPtr ln;
1612 gaiaPolygonPtr pg;
1613 GEOSGeometry *g;
1614 GEOSGeometry *g_pt;
1615 double length;
1616 double current_length = 0.0;
1617 GEOSContextHandle_t handle = NULL;
1618
1619 if (cache != NULL)
1620 {
1621 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1622 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1623 return NULL;
1624 handle = cache->GEOS_handle;
1625 if (handle == NULL)
1626 return NULL;
1627 }
1628 #ifdef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1629 if (handle == NULL)
1630 return NULL;
1631 #endif
1632 if (!geom)
1633 return NULL;
1634 if (distance <= 0.0)
1635 return NULL;
1636
1637 /* checking if a single Linestring has been passed */
1638 pt = geom->FirstPoint;
1639 while (pt)
1640 {
1641 pts++;
1642 pt = pt->Next;
1643 }
1644 ln = geom->FirstLinestring;
1645 while (ln)
1646 {
1647 lns++;
1648 ln = ln->Next;
1649 }
1650 pg = geom->FirstPolygon;
1651 while (pg)
1652 {
1653 pgs++;
1654 pg = pg->Next;
1655 }
1656 if (pts == 0 && lns == 1 && pgs == 0)
1657 ;
1658 else
1659 return NULL;
1660
1661 if (cache != NULL)
1662 {
1663 g = gaiaToGeos_r (cache, geom);
1664 if (GEOSLength_r (handle, g, &length))
1665 {
1666 if (length <= distance)
1667 {
1668 /* the line is too short to apply interpolation */
1669 GEOSGeom_destroy_r (handle, g);
1670 return NULL;
1671 }
1672 }
1673 else
1674 {
1675 GEOSGeom_destroy_r (handle, g);
1676 return NULL;
1677 }
1678 }
1679 else
1680 {
1681 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1682 g = gaiaToGeos (geom);
1683 if (GEOSLength (g, &length))
1684 {
1685 if (length <= distance)
1686 {
1687 /* the line is too short to apply interpolation */
1688 GEOSGeom_destroy (g);
1689 return NULL;
1690 }
1691 }
1692 else
1693 {
1694 GEOSGeom_destroy (g);
1695 return NULL;
1696 }
1697 #endif
1698 }
1699
1700 /* creating the MultiPoint [always supporting M] */
1701 if (geom->DimensionModel == GAIA_XY_Z
1702 || geom->DimensionModel == GAIA_XY_Z_M)
1703 result = gaiaAllocGeomCollXYZM ();
1704 else
1705 result = gaiaAllocGeomCollXYM ();
1706 if (result == NULL)
1707 {
1708 if (cache != NULL)
1709 {
1710 GEOSGeom_destroy_r (handle, g);
1711 return NULL;
1712 }
1713 else
1714 {
1715 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1716 GEOSGeom_destroy (g);
1717 return NULL;
1718 #endif
1719 }
1720 }
1721
1722 while (1)
1723 {
1724 /* increasing the current distance */
1725 current_length += distance;
1726 if (current_length >= length)
1727 break;
1728 /* interpolating a point */
1729 if (handle != NULL)
1730 g_pt = GEOSInterpolate_r (handle, g, current_length);
1731 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1732 else
1733 g_pt = GEOSInterpolate (g, current_length);
1734 #endif
1735 if (!g_pt)
1736 goto error;
1737 if (geom->DimensionModel == GAIA_XY_Z)
1738 {
1739 if (cache != NULL)
1740 xpt = gaiaFromGeos_XYZ_r (cache, g_pt);
1741 else
1742 xpt = gaiaFromGeos_XYZ (g_pt);
1743 if (!xpt)
1744 goto error;
1745 pt = xpt->FirstPoint;
1746 if (!pt)
1747 goto error;
1748 gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z,
1749 current_length);
1750 }
1751 else if (geom->DimensionModel == GAIA_XY_M)
1752 {
1753 if (cache != NULL)
1754 xpt = gaiaFromGeos_XYM_r (cache, g_pt);
1755 else
1756 xpt = gaiaFromGeos_XYM (g_pt);
1757 if (!xpt)
1758 goto error;
1759 pt = xpt->FirstPoint;
1760 if (!pt)
1761 goto error;
1762 gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y,
1763 current_length);
1764 }
1765 else if (geom->DimensionModel == GAIA_XY_Z_M)
1766 {
1767 if (cache != NULL)
1768 xpt = gaiaFromGeos_XYZM_r (cache, g_pt);
1769 else
1770 xpt = gaiaFromGeos_XYZM (g_pt);
1771 if (!xpt)
1772 goto error;
1773 pt = xpt->FirstPoint;
1774 if (!pt)
1775 goto error;
1776 gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z,
1777 current_length);
1778 }
1779 else
1780 {
1781 if (cache != NULL)
1782 xpt = gaiaFromGeos_XY_r (cache, g_pt);
1783 else
1784 xpt = gaiaFromGeos_XY (g_pt);
1785 if (!xpt)
1786 goto error;
1787 pt = xpt->FirstPoint;
1788 if (!pt)
1789 goto error;
1790 gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y,
1791 current_length);
1792 }
1793 if (handle != NULL)
1794 GEOSGeom_destroy_r (handle, g_pt);
1795 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1796 else
1797 GEOSGeom_destroy (g_pt);
1798 #endif
1799 gaiaFreeGeomColl (xpt);
1800 }
1801 if (handle != NULL)
1802 GEOSGeom_destroy_r (handle, g);
1803 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1804 else
1805 GEOSGeom_destroy (g);
1806 #endif
1807 result->Srid = geom->Srid;
1808 result->DeclaredType = GAIA_MULTIPOINT;
1809 return result;
1810
1811 error:
1812 if (handle != NULL)
1813 {
1814 if (g_pt)
1815 GEOSGeom_destroy_r (handle, g_pt);
1816 GEOSGeom_destroy_r (handle, g);
1817 }
1818 else
1819 {
1820 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1821 if (g_pt)
1822 GEOSGeom_destroy (g_pt);
1823 GEOSGeom_destroy (g);
1824 #endif
1825 }
1826 gaiaFreeGeomColl (result);
1827 return NULL;
1828 }
1829
1830 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineInterpolateEquidistantPoints(gaiaGeomCollPtr geom,double distance)1831 gaiaLineInterpolateEquidistantPoints (gaiaGeomCollPtr geom, double distance)
1832 {
1833 gaiaResetGeosMsg ();
1834 return gaiaLineInterpolateEquidistantPointsCommon (NULL, geom, distance);
1835 }
1836
1837 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineInterpolateEquidistantPoints_r(const void * p_cache,gaiaGeomCollPtr geom,double distance)1838 gaiaLineInterpolateEquidistantPoints_r (const void *p_cache,
1839 gaiaGeomCollPtr geom, double distance)
1840 {
1841 struct splite_internal_cache *cache =
1842 (struct splite_internal_cache *) p_cache;
1843 GEOSContextHandle_t handle = NULL;
1844 if (cache == NULL)
1845 return NULL;
1846 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1847 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1848 return NULL;
1849 handle = cache->GEOS_handle;
1850 if (handle == NULL)
1851 return NULL;
1852 gaiaResetGeosMsg_r (cache);
1853 return gaiaLineInterpolateEquidistantPointsCommon (cache, geom, distance);
1854 }
1855
1856 GAIAGEO_DECLARE double
gaiaLineLocatePoint(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)1857 gaiaLineLocatePoint (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1858 {
1859 /*
1860 * attempts to compute the location of the closest point on LineString
1861 * to the given Point, as a fraction of total 2d line length
1862 *
1863 * the fraction is expressed into the range from 0.0 to 1.0
1864 */
1865 double result = -1.0;
1866 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
1867 int pts1 = 0;
1868 int lns1 = 0;
1869 int pgs1 = 0;
1870 int pts2 = 0;
1871 int lns2 = 0;
1872 int pgs2 = 0;
1873 double length;
1874 double projection;
1875 gaiaPointPtr pt;
1876 gaiaLinestringPtr ln;
1877 gaiaPolygonPtr pg;
1878 GEOSGeometry *g1;
1879 GEOSGeometry *g2;
1880 gaiaResetGeosMsg ();
1881 if (!geom1 || !geom2)
1882 return -1.0;
1883
1884 /* checking if a single Linestring has been passed */
1885 pt = geom1->FirstPoint;
1886 while (pt)
1887 {
1888 pts1++;
1889 pt = pt->Next;
1890 }
1891 ln = geom1->FirstLinestring;
1892 while (ln)
1893 {
1894 lns1++;
1895 ln = ln->Next;
1896 }
1897 pg = geom1->FirstPolygon;
1898 while (pg)
1899 {
1900 pgs1++;
1901 pg = pg->Next;
1902 }
1903 if (pts1 == 0 && lns1 >= 1 && pgs1 == 0)
1904 ;
1905 else
1906 return -1.0;
1907
1908 /* checking if a single Point has been passed */
1909 pt = geom2->FirstPoint;
1910 while (pt)
1911 {
1912 pts2++;
1913 pt = pt->Next;
1914 }
1915 ln = geom2->FirstLinestring;
1916 while (ln)
1917 {
1918 lns2++;
1919 ln = ln->Next;
1920 }
1921 pg = geom2->FirstPolygon;
1922 while (pg)
1923 {
1924 pgs2++;
1925 pg = pg->Next;
1926 }
1927 if (pts2 == 1 && lns2 == 0 && pgs2 == 0)
1928 ;
1929 else
1930 return -1.0;
1931
1932 g1 = gaiaToGeos (geom1);
1933 g2 = gaiaToGeos (geom2);
1934 projection = GEOSProject (g1, g2);
1935 if (GEOSLength (g1, &length))
1936 {
1937 /* normalizing as a fraction between 0.0 and 1.0 */
1938 result = projection / length;
1939 }
1940 else
1941 result = -1.0;
1942 GEOSGeom_destroy (g1);
1943 GEOSGeom_destroy (g2);
1944 #else
1945 if (geom1 == NULL || geom2 == NULL)
1946 geom1 = NULL; /* silencing stupid compiler warnings */
1947 #endif
1948 return result;
1949 }
1950
1951 GAIAGEO_DECLARE double
gaiaLineLocatePoint_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)1952 gaiaLineLocatePoint_r (const void *p_cache, gaiaGeomCollPtr geom1,
1953 gaiaGeomCollPtr geom2)
1954 {
1955 /*
1956 * attempts to compute the location of the closest point on LineString
1957 * to the given Point, as a fraction of total 2d line length
1958 *
1959 * the fraction is expressed into the range from 0.0 to 1.0
1960 */
1961 int pts1 = 0;
1962 int lns1 = 0;
1963 int pgs1 = 0;
1964 int pts2 = 0;
1965 int lns2 = 0;
1966 int pgs2 = 0;
1967 double length;
1968 double projection;
1969 double result;
1970 gaiaPointPtr pt;
1971 gaiaLinestringPtr ln;
1972 gaiaPolygonPtr pg;
1973 GEOSGeometry *g1;
1974 GEOSGeometry *g2;
1975 struct splite_internal_cache *cache =
1976 (struct splite_internal_cache *) p_cache;
1977 GEOSContextHandle_t handle = NULL;
1978 if (cache == NULL)
1979 return -1.0;
1980 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1981 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1982 return -1.0;
1983 handle = cache->GEOS_handle;
1984 if (handle == NULL)
1985 return -1.0;
1986 gaiaResetGeosMsg_r (cache);
1987 if (!geom1 || !geom2)
1988 return -1.0;
1989
1990 /* checking if a single Linestring has been passed */
1991 pt = geom1->FirstPoint;
1992 while (pt)
1993 {
1994 pts1++;
1995 pt = pt->Next;
1996 }
1997 ln = geom1->FirstLinestring;
1998 while (ln)
1999 {
2000 lns1++;
2001 ln = ln->Next;
2002 }
2003 pg = geom1->FirstPolygon;
2004 while (pg)
2005 {
2006 pgs1++;
2007 pg = pg->Next;
2008 }
2009 if (pts1 == 0 && lns1 >= 1 && pgs1 == 0)
2010 ;
2011 else
2012 return -1.0;
2013
2014 /* checking if a single Point has been passed */
2015 pt = geom2->FirstPoint;
2016 while (pt)
2017 {
2018 pts2++;
2019 pt = pt->Next;
2020 }
2021 ln = geom2->FirstLinestring;
2022 while (ln)
2023 {
2024 lns2++;
2025 ln = ln->Next;
2026 }
2027 pg = geom2->FirstPolygon;
2028 while (pg)
2029 {
2030 pgs2++;
2031 pg = pg->Next;
2032 }
2033 if (pts2 == 1 && lns2 == 0 && pgs2 == 0)
2034 ;
2035 else
2036 return -1.0;
2037
2038 g1 = gaiaToGeos_r (cache, geom1);
2039 g2 = gaiaToGeos_r (cache, geom2);
2040 projection = GEOSProject_r (handle, g1, g2);
2041 if (GEOSLength_r (handle, g1, &length))
2042 {
2043 /* normalizing as a fraction between 0.0 and 1.0 */
2044 result = projection / length;
2045 }
2046 else
2047 result = -1.0;
2048 GEOSGeom_destroy_r (handle, g1);
2049 GEOSGeom_destroy_r (handle, g2);
2050 return result;
2051 }
2052
2053 static gaiaGeomCollPtr
gaiaLineSubstringCommon(struct splite_internal_cache * cache,gaiaGeomCollPtr geom,double start_fraction,double end_fraction)2054 gaiaLineSubstringCommon (struct splite_internal_cache *cache,
2055 gaiaGeomCollPtr geom, double start_fraction,
2056 double end_fraction)
2057 {
2058 /*
2059 * attempts to build a new Linestring being a substring of the input one starting
2060 * and ending at the given fractions of total 2d length
2061 */
2062 int pts = 0;
2063 int lns = 0;
2064 int pgs = 0;
2065 gaiaGeomCollPtr result;
2066 gaiaPointPtr pt;
2067 gaiaLinestringPtr ln;
2068 gaiaLinestringPtr out;
2069 gaiaPolygonPtr pg;
2070 GEOSGeometry *g;
2071 GEOSGeometry *g_start;
2072 GEOSGeometry *g_end;
2073 GEOSCoordSequence *cs;
2074 const GEOSCoordSequence *in_cs;
2075 GEOSGeometry *segm;
2076 double length;
2077 double total = 0.0;
2078 double start;
2079 double end;
2080 int iv;
2081 int i_start = -1;
2082 int i_end = -1;
2083 int points;
2084 double x;
2085 double y;
2086 double z;
2087 double m;
2088 double x0;
2089 double y0;
2090 unsigned int dims;
2091 GEOSContextHandle_t handle = NULL;
2092
2093 if (cache != NULL)
2094 {
2095 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
2096 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
2097 return NULL;
2098 handle = cache->GEOS_handle;
2099 if (handle == NULL)
2100 return NULL;
2101 }
2102 #ifdef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2103 if (handle == NULL)
2104 return NULL;
2105 #endif
2106 if (!geom)
2107 return NULL;
2108
2109 /* checking if a single Linestring has been passed */
2110 pt = geom->FirstPoint;
2111 while (pt)
2112 {
2113 pts++;
2114 pt = pt->Next;
2115 }
2116 ln = geom->FirstLinestring;
2117 while (ln)
2118 {
2119 lns++;
2120 ln = ln->Next;
2121 }
2122 pg = geom->FirstPolygon;
2123 while (pg)
2124 {
2125 pgs++;
2126 pg = pg->Next;
2127 }
2128 if (pts == 0 && lns == 1 && pgs == 0)
2129 ;
2130 else
2131 return NULL;
2132
2133 if (start_fraction < 0.0)
2134 start_fraction = 0.0;
2135 if (start_fraction > 1.0)
2136 start_fraction = 1.0;
2137 if (end_fraction < 0.0)
2138 end_fraction = 0.0;
2139 if (end_fraction > 1.0)
2140 end_fraction = 1.0;
2141 if (start_fraction >= end_fraction)
2142 return NULL;
2143 if (cache != NULL)
2144 {
2145 g = gaiaToGeos_r (cache, geom);
2146 if (GEOSLength_r (handle, g, &length))
2147 {
2148 start = length * start_fraction;
2149 end = length * end_fraction;
2150 }
2151 else
2152 {
2153 GEOSGeom_destroy_r (handle, g);
2154 return NULL;
2155 }
2156 g_start = GEOSInterpolate_r (handle, g, start);
2157 g_end = GEOSInterpolate_r (handle, g, end);
2158 GEOSGeom_destroy_r (handle, g);
2159 }
2160 else
2161 {
2162 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2163 g = gaiaToGeos (geom);
2164 if (GEOSLength (g, &length))
2165 {
2166 start = length * start_fraction;
2167 end = length * end_fraction;
2168 }
2169 else
2170 {
2171 GEOSGeom_destroy (g);
2172 return NULL;
2173 }
2174 g_start = GEOSInterpolate (g, start);
2175 g_end = GEOSInterpolate (g, end);
2176 GEOSGeom_destroy (g);
2177 #endif
2178 }
2179 if (!g_start || !g_end)
2180 return NULL;
2181
2182 /* identifying first and last valid vertex */
2183 x0 = 0.0;
2184 y0 = 0.0;
2185 ln = geom->FirstLinestring;
2186 for (iv = 0; iv < ln->Points; iv++)
2187 {
2188 switch (ln->DimensionModel)
2189 {
2190 case GAIA_XY_Z:
2191 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2192 break;
2193 case GAIA_XY_M:
2194 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2195 break;
2196 case GAIA_XY_Z_M:
2197 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2198 break;
2199 default:
2200 gaiaGetPoint (ln->Coords, iv, &x, &y);
2201 break;
2202 };
2203
2204 if (iv > 0)
2205 {
2206 if (handle != NULL)
2207 {
2208 cs = GEOSCoordSeq_create_r (handle, 2, 2);
2209 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
2210 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
2211 GEOSCoordSeq_setX_r (handle, cs, 1, x);
2212 GEOSCoordSeq_setY_r (handle, cs, 1, y);
2213 segm = GEOSGeom_createLineString_r (handle, cs);
2214 GEOSLength_r (handle, segm, &length);
2215 total += length;
2216 GEOSGeom_destroy_r (handle, segm);
2217 }
2218 else
2219 {
2220 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2221 cs = GEOSCoordSeq_create (2, 2);
2222 GEOSCoordSeq_setX (cs, 0, x0);
2223 GEOSCoordSeq_setY (cs, 0, y0);
2224 GEOSCoordSeq_setX (cs, 1, x);
2225 GEOSCoordSeq_setY (cs, 1, y);
2226 segm = GEOSGeom_createLineString (cs);
2227 GEOSLength (segm, &length);
2228 total += length;
2229 GEOSGeom_destroy (segm);
2230 #endif
2231 }
2232 if (total > start && i_start < 0)
2233 i_start = iv;
2234 if (total < end)
2235 i_end = iv;
2236 }
2237 x0 = x;
2238 y0 = y;
2239 }
2240 if (i_start < 0 || i_end < 0)
2241 {
2242 i_start = -1;
2243 i_end = -1;
2244 points = 2;
2245 }
2246 else
2247 points = i_end - i_start + 3;
2248
2249 /* creating the output geometry */
2250 switch (ln->DimensionModel)
2251 {
2252 case GAIA_XY_Z:
2253 result = gaiaAllocGeomCollXYZ ();
2254 break;
2255 case GAIA_XY_M:
2256 result = gaiaAllocGeomCollXYM ();
2257 break;
2258 case GAIA_XY_Z_M:
2259 result = gaiaAllocGeomCollXYZM ();
2260 break;
2261 default:
2262 result = gaiaAllocGeomColl ();
2263 break;
2264 };
2265 result->Srid = geom->Srid;
2266 out = gaiaAddLinestringToGeomColl (result, points);
2267
2268 /* start vertex */
2269 points = 0;
2270 if (handle)
2271 {
2272 in_cs = GEOSGeom_getCoordSeq_r (handle, g_start);
2273 GEOSCoordSeq_getDimensions_r (handle, in_cs, &dims);
2274 if (dims == 3)
2275 {
2276 GEOSCoordSeq_getX_r (handle, in_cs, 0, &x);
2277 GEOSCoordSeq_getY_r (handle, in_cs, 0, &y);
2278 GEOSCoordSeq_getZ_r (handle, in_cs, 0, &z);
2279 m = 0.0;
2280 }
2281 else
2282 {
2283 GEOSCoordSeq_getX_r (handle, in_cs, 0, &x);
2284 GEOSCoordSeq_getY_r (handle, in_cs, 0, &y);
2285 z = 0.0;
2286 m = 0.0;
2287 }
2288 GEOSGeom_destroy_r (handle, g_start);
2289 }
2290 else
2291 {
2292 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2293 in_cs = GEOSGeom_getCoordSeq (g_start);
2294 GEOSCoordSeq_getDimensions (in_cs, &dims);
2295 if (dims == 3)
2296 {
2297 GEOSCoordSeq_getX (in_cs, 0, &x);
2298 GEOSCoordSeq_getY (in_cs, 0, &y);
2299 GEOSCoordSeq_getZ (in_cs, 0, &z);
2300 m = 0.0;
2301 }
2302 else
2303 {
2304 GEOSCoordSeq_getX (in_cs, 0, &x);
2305 GEOSCoordSeq_getY (in_cs, 0, &y);
2306 z = 0.0;
2307 m = 0.0;
2308 }
2309 GEOSGeom_destroy (g_start);
2310 #endif
2311 }
2312 switch (out->DimensionModel)
2313 {
2314 case GAIA_XY_Z:
2315 gaiaSetPointXYZ (out->Coords, points, x, y, z);
2316 break;
2317 case GAIA_XY_M:
2318 gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
2319 break;
2320 case GAIA_XY_Z_M:
2321 gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
2322 break;
2323 default:
2324 gaiaSetPoint (out->Coords, points, x, y);
2325 break;
2326 };
2327 points++;
2328
2329 if (i_start < 0 || i_end < 0)
2330 ;
2331 else
2332 {
2333 for (iv = i_start; iv <= i_end; iv++)
2334 {
2335 z = 0.0;
2336 m = 0.0;
2337 switch (ln->DimensionModel)
2338 {
2339 case GAIA_XY_Z:
2340 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2341 break;
2342 case GAIA_XY_M:
2343 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2344 break;
2345 case GAIA_XY_Z_M:
2346 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2347 break;
2348 default:
2349 gaiaGetPoint (ln->Coords, iv, &x, &y);
2350 break;
2351 };
2352 switch (out->DimensionModel)
2353 {
2354 case GAIA_XY_Z:
2355 gaiaSetPointXYZ (out->Coords, points, x, y, z);
2356 break;
2357 case GAIA_XY_M:
2358 gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
2359 break;
2360 case GAIA_XY_Z_M:
2361 gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
2362 break;
2363 default:
2364 gaiaSetPoint (out->Coords, points, x, y);
2365 break;
2366 };
2367 points++;
2368 }
2369 }
2370
2371 /* end vertex */
2372 if (handle != NULL)
2373 {
2374 in_cs = GEOSGeom_getCoordSeq_r (handle, g_end);
2375 GEOSCoordSeq_getDimensions_r (handle, in_cs, &dims);
2376 if (dims == 3)
2377 {
2378 GEOSCoordSeq_getX_r (handle, in_cs, 0, &x);
2379 GEOSCoordSeq_getY_r (handle, in_cs, 0, &y);
2380 GEOSCoordSeq_getZ_r (handle, in_cs, 0, &z);
2381 m = 0.0;
2382 }
2383 else
2384 {
2385 GEOSCoordSeq_getX_r (handle, in_cs, 0, &x);
2386 GEOSCoordSeq_getY_r (handle, in_cs, 0, &y);
2387 z = 0.0;
2388 m = 0.0;
2389 }
2390 GEOSGeom_destroy_r (handle, g_end);
2391 }
2392 else
2393 {
2394 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2395 in_cs = GEOSGeom_getCoordSeq (g_end);
2396 GEOSCoordSeq_getDimensions (in_cs, &dims);
2397 if (dims == 3)
2398 {
2399 GEOSCoordSeq_getX (in_cs, 0, &x);
2400 GEOSCoordSeq_getY (in_cs, 0, &y);
2401 GEOSCoordSeq_getZ (in_cs, 0, &z);
2402 m = 0.0;
2403 }
2404 else
2405 {
2406 GEOSCoordSeq_getX (in_cs, 0, &x);
2407 GEOSCoordSeq_getY (in_cs, 0, &y);
2408 z = 0.0;
2409 m = 0.0;
2410 }
2411 GEOSGeom_destroy (g_end);
2412 #endif
2413 }
2414 switch (out->DimensionModel)
2415 {
2416 case GAIA_XY_Z:
2417 gaiaSetPointXYZ (out->Coords, points, x, y, z);
2418 break;
2419 case GAIA_XY_M:
2420 gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
2421 break;
2422 case GAIA_XY_Z_M:
2423 gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
2424 break;
2425 default:
2426 gaiaSetPoint (out->Coords, points, x, y);
2427 break;
2428 };
2429 return result;
2430 }
2431
2432 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineSubstring(gaiaGeomCollPtr geom,double start_fraction,double end_fraction)2433 gaiaLineSubstring (gaiaGeomCollPtr geom, double start_fraction,
2434 double end_fraction)
2435 {
2436 gaiaResetGeosMsg ();
2437 return gaiaLineSubstringCommon (NULL, geom, start_fraction, end_fraction);
2438 }
2439
2440 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineSubstring_r(const void * p_cache,gaiaGeomCollPtr geom,double start_fraction,double end_fraction)2441 gaiaLineSubstring_r (const void *p_cache, gaiaGeomCollPtr geom,
2442 double start_fraction, double end_fraction)
2443 {
2444 struct splite_internal_cache *cache =
2445 (struct splite_internal_cache *) p_cache;
2446 GEOSContextHandle_t handle = NULL;
2447 if (cache == NULL)
2448 return NULL;
2449 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
2450 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
2451 return NULL;
2452 handle = cache->GEOS_handle;
2453 if (handle == NULL)
2454 return NULL;
2455 gaiaResetGeosMsg_r (cache);
2456 return gaiaLineSubstringCommon (cache, geom, start_fraction, end_fraction);
2457 }
2458
2459 static GEOSGeometry *
buildGeosPoints(GEOSContextHandle_t handle,const gaiaGeomCollPtr gaia)2460 buildGeosPoints (GEOSContextHandle_t handle, const gaiaGeomCollPtr gaia)
2461 {
2462 /* converting a GAIA Geometry into a GEOS Geometry of POINTS */
2463 int pts = 0;
2464 unsigned int dims;
2465 int iv;
2466 int ib;
2467 int nItem;
2468 double x;
2469 double y;
2470 double z;
2471 double m;
2472 gaiaPointPtr pt;
2473 gaiaLinestringPtr ln;
2474 gaiaPolygonPtr pg;
2475 gaiaRingPtr rng;
2476 GEOSGeometry *geos = NULL;
2477 GEOSGeometry *geos_item;
2478 GEOSGeometry **geos_coll;
2479 GEOSCoordSequence *cs;
2480
2481 #ifdef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2482 if (handle == NULL)
2483 return NULL;
2484 #endif
2485 if (!gaia)
2486 return NULL;
2487
2488 pt = gaia->FirstPoint;
2489 while (pt)
2490 {
2491 /* counting how many POINTs are there */
2492 pts++;
2493 pt = pt->Next;
2494 }
2495 ln = gaia->FirstLinestring;
2496 while (ln)
2497 {
2498 /* counting how many POINTs are there */
2499 pts += ln->Points;
2500 ln = ln->Next;
2501 }
2502 pg = gaia->FirstPolygon;
2503 while (pg)
2504 {
2505 /* counting how many POINTs are there */
2506 rng = pg->Exterior;
2507 pts += rng->Points - 1; /* exterior ring */
2508 for (ib = 0; ib < pg->NumInteriors; ib++)
2509 {
2510 /* interior ring */
2511 rng = pg->Interiors + ib;
2512 pts += rng->Points - 1;
2513 }
2514 pg = pg->Next;
2515 }
2516 if (pts == 0)
2517 return NULL;
2518 switch (gaia->DimensionModel)
2519 {
2520 case GAIA_XY_Z:
2521 case GAIA_XY_Z_M:
2522 dims = 3;
2523 break;
2524 default:
2525 dims = 2;
2526 break;
2527 };
2528 nItem = 0;
2529 geos_coll = malloc (sizeof (GEOSGeometry *) * (pts));
2530 pt = gaia->FirstPoint;
2531 while (pt)
2532 {
2533 if (handle != NULL)
2534 {
2535 cs = GEOSCoordSeq_create_r (handle, 1, dims);
2536 switch (pt->DimensionModel)
2537 {
2538 case GAIA_XY_Z:
2539 case GAIA_XY_Z_M:
2540 GEOSCoordSeq_setX_r (handle, cs, 0, pt->X);
2541 GEOSCoordSeq_setY_r (handle, cs, 0, pt->Y);
2542 GEOSCoordSeq_setZ_r (handle, cs, 0, pt->Z);
2543 break;
2544 default:
2545 GEOSCoordSeq_setX_r (handle, cs, 0, pt->X);
2546 GEOSCoordSeq_setY_r (handle, cs, 0, pt->Y);
2547 break;
2548 };
2549 geos_item = GEOSGeom_createPoint_r (handle, cs);
2550 }
2551 else
2552 {
2553 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2554 cs = GEOSCoordSeq_create (1, dims);
2555 switch (pt->DimensionModel)
2556 {
2557 case GAIA_XY_Z:
2558 case GAIA_XY_Z_M:
2559 GEOSCoordSeq_setX (cs, 0, pt->X);
2560 GEOSCoordSeq_setY (cs, 0, pt->Y);
2561 GEOSCoordSeq_setZ (cs, 0, pt->Z);
2562 break;
2563 default:
2564 GEOSCoordSeq_setX (cs, 0, pt->X);
2565 GEOSCoordSeq_setY (cs, 0, pt->Y);
2566 break;
2567 };
2568 geos_item = GEOSGeom_createPoint (cs);
2569 #endif
2570 }
2571 *(geos_coll + nItem++) = geos_item;
2572 pt = pt->Next;
2573 }
2574 ln = gaia->FirstLinestring;
2575 while (ln)
2576 {
2577 for (iv = 0; iv < ln->Points; iv++)
2578 {
2579 z = 0.0;
2580 m = 0.0;
2581 switch (ln->DimensionModel)
2582 {
2583 case GAIA_XY_Z:
2584 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2585 break;
2586 case GAIA_XY_M:
2587 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2588 break;
2589 case GAIA_XY_Z_M:
2590 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2591 break;
2592 default:
2593 gaiaGetPoint (ln->Coords, iv, &x, &y);
2594 break;
2595 };
2596 if (handle != NULL)
2597 {
2598 cs = GEOSCoordSeq_create_r (handle, 1, dims);
2599 if (dims == 3)
2600 {
2601 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2602 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2603 GEOSCoordSeq_setZ_r (handle, cs, 0, z);
2604 }
2605 else
2606 {
2607 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2608 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2609 }
2610 geos_item = GEOSGeom_createPoint_r (handle, cs);
2611 }
2612 else
2613 {
2614 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2615 cs = GEOSCoordSeq_create (1, dims);
2616 if (dims == 3)
2617 {
2618 GEOSCoordSeq_setX (cs, 0, x);
2619 GEOSCoordSeq_setY (cs, 0, y);
2620 GEOSCoordSeq_setZ (cs, 0, z);
2621 }
2622 else
2623 {
2624 GEOSCoordSeq_setX (cs, 0, x);
2625 GEOSCoordSeq_setY (cs, 0, y);
2626 }
2627 geos_item = GEOSGeom_createPoint (cs);
2628 #endif
2629 }
2630 *(geos_coll + nItem++) = geos_item;
2631 }
2632 ln = ln->Next;
2633 }
2634 pg = gaia->FirstPolygon;
2635 while (pg)
2636 {
2637 rng = pg->Exterior;
2638 for (iv = 1; iv < rng->Points; iv++)
2639 {
2640 /* exterior ring */
2641 z = 0.0;
2642 m = 0.0;
2643 switch (rng->DimensionModel)
2644 {
2645 case GAIA_XY_Z:
2646 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2647 break;
2648 case GAIA_XY_M:
2649 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2650 break;
2651 case GAIA_XY_Z_M:
2652 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2653 break;
2654 default:
2655 gaiaGetPoint (rng->Coords, iv, &x, &y);
2656 break;
2657 };
2658 if (handle != NULL)
2659 {
2660 cs = GEOSCoordSeq_create_r (handle, 1, dims);
2661 if (dims == 3)
2662 {
2663 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2664 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2665 GEOSCoordSeq_setZ_r (handle, cs, 0, z);
2666 }
2667 else
2668 {
2669 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2670 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2671 }
2672 geos_item = GEOSGeom_createPoint_r (handle, cs);
2673 }
2674 else
2675 {
2676 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2677 cs = GEOSCoordSeq_create (1, dims);
2678 if (dims == 3)
2679 {
2680 GEOSCoordSeq_setX (cs, 0, x);
2681 GEOSCoordSeq_setY (cs, 0, y);
2682 GEOSCoordSeq_setZ (cs, 0, z);
2683 }
2684 else
2685 {
2686 GEOSCoordSeq_setX (cs, 0, x);
2687 GEOSCoordSeq_setY (cs, 0, y);
2688 }
2689 geos_item = GEOSGeom_createPoint (cs);
2690 #endif
2691 }
2692 *(geos_coll + nItem++) = geos_item;
2693 }
2694 for (ib = 0; ib < pg->NumInteriors; ib++)
2695 {
2696 /* interior ring */
2697 rng = pg->Interiors + ib;
2698 for (iv = 1; iv < rng->Points; iv++)
2699 {
2700 /* exterior ring */
2701 z = 0.0;
2702 m = 0.0;
2703 switch (rng->DimensionModel)
2704 {
2705 case GAIA_XY_Z:
2706 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2707 break;
2708 case GAIA_XY_M:
2709 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2710 break;
2711 case GAIA_XY_Z_M:
2712 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2713 break;
2714 default:
2715 gaiaGetPoint (rng->Coords, iv, &x, &y);
2716 break;
2717 };
2718 if (handle != NULL)
2719 {
2720 cs = GEOSCoordSeq_create_r (handle, 1, dims);
2721 if (dims == 3)
2722 {
2723 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2724 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2725 GEOSCoordSeq_setZ_r (handle, cs, 0, z);
2726 }
2727 else
2728 {
2729 GEOSCoordSeq_setX_r (handle, cs, 0, x);
2730 GEOSCoordSeq_setY_r (handle, cs, 0, y);
2731 }
2732 geos_item = GEOSGeom_createPoint_r (handle, cs);
2733 }
2734 else
2735 {
2736 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2737 cs = GEOSCoordSeq_create (1, dims);
2738 if (dims == 3)
2739 {
2740 GEOSCoordSeq_setX (cs, 0, x);
2741 GEOSCoordSeq_setY (cs, 0, y);
2742 GEOSCoordSeq_setZ (cs, 0, z);
2743 }
2744 else
2745 {
2746 GEOSCoordSeq_setX (cs, 0, x);
2747 GEOSCoordSeq_setY (cs, 0, y);
2748 }
2749 geos_item = GEOSGeom_createPoint (cs);
2750 #endif
2751 }
2752 *(geos_coll + nItem++) = geos_item;
2753 }
2754 }
2755 pg = pg->Next;
2756 }
2757 if (handle != NULL)
2758 {
2759 geos =
2760 GEOSGeom_createCollection_r (handle, GEOS_MULTIPOINT, geos_coll,
2761 pts);
2762 free (geos_coll);
2763 GEOSSetSRID_r (handle, geos, gaia->Srid);
2764 }
2765 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2766 else
2767 {
2768 geos = GEOSGeom_createCollection (GEOS_MULTIPOINT, geos_coll, pts);
2769 free (geos_coll);
2770 GEOSSetSRID (geos, gaia->Srid);
2771 }
2772 #endif
2773 return geos;
2774 }
2775
2776 static GEOSGeometry *
buildGeosSegments(GEOSContextHandle_t handle,const gaiaGeomCollPtr gaia)2777 buildGeosSegments (GEOSContextHandle_t handle, const gaiaGeomCollPtr gaia)
2778 {
2779 /* converting a GAIA Geometry into a GEOS Geometry of SEGMENTS */
2780 int segms = 0;
2781 unsigned int dims;
2782 int iv;
2783 int ib;
2784 int nItem;
2785 double x;
2786 double y;
2787 double z;
2788 double m;
2789 double x0 = 0.0;
2790 double y0 = 0.0;
2791 double z0 = 0.0;
2792 gaiaLinestringPtr ln;
2793 gaiaPolygonPtr pg;
2794 gaiaRingPtr rng;
2795 GEOSGeometry *geos = NULL;
2796 GEOSGeometry *geos_item;
2797 GEOSGeometry **geos_coll;
2798 GEOSCoordSequence *cs;
2799
2800 #ifdef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2801 if (handle == NULL)
2802 return NULL;
2803 #endif
2804 if (!gaia)
2805 return NULL;
2806
2807 ln = gaia->FirstLinestring;
2808 while (ln)
2809 {
2810 /* counting how many SEGMENTs are there */
2811 segms += ln->Points - 1;
2812 ln = ln->Next;
2813 }
2814 pg = gaia->FirstPolygon;
2815 while (pg)
2816 {
2817 /* counting how many SEGMENTs are there */
2818 rng = pg->Exterior;
2819 segms += rng->Points - 1; /* exterior ring */
2820 for (ib = 0; ib < pg->NumInteriors; ib++)
2821 {
2822 /* interior ring */
2823 rng = pg->Interiors + ib;
2824 segms += rng->Points - 1;
2825 }
2826 pg = pg->Next;
2827 }
2828 if (segms == 0)
2829 return NULL;
2830 switch (gaia->DimensionModel)
2831 {
2832 case GAIA_XY_Z:
2833 case GAIA_XY_Z_M:
2834 dims = 3;
2835 break;
2836 default:
2837 dims = 2;
2838 break;
2839 };
2840 nItem = 0;
2841 geos_coll = malloc (sizeof (GEOSGeometry *) * (segms));
2842 ln = gaia->FirstLinestring;
2843 while (ln)
2844 {
2845 for (iv = 0; iv < ln->Points; iv++)
2846 {
2847 z = 0.0;
2848 m = 0.0;
2849 switch (ln->DimensionModel)
2850 {
2851 case GAIA_XY_Z:
2852 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2853 break;
2854 case GAIA_XY_M:
2855 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2856 break;
2857 case GAIA_XY_Z_M:
2858 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2859 break;
2860 default:
2861 gaiaGetPoint (ln->Coords, iv, &x, &y);
2862 break;
2863 };
2864 if (iv > 0)
2865 {
2866 if (handle != NULL)
2867 {
2868 cs = GEOSCoordSeq_create_r (handle, 2, dims);
2869 if (dims == 3)
2870 {
2871 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
2872 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
2873 GEOSCoordSeq_setZ_r (handle, cs, 0, z0);
2874 GEOSCoordSeq_setX_r (handle, cs, 1, x);
2875 GEOSCoordSeq_setY_r (handle, cs, 1, y);
2876 GEOSCoordSeq_setZ_r (handle, cs, 1, z);
2877 }
2878 else
2879 {
2880 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
2881 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
2882 GEOSCoordSeq_setX_r (handle, cs, 1, x);
2883 GEOSCoordSeq_setY_r (handle, cs, 1, y);
2884 }
2885 geos_item =
2886 GEOSGeom_createLineString_r (handle, cs);
2887 }
2888 else
2889 {
2890 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2891 cs = GEOSCoordSeq_create (2, dims);
2892 if (dims == 3)
2893 {
2894 GEOSCoordSeq_setX (cs, 0, x0);
2895 GEOSCoordSeq_setY (cs, 0, y0);
2896 GEOSCoordSeq_setZ (cs, 0, z0);
2897 GEOSCoordSeq_setX (cs, 1, x);
2898 GEOSCoordSeq_setY (cs, 1, y);
2899 GEOSCoordSeq_setZ (cs, 1, z);
2900 }
2901 else
2902 {
2903 GEOSCoordSeq_setX (cs, 0, x0);
2904 GEOSCoordSeq_setY (cs, 0, y0);
2905 GEOSCoordSeq_setX (cs, 1, x);
2906 GEOSCoordSeq_setY (cs, 1, y);
2907 }
2908 geos_item = GEOSGeom_createLineString (cs);
2909 #endif
2910 }
2911 *(geos_coll + nItem++) = geos_item;
2912 }
2913 x0 = x;
2914 y0 = y;
2915 z0 = z;
2916 }
2917 ln = ln->Next;
2918 }
2919 pg = gaia->FirstPolygon;
2920 while (pg)
2921 {
2922 rng = pg->Exterior;
2923 for (iv = 0; iv < rng->Points; iv++)
2924 {
2925 /* exterior ring */
2926 z = 0.0;
2927 m = 0.0;
2928 switch (rng->DimensionModel)
2929 {
2930 case GAIA_XY_Z:
2931 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2932 break;
2933 case GAIA_XY_M:
2934 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2935 break;
2936 case GAIA_XY_Z_M:
2937 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2938 break;
2939 default:
2940 gaiaGetPoint (rng->Coords, iv, &x, &y);
2941 break;
2942 };
2943 if (iv > 0)
2944 {
2945 if (handle != NULL)
2946 {
2947 cs = GEOSCoordSeq_create_r (handle, 2, dims);
2948 if (dims == 3)
2949 {
2950 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
2951 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
2952 GEOSCoordSeq_setZ_r (handle, cs, 0, z0);
2953 GEOSCoordSeq_setX_r (handle, cs, 1, x);
2954 GEOSCoordSeq_setY_r (handle, cs, 1, y);
2955 GEOSCoordSeq_setZ_r (handle, cs, 1, z);
2956 }
2957 else
2958 {
2959 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
2960 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
2961 GEOSCoordSeq_setX_r (handle, cs, 1, x);
2962 GEOSCoordSeq_setY_r (handle, cs, 1, y);
2963 }
2964 geos_item =
2965 GEOSGeom_createLineString_r (handle, cs);
2966 }
2967 else
2968 {
2969 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
2970 cs = GEOSCoordSeq_create (2, dims);
2971 if (dims == 3)
2972 {
2973 GEOSCoordSeq_setX (cs, 0, x0);
2974 GEOSCoordSeq_setY (cs, 0, y0);
2975 GEOSCoordSeq_setZ (cs, 0, z0);
2976 GEOSCoordSeq_setX (cs, 1, x);
2977 GEOSCoordSeq_setY (cs, 1, y);
2978 GEOSCoordSeq_setZ (cs, 1, z);
2979 }
2980 else
2981 {
2982 GEOSCoordSeq_setX (cs, 0, x0);
2983 GEOSCoordSeq_setY (cs, 0, y0);
2984 GEOSCoordSeq_setX (cs, 1, x);
2985 GEOSCoordSeq_setY (cs, 1, y);
2986 }
2987 geos_item = GEOSGeom_createLineString (cs);
2988 #endif
2989 }
2990 *(geos_coll + nItem++) = geos_item;
2991 }
2992 x0 = x;
2993 y0 = y;
2994 z0 = z;
2995 }
2996 for (ib = 0; ib < pg->NumInteriors; ib++)
2997 {
2998 /* interior ring */
2999 rng = pg->Interiors + ib;
3000 for (iv = 0; iv < rng->Points; iv++)
3001 {
3002 /* exterior ring */
3003 z = 0.0;
3004 m = 0.0;
3005 switch (rng->DimensionModel)
3006 {
3007 case GAIA_XY_Z:
3008 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3009 break;
3010 case GAIA_XY_M:
3011 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3012 break;
3013 case GAIA_XY_Z_M:
3014 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3015 break;
3016 default:
3017 gaiaGetPoint (rng->Coords, iv, &x, &y);
3018 break;
3019 };
3020 if (iv > 0)
3021 {
3022 if (handle != NULL)
3023 {
3024 cs = GEOSCoordSeq_create_r (handle, 2, dims);
3025 if (dims == 3)
3026 {
3027 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
3028 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
3029 GEOSCoordSeq_setZ_r (handle, cs, 0, z0);
3030 GEOSCoordSeq_setX_r (handle, cs, 1, x);
3031 GEOSCoordSeq_setY_r (handle, cs, 1, y);
3032 GEOSCoordSeq_setZ_r (handle, cs, 1, z);
3033 }
3034 else
3035 {
3036 GEOSCoordSeq_setX_r (handle, cs, 0, x0);
3037 GEOSCoordSeq_setY_r (handle, cs, 0, y0);
3038 GEOSCoordSeq_setX_r (handle, cs, 1, x);
3039 GEOSCoordSeq_setY_r (handle, cs, 1, y);
3040 }
3041 geos_item =
3042 GEOSGeom_createLineString_r (handle, cs);
3043 }
3044 else
3045 {
3046 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3047 cs = GEOSCoordSeq_create (2, dims);
3048 if (dims == 3)
3049 {
3050 GEOSCoordSeq_setX (cs, 0, x0);
3051 GEOSCoordSeq_setY (cs, 0, y0);
3052 GEOSCoordSeq_setZ (cs, 0, z0);
3053 GEOSCoordSeq_setX (cs, 1, x);
3054 GEOSCoordSeq_setY (cs, 1, y);
3055 GEOSCoordSeq_setZ (cs, 1, z);
3056 }
3057 else
3058 {
3059 GEOSCoordSeq_setX (cs, 0, x0);
3060 GEOSCoordSeq_setY (cs, 0, y0);
3061 GEOSCoordSeq_setX (cs, 1, x);
3062 GEOSCoordSeq_setY (cs, 1, y);
3063 }
3064 geos_item = GEOSGeom_createLineString (cs);
3065 #endif
3066 }
3067 *(geos_coll + nItem++) = geos_item;
3068 }
3069 x0 = x;
3070 y0 = y;
3071 z0 = z;
3072 }
3073 }
3074 pg = pg->Next;
3075 }
3076 if (handle != NULL)
3077 {
3078 geos =
3079 GEOSGeom_createCollection_r (handle, GEOS_MULTILINESTRING,
3080 geos_coll, segms);
3081 free (geos_coll);
3082 GEOSSetSRID_r (handle, geos, gaia->Srid);
3083 }
3084 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3085 else
3086 {
3087 geos =
3088 GEOSGeom_createCollection (GEOS_MULTILINESTRING, geos_coll,
3089 segms);
3090 free (geos_coll);
3091 GEOSSetSRID (geos, gaia->Srid);
3092 }
3093 #endif
3094 return geos;
3095 }
3096
3097 static gaiaGeomCollPtr
gaiaShortestLineCommon(struct splite_internal_cache * cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)3098 gaiaShortestLineCommon (struct splite_internal_cache *cache,
3099 gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
3100 {
3101 /* attempts to compute the shortest line between two geometries */
3102 GEOSGeometry *g1_points;
3103 GEOSGeometry *g1_segments;
3104 const GEOSGeometry *g1_item;
3105 GEOSGeometry *g2_points;
3106 GEOSGeometry *g2_segments;
3107 const GEOSGeometry *g2_item;
3108 const GEOSCoordSequence *cs;
3109 GEOSGeometry *g_pt;
3110 gaiaGeomCollPtr result = NULL;
3111 gaiaLinestringPtr ln;
3112 int nItems1;
3113 int nItems2;
3114 int it1;
3115 int it2;
3116 unsigned int dims;
3117 double x_ini = 0.0;
3118 double y_ini = 0.0;
3119 double z_ini = 0.0;
3120 double x_fin = 0.0;
3121 double y_fin = 0.0;
3122 double z_fin = 0.0;
3123 double dist;
3124 double min_dist = DBL_MAX;
3125 double projection;
3126 GEOSContextHandle_t handle = NULL;
3127
3128 if (cache != NULL)
3129 {
3130 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
3131 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
3132 return NULL;
3133 handle = cache->GEOS_handle;
3134 if (handle == NULL)
3135 return NULL;
3136 }
3137 #ifdef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3138 if (handle == NULL)
3139 return NULL;
3140 #endif
3141 if (!geom1 || !geom2)
3142 return NULL;
3143
3144 g1_points = buildGeosPoints (handle, geom1);
3145 g1_segments = buildGeosSegments (handle, geom1);
3146 g2_points = buildGeosPoints (handle, geom2);
3147 g2_segments = buildGeosSegments (handle, geom2);
3148
3149 if (g1_points && g2_points)
3150 {
3151 /* computing distances between POINTs */
3152 if (handle != NULL)
3153 {
3154 nItems1 = GEOSGetNumGeometries_r (handle, g1_points);
3155 nItems2 = GEOSGetNumGeometries_r (handle, g2_points);
3156 }
3157 else
3158 {
3159 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3160 nItems1 = GEOSGetNumGeometries (g1_points);
3161 nItems2 = GEOSGetNumGeometries (g2_points);
3162 #endif
3163 }
3164 for (it1 = 0; it1 < nItems1; it1++)
3165 {
3166 if (handle != NULL)
3167 g1_item = GEOSGetGeometryN_r (handle, g1_points, it1);
3168 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3169 else
3170 g1_item = GEOSGetGeometryN (g1_points, it1);
3171 #endif
3172 for (it2 = 0; it2 < nItems2; it2++)
3173 {
3174 int distret;
3175 if (handle != NULL)
3176 {
3177 g2_item =
3178 GEOSGetGeometryN_r (handle, g2_points, it2);
3179 distret =
3180 GEOSDistance_r (handle, g1_item, g2_item,
3181 &dist);
3182 }
3183 else
3184 {
3185 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3186 g2_item = GEOSGetGeometryN (g2_points, it2);
3187 distret = GEOSDistance (g1_item, g2_item, &dist);
3188 #endif
3189 }
3190 if (distret)
3191 {
3192 if (dist < min_dist)
3193 {
3194 /* saving min-dist points */
3195 min_dist = dist;
3196 if (handle != NULL)
3197 {
3198 cs = GEOSGeom_getCoordSeq_r (handle,
3199 g1_item);
3200 GEOSCoordSeq_getDimensions_r (handle,
3201 cs,
3202 &dims);
3203 if (dims == 3)
3204 {
3205 GEOSCoordSeq_getX_r (handle, cs,
3206 0, &x_ini);
3207 GEOSCoordSeq_getY_r (handle, cs,
3208 0, &y_ini);
3209 GEOSCoordSeq_getZ_r (handle, cs,
3210 0, &z_ini);
3211 }
3212 else
3213 {
3214 GEOSCoordSeq_getX_r (handle, cs,
3215 0, &x_ini);
3216 GEOSCoordSeq_getY_r (handle, cs,
3217 0, &y_ini);
3218 z_ini = 0.0;
3219 }
3220 cs = GEOSGeom_getCoordSeq_r (handle,
3221 g2_item);
3222 GEOSCoordSeq_getDimensions_r (handle,
3223 cs,
3224 &dims);
3225 if (dims == 3)
3226 {
3227 GEOSCoordSeq_getX_r (handle, cs,
3228 0, &x_fin);
3229 GEOSCoordSeq_getY_r (handle, cs,
3230 0, &y_fin);
3231 GEOSCoordSeq_getZ_r (handle, cs,
3232 0, &z_fin);
3233 }
3234 else
3235 {
3236 GEOSCoordSeq_getX_r (handle, cs,
3237 0, &x_fin);
3238 GEOSCoordSeq_getY_r (handle, cs,
3239 0, &y_fin);
3240 z_fin = 0.0;
3241 }
3242 }
3243 else
3244 {
3245 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3246 cs = GEOSGeom_getCoordSeq (g1_item);
3247 GEOSCoordSeq_getDimensions (cs, &dims);
3248 if (dims == 3)
3249 {
3250 GEOSCoordSeq_getX (cs, 0, &x_ini);
3251 GEOSCoordSeq_getY (cs, 0, &y_ini);
3252 GEOSCoordSeq_getZ (cs, 0, &z_ini);
3253 }
3254 else
3255 {
3256 GEOSCoordSeq_getX (cs, 0, &x_ini);
3257 GEOSCoordSeq_getY (cs, 0, &y_ini);
3258 z_ini = 0.0;
3259 }
3260 cs = GEOSGeom_getCoordSeq (g2_item);
3261 GEOSCoordSeq_getDimensions (cs, &dims);
3262 if (dims == 3)
3263 {
3264 GEOSCoordSeq_getX (cs, 0, &x_fin);
3265 GEOSCoordSeq_getY (cs, 0, &y_fin);
3266 GEOSCoordSeq_getZ (cs, 0, &z_fin);
3267 }
3268 else
3269 {
3270 GEOSCoordSeq_getX (cs, 0, &x_fin);
3271 GEOSCoordSeq_getY (cs, 0, &y_fin);
3272 z_fin = 0.0;
3273 }
3274 #endif
3275 }
3276 }
3277 }
3278 }
3279 }
3280 }
3281
3282 if (g1_points && g2_segments)
3283 {
3284 /* computing distances between POINTs (g1) and SEGMENTs (g2) */
3285 if (handle != NULL)
3286 {
3287 nItems1 = GEOSGetNumGeometries_r (handle, g1_points);
3288 nItems2 = GEOSGetNumGeometries_r (handle, g2_segments);
3289 }
3290 else
3291 {
3292 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3293 nItems1 = GEOSGetNumGeometries (g1_points);
3294 nItems2 = GEOSGetNumGeometries (g2_segments);
3295 #endif
3296 }
3297 for (it1 = 0; it1 < nItems1; it1++)
3298 {
3299 if (handle != NULL)
3300 g1_item = GEOSGetGeometryN_r (handle, g1_points, it1);
3301 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3302 else
3303 g1_item = GEOSGetGeometryN (g1_points, it1);
3304 #endif
3305 for (it2 = 0; it2 < nItems2; it2++)
3306 {
3307 int distret;
3308 if (handle != NULL)
3309 {
3310 g2_item =
3311 GEOSGetGeometryN_r (handle, g2_segments, it2);
3312 distret =
3313 GEOSDistance_r (handle, g1_item, g2_item,
3314 &dist);
3315 }
3316 else
3317 {
3318 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3319 g2_item = GEOSGetGeometryN (g2_segments, it2);
3320 distret = GEOSDistance (g1_item, g2_item, &dist);
3321 #endif
3322 }
3323 if (distret)
3324 {
3325 if (dist < min_dist)
3326 {
3327 /* saving min-dist points */
3328 if (handle != NULL)
3329 {
3330 projection =
3331 GEOSProject_r (handle, g2_item,
3332 g1_item);
3333 g_pt =
3334 GEOSInterpolate_r (handle, g2_item,
3335 projection);
3336 if (g_pt)
3337 {
3338 min_dist = dist;
3339 cs = GEOSGeom_getCoordSeq_r
3340 (handle, g1_item);
3341 GEOSCoordSeq_getDimensions_r
3342 (handle, cs, &dims);
3343 if (dims == 3)
3344 {
3345 GEOSCoordSeq_getX_r (handle,
3346 cs, 0,
3347 &x_ini);
3348 GEOSCoordSeq_getY_r (handle,
3349 cs, 0,
3350 &y_ini);
3351 GEOSCoordSeq_getZ_r (handle,
3352 cs, 0,
3353 &z_ini);
3354 }
3355 else
3356 {
3357 GEOSCoordSeq_getX_r (handle,
3358 cs, 0,
3359 &x_ini);
3360 GEOSCoordSeq_getY_r (handle,
3361 cs, 0,
3362 &y_ini);
3363 z_ini = 0.0;
3364 }
3365 cs = GEOSGeom_getCoordSeq_r
3366 (handle, g_pt);
3367 GEOSCoordSeq_getDimensions_r
3368 (handle, cs, &dims);
3369 if (dims == 3)
3370 {
3371 GEOSCoordSeq_getX_r (handle,
3372 cs, 0,
3373 &x_fin);
3374 GEOSCoordSeq_getY_r (handle,
3375 cs, 0,
3376 &y_fin);
3377 GEOSCoordSeq_getZ_r (handle,
3378 cs, 0,
3379 &z_fin);
3380 }
3381 else
3382 {
3383 GEOSCoordSeq_getX_r (handle,
3384 cs, 0,
3385 &x_fin);
3386 GEOSCoordSeq_getY_r (handle,
3387 cs, 0,
3388 &y_fin);
3389 z_fin = 0.0;
3390 }
3391 GEOSGeom_destroy_r (handle, g_pt);
3392 }
3393 }
3394 else
3395 {
3396 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3397 projection =
3398 GEOSProject (g2_item, g1_item);
3399 g_pt =
3400 GEOSInterpolate (g2_item,
3401 projection);
3402 if (g_pt)
3403 {
3404 min_dist = dist;
3405 cs = GEOSGeom_getCoordSeq
3406 (g1_item);
3407 GEOSCoordSeq_getDimensions (cs,
3408 &dims);
3409 if (dims == 3)
3410 {
3411 GEOSCoordSeq_getX (cs, 0,
3412 &x_ini);
3413 GEOSCoordSeq_getY (cs, 0,
3414 &y_ini);
3415 GEOSCoordSeq_getZ (cs, 0,
3416 &z_ini);
3417 }
3418 else
3419 {
3420 GEOSCoordSeq_getX (cs, 0,
3421 &x_ini);
3422 GEOSCoordSeq_getY (cs, 0,
3423 &y_ini);
3424 z_ini = 0.0;
3425 }
3426 cs = GEOSGeom_getCoordSeq (g_pt);
3427 GEOSCoordSeq_getDimensions (cs,
3428 &dims);
3429 if (dims == 3)
3430 {
3431 GEOSCoordSeq_getX (cs, 0,
3432 &x_fin);
3433 GEOSCoordSeq_getY (cs, 0,
3434 &y_fin);
3435 GEOSCoordSeq_getZ (cs, 0,
3436 &z_fin);
3437 }
3438 else
3439 {
3440 GEOSCoordSeq_getX (cs, 0,
3441 &x_fin);
3442 GEOSCoordSeq_getY (cs, 0,
3443 &y_fin);
3444 z_fin = 0.0;
3445 }
3446 GEOSGeom_destroy (g_pt);
3447 }
3448 #endif
3449 }
3450 }
3451 }
3452 }
3453 }
3454 }
3455
3456 if (g1_segments && g2_points)
3457 {
3458 /* computing distances between SEGMENTs (g1) and POINTs (g2) */
3459 if (handle != NULL)
3460 {
3461 nItems1 = GEOSGetNumGeometries_r (handle, g1_segments);
3462 nItems2 = GEOSGetNumGeometries_r (handle, g2_points);
3463 }
3464 else
3465 {
3466 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3467 nItems1 = GEOSGetNumGeometries (g1_segments);
3468 nItems2 = GEOSGetNumGeometries (g2_points);
3469 #endif
3470 }
3471 for (it1 = 0; it1 < nItems1; it1++)
3472 {
3473 if (handle != NULL)
3474 g1_item = GEOSGetGeometryN_r (handle, g1_segments, it1);
3475 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3476 else
3477 g1_item = GEOSGetGeometryN (g1_segments, it1);
3478 #endif
3479 for (it2 = 0; it2 < nItems2; it2++)
3480 {
3481 int distret;
3482 if (handle != NULL)
3483 {
3484 g2_item =
3485 GEOSGetGeometryN_r (handle, g2_points, it2);
3486 distret =
3487 GEOSDistance_r (handle, g1_item, g2_item,
3488 &dist);
3489 }
3490 else
3491 {
3492 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3493 g2_item = GEOSGetGeometryN (g2_points, it2);
3494 distret = GEOSDistance (g1_item, g2_item, &dist);
3495 #endif
3496 }
3497 if (distret)
3498 {
3499 if (dist < min_dist)
3500 {
3501 /* saving min-dist points */
3502 if (handle != NULL)
3503 {
3504 projection =
3505 GEOSProject_r (handle, g1_item,
3506 g2_item);
3507 g_pt =
3508 GEOSInterpolate_r (handle, g1_item,
3509 projection);
3510 if (g_pt)
3511 {
3512 min_dist = dist;
3513 cs = GEOSGeom_getCoordSeq_r
3514 (handle, g_pt);
3515 GEOSCoordSeq_getDimensions_r
3516 (handle, cs, &dims);
3517 if (dims == 3)
3518 {
3519 GEOSCoordSeq_getX_r (handle,
3520 cs, 0,
3521 &x_ini);
3522 GEOSCoordSeq_getY_r (handle,
3523 cs, 0,
3524 &y_ini);
3525 GEOSCoordSeq_getZ_r (handle,
3526 cs, 0,
3527 &z_ini);
3528 }
3529 else
3530 {
3531 GEOSCoordSeq_getX_r (handle,
3532 cs, 0,
3533 &x_ini);
3534 GEOSCoordSeq_getY_r (handle,
3535 cs, 0,
3536 &y_ini);
3537 z_ini = 0.0;
3538 }
3539 cs = GEOSGeom_getCoordSeq_r
3540 (handle, g2_item);
3541 GEOSCoordSeq_getDimensions_r
3542 (handle, cs, &dims);
3543 if (dims == 3)
3544 {
3545 GEOSCoordSeq_getX_r (handle,
3546 cs, 0,
3547 &x_fin);
3548 GEOSCoordSeq_getY_r (handle,
3549 cs, 0,
3550 &y_fin);
3551 GEOSCoordSeq_getZ_r (handle,
3552 cs, 0,
3553 &z_fin);
3554 }
3555 else
3556 {
3557 GEOSCoordSeq_getX_r (handle,
3558 cs, 0,
3559 &x_fin);
3560 GEOSCoordSeq_getY_r (handle,
3561 cs, 0,
3562 &y_fin);
3563 z_fin = 0.0;
3564 }
3565 GEOSGeom_destroy_r (handle, g_pt);
3566 }
3567 }
3568 else
3569 {
3570 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3571 projection =
3572 GEOSProject (g1_item, g2_item);
3573 g_pt =
3574 GEOSInterpolate (g1_item,
3575 projection);
3576 if (g_pt)
3577 {
3578 min_dist = dist;
3579 cs = GEOSGeom_getCoordSeq (g_pt);
3580 GEOSCoordSeq_getDimensions (cs,
3581 &dims);
3582 if (dims == 3)
3583 {
3584 GEOSCoordSeq_getX (cs, 0,
3585 &x_ini);
3586 GEOSCoordSeq_getY (cs, 0,
3587 &y_ini);
3588 GEOSCoordSeq_getZ (cs, 0,
3589 &z_ini);
3590 }
3591 else
3592 {
3593 GEOSCoordSeq_getX (cs, 0,
3594 &x_ini);
3595 GEOSCoordSeq_getY (cs, 0,
3596 &y_ini);
3597 z_ini = 0.0;
3598 }
3599 cs = GEOSGeom_getCoordSeq
3600 (g2_item);
3601 GEOSCoordSeq_getDimensions (cs,
3602 &dims);
3603 if (dims == 3)
3604 {
3605 GEOSCoordSeq_getX (cs, 0,
3606 &x_fin);
3607 GEOSCoordSeq_getY (cs, 0,
3608 &y_fin);
3609 GEOSCoordSeq_getZ (cs, 0,
3610 &z_fin);
3611 }
3612 else
3613 {
3614 GEOSCoordSeq_getX (cs, 0,
3615 &x_fin);
3616 GEOSCoordSeq_getY (cs, 0,
3617 &y_fin);
3618 z_fin = 0.0;
3619 }
3620 GEOSGeom_destroy (g_pt);
3621 }
3622 #endif
3623 }
3624 }
3625 }
3626 }
3627 }
3628 }
3629 if (handle != NULL)
3630 {
3631 if (g1_points)
3632 GEOSGeom_destroy_r (handle, g1_points);
3633 if (g1_segments)
3634 GEOSGeom_destroy_r (handle, g1_segments);
3635 if (g2_points)
3636 GEOSGeom_destroy_r (handle, g2_points);
3637 if (g2_segments)
3638 GEOSGeom_destroy_r (handle, g2_segments);
3639 }
3640 else
3641 {
3642 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3643 if (g1_points)
3644 GEOSGeom_destroy (g1_points);
3645 if (g1_segments)
3646 GEOSGeom_destroy (g1_segments);
3647 if (g2_points)
3648 GEOSGeom_destroy (g2_points);
3649 if (g2_segments)
3650 GEOSGeom_destroy (g2_segments);
3651 #endif
3652 }
3653 if (min_dist == DBL_MAX || min_dist <= 0.0)
3654 return NULL;
3655
3656 /* building the shortest line */
3657 switch (geom1->DimensionModel)
3658 {
3659 case GAIA_XY_Z:
3660 result = gaiaAllocGeomCollXYZ ();
3661 break;
3662 case GAIA_XY_M:
3663 result = gaiaAllocGeomCollXYM ();
3664 break;
3665 case GAIA_XY_Z_M:
3666 result = gaiaAllocGeomCollXYZM ();
3667 break;
3668 default:
3669 result = gaiaAllocGeomColl ();
3670 break;
3671 };
3672 result->Srid = geom1->Srid;
3673 ln = gaiaAddLinestringToGeomColl (result, 2);
3674 switch (ln->DimensionModel)
3675 {
3676 case GAIA_XY_Z:
3677 gaiaSetPointXYZ (ln->Coords, 0, x_ini, y_ini, z_ini);
3678 gaiaSetPointXYZ (ln->Coords, 1, x_fin, y_fin, z_fin);
3679 break;
3680 case GAIA_XY_M:
3681 gaiaSetPointXYM (ln->Coords, 0, x_ini, y_ini, 0.0);
3682 gaiaSetPointXYM (ln->Coords, 1, x_fin, y_fin, 0.0);
3683 break;
3684 case GAIA_XY_Z_M:
3685 gaiaSetPointXYZM (ln->Coords, 0, x_ini, y_ini, z_ini, 0.0);
3686 gaiaSetPointXYZM (ln->Coords, 1, x_fin, y_fin, z_fin, 0.0);
3687 break;
3688 default:
3689 gaiaSetPoint (ln->Coords, 0, x_ini, y_ini);
3690 gaiaSetPoint (ln->Coords, 1, x_fin, y_fin);
3691 break;
3692 };
3693 return result;
3694 }
3695
3696 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaShortestLine(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)3697 gaiaShortestLine (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
3698 {
3699 gaiaResetGeosMsg ();
3700 return gaiaShortestLineCommon (NULL, geom1, geom2);
3701 }
3702
3703 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaShortestLine_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)3704 gaiaShortestLine_r (const void *p_cache, gaiaGeomCollPtr geom1,
3705 gaiaGeomCollPtr geom2)
3706 {
3707 struct splite_internal_cache *cache =
3708 (struct splite_internal_cache *) p_cache;
3709 GEOSContextHandle_t handle = NULL;
3710 if (cache == NULL)
3711 return NULL;
3712 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
3713 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
3714 return NULL;
3715 handle = cache->GEOS_handle;
3716 if (handle == NULL)
3717 return NULL;
3718 gaiaResetGeosMsg_r (cache);
3719 return gaiaShortestLineCommon (cache, geom1, geom2);
3720 }
3721
3722 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSnap(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double tolerance)3723 gaiaSnap (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2, double tolerance)
3724 {
3725 /* attempts to "snap" geom1 on geom2 using the given tolerance */
3726 gaiaGeomCollPtr result = NULL;
3727 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3728 GEOSGeometry *g1;
3729 GEOSGeometry *g2;
3730 GEOSGeometry *g3;
3731 gaiaResetGeosMsg ();
3732 if (!geom1 || !geom2)
3733 return NULL;
3734
3735 g1 = gaiaToGeos (geom1);
3736 g2 = gaiaToGeos (geom2);
3737 g3 = GEOSSnap (g1, g2, tolerance);
3738 GEOSGeom_destroy (g1);
3739 GEOSGeom_destroy (g2);
3740 if (!g3)
3741 return NULL;
3742 if (geom1->DimensionModel == GAIA_XY_Z)
3743 result = gaiaFromGeos_XYZ (g3);
3744 else if (geom1->DimensionModel == GAIA_XY_M)
3745 result = gaiaFromGeos_XYM (g3);
3746 else if (geom1->DimensionModel == GAIA_XY_Z_M)
3747 result = gaiaFromGeos_XYZM (g3);
3748 else
3749 result = gaiaFromGeos_XY (g3);
3750 GEOSGeom_destroy (g3);
3751 if (result == NULL)
3752 return NULL;
3753 result->Srid = geom1->Srid;
3754 #else
3755 if (geom1 == NULL || geom2 == NULL || tolerance == 0.0)
3756 geom1 = NULL; /* silencing stupid compiler warnings */
3757 #endif
3758 return result;
3759 }
3760
3761 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSnap_r(const void * p_cache,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double tolerance)3762 gaiaSnap_r (const void *p_cache, gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
3763 double tolerance)
3764 {
3765 /* attempts to "snap" geom1 on geom2 using the given tolerance */
3766 GEOSGeometry *g1;
3767 GEOSGeometry *g2;
3768 GEOSGeometry *g3;
3769 gaiaGeomCollPtr result;
3770 struct splite_internal_cache *cache =
3771 (struct splite_internal_cache *) p_cache;
3772 GEOSContextHandle_t handle = NULL;
3773 if (cache == NULL)
3774 return NULL;
3775 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
3776 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
3777 return NULL;
3778 handle = cache->GEOS_handle;
3779 if (handle == NULL)
3780 return NULL;
3781 gaiaResetGeosMsg_r (cache);
3782 if (!geom1 || !geom2)
3783 return NULL;
3784
3785 g1 = gaiaToGeos_r (cache, geom1);
3786 g2 = gaiaToGeos_r (cache, geom2);
3787 g3 = GEOSSnap_r (handle, g1, g2, tolerance);
3788 GEOSGeom_destroy_r (handle, g1);
3789 GEOSGeom_destroy_r (handle, g2);
3790 if (!g3)
3791 return NULL;
3792 if (geom1->DimensionModel == GAIA_XY_Z)
3793 result = gaiaFromGeos_XYZ_r (cache, g3);
3794 else if (geom1->DimensionModel == GAIA_XY_M)
3795 result = gaiaFromGeos_XYM_r (cache, g3);
3796 else if (geom1->DimensionModel == GAIA_XY_Z_M)
3797 result = gaiaFromGeos_XYZM_r (cache, g3);
3798 else
3799 result = gaiaFromGeos_XY_r (cache, g3);
3800 GEOSGeom_destroy_r (handle, g3);
3801 if (result == NULL)
3802 return NULL;
3803 result->Srid = geom1->Srid;
3804 return result;
3805 }
3806
3807 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineMerge(gaiaGeomCollPtr geom)3808 gaiaLineMerge (gaiaGeomCollPtr geom)
3809 {
3810 /* attempts to reassemble lines from a collection of sparse fragments */
3811 gaiaGeomCollPtr result = NULL;
3812 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3813 GEOSGeometry *g1;
3814 GEOSGeometry *g2;
3815 gaiaResetGeosMsg ();
3816 if (!geom)
3817 return NULL;
3818 if (gaiaIsToxic (geom))
3819 return NULL;
3820
3821 g1 = gaiaToGeos (geom);
3822 g2 = GEOSLineMerge (g1);
3823 GEOSGeom_destroy (g1);
3824 if (!g2)
3825 return NULL;
3826 if (geom->DimensionModel == GAIA_XY_Z)
3827 result = gaiaFromGeos_XYZ (g2);
3828 else if (geom->DimensionModel == GAIA_XY_M)
3829 result = gaiaFromGeos_XYM (g2);
3830 else if (geom->DimensionModel == GAIA_XY_Z_M)
3831 result = gaiaFromGeos_XYZM (g2);
3832 else
3833 result = gaiaFromGeos_XY (g2);
3834 GEOSGeom_destroy (g2);
3835 if (result == NULL)
3836 return NULL;
3837 result->Srid = geom->Srid;
3838 #else
3839 if (geom == NULL)
3840 geom = NULL; /* silencing stupid compiler warnings */
3841 #endif
3842 return result;
3843 }
3844
3845 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLineMerge_r(const void * p_cache,gaiaGeomCollPtr geom)3846 gaiaLineMerge_r (const void *p_cache, gaiaGeomCollPtr geom)
3847 {
3848 /* attempts to reassemble lines from a collection of sparse fragments */
3849 GEOSGeometry *g1;
3850 GEOSGeometry *g2;
3851 gaiaGeomCollPtr result;
3852 struct splite_internal_cache *cache =
3853 (struct splite_internal_cache *) p_cache;
3854 GEOSContextHandle_t handle = NULL;
3855 if (cache == NULL)
3856 return NULL;
3857 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
3858 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
3859 return NULL;
3860 handle = cache->GEOS_handle;
3861 if (handle == NULL)
3862 return NULL;
3863 gaiaResetGeosMsg_r (cache);
3864 if (!geom)
3865 return NULL;
3866 if (gaiaIsToxic_r (cache, geom))
3867 return NULL;
3868
3869 g1 = gaiaToGeos_r (cache, geom);
3870 g2 = GEOSLineMerge_r (handle, g1);
3871 GEOSGeom_destroy_r (handle, g1);
3872 if (!g2)
3873 return NULL;
3874 if (geom->DimensionModel == GAIA_XY_Z)
3875 result = gaiaFromGeos_XYZ_r (cache, g2);
3876 else if (geom->DimensionModel == GAIA_XY_M)
3877 result = gaiaFromGeos_XYM_r (cache, g2);
3878 else if (geom->DimensionModel == GAIA_XY_Z_M)
3879 result = gaiaFromGeos_XYZM_r (cache, g2);
3880 else
3881 result = gaiaFromGeos_XY_r (cache, g2);
3882 GEOSGeom_destroy_r (handle, g2);
3883 if (result == NULL)
3884 return NULL;
3885 result->Srid = geom->Srid;
3886 return result;
3887 }
3888
3889 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaUnaryUnion(gaiaGeomCollPtr geom)3890 gaiaUnaryUnion (gaiaGeomCollPtr geom)
3891 {
3892 /* Unary Union (single Collection) */
3893 gaiaGeomCollPtr result = NULL;
3894 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
3895 GEOSGeometry *g1;
3896 GEOSGeometry *g2;
3897 gaiaResetGeosMsg ();
3898 if (!geom)
3899 return NULL;
3900 if (gaiaIsToxic (geom))
3901 return NULL;
3902 g1 = gaiaToGeos (geom);
3903 g2 = GEOSUnaryUnion (g1);
3904 GEOSGeom_destroy (g1);
3905 if (!g2)
3906 return NULL;
3907 if (geom->DimensionModel == GAIA_XY_Z)
3908 result = gaiaFromGeos_XYZ (g2);
3909 else if (geom->DimensionModel == GAIA_XY_M)
3910 result = gaiaFromGeos_XYM (g2);
3911 else if (geom->DimensionModel == GAIA_XY_Z_M)
3912 result = gaiaFromGeos_XYZM (g2);
3913 else
3914 result = gaiaFromGeos_XY (g2);
3915 GEOSGeom_destroy (g2);
3916 if (result == NULL)
3917 return NULL;
3918 result->Srid = geom->Srid;
3919 #else
3920 if (geom == NULL)
3921 geom = NULL; /* silencing stupid compiler warnings */
3922 #endif
3923 return result;
3924 }
3925
3926 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaUnaryUnion_r(const void * p_cache,gaiaGeomCollPtr geom)3927 gaiaUnaryUnion_r (const void *p_cache, gaiaGeomCollPtr geom)
3928 {
3929 /* Unary Union (single Collection) */
3930 GEOSGeometry *g1;
3931 GEOSGeometry *g2;
3932 gaiaGeomCollPtr result;
3933 struct splite_internal_cache *cache =
3934 (struct splite_internal_cache *) p_cache;
3935 GEOSContextHandle_t handle = NULL;
3936 if (cache == NULL)
3937 return NULL;
3938 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
3939 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
3940 return NULL;
3941 handle = cache->GEOS_handle;
3942 if (handle == NULL)
3943 return NULL;
3944 gaiaResetGeosMsg_r (cache);
3945 if (!geom)
3946 return NULL;
3947 if (gaiaIsToxic_r (cache, geom))
3948 return NULL;
3949 g1 = gaiaToGeos_r (cache, geom);
3950 g2 = GEOSUnaryUnion_r (handle, g1);
3951 GEOSGeom_destroy_r (handle, g1);
3952 if (!g2)
3953 return NULL;
3954 if (geom->DimensionModel == GAIA_XY_Z)
3955 result = gaiaFromGeos_XYZ_r (cache, g2);
3956 else if (geom->DimensionModel == GAIA_XY_M)
3957 result = gaiaFromGeos_XYM_r (cache, g2);
3958 else if (geom->DimensionModel == GAIA_XY_Z_M)
3959 result = gaiaFromGeos_XYZM_r (cache, g2);
3960 else
3961 result = gaiaFromGeos_XY_r (cache, g2);
3962 GEOSGeom_destroy_r (handle, g2);
3963 if (result == NULL)
3964 return NULL;
3965 result->Srid = geom->Srid;
3966 return result;
3967 }
3968
3969 static void
rotateRingBeforeCut(gaiaLinestringPtr ln,gaiaPointPtr node)3970 rotateRingBeforeCut (gaiaLinestringPtr ln, gaiaPointPtr node)
3971 {
3972 /* rotating a Ring, so to ensure that Start/End points match the node */
3973 int io = 0;
3974 int iv;
3975 int copy = 0;
3976 int base_idx = -1;
3977 double x = 0.0;
3978 double y = 0.0;
3979 double z = 0.0;
3980 double m = 0.0;
3981 gaiaLinestringPtr new_ln = NULL;
3982
3983 if (ln->DimensionModel == GAIA_XY_Z)
3984 new_ln = gaiaAllocLinestringXYZ (ln->Points);
3985 else if (ln->DimensionModel == GAIA_XY_M)
3986 new_ln = gaiaAllocLinestringXYM (ln->Points);
3987 else if (ln->DimensionModel == GAIA_XY_Z_M)
3988 new_ln = gaiaAllocLinestringXYZM (ln->Points);
3989 else
3990 new_ln = gaiaAllocLinestring (ln->Points);
3991
3992 /* first pass */
3993 for (iv = 0; iv < ln->Points; iv++)
3994 {
3995 if (ln->DimensionModel == GAIA_XY_Z)
3996 {
3997 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3998 }
3999 else if (ln->DimensionModel == GAIA_XY_M)
4000 {
4001 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4002 }
4003 else if (ln->DimensionModel == GAIA_XY_Z_M)
4004 {
4005 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4006 }
4007 else
4008 {
4009 gaiaGetPoint (ln->Coords, iv, &x, &y);
4010 }
4011 if (!copy)
4012 {
4013 if (ln->DimensionModel == GAIA_XY_Z
4014 || ln->DimensionModel == GAIA_XY_Z_M)
4015 {
4016 if (node->X == x && node->Y == y && node->Z == z)
4017 {
4018 base_idx = iv;
4019 copy = 1;
4020 }
4021 }
4022 else if (node->X == x && node->Y == y)
4023 {
4024 base_idx = iv;
4025 copy = 1;
4026 }
4027 }
4028 if (copy)
4029 {
4030 /* copying points */
4031 if (ln->DimensionModel == GAIA_XY_Z)
4032 {
4033 gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4034 }
4035 else if (ln->DimensionModel == GAIA_XY_M)
4036 {
4037 gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4038 }
4039 else if (ln->DimensionModel == GAIA_XY_Z_M)
4040 {
4041 gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4042 }
4043 else
4044 {
4045 gaiaSetPoint (new_ln->Coords, io, x, y);
4046 }
4047 io++;
4048 }
4049 }
4050 if (base_idx <= 0)
4051 {
4052 gaiaFreeLinestring (new_ln);
4053 return;
4054 }
4055
4056 /* second pass */
4057 for (iv = 1; iv <= base_idx; iv++)
4058 {
4059 if (ln->DimensionModel == GAIA_XY_Z)
4060 {
4061 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4062 }
4063 else if (ln->DimensionModel == GAIA_XY_M)
4064 {
4065 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4066 }
4067 else if (ln->DimensionModel == GAIA_XY_Z_M)
4068 {
4069 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4070 }
4071 else
4072 {
4073 gaiaGetPoint (ln->Coords, iv, &x, &y);
4074 }
4075 if (ln->DimensionModel == GAIA_XY_Z)
4076 {
4077 gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4078 }
4079 else if (ln->DimensionModel == GAIA_XY_M)
4080 {
4081 gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4082 }
4083 else if (ln->DimensionModel == GAIA_XY_Z_M)
4084 {
4085 gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4086 }
4087 else
4088 {
4089 gaiaSetPoint (new_ln->Coords, io, x, y);
4090 }
4091 io++;
4092 }
4093
4094 /* copying back */
4095 for (iv = 0; iv < new_ln->Points; iv++)
4096 {
4097 if (ln->DimensionModel == GAIA_XY_Z)
4098 {
4099 gaiaGetPointXYZ (new_ln->Coords, iv, &x, &y, &z);
4100 }
4101 else if (ln->DimensionModel == GAIA_XY_M)
4102 {
4103 gaiaGetPointXYM (new_ln->Coords, iv, &x, &y, &m);
4104 }
4105 else if (ln->DimensionModel == GAIA_XY_Z_M)
4106 {
4107 gaiaGetPointXYZM (new_ln->Coords, iv, &x, &y, &z, &m);
4108 }
4109 else
4110 {
4111 gaiaGetPoint (new_ln->Coords, iv, &x, &y);
4112 }
4113 if (ln->DimensionModel == GAIA_XY_Z)
4114 {
4115 gaiaSetPointXYZ (ln->Coords, iv, x, y, z);
4116 }
4117 else if (ln->DimensionModel == GAIA_XY_M)
4118 {
4119 gaiaSetPointXYM (ln->Coords, iv, x, y, m);
4120 }
4121 else if (ln->DimensionModel == GAIA_XY_Z_M)
4122 {
4123 gaiaSetPointXYZM (ln->Coords, iv, x, y, z, m);
4124 }
4125 else
4126 {
4127 gaiaSetPoint (ln->Coords, iv, x, y);
4128 }
4129 }
4130 gaiaFreeLinestring (new_ln);
4131 }
4132
4133 static void
extractSubLine(gaiaGeomCollPtr result,gaiaLinestringPtr ln,int i_start,int i_end)4134 extractSubLine (gaiaGeomCollPtr result, gaiaLinestringPtr ln, int i_start,
4135 int i_end)
4136 {
4137 /* extracting s SubLine */
4138 int iv;
4139 int io = 0;
4140 int pts = i_end - i_start + 1;
4141 gaiaLinestringPtr new_ln = NULL;
4142 double x;
4143 double y;
4144 double z;
4145 double m;
4146
4147 new_ln = gaiaAddLinestringToGeomColl (result, pts);
4148
4149 for (iv = i_start; iv <= i_end; iv++)
4150 {
4151 m = 0.0;
4152 z = 0.0;
4153 if (ln->DimensionModel == GAIA_XY_Z)
4154 {
4155 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4156 }
4157 else if (ln->DimensionModel == GAIA_XY_M)
4158 {
4159 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4160 }
4161 else if (ln->DimensionModel == GAIA_XY_Z_M)
4162 {
4163 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4164 }
4165 else
4166 {
4167 gaiaGetPoint (ln->Coords, iv, &x, &y);
4168 }
4169 if (ln->DimensionModel == GAIA_XY_Z)
4170 {
4171 gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4172 }
4173 else if (ln->DimensionModel == GAIA_XY_M)
4174 {
4175 gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4176 }
4177 else if (ln->DimensionModel == GAIA_XY_Z_M)
4178 {
4179 gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4180 }
4181 else
4182 {
4183 gaiaSetPoint (new_ln->Coords, io, x, y);
4184 }
4185 io++;
4186 }
4187 }
4188
4189 static void
cutLineAtNodes(gaiaLinestringPtr ln,gaiaPointPtr pt_base,gaiaGeomCollPtr result)4190 cutLineAtNodes (gaiaLinestringPtr ln, gaiaPointPtr pt_base,
4191 gaiaGeomCollPtr result)
4192 {
4193 /* attempts to cut a single Line accordingly to given nodes */
4194 int closed = 0;
4195 int match = 0;
4196 int iv;
4197 int i_start;
4198 double x = 0.0;
4199 double y = 0.0;
4200 double z = 0.0;
4201 double m = 0.0;
4202 gaiaPointPtr pt;
4203 gaiaPointPtr node = NULL;
4204
4205 if (gaiaIsClosed (ln))
4206 closed = 1;
4207 /* pre-check */
4208 for (iv = 0; iv < ln->Points; iv++)
4209 {
4210 if (ln->DimensionModel == GAIA_XY_Z)
4211 {
4212 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4213 }
4214 else if (ln->DimensionModel == GAIA_XY_M)
4215 {
4216 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4217 }
4218 else if (ln->DimensionModel == GAIA_XY_Z_M)
4219 {
4220 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4221 }
4222 else
4223 {
4224 gaiaGetPoint (ln->Coords, iv, &x, &y);
4225 }
4226 pt = pt_base;
4227 while (pt)
4228 {
4229 if (ln->DimensionModel == GAIA_XY_Z
4230 || ln->DimensionModel == GAIA_XY_Z_M)
4231 {
4232 if (pt->X == x && pt->Y == y && pt->Z == z)
4233 {
4234 node = pt;
4235 match++;
4236 }
4237 }
4238 else if (pt->X == x && pt->Y == y)
4239 {
4240 node = pt;
4241 match++;
4242 }
4243 pt = pt->Next;
4244 }
4245 }
4246
4247 if (closed && node)
4248 rotateRingBeforeCut (ln, node);
4249
4250 i_start = 0;
4251 for (iv = 1; iv < ln->Points - 1; iv++)
4252 {
4253 /* identifying sub-linestrings */
4254 if (ln->DimensionModel == GAIA_XY_Z)
4255 {
4256 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4257 }
4258 else if (ln->DimensionModel == GAIA_XY_M)
4259 {
4260 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4261 }
4262 else if (ln->DimensionModel == GAIA_XY_Z_M)
4263 {
4264 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4265 }
4266 else
4267 {
4268 gaiaGetPoint (ln->Coords, iv, &x, &y);
4269 }
4270 match = 0;
4271 pt = pt_base;
4272 while (pt)
4273 {
4274 if (ln->DimensionModel == GAIA_XY_Z
4275 || ln->DimensionModel == GAIA_XY_Z_M)
4276 {
4277 if (pt->X == x && pt->Y == y && pt->Z == z)
4278 {
4279 match = 1;
4280 break;
4281 }
4282 }
4283 else if (pt->X == x && pt->Y == y)
4284 {
4285 match = 1;
4286 break;
4287 }
4288 pt = pt->Next;
4289 }
4290 if (match)
4291 {
4292 /* cutting the line */
4293 extractSubLine (result, ln, i_start, iv);
4294 i_start = iv;
4295 }
4296 }
4297 if (i_start != 0 && i_start != ln->Points - 1)
4298 {
4299 /* extracting the last SubLine */
4300 extractSubLine (result, ln, i_start, ln->Points - 1);
4301 }
4302 else
4303 {
4304 /* cloning the untouched Line */
4305 extractSubLine (result, ln, 0, ln->Points - 1);
4306 }
4307 }
4308
4309 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLinesCutAtNodes(gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2)4310 gaiaLinesCutAtNodes (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
4311 {
4312 /* attempts to cut lines accordingly to nodes */
4313 int pts1 = 0;
4314 int lns1 = 0;
4315 int pgs1 = 0;
4316 int pts2 = 0;
4317 int lns2 = 0;
4318 int pgs2 = 0;
4319 gaiaPointPtr pt;
4320 gaiaLinestringPtr ln;
4321 gaiaPolygonPtr pg;
4322 gaiaGeomCollPtr result = NULL;
4323
4324 if (!geom1)
4325 return NULL;
4326 if (!geom2)
4327 return NULL;
4328
4329 /* both Geometryes should have identical Dimensions */
4330 if (geom1->DimensionModel != geom2->DimensionModel)
4331 return NULL;
4332
4333 pt = geom1->FirstPoint;
4334 while (pt)
4335 {
4336 pts1++;
4337 pt = pt->Next;
4338 }
4339 ln = geom1->FirstLinestring;
4340 while (ln)
4341 {
4342 lns1++;
4343 ln = ln->Next;
4344 }
4345 pg = geom1->FirstPolygon;
4346 while (pg)
4347 {
4348 pgs1++;
4349 pg = pg->Next;
4350 }
4351 pt = geom2->FirstPoint;
4352 while (pt)
4353 {
4354 pts2++;
4355 pt = pt->Next;
4356 }
4357 ln = geom2->FirstLinestring;
4358 while (ln)
4359 {
4360 lns2++;
4361 ln = ln->Next;
4362 }
4363 pg = geom2->FirstPolygon;
4364 while (pg)
4365 {
4366 pgs2++;
4367 pg = pg->Next;
4368 }
4369
4370 /* the first Geometry is expected to contain one or more Linestring(s) */
4371 if (pts1 == 0 && lns1 > 0 && pgs1 == 0)
4372 ;
4373 else
4374 return NULL;
4375 /* the second Geometry is expected to contain one or more Point(s) */
4376 if (pts2 > 0 && lns2 == 0 && pgs2 == 0)
4377 ;
4378 else
4379 return NULL;
4380
4381 /* attempting to cut Lines accordingly to Nodes */
4382 if (geom1->DimensionModel == GAIA_XY_Z)
4383 result = gaiaAllocGeomCollXYZ ();
4384 else if (geom1->DimensionModel == GAIA_XY_M)
4385 result = gaiaAllocGeomCollXYM ();
4386 else if (geom1->DimensionModel == GAIA_XY_Z_M)
4387 result = gaiaAllocGeomCollXYZM ();
4388 else
4389 result = gaiaAllocGeomColl ();
4390 ln = geom1->FirstLinestring;
4391 while (ln)
4392 {
4393 cutLineAtNodes (ln, geom2->FirstPoint, result);
4394 ln = ln->Next;
4395 }
4396 if (result->FirstLinestring == NULL)
4397 {
4398 gaiaFreeGeomColl (result);
4399 return NULL;
4400 }
4401 result->Srid = geom1->Srid;
4402 return result;
4403 }
4404
4405 #ifdef GEOS_ADVANCED /* GEOS advanced features - 3.4.0 */
4406
4407 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDelaunayTriangulation(gaiaGeomCollPtr geom,double tolerance,int only_edges)4408 gaiaDelaunayTriangulation (gaiaGeomCollPtr geom, double tolerance,
4409 int only_edges)
4410 {
4411 /* Delaunay Triangulation */
4412 gaiaGeomCollPtr result = NULL;
4413 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
4414 GEOSGeometry *g1;
4415 GEOSGeometry *g2;
4416 gaiaResetGeosMsg ();
4417 if (!geom)
4418 return NULL;
4419 g1 = gaiaToGeos (geom);
4420 g2 = GEOSDelaunayTriangulation (g1, tolerance, only_edges);
4421 GEOSGeom_destroy (g1);
4422 if (!g2)
4423 return NULL;
4424 if (geom->DimensionModel == GAIA_XY_Z)
4425 result = gaiaFromGeos_XYZ (g2);
4426 else if (geom->DimensionModel == GAIA_XY_M)
4427 result = gaiaFromGeos_XYM (g2);
4428 else if (geom->DimensionModel == GAIA_XY_Z_M)
4429 result = gaiaFromGeos_XYZM (g2);
4430 else
4431 result = gaiaFromGeos_XY (g2);
4432 GEOSGeom_destroy (g2);
4433 if (result == NULL)
4434 return NULL;
4435 result->Srid = geom->Srid;
4436 if (only_edges)
4437 result->DeclaredType = GAIA_MULTILINESTRING;
4438 else
4439 result->DeclaredType = GAIA_MULTIPOLYGON;
4440 #else
4441 if (geom == NULL || tolerance == 0.0 || only_edges == 0)
4442 geom = NULL; /* silencing stupid compiler warnings */
4443 #endif
4444 return result;
4445 }
4446
4447 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDelaunayTriangulation_r(const void * p_cache,gaiaGeomCollPtr geom,double tolerance,int only_edges)4448 gaiaDelaunayTriangulation_r (const void *p_cache, gaiaGeomCollPtr geom,
4449 double tolerance, int only_edges)
4450 {
4451 /* Delaunay Triangulation */
4452 GEOSGeometry *g1;
4453 GEOSGeometry *g2;
4454 gaiaGeomCollPtr result;
4455 struct splite_internal_cache *cache =
4456 (struct splite_internal_cache *) p_cache;
4457 GEOSContextHandle_t handle = NULL;
4458 if (cache == NULL)
4459 return NULL;
4460 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
4461 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
4462 return NULL;
4463 handle = cache->GEOS_handle;
4464 if (handle == NULL)
4465 return NULL;
4466 gaiaResetGeosMsg_r (cache);
4467 if (!geom)
4468 return NULL;
4469 g1 = gaiaToGeos_r (cache, geom);
4470 g2 = GEOSDelaunayTriangulation_r (handle, g1, tolerance, only_edges);
4471 GEOSGeom_destroy_r (handle, g1);
4472 if (!g2)
4473 return NULL;
4474 if (geom->DimensionModel == GAIA_XY_Z)
4475 result = gaiaFromGeos_XYZ_r (cache, g2);
4476 else if (geom->DimensionModel == GAIA_XY_M)
4477 result = gaiaFromGeos_XYM_r (cache, g2);
4478 else if (geom->DimensionModel == GAIA_XY_Z_M)
4479 result = gaiaFromGeos_XYZM_r (cache, g2);
4480 else
4481 result = gaiaFromGeos_XY_r (cache, g2);
4482 GEOSGeom_destroy_r (handle, g2);
4483 if (result == NULL)
4484 return NULL;
4485 result->Srid = geom->Srid;
4486 if (only_edges)
4487 result->DeclaredType = GAIA_MULTILINESTRING;
4488 else
4489 result->DeclaredType = GAIA_MULTIPOLYGON;
4490 return result;
4491 }
4492
4493 #ifdef GEOS_REENTRANT /* only if GEOS >= 3.5.0 directly supporting Voronoj */
4494 static gaiaGeomCollPtr
voronoj_envelope(gaiaGeomCollPtr geom,double extra_frame_size)4495 voronoj_envelope (gaiaGeomCollPtr geom, double extra_frame_size)
4496 {
4497 /* building the extended envelope for Voronoj */
4498 gaiaGeomCollPtr bbox;
4499 gaiaPolygonPtr pg;
4500 gaiaRingPtr rect;
4501 double minx;
4502 double miny;
4503 double maxx;
4504 double maxy;
4505 double ext_x;
4506 double ext_y;
4507 double delta;
4508 double delta2;
4509
4510 gaiaMbrGeometry (geom);
4511 /* setting the frame extent */
4512 if (extra_frame_size < 0.0)
4513 extra_frame_size = 5.0;
4514 ext_x = geom->MaxX - geom->MinX;
4515 ext_y = geom->MaxY - geom->MinY;
4516 delta = (ext_x * extra_frame_size) / 100.0;
4517 delta2 = (ext_y * extra_frame_size) / 100.0;
4518 if (delta2 > delta)
4519 delta = delta2;
4520 minx = geom->MinX - delta;
4521 miny = geom->MinY - delta;
4522 maxx = geom->MaxX + delta;
4523 maxy = geom->MaxY + delta;
4524
4525 /* building the frame */
4526 if (geom->DimensionModel == GAIA_XY_Z)
4527 bbox = gaiaAllocGeomCollXYZ ();
4528 else if (geom->DimensionModel == GAIA_XY_M)
4529 bbox = gaiaAllocGeomCollXYM ();
4530 else if (geom->DimensionModel == GAIA_XY_Z_M)
4531 bbox = gaiaAllocGeomCollXYZM ();
4532 else
4533 bbox = gaiaAllocGeomColl ();
4534 bbox->Srid = geom->Srid;
4535 bbox->DeclaredType = GAIA_POLYGON;
4536 pg = gaiaAddPolygonToGeomColl (bbox, 5, 0);
4537 rect = pg->Exterior;
4538 if (geom->DimensionModel == GAIA_XY_Z)
4539 {
4540 gaiaSetPointXYZ (rect->Coords, 0, minx, miny, 0.0);
4541 gaiaSetPointXYZ (rect->Coords, 1, maxx, miny, 0.0);
4542 gaiaSetPointXYZ (rect->Coords, 2, maxx, maxy, 0.0);
4543 gaiaSetPointXYZ (rect->Coords, 3, minx, maxy, 0.0);
4544 gaiaSetPointXYZ (rect->Coords, 4, minx, miny, 0.0);
4545 }
4546 else if (geom->DimensionModel == GAIA_XY_M)
4547 {
4548 gaiaSetPointXYM (rect->Coords, 0, minx, miny, 0.0);
4549 gaiaSetPointXYM (rect->Coords, 1, maxx, miny, 0.0);
4550 gaiaSetPointXYM (rect->Coords, 2, maxx, maxy, 0.0);
4551 gaiaSetPointXYM (rect->Coords, 3, minx, maxy, 0.0);
4552 gaiaSetPointXYM (rect->Coords, 4, minx, miny, 0.0);
4553 }
4554 else if (geom->DimensionModel == GAIA_XY_Z_M)
4555 {
4556 gaiaSetPointXYZM (rect->Coords, 0, minx, miny, 0.0, 0.0);
4557 gaiaSetPointXYZM (rect->Coords, 1, maxx, miny, 0.0, 0.0);
4558 gaiaSetPointXYZM (rect->Coords, 2, maxx, maxy, 0.0, 0.0);
4559 gaiaSetPointXYZM (rect->Coords, 3, minx, maxy, 0.0, 0.0);
4560 gaiaSetPointXYZM (rect->Coords, 4, minx, miny, 0.0, 0.0);
4561 }
4562 else
4563 {
4564 gaiaSetPoint (rect->Coords, 0, minx, miny); /* vertex # 1 */
4565 gaiaSetPoint (rect->Coords, 1, maxx, miny); /* vertex # 2 */
4566 gaiaSetPoint (rect->Coords, 2, maxx, maxy); /* vertex # 3 */
4567 gaiaSetPoint (rect->Coords, 3, minx, maxy); /* vertex # 4 */
4568 gaiaSetPoint (rect->Coords, 4, minx, miny); /* vertex # 5 [same as vertex # 1 to close the polygon] */
4569 }
4570
4571 return bbox;
4572 }
4573
4574 static int
voronoj_mbr_contains(gaiaGeomCollPtr g1,gaiaGeomCollPtr g2)4575 voronoj_mbr_contains (gaiaGeomCollPtr g1, gaiaGeomCollPtr g2)
4576 {
4577 /* checks if MBR#1 fully contains MBR#2 */
4578 if (g2->MinX < g1->MinX)
4579 return 0;
4580 if (g2->MaxX > g1->MaxX)
4581 return 0;
4582 if (g2->MinY < g1->MinY)
4583 return 0;
4584 if (g2->MaxY > g1->MaxY)
4585 return 0;
4586 return 1;
4587 }
4588
4589 static int
voronoj_mbr_overlaps(gaiaGeomCollPtr g1,gaiaGeomCollPtr g2)4590 voronoj_mbr_overlaps (gaiaGeomCollPtr g1, gaiaGeomCollPtr g2)
4591 {
4592 /* checks if two MBRs do overlap */
4593 if (g1->MaxX < g2->MinX)
4594 return 0;
4595 if (g1->MinX > g2->MaxX)
4596 return 0;
4597 if (g1->MaxY < g2->MinY)
4598 return 0;
4599 if (g1->MinY > g2->MaxY)
4600 return 0;
4601 return 1;
4602 }
4603
4604 static gaiaGeomCollPtr
voronoj_postprocess(struct splite_internal_cache * cache,gaiaGeomCollPtr raw,gaiaGeomCollPtr envelope,int only_edges)4605 voronoj_postprocess (struct splite_internal_cache *cache, gaiaGeomCollPtr raw,
4606 gaiaGeomCollPtr envelope, int only_edges)
4607 {
4608 /* postprocessing the result returned by GEOS Voronoj */
4609 gaiaGeomCollPtr candidate;
4610 gaiaGeomCollPtr result;
4611 gaiaGeomCollPtr framed;
4612 gaiaPolygonPtr pg;
4613 gaiaPolygonPtr new_pg;
4614
4615 if (raw->DimensionModel == GAIA_XY_Z)
4616 result = gaiaAllocGeomCollXYZ ();
4617 else if (raw->DimensionModel == GAIA_XY_M)
4618 result = gaiaAllocGeomCollXYM ();
4619 else if (raw->DimensionModel == GAIA_XY_Z_M)
4620 result = gaiaAllocGeomCollXYZM ();
4621 else
4622 result = gaiaAllocGeomColl ();
4623 result->Srid = raw->Srid;
4624 result->DeclaredType = GAIA_MULTIPOLYGON;
4625
4626 if (raw->DimensionModel == GAIA_XY_Z)
4627 candidate = gaiaAllocGeomCollXYZ ();
4628 else if (raw->DimensionModel == GAIA_XY_M)
4629 candidate = gaiaAllocGeomCollXYM ();
4630 else if (raw->DimensionModel == GAIA_XY_Z_M)
4631 candidate = gaiaAllocGeomCollXYZM ();
4632 else
4633 candidate = gaiaAllocGeomColl ();
4634 candidate->Srid = raw->Srid;
4635 candidate->DeclaredType = GAIA_POLYGON;
4636
4637 gaiaMbrGeometry (raw);
4638 gaiaMbrGeometry (envelope);
4639 pg = raw->FirstPolygon;
4640 while (pg != NULL)
4641 {
4642 candidate->FirstPolygon = pg;
4643 candidate->LastPolygon = pg;
4644 candidate->MinX = pg->MinX;
4645 candidate->MinY = pg->MinY;
4646 candidate->MaxX = pg->MaxX;
4647 candidate->MaxY = pg->MaxY;
4648 if (voronoj_mbr_contains (envelope, candidate))
4649 {
4650 /* copying a Polygon fully contained within the frame */
4651 new_pg = gaiaClonePolygon (pg);
4652 if (result->FirstPolygon == NULL)
4653 result->FirstPolygon = new_pg;
4654 if (result->LastPolygon != NULL)
4655 result->LastPolygon->Next = new_pg;
4656 result->LastPolygon = new_pg;
4657 }
4658 else if (voronoj_mbr_overlaps (envelope, candidate))
4659 {
4660 /* found a polygon only partially contained within the frame */
4661 new_pg = gaiaClonePolygon (pg);
4662 candidate->FirstPolygon = new_pg;
4663 candidate->LastPolygon = new_pg;
4664 if (cache == NULL)
4665 framed = gaiaGeometryIntersection (envelope, candidate);
4666 else
4667 framed =
4668 gaiaGeometryIntersection_r (cache, envelope, candidate);
4669 candidate->FirstPolygon = NULL;
4670 candidate->LastPolygon = NULL;
4671 gaiaFreePolygon (new_pg);
4672 if (framed)
4673 {
4674 /* copying all framed polygons into the result */
4675 gaiaPolygonPtr pg2 = framed->FirstPolygon;
4676 while (pg2 != NULL)
4677 {
4678 if (result->FirstPolygon == NULL)
4679 result->FirstPolygon = pg2;
4680 if (result->LastPolygon != NULL)
4681 result->LastPolygon->Next = pg2;
4682 result->LastPolygon = pg2;
4683 pg2 = pg2->Next;
4684 }
4685 framed->FirstPolygon = NULL;
4686 framed->LastPolygon = NULL;
4687 gaiaFreeGeomColl (framed);
4688 }
4689 }
4690 pg = pg->Next;
4691 }
4692
4693 candidate->FirstPolygon = NULL;
4694 candidate->LastPolygon = NULL;
4695 gaiaFreeGeomColl (candidate);
4696 gaiaFreeGeomColl (raw);
4697 if (only_edges)
4698 {
4699 gaiaGeomCollPtr lines = gaiaLinearize (result, 1);
4700 gaiaFreeGeomColl (result);
4701 return lines;
4702 }
4703 return result;
4704 }
4705 #endif
4706
4707 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaVoronojDiagram(gaiaGeomCollPtr geom,double extra_frame_size,double tolerance,int only_edges)4708 gaiaVoronojDiagram (gaiaGeomCollPtr geom, double extra_frame_size,
4709 double tolerance, int only_edges)
4710 {
4711 /* Voronoj Diagram */
4712 gaiaGeomCollPtr result = NULL;
4713 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
4714 GEOSGeometry *g1;
4715 GEOSGeometry *g2;
4716 #ifdef GEOS_REENTRANT
4717 GEOSGeometry *env;
4718 gaiaGeomCollPtr envelope;
4719 #else
4720 void *voronoj;
4721 gaiaPolygonPtr pg;
4722 int pgs = 0;
4723 int errs = 0;
4724 #endif
4725 gaiaResetGeosMsg ();
4726 if (!geom)
4727 return NULL;
4728 g1 = gaiaToGeos (geom);
4729 #ifdef GEOS_REENTRANT /* GEOS >= 3.5.0 directly supports Voronoj */
4730 envelope = voronoj_envelope (geom, extra_frame_size);
4731 env = gaiaToGeos (envelope);
4732 g2 = GEOSVoronoiDiagram (g1, env, tolerance, 0);
4733 GEOSGeom_destroy (g1);
4734 GEOSGeom_destroy (env);
4735 if (!g2)
4736 {
4737 gaiaFreeGeomColl (envelope);
4738 return NULL;
4739 }
4740 if (geom->DimensionModel == GAIA_XY_Z)
4741 result = gaiaFromGeos_XYZ (g2);
4742 else if (geom->DimensionModel == GAIA_XY_M)
4743 result = gaiaFromGeos_XYM (g2);
4744 else if (geom->DimensionModel == GAIA_XY_Z_M)
4745 result = gaiaFromGeos_XYZM (g2);
4746 else
4747 result = gaiaFromGeos_XY (g2);
4748 GEOSGeom_destroy (g2);
4749 result = voronoj_postprocess (NULL, result, envelope, only_edges);
4750 gaiaFreeGeomColl (envelope);
4751 if (result == NULL)
4752 return NULL;
4753 #else
4754 g2 = GEOSDelaunayTriangulation (g1, tolerance, 0);
4755 GEOSGeom_destroy (g1);
4756 if (!g2)
4757 return NULL;
4758 if (geom->DimensionModel == GAIA_XY_Z)
4759 result = gaiaFromGeos_XYZ (g2);
4760 else if (geom->DimensionModel == GAIA_XY_M)
4761 result = gaiaFromGeos_XYM (g2);
4762 else if (geom->DimensionModel == GAIA_XY_Z_M)
4763 result = gaiaFromGeos_XYZM (g2);
4764 else
4765 result = gaiaFromGeos_XY (g2);
4766 GEOSGeom_destroy (g2);
4767 if (result == NULL)
4768 return NULL;
4769 pg = result->FirstPolygon;
4770 while (pg)
4771 {
4772 /* counting how many triangles are in Delaunay */
4773 if (delaunay_triangle_check (pg))
4774 pgs++;
4775 else
4776 errs++;
4777 pg = pg->Next;
4778 }
4779 if (pgs == 0 || errs)
4780 {
4781 gaiaFreeGeomColl (result);
4782 return NULL;
4783 }
4784
4785 /* building the Voronoj Diagram from Delaunay */
4786 voronoj = voronoj_build (pgs, result->FirstPolygon, extra_frame_size);
4787 gaiaFreeGeomColl (result);
4788
4789 /* creating the Geometry representing Voronoj */
4790 if (geom->DimensionModel == GAIA_XY_Z)
4791 result = gaiaAllocGeomCollXYZ ();
4792 else if (geom->DimensionModel == GAIA_XY_M)
4793 result = gaiaAllocGeomCollXYM ();
4794 else if (geom->DimensionModel == GAIA_XY_Z_M)
4795 result = gaiaAllocGeomCollXYZM ();
4796 else
4797 result = gaiaAllocGeomColl ();
4798 result = voronoj_export (voronoj, result, only_edges);
4799 voronoj_free (voronoj);
4800
4801 result->Srid = geom->Srid;
4802 if (only_edges)
4803 result->DeclaredType = GAIA_MULTILINESTRING;
4804 else
4805 result->DeclaredType = GAIA_MULTIPOLYGON;
4806 #endif /* end GEOS_REENTRANT */
4807 #else
4808 if (geom == NULL || extra_frame_size == 0.0 || tolerance == 0.0
4809 || only_edges == 0)
4810 geom = NULL; /* silencing stupid compiler warnings */
4811 #endif /* end GEOS_USE_ONLY_R_API */
4812 return result;
4813 }
4814
4815 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaVoronojDiagram_r(const void * p_cache,gaiaGeomCollPtr geom,double extra_frame_size,double tolerance,int only_edges)4816 gaiaVoronojDiagram_r (const void *p_cache, gaiaGeomCollPtr geom,
4817 double extra_frame_size, double tolerance, int only_edges)
4818 {
4819 /* Voronoj Diagram */
4820 GEOSGeometry *g1;
4821 GEOSGeometry *g2;
4822 gaiaGeomCollPtr result;
4823 #ifdef GEOS_REENTRANT
4824 GEOSGeometry *env;
4825 gaiaGeomCollPtr envelope;
4826 #else
4827 gaiaPolygonPtr pg;
4828 int pgs = 0;
4829 int errs = 0;
4830 void *voronoj;
4831 #endif
4832 struct splite_internal_cache *cache =
4833 (struct splite_internal_cache *) p_cache;
4834 GEOSContextHandle_t handle = NULL;
4835 if (cache == NULL)
4836 return NULL;
4837 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
4838 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
4839 return NULL;
4840 handle = cache->GEOS_handle;
4841 if (handle == NULL)
4842 return NULL;
4843 gaiaResetGeosMsg_r (cache);
4844 if (!geom)
4845 return NULL;
4846 g1 = gaiaToGeos_r (cache, geom);
4847 #ifdef GEOS_REENTRANT /* GEOS >= 3.5.0 directly supports Voronoj */
4848 envelope = voronoj_envelope (geom, extra_frame_size);
4849 env = gaiaToGeos_r (cache, envelope);
4850 g2 = GEOSVoronoiDiagram_r (handle, g1, env, tolerance, 0);
4851 GEOSGeom_destroy_r (handle, g1);
4852 GEOSGeom_destroy_r (handle, env);
4853 if (!g2)
4854 {
4855 gaiaFreeGeomColl (envelope);
4856 return NULL;
4857 }
4858 if (geom->DimensionModel == GAIA_XY_Z)
4859 result = gaiaFromGeos_XYZ_r (cache, g2);
4860 else if (geom->DimensionModel == GAIA_XY_M)
4861 result = gaiaFromGeos_XYM_r (cache, g2);
4862 else if (geom->DimensionModel == GAIA_XY_Z_M)
4863 result = gaiaFromGeos_XYZM_r (cache, g2);
4864 else
4865 result = gaiaFromGeos_XY_r (cache, g2);
4866 GEOSGeom_destroy_r (handle, g2);
4867 result = voronoj_postprocess (cache, result, envelope, only_edges);
4868 gaiaFreeGeomColl (envelope);
4869 if (result == NULL)
4870 return NULL;
4871 #else
4872 g2 = GEOSDelaunayTriangulation_r (handle, g1, tolerance, 0);
4873 GEOSGeom_destroy_r (handle, g1);
4874 if (!g2)
4875 return NULL;
4876 if (geom->DimensionModel == GAIA_XY_Z)
4877 result = gaiaFromGeos_XYZ_r (cache, g2);
4878 else if (geom->DimensionModel == GAIA_XY_M)
4879 result = gaiaFromGeos_XYM_r (cache, g2);
4880 else if (geom->DimensionModel == GAIA_XY_Z_M)
4881 result = gaiaFromGeos_XYZM_r (cache, g2);
4882 else
4883 result = gaiaFromGeos_XY_r (cache, g2);
4884 GEOSGeom_destroy_r (handle, g2);
4885 if (result == NULL)
4886 return NULL;
4887 pg = result->FirstPolygon;
4888 while (pg)
4889 {
4890 /* counting how many triangles are in Delaunay */
4891 if (delaunay_triangle_check (pg))
4892 pgs++;
4893 else
4894 errs++;
4895 pg = pg->Next;
4896 }
4897 if (pgs == 0 || errs)
4898 {
4899 gaiaFreeGeomColl (result);
4900 return NULL;
4901 }
4902
4903 /* building the Voronoj Diagram from Delaunay */
4904 voronoj =
4905 voronoj_build_r (cache, pgs, result->FirstPolygon, extra_frame_size);
4906 gaiaFreeGeomColl (result);
4907
4908 /* creating the Geometry representing Voronoj */
4909 if (geom->DimensionModel == GAIA_XY_Z)
4910 result = gaiaAllocGeomCollXYZ ();
4911 else if (geom->DimensionModel == GAIA_XY_M)
4912 result = gaiaAllocGeomCollXYM ();
4913 else if (geom->DimensionModel == GAIA_XY_Z_M)
4914 result = gaiaAllocGeomCollXYZM ();
4915 else
4916 result = gaiaAllocGeomColl ();
4917 result = voronoj_export_r (cache, voronoj, result, only_edges);
4918 voronoj_free (voronoj);
4919
4920 result->Srid = geom->Srid;
4921 if (only_edges)
4922 result->DeclaredType = GAIA_MULTILINESTRING;
4923 else
4924 result->DeclaredType = GAIA_MULTIPOLYGON;
4925 #endif /* end GEOS_REENTRANT */
4926 return result;
4927 }
4928
4929 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaConcaveHull(gaiaGeomCollPtr geom,double factor,double tolerance,int allow_holes)4930 gaiaConcaveHull (gaiaGeomCollPtr geom, double factor, double tolerance,
4931 int allow_holes)
4932 {
4933 /* Concave Hull */
4934 gaiaGeomCollPtr result = NULL;
4935 #ifndef GEOS_USE_ONLY_R_API /* obsolete versions non fully thread-safe */
4936 GEOSGeometry *g1;
4937 GEOSGeometry *g2;
4938 gaiaGeomCollPtr concave_hull;
4939 gaiaPolygonPtr pg;
4940 int pgs = 0;
4941 int errs = 0;
4942 gaiaResetGeosMsg ();
4943 if (!geom)
4944 return NULL;
4945 g1 = gaiaToGeos (geom);
4946 g2 = GEOSDelaunayTriangulation (g1, tolerance, 0);
4947 GEOSGeom_destroy (g1);
4948 if (!g2)
4949 return NULL;
4950 if (geom->DimensionModel == GAIA_XY_Z)
4951 result = gaiaFromGeos_XYZ (g2);
4952 else if (geom->DimensionModel == GAIA_XY_M)
4953 result = gaiaFromGeos_XYM (g2);
4954 else if (geom->DimensionModel == GAIA_XY_Z_M)
4955 result = gaiaFromGeos_XYZM (g2);
4956 else
4957 result = gaiaFromGeos_XY (g2);
4958 GEOSGeom_destroy (g2);
4959 if (result == NULL)
4960 return NULL;
4961 pg = result->FirstPolygon;
4962 while (pg)
4963 {
4964 /* counting how many triangles are in Delaunay */
4965 if (delaunay_triangle_check (pg))
4966 pgs++;
4967 else
4968 errs++;
4969 pg = pg->Next;
4970 }
4971 if (pgs == 0 || errs)
4972 {
4973 gaiaFreeGeomColl (result);
4974 return NULL;
4975 }
4976
4977 /* building the Concave Hull from Delaunay */
4978 concave_hull =
4979 concave_hull_build (result->FirstPolygon, geom->DimensionModel, factor,
4980 allow_holes);
4981 gaiaFreeGeomColl (result);
4982 if (!concave_hull)
4983 return NULL;
4984 result = concave_hull;
4985
4986 result->Srid = geom->Srid;
4987 #else
4988 if (geom == NULL || factor == 0.0 || tolerance == 0.0 || allow_holes == 0)
4989 geom = NULL; /* silencing stupid compiler warnings */
4990 #endif
4991 return result;
4992 }
4993
4994 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaConcaveHull_r(const void * p_cache,gaiaGeomCollPtr geom,double factor,double tolerance,int allow_holes)4995 gaiaConcaveHull_r (const void *p_cache, gaiaGeomCollPtr geom, double factor,
4996 double tolerance, int allow_holes)
4997 {
4998 /* Concave Hull */
4999 GEOSGeometry *g1;
5000 GEOSGeometry *g2;
5001 gaiaGeomCollPtr result;
5002 gaiaGeomCollPtr concave_hull;
5003 gaiaPolygonPtr pg;
5004 int pgs = 0;
5005 int errs = 0;
5006 struct splite_internal_cache *cache =
5007 (struct splite_internal_cache *) p_cache;
5008 GEOSContextHandle_t handle = NULL;
5009 if (cache == NULL)
5010 return NULL;
5011 if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
5012 || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
5013 return NULL;
5014 handle = cache->GEOS_handle;
5015 if (handle == NULL)
5016 return NULL;
5017 gaiaResetGeosMsg_r (cache);
5018 if (!geom)
5019 return NULL;
5020 g1 = gaiaToGeos_r (cache, geom);
5021 g2 = GEOSDelaunayTriangulation_r (handle, g1, tolerance, 0);
5022 GEOSGeom_destroy_r (handle, g1);
5023 if (!g2)
5024 return NULL;
5025 if (geom->DimensionModel == GAIA_XY_Z)
5026 result = gaiaFromGeos_XYZ_r (cache, g2);
5027 else if (geom->DimensionModel == GAIA_XY_M)
5028 result = gaiaFromGeos_XYM_r (cache, g2);
5029 else if (geom->DimensionModel == GAIA_XY_Z_M)
5030 result = gaiaFromGeos_XYZM_r (cache, g2);
5031 else
5032 result = gaiaFromGeos_XY_r (cache, g2);
5033 GEOSGeom_destroy_r (handle, g2);
5034 if (result == NULL)
5035 return NULL;
5036 pg = result->FirstPolygon;
5037 while (pg)
5038 {
5039 /* counting how many triangles are in Delaunay */
5040 if (delaunay_triangle_check (pg))
5041 pgs++;
5042 else
5043 errs++;
5044 pg = pg->Next;
5045 }
5046 if (pgs == 0 || errs)
5047 {
5048 gaiaFreeGeomColl (result);
5049 return NULL;
5050 }
5051
5052 /* building the Concave Hull from Delaunay */
5053 concave_hull =
5054 concave_hull_build_r (p_cache, result->FirstPolygon,
5055 geom->DimensionModel, factor, allow_holes);
5056 gaiaFreeGeomColl (result);
5057 if (!concave_hull)
5058 return NULL;
5059 result = concave_hull;
5060
5061 result->Srid = geom->Srid;
5062 return result;
5063 }
5064
5065 #endif /* end GEOS advanced features */
5066
5067 #endif /* end including GEOS */
5068