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