5 DELABELLA - Delaunay triangulation library
6 Copyright (C) 2018 GUMIX - Marcin Sokalski
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:
15 The above copyright notice and this permission notice shall be included in all
16 copies or substantial portions of the Software.
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
31 // gcc 4.9 for Android doesn't have search.h
32 #if !defined(__ANDROID__) || defined(__clang__)
36 #if (defined(__APPLE__))
37 #include <malloc/malloc.h>
43 #include "delabella.pxx" // we just need LOG() macro
45 // assuming BITS is max(X_BITS,Y_BITS)
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)
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:
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)
68 // alternatively, one could use some BigInt implementation
69 // in order to expand xy range
73 static Unsigned28 s14sqr(const Signed14& s)
75 return (Unsigned28)((Signed29)s*s);
90 Norm cross (const Vect& v) const // cross prod
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;
100 struct CDelaBella : IDelaBella
104 struct Vert : DelaBella_Vertex
109 Vect operator - (const Vert& v) const // diff
112 d.x = (Signed15)x - (Signed15)v.x;
113 d.y = (Signed15)y - (Signed15)v.y;
114 d.z = (Signed29)z - (Signed29)v.z;
118 static bool overlap(const Vert* v1, const Vert* v2)
120 return v1->x == v2->x && v1->y == v2->y;
123 bool operator < (const Vert& v) const
125 return u28cmp(this, &v) < 0;
128 static int u28cmp(const void* a, const void* b)
152 struct Face : DelaBella_Triangle
156 static Face* Alloc(Face** from)
159 *from = (Face*)f->next;
170 Face* Next(const Vert* p) const
172 // return next face in same direction as face vertices are (cw/ccw)
183 Signed62 dot(const Vert& p) const // dot
185 Vect d = p - *(Vert*)v[0];
186 return (Signed62)n.x * d.x + (Signed62)n.y * d.y + (Signed62)n.z * d.z;
189 Norm cross() const // cross of diffs
191 return (*(Vert*)v[1] - *(Vert*)v[0]).cross(*(Vert*)v[2] - *(Vert*)v[0]);
200 Face* first_dela_face;
201 Face* first_hull_face;
202 Vert* first_hull_vert;
207 int(*errlog_proc)(void* file, const char* fmt, ...);
210 virtual ~CDelaBella ()
216 int points = inp_verts;
217 std::sort(vert_alloc, vert_alloc + points);
221 int w = 0, r = 1; // skip initial no-dups block
222 while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + w))
235 while (r < points && Vert::overlap(vert_alloc + r, vert_alloc + r - 1))
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++];
246 errlog_proc(errlog_file, "[WRN] detected %d dups in xy array!\n", points - w);
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;
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;
272 int hull_faces = 2 * points - 4;
274 if (max_faces < hull_faces)
279 face_alloc = (Face*)malloc(sizeof(Face) * hull_faces);
281 max_faces = hull_faces;
285 errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
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;
294 Face* cache = face_alloc;
298 f.v[0] = vert_alloc + 0;
299 f.v[1] = vert_alloc + 1;
300 f.v[2] = vert_alloc + 2;
303 bool colinear = f.n.z == 0;
306 /////////////////////////////////////////////////////////////////////////
307 // UNTIL INPUT IS COPLANAR, GROW IT IN FORM OF A 2D CONTOUR
309 . | | after adding . | ________* L
310 . \ Last points to / Head next point . \ ______/ /
311 . *____ | -----> .H *____ |
312 . |\_ \_____ | . |\_ \_____ |
313 . \ \_ \__* - Tail points to Last . \ \_ \__* T
315 . \__ \_ __/ . \__ \_ __/
316 . \__* - Head points to Tail . \__/
319 Vert* head = (Vert*)f.v[0];
320 Vert* tail = (Vert*)f.v[1];
321 Vert* last = (Vert*)f.v[2];
327 while (i < points && f.dot(vert_alloc[i]) == 0)
329 Vert* v = vert_alloc + i;
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...
335 // why not testing sign of dot prod of 2 normals?
336 // that way we'd fall into precision problems
338 Norm LvH = (*v - *last).cross(*head - *last);
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);
344 Norm TvL = (*v - *tail).cross(*last - *tail);
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);
350 if (lvh && !tvl) // insert new f on top of e(2,0) = (last,head)
361 if (tvl && !lvh) // insert new f on top of e(1,2) = (tail,last)
373 // wtf? dilithium crystals are fucked.
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
396 errlog_proc(errlog_file, "[NFO] all input points are cocircular.\n");
401 errlog_proc(errlog_file, "[NFO] trivial case of 3 points, thank you.\n");
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;
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;
411 head->sew = tail->sew = last->sew = first_hull_face;
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;
422 // reverse silhouette
427 first_hull_vert = last;
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;
438 first_hull_vert = head;
441 first_dela_face->n = first_dela_face->cross();
442 first_hull_face->n = first_hull_face->cross();
447 // retract last point it will be added as a cone's top later
449 head = (Vert*)head->next;
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
458 Vert* q = vert_alloc + i;
463 Vert* n = (Vert*)p->next;
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();
475 Face* prev_side = first_side;
477 Face* first_base = 0;
481 Face* base = Face::Alloc(&cache);
485 base->n = base->cross();
487 Face* side = Face::Alloc(&cache);
491 side->n = side->cross();
496 side->f[1] = prev_side;
497 prev_side->f[0] = side;
499 base->f[0] = prev_base;
501 prev_base->f[1] = base;
512 Face* last_side = Face::Alloc(&cache);
516 last_side->n = last_side->cross();
518 last_side->f[1] = prev_side;
519 prev_side->f[0] = last_side;
521 last_side->f[0] = first_side;
522 first_side->f[1] = last_side;
524 first_base->f[0] = first_side;
525 first_side->f[2] = first_base;
527 last_side->f[2] = prev_base;
528 prev_base->f[1] = last_side;
535 Vert* n = (Vert*)p->next;
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();
547 Face* prev_side = first_side;
549 Face* first_base = 0;
553 Face* base = Face::Alloc(&cache);
557 base->n = base->cross();
559 Face* side = Face::Alloc(&cache);
563 side->n = side->cross();
568 side->f[0] = prev_side;
569 prev_side->f[1] = side;
571 base->f[1] = prev_base;
573 prev_base->f[0] = base;
584 Face* last_side = Face::Alloc(&cache);
588 last_side->n = last_side->cross();
590 last_side->f[0] = prev_side;
591 prev_side->f[1] = last_side;
593 last_side->f[1] = first_side;
594 first_side->f[0] = last_side;
596 first_base->f[1] = first_side;
597 first_side->f[2] = first_base;
599 last_side->f[2] = prev_base;
600 prev_base->f[0] = last_side;
605 /////////////////////////////////////////////////////////////////////////
608 for (; i < points; i++)
610 //ValidateHull(alloc, 2 * i - 4);
612 Vert* _q = vert_alloc + i;
613 Vert* _p = vert_alloc + i - 1;
616 // 1. FIND FIRST VISIBLE FACE
617 // simply iterate around last vertex using last added triangle adjecency info
618 while (_f->dot(*_q) <= 0)
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)
629 assert(_f != face_alloc); // no face is visible? you must be kidding!
635 // 2. DELETE VISIBLE FACES & ADD NEW ONES
636 // (we also build silhouette (vertex loop) between visible & invisible faces)
641 // push first visible face onto stack (of visible faces)
643 _f->next = _f; // old trick to use list pointers as 'on-stack' markers
646 // pop, take care of last item ptr (it's not null!)
648 stack = (Face*)_f->next;
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] };
657 // delete visible face
661 // check all 3 neighbors
662 for (int e = 0; e < 3; e++)
665 if (n && !n->next) // ensure neighbor is not processed yet & isn't on stack
667 if (n->dot(*_q) <= 0) // if neighbor is not visible we have slihouette edge
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 } };
676 Vert* a = fv[ab[e][0]];
677 Vert* b = fv[ab[e][1]];
679 Face* s = Face::Alloc(&cache);
687 // change neighbour's adjacency from old visible face to cone side
699 // build silhouette needed for sewing sides in the second pass
705 // disjoin visible faces
706 // so they won't be processed more than once
719 // push neighbor face, it's visible and requires processing
720 n->next = stack ? stack : n;
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);
732 // 3. SEW SIDES OF CONE BUILT ON SLIHOUTTE SEGMENTS
734 hull = face_alloc + 2 * i - 4 + 1; // last added face
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];
744 Vert* nx = (Vert*)pr->next;
745 pr->sew->f[0] = nx->sew;
746 nx->sew->f[1] = pr->sew;
748 } while (pr != entry);
751 assert(2 * i - 4 == hull_faces);
752 //ValidateHull(alloc, hull_faces);
755 for (i = 0; i < points; i++)
757 vert_alloc[i].next = 0;
758 vert_alloc[i].sew = 0;
762 Face** prev_dela = &first_dela_face;
763 Face** prev_hull = &first_hull_face;
764 for (int j = 0; j < hull_faces; j++)
766 Face* _f = face_alloc + j;
770 prev_dela = (Face**)&_f->next;
776 prev_hull = (Face**)&_f->next;
777 if (((Face*)_f->f[0])->n.z < 0)
779 _f->v[1]->next = _f->v[2];
780 ((Vert*)_f->v[1])->sew = _f;
782 if (((Face*)_f->f[1])->n.z < 0)
784 _f->v[2]->next = _f->v[0];
785 ((Vert*)_f->v[2])->sew = _f;
787 if (((Face*)_f->f[2])->n.z < 0)
789 _f->v[0]->next = _f->v[1];
790 ((Vert*)_f->v[0])->sew = _f;
798 first_hull_vert = (Vert*)first_hull_face->v[0];
800 // todo link slihouette verts into contour
801 // and sew its edges with hull faces
806 bool ReallocVerts(int points)
815 if (max_verts < points)
824 vert_alloc = (Vert*)malloc(sizeof(Vert)*points);
830 errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
838 virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0)
846 if (advance_bytes < static_cast<int>(sizeof(float) * 2))
847 advance_bytes = static_cast<int>(sizeof(float) * 2);
849 if (!ReallocVerts(points))
852 for (int i = 0; i < points; i++)
854 Vert* v = vert_alloc + 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);
861 out_verts = Triangulate();
865 virtual int Triangulate(int points, const double* x, const double* y, int advance_bytes)
873 if (advance_bytes < static_cast<int>(sizeof(double) * 2))
874 advance_bytes = static_cast<int>(sizeof(double) * 2);
876 if (!ReallocVerts(points))
879 for (int i = 0; i < points; i++)
881 Vert* v = vert_alloc + 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);
888 out_verts = Triangulate();
892 virtual void Destroy()
901 // num of points passed to last call to Triangulate()
902 virtual int GetNumInputPoints() const
907 // num of verts returned from last call to Triangulate()
908 virtual int GetNumOutputVerts() const
913 virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const
915 return first_dela_face;
918 virtual const DelaBella_Triangle* GetFirstHullTriangle() const
920 return first_hull_face;
923 virtual const DelaBella_Vertex* GetFirstHullVertex() const
925 return first_hull_vert;
928 virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream)
931 errlog_file = stream;
935 IDelaBella* IDelaBella::Create()
937 CDelaBella* db = new CDelaBella;
946 db->first_dela_face = 0;
947 db->first_hull_face = 0;
948 db->first_hull_vert = 0;
959 void* DelaBella_Create()
961 return IDelaBella::Create();
964 void DelaBella_Destroy(void* db)
966 ((IDelaBella*)db)->Destroy();
969 void DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream)
971 ((IDelaBella*)db)->SetErrLog(proc, stream);
974 int DelaBella_TriangulateFloat(void* db, int points, float* x, float* y, int advance_bytes)
976 return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
979 int DelaBella_TriangulateDouble(void* db, int points, double* x, double* y, int advance_bytes)
981 return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
984 int DelaBella_GetNumInputPoints(void* db)
986 return ((IDelaBella*)db)->GetNumInputPoints();
989 int DelaBella_GetNumOutputVerts(void* db)
991 return ((IDelaBella*)db)->GetNumOutputVerts();
994 const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db)
996 return ((IDelaBella*)db)->GetFirstDelaunayTriangle();
999 const DelaBella_Triangle* GetFirstHullTriangle(void* db)
1001 return ((IDelaBella*)db)->GetFirstHullTriangle();
1004 const DelaBella_Vertex* GetFirstHullVertex(void* db)
1006 return ((IDelaBella*)db)->GetFirstHullVertex();
1010 int DelaBella(int points, const double* xy, int* abc, int(*errlog)(const char* fmt, ...))
1013 errlog("[WRN] Depreciated interface! errlog disabled.\n");
1015 if (!xy || points <= 0)
1018 IDelaBella* db = IDelaBella::Create();
1019 int verts = db->Triangulate(points, xy, 0, 0);
1026 int tris = verts / 3;
1027 const DelaBella_Triangle* dela = db->GetFirstDelaunayTriangle();
1028 for (int i = 0; i < tris; i++)
1030 for (int j=0; j<3; j++)
1031 abc[3 * i + j] = dela->v[j]->i;
1038 const DelaBella_Vertex* line = db->GetFirstHullVertex();
1039 for (int i = 0; i < pnts; i++)