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