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