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