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