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