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