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