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