4041fc2c755a734f9651b21d3e343de634731905
[blender.git] / source / blender / bmesh / intern / bmesh_polygon.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  * Contributor(s): Joseph Eagar, Geoffrey Bantle, Campbell Barton
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/intern/bmesh_polygon.c
24  *  \ingroup bmesh
25  *
26  * This file contains code for dealing
27  * with polygons (normal/area calculation,
28  * tessellation, etc)
29  *
30  * BMESH_TODO:
31  *  - Add in Tessellator frontend that creates
32  *    BMTriangles from copied faces
33  *
34  *  - Add in Function that checks for and flags
35  *    degenerate faces.
36  */
37
38 #include "DNA_listBase.h"
39
40 #include "MEM_guardedalloc.h"
41
42 #include "BLI_math.h"
43 #include "BLI_array.h"
44 #include "BLI_scanfill.h"
45 #include "BLI_listbase.h"
46
47 #include "bmesh.h"
48
49 #include "intern/bmesh_private.h"
50
51 /**
52  * \brief TEST EDGE SIDE and POINT IN TRIANGLE
53  *
54  * Point in triangle tests stolen from scanfill.c.
55  * Used for tessellator
56  */
57
58 static bool testedgesidef(const float v1[2], const float v2[2], const float v3[2])
59 {
60         /* is v3 to the right of v1 - v2 ? With exception: v3 == v1 || v3 == v2 */
61         double inp;
62
63         //inp = (v2[cox] - v1[cox]) * (v1[coy] - v3[coy]) + (v1[coy] - v2[coy]) * (v1[cox] - v3[cox]);
64         inp = (v2[0] - v1[0]) * (v1[1] - v3[1]) + (v1[1] - v2[1]) * (v1[0] - v3[0]);
65
66         if (inp < 0.0) {
67                 return false;
68         }
69         else if (inp == 0) {
70                 if (v1[0] == v3[0] && v1[1] == v3[1]) return false;
71                 if (v2[0] == v3[0] && v2[1] == v3[1]) return false;
72         }
73         return true;
74 }
75
76 /**
77  * \brief COMPUTE POLY NORMAL
78  *
79  * Computes the normal of a planar
80  * polygon See Graphics Gems for
81  * computing newell normal.
82  */
83 static void calc_poly_normal(float normal[3], float verts[][3], int nverts)
84 {
85         float const *v_prev = verts[nverts - 1];
86         float const *v_curr = verts[0];
87         float n[3] = {0.0f};
88         int i;
89
90         /* Newell's Method */
91         for (i = 0; i < nverts; v_prev = v_curr, v_curr = verts[++i]) {
92                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
93         }
94
95         if (UNLIKELY(normalize_v3_v3(normal, n) == 0.0f)) {
96                 normal[2] = 1.0f; /* other axis set to 0.0 */
97         }
98 }
99
100 /**
101  * \brief COMPUTE POLY NORMAL (BMFace)
102  *
103  * Same as #calc_poly_normal but operates directly on a bmesh face.
104  */
105 static void bm_face_calc_poly_normal(BMFace *f)
106 {
107         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
108         BMLoop *l_iter  = l_first;
109         float const *v_prev = l_first->prev->v->co;
110         float const *v_curr = l_first->v->co;
111         float n[3] = {0.0f};
112
113         /* Newell's Method */
114         do {
115                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
116
117                 l_iter = l_iter->next;
118                 v_prev = v_curr;
119                 v_curr = l_iter->v->co;
120
121         } while (l_iter != l_first);
122
123         if (UNLIKELY(normalize_v3_v3(f->no, n) == 0.0f)) {
124                 f->no[2] = 1.0f; /* other axis set to 0.0 */
125         }
126 }
127
128 /**
129  * \brief COMPUTE POLY NORMAL (BMFace)
130  *
131  * Same as #calc_poly_normal and #bm_face_calc_poly_normal
132  * but takes an array of vertex locations.
133  */
134 static void bm_face_calc_poly_normal_vertex_cos(BMFace *f, float n[3],
135                                                 float const (*vertexCos)[3])
136 {
137         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
138         BMLoop *l_iter  = l_first;
139         float const *v_prev = vertexCos[BM_elem_index_get(l_first->prev->v)];
140         float const *v_curr = vertexCos[BM_elem_index_get(l_first->v)];
141
142         zero_v3(n);
143
144         /* Newell's Method */
145         do {
146                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
147
148                 l_iter = l_iter->next;
149                 v_prev = v_curr;
150                 v_curr = vertexCos[BM_elem_index_get(l_iter->v)];
151         } while (l_iter != l_first);
152
153         if (UNLIKELY(normalize_v3(n) == 0.0f)) {
154                 n[2] = 1.0f; /* other axis set to 0.0 */
155         }
156 }
157
158 /**
159  * For tools that insist on using triangles, ideally we would cache this data.
160  *
161  * \param r_loops  Store face loop pointers, (f->len)
162  * \param r_index  Store triangle triples, indicies into \a r_loops,  ((f->len - 2) * 3)
163  */
164 int BM_face_calc_tessellation(BMFace *f, BMLoop **r_loops, int (*_r_index)[3])
165 {
166         int *r_index = (int *)_r_index;
167         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
168         BMLoop *l_iter;
169         int totfilltri;
170
171         if (f->len == 3) {
172                 *r_loops++ = (l_iter = l_first);
173                 *r_loops++ = (l_iter = l_iter->next);
174                 *r_loops++ = (         l_iter->next);
175
176                 r_index[0] = 0;
177                 r_index[1] = 1;
178                 r_index[2] = 2;
179                 totfilltri = 1;
180         }
181         else if (f->len == 4) {
182                 *r_loops++ = (l_iter = l_first);
183                 *r_loops++ = (l_iter = l_iter->next);
184                 *r_loops++ = (l_iter = l_iter->next);
185                 *r_loops++ = (         l_iter->next);
186
187                 r_index[0] = 0;
188                 r_index[1] = 1;
189                 r_index[2] = 2;
190
191                 r_index[3] = 0;
192                 r_index[4] = 2;
193                 r_index[5] = 3;
194                 totfilltri = 2;
195         }
196         else {
197                 int j;
198
199                 ScanFillContext sf_ctx;
200                 ScanFillVert *sf_vert, *sf_vert_last = NULL, *sf_vert_first = NULL;
201                 /* ScanFillEdge *e; */ /* UNUSED */
202                 ScanFillFace *sf_tri;
203
204                 BLI_scanfill_begin(&sf_ctx);
205
206                 j = 0;
207                 l_iter = l_first;
208                 do {
209                         sf_vert = BLI_scanfill_vert_add(&sf_ctx, l_iter->v->co);
210                         sf_vert->tmp.p = l_iter;
211
212                         if (sf_vert_last) {
213                                 /* e = */ BLI_scanfill_edge_add(&sf_ctx, sf_vert_last, sf_vert);
214                         }
215
216                         sf_vert_last = sf_vert;
217                         if (sf_vert_first == NULL) {
218                                 sf_vert_first = sf_vert;
219                         }
220
221                         r_loops[j] = l_iter;
222
223                         /* mark order */
224                         BM_elem_index_set(l_iter, j++); /* set_loop */
225
226                 } while ((l_iter = l_iter->next) != l_first);
227
228                 /* complete the loop */
229                 BLI_scanfill_edge_add(&sf_ctx, sf_vert_first, sf_vert);
230
231                 totfilltri = BLI_scanfill_calc_ex(&sf_ctx, 0, f->no);
232                 BLI_assert(totfilltri <= f->len - 2);
233                 BLI_assert(totfilltri == BLI_countlist(&sf_ctx.fillfacebase));
234
235                 for (sf_tri = sf_ctx.fillfacebase.first; sf_tri; sf_tri = sf_tri->next) {
236                         int i1 = BM_elem_index_get((BMLoop *)sf_tri->v1->tmp.p);
237                         int i2 = BM_elem_index_get((BMLoop *)sf_tri->v2->tmp.p);
238                         int i3 = BM_elem_index_get((BMLoop *)sf_tri->v3->tmp.p);
239
240                         if (i1 > i2) { SWAP(int, i1, i2); }
241                         if (i2 > i3) { SWAP(int, i2, i3); }
242                         if (i1 > i2) { SWAP(int, i1, i2); }
243
244                         *r_index++ = i1;
245                         *r_index++ = i2;
246                         *r_index++ = i3;
247                 }
248
249                 BLI_scanfill_end(&sf_ctx);
250         }
251
252         return totfilltri;
253 }
254
255 /**
256  * get the area of the face
257  */
258 float BM_face_calc_area(BMFace *f)
259 {
260         BMLoop *l;
261         BMIter iter;
262         float (*verts)[3] = BLI_array_alloca(verts, f->len);
263         float area;
264         int i;
265
266         BM_ITER_ELEM_INDEX (l, &iter, f, BM_LOOPS_OF_FACE, i) {
267                 copy_v3_v3(verts[i], l->v->co);
268         }
269
270         if (f->len == 3) {
271                 area = area_tri_v3(verts[0], verts[1], verts[2]);
272         }
273         else if (f->len == 4) {
274                 area = area_quad_v3(verts[0], verts[1], verts[2], verts[3]);
275         }
276         else {
277                 float normal[3];
278                 calc_poly_normal(normal, verts, f->len);
279                 area = area_poly_v3(f->len, verts, normal);
280         }
281
282         return area;
283 }
284
285 /**
286  * compute the perimeter of an ngon
287  */
288 float BM_face_calc_perimeter(BMFace *f)
289 {
290         BMLoop *l_iter, *l_first;
291         float perimeter = 0.0f;
292
293         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
294         do {
295                 perimeter += len_v3v3(l_iter->v->co, l_iter->next->v->co);
296         } while ((l_iter = l_iter->next) != l_first);
297
298         return perimeter;
299 }
300
301 /**
302  * computes center of face in 3d.  uses center of bounding box.
303  */
304 void BM_face_calc_center_bounds(BMFace *f, float r_cent[3])
305 {
306         BMLoop *l_iter;
307         BMLoop *l_first;
308         float min[3], max[3];
309
310         INIT_MINMAX(min, max);
311
312         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
313         do {
314                 minmax_v3v3_v3(min, max, l_iter->v->co);
315         } while ((l_iter = l_iter->next) != l_first);
316
317         mid_v3_v3v3(r_cent, min, max);
318 }
319
320 /**
321  * computes the center of a face, using the mean average
322  */
323 void BM_face_calc_center_mean(BMFace *f, float r_cent[3])
324 {
325         BMLoop *l_iter;
326         BMLoop *l_first;
327
328         zero_v3(r_cent);
329
330         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
331         do {
332                 add_v3_v3(r_cent, l_iter->v->co);
333         } while ((l_iter = l_iter->next) != l_first);
334
335         if (f->len)
336                 mul_v3_fl(r_cent, 1.0f / (float) f->len);
337 }
338
339 /**
340  * COMPUTE POLY PLANE
341  *
342  * Projects a set polygon's vertices to
343  * a plane defined by the average
344  * of its edges cross products
345  */
346 void calc_poly_plane(float (*verts)[3], const int nverts)
347 {
348         
349         float avgc[3], norm[3], mag, avgn[3];
350         float *v1, *v2, *v3;
351         int i;
352         
353         if (nverts < 3)
354                 return;
355
356         zero_v3(avgn);
357         zero_v3(avgc);
358
359         for (i = 0; i < nverts; i++) {
360                 v1 = verts[i];
361                 v2 = verts[(i + 1) % nverts];
362                 v3 = verts[(i + 2) % nverts];
363                 normal_tri_v3(norm, v1, v2, v3);
364
365                 add_v3_v3(avgn, norm);
366         }
367
368         if (UNLIKELY(normalize_v3(avgn) == 0.0f)) {
369                 avgn[2] = 1.0f;
370         }
371         
372         for (i = 0; i < nverts; i++) {
373                 v1 = verts[i];
374                 mag = dot_v3v3(v1, avgn);
375                 madd_v3_v3fl(v1, avgn, -mag);
376         }
377 }
378
379 /**
380  * \brief BM LEGAL EDGES
381  *
382  * takes in a face and a list of edges, and sets to NULL any edge in
383  * the list that bridges a concave region of the face or intersects
384  * any of the faces's edges.
385  */
386 static void scale_edge_v3f(float v1[3], float v2[3], const float fac)
387 {
388         float mid[3];
389
390         mid_v3_v3v3(mid, v1, v2);
391
392         sub_v3_v3v3(v1, v1, mid);
393         sub_v3_v3v3(v2, v2, mid);
394
395         mul_v3_fl(v1, fac);
396         mul_v3_fl(v2, fac);
397
398         add_v3_v3v3(v1, v1, mid);
399         add_v3_v3v3(v2, v2, mid);
400 }
401
402 /**
403  * \brief POLY ROTATE PLANE
404  *
405  * Rotates a polygon so that it's
406  * normal is pointing towards the mesh Z axis
407  */
408 void poly_rotate_plane(const float normal[3], float (*verts)[3], const int nverts)
409 {
410         float mat[3][3];
411
412         if (axis_dominant_v3_to_m3(mat, normal)) {
413                 int i;
414                 for (i = 0; i < nverts; i++) {
415                         mul_m3_v3(mat, verts[i]);
416                 }
417         }
418 }
419
420 /**
421  * updates face and vertex normals incident on an edge
422  */
423 void BM_edge_normals_update(BMEdge *e)
424 {
425         BMIter iter;
426         BMFace *f;
427         
428         BM_ITER_ELEM (f, &iter, e, BM_FACES_OF_EDGE) {
429                 BM_face_normal_update(f);
430         }
431
432         BM_vert_normal_update(e->v1);
433         BM_vert_normal_update(e->v2);
434 }
435
436 /**
437  * update a vert normal (but not the faces incident on it)
438  */
439 void BM_vert_normal_update(BMVert *v)
440 {
441         /* TODO, we can normalize each edge only once, then compare with previous edge */
442
443         BMIter liter;
444         BMLoop *l;
445         float vec1[3], vec2[3], fac;
446         int len = 0;
447
448         zero_v3(v->no);
449
450         BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
451                 /* Same calculation used in BM_mesh_normals_update */
452                 sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
453                 sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
454                 normalize_v3(vec1);
455                 normalize_v3(vec2);
456
457                 fac = saacos(-dot_v3v3(vec1, vec2));
458
459                 madd_v3_v3fl(v->no, l->f->no, fac);
460
461                 len++;
462         }
463
464         if (len) {
465                 normalize_v3(v->no);
466         }
467 }
468
469 void BM_vert_normal_update_all(BMVert *v)
470 {
471         BMIter iter;
472         BMFace *f;
473
474         BM_ITER_ELEM (f, &iter, v, BM_FACES_OF_VERT) {
475                 BM_face_normal_update(f);
476         }
477
478         BM_vert_normal_update(v);
479 }
480
481 /**
482  * \brief BMESH UPDATE FACE NORMAL
483  *
484  * Updates the stored normal for the
485  * given face. Requires that a buffer
486  * of sufficient length to store projected
487  * coordinates for all of the face's vertices
488  * is passed in as well.
489  */
490
491 void BM_face_normal_update(BMFace *f)
492 {
493         BMLoop *l;
494
495         /* common cases first */
496         switch (f->len) {
497                 case 4:
498                 {
499                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
500                         const float *co2 = (l = l->next)->v->co;
501                         const float *co3 = (l = l->next)->v->co;
502                         const float *co4 = (l->next)->v->co;
503
504                         normal_quad_v3(f->no, co1, co2, co3, co4);
505                         break;
506                 }
507                 case 3:
508                 {
509                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
510                         const float *co2 = (l = l->next)->v->co;
511                         const float *co3 = (l->next)->v->co;
512
513                         normal_tri_v3(f->no, co1, co2, co3);
514                         break;
515                 }
516                 case 0:
517                 {
518                         zero_v3(f->no);
519                         break;
520                 }
521                 default:
522                 {
523                         bm_face_calc_poly_normal(f);
524                         break;
525                 }
526         }
527 }
528 /* exact same as 'bmesh_face_normal_update' but accepts vertex coords */
529 void BM_face_normal_update_vcos(BMesh *bm, BMFace *f, float no[3],
530                                 float const (*vertexCos)[3])
531 {
532         BMLoop *l;
533
534         /* must have valid index data */
535         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
536         (void)bm;
537
538         /* common cases first */
539         switch (f->len) {
540                 case 4:
541                 {
542                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
543                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
544                         const float *co3 = vertexCos[BM_elem_index_get((l = l->next)->v)];
545                         const float *co4 = vertexCos[BM_elem_index_get((l->next)->v)];
546
547                         normal_quad_v3(no, co1, co2, co3, co4);
548                         break;
549                 }
550                 case 3:
551                 {
552                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
553                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
554                         const float *co3 = vertexCos[BM_elem_index_get((l->next)->v)];
555
556                         normal_tri_v3(no, co1, co2, co3);
557                         break;
558                 }
559                 case 0:
560                 {
561                         zero_v3(no);
562                         break;
563                 }
564                 default:
565                 {
566                         bm_face_calc_poly_normal_vertex_cos(f, no, vertexCos);
567                         break;
568                 }
569         }
570 }
571
572 /**
573  * \brief Face Flip Normal
574  *
575  * Reverses the winding of a face.
576  * \note This updates the calculated normal.
577  */
578 void BM_face_normal_flip(BMesh *bm, BMFace *f)
579 {
580         bmesh_loop_reverse(bm, f);
581         negate_v3(f->no);
582 }
583
584 /* detects if two line segments cross each other (intersects).
585  * note, there could be more winding cases then there needs to be. */
586 static bool line_crosses_v2f(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
587 {
588
589 #define GETMIN2_AXIS(a, b, ma, mb, axis)   \
590         {                                      \
591                 ma[axis] = min_ff(a[axis], b[axis]); \
592                 mb[axis] = max_ff(a[axis], b[axis]); \
593         } (void)0
594
595 #define GETMIN2(a, b, ma, mb)          \
596         {                                  \
597                 GETMIN2_AXIS(a, b, ma, mb, 0); \
598                 GETMIN2_AXIS(a, b, ma, mb, 1); \
599         } (void)0
600
601 #define EPS (FLT_EPSILON * 15)
602
603         int w1, w2, w3, w4, w5 /*, re */;
604         float mv1[2], mv2[2], mv3[2], mv4[2];
605         
606         /* now test winding */
607         w1 = testedgesidef(v1, v3, v2);
608         w2 = testedgesidef(v2, v4, v1);
609         w3 = !testedgesidef(v1, v2, v3);
610         w4 = testedgesidef(v3, v2, v4);
611         w5 = !testedgesidef(v3, v1, v4);
612         
613         if (w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5) {
614                 return true;
615         }
616         
617         GETMIN2(v1, v2, mv1, mv2);
618         GETMIN2(v3, v4, mv3, mv4);
619         
620         /* do an interval test on the x and y axes */
621         /* first do x axis */
622         if (fabsf(v1[1] - v2[1]) < EPS &&
623             fabsf(v3[1] - v4[1]) < EPS &&
624             fabsf(v1[1] - v3[1]) < EPS)
625         {
626                 return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
627         }
628
629         /* now do y axis */
630         if (fabsf(v1[0] - v2[0]) < EPS &&
631             fabsf(v3[0] - v4[0]) < EPS &&
632             fabsf(v1[0] - v3[0]) < EPS)
633         {
634                 return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
635         }
636
637         return false;
638
639 #undef GETMIN2_AXIS
640 #undef GETMIN2
641 #undef EPS
642
643 }
644
645 /**
646  *  BM POINT IN FACE
647  *
648  * Projects co onto face f, and returns true if it is inside
649  * the face bounds.
650  *
651  * \note this uses a best-axis projection test,
652  * instead of projecting co directly into f's orientation space,
653  * so there might be accuracy issues.
654  */
655 bool BM_face_point_inside_test(BMFace *f, const float co[3])
656 {
657         int ax, ay;
658         float co2[2], cent[2] = {0.0f, 0.0f}, out[2] = {FLT_MAX * 0.5f, FLT_MAX * 0.5f};
659         BMLoop *l_iter;
660         BMLoop *l_first;
661         int crosses = 0;
662         float onepluseps = 1.0f + (float)FLT_EPSILON * 150.0f;
663         
664         if (dot_v3v3(f->no, f->no) <= FLT_EPSILON * 10)
665                 BM_face_normal_update(f);
666         
667         /* find best projection of face XY, XZ or YZ: barycentric weights of
668          * the 2d projected coords are the same and faster to compute
669          *
670          * this probably isn't all that accurate, but it has the advantage of
671          * being fast (especially compared to projecting into the face orientation)
672          */
673         axis_dominant_v3(&ax, &ay, f->no);
674
675         co2[0] = co[ax];
676         co2[1] = co[ay];
677         
678         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
679         do {
680                 cent[0] += l_iter->v->co[ax];
681                 cent[1] += l_iter->v->co[ay];
682         } while ((l_iter = l_iter->next) != l_first);
683         
684         mul_v2_fl(cent, 1.0f / (float)f->len);
685         
686         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
687         do {
688                 float v1[2], v2[2];
689                 
690                 v1[0] = (l_iter->prev->v->co[ax] - cent[0]) * onepluseps + cent[0];
691                 v1[1] = (l_iter->prev->v->co[ay] - cent[1]) * onepluseps + cent[1];
692                 
693                 v2[0] = (l_iter->v->co[ax] - cent[0]) * onepluseps + cent[0];
694                 v2[1] = (l_iter->v->co[ay] - cent[1]) * onepluseps + cent[1];
695                 
696                 crosses += line_crosses_v2f(v1, v2, co2, out) != 0;
697         } while ((l_iter = l_iter->next) != l_first);
698         
699         return crosses % 2 != 0;
700 }
701
702 static bool bm_face_goodline(float const (*projectverts)[2], BMFace *f, int v1i, int v2i, int v3i)
703 {
704         BMLoop *l_iter;
705         BMLoop *l_first;
706
707         float pv1[2];
708         const float *v1 = projectverts[v1i];
709         const float *v2 = projectverts[v2i];
710         const float *v3 = projectverts[v3i];
711         int i;
712
713         /* v3 must be on the left side of [v1, v2] line, else we know [v1, v3] is outside of f! */
714         if (testedgesidef(v1, v2, v3)) {
715                 return false;
716         }
717
718         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
719         do {
720                 i = BM_elem_index_get(l_iter->v);
721                 copy_v2_v2(pv1, projectverts[i]);
722
723                 if (ELEM3(i, v1i, v2i, v3i)) {
724 #if 0
725                         printf("%d in (%d, %d, %d) tri (from indices!), continuing\n", i, v1i, v2i, v3i);
726 #endif
727                         continue;
728                 }
729
730                 if (isect_point_tri_v2(pv1, v1, v2, v3) || isect_point_tri_v2(pv1, v3, v2, v1)) {
731 #if 0
732                         if (isect_point_tri_v2(pv1, v1, v2, v3))
733                                 printf("%d in (%d, %d, %d)\n", v3i, i, v1i, v2i);
734                         else
735                                 printf("%d in (%d, %d, %d)\n", v1i, i, v3i, v2i);
736 #endif
737                         return false;
738                 }
739         } while ((l_iter = l_iter->next) != l_first);
740         return true;
741 }
742
743 /**
744  * \brief Find Ear
745  *
746  * Used by tessellator to find the next triangle to 'clip off' of a polygon while tessellating.
747  *
748  * \param f The face to search.
749  * \param projectverts an array of face vert coords.
750  * \param use_beauty Currently only applies to quads, can be extended later on.
751  * \param abscoss Must be allocated by caller, and at least f->len length
752  *        (allow to avoid allocating a new one for each tri!).
753  */
754 static BMLoop *poly_find_ear(BMFace *f, float (*projectverts)[2], const bool use_beauty, float *abscoss)
755 {
756         BMLoop *bestear = NULL;
757
758         BMLoop *l_iter;
759         BMLoop *l_first;
760
761         const float cos_threshold = 0.9f;
762         const float bias = 1.0f + 1e-6f;
763
764         if (f->len == 4) {
765                 BMLoop *larr[4];
766                 int i = 0, i4;
767                 float cos1, cos2;
768                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
769                 do {
770                         larr[i] = l_iter;
771                         i++;
772                 } while ((l_iter = l_iter->next) != l_first);
773
774                 /* pick 0/1 based on best lenth */
775                 /* XXX Can't only rely on such test, also must check we do not get (too much) degenerated triangles!!! */
776                 i = (((len_squared_v3v3(larr[0]->v->co, larr[2]->v->co) >
777                      len_squared_v3v3(larr[1]->v->co, larr[3]->v->co) * bias)) != use_beauty);
778                 i4 = (i + 3) % 4;
779                 /* Check produced tris aren't too flat/narrow...
780                  * Probably not the best test, but is quite efficient and should at least avoid null-area faces! */
781                 cos1 = fabsf(cos_v3v3v3(larr[i4]->v->co, larr[i]->v->co, larr[i + 1]->v->co));
782                 cos2 = fabsf(cos_v3v3v3(larr[i4]->v->co, larr[i + 2]->v->co, larr[i + 1]->v->co));
783 #if 0
784                 printf("%d, (%f, %f), (%f, %f)\n", i, cos1, cos2,
785                        fabsf(cos_v3v3v3(larr[i]->v->co, larr[i4]->v->co, larr[i + 2]->v->co)),
786                        fabsf(cos_v3v3v3(larr[i]->v->co, larr[i + 1]->v->co, larr[i + 2]->v->co)));
787 #endif
788                 if (cos1 < cos2)
789                         cos1 = cos2;
790                 if (cos1 > cos_threshold) {
791                         if (cos1 > fabsf(cos_v3v3v3(larr[i]->v->co, larr[i4]->v->co, larr[i + 2]->v->co)) &&
792                             cos1 > fabsf(cos_v3v3v3(larr[i]->v->co, larr[i + 1]->v->co, larr[i + 2]->v->co)))
793                         {
794                                 i = !i;
795                         }
796                 }
797                 /* Last check we do not get overlapping triangles
798                  * (as much as possible, there are some cases with no good solution!) */
799                 i4 = (i + 3) % 4;
800                 if (!bm_face_goodline((float const (*)[2])projectverts, f, BM_elem_index_get(larr[i4]->v),
801                                       BM_elem_index_get(larr[i]->v), BM_elem_index_get(larr[i + 1]->v)))
802                 {
803                         i = !i;
804                 }
805 /*              printf("%d\n", i);*/
806                 bestear = larr[i];
807
808         }
809         else {
810                 /* float angle, bestangle = 180.0f; */
811                 float cos, bestcos = 1.0f;
812                 int i, j, len;
813
814                 /* Compute cos of all corners! */
815                 i = 0;
816                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
817                 len = l_iter->f->len;
818                 do {
819                         const BMVert *v1 = l_iter->prev->v;
820                         const BMVert *v2 = l_iter->v;
821                         const BMVert *v3 = l_iter->next->v;
822
823                         abscoss[i] = fabsf(cos_v3v3v3(v1->co, v2->co, v3->co));
824 /*                      printf("tcoss: %f\n", *tcoss);*/
825                         i++;
826                 } while ((l_iter = l_iter->next) != l_first);
827
828                 i = 0;
829                 l_iter = l_first;
830                 do {
831                         const BMVert *v1 = l_iter->prev->v;
832                         const BMVert *v2 = l_iter->v;
833                         const BMVert *v3 = l_iter->next->v;
834
835                         if (bm_face_goodline((float const (*)[2])projectverts, f,
836                                              BM_elem_index_get(v1), BM_elem_index_get(v2), BM_elem_index_get(v3)))
837                         {
838                                 /* Compute highest cos (i.e. narrowest angle) of this tri. */
839                                 cos = max_fff(abscoss[i],
840                                               fabsf(cos_v3v3v3(v2->co, v3->co, v1->co)),
841                                               fabsf(cos_v3v3v3(v3->co, v1->co, v2->co)));
842
843                                 /* Compare to prev best (i.e. lowest) cos. */
844                                 if (cos < bestcos) {
845                                         /* We must check this tri would not leave a (too much) degenerated remaining face! */
846                                         /* For now just assume if the average of cos of all
847                                          * "remaining face"'s corners is below a given threshold, it's OK. */
848                                         float avgcos = fabsf(cos_v3v3v3(v1->co, v3->co, l_iter->next->next->v->co));
849                                         const int i_limit = (i - 1 + len) % len;
850                                         avgcos += fabsf(cos_v3v3v3(l_iter->prev->prev->v->co, v1->co, v3->co));
851                                         j = (i + 2) % len;
852                                         do {
853                                                 avgcos += abscoss[j];
854                                         } while ((j = (j + 1) % len) != i_limit);
855                                         avgcos /= len - 1;
856
857                                         /* We need a best ear in any case... */
858                                         if (avgcos < cos_threshold || (!bestear && avgcos < 1.0f)) {
859                                                 /* OKI, keep this ear (corner...) as a potential best one! */
860                                                 bestear = l_iter;
861                                                 bestcos = cos;
862                                         }
863 #if 0
864                                         else
865                                                 printf("Had a nice tri (higest cos of %f, current bestcos is %f), "
866                                                        "but average cos of all \"remaining face\"'s corners is too high (%f)!\n",
867                                                        cos, bestcos, avgcos);
868 #endif
869                                 }
870                         }
871                         i++;
872                 } while ((l_iter = l_iter->next) != l_first);
873         }
874
875         return bestear;
876 }
877
878 /**
879  * \brief BMESH TRIANGULATE FACE
880  *
881  * Currently repeatedly find the best triangle (i.e. the most "open" one), provided it does not
882  * produces a "remaining" face with too much wide/narrow angles
883  * (using cos (i.e. dot product of normalized vectors) of angles).
884  *
885  * \param r_faces_new if non-null, must be an array of BMFace pointers,
886  * with a length equal to (f->len - 2). It will be filled with the new
887  * triangles.
888  *
889  * \note use_tag tags new flags and edges.
890  */
891 void BM_face_triangulate(BMesh *bm, BMFace *f,
892                          BMFace **r_faces_new,
893                          const bool use_beauty, const bool use_tag)
894 {
895         const float f_len_orig = f->len;
896         int i, nf_i = 0;
897         BMLoop *l_new;
898         BMLoop *l_iter;
899         BMLoop *l_first;
900         /* BM_face_triangulate: temp absolute cosines of face corners */
901         float (*projectverts)[2] = BLI_array_alloca(projectverts, f_len_orig);
902         float *abscoss = BLI_array_alloca(abscoss, f_len_orig);
903         float mat[3][3];
904
905         axis_dominant_v3_to_m3(mat, f->no);
906
907         /* copy vertex coordinates to vertspace area */
908         i = 0;
909         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
910         do {
911                 mul_v2_m3v3(projectverts[i], mat, l_iter->v->co);
912                 BM_elem_index_set(l_iter->v, i++); /* set dirty! */
913         } while ((l_iter = l_iter->next) != l_first);
914
915         bm->elem_index_dirty |= BM_VERT; /* see above */
916
917         while (f->len > 3) {
918                 l_iter = poly_find_ear(f, projectverts, use_beauty, abscoss);
919
920                 /* force triangulation - if we can't find an ear the face is degenerate */
921                 if (l_iter == NULL) {
922                         l_iter = BM_FACE_FIRST_LOOP(f);
923                 }
924
925 /*              printf("Subdividing face...\n");*/
926                 f = BM_face_split(bm, l_iter->f, l_iter->prev->v, l_iter->next->v, &l_new, NULL, true);
927
928                 if (UNLIKELY(!f)) {
929                         fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
930                         break;
931                 }
932
933                 copy_v3_v3(f->no, l_iter->f->no);
934
935                 if (use_tag) {
936                         BM_elem_flag_enable(l_new->e, BM_ELEM_TAG);
937                         BM_elem_flag_enable(f, BM_ELEM_TAG);
938                 }
939
940                 if (r_faces_new) {
941                         r_faces_new[nf_i++] = f;
942                 }
943         }
944
945         BLI_assert(f->len == 3);
946 }
947
948 /**
949  * each pair of loops defines a new edge, a split.  this function goes
950  * through and sets pairs that are geometrically invalid to null.  a
951  * split is invalid, if it forms a concave angle or it intersects other
952  * edges in the face, or it intersects another split.  in the case of
953  * intersecting splits, only the first of the set of intersecting
954  * splits survives
955  */
956 void BM_face_legal_splits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
957 {
958         BMIter iter;
959         BMLoop *l;
960         float v1[3], v2[3], v3[3] /*, v4[3 */, no[3], mid[3], *p1, *p2, *p3, *p4;
961         float out[3] = {-FLT_MAX, -FLT_MAX, 0.0f};
962         float (*projverts)[3] = BLI_array_alloca(projverts, f->len);
963         float (*edgeverts)[3] = BLI_array_alloca(edgeverts, len * 2);
964         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
965         int i, j, a = 0, clen;
966         
967         i = 0;
968         l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
969         for ( ; l; l = BM_iter_step(&iter)) {
970                 BM_elem_index_set(l, i); /* set_loop */
971                 copy_v3_v3(projverts[i], l->v->co);
972                 i++;
973         }
974         
975         for (i = 0; i < len; i++) {
976                 copy_v3_v3(v1, loops[i][0]->v->co);
977                 copy_v3_v3(v2, loops[i][1]->v->co);
978
979                 scale_edge_v3f(v1, v2, fac2);
980                 
981                 copy_v3_v3(edgeverts[a], v1);
982                 a++;
983                 copy_v3_v3(edgeverts[a], v2);
984                 a++;
985         }
986         
987         calc_poly_normal(no, projverts, f->len);
988         poly_rotate_plane(no, projverts, f->len);
989         poly_rotate_plane(no, edgeverts, len * 2);
990
991         for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
992                 p1 = projverts[i];
993                 out[0] = max_ff(out[0], p1[0]);
994                 out[1] = max_ff(out[1], p1[1]);
995                 /* out[2] = 0.0f; */ /* keep at zero */
996
997                 p1[2] = 0.0f;
998         }
999         
1000         /* ensure we are well outside the face bounds (value is arbitrary) */
1001         add_v2_fl(out, 1.0f);
1002
1003         for (i = 0; i < len; i++) {
1004                 edgeverts[i * 2][2] = 0.0f;
1005                 edgeverts[i * 2 + 1][2] = 0.0f;
1006         }
1007
1008         /* do convexity test */
1009         for (i = 0; i < len; i++) {
1010                 copy_v3_v3(v2, edgeverts[i * 2]);
1011                 copy_v3_v3(v3, edgeverts[i * 2 + 1]);
1012
1013                 mid_v3_v3v3(mid, v2, v3);
1014                 
1015                 clen = 0;
1016                 for (j = 0; j < f->len; j++) {
1017                         p1 = projverts[j];
1018                         p2 = projverts[(j + 1) % f->len];
1019                         
1020 #if 0
1021                         copy_v3_v3(v1, p1);
1022                         copy_v3_v3(v2, p2);
1023
1024                         scale_edge_v3f(v1, v2, fac1);
1025                         if (line_crosses_v2f(v1, v2, mid, out)) {
1026                                 clen++;
1027                         }
1028 #else
1029                         if (line_crosses_v2f(p1, p2, mid, out)) {
1030                                 clen++;
1031                         }
1032 #endif
1033                 }
1034
1035                 if (clen % 2 == 0) {
1036                         loops[i][0] = NULL;
1037                 }
1038         }
1039
1040         /* do line crossing tests */
1041         for (i = 0; i < f->len; i++) {
1042                 p1 = projverts[i];
1043                 p2 = projverts[(i + 1) % f->len];
1044                 
1045                 copy_v3_v3(v1, p1);
1046                 copy_v3_v3(v2, p2);
1047
1048                 scale_edge_v3f(v1, v2, fac1);
1049
1050                 for (j = 0; j < len; j++) {
1051                         if (!loops[j][0]) {
1052                                 continue;
1053                         }
1054
1055                         p3 = edgeverts[j * 2];
1056                         p4 = edgeverts[j * 2 + 1];
1057
1058                         if (line_crosses_v2f(v1, v2, p3, p4)) {
1059                                 loops[j][0] = NULL;
1060                         }
1061                 }
1062         }
1063
1064         for (i = 0; i < len; i++) {
1065                 for (j = 0; j < len; j++) {
1066                         if (j != i && loops[i][0] && loops[j][0]) {
1067                                 p1 = edgeverts[i * 2];
1068                                 p2 = edgeverts[i * 2 + 1];
1069                                 p3 = edgeverts[j * 2];
1070                                 p4 = edgeverts[j * 2 + 1];
1071
1072                                 copy_v3_v3(v1, p1);
1073                                 copy_v3_v3(v2, p2);
1074
1075                                 scale_edge_v3f(v1, v2, fac1);
1076
1077                                 if (line_crosses_v2f(v1, v2, p3, p4)) {
1078                                         loops[i][0] = NULL;
1079                                 }
1080                         }
1081                 }
1082         }
1083 }
1084
1085
1086 /**
1087  * Small utility functions for fast access
1088  *
1089  * faster alternative to:
1090  *  BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void**)v, 3);
1091  */
1092 void BM_face_as_array_vert_tri(BMFace *f, BMVert *r_verts[3])
1093 {
1094         BMLoop *l = BM_FACE_FIRST_LOOP(f);
1095
1096         BLI_assert(f->len == 3);
1097
1098         r_verts[0] = l->v; l = l->next;
1099         r_verts[1] = l->v; l = l->next;
1100         r_verts[2] = l->v;
1101 }
1102
1103 /**
1104  * faster alternative to:
1105  *  BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void**)v, 4);
1106  */
1107 void BM_face_as_array_vert_quad(BMFace *f, BMVert *r_verts[4])
1108 {
1109         BMLoop *l = BM_FACE_FIRST_LOOP(f);
1110
1111         BLI_assert(f->len == 4);
1112
1113         r_verts[0] = l->v; l = l->next;
1114         r_verts[1] = l->v; l = l->next;
1115         r_verts[2] = l->v; l = l->next;
1116         r_verts[3] = l->v;
1117 }