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