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