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