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