source/blender/bmesh/tools/bmesh_decimate_collapse.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): Campbell Barton
  19  *
  20  * ***** END GPL LICENSE BLOCK *****
  21  */
  22
  23 /** \file blender/bmesh/tools/bmesh_decimate_collapse.c
  24  *  \ingroup bmesh
  25  *
  26  * BMesh decimator that uses an edge collapse method.
  27  */
  28
  29 #include <stddef.h>
  30
  31 #include "MEM_guardedalloc.h"
  32
  33
  34 #include "BLI_math.h"
  35 #include "BLI_quadric.h"
  36 #include "BLI_heap.h"
  37 #include "BLI_linklist.h"
  38 #include "BLI_alloca.h"
  39 #include "BLI_memarena.h"
  40 #include "BLI_edgehash.h"
  41 #include "BLI_polyfill2d.h"
  42 #include "BLI_polyfill2d_beautify.h"
  43 #include "BLI_utildefines_stack.h"
  44
  45
  46 #include "BKE_customdata.h"
  47
  48 #include "bmesh.h"
  49 #include "bmesh_decimate.h"  /* own include */
  50
  51 #include "../intern/bmesh_structure.h"
  52
  53 #define USE_SYMMETRY
  54 #ifdef USE_SYMMETRY
  55 #include "BLI_kdtree.h"
  56 #endif
  57
  58 /* defines for testing */
  59 #define USE_CUSTOMDATA
  60 #define USE_TRIANGULATE
  61 #define USE_VERT_NORMAL_INTERP  /* has the advantage that flipped faces don't mess up vertex normals */
  62
  63 /* if the cost from #BLI_quadric_evaluate is 'noise', fallback to topology */
  64 #define USE_TOPOLOGY_FALLBACK
  65 #ifdef  USE_TOPOLOGY_FALLBACK
  66 /* cost is calculated with double precision, it's ok to use a very small epsilon, see T48154. */
  67 #  define   TOPOLOGY_FALLBACK_EPS  1e-12f
  68 #endif
  69
  70 #define BOUNDARY_PRESERVE_WEIGHT 100.0f
  71 /* Uses double precision, impacts behavior on near-flat surfaces,
  72  * cane give issues with very small faces. 1e-2 is too big, see: T48154. */
  73 #define OPTIMIZE_EPS 1e-8
  74 #define COST_INVALID FLT_MAX
  75
  76 typedef enum CD_UseFlag {
  77         CD_DO_VERT = (1 << 0),
  78         CD_DO_EDGE = (1 << 1),
  79         CD_DO_LOOP = (1 << 2)
  80 } CD_UseFlag;
  81
  82
  83 /* BMesh Helper Functions
  84  * ********************** */
  85
  86 /**
  87  * \param vquadrics must be calloc'd
  88  */
  89 static void bm_decim_build_quadrics(BMesh *bm, Quadric *vquadrics)
  90 {
  91         BMIter iter;
  92         BMFace *f;
  93         BMEdge *e;
  94
  95         BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
  96                 BMLoop *l_first;
  97                 BMLoop *l_iter;
  98
  99                 float center[3];
 100                 double plane_db[4];
 101                 Quadric q;
 102
 103                 BM_face_calc_center_mean(f, center);
 104                 copy_v3db_v3fl(plane_db, f->no);
 105                 plane_db[3] = -dot_v3db_v3fl(plane_db, center);
 106
 107                 BLI_quadric_from_plane(&q, plane_db);
 108
 109                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
 110                 do {
 111                         BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(l_iter->v)], &q);
 112                 } while ((l_iter = l_iter->next) != l_first);
 113         }
 114
 115         /* boundary edges */
 116         BM_ITER_MESH (e, &iter, bm, BM_EDGES_OF_MESH) {
 117                 if (UNLIKELY(BM_edge_is_boundary(e))) {
 118                         float edge_vector[3];
 119                         float edge_plane[3];
 120                         double edge_plane_db[4];
 121                         sub_v3_v3v3(edge_vector, e->v2->co, e->v1->co);
 122                         f = e->l->f;
 123
 124                         cross_v3_v3v3(edge_plane, edge_vector, f->no);
 125                         copy_v3db_v3fl(edge_plane_db, edge_plane);
 126
 127                         if (normalize_v3_d(edge_plane_db) > (double)FLT_EPSILON) {
 128                                 Quadric q;
 129                                 float center[3];
 130
 131                                 mid_v3_v3v3(center, e->v1->co, e->v2->co);
 132
 133                                 edge_plane_db[3] = -dot_v3db_v3fl(edge_plane_db, center);
 134                                 BLI_quadric_from_plane(&q, edge_plane_db);
 135                                 BLI_quadric_mul(&q, BOUNDARY_PRESERVE_WEIGHT);
 136
 137                                 BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(e->v1)], &q);
 138                                 BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(e->v2)], &q);
 139                         }
 140                 }
 141         }
 142 }
 143
 144
 145 static void bm_decim_calc_target_co_db(
 146         BMEdge *e, double optimize_co[3],
 147         const Quadric *vquadrics)
 148 {
 149         /* compute an edge contraction target for edge 'e'
 150          * this is computed by summing it's vertices quadrics and
 151          * optimizing the result. */
 152         Quadric q;
 153
 154         BLI_quadric_add_qu_ququ(
 155                 &q,
 156                 &vquadrics[BM_elem_index_get(e->v1)],
 157                 &vquadrics[BM_elem_index_get(e->v2)]);
 158
 159         if (BLI_quadric_optimize(&q, optimize_co, OPTIMIZE_EPS)) {
 160                 return;  /* all is good */
 161         }
 162         else {
 163                 optimize_co[0] = 0.5 * ((double)e->v1->co[0] + (double)e->v2->co[0]);
 164                 optimize_co[1] = 0.5 * ((double)e->v1->co[1] + (double)e->v2->co[1]);
 165                 optimize_co[2] = 0.5 * ((double)e->v1->co[2] + (double)e->v2->co[2]);
 166         }
 167 }
 168
 169 static void bm_decim_calc_target_co_fl(
 170         BMEdge *e, float optimize_co[3],
 171         const Quadric *vquadrics)
 172 {
 173         double optimize_co_db[3];
 174         bm_decim_calc_target_co_db(e, optimize_co_db, vquadrics);
 175         copy_v3fl_v3db(optimize_co, optimize_co_db);
 176 }
 177
 178
 179 static bool bm_edge_collapse_is_degenerate_flip(BMEdge *e, const float optimize_co[3])
 180 {
 181         BMIter liter;
 182         BMLoop *l;
 183         uint i;
 184
 185         for (i = 0; i < 2; i++) {
 186                 /* loop over both verts */
 187                 BMVert *v = *((&e->v1) + i);
 188
 189                 BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
 190                         if (l->e != e && l->prev->e != e) {
 191                                 const float *co_prev = l->prev->v->co;
 192                                 const float *co_next = l->next->v->co;
 193                                 float cross_exist[3];
 194                                 float cross_optim[3];
 195
 196 #if 1
 197                                 float vec_other[3];  /* line between the two outer verts, re-use for both cross products */
 198                                 float vec_exist[3];  /* before collapse */
 199                                 float vec_optim[3];  /* after collapse */
 200
 201                                 sub_v3_v3v3(vec_other, co_prev, co_next);
 202                                 sub_v3_v3v3(vec_exist, co_prev, v->co);
 203                                 sub_v3_v3v3(vec_optim, co_prev, optimize_co);
 204
 205                                 cross_v3_v3v3(cross_exist, vec_other, vec_exist);
 206                                 cross_v3_v3v3(cross_optim, vec_other, vec_optim);
 207
 208                                 /* avoid normalize */
 209                                 if (dot_v3v3(cross_exist, cross_optim) <=
 210                                     (len_squared_v3(cross_exist) + len_squared_v3(cross_optim)) * 0.01f)
 211                                 {
 212                                         return true;
 213                                 }
 214 #else
 215                                 normal_tri_v3(cross_exist, v->co,       co_prev, co_next);
 216                                 normal_tri_v3(cross_optim, optimize_co, co_prev, co_next);
 217
 218                                 /* use a small value rather then zero so we don't flip a face in multiple steps
 219                                  * (first making it zero area, then flipping again) */
 220                                 if (dot_v3v3(cross_exist, cross_optim) <= FLT_EPSILON) {
 221                                         //printf("no flip\n");
 222                                         return true;
 223                                 }
 224 #endif
 225
 226                         }
 227                 }
 228         }
 229
 230         return false;
 231 }
 232
 233 #ifdef USE_TOPOLOGY_FALLBACK
 234 /**
 235  * when the cost is so small that its not useful (flat surfaces),
 236  * fallback to using a 'topology' cost.
 237  *
 238  * This avoids cases where a flat (or near flat) areas get very un-even geometry.
 239  */
 240 static float bm_decim_build_edge_cost_single_squared__topology(BMEdge *e)
 241 {
 242         return fabsf(dot_v3v3(e->v1->no, e->v2->no)) / min_ff(-len_squared_v3v3(e->v1->co, e->v2->co), -FLT_EPSILON);
 243 }
 244 static float bm_decim_build_edge_cost_single__topology(BMEdge *e)
 245 {
 246         return fabsf(dot_v3v3(e->v1->no, e->v2->no)) / min_ff(-len_v3v3(e->v1->co, e->v2->co), -FLT_EPSILON);
 247 }
 248
 249 #endif  /* USE_TOPOLOGY_FALLBACK */
 250
 251 static void bm_decim_build_edge_cost_single(
 252         BMEdge *e,
 253         const Quadric *vquadrics,
 254         const float *vweights, const float vweight_factor,
 255         Heap *eheap, HeapNode **eheap_table)
 256 {
 257         float cost;
 258
 259         if (UNLIKELY(vweights &&
 260                      ((vweights[BM_elem_index_get(e->v1)] == 0.0f) ||
 261                       (vweights[BM_elem_index_get(e->v2)] == 0.0f))))
 262         {
 263                 goto clear;
 264         }
 265
 266         /* check we can collapse, some edges we better not touch */
 267         if (BM_edge_is_boundary(e)) {
 268                 if (e->l->f->len == 3) {
 269                         /* pass */
 270                 }
 271                 else {
 272                         /* only collapse tri's */
 273                         goto clear;
 274                 }
 275         }
 276         else if (BM_edge_is_manifold(e)) {
 277                 if ((e->l->f->len == 3) && (e->l->radial_next->f->len == 3)) {
 278                         /* pass */
 279                 }
 280                 else {
 281                         /* only collapse tri's */
 282                         goto clear;
 283                 }
 284         }
 285         else {
 286                 goto clear;
 287         }
 288         /* end sanity check */
 289
 290         {
 291                 double optimize_co[3];
 292                 bm_decim_calc_target_co_db(e, optimize_co, vquadrics);
 293
 294                 const Quadric *q1, *q2;
 295                 q1 = &vquadrics[BM_elem_index_get(e->v1)];
 296                 q2 = &vquadrics[BM_elem_index_get(e->v2)];
 297
 298                 cost = (BLI_quadric_evaluate(q1, optimize_co) +
 299                         BLI_quadric_evaluate(q2, optimize_co));
 300         }
 301
 302         /* note, 'cost' shouldn't be negative but happens sometimes with small values.
 303          * this can cause faces that make up a flat surface to over-collapse, see [#37121] */
 304         cost = fabsf(cost);
 305
 306 #ifdef USE_TOPOLOGY_FALLBACK
 307         if (UNLIKELY(cost < TOPOLOGY_FALLBACK_EPS)) {
 308                 /* subtract existing cost to further differentiate edges from one another
 309                  *
 310                  * keep topology cost below 0.0 so their values don't interfere with quadric cost,
 311                  * (and they get handled first).
 312                  * */
 313                 if (vweights == NULL) {
 314                         cost = bm_decim_build_edge_cost_single_squared__topology(e) - cost;
 315                 }
 316                 else {
 317                         /* with weights we need to use the real length so we can scale them properly */
 318                         const float e_weight = (vweights[BM_elem_index_get(e->v1)] +
 319                                                 vweights[BM_elem_index_get(e->v2)]);
 320                         cost = bm_decim_build_edge_cost_single__topology(e) - cost;
 321                         /* note, this is rather arbitrary max weight is 2 here,
 322                          * allow for skipping edges 4x the length, based on weights */
 323                         if (e_weight) {
 324                                 cost *= 1.0f + (e_weight * vweight_factor);
 325                         }
 326
 327                         BLI_assert(cost <= 0.0f);
 328                 }
 329         }
 330         else
 331 #endif
 332         if (vweights) {
 333                 const float e_weight = 2.0f - (vweights[BM_elem_index_get(e->v1)] +
 334                                                vweights[BM_elem_index_get(e->v2)]);
 335                 if (e_weight) {
 336                         cost += (BM_edge_calc_length(e) * ((e_weight * vweight_factor)));
 337                 }
 338         }
 339
 340         BLI_heap_insert_or_update(eheap, &eheap_table[BM_elem_index_get(e)], cost, e);
 341         return;
 342
 343 clear:
 344         if (eheap_table[BM_elem_index_get(e)]) {
 345                 BLI_heap_remove(eheap, eheap_table[BM_elem_index_get(e)]);
 346         }
 347         eheap_table[BM_elem_index_get(e)] = NULL;
 348 }
 349
 350
 351 /* use this for degenerate cases - add back to the heap with an invalid cost,
 352  * this way it may be calculated again if surrounding geometry changes */
 353 static void bm_decim_invalid_edge_cost_single(
 354         BMEdge *e,
 355         Heap *eheap, HeapNode **eheap_table)
 356 {
 357         BLI_assert(eheap_table[BM_elem_index_get(e)] == NULL);
 358         eheap_table[BM_elem_index_get(e)] = BLI_heap_insert(eheap, COST_INVALID, e);
 359 }
 360
 361 static void bm_decim_build_edge_cost(
 362         BMesh *bm,
 363         const Quadric *vquadrics,
 364         const float *vweights, const float vweight_factor,
 365         Heap *eheap, HeapNode **eheap_table)
 366 {
 367         BMIter iter;
 368         BMEdge *e;
 369         uint i;
 370
 371         BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
 372                 eheap_table[i] = NULL;  /* keep sanity check happy */
 373                 bm_decim_build_edge_cost_single(e, vquadrics, vweights, vweight_factor, eheap, eheap_table);
 374         }
 375 }
 376
 377 #ifdef USE_SYMMETRY
 378
 379 struct KD_Symmetry_Data {
 380         /* pre-flipped coords */
 381         float e_v1_co[3], e_v2_co[3];
 382         /* Use to compare the correct endpoints incase v1/v2 are swapped */
 383         float e_dir[3];
 384
 385         int e_found_index;
 386
 387         /* same for all */
 388         BMEdge **etable;
 389         float limit_sq;
 390 };
 391
 392 static bool bm_edge_symmetry_check_cb(void *user_data, int index, const float UNUSED(co[3]), float UNUSED(dist_sq))
 393 {
 394         struct KD_Symmetry_Data *sym_data = user_data;
 395         BMEdge *e_other = sym_data->etable[index];
 396         float e_other_dir[3];
 397
 398         sub_v3_v3v3(e_other_dir, e_other->v2->co, e_other->v1->co);
 399
 400         if (dot_v3v3(e_other_dir, sym_data->e_dir) > 0.0f) {
 401                 if ((len_squared_v3v3(sym_data->e_v1_co, e_other->v1->co) > sym_data->limit_sq) ||
 402                     (len_squared_v3v3(sym_data->e_v2_co, e_other->v2->co) > sym_data->limit_sq))
 403                 {
 404                         return true;
 405                 }
 406         }
 407         else {
 408                 if ((len_squared_v3v3(sym_data->e_v1_co, e_other->v2->co) > sym_data->limit_sq) ||
 409                     (len_squared_v3v3(sym_data->e_v2_co, e_other->v1->co) > sym_data->limit_sq))
 410                 {
 411                         return true;
 412                 }
 413         }
 414
 415         /* exit on first-hit, this is OK since the search range is very small */
 416         sym_data->e_found_index = index;
 417         return false;
 418 }
 419
 420 static int *bm_edge_symmetry_map(BMesh *bm, uint symmetry_axis, float limit)
 421 {
 422         struct KD_Symmetry_Data sym_data;
 423         BMIter iter;
 424         BMEdge *e, **etable;
 425         uint i;
 426         int *edge_symmetry_map;
 427         const float limit_sq = SQUARE(limit);
 428         KDTree *tree;
 429
 430         tree = BLI_kdtree_new(bm->totedge);
 431
 432         etable = MEM_mallocN(sizeof(*etable) * bm->totedge, __func__);
 433         edge_symmetry_map = MEM_mallocN(sizeof(*edge_symmetry_map) * bm->totedge, __func__);
 434
 435         BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
 436                 float co[3];
 437                 mid_v3_v3v3(co, e->v1->co, e->v2->co);
 438                 BLI_kdtree_insert(tree, i, co);
 439                 etable[i] = e;
 440                 edge_symmetry_map[i] = -1;
 441         }
 442
 443         BLI_kdtree_balance(tree);
 444
 445         sym_data.etable = etable;
 446         sym_data.limit_sq = limit_sq;
 447
 448         BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
 449                 if (edge_symmetry_map[i] == -1) {
 450                         float co[3];
 451                         mid_v3_v3v3(co, e->v1->co, e->v2->co);
 452                         co[symmetry_axis] *= -1.0f;
 453
 454                         copy_v3_v3(sym_data.e_v1_co, e->v1->co);
 455                         copy_v3_v3(sym_data.e_v2_co, e->v2->co);
 456                         sym_data.e_v1_co[symmetry_axis] *= -1.0f;
 457                         sym_data.e_v2_co[symmetry_axis] *= -1.0f;
 458                         sub_v3_v3v3(sym_data.e_dir, sym_data.e_v2_co, sym_data.e_v1_co);
 459                         sym_data.e_found_index = -1;
 460
 461                         BLI_kdtree_range_search_cb(tree, co, limit, bm_edge_symmetry_check_cb, &sym_data);
 462
 463                         if (sym_data.e_found_index != -1) {
 464                                 const int i_other = sym_data.e_found_index;
 465                                 edge_symmetry_map[i] = i_other;
 466                                 edge_symmetry_map[i_other] = i;
 467                         }
 468                 }
 469         }
 470
 471         MEM_freeN(etable);
 472         BLI_kdtree_free(tree);
 473
 474         return edge_symmetry_map;
 475 }
 476 #endif  /* USE_SYMMETRY */
 477
 478 #ifdef USE_TRIANGULATE
 479 /* Temp Triangulation
 480  * ****************** */
 481
 482 /**
 483  * To keep things simple we can only collapse edges on triangulated data
 484  * (limitation with edge collapse and error calculation functions).
 485  *
 486  * But to avoid annoying users by only giving triangle results, we can
 487  * triangulate, keeping a reference between the faces, then join after
 488  * if the edges don't collapse, this will also allow more choices when
 489  * collapsing edges so even has some advantage over decimating quads
 490  * directly.
 491  *
 492  * \return true if any faces were triangulated.
 493  */
 494 static bool bm_face_triangulate(
 495         BMesh *bm, BMFace *f_base, LinkNode **r_faces_double, int *r_edges_tri_tot,
 496
 497         MemArena *pf_arena,
 498         /* use for MOD_TRIANGULATE_NGON_BEAUTY only! */
 499         struct Heap *pf_heap)
 500 {
 501         const int f_base_len = f_base->len;
 502         int faces_array_tot = f_base_len - 3;
 503         int edges_array_tot = f_base_len - 3;
 504         BMFace  **faces_array = BLI_array_alloca(faces_array, faces_array_tot);
 505         BMEdge  **edges_array = BLI_array_alloca(edges_array, edges_array_tot);
 506         const int quad_method = 0, ngon_method = 0;  /* beauty */
 507
 508         bool has_cut = false;
 509
 510         const int f_index = BM_elem_index_get(f_base);
 511
 512         BM_face_triangulate(
 513                 bm, f_base,
 514                 faces_array, &faces_array_tot,
 515                 edges_array, &edges_array_tot,
 516                 r_faces_double,
 517                 quad_method, ngon_method, false,
 518                 pf_arena, pf_heap);
 519
 520         for (int i = 0; i < edges_array_tot; i++) {
 521                 BMLoop *l_iter, *l_first;
 522                 l_iter = l_first = edges_array[i]->l;
 523                 do {
 524                         BM_elem_index_set(l_iter, f_index);  /* set_dirty */
 525                         has_cut = true;
 526                 } while ((l_iter = l_iter->radial_next) != l_first);
 527         }
 528
 529         for (int i = 0; i < faces_array_tot; i++) {
 530                 BM_face_normal_update(faces_array[i]);
 531         }
 532
 533         *r_edges_tri_tot += edges_array_tot;
 534
 535         return has_cut;
 536 }
 537
 538
 539 static bool bm_decim_triangulate_begin(BMesh *bm, int *r_edges_tri_tot)
 540 {
 541         BMIter iter;
 542         BMFace *f;
 543         bool has_quad = false;
 544         bool has_ngon = false;
 545         bool has_cut = false;
 546
 547         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
 548
 549         /* first clear loop index values */
 550         BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
 551                 BMLoop *l_iter;
 552                 BMLoop *l_first;
 553
 554                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
 555                 do {
 556                         BM_elem_index_set(l_iter, -1);  /* set_dirty */
 557                 } while ((l_iter = l_iter->next) != l_first);
 558
 559                 has_quad |= (f->len > 3);
 560                 has_ngon |= (f->len > 4);
 561         }
 562
 563         bm->elem_index_dirty |= BM_LOOP;
 564
 565         {
 566                 MemArena *pf_arena;
 567                 Heap *pf_heap;
 568
 569                 LinkNode *faces_double = NULL;
 570
 571                 if (has_ngon) {
 572                         pf_arena = BLI_memarena_new(BLI_POLYFILL_ARENA_SIZE, __func__);
 573                         pf_heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
 574                 }
 575                 else {
 576                         pf_arena = NULL;
 577                         pf_heap = NULL;
 578                 }
 579
 580                 /* adding new faces as we loop over faces
 581                  * is normally best avoided, however in this case its not so bad because any face touched twice
 582                  * will already be triangulated*/
 583                 BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
 584                         if (f->len > 3) {
 585                                 has_cut |= bm_face_triangulate(
 586                                         bm, f, &faces_double,
 587                                         r_edges_tri_tot,
 588
 589                                         pf_arena, pf_heap);
 590                         }
 591                 }
 592
 593                 while (faces_double) {
 594                         LinkNode *next = faces_double->next;
 595                         BM_face_kill(bm, faces_double->link);
 596                         MEM_freeN(faces_double);
 597                         faces_double = next;
 598                 }
 599
 600                 if (has_ngon) {
 601                         BLI_memarena_free(pf_arena);
 602                         BLI_heap_free(pf_heap, NULL);
 603                 }
 604
 605                 BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
 606
 607                 if (has_cut) {
 608                         /* now triangulation is done we need to correct index values */
 609                         BM_mesh_elem_index_ensure(bm, BM_EDGE | BM_FACE);
 610                 }
 611         }
 612
 613         return has_cut;
 614 }
 615
 616
 617 static void bm_decim_triangulate_end(BMesh *bm, const int edges_tri_tot)
 618 {
 619         /* decimation finished, now re-join */
 620         BMIter iter;
 621         BMEdge *e;
 622
 623         /* we need to collect before merging for ngons since the loops indices will be lost */
 624         BMEdge **edges_tri = MEM_mallocN(MIN2(edges_tri_tot, bm->totedge) * sizeof(*edges_tri), __func__);
 625         STACK_DECLARE(edges_tri);
 626
 627         STACK_INIT(edges_tri, MIN2(edges_tri_tot, bm->totedge));
 628
 629         /* boundary edges */
 630         BM_ITER_MESH (e, &iter, bm, BM_EDGES_OF_MESH) {
 631                 BMLoop *l_a, *l_b;
 632                 if (BM_edge_loop_pair(e, &l_a, &l_b)) {
 633                         const int l_a_index = BM_elem_index_get(l_a);
 634                         if (l_a_index != -1) {
 635                                 const int l_b_index = BM_elem_index_get(l_b);
 636                                 if (l_a_index == l_b_index) {
 637                                         if (l_a->v != l_b->v) {  /* if this is the case, faces have become flipped */
 638                                                 /* check we are not making a degenerate quad */
 639
 640 #define CAN_LOOP_MERGE(l) \
 641         (BM_loop_is_manifold(l) && \
 642          ((l)->v != (l)->radial_next->v) && \
 643          (l_a_index == BM_elem_index_get(l)) && \
 644          (l_a_index == BM_elem_index_get((l)->radial_next)))
 645
 646                                                 if ((l_a->f->len == 3 && l_b->f->len == 3) &&
 647                                                     (!CAN_LOOP_MERGE(l_a->next)) &&
 648                                                     (!CAN_LOOP_MERGE(l_a->prev)) &&
 649                                                     (!CAN_LOOP_MERGE(l_b->next)) &&
 650                                                     (!CAN_LOOP_MERGE(l_b->prev)))
 651                                                 {
 652                                                         BMVert *vquad[4] = {
 653                                                                 e->v1,
 654                                                                 BM_vert_in_edge(e, l_a->next->v) ? l_a->prev->v : l_a->next->v,
 655                                                                 e->v2,
 656                                                                 BM_vert_in_edge(e, l_b->next->v) ? l_b->prev->v : l_b->next->v,
 657                                                         };
 658
 659                                                         BLI_assert(ELEM(vquad[0], vquad[1], vquad[2], vquad[3]) == false);
 660                                                         BLI_assert(ELEM(vquad[1], vquad[0], vquad[2], vquad[3]) == false);
 661                                                         BLI_assert(ELEM(vquad[2], vquad[1], vquad[0], vquad[3]) == false);
 662                                                         BLI_assert(ELEM(vquad[3], vquad[1], vquad[2], vquad[0]) == false);
 663
 664                                                         if (!is_quad_convex_v3(vquad[0]->co, vquad[1]->co, vquad[2]->co, vquad[3]->co)) {
 665                                                                 continue;
 666                                                         }
 667                                                 }
 668 #undef CAN_LOOP_MERGE
 669
 670                                                 /* highly unlikely to fail, but prevents possible double-ups */
 671                                                 STACK_PUSH(edges_tri, e);
 672                                         }
 673                                 }
 674                         }
 675                 }
 676         }
 677
 678         for (int i = 0; i < STACK_SIZE(edges_tri); i++) {
 679                 BMLoop *l_a, *l_b;
 680                 e = edges_tri[i];
 681                 if (BM_edge_loop_pair(e, &l_a, &l_b)) {
 682                         BMFace *f_array[2] = {l_a->f, l_b->f};
 683                         BM_faces_join(bm, f_array, 2, false);
 684                         if (e->l == NULL) {
 685                                 BM_edge_kill(bm, e);
 686                         }
 687                 }
 688         }
 689         MEM_freeN(edges_tri);
 690 }
 691
 692 #endif  /* USE_TRIANGULATE */
 693
 694 /* Edge Collapse Functions
 695  * *********************** */
 696
 697 #ifdef USE_CUSTOMDATA
 698
 699 /**
 700  * \param l: defines the vert to collapse into.
 701  */
 702 static void bm_edge_collapse_loop_customdata(
 703         BMesh *bm, BMLoop *l, BMVert *v_clear, BMVert *v_other,
 704         const float customdata_fac)
 705 {
 706         /* disable seam check - the seam check would have to be done per layer, its not really that important */
 707 //#define USE_SEAM
 708         /* these don't need to be updated, since they will get removed when the edge collapses */
 709         BMLoop *l_clear, *l_other;
 710         const bool is_manifold = BM_edge_is_manifold(l->e);
 711         int side;
 712
 713         /* first find the loop of 'v_other' thats attached to the face of 'l' */
 714         if (l->v == v_clear) {
 715                 l_clear = l;
 716                 l_other = l->next;
 717         }
 718         else {
 719                 l_clear = l->next;
 720                 l_other = l;
 721         }
 722
 723         BLI_assert(l_clear->v == v_clear);
 724         BLI_assert(l_other->v == v_other);
 725         (void)v_other;  /* quiet warnings for release */
 726
 727         /* now we have both corners of the face 'l->f' */
 728         for (side = 0; side < 2; side++) {
 729 #ifdef USE_SEAM
 730                 bool is_seam = false;
 731 #endif
 732                 void *src[2];
 733                 BMFace *f_exit = is_manifold ? l->radial_next->f : NULL;
 734                 BMEdge *e_prev = l->e;
 735                 BMLoop *l_first;
 736                 BMLoop *l_iter;
 737                 float w[2];
 738
 739                 if (side == 0) {
 740                         l_iter = l_first = l_clear;
 741                         src[0] = l_clear->head.data;
 742                         src[1] = l_other->head.data;
 743
 744                         w[0] = customdata_fac;
 745                         w[1] = 1.0f - customdata_fac;
 746                 }
 747                 else {
 748                         l_iter = l_first = l_other;
 749                         src[0] = l_other->head.data;
 750                         src[1] = l_clear->head.data;
 751
 752                         w[0] = 1.0f - customdata_fac;
 753                         w[1] = customdata_fac;
 754                 }
 755
 756                 // print_v2("weights", w);
 757
 758                 /* WATCH IT! - should NOT reference (_clear or _other) vars for this while loop */
 759
 760                 /* walk around the fan using 'e_prev' */
 761                 while (((l_iter = BM_vert_step_fan_loop(l_iter, &e_prev)) != l_first) && (l_iter != NULL)) {
 762                         int i;
 763                         /* quit once we hit the opposite face, if we have one */
 764                         if (f_exit && UNLIKELY(f_exit == l_iter->f)) {
 765                                 break;
 766                         }
 767
 768 #ifdef USE_SEAM
 769                         /* break out unless we find a match */
 770                         is_seam = true;
 771 #endif
 772
 773                         /* ok. we have a loop. now be smart with it! */
 774                         for (i = 0; i < bm->ldata.totlayer; i++) {
 775                                 if (CustomData_layer_has_math(&bm->ldata, i)) {
 776                                         const int offset = bm->ldata.layers[i].offset;
 777                                         const int type = bm->ldata.layers[i].type;
 778                                         const void *cd_src[2] = {
 779                                             POINTER_OFFSET(src[0], offset),
 780                                             POINTER_OFFSET(src[1], offset),
 781                                         };
 782                                         void *cd_iter = POINTER_OFFSET(l_iter->head.data, offset);
 783
 784                                         /* detect seams */
 785                                         if (CustomData_data_equals(type, cd_src[0], cd_iter)) {
 786                                                 CustomData_bmesh_interp_n(
 787                                                         &bm->ldata, cd_src, w, NULL, ARRAY_SIZE(cd_src),
 788                                                         POINTER_OFFSET(l_iter->head.data, offset), i);
 789 #ifdef USE_SEAM
 790                                                 is_seam = false;
 791 #endif
 792                                         }
 793                                 }
 794                         }
 795
 796 #ifdef USE_SEAM
 797                         if (is_seam) {
 798                                 break;
 799                         }
 800 #endif
 801                 }
 802         }
 803
 804 //#undef USE_SEAM
 805
 806 }
 807 #endif  /* USE_CUSTOMDATA */
 808
 809 /**
 810  * Check if the collapse will result in a degenerate mesh,
 811  * that is - duplicate edges or faces.
 812  *
 813  * This situation could be checked for when calculating collapse cost
 814  * however its quite slow and a degenerate collapse could eventuate
 815  * after the cost is calculated, so instead, check just before collapsing.
 816  */
 817
 818 static void bm_edge_tag_enable(BMEdge *e)
 819 {
 820         BM_elem_flag_enable(e->v1, BM_ELEM_TAG);
 821         BM_elem_flag_enable(e->v2, BM_ELEM_TAG);
 822         if (e->l) {
 823                 BM_elem_flag_enable(e->l->f, BM_ELEM_TAG);
 824                 if (e->l != e->l->radial_next) {
 825                         BM_elem_flag_enable(e->l->radial_next->f, BM_ELEM_TAG);
 826                 }
 827         }
 828 }
 829
 830 static void bm_edge_tag_disable(BMEdge *e)
 831 {
 832         BM_elem_flag_disable(e->v1, BM_ELEM_TAG);
 833         BM_elem_flag_disable(e->v2, BM_ELEM_TAG);
 834         if (e->l) {
 835                 BM_elem_flag_disable(e->l->f, BM_ELEM_TAG);
 836                 if (e->l != e->l->radial_next) {
 837                         BM_elem_flag_disable(e->l->radial_next->f, BM_ELEM_TAG);
 838                 }
 839         }
 840 }
 841
 842 static bool bm_edge_tag_test(BMEdge *e)
 843 {
 844         /* is the edge or one of its faces tagged? */
 845         return (BM_elem_flag_test(e->v1, BM_ELEM_TAG) ||
 846                 BM_elem_flag_test(e->v2, BM_ELEM_TAG) ||
 847                 (e->l && (BM_elem_flag_test(e->l->f, BM_ELEM_TAG) ||
 848                           (e->l != e->l->radial_next &&
 849                           BM_elem_flag_test(e->l->radial_next->f, BM_ELEM_TAG))))
 850                 );
 851 }
 852
 853 /* takes the edges loop */
 854 BLI_INLINE int bm_edge_is_manifold_or_boundary(BMLoop *l)
 855 {
 856 #if 0
 857         /* less optimized version of check below */
 858         return (BM_edge_is_manifold(l->e) || BM_edge_is_boundary(l->e);
 859 #else
 860         /* if the edge is a boundary it points to its self, else this must be a manifold */
 861         return LIKELY(l) && LIKELY(l->radial_next->radial_next == l);
 862 #endif
 863 }
 864
 865 static bool bm_edge_collapse_is_degenerate_topology(BMEdge *e_first)
 866 {
 867         /* simply check that there is no overlap between faces and edges of each vert,
 868          * (excluding the 2 faces attached to 'e' and 'e' its self) */
 869
 870         BMEdge *e_iter;
 871
 872         /* clear flags on both disks */
 873         e_iter = e_first;
 874         do {
 875                 if (!bm_edge_is_manifold_or_boundary(e_iter->l)) {
 876                         return true;
 877                 }
 878                 bm_edge_tag_disable(e_iter);
 879         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v1)) != e_first);
 880
 881         e_iter = e_first;
 882         do {
 883                 if (!bm_edge_is_manifold_or_boundary(e_iter->l)) {
 884                         return true;
 885                 }
 886                 bm_edge_tag_disable(e_iter);
 887         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v2)) != e_first);
 888
 889         /* now enable one side... */
 890         e_iter = e_first;
 891         do {
 892                 bm_edge_tag_enable(e_iter);
 893         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v1)) != e_first);
 894
 895         /* ... except for the edge we will collapse, we know thats shared,
 896          * disable this to avoid false positive. We could be smart and never enable these
 897          * face/edge tags in the first place but easier to do this */
 898         // bm_edge_tag_disable(e_first);
 899         /* do inline... */
 900         {
 901 #if 0
 902                 BMIter iter;
 903                 BMIter liter;
 904                 BMLoop *l;
 905                 BMVert *v;
 906                 BM_ITER_ELEM (l, &liter, e_first, BM_LOOPS_OF_EDGE) {
 907                         BM_elem_flag_disable(l->f, BM_ELEM_TAG);
 908                         BM_ITER_ELEM (v, &iter, l->f, BM_VERTS_OF_FACE) {
 909                                 BM_elem_flag_disable(v, BM_ELEM_TAG);
 910                         }
 911                 }
 912 #else
 913                 /* we know each face is a triangle, no looping/iterators needed here */
 914
 915                 BMLoop *l_radial;
 916                 BMLoop *l_face;
 917
 918                 l_radial = e_first->l;
 919                 l_face = l_radial;
 920                 BLI_assert(l_face->f->len == 3);
 921                 BM_elem_flag_disable(l_face->f, BM_ELEM_TAG);
 922                 BM_elem_flag_disable((l_face = l_radial)->v,     BM_ELEM_TAG);
 923                 BM_elem_flag_disable((l_face = l_face->next)->v, BM_ELEM_TAG);
 924                 BM_elem_flag_disable((         l_face->next)->v, BM_ELEM_TAG);
 925                 l_face = l_radial->radial_next;
 926                 if (l_radial != l_face) {
 927                         BLI_assert(l_face->f->len == 3);
 928                         BM_elem_flag_disable(l_face->f, BM_ELEM_TAG);
 929                         BM_elem_flag_disable((l_face = l_radial->radial_next)->v, BM_ELEM_TAG);
 930                         BM_elem_flag_disable((l_face = l_face->next)->v,          BM_ELEM_TAG);
 931                         BM_elem_flag_disable((         l_face->next)->v,          BM_ELEM_TAG);
 932                 }
 933 #endif
 934         }
 935
 936         /* and check for overlap */
 937         e_iter = e_first;
 938         do {
 939                 if (bm_edge_tag_test(e_iter)) {
 940                         return true;
 941                 }
 942         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v2)) != e_first);
 943
 944         return false;
 945 }
 946
 947 /**
 948  * special, highly limited edge collapse function
 949  * intended for speed over flexibility.
 950  * can only collapse edges connected to (1, 2) tris.
 951  *
 952  * Important - dont add vert/edge/face data on collapsing!
 953  *
 954  * \param r_e_clear_other: Let caller know what edges we remove besides \a e_clear
 955  * \param customdata_flag: Merge factor, scales from 0 - 1 ('v_clear' -> 'v_other')
 956  */
 957 static bool bm_edge_collapse(
 958         BMesh *bm, BMEdge *e_clear, BMVert *v_clear, int r_e_clear_other[2],
 959 #ifdef USE_SYMMETRY
 960         int *edge_symmetry_map,
 961 #endif
 962 #ifdef USE_CUSTOMDATA
 963         const CD_UseFlag customdata_flag,
 964         const float customdata_fac
 965 #else
 966         const CD_UseFlag UNUSED(customdata_flag),
 967         const float UNUSED(customdata_fac)
 968 #endif
 969         )
 970 {
 971         BMVert *v_other;
 972
 973         v_other = BM_edge_other_vert(e_clear, v_clear);
 974         BLI_assert(v_other != NULL);
 975
 976         if (BM_edge_is_manifold(e_clear)) {
 977                 BMLoop *l_a, *l_b;
 978                 BMEdge *e_a_other[2], *e_b_other[2];
 979                 bool ok;
 980
 981                 ok = BM_edge_loop_pair(e_clear, &l_a, &l_b);
 982
 983                 BLI_assert(ok == true);
 984                 BLI_assert(l_a->f->len == 3);
 985                 BLI_assert(l_b->f->len == 3);
 986                 UNUSED_VARS_NDEBUG(ok);
 987
 988                 /* keep 'v_clear' 0th */
 989                 if (BM_vert_in_edge(l_a->prev->e, v_clear)) {
 990                         e_a_other[0] = l_a->prev->e;
 991                         e_a_other[1] = l_a->next->e;
 992                 }
 993                 else {
 994                         e_a_other[1] = l_a->prev->e;
 995                         e_a_other[0] = l_a->next->e;
 996                 }
 997
 998                 if (BM_vert_in_edge(l_b->prev->e, v_clear)) {
 999                         e_b_other[0] = l_b->prev->e;
1000                         e_b_other[1] = l_b->next->e;
1001                 }
1002                 else {
1003                         e_b_other[1] = l_b->prev->e;
1004                         e_b_other[0] = l_b->next->e;
1005                 }
1006
1007                 /* we could assert this case, but better just bail out */
1008 #if 0
1009                 BLI_assert(e_a_other[0] != e_b_other[0]);
1010                 BLI_assert(e_a_other[0] != e_b_other[1]);
1011                 BLI_assert(e_b_other[0] != e_a_other[0]);
1012                 BLI_assert(e_b_other[0] != e_a_other[1]);
1013 #endif
1014                 /* not totally common but we want to avoid */
1015                 if (ELEM(e_a_other[0], e_b_other[0], e_b_other[1]) ||
1016                     ELEM(e_a_other[1], e_b_other[0], e_b_other[1]))
1017                 {
1018                         return false;
1019                 }
1020
1021                 BLI_assert(BM_edge_share_vert(e_a_other[0], e_b_other[0]));
1022                 BLI_assert(BM_edge_share_vert(e_a_other[1], e_b_other[1]));
1023
1024                 r_e_clear_other[0] = BM_elem_index_get(e_a_other[0]);
1025                 r_e_clear_other[1] = BM_elem_index_get(e_b_other[0]);
1026
1027 #ifdef USE_CUSTOMDATA
1028                 /* before killing, do customdata */
1029                 if (customdata_flag & CD_DO_VERT) {
1030                         BM_data_interp_from_verts(bm, v_other, v_clear, v_other, customdata_fac);
1031                 }
1032                 if (customdata_flag & CD_DO_EDGE) {
1033                         BM_data_interp_from_edges(bm, e_a_other[1], e_a_other[0], e_a_other[1], customdata_fac);
1034                         BM_data_interp_from_edges(bm, e_b_other[1], e_b_other[0], e_b_other[1], customdata_fac);
1035                 }
1036                 if (customdata_flag & CD_DO_LOOP) {
1037                         bm_edge_collapse_loop_customdata(bm, e_clear->l,              v_clear, v_other, customdata_fac);
1038                         bm_edge_collapse_loop_customdata(bm, e_clear->l->radial_next, v_clear, v_other, customdata_fac);
1039                 }
1040 #endif
1041
1042                 BM_edge_kill(bm, e_clear);
1043
1044                 v_other->head.hflag |= v_clear->head.hflag;
1045                 BM_vert_splice(bm, v_other, v_clear);
1046
1047                 e_a_other[1]->head.hflag |= e_a_other[0]->head.hflag;
1048                 e_b_other[1]->head.hflag |= e_b_other[0]->head.hflag;
1049                 BM_edge_splice(bm, e_a_other[1], e_a_other[0]);
1050                 BM_edge_splice(bm, e_b_other[1], e_b_other[0]);
1051
1052 #ifdef USE_SYMMETRY
1053                 /* update mirror map */
1054                 if (edge_symmetry_map) {
1055                         if (edge_symmetry_map[r_e_clear_other[0]] != -1) {
1056                                 edge_symmetry_map[edge_symmetry_map[r_e_clear_other[0]]] = BM_elem_index_get(e_a_other[1]);
1057                         }
1058                         if (edge_symmetry_map[r_e_clear_other[1]] != -1) {
1059                                 edge_symmetry_map[edge_symmetry_map[r_e_clear_other[1]]] = BM_elem_index_get(e_b_other[1]);
1060                         }
1061                 }
1062 #endif
1063
1064                 // BM_mesh_validate(bm);
1065
1066                 return true;
1067         }
1068         else if (BM_edge_is_boundary(e_clear)) {
1069                 /* same as above but only one triangle */
1070                 BMLoop *l_a;
1071                 BMEdge *e_a_other[2];
1072
1073                 l_a = e_clear->l;
1074
1075                 BLI_assert(l_a->f->len == 3);
1076
1077                 /* keep 'v_clear' 0th */
1078                 if (BM_vert_in_edge(l_a->prev->e, v_clear)) {
1079                         e_a_other[0] = l_a->prev->e;
1080                         e_a_other[1] = l_a->next->e;
1081                 }
1082                 else {
1083                         e_a_other[1] = l_a->prev->e;
1084                         e_a_other[0] = l_a->next->e;
1085                 }
1086
1087                 r_e_clear_other[0] = BM_elem_index_get(e_a_other[0]);
1088                 r_e_clear_other[1] = -1;
1089
1090 #ifdef USE_CUSTOMDATA
1091                 /* before killing, do customdata */
1092                 if (customdata_flag & CD_DO_VERT) {
1093                         BM_data_interp_from_verts(bm, v_other, v_clear, v_other, customdata_fac);
1094                 }
1095                 if (customdata_flag & CD_DO_EDGE) {
1096                         BM_data_interp_from_edges(bm, e_a_other[1], e_a_other[0], e_a_other[1], customdata_fac);
1097                 }
1098                 if (customdata_flag & CD_DO_LOOP) {
1099                         bm_edge_collapse_loop_customdata(bm, e_clear->l, v_clear, v_other, customdata_fac);
1100                 }
1101 #endif
1102
1103                 BM_edge_kill(bm, e_clear);
1104
1105                 v_other->head.hflag |= v_clear->head.hflag;
1106                 BM_vert_splice(bm, v_other, v_clear);
1107
1108                 e_a_other[1]->head.hflag |= e_a_other[0]->head.hflag;
1109                 BM_edge_splice(bm, e_a_other[1], e_a_other[0]);
1110
1111 #ifdef USE_SYMMETRY
1112                 /* update mirror map */
1113                 if (edge_symmetry_map) {
1114                         if (edge_symmetry_map[r_e_clear_other[0]] != -1) {
1115                                 edge_symmetry_map[edge_symmetry_map[r_e_clear_other[0]]] = BM_elem_index_get(e_a_other[1]);
1116                         }
1117                 }
1118 #endif
1119
1120                 // BM_mesh_validate(bm);
1121
1122                 return true;
1123         }
1124         else {
1125                 return false;
1126         }
1127 }
1128
1129
1130 /**
1131  * Collapse e the edge, removing e->v2
1132  *
1133  * \return true when the edge was collapsed.
1134  */
1135 static bool bm_decim_edge_collapse(
1136         BMesh *bm, BMEdge *e,
1137         Quadric *vquadrics,
1138         float *vweights, const float vweight_factor,
1139         Heap *eheap, HeapNode **eheap_table,
1140 #ifdef USE_SYMMETRY
1141         int *edge_symmetry_map,
1142 #endif
1143         const CD_UseFlag customdata_flag,
1144         float optimize_co[3], bool optimize_co_calc
1145         )
1146 {
1147         int e_clear_other[2];
1148         BMVert *v_other = e->v1;
1149         const int v_other_index = BM_elem_index_get(e->v1);
1150         const int v_clear_index = BM_elem_index_get(e->v2);  /* the vert is removed so only store the index */
1151         float customdata_fac;
1152
1153 #ifdef USE_VERT_NORMAL_INTERP
1154         float v_clear_no[3];
1155         copy_v3_v3(v_clear_no, e->v2->no);
1156 #endif
1157
1158         /* when false, use without degenerate checks */
1159         if (optimize_co_calc) {
1160                 /* disallow collapsing which results in degenerate cases */
1161                 if (UNLIKELY(bm_edge_collapse_is_degenerate_topology(e))) {
1162                         bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);  /* add back with a high cost */
1163                         return false;
1164                 }
1165
1166                 bm_decim_calc_target_co_fl(e, optimize_co, vquadrics);
1167
1168                 /* check if this would result in an overlapping face */
1169                 if (UNLIKELY(bm_edge_collapse_is_degenerate_flip(e, optimize_co))) {
1170                         bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);  /* add back with a high cost */
1171                         return false;
1172                 }
1173         }
1174
1175         /* use for customdata merging */
1176         if (LIKELY(compare_v3v3(e->v1->co, e->v2->co, FLT_EPSILON) == false)) {
1177                 customdata_fac = line_point_factor_v3(optimize_co, e->v1->co, e->v2->co);
1178 #if 0
1179                 /* simple test for stupid collapse */
1180                 if (customdata_fac < 0.0 - FLT_EPSILON || customdata_fac > 1.0f + FLT_EPSILON) {
1181                         return false;
1182                 }
1183 #endif
1184         }
1185         else {
1186                 /* avoid divide by zero */
1187                 customdata_fac = 0.5f;
1188         }
1189
1190         if (bm_edge_collapse(
1191                 bm, e, e->v2, e_clear_other,
1192 #ifdef USE_SYMMETRY
1193                 edge_symmetry_map,
1194 #endif
1195                 customdata_flag, customdata_fac))
1196         {
1197                 /* update collapse info */
1198                 int i;
1199
1200                 if (vweights) {
1201                         float v_other_weight = interpf(vweights[v_other_index], vweights[v_clear_index], customdata_fac);
1202                         CLAMP(v_other_weight, 0.0f, 1.0f);
1203                         vweights[v_other_index] = v_other_weight;
1204                 }
1205
1206                 e = NULL;  /* paranoid safety check */
1207
1208                 copy_v3_v3(v_other->co, optimize_co);
1209
1210                 /* remove eheap */
1211                 for (i = 0; i < 2; i++) {
1212                         /* highly unlikely 'eheap_table[ke_other[i]]' would be NULL, but do for sanity sake */
1213                         if ((e_clear_other[i] != -1) && (eheap_table[e_clear_other[i]] != NULL)) {
1214                                 BLI_heap_remove(eheap, eheap_table[e_clear_other[i]]);
1215                                 eheap_table[e_clear_other[i]] = NULL;
1216                         }
1217                 }
1218
1219                 /* update vertex quadric, add kept vertex from killed vertex */
1220                 BLI_quadric_add_qu_qu(&vquadrics[v_other_index], &vquadrics[v_clear_index]);
1221
1222                 /* update connected normals */
1223
1224                 /* in fact face normals are not used for progressive updates, no need to update them */
1225                 // BM_vert_normal_update_all(v);
1226 #ifdef USE_VERT_NORMAL_INTERP
1227                 interp_v3_v3v3(v_other->no, v_other->no, v_clear_no, customdata_fac);
1228                 normalize_v3(v_other->no);
1229 #else
1230                 BM_vert_normal_update(v_other);
1231 #endif
1232
1233
1234                 /* update error costs and the eheap */
1235                 if (LIKELY(v_other->e)) {
1236                         BMEdge *e_iter;
1237                         BMEdge *e_first;
1238                         e_iter = e_first = v_other->e;
1239                         do {
1240                                 BLI_assert(BM_edge_find_double(e_iter) == NULL);
1241                                 bm_decim_build_edge_cost_single(e_iter, vquadrics, vweights, vweight_factor, eheap, eheap_table);
1242                         } while ((e_iter = bmesh_disk_edge_next(e_iter, v_other)) != e_first);
1243                 }
1244
1245                 /* this block used to be disabled,
1246                  * but enable now since surrounding faces may have been
1247                  * set to COST_INVALID because of a face overlap that no longer occurs */
1248 #if 1
1249                 /* optional, update edges around the vertex face fan */
1250                 {
1251                         BMIter liter;
1252                         BMLoop *l;
1253                         BM_ITER_ELEM (l, &liter, v_other, BM_LOOPS_OF_VERT) {
1254                                 if (l->f->len == 3) {
1255                                         BMEdge *e_outer;
1256                                         if (BM_vert_in_edge(l->prev->e, l->v))
1257                                                 e_outer = l->next->e;
1258                                         else
1259                                                 e_outer = l->prev->e;
1260
1261                                         BLI_assert(BM_vert_in_edge(e_outer, l->v) == false);
1262
1263                                         bm_decim_build_edge_cost_single(e_outer, vquadrics, vweights, vweight_factor, eheap, eheap_table);
1264                                 }
1265                         }
1266                 }
1267                 /* end optional update */
1268                 return true;
1269 #endif
1270         }
1271         else {
1272                 /* add back with a high cost */
1273                 bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);
1274                 return false;
1275         }
1276 }
1277
1278
1279 /* Main Decimate Function
1280  * ********************** */
1281
1282 /**
1283  * \brief BM_mesh_decimate
1284  * \param bm The mesh
1285  * \param factor face count multiplier [0 - 1]
1286  * \param vweights Optional array of vertex  aligned weights [0 - 1],
1287  *        a vertex group is the usual source for this.
1288  * \param symmetry_axis: Axis of symmetry, -1 to disable mirror decimate.
1289  * \param symmetry_eps: Threshold when matching mirror verts.
1290  */
1291 void BM_mesh_decimate_collapse(
1292         BMesh *bm,
1293         const float factor,
1294         float *vweights, float vweight_factor,
1295         const bool do_triangulate,
1296         const int symmetry_axis, const float symmetry_eps)
1297 {
1298         Heap *eheap;             /* edge heap */
1299         HeapNode **eheap_table;  /* edge index aligned table pointing to the eheap */
1300         Quadric *vquadrics;      /* vert index aligned quadrics */
1301         int tot_edge_orig;
1302         int face_tot_target;
1303
1304         CD_UseFlag customdata_flag = 0;
1305
1306 #ifdef USE_SYMMETRY
1307         bool use_symmetry = (symmetry_axis != -1);
1308         int *edge_symmetry_map;
1309 #endif
1310
1311 #ifdef USE_TRIANGULATE
1312         int edges_tri_tot = 0;
1313         /* temp convert quads to triangles */
1314         bool use_triangulate = bm_decim_triangulate_begin(bm, &edges_tri_tot);
1315 #else
1316         UNUSED_VARS(do_triangulate);
1317 #endif
1318
1319
1320         /* alloc vars */
1321         vquadrics = MEM_callocN(sizeof(Quadric) * bm->totvert, __func__);
1322         /* since some edges may be degenerate, we might be over allocing a little here */
1323         eheap = BLI_heap_new_ex(bm->totedge);
1324         eheap_table = MEM_mallocN(sizeof(HeapNode *) * bm->totedge, __func__);
1325         tot_edge_orig = bm->totedge;
1326
1327
1328         /* build initial edge collapse cost data */
1329         bm_decim_build_quadrics(bm, vquadrics);
1330
1331         bm_decim_build_edge_cost(bm, vquadrics, vweights, vweight_factor, eheap, eheap_table);
1332
1333         face_tot_target = bm->totface * factor;
1334         bm->elem_index_dirty |= BM_ALL;
1335
1336 #ifdef USE_SYMMETRY
1337         edge_symmetry_map = (use_symmetry) ? bm_edge_symmetry_map(bm, symmetry_axis, symmetry_eps) : NULL;
1338 #else
1339         UNUSED_VARS(symmetry_axis, symmetry_eps);
1340 #endif
1341
1342 #ifdef USE_CUSTOMDATA
1343         /* initialize customdata flag, we only need math for loops */
1344         if (CustomData_has_interp(&bm->vdata))  customdata_flag |= CD_DO_VERT;
1345         if (CustomData_has_interp(&bm->edata))  customdata_flag |= CD_DO_EDGE;
1346         if (CustomData_has_math(&bm->ldata))    customdata_flag |= CD_DO_LOOP;
1347 #endif
1348
1349         /* iterative edge collapse and maintain the eheap */
1350 #ifdef USE_SYMMETRY
1351         if (use_symmetry == false)
1352 #endif
1353         {
1354                 /* simple non-mirror case */
1355                 while ((bm->totface > face_tot_target) &&
1356                        (BLI_heap_is_empty(eheap) == false) &&
1357                        (BLI_heap_node_value(BLI_heap_top(eheap)) != COST_INVALID))
1358                 {
1359                         // const float value = BLI_heap_node_value(BLI_heap_top(eheap));
1360                         BMEdge *e = BLI_heap_pop_min(eheap);
1361                         float optimize_co[3];
1362                         BLI_assert(BM_elem_index_get(e) < tot_edge_orig);  /* handy to detect corruptions elsewhere */
1363
1364                         /* under normal conditions wont be accessed again,
1365                          * but NULL just incase so we don't use freed node */
1366                         eheap_table[BM_elem_index_get(e)] = NULL;
1367
1368                         bm_decim_edge_collapse(
1369                                 bm, e, vquadrics, vweights, vweight_factor, eheap, eheap_table,
1370 #ifdef USE_SYMMETRY
1371                                 edge_symmetry_map,
1372 #endif
1373                                 customdata_flag,
1374                                 optimize_co, true
1375                                 );
1376                 }
1377         }
1378 #ifdef USE_SYMMETRY
1379         else {
1380                 while ((bm->totface > face_tot_target) &&
1381                        (BLI_heap_is_empty(eheap) == false) &&
1382                        (BLI_heap_node_value(BLI_heap_top(eheap)) != COST_INVALID))
1383                 {
1384                         /**
1385                          * \note
1386                          * - `eheap_table[e_index_mirr]` is only removed from the heap at the last moment
1387                          *   since its possible (in theory) for collapsing `e` to remove `e_mirr`.
1388                          * - edges sharing a vertex are ignored, so the pivot vertex isnt moved to one side.
1389                          */
1390
1391                         BMEdge *e = BLI_heap_pop_min(eheap);
1392                         const int e_index = BM_elem_index_get(e);
1393                         const int e_index_mirr = edge_symmetry_map[e_index];
1394                         BMEdge *e_mirr = NULL;
1395                         float optimize_co[3];
1396                         char e_invalidate = 0;
1397
1398                         BLI_assert(e_index < tot_edge_orig);
1399
1400                         eheap_table[e_index] = NULL;
1401
1402                         if (e_index_mirr != -1) {
1403                                 if (e_index_mirr == e_index) {
1404                                         /* pass */
1405                                 }
1406                                 else if (eheap_table[e_index_mirr]) {
1407                                         e_mirr = BLI_heap_node_ptr(eheap_table[e_index_mirr]);
1408                                         /* for now ignore edges with a shared vertex */
1409                                         if (BM_edge_share_vert_check(e, e_mirr)) {
1410                                                 /* ignore permanently!
1411                                                  * Otherwise we would keep re-evaluating and attempting to collapse. */
1412                                                 // e_invalidate |= (1 | 2);
1413                                                 goto invalidate;
1414                                         }
1415                                 }
1416                                 else {
1417                                         /* mirror edge can't be operated on (happens with asymmetrical meshes) */
1418                                         e_invalidate |= 1;
1419                                         goto invalidate;
1420                                 }
1421                         }
1422
1423                         /* when false, use without degenerate checks */
1424                         {
1425                                 /* run both before checking (since they invalidate surrounding geometry) */
1426                                 bool ok_a, ok_b;
1427
1428                                 ok_a = !bm_edge_collapse_is_degenerate_topology(e);
1429                                 ok_b = e_mirr ? !bm_edge_collapse_is_degenerate_topology(e_mirr) : true;
1430
1431                                 /* disallow collapsing which results in degenerate cases */
1432
1433                                 if (UNLIKELY(!ok_a || !ok_b)) {
1434                                         e_invalidate |= (1 | (e_mirr ? 2 : 0));
1435                                         goto invalidate;
1436                                 }
1437
1438                                 bm_decim_calc_target_co_fl(e, optimize_co, vquadrics);
1439
1440                                 if (e_index_mirr == e_index) {
1441                                         optimize_co[symmetry_axis] = 0.0f;
1442                                 }
1443
1444                                 /* check if this would result in an overlapping face */
1445                                 if (UNLIKELY(bm_edge_collapse_is_degenerate_flip(e, optimize_co))) {
1446                                         e_invalidate |= (1 | (e_mirr ? 2 : 0));
1447                                         goto invalidate;
1448                                 }
1449                         }
1450
1451                         if (bm_decim_edge_collapse(
1452                                 bm, e, vquadrics, vweights, vweight_factor, eheap, eheap_table,
1453                                 edge_symmetry_map,
1454                                 customdata_flag,
1455                                 optimize_co, false))
1456                         {
1457                                 if (e_mirr && (eheap_table[e_index_mirr])) {
1458                                         BLI_assert(e_index_mirr != e_index);
1459                                         BLI_heap_remove(eheap, eheap_table[e_index_mirr]);
1460                                         eheap_table[e_index_mirr] = NULL;
1461                                         optimize_co[symmetry_axis] *= -1.0f;
1462                                         bm_decim_edge_collapse(
1463                                                 bm, e_mirr, vquadrics, vweights, vweight_factor, eheap, eheap_table,
1464                                                 edge_symmetry_map,
1465                                                 customdata_flag,
1466                                                 optimize_co, false);
1467                                 }
1468                         }
1469                         else {
1470                                 if (e_mirr && (eheap_table[e_index_mirr])) {
1471                                         e_invalidate |= 2;
1472                                         goto invalidate;
1473                                 }
1474                         }
1475
1476                         BLI_assert(e_invalidate == 0);
1477                         continue;
1478
1479 invalidate:
1480                         if (e_invalidate & 1) {
1481                                 bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);
1482                         }
1483
1484                         if (e_invalidate & 2) {
1485                                 BLI_assert(eheap_table[e_index_mirr] != NULL);
1486                                 BLI_heap_remove(eheap, eheap_table[e_index_mirr]);
1487                                 eheap_table[e_index_mirr] = NULL;
1488                                 bm_decim_invalid_edge_cost_single(e_mirr, eheap, eheap_table);
1489                         }
1490                 }
1491
1492                 MEM_freeN((void *)edge_symmetry_map);
1493         }
1494 #endif  /* USE_SYMMETRY */
1495
1496 #ifdef USE_TRIANGULATE
1497         if (do_triangulate == false) {
1498                 /* its possible we only had triangles, skip this step in that case */
1499                 if (LIKELY(use_triangulate)) {
1500                         /* temp convert quads to triangles */
1501                         bm_decim_triangulate_end(bm, edges_tri_tot);
1502                 }
1503         }
1504 #endif
1505
1506         /* free vars */
1507         MEM_freeN(vquadrics);
1508         MEM_freeN(eheap_table);
1509         BLI_heap_free(eheap, NULL);
1510
1511         /* testing only */
1512         // BM_mesh_validate(bm);
1513
1514         (void)tot_edge_orig;  /* quiet release build warning */
1515 }