svn merge ^/trunk/blender -r42333:42361
[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 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         axis_dominant_v3(&ax, &ay, f->no);
628
629         co2[0] = co[ax];
630         co2[1] = co[ay];
631         co2[2] = 0;
632         
633         l = bm_firstfaceloop(f);
634         do {
635                 cent[0] += l->v->co[ax];
636                 cent[1] += l->v->co[ay];
637                 l = l->next;
638         } while (l != bm_firstfaceloop(f));
639         
640         mul_v2_fl(cent, 1.0f/(float)f->len);
641         
642         l = bm_firstfaceloop(f);
643         do {
644                 float v1[3], v2[3];
645                 
646                 v1[0] = (l->prev->v->co[ax] - cent[ax])*eps + cent[ax];
647                 v1[1] = (l->prev->v->co[ay] - cent[ay])*eps + cent[ay];
648                 v1[2] = 0.0f;
649                 
650                 v2[0] = (l->v->co[ax] - cent[ax])*eps + cent[ax];
651                 v2[1] = (l->v->co[ay] - cent[ay])*eps + cent[ay];
652                 v2[2] = 0.0f;
653                 
654                 crosses += linecrossesf(v1, v2, co2, out) != 0;
655                 
656                 l = l->next;
657         } while (l != bm_firstfaceloop(f));
658         
659         return crosses%2 != 0;
660 }
661
662 static int goodline(float (*projectverts)[3], BMFace *f, int v1i,
663                     int v2i, int v3i, int UNUSED(nvert))
664 {
665         BMLoop *l = bm_firstfaceloop(f);
666         double v1[3], v2[3], v3[3], pv1[3], pv2[3];
667         int i;
668
669         VECCOPY(v1, projectverts[v1i]);
670         VECCOPY(v2, projectverts[v2i]);
671         VECCOPY(v3, projectverts[v3i]);
672         
673         if (testedgeside(v1, v2, v3)) return 0;
674         
675         //for (i=0; i<nvert; i++) {
676         do {
677                 i = BM_GetIndex(l->v);
678                 if (i == v1i || i == v2i || i == v3i) {
679                         l = l->next;
680                         continue;
681                 }
682                 
683                 VECCOPY(pv1, projectverts[BM_GetIndex(l->v)]);
684                 VECCOPY(pv2, projectverts[BM_GetIndex(l->next->v)]);
685                 
686                 //if (linecrosses(pv1, pv2, v1, v3)) return 0;
687                 if (point_in_triangle(v1, v2, v3, pv1)) return 0;
688                 if (point_in_triangle(v3, v2, v1, pv1)) return 0;
689
690                 l = l->next;
691         } while (l != bm_firstfaceloop(f));
692         return 1;
693 }
694 /*
695  * FIND EAR
696  *
697  * Used by tesselator to find
698  * the next triangle to 'clip off'
699  * of a polygon while tesselating.
700  *
701 */
702
703 static BMLoop *find_ear(BMesh *UNUSED(bm), BMFace *f, float (*verts)[3],
704                         int nvert)
705 {
706         BMVert *v1, *v2, *v3;
707         BMLoop *bestear = NULL, *l;
708         /*float angle, bestangle = 180.0f;*/
709         int isear /*, i=0*/;
710         
711         l = bm_firstfaceloop(f);
712         do {
713                 isear = 1;
714                 
715                 v1 = l->prev->v;
716                 v2 = l->v;
717                 v3 = l->next->v;
718
719                 if (BM_Edge_Exist(v1, v3)) isear = 0;
720
721                 if (isear && !goodline(verts, f, BM_GetIndex(v1), BM_GetIndex(v2),
722                                        BM_GetIndex(v3), nvert))
723                         isear = 0;
724                 
725                 if(isear) {
726                         /*angle = angle_v3v3v3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
727                         if(!bestear || ABS(angle-45.0f) < bestangle) {
728                                 bestear = l;
729                                 bestangle = ABS(45.0f-angle);
730                         }
731                         
732                         if (angle > 20 && angle < 90) break;
733                         if (angle < 100 && i > 5) break;
734                         i += 1;*/
735                         bestear = l;
736                         break;
737                 }
738                 l = l->next;
739         }
740         while(l != bm_firstfaceloop(f));
741
742         return bestear;
743 }
744
745 /*
746  * BMESH TRIANGULATE FACE
747  *
748  * Triangulates a face using a 
749  * simple 'ear clipping' algorithm
750  * that tries to favor non-skinny
751  * triangles (angles less than 
752  * 90 degrees). If the triangulator
753  * has bits left over (or cannot
754  * triangulate at all) it uses a
755  * simple fan triangulation
756  *
757  * newfaces, if non-null, must be an array of BMFace pointers,
758  * with a length equal to f->len.  it will be filled with the new
759  * triangles, and will be NULL-terminated.
760 */
761 void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], 
762                          int newedgeflag, int newfaceflag, BMFace **newfaces)
763 {
764         int i, done, nvert, nf_i = 0;
765         BMLoop *l, *newl, *nextloop;
766         /* BMVert *v; */ /* UNUSED */
767
768         /*copy vertex coordinates to vertspace array*/
769         i = 0;
770         l = bm_firstfaceloop(f);
771         do{
772                 copy_v3_v3(projectverts[i], l->v->co);
773                 BM_SetIndex(l->v, i); /* set dirty! */
774                 i++;
775                 l = l->next;
776         }while(l != bm_firstfaceloop(f));
777
778         bm->elem_index_dirty |= BM_VERT; /* see above */
779
780         ///bmesh_update_face_normal(bm, f, projectverts);
781
782         compute_poly_normal(f->no, projectverts, f->len);
783         poly_rotate_plane(f->no, projectverts, i);
784
785         nvert = f->len;
786
787         //compute_poly_plane(projectverts, i);
788         for (i=0; i<nvert; i++) {
789                 projectverts[i][2] = 0.0f;
790         }
791
792         done = 0;
793         while(!done && f->len > 3){
794                 done = 1;
795                 l = find_ear(bm, f, projectverts, nvert);
796                 if(l) {
797                         done = 0;
798                         /* v = l->v; */ /* UNUSED */
799                         f = BM_Split_Face(bm, l->f, l->prev->v,
800                                           l->next->v,
801                                           &newl, NULL);
802                         copy_v3_v3(f->no, l->f->no);
803
804                         if (!f) {
805                                 fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
806                                 break;
807                         }
808
809                         BMO_SetFlag(bm, newl->e, newedgeflag);
810                         BMO_SetFlag(bm, f, newfaceflag);
811                         
812                         if (newfaces) newfaces[nf_i++] = f;
813
814                         /*l = f->loopbase;
815                         do {
816                                 if (l->v == v) {
817                                         f->loopbase = l;
818                                         break;
819                                 }
820                                 l = l->next;
821                         } while (l != f->loopbase);*/
822                 }
823         }
824
825         if (f->len > 3){
826                 l = bm_firstfaceloop(f);
827                 while (l->f->len > 3){
828                         nextloop = l->next->next;
829                         f = BM_Split_Face(bm, l->f, l->v, nextloop->v, 
830                                           &newl, NULL);
831                         if (!f) {
832                                 printf("triangle fan step of triangulator failed.\n");
833
834                                 /*NULL-terminate*/
835                                 if (newfaces) newfaces[nf_i] = NULL;
836                                 return;
837                         }
838
839                         if (newfaces) newfaces[nf_i++] = f;
840                         
841                         BMO_SetFlag(bm, newl->e, newedgeflag);
842                         BMO_SetFlag(bm, f, newfaceflag);
843                         l = nextloop;
844                 }
845         }
846         
847         /*NULL-terminate*/
848         if (newfaces) newfaces[nf_i] = NULL;
849 }
850
851 /*each pair of loops defines a new edge, a split.  this function goes
852   through and sets pairs that are geometrically invalid to null.  a
853   split is invalid, if it forms a concave angle or it intersects other
854   edges in the face, or it intersects another split.  in the case of
855   intersecting splits, only the first of the set of intersecting
856   splits survives.*/
857 void BM_LegalSplits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
858 {
859         BMIter iter;
860         BMLoop *l;
861         float v1[3], v2[3], v3[3]/*, v4[3]*/, no[3], mid[3], *p1, *p2, *p3, *p4;
862         float out[3] = {-234324.0f, -234324.0f, 0.0f};
863         float (*projverts)[3];
864         float (*edgeverts)[3];
865         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
866         int i, j, a=0, clen;
867
868         BLI_array_fixedstack_declare(projverts, BM_NGON_STACK_SIZE, f->len,         "projvertsb");
869         BLI_array_fixedstack_declare(edgeverts, BM_NGON_STACK_SIZE * 2, len * 2, "edgevertsb");
870         
871         i = 0;
872         l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
873         for (; l; l=BMIter_Step(&iter)) {
874                 BM_SetIndex(l, i); /* set_loop */
875                 copy_v3_v3(projverts[i], l->v->co);
876                 i++;
877         }
878         
879         for (i=0; i<len; i++) {
880                 copy_v3_v3(v1, loops[i][0]->v->co);
881                 copy_v3_v3(v2, loops[i][1]->v->co);
882
883                 shrink_edgef(v1, v2, fac2);
884                 
885                 copy_v3_v3(edgeverts[a], v1);
886                 a++;
887                 copy_v3_v3(edgeverts[a], v2);
888                 a++;
889         }
890         
891         compute_poly_normal(no, projverts, f->len);
892         poly_rotate_plane(no, projverts, f->len);
893         poly_rotate_plane(no, edgeverts, len*2);
894         
895         l = bm_firstfaceloop(f);
896         for (i=0; i<f->len; i++) {
897                 p1 = projverts[i];
898                 out[0] = MAX2(out[0], p1[0]) + 0.01f;
899                 out[1] = MAX2(out[1], p1[1]) + 0.01f;
900                 out[2] = 0.0f;
901                 p1[2] = 0.0f;
902
903                 //copy_v3_v3(l->v->co, p1);
904
905                 l = l->next;
906         }
907         
908         for (i=0; i<len; i++) {
909                 edgeverts[i*2][2] = 0.0f;
910                 edgeverts[i*2+1][2] = 0.0f;
911         }
912
913         /*do convexity test*/
914         for (i=0; i<len; i++) {
915                 copy_v3_v3(v2, edgeverts[i*2]);
916                 copy_v3_v3(v3, edgeverts[i*2+1]);
917
918                 mid_v3_v3v3(mid, v2, v3);
919                 
920                 clen = 0;
921                 for (j=0; j<f->len; j++) {
922                         p1 = projverts[j];
923                         p2 = projverts[(j+1)%f->len];
924                         
925                         copy_v3_v3(v1, p1);
926                         copy_v3_v3(v2, p2);
927
928                         shrink_edgef(v1, v2, fac1);
929
930                         if (linecrossesf(p1, p2, mid, out)) clen++;
931                 }
932                 
933                 if (clen%2 == 0) {
934                         loops[i][0] = NULL;
935                 }
936         }
937         
938         /*do line crossing tests*/
939         for (i=0; i<f->len; i++) {
940                 p1 = projverts[i];
941                 p2 = projverts[(i+1)%f->len];
942                 
943                 copy_v3_v3(v1, p1);
944                 copy_v3_v3(v2, p2);
945
946                 shrink_edgef(v1, v2, fac1);
947
948                 for (j=0; j<len; j++) {
949                         if (!loops[j][0]) continue;
950
951                         p3 = edgeverts[j*2];
952                         p4 = edgeverts[j*2+1];
953
954                         if (linecrossesf(v1, v2, p3, p4))
955                         {
956                                 loops[j][0] = NULL;
957                         }
958                 }
959         }
960
961         for (i=0; i<len; i++) {
962                 for (j=0; j<len; j++) {
963                         if (j == i) continue;
964                         if (!loops[i][0]) continue;
965                         if (!loops[j][0]) continue;
966
967                         p1 = edgeverts[i*2];
968                         p2 = edgeverts[i*2+1];
969                         p3 = edgeverts[j*2];
970                         p4 = edgeverts[j*2+1];
971
972                         copy_v3_v3(v1, p1);
973                         copy_v3_v3(v2, p2);
974
975                         shrink_edgef(v1, v2, fac1);
976
977                         if (linecrossesf(v1, v2, p3, p4)) {
978                                 loops[i][0]=NULL;
979                         }
980                 }
981         }
982
983         BLI_array_fixedstack_free(projverts);
984         BLI_array_fixedstack_free(edgeverts);
985 }