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