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