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