1 /*
2
3 gg_voronoj.c -- Voronoj Diagram implementation
4
5 version 5.0, 2020 August 1
6
7 Author: Sandro Furieri a.furieri@lqt.it
8
9 ------------------------------------------------------------------------------
10
11 Version: MPL 1.1/GPL 2.0/LGPL 2.1
12
13 The contents of this file are subject to the Mozilla Public License Version
14 1.1 (the "License"); you may not use this file except in compliance with
15 the License. You may obtain a copy of the License at
16 http://www.mozilla.org/MPL/
17
18 Software distributed under the License is distributed on an "AS IS" basis,
19 WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
20 for the specific language governing rights and limitations under the
21 License.
22
23 The Original Code is the SpatiaLite library
24
25 The Initial Developer of the Original Code is Alessandro Furieri
26
27 Portions created by the Initial Developer are Copyright (C) 2008-2021
28 the Initial Developer. All Rights Reserved.
29
30 Contributor(s):
31
32 Alternatively, the contents of this file may be used under the terms of
33 either the GNU General Public License Version 2 or later (the "GPL"), or
34 the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
35 in which case the provisions of the GPL or the LGPL are applicable instead
36 of those above. If you wish to allow use of your version of this file only
37 under the terms of either the GPL or the LGPL, and not to allow others to
38 use your version of this file under the terms of the MPL, indicate your
39 decision by deleting the provisions above and replace them with the notice
40 and other provisions required by the GPL or the LGPL. If you do not delete
41 the provisions above, a recipient may use your version of this file under
42 the terms of any one of the MPL, the GPL or the LGPL.
43
44 */
45
46 #include <sys/types.h>
47 #include <stdlib.h>
48 #include <stdio.h>
49 #include <string.h>
50 #include <float.h>
51 #include <math.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_private.h>
60 #include <spatialite/sqlite.h>
61
62 #include <spatialite/gaiageo.h>
63
64 #ifdef GEOS_ADVANCED /* GEOS advanced features */
65
66 struct voronoj_triangle
67 {
68 /* a struct representing a Delaunay triangle */
69 double x1; /* vertex #1 */
70 double y1;
71 double x2; /* vertex #2 */
72 double y2;
73 double x3; /* vertex #3 */
74 double y3;
75 double cx; /* circumcenter */
76 double cy;
77 double x_1_2; /* vertex on the frame - edge #1-#2 */
78 double y_1_2;
79 double x_2_3; /* vertex on the frame - edge #2-#3 */
80 double y_2_3;
81 double x_3_1; /* vertex on the frame - edge #3-#1 */
82 double y_3_1;
83
84 struct voronoj_triangle *tri_1_2; /* triangle sharing edge #1-#2 */
85 struct voronoj_triangle *tri_2_3; /* triangle sharing edge #2-#3 */
86 struct voronoj_triangle *tri_3_1; /* triangle sharing edge #3-#1 */
87 char trace_1_2; /* flags: to be traced */
88 char trace_2_3;
89 char trace_3_1;
90 };
91
92 struct voronoj_point
93 {
94 /* an auxiliary struct - point on Voronoj's frame */
95 double coord;
96 struct voronoj_point *next;
97 };
98
99 struct voronoj_aux
100 {
101 /* an auxiliary struct supporting Voronoj */
102 struct voronoj_triangle *array; /* array of Triangles */
103 int count; /* number of Triangles */
104 double minx; /* the frame extent */
105 double miny;
106 double maxx;
107 double maxy;
108 struct voronoj_point *first_up; /* linked list: horz-up frame's edge */
109 struct voronoj_point *last_up;
110 struct voronoj_point *first_low; /* linked list: horz-down frame's edge */
111 struct voronoj_point *last_low;
112 struct voronoj_point *first_left; /* linked list: vert-left frame's edge */
113 struct voronoj_point *last_left;
114 struct voronoj_point *first_right; /* linked list: vert-right frame's edge */
115 struct voronoj_point *last_right;
116 };
117
118 struct concave_hull_str
119 {
120 /* a struct to implement StandardVariation and Variance for Concave Hull */
121 double mean;
122 double quot;
123 double count;
124 };
125
126 #ifndef GEOS_REENTRANT /* GEOS >= 3.5.0 directly supports Voronoj */
127
128 static double *
voronoj_sorted_up(struct voronoj_aux * voronoj,int * count)129 voronoj_sorted_up (struct voronoj_aux *voronoj, int *count)
130 {
131 /* returning a sorted array of coordinates */
132 double *array = NULL;
133 int cnt = 0;
134 int ok = 1;
135 int i;
136 struct voronoj_point *pt = voronoj->first_up;
137 while (pt)
138 {
139 /* counting how many points are there */
140 cnt++;
141 pt = pt->next;
142 }
143 *count = cnt;
144 if (cnt == 0)
145 return NULL;
146
147 /* allocating and populating the array */
148 array = malloc (sizeof (double) * *count);
149 cnt = 0;
150 pt = voronoj->first_up;
151 while (pt)
152 {
153 *(array + cnt++) = pt->coord;
154 pt = pt->next;
155 }
156
157 /* sorting the array */
158 while (ok)
159 {
160 ok = 0;
161 for (i = 1; i < *count; i++)
162 {
163 if (*(array + i - 1) > *(array + i))
164 {
165 /* swapping two values */
166 double save = *(array + i - 1);
167 *(array + i - 1) = *(array + i);
168 *(array + i) = save;
169 ok = 1;
170 }
171 }
172 }
173 return array;
174 }
175
176 static double *
voronoj_sorted_low(struct voronoj_aux * voronoj,int * count)177 voronoj_sorted_low (struct voronoj_aux *voronoj, int *count)
178 {
179 /* returning a sorted array of coordinates */
180 double *array = NULL;
181 int cnt = 0;
182 int ok = 1;
183 int i;
184 struct voronoj_point *pt = voronoj->first_low;
185 while (pt)
186 {
187 /* counting how many points are there */
188 cnt++;
189 pt = pt->next;
190 }
191 *count = cnt;
192 if (cnt == 0)
193 return NULL;
194
195 /* allocating and populating the array */
196 array = malloc (sizeof (double) * *count);
197 cnt = 0;
198 pt = voronoj->first_low;
199 while (pt)
200 {
201 *(array + cnt++) = pt->coord;
202 pt = pt->next;
203 }
204
205 /* sorting the array */
206 while (ok)
207 {
208 ok = 0;
209 for (i = 1; i < *count; i++)
210 {
211 if (*(array + i - 1) > *(array + i))
212 {
213 /* swapping two values */
214 double save = *(array + i - 1);
215 *(array + i - 1) = *(array + i);
216 *(array + i) = save;
217 ok = 1;
218 }
219 }
220 }
221 return array;
222 }
223
224 static double *
voronoj_sorted_left(struct voronoj_aux * voronoj,int * count)225 voronoj_sorted_left (struct voronoj_aux *voronoj, int *count)
226 {
227 /* returning a sorted array of coordinates */
228 double *array = NULL;
229 int cnt = 0;
230 int ok = 1;
231 int i;
232 struct voronoj_point *pt = voronoj->first_left;
233 while (pt)
234 {
235 /* counting how many points are there */
236 cnt++;
237 pt = pt->next;
238 }
239 *count = cnt;
240 if (cnt == 0)
241 return NULL;
242
243 /* allocating and populating the array */
244 array = malloc (sizeof (double) * *count);
245 cnt = 0;
246 pt = voronoj->first_left;
247 while (pt)
248 {
249 *(array + cnt++) = pt->coord;
250 pt = pt->next;
251 }
252
253 /* sorting the array */
254 while (ok)
255 {
256 ok = 0;
257 for (i = 1; i < *count; i++)
258 {
259 if (*(array + i - 1) > *(array + i))
260 {
261 /* swapping two values */
262 double save = *(array + i - 1);
263 *(array + i - 1) = *(array + i);
264 *(array + i) = save;
265 ok = 1;
266 }
267 }
268 }
269 return array;
270 }
271
272 static double *
voronoj_sorted_right(struct voronoj_aux * voronoj,int * count)273 voronoj_sorted_right (struct voronoj_aux *voronoj, int *count)
274 {
275 /* returning a sorted array of coordinates */
276 double *array = NULL;
277 int cnt = 0;
278 int ok = 1;
279 int i;
280 struct voronoj_point *pt = voronoj->first_right;
281 while (pt)
282 {
283 /* counting how many points are there */
284 cnt++;
285 pt = pt->next;
286 }
287 *count = cnt;
288 if (cnt == 0)
289 return NULL;
290
291 /* allocating and populating the array */
292 array = malloc (sizeof (double) * *count);
293 cnt = 0;
294 pt = voronoj->first_right;
295 while (pt)
296 {
297 *(array + cnt++) = pt->coord;
298 pt = pt->next;
299 }
300
301 /* sorting the array */
302 while (ok)
303 {
304 ok = 0;
305 for (i = 1; i < *count; i++)
306 {
307 if (*(array + i - 1) > *(array + i))
308 {
309 /* swapping two values */
310 double save = *(array + i - 1);
311 *(array + i - 1) = *(array + i);
312 *(array + i) = save;
313 ok = 1;
314 }
315 }
316 }
317 return array;
318 }
319
320 static void
voronoj_add_frame_point(struct voronoj_aux * voronoj,double x,double y)321 voronoj_add_frame_point (struct voronoj_aux *voronoj, double x, double y)
322 {
323 /* adding some frame point */
324 struct voronoj_point *pt = NULL;
325
326 /* skipping any corner */
327 if (x == voronoj->minx && y == voronoj->miny)
328 return;
329 if (x == voronoj->minx && y == voronoj->maxy)
330 return;
331 if (x == voronoj->maxx && y == voronoj->miny)
332 return;
333 if (x == voronoj->maxx && y == voronoj->maxy)
334 return;
335
336 if (x == voronoj->minx)
337 {
338 pt = malloc (sizeof (struct voronoj_point));
339 pt->coord = y;
340 pt->next = NULL;
341 if (voronoj->first_left == NULL)
342 voronoj->first_left = pt;
343 if (voronoj->last_left != NULL)
344 voronoj->last_left->next = pt;
345 voronoj->last_left = pt;
346 }
347 if (x == voronoj->maxx)
348 {
349 pt = malloc (sizeof (struct voronoj_point));
350 pt->coord = y;
351 pt->next = NULL;
352 if (voronoj->first_right == NULL)
353 voronoj->first_right = pt;
354 if (voronoj->last_right != NULL)
355 voronoj->last_right->next = pt;
356 voronoj->last_right = pt;
357 }
358 if (y == voronoj->miny)
359 {
360 pt = malloc (sizeof (struct voronoj_point));
361 pt->coord = x;
362 pt->next = NULL;
363 if (voronoj->first_low == NULL)
364 voronoj->first_low = pt;
365 if (voronoj->last_low != NULL)
366 voronoj->last_low->next = pt;
367 voronoj->last_low = pt;
368 }
369 if (y == voronoj->maxy)
370 {
371 pt = malloc (sizeof (struct voronoj_point));
372 pt->coord = x;
373 pt->next = NULL;
374 if (voronoj->first_up == NULL)
375 voronoj->first_up = pt;
376 if (voronoj->last_up != NULL)
377 voronoj->last_up->next = pt;
378 voronoj->last_up = pt;
379 }
380 }
381
382 static int
voronoj_same_edge(double ax1,double ay1,double ax2,double ay2,double bx1,double by1,double bx2,double by2)383 voronoj_same_edge (double ax1, double ay1, double ax2, double ay2, double bx1,
384 double by1, double bx2, double by2)
385 {
386 /* testing if two segments are the same */
387 if (ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2)
388 return 1;
389 if (ax1 == bx2 && ay1 == by2 && ax2 == bx1 && ay2 == by1)
390 return 1;
391 return 0;
392 }
393
394 static int
voronoj_internal(const void * p_cache,struct voronoj_triangle * triangle)395 voronoj_internal (const void *p_cache, struct voronoj_triangle *triangle)
396 {
397 /* checking if the circumcenter falls inside the triangle */
398 int ret;
399 gaiaGeomCollPtr pt = gaiaAllocGeomColl ();
400 gaiaGeomCollPtr tri = gaiaAllocGeomColl ();
401 gaiaPolygonPtr pg = gaiaAddPolygonToGeomColl (tri, 4, 0);
402 gaiaRingPtr rng = pg->Exterior;
403 gaiaSetPoint (rng->Coords, 0, triangle->x1, triangle->y1);
404 gaiaSetPoint (rng->Coords, 1, triangle->x2, triangle->y2);
405 gaiaSetPoint (rng->Coords, 2, triangle->x3, triangle->y3);
406 gaiaSetPoint (rng->Coords, 3, triangle->x1, triangle->y1);
407 gaiaAddPointToGeomColl (pt, triangle->cx, triangle->cy);
408 gaiaMbrGeometry (pt);
409 gaiaMbrGeometry (tri);
410 if (p_cache != NULL)
411 ret = gaiaGeomCollIntersects_r (p_cache, tri, pt);
412 else
413 ret = gaiaGeomCollIntersects (tri, pt);
414 gaiaFreeGeomColl (pt);
415 gaiaFreeGeomColl (tri);
416 return ret;
417 }
418
419 static double
voronoj_test_point(const void * p_cache,double x1,double y1,double x2,double y2,double x,double y)420 voronoj_test_point (const void *p_cache, double x1, double y1, double x2,
421 double y2, double x, double y)
422 {
423 /* point-segment distance */
424 double dist;
425 gaiaGeomCollPtr pt = gaiaAllocGeomColl ();
426 gaiaGeomCollPtr segm = gaiaAllocGeomColl ();
427 gaiaLinestringPtr ln = gaiaAddLinestringToGeomColl (segm, 2);
428 gaiaSetPoint (ln->Coords, 0, x1, y1);
429 gaiaSetPoint (ln->Coords, 1, x2, y2);
430 gaiaAddPointToGeomColl (pt, x, y);
431 if (p_cache != NULL)
432 gaiaGeomCollDistance_r (p_cache, segm, pt, &dist);
433 else
434 gaiaGeomCollDistance (segm, pt, &dist);
435 gaiaFreeGeomColl (pt);
436 gaiaFreeGeomColl (segm);
437 return dist;
438 }
439
440 static int
voronoj_check_nearest_edge(const void * p_cache,struct voronoj_triangle * tri,int which)441 voronoj_check_nearest_edge (const void *p_cache, struct voronoj_triangle *tri,
442 int which)
443 {
444 /* testing if direction outside */
445 double d_1_2;
446 double d_2_3;
447 double d_3_1;
448 gaiaGeomCollPtr pt = gaiaAllocGeomColl ();
449 gaiaGeomCollPtr segm = gaiaAllocGeomColl ();
450 gaiaLinestringPtr ln = gaiaAddLinestringToGeomColl (segm, 2);
451 gaiaSetPoint (ln->Coords, 0, tri->x1, tri->y1);
452 gaiaSetPoint (ln->Coords, 1, tri->x2, tri->y2);
453 gaiaAddPointToGeomColl (pt, tri->cx, tri->cy);
454 if (p_cache != NULL)
455 gaiaGeomCollDistance_r (p_cache, segm, pt, &d_1_2);
456 else
457 gaiaGeomCollDistance (segm, pt, &d_1_2);
458 gaiaFreeGeomColl (segm);
459 segm = gaiaAllocGeomColl ();
460 ln = gaiaAddLinestringToGeomColl (segm, 2);
461 gaiaSetPoint (ln->Coords, 0, tri->x2, tri->y2);
462 gaiaSetPoint (ln->Coords, 1, tri->x3, tri->y3);
463 if (p_cache != NULL)
464 gaiaGeomCollDistance_r (p_cache, segm, pt, &d_2_3);
465 else
466 gaiaGeomCollDistance (segm, pt, &d_2_3);
467 gaiaFreeGeomColl (segm);
468 segm = gaiaAllocGeomColl ();
469 ln = gaiaAddLinestringToGeomColl (segm, 2);
470 gaiaSetPoint (ln->Coords, 0, tri->x3, tri->y3);
471 gaiaSetPoint (ln->Coords, 1, tri->x1, tri->y1);
472 if (p_cache != NULL)
473 gaiaGeomCollDistance_r (p_cache, segm, pt, &d_3_1);
474 else
475 gaiaGeomCollDistance (segm, pt, &d_3_1);
476 gaiaFreeGeomColl (segm);
477 gaiaFreeGeomColl (pt);
478
479 if (which == 12 && d_1_2 < d_2_3 && d_1_2 < d_3_1)
480 return 0;
481 if (which == 23 && d_2_3 < d_1_2 && d_2_3 < d_3_1)
482 return 0;
483 if (which == 31 && d_3_1 < d_1_2 && d_3_1 < d_2_3)
484 return 0;
485 return 1;
486 }
487
488 static void
voronoj_minmax(double x,double y,double * minx,double * miny,double * maxx,double * maxy)489 voronoj_minmax (double x, double y, double *minx, double *miny, double *maxx,
490 double *maxy)
491 {
492 /* updating the frame extent */
493 if (x < *minx)
494 *minx = x;
495 if (y < *miny)
496 *miny = y;
497 if (x > *maxx)
498 *maxx = x;
499 if (y > *maxy)
500 *maxy = y;
501 }
502
503 static void
voronoj_frame_point(const void * p_cache,double intercept,double slope,struct voronoj_aux * voronoj,double cx,double cy,double mx,double my,int direct,double * x,double * y)504 voronoj_frame_point (const void *p_cache, double intercept, double slope,
505 struct voronoj_aux *voronoj, double cx, double cy,
506 double mx, double my, int direct, double *x, double *y)
507 {
508 /* determining a vertex on the frame */
509 double x_up;
510 double x_low;
511 double y_left;
512 double y_right;
513 double pre_x1 = DBL_MAX;
514 double pre_y1 = DBL_MAX;
515 double pre_x2 = DBL_MAX;
516 double pre_y2 = DBL_MAX;
517 double d1;
518 double d2;
519
520 if (slope == DBL_MAX)
521 {
522 x_up = cx;
523 x_low = cx;
524 y_left = DBL_MAX;
525 y_right = DBL_MAX;
526 }
527 else
528 {
529 x_up = (voronoj->maxy - intercept) / slope;
530 x_low = (voronoj->miny - intercept) / slope;
531 y_left = (slope * voronoj->minx) + intercept;
532 y_right = (slope * voronoj->maxx) + intercept;
533 }
534
535 if (x_up >= voronoj->minx && x_up <= voronoj->maxx)
536 {
537 if (pre_x1 == DBL_MAX && pre_y1 == DBL_MAX)
538 {
539 pre_x1 = x_up;
540 pre_y1 = voronoj->maxy;
541 }
542 else if (pre_x2 == DBL_MAX && pre_y2 == DBL_MAX)
543 {
544 pre_x2 = x_up;
545 pre_y2 = voronoj->maxy;
546 }
547 }
548 if (x_low >= voronoj->minx && x_low <= voronoj->maxx)
549 {
550 if (pre_x1 == DBL_MAX && pre_y1 == DBL_MAX)
551 {
552 pre_x1 = x_low;
553 pre_y1 = voronoj->miny;
554 }
555 else if (pre_x2 == DBL_MAX && pre_y2 == DBL_MAX)
556 {
557 pre_x2 = x_low;
558 pre_y2 = voronoj->miny;
559 }
560 }
561 if (y_left >= voronoj->miny && y_left <= voronoj->maxy)
562 {
563 if (pre_x1 == DBL_MAX && pre_y1 == DBL_MAX)
564 {
565 pre_x1 = voronoj->minx;
566 pre_y1 = y_left;
567 }
568 else if (pre_x2 == DBL_MAX && pre_y2 == DBL_MAX)
569 {
570 pre_x2 = voronoj->minx;
571 pre_y2 = y_left;
572 }
573 }
574 if (y_right >= voronoj->miny && y_right <= voronoj->maxy)
575 {
576 if (pre_x1 == DBL_MAX && pre_y1 == DBL_MAX)
577 {
578 pre_x1 = voronoj->maxx;
579 pre_y1 = y_right;
580 }
581 else if (pre_x2 == DBL_MAX && pre_y2 == DBL_MAX)
582 {
583 pre_x2 = voronoj->maxx;
584 pre_y2 = y_right;
585 }
586 }
587
588 /* choosing wich point has to be returned */
589 if (direct)
590 {
591 /* cutting the edge in two */
592 d1 = voronoj_test_point (p_cache, cx, cy, pre_x1, pre_y1, mx, my);
593 d2 = voronoj_test_point (p_cache, cx, cy, pre_x2, pre_y2, mx, my);
594 if (d1 < d2)
595 {
596 *x = pre_x1;
597 *y = pre_y1;
598 }
599 else
600 {
601 *x = pre_x2;
602 *y = pre_y2;
603 }
604 }
605 else
606 {
607 /* going outside */
608 d1 = voronoj_test_point (p_cache, cx, cy, pre_x1, pre_y1, mx, my);
609 d2 = voronoj_test_point (p_cache, cx, cy, pre_x2, pre_y2, mx, my);
610 if (d1 > d2)
611 {
612 *x = pre_x1;
613 *y = pre_y1;
614 }
615 else
616 {
617 *x = pre_x2;
618 *y = pre_y2;
619 }
620 }
621 }
622
623 SPATIALITE_PRIVATE void *
voronoj_build(int count,void * p_first,double extra_frame_size)624 voronoj_build (int count, void *p_first, double extra_frame_size)
625 {
626 return voronoj_build_r (NULL, count, p_first, extra_frame_size);
627 }
628
629 SPATIALITE_PRIVATE void *
voronoj_build_r(const void * p_cache,int count,void * p_first,double extra_frame_size)630 voronoj_build_r (const void *p_cache, int count, void *p_first,
631 double extra_frame_size)
632 {
633 /* building the Voronoj auxiliary struct */
634 gaiaPolygonPtr first = (gaiaPolygonPtr) p_first;
635 struct voronoj_aux *voronoj = NULL;
636 struct voronoj_triangle *triangle;
637 struct voronoj_triangle *tri2;
638 gaiaPolygonPtr pg;
639 gaiaRingPtr rng;
640 int ind = 0;
641 int i2;
642 int direct;
643 double x;
644 double y;
645 double z;
646 double m;
647 double xba;
648 double yba;
649 double xca;
650 double yca;
651 double bl;
652 double cl;
653 double d;
654 double mx;
655 double my;
656 double slope;
657 double intercept = 0.0;
658 double minx = DBL_MAX;
659 double miny = DBL_MAX;
660 double maxx = -DBL_MAX;
661 double maxy = -DBL_MAX;
662 double ext_x;
663 double ext_y;
664 double delta;
665 double delta2;
666
667 /* allocating the Voronoj struct */
668 voronoj = malloc (sizeof (struct voronoj_aux));
669 voronoj->count = count;
670 voronoj->first_up = NULL;
671 voronoj->last_up = NULL;
672 voronoj->first_low = NULL;
673 voronoj->last_low = NULL;
674 voronoj->first_left = NULL;
675 voronoj->last_left = NULL;
676 voronoj->first_right = NULL;
677 voronoj->last_right = NULL;
678 voronoj->array = malloc (sizeof (struct voronoj_triangle) * count);
679
680 /* initializing the Voronoj struct */
681 pg = first;
682 while (pg)
683 {
684 rng = pg->Exterior;
685 triangle = voronoj->array + ind;
686 if (pg->DimensionModel == GAIA_XY_Z)
687 {
688 gaiaGetPointXYZ (rng->Coords, 0, &x, &y, &z);
689 triangle->x1 = x;
690 triangle->y1 = y;
691 gaiaGetPointXYZ (rng->Coords, 1, &x, &y, &z);
692 triangle->x2 = x;
693 triangle->y2 = y;
694 gaiaGetPointXYZ (rng->Coords, 2, &x, &y, &z);
695 triangle->x3 = x;
696 triangle->y3 = y;
697 }
698 else if (pg->DimensionModel == GAIA_XY_M)
699 {
700 gaiaGetPointXYM (rng->Coords, 0, &x, &y, &m);
701 triangle->x1 = x;
702 triangle->y1 = y;
703 gaiaGetPointXYM (rng->Coords, 1, &x, &y, &m);
704 triangle->x2 = x;
705 triangle->y2 = y;
706 gaiaGetPointXYM (rng->Coords, 2, &x, &y, &m);
707 triangle->x3 = x;
708 triangle->y3 = y;
709 }
710 else if (pg->DimensionModel == GAIA_XY_Z_M)
711 {
712 gaiaGetPointXYZM (rng->Coords, 0, &x, &y, &z, &m);
713 triangle->x1 = x;
714 triangle->y1 = y;
715 gaiaGetPointXYZM (rng->Coords, 1, &x, &y, &z, &m);
716 triangle->x2 = x;
717 triangle->y2 = y;
718 gaiaGetPointXYZM (rng->Coords, 2, &x, &y, &z, &m);
719 triangle->x3 = x;
720 triangle->y3 = y;
721 }
722 else
723 {
724 gaiaGetPoint (rng->Coords, 0, &x, &y);
725 triangle->x1 = x;
726 triangle->y1 = y;
727 gaiaGetPoint (rng->Coords, 1, &x, &y);
728 triangle->x2 = x;
729 triangle->y2 = y;
730 gaiaGetPoint (rng->Coords, 2, &x, &y);
731 triangle->x3 = x;
732 triangle->y3 = y;
733 }
734
735 /* computing the triangle circumcenter */
736 xba = triangle->x2 - triangle->x1;
737 yba = triangle->y2 - triangle->y1;
738 xca = triangle->x3 - triangle->x1;
739 yca = triangle->y3 - triangle->y1;
740 bl = xba * xba + yba * yba;
741 cl = xca * xca + yca * yca;
742 d = 0.5 / (xba * yca - yba * xca);
743 triangle->cx = triangle->x1 + ((yca * bl - yba * cl) * d);
744 triangle->cy = triangle->y1 + ((xba * cl - xca * bl) * d);
745 /* adjusting the frame extent */
746 voronoj_minmax (triangle->x1, triangle->y1, &minx, &miny, &maxx,
747 &maxy);
748 voronoj_minmax (triangle->x2, triangle->y2, &minx, &miny, &maxx,
749 &maxy);
750 voronoj_minmax (triangle->x3, triangle->y3, &minx, &miny, &maxx,
751 &maxy);
752 voronoj_minmax (triangle->cx, triangle->cy, &minx, &miny, &maxx,
753 &maxy);
754
755 triangle->tri_1_2 = NULL;
756 triangle->tri_2_3 = NULL;
757 triangle->tri_3_1 = NULL;
758 triangle->trace_1_2 = 0;
759 triangle->trace_2_3 = 0;
760 triangle->trace_3_1 = 0;
761 ind++;
762 pg = pg->Next;
763 }
764
765 /* setting the frame extent */
766 if (extra_frame_size < 0.0)
767 extra_frame_size = 5.0;
768 ext_x = maxx - minx;
769 ext_y = maxy - miny;
770 delta = (ext_x * extra_frame_size) / 100.0;
771 delta2 = (ext_y * extra_frame_size) / 100.0;
772 if (delta2 > delta)
773 delta = delta2;
774 voronoj->minx = minx - delta;
775 voronoj->miny = miny - delta;
776 voronoj->maxx = maxx + delta;
777 voronoj->maxy = maxy + delta;
778
779 /* identifying triangles sharing the same edge */
780 for (ind = 0; ind < voronoj->count; ind++)
781 {
782 triangle = voronoj->array + ind;
783 for (i2 = ind + 1; i2 < voronoj->count; i2++)
784 {
785 tri2 = voronoj->array + i2;
786 if (triangle->tri_1_2 == NULL && tri2->tri_1_2 == NULL)
787 {
788 if (voronoj_same_edge
789 (triangle->x1, triangle->y1, triangle->x2,
790 triangle->y2, tri2->x1, tri2->y1, tri2->x2,
791 tri2->y2))
792 {
793 triangle->tri_1_2 = tri2;
794 triangle->trace_1_2 = 1;
795 tri2->tri_1_2 = triangle;
796 }
797 }
798 if (triangle->tri_1_2 == NULL && tri2->tri_2_3 == NULL)
799 {
800 if (voronoj_same_edge
801 (triangle->x1, triangle->y1, triangle->x2,
802 triangle->y2, tri2->x2, tri2->y2, tri2->x3,
803 tri2->y3))
804 {
805 triangle->tri_1_2 = tri2;
806 triangle->trace_1_2 = 1;
807 tri2->tri_2_3 = triangle;
808 }
809 }
810 if (triangle->tri_1_2 == NULL && tri2->tri_3_1 == NULL)
811 {
812 if (voronoj_same_edge
813 (triangle->x1, triangle->y1, triangle->x2,
814 triangle->y2, tri2->x3, tri2->y3, tri2->x1,
815 tri2->y1))
816 {
817 triangle->tri_1_2 = tri2;
818 triangle->trace_1_2 = 1;
819 tri2->tri_3_1 = triangle;
820 }
821 }
822 if (triangle->tri_2_3 == NULL && tri2->tri_1_2 == NULL)
823 {
824 if (voronoj_same_edge
825 (triangle->x2, triangle->y2, triangle->x3,
826 triangle->y3, tri2->x1, tri2->y1, tri2->x2,
827 tri2->y2))
828 {
829 triangle->tri_2_3 = tri2;
830 triangle->trace_2_3 = 1;
831 tri2->tri_1_2 = triangle;
832 }
833 }
834 if (triangle->tri_2_3 == NULL && tri2->tri_2_3 == NULL)
835 {
836 if (voronoj_same_edge
837 (triangle->x2, triangle->y2, triangle->x3,
838 triangle->y3, tri2->x2, tri2->y2, tri2->x3,
839 tri2->y3))
840 {
841 triangle->tri_2_3 = tri2;
842 triangle->trace_2_3 = 1;
843 tri2->tri_2_3 = triangle;
844 }
845 }
846 if (triangle->tri_2_3 == NULL && tri2->tri_3_1 == NULL)
847 {
848 if (voronoj_same_edge
849 (triangle->x2, triangle->y2, triangle->x3,
850 triangle->y3, tri2->x3, tri2->y3, tri2->x1,
851 tri2->y1))
852 {
853 triangle->tri_2_3 = tri2;
854 triangle->trace_2_3 = 1;
855 tri2->tri_3_1 = triangle;
856 }
857 }
858 if (triangle->tri_3_1 == NULL && tri2->tri_1_2 == NULL)
859 {
860 if (voronoj_same_edge
861 (triangle->x3, triangle->y3, triangle->x1,
862 triangle->y1, tri2->x1, tri2->y1, tri2->x2,
863 tri2->y2))
864 {
865 triangle->tri_3_1 = tri2;
866 triangle->trace_3_1 = 1;
867 tri2->tri_1_2 = triangle;
868 }
869 }
870 if (triangle->tri_3_1 == NULL && tri2->tri_2_3 == NULL)
871 {
872 if (voronoj_same_edge
873 (triangle->x3, triangle->y3, triangle->x1,
874 triangle->y1, tri2->x2, tri2->y2, tri2->x3,
875 tri2->y3))
876 {
877 triangle->tri_3_1 = tri2;
878 triangle->trace_3_1 = 1;
879 tri2->tri_2_3 = triangle;
880 }
881 }
882 if (triangle->tri_3_1 == NULL && tri2->tri_3_1 == NULL)
883 {
884 if (voronoj_same_edge
885 (triangle->x3, triangle->y3, triangle->x1,
886 triangle->y1, tri2->x3, tri2->y3, tri2->x1,
887 tri2->y1))
888 {
889 triangle->tri_3_1 = tri2;
890 triangle->trace_3_1 = 1;
891 tri2->tri_3_1 = triangle;
892 }
893 }
894 }
895
896 /* identifying vertices on the frame */
897 if (triangle->tri_1_2 == NULL)
898 {
899 mx = (triangle->x1 + triangle->x2) / 2.0;
900 my = (triangle->y1 + triangle->y2) / 2.0;
901 if (mx == triangle->cx)
902 slope = DBL_MAX;
903 else
904 {
905 slope = (my - triangle->cy) / (mx - triangle->cx);
906 intercept = my - (slope * mx);
907 }
908 direct = 1;
909 if (!voronoj_internal (p_cache, triangle))
910 direct = voronoj_check_nearest_edge (p_cache, triangle, 12);
911 voronoj_frame_point (p_cache, intercept, slope, voronoj,
912 triangle->cx, triangle->cy, mx, my, direct,
913 &x, &y);
914 triangle->x_1_2 = x;
915 triangle->y_1_2 = y;
916 }
917 if (triangle->tri_2_3 == NULL)
918 {
919 mx = (triangle->x2 + triangle->x3) / 2.0;
920 my = (triangle->y2 + triangle->y3) / 2.0;
921 if (mx == triangle->cx)
922 slope = DBL_MAX;
923 else
924 {
925 slope = (my - triangle->cy) / (mx - triangle->cx);
926 intercept = my - (slope * mx);
927 }
928 direct = 1;
929 if (!voronoj_internal (p_cache, triangle))
930 direct = voronoj_check_nearest_edge (p_cache, triangle, 23);
931 voronoj_frame_point (p_cache, intercept, slope, voronoj,
932 triangle->cx, triangle->cy, mx, my, direct,
933 &x, &y);
934 triangle->x_2_3 = x;
935 triangle->y_2_3 = y;
936 }
937 if (triangle->tri_3_1 == NULL)
938 {
939 mx = (triangle->x3 + triangle->x1) / 2.0;
940 my = (triangle->y3 + triangle->y1) / 2.0;
941 if (mx == triangle->cx)
942 slope = DBL_MAX;
943 else
944 {
945 slope = (my - triangle->cy) / (mx - triangle->cx);
946 intercept = my - (slope * mx);
947 }
948 direct = 1;
949 if (!voronoj_internal (p_cache, triangle))
950 direct = voronoj_check_nearest_edge (p_cache, triangle, 31);
951 voronoj_frame_point (p_cache, intercept, slope, voronoj,
952 triangle->cx, triangle->cy, mx, my, direct,
953 &x, &y);
954 triangle->x_3_1 = x;
955 triangle->y_3_1 = y;
956 }
957 }
958 return voronoj;
959 }
960
961 static void *
voronoj_export_common(const void * p_cache,void * p_voronoj,void * p_result,int only_edges)962 voronoj_export_common (const void *p_cache, void *p_voronoj, void *p_result,
963 int only_edges)
964 {
965 /* building the Geometry representing Voronoj */
966 gaiaGeomCollPtr result = (gaiaGeomCollPtr) p_result;
967 gaiaGeomCollPtr lines;
968 struct voronoj_aux *voronoj = (struct voronoj_aux *) p_voronoj;
969 int ind;
970 gaiaLinestringPtr ln;
971 struct voronoj_triangle *triangle;
972 struct voronoj_triangle *tri2;
973 int i;
974 int count;
975 double last;
976 double *array;
977
978 for (ind = 0; ind < voronoj->count; ind++)
979 {
980 triangle = voronoj->array + ind;
981
982 /* segments connecting two circumcenters */
983 if (triangle->tri_1_2 && triangle->trace_1_2)
984 {
985 tri2 = triangle->tri_1_2;
986 ln = gaiaAddLinestringToGeomColl (result, 2);
987 if (result->DimensionModel == GAIA_XY_Z)
988 {
989 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
990 triangle->cy, 0.0);
991 gaiaSetPointXYZ (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
992 }
993 else if (result->DimensionModel == GAIA_XY_M)
994 {
995 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
996 triangle->cy, 0.0);
997 gaiaSetPointXYM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
998 }
999 else if (result->DimensionModel == GAIA_XY_Z_M)
1000 {
1001 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1002 triangle->cy, 0.0, 0.0);
1003 gaiaSetPointXYZM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0,
1004 0.0);
1005 }
1006 else
1007 {
1008 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1009 gaiaSetPoint (ln->Coords, 1, tri2->cx, tri2->cy);
1010 }
1011 }
1012 if (triangle->tri_2_3 && triangle->trace_2_3)
1013 {
1014 tri2 = triangle->tri_2_3;
1015 ln = gaiaAddLinestringToGeomColl (result, 2);
1016 if (result->DimensionModel == GAIA_XY_Z)
1017 {
1018 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
1019 triangle->cy, 0.0);
1020 gaiaSetPointXYZ (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
1021 }
1022 else if (result->DimensionModel == GAIA_XY_M)
1023 {
1024 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
1025 triangle->cy, 0.0);
1026 gaiaSetPointXYM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
1027 }
1028 else if (result->DimensionModel == GAIA_XY_Z_M)
1029 {
1030 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1031 triangle->cy, 0.0, 0.0);
1032 gaiaSetPointXYZM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0,
1033 0.0);
1034 }
1035 else
1036 {
1037 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1038 gaiaSetPoint (ln->Coords, 1, tri2->cx, tri2->cy);
1039 }
1040 }
1041 if (triangle->tri_3_1 && triangle->trace_3_1)
1042 {
1043 tri2 = triangle->tri_3_1;
1044 ln = gaiaAddLinestringToGeomColl (result, 2);
1045 if (result->DimensionModel == GAIA_XY_Z)
1046 {
1047 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
1048 triangle->cy, 0.0);
1049 gaiaSetPointXYZ (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
1050 }
1051 else if (result->DimensionModel == GAIA_XY_M)
1052 {
1053 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
1054 triangle->cy, 0.0);
1055 gaiaSetPointXYM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0);
1056 }
1057 else if (result->DimensionModel == GAIA_XY_Z_M)
1058 {
1059 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1060 triangle->cy, 0.0, 0.0);
1061 gaiaSetPointXYZM (ln->Coords, 1, tri2->cx, tri2->cy, 0.0,
1062 0.0);
1063 }
1064 else
1065 {
1066 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1067 gaiaSetPoint (ln->Coords, 1, tri2->cx, tri2->cy);
1068 }
1069 }
1070
1071 /* segments connecting a circumcenter to the frame */
1072 if (triangle->tri_1_2 == NULL)
1073 {
1074 ln = gaiaAddLinestringToGeomColl (result, 2);
1075 if (result->DimensionModel == GAIA_XY_Z)
1076 {
1077 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
1078 triangle->cy, 0.0);
1079 gaiaSetPointXYZ (ln->Coords, 1, triangle->x_1_2,
1080 triangle->y_1_2, 0.0);
1081 }
1082 else if (result->DimensionModel == GAIA_XY_M)
1083 {
1084 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
1085 triangle->cy, 0.0);
1086 gaiaSetPointXYM (ln->Coords, 1, triangle->x_1_2,
1087 triangle->y_1_2, 0.0);
1088 }
1089 else if (result->DimensionModel == GAIA_XY_Z_M)
1090 {
1091 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1092 triangle->cy, 0.0, 0.0);
1093 gaiaSetPointXYZM (ln->Coords, 1, triangle->x_1_2,
1094 triangle->y_1_2, 0.0, 0.0);
1095 }
1096 else
1097 {
1098 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1099 gaiaSetPoint (ln->Coords, 1, triangle->x_1_2,
1100 triangle->y_1_2);
1101 }
1102 if (!only_edges)
1103 voronoj_add_frame_point (voronoj, triangle->x_1_2,
1104 triangle->y_1_2);
1105 }
1106 if (triangle->tri_2_3 == NULL)
1107 {
1108 ln = gaiaAddLinestringToGeomColl (result, 2);
1109 if (result->DimensionModel == GAIA_XY_Z)
1110 {
1111 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
1112 triangle->cy, 0.0);
1113 gaiaSetPointXYZ (ln->Coords, 1, triangle->x_2_3,
1114 triangle->y_2_3, 0.0);
1115 }
1116 else if (result->DimensionModel == GAIA_XY_M)
1117 {
1118 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
1119 triangle->cy, 0.0);
1120 gaiaSetPointXYM (ln->Coords, 1, triangle->x_2_3,
1121 triangle->y_2_3, 0.0);
1122 }
1123 else if (result->DimensionModel == GAIA_XY_Z_M)
1124 {
1125 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1126 triangle->cy, 0.0, 0.0);
1127 gaiaSetPointXYZM (ln->Coords, 1, triangle->x_2_3,
1128 triangle->y_2_3, 0.0, 0.0);
1129 }
1130 else
1131 {
1132 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1133 gaiaSetPoint (ln->Coords, 1, triangle->x_2_3,
1134 triangle->y_2_3);
1135 }
1136 if (!only_edges)
1137 voronoj_add_frame_point (voronoj, triangle->x_2_3,
1138 triangle->y_2_3);
1139 }
1140 if (triangle->tri_3_1 == NULL)
1141 {
1142 ln = gaiaAddLinestringToGeomColl (result, 2);
1143 if (result->DimensionModel == GAIA_XY_Z)
1144 {
1145 gaiaSetPointXYZ (ln->Coords, 0, triangle->cx,
1146 triangle->cy, 0.0);
1147 gaiaSetPointXYZ (ln->Coords, 1, triangle->x_3_1,
1148 triangle->y_3_1, 0.0);
1149 }
1150 else if (result->DimensionModel == GAIA_XY_M)
1151 {
1152 gaiaSetPointXYM (ln->Coords, 0, triangle->cx,
1153 triangle->cy, 0.0);
1154 gaiaSetPointXYM (ln->Coords, 1, triangle->x_3_1,
1155 triangle->y_3_1, 0.0);
1156 }
1157 else if (result->DimensionModel == GAIA_XY_Z_M)
1158 {
1159 gaiaSetPointXYZM (ln->Coords, 0, triangle->cx,
1160 triangle->cy, 0.0, 0.0);
1161 gaiaSetPointXYZM (ln->Coords, 1, triangle->x_3_1,
1162 triangle->y_3_1, 0.0, 0.0);
1163 }
1164 else
1165 {
1166 gaiaSetPoint (ln->Coords, 0, triangle->cx, triangle->cy);
1167 gaiaSetPoint (ln->Coords, 1, triangle->x_3_1,
1168 triangle->y_3_1);
1169 }
1170 if (!only_edges)
1171 voronoj_add_frame_point (voronoj, triangle->x_3_1,
1172 triangle->y_3_1);
1173 }
1174 }
1175
1176 if (only_edges)
1177 return result;
1178
1179 /* setting up the frame's upper edge */
1180 last = voronoj->minx;
1181 array = voronoj_sorted_up (voronoj, &count);
1182 if (array)
1183 {
1184 for (i = 0; i < count; i++)
1185 {
1186 ln = gaiaAddLinestringToGeomColl (result, 2);
1187 if (result->DimensionModel == GAIA_XY_Z)
1188 {
1189 gaiaSetPointXYZ (ln->Coords, 0, last, voronoj->maxy, 0.0);
1190 gaiaSetPointXYZ (ln->Coords, 1, *(array + i),
1191 voronoj->maxy, 0.0);
1192 }
1193 else if (result->DimensionModel == GAIA_XY_M)
1194 {
1195 gaiaSetPointXYM (ln->Coords, 0, last, voronoj->maxy, 0.0);
1196 gaiaSetPointXYM (ln->Coords, 1, *(array + i),
1197 voronoj->maxy, 0.0);
1198 }
1199 else if (result->DimensionModel == GAIA_XY_Z_M)
1200 {
1201 gaiaSetPointXYZM (ln->Coords, 0, last,
1202 voronoj->maxy, 0.0, 0.0);
1203 gaiaSetPointXYZM (ln->Coords, 1, *(array + i),
1204 voronoj->maxy, 0.0, 0.0);
1205 }
1206 else
1207 {
1208 gaiaSetPoint (ln->Coords, 0, last, voronoj->maxy);
1209 gaiaSetPoint (ln->Coords, 1, *(array + i), voronoj->maxy);
1210 }
1211 last = *(array + i);
1212 }
1213 free (array);
1214 }
1215
1216 /* closing the frame's upper edge */
1217 ln = gaiaAddLinestringToGeomColl (result, 2);
1218 if (result->DimensionModel == GAIA_XY_Z)
1219 {
1220 gaiaSetPointXYZ (ln->Coords, 0, last, voronoj->maxy, 0.0);
1221 gaiaSetPointXYZ (ln->Coords, 1, voronoj->maxx, voronoj->maxy, 0.0);
1222 }
1223 else if (result->DimensionModel == GAIA_XY_M)
1224 {
1225 gaiaSetPointXYM (ln->Coords, 0, last, voronoj->maxy, 0.0);
1226 gaiaSetPointXYM (ln->Coords, 1, voronoj->maxx, voronoj->maxy, 0.0);
1227 }
1228 else if (result->DimensionModel == GAIA_XY_Z_M)
1229 {
1230 gaiaSetPointXYZM (ln->Coords, 0, last, voronoj->maxy, 0.0, 0.0);
1231 gaiaSetPointXYZM (ln->Coords, 1, voronoj->maxx,
1232 voronoj->maxy, 0.0, 0.0);
1233 }
1234 else
1235 {
1236 gaiaSetPoint (ln->Coords, 0, last, voronoj->maxy);
1237 gaiaSetPoint (ln->Coords, 1, voronoj->maxx, voronoj->maxy);
1238 }
1239
1240 /* setting up the frame's lower edge */
1241 last = voronoj->minx;
1242 array = voronoj_sorted_low (voronoj, &count);
1243 if (array)
1244 {
1245 for (i = 0; i < count; i++)
1246 {
1247 ln = gaiaAddLinestringToGeomColl (result, 2);
1248 if (result->DimensionModel == GAIA_XY_Z)
1249 {
1250 gaiaSetPointXYZ (ln->Coords, 0, last, voronoj->miny, 0.0);
1251 gaiaSetPointXYZ (ln->Coords, 1, *(array + i),
1252 voronoj->miny, 0.0);
1253 }
1254 else if (result->DimensionModel == GAIA_XY_M)
1255 {
1256 gaiaSetPointXYM (ln->Coords, 0, last, voronoj->miny, 0.0);
1257 gaiaSetPointXYM (ln->Coords, 1, *(array + i),
1258 voronoj->miny, 0.0);
1259 }
1260 else if (result->DimensionModel == GAIA_XY_Z_M)
1261 {
1262 gaiaSetPointXYZM (ln->Coords, 0, last,
1263 voronoj->miny, 0.0, 0.0);
1264 gaiaSetPointXYZM (ln->Coords, 1, *(array + i),
1265 voronoj->miny, 0.0, 0.0);
1266 }
1267 else
1268 {
1269 gaiaSetPoint (ln->Coords, 0, last, voronoj->miny);
1270 gaiaSetPoint (ln->Coords, 1, *(array + i), voronoj->miny);
1271 }
1272 last = *(array + i);
1273 }
1274 free (array);
1275 }
1276
1277 /* closing the frame's lower edge */
1278 ln = gaiaAddLinestringToGeomColl (result, 2);
1279 if (result->DimensionModel == GAIA_XY_Z)
1280 {
1281 gaiaSetPointXYZ (ln->Coords, 0, last, voronoj->miny, 0.0);
1282 gaiaSetPointXYZ (ln->Coords, 1, voronoj->maxx, voronoj->miny, 0.0);
1283 }
1284 else if (result->DimensionModel == GAIA_XY_M)
1285 {
1286 gaiaSetPointXYM (ln->Coords, 0, last, voronoj->miny, 0.0);
1287 gaiaSetPointXYM (ln->Coords, 1, voronoj->maxx, voronoj->miny, 0.0);
1288 }
1289 else if (result->DimensionModel == GAIA_XY_Z_M)
1290 {
1291 gaiaSetPointXYZM (ln->Coords, 0, last, voronoj->miny, 0.0, 0.0);
1292 gaiaSetPointXYZM (ln->Coords, 1, voronoj->maxx,
1293 voronoj->miny, 0.0, 0.0);
1294 }
1295 else
1296 {
1297 gaiaSetPoint (ln->Coords, 0, last, voronoj->miny);
1298 gaiaSetPoint (ln->Coords, 1, voronoj->maxx, voronoj->miny);
1299 }
1300
1301 /* setting up the frame's left edge */
1302 last = voronoj->miny;
1303 array = voronoj_sorted_left (voronoj, &count);
1304 if (array)
1305 {
1306 for (i = 0; i < count; i++)
1307 {
1308 ln = gaiaAddLinestringToGeomColl (result, 2);
1309 if (result->DimensionModel == GAIA_XY_Z)
1310 {
1311 gaiaSetPointXYZ (ln->Coords, 0, voronoj->minx, last, 0.0);
1312 gaiaSetPointXYZ (ln->Coords, 1, voronoj->minx,
1313 *(array + i), 0.0);
1314 }
1315 else if (result->DimensionModel == GAIA_XY_M)
1316 {
1317 gaiaSetPointXYM (ln->Coords, 0, voronoj->minx, last, 0.0);
1318 gaiaSetPointXYM (ln->Coords, 1, voronoj->minx,
1319 *(array + i), 0.0);
1320 }
1321 else if (result->DimensionModel == GAIA_XY_Z_M)
1322 {
1323 gaiaSetPointXYZM (ln->Coords, 0, voronoj->minx,
1324 last, 0.0, 0.0);
1325 gaiaSetPointXYZM (ln->Coords, 1, voronoj->minx,
1326 *(array + i), 0.0, 0.0);
1327 }
1328 else
1329 {
1330 gaiaSetPoint (ln->Coords, 0, voronoj->minx, last);
1331 gaiaSetPoint (ln->Coords, 1, voronoj->minx, *(array + i));
1332 }
1333 last = *(array + i);
1334 }
1335 free (array);
1336 }
1337
1338 /* closing the frame's left edge */
1339 ln = gaiaAddLinestringToGeomColl (result, 2);
1340 if (result->DimensionModel == GAIA_XY_Z)
1341 {
1342 gaiaSetPointXYZ (ln->Coords, 0, voronoj->minx, last, 0.0);
1343 gaiaSetPointXYZ (ln->Coords, 1, voronoj->minx, voronoj->maxy, 0.0);
1344 }
1345 else if (result->DimensionModel == GAIA_XY_M)
1346 {
1347 gaiaSetPointXYM (ln->Coords, 0, voronoj->minx, last, 0.0);
1348 gaiaSetPointXYM (ln->Coords, 1, voronoj->minx, voronoj->maxy, 0.0);
1349 }
1350 else if (result->DimensionModel == GAIA_XY_Z_M)
1351 {
1352 gaiaSetPointXYZM (ln->Coords, 0, voronoj->minx, last, 0.0, 0.0);
1353 gaiaSetPointXYZM (ln->Coords, 1, voronoj->minx,
1354 voronoj->maxy, 0.0, 0.0);
1355 }
1356 else
1357 {
1358 gaiaSetPoint (ln->Coords, 0, voronoj->minx, last);
1359 gaiaSetPoint (ln->Coords, 1, voronoj->minx, voronoj->maxy);
1360 }
1361
1362 /* setting up the frame's right edge */
1363 last = voronoj->miny;
1364 array = voronoj_sorted_right (voronoj, &count);
1365 if (array)
1366 {
1367 for (i = 0; i < count; i++)
1368 {
1369 ln = gaiaAddLinestringToGeomColl (result, 2);
1370 if (result->DimensionModel == GAIA_XY_Z)
1371 {
1372 gaiaSetPointXYZ (ln->Coords, 0, voronoj->maxx, last, 0.0);
1373 gaiaSetPointXYZ (ln->Coords, 1, voronoj->maxx,
1374 *(array + i), 0.0);
1375 }
1376 else if (result->DimensionModel == GAIA_XY_M)
1377 {
1378 gaiaSetPointXYM (ln->Coords, 0, voronoj->maxx, last, 0.0);
1379 gaiaSetPointXYM (ln->Coords, 1, voronoj->maxx,
1380 *(array + i), 0.0);
1381 }
1382 else if (result->DimensionModel == GAIA_XY_Z_M)
1383 {
1384 gaiaSetPointXYZM (ln->Coords, 0, voronoj->maxx,
1385 last, 0.0, 0.0);
1386 gaiaSetPointXYZM (ln->Coords, 1, voronoj->maxx,
1387 *(array + i), 0.0, 0.0);
1388 }
1389 else
1390 {
1391 gaiaSetPoint (ln->Coords, 0, voronoj->maxx, last);
1392 gaiaSetPoint (ln->Coords, 1, voronoj->maxx, *(array + i));
1393 }
1394 last = *(array + i);
1395 }
1396 free (array);
1397 }
1398
1399 /* closing the frame's right edge */
1400 ln = gaiaAddLinestringToGeomColl (result, 2);
1401 if (result->DimensionModel == GAIA_XY_Z)
1402 {
1403 gaiaSetPointXYZ (ln->Coords, 0, voronoj->maxx, last, 0.0);
1404 gaiaSetPointXYZ (ln->Coords, 1, voronoj->maxx, voronoj->maxy, 0.0);
1405 }
1406 else if (result->DimensionModel == GAIA_XY_M)
1407 {
1408 gaiaSetPointXYM (ln->Coords, 0, voronoj->maxx, last, 0.0);
1409 gaiaSetPointXYM (ln->Coords, 1, voronoj->maxx, voronoj->maxy, 0.0);
1410 }
1411 else if (result->DimensionModel == GAIA_XY_Z_M)
1412 {
1413 gaiaSetPointXYZM (ln->Coords, 0, voronoj->maxx, last, 0.0, 0.0);
1414 gaiaSetPointXYZM (ln->Coords, 1, voronoj->maxx,
1415 voronoj->maxy, 0.0, 0.0);
1416 }
1417 else
1418 {
1419 gaiaSetPoint (ln->Coords, 0, voronoj->maxx, last);
1420 gaiaSetPoint (ln->Coords, 1, voronoj->maxx, voronoj->maxy);
1421 }
1422
1423 /* building Polygons */
1424 lines = result;
1425 if (p_cache != NULL)
1426 result = gaiaPolygonize_r (p_cache, lines, 1);
1427 else
1428 result = gaiaPolygonize (lines, 1);
1429 gaiaFreeGeomColl (lines);
1430 return result;
1431 }
1432
1433 SPATIALITE_PRIVATE void *
voronoj_export(void * p_voronoj,void * p_result,int only_edges)1434 voronoj_export (void *p_voronoj, void *p_result, int only_edges)
1435 {
1436 return voronoj_export_common (NULL, p_voronoj, p_result, only_edges);
1437 }
1438
1439 SPATIALITE_PRIVATE void *
voronoj_export_r(const void * p_cache,void * p_voronoj,void * p_result,int only_edges)1440 voronoj_export_r (const void *p_cache, void *p_voronoj, void *p_result,
1441 int only_edges)
1442 {
1443 return voronoj_export_common (p_cache, p_voronoj, p_result, only_edges);
1444 }
1445
1446 SPATIALITE_PRIVATE void
voronoj_free(void * p_voronoj)1447 voronoj_free (void *p_voronoj)
1448 {
1449 /* memory cleanup: destroying the Voronoj auxiliary struct */
1450 struct voronoj_aux *voronoj = (struct voronoj_aux *) p_voronoj;
1451 struct voronoj_point *pt;
1452 struct voronoj_point *ptn;
1453 free (voronoj->array);
1454 pt = voronoj->first_up;
1455 while (pt)
1456 {
1457 ptn = pt->next;
1458 free (pt);
1459 pt = ptn;
1460 }
1461 pt = voronoj->first_low;
1462 while (pt)
1463 {
1464 ptn = pt->next;
1465 free (pt);
1466 pt = ptn;
1467 }
1468 pt = voronoj->first_left;
1469 while (pt)
1470 {
1471 ptn = pt->next;
1472 free (pt);
1473 pt = ptn;
1474 }
1475 pt = voronoj->first_right;
1476 while (pt)
1477 {
1478 ptn = pt->next;
1479 free (pt);
1480 pt = ptn;
1481 }
1482 free (voronoj);
1483 }
1484
1485 #endif /* END GEOS_REENTRANT - GEOS >= 3.5.0 */
1486
1487 SPATIALITE_PRIVATE int
delaunay_triangle_check(void * ppg)1488 delaunay_triangle_check (void *ppg)
1489 {
1490 /* test if it's really a triangle */
1491 gaiaPolygonPtr pg = (gaiaPolygonPtr) ppg;
1492 gaiaRingPtr rng = pg->Exterior;
1493 if (rng->Points == 4 && pg->NumInteriors == 0)
1494 return 1;
1495 return 0;
1496 }
1497
1498 static void
concave_hull_stats(struct concave_hull_str * concave,double length)1499 concave_hull_stats (struct concave_hull_str *concave, double length)
1500 {
1501 /* update concave hull statistics */
1502 if (concave->count == 0)
1503 {
1504 concave->count = 1.0;
1505 concave->mean = length;
1506 return;
1507 }
1508
1509 concave->count += 1.0;
1510 concave->quot = concave->quot + (((concave->count - 1.0) *
1511 ((length - concave->mean) *
1512 (length - concave->mean))) /
1513 concave->count);
1514 concave->mean = concave->mean + ((length - concave->mean) / concave->count);
1515 }
1516
1517 static int
concave_hull_filter(const void * p_cache,double x1,double y1,double x2,double y2,double x3,double y3,double limit)1518 concave_hull_filter (const void *p_cache, double x1, double y1,
1519 double x2, double y2, double x3, double y3, double limit)
1520 {
1521 /* filtering triangles to be inserted into the Concave Hull */
1522 gaiaGeomCollPtr segm;
1523 gaiaLinestringPtr ln;
1524 double length;
1525
1526 segm = gaiaAllocGeomColl ();
1527 ln = gaiaAddLinestringToGeomColl (segm, 2);
1528 gaiaSetPoint (ln->Coords, 0, x1, y1);
1529 gaiaSetPoint (ln->Coords, 1, x2, y2);
1530 if (p_cache != NULL)
1531 gaiaGeomCollLength_r (p_cache, segm, &length);
1532 else
1533 gaiaGeomCollLength (segm, &length);
1534 gaiaFreeGeomColl (segm);
1535 if (length >= limit)
1536 return 0;
1537
1538 segm = gaiaAllocGeomColl ();
1539 ln = gaiaAddLinestringToGeomColl (segm, 2);
1540 gaiaSetPoint (ln->Coords, 0, x2, y2);
1541 gaiaSetPoint (ln->Coords, 1, x3, y3);
1542 if (p_cache != NULL)
1543 gaiaGeomCollLength_r (p_cache, segm, &length);
1544 else
1545 gaiaGeomCollLength (segm, &length);
1546 gaiaFreeGeomColl (segm);
1547 if (length >= limit)
1548 return 0;
1549
1550 segm = gaiaAllocGeomColl ();
1551 ln = gaiaAddLinestringToGeomColl (segm, 2);
1552 gaiaSetPoint (ln->Coords, 0, x3, y3);
1553 gaiaSetPoint (ln->Coords, 1, x1, y1);
1554 if (p_cache != NULL)
1555 gaiaGeomCollLength_r (p_cache, segm, &length);
1556 else
1557 gaiaGeomCollLength (segm, &length);
1558 gaiaFreeGeomColl (segm);
1559 if (length >= limit)
1560 return 0;
1561
1562 return 1;
1563 }
1564
1565 static gaiaGeomCollPtr
concave_hull_no_holes(gaiaGeomCollPtr in)1566 concave_hull_no_holes (gaiaGeomCollPtr in)
1567 {
1568 /* returning a Polygon surely not containing any hole */
1569 gaiaGeomCollPtr out = NULL;
1570 gaiaPolygonPtr pg_in;
1571 gaiaPolygonPtr pg_out;
1572 gaiaRingPtr rng_in;
1573 gaiaRingPtr rng_out;
1574 int iv;
1575 double x;
1576 double y;
1577 double z;
1578 double m;
1579
1580 if (in->DimensionModel == GAIA_XY_Z)
1581 out = gaiaAllocGeomCollXYZ ();
1582 else if (in->DimensionModel == GAIA_XY_M)
1583 out = gaiaAllocGeomCollXYM ();
1584 else if (in->DimensionModel == GAIA_XY_Z_M)
1585 out = gaiaAllocGeomCollXYZM ();
1586 else
1587 out = gaiaAllocGeomColl ();
1588 out->Srid = in->Srid;
1589
1590 pg_in = in->FirstPolygon;
1591 while (pg_in)
1592 {
1593 rng_in = pg_in->Exterior;
1594 pg_out = gaiaAddPolygonToGeomColl (out, rng_in->Points, 0);
1595 rng_out = pg_out->Exterior;
1596 for (iv = 0; iv < rng_in->Points; iv++)
1597 {
1598 /* copying Exterior Ring vertices */
1599 if (in->DimensionModel == GAIA_XY_Z)
1600 {
1601 gaiaGetPointXYZ (rng_in->Coords, iv, &x, &y, &z);
1602 gaiaSetPointXYZ (rng_out->Coords, iv, x, y, z);
1603 }
1604 else if (in->DimensionModel == GAIA_XY_M)
1605 {
1606 gaiaGetPointXYM (rng_in->Coords, iv, &x, &y, &m);
1607 gaiaSetPointXYM (rng_out->Coords, iv, x, y, m);
1608 }
1609 else if (in->DimensionModel == GAIA_XY_Z_M)
1610 {
1611 gaiaGetPointXYZM (rng_in->Coords, iv, &x, &y, &z, &m);
1612 gaiaSetPointXYZM (rng_out->Coords, iv, x, y, z, m);
1613 }
1614 else
1615 {
1616 gaiaGetPoint (rng_in->Coords, iv, &x, &y);
1617 gaiaSetPoint (rng_out->Coords, iv, x, y);
1618 }
1619 }
1620 pg_in = pg_in->Next;
1621 }
1622 return out;
1623 }
1624
1625 static void *
concave_hull_build_common(const void * p_cache,void * p_first,int dimension_model,double factor,int allow_holes)1626 concave_hull_build_common (const void *p_cache, void *p_first,
1627 int dimension_model, double factor, int allow_holes)
1628 {
1629 /* building the Concave Hull */
1630 struct concave_hull_str concave;
1631 gaiaPolygonPtr first = (gaiaPolygonPtr) p_first;
1632 gaiaPolygonPtr pg;
1633 gaiaRingPtr rng;
1634 gaiaPolygonPtr pg_out;
1635 gaiaRingPtr rng_out;
1636 gaiaGeomCollPtr segm;
1637 gaiaGeomCollPtr result;
1638 gaiaLinestringPtr ln;
1639 double x;
1640 double y;
1641 double z;
1642 double m;
1643 double x1;
1644 double y1;
1645 double x2;
1646 double y2;
1647 double x3;
1648 double y3;
1649 double length;
1650 double std_dev;
1651 int count;
1652
1653 /* initializing the struct for mean and standard deviation */
1654 concave.mean = 0.0;
1655 concave.quot = 0.0;
1656 concave.count = 0.0;
1657
1658 pg = first;
1659 while (pg)
1660 {
1661 /* examining each triangle / computing statistics distribution */
1662 rng = pg->Exterior;
1663 if (pg->DimensionModel == GAIA_XY_Z)
1664 {
1665 gaiaGetPointXYZ (rng->Coords, 0, &x, &y, &z);
1666 x1 = x;
1667 y1 = y;
1668 gaiaGetPointXYZ (rng->Coords, 1, &x, &y, &z);
1669 x2 = x;
1670 y2 = y;
1671 gaiaGetPointXYZ (rng->Coords, 2, &x, &y, &z);
1672 x3 = x;
1673 y3 = y;
1674 }
1675 else if (pg->DimensionModel == GAIA_XY_M)
1676 {
1677 gaiaGetPointXYM (rng->Coords, 0, &x, &y, &m);
1678 x1 = x;
1679 y1 = y;
1680 gaiaGetPointXYM (rng->Coords, 1, &x, &y, &m);
1681 x2 = x;
1682 y2 = y;
1683 gaiaGetPointXYM (rng->Coords, 2, &x, &y, &m);
1684 x3 = x;
1685 y3 = y;
1686 }
1687 else if (pg->DimensionModel == GAIA_XY_Z_M)
1688 {
1689 gaiaGetPointXYZM (rng->Coords, 0, &x, &y, &z, &m);
1690 x1 = x;
1691 y1 = y;
1692 gaiaGetPointXYZM (rng->Coords, 1, &x, &y, &z, &m);
1693 x2 = x;
1694 y2 = y;
1695 gaiaGetPointXYZM (rng->Coords, 2, &x, &y, &z, &m);
1696 x3 = x;
1697 y3 = y;
1698 }
1699 else
1700 {
1701 gaiaGetPoint (rng->Coords, 0, &x, &y);
1702 x1 = x;
1703 y1 = y;
1704 gaiaGetPoint (rng->Coords, 1, &x, &y);
1705 x2 = x;
1706 y2 = y;
1707 gaiaGetPoint (rng->Coords, 2, &x, &y);
1708 x3 = x;
1709 y3 = y;
1710 }
1711
1712 segm = gaiaAllocGeomColl ();
1713 ln = gaiaAddLinestringToGeomColl (segm, 2);
1714 gaiaSetPoint (ln->Coords, 0, x1, y1);
1715 gaiaSetPoint (ln->Coords, 1, x2, y2);
1716 if (p_cache != NULL)
1717 gaiaGeomCollLength_r (p_cache, segm, &length);
1718 else
1719 gaiaGeomCollLength (segm, &length);
1720 gaiaFreeGeomColl (segm);
1721 concave_hull_stats (&concave, length);
1722
1723 segm = gaiaAllocGeomColl ();
1724 ln = gaiaAddLinestringToGeomColl (segm, 2);
1725 gaiaSetPoint (ln->Coords, 0, x2, y2);
1726 gaiaSetPoint (ln->Coords, 1, x3, y3);
1727 if (p_cache != NULL)
1728 gaiaGeomCollLength_r (p_cache, segm, &length);
1729 else
1730 gaiaGeomCollLength (segm, &length);
1731 gaiaFreeGeomColl (segm);
1732 concave_hull_stats (&concave, length);
1733
1734 segm = gaiaAllocGeomColl ();
1735 ln = gaiaAddLinestringToGeomColl (segm, 2);
1736 gaiaSetPoint (ln->Coords, 0, x3, y3);
1737 gaiaSetPoint (ln->Coords, 1, x1, y1);
1738 if (p_cache != NULL)
1739 gaiaGeomCollLength_r (p_cache, segm, &length);
1740 else
1741 gaiaGeomCollLength (segm, &length);
1742 gaiaFreeGeomColl (segm);
1743 concave_hull_stats (&concave, length);
1744
1745 pg = pg->Next;
1746 }
1747
1748 std_dev = sqrt (concave.quot / concave.count);
1749
1750 /* creating the Geometry representing the Concave Hull */
1751 if (dimension_model == GAIA_XY_Z)
1752 result = gaiaAllocGeomCollXYZ ();
1753 else if (dimension_model == GAIA_XY_M)
1754 result = gaiaAllocGeomCollXYM ();
1755 else if (dimension_model == GAIA_XY_Z_M)
1756 result = gaiaAllocGeomCollXYZM ();
1757 else
1758 result = gaiaAllocGeomColl ();
1759
1760 count = 0;
1761 pg = first;
1762 while (pg)
1763 {
1764 /* selecting triangles to be inserted */
1765 rng = pg->Exterior;
1766 if (pg->DimensionModel == GAIA_XY_Z)
1767 {
1768 gaiaGetPointXYZ (rng->Coords, 0, &x, &y, &z);
1769 x1 = x;
1770 y1 = y;
1771 gaiaGetPointXYZ (rng->Coords, 1, &x, &y, &z);
1772 x2 = x;
1773 y2 = y;
1774 gaiaGetPointXYZ (rng->Coords, 2, &x, &y, &z);
1775 x3 = x;
1776 y3 = y;
1777 }
1778 else if (pg->DimensionModel == GAIA_XY_M)
1779 {
1780 gaiaGetPointXYM (rng->Coords, 0, &x, &y, &m);
1781 x1 = x;
1782 y1 = y;
1783 gaiaGetPointXYM (rng->Coords, 1, &x, &y, &m);
1784 x2 = x;
1785 y2 = y;
1786 gaiaGetPointXYM (rng->Coords, 2, &x, &y, &m);
1787 x3 = x;
1788 y3 = y;
1789 }
1790 else if (pg->DimensionModel == GAIA_XY_Z_M)
1791 {
1792 gaiaGetPointXYZM (rng->Coords, 0, &x, &y, &z, &m);
1793 x1 = x;
1794 y1 = y;
1795 gaiaGetPointXYZM (rng->Coords, 1, &x, &y, &z, &m);
1796 x2 = x;
1797 y2 = y;
1798 gaiaGetPointXYZM (rng->Coords, 2, &x, &y, &z, &m);
1799 x3 = x;
1800 y3 = y;
1801 }
1802 else
1803 {
1804 gaiaGetPoint (rng->Coords, 0, &x, &y);
1805 x1 = x;
1806 y1 = y;
1807 gaiaGetPoint (rng->Coords, 1, &x, &y);
1808 x2 = x;
1809 y2 = y;
1810 gaiaGetPoint (rng->Coords, 2, &x, &y);
1811 x3 = x;
1812 y3 = y;
1813 }
1814
1815 if (concave_hull_filter
1816 (p_cache, x1, y1, x2, y2, x3, y3, std_dev * factor))
1817 {
1818 /* inserting this triangle into the Concave Hull */
1819 pg_out = gaiaAddPolygonToGeomColl (result, 4, 0);
1820 rng_out = pg_out->Exterior;
1821 if (pg->DimensionModel == GAIA_XY_Z)
1822 {
1823 gaiaGetPointXYZ (rng->Coords, 0, &x, &y, &z);
1824 gaiaSetPointXYZ (rng_out->Coords, 0, x, y, z);
1825 gaiaGetPointXYZ (rng->Coords, 1, &x, &y, &z);
1826 gaiaSetPointXYZ (rng_out->Coords, 1, x, y, z);
1827 gaiaGetPointXYZ (rng->Coords, 2, &x, &y, &z);
1828 gaiaSetPointXYZ (rng_out->Coords, 2, x, y, z);
1829 gaiaGetPointXYZ (rng->Coords, 3, &x, &y, &z);
1830 gaiaSetPointXYZ (rng_out->Coords, 3, x, y, z);
1831 }
1832 else if (pg->DimensionModel == GAIA_XY_M)
1833 {
1834 gaiaGetPointXYM (rng->Coords, 0, &x, &y, &m);
1835 gaiaSetPointXYM (rng_out->Coords, 0, x, y, m);
1836 gaiaGetPointXYM (rng->Coords, 1, &x, &y, &m);
1837 gaiaSetPointXYM (rng_out->Coords, 1, x, y, m);
1838 gaiaGetPointXYM (rng->Coords, 2, &x, &y, &m);
1839 gaiaSetPointXYM (rng_out->Coords, 2, x, y, m);
1840 gaiaGetPointXYM (rng->Coords, 3, &x, &y, &m);
1841 gaiaSetPointXYM (rng_out->Coords, 3, x, y, m);
1842 }
1843 else if (pg->DimensionModel == GAIA_XY_Z_M)
1844 {
1845 gaiaGetPointXYZM (rng->Coords, 0, &x, &y, &z, &m);
1846 gaiaSetPointXYZM (rng_out->Coords, 0, x, y, z, m);
1847 gaiaGetPointXYZM (rng->Coords, 1, &x, &y, &z, &m);
1848 gaiaSetPointXYZM (rng_out->Coords, 1, x, y, z, m);
1849 gaiaGetPointXYZM (rng->Coords, 2, &x, &y, &z, &m);
1850 gaiaSetPointXYZM (rng_out->Coords, 2, x, y, z, m);
1851 gaiaGetPointXYZM (rng->Coords, 3, &x, &y, &z, &m);
1852 gaiaSetPointXYZM (rng_out->Coords, 3, x, y, z, m);
1853 }
1854 else
1855 {
1856 gaiaGetPoint (rng->Coords, 0, &x, &y);
1857 gaiaSetPoint (rng_out->Coords, 0, x, y);
1858 gaiaGetPoint (rng->Coords, 1, &x, &y);
1859 gaiaSetPoint (rng_out->Coords, 1, x, y);
1860 gaiaGetPoint (rng->Coords, 2, &x, &y);
1861 gaiaSetPoint (rng_out->Coords, 2, x, y);
1862 gaiaGetPoint (rng->Coords, 3, &x, &y);
1863 gaiaSetPoint (rng_out->Coords, 3, x, y);
1864 }
1865 count++;
1866 }
1867
1868 pg = pg->Next;
1869 }
1870
1871 if (count == 0)
1872 {
1873 gaiaFreeGeomColl (result);
1874 return NULL;
1875 }
1876
1877 /* merging all triangles into the Concave Hull */
1878 segm = result;
1879 if (p_cache != NULL)
1880 result = gaiaUnaryUnion_r (p_cache, segm);
1881 else
1882 result = gaiaUnaryUnion (segm);
1883 gaiaFreeGeomColl (segm);
1884 if (!result)
1885 return NULL;
1886 if (result->FirstPolygon == NULL)
1887 {
1888 gaiaFreeGeomColl (result);
1889 return NULL;
1890 }
1891 if (allow_holes)
1892 return result;
1893
1894 /* suppressing any interior hole */
1895 segm = result;
1896 result = concave_hull_no_holes (segm);
1897 gaiaFreeGeomColl (segm);
1898 if (!result)
1899 return NULL;
1900 if (result->FirstPolygon == NULL)
1901 {
1902 gaiaFreeGeomColl (result);
1903 return NULL;
1904 }
1905 return result;
1906 }
1907
1908 SPATIALITE_PRIVATE void *
concave_hull_build(void * p_first,int dimension_model,double factor,int allow_holes)1909 concave_hull_build (void *p_first, int dimension_model, double factor,
1910 int allow_holes)
1911 {
1912 return concave_hull_build_common (NULL, p_first, dimension_model, factor,
1913 allow_holes);
1914 }
1915
1916 SPATIALITE_PRIVATE void *
concave_hull_build_r(const void * p_cache,void * p_first,int dimension_model,double factor,int allow_holes)1917 concave_hull_build_r (const void *p_cache, void *p_first, int dimension_model,
1918 double factor, int allow_holes)
1919 {
1920 return concave_hull_build_common (p_cache, p_first, dimension_model, factor,
1921 allow_holes);
1922 }
1923
1924 #endif /* end GEOS advanced features */
1925