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