remove BM_ITER, BM_ITER_INDEX macros, use ELEM or MESH variants only (the maceros...
[blender-staging.git] / source / blender / bmesh / intern / bmesh_polygon.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): Joseph Eagar, Geoffrey Bantle, Campbell Barton
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/intern/bmesh_polygon.c
24  *  \ingroup bmesh
25  *
26  * This file contains code for dealing
27  * with polygons (normal/area calculation,
28  * tessellation, etc)
29  *
30  * BMESH_TODO:
31  *  - Add in Tessellator frontend that creates
32  *    BMTriangles from copied faces
33  *
34  *  - Add in Function that checks for and flags
35  *    degenerate faces.
36  */
37
38 #include "BLI_math.h"
39 #include "BLI_array.h"
40
41 #include "MEM_guardedalloc.h"
42
43 #include "bmesh.h"
44 #include "intern/bmesh_private.h"
45
46 /**
47  * \brief TEST EDGE SIDE and POINT IN TRIANGLE
48  *
49  * Point in triangle tests stolen from scanfill.c.
50  * Used for tessellator
51  */
52
53 static short testedgesidef(const float v1[2], const float v2[2], const float v3[2])
54 {
55         /* is v3 to the right of v1 - v2 ? With exception: v3 == v1 || v3 == v2 */
56         double inp;
57
58         //inp = (v2[cox] - v1[cox]) * (v1[coy] - v3[coy]) + (v1[coy] - v2[coy]) * (v1[cox] - v3[cox]);
59         inp = (v2[0] - v1[0]) * (v1[1] - v3[1]) + (v1[1] - v2[1]) * (v1[0] - v3[0]);
60
61         if (inp < 0.0) {
62                 return FALSE;
63         }
64         else if (inp == 0) {
65                 if (v1[0] == v3[0] && v1[1] == v3[1]) return FALSE;
66                 if (v2[0] == v3[0] && v2[1] == v3[1]) return FALSE;
67         }
68         return TRUE;
69 }
70
71 /**
72  * \brief COMPUTE POLY NORMAL
73  *
74  * Computes the normal of a planar
75  * polygon See Graphics Gems for
76  * computing newell normal.
77  */
78 static void compute_poly_normal(float normal[3], float verts[][3], int nverts)
79 {
80         float const *v_prev = verts[nverts - 1];
81         float const *v_curr = verts[0];
82         float n[3] = {0.0f};
83         int i;
84
85         /* Newell's Method */
86         for (i = 0; i < nverts; v_prev = v_curr, v_curr = verts[++i]) {
87                 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 #compute_poly_normal but operates directly on a bmesh face.
99  */
100 static void bm_face_compute_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 #compute_poly_normal and #bm_face_compute_poly_normal
127  * but takes an array of vertex locations.
128  */
129 static void bm_face_compute_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_area_calc(BMesh *bm, 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                 compute_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_perimeter_calc(BMesh *UNUSED(bm), 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_center_bounds_calc(BMesh *UNUSED(bm), 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_center_mean_calc(BMesh *UNUSED(bm), 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 compute_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(BMesh *bm, BMEdge *e)
336 {
337         BMIter iter;
338         BMFace *f;
339         
340         f = BM_iter_new(&iter, bm, BM_FACES_OF_EDGE, e);
341         for (; f; f = BM_iter_step(&iter)) {
342                 BM_face_normal_update(bm, f);
343         }
344
345         BM_vert_normal_update(bm, e->v1);
346         BM_vert_normal_update(bm, e->v2);
347 }
348
349 /**
350  * update a vert normal (but not the faces incident on it)
351  */
352 void BM_vert_normal_update(BMesh *bm, BMVert *v)
353 {
354         /* TODO, we can normalize each edge only once, then compare with previous edge */
355
356         BMIter liter;
357         BMLoop *l;
358         float vec1[3], vec2[3], fac;
359         int len = 0;
360
361         zero_v3(v->no);
362
363         BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
364                 /* Same calculation used in BM_mesh_normals_update */
365                 sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
366                 sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
367                 normalize_v3(vec1);
368                 normalize_v3(vec2);
369
370                 fac = saacos(-dot_v3v3(vec1, vec2));
371
372                 madd_v3_v3fl(v->no, l->f->no, fac);
373
374                 len++;
375         }
376
377         if (len) {
378                 normalize_v3(v->no);
379         }
380 }
381
382 void BM_vert_normal_update_all(BMesh *bm, BMVert *v)
383 {
384         BMIter iter;
385         BMFace *f;
386
387         BM_ITER_ELEM (f, &iter, v, BM_FACES_OF_VERT) {
388                 BM_face_normal_update(bm, f);
389         }
390
391         BM_vert_normal_update(bm, v);
392 }
393
394 /**
395  * \brief BMESH UPDATE FACE NORMAL
396  *
397  * Updates the stored normal for the
398  * given face. Requires that a buffer
399  * of sufficient length to store projected
400  * coordinates for all of the face's vertices
401  * is passed in as well.
402  */
403
404 void BM_face_normal_update(BMesh *UNUSED(bm), BMFace *f)
405 {
406         BMLoop *l;
407
408         /* common cases first */
409         switch (f->len) {
410                 case 4:
411                 {
412                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
413                         const float *co2 = (l = l->next)->v->co;
414                         const float *co3 = (l = l->next)->v->co;
415                         const float *co4 = (l->next)->v->co;
416
417                         normal_quad_v3(f->no, co1, co2, co3, co4);
418                         break;
419                 }
420                 case 3:
421                 {
422                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
423                         const float *co2 = (l = l->next)->v->co;
424                         const float *co3 = (l->next)->v->co;
425
426                         normal_tri_v3(f->no, co1, co2, co3);
427                         break;
428                 }
429                 case 0:
430                 {
431                         zero_v3(f->no);
432                         break;
433                 }
434                 default:
435                 {
436                         bm_face_compute_poly_normal(f);
437                         break;
438                 }
439         }
440 }
441 /* exact same as 'bmesh_face_normal_update' but accepts vertex coords */
442 void BM_face_normal_update_vcos(BMesh *bm, BMFace *f, float no[3],
443                                 float const (*vertexCos)[3])
444 {
445         BMLoop *l;
446
447         /* must have valid index data */
448         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
449
450         /* common cases first */
451         switch (f->len) {
452                 case 4:
453                 {
454                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
455                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
456                         const float *co3 = vertexCos[BM_elem_index_get((l = l->next)->v)];
457                         const float *co4 = vertexCos[BM_elem_index_get((l->next)->v)];
458
459                         normal_quad_v3(no, co1, co2, co3, co4);
460                         break;
461                 }
462                 case 3:
463                 {
464                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
465                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
466                         const float *co3 = vertexCos[BM_elem_index_get((l->next)->v)];
467
468                         normal_tri_v3(no, co1, co2, co3);
469                         break;
470                 }
471                 case 0:
472                 {
473                         zero_v3(no);
474                         break;
475                 }
476                 default:
477                 {
478                         bm_face_compute_poly_normal_vertex_cos(f, no, vertexCos);
479                         break;
480                 }
481         }
482 }
483
484 /**
485  * \brief Face Flip Normal
486  *
487  * Reverses the winding of a face.
488  * \note This updates the calculated normal.
489  */
490 void BM_face_normal_flip(BMesh *bm, BMFace *f)
491 {
492         bmesh_loop_reverse(bm, f);
493         negate_v3(f->no);
494 }
495
496 /* detects if two line segments cross each other (intersects).
497  * note, there could be more winding cases then there needs to be. */
498 static int linecrossesf(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
499 {
500
501 #define GETMIN2_AXIS(a, b, ma, mb, axis)   \
502         {                                      \
503                 ma[axis] = MIN2(a[axis], b[axis]); \
504                 mb[axis] = MAX2(a[axis], b[axis]); \
505         } (void)0
506
507 #define GETMIN2(a, b, ma, mb)          \
508         {                                  \
509                 GETMIN2_AXIS(a, b, ma, mb, 0); \
510                 GETMIN2_AXIS(a, b, ma, mb, 1); \
511         } (void)0
512
513 #define EPS (FLT_EPSILON * 15)
514
515         int w1, w2, w3, w4, w5 /*, re */;
516         float mv1[2], mv2[2], mv3[2], mv4[2];
517         
518         /* now test winding */
519         w1 = testedgesidef(v1, v3, v2);
520         w2 = testedgesidef(v2, v4, v1);
521         w3 = !testedgesidef(v1, v2, v3);
522         w4 = testedgesidef(v3, v2, v4);
523         w5 = !testedgesidef(v3, v1, v4);
524         
525         if (w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5) {
526                 return TRUE;
527         }
528         
529         GETMIN2(v1, v2, mv1, mv2);
530         GETMIN2(v3, v4, mv3, mv4);
531         
532         /* do an interval test on the x and y axes */
533         /* first do x axis */
534         if (ABS(v1[1] - v2[1]) < EPS &&
535             ABS(v3[1] - v4[1]) < EPS &&
536             ABS(v1[1] - v3[1]) < EPS)
537         {
538                 return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
539         }
540
541         /* now do y axis */
542         if (ABS(v1[0] - v2[0]) < EPS &&
543             ABS(v3[0] - v4[0]) < EPS &&
544             ABS(v1[0] - v3[0]) < EPS)
545         {
546                 return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
547         }
548
549         return FALSE;
550
551 #undef GETMIN2_AXIS
552 #undef GETMIN2
553 #undef EPS
554
555 }
556
557 /**
558  *  BM POINT IN FACE
559  *
560  * Projects co onto face f, and returns true if it is inside
561  * the face bounds.
562  *
563  * \note this uses a best-axis projection test,
564  * instead of projecting co directly into f's orientation space,
565  * so there might be accuracy issues.
566  */
567 int BM_face_point_inside_test(BMesh *bm, BMFace *f, const float co[3])
568 {
569         int ax, ay;
570         float co2[2], cent[2] = {0.0f, 0.0f}, out[2] = {FLT_MAX * 0.5f, FLT_MAX * 0.5f};
571         BMLoop *l_iter;
572         BMLoop *l_first;
573         int crosses = 0;
574         float onepluseps = 1.0f + (float)FLT_EPSILON * 150.0f;
575         
576         if (dot_v3v3(f->no, f->no) <= FLT_EPSILON * 10)
577                 BM_face_normal_update(bm, f);
578         
579         /* find best projection of face XY, XZ or YZ: barycentric weights of
580          * the 2d projected coords are the same and faster to compute
581          *
582          * this probably isn't all that accurate, but it has the advantage of
583          * being fast (especially compared to projecting into the face orientation)
584          */
585         axis_dominant_v3(&ax, &ay, f->no);
586
587         co2[0] = co[ax];
588         co2[1] = co[ay];
589         
590         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
591         do {
592                 cent[0] += l_iter->v->co[ax];
593                 cent[1] += l_iter->v->co[ay];
594         } while ((l_iter = l_iter->next) != l_first);
595         
596         mul_v2_fl(cent, 1.0f / (float)f->len);
597         
598         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
599         do {
600                 float v1[2], v2[2];
601                 
602                 v1[0] = (l_iter->prev->v->co[ax] - cent[ax]) * onepluseps + cent[ax];
603                 v1[1] = (l_iter->prev->v->co[ay] - cent[ay]) * onepluseps + cent[ay];
604                 
605                 v2[0] = (l_iter->v->co[ax] - cent[ax]) * onepluseps + cent[ax];
606                 v2[1] = (l_iter->v->co[ay] - cent[ay]) * onepluseps + cent[ay];
607                 
608                 crosses += linecrossesf(v1, v2, co2, out) != 0;
609         } while ((l_iter = l_iter->next) != l_first);
610         
611         return crosses % 2 != 0;
612 }
613
614 static int goodline(float const (*projectverts)[3], BMFace *f,
615                     int v1i, int v2i, int v3i,
616                     int UNUSED(nvert))
617 {
618         BMLoop *l_iter;
619         BMLoop *l_first;
620         float v1[3], v2[3], v3[3], pv1[3], pv2[3];
621         int i;
622
623         copy_v3_v3(v1, projectverts[v1i]);
624         copy_v3_v3(v2, projectverts[v2i]);
625         copy_v3_v3(v3, projectverts[v3i]);
626         
627         if (testedgesidef(v1, v2, v3)) {
628                 return FALSE;
629         }
630
631         //for (i = 0; i < nvert; i++) {
632         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
633         do {
634                 i = BM_elem_index_get(l_iter->v);
635                 if (i == v1i || i == v2i || i == v3i) {
636                         continue;
637                 }
638                 
639                 copy_v3_v3(pv1, projectverts[BM_elem_index_get(l_iter->v)]);
640                 copy_v3_v3(pv2, projectverts[BM_elem_index_get(l_iter->next->v)]);
641                 
642                 //if (linecrossesf(pv1, pv2, v1, v3)) return FALSE;
643
644                 if (isect_point_tri_v2(pv1, v1, v2, v3) ||
645                     isect_point_tri_v2(pv1, v3, v2, v1))
646                 {
647                         return FALSE;
648                 }
649         } while ((l_iter = l_iter->next) != l_first);
650         return TRUE;
651 }
652
653 /**
654  * \brief Find Ear
655  *
656  * Used by tessellator to find
657  * the next triangle to 'clip off'
658  * of a polygon while tessellating.
659  *
660  * \param use_beauty Currently only applies to quads, can be extended later on.
661  */
662 static BMLoop *find_ear(BMesh *UNUSED(bm), BMFace *f, float (*verts)[3], const int nvert, const int use_beauty)
663 {
664         BMLoop *bestear = NULL;
665
666         BMLoop *l_iter;
667         BMLoop *l_first;
668
669         if (f->len == 4) {
670                 BMLoop *larr[4];
671                 int i = 0;
672
673                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
674                 do {
675                         larr[i] = l_iter;
676                         i++;
677                 } while ((l_iter = l_iter->next) != l_first);
678
679                 /* pick 0/1 based on best lenth */
680                 bestear = larr[(((len_squared_v3v3(larr[0]->v->co, larr[2]->v->co) >
681                                   len_squared_v3v3(larr[1]->v->co, larr[3]->v->co))) != use_beauty)];
682
683         }
684         else {
685                 BMVert *v1, *v2, *v3;
686
687                 /* float angle, bestangle = 180.0f; */
688                 int isear /*, i = 0 */;
689
690                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
691                 do {
692                         isear = TRUE;
693
694                         v1 = l_iter->prev->v;
695                         v2 = l_iter->v;
696                         v3 = l_iter->next->v;
697
698                         if (BM_edge_exists(v1, v3)) {
699                                 isear = FALSE;
700                         }
701                         else if (!goodline((float const (*)[3])verts, f,
702                                            BM_elem_index_get(v1), BM_elem_index_get(v2), BM_elem_index_get(v3),
703                                            nvert))
704                         {
705                                 isear = FALSE;
706                         }
707
708                         if (isear) {
709         #if 0
710                                 /* if this code comes back, it needs to be converted to radians */
711                                 angle = angle_v3v3v3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
712                                 if (!bestear || ABS(angle - 45.0f) < bestangle) {
713                                         bestear = l;
714                                         bestangle = ABS(45.0f - angle);
715                                 }
716
717                                 if (angle > 20 && angle < 90) break;
718                                 if (angle < 100 && i > 5) break;
719                                 i += 1;
720         #endif
721
722                                 bestear = l_iter;
723                                 break;
724                         }
725                 } while ((l_iter = l_iter->next) != l_first);
726         }
727
728         return bestear;
729 }
730
731 /**
732  * \brief BMESH TRIANGULATE FACE
733  *
734  * Triangulates a face using a simple 'ear clipping' algorithm that tries to
735  * favor non-skinny triangles (angles less than 90 degrees).
736  *
737  * If the triangulator has bits left over (or cannot triangulate at all)
738  * it uses a simple fan triangulation,
739  *
740  * newfaces, if non-null, must be an array of BMFace pointers,
741  * with a length equal to f->len.  it will be filled with the new
742  * triangles, and will be NULL-terminated.
743  *
744  * \note newedgeflag sets a flag layer flag, obviously not the header flag.
745  */
746 void BM_face_triangulate(BMesh *bm, BMFace *f, float (*projectverts)[3],
747                          const short newedge_oflag, const short newface_oflag, BMFace **newfaces,
748                          const short use_beauty)
749 {
750         int i, done, nvert, nf_i = 0;
751         BMLoop *newl, *nextloop;
752         BMLoop *l_iter;
753         BMLoop *l_first;
754         /* BMVert *v; */ /* UNUSED */
755
756         /* copy vertex coordinates to vertspace arra */
757         i = 0;
758         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
759         do {
760                 copy_v3_v3(projectverts[i], l_iter->v->co);
761                 BM_elem_index_set(l_iter->v, i); /* set dirty! */
762                 i++;
763         } while ((l_iter = l_iter->next) != l_first);
764
765         bm->elem_index_dirty |= BM_VERT; /* see above */
766
767         ///bmesh_face_normal_update(bm, f, f->no, projectverts);
768
769         compute_poly_normal(f->no, projectverts, f->len);
770         poly_rotate_plane(f->no, projectverts, i);
771
772         nvert = f->len;
773
774         //compute_poly_plane(projectverts, i);
775         for (i = 0; i < nvert; i++) {
776                 projectverts[i][2] = 0.0f;
777         }
778
779         done = 0;
780         while (!done && f->len > 3) {
781                 done = 1;
782                 l_iter = find_ear(bm, f, projectverts, nvert, use_beauty);
783                 if (l_iter) {
784                         done = 0;
785                         /* v = l->v; */ /* UNUSED */
786                         f = BM_face_split(bm, l_iter->f, l_iter->prev->v,
787                                           l_iter->next->v,
788                                           &newl, NULL, TRUE);
789
790                         if (UNLIKELY(!f)) {
791                                 fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
792                                 break;
793                         }
794
795                         copy_v3_v3(f->no, l_iter->f->no);
796                         BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
797                         BMO_elem_flag_enable(bm, f, newface_oflag);
798                         
799                         if (newfaces) newfaces[nf_i++] = f;
800
801 #if 0
802                         l = f->loopbase;
803                         do {
804                                 if (l->v == v) {
805                                         f->loopbase = l;
806                                         break;
807                                 }
808                                 l = l->next;
809                         } while (l != f->loopbase);
810 #endif
811
812                 }
813         }
814
815         if (f->len > 3) {
816                 l_iter = BM_FACE_FIRST_LOOP(f);
817                 while (l_iter->f->len > 3) {
818                         nextloop = l_iter->next->next;
819                         f = BM_face_split(bm, l_iter->f, l_iter->v, nextloop->v,
820                                           &newl, NULL, TRUE);
821                         if (!f) {
822                                 printf("triangle fan step of triangulator failed.\n");
823
824                                 /* NULL-terminate */
825                                 if (newfaces) newfaces[nf_i] = NULL;
826                                 return;
827                         }
828
829                         if (newfaces) newfaces[nf_i++] = f;
830                         
831                         BMO_elem_flag_enable(bm, newl->e, newedge_oflag);
832                         BMO_elem_flag_enable(bm, f, newface_oflag);
833                         l_iter = nextloop;
834                 }
835         }
836         
837         /* NULL-terminate */
838         if (newfaces) newfaces[nf_i] = NULL;
839 }
840
841 /**
842  * each pair of loops defines a new edge, a split.  this function goes
843  * through and sets pairs that are geometrically invalid to null.  a
844  * split is invalid, if it forms a concave angle or it intersects other
845  * edges in the face, or it intersects another split.  in the case of
846  * intersecting splits, only the first of the set of intersecting
847  * splits survives
848  */
849 void BM_face_legal_splits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
850 {
851         BMIter iter;
852         BMLoop *l;
853         float v1[3], v2[3], v3[3] /*, v4[3 */, no[3], mid[3], *p1, *p2, *p3, *p4;
854         float out[3] = {-234324.0f, -234324.0f, 0.0f};
855         float (*projverts)[3];
856         float (*edgeverts)[3];
857         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
858         int i, j, a = 0, clen;
859
860         BLI_array_fixedstack_declare(projverts, BM_NGON_STACK_SIZE, f->len,      "projvertsb");
861         BLI_array_fixedstack_declare(edgeverts, BM_NGON_STACK_SIZE * 2, len * 2, "edgevertsb");
862         
863         i = 0;
864         l = BM_iter_new(&iter, bm, BM_LOOPS_OF_FACE, f);
865         for ( ; l; l = BM_iter_step(&iter)) {
866                 BM_elem_index_set(l, i); /* set_loop */
867                 copy_v3_v3(projverts[i], l->v->co);
868                 i++;
869         }
870         
871         for (i = 0; i < len; i++) {
872                 copy_v3_v3(v1, loops[i][0]->v->co);
873                 copy_v3_v3(v2, loops[i][1]->v->co);
874
875                 shrink_edgef(v1, v2, fac2);
876                 
877                 copy_v3_v3(edgeverts[a], v1);
878                 a++;
879                 copy_v3_v3(edgeverts[a], v2);
880                 a++;
881         }
882         
883         compute_poly_normal(no, projverts, f->len);
884         poly_rotate_plane(no, projverts, f->len);
885         poly_rotate_plane(no, edgeverts, len * 2);
886
887         for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
888                 p1 = projverts[i];
889                 out[0] = MAX2(out[0], p1[0]) + 0.01f;
890                 out[1] = MAX2(out[1], p1[1]) + 0.01f;
891                 out[2] = 0.0f;
892                 p1[2] = 0.0f;
893
894                 //copy_v3_v3(l->v->co, p1);
895         }
896         
897         for (i = 0; i < len; i++) {
898                 edgeverts[i * 2][2] = 0.0f;
899                 edgeverts[i * 2 + 1][2] = 0.0f;
900         }
901
902         /* do convexity test */
903         for (i = 0; i < len; i++) {
904                 copy_v3_v3(v2, edgeverts[i * 2]);
905                 copy_v3_v3(v3, edgeverts[i * 2 + 1]);
906
907                 mid_v3_v3v3(mid, v2, v3);
908                 
909                 clen = 0;
910                 for (j = 0; j < f->len; j++) {
911                         p1 = projverts[j];
912                         p2 = projverts[(j + 1) % f->len];
913                         
914                         copy_v3_v3(v1, p1);
915                         copy_v3_v3(v2, p2);
916
917                         shrink_edgef(v1, v2, fac1);
918
919                         if (linecrossesf(p1, p2, mid, out)) clen++;
920                 }
921                 
922                 if (clen % 2 == 0) {
923                         loops[i][0] = NULL;
924                 }
925         }
926         
927         /* do line crossing test */
928         for (i = 0; i < f->len; i++) {
929                 p1 = projverts[i];
930                 p2 = projverts[(i + 1) % f->len];
931                 
932                 copy_v3_v3(v1, p1);
933                 copy_v3_v3(v2, p2);
934
935                 shrink_edgef(v1, v2, fac1);
936
937                 for (j = 0; j < len; j++) {
938                         if (!loops[j][0]) {
939                                 continue;
940                         }
941
942                         p3 = edgeverts[j * 2];
943                         p4 = edgeverts[j * 2 + 1];
944
945                         if (linecrossesf(v1, v2, p3, p4)) {
946                                 loops[j][0] = NULL;
947                         }
948                 }
949         }
950
951         for (i = 0; i < len; i++) {
952                 for (j = 0; j < len; j++) {
953                         if (j != i && loops[i][0] && loops[j][0]) {
954                                 p1 = edgeverts[i * 2];
955                                 p2 = edgeverts[i * 2 + 1];
956                                 p3 = edgeverts[j * 2];
957                                 p4 = edgeverts[j * 2 + 1];
958
959                                 copy_v3_v3(v1, p1);
960                                 copy_v3_v3(v2, p2);
961
962                                 shrink_edgef(v1, v2, fac1);
963
964                                 if (linecrossesf(v1, v2, p3, p4)) {
965                                         loops[i][0] = NULL;
966                                 }
967                         }
968                 }
969         }
970
971         BLI_array_fixedstack_free(projverts);
972         BLI_array_fixedstack_free(edgeverts);
973 }