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