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