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