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