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