1 /*
2
3 MIT License
4
5 DELABELLA - Delaunay triangulation library
6 Copyright (C) 2018 GUMIX - Marcin Sokalski
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this software and associated documentation files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in all
16 copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24 SOFTWARE.
25
26 */
27
28 #include <assert.h>
29 #include <stdio.h>
30
31 // gcc 4.9 for Android doesn't have search.h
32 #if !defined(__ANDROID__) || defined(__clang__)
33 #include <search.h>
34 #endif
35
36 #if (defined(__APPLE__))
37 #include <malloc/malloc.h>
38 #else
39 #include <stdlib.h>
40 #endif
41
42 #include <algorithm>
43 #include "delabella.pxx" // we just need LOG() macro
44
45 // assuming BITS is max(X_BITS,Y_BITS)
46
47 typedef double Signed14; // BITS xy coords
48 typedef double Signed15; // BITS + 1 vect::xy
49 typedef long double Unsigned28; // 2xBITS z coord
50 typedef long double Signed29; // 2xBITS + 1 vect::z
51 typedef long double Signed31; // 2xBITS + 3 norm::z
52 typedef long double Signed45; // 3xBITS + 3 norm::xy
53 typedef long double Signed62; // 4xBITS + 6 dot(vect,norm)
54
55 /*
56 // above typedefs can be used to configure delabella arithmetic types
57 // in example, EXACT SOLVER (with xy coords limited to 14bits integers in range: +/-8192)
58 // could be achieved with following configuration:
59
60 typedef int16_t Signed14; // BITS xy coords
61 typedef int16_t Signed15; // BITS + 1 vect::xy
62 typedef uint32_t Unsigned28; // 2xBITS z coord
63 typedef int32_t Signed29; // 2xBITS + 1 vect::z
64 typedef int32_t Signed31; // 2xBITS + 3 norm::z
65 typedef int64_t Signed45; // 3xBITS + 3 norm::xy
66 typedef int64_t Signed62; // 4xBITS + 6 dot(vect,norm)
67
68 // alternatively, one could use some BigInt implementation
69 // in order to expand xy range
70 */
71
72
s14sqr(const Signed14 & s)73 static Unsigned28 s14sqr(const Signed14& s)
74 {
75 return (Unsigned28)((Signed29)s*s);
76 }
77
78 struct Norm
79 {
80 Signed45 x;
81 Signed45 y;
82 Signed31 z;
83 };
84
85 struct Vect
86 {
87 Signed15 x, y;
88 Signed29 z;
89
crossVect90 Norm cross (const Vect& v) const // cross prod
91 {
92 Norm n;
93 n.x = (Signed45)y*v.z - (Signed45)z*v.y;
94 n.y = (Signed45)z*v.x - (Signed45)x*v.z;
95 n.z = (Signed29)x*v.y - (Signed29)y*v.x;
96 return n;
97 }
98 };
99
100 struct CDelaBella : IDelaBella
101 {
102 struct Face;
103
104 struct Vert : DelaBella_Vertex
105 {
106 Unsigned28 z;
107 Face* sew;
108
operator -CDelaBella::Vert109 Vect operator - (const Vert& v) const // diff
110 {
111 Vect d;
112 d.x = (Signed15)x - (Signed15)v.x;
113 d.y = (Signed15)y - (Signed15)v.y;
114 d.z = (Signed29)z - (Signed29)v.z;
115 return d;
116 }
117
overlapCDelaBella::Vert118 static bool overlap(const Vert* v1, const Vert* v2)
119 {
120 return v1->x == v2->x && v1->y == v2->y;
121 }
122
operator <CDelaBella::Vert123 bool operator < (const Vert& v) const
124 {
125 return u28cmp(this, &v) < 0;
126 }
127
u28cmpCDelaBella::Vert128 static int u28cmp(const void* a, const void* b)
129 {
130 Vert* va = (Vert*)a;
131 Vert* vb = (Vert*)b;
132 if (va->z < vb->z)
133 return -1;
134 if (va->z > vb->z)
135 return 1;
136 if (va->y < vb->y)
137 return -1;
138 if (va->y > vb->y)
139 return 1;
140 if (va->x < vb->x)
141 return -1;
142 if (va->x > vb->x)
143 return 1;
144 if (va->i < vb->i)
145 return -1;
146 if (va->i > vb->i)
147 return 1;
148 return 0;
149 }
150 };
151
152 struct Face : DelaBella_Triangle
153 {
154 Norm n;
155
AllocCDelaBella::Face156 static Face* Alloc(Face** from)
157 {
158 Face* f = *from;
159 *from = (Face*)f->next;
160 f->next = 0;
161 return f;
162 }
163
FreeCDelaBella::Face164 void Free(Face** to)
165 {
166 next = *to;
167 *to = this;
168 }
169
NextCDelaBella::Face170 Face* Next(const Vert* p) const
171 {
172 // return next face in same direction as face vertices are (cw/ccw)
173
174 if (v[0] == p)
175 return (Face*)f[1];
176 if (v[1] == p)
177 return (Face*)f[2];
178 if (v[2] == p)
179 return (Face*)f[0];
180 return 0;
181 }
182
dotCDelaBella::Face183 Signed62 dot(const Vert& p) const // dot
184 {
185 Vect d = p - *(Vert*)v[0];
186 return (Signed62)n.x * d.x + (Signed62)n.y * d.y + (Signed62)n.z * d.z;
187 }
188
crossCDelaBella::Face189 Norm cross() const // cross of diffs
190 {
191 return (*(Vert*)v[1] - *(Vert*)v[0]).cross(*(Vert*)v[2] - *(Vert*)v[0]);
192 }
193 };
194
195 Vert* vert_alloc;
196 Face* face_alloc;
197 int max_verts;
198 int max_faces;
199
200 Face* first_dela_face;
201 Face* first_hull_face;
202 Vert* first_hull_vert;
203
204 int inp_verts;
205 int out_verts;
206
207 int(*errlog_proc)(void* file, const char* fmt, ...);
208 void* errlog_file;
209
~CDelaBellaCDelaBella210 virtual ~CDelaBella ()
211 {
212 }
213
TriangulateCDelaBella214 int Triangulate()
215 {
216 int points = inp_verts;
217 std::sort(vert_alloc, vert_alloc + points);
218
219 // remove dups
220 {
221 int w = 0, r = 1; // skip initial no-dups block
222 while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + w))
223 {
224 w++;
225 r++;
226 }
227
228 w++;
229
230 while (r < points)
231 {
232 r++;
233
234 // skip dups
235 while (r < points && Vert::overlap(vert_alloc + r, vert_alloc + r - 1))
236 r++;
237
238 // copy next no-dups block
239 while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + r - 1))
240 vert_alloc[w++] = vert_alloc[r++];
241 }
242
243 if (points - w)
244 {
245 if (errlog_proc)
246 errlog_proc(errlog_file, "[WRN] detected %d dups in xy array!\n", points - w);
247 points = w;
248 }
249 }
250
251 if (points < 3)
252 {
253 if (points == 2)
254 {
255 if (errlog_proc)
256 errlog_proc(errlog_file, "[WRN] all input points are colinear, returning single segment!\n");
257 first_hull_vert = vert_alloc + 0;
258 vert_alloc[0].next = (DelaBella_Vertex*)vert_alloc + 1;
259 vert_alloc[1].next = 0;
260 }
261 else
262 {
263 if (errlog_proc)
264 errlog_proc(errlog_file, "[WRN] all input points are identical, returning single point!\n");
265 first_hull_vert = vert_alloc + 0;
266 vert_alloc[0].next = 0;
267 }
268
269 return -points;
270 }
271
272 int hull_faces = 2 * points - 4;
273
274 if (max_faces < hull_faces)
275 {
276 if (max_faces)
277 free(face_alloc);
278 max_faces = 0;
279 face_alloc = (Face*)malloc(sizeof(Face) * hull_faces);
280 if (face_alloc)
281 max_faces = hull_faces;
282 else
283 {
284 if (errlog_proc)
285 errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
286 return 0;
287 }
288 }
289
290 for (int i = 1; i < hull_faces; i++)
291 face_alloc[i - 1].next = face_alloc + i;
292 face_alloc[hull_faces - 1].next = 0;
293
294 Face* cache = face_alloc;
295 Face* hull = 0;
296
297 Face f; // tmp
298 f.v[0] = vert_alloc + 0;
299 f.v[1] = vert_alloc + 1;
300 f.v[2] = vert_alloc + 2;
301 f.n = f.cross();
302
303 bool colinear = f.n.z == 0;
304 int i = 3;
305
306 /////////////////////////////////////////////////////////////////////////
307 // UNTIL INPUT IS COPLANAR, GROW IT IN FORM OF A 2D CONTOUR
308 /*
309 . | | after adding . | ________* L
310 . \ Last points to / Head next point . \ ______/ /
311 . *____ | -----> .H *____ |
312 . |\_ \_____ | . |\_ \_____ |
313 . \ \_ \__* - Tail points to Last . \ \_ \__* T
314 . \ \_ / . \ \_ /
315 . \__ \_ __/ . \__ \_ __/
316 . \__* - Head points to Tail . \__/
317 */
318
319 Vert* head = (Vert*)f.v[0];
320 Vert* tail = (Vert*)f.v[1];
321 Vert* last = (Vert*)f.v[2];
322
323 head->next = tail;
324 tail->next = last;
325 last->next = head;
326
327 while (i < points && f.dot(vert_alloc[i]) == 0)
328 {
329 Vert* v = vert_alloc + i;
330
331 // it is enough to test just 1 non-zero coord
332 // but we want also to test stability (assert)
333 // so we calc all signs...
334
335 // why not testing sign of dot prod of 2 normals?
336 // that way we'd fall into precision problems
337
338 Norm LvH = (*v - *last).cross(*head - *last);
339 bool lvh =
340 (f.n.x > 0 && LvH.x > 0) || (f.n.x < 0 && LvH.x < 0) ||
341 (f.n.y > 0 && LvH.y > 0) || (f.n.y < 0 && LvH.y < 0) ||
342 (f.n.z > 0 && LvH.z > 0) || (f.n.z < 0 && LvH.z < 0);
343
344 Norm TvL = (*v - *tail).cross(*last - *tail);
345 bool tvl =
346 (f.n.x > 0 && TvL.x > 0) || (f.n.x < 0 && TvL.x < 0) ||
347 (f.n.y > 0 && TvL.y > 0) || (f.n.y < 0 && TvL.y < 0) ||
348 (f.n.z > 0 && TvL.z > 0) || (f.n.z < 0 && TvL.z < 0);
349
350 if (lvh && !tvl) // insert new f on top of e(2,0) = (last,head)
351 {
352 // f.v[0] = head;
353 f.v[1] = last;
354 f.v[2] = v;
355
356 last->next = v;
357 v->next = head;
358 tail = last;
359 }
360 else
361 if (tvl && !lvh) // insert new f on top of e(1,2) = (tail,last)
362 {
363 f.v[0] = last;
364 //f.v[1] = tail;
365 f.v[2] = v;
366
367 tail->next = v;
368 v->next = last;
369 head = last;
370 }
371 else
372 {
373 // wtf? dilithium crystals are fucked.
374 assert(0);
375 }
376
377 last = v;
378 i++;
379 }
380
381 if (i == points)
382 {
383 if (colinear)
384 {
385 if (errlog_proc)
386 errlog_proc(errlog_file, "[WRN] all input points are colinear, returning segment list!\n");
387 first_hull_vert = head;
388 last->next = 0; // break contour, make it a list
389 return -points;
390 }
391 else
392 {
393 if (points > 3)
394 {
395 if (errlog_proc)
396 errlog_proc(errlog_file, "[NFO] all input points are cocircular.\n");
397 }
398 else
399 {
400 if (errlog_proc)
401 errlog_proc(errlog_file, "[NFO] trivial case of 3 points, thank you.\n");
402
403 first_dela_face = Face::Alloc(&cache);
404 first_dela_face->next = 0;
405 first_hull_face = Face::Alloc(&cache);
406 first_hull_face->next = 0;
407
408 first_dela_face->f[0] = first_dela_face->f[1] = first_dela_face->f[2] = first_hull_face;
409 first_hull_face->f[0] = first_hull_face->f[1] = first_hull_face->f[2] = first_dela_face;
410
411 head->sew = tail->sew = last->sew = first_hull_face;
412
413 if (f.n.z < 0)
414 {
415 first_dela_face->v[0] = head;
416 first_dela_face->v[1] = tail;
417 first_dela_face->v[2] = last;
418 first_hull_face->v[0] = last;
419 first_hull_face->v[1] = tail;
420 first_hull_face->v[2] = head;
421
422 // reverse silhouette
423 head->next = last;
424 last->next = tail;
425 tail->next = head;
426
427 first_hull_vert = last;
428 }
429 else
430 {
431 first_dela_face->v[0] = last;
432 first_dela_face->v[1] = tail;
433 first_dela_face->v[2] = head;
434 first_hull_face->v[0] = head;
435 first_hull_face->v[1] = tail;
436 first_hull_face->v[2] = last;
437
438 first_hull_vert = head;
439 }
440
441 first_dela_face->n = first_dela_face->cross();
442 first_hull_face->n = first_hull_face->cross();
443
444 return 3;
445 }
446
447 // retract last point it will be added as a cone's top later
448 last = head;
449 head = (Vert*)head->next;
450 i--;
451 }
452 }
453
454 /////////////////////////////////////////////////////////////////////////
455 // CREATE CONE HULL WITH TOP AT cloud[i] AND BASE MADE OF CONTOUR LIST
456 // in 2 ways :( - depending on at which side of the contour a top vertex appears
457
458 Vert* q = vert_alloc + i;
459
460 if (f.dot(*q) > 0)
461 {
462 Vert* p = last;
463 Vert* n = (Vert*)p->next;
464
465 Face* first_side = Face::Alloc(&cache);
466 first_side->v[0] = p;
467 first_side->v[1] = n;
468 first_side->v[2] = q;
469 first_side->n = first_side->cross();
470 hull = first_side;
471
472 p = n;
473 n = (Vert*)n->next;
474
475 Face* prev_side = first_side;
476 Face* prev_base = 0;
477 Face* first_base = 0;
478
479 do
480 {
481 Face* base = Face::Alloc(&cache);
482 base->v[0] = n;
483 base->v[1] = p;
484 base->v[2] = last;
485 base->n = base->cross();
486
487 Face* side = Face::Alloc(&cache);
488 side->v[0] = p;
489 side->v[1] = n;
490 side->v[2] = q;
491 side->n = side->cross();
492
493 side->f[2] = base;
494 base->f[2] = side;
495
496 side->f[1] = prev_side;
497 prev_side->f[0] = side;
498
499 base->f[0] = prev_base;
500 if (prev_base)
501 prev_base->f[1] = base;
502 else
503 first_base = base;
504
505 prev_base = base;
506 prev_side = side;
507
508 p = n;
509 n = (Vert*)n->next;
510 } while (n != last);
511
512 Face* last_side = Face::Alloc(&cache);
513 last_side->v[0] = p;
514 last_side->v[1] = n;
515 last_side->v[2] = q;
516 last_side->n = last_side->cross();
517
518 last_side->f[1] = prev_side;
519 prev_side->f[0] = last_side;
520
521 last_side->f[0] = first_side;
522 first_side->f[1] = last_side;
523
524 first_base->f[0] = first_side;
525 first_side->f[2] = first_base;
526
527 last_side->f[2] = prev_base;
528 prev_base->f[1] = last_side;
529
530 i++;
531 }
532 else
533 {
534 Vert* p = last;
535 Vert* n = (Vert*)p->next;
536
537 Face* first_side = Face::Alloc(&cache);
538 first_side->v[0] = n;
539 first_side->v[1] = p;
540 first_side->v[2] = q;
541 first_side->n = first_side->cross();
542 hull = first_side;
543
544 p = n;
545 n = (Vert*)n->next;
546
547 Face* prev_side = first_side;
548 Face* prev_base = 0;
549 Face* first_base = 0;
550
551 do
552 {
553 Face* base = Face::Alloc(&cache);
554 base->v[0] = p;
555 base->v[1] = n;
556 base->v[2] = last;
557 base->n = base->cross();
558
559 Face* side = Face::Alloc(&cache);
560 side->v[0] = n;
561 side->v[1] = p;
562 side->v[2] = q;
563 side->n = side->cross();
564
565 side->f[2] = base;
566 base->f[2] = side;
567
568 side->f[0] = prev_side;
569 prev_side->f[1] = side;
570
571 base->f[1] = prev_base;
572 if (prev_base)
573 prev_base->f[0] = base;
574 else
575 first_base = base;
576
577 prev_base = base;
578 prev_side = side;
579
580 p = n;
581 n = (Vert*)n->next;
582 } while (n != last);
583
584 Face* last_side = Face::Alloc(&cache);
585 last_side->v[0] = n;
586 last_side->v[1] = p;
587 last_side->v[2] = q;
588 last_side->n = last_side->cross();
589
590 last_side->f[0] = prev_side;
591 prev_side->f[1] = last_side;
592
593 last_side->f[1] = first_side;
594 first_side->f[0] = last_side;
595
596 first_base->f[1] = first_side;
597 first_side->f[2] = first_base;
598
599 last_side->f[2] = prev_base;
600 prev_base->f[0] = last_side;
601
602 i++;
603 }
604
605 /////////////////////////////////////////////////////////////////////////
606 // ACTUAL ALGORITHM
607
608 for (; i < points; i++)
609 {
610 //ValidateHull(alloc, 2 * i - 4);
611
612 Vert* _q = vert_alloc + i;
613 Vert* _p = vert_alloc + i - 1;
614 Face* _f = hull;
615
616 // 1. FIND FIRST VISIBLE FACE
617 // simply iterate around last vertex using last added triangle adjecency info
618 while (_f->dot(*_q) <= 0)
619 {
620 _f = _f->Next(_p);
621 if (_f == hull)
622 {
623 // if no visible face can be located at last vertex,
624 // let's run through all faces (approximately last to first),
625 // yes this is emergency fallback and should not ever happen.
626 _f = face_alloc + 2 * i - 4 - 1;
627 while (_f->dot(*_q) <= 0)
628 {
629 assert(_f != face_alloc); // no face is visible? you must be kidding!
630 _f--;
631 }
632 }
633 }
634
635 // 2. DELETE VISIBLE FACES & ADD NEW ONES
636 // (we also build silhouette (vertex loop) between visible & invisible faces)
637
638 int del = 0;
639 int add = 0;
640
641 // push first visible face onto stack (of visible faces)
642 Face* stack = _f;
643 _f->next = _f; // old trick to use list pointers as 'on-stack' markers
644 while (stack)
645 {
646 // pop, take care of last item ptr (it's not null!)
647 _f = stack;
648 stack = (Face*)_f->next;
649 if (stack == _f)
650 stack = 0;
651 _f->next = 0;
652
653 // copy parts of old face that we still need after removal
654 Vert* fv[3] = { (Vert*)_f->v[0],(Vert*)_f->v[1],(Vert*)_f->v[2] };
655 Face* ff[3] = { (Face*)_f->f[0],(Face*)_f->f[1],(Face*)_f->f[2] };
656
657 // delete visible face
658 _f->Free(&cache);
659 del++;
660
661 // check all 3 neighbors
662 for (int e = 0; e < 3; e++)
663 {
664 Face* n = ff[e];
665 if (n && !n->next) // ensure neighbor is not processed yet & isn't on stack
666 {
667 if (n->dot(*_q) <= 0) // if neighbor is not visible we have slihouette edge
668 {
669 // build face
670 add++;
671
672 // ab: given face adjacency [index][],
673 // it provides [][2] vertex indices on shared edge (CCW order)
674 const static int ab[3][2] = { { 1,2 },{ 2,0 },{ 0,1 } };
675
676 Vert* a = fv[ab[e][0]];
677 Vert* b = fv[ab[e][1]];
678
679 Face* s = Face::Alloc(&cache);
680 s->v[0] = a;
681 s->v[1] = b;
682 s->v[2] = _q;
683
684 s->n = s->cross();
685 s->f[2] = n;
686
687 // change neighbour's adjacency from old visible face to cone side
688 if (n->f[0] == _f)
689 n->f[0] = s;
690 else
691 if (n->f[1] == _f)
692 n->f[1] = s;
693 else
694 if (n->f[2] == _f)
695 n->f[2] = s;
696 else
697 assert(0);
698
699 // build silhouette needed for sewing sides in the second pass
700 a->sew = s;
701 a->next = b;
702 }
703 else
704 {
705 // disjoin visible faces
706 // so they won't be processed more than once
707
708 if (n->f[0] == _f)
709 n->f[0] = 0;
710 else
711 if (n->f[1] == _f)
712 n->f[1] = 0;
713 else
714 if (n->f[2] == _f)
715 n->f[2] = 0;
716 else
717 assert(0);
718
719 // push neighbor face, it's visible and requires processing
720 n->next = stack ? stack : n;
721 stack = n;
722 }
723 }
724 }
725 }
726
727 // if add<del+2 hungry hull has consumed some point
728 // that means we can't do delaunay for some under precision reasons
729 // although convex hull would be fine with it
730 assert(add == del + 2);
731
732 // 3. SEW SIDES OF CONE BUILT ON SLIHOUTTE SEGMENTS
733
734 hull = face_alloc + 2 * i - 4 + 1; // last added face
735
736 // last face must contain part of the silhouette
737 // (edge between its v[0] and v[1])
738 Vert* entry = (Vert*)hull->v[0];
739
740 Vert* pr = entry;
741 do
742 {
743 // sew pr<->nx
744 Vert* nx = (Vert*)pr->next;
745 pr->sew->f[0] = nx->sew;
746 nx->sew->f[1] = pr->sew;
747 pr = nx;
748 } while (pr != entry);
749 }
750
751 assert(2 * i - 4 == hull_faces);
752 //ValidateHull(alloc, hull_faces);
753
754 // needed?
755 for (i = 0; i < points; i++)
756 {
757 vert_alloc[i].next = 0;
758 vert_alloc[i].sew = 0;
759 }
760
761 i = 0;
762 Face** prev_dela = &first_dela_face;
763 Face** prev_hull = &first_hull_face;
764 for (int j = 0; j < hull_faces; j++)
765 {
766 Face* _f = face_alloc + j;
767 if (_f->n.z < 0)
768 {
769 *prev_dela = _f;
770 prev_dela = (Face**)&_f->next;
771 i++;
772 }
773 else
774 {
775 *prev_hull = _f;
776 prev_hull = (Face**)&_f->next;
777 if (((Face*)_f->f[0])->n.z < 0)
778 {
779 _f->v[1]->next = _f->v[2];
780 ((Vert*)_f->v[1])->sew = _f;
781 }
782 if (((Face*)_f->f[1])->n.z < 0)
783 {
784 _f->v[2]->next = _f->v[0];
785 ((Vert*)_f->v[2])->sew = _f;
786 }
787 if (((Face*)_f->f[2])->n.z < 0)
788 {
789 _f->v[0]->next = _f->v[1];
790 ((Vert*)_f->v[0])->sew = _f;
791 }
792 }
793 }
794
795 *prev_dela = 0;
796 *prev_hull = 0;
797
798 first_hull_vert = (Vert*)first_hull_face->v[0];
799
800 // todo link slihouette verts into contour
801 // and sew its edges with hull faces
802
803 return 3*i;
804 }
805
ReallocVertsCDelaBella806 bool ReallocVerts(int points)
807 {
808 inp_verts = points;
809 out_verts = 0;
810
811 first_dela_face = 0;
812 first_hull_face = 0;
813 first_hull_vert = 0;
814
815 if (max_verts < points)
816 {
817 if (max_verts)
818 {
819 free(vert_alloc);
820 vert_alloc = 0;
821 max_verts = 0;
822 }
823
824 vert_alloc = (Vert*)malloc(sizeof(Vert)*points);
825 if (vert_alloc)
826 max_verts = points;
827 else
828 {
829 if (errlog_proc)
830 errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
831 return false;
832 }
833 }
834
835 return true;
836 }
837
TriangulateCDelaBella838 virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0)
839 {
840 if (!x)
841 return 0;
842
843 if (!y)
844 y = x + 1;
845
846 if (advance_bytes < static_cast<int>(sizeof(float) * 2))
847 advance_bytes = static_cast<int>(sizeof(float) * 2);
848
849 if (!ReallocVerts(points))
850 return 0;
851
852 for (int i = 0; i < points; i++)
853 {
854 Vert* v = vert_alloc + i;
855 v->i = i;
856 v->x = (Signed14)*(const float*)((const char*)x + i*advance_bytes);
857 v->y = (Signed14)*(const float*)((const char*)y + i*advance_bytes);
858 v->z = s14sqr(v->x) + s14sqr(v->y);
859 }
860
861 out_verts = Triangulate();
862 return out_verts;
863 }
864
TriangulateCDelaBella865 virtual int Triangulate(int points, const double* x, const double* y, int advance_bytes)
866 {
867 if (!x)
868 return 0;
869
870 if (!y)
871 y = x + 1;
872
873 if (advance_bytes < static_cast<int>(sizeof(double) * 2))
874 advance_bytes = static_cast<int>(sizeof(double) * 2);
875
876 if (!ReallocVerts(points))
877 return 0;
878
879 for (int i = 0; i < points; i++)
880 {
881 Vert* v = vert_alloc + i;
882 v->i = i;
883 v->x = (Signed14)*(const double*)((const char*)x + i*advance_bytes);
884 v->y = (Signed14)*(const double*)((const char*)y + i*advance_bytes);
885 v->z = s14sqr (v->x) + s14sqr (v->y);
886 }
887
888 out_verts = Triangulate();
889 return out_verts;
890 }
891
DestroyCDelaBella892 virtual void Destroy()
893 {
894 if (face_alloc)
895 free(face_alloc);
896 if (vert_alloc)
897 free(vert_alloc);
898 delete this;
899 }
900
901 // num of points passed to last call to Triangulate()
GetNumInputPointsCDelaBella902 virtual int GetNumInputPoints() const
903 {
904 return inp_verts;
905 }
906
907 // num of verts returned from last call to Triangulate()
GetNumOutputVertsCDelaBella908 virtual int GetNumOutputVerts() const
909 {
910 return out_verts;
911 }
912
GetFirstDelaunayTriangleCDelaBella913 virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const
914 {
915 return first_dela_face;
916 }
917
GetFirstHullTriangleCDelaBella918 virtual const DelaBella_Triangle* GetFirstHullTriangle() const
919 {
920 return first_hull_face;
921 }
922
GetFirstHullVertexCDelaBella923 virtual const DelaBella_Vertex* GetFirstHullVertex() const
924 {
925 return first_hull_vert;
926 }
927
SetErrLogCDelaBella928 virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream)
929 {
930 errlog_proc = proc;
931 errlog_file = stream;
932 }
933 };
934
Create()935 IDelaBella* IDelaBella::Create()
936 {
937 CDelaBella* db = new CDelaBella;
938 if (!db)
939 return 0;
940
941 db->vert_alloc = 0;
942 db->face_alloc = 0;
943 db->max_verts = 0;
944 db->max_faces = 0;
945
946 db->first_dela_face = 0;
947 db->first_hull_face = 0;
948 db->first_hull_vert = 0;
949
950 db->inp_verts = 0;
951 db->out_verts = 0;
952
953 db->errlog_proc = 0;
954 db->errlog_file = 0;
955
956 return db;
957 }
958
DelaBella_Create()959 void* DelaBella_Create()
960 {
961 return IDelaBella::Create();
962 }
963
DelaBella_Destroy(void * db)964 void DelaBella_Destroy(void* db)
965 {
966 ((IDelaBella*)db)->Destroy();
967 }
968
DelaBella_SetErrLog(void * db,int (* proc)(void * stream,const char * fmt,...),void * stream)969 void DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream)
970 {
971 ((IDelaBella*)db)->SetErrLog(proc, stream);
972 }
973
DelaBella_TriangulateFloat(void * db,int points,float * x,float * y,int advance_bytes)974 int DelaBella_TriangulateFloat(void* db, int points, float* x, float* y, int advance_bytes)
975 {
976 return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
977 }
978
DelaBella_TriangulateDouble(void * db,int points,double * x,double * y,int advance_bytes)979 int DelaBella_TriangulateDouble(void* db, int points, double* x, double* y, int advance_bytes)
980 {
981 return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
982 }
983
DelaBella_GetNumInputPoints(void * db)984 int DelaBella_GetNumInputPoints(void* db)
985 {
986 return ((IDelaBella*)db)->GetNumInputPoints();
987 }
988
DelaBella_GetNumOutputVerts(void * db)989 int DelaBella_GetNumOutputVerts(void* db)
990 {
991 return ((IDelaBella*)db)->GetNumOutputVerts();
992 }
993
GetFirstDelaunayTriangle(void * db)994 const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db)
995 {
996 return ((IDelaBella*)db)->GetFirstDelaunayTriangle();
997 }
998
GetFirstHullTriangle(void * db)999 const DelaBella_Triangle* GetFirstHullTriangle(void* db)
1000 {
1001 return ((IDelaBella*)db)->GetFirstHullTriangle();
1002 }
1003
GetFirstHullVertex(void * db)1004 const DelaBella_Vertex* GetFirstHullVertex(void* db)
1005 {
1006 return ((IDelaBella*)db)->GetFirstHullVertex();
1007 }
1008
1009 // depreciated!
DelaBella(int points,const double * xy,int * abc,int (* errlog)(const char * fmt,...))1010 int DelaBella(int points, const double* xy, int* abc, int(*errlog)(const char* fmt, ...))
1011 {
1012 if (errlog)
1013 errlog("[WRN] Depreciated interface! errlog disabled.\n");
1014
1015 if (!xy || points <= 0)
1016 return 0;
1017
1018 IDelaBella* db = IDelaBella::Create();
1019 int verts = db->Triangulate(points, xy, 0, 0);
1020
1021 if (!abc)
1022 return verts;
1023
1024 if (verts > 0)
1025 {
1026 int tris = verts / 3;
1027 const DelaBella_Triangle* dela = db->GetFirstDelaunayTriangle();
1028 for (int i = 0; i < tris; i++)
1029 {
1030 for (int j=0; j<3; j++)
1031 abc[3 * i + j] = dela->v[j]->i;
1032 dela = dela->next;
1033 }
1034 }
1035 else
1036 {
1037 int pnts = -verts;
1038 const DelaBella_Vertex* line = db->GetFirstHullVertex();
1039 for (int i = 0; i < pnts; i++)
1040 {
1041 abc[i] = line->i;
1042 line = line->next;
1043 }
1044 }
1045
1046 return verts;
1047 }
1048