1 /*
2
3 gg_advanced.c -- Gaia advanced geometric operations
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 <math.h>
50 #include <float.h>
51 #include <string.h>
52
53 #if defined(_WIN32) && !defined(__MINGW32__)
54 #include "config-msvc.h"
55 #else
56 #include "config.h"
57 #endif
58
59 #include <spatialite/sqlite.h>
60
61 #include <spatialite.h>
62 #include <spatialite_private.h>
63 #include <spatialite/gaiageo.h>
64 #include <spatialite/debug.h>
65
66 GAIAGEO_DECLARE double
gaiaMeasureLength(int dims,double * coords,int vert)67 gaiaMeasureLength (int dims, double *coords, int vert)
68 {
69 /* computes the total length */
70 double lung = 0.0;
71 double xx1;
72 double xx2;
73 double yy1;
74 double yy2;
75 double x;
76 double y;
77 double z;
78 double m;
79 double dist;
80 int ind;
81 if (vert <= 0)
82 return lung;
83 if (dims == GAIA_XY_Z)
84 {
85 gaiaGetPointXYZ (coords, 0, &xx1, &yy1, &z);
86 }
87 else if (dims == GAIA_XY_M)
88 {
89 gaiaGetPointXYM (coords, 0, &xx1, &yy1, &m);
90 }
91 else if (dims == GAIA_XY_Z_M)
92 {
93 gaiaGetPointXYZM (coords, 0, &xx1, &yy1, &z, &m);
94 }
95 else
96 {
97 gaiaGetPoint (coords, 0, &xx1, &yy1);
98 }
99 for (ind = 1; ind < vert; ind++)
100 {
101 if (dims == GAIA_XY_Z)
102 {
103 gaiaGetPointXYZ (coords, ind, &xx2, &yy2, &z);
104 }
105 else if (dims == GAIA_XY_M)
106 {
107 gaiaGetPointXYM (coords, ind, &xx2, &yy2, &m);
108 }
109 else if (dims == GAIA_XY_Z_M)
110 {
111 gaiaGetPointXYZM (coords, ind, &xx2, &yy2, &z, &m);
112 }
113 else
114 {
115 gaiaGetPoint (coords, ind, &xx2, &yy2);
116 }
117 x = xx1 - xx2;
118 y = yy1 - yy2;
119 dist = sqrt ((x * x) + (y * y));
120 lung += dist;
121 xx1 = xx2;
122 yy1 = yy2;
123 }
124 return lung;
125 }
126
127 GAIAGEO_DECLARE double
gaiaMeasureArea(gaiaRingPtr ring)128 gaiaMeasureArea (gaiaRingPtr ring)
129 {
130 /* computes the area */
131 int iv;
132 double xx;
133 double yy;
134 double x;
135 double y;
136 double z;
137 double m;
138 double area = 0.0;
139 if (!ring)
140 return 0.0;
141 if (ring->DimensionModel == GAIA_XY_Z)
142 {
143 gaiaGetPointXYZ (ring->Coords, 0, &xx, &yy, &z);
144 }
145 else if (ring->DimensionModel == GAIA_XY_M)
146 {
147 gaiaGetPointXYM (ring->Coords, 0, &xx, &yy, &m);
148 }
149 else if (ring->DimensionModel == GAIA_XY_Z_M)
150 {
151 gaiaGetPointXYZM (ring->Coords, 0, &xx, &yy, &z, &m);
152 }
153 else
154 {
155 gaiaGetPoint (ring->Coords, 0, &xx, &yy);
156 }
157 for (iv = 1; iv < ring->Points; iv++)
158 {
159 if (ring->DimensionModel == GAIA_XY_Z)
160 {
161 gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
162 }
163 else if (ring->DimensionModel == GAIA_XY_M)
164 {
165 gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
166 }
167 else if (ring->DimensionModel == GAIA_XY_Z_M)
168 {
169 gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
170 }
171 else
172 {
173 gaiaGetPoint (ring->Coords, iv, &x, &y);
174 }
175 area += ((xx * y) - (x * yy));
176 xx = x;
177 yy = y;
178 }
179 area /= 2.0;
180 return fabs (area);
181 }
182
183 GAIAGEO_DECLARE void
gaiaRingCentroid(gaiaRingPtr ring,double * rx,double * ry)184 gaiaRingCentroid (gaiaRingPtr ring, double *rx, double *ry)
185 {
186 /* computes the simple ring centroid */
187 double cx = 0.0;
188 double cy = 0.0;
189 double xx;
190 double yy;
191 double x;
192 double y;
193 double z;
194 double m;
195 double coeff;
196 double area;
197 double term;
198 int iv;
199 if (!ring)
200 {
201 *rx = -DBL_MAX;
202 *ry = -DBL_MAX;
203 return;
204 }
205 area = gaiaMeasureArea (ring);
206 coeff = 1.0 / (area * 6.0);
207 if (ring->DimensionModel == GAIA_XY_Z)
208 {
209 gaiaGetPointXYZ (ring->Coords, 0, &xx, &yy, &z);
210 }
211 else if (ring->DimensionModel == GAIA_XY_M)
212 {
213 gaiaGetPointXYM (ring->Coords, 0, &xx, &yy, &m);
214 }
215 else if (ring->DimensionModel == GAIA_XY_Z_M)
216 {
217 gaiaGetPointXYZM (ring->Coords, 0, &xx, &yy, &z, &m);
218 }
219 else
220 {
221 gaiaGetPoint (ring->Coords, 0, &xx, &yy);
222 }
223 for (iv = 1; iv < ring->Points; iv++)
224 {
225 if (ring->DimensionModel == GAIA_XY_Z)
226 {
227 gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
228 }
229 else if (ring->DimensionModel == GAIA_XY_M)
230 {
231 gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
232 }
233 else if (ring->DimensionModel == GAIA_XY_Z_M)
234 {
235 gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
236 }
237 else
238 {
239 gaiaGetPoint (ring->Coords, iv, &x, &y);
240 }
241 term = (xx * y) - (x * yy);
242 cx += (xx + x) * term;
243 cy += (yy + y) * term;
244 xx = x;
245 yy = y;
246 }
247 *rx = fabs (cx * coeff);
248 *ry = fabs (cy * coeff);
249 }
250
251 GAIAGEO_DECLARE void
gaiaClockwise(gaiaRingPtr p)252 gaiaClockwise (gaiaRingPtr p)
253 {
254 /* determines clockwise or anticlockwise direction */
255 int ind;
256 int ix;
257 double xx;
258 double yy;
259 double x;
260 double y;
261 double z;
262 double m;
263 double area = 0.0;
264 for (ind = 0; ind < p->Points; ind++)
265 {
266 if (p->DimensionModel == GAIA_XY_Z)
267 {
268 gaiaGetPointXYZ (p->Coords, ind, &xx, &yy, &z);
269 }
270 else if (p->DimensionModel == GAIA_XY_M)
271 {
272 gaiaGetPointXYM (p->Coords, ind, &xx, &yy, &m);
273 }
274 else if (p->DimensionModel == GAIA_XY_Z_M)
275 {
276 gaiaGetPointXYZM (p->Coords, ind, &xx, &yy, &z, &m);
277 }
278 else
279 {
280 gaiaGetPoint (p->Coords, ind, &xx, &yy);
281 }
282 ix = (ind + 1) % p->Points;
283 if (p->DimensionModel == GAIA_XY_Z)
284 {
285 gaiaGetPointXYZ (p->Coords, ix, &x, &y, &z);
286 }
287 else if (p->DimensionModel == GAIA_XY_M)
288 {
289 gaiaGetPointXYM (p->Coords, ix, &x, &y, &m);
290 }
291 else if (p->DimensionModel == GAIA_XY_Z_M)
292 {
293 gaiaGetPointXYZM (p->Coords, ix, &x, &y, &z, &m);
294 }
295 else
296 {
297 gaiaGetPoint (p->Coords, ix, &x, &y);
298 }
299 area += ((xx * y) - (x * yy));
300 }
301 area /= 2.0;
302 if (area >= 0.0)
303 p->Clockwise = 0;
304 else
305 p->Clockwise = 1;
306 }
307
308 GAIAGEO_DECLARE double
gaiaCurvosityIndex(const void * p_cache,gaiaLinestringPtr ln,int extra_points)309 gaiaCurvosityIndex (const void *p_cache, gaiaLinestringPtr ln, int extra_points)
310 {
311 /* calculates the Curvosity Index of some Linestring */
312 #ifndef OMIT_GEOS /* including GEOS */
313 double x;
314 double y;
315 double z;
316 double m;
317 int iv;
318 int ov;
319 gaiaGeomCollPtr geo = NULL;
320 gaiaLinestringPtr refln;
321 double refline_length;
322 double line_length =
323 gaiaMeasureLength (ln->DimensionModel, ln->Coords, ln->Points);
324
325 /* building the Reference Line */
326 refln = gaiaAllocLinestring (2 + extra_points);
327 /* inserting the first vertex of Line into the Reference Line */
328 iv = 0;
329 ov = 0;
330 if (ln->DimensionModel == GAIA_XY_Z)
331 {
332 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
333 }
334 else if (ln->DimensionModel == GAIA_XY_M)
335 {
336 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
337 }
338 else if (ln->DimensionModel == GAIA_XY_Z_M)
339 {
340 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
341 }
342 else
343 {
344 gaiaGetPoint (ln->Coords, iv, &x, &y);
345 }
346 gaiaSetPoint (refln->Coords, 0, x, y);
347 ov++;
348 if (extra_points > 0)
349 {
350 /* inserting all interpolated points into the Reference Line */
351 int i;
352 double perc = 1.0 / (double) (extra_points + 1);
353 double progr = perc;
354 gaiaGeomCollPtr geo;
355 gaiaGeomCollPtr result;
356 gaiaPointPtr pt;
357 if (ln->DimensionModel == GAIA_XY_Z)
358 geo = gaiaAllocGeomCollXYZ ();
359 else if (ln->DimensionModel == GAIA_XY_M)
360 geo = gaiaAllocGeomCollXYM ();
361 else if (ln->DimensionModel == GAIA_XY_Z_M)
362 geo = gaiaAllocGeomCollXYZM ();
363 else
364 geo = gaiaAllocGeomColl ();
365 gaiaInsertLinestringInGeomColl (geo, ln);
366 for (i = 0; i < extra_points; i++)
367 {
368 if (p_cache != NULL)
369 result = gaiaLineInterpolatePoint_r (p_cache, geo, progr);
370 else
371 result = gaiaLineInterpolatePoint (geo, progr);
372 if (result == NULL)
373 goto error;
374 pt = result->FirstPoint;
375 if (pt == NULL)
376 goto error;
377 x = pt->X;
378 y = pt->Y;
379 gaiaFreeGeomColl (result);
380 gaiaSetPoint (refln->Coords, ov, x, y);
381 ov++;
382 progr += perc;
383 }
384 geo->FirstLinestring = NULL; /* releasing ownership on Line */
385 geo->LastLinestring = NULL;
386 gaiaFreeGeomColl (geo);
387 }
388 /* inserting the last vertex of Line into the Reference Line */
389 iv = ln->Points - 1;
390 if (ln->DimensionModel == GAIA_XY_Z)
391 {
392 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
393 }
394 else if (ln->DimensionModel == GAIA_XY_M)
395 {
396 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
397 }
398 else if (ln->DimensionModel == GAIA_XY_Z_M)
399 {
400 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
401 }
402 else
403 {
404 gaiaGetPoint (ln->Coords, iv, &x, &y);
405 }
406 gaiaSetPoint (refln->Coords, ov, x, y);
407 refline_length =
408 gaiaMeasureLength (refln->DimensionModel, refln->Coords, refln->Points);
409 return (refline_length / line_length);
410
411 error:
412 if (geo != NULL)
413 {
414 geo->FirstLinestring = NULL; /* releasing ownership on Line */
415 geo->LastLinestring = NULL;
416 gaiaFreeGeomColl (geo);
417 }
418 #endif /* end including GEOS */
419 return -1.0;
420 }
421
422 GAIAGEO_DECLARE void
gaiaUpDownHeight(gaiaLinestringPtr ln,double * up,double * down)423 gaiaUpDownHeight (gaiaLinestringPtr ln, double *up, double *down)
424 {
425 /* calculates the Uphill and Downhill Height of some Linestring */
426 double x;
427 double y;
428 double z;
429 double m;
430 int iv;
431 double prev_z;
432 double tot_up = 0.0;
433 double tot_down = 0.0;
434
435 /* inserting the last vertex of Line into the Reference Line */
436 if (ln->DimensionModel == GAIA_XY || ln->DimensionModel == GAIA_XY_M)
437 {
438 *up = 0.0;
439 *down = 0.0;
440 goto stop;
441 }
442 for (iv = 0; iv < ln->Points; iv++)
443 {
444 if (ln->DimensionModel == GAIA_XY_Z)
445 {
446 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
447 }
448 else if (ln->DimensionModel == GAIA_XY_Z_M)
449 {
450 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
451 }
452 if (iv == 0)
453 {
454 prev_z = z;
455 continue;
456 }
457 if (z > prev_z)
458 tot_up += z - prev_z;
459 else
460 tot_down += prev_z - z;
461 prev_z = z;
462 }
463 stop:
464 *up = tot_up;
465 *down = tot_down;
466 }
467
468 GAIAGEO_DECLARE int
gaiaIsPointOnRingSurface(gaiaRingPtr ring,double pt_x,double pt_y)469 gaiaIsPointOnRingSurface (gaiaRingPtr ring, double pt_x, double pt_y)
470 {
471 /* tests if a POINT falls inside a RING */
472 int isInternal = 0;
473 int cnt;
474 int i;
475 int j;
476 double x;
477 double y;
478 double z;
479 double m;
480 double *vert_x;
481 double *vert_y;
482 double minx = DBL_MAX;
483 double miny = DBL_MAX;
484 double maxx = -DBL_MAX;
485 double maxy = -DBL_MAX;
486 cnt = ring->Points;
487 cnt--; /* ignoring last vertex because surely identical to the first one */
488 if (cnt < 2)
489 return 0;
490 /* allocating and loading an array of vertices */
491 vert_x = malloc (sizeof (double) * (cnt));
492 vert_y = malloc (sizeof (double) * (cnt));
493 for (i = 0; i < cnt; i++)
494 {
495 if (ring->DimensionModel == GAIA_XY_Z)
496 {
497 gaiaGetPointXYZ (ring->Coords, i, &x, &y, &z);
498 }
499 else if (ring->DimensionModel == GAIA_XY_M)
500 {
501 gaiaGetPointXYM (ring->Coords, i, &x, &y, &m);
502 }
503 else if (ring->DimensionModel == GAIA_XY_Z_M)
504 {
505 gaiaGetPointXYZM (ring->Coords, i, &x, &y, &z, &m);
506 }
507 else
508 {
509 gaiaGetPoint (ring->Coords, i, &x, &y);
510 }
511 vert_x[i] = x;
512 vert_y[i] = y;
513 if (x < minx)
514 minx = x;
515 if (x > maxx)
516 maxx = x;
517 if (y < miny)
518 miny = y;
519 if (y > maxy)
520 maxy = y;
521 }
522 if (pt_x < minx || pt_x > maxx)
523 goto end; /* outside the bounding box (x axis) */
524 if (pt_y < miny || pt_y > maxy)
525 goto end; /* outside the bounding box (y axis) */
526 for (i = 0, j = cnt - 1; i < cnt; j = i++)
527 {
528 /* The definitive reference is "Point in Polyon Strategies" by
529 / Eric Haines [Gems IV] pp. 24-46.
530 / The code in the Sedgewick book Algorithms (2nd Edition, p.354) is
531 / incorrect.
532 */
533 if ((((vert_y[i] <= pt_y) && (pt_y < vert_y[j]))
534 || ((vert_y[j] <= pt_y) && (pt_y < vert_y[i])))
535 && (pt_x <
536 (vert_x[j] - vert_x[i]) * (pt_y - vert_y[i]) / (vert_y[j] -
537 vert_y[i]) +
538 vert_x[i]))
539 isInternal = !isInternal;
540 }
541 end:
542 free (vert_x);
543 free (vert_y);
544 return isInternal;
545 }
546
547 GAIAGEO_DECLARE double
gaiaMinDistance(double x0,double y0,int dims,double * coords,int n_vert)548 gaiaMinDistance (double x0, double y0, int dims, double *coords, int n_vert)
549 {
550 /* computing minimal distance between a POINT and a linestring/ring */
551 double x;
552 double y;
553 double z;
554 double m;
555 double ox;
556 double oy;
557 double lineMag;
558 double u;
559 double px;
560 double py;
561 double dist;
562 double min_dist = DBL_MAX;
563 int iv;
564 if (n_vert < 2)
565 return min_dist; /* not a valid linestring */
566 /* computing distance from first vertex */
567 ox = *(coords + 0);
568 oy = *(coords + 1);
569 min_dist = sqrt (((x0 - ox) * (x0 - ox)) + ((y0 - oy) * (y0 - oy)));
570 for (iv = 1; iv < n_vert; iv++)
571 {
572 /* segment start-end coordinates */
573 if (dims == GAIA_XY_Z)
574 {
575 gaiaGetPointXYZ (coords, iv - 1, &ox, &oy, &z);
576 gaiaGetPointXYZ (coords, iv, &x, &y, &z);
577 }
578 else if (dims == GAIA_XY_M)
579 {
580 gaiaGetPointXYM (coords, iv - 1, &ox, &oy, &m);
581 gaiaGetPointXYM (coords, iv, &x, &y, &m);
582 }
583 else if (dims == GAIA_XY_Z_M)
584 {
585 gaiaGetPointXYZM (coords, iv - 1, &ox, &oy, &z, &m);
586 gaiaGetPointXYZM (coords, iv, &x, &y, &z, &m);
587 }
588 else
589 {
590 gaiaGetPoint (coords, iv - 1, &ox, &oy);
591 gaiaGetPoint (coords, iv, &x, &y);
592 }
593 /* computing distance from vertex */
594 dist = sqrt (((x0 - x) * (x0 - x)) + ((y0 - y) * (y0 - y)));
595 if (dist < min_dist)
596 min_dist = dist;
597 /* computing a projection */
598 lineMag = ((x - ox) * (x - ox)) + ((y - oy) * (y - oy));
599 u = (((x0 - ox) * (x - ox)) + ((y0 - oy) * (y - oy))) / lineMag;
600 if (u < 0.0 || u > 1.0)
601 ; /* closest point does not fall within the line segment */
602 else
603 {
604 px = ox + u * (x - ox);
605 py = oy + u * (y - oy);
606 dist = sqrt (((x0 - px) * (x0 - px)) + ((y0 - py) * (y0 - py)));
607 if (dist < min_dist)
608 min_dist = dist;
609 }
610 }
611 return min_dist;
612 }
613
614 GAIAGEO_DECLARE int
gaiaIsPointOnPolygonSurface(gaiaPolygonPtr polyg,double x,double y)615 gaiaIsPointOnPolygonSurface (gaiaPolygonPtr polyg, double x, double y)
616 {
617 /* tests if a POINT falls inside a POLYGON */
618 int ib;
619 gaiaRingPtr ring = polyg->Exterior;
620 if (gaiaIsPointOnRingSurface (ring, x, y))
621 {
622 /* ok, the POINT falls inside the polygon */
623 for (ib = 0; ib < polyg->NumInteriors; ib++)
624 {
625 ring = polyg->Interiors + ib;
626 if (gaiaIsPointOnRingSurface (ring, x, y))
627 {
628 /* no, the POINT fall inside some hole */
629 return 0;
630 }
631 }
632 return 1;
633 }
634 return 0;
635 }
636
637 GAIAGEO_DECLARE int
gaiaIntersect(double * x0,double * y0,double x1,double y1,double x2,double y2,double x3,double y3,double x4,double y4)638 gaiaIntersect (double *x0, double *y0, double x1, double y1, double x2,
639 double y2, double x3, double y3, double x4, double y4)
640 {
641 /* computes intersection [if any] between two line segments
642 / the intersection POINT has coordinates (x0, y0)
643 / first line is identified by(x1, y1) and (x2, y2)
644 / second line is identified by (x3, y3) and (x4, y4)
645 */
646 double x;
647 double y;
648 double a1;
649 double b1;
650 double c1;
651 double a2;
652 double b2;
653 double c2;
654 double m1;
655 double m2;
656 double p;
657 double det_inv;
658 double minx1;
659 double miny1;
660 double maxx1;
661 double maxy1;
662 double minx2;
663 double miny2;
664 double maxx2;
665 double maxy2;
666 int ok1 = 0;
667 int ok2 = 0;
668 /* building line segment's MBRs */
669 if (x2 < x1)
670 {
671 minx1 = x2;
672 maxx1 = x1;
673 }
674 else
675 {
676 minx1 = x1;
677 maxx1 = x2;
678 }
679 if (y2 < y1)
680 {
681 miny1 = y2;
682 maxy1 = y1;
683 }
684 else
685 {
686 miny1 = y1;
687 maxy1 = y2;
688 }
689 if (x4 < x3)
690 {
691 minx2 = x4;
692 maxx2 = x3;
693 }
694 else
695 {
696 minx2 = x3;
697 maxx2 = x4;
698 }
699 if (y4 < y3)
700 {
701 miny2 = y4;
702 maxy2 = y3;
703 }
704 else
705 {
706 miny2 = y3;
707 maxy2 = y4;
708 }
709 /* checking MBRs first */
710 if (minx1 > maxx2)
711 return 0;
712 if (miny1 > maxy2)
713 return 0;
714 if (maxx1 < minx2)
715 return 0;
716 if (maxy1 < miny2)
717 return 0;
718 if (minx2 > maxx1)
719 return 0;
720 if (miny2 > maxy1)
721 return 0;
722 if (maxx2 < minx1)
723 return 0;
724 if (maxy2 < miny1)
725 return 0;
726 /* there is an MBRs intersection - proceeding */
727 if ((x2 - x1) != 0.0)
728 m1 = (y2 - y1) / (x2 - x1);
729 else
730 m1 = DBL_MAX;
731 if ((x4 - x3) != 0)
732 m2 = (y4 - y3) / (x4 - x3);
733 else
734 m2 = DBL_MAX;
735 if (m1 == m2) /* parallel lines */
736 return 0;
737 if (m1 == DBL_MAX)
738 c1 = y1;
739 else
740 c1 = (y1 - m1 * x1);
741 if (m2 == DBL_MAX)
742 c2 = y3;
743 else
744 c2 = (y3 - m2 * x3);
745 if (m1 == DBL_MAX)
746 {
747 x = x1;
748 p = m2 * x1;
749 y = p + c2; /* first line is vertical */
750 goto check_bbox;
751 }
752 if (m2 == DBL_MAX)
753 {
754 x = x3;
755 p = m1 * x3;
756 y = p + c1; /* second line is vertical */
757 goto check_bbox;
758 }
759 a1 = m1;
760 a2 = m2;
761 b1 = -1;
762 b2 = -1;
763 det_inv = 1 / (a1 * b2 - a2 * b1);
764 x = ((b1 * c2 - b2 * c1) * det_inv);
765 y = ((a2 * c1 - a1 * c2) * det_inv);
766 /* now checking if intersection falls within both segment boundaries */
767 check_bbox:
768 if (x >= minx1 && x <= maxx1 && y >= miny1 && y <= maxy1)
769 ok1 = 1;
770 if (x >= minx2 && x <= maxx2 && y >= miny2 && y <= maxy2)
771 ok2 = 1;
772 if (ok1 && ok2)
773 {
774 /* intersection point falls within the segments */
775 *x0 = x;
776 *y0 = y;
777 return 1;
778 }
779 return 0;
780 }
781
782 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaSanitize(gaiaGeomCollPtr geom)783 gaiaSanitize (gaiaGeomCollPtr geom)
784 {
785 /*
786 / sanitizes a GEOMETRYCOLLECTION:
787 / - repeated vertices are omitted
788 / - ring closure is enforced anyway
789 */
790 int iv;
791 int ib;
792 double x = 0.0;
793 double y = 0.0;
794 double z = 0.0;
795 double m = 0.0;
796 double last_x = 0.0;
797 double last_y = 0.0;
798 double last_z = 0.0;
799 int points;
800 gaiaPointPtr point;
801 gaiaLinestringPtr line;
802 gaiaLinestringPtr new_line;
803 gaiaPolygonPtr polyg;
804 gaiaPolygonPtr new_polyg;
805 gaiaGeomCollPtr new_geom;
806 gaiaRingPtr i_ring;
807 gaiaRingPtr o_ring;
808 if (!geom)
809 return NULL;
810 if (geom->DimensionModel == GAIA_XY_Z)
811 new_geom = gaiaAllocGeomCollXYZ ();
812 else if (geom->DimensionModel == GAIA_XY_M)
813 new_geom = gaiaAllocGeomCollXYM ();
814 else if (geom->DimensionModel == GAIA_XY_Z_M)
815 new_geom = gaiaAllocGeomCollXYZM ();
816 else
817 new_geom = gaiaAllocGeomColl ();
818 new_geom->Srid = geom->Srid;
819 new_geom->DeclaredType = geom->DeclaredType;
820 point = geom->FirstPoint;
821 while (point)
822 {
823 /* copying POINTs */
824 gaiaAddPointToGeomCollXYZM (new_geom, point->X, point->Y, point->Z,
825 point->M);
826 point = point->Next;
827 }
828 line = geom->FirstLinestring;
829 while (line)
830 {
831 /* sanitizing LINESTRINGs */
832 points = 0;
833 for (iv = 0; iv < line->Points; iv++)
834 {
835 /* PASS I - checking points */
836 z = 0.0;
837 m = 0.0;
838 if (line->DimensionModel == GAIA_XY_Z)
839 {
840 gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
841 }
842 else if (line->DimensionModel == GAIA_XY_M)
843 {
844 gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
845 }
846 else if (line->DimensionModel == GAIA_XY_Z_M)
847 {
848 gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
849 }
850 else
851 {
852 gaiaGetPoint (line->Coords, iv, &x, &y);
853 }
854 if (iv > 0)
855 {
856 if (last_x == x && last_y == y && last_z == z)
857 ;
858 else
859 points++;
860 }
861 else
862 points++;
863 last_x = x;
864 last_y = y;
865 last_z = z;
866 }
867 if (points < 2)
868 {
869 /* illegal LINESTRING - copying the original one */
870 new_line = gaiaAddLinestringToGeomColl (new_geom, line->Points);
871 gaiaCopyLinestringCoords (new_line, line);
872 }
873 else
874 {
875 /* valid LINESTRING - sanitizing */
876 new_line = gaiaAddLinestringToGeomColl (new_geom, points);
877 points = 0;
878 for (iv = 0; iv < line->Points; iv++)
879 {
880 /* PASS II - inserting points */
881 z = 0.0;
882 m = 0.0;
883 if (line->DimensionModel == GAIA_XY_Z)
884 {
885 gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
886 }
887 else if (line->DimensionModel == GAIA_XY_M)
888 {
889 gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
890 }
891 else if (line->DimensionModel == GAIA_XY_Z_M)
892 {
893 gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
894 }
895 else
896 {
897 gaiaGetPoint (line->Coords, iv, &x, &y);
898 }
899 if (iv > 0)
900 {
901 if (last_x == x && last_y == y && last_z == z)
902 ;
903 else
904 {
905 if (new_line->DimensionModel == GAIA_XY_Z)
906 {
907 gaiaSetPointXYZ (new_line->Coords,
908 points, x, y, z);
909 }
910 else if (new_line->DimensionModel ==
911 GAIA_XY_M)
912 {
913 gaiaSetPointXYM (new_line->Coords,
914 points, x, y, m);
915 }
916 else if (new_line->DimensionModel ==
917 GAIA_XY_Z_M)
918 {
919 gaiaSetPointXYZM (new_line->Coords,
920 points, x, y, z, m);
921 }
922 else
923 {
924 gaiaSetPoint (new_line->Coords, points,
925 x, y);
926 }
927 points++;
928 }
929 }
930 else
931 {
932 if (new_line->DimensionModel == GAIA_XY_Z)
933 {
934 gaiaSetPointXYZ (new_line->Coords, points, x,
935 y, z);
936 }
937 else if (new_line->DimensionModel == GAIA_XY_M)
938 {
939 gaiaSetPointXYM (new_line->Coords, points, x,
940 y, m);
941 }
942 else if (new_line->DimensionModel == GAIA_XY_Z_M)
943 {
944 gaiaSetPointXYZM (new_line->Coords, points, x,
945 y, z, m);
946 }
947 else
948 {
949 gaiaSetPoint (new_line->Coords, points, x, y);
950 }
951 points++;
952 }
953 last_x = x;
954 last_y = y;
955 last_z = z;
956 }
957 }
958 line = line->Next;
959 }
960 polyg = geom->FirstPolygon;
961 while (polyg)
962 {
963 /* copying POLYGONs */
964 i_ring = polyg->Exterior;
965 /* sanitizing EXTERIOR RING */
966 points = 0;
967 for (iv = 0; iv < i_ring->Points; iv++)
968 {
969 /* PASS I - checking points */
970 z = 0.0;
971 m = 0.0;
972 if (i_ring->DimensionModel == GAIA_XY_Z)
973 {
974 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
975 }
976 else if (i_ring->DimensionModel == GAIA_XY_M)
977 {
978 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
979 }
980 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
981 {
982 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z, &m);
983 }
984 else
985 {
986 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
987 }
988 if (iv > 0)
989 {
990 if (last_x == x && last_y == y && last_z == z)
991 ;
992 else
993 points++;
994 }
995 else
996 points++;
997 last_x = x;
998 last_y = y;
999 last_z = z;
1000 }
1001 if (last_x == x && last_y == y && last_z == z)
1002 ;
1003 else
1004 {
1005 /* forcing RING closure */
1006 points++;
1007 }
1008 if (points < 4)
1009 {
1010 /* illegal RING - copying the original one */
1011 new_polyg =
1012 gaiaAddPolygonToGeomColl (new_geom, i_ring->Points,
1013 polyg->NumInteriors);
1014 o_ring = new_polyg->Exterior;
1015 gaiaCopyRingCoords (o_ring, i_ring);
1016 }
1017 else
1018 {
1019 /* valid RING - sanitizing */
1020 new_polyg =
1021 gaiaAddPolygonToGeomColl (new_geom, points,
1022 polyg->NumInteriors);
1023 o_ring = new_polyg->Exterior;
1024 points = 0;
1025 for (iv = 0; iv < i_ring->Points; iv++)
1026 {
1027 /* PASS II - inserting points */
1028 z = 0.0;
1029 m = 0.0;
1030 if (i_ring->DimensionModel == GAIA_XY_Z)
1031 {
1032 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
1033 }
1034 else if (i_ring->DimensionModel == GAIA_XY_M)
1035 {
1036 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
1037 }
1038 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1039 {
1040 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z,
1041 &m);
1042 }
1043 else
1044 {
1045 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1046 }
1047 if (iv > 0)
1048 {
1049 if (last_x == x && last_y == y && last_z == z)
1050 ;
1051 else
1052 {
1053 if (o_ring->DimensionModel == GAIA_XY_Z)
1054 {
1055 gaiaSetPointXYZ (o_ring->Coords, points,
1056 x, y, z);
1057 }
1058 else if (o_ring->DimensionModel == GAIA_XY_M)
1059 {
1060 gaiaSetPointXYM (o_ring->Coords, points,
1061 x, y, m);
1062 }
1063 else if (o_ring->DimensionModel ==
1064 GAIA_XY_Z_M)
1065 {
1066 gaiaSetPointXYZM (o_ring->Coords,
1067 points, x, y, z, m);
1068 }
1069 else
1070 {
1071 gaiaSetPoint (o_ring->Coords, points, x,
1072 y);
1073 }
1074 points++;
1075 }
1076 }
1077 else
1078 {
1079 if (o_ring->DimensionModel == GAIA_XY_Z)
1080 {
1081 gaiaSetPointXYZ (o_ring->Coords, points, x,
1082 y, z);
1083 }
1084 else if (o_ring->DimensionModel == GAIA_XY_M)
1085 {
1086 gaiaSetPointXYM (o_ring->Coords, points, x,
1087 y, m);
1088 }
1089 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1090 {
1091 gaiaSetPointXYZM (o_ring->Coords, points, x,
1092 y, z, m);
1093 }
1094 else
1095 {
1096 gaiaSetPoint (o_ring->Coords, points, x, y);
1097 }
1098 points++;
1099 }
1100 last_x = x;
1101 last_y = y;
1102 last_z = z;
1103 }
1104 }
1105 /* PASS III - forcing RING closure */
1106 z = 0.0;
1107 m = 0.0;
1108 if (i_ring->DimensionModel == GAIA_XY_Z)
1109 {
1110 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1111 }
1112 else if (i_ring->DimensionModel == GAIA_XY_M)
1113 {
1114 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1115 }
1116 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1117 {
1118 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z, &m);
1119 }
1120 else
1121 {
1122 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1123 }
1124 points = o_ring->Points - 1;
1125 if (o_ring->DimensionModel == GAIA_XY_Z)
1126 {
1127 gaiaSetPointXYZ (o_ring->Coords, points, x, y, z);
1128 }
1129 else if (o_ring->DimensionModel == GAIA_XY_M)
1130 {
1131 gaiaSetPointXYM (o_ring->Coords, points, x, y, m);
1132 }
1133 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1134 {
1135 gaiaSetPointXYZM (o_ring->Coords, points, x, y, z, m);
1136 }
1137 else
1138 {
1139 gaiaSetPoint (o_ring->Coords, points, x, y);
1140 }
1141 for (ib = 0; ib < new_polyg->NumInteriors; ib++)
1142 {
1143 /* copying each INTERIOR RING [if any] */
1144 i_ring = polyg->Interiors + ib;
1145 /* sanitizing an INTERIOR RING */
1146 points = 0;
1147 for (iv = 0; iv < i_ring->Points; iv++)
1148 {
1149 /* PASS I - checking points */
1150 z = 0.0;
1151 m = 0.0;
1152 if (i_ring->DimensionModel == GAIA_XY_Z)
1153 {
1154 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
1155 }
1156 else if (i_ring->DimensionModel == GAIA_XY_M)
1157 {
1158 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
1159 }
1160 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1161 {
1162 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z,
1163 &m);
1164 }
1165 else
1166 {
1167 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1168 }
1169 if (iv > 0)
1170 {
1171 if (last_x == x && last_y == y && last_z == z)
1172 ;
1173 else
1174 points++;
1175 }
1176 else
1177 points++;
1178 last_x = x;
1179 last_y = y;
1180 last_z = z;
1181 }
1182 if (last_x == x && last_y == y && last_z == z)
1183 ;
1184 else
1185 {
1186 /* forcing RING closure */
1187 points++;
1188 }
1189 if (points < 4)
1190 {
1191 /* illegal RING - copying the original one */
1192 o_ring =
1193 gaiaAddInteriorRing (new_polyg, ib, i_ring->Points);
1194 gaiaCopyRingCoords (o_ring, i_ring);
1195 }
1196 else
1197 {
1198 /* valid RING - sanitizing */
1199 o_ring = gaiaAddInteriorRing (new_polyg, ib, points);
1200 points = 0;
1201 for (iv = 0; iv < i_ring->Points; iv++)
1202 {
1203 /* PASS II - inserting points */
1204 z = 0.0;
1205 m = 0.0;
1206 if (i_ring->DimensionModel == GAIA_XY_Z)
1207 {
1208 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y,
1209 &z);
1210 }
1211 else if (i_ring->DimensionModel == GAIA_XY_M)
1212 {
1213 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y,
1214 &m);
1215 }
1216 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1217 {
1218 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y,
1219 &z, &m);
1220 }
1221 else
1222 {
1223 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1224 }
1225 if (iv > 0)
1226 {
1227 if (last_x == x && last_y == y && last_z == z)
1228 ;
1229 else
1230 {
1231 if (o_ring->DimensionModel == GAIA_XY_Z)
1232 {
1233 gaiaSetPointXYZ (o_ring->Coords,
1234 points, x, y, z);
1235 }
1236 else if (o_ring->DimensionModel ==
1237 GAIA_XY_M)
1238 {
1239 gaiaSetPointXYM (o_ring->Coords,
1240 points, x, y, m);
1241 }
1242 else if (o_ring->DimensionModel ==
1243 GAIA_XY_Z_M)
1244 {
1245 gaiaSetPointXYZM (o_ring->Coords,
1246 points, x, y, z,
1247 m);
1248 }
1249 else
1250 {
1251 gaiaSetPoint (o_ring->Coords,
1252 points, x, y);
1253 }
1254 points++;
1255 }
1256 }
1257 else
1258 {
1259 if (o_ring->DimensionModel == GAIA_XY_Z)
1260 {
1261 gaiaSetPointXYZ (o_ring->Coords, points,
1262 x, y, z);
1263 }
1264 else if (o_ring->DimensionModel == GAIA_XY_M)
1265 {
1266 gaiaSetPointXYM (o_ring->Coords, points,
1267 x, y, m);
1268 }
1269 else if (o_ring->DimensionModel ==
1270 GAIA_XY_Z_M)
1271 {
1272 gaiaSetPointXYZM (o_ring->Coords,
1273 points, x, y, z, m);
1274 }
1275 else
1276 {
1277 gaiaSetPoint (o_ring->Coords, points, x,
1278 y);
1279 }
1280 points++;
1281 }
1282 last_x = x;
1283 last_y = y;
1284 last_z = z;
1285 }
1286 /* PASS III - forcing RING closure */
1287 z = 0.0;
1288 m = 0.0;
1289 if (i_ring->DimensionModel == GAIA_XY_Z)
1290 {
1291 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1292 }
1293 else if (i_ring->DimensionModel == GAIA_XY_M)
1294 {
1295 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1296 }
1297 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1298 {
1299 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z,
1300 &m);
1301 }
1302 else
1303 {
1304 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1305 }
1306 points = o_ring->Points - 1;
1307 if (o_ring->DimensionModel == GAIA_XY_Z)
1308 {
1309 gaiaSetPointXYZ (o_ring->Coords, points, x, y, z);
1310 }
1311 else if (o_ring->DimensionModel == GAIA_XY_M)
1312 {
1313 gaiaSetPointXYM (o_ring->Coords, points, x, y, m);
1314 }
1315 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1316 {
1317 gaiaSetPointXYZM (o_ring->Coords, points, x, y, z,
1318 m);
1319 }
1320 else
1321 {
1322 gaiaSetPoint (o_ring->Coords, points, x, y);
1323 }
1324 }
1325 }
1326 polyg = polyg->Next;
1327 }
1328 return new_geom;
1329 }
1330
1331 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaEnsureClosedRings(gaiaGeomCollPtr geom)1332 gaiaEnsureClosedRings (gaiaGeomCollPtr geom)
1333 {
1334 /*
1335 / sanitizes a GEOMETRYCOLLECTION:
1336 / - ring closure is enforced anyway
1337 */
1338 int iv;
1339 int ib;
1340 double x = 0.0;
1341 double y = 0.0;
1342 double z = 0.0;
1343 double m = 0.0;
1344 double last_x = 0.0;
1345 double last_y = 0.0;
1346 double last_z = 0.0;
1347 int points;
1348 int end_point;
1349 int enforce;
1350 gaiaPointPtr point;
1351 gaiaLinestringPtr line;
1352 gaiaLinestringPtr new_line;
1353 gaiaPolygonPtr polyg;
1354 gaiaPolygonPtr new_polyg;
1355 gaiaGeomCollPtr new_geom;
1356 gaiaRingPtr i_ring;
1357 gaiaRingPtr o_ring;
1358 if (!geom)
1359 return NULL;
1360 if (geom->DimensionModel == GAIA_XY_Z)
1361 new_geom = gaiaAllocGeomCollXYZ ();
1362 else if (geom->DimensionModel == GAIA_XY_M)
1363 new_geom = gaiaAllocGeomCollXYM ();
1364 else if (geom->DimensionModel == GAIA_XY_Z_M)
1365 new_geom = gaiaAllocGeomCollXYZM ();
1366 else
1367 new_geom = gaiaAllocGeomColl ();
1368 new_geom->Srid = geom->Srid;
1369 new_geom->DeclaredType = geom->DeclaredType;
1370 point = geom->FirstPoint;
1371 while (point)
1372 {
1373 /* copying POINTs */
1374 gaiaAddPointToGeomCollXYZM (new_geom, point->X, point->Y, point->Z,
1375 point->M);
1376 point = point->Next;
1377 }
1378 line = geom->FirstLinestring;
1379 while (line)
1380 {
1381 /* copying LINESTRINGs */
1382 new_line = gaiaAddLinestringToGeomColl (new_geom, line->Points);
1383 gaiaCopyLinestringCoords (new_line, line);
1384 line = line->Next;
1385 }
1386 polyg = geom->FirstPolygon;
1387 while (polyg)
1388 {
1389 /* copying POLYGONs */
1390 i_ring = polyg->Exterior;
1391 /* sanitizing EXTERIOR RING */
1392 points = i_ring->Points;
1393 z = 0.0;
1394 m = 0.0;
1395 if (i_ring->DimensionModel == GAIA_XY_Z)
1396 {
1397 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1398 }
1399 else if (i_ring->DimensionModel == GAIA_XY_M)
1400 {
1401 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1402 }
1403 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1404 {
1405 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z, &m);
1406 }
1407 else
1408 {
1409 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1410 }
1411 last_z = 0.0;
1412 m = 0.0;
1413 end_point = points - 1;
1414 if (i_ring->DimensionModel == GAIA_XY_Z)
1415 {
1416 gaiaGetPointXYZ (i_ring->Coords, end_point, &last_x, &last_y,
1417 &last_z);
1418 }
1419 else if (i_ring->DimensionModel == GAIA_XY_M)
1420 {
1421 gaiaGetPointXYM (i_ring->Coords, end_point, &last_x, &last_y,
1422 &m);
1423 }
1424 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1425 {
1426 gaiaGetPointXYZM (i_ring->Coords, end_point, &last_x, &last_y,
1427 &last_z, &m);
1428 }
1429 else
1430 {
1431 gaiaGetPoint (i_ring->Coords, end_point, &last_x, &last_y);
1432 }
1433 if (last_x == x && last_y == y && last_z == z)
1434 enforce = 0;
1435 else
1436 {
1437 /* forcing RING closure */
1438 enforce = 1;
1439 }
1440 if (!enforce)
1441 {
1442 /* already closed */
1443 new_polyg =
1444 gaiaAddPolygonToGeomColl (new_geom, i_ring->Points,
1445 polyg->NumInteriors);
1446 o_ring = new_polyg->Exterior;
1447 gaiaCopyRingCoords (o_ring, i_ring);
1448 }
1449 else
1450 {
1451 /* forcing closure */
1452 new_polyg =
1453 gaiaAddPolygonToGeomColl (new_geom, i_ring->Points + 1,
1454 polyg->NumInteriors);
1455 o_ring = new_polyg->Exterior;
1456 for (iv = 0; iv < i_ring->Points; iv++)
1457 {
1458 /* inserting points */
1459 z = 0.0;
1460 m = 0.0;
1461 if (i_ring->DimensionModel == GAIA_XY_Z)
1462 {
1463 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
1464 }
1465 else if (i_ring->DimensionModel == GAIA_XY_M)
1466 {
1467 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
1468 }
1469 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1470 {
1471 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z,
1472 &m);
1473 }
1474 else
1475 {
1476 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1477 }
1478 if (o_ring->DimensionModel == GAIA_XY_Z)
1479 {
1480 gaiaSetPointXYZ (o_ring->Coords, iv, x, y, z);
1481 }
1482 else if (o_ring->DimensionModel == GAIA_XY_M)
1483 {
1484 gaiaSetPointXYM (o_ring->Coords, iv, x, y, m);
1485 }
1486 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1487 {
1488 gaiaSetPointXYZM (o_ring->Coords, iv, x, y, z, m);
1489 }
1490 else
1491 {
1492 gaiaSetPoint (o_ring->Coords, iv, x, y);
1493 }
1494 }
1495 /* adding last vertex so to ensure closure */
1496 z = 0.0;
1497 m = 0.0;
1498 if (i_ring->DimensionModel == GAIA_XY_Z)
1499 {
1500 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1501 }
1502 else if (i_ring->DimensionModel == GAIA_XY_M)
1503 {
1504 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1505 }
1506 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1507 {
1508 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z, &m);
1509 }
1510 else
1511 {
1512 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1513 }
1514 end_point = o_ring->Points - 1;
1515 if (o_ring->DimensionModel == GAIA_XY_Z)
1516 {
1517 gaiaSetPointXYZ (o_ring->Coords, end_point, x, y, z);
1518 }
1519 else if (o_ring->DimensionModel == GAIA_XY_M)
1520 {
1521 gaiaSetPointXYM (o_ring->Coords, end_point, x, y, m);
1522 }
1523 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1524 {
1525 gaiaSetPointXYZM (o_ring->Coords, end_point, x, y, z, m);
1526 }
1527 else
1528 {
1529 gaiaSetPoint (o_ring->Coords, end_point, x, y);
1530 }
1531 }
1532
1533 for (ib = 0; ib < new_polyg->NumInteriors; ib++)
1534 {
1535 /* copying each INTERIOR RING [if any] */
1536 i_ring = polyg->Interiors + ib;
1537 points = i_ring->Points;
1538 z = 0.0;
1539 m = 0.0;
1540 if (i_ring->DimensionModel == GAIA_XY_Z)
1541 {
1542 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1543 }
1544 else if (i_ring->DimensionModel == GAIA_XY_M)
1545 {
1546 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1547 }
1548 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1549 {
1550 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z, &m);
1551 }
1552 else
1553 {
1554 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1555 }
1556 last_z = 0.0;
1557 m = 0.0;
1558 end_point = points - 1;
1559 if (i_ring->DimensionModel == GAIA_XY_Z)
1560 {
1561 gaiaGetPointXYZ (i_ring->Coords, end_point, &last_x,
1562 &last_y, &last_z);
1563 }
1564 else if (i_ring->DimensionModel == GAIA_XY_M)
1565 {
1566 gaiaGetPointXYM (i_ring->Coords, end_point, &last_x,
1567 &last_y, &m);
1568 }
1569 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1570 {
1571 gaiaGetPointXYZM (i_ring->Coords, end_point, &last_x,
1572 &last_y, &last_z, &m);
1573 }
1574 else
1575 {
1576 gaiaGetPoint (i_ring->Coords, end_point, &last_x,
1577 &last_y);
1578 }
1579 last_z = z;
1580 if (last_x == x && last_y == y && last_z == z)
1581 enforce = 0;
1582 else
1583 {
1584 /* forcing RING closure */
1585 enforce = 1;
1586 }
1587 if (!enforce)
1588 {
1589 /* already closed */
1590 o_ring =
1591 gaiaAddInteriorRing (new_polyg, ib, i_ring->Points);
1592 gaiaCopyRingCoords (o_ring, i_ring);
1593 }
1594 else
1595 {
1596 /* forcing closure */
1597 o_ring =
1598 gaiaAddInteriorRing (new_polyg, ib,
1599 i_ring->Points + 1);
1600 for (iv = 0; iv < i_ring->Points; iv++)
1601 {
1602 /* inserting points */
1603 z = 0.0;
1604 m = 0.0;
1605 if (i_ring->DimensionModel == GAIA_XY_Z)
1606 {
1607 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y,
1608 &z);
1609 }
1610 else if (i_ring->DimensionModel == GAIA_XY_M)
1611 {
1612 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y,
1613 &m);
1614 }
1615 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1616 {
1617 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y,
1618 &z, &m);
1619 }
1620 else
1621 {
1622 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1623 }
1624 if (o_ring->DimensionModel == GAIA_XY_Z)
1625 {
1626 gaiaSetPointXYZ (o_ring->Coords, iv, x, y, z);
1627 }
1628 else if (o_ring->DimensionModel == GAIA_XY_M)
1629 {
1630 gaiaSetPointXYM (o_ring->Coords, iv, x, y, m);
1631 }
1632 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1633 {
1634 gaiaSetPointXYZM (o_ring->Coords,
1635 iv, x, y, z, m);
1636 }
1637 else
1638 {
1639 gaiaSetPoint (o_ring->Coords, iv, x, y);
1640 }
1641 }
1642 /* adding last vertex so to ensure closure */
1643 z = 0.0;
1644 m = 0.0;
1645 if (i_ring->DimensionModel == GAIA_XY_Z)
1646 {
1647 gaiaGetPointXYZ (i_ring->Coords, 0, &x, &y, &z);
1648 }
1649 else if (i_ring->DimensionModel == GAIA_XY_M)
1650 {
1651 gaiaGetPointXYM (i_ring->Coords, 0, &x, &y, &m);
1652 }
1653 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1654 {
1655 gaiaGetPointXYZM (i_ring->Coords, 0, &x, &y, &z,
1656 &m);
1657 }
1658 else
1659 {
1660 gaiaGetPoint (i_ring->Coords, 0, &x, &y);
1661 }
1662 end_point = o_ring->Points - 1;
1663 if (o_ring->DimensionModel == GAIA_XY_Z)
1664 {
1665 gaiaSetPointXYZ (o_ring->Coords, end_point, x, y,
1666 z);
1667 }
1668 else if (o_ring->DimensionModel == GAIA_XY_M)
1669 {
1670 gaiaSetPointXYM (o_ring->Coords, end_point, x, y,
1671 m);
1672 }
1673 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
1674 {
1675 gaiaSetPointXYZM (o_ring->Coords, end_point, x, y,
1676 z, m);
1677 }
1678 else
1679 {
1680 gaiaSetPoint (o_ring->Coords, end_point, x, y);
1681 }
1682 }
1683 }
1684 polyg = polyg->Next;
1685 }
1686 return new_geom;
1687 }
1688
1689 static double
point_point_distance(double x1,double y1,double x2,double y2)1690 point_point_distance (double x1, double y1, double x2, double y2)
1691 {
1692 return sqrt (((x1 - x2) * (x1 - x2)) + ((y1 - y2) * (y1 - y2)));
1693 }
1694
1695 static int
repeated_matching_point(gaiaGeomCollPtr geom,double x,double y,double z,double tolerance)1696 repeated_matching_point (gaiaGeomCollPtr geom, double x, double y, double z,
1697 double tolerance)
1698 {
1699 /* checking for duplicate points - multipoint */
1700 gaiaPointPtr point;
1701 point = geom->FirstPoint;
1702 while (point)
1703 {
1704 if (tolerance <= 0.0)
1705 {
1706 if (point->X == x && point->Y == y && point->Z == z)
1707 return 1;
1708 }
1709 else
1710 {
1711 if (point_point_distance (x, y, point->X, point->Y) <=
1712 tolerance)
1713 return 1;
1714 }
1715 point = point->Next;
1716 }
1717 return 0;
1718 }
1719
1720 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaRemoveRepeatedPoints(gaiaGeomCollPtr geom,double tolerance)1721 gaiaRemoveRepeatedPoints (gaiaGeomCollPtr geom, double tolerance)
1722 {
1723 /*
1724 / sanitizes a GEOMETRYCOLLECTION:
1725 / - repeated vertices are omitted
1726 */
1727 int iv;
1728 int ib;
1729 double x = 0.0;
1730 double y = 0.0;
1731 double z = 0.0;
1732 double m = 0.0;
1733 double last_x = 0.0;
1734 double last_y = 0.0;
1735 double last_z = 0.0;
1736 int points;
1737 gaiaPointPtr point;
1738 gaiaLinestringPtr line;
1739 gaiaLinestringPtr new_line;
1740 gaiaPolygonPtr polyg;
1741 gaiaPolygonPtr new_polyg;
1742 gaiaGeomCollPtr new_geom;
1743 gaiaRingPtr i_ring;
1744 gaiaRingPtr o_ring;
1745 if (!geom)
1746 return NULL;
1747 if (geom->DimensionModel == GAIA_XY_Z)
1748 new_geom = gaiaAllocGeomCollXYZ ();
1749 else if (geom->DimensionModel == GAIA_XY_M)
1750 new_geom = gaiaAllocGeomCollXYM ();
1751 else if (geom->DimensionModel == GAIA_XY_Z_M)
1752 new_geom = gaiaAllocGeomCollXYZM ();
1753 else
1754 new_geom = gaiaAllocGeomColl ();
1755 new_geom->Srid = geom->Srid;
1756 new_geom->DeclaredType = geom->DeclaredType;
1757 point = geom->FirstPoint;
1758 while (point)
1759 {
1760 /* copying POINTs */
1761 if (!repeated_matching_point
1762 (new_geom, point->X, point->Y, point->Z, tolerance))
1763 gaiaAddPointToGeomCollXYZM (new_geom, point->X, point->Y,
1764 point->Z, point->M);
1765 point = point->Next;
1766 }
1767 line = geom->FirstLinestring;
1768 while (line)
1769 {
1770 /* sanitizing LINESTRINGs */
1771 points = 0;
1772 for (iv = 0; iv < line->Points; iv++)
1773 {
1774 /* PASS I - checking points */
1775 z = 0.0;
1776 m = 0.0;
1777 if (line->DimensionModel == GAIA_XY_Z)
1778 {
1779 gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
1780 }
1781 else if (line->DimensionModel == GAIA_XY_M)
1782 {
1783 gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
1784 }
1785 else if (line->DimensionModel == GAIA_XY_Z_M)
1786 {
1787 gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
1788 }
1789 else
1790 {
1791 gaiaGetPoint (line->Coords, iv, &x, &y);
1792 }
1793 if (iv > 0)
1794 {
1795 if (tolerance <= 0.0)
1796 {
1797 if (last_x == x && last_y == y && last_z == z)
1798 ;
1799 else
1800 points++;
1801 }
1802 else
1803 {
1804 if (point_point_distance (x, y, last_x, last_y) <=
1805 tolerance)
1806 ;
1807 else
1808 points++;
1809 }
1810 }
1811 else
1812 points++;
1813 last_x = x;
1814 last_y = y;
1815 last_z = z;
1816 }
1817 if (points < 2)
1818 {
1819 /* illegal LINESTRING - copying the original one */
1820 new_line = gaiaAddLinestringToGeomColl (new_geom, line->Points);
1821 gaiaCopyLinestringCoords (new_line, line);
1822 }
1823 else
1824 {
1825 /* valid LINESTRING - sanitizing */
1826 new_line = gaiaAddLinestringToGeomColl (new_geom, points);
1827 points = 0;
1828 for (iv = 0; iv < line->Points; iv++)
1829 {
1830 /* PASS II - inserting points */
1831 z = 0.0;
1832 m = 0.0;
1833 if (line->DimensionModel == GAIA_XY_Z)
1834 {
1835 gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
1836 }
1837 else if (line->DimensionModel == GAIA_XY_M)
1838 {
1839 gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
1840 }
1841 else if (line->DimensionModel == GAIA_XY_Z_M)
1842 {
1843 gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
1844 }
1845 else
1846 {
1847 gaiaGetPoint (line->Coords, iv, &x, &y);
1848 }
1849 if (iv > 0)
1850 {
1851 int skip = 0;
1852 if (tolerance <= 0.0)
1853 {
1854 if (last_x == x && last_y == y && last_z == z)
1855 skip = 1;
1856 }
1857 else
1858 {
1859 if (point_point_distance
1860 (x, y, last_x, last_y) <= tolerance)
1861 skip = 1;
1862 }
1863 if (skip)
1864 ;
1865 else
1866 {
1867 if (new_line->DimensionModel == GAIA_XY_Z)
1868 {
1869 gaiaSetPointXYZ (new_line->Coords,
1870 points, x, y, z);
1871 }
1872 else if (new_line->DimensionModel ==
1873 GAIA_XY_M)
1874 {
1875 gaiaSetPointXYM (new_line->Coords,
1876 points, x, y, m);
1877 }
1878 else if (new_line->DimensionModel ==
1879 GAIA_XY_Z_M)
1880 {
1881 gaiaSetPointXYZM (new_line->Coords,
1882 points, x, y, z, m);
1883 }
1884 else
1885 {
1886 gaiaSetPoint (new_line->Coords, points,
1887 x, y);
1888 }
1889 points++;
1890 }
1891 }
1892 else
1893 {
1894 if (new_line->DimensionModel == GAIA_XY_Z)
1895 {
1896 gaiaSetPointXYZ (new_line->Coords, points, x,
1897 y, z);
1898 }
1899 else if (new_line->DimensionModel == GAIA_XY_M)
1900 {
1901 gaiaSetPointXYM (new_line->Coords, points, x,
1902 y, m);
1903 }
1904 else if (new_line->DimensionModel == GAIA_XY_Z_M)
1905 {
1906 gaiaSetPointXYZM (new_line->Coords, points, x,
1907 y, z, m);
1908 }
1909 else
1910 {
1911 gaiaSetPoint (new_line->Coords, points, x, y);
1912 }
1913 points++;
1914 }
1915 last_x = x;
1916 last_y = y;
1917 last_z = z;
1918 }
1919 }
1920 line = line->Next;
1921 }
1922 polyg = geom->FirstPolygon;
1923 while (polyg)
1924 {
1925 /* copying POLYGONs */
1926 i_ring = polyg->Exterior;
1927 /* sanitizing EXTERIOR RING */
1928 points = 0;
1929 for (iv = 0; iv < i_ring->Points; iv++)
1930 {
1931 /* PASS I - checking points */
1932 z = 0.0;
1933 m = 0.0;
1934 if (i_ring->DimensionModel == GAIA_XY_Z)
1935 {
1936 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
1937 }
1938 else if (i_ring->DimensionModel == GAIA_XY_M)
1939 {
1940 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
1941 }
1942 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
1943 {
1944 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z, &m);
1945 }
1946 else
1947 {
1948 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
1949 }
1950 if (iv > 0)
1951 {
1952 if (tolerance <= 0.0)
1953 {
1954 if (last_x == x && last_y == y && last_z == z)
1955 ;
1956 else
1957 points++;
1958 }
1959 else
1960 {
1961 if (point_point_distance (x, y, last_x, last_y) <=
1962 tolerance)
1963 ;
1964 else
1965 points++;
1966 }
1967 }
1968 else
1969 points++;
1970 last_x = x;
1971 last_y = y;
1972 last_z = z;
1973 }
1974 if (points < 4)
1975 {
1976 /* illegal RING - copying the original one */
1977 new_polyg =
1978 gaiaAddPolygonToGeomColl (new_geom, i_ring->Points,
1979 polyg->NumInteriors);
1980 o_ring = new_polyg->Exterior;
1981 gaiaCopyRingCoords (o_ring, i_ring);
1982 }
1983 else
1984 {
1985 /* valid RING - sanitizing */
1986 new_polyg =
1987 gaiaAddPolygonToGeomColl (new_geom, points,
1988 polyg->NumInteriors);
1989 o_ring = new_polyg->Exterior;
1990 points = 0;
1991 for (iv = 0; iv < i_ring->Points; iv++)
1992 {
1993 /* PASS II - inserting points */
1994 z = 0.0;
1995 m = 0.0;
1996 if (i_ring->DimensionModel == GAIA_XY_Z)
1997 {
1998 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
1999 }
2000 else if (i_ring->DimensionModel == GAIA_XY_M)
2001 {
2002 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
2003 }
2004 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
2005 {
2006 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z,
2007 &m);
2008 }
2009 else
2010 {
2011 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
2012 }
2013 if (iv > 0)
2014 {
2015 int skip = 0;
2016 if (tolerance <= 0.0)
2017 {
2018 if (last_x == x && last_y == y && last_z == z)
2019 skip = 1;
2020 }
2021 else
2022 {
2023 if (point_point_distance
2024 (x, y, last_x, last_y) <= tolerance)
2025 skip = 1;
2026 }
2027 if (skip)
2028 ;
2029 else
2030 {
2031 if (o_ring->DimensionModel == GAIA_XY_Z)
2032 {
2033 gaiaSetPointXYZ (o_ring->Coords, points,
2034 x, y, z);
2035 }
2036 else if (o_ring->DimensionModel == GAIA_XY_M)
2037 {
2038 gaiaSetPointXYM (o_ring->Coords, points,
2039 x, y, m);
2040 }
2041 else if (o_ring->DimensionModel ==
2042 GAIA_XY_Z_M)
2043 {
2044 gaiaSetPointXYZM (o_ring->Coords,
2045 points, x, y, z, m);
2046 }
2047 else
2048 {
2049 gaiaSetPoint (o_ring->Coords, points, x,
2050 y);
2051 }
2052 points++;
2053 }
2054 }
2055 else
2056 {
2057 if (o_ring->DimensionModel == GAIA_XY_Z)
2058 {
2059 gaiaSetPointXYZ (o_ring->Coords, points, x,
2060 y, z);
2061 }
2062 else if (o_ring->DimensionModel == GAIA_XY_M)
2063 {
2064 gaiaSetPointXYM (o_ring->Coords, points, x,
2065 y, m);
2066 }
2067 else if (o_ring->DimensionModel == GAIA_XY_Z_M)
2068 {
2069 gaiaSetPointXYZM (o_ring->Coords, points, x,
2070 y, z, m);
2071 }
2072 else
2073 {
2074 gaiaSetPoint (o_ring->Coords, points, x, y);
2075 }
2076 points++;
2077 }
2078 last_x = x;
2079 last_y = y;
2080 last_z = z;
2081 }
2082 }
2083
2084 for (ib = 0; ib < new_polyg->NumInteriors; ib++)
2085 {
2086 /* copying each INTERIOR RING [if any] */
2087 i_ring = polyg->Interiors + ib;
2088 /* sanitizing an INTERIOR RING */
2089 points = 0;
2090 for (iv = 0; iv < i_ring->Points; iv++)
2091 {
2092 /* PASS I - checking points */
2093 z = 0.0;
2094 m = 0.0;
2095 if (i_ring->DimensionModel == GAIA_XY_Z)
2096 {
2097 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y, &z);
2098 }
2099 else if (i_ring->DimensionModel == GAIA_XY_M)
2100 {
2101 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y, &m);
2102 }
2103 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
2104 {
2105 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y, &z,
2106 &m);
2107 }
2108 else
2109 {
2110 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
2111 }
2112 if (iv > 0)
2113 {
2114 if (tolerance <= 0.0)
2115 {
2116 if (last_x == x && last_y == y && last_z == z)
2117 ;
2118 else
2119 points++;
2120 }
2121 else
2122 {
2123 if (point_point_distance
2124 (x, y, last_x, last_y) <= tolerance)
2125 ;
2126 else
2127 points++;
2128 }
2129 }
2130 else
2131 points++;
2132 last_x = x;
2133 last_y = y;
2134 last_z = z;
2135 }
2136 if (points < 4)
2137 {
2138 /* illegal RING - copying the original one */
2139 o_ring =
2140 gaiaAddInteriorRing (new_polyg, ib, i_ring->Points);
2141 gaiaCopyRingCoords (o_ring, i_ring);
2142 }
2143 else
2144 {
2145 /* valid RING - sanitizing */
2146 o_ring = gaiaAddInteriorRing (new_polyg, ib, points);
2147 points = 0;
2148 for (iv = 0; iv < i_ring->Points; iv++)
2149 {
2150 /* PASS II - inserting points */
2151 z = 0.0;
2152 m = 0.0;
2153 if (i_ring->DimensionModel == GAIA_XY_Z)
2154 {
2155 gaiaGetPointXYZ (i_ring->Coords, iv, &x, &y,
2156 &z);
2157 }
2158 else if (i_ring->DimensionModel == GAIA_XY_M)
2159 {
2160 gaiaGetPointXYM (i_ring->Coords, iv, &x, &y,
2161 &m);
2162 }
2163 else if (i_ring->DimensionModel == GAIA_XY_Z_M)
2164 {
2165 gaiaGetPointXYZM (i_ring->Coords, iv, &x, &y,
2166 &z, &m);
2167 }
2168 else
2169 {
2170 gaiaGetPoint (i_ring->Coords, iv, &x, &y);
2171 }
2172 if (iv > 0)
2173 {
2174 int skip = 0;
2175 if (tolerance <= 0.0)
2176 {
2177 if (last_x == x && last_y == y
2178 && last_z == z)
2179 skip = 1;
2180 }
2181 else
2182 {
2183 if (point_point_distance
2184 (x, y, last_x, last_y) <= tolerance)
2185 skip = 1;
2186 }
2187 if (skip)
2188 ;
2189 else
2190 {
2191 if (o_ring->DimensionModel == GAIA_XY_Z)
2192 {
2193 gaiaSetPointXYZ (o_ring->Coords,
2194 points, x, y, z);
2195 }
2196 else if (o_ring->DimensionModel ==
2197 GAIA_XY_M)
2198 {
2199 gaiaSetPointXYM (o_ring->Coords,
2200 points, x, y, m);
2201 }
2202 else if (o_ring->DimensionModel ==
2203 GAIA_XY_Z_M)
2204 {
2205 gaiaSetPointXYZM (o_ring->Coords,
2206 points, x, y, z,
2207 m);
2208 }
2209 else
2210 {
2211 gaiaSetPoint (o_ring->Coords,
2212 points, x, y);
2213 }
2214 points++;
2215 }
2216 }
2217 else
2218 {
2219 if (o_ring->DimensionModel == GAIA_XY_Z)
2220 {
2221 gaiaSetPointXYZ (o_ring->Coords, points,
2222 x, y, z);
2223 }
2224 else if (o_ring->DimensionModel == GAIA_XY_M)
2225 {
2226 gaiaSetPointXYM (o_ring->Coords, points,
2227 x, y, m);
2228 }
2229 else if (o_ring->DimensionModel ==
2230 GAIA_XY_Z_M)
2231 {
2232 gaiaSetPointXYZM (o_ring->Coords,
2233 points, x, y, z, m);
2234 }
2235 else
2236 {
2237 gaiaSetPoint (o_ring->Coords, points, x,
2238 y);
2239 }
2240 points++;
2241 }
2242 last_x = x;
2243 last_y = y;
2244 last_z = z;
2245 }
2246 }
2247 }
2248 polyg = polyg->Next;
2249 }
2250 return new_geom;
2251 }
2252
2253 static int
gaiaIsToxicLinestring(gaiaLinestringPtr line)2254 gaiaIsToxicLinestring (gaiaLinestringPtr line)
2255 {
2256 /* checking a Linestring */
2257 if (line->Points < 2)
2258 return 1;
2259
2260 /* not a valid Linestring, simply a degenerated Point */
2261 return 0;
2262 }
2263
2264 static int
gaiaIsToxicRing(gaiaRingPtr ring)2265 gaiaIsToxicRing (gaiaRingPtr ring)
2266 {
2267 /* checking a Ring */
2268 if (ring->Points < 4)
2269 return 1;
2270 return 0;
2271 }
2272
2273 GAIAGEO_DECLARE int
gaiaIsToxic(gaiaGeomCollPtr geom)2274 gaiaIsToxic (gaiaGeomCollPtr geom)
2275 {
2276 return gaiaIsToxic_r (NULL, geom);
2277 }
2278
2279 GAIAGEO_DECLARE int
gaiaIsToxic_r(const void * cache,gaiaGeomCollPtr geom)2280 gaiaIsToxic_r (const void *cache, gaiaGeomCollPtr geom)
2281 {
2282 /*
2283 / identifying toxic geometries
2284 / i.e. geoms making GEOS to crash !!!
2285 */
2286 int ib;
2287 gaiaPointPtr point;
2288 gaiaLinestringPtr line;
2289 gaiaPolygonPtr polyg;
2290 gaiaRingPtr ring;
2291 if (!geom)
2292 return 0;
2293 if (gaiaIsEmpty (geom))
2294 return 1;
2295 point = geom->FirstPoint;
2296 while (point)
2297 {
2298 /* checking POINTs */
2299 point = point->Next;
2300 }
2301 line = geom->FirstLinestring;
2302 while (line)
2303 {
2304 /* checking LINESTRINGs */
2305 if (gaiaIsToxicLinestring (line))
2306 {
2307 if (cache != NULL)
2308 gaiaSetGeosAuxErrorMsg_r
2309 (cache,
2310 "gaiaIsToxic detected a toxic Linestring: < 2 pts");
2311 else
2312 gaiaSetGeosAuxErrorMsg
2313 ("gaiaIsToxic detected a toxic Linestring: < 2 pts");
2314 return 1;
2315 }
2316 line = line->Next;
2317 }
2318 polyg = geom->FirstPolygon;
2319 while (polyg)
2320 {
2321 /* checking POLYGONs */
2322 ring = polyg->Exterior;
2323 if (gaiaIsToxicRing (ring))
2324 {
2325 if (cache != NULL)
2326 gaiaSetGeosAuxErrorMsg_r
2327 (cache, "gaiaIsToxic detected a toxic Ring: < 4 pts");
2328 else
2329 gaiaSetGeosAuxErrorMsg
2330 ("gaiaIsToxic detected a toxic Ring: < 4 pts");
2331 return 1;
2332 }
2333 for (ib = 0; ib < polyg->NumInteriors; ib++)
2334 {
2335 ring = polyg->Interiors + ib;
2336 if (gaiaIsToxicRing (ring))
2337 {
2338 if (cache != NULL)
2339 gaiaSetGeosAuxErrorMsg_r
2340 (cache,
2341 "gaiaIsToxic detected a toxic Ring: < 4 pts");
2342 else
2343 gaiaSetGeosAuxErrorMsg
2344 ("gaiaIsToxic detected a toxic Ring: < 4 pts");
2345 return 1;
2346 }
2347 }
2348 polyg = polyg->Next;
2349 }
2350 return 0;
2351 }
2352
2353 GAIAGEO_DECLARE int
gaiaIsNotClosedRing(gaiaRingPtr ring)2354 gaiaIsNotClosedRing (gaiaRingPtr ring)
2355 {
2356 return gaiaIsNotClosedRing_r (NULL, ring);
2357 }
2358
2359 GAIAGEO_DECLARE int
gaiaIsNotClosedRing_r(const void * cache,gaiaRingPtr ring)2360 gaiaIsNotClosedRing_r (const void *cache, gaiaRingPtr ring)
2361 {
2362 /* checking a Ring for closure */
2363 double x0;
2364 double y0;
2365 double z0;
2366 double m0;
2367 double x1;
2368 double y1;
2369 double z1;
2370 double m1;
2371 gaiaRingGetPoint (ring, 0, &x0, &y0, &z0, &m0);
2372 gaiaRingGetPoint (ring, ring->Points - 1, &x1, &y1, &z1, &m1);
2373 if (x0 == x1 && y0 == y1 && z0 == z1 && m0 == m1)
2374 return 0;
2375 else
2376 {
2377 if (cache != NULL)
2378 gaiaSetGeosAuxErrorMsg_r (cache,
2379 "gaia detected a not-closed Ring");
2380 else
2381 gaiaSetGeosAuxErrorMsg ("gaia detected a not-closed Ring");
2382 return 1;
2383 }
2384 }
2385
2386 GAIAGEO_DECLARE int
gaiaIsNotClosedGeomColl(gaiaGeomCollPtr geom)2387 gaiaIsNotClosedGeomColl (gaiaGeomCollPtr geom)
2388 {
2389 return gaiaIsNotClosedGeomColl_r (NULL, geom);
2390 }
2391
2392 GAIAGEO_DECLARE int
gaiaIsNotClosedGeomColl_r(const void * cache,gaiaGeomCollPtr geom)2393 gaiaIsNotClosedGeomColl_r (const void *cache, gaiaGeomCollPtr geom)
2394 {
2395 /*
2396 / identifying not properly closed Rings
2397 / i.e. geoms making GEOS to crash !!!
2398 */
2399 int ret;
2400 int ib;
2401 gaiaPolygonPtr polyg;
2402 gaiaRingPtr ring;
2403 if (!geom)
2404 return 0;
2405 polyg = geom->FirstPolygon;
2406 while (polyg)
2407 {
2408 /* checking POLYGONs */
2409 ring = polyg->Exterior;
2410 if (cache != NULL)
2411 ret = gaiaIsNotClosedRing_r (cache, ring);
2412 else
2413 ret = gaiaIsNotClosedRing (ring);
2414 if (ret)
2415 return 1;
2416 for (ib = 0; ib < polyg->NumInteriors; ib++)
2417 {
2418 ring = polyg->Interiors + ib;
2419 if (cache != NULL)
2420 ret = gaiaIsNotClosedRing_r (cache, ring);
2421 else
2422 ret = gaiaIsNotClosedRing (ring);
2423 if (ret)
2424 return 1;
2425 }
2426 polyg = polyg->Next;
2427 }
2428 return 0;
2429 }
2430
2431 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaLinearize(gaiaGeomCollPtr geom,int force_multi)2432 gaiaLinearize (gaiaGeomCollPtr geom, int force_multi)
2433 {
2434 /* attempts to rearrange a generic Geometry into a (multi)linestring */
2435 int pts = 0;
2436 int lns = 0;
2437 gaiaGeomCollPtr result;
2438 gaiaPointPtr pt;
2439 gaiaLinestringPtr ln;
2440 gaiaLinestringPtr new_ln;
2441 gaiaPolygonPtr pg;
2442 gaiaRingPtr rng;
2443 int iv;
2444 int ib;
2445 double x;
2446 double y;
2447 double m;
2448 double z;
2449 if (!geom)
2450 return NULL;
2451 pt = geom->FirstPoint;
2452 while (pt)
2453 {
2454 pts++;
2455 pt = pt->Next;
2456 }
2457 ln = geom->FirstLinestring;
2458 while (ln)
2459 {
2460 lns++;
2461 ln = ln->Next;
2462 }
2463 if (pts || lns)
2464 return NULL;
2465
2466 if (geom->DimensionModel == GAIA_XY_Z_M)
2467 result = gaiaAllocGeomCollXYZM ();
2468 else if (geom->DimensionModel == GAIA_XY_Z)
2469 result = gaiaAllocGeomCollXYZ ();
2470 else if (geom->DimensionModel == GAIA_XY_M)
2471 result = gaiaAllocGeomCollXYM ();
2472 else
2473 result = gaiaAllocGeomColl ();
2474 result->Srid = geom->Srid;
2475 if (force_multi)
2476 result->DeclaredType = GAIA_MULTILINESTRING;
2477
2478 pg = geom->FirstPolygon;
2479 while (pg)
2480 {
2481 /* dissolving any POLYGON as simple LINESTRINGs (rings) */
2482 rng = pg->Exterior;
2483 new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
2484 for (iv = 0; iv < rng->Points; iv++)
2485 {
2486 /* copying the EXTERIOR RING as LINESTRING */
2487 if (geom->DimensionModel == GAIA_XY_Z_M)
2488 {
2489 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2490 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
2491 }
2492 else if (geom->DimensionModel == GAIA_XY_Z)
2493 {
2494 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2495 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
2496 }
2497 else if (geom->DimensionModel == GAIA_XY_M)
2498 {
2499 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2500 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
2501 }
2502 else
2503 {
2504 gaiaGetPoint (rng->Coords, iv, &x, &y);
2505 gaiaSetPoint (new_ln->Coords, iv, x, y);
2506 }
2507 }
2508 for (ib = 0; ib < pg->NumInteriors; ib++)
2509 {
2510 rng = pg->Interiors + ib;
2511 new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
2512 for (iv = 0; iv < rng->Points; iv++)
2513 {
2514 /* copying an INTERIOR RING as LINESTRING */
2515 if (geom->DimensionModel == GAIA_XY_Z_M)
2516 {
2517 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2518 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
2519 }
2520 else if (geom->DimensionModel == GAIA_XY_Z)
2521 {
2522 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2523 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
2524 }
2525 else if (geom->DimensionModel == GAIA_XY_M)
2526 {
2527 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2528 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
2529 }
2530 else
2531 {
2532 gaiaGetPoint (rng->Coords, iv, &x, &y);
2533 gaiaSetPoint (new_ln->Coords, iv, x, y);
2534 }
2535 }
2536 }
2537 pg = pg->Next;
2538 }
2539 if (result->FirstLinestring == NULL)
2540 {
2541 gaiaFreeGeomColl (result);
2542 return NULL;
2543 }
2544 return result;
2545 }
2546
2547 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDissolveSegments(gaiaGeomCollPtr geom)2548 gaiaDissolveSegments (gaiaGeomCollPtr geom)
2549 {
2550 /* attempts to dissolve a Geometry into segments */
2551 gaiaGeomCollPtr result;
2552 gaiaPointPtr pt;
2553 gaiaLinestringPtr ln;
2554 gaiaLinestringPtr segment;
2555 gaiaPolygonPtr pg;
2556 gaiaRingPtr rng;
2557 int iv;
2558 int ib;
2559 double x = 0.0;
2560 double y = 0.0;
2561 double z = 0.0;
2562 double m = 0.0;
2563 double x0 = 0.0;
2564 double y0 = 0.0;
2565 double z0 = 0.0;
2566 double m0 = 0.0;
2567 if (!geom)
2568 return NULL;
2569
2570 if (geom->DimensionModel == GAIA_XY_Z_M)
2571 result = gaiaAllocGeomCollXYZM ();
2572 else if (geom->DimensionModel == GAIA_XY_Z)
2573 result = gaiaAllocGeomCollXYZ ();
2574 else if (geom->DimensionModel == GAIA_XY_M)
2575 result = gaiaAllocGeomCollXYM ();
2576 else
2577 result = gaiaAllocGeomColl ();
2578 pt = geom->FirstPoint;
2579 while (pt)
2580 {
2581 if (geom->DimensionModel == GAIA_XY_Z_M)
2582 gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z, pt->M);
2583 else if (geom->DimensionModel == GAIA_XY_Z)
2584 gaiaAddPointToGeomCollXYZ (result, pt->X, pt->Y, pt->Z);
2585 else if (geom->DimensionModel == GAIA_XY_M)
2586 gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y, pt->M);
2587 else
2588 gaiaAddPointToGeomColl (result, pt->X, pt->Y);
2589 pt = pt->Next;
2590 }
2591 ln = geom->FirstLinestring;
2592 while (ln)
2593 {
2594 for (iv = 0; iv < ln->Points; iv++)
2595 {
2596 if (ln->DimensionModel == GAIA_XY_Z_M)
2597 {
2598 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2599 }
2600 else if (ln->DimensionModel == GAIA_XY_Z)
2601 {
2602 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2603 }
2604 else if (ln->DimensionModel == GAIA_XY_M)
2605 {
2606 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2607 }
2608 else
2609 {
2610 gaiaGetPoint (ln->Coords, iv, &x, &y);
2611 }
2612 if (iv > 0)
2613 {
2614 if (geom->DimensionModel == GAIA_XY_Z_M)
2615 {
2616 if (x != x0 || y != y0 || z != z0 || m != m0)
2617 {
2618 segment =
2619 gaiaAddLinestringToGeomColl (result, 2);
2620 gaiaSetPointXYZM (segment->Coords, 0, x0, y0,
2621 z0, m0);
2622 gaiaSetPointXYZM (segment->Coords, 1, x, y, z,
2623 m);
2624 }
2625 }
2626 else if (geom->DimensionModel == GAIA_XY_Z)
2627 {
2628 if (x != x0 || y != y0 || z != z0)
2629 {
2630 segment =
2631 gaiaAddLinestringToGeomColl (result, 2);
2632 gaiaSetPointXYZ (segment->Coords, 0, x0, y0,
2633 z0);
2634 gaiaSetPointXYZ (segment->Coords, 1, x, y, z);
2635 }
2636 }
2637 else if (geom->DimensionModel == GAIA_XY_M)
2638 {
2639 if (x != x0 || y != y0 || m != m0)
2640 {
2641 segment =
2642 gaiaAddLinestringToGeomColl (result, 2);
2643 gaiaSetPointXYM (segment->Coords, 0, x0, y0,
2644 m0);
2645 gaiaSetPointXYM (segment->Coords, 1, x, y, m);
2646 }
2647 }
2648 else
2649 {
2650 if (x != x0 || y != y0)
2651 {
2652 segment =
2653 gaiaAddLinestringToGeomColl (result, 2);
2654 gaiaSetPoint (segment->Coords, 0, x0, y0);
2655 gaiaSetPoint (segment->Coords, 1, x, y);
2656 }
2657 }
2658 }
2659 x0 = x;
2660 y0 = y;
2661 z0 = z;
2662 m0 = m;
2663 }
2664 ln = ln->Next;
2665 }
2666 pg = geom->FirstPolygon;
2667 while (pg)
2668 {
2669 rng = pg->Exterior;
2670 for (iv = 0; iv < rng->Points; iv++)
2671 {
2672 /* exterior Ring */
2673 if (rng->DimensionModel == GAIA_XY_Z_M)
2674 {
2675 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2676 }
2677 else if (rng->DimensionModel == GAIA_XY_Z)
2678 {
2679 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2680 }
2681 else if (rng->DimensionModel == GAIA_XY_M)
2682 {
2683 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2684 }
2685 else
2686 {
2687 gaiaGetPoint (rng->Coords, iv, &x, &y);
2688 }
2689 if (iv > 0)
2690 {
2691 if (geom->DimensionModel == GAIA_XY_Z_M)
2692 {
2693 if (x != x0 || y != y0 || z != z0 || m != m0)
2694 {
2695 segment =
2696 gaiaAddLinestringToGeomColl (result, 2);
2697 gaiaSetPointXYZM (segment->Coords, 0, x0, y0,
2698 z0, m0);
2699 gaiaSetPointXYZM (segment->Coords, 1, x, y, z,
2700 m);
2701 }
2702 }
2703 else if (geom->DimensionModel == GAIA_XY_Z)
2704 {
2705 if (x != x0 || y != y0 || z != z0)
2706 {
2707 segment =
2708 gaiaAddLinestringToGeomColl (result, 2);
2709 gaiaSetPointXYZ (segment->Coords, 0, x0, y0,
2710 z0);
2711 gaiaSetPointXYZ (segment->Coords, 1, x, y, z);
2712 }
2713 }
2714 else if (geom->DimensionModel == GAIA_XY_M)
2715 {
2716 if (x != x0 || y != y0 || m != m0)
2717 {
2718 segment =
2719 gaiaAddLinestringToGeomColl (result, 2);
2720 gaiaSetPointXYM (segment->Coords, 0, x0, y0,
2721 m0);
2722 gaiaSetPointXYM (segment->Coords, 1, x, y, m);
2723 }
2724 }
2725 else
2726 {
2727 if (x != x0 || y != y0)
2728 {
2729 segment =
2730 gaiaAddLinestringToGeomColl (result, 2);
2731 gaiaSetPoint (segment->Coords, 0, x0, y0);
2732 gaiaSetPoint (segment->Coords, 1, x, y);
2733 }
2734 }
2735 }
2736 x0 = x;
2737 y0 = y;
2738 z0 = z;
2739 m0 = m;
2740 }
2741 for (ib = 0; ib < pg->NumInteriors; ib++)
2742 {
2743 rng = pg->Interiors + ib;
2744 for (iv = 0; iv < rng->Points; iv++)
2745 {
2746 /* interior Ring */
2747 if (rng->DimensionModel == GAIA_XY_Z_M)
2748 {
2749 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2750 }
2751 else if (rng->DimensionModel == GAIA_XY_Z)
2752 {
2753 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2754 }
2755 else if (rng->DimensionModel == GAIA_XY_M)
2756 {
2757 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2758 }
2759 else
2760 {
2761 gaiaGetPoint (rng->Coords, iv, &x, &y);
2762 }
2763 if (iv > 0)
2764 {
2765 if (geom->DimensionModel == GAIA_XY_Z_M)
2766 {
2767 if (x != x0 || y != y0 || z != z0 || m != m0)
2768 {
2769 segment =
2770 gaiaAddLinestringToGeomColl (result,
2771 2);
2772 gaiaSetPointXYZM (segment->Coords, 0,
2773 x0, y0, z0, m0);
2774 gaiaSetPointXYZM (segment->Coords, 1, x,
2775 y, z, m);
2776 }
2777 }
2778 else if (geom->DimensionModel == GAIA_XY_Z)
2779 {
2780 if (x != x0 || y != y0 || z != z0)
2781 {
2782 segment =
2783 gaiaAddLinestringToGeomColl (result,
2784 2);
2785 gaiaSetPointXYZ (segment->Coords, 0, x0,
2786 y0, z0);
2787 gaiaSetPointXYZ (segment->Coords, 1, x,
2788 y, z);
2789 }
2790 }
2791 else if (geom->DimensionModel == GAIA_XY_M)
2792 {
2793 if (x != x0 || y != y0 || m != m0)
2794 {
2795 segment =
2796 gaiaAddLinestringToGeomColl (result,
2797 2);
2798 gaiaSetPointXYM (segment->Coords, 0, x0,
2799 y0, m0);
2800 gaiaSetPointXYM (segment->Coords, 1, x,
2801 y, m);
2802 }
2803 }
2804 else
2805 {
2806 if (x != x0 || y != y0)
2807 {
2808 segment =
2809 gaiaAddLinestringToGeomColl (result,
2810 2);
2811 gaiaSetPoint (segment->Coords, 0, x0,
2812 y0);
2813 gaiaSetPoint (segment->Coords, 1, x, y);
2814 }
2815 }
2816 }
2817 x0 = x;
2818 y0 = y;
2819 z0 = z;
2820 m0 = m;
2821 }
2822 }
2823 pg = pg->Next;
2824 }
2825 result->Srid = geom->Srid;
2826 return result;
2827 }
2828
2829 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDissolvePoints(gaiaGeomCollPtr geom)2830 gaiaDissolvePoints (gaiaGeomCollPtr geom)
2831 {
2832 /* attempts to dissolve a Geometry into points */
2833 gaiaGeomCollPtr result;
2834 gaiaPointPtr pt;
2835 gaiaLinestringPtr ln;
2836 gaiaPolygonPtr pg;
2837 gaiaRingPtr rng;
2838 int iv;
2839 int ib;
2840 double x = 0.0;
2841 double y = 0.0;
2842 double z = 0.0;
2843 double m = 0.0;
2844 if (!geom)
2845 return NULL;
2846
2847 if (geom->DimensionModel == GAIA_XY_Z_M)
2848 result = gaiaAllocGeomCollXYZM ();
2849 else if (geom->DimensionModel == GAIA_XY_Z)
2850 result = gaiaAllocGeomCollXYZ ();
2851 else if (geom->DimensionModel == GAIA_XY_M)
2852 result = gaiaAllocGeomCollXYM ();
2853 else
2854 result = gaiaAllocGeomColl ();
2855 pt = geom->FirstPoint;
2856 while (pt)
2857 {
2858 if (geom->DimensionModel == GAIA_XY_Z_M)
2859 gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z, pt->M);
2860 else if (geom->DimensionModel == GAIA_XY_Z)
2861 gaiaAddPointToGeomCollXYZ (result, pt->X, pt->Y, pt->Z);
2862 else if (geom->DimensionModel == GAIA_XY_M)
2863 gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y, pt->M);
2864 else
2865 gaiaAddPointToGeomColl (result, pt->X, pt->Y);
2866 pt = pt->Next;
2867 }
2868 ln = geom->FirstLinestring;
2869 while (ln)
2870 {
2871 for (iv = 0; iv < ln->Points; iv++)
2872 {
2873 if (ln->DimensionModel == GAIA_XY_Z_M)
2874 {
2875 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2876 }
2877 else if (ln->DimensionModel == GAIA_XY_Z)
2878 {
2879 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2880 }
2881 else if (ln->DimensionModel == GAIA_XY_M)
2882 {
2883 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2884 }
2885 else
2886 {
2887 gaiaGetPoint (ln->Coords, iv, &x, &y);
2888 }
2889 if (geom->DimensionModel == GAIA_XY_Z_M)
2890 gaiaAddPointToGeomCollXYZM (result, x, y, z, m);
2891 else if (geom->DimensionModel == GAIA_XY_Z)
2892 gaiaAddPointToGeomCollXYZ (result, x, y, z);
2893 else if (geom->DimensionModel == GAIA_XY_M)
2894 gaiaAddPointToGeomCollXYM (result, x, y, m);
2895 else
2896 gaiaAddPointToGeomColl (result, x, y);
2897 }
2898 ln = ln->Next;
2899 }
2900 pg = geom->FirstPolygon;
2901 while (pg)
2902 {
2903 rng = pg->Exterior;
2904 for (iv = 0; iv < rng->Points; iv++)
2905 {
2906 /* exterior Ring */
2907 if (rng->DimensionModel == GAIA_XY_Z_M)
2908 {
2909 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2910 }
2911 else if (rng->DimensionModel == GAIA_XY_Z)
2912 {
2913 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2914 }
2915 else if (rng->DimensionModel == GAIA_XY_M)
2916 {
2917 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2918 }
2919 else
2920 {
2921 gaiaGetPoint (rng->Coords, iv, &x, &y);
2922 }
2923 if (geom->DimensionModel == GAIA_XY_Z_M)
2924 gaiaAddPointToGeomCollXYZM (result, x, y, z, m);
2925 else if (geom->DimensionModel == GAIA_XY_Z)
2926 gaiaAddPointToGeomCollXYZ (result, x, y, z);
2927 else if (geom->DimensionModel == GAIA_XY_M)
2928 gaiaAddPointToGeomCollXYM (result, x, y, m);
2929 else
2930 gaiaAddPointToGeomColl (result, x, y);
2931 }
2932 for (ib = 0; ib < pg->NumInteriors; ib++)
2933 {
2934 rng = pg->Interiors + ib;
2935 for (iv = 0; iv < rng->Points; iv++)
2936 {
2937 /* interior Ring */
2938 if (rng->DimensionModel == GAIA_XY_Z_M)
2939 {
2940 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2941 }
2942 else if (rng->DimensionModel == GAIA_XY_Z)
2943 {
2944 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2945 }
2946 else if (rng->DimensionModel == GAIA_XY_M)
2947 {
2948 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2949 }
2950 else
2951 {
2952 gaiaGetPoint (rng->Coords, iv, &x, &y);
2953 }
2954 if (geom->DimensionModel == GAIA_XY_Z_M)
2955 gaiaAddPointToGeomCollXYZM (result, x, y, z, m);
2956 else if (geom->DimensionModel == GAIA_XY_Z)
2957 gaiaAddPointToGeomCollXYZ (result, x, y, z);
2958 else if (geom->DimensionModel == GAIA_XY_M)
2959 gaiaAddPointToGeomCollXYM (result, x, y, m);
2960 else
2961 gaiaAddPointToGeomColl (result, x, y);
2962 }
2963 }
2964 pg = pg->Next;
2965 }
2966 result->Srid = geom->Srid;
2967 return result;
2968 }
2969
2970 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaExtractPointsFromGeomColl(gaiaGeomCollPtr geom)2971 gaiaExtractPointsFromGeomColl (gaiaGeomCollPtr geom)
2972 {
2973 /* extracts any POINT from a GeometryCollection */
2974 gaiaGeomCollPtr result;
2975 gaiaPointPtr pt;
2976 int pts = 0;
2977 if (!geom)
2978 return NULL;
2979
2980 pt = geom->FirstPoint;
2981 while (pt)
2982 {
2983 pts++;
2984 pt = pt->Next;
2985 }
2986 if (!pts)
2987 return NULL;
2988
2989 if (geom->DimensionModel == GAIA_XY_Z_M)
2990 result = gaiaAllocGeomCollXYZM ();
2991 else if (geom->DimensionModel == GAIA_XY_Z)
2992 result = gaiaAllocGeomCollXYZ ();
2993 else if (geom->DimensionModel == GAIA_XY_M)
2994 result = gaiaAllocGeomCollXYM ();
2995 else
2996 result = gaiaAllocGeomColl ();
2997 pt = geom->FirstPoint;
2998 while (pt)
2999 {
3000 if (geom->DimensionModel == GAIA_XY_Z_M)
3001 gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z, pt->M);
3002 else if (geom->DimensionModel == GAIA_XY_Z)
3003 gaiaAddPointToGeomCollXYZ (result, pt->X, pt->Y, pt->Z);
3004 else if (geom->DimensionModel == GAIA_XY_M)
3005 gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y, pt->M);
3006 else
3007 gaiaAddPointToGeomColl (result, pt->X, pt->Y);
3008 pt = pt->Next;
3009 }
3010 result->Srid = geom->Srid;
3011 if (pts == 1)
3012 result->DeclaredType = GAIA_POINT;
3013 else
3014 result->DeclaredType = GAIA_MULTIPOINT;
3015 return result;
3016 }
3017
3018 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaExtractLinestringsFromGeomColl(gaiaGeomCollPtr geom)3019 gaiaExtractLinestringsFromGeomColl (gaiaGeomCollPtr geom)
3020 {
3021 /* extracts any LINESTRING from a GeometryCollection */
3022 gaiaGeomCollPtr result;
3023 gaiaLinestringPtr ln;
3024 gaiaLinestringPtr new_ln;
3025 int lns = 0;
3026 int iv;
3027 double x;
3028 double y;
3029 double z;
3030 double m;
3031 if (!geom)
3032 return NULL;
3033
3034 ln = geom->FirstLinestring;
3035 while (ln)
3036 {
3037 lns++;
3038 ln = ln->Next;
3039 }
3040 if (!lns)
3041 return NULL;
3042
3043 if (geom->DimensionModel == GAIA_XY_Z_M)
3044 result = gaiaAllocGeomCollXYZM ();
3045 else if (geom->DimensionModel == GAIA_XY_Z)
3046 result = gaiaAllocGeomCollXYZ ();
3047 else if (geom->DimensionModel == GAIA_XY_M)
3048 result = gaiaAllocGeomCollXYM ();
3049 else
3050 result = gaiaAllocGeomColl ();
3051 ln = geom->FirstLinestring;
3052 while (ln)
3053 {
3054 new_ln = gaiaAddLinestringToGeomColl (result, ln->Points);
3055 for (iv = 0; iv < ln->Points; iv++)
3056 {
3057 if (ln->DimensionModel == GAIA_XY_Z_M)
3058 {
3059 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3060 gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
3061 }
3062 else if (ln->DimensionModel == GAIA_XY_Z)
3063 {
3064 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3065 gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
3066 }
3067 else if (ln->DimensionModel == GAIA_XY_M)
3068 {
3069 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3070 gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
3071 }
3072 else
3073 {
3074 gaiaGetPoint (ln->Coords, iv, &x, &y);
3075 gaiaSetPoint (new_ln->Coords, iv, x, y);
3076 }
3077 }
3078 ln = ln->Next;
3079 }
3080 result->Srid = geom->Srid;
3081 if (lns == 1)
3082 result->DeclaredType = GAIA_LINESTRING;
3083 else
3084 result->DeclaredType = GAIA_MULTILINESTRING;
3085 return result;
3086 }
3087
3088 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaExtractPolygonsFromGeomColl(gaiaGeomCollPtr geom)3089 gaiaExtractPolygonsFromGeomColl (gaiaGeomCollPtr geom)
3090 {
3091 /* extracts any POLYGON from a GeometryCollection */
3092 gaiaGeomCollPtr result;
3093 gaiaPolygonPtr pg;
3094 gaiaPolygonPtr new_pg;
3095 gaiaRingPtr rng;
3096 gaiaRingPtr new_rng;
3097 int pgs = 0;
3098 int iv;
3099 int ib;
3100 double x;
3101 double y;
3102 double z;
3103 double m;
3104 if (!geom)
3105 return NULL;
3106
3107 pg = geom->FirstPolygon;
3108 while (pg)
3109 {
3110 pgs++;
3111 pg = pg->Next;
3112 }
3113 if (!pgs)
3114 return NULL;
3115
3116 if (geom->DimensionModel == GAIA_XY_Z_M)
3117 result = gaiaAllocGeomCollXYZM ();
3118 else if (geom->DimensionModel == GAIA_XY_Z)
3119 result = gaiaAllocGeomCollXYZ ();
3120 else if (geom->DimensionModel == GAIA_XY_M)
3121 result = gaiaAllocGeomCollXYM ();
3122 else
3123 result = gaiaAllocGeomColl ();
3124 pg = geom->FirstPolygon;
3125 while (pg)
3126 {
3127 rng = pg->Exterior;
3128 new_pg =
3129 gaiaAddPolygonToGeomColl (result, rng->Points, pg->NumInteriors);
3130 new_rng = new_pg->Exterior;
3131 for (iv = 0; iv < rng->Points; iv++)
3132 {
3133 if (rng->DimensionModel == GAIA_XY_Z_M)
3134 {
3135 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3136 gaiaSetPointXYZM (new_rng->Coords, iv, x, y, z, m);
3137 }
3138 else if (rng->DimensionModel == GAIA_XY_Z)
3139 {
3140 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3141 gaiaSetPointXYZ (new_rng->Coords, iv, x, y, z);
3142 }
3143 else if (rng->DimensionModel == GAIA_XY_M)
3144 {
3145 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3146 gaiaSetPointXYM (new_rng->Coords, iv, x, y, m);
3147 }
3148 else
3149 {
3150 gaiaGetPoint (rng->Coords, iv, &x, &y);
3151 gaiaSetPoint (new_rng->Coords, iv, x, y);
3152 }
3153 }
3154 for (ib = 0; ib < pg->NumInteriors; ib++)
3155 {
3156 rng = pg->Interiors + ib;
3157 new_rng = gaiaAddInteriorRing (new_pg, ib, rng->Points);
3158 for (iv = 0; iv < rng->Points; iv++)
3159 {
3160 if (rng->DimensionModel == GAIA_XY_Z_M)
3161 {
3162 gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3163 gaiaSetPointXYZM (new_rng->Coords, iv, x, y, z, m);
3164 }
3165 else if (rng->DimensionModel == GAIA_XY_Z)
3166 {
3167 gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3168 gaiaSetPointXYZ (new_rng->Coords, iv, x, y, z);
3169 }
3170 else if (rng->DimensionModel == GAIA_XY_M)
3171 {
3172 gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3173 gaiaSetPointXYM (new_rng->Coords, iv, x, y, m);
3174 }
3175 else
3176 {
3177 gaiaGetPoint (rng->Coords, iv, &x, &y);
3178 gaiaSetPoint (new_rng->Coords, iv, x, y);
3179 }
3180 }
3181 }
3182 pg = pg->Next;
3183 }
3184 result->Srid = geom->Srid;
3185 if (pgs == 1)
3186 result->DeclaredType = GAIA_POLYGON;
3187 else
3188 result->DeclaredType = GAIA_MULTIPOLYGON;
3189 return result;
3190 }
3191
3192 SPATIALITE_PRIVATE int
gaia_do_check_linestring(const void * g)3193 gaia_do_check_linestring (const void *g)
3194 {
3195 /* testing if the Geometry is a simple Linestring */
3196 gaiaGeomCollPtr geom = (gaiaGeomCollPtr) g;
3197 gaiaPointPtr pt;
3198 gaiaLinestringPtr ln;
3199 gaiaPolygonPtr pg;
3200 int pts = 0;
3201 int lns = 0;
3202 int pgs = 0;
3203 pt = geom->FirstPoint;
3204 while (pt != NULL)
3205 {
3206 /* counting how many Points are there */
3207 pts++;
3208 pt = pt->Next;
3209 }
3210 ln = geom->FirstLinestring;
3211 while (ln != NULL)
3212 {
3213 /* counting how many Linestrings are there */
3214 lns++;
3215 ln = ln->Next;
3216 }
3217 pg = geom->FirstPolygon;
3218 while (pg != NULL)
3219 {
3220 /* counting how many Polygons are there */
3221 pgs++;
3222 pg = pg->Next;
3223 }
3224 if (pts == 0 && lns == 1 && pgs == 0)
3225 return 1;
3226 return 0;
3227 }
3228
3229 static int
do_create_points(sqlite3 * mem_db,const char * table)3230 do_create_points (sqlite3 * mem_db, const char *table)
3231 {
3232 /* creating a table for storing Points */
3233 int ret;
3234 char *sql;
3235 char *err_msg = NULL;
3236
3237 /* creating the main table */
3238 if (strcmp (table, "points1") == 0)
3239 sql = sqlite3_mprintf ("CREATE TABLE %s "
3240 "(id INTEGER PRIMARY KEY AUTOINCREMENT, "
3241 "geom BLOB NOT NULL, "
3242 "needs_interpolation INTEGER NOT NULL)", table);
3243 else
3244 sql = sqlite3_mprintf ("CREATE TABLE %s "
3245 "(id INTEGER PRIMARY KEY AUTOINCREMENT, "
3246 "geom BLOB NOT NULL)", table);
3247 ret = sqlite3_exec (mem_db, sql, NULL, NULL, &err_msg);
3248 sqlite3_free (sql);
3249 if (ret != SQLITE_OK)
3250 {
3251 spatialite_e ("gaiaDrapeLine: CREATE TABLE \"%s\" error: %s\n",
3252 table, err_msg);
3253 sqlite3_free (err_msg);
3254 return 0;
3255 }
3256
3257 if (strcmp (table, "points1") == 0)
3258 return 1;
3259
3260 /* creating the companion R*Tree table */
3261 sql = sqlite3_mprintf ("CREATE VIRTUAL TABLE rtree_%s "
3262 "USING rtree(pkid, xmin, xmax, ymin, ymax)", table);
3263 ret = sqlite3_exec (mem_db, sql, NULL, NULL, &err_msg);
3264 sqlite3_free (sql);
3265 if (ret != SQLITE_OK)
3266 {
3267 spatialite_e ("gaiaDrapeLine: CREATE TABLE \"rtree_%s\" error: %s\n",
3268 table, err_msg);
3269 sqlite3_free (err_msg);
3270 return 0;
3271 }
3272
3273 return 1;
3274 }
3275
3276 static int
do_insert_point(sqlite3 * mem_db,sqlite3_stmt * stmt_pts,sqlite3_stmt * stmt_rtree_pts,double x,double y,double z,double m)3277 do_insert_point (sqlite3 * mem_db, sqlite3_stmt * stmt_pts,
3278 sqlite3_stmt * stmt_rtree_pts, double x,
3279 double y, double z, double m)
3280 {
3281 /* inserting into the Points2 helper table */
3282 int ret;
3283 sqlite3_int64 rowid;
3284
3285 sqlite3_reset (stmt_pts);
3286 sqlite3_clear_bindings (stmt_pts);
3287 sqlite3_bind_double (stmt_pts, 1, x);
3288 sqlite3_bind_double (stmt_pts, 2, y);
3289 sqlite3_bind_double (stmt_pts, 3, z);
3290 sqlite3_bind_double (stmt_pts, 4, m);
3291 ret = sqlite3_step (stmt_pts);
3292 if (ret == SQLITE_DONE || ret == SQLITE_ROW)
3293 rowid = sqlite3_last_insert_rowid (mem_db);
3294 else
3295 {
3296 spatialite_e ("INSERT INTO \"Points\" error: \"%s\"\n",
3297 sqlite3_errmsg (mem_db));
3298 goto error;
3299 }
3300
3301 /* inserting into the companion R*Tree */
3302 sqlite3_reset (stmt_rtree_pts);
3303 sqlite3_clear_bindings (stmt_rtree_pts);
3304 sqlite3_bind_int64 (stmt_rtree_pts, 1, rowid);
3305 sqlite3_bind_double (stmt_rtree_pts, 2, x);
3306 sqlite3_bind_double (stmt_rtree_pts, 3, x);
3307 sqlite3_bind_double (stmt_rtree_pts, 4, y);
3308 sqlite3_bind_double (stmt_rtree_pts, 5, y);
3309 ret = sqlite3_step (stmt_rtree_pts);
3310 if (ret == SQLITE_DONE || ret == SQLITE_ROW)
3311 ;
3312 else
3313 {
3314 spatialite_e ("INSERT INTO \"RTree_Points\" error: \"%s\"\n",
3315 sqlite3_errmsg (mem_db));
3316 goto error;
3317 }
3318
3319 return 1;
3320
3321 error:
3322 return 0;
3323 }
3324
3325 static int
do_populate_points2(sqlite3 * mem_db,gaiaGeomCollPtr geom)3326 do_populate_points2 (sqlite3 * mem_db, gaiaGeomCollPtr geom)
3327 {
3328 /* populating Points-2 */
3329 int ret;
3330 const char *sql;
3331 int iv;
3332 gaiaLinestringPtr ln;
3333 sqlite3_stmt *stmt_pts = NULL;
3334 sqlite3_stmt *stmt_rtree_pts = NULL;
3335 double ox;
3336 double oy;
3337 double oz;
3338 double om;
3339 double fx;
3340 double fy;
3341 double fz;
3342 double fm;
3343
3344 /* creating an SQL statement for inserting rows into the Points2 table */
3345 sql = "INSERT INTO points2 (id, geom) VALUES "
3346 "(NULL, MakePointZM(?, ?, ?, ?))";
3347 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt_pts, NULL);
3348 if (ret != SQLITE_OK)
3349 {
3350 spatialite_e ("INSERT INTO Points2: error %d \"%s\"\n",
3351 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3352 goto error;
3353 }
3354
3355 /* creating an SQL statement for inserting rows into the R*TREE table */
3356 sql = "INSERT INTO rtree_points2 (pkid, xmin, xmax, ymin, ymax) "
3357 "VALUES (?, ?, ?, ?, ?)";
3358 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt_rtree_pts, NULL);
3359 if (ret != SQLITE_OK)
3360 {
3361 spatialite_e ("INSERT INTO RTree_Points2: error %d \"%s\"\n",
3362 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3363 goto error;
3364 }
3365
3366 /* starting a Transaction */
3367 sql = "BEGIN";
3368 ret = sqlite3_exec (mem_db, sql, NULL, NULL, NULL);
3369 if (ret != SQLITE_OK)
3370 {
3371 spatialite_e ("BEGIN: error: %d \"%s\"\n",
3372 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3373 goto error;
3374 }
3375
3376 ln = geom->FirstLinestring;
3377 for (iv = 0; iv < ln->Points; iv++)
3378 {
3379 /* processing all Vertices from the input Linestring */
3380 double x;
3381 double y;
3382 double z = 0.0;
3383 double m = 0.0;
3384 int repeated;
3385 if (ln->DimensionModel == GAIA_XY_Z)
3386 {
3387 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3388 }
3389 else if (ln->DimensionModel == GAIA_XY_M)
3390 {
3391 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3392 }
3393 else if (ln->DimensionModel == GAIA_XY_Z_M)
3394 {
3395 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3396 }
3397 else
3398 {
3399 gaiaGetPoint (ln->Coords, iv, &x, &y);
3400 }
3401 repeated = 0;
3402 if (iv != 0)
3403 {
3404 /* checking for repeated points */
3405 if (x == ox && y == oy && z == oz && m == om)
3406 repeated = 1;
3407 }
3408 if (iv == ln->Points - 1)
3409 {
3410 /* checkink for a closed Linestring */
3411 if (x == fx && y == fy && z == fz && m == fm)
3412 repeated = 1;
3413 }
3414 if (!repeated)
3415 {
3416 /* inserting a Point */
3417 if (!do_insert_point
3418 (mem_db, stmt_pts, stmt_rtree_pts, x, y, z, m))
3419 goto error;
3420 }
3421 /* saving the current coords */
3422 ox = x;
3423 oy = y;
3424 oz = z;
3425 om = m;
3426 if (iv == 0)
3427 {
3428 /* saving the Start Point coords */
3429 fx = x;
3430 fy = y;
3431 fz = z;
3432 fm = m;
3433 }
3434 }
3435
3436 /* committing the pending Transaction */
3437 sql = "COMMIT";
3438 ret = sqlite3_exec (mem_db, sql, NULL, NULL, NULL);
3439 if (ret != SQLITE_OK)
3440 {
3441 spatialite_e ("COMMIT: error: %d \"%s\"\n",
3442 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3443 goto error;
3444 }
3445
3446 sqlite3_finalize (stmt_pts);
3447 sqlite3_finalize (stmt_rtree_pts);
3448 return 1;
3449
3450 error:
3451 if (stmt_pts != NULL)
3452 sqlite3_finalize (stmt_pts);
3453 if (stmt_rtree_pts != NULL)
3454 sqlite3_finalize (stmt_rtree_pts);
3455 return 0;
3456 }
3457
3458 static int
do_insert_draped_point(sqlite3 * mem_db,sqlite3_stmt * stmt_pts,int needs_interpolation,gaiaGeomCollPtr geom)3459 do_insert_draped_point (sqlite3 * mem_db, sqlite3_stmt * stmt_pts,
3460 int needs_interpolation, gaiaGeomCollPtr geom)
3461 {
3462 /* inserting into the Points helper table */
3463 int ret;
3464 gaiaPointPtr pt = geom->FirstPoint;
3465
3466 if (pt == NULL)
3467 return 0;
3468
3469 sqlite3_reset (stmt_pts);
3470 sqlite3_clear_bindings (stmt_pts);
3471 sqlite3_bind_double (stmt_pts, 1, pt->X);
3472 sqlite3_bind_double (stmt_pts, 2, pt->Y);
3473 sqlite3_bind_double (stmt_pts, 3, pt->Z);
3474 sqlite3_bind_double (stmt_pts, 4, pt->M);
3475 sqlite3_bind_int (stmt_pts, 5, needs_interpolation);
3476 ret = sqlite3_step (stmt_pts);
3477 if (ret == SQLITE_DONE || ret == SQLITE_ROW)
3478 return 1;
3479
3480 spatialite_e ("INSERT INTO \"Points1\" error: \"%s\"\n",
3481 sqlite3_errmsg (mem_db));
3482 return 0;
3483 }
3484
3485 static gaiaGeomCollPtr
do_point_drape_coords(int srid,double x,double y,gaiaGeomCollPtr geom_3d)3486 do_point_drape_coords (int srid, double x, double y, gaiaGeomCollPtr geom_3d)
3487 {
3488 /* draping a Point */
3489 gaiaPointPtr pt_3d = geom_3d->FirstPoint;
3490 gaiaGeomCollPtr geom = gaiaAllocGeomCollXYZM ();
3491 geom->Srid = srid;
3492 gaiaAddPointToGeomCollXYZM (geom, x, y, pt_3d->Z, pt_3d->M);
3493 return geom;
3494 }
3495
3496 static gaiaGeomCollPtr
do_point_same_coords(int srid,double x,double y,double z,double m)3497 do_point_same_coords (int srid, double x, double y, double z, double m)
3498 {
3499 /* creating a Point (unchanged coords) */
3500 gaiaGeomCollPtr geom;
3501 geom = gaiaAllocGeomCollXYZM ();
3502 geom->Srid = srid;
3503 gaiaAddPointToGeomCollXYZM (geom, x, y, z, m);
3504 return geom;
3505 }
3506
3507 static int
do_drape_vertex(sqlite3 * mem_db,sqlite3_stmt * stmt,sqlite3_stmt * stmt_pts,int srid,double tolerance,double x,double y,double z,double m)3508 do_drape_vertex (sqlite3 * mem_db, sqlite3_stmt * stmt,
3509 sqlite3_stmt * stmt_pts, int srid, double tolerance,
3510 double x, double y, double z, double m)
3511 {
3512 /* draping the segment */
3513 double minx = x;
3514 double miny = y;
3515 double maxx = x;
3516 double maxy = y;
3517 int ret;
3518 int count = 0;
3519 gaiaGeomCollPtr g2;
3520
3521 /* preparing an extendend BBOX */
3522 minx = x - (tolerance * 2.0);
3523 miny = y - (tolerance * 2.0);
3524 maxx = x + (tolerance * 2.0);
3525 maxy = y + (tolerance * 2.0);
3526
3527 /* querying Points-2 by distance from Segment */
3528 sqlite3_reset (stmt);
3529 sqlite3_clear_bindings (stmt);
3530 sqlite3_bind_double (stmt, 1, minx); /* R*Tree BBOX */
3531 sqlite3_bind_double (stmt, 2, miny);
3532 sqlite3_bind_double (stmt, 3, maxx);
3533 sqlite3_bind_double (stmt, 4, maxy);
3534 sqlite3_bind_double (stmt, 5, x); /* point (distance) */
3535 sqlite3_bind_double (stmt, 6, y);
3536 sqlite3_bind_double (stmt, 7, tolerance);
3537 while (1)
3538 {
3539 /* scrolling the result set rows */
3540 ret = sqlite3_step (stmt);
3541 if (ret == SQLITE_DONE)
3542 break; /* end of result set */
3543 if (ret == SQLITE_ROW)
3544 {
3545 gaiaGeomCollPtr g2 = NULL;
3546 gaiaGeomCollPtr geom = NULL;
3547 if (sqlite3_column_type (stmt, 0) == SQLITE_BLOB)
3548 {
3549 unsigned char *p_blob =
3550 (unsigned char *) sqlite3_column_blob (stmt, 0);
3551 int n_bytes = sqlite3_column_bytes (stmt, 0);
3552 geom = gaiaFromSpatiaLiteBlobWkb (p_blob, n_bytes);
3553 }
3554 if (geom != NULL)
3555 {
3556 /* found a valid Point */
3557 g2 = do_point_drape_coords (srid, x, y, geom);
3558 gaiaFreeGeomColl (geom);
3559 if (!do_insert_draped_point (mem_db, stmt_pts, 0, g2))
3560 goto error;
3561 gaiaFreeGeomColl (g2);
3562 count++;
3563 }
3564 }
3565 }
3566 if (count == 0)
3567 {
3568 /* not found - inserting the Point itself */
3569 g2 = do_point_same_coords (srid, x, y, z, m);
3570 if (!do_insert_draped_point (mem_db, stmt_pts, 1, g2))
3571 goto error;
3572 gaiaFreeGeomColl (g2);
3573 }
3574 return 1;
3575
3576 error:
3577 return 0;
3578 }
3579
3580 static int
do_drape_line(sqlite3 * mem_db,gaiaGeomCollPtr geom,double tolerance)3581 do_drape_line (sqlite3 * mem_db, gaiaGeomCollPtr geom, double tolerance)
3582 {
3583 /* draping the line */
3584 int ret;
3585 sqlite3_stmt *stmt = NULL;
3586 sqlite3_stmt *stmt_pts = NULL;
3587 const char *sql;
3588 gaiaLinestringPtr ln;
3589 int iv;
3590
3591 /* creating an SQL statement for querying Points-2 */
3592 sql = "SELECT geom FROM points2 "
3593 "WHERE ROWID IN (SELECT pkid FROM rtree_points2 WHERE "
3594 "MbrIntersects(geom, BuildMbr(?, ?, ?, ?)) = 1) "
3595 "AND ST_Distance(geom, MakePoint(?, ?)) <= ? " "ORDER BY id";
3596 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt, NULL);
3597 if (ret != SQLITE_OK)
3598 {
3599 spatialite_e ("SELECT Points2: error %d \"%s\"\n",
3600 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3601 goto error;
3602 }
3603
3604 /* creating an SQL statement for inserting rows into the Points-1 table */
3605 sql = "INSERT INTO points1 (id, geom, needs_interpolation) VALUES "
3606 "(NULL, MakePointZM(?, ?, ?, ?), ?)";
3607 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt_pts, NULL);
3608 if (ret != SQLITE_OK)
3609 {
3610 spatialite_e ("INSERT INTO Points1: error %d \"%s\"\n",
3611 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3612 goto error;
3613 }
3614
3615 /* starting a Transaction */
3616 sql = "BEGIN";
3617 ret = sqlite3_exec (mem_db, sql, NULL, NULL, NULL);
3618 if (ret != SQLITE_OK)
3619 {
3620 spatialite_e ("BEGIN: error: %d \"%s\"\n",
3621 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3622 goto error;
3623 }
3624
3625 ln = geom->FirstLinestring;
3626 for (iv = 0; iv < ln->Points; iv++)
3627 {
3628 /* processing all Vertices from the input Linestring */
3629 double x;
3630 double y;
3631 double z = 0.0;
3632 double m = 0.0;
3633 if (ln->DimensionModel == GAIA_XY_Z)
3634 {
3635 gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3636 }
3637 else if (ln->DimensionModel == GAIA_XY_M)
3638 {
3639 gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3640 }
3641 else if (ln->DimensionModel == GAIA_XY_Z_M)
3642 {
3643 gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3644 }
3645 else
3646 {
3647 gaiaGetPoint (ln->Coords, iv, &x, &y);
3648 }
3649 /* processing a Vertex */
3650 if (!do_drape_vertex
3651 (mem_db, stmt, stmt_pts, geom->Srid, tolerance, x, y, z, m))
3652 goto error;
3653 }
3654
3655 /* committing the pending Transaction */
3656 sql = "COMMIT";
3657 ret = sqlite3_exec (mem_db, sql, NULL, NULL, NULL);
3658 if (ret != SQLITE_OK)
3659 {
3660 spatialite_e ("COMMIT: error: %d \"%s\"\n",
3661 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3662 goto error;
3663 }
3664
3665 sqlite3_finalize (stmt);
3666 sqlite3_finalize (stmt_pts);
3667 return 1;
3668
3669 error:
3670 if (stmt != NULL)
3671 sqlite3_finalize (stmt);
3672 if (stmt_pts != NULL)
3673 sqlite3_finalize (stmt_pts);
3674 return 0;
3675 }
3676
3677 static int
get_prev_coords(int index,gaiaDynamicLinePtr dyn,double * z,double * m,double * dist)3678 get_prev_coords (int index, gaiaDynamicLinePtr dyn, double *z, double *m,
3679 double *dist)
3680 {
3681 /* retrieving the previous Point */
3682 double ox;
3683 double oy;
3684 double x;
3685 double y;
3686 double zz;
3687 double mm;
3688 int ok = 0;
3689 int count = 0;
3690 gaiaPointPtr pt = dyn->First;
3691 while (pt != NULL)
3692 {
3693 if (index - 1 == count)
3694 {
3695 /* this is the previous point */
3696 ox = pt->X;
3697 oy = pt->Y;
3698 zz = pt->Z;
3699 mm = pt->M;
3700 ok = 1;
3701 }
3702 if (index == count)
3703 {
3704 /* this is the current point */
3705 x = pt->X;
3706 y = pt->Y;
3707 if (ok)
3708 {
3709 *z = zz;
3710 *m = mm;
3711 *dist =
3712 sqrt (((ox - x) * (ox - x)) + ((oy - y) * (oy - y)));
3713 return 1;
3714 }
3715 return 0;
3716 }
3717 count++;
3718 pt = pt->Next;
3719 }
3720 return 0;
3721 }
3722
3723 static int
get_next_coords(int index,gaiaDynamicLinePtr dyn,const char * interpolate,double * z,double * m,double * dist)3724 get_next_coords (int index, gaiaDynamicLinePtr dyn, const char *interpolate,
3725 double *z, double *m, double *dist)
3726 {
3727 /* retrieving the next Point */
3728 double ox;
3729 double oy;
3730 double x;
3731 double y;
3732 double d = 0.0;
3733 int ok = 0;
3734 int count = 0;
3735 gaiaPointPtr pt = dyn->First;
3736 while (pt != NULL)
3737 {
3738 if (index == count)
3739 {
3740 /* this is the current Point */
3741 ox = pt->X;
3742 oy = pt->Y;
3743 ok = 1;
3744 }
3745 if (index < count)
3746 {
3747 /* this is a following Point */
3748 if (!ok)
3749 return 0;
3750 x = pt->X;
3751 y = pt->Y;
3752 d += sqrt (((ox - x) * (ox - x)) + ((oy - y) * (oy - y)));
3753 if (*(interpolate + count) == 'N')
3754 {
3755 /* found a valid Point */
3756 *z = pt->Z;
3757 *m = pt->M;
3758 *dist = d;
3759 return 1;
3760 }
3761 }
3762 count++;
3763 pt = pt->Next;
3764 }
3765 return 0;
3766 }
3767
3768 static int
do_update_coords(int index,gaiaDynamicLinePtr dyn,double z,double m)3769 do_update_coords (int index, gaiaDynamicLinePtr dyn, double z, double m)
3770 {
3771 /* updating the coords */
3772 int count = 0;
3773 gaiaPointPtr pt = dyn->First;
3774 while (pt != NULL)
3775 {
3776 if (index == count)
3777 {
3778 pt->Z = z;
3779 pt->M = m;
3780 return 1;
3781 }
3782 count++;
3783 pt = pt->Next;
3784 }
3785 return 0;
3786 }
3787
3788 static void
do_interpolate_coords(int index,gaiaDynamicLinePtr dyn,char * interpolate)3789 do_interpolate_coords (int index, gaiaDynamicLinePtr dyn, char *interpolate)
3790 {
3791 /* attempting to interpolate coords */
3792 double pz;
3793 double pm;
3794 double nz;
3795 double nm;
3796 double z;
3797 double m;
3798 double pdist;
3799 double ndist;
3800 double perc;
3801 if (!get_prev_coords (index, dyn, &pz, &pm, &pdist))
3802 return;
3803 if (!get_next_coords (index, dyn, interpolate, &nz, &nm, &ndist))
3804 return;
3805 perc = pdist / (pdist + ndist);
3806 z = pz + ((nz - pz) * perc);
3807 m = pm + ((nm - pm) * perc);
3808 if (do_update_coords (index, dyn, z, m))
3809 *(interpolate + index) = 'I';
3810 }
3811
3812 static gaiaGeomCollPtr
do_reassemble_line(sqlite3 * mem_db,int dims,int srid)3813 do_reassemble_line (sqlite3 * mem_db, int dims, int srid)
3814 {
3815 /* reassembling the Linestring to be returned */
3816 gaiaGeomCollPtr geom = NULL;
3817 gaiaPointPtr pt;
3818 gaiaLinestringPtr ln;
3819 gaiaDynamicLinePtr dyn = gaiaAllocDynamicLine ();
3820 const char *sql;
3821 int ret;
3822 sqlite3_stmt *stmt = NULL;
3823 int count = 0;
3824 int needs_interpolation = 0;
3825
3826 /* creating an SQL statement for querying Points-1 */
3827 sql = "SELECT geom, needs_interpolation FROM points1 ORDER BY id";
3828 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt, NULL);
3829 if (ret != SQLITE_OK)
3830 {
3831 spatialite_e ("SELECT Points1: error %d \"%s\"\n",
3832 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3833 goto end;
3834 }
3835 while (1)
3836 {
3837 /* scrolling the result set rows */
3838 ret = sqlite3_step (stmt);
3839 if (ret == SQLITE_DONE)
3840 break; /* end of result set */
3841 if (ret == SQLITE_ROW)
3842 {
3843 gaiaGeomCollPtr g = NULL;
3844 if (sqlite3_column_type (stmt, 0) == SQLITE_BLOB)
3845 {
3846 unsigned char *p_blob =
3847 (unsigned char *) sqlite3_column_blob (stmt, 0);
3848 int n_bytes = sqlite3_column_bytes (stmt, 0);
3849 g = gaiaFromSpatiaLiteBlobWkb (p_blob, n_bytes);
3850 }
3851 if (g != NULL)
3852 {
3853 /* found a valid Point */
3854 pt = g->FirstPoint;
3855 if (dims == GAIA_XY_Z_M)
3856 gaiaAppendPointZMToDynamicLine (dyn, pt->X, pt->Y,
3857 pt->Z, pt->M);
3858 else if (dims == GAIA_XY_Z)
3859 gaiaAppendPointZToDynamicLine (dyn, pt->X, pt->Y,
3860 pt->Z);
3861 else if (dims == GAIA_XY_M)
3862 gaiaAppendPointMToDynamicLine (dyn, pt->X, pt->Y,
3863 pt->M);
3864 else
3865 gaiaAppendPointToDynamicLine (dyn, pt->X, pt->Y);
3866 gaiaFreeGeomColl (g);
3867 }
3868 if (sqlite3_column_int (stmt, 1) == 1)
3869 needs_interpolation = 1;
3870 }
3871 }
3872
3873 pt = dyn->First;
3874 while (pt)
3875 {
3876 /* counting how many points are there */
3877 count++;
3878 pt = pt->Next;
3879 }
3880 if (count < 2)
3881 goto end;
3882
3883 if (needs_interpolation)
3884 {
3885 /* attempting to interpolate missing coords */
3886 int i;
3887 int max = count;
3888 char *interpolate = malloc (max + 1);
3889 memset (interpolate, '\0', max + 1);
3890 sqlite3_reset (stmt); /* rewinding the resultset */
3891 count = 0;
3892 while (1)
3893 {
3894 /* scrolling the result set rows */
3895 ret = sqlite3_step (stmt);
3896 if (ret == SQLITE_DONE)
3897 break; /* end of result set */
3898 if (ret == SQLITE_ROW)
3899 {
3900 if (sqlite3_column_int (stmt, 1) == 0)
3901 *(interpolate + count) = 'N';
3902 else
3903 *(interpolate + count) = 'Y';
3904 count++;
3905 }
3906 }
3907 for (i = 0; i < max; i++)
3908 {
3909 if (*(interpolate + i) == 'Y')
3910 do_interpolate_coords (i, dyn, interpolate);
3911 }
3912 free (interpolate);
3913 }
3914
3915 sqlite3_finalize (stmt);
3916 stmt = NULL;
3917
3918 /* creating the final Linestring */
3919 if (dims == GAIA_XY_Z_M)
3920 geom = gaiaAllocGeomCollXYZM ();
3921 else if (dims == GAIA_XY_Z)
3922 geom = gaiaAllocGeomCollXYZ ();
3923 else if (dims == GAIA_XY_M)
3924 geom = gaiaAllocGeomCollXYM ();
3925 else
3926 geom = gaiaAllocGeomColl ();
3927 geom->Srid = srid;
3928
3929 ln = gaiaAddLinestringToGeomColl (geom, count);
3930 count = 0;
3931 pt = dyn->First;
3932 while (pt != NULL)
3933 {
3934 if (dims == GAIA_XY_Z_M)
3935 {
3936 gaiaSetPointXYZM (ln->Coords, count, pt->X, pt->Y,
3937 pt->Z, pt->M);
3938 }
3939 else if (dims == GAIA_XY_Z)
3940 {
3941 gaiaSetPointXYZ (ln->Coords, count, pt->X, pt->Y, pt->Z);
3942 }
3943 else if (dims == GAIA_XY_M)
3944 {
3945 gaiaSetPointXYM (ln->Coords, count, pt->X, pt->Y, pt->M);
3946 }
3947 else
3948 {
3949 gaiaSetPoint (ln->Coords, count, pt->X, pt->Y);
3950 }
3951 count++;
3952 pt = pt->Next;
3953 }
3954
3955 end:
3956 gaiaFreeDynamicLine (dyn);
3957 if (stmt != NULL)
3958 sqlite3_finalize (stmt);
3959 return geom;
3960 }
3961
3962 static gaiaGeomCollPtr
do_reassemble_multi_point(sqlite3 * mem_db,int dims,int srid,int interpolated)3963 do_reassemble_multi_point (sqlite3 * mem_db, int dims, int srid,
3964 int interpolated)
3965 {
3966 /* reassembling the MultiPoint to be returned */
3967 gaiaGeomCollPtr geom = NULL;
3968 gaiaPointPtr pt;
3969 gaiaDynamicLinePtr dyn = gaiaAllocDynamicLine ();
3970 const char *sql;
3971 int ret;
3972 sqlite3_stmt *stmt = NULL;
3973 int count = 0;
3974 int needs_interpolation = 0;
3975 char *interpolate = NULL;
3976
3977 /* creating an SQL statement for querying Points-1 */
3978 sql = "SELECT geom, needs_interpolation FROM points1 ORDER BY id";
3979 ret = sqlite3_prepare_v2 (mem_db, sql, strlen (sql), &stmt, NULL);
3980 if (ret != SQLITE_OK)
3981 {
3982 spatialite_e ("SELECT Points1: error %d \"%s\"\n",
3983 sqlite3_errcode (mem_db), sqlite3_errmsg (mem_db));
3984 goto end;
3985 }
3986 while (1)
3987 {
3988 /* scrolling the result set rows */
3989 ret = sqlite3_step (stmt);
3990 if (ret == SQLITE_DONE)
3991 break; /* end of result set */
3992 if (ret == SQLITE_ROW)
3993 {
3994 gaiaGeomCollPtr g = NULL;
3995 if (sqlite3_column_type (stmt, 0) == SQLITE_BLOB)
3996 {
3997 unsigned char *p_blob =
3998 (unsigned char *) sqlite3_column_blob (stmt, 0);
3999 int n_bytes = sqlite3_column_bytes (stmt, 0);
4000 g = gaiaFromSpatiaLiteBlobWkb (p_blob, n_bytes);
4001 }
4002 if (g != NULL)
4003 {
4004 /* found a valid Point */
4005 pt = g->FirstPoint;
4006 if (dims == GAIA_XY_Z_M)
4007 gaiaAppendPointZMToDynamicLine (dyn, pt->X, pt->Y,
4008 pt->Z, pt->M);
4009 else if (dims == GAIA_XY_Z)
4010 gaiaAppendPointZToDynamicLine (dyn, pt->X, pt->Y,
4011 pt->Z);
4012 else if (dims == GAIA_XY_M)
4013 gaiaAppendPointMToDynamicLine (dyn, pt->X, pt->Y,
4014 pt->M);
4015 else
4016 gaiaAppendPointToDynamicLine (dyn, pt->X, pt->Y);
4017 gaiaFreeGeomColl (g);
4018 }
4019 if (sqlite3_column_int (stmt, 1) == 1)
4020 needs_interpolation = 1;
4021 }
4022 }
4023
4024 pt = dyn->First;
4025 while (pt)
4026 {
4027 /* counting how many points are there */
4028 count++;
4029 pt = pt->Next;
4030 }
4031 if (count < 2)
4032 goto end;
4033
4034 if (needs_interpolation)
4035 {
4036 /* attempting to interpolate missing coords */
4037 int i;
4038 int max = count;
4039 interpolate = malloc (max + 1);
4040 memset (interpolate, '\0', max + 1);
4041 sqlite3_reset (stmt); /* rewinding the resultset */
4042 count = 0;
4043 while (1)
4044 {
4045 /* scrolling the result set rows */
4046 ret = sqlite3_step (stmt);
4047 if (ret == SQLITE_DONE)
4048 break; /* end of result set */
4049 if (ret == SQLITE_ROW)
4050 {
4051 if (sqlite3_column_int (stmt, 1) == 0)
4052 *(interpolate + count) = 'N';
4053 else
4054 *(interpolate + count) = 'Y';
4055 count++;
4056 }
4057 }
4058 for (i = 0; i < max; i++)
4059 {
4060 if (*(interpolate + i) == 'Y')
4061 do_interpolate_coords (i, dyn, interpolate);
4062 }
4063 }
4064
4065 sqlite3_finalize (stmt);
4066 stmt = NULL;
4067
4068 /* creating the final MultiPoint */
4069 if (dims == GAIA_XY_Z_M)
4070 geom = gaiaAllocGeomCollXYZM ();
4071 else if (dims == GAIA_XY_Z)
4072 geom = gaiaAllocGeomCollXYZ ();
4073 else if (dims == GAIA_XY_M)
4074 geom = gaiaAllocGeomCollXYM ();
4075 else
4076 geom = gaiaAllocGeomColl ();
4077 geom->Srid = srid;
4078 geom->DeclaredType = GAIA_MULTIPOINT;
4079
4080 pt = dyn->First;
4081 count = 0;
4082 while (pt != NULL)
4083 {
4084 int ok = 0;
4085 if (*(interpolate + count) == 'Y')
4086 ok = 1;
4087 if (!interpolated && *(interpolate + count) == 'I')
4088 ok = 1;
4089 if (ok)
4090 {
4091 if (dims == GAIA_XY_Z_M)
4092 gaiaAddPointToGeomCollXYZM (geom, pt->X, pt->Y, pt->Z,
4093 pt->M);
4094 else if (dims == GAIA_XY_Z)
4095 gaiaAddPointToGeomCollXYZ (geom, pt->X, pt->Y, pt->Z);
4096 else if (dims == GAIA_XY_M)
4097 gaiaAddPointToGeomCollXYM (geom, pt->X, pt->Y, pt->M);
4098 else
4099 gaiaAddPointToGeomColl (geom, pt->X, pt->Y);
4100 }
4101 count++;
4102 pt = pt->Next;
4103 }
4104
4105 end:
4106 if (interpolate != NULL)
4107 free (interpolate);
4108 gaiaFreeDynamicLine (dyn);
4109 if (stmt != NULL)
4110 sqlite3_finalize (stmt);
4111 return geom;
4112 }
4113
4114 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDrapeLine(sqlite3 * db_handle,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double tolerance)4115 gaiaDrapeLine (sqlite3 * db_handle, gaiaGeomCollPtr geom1,
4116 gaiaGeomCollPtr geom2, double tolerance)
4117 {
4118 /* will return a 3D Linestring by draping line-1 (2D) over line-2 (3D) */
4119 int ret;
4120 sqlite3 *mem_db = NULL;
4121 void *cache;
4122 char *sql;
4123 char *err_msg = NULL;
4124 gaiaGeomCollPtr geom3 = NULL;
4125
4126 /* arguments validation */
4127 if (db_handle == NULL)
4128 return NULL;
4129 if (geom1 == NULL || geom2 == NULL)
4130 return NULL;
4131 if (tolerance < 0.0)
4132 return NULL;
4133 if (geom1->Srid != geom2->Srid)
4134 return NULL;
4135 if (geom1->DimensionModel != GAIA_XY)
4136 return NULL;
4137 if (geom2->DimensionModel != GAIA_XY_Z)
4138 return NULL;
4139 if (!gaia_do_check_linestring (geom1))
4140 return NULL;
4141 if (!gaia_do_check_linestring (geom2))
4142 return NULL;
4143
4144 /* creating a Temporary MemoryDB */
4145 ret =
4146 sqlite3_open_v2 (":memory:", &mem_db,
4147 SQLITE_OPEN_READWRITE | SQLITE_OPEN_CREATE, NULL);
4148 if (ret != SQLITE_OK)
4149 {
4150 spatialite_e ("gaiaDrapeLine: sqlite3_open_v2 error: %s\n",
4151 sqlite3_errmsg (mem_db));
4152 sqlite3_close (mem_db);
4153 return NULL;
4154 }
4155 cache = spatialite_alloc_connection ();
4156 spatialite_internal_init (mem_db, cache);
4157
4158 /* initializing a minimal SpatiaLite DB */
4159 sql = "SELECT InitSpatialMetadata(1, 'NONE')";
4160 ret = sqlite3_exec (mem_db, sql, NULL, NULL, &err_msg);
4161 if (ret != SQLITE_OK)
4162 {
4163 spatialite_e ("gaiaDrapeLine: InitSpatialMetadata() error: %s\n",
4164 err_msg);
4165 sqlite3_free (err_msg);
4166 goto end;
4167 }
4168
4169 /* creating the helper tables on the Temporary MemoryDB */
4170 if (!do_create_points (mem_db, "points1"))
4171 goto end;
4172 if (!do_create_points (mem_db, "points2"))
4173 goto end;
4174
4175 /* populating the Points-2 helper table */
4176 if (!do_populate_points2 (mem_db, geom2))
4177 goto end;
4178
4179 /* draping the first line over Points-2 */
4180 if (!do_drape_line (mem_db, geom1, tolerance))
4181 goto end;
4182
4183 /* building the final linestring to be returned */
4184 geom3 = do_reassemble_line (mem_db, geom2->DimensionModel, geom2->Srid);
4185
4186 end:
4187 /* releasing the Temporary MemoryDB */
4188 ret = sqlite3_close (mem_db);
4189 if (ret != SQLITE_OK)
4190 spatialite_e ("gaiaDrapeLine: sqlite3_close() error: %s\n",
4191 sqlite3_errmsg (mem_db));
4192 spatialite_internal_cleanup (cache);
4193
4194 return geom3;
4195 }
4196
4197 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaDrapeLineExceptions(sqlite3 * db_handle,gaiaGeomCollPtr geom1,gaiaGeomCollPtr geom2,double tolerance,int interpolated)4198 gaiaDrapeLineExceptions (sqlite3 * db_handle, gaiaGeomCollPtr geom1,
4199 gaiaGeomCollPtr geom2, double tolerance,
4200 int interpolated)
4201 {
4202 /*
4203 * will return a 3D MultiPoint containing all Vertices from geom1
4204 * not being correctly draped over geom2
4205 */
4206 int ret;
4207 sqlite3 *mem_db = NULL;
4208 void *cache;
4209 char *sql;
4210 char *err_msg = NULL;
4211 gaiaGeomCollPtr geom3 = NULL;
4212
4213 /* arguments validation */
4214 if (db_handle == NULL)
4215 return NULL;
4216 if (geom1 == NULL || geom2 == NULL)
4217 return NULL;
4218 if (tolerance < 0.0)
4219 return NULL;
4220 if (geom1->Srid != geom2->Srid)
4221 return NULL;
4222 if (geom1->DimensionModel != GAIA_XY)
4223 return NULL;
4224 if (geom2->DimensionModel != GAIA_XY_Z)
4225 return NULL;
4226 if (!gaia_do_check_linestring (geom1))
4227 return NULL;
4228 if (!gaia_do_check_linestring (geom2))
4229 return NULL;
4230
4231 /* creating a Temporary MemoryDB */
4232 ret =
4233 sqlite3_open_v2 (":memory:", &mem_db,
4234 SQLITE_OPEN_READWRITE | SQLITE_OPEN_CREATE, NULL);
4235 if (ret != SQLITE_OK)
4236 {
4237 spatialite_e ("gaiaDrapeLine: sqlite3_open_v2 error: %s\n",
4238 sqlite3_errmsg (mem_db));
4239 sqlite3_close (mem_db);
4240 return NULL;
4241 }
4242 cache = spatialite_alloc_connection ();
4243 spatialite_internal_init (mem_db, cache);
4244
4245 /* initializing a minimal SpatiaLite DB */
4246 sql = "SELECT InitSpatialMetadata(1, 'NONE')";
4247 ret = sqlite3_exec (mem_db, sql, NULL, NULL, &err_msg);
4248 if (ret != SQLITE_OK)
4249 {
4250 spatialite_e
4251 ("gaiaDrapeLineExceptions: InitSpatialMetadata() error: %s\n",
4252 err_msg);
4253 sqlite3_free (err_msg);
4254 goto end;
4255 }
4256
4257 /* creating the helper tables on the Temporary MemoryDB */
4258 if (!do_create_points (mem_db, "points1"))
4259 goto end;
4260 if (!do_create_points (mem_db, "points2"))
4261 goto end;
4262
4263 /* populating the Points-2 helper table */
4264 if (!do_populate_points2 (mem_db, geom2))
4265 goto end;
4266
4267 /* draping the first line over Points-2 */
4268 if (!do_drape_line (mem_db, geom1, tolerance))
4269 goto end;
4270
4271 /* building the final MultiPoint to be returned */
4272 geom3 =
4273 do_reassemble_multi_point (mem_db, geom2->DimensionModel, geom2->Srid,
4274 interpolated);
4275
4276 end:
4277 /* releasing the Temporary MemoryDB */
4278 ret = sqlite3_close (mem_db);
4279 if (ret != SQLITE_OK)
4280 spatialite_e ("gaiaDrapeLineExceptions: sqlite3_close() error: %s\n",
4281 sqlite3_errmsg (mem_db));
4282 spatialite_internal_cleanup (cache);
4283
4284 return geom3;
4285 }
4286