part 1 of cleaning up my little array macro library to be a formal API. also removed...
[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_arithb.h"
8 #include "BLI_blenlib.h"
9 #include "BLI_array.h"
10
11 #include "MEM_guardedalloc.h"
12
13 #include "bmesh.h"
14 #include "bmesh_private.h"
15
16 /*
17  *
18  * BME POLYGON.C
19  *
20  * This file contains code for dealing
21  * with polygons (normal/area calculation,
22  * tesselation, ect)
23  *
24  * TODO:
25  *   -Add in Tesselator frontend that creates
26  *     BMTriangles from copied faces
27  *  -Add in Function that checks for and flags
28  *   degenerate faces.
29  *
30 */
31
32 /*
33  * TEST EDGE SIDE and POINT IN TRIANGLE
34  *
35  * Point in triangle tests stolen from scanfill.c.
36  * Used for tesselator
37  *
38 */
39
40 static short testedgeside(double *v1, double *v2, double *v3)
41 /* is v3 to the right of v1-v2 ? With exception: v3==v1 || v3==v2 */
42 {
43         double inp;
44
45         //inp= (v2[cox]-v1[cox])*(v1[coy]-v3[coy]) +(v1[coy]-v2[coy])*(v1[cox]-v3[cox]);
46         inp= (v2[0]-v1[0])*(v1[1]-v3[1]) +(v1[1]-v2[1])*(v1[0]-v3[0]);
47
48         if(inp<0.0) return 0;
49         else if(inp==0) {
50                 if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
51                 if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
52         }
53         return 1;
54 }
55
56 static short testedgesidef(float *v1, float *v2, float *v3)
57 /* is v3 to the right of v1-v2 ? With exception: v3==v1 || v3==v2 */
58 {
59         double inp;
60
61         //inp= (v2[cox]-v1[cox])*(v1[coy]-v3[coy]) +(v1[coy]-v2[coy])*(v1[cox]-v3[cox]);
62         inp= (v2[0]-v1[0])*(v1[1]-v3[1]) +(v1[1]-v2[1])*(v1[0]-v3[0]);
63
64         if(inp<0.0) return 0;
65         else if(inp==0) {
66                 if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
67                 if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
68         }
69         return 1;
70 }
71
72 static int point_in_triangle(double *v1, double *v2, double *v3, double *pt)
73 {
74         if(testedgeside(v1,v2,pt) && testedgeside(v2,v3,pt) && testedgeside(v3,v1,pt))
75                 return 1;
76         return 0;
77 }
78
79 /*
80  * COMPUTE POLY NORMAL
81  *
82  * Computes the normal of a planar 
83  * polygon See Graphics Gems for 
84  * computing newell normal.
85  *
86 */
87 #define FEQ(f1, f2) (ABS((double)(f1)-(double)(f2)) < 0.1)
88
89 static void compute_poly_normal(float normal[3], float (*verts)[3], int nverts)
90 {
91
92         double u[3],  v[3], w[3];/*, *w, v1[3], v2[3];*/
93         double n[3] = {0.0, 0.0, 0.0}, l, v1[3], v2[3];
94         double l2;
95         int i, s=0;
96
97         /*this fixes some weird numerical error*/
98         verts[0][0] += 0.0001f;
99         verts[0][1] += 0.0001f;
100         verts[0][2] += 0.0001f;
101
102         for(i = 0; i < nverts; i++){
103                 VECCOPY(u, verts[i]);
104                 VECCOPY(v, verts[(i+1) % nverts]);
105                 VECCOPY(w, verts[(i+2) % nverts]);
106                 
107 #if 0
108                 VECSUB(v1, w, v);
109                 VECSUB(v2, v, u);
110                 Normalize_d(v1);
111                 Normalize_d(v2);
112
113                 l = INPR(v1, v2);
114                 if (ABS(l-1.0) < 0.1)
115                         continue;
116 #endif
117                 /* newell's method
118                 
119                 so thats?:
120                 (a[1] - b[1]) * (a[2] + b[2]);
121                 a[1]*b[2] - b[1]*a[2] - b[1]*b[2] + a[1]*a[2]
122
123                 odd.  half of that is the cross product. . .what's the
124                 other half?
125
126                 also could be like a[1]*(b[2] + a[2]) - b[1]*(a[2] - b[2])
127                 */
128
129                 n[0] += (u[1] - v[1]) * (u[2] + v[2]);
130                 n[1] += (u[2] - v[2]) * (u[0] + v[0]);
131                 n[2] += (u[0] - v[0]) * (u[1] + v[1]);
132         }
133         
134         l = n[0]*n[0]+n[1]*n[1]+n[2]*n[2];
135         l = sqrt(l);
136         /*fast square root, newton/babylonian method:
137         l2 = l*0.1;
138
139         l2 = (l/l2 + l2)*0.5;
140         l2 = (l/l2 + l2)*0.5;
141         l2 = (l/l2 + l2)*0.5;
142         */
143
144         if (l == 0.0) {
145                 normal[0] = 0.0f;
146                 normal[1] = 0.0f;
147                 normal[2] = 1.0f;
148
149                 return;
150         } else l = 1.0f / l;
151
152         n[0] *= l;
153         n[1] *= l;
154         n[2] *= l;
155         
156         normal[0] = (float) n[0];
157         normal[1] = (float) n[1];
158         normal[2] = (float) n[2];
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         center[0] = center[1] = center[2] = 0.0;        
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                 VECCOPY(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         VECADD(center, min, max);
246         VECMUL(center, 0.5f);
247 }
248
249 /*
250  * COMPUTE POLY PLANE
251  *
252  * Projects a set polygon's vertices to 
253  * a plane defined by the average
254  * of its edges cross products
255  *
256 */
257
258 void compute_poly_plane(float (*verts)[3], int nverts)
259 {
260         
261         float avgc[3], norm[3], temp[3], mag, avgn[3];
262         float *v1, *v2, *v3;
263         int i;
264         
265         if(nverts < 3) 
266                 return;
267
268         avgn[0] = avgn[1] = avgn[2] = 0.0;
269         avgc[0] = avgc[1] = avgc[2] = 0.0;
270
271         for(i = 0; i < nverts; i++){
272                 v1 = verts[i];
273                 v2 = verts[(i+1) % nverts];
274                 v3 = verts[(i+2) % nverts];
275                 CalcNormFloat(v1, v2, v3, norm);        
276         
277                 avgn[0] += norm[0];
278                 avgn[1] += norm[1];
279                 avgn[2] += norm[2];
280         }
281
282         /*what was this bit for?*/
283         if(avgn[0] == 0.0 && avgn[1] == 0.0 && avgn[2] == 0.0){
284                 avgn[0] = 0.0;
285                 avgn[1] = 0.0;
286                 avgn[2] = 1.0;
287         } else {
288                 avgn[0] /= nverts;
289                 avgn[1] /= nverts;
290                 avgn[2] /= nverts;
291                 Normalize(avgn);
292         }
293         
294         for(i = 0; i < nverts; i++){
295                 v1 = verts[i];
296                 VECCOPY(temp, v1);
297                 mag = 0.0;
298                 mag += (temp[0] * avgn[0]);
299                 mag += (temp[1] * avgn[1]);
300                 mag += (temp[2] * avgn[2]);
301                 
302                 temp[0] = (avgn[0] * mag);
303                 temp[1] = (avgn[1] * mag);
304                 temp[2] = (avgn[2] * mag);
305
306                 VecSubf(v1, v1, temp);
307         }       
308 }
309
310 /*
311   BM LEGAL EDGES
312
313   takes in a face and a list of edges, and sets to NULL any edge in
314   the list that bridges a concave region of the face or intersects
315   any of the faces's edges.
316 */
317 static void shrink_edged(double *v1, double *v2, double fac)
318 {
319         double mid[3];
320
321         VECADD(mid, v1, v2);
322         VECMUL(mid, 0.5);
323
324         VECSUB(v1, v1, mid);
325         VECSUB(v2, v2, mid);
326
327         VECMUL(v1, fac);
328         VECMUL(v2, fac);
329
330         VECADD(v1, v1, mid);
331         VECADD(v2, v2, mid);
332 }
333
334 static void shrink_edgef(float *v1, float *v2, float fac)
335 {
336         float mid[3];
337
338         VECADD(mid, v1, v2);
339         VECMUL(mid, 0.5);
340
341         VECSUB(v1, v1, mid);
342         VECSUB(v2, v2, mid);
343
344         VECMUL(v1, fac);
345         VECMUL(v2, fac);
346
347         VECADD(v1, v1, mid);
348         VECADD(v2, v2, mid);
349 }
350
351 /*
352  * POLY ROTATE PLANE
353  *
354  * Rotates a polygon so that it's
355  * normal is pointing towards the mesh Z axis
356  *
357 */
358
359 void poly_rotate_plane(float normal[3], float (*verts)[3], int nverts)
360 {
361
362         float up[3] = {0.0f,0.0f,1.0f}, axis[3], q[4];
363         float mat[3][3];
364         double angle;
365         int i;
366
367         Crossf(axis, up, normal);
368         axis[0] *= -1;
369         axis[1] *= -1;
370         axis[2] *= -1;
371
372         angle = saacos(normal[0]*up[0]+normal[1]*up[1] + normal[2]*up[2]);
373
374         if (angle == 0.0f) return;
375
376         AxisAngleToQuatd(q, axis, angle);
377         QuatToMat3(q, mat);
378
379         for(i = 0;  i < nverts;  i++)
380                 Mat3MulVecfl(mat, verts[i]);
381 }
382
383 /*
384  * BMESH UPDATE FACE NORMAL
385  *
386  * Updates the stored normal for the
387  * given face. Requires that a buffer
388  * of sufficient length to store projected
389  * coordinates for all of the face's vertices
390  * is passed in as well.
391  *
392 */
393
394 void BM_Face_UpdateNormal(BMesh *bm, BMFace *f)
395 {
396         float projverts[200][3];
397         float (*proj)[3] = f->len < 200 ? projverts : MEM_mallocN(sizeof(float)*f->len*3, "projvertsn");
398         BMLoop *l = f->loopbase;
399         int i=0;
400
401         if (f->len < 3) return;
402         
403         do {
404                 VECCOPY(proj[i], l->v->co);
405                 i += 1;
406         } while (l != f->loopbase);
407
408         bmesh_update_face_normal(bm, f, proj);
409
410         if (projverts != proj) MEM_freeN(proj);
411 }
412
413 void BM_Edge_UpdateNormals(BMesh *bm, BMEdge *e)
414 {
415         BMIter iter;
416         BMFace *f;
417         
418         f = BMIter_New(&iter, bm, BM_FACES_OF_EDGE, e);
419         for (; f; f=BMIter_Step(&iter)) {
420                 BM_Face_UpdateNormal(bm, f);
421         }
422
423         BM_Vert_UpdateNormal(bm, e->v1);
424         BM_Vert_UpdateNormal(bm, e->v2);
425 }
426
427 void BM_Vert_UpdateNormal(BMesh *bm, BMVert *v)
428 {
429         BMIter iter;
430         BMFace *f;
431         int len=0;
432
433         v->no[0] = v->no[1] = v->no[2] = 0.0f;
434
435         f = BMIter_New(&iter, bm, BM_FACES_OF_VERT, v);
436         for (; f; f=BMIter_Step(&iter), len++) {
437                 VecAddf(v->no, f->no, v->no);
438         }
439
440         if (!len) return;
441
442         VecMulf(v->no, 1.0f/(int)len);
443 }
444
445 void bmesh_update_face_normal(BMesh *bm, BMFace *f, float (*projectverts)[3])
446 {
447         BMLoop *l;
448         int i;
449
450         if(f->len > 4) {
451                 i = 0;
452                 l = f->loopbase;
453                 do{
454                         VECCOPY(projectverts[i], l->v->co);
455                         l = (BMLoop*)(l->head.next);
456                         i += 1;
457                 }while(l!=f->loopbase);
458
459                 compute_poly_normal(f->no, projectverts, f->len);       
460         }
461         else if(f->len == 3){
462                 BMVert *v1, *v2, *v3;
463                 v1 = f->loopbase->v;
464                 v2 = ((BMLoop*)(f->loopbase->head.next))->v;
465                 v3 = ((BMLoop*)(f->loopbase->head.next->next))->v;
466                 CalcNormFloat(v1->co, v2->co, v3->co, f->no);
467         }
468         else if(f->len == 4){
469                 BMVert *v1, *v2, *v3, *v4;
470                 v1 = f->loopbase->v;
471                 v2 = ((BMLoop*)(f->loopbase->head.next))->v;
472                 v3 = ((BMLoop*)(f->loopbase->head.next->next))->v;
473                 v4 = ((BMLoop*)(f->loopbase->head.prev))->v;
474                 CalcNormFloat4(v1->co, v2->co, v3->co, v4->co, f->no);
475         }
476         else{ /*horrible, two sided face!*/
477                 f->no[0] = 0.0;
478                 f->no[1] = 0.0;
479                 f->no[2] = 1.0;
480         }
481
482 }
483
484
485 /*
486  * BMESH FLIP NORMAL
487  * 
488  *  Reverses the winding of a face.
489  *  Note that this updates the calculated 
490  *  normal.
491 */
492 void BM_flip_normal(BMesh *bm, BMFace *f)
493 {       
494         bmesh_loop_reverse(bm, f);
495         BM_Face_UpdateNormal(bm, f);
496 }
497
498 /* detects if two line segments cross each other (intersects).
499    note, there could be more winding cases then there needs to be. */
500 int linecrosses(double *v1, double *v2, double *v3, double *v4)
501 {
502         int w1, w2, w3, w4, w5;
503         
504         /*w1 = winding(v1, v3, v4);
505         w2 = winding(v2, v3, v4);
506         w3 = winding(v3, v1, v2);
507         w4 = winding(v4, v1, v2);
508         
509         return (w1 == w2) && (w3 == w4);*/
510
511         w1 = testedgeside(v1, v3, v2);
512         w2 = testedgeside(v2, v4, v1);
513         w3 = !testedgeside(v1, v2, v3);
514         w4 = testedgeside(v3, v2, v4);
515         w5 = !testedgeside(v3, v1, v4);
516         return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5;
517 }
518
519 /* detects if two line segments cross each other (intersects).
520    note, there could be more winding cases then there needs to be. */
521 int linecrossesf(float *v1, float *v2, float *v3, float *v4)
522 {
523         int w1, w2, w3, w4, w5, ret;
524         
525 /*         int test1_a, test1_a, test2_a, test2_a;
526
527    test1_a = check_tri_clock_dir(l1p1, l1p2, l2p1);
528    test1_b = check_tri_clock_dir(l1p1, l1p2, l2p2);
529    if (test1_a != test1_b)
530    {
531       test2_a = check_tri_clock_dir(l2p1, l2p2, l1p1);
532       test2_b = check_tri_clock_dir(l2p1, l2p2, l1p2);
533       if (test2_a != test2_b)
534       {
535          return 1;
536       }
537    }*/
538         /*w1 = testedgesidef(v1, v2, v3);
539         w2 = testedgesidef(v1, v2, v4);
540         if(w1 != w2) {
541                 w3 = testedgesidef(v3, v4, v1);
542                 w4 = testedgesidef(v3, v4, v2);
543                 if (w3 != w4) return 1;
544         }
545                 
546         return 0;*/
547
548         /*w1 = testedgesidef(v1, v3, v4);
549         w2 = testedgesidef(v2, v3, v4);
550         w3 = testedgesidef(v3, v1, v2);
551         w4 = testedgesidef(v4, v1, v2);
552         
553         return (w1 == w2) && (w2 == w3) && (w3 == w4);*/
554
555         /*do an interval test on the x and y axes*/
556         /*first do x axis*/
557         #define T 0.01
558         if (ABS(v1[1]-v2[1]) < T && ABS(v3[1]-v4[1]) < T &&
559             ABS(v1[1]-v3[1]) < T) {
560                 if (v3[0] >= v1[0] && v3[0] <= v2[0])
561                         return 1;
562                 if (v4[0] >= v1[0] && v4[0] <= v2[0])
563                         return 1;
564                 if (v3[0] <= v1[0] && v4[0] >= v2[0])
565                         return 1;
566         }
567
568         /*now do y axis*/
569         if (ABS(v1[0]-v2[0]) < T && ABS(v3[0]-v4[0]) < T &&
570             ABS(v1[0]-v3[0]) < T) {
571                 if (v3[1] >= v1[1] && v3[1] <= v2[1])
572                         return 1;
573                 if (v4[1] >= v1[1] && v4[1] <= v2[1])
574                         return 1;
575                 if (v3[1] <= v1[1] && v4[1] >= v2[1])
576                         return 1;
577         }
578
579         /*now test winding*/
580         w1 = testedgesidef(v1, v3, v2);
581         w2 = testedgesidef(v2, v4, v1);
582         w3 = !testedgesidef(v1, v2, v3);
583         w4 = testedgesidef(v3, v2, v4);
584         w5 = !testedgesidef(v3, v1, v4);
585
586         return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5; 
587 }
588
589 int goodline(float (*projectverts)[3], BMFace *f, int v1i,
590              int v2i, int v3i, int nvert) {
591         BMLoop *l = f->loopbase;
592         double v1[3], v2[3], v3[3], pv1[3], pv2[3];
593         int i;
594
595         VECCOPY(v1, projectverts[v1i]);
596         VECCOPY(v2, projectverts[v2i]);
597         VECCOPY(v3, projectverts[v3i]);
598         
599         if (testedgeside(v1, v2, v3)) return 0;
600         
601         //for (i=0; i<nvert; i++) {
602         do {
603                 i = l->v->head.eflag2;
604                 if (i == v1i || i == v2i || i == v3i) {
605                         l = (BMLoop*)l->head.next;
606                         continue;
607                 }
608                 
609                 VECCOPY(pv1, projectverts[l->v->head.eflag2]);
610                 VECCOPY(pv2, projectverts[((BMLoop*)l->head.next)->v->head.eflag2]);
611                 
612                 //if (linecrosses(pv1, pv2, v1, v3)) return 0;
613                 if (point_in_triangle(v1, v2, v3, pv1)) return 0;
614                 if (point_in_triangle(v3, v2, v1, pv1)) return 0;
615
616                 l = (BMLoop*)l->head.next;
617         } while (l != f->loopbase);
618         return 1;
619 }
620 /*
621  * FIND EAR
622  *
623  * Used by tesselator to find
624  * the next triangle to 'clip off'
625  * of a polygon while tesselating.
626  *
627 */
628
629 static BMLoop *find_ear(BMesh *bm, BMFace *f, float (*verts)[3], 
630                         int nvert)
631 {
632         BMVert *v1, *v2, *v3;
633         BMLoop *bestear = NULL, *l;
634         float angle, bestangle = 180.0f;
635         int isear, i=0;
636         
637         l = f->loopbase;
638         do {
639                 isear = 1;
640                 
641                 v1 = ((BMLoop*)(l->head.prev))->v;
642                 v2 = l->v;
643                 v3 = ((BMLoop*)(l->head.next))->v;
644
645                 if (BM_Edge_Exist(v1, v3)) isear = 0;
646
647                 if (isear && !goodline(verts, f, v1->head.eflag2, v2->head.eflag2,
648                                        v3->head.eflag2, nvert))
649                         isear = 0;
650                 
651                 if(isear) {
652                         /*angle = VecAngle3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
653                         if(!bestear || ABS(angle-45.0f) < bestangle) {
654                                 bestear = l;
655                                 bestangle = ABS(45.0f-angle);
656                         }
657                         
658                         if (angle > 20 && angle < 90) break;
659                         if (angle < 100 && i > 5) break;
660                         i += 1;*/
661                         bestear = l;
662                         break;
663                 }
664                 l = (BMLoop*)(l->head.next);
665         }
666         while(l != f->loopbase);
667
668         return bestear;
669 }
670
671 /*
672  * BMESH TRIANGULATE FACE
673  *
674  * Triangulates a face using a 
675  * simple 'ear clipping' algorithm
676  * that tries to favor non-skinny
677  * triangles (angles less than 
678  * 90 degrees). If the triangulator
679  * has bits left over (or cannot
680  * triangulate at all) it uses a
681  * simple fan triangulation
682  *
683  * newfaces, if non-null, must be an array of BMFace pointers,
684  * with a length equal to f->len.  it will be filled with the new
685  * triangles, and will be NULL-terminated.
686 */
687 void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], 
688                          int newedgeflag, int newfaceflag, BMFace **newfaces)
689 {
690         int i, done, nvert, nf_i = 0;
691         BMLoop *l, *newl, *nextloop;
692         BMVert *v;
693
694         /*copy vertex coordinates to vertspace array*/
695         i = 0;
696         l = f->loopbase;
697         do{
698                 VECCOPY(projectverts[i], l->v->co);
699                 l->v->head.eflag2 = i; /*warning, abuse! never duplicate in tools code! never you hear?*/ /*actually, get rid of this completely, use a new structure for this....*/
700                 i++;
701                 l = (BMLoop*)(l->head.next);
702         }while(l != f->loopbase);
703         
704         ///bmesh_update_face_normal(bm, f, projectverts);
705
706         compute_poly_normal(f->no, projectverts, f->len);
707         poly_rotate_plane(f->no, projectverts, i);
708
709         nvert = f->len;
710
711         //compute_poly_plane(projectverts, i);
712         for (i=0; i<nvert; i++) {
713                 projectverts[i][2] = 0.0f;
714         }
715
716         done = 0;
717         while(!done && f->len > 3){
718                 done = 1;
719                 l = find_ear(bm, f, projectverts, nvert);
720                 if(l) {
721                         done = 0;
722                         v = l->v;
723                         f = BM_Split_Face(bm, l->f, ((BMLoop*)(l->head.prev))->v, 
724                                           ((BMLoop*)(l->head.next))->v, 
725                                           &newl, NULL);
726                         VECCOPY(f->no, l->f->no);
727
728                         if (!f) {
729                                 printf("yeek! triangulator failed to split face!\n");
730                                 break;
731                         }
732
733                         BMO_SetFlag(bm, newl->e, newedgeflag);
734                         BMO_SetFlag(bm, f, newfaceflag);
735                         
736                         if (newfaces) newfaces[nf_i++] = f;
737
738                         /*l = f->loopbase;
739                         do {
740                                 if (l->v == v) {
741                                         f->loopbase = l;
742                                         break;
743                                 }
744                                 l = l->head.next;
745                         } while (l != f->loopbase);*/
746                 }
747         }
748
749         if (f->len > 3){
750                 l = f->loopbase;
751                 while (l->f->len > 3){
752                         nextloop = ((BMLoop*)(l->head.next->next));
753                         f = BM_Split_Face(bm, l->f, l->v, nextloop->v, 
754                                           &newl, NULL);
755                         if (!f) {
756                                 printf("triangle fan step of triangulator failed.\n");
757
758                                 /*NULL-terminate*/
759                                 if (newfaces) newfaces[nf_i] = NULL;
760                                 return;
761                         }
762
763                         if (newfaces) newfaces[nf_i++] = f;
764                         
765                         BMO_SetFlag(bm, newl->e, newedgeflag);
766                         BMO_SetFlag(bm, f, newfaceflag);
767                         l = nextloop;
768                 }
769         }
770         
771         /*NULL-terminate*/
772         if (newfaces) newfaces[nf_i] = NULL;
773 }
774
775 /*each pair of loops defines a new edge, a split.  this function goes
776   through and sets pairs that are geometrically invalid to null.  a
777   split is invalid, if it forms a concave angle or it intersects other
778   edges in the face, or it intersects another split.  in the case of
779   intersecting splits, only the first of the set of intersecting
780   splits survives.*/
781 void BM_LegalSplits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
782 {
783         BMIter iter;
784         BMLoop *l;
785         float v1[3], v2[3], v3[3], v4[3], no[3], mid[3], *p1, *p2, *p3, *p4;
786         float out[3] = {-234324.0f, -234324.0f, 0.0f};
787         float projectverts[100][3];
788         float edgevertsstack[200][3];
789         float (*projverts)[3] = projectverts;
790         float (*edgeverts)[3] = edgevertsstack;
791         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
792         int i, j, a=0, clen;
793
794         if (f->len > 100) projverts = MEM_mallocN(sizeof(float)*3*f->len, "projvertsb");
795         if (len > 100) edgeverts = MEM_mallocN(sizeof(float)*3*2*len, "edgevertsb");
796         
797         i = 0;
798         l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
799         for (; l; l=BMIter_Step(&iter)) {
800                 l->head.eflag2 = i;
801                 VECCOPY(projverts[i], l->v->co);
802                 i++;
803         }
804         
805         for (i=0; i<len; i++) {
806                 VECCOPY(v1, loops[i][0]->v->co);
807                 VECCOPY(v2, loops[i][1]->v->co);
808
809                 shrink_edgef(v1, v2, fac2);
810                 
811                 VECCOPY(edgeverts[a], v1);
812                 a++;
813                 VECCOPY(edgeverts[a], v2);
814                 a++;
815         }
816         
817         compute_poly_normal(no, projverts, f->len);
818         poly_rotate_plane(no, projverts, f->len);
819         poly_rotate_plane(no, edgeverts, len*2);
820         
821         l = f->loopbase;
822         for (i=0; i<f->len; i++) {
823                 p1 = projverts[i];
824                 out[0] = MAX2(out[0], p1[0]) + 0.01f;
825                 out[1] = MAX2(out[1], p1[1]) + 0.01f;
826                 out[2] = 0.0f;
827                 p1[2] = 0.0f;
828
829                 //VECCOPY(l->v->co, p1);
830
831                 l = (BMLoop*) l->head.next;
832         }
833         
834         for (i=0; i<len; i++) {
835                 edgeverts[i*2][2] = 0.0f;
836                 edgeverts[i*2+1][2] = 0.0f;
837         }
838
839         /*do convexity test*/
840         for (i=0; i<len; i++) {
841                 VECCOPY(v2, edgeverts[i*2]);
842                 VECCOPY(v3, edgeverts[i*2+1]);
843
844                 VecAddf(mid, v2, v3);
845                 VecMulf(mid, 0.5f);
846                 
847                 clen = 0;
848                 for (j=0; j<f->len; j++) {
849                         p1 = projverts[j];
850                         p2 = projverts[(j+1)%f->len];
851                         
852                         VECCOPY(v1, p1);
853                         VECCOPY(v2, p2);
854
855                         shrink_edgef(v1, v2, fac1);
856
857                         if (linecrossesf(p1, p2, mid, out)) clen++;
858                 }
859                 
860                 if (clen%2 == 0) {
861                         loops[i][0] = NULL;
862                 }
863         }
864         
865         /*do line crossing tests*/
866         for (i=0; i<f->len; i++) {
867                 p1 = projverts[i];
868                 p2 = projverts[(i+1)%f->len];
869                 
870                 VECCOPY(v1, p1);
871                 VECCOPY(v2, p2);
872
873                 shrink_edgef(v1, v2, fac1);
874
875                 for (j=0; j<len; j++) {
876                         if (!loops[j][0]) continue;
877
878                         p3 = edgeverts[j*2];
879                         p4 = edgeverts[j*2+1];
880
881                         if (linecrossesf(v1, v2, p3, p4))
882                         {
883                                 loops[j][0] = NULL;
884                         }
885                 }
886         }
887
888         for (i=0; i<len; i++) {
889                 for (j=0; j<len; j++) {
890                         if (j == i) continue;
891                         if (!loops[i][0]) continue;
892                         if (!loops[j][0]) continue;
893
894                         p1 = edgeverts[i*2];
895                         p2 = edgeverts[i*2+1];
896                         p3 = edgeverts[j*2];
897                         p4 = edgeverts[j*2+1];
898
899                         VECCOPY(v1, p1);
900                         VECCOPY(v2, p2);
901
902                         shrink_edgef(v1, v2, fac1);
903
904                         if (linecrossesf(v1, v2, p3, p4)) {
905                                 loops[i][0]=NULL;
906                         }
907                 }
908         }
909         
910         if (projverts != projectverts) MEM_freeN(projverts);
911         if (edgeverts != edgevertsstack) MEM_freeN(edgeverts);
912 }