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