Merging r58464 through r58474 from trunk into soc-2013-depsgraph_mt
[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
31 #include "DNA_listBase.h"
32
33 #include "MEM_guardedalloc.h"
34
35 #include "BLI_math.h"
36 #include "BLI_array.h"
37 #include "BLI_scanfill.h"
38 #include "BLI_listbase.h"
39
40 #include "bmesh.h"
41
42 #include "intern/bmesh_private.h"
43
44 /**
45  * \brief TEST EDGE SIDE and POINT IN TRIANGLE
46  *
47  * Point in triangle tests stolen from scanfill.c.
48  * Used for tessellator
49  */
50
51 static bool testedgesidef(const float v1[2], const float v2[2], const float v3[2])
52 {
53         /* is v3 to the right of v1 - v2 ? With exception: v3 == v1 || v3 == v2 */
54         double inp;
55
56         //inp = (v2[cox] - v1[cox]) * (v1[coy] - v3[coy]) + (v1[coy] - v2[coy]) * (v1[cox] - v3[cox]);
57         inp = (v2[0] - v1[0]) * (v1[1] - v3[1]) + (v1[1] - v2[1]) * (v1[0] - v3[0]);
58
59         if (inp < 0.0) {
60                 return false;
61         }
62         else if (inp == 0) {
63                 if (v1[0] == v3[0] && v1[1] == v3[1]) return false;
64                 if (v2[0] == v3[0] && v2[1] == v3[1]) return false;
65         }
66         return true;
67 }
68
69 /**
70  * \brief COMPUTE POLY NORMAL
71  *
72  * Computes the normal of a planar
73  * polygon See Graphics Gems for
74  * computing newell normal.
75  */
76 static void calc_poly_normal(float normal[3], float verts[][3], int nverts)
77 {
78         float const *v_prev = verts[nverts - 1];
79         float const *v_curr = verts[0];
80         float n[3] = {0.0f};
81         int i;
82
83         /* Newell's Method */
84         for (i = 0; i < nverts; v_prev = v_curr, v_curr = verts[++i]) {
85                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
86         }
87
88         if (UNLIKELY(normalize_v3_v3(normal, n) == 0.0f)) {
89                 normal[2] = 1.0f; /* other axis set to 0.0 */
90         }
91 }
92
93 /**
94  * \brief COMPUTE POLY NORMAL (BMFace)
95  *
96  * Same as #calc_poly_normal but operates directly on a bmesh face.
97  */
98 static void bm_face_calc_poly_normal(BMFace *f, float n[3])
99 {
100         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
101         BMLoop *l_iter  = l_first;
102         float const *v_prev = l_first->prev->v->co;
103         float const *v_curr = l_first->v->co;
104
105         zero_v3(n);
106
107         /* Newell's Method */
108         do {
109                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
110
111                 l_iter = l_iter->next;
112                 v_prev = v_curr;
113                 v_curr = l_iter->v->co;
114
115         } while (l_iter != l_first);
116
117         if (UNLIKELY(normalize_v3(n) == 0.0f)) {
118                 n[2] = 1.0f;
119         }
120 }
121
122 /**
123  * \brief COMPUTE POLY NORMAL (BMFace)
124  *
125  * Same as #calc_poly_normal and #bm_face_calc_poly_normal
126  * but takes an array of vertex locations.
127  */
128 static void bm_face_calc_poly_normal_vertex_cos(BMFace *f, float n[3],
129                                                 float const (*vertexCos)[3])
130 {
131         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
132         BMLoop *l_iter  = l_first;
133         float const *v_prev = vertexCos[BM_elem_index_get(l_first->prev->v)];
134         float const *v_curr = vertexCos[BM_elem_index_get(l_first->v)];
135
136         zero_v3(n);
137
138         /* Newell's Method */
139         do {
140                 add_newell_cross_v3_v3v3(n, v_prev, v_curr);
141
142                 l_iter = l_iter->next;
143                 v_prev = v_curr;
144                 v_curr = vertexCos[BM_elem_index_get(l_iter->v)];
145         } while (l_iter != l_first);
146
147         if (UNLIKELY(normalize_v3(n) == 0.0f)) {
148                 n[2] = 1.0f; /* other axis set to 0.0 */
149         }
150 }
151
152 /**
153  * For tools that insist on using triangles, ideally we would cache this data.
154  *
155  * \param r_loops  Store face loop pointers, (f->len)
156  * \param r_index  Store triangle triples, indicies into \a r_loops,  ((f->len - 2) * 3)
157  */
158 int BM_face_calc_tessellation(BMFace *f, BMLoop **r_loops, int (*_r_index)[3])
159 {
160         int *r_index = (int *)_r_index;
161         BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
162         BMLoop *l_iter;
163         int totfilltri;
164
165         if (f->len == 3) {
166                 *r_loops++ = (l_iter = l_first);
167                 *r_loops++ = (l_iter = l_iter->next);
168                 *r_loops++ = (         l_iter->next);
169
170                 r_index[0] = 0;
171                 r_index[1] = 1;
172                 r_index[2] = 2;
173                 totfilltri = 1;
174         }
175         else if (f->len == 4) {
176                 *r_loops++ = (l_iter = l_first);
177                 *r_loops++ = (l_iter = l_iter->next);
178                 *r_loops++ = (l_iter = l_iter->next);
179                 *r_loops++ = (         l_iter->next);
180
181                 r_index[0] = 0;
182                 r_index[1] = 1;
183                 r_index[2] = 2;
184
185                 r_index[3] = 0;
186                 r_index[4] = 2;
187                 r_index[5] = 3;
188                 totfilltri = 2;
189         }
190         else {
191                 int j;
192
193                 ScanFillContext sf_ctx;
194                 ScanFillVert *sf_vert, *sf_vert_last = NULL, *sf_vert_first = NULL;
195                 /* ScanFillEdge *e; */ /* UNUSED */
196                 ScanFillFace *sf_tri;
197
198                 BLI_scanfill_begin(&sf_ctx);
199
200                 j = 0;
201                 l_iter = l_first;
202                 do {
203                         sf_vert = BLI_scanfill_vert_add(&sf_ctx, l_iter->v->co);
204                         sf_vert->tmp.p = l_iter;
205
206                         if (sf_vert_last) {
207                                 /* e = */ BLI_scanfill_edge_add(&sf_ctx, sf_vert_last, sf_vert);
208                         }
209
210                         sf_vert_last = sf_vert;
211                         if (sf_vert_first == NULL) {
212                                 sf_vert_first = sf_vert;
213                         }
214
215                         r_loops[j] = l_iter;
216
217                         /* mark order */
218                         BM_elem_index_set(l_iter, j++); /* set_loop */
219
220                 } while ((l_iter = l_iter->next) != l_first);
221
222                 /* complete the loop */
223                 BLI_scanfill_edge_add(&sf_ctx, sf_vert_first, sf_vert);
224
225                 totfilltri = BLI_scanfill_calc_ex(&sf_ctx, 0, f->no);
226                 BLI_assert(totfilltri <= f->len - 2);
227                 BLI_assert(totfilltri == BLI_countlist(&sf_ctx.fillfacebase));
228
229                 for (sf_tri = sf_ctx.fillfacebase.first; sf_tri; sf_tri = sf_tri->next) {
230                         int i1 = BM_elem_index_get((BMLoop *)sf_tri->v1->tmp.p);
231                         int i2 = BM_elem_index_get((BMLoop *)sf_tri->v2->tmp.p);
232                         int i3 = BM_elem_index_get((BMLoop *)sf_tri->v3->tmp.p);
233
234                         if (i1 > i2) { SWAP(int, i1, i2); }
235                         if (i2 > i3) { SWAP(int, i2, i3); }
236                         if (i1 > i2) { SWAP(int, i1, i2); }
237
238                         *r_index++ = i1;
239                         *r_index++ = i2;
240                         *r_index++ = i3;
241                 }
242
243                 BLI_scanfill_end(&sf_ctx);
244         }
245
246         return totfilltri;
247 }
248
249 /**
250  * get the area of the face
251  */
252 float BM_face_calc_area(BMFace *f)
253 {
254         BMLoop *l;
255         BMIter iter;
256         float (*verts)[3] = BLI_array_alloca(verts, f->len);
257         float area;
258         int i;
259
260         BM_ITER_ELEM_INDEX (l, &iter, f, BM_LOOPS_OF_FACE, i) {
261                 copy_v3_v3(verts[i], l->v->co);
262         }
263
264         if (f->len == 3) {
265                 area = area_tri_v3(verts[0], verts[1], verts[2]);
266         }
267         else if (f->len == 4) {
268                 area = area_quad_v3(verts[0], verts[1], verts[2], verts[3]);
269         }
270         else {
271                 float normal[3];
272                 calc_poly_normal(normal, verts, f->len);
273                 area = area_poly_v3(f->len, verts, normal);
274         }
275
276         return area;
277 }
278
279 /**
280  * compute the perimeter of an ngon
281  */
282 float BM_face_calc_perimeter(BMFace *f)
283 {
284         BMLoop *l_iter, *l_first;
285         float perimeter = 0.0f;
286
287         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
288         do {
289                 perimeter += len_v3v3(l_iter->v->co, l_iter->next->v->co);
290         } while ((l_iter = l_iter->next) != l_first);
291
292         return perimeter;
293 }
294
295 /**
296  * Compute a meaningful direction along the face (use for manipulator axis).
297  * \note result isnt normalized.
298  */
299 void BM_face_calc_plane(BMFace *f, float r_plane[3])
300 {
301         if (f->len == 3) {
302                 BMVert *verts[3];
303                 float lens[3];
304                 float difs[3];
305                 int  order[3] = {0, 1, 2};
306
307                 BM_face_as_array_vert_tri(f, verts);
308
309                 lens[0] = len_v3v3(verts[0]->co, verts[1]->co);
310                 lens[1] = len_v3v3(verts[1]->co, verts[2]->co);
311                 lens[2] = len_v3v3(verts[2]->co, verts[0]->co);
312
313                 /* find the shortest or the longest loop */
314                 difs[0] = fabsf(lens[1] - lens[2]);
315                 difs[1] = fabsf(lens[2] - lens[0]);
316                 difs[2] = fabsf(lens[0] - lens[1]);
317
318                 axis_sort_v3(difs, order);
319                 sub_v3_v3v3(r_plane, verts[order[0]]->co, verts[(order[0] + 1) % 3]->co);
320         }
321         else if (f->len == 4) {
322                 BMVert *verts[4];
323                 float vec[3], vec_a[3], vec_b[3];
324
325                 // BM_iter_as_array(NULL, BM_VERTS_OF_FACE, efa, (void **)verts, 4);
326                 BM_face_as_array_vert_quad(f, verts);
327
328                 sub_v3_v3v3(vec_a, verts[3]->co, verts[2]->co);
329                 sub_v3_v3v3(vec_b, verts[0]->co, verts[1]->co);
330                 add_v3_v3v3(r_plane, vec_a, vec_b);
331
332                 sub_v3_v3v3(vec_a, verts[0]->co, verts[3]->co);
333                 sub_v3_v3v3(vec_b, verts[1]->co, verts[2]->co);
334                 add_v3_v3v3(vec, vec_a, vec_b);
335                 /* use the biggest edge length */
336                 if (dot_v3v3(r_plane, r_plane) < dot_v3v3(vec, vec)) {
337                         copy_v3_v3(r_plane, vec);
338                 }
339         }
340         else {
341                 BMLoop *l_long  = BM_face_find_longest_loop(f);
342
343                 sub_v3_v3v3(r_plane, l_long->v->co, l_long->next->v->co);
344         }
345
346         normalize_v3(r_plane);
347 }
348
349 /**
350  * computes center of face in 3d.  uses center of bounding box.
351  */
352 void BM_face_calc_center_bounds(BMFace *f, float r_cent[3])
353 {
354         BMLoop *l_iter;
355         BMLoop *l_first;
356         float min[3], max[3];
357
358         INIT_MINMAX(min, max);
359
360         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
361         do {
362                 minmax_v3v3_v3(min, max, l_iter->v->co);
363         } while ((l_iter = l_iter->next) != l_first);
364
365         mid_v3_v3v3(r_cent, min, max);
366 }
367
368 /**
369  * computes the center of a face, using the mean average
370  */
371 void BM_face_calc_center_mean(BMFace *f, float r_cent[3])
372 {
373         BMLoop *l_iter;
374         BMLoop *l_first;
375
376         zero_v3(r_cent);
377
378         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
379         do {
380                 add_v3_v3(r_cent, l_iter->v->co);
381         } while ((l_iter = l_iter->next) != l_first);
382
383         if (f->len)
384                 mul_v3_fl(r_cent, 1.0f / (float) f->len);
385 }
386
387 /**
388  * computes the center of a face, using the mean average
389  * weighted by edge length
390  */
391 void BM_face_calc_center_mean_weighted(BMFace *f, float r_cent[3])
392 {
393         BMLoop *l_iter;
394         BMLoop *l_first;
395         float totw = 0.0f;
396         float w_prev;
397
398         zero_v3(r_cent);
399
400
401         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
402         w_prev = BM_edge_calc_length(l_iter->prev->e);
403         do {
404                 const float w_curr = BM_edge_calc_length(l_iter->e);
405                 const float w = (w_curr + w_prev);
406                 madd_v3_v3fl(r_cent, l_iter->v->co, w);
407                 totw += w;
408                 w_prev = w_curr;
409         } while ((l_iter = l_iter->next) != l_first);
410
411         if (totw != 0.0f)
412                 mul_v3_fl(r_cent, 1.0f / (float) totw);
413 }
414
415 /**
416  * COMPUTE POLY PLANE
417  *
418  * Projects a set polygon's vertices to
419  * a plane defined by the average
420  * of its edges cross products
421  */
422 void calc_poly_plane(float (*verts)[3], const int nverts)
423 {
424         
425         float avgc[3], norm[3], mag, avgn[3];
426         float *v1, *v2, *v3;
427         int i;
428         
429         if (nverts < 3)
430                 return;
431
432         zero_v3(avgn);
433         zero_v3(avgc);
434
435         for (i = 0; i < nverts; i++) {
436                 v1 = verts[i];
437                 v2 = verts[(i + 1) % nverts];
438                 v3 = verts[(i + 2) % nverts];
439                 normal_tri_v3(norm, v1, v2, v3);
440
441                 add_v3_v3(avgn, norm);
442         }
443
444         if (UNLIKELY(normalize_v3(avgn) == 0.0f)) {
445                 avgn[2] = 1.0f;
446         }
447         
448         for (i = 0; i < nverts; i++) {
449                 v1 = verts[i];
450                 mag = dot_v3v3(v1, avgn);
451                 madd_v3_v3fl(v1, avgn, -mag);
452         }
453 }
454
455 /**
456  * \brief BM LEGAL EDGES
457  *
458  * takes in a face and a list of edges, and sets to NULL any edge in
459  * the list that bridges a concave region of the face or intersects
460  * any of the faces's edges.
461  */
462 static void scale_edge_v2f(float v1[2], float v2[2], const float fac)
463 {
464         float mid[2];
465
466         mid_v2_v2v2(mid, v1, v2);
467
468         sub_v2_v2v2(v1, v1, mid);
469         sub_v2_v2v2(v2, v2, mid);
470
471         mul_v2_fl(v1, fac);
472         mul_v2_fl(v2, fac);
473
474         add_v2_v2v2(v1, v1, mid);
475         add_v2_v2v2(v2, v2, mid);
476 }
477
478 /**
479  * \brief POLY ROTATE PLANE
480  *
481  * Rotates a polygon so that it's
482  * normal is pointing towards the mesh Z axis
483  */
484 void poly_rotate_plane(const float normal[3], float (*verts)[3], const int nverts)
485 {
486         float mat[3][3];
487
488         if (axis_dominant_v3_to_m3(mat, normal)) {
489                 int i;
490                 for (i = 0; i < nverts; i++) {
491                         mul_m3_v3(mat, verts[i]);
492                 }
493         }
494 }
495
496 /**
497  * updates face and vertex normals incident on an edge
498  */
499 void BM_edge_normals_update(BMEdge *e)
500 {
501         BMIter iter;
502         BMFace *f;
503         
504         BM_ITER_ELEM (f, &iter, e, BM_FACES_OF_EDGE) {
505                 BM_face_normal_update(f);
506         }
507
508         BM_vert_normal_update(e->v1);
509         BM_vert_normal_update(e->v2);
510 }
511
512 /**
513  * update a vert normal (but not the faces incident on it)
514  */
515 void BM_vert_normal_update(BMVert *v)
516 {
517         /* TODO, we can normalize each edge only once, then compare with previous edge */
518
519         BMIter liter;
520         BMLoop *l;
521         float vec1[3], vec2[3], fac;
522         int len = 0;
523
524         zero_v3(v->no);
525
526         BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
527                 /* Same calculation used in BM_mesh_normals_update */
528                 sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
529                 sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
530                 normalize_v3(vec1);
531                 normalize_v3(vec2);
532
533                 fac = saacos(-dot_v3v3(vec1, vec2));
534
535                 madd_v3_v3fl(v->no, l->f->no, fac);
536
537                 len++;
538         }
539
540         if (len) {
541                 normalize_v3(v->no);
542         }
543 }
544
545 void BM_vert_normal_update_all(BMVert *v)
546 {
547         BMIter iter;
548         BMFace *f;
549
550         BM_ITER_ELEM (f, &iter, v, BM_FACES_OF_VERT) {
551                 BM_face_normal_update(f);
552         }
553
554         BM_vert_normal_update(v);
555 }
556
557 /**
558  * \brief BMESH UPDATE FACE NORMAL
559  *
560  * Updates the stored normal for the
561  * given face. Requires that a buffer
562  * of sufficient length to store projected
563  * coordinates for all of the face's vertices
564  * is passed in as well.
565  */
566
567 void BM_face_calc_normal(BMFace *f, float r_no[3])
568 {
569         BMLoop *l;
570
571         /* common cases first */
572         switch (f->len) {
573                 case 4:
574                 {
575                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
576                         const float *co2 = (l = l->next)->v->co;
577                         const float *co3 = (l = l->next)->v->co;
578                         const float *co4 = (l->next)->v->co;
579
580                         normal_quad_v3(r_no, co1, co2, co3, co4);
581                         break;
582                 }
583                 case 3:
584                 {
585                         const float *co1 = (l = BM_FACE_FIRST_LOOP(f))->v->co;
586                         const float *co2 = (l = l->next)->v->co;
587                         const float *co3 = (l->next)->v->co;
588
589                         normal_tri_v3(r_no, co1, co2, co3);
590                         break;
591                 }
592                 default:
593                 {
594                         bm_face_calc_poly_normal(f, r_no);
595                         break;
596                 }
597         }
598 }
599 void BM_face_normal_update(BMFace *f)
600 {
601         BM_face_calc_normal(f, f->no);
602 }
603
604 /* exact same as 'bmesh_face_normal_update' but accepts vertex coords */
605 void BM_face_normal_update_vcos(BMesh *bm, BMFace *f, float no[3],
606                                 float const (*vertexCos)[3])
607 {
608         BMLoop *l;
609
610         /* must have valid index data */
611         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
612         (void)bm;
613
614         /* common cases first */
615         switch (f->len) {
616                 case 4:
617                 {
618                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
619                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
620                         const float *co3 = vertexCos[BM_elem_index_get((l = l->next)->v)];
621                         const float *co4 = vertexCos[BM_elem_index_get((l->next)->v)];
622
623                         normal_quad_v3(no, co1, co2, co3, co4);
624                         break;
625                 }
626                 case 3:
627                 {
628                         const float *co1 = vertexCos[BM_elem_index_get((l = BM_FACE_FIRST_LOOP(f))->v)];
629                         const float *co2 = vertexCos[BM_elem_index_get((l = l->next)->v)];
630                         const float *co3 = vertexCos[BM_elem_index_get((l->next)->v)];
631
632                         normal_tri_v3(no, co1, co2, co3);
633                         break;
634                 }
635                 case 0:
636                 {
637                         zero_v3(no);
638                         break;
639                 }
640                 default:
641                 {
642                         bm_face_calc_poly_normal_vertex_cos(f, no, vertexCos);
643                         break;
644                 }
645         }
646 }
647
648 /**
649  * \brief Face Flip Normal
650  *
651  * Reverses the winding of a face.
652  * \note This updates the calculated normal.
653  */
654 void BM_face_normal_flip(BMesh *bm, BMFace *f)
655 {
656         bmesh_loop_reverse(bm, f);
657         negate_v3(f->no);
658 }
659
660 /* detects if two line segments cross each other (intersects).
661  * note, there could be more winding cases then there needs to be. */
662 static bool line_crosses_v2f(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
663 {
664
665 #define GETMIN2_AXIS(a, b, ma, mb, axis)   \
666         {                                      \
667                 ma[axis] = min_ff(a[axis], b[axis]); \
668                 mb[axis] = max_ff(a[axis], b[axis]); \
669         } (void)0
670
671 #define GETMIN2(a, b, ma, mb)          \
672         {                                  \
673                 GETMIN2_AXIS(a, b, ma, mb, 0); \
674                 GETMIN2_AXIS(a, b, ma, mb, 1); \
675         } (void)0
676
677 #define EPS (FLT_EPSILON * 15)
678
679         int w1, w2, w3, w4, w5 /*, re */;
680         float mv1[2], mv2[2], mv3[2], mv4[2];
681         
682         /* now test winding */
683         w1 = testedgesidef(v1, v3, v2);
684         w2 = testedgesidef(v2, v4, v1);
685         w3 = !testedgesidef(v1, v2, v3);
686         w4 = testedgesidef(v3, v2, v4);
687         w5 = !testedgesidef(v3, v1, v4);
688         
689         if (w1 == w2 && w2 == w3 && w3 == w4 && w4 == w5) {
690                 return true;
691         }
692         
693         GETMIN2(v1, v2, mv1, mv2);
694         GETMIN2(v3, v4, mv3, mv4);
695         
696         /* do an interval test on the x and y axes */
697         /* first do x axis */
698         if (fabsf(v1[1] - v2[1]) < EPS &&
699             fabsf(v3[1] - v4[1]) < EPS &&
700             fabsf(v1[1] - v3[1]) < EPS)
701         {
702                 return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
703         }
704
705         /* now do y axis */
706         if (fabsf(v1[0] - v2[0]) < EPS &&
707             fabsf(v3[0] - v4[0]) < EPS &&
708             fabsf(v1[0] - v3[0]) < EPS)
709         {
710                 return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
711         }
712
713         return false;
714
715 #undef GETMIN2_AXIS
716 #undef GETMIN2
717 #undef EPS
718
719 }
720
721 /**
722  *  BM POINT IN FACE
723  *
724  * Projects co onto face f, and returns true if it is inside
725  * the face bounds.
726  *
727  * \note this uses a best-axis projection test,
728  * instead of projecting co directly into f's orientation space,
729  * so there might be accuracy issues.
730  */
731 bool BM_face_point_inside_test(BMFace *f, const float co[3])
732 {
733         int ax, ay;
734         float co2[2], cent[2] = {0.0f, 0.0f}, out[2] = {FLT_MAX * 0.5f, FLT_MAX * 0.5f};
735         BMLoop *l_iter;
736         BMLoop *l_first;
737         int crosses = 0;
738         float onepluseps = 1.0f + (float)FLT_EPSILON * 150.0f;
739         
740         if (dot_v3v3(f->no, f->no) <= FLT_EPSILON * 10)
741                 BM_face_normal_update(f);
742         
743         /* find best projection of face XY, XZ or YZ: barycentric weights of
744          * the 2d projected coords are the same and faster to compute
745          *
746          * this probably isn't all that accurate, but it has the advantage of
747          * being fast (especially compared to projecting into the face orientation)
748          */
749         axis_dominant_v3(&ax, &ay, f->no);
750
751         co2[0] = co[ax];
752         co2[1] = co[ay];
753         
754         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
755         do {
756                 cent[0] += l_iter->v->co[ax];
757                 cent[1] += l_iter->v->co[ay];
758         } while ((l_iter = l_iter->next) != l_first);
759         
760         mul_v2_fl(cent, 1.0f / (float)f->len);
761         
762         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
763         do {
764                 float v1[2], v2[2];
765                 
766                 v1[0] = (l_iter->prev->v->co[ax] - cent[0]) * onepluseps + cent[0];
767                 v1[1] = (l_iter->prev->v->co[ay] - cent[1]) * onepluseps + cent[1];
768                 
769                 v2[0] = (l_iter->v->co[ax] - cent[0]) * onepluseps + cent[0];
770                 v2[1] = (l_iter->v->co[ay] - cent[1]) * onepluseps + cent[1];
771                 
772                 crosses += line_crosses_v2f(v1, v2, co2, out) != 0;
773         } while ((l_iter = l_iter->next) != l_first);
774         
775         return crosses % 2 != 0;
776 }
777
778 static bool bm_face_goodline(float const (*projectverts)[2], BMFace *f, int v1i, int v2i, int v3i)
779 {
780         BMLoop *l_iter;
781         BMLoop *l_first;
782
783         float pv1[2];
784         const float *v1 = projectverts[v1i];
785         const float *v2 = projectverts[v2i];
786         const float *v3 = projectverts[v3i];
787         int i;
788
789         /* v3 must be on the left side of [v1, v2] line, else we know [v1, v3] is outside of f! */
790         if (testedgesidef(v1, v2, v3)) {
791                 return false;
792         }
793
794         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
795         do {
796                 i = BM_elem_index_get(l_iter->v);
797                 copy_v2_v2(pv1, projectverts[i]);
798
799                 if (ELEM3(i, v1i, v2i, v3i)) {
800 #if 0
801                         printf("%d in (%d, %d, %d) tri (from indices!), continuing\n", i, v1i, v2i, v3i);
802 #endif
803                         continue;
804                 }
805
806                 if (isect_point_tri_v2(pv1, v1, v2, v3) || isect_point_tri_v2(pv1, v3, v2, v1)) {
807 #if 0
808                         if (isect_point_tri_v2(pv1, v1, v2, v3))
809                                 printf("%d in (%d, %d, %d)\n", v3i, i, v1i, v2i);
810                         else
811                                 printf("%d in (%d, %d, %d)\n", v1i, i, v3i, v2i);
812 #endif
813                         return false;
814                 }
815         } while ((l_iter = l_iter->next) != l_first);
816         return true;
817 }
818
819 /**
820  * \brief Find Ear
821  *
822  * Used by tessellator to find the next triangle to 'clip off' of a polygon while tessellating.
823  *
824  * \param f The face to search.
825  * \param projectverts an array of face vert coords.
826  * \param use_beauty Currently only applies to quads, can be extended later on.
827  * \param abscoss Must be allocated by caller, and at least f->len length
828  *        (allow to avoid allocating a new one for each tri!).
829  */
830 static BMLoop *poly_find_ear(BMFace *f, float (*projectverts)[2], const bool use_beauty, float *abscoss)
831 {
832         BMLoop *bestear = NULL;
833
834         BMLoop *l_iter;
835         BMLoop *l_first;
836
837         const float cos_threshold = 0.9f;
838         const float bias = 1.0f + 1e-6f;
839
840         BLI_assert(len_squared_v3(f->no) > FLT_EPSILON);
841
842         if (f->len == 4) {
843                 BMLoop *larr[4];
844                 int i = 0, i4;
845                 float cos1, cos2;
846                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
847                 do {
848                         larr[i] = l_iter;
849                         i++;
850                 } while ((l_iter = l_iter->next) != l_first);
851
852                 /* pick 0/1 based on best lenth */
853                 /* XXX Can't only rely on such test, also must check we do not get (too much) degenerated triangles!!! */
854                 i = (((len_squared_v3v3(larr[0]->v->co, larr[2]->v->co) >
855                      len_squared_v3v3(larr[1]->v->co, larr[3]->v->co) * bias)) != use_beauty);
856                 i4 = (i + 3) % 4;
857                 /* Check produced tris aren't too flat/narrow...
858                  * Probably not the best test, but is quite efficient and should at least avoid null-area faces! */
859                 cos1 = fabsf(cos_v3v3v3(larr[i4]->v->co, larr[i]->v->co, larr[i + 1]->v->co));
860                 cos2 = fabsf(cos_v3v3v3(larr[i4]->v->co, larr[i + 2]->v->co, larr[i + 1]->v->co));
861 #if 0
862                 printf("%d, (%f, %f), (%f, %f)\n", i, cos1, cos2,
863                        fabsf(cos_v3v3v3(larr[i]->v->co, larr[i4]->v->co, larr[i + 2]->v->co)),
864                        fabsf(cos_v3v3v3(larr[i]->v->co, larr[i + 1]->v->co, larr[i + 2]->v->co)));
865 #endif
866                 if (cos1 < cos2)
867                         cos1 = cos2;
868                 if (cos1 > cos_threshold) {
869                         if (cos1 > fabsf(cos_v3v3v3(larr[i]->v->co, larr[i4]->v->co, larr[i + 2]->v->co)) &&
870                             cos1 > fabsf(cos_v3v3v3(larr[i]->v->co, larr[i + 1]->v->co, larr[i + 2]->v->co)))
871                         {
872                                 i = !i;
873                         }
874                 }
875                 /* Last check we do not get overlapping triangles
876                  * (as much as possible, there are some cases with no good solution!) */
877                 i4 = (i + 3) % 4;
878                 if (!bm_face_goodline((float const (*)[2])projectverts, f, BM_elem_index_get(larr[i4]->v),
879                                       BM_elem_index_get(larr[i]->v), BM_elem_index_get(larr[i + 1]->v)))
880                 {
881                         i = !i;
882                 }
883 /*              printf("%d\n", i);*/
884                 bestear = larr[i];
885
886         }
887         else {
888                 /* float angle, bestangle = 180.0f; */
889                 float cos, bestcos = 1.0f;
890                 int i, j, len;
891
892                 /* Compute cos of all corners! */
893                 i = 0;
894                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
895                 len = l_iter->f->len;
896                 do {
897                         const BMVert *v1 = l_iter->prev->v;
898                         const BMVert *v2 = l_iter->v;
899                         const BMVert *v3 = l_iter->next->v;
900
901                         abscoss[i] = fabsf(cos_v3v3v3(v1->co, v2->co, v3->co));
902 /*                      printf("tcoss: %f\n", *tcoss);*/
903                         i++;
904                 } while ((l_iter = l_iter->next) != l_first);
905
906                 i = 0;
907                 l_iter = l_first;
908                 do {
909                         const BMVert *v1 = l_iter->prev->v;
910                         const BMVert *v2 = l_iter->v;
911                         const BMVert *v3 = l_iter->next->v;
912
913                         if (bm_face_goodline((float const (*)[2])projectverts, f,
914                                              BM_elem_index_get(v1), BM_elem_index_get(v2), BM_elem_index_get(v3)))
915                         {
916                                 /* Compute highest cos (i.e. narrowest angle) of this tri. */
917                                 cos = max_fff(abscoss[i],
918                                               fabsf(cos_v3v3v3(v2->co, v3->co, v1->co)),
919                                               fabsf(cos_v3v3v3(v3->co, v1->co, v2->co)));
920
921                                 /* Compare to prev best (i.e. lowest) cos. */
922                                 if (cos < bestcos) {
923                                         /* We must check this tri would not leave a (too much) degenerated remaining face! */
924                                         /* For now just assume if the average of cos of all
925                                          * "remaining face"'s corners is below a given threshold, it's OK. */
926                                         float avgcos = fabsf(cos_v3v3v3(v1->co, v3->co, l_iter->next->next->v->co));
927                                         const int i_limit = (i - 1 + len) % len;
928                                         avgcos += fabsf(cos_v3v3v3(l_iter->prev->prev->v->co, v1->co, v3->co));
929                                         j = (i + 2) % len;
930                                         do {
931                                                 avgcos += abscoss[j];
932                                         } while ((j = (j + 1) % len) != i_limit);
933                                         avgcos /= len - 1;
934
935                                         /* We need a best ear in any case... */
936                                         if (avgcos < cos_threshold || (!bestear && avgcos < 1.0f)) {
937                                                 /* OKI, keep this ear (corner...) as a potential best one! */
938                                                 bestear = l_iter;
939                                                 bestcos = cos;
940                                         }
941 #if 0
942                                         else
943                                                 printf("Had a nice tri (higest cos of %f, current bestcos is %f), "
944                                                        "but average cos of all \"remaining face\"'s corners is too high (%f)!\n",
945                                                        cos, bestcos, avgcos);
946 #endif
947                                 }
948                         }
949                         i++;
950                 } while ((l_iter = l_iter->next) != l_first);
951         }
952
953         return bestear;
954 }
955
956 /**
957  * \brief BMESH TRIANGULATE FACE
958  *
959  * Currently repeatedly find the best triangle (i.e. the most "open" one), provided it does not
960  * produces a "remaining" face with too much wide/narrow angles
961  * (using cos (i.e. dot product of normalized vectors) of angles).
962  *
963  * \param r_faces_new if non-null, must be an array of BMFace pointers,
964  * with a length equal to (f->len - 2). It will be filled with the new
965  * triangles.
966  *
967  * \note use_tag tags new flags and edges.
968  */
969 void BM_face_triangulate(BMesh *bm, BMFace *f,
970                          BMFace **r_faces_new,
971                          const bool use_beauty, const bool use_tag)
972 {
973         const float f_len_orig = f->len;
974         int i, nf_i = 0;
975         BMLoop *l_new;
976         BMLoop *l_iter;
977         BMLoop *l_first;
978         /* BM_face_triangulate: temp absolute cosines of face corners */
979         float (*projectverts)[2] = BLI_array_alloca(projectverts, f_len_orig);
980         float *abscoss = BLI_array_alloca(abscoss, f_len_orig);
981         float mat[3][3];
982
983         axis_dominant_v3_to_m3(mat, f->no);
984
985         /* copy vertex coordinates to vertspace area */
986         i = 0;
987         l_iter = l_first = BM_FACE_FIRST_LOOP(f);
988         do {
989                 mul_v2_m3v3(projectverts[i], mat, l_iter->v->co);
990                 BM_elem_index_set(l_iter->v, i++); /* set dirty! */
991         } while ((l_iter = l_iter->next) != l_first);
992
993         bm->elem_index_dirty |= BM_VERT; /* see above */
994
995         while (f->len > 3) {
996                 l_iter = poly_find_ear(f, projectverts, use_beauty, abscoss);
997
998                 /* force triangulation - if we can't find an ear the face is degenerate */
999                 if (l_iter == NULL) {
1000                         l_iter = BM_FACE_FIRST_LOOP(f);
1001                 }
1002
1003 /*              printf("Subdividing face...\n");*/
1004                 f = BM_face_split(bm, l_iter->f, l_iter->prev->v, l_iter->next->v, &l_new, NULL, true);
1005
1006                 if (UNLIKELY(!f)) {
1007                         fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
1008                         break;
1009                 }
1010
1011                 copy_v3_v3(f->no, l_iter->f->no);
1012
1013                 if (use_tag) {
1014                         BM_elem_flag_enable(l_new->e, BM_ELEM_TAG);
1015                         BM_elem_flag_enable(f, BM_ELEM_TAG);
1016                 }
1017
1018                 if (r_faces_new) {
1019                         r_faces_new[nf_i++] = f;
1020                 }
1021         }
1022
1023         BLI_assert(f->len == 3);
1024 }
1025
1026 /**
1027  * each pair of loops defines a new edge, a split.  this function goes
1028  * through and sets pairs that are geometrically invalid to null.  a
1029  * split is invalid, if it forms a concave angle or it intersects other
1030  * edges in the face, or it intersects another split.  in the case of
1031  * intersecting splits, only the first of the set of intersecting
1032  * splits survives
1033  */
1034 void BM_face_legal_splits(BMFace *f, BMLoop *(*loops)[2], int len)
1035 {
1036         const int len2 = len * 2;
1037         BMLoop *l;
1038         float v1[2], v2[2], v3[2] /*, v4[3 */, no[3], mid[2], *p1, *p2, *p3, *p4;
1039         float out[2] = {-FLT_MAX, -FLT_MAX};
1040         float axis_mat[3][3];
1041         float (*projverts)[2] = BLI_array_alloca(projverts, f->len);
1042         float (*edgeverts)[2] = BLI_array_alloca(edgeverts, len2);
1043         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
1044         int i, j, a = 0, clen;
1045
1046         /* TODO, the face normal may already be correct */
1047         BM_face_calc_normal(f, no);
1048
1049         axis_dominant_v3_to_m3(axis_mat, no);
1050
1051         for (i = 0, l = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l = l->next) {
1052                 BM_elem_index_set(l, i); /* set_loop */
1053                 mul_v2_m3v3(projverts[i], axis_mat, l->v->co);
1054
1055                 out[0] = max_ff(out[0], projverts[i][0]);
1056                 out[1] = max_ff(out[1], projverts[i][1]);
1057         }
1058         
1059         /* ensure we are well outside the face bounds (value is arbitrary) */
1060         add_v2_fl(out, 1.0f);
1061
1062         for (i = 0; i < len; i++) {
1063                 copy_v2_v2(edgeverts[a + 0], projverts[BM_elem_index_get(loops[i][0])]);
1064                 copy_v2_v2(edgeverts[a + 1], projverts[BM_elem_index_get(loops[i][1])]);
1065                 scale_edge_v2f(edgeverts[a + 0], edgeverts[a + 1], fac2);
1066                 a += 2;
1067         }
1068
1069         /* do convexity test */
1070         for (i = 0; i < len; i++) {
1071                 copy_v2_v2(v2, edgeverts[i * 2 + 0]);
1072                 copy_v2_v2(v3, edgeverts[i * 2 + 1]);
1073
1074                 mid_v2_v2v2(mid, v2, v3);
1075                 
1076                 clen = 0;
1077                 for (j = 0; j < f->len; j++) {
1078                         p1 = projverts[j];
1079                         p2 = projverts[(j + 1) % f->len];
1080                         
1081 #if 0
1082                         copy_v2_v2(v1, p1);
1083                         copy_v2_v2(v2, p2);
1084
1085                         scale_edge_v2f(v1, v2, fac1);
1086                         if (line_crosses_v2f(v1, v2, mid, out)) {
1087                                 clen++;
1088                         }
1089 #else
1090                         if (line_crosses_v2f(p1, p2, mid, out)) {
1091                                 clen++;
1092                         }
1093 #endif
1094                 }
1095
1096                 if (clen % 2 == 0) {
1097                         loops[i][0] = NULL;
1098                 }
1099         }
1100
1101         /* do line crossing tests */
1102         for (i = 0; i < f->len; i++) {
1103                 p1 = projverts[i];
1104                 p2 = projverts[(i + 1) % f->len];
1105                 
1106                 copy_v2_v2(v1, p1);
1107                 copy_v2_v2(v2, p2);
1108
1109                 scale_edge_v2f(v1, v2, fac1);
1110
1111                 for (j = 0; j < len; j++) {
1112                         if (!loops[j][0]) {
1113                                 continue;
1114                         }
1115
1116                         p3 = edgeverts[j * 2];
1117                         p4 = edgeverts[j * 2 + 1];
1118
1119                         if (line_crosses_v2f(v1, v2, p3, p4)) {
1120                                 loops[j][0] = NULL;
1121                         }
1122                 }
1123         }
1124
1125         for (i = 0; i < len; i++) {
1126                 for (j = 0; j < len; j++) {
1127                         if (j != i && loops[i][0] && loops[j][0]) {
1128                                 p1 = edgeverts[i * 2];
1129                                 p2 = edgeverts[i * 2 + 1];
1130                                 p3 = edgeverts[j * 2];
1131                                 p4 = edgeverts[j * 2 + 1];
1132
1133                                 copy_v2_v2(v1, p1);
1134                                 copy_v2_v2(v2, p2);
1135
1136                                 scale_edge_v2f(v1, v2, fac1);
1137
1138                                 if (line_crosses_v2f(v1, v2, p3, p4)) {
1139                                         loops[i][0] = NULL;
1140                                 }
1141                         }
1142                 }
1143         }
1144 }
1145
1146
1147 /**
1148  * Small utility functions for fast access
1149  *
1150  * faster alternative to:
1151  *  BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void **)v, 3);
1152  */
1153 void BM_face_as_array_vert_tri(BMFace *f, BMVert *r_verts[3])
1154 {
1155         BMLoop *l = BM_FACE_FIRST_LOOP(f);
1156
1157         BLI_assert(f->len == 3);
1158
1159         r_verts[0] = l->v; l = l->next;
1160         r_verts[1] = l->v; l = l->next;
1161         r_verts[2] = l->v;
1162 }
1163
1164 /**
1165  * faster alternative to:
1166  *  BM_iter_as_array(bm, BM_VERTS_OF_FACE, f, (void **)v, 4);
1167  */
1168 void BM_face_as_array_vert_quad(BMFace *f, BMVert *r_verts[4])
1169 {
1170         BMLoop *l = BM_FACE_FIRST_LOOP(f);
1171
1172         BLI_assert(f->len == 4);
1173
1174         r_verts[0] = l->v; l = l->next;
1175         r_verts[1] = l->v; l = l->next;
1176         r_verts[2] = l->v; l = l->next;
1177         r_verts[3] = l->v;
1178 }