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