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