62755d6b55845c32c8ceb9661af103adf9e744a2
[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         } else l = 1.0f / l;
174
175         mul_v3_fl(n, l);
176
177         copy_v3_v3(normal, n);
178 #endif
179 }
180
181 /*
182  * COMPUTE POLY CENTER
183  *
184  * Computes the centroid and
185  * area of a polygon in the X/Y
186  * plane.
187  *
188 */
189
190 static int compute_poly_center(float center[3], float *r_area, float (*verts)[3], int nverts)
191 {
192         int i, j;
193         float atmp = 0.0, xtmp = 0.0, ytmp = 0.0, ai;
194
195         zero_v3(center);
196
197         if(nverts < 3) 
198                 return 0;
199
200         i = nverts-1;
201         j = 0;
202         
203         while(j < nverts) {
204                 ai = verts[i][0] * verts[j][1] - verts[j][0] * verts[i][1];                             
205                 atmp += ai;
206                 xtmp += (verts[j][0] + verts[i][0]) * ai;
207                 ytmp += (verts[j][1] + verts[i][1]) * ai;
208                 i = j;
209                 j += 1;
210         }
211
212         if(r_area)
213                 *r_area = atmp / 2.0f;
214         
215         if (atmp != 0) {
216                 center[0] = xtmp /  (3.0f * atmp);
217                 center[1] = xtmp /  (3.0f * atmp);
218                 return 1;
219         }
220         return 0;
221 }
222
223 float BM_Compute_Face_Area(BMesh *bm, BMFace *f)
224 {
225         BMLoop *l;
226         BMIter iter;
227         float (*verts)[3];
228         float center[3];
229         float area = 0.0f;
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, bm, 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         if (f->len) 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[3], float v2[3], const 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, f->no, proj);
424
425                 BLI_array_fixedstack_free(proj);
426         }
427 }
428 /* same as BM_Face_UpdateNormal but takes vertex coords */
429 void BM_Face_UpdateNormal_VertexCos(BMesh *bm, BMFace *f, float no[3], float (*vertexCos)[3])
430 {
431         if (f->len >= 3) {
432                 float (*proj)[3];
433
434                 BLI_array_fixedstack_declare(proj, BM_NGON_STACK_SIZE, f->len, __func__);
435
436                 bmesh_update_face_normal_vertex_cos(bm, f, no, proj, vertexCos);
437
438                 BLI_array_fixedstack_free(proj);
439         }
440 }
441
442 void BM_Edge_UpdateNormals(BMesh *bm, BMEdge *e)
443 {
444         BMIter iter;
445         BMFace *f;
446         
447         f = BMIter_New(&iter, bm, BM_FACES_OF_EDGE, e);
448         for (; f; f=BMIter_Step(&iter)) {
449                 BM_Face_UpdateNormal(bm, f);
450         }
451
452         BM_Vert_UpdateNormal(bm, e->v1);
453         BM_Vert_UpdateNormal(bm, e->v2);
454 }
455
456 void BM_Vert_UpdateNormal(BMesh *bm, BMVert *v)
457 {
458         BMIter eiter, liter;
459         BMEdge *e;
460         BMLoop *l;
461         float vec1[3], vec2[3], fac;
462         int len=0;
463
464         zero_v3(v->no);
465
466         BM_ITER(e, &eiter, bm, BM_EDGES_OF_VERT, v) {
467                 BM_ITER(l, &liter, bm, BM_LOOPS_OF_EDGE, e) {
468                         if (l->v == v) {
469                                 /* Same calculation used in BM_Compute_Normals */
470                                 sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
471                                 sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
472                                 normalize_v3(vec1);
473                                 normalize_v3(vec2);
474
475                                 fac = saacos(-dot_v3v3(vec1, vec2));
476                                 
477                                 madd_v3_v3fl(v->no, l->f->no, fac);
478
479                                 len++;
480                         }
481                 }
482         }
483
484         if (len) {
485                 normalize_v3(v->no);
486         }
487 }
488
489 void BM_Vert_UpdateAllNormals(BMesh *bm, BMVert *v)
490 {
491         BMIter iter;
492         BMFace *f;
493         int len=0;
494
495         f = BMIter_New(&iter, bm, BM_FACES_OF_VERT, v);
496         for (; f; f=BMIter_Step(&iter), len++) {
497                 BM_Face_UpdateNormal(bm, f);
498         }
499
500         BM_Vert_UpdateNormal(bm, v);
501 }
502
503 void bmesh_update_face_normal(BMesh *bm, BMFace *f, float no[3],
504                               float (*projectverts)[3])
505 {
506         BMLoop *l;
507
508         /* common cases first */
509         switch (f->len) {
510                 case 4:
511                 {
512                         BMVert *v1 = (l = bm_firstfaceloop(f))->v;
513                         BMVert *v2 = (l = l->next)->v;
514                         BMVert *v3 = (l = l->next)->v;
515                         BMVert *v4 = (l->next)->v;
516                         normal_quad_v3(no, v1->co, v2->co, v3->co, v4->co);
517                         break;
518                 }
519                 case 3:
520                 {
521                         BMVert *v1 = (l = bm_firstfaceloop(f))->v;
522                         BMVert *v2 = (l = l->next)->v;
523                         BMVert *v3 = (l->next)->v;
524                         normal_tri_v3(no, v1->co, v2->co, v3->co);
525                         break;
526                 }
527                 case 0:
528                 {
529                         zero_v3(no);
530                         break;
531                 }
532                 default:
533                 {
534                         BMIter iter;
535                         int i = 0;
536                         BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
537                                 copy_v3_v3(projectverts[i], l->v->co);
538                                 i += 1;
539                         }
540                         compute_poly_normal(no, projectverts, f->len);
541                         break;
542                 }
543         }
544 }
545 /* exact same as 'bmesh_update_face_normal' but accepts vertex coords */
546 void bmesh_update_face_normal_vertex_cos(BMesh *bm, BMFace *f, float no[3],
547                                    float (*projectverts)[3], float (*vertexCos)[3])
548 {
549         BMLoop *l;
550
551         /* must have valid index data */
552         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
553
554         /* common cases first */
555         switch (f->len) {
556                 case 4:
557                 {
558                         BMVert *v1 = (l = bm_firstfaceloop(f))->v;
559                         BMVert *v2 = (l = l->next)->v;
560                         BMVert *v3 = (l = l->next)->v;
561                         BMVert *v4 = (l->next)->v;
562                         normal_quad_v3(no,
563                                        vertexCos[BM_GetIndex(v1)],
564                                        vertexCos[BM_GetIndex(v2)],
565                                        vertexCos[BM_GetIndex(v3)],
566                                        vertexCos[BM_GetIndex(v4)]);
567                         break;
568                 }
569                 case 3:
570                 {
571                         BMVert *v1 = (l = bm_firstfaceloop(f))->v;
572                         BMVert *v2 = (l = l->next)->v;
573                         BMVert *v3 = (l->next)->v;
574                         normal_tri_v3(no,
575                                       vertexCos[BM_GetIndex(v1)],
576                                       vertexCos[BM_GetIndex(v2)],
577                                       vertexCos[BM_GetIndex(v3)]);
578                         break;
579                 }
580                 case 0:
581                 {
582                         zero_v3(no);
583                         break;
584                 }
585                 default:
586                 {
587                         BMIter iter;
588                         int i = 0;
589                         BM_ITER(l, &iter, bm, BM_LOOPS_OF_FACE, f) {
590                                 copy_v3_v3(projectverts[i], vertexCos[BM_GetIndex(l->v)]);
591                                 i += 1;
592                         }
593                         compute_poly_normal(no, projectverts, f->len);
594                         break;
595                 }
596         }
597 }
598
599 /*
600  * BMESH FLIP NORMAL
601  * 
602  *  Reverses the winding of a face.
603  *  Note that this updates the calculated 
604  *  normal.
605 */
606 void BM_flip_normal(BMesh *bm, BMFace *f)
607 {       
608         bmesh_loop_reverse(bm, f);
609         negate_v3(f->no);
610 }
611
612 /* detects if two line segments cross each other (intersects).
613    note, there could be more winding cases then there needs to be. */
614 static int UNUSED_FUNCTION(linecrosses)(const double v1[2], const double v2[2], const double v3[2], const double v4[2])
615 {
616         int w1, w2, w3, w4, w5;
617         
618         /*w1 = winding(v1, v3, v4);
619         w2 = winding(v2, v3, v4);
620         w3 = winding(v3, v1, v2);
621         w4 = winding(v4, v1, v2);
622         
623         return (w1 == w2) && (w3 == w4);*/
624
625         w1 = testedgeside(v1, v3, v2);
626         w2 = testedgeside(v2, v4, v1);
627         w3 = !testedgeside(v1, v2, v3);
628         w4 = testedgeside(v3, v2, v4);
629         w5 = !testedgeside(v3, v1, v4);
630         return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5;
631 }
632
633 /* detects if two line segments cross each other (intersects).
634    note, there could be more winding cases then there needs to be. */
635 static int linecrossesf(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
636 {
637         int w1, w2, w3, w4, w5 /*, ret*/;
638         float mv1[2], mv2[2], mv3[2], mv4[2];
639         
640         /*now test winding*/
641         w1 = testedgesidef(v1, v3, v2);
642         w2 = testedgesidef(v2, v4, v1);
643         w3 = !testedgesidef(v1, v2, v3);
644         w4 = testedgesidef(v3, v2, v4);
645         w5 = !testedgesidef(v3, v1, v4);
646         
647         if (w1 == w2 && w2 == w3 && w3 == w4 && w4==w5)
648                 return 1;
649         
650 #define GETMIN2_AXIS(a, b, ma, mb, axis) ma[axis] = MIN2(a[axis], b[axis]), mb[axis] = MAX2(a[axis], b[axis])
651 #define GETMIN2(a, b, ma, mb) GETMIN2_AXIS(a, b, ma, mb, 0); GETMIN2_AXIS(a, b, ma, mb, 1);
652         
653         GETMIN2(v1, v2, mv1, mv2);
654         GETMIN2(v3, v4, mv3, mv4);
655         
656         /*do an interval test on the x and y axes*/
657         /*first do x axis*/
658         #define T FLT_EPSILON*15
659         if (ABS(v1[1]-v2[1]) < T && ABS(v3[1]-v4[1]) < T &&
660             ABS(v1[1]-v3[1]) < T) 
661         {
662                 return (mv4[0] >= mv1[0] && mv3[0] <= mv2[0]);
663         }
664
665         /*now do y axis*/
666         if (ABS(v1[0]-v2[0]) < T && ABS(v3[0]-v4[0]) < T &&
667             ABS(v1[0]-v3[0]) < T)
668         {
669                 return (mv4[1] >= mv1[1] && mv3[1] <= mv2[1]);
670         }
671
672         return 0; 
673 }
674
675 /*
676    BM POINT IN FACE
677    
678   Projects co onto face f, and returns true if it is inside
679   the face bounds.  Note that this uses a best-axis projection
680   test, instead of projecting co directly into f's orientation
681   space, so there might be accuracy issues.
682 */
683 int BM_Point_In_Face(BMesh *bm, BMFace *f, const float co[3])
684 {
685         int ax, ay;
686         float co2[3], cent[3] = {0.0f, 0.0f, 0.0f}, out[3] = {FLT_MAX*0.5f, FLT_MAX*0.5f, 0};
687         BMLoop *l;
688         int crosses = 0;
689         float eps = 1.0f+(float)FLT_EPSILON*150.0f;
690         
691         if (dot_v3v3(f->no, f->no) <= FLT_EPSILON*10)
692                 BM_Face_UpdateNormal(bm, f);
693         
694         /* find best projection of face XY, XZ or YZ: barycentric weights of
695          * the 2d projected coords are the same and faster to compute
696          *
697          * this probably isn't all that accurate, but it has the advantage of
698          * being fast (especially compared to projecting into the face orientation)
699          */
700         axis_dominant_v3(&ax, &ay, f->no);
701
702         co2[0] = co[ax];
703         co2[1] = co[ay];
704         co2[2] = 0;
705         
706         l = bm_firstfaceloop(f);
707         do {
708                 cent[0] += l->v->co[ax];
709                 cent[1] += l->v->co[ay];
710                 l = l->next;
711         } while (l != bm_firstfaceloop(f));
712         
713         mul_v2_fl(cent, 1.0f/(float)f->len);
714         
715         l = bm_firstfaceloop(f);
716         do {
717                 float v1[3], v2[3];
718                 
719                 v1[0] = (l->prev->v->co[ax] - cent[ax])*eps + cent[ax];
720                 v1[1] = (l->prev->v->co[ay] - cent[ay])*eps + cent[ay];
721                 v1[2] = 0.0f;
722                 
723                 v2[0] = (l->v->co[ax] - cent[ax])*eps + cent[ax];
724                 v2[1] = (l->v->co[ay] - cent[ay])*eps + cent[ay];
725                 v2[2] = 0.0f;
726                 
727                 crosses += linecrossesf(v1, v2, co2, out) != 0;
728                 
729                 l = l->next;
730         } while (l != bm_firstfaceloop(f));
731         
732         return crosses%2 != 0;
733 }
734
735 static int goodline(float (*projectverts)[3], BMFace *f, int v1i,
736                     int v2i, int v3i, int UNUSED(nvert))
737 {
738         BMLoop *l = bm_firstfaceloop(f);
739         double v1[3], v2[3], v3[3], pv1[3], pv2[3];
740         int i;
741
742         VECCOPY(v1, projectverts[v1i]);
743         VECCOPY(v2, projectverts[v2i]);
744         VECCOPY(v3, projectverts[v3i]);
745         
746         if (testedgeside(v1, v2, v3)) return 0;
747         
748         //for (i=0; i<nvert; i++) {
749         do {
750                 i = BM_GetIndex(l->v);
751                 if (i == v1i || i == v2i || i == v3i) {
752                         l = l->next;
753                         continue;
754                 }
755                 
756                 VECCOPY(pv1, projectverts[BM_GetIndex(l->v)]);
757                 VECCOPY(pv2, projectverts[BM_GetIndex(l->next->v)]);
758                 
759                 //if (linecrosses(pv1, pv2, v1, v3)) return 0;
760                 if (point_in_triangle(v1, v2, v3, pv1)) return 0;
761                 if (point_in_triangle(v3, v2, v1, pv1)) return 0;
762
763                 l = l->next;
764         } while (l != bm_firstfaceloop(f));
765         return 1;
766 }
767 /*
768  * FIND EAR
769  *
770  * Used by tesselator to find
771  * the next triangle to 'clip off'
772  * of a polygon while tesselating.
773  *
774 */
775
776 static BMLoop *find_ear(BMesh *UNUSED(bm), BMFace *f, float (*verts)[3],
777                         int nvert)
778 {
779         BMVert *v1, *v2, *v3;
780         BMLoop *bestear = NULL, *l;
781         /*float angle, bestangle = 180.0f;*/
782         int isear /*, i=0*/;
783         
784         l = bm_firstfaceloop(f);
785         do {
786                 isear = 1;
787                 
788                 v1 = l->prev->v;
789                 v2 = l->v;
790                 v3 = l->next->v;
791
792                 if (BM_Edge_Exist(v1, v3)) isear = 0;
793
794                 if (isear && !goodline(verts, f, BM_GetIndex(v1), BM_GetIndex(v2),
795                                        BM_GetIndex(v3), nvert))
796                         isear = 0;
797                 
798                 if(isear) {
799                         /*angle = angle_v3v3v3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
800                         if(!bestear || ABS(angle-45.0f) < bestangle) {
801                                 bestear = l;
802                                 bestangle = ABS(45.0f-angle);
803                         }
804                         
805                         if (angle > 20 && angle < 90) break;
806                         if (angle < 100 && i > 5) break;
807                         i += 1;*/
808                         bestear = l;
809                         break;
810                 }
811                 l = l->next;
812         }
813         while(l != bm_firstfaceloop(f));
814
815         return bestear;
816 }
817
818 /*
819  * BMESH TRIANGULATE FACE
820  *
821  * Triangulates a face using a 
822  * simple 'ear clipping' algorithm
823  * that tries to favor non-skinny
824  * triangles (angles less than 
825  * 90 degrees). If the triangulator
826  * has bits left over (or cannot
827  * triangulate at all) it uses a
828  * simple fan triangulation
829  *
830  * newfaces, if non-null, must be an array of BMFace pointers,
831  * with a length equal to f->len.  it will be filled with the new
832  * triangles, and will be NULL-terminated.
833 */
834 void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], 
835                          int newedgeflag, int newfaceflag, BMFace **newfaces)
836 {
837         int i, done, nvert, nf_i = 0;
838         BMLoop *l, *newl, *nextloop;
839         /* BMVert *v; */ /* UNUSED */
840
841         /*copy vertex coordinates to vertspace array*/
842         i = 0;
843         l = bm_firstfaceloop(f);
844         do {
845                 copy_v3_v3(projectverts[i], l->v->co);
846                 BM_SetIndex(l->v, i); /* set dirty! */
847                 i++;
848                 l = l->next;
849         } while(l != bm_firstfaceloop(f));
850
851         bm->elem_index_dirty |= BM_VERT; /* see above */
852
853         ///bmesh_update_face_normal(bm, f, f->no, projectverts);
854
855         compute_poly_normal(f->no, projectverts, f->len);
856         poly_rotate_plane(f->no, projectverts, i);
857
858         nvert = f->len;
859
860         //compute_poly_plane(projectverts, i);
861         for (i=0; i<nvert; i++) {
862                 projectverts[i][2] = 0.0f;
863         }
864
865         done = 0;
866         while(!done && f->len > 3) {
867                 done = 1;
868                 l = find_ear(bm, f, projectverts, nvert);
869                 if(l) {
870                         done = 0;
871                         /* v = l->v; */ /* UNUSED */
872                         f = BM_Split_Face(bm, l->f, l->prev->v,
873                                           l->next->v,
874                                           &newl, NULL);
875                         copy_v3_v3(f->no, l->f->no);
876
877                         if (!f) {
878                                 fprintf(stderr, "%s: triangulator failed to split face! (bmesh internal error)\n", __func__);
879                                 break;
880                         }
881
882                         BMO_SetFlag(bm, newl->e, newedgeflag);
883                         BMO_SetFlag(bm, f, newfaceflag);
884                         
885                         if (newfaces) newfaces[nf_i++] = f;
886
887                         /*l = f->loopbase;
888                         do {
889                                 if (l->v == v) {
890                                         f->loopbase = l;
891                                         break;
892                                 }
893                                 l = l->next;
894                         } while (l != f->loopbase);*/
895                 }
896         }
897
898         if (f->len > 3) {
899                 l = bm_firstfaceloop(f);
900                 while (l->f->len > 3) {
901                         nextloop = l->next->next;
902                         f = BM_Split_Face(bm, l->f, l->v, nextloop->v, 
903                                           &newl, NULL);
904                         if (!f) {
905                                 printf("triangle fan step of triangulator failed.\n");
906
907                                 /*NULL-terminate*/
908                                 if (newfaces) newfaces[nf_i] = NULL;
909                                 return;
910                         }
911
912                         if (newfaces) newfaces[nf_i++] = f;
913                         
914                         BMO_SetFlag(bm, newl->e, newedgeflag);
915                         BMO_SetFlag(bm, f, newfaceflag);
916                         l = nextloop;
917                 }
918         }
919         
920         /*NULL-terminate*/
921         if (newfaces) newfaces[nf_i] = NULL;
922 }
923
924 /*each pair of loops defines a new edge, a split.  this function goes
925   through and sets pairs that are geometrically invalid to null.  a
926   split is invalid, if it forms a concave angle or it intersects other
927   edges in the face, or it intersects another split.  in the case of
928   intersecting splits, only the first of the set of intersecting
929   splits survives.*/
930 void BM_LegalSplits(BMesh *bm, BMFace *f, BMLoop *(*loops)[2], int len)
931 {
932         BMIter iter;
933         BMLoop *l;
934         float v1[3], v2[3], v3[3]/*, v4[3]*/, no[3], mid[3], *p1, *p2, *p3, *p4;
935         float out[3] = {-234324.0f, -234324.0f, 0.0f};
936         float (*projverts)[3];
937         float (*edgeverts)[3];
938         float fac1 = 1.0000001f, fac2 = 0.9f; //9999f; //0.999f;
939         int i, j, a=0, clen;
940
941         BLI_array_fixedstack_declare(projverts, BM_NGON_STACK_SIZE, f->len,         "projvertsb");
942         BLI_array_fixedstack_declare(edgeverts, BM_NGON_STACK_SIZE * 2, len * 2, "edgevertsb");
943         
944         i = 0;
945         l = BMIter_New(&iter, bm, BM_LOOPS_OF_FACE, f);
946         for (; l; l=BMIter_Step(&iter)) {
947                 BM_SetIndex(l, i); /* set_loop */
948                 copy_v3_v3(projverts[i], l->v->co);
949                 i++;
950         }
951         
952         for (i=0; i<len; i++) {
953                 copy_v3_v3(v1, loops[i][0]->v->co);
954                 copy_v3_v3(v2, loops[i][1]->v->co);
955
956                 shrink_edgef(v1, v2, fac2);
957                 
958                 copy_v3_v3(edgeverts[a], v1);
959                 a++;
960                 copy_v3_v3(edgeverts[a], v2);
961                 a++;
962         }
963         
964         compute_poly_normal(no, projverts, f->len);
965         poly_rotate_plane(no, projverts, f->len);
966         poly_rotate_plane(no, edgeverts, len*2);
967         
968         l = bm_firstfaceloop(f);
969         for (i=0; i<f->len; i++) {
970                 p1 = projverts[i];
971                 out[0] = MAX2(out[0], p1[0]) + 0.01f;
972                 out[1] = MAX2(out[1], p1[1]) + 0.01f;
973                 out[2] = 0.0f;
974                 p1[2] = 0.0f;
975
976                 //copy_v3_v3(l->v->co, p1);
977
978                 l = l->next;
979         }
980         
981         for (i=0; i<len; i++) {
982                 edgeverts[i*2][2] = 0.0f;
983                 edgeverts[i*2+1][2] = 0.0f;
984         }
985
986         /*do convexity test*/
987         for (i=0; i<len; i++) {
988                 copy_v3_v3(v2, edgeverts[i*2]);
989                 copy_v3_v3(v3, edgeverts[i*2+1]);
990
991                 mid_v3_v3v3(mid, v2, v3);
992                 
993                 clen = 0;
994                 for (j=0; j<f->len; j++) {
995                         p1 = projverts[j];
996                         p2 = projverts[(j+1)%f->len];
997                         
998                         copy_v3_v3(v1, p1);
999                         copy_v3_v3(v2, p2);
1000
1001                         shrink_edgef(v1, v2, fac1);
1002
1003                         if (linecrossesf(p1, p2, mid, out)) clen++;
1004                 }
1005                 
1006                 if (clen%2 == 0) {
1007                         loops[i][0] = NULL;
1008                 }
1009         }
1010         
1011         /*do line crossing tests*/
1012         for (i=0; i<f->len; i++) {
1013                 p1 = projverts[i];
1014                 p2 = projverts[(i+1)%f->len];
1015                 
1016                 copy_v3_v3(v1, p1);
1017                 copy_v3_v3(v2, p2);
1018
1019                 shrink_edgef(v1, v2, fac1);
1020
1021                 for (j=0; j<len; j++) {
1022                         if (!loops[j][0]) continue;
1023
1024                         p3 = edgeverts[j*2];
1025                         p4 = edgeverts[j*2+1];
1026
1027                         if (linecrossesf(v1, v2, p3, p4))
1028                         {
1029                                 loops[j][0] = NULL;
1030                         }
1031                 }
1032         }
1033
1034         for (i=0; i<len; i++) {
1035                 for (j=0; j<len; j++) {
1036                         if (j == i) continue;
1037                         if (!loops[i][0]) continue;
1038                         if (!loops[j][0]) continue;
1039
1040                         p1 = edgeverts[i*2];
1041                         p2 = edgeverts[i*2+1];
1042                         p3 = edgeverts[j*2];
1043                         p4 = edgeverts[j*2+1];
1044
1045                         copy_v3_v3(v1, p1);
1046                         copy_v3_v3(v2, p2);
1047
1048                         shrink_edgef(v1, v2, fac1);
1049
1050                         if (linecrossesf(v1, v2, p3, p4)) {
1051                                 loops[i][0]=NULL;
1052                         }
1053                 }
1054         }
1055
1056         BLI_array_fixedstack_free(projverts);
1057         BLI_array_fixedstack_free(edgeverts);
1058 }