BLI_math_rotation: properly name the quaternion power function.
[blender.git] / source / blender / blenlib / intern / polyfill_2d.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * ***** END GPL LICENSE BLOCK *****
19  */
20
21 /** \file blender/blenlib/intern/polyfill_2d.c
22  *  \ingroup bli
23  *
24  * An ear clipping algorithm to triangulate single boundary polygons.
25  *
26  * Details:
27  *
28  * - The algorithm guarantees all triangles are assigned (number of coords - 2)
29  *   and that triangles will have non-overlapping indices (even for degenerate geometry).
30  * - Self-intersections are considered degenerate (resulting triangles will overlap).
31  * - While multiple polygons aren't supported, holes can still be defined using *key-holes*
32  *   (where the polygon doubles back on its self with *exactly* matching coordinates).
33  *
34  * \note
35  *
36  * Changes made for Blender.
37  *
38  * - loop the array to clip last verts first (less array resizing)
39  *
40  * - advance the ear to clip each iteration
41  *   to avoid fan-filling convex shapes (USE_CLIP_EVEN).
42  *
43  * - avoid intersection tests when there are no convex points (USE_CONVEX_SKIP).
44  *
45  * \note
46  *
47  * No globals - keep threadsafe.
48  */
49
50 #include "BLI_utildefines.h"
51 #include "BLI_math.h"
52
53 #include "BLI_memarena.h"
54 #include "BLI_alloca.h"
55
56 #include "BLI_polyfill_2d.h"  /* own include */
57
58 #include "BLI_strict_flags.h"
59
60 /* avoid fan-fill topology */
61 #define USE_CLIP_EVEN
62 #define USE_CONVEX_SKIP
63 /* sweep back-and-forth about convex ears (avoids lop-sided fans) */
64 #define USE_CLIP_SWEEP
65 // #define USE_CONVEX_SKIP_TEST
66
67 #ifdef USE_CONVEX_SKIP
68 #  define USE_KDTREE
69 #endif
70
71 /* disable in production, it can fail on near zero area ngons */
72 // #define USE_STRICT_ASSERT
73
74 // #define DEBUG_TIME
75 #ifdef DEBUG_TIME
76 #  include "PIL_time_utildefines.h"
77 #endif
78
79
80 typedef signed char eSign;
81
82 #ifdef USE_KDTREE
83 /**
84  * Spatial optimization for point-in-triangle intersection checks.
85  * The simple version of this algorithm is ``O(n^2)`` complexity
86  * (every point needing to check the triangle defined by every other point),
87  * Using a binary-tree reduces the complexity to ``O(n log n)``
88  * plus some overhead of creating the tree.
89  *
90  * This is a single purpose KDTree based on BLI_kdtree with some modifications
91  * to better suit polyfill2d.
92  *
93  *
94  * - #KDTreeNode2D is kept small (only 16 bytes),
95  *   by not storing coords in the nodes and using index values rather then pointers
96  *   to reference neg/pos values.
97  *
98  * - #kdtree2d_isect_tri is the only function currently used.
99  *   This simply intersects a triangle with the kdtree points.
100  *
101  * - the KDTree is only built & used when the polygon is concave.
102  */
103
104 typedef bool axis_t;
105
106 /* use for sorting */
107 typedef struct KDTreeNode2D_head {
108         uint neg, pos;
109         uint index;
110 } KDTreeNode2D_head;
111
112 typedef struct KDTreeNode2D {
113         uint neg, pos;
114         uint index;
115         axis_t axis;  /* range is only (0-1) */
116         ushort flag;
117         uint parent;
118 } KDTreeNode2D;
119
120 struct KDTree2D {
121         KDTreeNode2D *nodes;
122         const float (*coords)[2];
123         uint root;
124         uint totnode;
125         uint *nodes_map;  /* index -> node lookup */
126 };
127
128 struct KDRange2D {
129         float min, max;
130 };
131 #endif  /* USE_KDTREE */
132
133 enum {
134         CONCAVE = -1,
135         TANGENTIAL = 0,
136         CONVEX = 1,
137 };
138
139 typedef struct PolyFill {
140         struct PolyIndex *indices;  /* vertex aligned */
141
142         const float (*coords)[2];
143         uint  coords_tot;
144 #ifdef USE_CONVEX_SKIP
145         uint  coords_tot_concave;
146 #endif
147
148         /* A polygon with n vertices has a triangulation of n-2 triangles. */
149         uint (*tris)[3];
150         uint   tris_tot;
151
152 #ifdef USE_KDTREE
153         struct KDTree2D kdtree;
154 #endif
155 } PolyFill;
156
157
158 /* circular linklist */
159 typedef struct PolyIndex {
160         struct PolyIndex *next, *prev;
161         uint index;
162         eSign sign;
163 } PolyIndex;
164
165
166 /* based on libgdx 2013-11-28, apache 2.0 licensed */
167
168 static void pf_coord_sign_calc(PolyFill *pf, PolyIndex *pi);
169
170 static PolyIndex *pf_ear_tip_find(
171         PolyFill *pf
172 #ifdef USE_CLIP_EVEN
173         , PolyIndex *pi_ear_init
174 #endif
175 #ifdef USE_CLIP_SWEEP
176         , bool reverse
177 #endif
178         );
179
180 static bool       pf_ear_tip_check(PolyFill *pf, PolyIndex *pi_ear_tip);
181 static void       pf_ear_tip_cut(PolyFill *pf, PolyIndex *pi_ear_tip);
182
183
184 BLI_INLINE eSign signum_enum(float a)
185 {
186         if (UNLIKELY(a == 0.0f))
187                 return  0;
188         else if (a > 0.0f)
189                 return  1;
190         else
191                 return -1;
192 }
193
194 /**
195  * alternative version of #area_tri_signed_v2
196  * needed because of float precision issues
197  *
198  * \note removes / 2 since its not needed since we only need the sign.
199  */
200 BLI_INLINE float area_tri_signed_v2_alt_2x(const float v1[2], const float v2[2], const float v3[2])
201 {
202         return ((v1[0] * (v2[1] - v3[1])) +
203                 (v2[0] * (v3[1] - v1[1])) +
204                 (v3[0] * (v1[1] - v2[1])));
205 }
206
207
208 static eSign span_tri_v2_sign(const float v1[2], const float v2[2], const float v3[2])
209 {
210         return signum_enum(area_tri_signed_v2_alt_2x(v3, v2, v1));
211 }
212
213
214 #ifdef USE_KDTREE
215 #define KDNODE_UNSET ((uint)-1)
216
217 enum {
218         KDNODE_FLAG_REMOVED = (1 << 0),
219 };
220
221 static void kdtree2d_new(
222         struct KDTree2D *tree,
223         uint tot,
224         const float (*coords)[2])
225 {
226         /* set by caller */
227         // tree->nodes = nodes;
228         tree->coords = coords;
229         tree->root = KDNODE_UNSET;
230         tree->totnode = tot;
231 }
232
233 /**
234  * no need for kdtree2d_insert, since we know the coords array.
235  */
236 static void kdtree2d_init(
237         struct KDTree2D *tree,
238         const uint coords_tot,
239         const PolyIndex *indices)
240 {
241         KDTreeNode2D *node;
242         uint i;
243
244         for (i = 0, node = tree->nodes; i < coords_tot; i++) {
245                 if (indices[i].sign != CONVEX) {
246                         node->neg = node->pos = KDNODE_UNSET;
247                         node->index = indices[i].index;
248                         node->axis = 0;
249                         node->flag = 0;
250                         node++;
251                 }
252         }
253
254         BLI_assert(tree->totnode == (uint)(node - tree->nodes));
255 }
256
257 static uint kdtree2d_balance_recursive(
258         KDTreeNode2D *nodes, uint totnode, axis_t axis,
259         const float (*coords)[2], const uint ofs)
260 {
261         KDTreeNode2D *node;
262         uint neg, pos, median, i, j;
263
264         if (totnode <= 0) {
265                 return KDNODE_UNSET;
266         }
267         else if (totnode == 1) {
268                 return 0 + ofs;
269         }
270
271         /* quicksort style sorting around median */
272         neg = 0;
273         pos = totnode - 1;
274         median = totnode / 2;
275
276         while (pos > neg) {
277                 const float co = coords[nodes[pos].index][axis];
278                 i = neg - 1;
279                 j = pos;
280
281                 while (1) {
282                         while (coords[nodes[++i].index][axis] < co) ;
283                         while (coords[nodes[--j].index][axis] > co && j > neg) ;
284
285                         if (i >= j) {
286                                 break;
287                         }
288                         SWAP(KDTreeNode2D_head, *(KDTreeNode2D_head *)&nodes[i], *(KDTreeNode2D_head *)&nodes[j]);
289                 }
290
291                 SWAP(KDTreeNode2D_head, *(KDTreeNode2D_head *)&nodes[i], *(KDTreeNode2D_head *)&nodes[pos]);
292                 if (i >= median) {
293                         pos = i - 1;
294                 }
295                 if (i <= median) {
296                         neg = i + 1;
297                 }
298         }
299
300         /* set node and sort subnodes */
301         node = &nodes[median];
302         node->axis = axis;
303         axis = !axis;
304         node->neg = kdtree2d_balance_recursive(nodes, median, axis, coords, ofs);
305         node->pos = kdtree2d_balance_recursive(&nodes[median + 1], (totnode - (median + 1)), axis, coords, (median + 1) + ofs);
306
307         return median + ofs;
308 }
309
310 static void kdtree2d_balance(
311         struct KDTree2D *tree)
312 {
313         tree->root = kdtree2d_balance_recursive(tree->nodes, tree->totnode, 0, tree->coords, 0);
314 }
315
316
317 static void kdtree2d_init_mapping(
318         struct KDTree2D *tree)
319 {
320         uint i;
321         KDTreeNode2D *node;
322
323         for (i = 0, node = tree->nodes; i < tree->totnode; i++, node++) {
324                 if (node->neg != KDNODE_UNSET) {
325                         tree->nodes[node->neg].parent = i;
326                 }
327                 if (node->pos != KDNODE_UNSET) {
328                         tree->nodes[node->pos].parent = i;
329                 }
330
331                 /* build map */
332                 BLI_assert(tree->nodes_map[node->index] == KDNODE_UNSET);
333                 tree->nodes_map[node->index] = i;
334         }
335
336         tree->nodes[tree->root].parent = KDNODE_UNSET;
337 }
338
339 static void kdtree2d_node_remove(
340         struct KDTree2D *tree,
341         uint index)
342 {
343         uint node_index = tree->nodes_map[index];
344         KDTreeNode2D *node;
345
346         if (node_index == KDNODE_UNSET) {
347                 return;
348         }
349         else {
350                 tree->nodes_map[index] = KDNODE_UNSET;
351         }
352
353         node = &tree->nodes[node_index];
354         tree->totnode -= 1;
355
356         BLI_assert((node->flag & KDNODE_FLAG_REMOVED) == 0);
357         node->flag |= KDNODE_FLAG_REMOVED;
358
359         while ((node->neg == KDNODE_UNSET) &&
360                (node->pos == KDNODE_UNSET) &&
361                (node->parent != KDNODE_UNSET))
362         {
363                 KDTreeNode2D *node_parent = &tree->nodes[node->parent];
364
365                 BLI_assert((uint)(node - tree->nodes) == node_index);
366                 if (node_parent->neg == node_index) {
367                         node_parent->neg = KDNODE_UNSET;
368                 }
369                 else {
370                         BLI_assert(node_parent->pos == node_index);
371                         node_parent->pos = KDNODE_UNSET;
372                 }
373
374                 if (node_parent->flag & KDNODE_FLAG_REMOVED) {
375                         node_index = node->parent;
376                         node = node_parent;
377                 }
378                 else {
379                         break;
380                 }
381         }
382 }
383
384 static bool kdtree2d_isect_tri_recursive(
385         const struct KDTree2D *tree,
386         const uint         tri_index[3],
387         const float       *tri_coords[3],
388         const float        tri_center[2],
389         const struct KDRange2D bounds[2],
390         const KDTreeNode2D *node)
391 {
392         const float *co = tree->coords[node->index];
393
394         /* bounds then triangle intersect */
395         if ((node->flag & KDNODE_FLAG_REMOVED) == 0) {
396                 /* bounding box test first */
397                 if ((co[0] >= bounds[0].min) &&
398                     (co[0] <= bounds[0].max) &&
399                     (co[1] >= bounds[1].min) &&
400                     (co[1] <= bounds[1].max))
401                 {
402                         if ((span_tri_v2_sign(tri_coords[0], tri_coords[1], co) != CONCAVE) &&
403                             (span_tri_v2_sign(tri_coords[1], tri_coords[2], co) != CONCAVE) &&
404                             (span_tri_v2_sign(tri_coords[2], tri_coords[0], co) != CONCAVE))
405                         {
406                                 if (!ELEM(node->index, tri_index[0], tri_index[1], tri_index[2])) {
407                                         return true;
408                                 }
409                         }
410                 }
411         }
412
413 #define KDTREE2D_ISECT_TRI_RECURSE_NEG \
414         (((node->neg != KDNODE_UNSET) && (co[node->axis] >= bounds[node->axis].min)) &&  \
415           (kdtree2d_isect_tri_recursive(tree, tri_index, tri_coords, tri_center, bounds, \
416                                         &tree->nodes[node->neg])))
417 #define KDTREE2D_ISECT_TRI_RECURSE_POS \
418         (((node->pos != KDNODE_UNSET) && (co[node->axis] <= bounds[node->axis].max)) &&  \
419           (kdtree2d_isect_tri_recursive(tree, tri_index, tri_coords, tri_center, bounds, \
420                                         &tree->nodes[node->pos])))
421
422         if (tri_center[node->axis] > co[node->axis]) {
423                 if (KDTREE2D_ISECT_TRI_RECURSE_POS) {
424                         return true;
425                 }
426                 if (KDTREE2D_ISECT_TRI_RECURSE_NEG) {
427                         return true;
428                 }
429         }
430         else {
431                 if (KDTREE2D_ISECT_TRI_RECURSE_NEG) {
432                         return true;
433                 }
434                 if (KDTREE2D_ISECT_TRI_RECURSE_POS) {
435                         return true;
436                 }
437         }
438
439 #undef KDTREE2D_ISECT_TRI_RECURSE_NEG
440 #undef KDTREE2D_ISECT_TRI_RECURSE_POS
441
442         BLI_assert(node->index != KDNODE_UNSET);
443
444         return false;
445 }
446
447 static bool kdtree2d_isect_tri(
448         struct KDTree2D *tree,
449         const uint ind[3])
450 {
451         const float *vs[3];
452         uint i;
453         struct KDRange2D bounds[2] = {
454             {FLT_MAX, -FLT_MAX},
455             {FLT_MAX, -FLT_MAX},
456         };
457         float tri_center[2] = {0.0f, 0.0f};
458
459         for (i = 0; i < 3; i++) {
460                 vs[i] = tree->coords[ind[i]];
461
462                 add_v2_v2(tri_center, vs[i]);
463
464                 CLAMP_MAX(bounds[0].min, vs[i][0]);
465                 CLAMP_MIN(bounds[0].max, vs[i][0]);
466                 CLAMP_MAX(bounds[1].min, vs[i][1]);
467                 CLAMP_MIN(bounds[1].max, vs[i][1]);
468         }
469
470         mul_v2_fl(tri_center, 1.0f / 3.0f);
471
472         return kdtree2d_isect_tri_recursive(tree, ind, vs, tri_center, bounds, &tree->nodes[tree->root]);
473 }
474
475 #endif  /* USE_KDTREE */
476
477
478 static uint *pf_tri_add(PolyFill *pf)
479 {
480         return pf->tris[pf->tris_tot++];
481 }
482
483 static void pf_coord_remove(PolyFill *pf, PolyIndex *pi)
484 {
485 #ifdef USE_KDTREE
486         /* avoid double lookups, since convex coords are ignored when testing intersections */
487         if (pf->kdtree.totnode) {
488                 kdtree2d_node_remove(&pf->kdtree, pi->index);
489         }
490 #endif
491
492         pi->next->prev = pi->prev;
493         pi->prev->next = pi->next;
494
495         if (UNLIKELY(pf->indices == pi)) {
496                 pf->indices = pi->next;
497         }
498 #ifdef DEBUG
499         pi->index = (uint)-1;
500         pi->next = pi->prev = NULL;
501 #endif
502
503         pf->coords_tot -= 1;
504 }
505
506 static void pf_triangulate(PolyFill *pf)
507 {
508         /* localize */
509         PolyIndex *pi_ear;
510
511 #ifdef USE_CLIP_EVEN
512         PolyIndex *pi_ear_init = pf->indices;
513 #endif
514 #ifdef USE_CLIP_SWEEP
515         bool reverse = false;
516 #endif
517
518         while (pf->coords_tot > 3) {
519                 PolyIndex *pi_prev, *pi_next;
520                 eSign sign_orig_prev, sign_orig_next;
521
522                 pi_ear = pf_ear_tip_find(
523                         pf
524 #ifdef USE_CLIP_EVEN
525                         , pi_ear_init
526 #endif
527 #ifdef USE_CLIP_SWEEP
528                         , reverse
529 #endif
530                         );
531
532 #ifdef USE_CONVEX_SKIP
533                 if (pi_ear->sign != CONVEX) {
534                         pf->coords_tot_concave -= 1;
535                 }
536 #endif
537
538                 pi_prev = pi_ear->prev;
539                 pi_next = pi_ear->next;
540
541                 pf_ear_tip_cut(pf, pi_ear);
542
543                 /* The type of the two vertices adjacent to the clipped vertex may have changed. */
544                 sign_orig_prev = pi_prev->sign;
545                 sign_orig_next = pi_next->sign;
546
547                 /* check if any verts became convex the (else if)
548                  * case is highly unlikely but may happen with degenerate polygons */
549                 if (sign_orig_prev != CONVEX) {
550                         pf_coord_sign_calc(pf, pi_prev);
551 #ifdef USE_CONVEX_SKIP
552                         if (pi_prev->sign == CONVEX) {
553                                 pf->coords_tot_concave -= 1;
554 #ifdef USE_KDTREE
555                                 kdtree2d_node_remove(&pf->kdtree, pi_prev->index);
556 #endif
557                         }
558 #endif
559                 }
560                 if (sign_orig_next != CONVEX) {
561                         pf_coord_sign_calc(pf, pi_next);
562 #ifdef USE_CONVEX_SKIP
563                         if (pi_next->sign == CONVEX) {
564                                 pf->coords_tot_concave -= 1;
565 #ifdef USE_KDTREE
566                                 kdtree2d_node_remove(&pf->kdtree, pi_next->index);
567 #endif
568                         }
569 #endif
570                 }
571
572 #ifdef USE_CLIP_EVEN
573 #ifdef USE_CLIP_SWEEP
574                 pi_ear_init = reverse ? pi_prev->prev : pi_next->next;
575 #else
576                 pi_ear_init = pi_next->next;
577 #endif
578 #endif
579
580 #ifdef USE_CLIP_EVEN
581 #ifdef USE_CLIP_SWEEP
582                 if (pi_ear_init->sign != CONVEX) {
583                         /* take the extra step since this ear isn't a good candidate */
584                         pi_ear_init = reverse ? pi_ear_init->prev : pi_ear_init->next;
585                         reverse = !reverse;
586                 }
587 #endif
588 #else
589                 if ((reverse ? pi_prev->prev : pi_next->next)->sign != CONVEX) {
590                         reverse = !reverse;
591                 }
592 #endif
593
594         }
595
596         if (pf->coords_tot == 3) {
597                 uint *tri = pf_tri_add(pf);
598                 pi_ear = pf->indices;
599                 tri[0] = pi_ear->index; pi_ear = pi_ear->next;
600                 tri[1] = pi_ear->index; pi_ear = pi_ear->next;
601                 tri[2] = pi_ear->index;
602         }
603 }
604
605 /**
606  * \return CONCAVE, TANGENTIAL or CONVEX
607  */
608 static void pf_coord_sign_calc(PolyFill *pf, PolyIndex *pi)
609 {
610         /* localize */
611         const float (*coords)[2] = pf->coords;
612
613         pi->sign = span_tri_v2_sign(
614                 coords[pi->prev->index],
615                 coords[pi->index],
616                 coords[pi->next->index]);
617 }
618
619 static PolyIndex *pf_ear_tip_find(
620         PolyFill *pf
621 #ifdef USE_CLIP_EVEN
622         , PolyIndex *pi_ear_init
623 #endif
624 #ifdef USE_CLIP_SWEEP
625         , bool reverse
626 #endif
627         )
628 {
629         /* localize */
630         const uint coords_tot = pf->coords_tot;
631         PolyIndex *pi_ear;
632
633         uint i;
634
635 #ifdef USE_CLIP_EVEN
636         pi_ear = pi_ear_init;
637 #else
638         pi_ear = pf->indices;
639 #endif
640
641         i = coords_tot;
642         while (i--) {
643                 if (pf_ear_tip_check(pf, pi_ear)) {
644                         return pi_ear;
645                 }
646 #ifdef USE_CLIP_SWEEP
647                 pi_ear = reverse ? pi_ear->prev : pi_ear->next;
648 #else
649                 pi_ear = pi_ear->next;
650 #endif
651         }
652
653         /* Desperate mode: if no vertex is an ear tip, we are dealing with a degenerate polygon (e.g. nearly collinear).
654          * Note that the input was not necessarily degenerate, but we could have made it so by clipping some valid ears.
655          *
656          * Idea taken from Martin Held, "FIST: Fast industrial-strength triangulation of polygons", Algorithmica (1998),
657          * http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.115.291
658          *
659          * Return a convex or tangential vertex if one exists.
660          */
661
662 #ifdef USE_CLIP_EVEN
663         pi_ear = pi_ear_init;
664 #else
665         pi_ear = pf->indices;
666 #endif
667
668         i = coords_tot;
669         while (i--) {
670                 if (pi_ear->sign != CONCAVE) {
671                         return pi_ear;
672                 }
673                 pi_ear = pi_ear->next;
674         }
675
676         /* If all vertices are concave, just return the last one. */
677         return pi_ear;
678 }
679
680 static bool pf_ear_tip_check(PolyFill *pf, PolyIndex *pi_ear_tip)
681 {
682 #ifndef USE_KDTREE
683         /* localize */
684         const float (*coords)[2] = pf->coords;
685         PolyIndex *pi_curr;
686
687         const float *v1, *v2, *v3;
688 #endif
689
690 #if defined(USE_CONVEX_SKIP) && !defined(USE_KDTREE)
691         uint coords_tot_concave_checked = 0;
692 #endif
693
694
695 #ifdef USE_CONVEX_SKIP
696
697 #ifdef USE_CONVEX_SKIP_TEST
698         /* check if counting is wrong */
699         {
700                 uint coords_tot_concave_test = 0;
701                 PolyIndex *pi_iter = pi_ear_tip;
702                 do {
703                         if (pi_iter->sign != CONVEX) {
704                                 coords_tot_concave_test += 1;
705                         }
706                 } while ((pi_iter = pi_iter->next) != pi_ear_tip);
707                 BLI_assert(coords_tot_concave_test == pf->coords_tot_concave);
708         }
709 #endif
710
711         /* fast-path for circles */
712         if (pf->coords_tot_concave == 0) {
713                 return true;
714         }
715 #endif
716
717         if (UNLIKELY(pi_ear_tip->sign == CONCAVE)) {
718                 return false;
719         }
720
721 #ifdef USE_KDTREE
722         {
723                 const uint ind[3] = {
724                     pi_ear_tip->index,
725                     pi_ear_tip->next->index,
726                     pi_ear_tip->prev->index};
727
728                 if (kdtree2d_isect_tri(&pf->kdtree, ind)) {
729                         return false;
730                 }
731         }
732 #else
733
734         v1 = coords[pi_ear_tip->prev->index];
735         v2 = coords[pi_ear_tip->index];
736         v3 = coords[pi_ear_tip->next->index];
737
738         /* Check if any point is inside the triangle formed by previous, current and next vertices.
739          * Only consider vertices that are not part of this triangle, or else we'll always find one inside. */
740
741         for (pi_curr = pi_ear_tip->next->next; pi_curr != pi_ear_tip->prev; pi_curr = pi_curr->next) {
742                 /* Concave vertices can obviously be inside the candidate ear, but so can tangential vertices
743                  * if they coincide with one of the triangle's vertices. */
744                 if (pi_curr->sign != CONVEX) {
745                         const float *v = coords[pi_curr->index];
746                         /* Because the polygon has clockwise winding order,
747                          * the area sign will be positive if the point is strictly inside.
748                          * It will be 0 on the edge, which we want to include as well. */
749
750                         /* note: check (v3, v1) first since it fails _far_ more often then the other 2 checks (those fail equally).
751                          * It's logical - the chance is low that points exist on the same side as the ear we're clipping off. */
752                         if ((span_tri_v2_sign(v3, v1, v) != CONCAVE) &&
753                             (span_tri_v2_sign(v1, v2, v) != CONCAVE) &&
754                             (span_tri_v2_sign(v2, v3, v) != CONCAVE))
755                         {
756                                 return false;
757                         }
758
759 #ifdef USE_CONVEX_SKIP
760                         coords_tot_concave_checked += 1;
761                         if (coords_tot_concave_checked == pf->coords_tot_concave) {
762                                 break;
763                         }
764 #endif
765                 }
766         }
767 #endif  /* USE_KDTREE */
768
769         return true;
770 }
771
772 static void pf_ear_tip_cut(PolyFill *pf, PolyIndex *pi_ear_tip)
773 {
774         uint *tri = pf_tri_add(pf);
775
776         tri[0] = pi_ear_tip->prev->index;
777         tri[1] = pi_ear_tip->index;
778         tri[2] = pi_ear_tip->next->index;
779
780         pf_coord_remove(pf, pi_ear_tip);
781 }
782
783 /**
784  * Initializes the #PolyFill structure before tessellating with #polyfill_calc.
785  */
786 static void polyfill_prepare(
787         PolyFill *pf,
788         const float (*coords)[2],
789         const uint coords_tot,
790         int coords_sign,
791         uint (*r_tris)[3],
792         PolyIndex *r_indices)
793 {
794         /* localize */
795         PolyIndex *indices = r_indices;
796
797         uint i;
798
799         /* assign all polyfill members here */
800         pf->indices = r_indices;
801         pf->coords = coords;
802         pf->coords_tot = coords_tot;
803 #ifdef USE_CONVEX_SKIP
804         pf->coords_tot_concave = 0;
805 #endif
806         pf->tris = r_tris;
807         pf->tris_tot = 0;
808
809         if (coords_sign == 0) {
810                 coords_sign = (cross_poly_v2(coords, coords_tot) >= 0.0f) ? 1 : -1;
811         }
812         else {
813                 /* check we're passing in correcty args */
814 #ifdef USE_STRICT_ASSERT
815 #ifndef NDEBUG
816                 if (coords_sign == 1) {
817                         BLI_assert(cross_poly_v2(coords, coords_tot) >= 0.0f);
818                 }
819                 else {
820                         BLI_assert(cross_poly_v2(coords, coords_tot) <= 0.0f);
821                 }
822 #endif
823 #endif
824         }
825
826         if (coords_sign == 1) {
827                 for (i = 0; i < coords_tot; i++) {
828                         indices[i].next = &indices[i + 1];
829                         indices[i].prev = &indices[i - 1];
830                         indices[i].index = i;
831                 }
832         }
833         else {
834                 /* reversed */
835                 uint n = coords_tot - 1;
836                 for (i = 0; i < coords_tot; i++) {
837                         indices[i].next = &indices[i + 1];
838                         indices[i].prev = &indices[i - 1];
839                         indices[i].index = (n - i);
840                 }
841         }
842         indices[0].prev = &indices[coords_tot - 1];
843         indices[coords_tot - 1].next = &indices[0];
844
845         for (i = 0; i < coords_tot; i++) {
846                 PolyIndex *pi = &indices[i];
847                 pf_coord_sign_calc(pf, pi);
848 #ifdef USE_CONVEX_SKIP
849                 if (pi->sign != CONVEX) {
850                         pf->coords_tot_concave += 1;
851                 }
852 #endif
853         }
854 }
855
856 static void polyfill_calc(
857         PolyFill *pf)
858 {
859 #ifdef USE_KDTREE
860 #ifdef USE_CONVEX_SKIP
861         if (pf->coords_tot_concave)
862 #endif
863         {
864                 kdtree2d_new(&pf->kdtree, pf->coords_tot_concave, pf->coords);
865                 kdtree2d_init(&pf->kdtree, pf->coords_tot, pf->indices);
866                 kdtree2d_balance(&pf->kdtree);
867                 kdtree2d_init_mapping(&pf->kdtree);
868         }
869 #endif
870
871         pf_triangulate(pf);
872 }
873
874 /**
875  * A version of #BLI_polyfill_calc that uses a memory arena to avoid re-allocations.
876  */
877 void BLI_polyfill_calc_arena(
878         const float (*coords)[2],
879         const uint coords_tot,
880         const int coords_sign,
881         uint (*r_tris)[3],
882
883         struct MemArena *arena)
884 {
885         PolyFill pf;
886         PolyIndex *indices = BLI_memarena_alloc(arena, sizeof(*indices) * coords_tot);
887
888 #ifdef DEBUG_TIME
889         TIMEIT_START(polyfill2d);
890 #endif
891
892         polyfill_prepare(
893                 &pf,
894                 coords, coords_tot, coords_sign,
895                 r_tris,
896                 /* cache */
897                 indices);
898
899 #ifdef USE_KDTREE
900         if (pf.coords_tot_concave) {
901                 pf.kdtree.nodes = BLI_memarena_alloc(arena, sizeof(*pf.kdtree.nodes) * pf.coords_tot_concave);
902                 pf.kdtree.nodes_map = memset(BLI_memarena_alloc(arena, sizeof(*pf.kdtree.nodes_map) * coords_tot),
903                                              0xff, sizeof(*pf.kdtree.nodes_map) * coords_tot);
904         }
905         else {
906                 pf.kdtree.totnode = 0;
907         }
908 #endif
909
910         polyfill_calc(&pf);
911
912         /* indices are no longer needed,
913          * caller can clear arena */
914
915 #ifdef DEBUG_TIME
916         TIMEIT_END(polyfill2d);
917 #endif
918 }
919
920 /**
921  * Triangulates the given (convex or concave) simple polygon to a list of triangle vertices.
922  *
923  * \param coords: 2D coordinates describing vertices of the polygon,
924  * in either clockwise or counterclockwise order.
925  * \param coords_tot: Total points in the array.
926  * \param coords_sign: Pass this when we know the sign in advance to avoid extra calculations.
927  *
928  * \param r_tris: This array is filled in with triangle indices in clockwise order.
929  * The length of the array must be ``coords_tot - 2``.
930  * Indices are guaranteed to be assigned to unique triangles, with valid indices,
931  * even in the case of degenerate input (self intersecting polygons, zero area ears... etc).
932  */
933 void BLI_polyfill_calc(
934         const float (*coords)[2],
935         const uint coords_tot,
936         const int coords_sign,
937         uint (*r_tris)[3])
938 {
939         PolyFill pf;
940         PolyIndex *indices = BLI_array_alloca(indices, coords_tot);
941
942 #ifdef DEBUG_TIME
943         TIMEIT_START(polyfill2d);
944 #endif
945
946         polyfill_prepare(
947                 &pf,
948                 coords, coords_tot, coords_sign,
949                 r_tris,
950                 /* cache */
951                 indices);
952
953 #ifdef USE_KDTREE
954         if (pf.coords_tot_concave) {
955                 pf.kdtree.nodes = BLI_array_alloca(pf.kdtree.nodes, pf.coords_tot_concave);
956                 pf.kdtree.nodes_map = memset(BLI_array_alloca(pf.kdtree.nodes_map, coords_tot),
957                                              0xff, sizeof(*pf.kdtree.nodes_map) * coords_tot);
958         }
959         else {
960                 pf.kdtree.totnode = 0;
961         }
962 #endif
963
964         polyfill_calc(&pf);
965
966 #ifdef DEBUG_TIME
967         TIMEIT_END(polyfill2d);
968 #endif
969 }