code cleanup: change C naming convention (so py and C api match), eg:
[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                 DO_MINMAX(l_iter->v->co, min, max);
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
449         /* common cases first */
450         switch (f->len) {
451                 case 4:
452                 {
453                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
454                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
455                         const float *co3 = vertexCos[BM_elem_index_get((l = l->next)->v)];
456                         const float *co4 = vertexCos[BM_elem_index_get((l->next)->v)];
457
458                         normal_quad_v3(no, co1, co2, co3, co4);
459                         break;
460                 }
461                 case 3:
462                 {
463                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
464                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
465                         const float *co3 = vertexCos[BM_elem_index_get((l->next)->v)];
466
467                         normal_tri_v3(no, co1, co2, co3);
468                         break;
469                 }
470                 case 0:
471                 {
472                         zero_v3(no);
473                         break;
474                 }
475                 default:
476                 {
477                         bm_face_calc_poly_normal_vertex_cos(f, no, vertexCos);
478                         break;
479                 }
480         }
481 }
482
483 /**
484  * \brief Face Flip Normal
485  *
486  * Reverses the winding of a face.
487  * \note This updates the calculated normal.
488  */
489 void BM_face_normal_flip(BMesh *bm, BMFace *f)
490 {
491         bmesh_loop_reverse(bm, f);
492         negate_v3(f->no);
493 }
494
495 /* detects if two line segments cross each other (intersects).
496  * note, there could be more winding cases then there needs to be. */
497 static int linecrossesf(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
498 {
499
500 #define GETMIN2_AXIS(a, b, ma, mb, axis)   \
501         {                                      \
502                 ma[axis] = MIN2(a[axis], b[axis]); \
503                 mb[axis] = MAX2(a[axis], b[axis]); \
504         } (void)0
505
506 #define GETMIN2(a, b, ma, mb)          \
507         {                                  \
508                 GETMIN2_AXIS(a, b, ma, mb, 0); \
509                 GETMIN2_AXIS(a, b, ma, mb, 1); \
510         } (void)0
511
512 #define EPS (FLT_EPSILON * 15)
513
514         int w1, w2, w3, w4, w5 /*, re */;
515         float mv1[2], mv2[2], mv3[2], mv4[2];
516         
517         /* now test winding */
518         w1 = testedgesidef(v1, v3, v2);
519         w2 = testedgesidef(v2, v4, v1);
520         w3 = !testedgesidef(v1, v2, v3);
521         w4 = testedgesidef(v3, v2, v4);
522         w5 = !testedgesidef(v3, v1, v4);
523         
524         if (w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5) {
525                 return TRUE;
526         }
527         
528         GETMIN2(v1, v2, mv1, mv2);
529         GETMIN2(v3, v4, mv3, mv4);
530         
531         /* do an interval test on the x and y axes */
532         /* first do x axis */
533         if (ABS(v1[1] - v2[1]) < EPS &&
534             ABS(v3[1] - v4[1]) < EPS &&
535             ABS(v1[1] - v3[1]) < EPS)
536         {
537                 return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
538         }
539
540         /* now do y axis */
541         if (ABS(v1[0] - v2[0]) < EPS &&
542             ABS(v3[0] - v4[0]) < EPS &&
543             ABS(v1[0] - v3[0]) < EPS)
544         {
545                 return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
546         }
547
548         return FALSE;
549
550 #undef GETMIN2_AXIS
551 #undef GETMIN2
552 #undef EPS
553
554 }
555
556 /**
557  *  BM POINT IN FACE
558  *
559  * Projects co onto face f, and returns true if it is inside
560  * the face bounds.
561  *
562  * \note this uses a best-axis projection test,
563  * instead of projecting co directly into f's orientation space,
564  * so there might be accuracy issues.
565  */
566 int BM_face_point_inside_test(BMFace *f, const float co[3])
567 {
568         int ax, ay;
569         float co2[2], cent[2] = {0.0f, 0.0f}, out[2] = {FLT_MAX * 0.5f, FLT_MAX * 0.5f};
570         BMLoop *l_iter;
571         BMLoop *l_first;
572         int crosses = 0;
573         float onepluseps = 1.0f + (float)FLT_EPSILON * 150.0f;
574         
575         if (dot_v3v3(f->no, f->no) <= FLT_EPSILON * 10)
576                 BM_face_normal_update(f);
577         
578         /* find best projection of face XY, XZ or YZ: barycentric weights of
579          * the 2d projected coords are the same and faster to compute
580          *
581          * this probably isn't all that accurate, but it has the advantage of
582          * being fast (especially compared to projecting into the face orientation)
583          */
584         axis_dominant_v3(&ax, &ay, f->no);
585
586         co2[0] = co[ax];
587         co2[1] = co[ay];
588         
589         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
590         do {
591                 cent[0] += l_iter->v->co[ax];
592                 cent[1] += l_iter->v->co[ay];
593         } while ((l_iter = l_iter->next) != l_first);
594         
595         mul_v2_fl(cent, 1.0f / (float)f->len);
596         
597         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
598         do {
599                 float v1[2], v2[2];
600                 
601                 v1[0] = (l_iter->prev->v->co[ax] - cent[ax]) * onepluseps + cent[ax];
602                 v1[1] = (l_iter->prev->v->co[ay] - cent[ay]) * onepluseps + cent[ay];
603                 
604                 v2[0] = (l_iter->v->co[ax] - cent[ax]) * onepluseps + cent[ax];
605                 v2[1] = (l_iter->v->co[ay] - cent[ay]) * onepluseps + cent[ay];
606                 
607                 crosses += linecrossesf(v1, v2, co2, out) != 0;
608         } while ((l_iter = l_iter->next) != l_first);
609         
610         return crosses % 2 != 0;
611 }
612
613 static int bm_face_goodline(float const (*projectverts)[3], BMFace *f,
614                             int v1i, int v2i, int v3i,
615                             int UNUSED(nvert))
616 {
617         BMLoop *l_iter;
618         BMLoop *l_first;
619         float v1[3], v2[3], v3[3], pv1[3], pv2[3];
620         int i;
621
622         copy_v3_v3(v1, projectverts[v1i]);
623         copy_v3_v3(v2, projectverts[v2i]);
624         copy_v3_v3(v3, projectverts[v3i]);
625         
626         if (testedgesidef(v1, v2, v3)) {
627                 return FALSE;
628         }
629
630         //for (i = 0; i < nvert; i++) {
631         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
632         do {
633                 i = BM_elem_index_get(l_iter->v);
634                 if (i == v1i || i == v2i || i == v3i) {
635                         continue;
636                 }
637                 
638                 copy_v3_v3(pv1, projectverts[BM_elem_index_get(l_iter->v)]);
639                 copy_v3_v3(pv2, projectverts[BM_elem_index_get(l_iter->next->v)]);
640                 
641                 //if (linecrossesf(pv1, pv2, v1, v3)) return FALSE;
642
643                 if (isect_point_tri_v2(pv1, v1, v2, v3) ||
644                     isect_point_tri_v2(pv1, v3, v2, v1))
645                 {
646                         return FALSE;
647                 }
648         } while ((l_iter = l_iter->next) != l_first);
649         return TRUE;
650 }
651
652 /**
653  * \brief Find Ear
654  *
655  * Used by tessellator to find
656  * the next triangle to 'clip off'
657  * of a polygon while tessellating.
658  *
659  * \param use_beauty Currently only applies to quads, can be extended later on.
660  */
661 static BMLoop *find_ear(BMFace *f, float (*verts)[3], const int nvert, const int use_beauty)
662 {
663         BMLoop *bestear = NULL;
664
665         BMLoop *l_iter;
666         BMLoop *l_first;
667
668         if (f->len == 4) {
669                 BMLoop *larr[4];
670                 int i = 0;
671
672                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
673                 do {
674                         larr[i] = l_iter;
675                         i++;
676                 } while ((l_iter = l_iter->next) != l_first);
677
678                 /* pick 0/1 based on best lenth */
679                 bestear = larr[(((len_squared_v3v3(larr[0]->v->co, larr[2]->v->co) >
680                                   len_squared_v3v3(larr[1]->v->co, larr[3]->v->co))) != use_beauty)];
681
682         }
683         else {
684                 BMVert *v1, *v2, *v3;
685
686                 /* float angle, bestangle = 180.0f; */
687                 int isear /*, i = 0 */;
688
689                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
690                 do {
691                         isear = TRUE;
692
693                         v1 = l_iter->prev->v;
694                         v2 = l_iter->v;
695                         v3 = l_iter->next->v;
696
697                         if (BM_edge_exists(v1, v3)) {
698                                 isear = FALSE;
699                         }
700                         else if (!bm_face_goodline((float const (*)[3])verts, f,
701                                                    BM_elem_index_get(v1), BM_elem_index_get(v2), BM_elem_index_get(v3),
702                                                    nvert))
703                         {
704                                 isear = FALSE;
705                         }
706
707                         if (isear) {
708         #if 0
709                                 /* if this code comes back, it needs to be converted to radians */
710                                 angle = angle_v3v3v3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
711                                 if (!bestear || ABS(angle - 45.0f) < bestangle) {
712                                         bestear = l;
713                                         bestangle = ABS(45.0f - angle);
714                                 }
715
716                                 if (angle > 20 && angle < 90) break;
717                                 if (angle < 100 && i > 5) break;
718                                 i += 1;
719         #endif
720
721                                 bestear = l_iter;
722                                 break;
723                         }
724                 } while ((l_iter = l_iter->next) != l_first);
725         }
726
727         return bestear;
728 }
729
730 /**
731  * \brief BMESH TRIANGULATE FACE
732  *
733  * Triangulates a face using a simple 'ear clipping' algorithm that tries to
734  * favor non-skinny triangles (angles less than 90 degrees).
735  *
736  * If the triangulator has bits left over (or cannot triangulate at all)
737  * it uses a simple fan triangulation,
738  *
739  * newfaces, if non-null, must be an array of BMFace pointers,
740  * with a length equal to f->len.  it will be filled with the new
741  * triangles, and will be NULL-terminated.
742  *
743  * \note newedgeflag sets a flag layer flag, obviously not the header flag.
744  */
745 void BM_face_triangulate(BMesh *bm, BMFace *f, float (*projectverts)[3],
746                          const short newedge_oflag, const short newface_oflag, BMFace **newfaces,
747                          const short use_beauty)
748 {
749         int i, done, nvert, nf_i = 0;
750         BMLoop *newl, *nextloop;
751         BMLoop *l_iter;
752         BMLoop *l_first;
753         /* BMVert *v; */ /* UNUSED */
754
755         /* copy vertex coordinates to vertspace arra */
756         i = 0;
757         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
758         do {
759                 copy_v3_v3(projectverts[i], l_iter->v->co);
760                 BM_elem_index_set(l_iter->v, i); /* set dirty! */
761                 i++;
762         } while ((l_iter = l_iter->next) != l_first);
763
764         bm->elem_index_dirty |= BM_VERT; /* see above */
765
766         ///bmesh_face_normal_update(bm, f, f->no, projectverts);
767
768         calc_poly_normal(f->no, projectverts, f->len);
769         poly_rotate_plane(f->no, projectverts, i);
770
771         nvert = f->len;
772
773         //calc_poly_plane(projectverts, i);
774         for (i = 0; i < nvert; i++) {
775                 projectverts[i][2] = 0.0f;
776         }
777
778         done = 0;
779         while (!done && f->len > 3) {
780                 done = 1;
781                 l_iter = find_ear(f, projectverts, nvert, use_beauty);
782                 if (l_iter) {
783                         done = 0;
784                         /* v = l->v; */ /* UNUSED */
785                         f = BM_face_split(bm, l_iter->f, l_iter->prev->v,
786                                           l_iter->next->v,
787                                           &newl, NULL, TRUE);
788
789                         if (UNLIKELY(!f)) {
790                                 fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
791                                 break;
792                         }
793
794                         copy_v3_v3(f->no, l_iter->f->no);
795                         BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
796                         BMO_elem_flag_enable(bm, f, newface_oflag);
797                         
798                         if (newfaces) newfaces[nf_i++] = f;
799
800 #if 0
801                         l = f->loopbase;
802                         do {
803                                 if (l->v == v) {
804                                         f->loopbase = l;
805                                         break;
806                                 }
807                                 l = l->next;
808                         } while (l != f->loopbase);
809 #endif
810
811                 }
812         }
813
814         if (f->len > 3) {
815                 l_iter = BM_FACE_FIRST_LOOP(f);
816                 while (l_iter->f->len > 3) {
817                         nextloop = l_iter->next->next;
818                         f = BM_face_split(bm, l_iter->f, l_iter->v, nextloop->v,
819                                           &newl, NULL, TRUE);
820                         if (!f) {
821                                 printf("triangle fan step of triangulator failed.\n");
822
823                                 /* NULL-terminate */
824                                 if (newfaces) newfaces[nf_i] = NULL;
825                                 return;
826                         }
827
828                         if (newfaces) newfaces[nf_i++] = f;
829                         
830                         BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
831                         BMO_elem_flag_enable(bm, f, newface_oflag);
832                         l_iter = nextloop;
833                 }
834         }
835         
836         /* NULL-terminate */
837         if (newfaces) newfaces[nf_i] = NULL;
838 }
839
840 /**
841  * each pair of loops defines a new edge, a split.  this function goes
842  * through and sets pairs that are geometrically invalid to null.  a
843  * split is invalid, if it forms a concave angle or it intersects other
844  * edges in the face, or it intersects another split.  in the case of
845  * intersecting splits, only the first of the set of intersecting
846  * splits survives
847  */
848 void BM_face_legal_splits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
849 {
850         BMIter iter;
851         BMLoop *l;
852         float v1[3], v2[3], v3[3] /*, v4[3 */, no[3], mid[3], *p1, *p2, *p3, *p4;
853         float out[3] = {-234324.0f, -234324.0f, 0.0f};
854         float (*projverts)[3];
855         float (*edgeverts)[3];
856         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
857         int i, j, a = 0, clen;
858
859         BLI_array_fixedstack_declare(projverts, BM_NGON_STACK_SIZE, f->len,      "projvertsb");
860         BLI_array_fixedstack_declare(edgeverts, BM_NGON_STACK_SIZE * 2, len * 2, "edgevertsb");
861         
862         i = 0;
863         l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
864         for ( ; l; l = BM_iter_step(&iter)) {
865                 BM_elem_index_set(l, i); /* set_loop */
866                 copy_v3_v3(projverts[i], l->v->co);
867                 i++;
868         }
869         
870         for (i = 0; i < len; i++) {
871                 copy_v3_v3(v1, loops[i][0]->v->co);
872                 copy_v3_v3(v2, loops[i][1]->v->co);
873
874                 shrink_edgef(v1, v2, fac2);
875                 
876                 copy_v3_v3(edgeverts[a], v1);
877                 a++;
878                 copy_v3_v3(edgeverts[a], v2);
879                 a++;
880         }
881         
882         calc_poly_normal(no, projverts, f->len);
883         poly_rotate_plane(no, projverts, f->len);
884         poly_rotate_plane(no, edgeverts, len * 2);
885
886         for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
887                 p1 = projverts[i];
888                 out[0] = MAX2(out[0], p1[0]) + 0.01f;
889                 out[1] = MAX2(out[1], p1[1]) + 0.01f;
890                 out[2] = 0.0f;
891                 p1[2] = 0.0f;
892
893                 //copy_v3_v3(l->v->co, p1);
894         }
895         
896         for (i = 0; i < len; i++) {
897                 edgeverts[i * 2][2] = 0.0f;
898                 edgeverts[i * 2 + 1][2] = 0.0f;
899         }
900
901         /* do convexity test */
902         for (i = 0; i < len; i++) {
903                 copy_v3_v3(v2, edgeverts[i * 2]);
904                 copy_v3_v3(v3, edgeverts[i * 2 + 1]);
905
906                 mid_v3_v3v3(mid, v2, v3);
907                 
908                 clen = 0;
909                 for (j = 0; j < f->len; j++) {
910                         p1 = projverts[j];
911                         p2 = projverts[(j + 1) % f->len];
912                         
913                         copy_v3_v3(v1, p1);
914                         copy_v3_v3(v2, p2);
915
916                         shrink_edgef(v1, v2, fac1);
917
918                         if (linecrossesf(p1, p2, mid, out)) clen++;
919                 }
920                 
921                 if (clen % 2 == 0) {
922                         loops[i][0] = NULL;
923                 }
924         }
925         
926         /* do line crossing test */
927         for (i = 0; i < f->len; i++) {
928                 p1 = projverts[i];
929                 p2 = projverts[(i + 1) % f->len];
930                 
931                 copy_v3_v3(v1, p1);
932                 copy_v3_v3(v2, p2);
933
934                 shrink_edgef(v1, v2, fac1);
935
936                 for (j = 0; j < len; j++) {
937                         if (!loops[j][0]) {
938                                 continue;
939                         }
940
941                         p3 = edgeverts[j * 2];
942                         p4 = edgeverts[j * 2 + 1];
943
944                         if (linecrossesf(v1, v2, p3, p4)) {
945                                 loops[j][0] = NULL;
946                         }
947                 }
948         }
949
950         for (i = 0; i < len; i++) {
951                 for (j = 0; j < len; j++) {
952                         if (j != i && loops[i][0] && loops[j][0]) {
953                                 p1 = edgeverts[i * 2];
954                                 p2 = edgeverts[i * 2 + 1];
955                                 p3 = edgeverts[j * 2];
956                                 p4 = edgeverts[j * 2 + 1];
957
958                                 copy_v3_v3(v1, p1);
959                                 copy_v3_v3(v2, p2);
960
961                                 shrink_edgef(v1, v2, fac1);
962
963                                 if (linecrossesf(v1, v2, p3, p4)) {
964                                         loops[i][0] = NULL;
965                                 }
966                         }
967                 }
968         }
969
970         BLI_array_fixedstack_free(projverts);
971         BLI_array_fixedstack_free(edgeverts);
972 }