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