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