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