73f713c55910f6116b6aeb29444b4e938bc0d0b7
[blender.git] / source / blender / bmesh / intern / bmesh_polygon.c
1 #include "bmesh.h"
2 #include "bmesh_private.h"
3
4 #include "BLI_arithb.h"
5 #include "BKE_utildefines.h"
6 #include <string.h>
7
8 /*
9  *
10  * BME POLYGON.C
11  *
12  * This file contains code for dealing
13  * with polygons (normal/area calculation,
14  * tesselation, ect)
15  *
16  * TODO:
17  *   -Add in Tesselator frontend that creates
18  *     BMTriangles from copied faces
19  *  -Add in Function that checks for and flags
20  *   degenerate faces.
21  *
22 */
23
24 /*
25  * TEST EDGE SIDE and POINT IN TRIANGLE
26  *
27  * Point in triangle tests stolen from scanfill.c.
28  * Used for tesselator
29  *
30 */
31
32 static short testedgeside(float *v1, float *v2, float *v3)
33 /* is v3 to the right of v1-v2 ? With exception: v3==v1 || v3==v2 */
34 {
35         float inp;
36
37         //inp= (v2[cox]-v1[cox])*(v1[coy]-v3[coy]) +(v1[coy]-v2[coy])*(v1[cox]-v3[cox]);
38         inp= (v2[0]-v1[0])*(v1[1]-v3[1]) +(v1[1]-v2[1])*(v1[0]-v3[0]);
39
40         if(inp<0.0) return 0;
41         else if(inp==0) {
42                 if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
43                 if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
44         }
45         return 1;
46 }
47
48 static int point_in_triangle(float *v1, float *v2, float *v3, float *pt)
49 {
50         if(testedgeside(v1,v2,pt) && testedgeside(v2,v3,pt) && testedgeside(v3,v1,pt))
51                 return 1;
52         return 0;
53 }
54
55 /*
56  * CONVEX ANGLE 
57  *
58  * Tests whether or not a given angle in
59  * a polygon is convex or not. Note that 
60  * this assumes that the polygon has been
61  * projected to the x/y plane
62  *
63 */
64 static int convexangle(float *v1t, float *v2t, float *v3t)
65 {
66         float v1[3], v3[3], n[3];
67         VecSubf(v1, v1t, v2t);
68         VecSubf(v3, v3t, v2t);
69
70         Normalize(v1);
71         Normalize(v3);
72         Crossf(n, v1, v3);
73         
74         if(n[2] < 0.0)
75                 return 0;
76
77         return 1;
78 }
79
80 /*
81  * COMPUTE POLY NORMAL
82  *
83  * Computes the normal of a planar 
84  * polygon See Graphics Gems for 
85  * computing newell normal.
86  *
87 */
88 static void compute_poly_normal(float normal[3], float (*verts)[3], int nverts)
89 {
90
91         float *u,  *v;
92         int i;
93
94         normal[0] = 0.0;
95         normal[1] = 0.0;
96         normal[2] = 0.0;
97         
98         for(i = 0; i < nverts; i++){
99                 u = verts[i];
100                 v = verts[(i+1) % nverts];
101                 
102                 normal[0] += (u[1] - v[1]) * (u[2] + v[2]);
103                 normal[1] += (u[2] - v[2]) * (u[0] + v[0]);
104                 normal[2] += (u[0] - v[0]) * (u[1] + v[1]);
105                 i++;
106         }
107
108         Normalize(normal);
109 }
110
111 /*
112  * COMPUTE POLY CENTER
113  *
114  * Computes the centroid and
115  * area of a polygon in the X/Y
116  * plane.
117  *
118 */
119
120 static int compute_poly_center(float center[3], float *area, float (*verts)[3], int nverts)
121 {
122         int i, j;
123         float atmp = 0.0, xtmp = 0.0, ytmp = 0.0, ai;
124         
125         center[0] = center[1] = center[2] = 0.0;        
126
127         if(nverts < 3) 
128                 return 0;
129
130         i = nverts-1;
131         j = 0;
132         
133         while(j < nverts){
134                 ai = verts[i][0] * verts[j][1] - verts[j][0] * verts[i][1];                             
135                 atmp += ai;
136                 xtmp += (verts[j][0] + verts[i][0]) * ai;
137                 ytmp += (verts[j][1] + verts[i][1]) * ai;
138                 i = j;
139                 j += 1;
140         }
141
142         if(area)
143                 *area = atmp / 2.0f;    
144         
145         if (atmp != 0){
146                 center[0] = xtmp /  (3.0f * atmp);
147                 center[1] = xtmp /  (3.0f * atmp);
148                 return 1;
149         }
150         return 0;
151 }
152
153
154 /*
155  * COMPUTE POLY PLANE
156  *
157  * Projects a set polygon's vertices to 
158  * a plane defined by the average
159  * of its edges cross products
160  *
161 */
162
163 void compute_poly_plane(float (*verts)[3], int nverts)
164 {
165         
166         float avgc[3], norm[3], temp[3], mag, avgn[3];
167         float *v1, *v2, *v3;
168         int i;
169         
170         if(nverts < 3) 
171                 return;
172
173         avgn[0] = avgn[1] = avgn[2] = 0.0;
174         avgc[0] = avgc[1] = avgc[2] = 0.0;
175
176         for(i = 0; i < nverts; i++){
177                 v1 = verts[i];
178                 v2 = verts[(i+1) % nverts];
179                 v3 = verts[(i+2) % nverts];
180                 CalcNormFloat(v1, v2, v3, norm);        
181         
182                 avgn[0] += norm[0];
183                 avgn[1] += norm[1];
184                 avgn[2] += norm[2];
185         }
186
187         /*what was this bit for?*/
188         if(avgn[0] == 0.0 && avgn[1] == 0.0 && avgn[2] == 0.0){
189                 avgn[0] = 0.0;
190                 avgn[1] = 0.0;
191                 avgn[2] = 1.0;
192         } else {
193                 avgn[0] /= nverts;
194                 avgn[1] /= nverts;
195                 avgn[2] /= nverts;
196                 Normalize(avgn);
197         }
198         
199         for(i = 0; i < nverts; i++){
200                 v1 = verts[i];
201                 VECCOPY(temp, v1);
202                 mag = 0.0;
203                 mag += (temp[0] * avgn[0]);
204                 mag += (temp[1] * avgn[1]);
205                 mag += (temp[2] * avgn[2]);
206                 
207                 temp[0] = (avgn[0] * mag);
208                 temp[1] = (avgn[1] * mag);
209                 temp[2] = (avgn[2] * mag);
210
211                 VecSubf(v1, v1, temp);
212         }       
213 }
214
215 /*
216  * POLY ROTATE PLANE
217  *
218  * Rotates a polygon so that it's
219  * normal is pointing towards the mesh Z axis
220  *
221 */
222
223 void poly_rotate_plane(float normal[3], float (*verts)[3], int nverts)
224 {
225
226         float up[3] = {0.0,0.0,1.0}, axis[3], angle, q[4];
227         float mat[3][3];
228         int i;
229
230         Crossf(axis, up, normal);
231         angle = VecAngle2(normal, up);
232
233         if (angle == 0.0) return;
234         
235         AxisAngleToQuat(q, axis, angle);
236         QuatToMat3(q, mat);
237
238         for(i = 0;  i < nverts;  i++)
239                 Mat3MulVecfl(mat, verts[i]);
240 }
241
242 /*
243  * BMESH UPDATE FACE NORMAL
244  *
245  * Updates the stored normal for the
246  * given face. Requires that a buffer
247  * of sufficient length to store projected
248  * coordinates for all of the face's vertices
249  * is passed in as well.
250  *
251 */
252
253 void bmesh_update_face_normal(BMesh *bm, BMFace *f, float (*projectverts)[3])
254 {
255         BMLoop *l;
256         int i;
257
258         if(f->len > 4){
259                 i = 0;
260                 l = f->loopbase;
261                 do{
262                         VECCOPY(projectverts[i], l->v->co);
263                         l = (BMLoop*)(l->head.next);
264                 }while(l!=f->loopbase);
265
266                 compute_poly_plane(projectverts, f->len);
267                 compute_poly_normal(f->no, projectverts, f->len);       
268         }
269         else if(f->len == 3){
270                 BMVert *v1, *v2, *v3;
271                 v1 = f->loopbase->v;
272                 v2 = ((BMLoop*)(f->loopbase->head.next))->v;
273                 v3 = ((BMLoop*)(f->loopbase->head.next->next))->v;
274                 CalcNormFloat(v1->co, v2->co, v3->co, f->no);
275         }
276         else if(f->len == 4){
277                 BMVert *v1, *v2, *v3, *v4;
278                 v1 = f->loopbase->v;
279                 v2 = ((BMLoop*)(f->loopbase->head.next))->v;
280                 v3 = ((BMLoop*)(f->loopbase->head.next->next))->v;
281                 v4 = ((BMLoop*)(f->loopbase->head.prev))->v;
282                 CalcNormFloat4(v1->co, v2->co, v3->co, v4->co, f->no);
283         }
284         else{ /*horrible, two sided face!*/
285                 f->no[0] = 0.0;
286                 f->no[1] = 0.0;
287                 f->no[2] = 1.0;
288         }
289
290 }
291
292
293 /*
294  * BMESH FLIP NORMAL
295  * 
296  *  Reverses the winding of a  faces
297  *  Note that this does *not* update the calculated 
298  *  Normal 
299 */
300 void BM_flip_normal(BMesh *bm, BMFace *f)
301 {       
302         bmesh_loop_reverse(bm, f);
303 }
304
305
306
307 int winding(float *a, float *b, float *c)
308 {
309         float v1[3], v2[3], v[3];
310         
311         VecSubf(v1, b, a);
312         VecSubf(v2, b, c);
313         
314         v1[2] = 0;
315         v2[2] = 0;
316         
317         Normalize(v1);
318         Normalize(v2);
319         
320         Crossf(v, v1, v2);
321
322         /*!! (turns nonzero into 1) is likely not necassary, 
323           since '>' I *think* should always
324           return 0 or 1, but I'm not totally sure. . .*/
325         return !!(v[2] > 0);
326 }
327
328 /* detects if two line segments cross each other (intersects).
329    note, there could be more winding cases then there needs to be. */
330 int linecrosses(float *v1, float *v2, float *v3, float *v4)
331 {
332         int w1, w2, w3, w4, w5;
333
334         w1 = winding(v1, v3, v2);
335         w2 = winding(v2, v4, v1);
336         w3 = !winding(v1, v2, v3);
337         w4 = winding(v3, v2, v4);
338         w5 = !winding(v3, v1, v4);
339         return w1 == w2 && w2 == w3 && w3 == w4 && w4==w5;
340 }
341
342 int goodline(float (*projectverts)[3], int v1i, int v2i, int nvert)
343 {
344         /*the hardcoded stuff here, 200000, 0.999, and 1.0001 may be problems
345           in the future, not sure. - joeedh*/
346         static float outv[3] = {2000.0f, 2000.0f, 0.0f};
347         float v1[3], v2[3], p[3], vv1[3], vv2[3], mid[3], a[3], b[3];
348         int i = 0, lcount=0;
349         
350         VECCOPY(v1, projectverts[v1i])
351         VECCOPY(v2, projectverts[v2i])
352         
353         VecAddf(p, v1, v2);
354         VecMulf(p, 0.5f);
355         
356         VecSubf(a, v1, p);
357         VecSubf(b, v2, p);
358         
359         VecMulf(a, 0.999f);
360         VecMulf(b, 0.999f);
361         
362         VecAddf(v1, p, a);
363         VecAddf(v2, p, b);
364         
365         while (i < nvert) {
366                 VECCOPY(vv1, projectverts[i]);
367                 VECCOPY(vv2, projectverts[(i+1)%nvert]);
368                 
369                 if (linecrosses(vv1, vv2, v1, v2)) return 0;
370
371                 VecAddf(mid, vv1, vv2);
372                 VecMulf(mid, 0.5f);
373                 
374                 VecSubf(a, vv1, mid);
375                 VecSubf(b, vv2, mid);
376                 
377                 VecMulf(a, 1.0001f);
378                 VecMulf(b, 1.0001f);
379                 
380                 VecAddf(vv1, mid, a);
381                 VecAddf(vv2, mid, b);
382                                 
383                 if (linecrosses(vv1, vv2, p, outv)) lcount += 1;
384
385                 i += 1;
386         }
387         if ((lcount % 2) == 0) return 0;
388
389         return 1;
390 }
391
392 /*
393  * FIND EAR
394  *
395  * Used by tesselator to find
396  * the next triangle to 'clip off'
397  * of a polygon while tesselating.
398  *
399 */
400
401 static BMLoop *find_ear(BMFace *f, float (*verts)[3], int nvert)
402 {
403         BMVert *v1, *v2, *v3;
404         BMLoop *bestear = NULL, *l;
405         float angle, bestangle = 180.0f;
406         int isear;
407         
408         l = f->loopbase;
409         do{
410                 isear = 1;
411                 
412                 v1 = ((BMLoop*)(l->head.prev))->v;
413                 v2 = l->v;
414                 v3 = ((BMLoop*)(l->head.next))->v;
415
416                 if (BM_Edge_Exist(v1, v3)) isear = 0;
417
418                 if (isear && !goodline(verts, v1->head.eflag2, v3->head.eflag2, nvert)) isear = 0;
419                 if(isear){
420                         angle = VecAngle3(verts[v1->head.eflag2], verts[v2->head.eflag2], verts[v3->head.eflag2]);
421                         if(!bestear || ABS(angle-40.0) < bestangle){
422                                 bestear = l;
423                                 bestangle = ABS(40.0-angle);
424                         }
425                 }
426                 l = (BMLoop*)(l->head.next);
427         }
428         while(l != f->loopbase);
429
430         return bestear;
431 }
432
433 /*
434  * BMESH TRIANGULATE FACE
435  *
436  * Triangulates a face using a 
437  * simple 'ear clipping' algorithm
438  * that tries to favor non-skinny
439  * triangles (angles less than 
440  * 90 degrees). If the triangulator
441  * has bits left over (or cannot
442  * triangulate at all) it uses an 
443  * arbitrary triangulation.
444  *
445  * TODO:
446  * -Modify this to try and find ears that will not create a non-manifold face after conversion back to editmesh
447  *
448 */
449 void BM_Triangulate_Face(BMesh *bm, BMFace *f, float (*projectverts)[3], int newedgeflag, int newfaceflag)
450 {
451         int i, done, nvert;
452         BMLoop *l, *newl, *nextloop;
453
454         /*copy vertex coordinates to vertspace array*/
455         i = 0;
456         l = f->loopbase;
457         do{
458                 VECCOPY(projectverts[i], l->v->co);
459                 l->v->head.eflag2 = i; /*warning, abuse! never duplicate in tools code! never you hear?*/ /*actually, get rid of this completely, use a new structure for this....*/
460                 i++;
461                 l = (BMLoop*)(l->head.next);
462         }while(l != f->loopbase);
463         
464         //bmesh_update_face_normal(bm, f, projectverts);
465
466         compute_poly_normal(f->no, projectverts, f->len);       
467         compute_poly_plane(projectverts, i);
468         poly_rotate_plane(f->no, projectverts, i);
469         
470         nvert = f->len;
471         
472         done = 0;
473         while(!done && f->len > 3){
474                 done = 1;
475                 l = find_ear(f, projectverts, nvert);
476                 if(l) {
477                         done = 0;
478                         f = bmesh_sfme(bm, f, ((BMLoop*)(l->head.prev))->v, ((BMLoop*)(l->head.next))->v, &newl);
479                         if (!f) {
480                                 printf("yeek! triangulator failed to split face!\n");
481                                 break;
482                         }
483
484                         BMO_SetFlag(bm, newl->e, newedgeflag);
485                         BMO_SetFlag(bm, f, newfaceflag);
486                 }
487         }
488
489         if (f->len > 3){
490                 l = f->loopbase;
491                 while (l->f->len > 3){
492                         nextloop = ((BMLoop*)(l->head.next->next));
493                         f = bmesh_sfme(bm, l->f, l->v,nextloop->v, &newl);
494                         if (!f) {
495                                 printf("triangle fan step of triangulator failed.\n");
496                                 return;
497                         }
498
499                         BMO_SetFlag(bm, newl->e, newedgeflag);
500                         BMO_SetFlag(bm, f, newfaceflag);
501                         l = nextloop;
502                 }
503         }
504 }