BMesh decimate, improve behavior with weights
[blender-staging.git] / source / blender / bmesh / tools / bmesh_decimate_collapse.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): Campbell Barton
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/tools/bmesh_decimate_collapse.c
24  *  \ingroup bmesh
25  *
26  * BMesh decimator that uses an edge collapse method.
27  */
28
29 #include <stddef.h>
30
31 #include "MEM_guardedalloc.h"
32
33
34 #include "BLI_math.h"
35 #include "BLI_quadric.h"
36 #include "BLI_heap.h"
37
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 /* defines for testing */
46 #define USE_CUSTOMDATA
47 #define USE_TRIANGULATE
48 #define USE_VERT_NORMAL_INTERP  /* has the advantage that flipped faces don't mess up vertex normals */
49
50 /* if the cost from #BLI_quadric_evaluate is 'noise', fallback to topology */
51 #define USE_TOPOLOGY_FALLBACK
52 #ifdef  USE_TOPOLOGY_FALLBACK
53 #  define   TOPOLOGY_FALLBACK_EPS  FLT_EPSILON
54 #endif
55
56 /* these checks are for rare cases that we can't avoid since they are valid meshes still */
57 #define USE_SAFETY_CHECKS
58
59 #define BOUNDARY_PRESERVE_WEIGHT 100.0f
60 #define OPTIMIZE_EPS 0.01f  /* FLT_EPSILON is too small, see [#33106] */
61 #define COST_INVALID FLT_MAX
62
63 typedef enum CD_UseFlag {
64         CD_DO_VERT = (1 << 0),
65         CD_DO_EDGE = (1 << 1),
66         CD_DO_LOOP = (1 << 2)
67 } CD_UseFlag;
68
69
70 /* BMesh Helper Functions
71  * ********************** */
72
73 /**
74  * \param vquadrics must be calloc'd
75  */
76 static void bm_decim_build_quadrics(BMesh *bm, Quadric *vquadrics)
77 {
78         BMIter iter;
79         BMFace *f;
80         BMEdge *e;
81
82         BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
83                 BMLoop *l_first;
84                 BMLoop *l_iter;
85
86                 float center[3];
87                 double plane_db[4];
88                 Quadric q;
89
90                 BM_face_calc_center_mean(f, center);
91                 copy_v3db_v3fl(plane_db, f->no);
92                 plane_db[3] = -dot_v3db_v3fl(plane_db, center);
93
94                 BLI_quadric_from_plane(&q, plane_db);
95
96                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
97                 do {
98                         BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(l_iter->v)], &q);
99                 } while ((l_iter = l_iter->next) != l_first);
100         }
101
102         /* boundary edges */
103         BM_ITER_MESH (e, &iter, bm, BM_EDGES_OF_MESH) {
104                 if (UNLIKELY(BM_edge_is_boundary(e))) {
105                         float edge_vector[3];
106                         float edge_plane[3];
107                         double edge_plane_db[4];
108                         sub_v3_v3v3(edge_vector, e->v2->co, e->v1->co);
109                         f = e->l->f;
110
111                         cross_v3_v3v3(edge_plane, edge_vector, f->no);
112                         copy_v3db_v3fl(edge_plane_db, edge_plane);
113
114                         if (normalize_v3_d(edge_plane_db) > FLT_EPSILON) {
115                                 Quadric q;
116                                 float center[3];
117
118                                 mid_v3_v3v3(center, e->v1->co, e->v2->co);
119
120                                 edge_plane_db[3] = -dot_v3db_v3fl(edge_plane_db, center);
121                                 BLI_quadric_from_plane(&q, edge_plane_db);
122                                 BLI_quadric_mul(&q, BOUNDARY_PRESERVE_WEIGHT);
123
124                                 BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(e->v1)], &q);
125                                 BLI_quadric_add_qu_qu(&vquadrics[BM_elem_index_get(e->v2)], &q);
126                         }
127                 }
128         }
129 }
130
131
132 static void bm_decim_calc_target_co(
133         BMEdge *e, float optimize_co[3],
134         const Quadric *vquadrics)
135 {
136         /* compute an edge contraction target for edge 'e'
137          * this is computed by summing it's vertices quadrics and
138          * optimizing the result. */
139         Quadric q;
140
141         BLI_quadric_add_qu_ququ(
142                 &q,
143                 &vquadrics[BM_elem_index_get(e->v1)],
144                 &vquadrics[BM_elem_index_get(e->v2)]);
145
146         if (BLI_quadric_optimize(&q, optimize_co, OPTIMIZE_EPS)) {
147                 return;  /* all is good */
148         }
149         else {
150                 mid_v3_v3v3(optimize_co, e->v1->co, e->v2->co);
151         }
152 }
153
154 static bool bm_edge_collapse_is_degenerate_flip(BMEdge *e, const float optimize_co[3])
155 {
156         BMIter liter;
157         BMLoop *l;
158         unsigned int i;
159
160         for (i = 0; i < 2; i++) {
161                 /* loop over both verts */
162                 BMVert *v = *((&e->v1) + i);
163
164                 BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
165                         if (l->e != e && l->prev->e != e) {
166                                 const float *co_prev = l->prev->v->co;
167                                 const float *co_next = l->next->v->co;
168                                 float cross_exist[3];
169                                 float cross_optim[3];
170
171 #if 1
172                                 float vec_other[3];  /* line between the two outer verts, re-use for both cross products */
173                                 float vec_exist[3];  /* before collapse */
174                                 float vec_optim[3];  /* after collapse */
175
176                                 sub_v3_v3v3(vec_other, co_prev, co_next);
177                                 sub_v3_v3v3(vec_exist, co_prev, v->co);
178                                 sub_v3_v3v3(vec_optim, co_prev, optimize_co);
179
180                                 cross_v3_v3v3(cross_exist, vec_other, vec_exist);
181                                 cross_v3_v3v3(cross_optim, vec_other, vec_optim);
182
183                                 /* avoid normalize */
184                                 if (dot_v3v3(cross_exist, cross_optim) <=
185                                     (len_squared_v3(cross_exist) + len_squared_v3(cross_optim)) * 0.01f)
186                                 {
187                                         return true;
188                                 }
189 #else
190                                 normal_tri_v3(cross_exist, v->co,       co_prev, co_next);
191                                 normal_tri_v3(cross_optim, optimize_co, co_prev, co_next);
192
193                                 /* use a small value rather then zero so we don't flip a face in multiple steps
194                                  * (first making it zero area, then flipping again) */
195                                 if (dot_v3v3(cross_exist, cross_optim) <= FLT_EPSILON) {
196                                         //printf("no flip\n");
197                                         return true;
198                                 }
199 #endif
200
201                         }
202                 }
203         }
204
205         return false;
206 }
207
208 #ifdef USE_TOPOLOGY_FALLBACK
209 /**
210  * when the cost is so small that its not useful (flat surfaces),
211  * fallback to using a 'topology' cost.
212  *
213  * This avoids cases where a flat (or near flat) areas get very un-even geometry.
214  */
215 static float bm_decim_build_edge_cost_single_squared__topology(BMEdge *e)
216 {
217         return fabsf(dot_v3v3(e->v1->no, e->v2->no)) / min_ff(-len_squared_v3v3(e->v1->co, e->v2->co), -FLT_EPSILON);
218 }
219 static float bm_decim_build_edge_cost_single__topology(BMEdge *e)
220 {
221         return fabsf(dot_v3v3(e->v1->no, e->v2->no)) / min_ff(-len_v3v3(e->v1->co, e->v2->co), -FLT_EPSILON);
222 }
223
224 #endif  /* USE_TOPOLOGY_FALLBACK */
225
226 static void bm_decim_build_edge_cost_single(
227         BMEdge *e,
228         const Quadric *vquadrics, const float *vweights,
229         Heap *eheap, HeapNode **eheap_table)
230 {
231         const Quadric *q1, *q2;
232         float optimize_co[3];
233         float cost;
234
235         if (eheap_table[BM_elem_index_get(e)]) {
236                 BLI_heap_remove(eheap, eheap_table[BM_elem_index_get(e)]);
237         }
238
239         /* check we can collapse, some edges we better not touch */
240         if (BM_edge_is_boundary(e)) {
241                 if (e->l->f->len == 3) {
242                         /* pass */
243                 }
244                 else {
245                         /* only collapse tri's */
246                         goto clear;
247                 }
248         }
249         else if (BM_edge_is_manifold(e)) {
250                 if ((e->l->f->len == 3) && (e->l->radial_next->f->len == 3)) {
251                         /* pass */
252                 }
253                 else {
254                         /* only collapse tri's */
255                         goto clear;
256                 }
257         }
258         else {
259                 goto clear;
260         }
261         /* end sanity check */
262
263
264         bm_decim_calc_target_co(e, optimize_co, vquadrics);
265
266         q1 = &vquadrics[BM_elem_index_get(e->v1)];
267         q2 = &vquadrics[BM_elem_index_get(e->v2)];
268
269         cost = (BLI_quadric_evaluate(q1, optimize_co) +
270                 BLI_quadric_evaluate(q2, optimize_co));
271
272
273         // print("COST %.12f\n");
274
275         /* note, 'cost' shouldn't be negative but happens sometimes with small values.
276          * this can cause faces that make up a flat surface to over-collapse, see [#37121] */
277         cost = fabsf(cost);
278
279         /* edge weights can be 2 max, ((2 * 2) ^ 2) == 16x max,
280          * a high weight means being inserted into the heap ahead of 16x better edges.
281          * but thats fine since its the purpose of weighting. */
282 #define VWEIGHT_SCALE 2
283
284 #ifdef USE_TOPOLOGY_FALLBACK
285         if (UNLIKELY(cost < TOPOLOGY_FALLBACK_EPS)) {
286                 /* subtract existing cost to further differentiate edges from one another
287                  *
288                  * keep topology cost below 0.0 so their values don't interfere with quadric cost,
289                  * (and they get handled first).
290                  * */
291                 if (vweights == NULL) {
292                         cost = bm_decim_build_edge_cost_single_squared__topology(e) - cost;
293                 }
294                 else {
295                         /* with weights we need to use the real length so we can scale them properly */
296                         float e_weight = (vweights[BM_elem_index_get(e->v1)] +
297                                           vweights[BM_elem_index_get(e->v2)]);
298                         cost = bm_decim_build_edge_cost_single__topology(e) - cost;
299                         /* note, this is rather arbitrary max weight is 2 here,
300                          * allow for skipping edges 4x the length, based on weights */
301                         if (e_weight) {
302                                 e_weight *= VWEIGHT_SCALE;
303                                 cost *= 1.0f + SQUARE(e_weight);
304                         }
305
306                         BLI_assert(cost <= 0.0f);
307                 }
308         }
309         else
310 #endif
311         if (vweights) {
312                 float e_weight = 2.0f - (vweights[BM_elem_index_get(e->v1)] +
313                                          vweights[BM_elem_index_get(e->v2)]);
314                 if (e_weight) {
315                         e_weight *= VWEIGHT_SCALE;
316                         cost += (BM_edge_calc_length(e) * SQUARE(e_weight));
317                 }
318         }
319 #undef VWEIGHT_SCALE
320
321         eheap_table[BM_elem_index_get(e)] = BLI_heap_insert(eheap, cost, e);
322         return;
323
324 clear:
325         eheap_table[BM_elem_index_get(e)] = NULL;
326 }
327
328
329 /* use this for degenerate cases - add back to the heap with an invalid cost,
330  * this way it may be calculated again if surrounding geometry changes */
331 static void bm_decim_invalid_edge_cost_single(
332         BMEdge *e,
333         Heap *eheap, HeapNode **eheap_table)
334 {
335         BLI_assert(eheap_table[BM_elem_index_get(e)] == NULL);
336         eheap_table[BM_elem_index_get(e)] = BLI_heap_insert(eheap, COST_INVALID, e);
337 }
338
339 static void bm_decim_build_edge_cost(
340         BMesh *bm,
341         const Quadric *vquadrics, const float *vweights,
342         Heap *eheap, HeapNode **eheap_table)
343 {
344         BMIter iter;
345         BMEdge *e;
346         unsigned int i;
347
348         BM_ITER_MESH_INDEX (e, &iter, bm, BM_EDGES_OF_MESH, i) {
349                 eheap_table[i] = NULL;  /* keep sanity check happy */
350                 bm_decim_build_edge_cost_single(e, vquadrics, vweights, eheap, eheap_table);
351         }
352 }
353
354 #ifdef USE_TRIANGULATE
355 /* Temp Triangulation
356  * ****************** */
357
358 /**
359  * To keep things simple we can only collapse edges on triangulated data
360  * (limitation with edge collapse and error calculation functions).
361  *
362  * But to avoid annoying users by only giving triangle results, we can
363  * triangulate, keeping a reference between the faces, then join after
364  * if the edges don't collapse, this will also allow more choices when
365  * collapsing edges so even has some advantage over decimating quads
366  * directly.
367  *
368  * \return true if any faces were triangulated.
369  */
370
371 static bool bm_decim_triangulate_begin(BMesh *bm)
372 {
373         BMIter iter;
374         BMFace *f;
375         // bool has_quad;  // could optimize this a little
376         bool has_cut = false;
377
378         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
379
380         /* first clear loop index values */
381         BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
382                 BMLoop *l_iter;
383                 BMLoop *l_first;
384
385                 l_iter = l_first = BM_FACE_FIRST_LOOP(f);
386                 do {
387                         BM_elem_index_set(l_iter, -1);  /* set_dirty */
388                 } while ((l_iter = l_iter->next) != l_first);
389
390                 // has_quad |= (f->len == 4)
391         }
392
393         bm->elem_index_dirty |= BM_LOOP;
394
395         /* adding new faces as we loop over faces
396          * is normally best avoided, however in this case its not so bad because any face touched twice
397          * will already be triangulated*/
398         BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
399                 if (f->len == 4) {
400                         BMLoop *f_l[4];
401                         BMLoop *l_a, *l_b;
402
403                         {
404                                 BMLoop *l_iter = BM_FACE_FIRST_LOOP(f);
405
406                                 f_l[0] = l_iter; l_iter = l_iter->next;
407                                 f_l[1] = l_iter; l_iter = l_iter->next;
408                                 f_l[2] = l_iter; l_iter = l_iter->next;
409                                 f_l[3] = l_iter;
410                         }
411
412                         if (len_squared_v3v3(f_l[0]->v->co, f_l[2]->v->co) <
413                             len_squared_v3v3(f_l[1]->v->co, f_l[3]->v->co))
414                         {
415                                 l_a = f_l[0];
416                                 l_b = f_l[2];
417                         }
418                         else {
419                                 l_a = f_l[1];
420                                 l_b = f_l[3];
421                         }
422
423 #ifdef USE_SAFETY_CHECKS
424                         if (BM_edge_exists(l_a->v, l_b->v) == NULL)
425 #endif
426                         {
427                                 BMFace *f_new;
428                                 BMLoop *l_new;
429
430                                 /* warning, NO_DOUBLE option here isn't handled as nice as it could be
431                                  * - if there is a quad that has a free standing edge joining it along
432                                  * where we want to split the face, there isnt a good way we can handle this.
433                                  * currently that edge will get removed when joining the tris back into a quad. */
434                                 f_new = BM_face_split(bm, f, l_a, l_b, &l_new, NULL, false);
435
436                                 if (f_new) {
437                                         /* the value of this doesn't matter, only that the 2 loops match and have unique values */
438                                         const int f_index = BM_elem_index_get(f);
439
440                                         /* since we just split theres only ever 2 loops */
441                                         BLI_assert(BM_edge_is_manifold(l_new->e));
442
443                                         BM_elem_index_set(l_new, f_index);  /* set_dirty */
444                                         BM_elem_index_set(l_new->radial_next, f_index);  /* set_dirty */
445
446                                         BM_face_normal_update(f);
447                                         BM_face_normal_update(f_new);
448
449                                         has_cut = true;
450                                 }
451                         }
452                 }
453         }
454
455         BLI_assert((bm->elem_index_dirty & BM_VERT) == 0);
456
457         if (has_cut) {
458                 /* now triangulation is done we need to correct index values */
459                 BM_mesh_elem_index_ensure(bm, BM_EDGE | BM_FACE);
460         }
461
462         return has_cut;
463 }
464
465 static void bm_decim_triangulate_end(BMesh *bm)
466 {
467         /* decimation finished, now re-join */
468         BMIter iter;
469         BMEdge *e, *e_next;
470
471         /* boundary edges */
472         BM_ITER_MESH_MUTABLE (e, e_next, &iter, bm, BM_EDGES_OF_MESH) {
473                 BMLoop *l_a, *l_b;
474                 if (BM_edge_loop_pair(e, &l_a, &l_b)) {
475                         const int l_a_index = BM_elem_index_get(l_a);
476                         if (l_a_index != -1) {
477                                 const int l_b_index = BM_elem_index_get(l_b);
478                                 if (l_a_index == l_b_index) {
479                                         if (LIKELY(l_a->f->len == 3 && l_b->f->len == 3)) {
480                                                 if (l_a->v != l_b->v) {  /* if this is the case, faces have become flipped */
481                                                         /* check we are not making a degenerate quad */
482                                                         BMVert *vquad[4] = {
483                                                                 e->v1,
484                                                                 BM_vert_in_edge(e, l_a->next->v) ? l_a->prev->v : l_a->next->v,
485                                                                 e->v2,
486                                                                 BM_vert_in_edge(e, l_b->next->v) ? l_b->prev->v : l_b->next->v,
487                                                         };
488
489                                                         BLI_assert(ELEM(vquad[0], vquad[1], vquad[2], vquad[3]) == false);
490                                                         BLI_assert(ELEM(vquad[1], vquad[0], vquad[2], vquad[3]) == false);
491                                                         BLI_assert(ELEM(vquad[2], vquad[1], vquad[0], vquad[3]) == false);
492                                                         BLI_assert(ELEM(vquad[3], vquad[1], vquad[2], vquad[0]) == false);
493
494                                                         if (is_quad_convex_v3(vquad[0]->co, vquad[1]->co, vquad[2]->co, vquad[3]->co)) {
495                                                                 /* highly unlikely to fail, but prevents possible double-ups */
496                                                                 BMFace *f[2] = {l_a->f, l_b->f};
497                                                                 BM_faces_join(bm, f, 2, true);
498                                                         }
499                                                 }
500                                         }
501                                 }
502                         }
503                 }
504         }
505 }
506
507 #endif  /* USE_TRIANGULATE */
508
509 /* Edge Collapse Functions
510  * *********************** */
511
512 #ifdef USE_CUSTOMDATA
513
514 /**
515  * \param l: defines the vert to collapse into.
516  */
517 static void bm_edge_collapse_loop_customdata(
518         BMesh *bm, BMLoop *l, BMVert *v_clear, BMVert *v_other,
519         const float customdata_fac)
520 {
521         /* disable seam check - the seam check would have to be done per layer, its not really that important */
522 //#define USE_SEAM
523         /* these don't need to be updated, since they will get removed when the edge collapses */
524         BMLoop *l_clear, *l_other;
525         const bool is_manifold = BM_edge_is_manifold(l->e);
526         int side;
527
528         /* first find the loop of 'v_other' thats attached to the face of 'l' */
529         if (l->v == v_clear) {
530                 l_clear = l;
531                 l_other = l->next;
532         }
533         else {
534                 l_clear = l->next;
535                 l_other = l;
536         }
537
538         BLI_assert(l_clear->v == v_clear);
539         BLI_assert(l_other->v == v_other);
540         (void)v_other;  /* quiet warnings for release */
541
542         /* now we have both corners of the face 'l->f' */
543         for (side = 0; side < 2; side++) {
544 #ifdef USE_SEAM
545                 bool is_seam = false;
546 #endif
547                 void *src[2];
548                 BMFace *f_exit = is_manifold ? l->radial_next->f : NULL;
549                 BMEdge *e_prev = l->e;
550                 BMLoop *l_first;
551                 BMLoop *l_iter;
552                 float w[2];
553
554                 if (side == 0) {
555                         l_iter = l_first = l_clear;
556                         src[0] = l_clear->head.data;
557                         src[1] = l_other->head.data;
558
559                         w[0] = customdata_fac;
560                         w[1] = 1.0f - customdata_fac;
561                 }
562                 else {
563                         l_iter = l_first = l_other;
564                         src[0] = l_other->head.data;
565                         src[1] = l_clear->head.data;
566
567                         w[0] = 1.0f - customdata_fac;
568                         w[1] = customdata_fac;
569                 }
570
571                 // print_v2("weights", w);
572
573                 /* WATCH IT! - should NOT reference (_clear or _other) vars for this while loop */
574
575                 /* walk around the fan using 'e_prev' */
576                 while (((l_iter = BM_vert_step_fan_loop(l_iter, &e_prev)) != l_first) && (l_iter != NULL)) {
577                         int i;
578                         /* quit once we hit the opposite face, if we have one */
579                         if (f_exit && UNLIKELY(f_exit == l_iter->f)) {
580                                 break;
581                         }
582
583 #ifdef USE_SEAM
584                         /* break out unless we find a match */
585                         is_seam = true;
586 #endif
587
588                         /* ok. we have a loop. now be smart with it! */
589                         for (i = 0; i < bm->ldata.totlayer; i++) {
590                                 if (CustomData_layer_has_math(&bm->ldata, i)) {
591                                         const int offset = bm->ldata.layers[i].offset;
592                                         const int type = bm->ldata.layers[i].type;
593                                         const void *cd_src[2] = {
594                                             POINTER_OFFSET(src[0], offset),
595                                             POINTER_OFFSET(src[1], offset),
596                                         };
597                                         void *cd_iter = POINTER_OFFSET(l_iter->head.data, offset);
598
599                                         /* detect seams */
600                                         if (CustomData_data_equals(type, cd_src[0], cd_iter)) {
601                                                 CustomData_bmesh_interp_n(
602                                                         &bm->ldata, cd_src, w, NULL, ARRAY_SIZE(cd_src),
603                                                         POINTER_OFFSET(l_iter->head.data, offset), i);
604 #ifdef USE_SEAM
605                                                 is_seam = false;
606 #endif
607                                         }
608                                 }
609                         }
610
611 #ifdef USE_SEAM
612                         if (is_seam) {
613                                 break;
614                         }
615 #endif
616                 }
617         }
618
619 //#undef USE_SEAM
620
621 }
622 #endif  /* USE_CUSTOMDATA */
623
624 /**
625  * Check if the collapse will result in a degenerate mesh,
626  * that is - duplicate edges or faces.
627  *
628  * This situation could be checked for when calculating collapse cost
629  * however its quite slow and a degenerate collapse could eventuate
630  * after the cost is calculated, so instead, check just before collapsing.
631  */
632
633 static void bm_edge_tag_enable(BMEdge *e)
634 {
635         BM_elem_flag_enable(e->v1, BM_ELEM_TAG);
636         BM_elem_flag_enable(e->v2, BM_ELEM_TAG);
637         if (e->l) {
638                 BM_elem_flag_enable(e->l->f, BM_ELEM_TAG);
639                 if (e->l != e->l->radial_next) {
640                         BM_elem_flag_enable(e->l->radial_next->f, BM_ELEM_TAG);
641                 }
642         }
643 }
644
645 static void bm_edge_tag_disable(BMEdge *e)
646 {
647         BM_elem_flag_disable(e->v1, BM_ELEM_TAG);
648         BM_elem_flag_disable(e->v2, BM_ELEM_TAG);
649         if (e->l) {
650                 BM_elem_flag_disable(e->l->f, BM_ELEM_TAG);
651                 if (e->l != e->l->radial_next) {
652                         BM_elem_flag_disable(e->l->radial_next->f, BM_ELEM_TAG);
653                 }
654         }
655 }
656
657 static bool bm_edge_tag_test(BMEdge *e)
658 {
659         /* is the edge or one of its faces tagged? */
660         return (BM_elem_flag_test(e->v1, BM_ELEM_TAG) ||
661                 BM_elem_flag_test(e->v2, BM_ELEM_TAG) ||
662                 (e->l && (BM_elem_flag_test(e->l->f, BM_ELEM_TAG) ||
663                           (e->l != e->l->radial_next &&
664                           BM_elem_flag_test(e->l->radial_next->f, BM_ELEM_TAG))))
665                 );
666 }
667
668 /* takes the edges loop */
669 BLI_INLINE int bm_edge_is_manifold_or_boundary(BMLoop *l)
670 {
671 #if 0
672         /* less optimized version of check below */
673         return (BM_edge_is_manifold(l->e) || BM_edge_is_boundary(l->e);
674 #else
675         /* if the edge is a boundary it points to its self, else this must be a manifold */
676         return LIKELY(l) && LIKELY(l->radial_next->radial_next == l);
677 #endif
678 }
679
680 static bool bm_edge_collapse_is_degenerate_topology(BMEdge *e_first)
681 {
682         /* simply check that there is no overlap between faces and edges of each vert,
683          * (excluding the 2 faces attached to 'e' and 'e' its self) */
684
685         BMEdge *e_iter;
686
687         /* clear flags on both disks */
688         e_iter = e_first;
689         do {
690                 if (!bm_edge_is_manifold_or_boundary(e_iter->l)) {
691                         return true;
692                 }
693                 bm_edge_tag_disable(e_iter);
694         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v1)) != e_first);
695
696         e_iter = e_first;
697         do {
698                 if (!bm_edge_is_manifold_or_boundary(e_iter->l)) {
699                         return true;
700                 }
701                 bm_edge_tag_disable(e_iter);
702         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v2)) != e_first);
703
704         /* now enable one side... */
705         e_iter = e_first;
706         do {
707                 bm_edge_tag_enable(e_iter);
708         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v1)) != e_first);
709
710         /* ... except for the edge we will collapse, we know thats shared,
711          * disable this to avoid false positive. We could be smart and never enable these
712          * face/edge tags in the first place but easier to do this */
713         // bm_edge_tag_disable(e_first);
714         /* do inline... */
715         {
716 #if 0
717                 BMIter iter;
718                 BMIter liter;
719                 BMLoop *l;
720                 BMVert *v;
721                 BM_ITER_ELEM (l, &liter, e_first, BM_LOOPS_OF_EDGE) {
722                         BM_elem_flag_disable(l->f, BM_ELEM_TAG);
723                         BM_ITER_ELEM (v, &iter, l->f, BM_VERTS_OF_FACE) {
724                                 BM_elem_flag_disable(v, BM_ELEM_TAG);
725                         }
726                 }
727 #else
728                 /* we know each face is a triangle, no looping/iterators needed here */
729
730                 BMLoop *l_radial;
731                 BMLoop *l_face;
732
733                 l_radial = e_first->l;
734                 l_face = l_radial;
735                 BLI_assert(l_face->f->len == 3);
736                 BM_elem_flag_disable(l_face->f, BM_ELEM_TAG);
737                 BM_elem_flag_disable((l_face = l_radial)->v,     BM_ELEM_TAG);
738                 BM_elem_flag_disable((l_face = l_face->next)->v, BM_ELEM_TAG);
739                 BM_elem_flag_disable((         l_face->next)->v, BM_ELEM_TAG);
740                 l_face = l_radial->radial_next;
741                 if (l_radial != l_face) {
742                         BLI_assert(l_face->f->len == 3);
743                         BM_elem_flag_disable(l_face->f, BM_ELEM_TAG);
744                         BM_elem_flag_disable((l_face = l_radial->radial_next)->v, BM_ELEM_TAG);
745                         BM_elem_flag_disable((l_face = l_face->next)->v,          BM_ELEM_TAG);
746                         BM_elem_flag_disable((         l_face->next)->v,          BM_ELEM_TAG);
747                 }
748 #endif
749         }
750
751         /* and check for overlap */
752         e_iter = e_first;
753         do {
754                 if (bm_edge_tag_test(e_iter)) {
755                         return true;
756                 }
757         } while ((e_iter = bmesh_disk_edge_next(e_iter, e_first->v2)) != e_first);
758
759         return false;
760 }
761
762 /**
763  * special, highly limited edge collapse function
764  * intended for speed over flexibility.
765  * can only collapse edges connected to (1, 2) tris.
766  *
767  * Important - dont add vert/edge/face data on collapsing!
768  *
769  * \param r_e_clear_other: Let caller know what edges we remove besides \a e_clear
770  * \param customdata_flag: Merge factor, scales from 0 - 1 ('v_clear' -> 'v_other')
771  */
772 static bool bm_edge_collapse(
773         BMesh *bm, BMEdge *e_clear, BMVert *v_clear, int r_e_clear_other[2],
774 #ifdef USE_CUSTOMDATA
775         const CD_UseFlag customdata_flag,
776         const float customdata_fac
777 #else
778         const CD_UseFlag UNUSED(customdata_flag),
779         const float UNUSED(customdata_fac)
780 #endif
781         )
782 {
783         BMVert *v_other;
784
785         v_other = BM_edge_other_vert(e_clear, v_clear);
786         BLI_assert(v_other != NULL);
787
788         if (BM_edge_is_manifold(e_clear)) {
789                 BMLoop *l_a, *l_b;
790                 BMEdge *e_a_other[2], *e_b_other[2];
791                 bool ok;
792
793                 ok = BM_edge_loop_pair(e_clear, &l_a, &l_b);
794
795                 BLI_assert(ok == true);
796                 BLI_assert(l_a->f->len == 3);
797                 BLI_assert(l_b->f->len == 3);
798                 UNUSED_VARS_NDEBUG(ok);
799
800                 /* keep 'v_clear' 0th */
801                 if (BM_vert_in_edge(l_a->prev->e, v_clear)) {
802                         e_a_other[0] = l_a->prev->e;
803                         e_a_other[1] = l_a->next->e;
804                 }
805                 else {
806                         e_a_other[1] = l_a->prev->e;
807                         e_a_other[0] = l_a->next->e;
808                 }
809
810                 if (BM_vert_in_edge(l_b->prev->e, v_clear)) {
811                         e_b_other[0] = l_b->prev->e;
812                         e_b_other[1] = l_b->next->e;
813                 }
814                 else {
815                         e_b_other[1] = l_b->prev->e;
816                         e_b_other[0] = l_b->next->e;
817                 }
818
819                 /* we could assert this case, but better just bail out */
820 #if 0
821                 BLI_assert(e_a_other[0] != e_b_other[0]);
822                 BLI_assert(e_a_other[0] != e_b_other[1]);
823                 BLI_assert(e_b_other[0] != e_a_other[0]);
824                 BLI_assert(e_b_other[0] != e_a_other[1]);
825 #endif
826                 /* not totally common but we want to avoid */
827                 if (ELEM(e_a_other[0], e_b_other[0], e_b_other[1]) ||
828                     ELEM(e_a_other[1], e_b_other[0], e_b_other[1]))
829                 {
830                         return false;
831                 }
832
833                 BLI_assert(BM_edge_share_vert(e_a_other[0], e_b_other[0]));
834                 BLI_assert(BM_edge_share_vert(e_a_other[1], e_b_other[1]));
835
836                 r_e_clear_other[0] = BM_elem_index_get(e_a_other[0]);
837                 r_e_clear_other[1] = BM_elem_index_get(e_b_other[0]);
838
839 #ifdef USE_CUSTOMDATA
840                 /* before killing, do customdata */
841                 if (customdata_flag & CD_DO_VERT) {
842                         BM_data_interp_from_verts(bm, v_other, v_clear, v_other, customdata_fac);
843                 }
844                 if (customdata_flag & CD_DO_EDGE) {
845                         BM_data_interp_from_edges(bm, e_a_other[1], e_a_other[0], e_a_other[1], customdata_fac);
846                         BM_data_interp_from_edges(bm, e_b_other[1], e_b_other[0], e_b_other[1], customdata_fac);
847                 }
848                 if (customdata_flag & CD_DO_LOOP) {
849                         bm_edge_collapse_loop_customdata(bm, e_clear->l,              v_clear, v_other, customdata_fac);
850                         bm_edge_collapse_loop_customdata(bm, e_clear->l->radial_next, v_clear, v_other, customdata_fac);
851                 }
852 #endif
853
854                 BM_edge_kill(bm, e_clear);
855
856                 v_other->head.hflag |= v_clear->head.hflag;
857                 BM_vert_splice(bm, v_other, v_clear);
858
859                 e_a_other[1]->head.hflag |= e_a_other[0]->head.hflag;
860                 e_b_other[1]->head.hflag |= e_b_other[0]->head.hflag;
861                 BM_edge_splice(bm, e_a_other[1], e_a_other[0]);
862                 BM_edge_splice(bm, e_b_other[1], e_b_other[0]);
863
864                 // BM_mesh_validate(bm);
865
866                 return true;
867         }
868         else if (BM_edge_is_boundary(e_clear)) {
869                 /* same as above but only one triangle */
870                 BMLoop *l_a;
871                 BMEdge *e_a_other[2];
872
873                 l_a = e_clear->l;
874
875                 BLI_assert(l_a->f->len == 3);
876
877                 /* keep 'v_clear' 0th */
878                 if (BM_vert_in_edge(l_a->prev->e, v_clear)) {
879                         e_a_other[0] = l_a->prev->e;
880                         e_a_other[1] = l_a->next->e;
881                 }
882                 else {
883                         e_a_other[1] = l_a->prev->e;
884                         e_a_other[0] = l_a->next->e;
885                 }
886
887                 r_e_clear_other[0] = BM_elem_index_get(e_a_other[0]);
888                 r_e_clear_other[1] = -1;
889
890 #ifdef USE_CUSTOMDATA
891                 /* before killing, do customdata */
892                 if (customdata_flag & CD_DO_VERT) {
893                         BM_data_interp_from_verts(bm, v_other, v_clear, v_other, customdata_fac);
894                 }
895                 if (customdata_flag & CD_DO_EDGE) {
896                         BM_data_interp_from_edges(bm, e_a_other[1], e_a_other[0], e_a_other[1], customdata_fac);
897                 }
898                 if (customdata_flag & CD_DO_LOOP) {
899                         bm_edge_collapse_loop_customdata(bm, e_clear->l, v_clear, v_other, customdata_fac);
900                 }
901 #endif
902
903                 BM_edge_kill(bm, e_clear);
904
905                 v_other->head.hflag |= v_clear->head.hflag;
906                 BM_vert_splice(bm, v_other, v_clear);
907
908                 e_a_other[1]->head.hflag |= e_a_other[0]->head.hflag;
909                 BM_edge_splice(bm, e_a_other[1], e_a_other[0]);
910
911                 // BM_mesh_validate(bm);
912
913                 return true;
914         }
915         else {
916                 return false;
917         }
918 }
919
920
921 /* collapse e the edge, removing e->v2 */
922 static void bm_decim_edge_collapse(
923         BMesh *bm, BMEdge *e,
924         Quadric *vquadrics, float *vweights,
925         Heap *eheap, HeapNode **eheap_table,
926         const CD_UseFlag customdata_flag)
927 {
928         int e_clear_other[2];
929         BMVert *v_other = e->v1;
930         const int v_other_index = BM_elem_index_get(e->v1);
931         const int v_clear_index = BM_elem_index_get(e->v2);  /* the vert is removed so only store the index */
932         float optimize_co[3];
933         float customdata_fac;
934
935 #ifdef USE_VERT_NORMAL_INTERP
936         float v_clear_no[3];
937         copy_v3_v3(v_clear_no, e->v2->no);
938 #endif
939
940         /* disallow collapsing which results in degenerate cases */
941         if (UNLIKELY(bm_edge_collapse_is_degenerate_topology(e))) {
942                 bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);  /* add back with a high cost */
943                 return;
944         }
945
946         bm_decim_calc_target_co(e, optimize_co, vquadrics);
947
948         /* check if this would result in an overlapping face */
949         if (UNLIKELY(bm_edge_collapse_is_degenerate_flip(e, optimize_co))) {
950                 bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);  /* add back with a high cost */
951                 return;
952         }
953
954         /* use for customdata merging */
955         if (LIKELY(compare_v3v3(e->v1->co, e->v2->co, FLT_EPSILON) == false)) {
956                 customdata_fac = line_point_factor_v3(optimize_co, e->v1->co, e->v2->co);
957 #if 0
958                 /* simple test for stupid collapse */
959                 if (customdata_fac < 0.0 - FLT_EPSILON || customdata_fac > 1.0f + FLT_EPSILON) {
960                         return;
961                 }
962 #endif
963         }
964         else {
965                 /* avoid divide by zero */
966                 customdata_fac = 0.5f;
967         }
968
969         if (bm_edge_collapse(bm, e, e->v2, e_clear_other, customdata_flag, customdata_fac)) {
970                 /* update collapse info */
971                 int i;
972
973                 if (vweights) {
974                         float v_other_weight = interpf(vweights[v_other_index], vweights[v_clear_index], customdata_fac);
975                         CLAMP(v_other_weight, 0.0f, 1.0f);
976                         vweights[v_other_index] = v_other_weight;
977                 }
978
979                 e = NULL;  /* paranoid safety check */
980
981                 copy_v3_v3(v_other->co, optimize_co);
982
983                 /* remove eheap */
984                 for (i = 0; i < 2; i++) {
985                         /* highly unlikely 'eheap_table[ke_other[i]]' would be NULL, but do for sanity sake */
986                         if ((e_clear_other[i] != -1) && (eheap_table[e_clear_other[i]] != NULL)) {
987                                 BLI_heap_remove(eheap, eheap_table[e_clear_other[i]]);
988                                 eheap_table[e_clear_other[i]] = NULL;
989                         }
990                 }
991
992                 /* update vertex quadric, add kept vertex from killed vertex */
993                 BLI_quadric_add_qu_qu(&vquadrics[v_other_index], &vquadrics[v_clear_index]);
994
995                 /* update connected normals */
996
997                 /* in fact face normals are not used for progressive updates, no need to update them */
998                 // BM_vert_normal_update_all(v);
999 #ifdef USE_VERT_NORMAL_INTERP
1000                 interp_v3_v3v3(v_other->no, v_other->no, v_clear_no, customdata_fac);
1001                 normalize_v3(v_other->no);
1002 #else
1003                 BM_vert_normal_update(v_other);
1004 #endif
1005
1006
1007                 /* update error costs and the eheap */
1008                 if (LIKELY(v_other->e)) {
1009                         BMEdge *e_iter;
1010                         BMEdge *e_first;
1011                         e_iter = e_first = v_other->e;
1012                         do {
1013                                 BLI_assert(BM_edge_find_double(e_iter) == NULL);
1014                                 bm_decim_build_edge_cost_single(e_iter, vquadrics, vweights, eheap, eheap_table);
1015                         } while ((e_iter = bmesh_disk_edge_next(e_iter, v_other)) != e_first);
1016                 }
1017
1018                 /* this block used to be disabled,
1019                  * but enable now since surrounding faces may have been
1020                  * set to COST_INVALID because of a face overlap that no longer occurs */
1021 #if 1
1022                 /* optional, update edges around the vertex face fan */
1023                 {
1024                         BMIter liter;
1025                         BMLoop *l;
1026                         BM_ITER_ELEM (l, &liter, v_other, BM_LOOPS_OF_VERT) {
1027                                 if (l->f->len == 3) {
1028                                         BMEdge *e_outer;
1029                                         if (BM_vert_in_edge(l->prev->e, l->v))
1030                                                 e_outer = l->next->e;
1031                                         else
1032                                                 e_outer = l->prev->e;
1033
1034                                         BLI_assert(BM_vert_in_edge(e_outer, l->v) == false);
1035
1036                                         bm_decim_build_edge_cost_single(e_outer, vquadrics, vweights, eheap, eheap_table);
1037                                 }
1038                         }
1039                 }
1040                 /* end optional update */
1041 #endif
1042         }
1043         else {
1044                 /* add back with a high cost */
1045                 bm_decim_invalid_edge_cost_single(e, eheap, eheap_table);
1046         }
1047 }
1048
1049
1050 /* Main Decimate Function
1051  * ********************** */
1052
1053 /**
1054  * \brief BM_mesh_decimate
1055  * \param bm The mesh
1056  * \param factor face count multiplier [0 - 1]
1057  * \param vweights Optional array of vertex  aligned weights [0 - 1],
1058  *        a vertex group is the usual source for this.
1059  */
1060 void BM_mesh_decimate_collapse(BMesh *bm, const float factor, float *vweights, const bool do_triangulate)
1061 {
1062         Heap *eheap;             /* edge heap */
1063         HeapNode **eheap_table;  /* edge index aligned table pointing to the eheap */
1064         Quadric *vquadrics;      /* vert index aligned quadrics */
1065         int tot_edge_orig;
1066         int face_tot_target;
1067         bool use_triangulate;
1068
1069         CD_UseFlag customdata_flag = 0;
1070
1071 #ifdef USE_TRIANGULATE
1072         /* temp convert quads to triangles */
1073         use_triangulate = bm_decim_triangulate_begin(bm);
1074 #endif
1075
1076
1077         /* alloc vars */
1078         vquadrics = MEM_callocN(sizeof(Quadric) * bm->totvert, __func__);
1079         /* since some edges may be degenerate, we might be over allocing a little here */
1080         eheap = BLI_heap_new_ex(bm->totedge);
1081         eheap_table = MEM_mallocN(sizeof(HeapNode *) * bm->totedge, __func__);
1082         tot_edge_orig = bm->totedge;
1083
1084
1085         /* build initial edge collapse cost data */
1086         bm_decim_build_quadrics(bm, vquadrics);
1087
1088         bm_decim_build_edge_cost(bm, vquadrics, vweights, eheap, eheap_table);
1089
1090         face_tot_target = bm->totface * factor;
1091         bm->elem_index_dirty |= BM_ALL;
1092
1093
1094 #ifdef USE_CUSTOMDATA
1095         /* initialize customdata flag, we only need math for loops */
1096         if (CustomData_has_interp(&bm->vdata))  customdata_flag |= CD_DO_VERT;
1097         if (CustomData_has_interp(&bm->edata))  customdata_flag |= CD_DO_EDGE;
1098         if (CustomData_has_math(&bm->ldata))    customdata_flag |= CD_DO_LOOP;
1099 #endif
1100
1101         /* iterative edge collapse and maintain the eheap */
1102         while ((bm->totface > face_tot_target) &&
1103                (BLI_heap_is_empty(eheap) == false) &&
1104                (BLI_heap_node_value(BLI_heap_top(eheap)) != COST_INVALID))
1105         {
1106                 // const float value = BLI_heap_node_value(BLI_heap_top(eheap));
1107                 BMEdge *e = BLI_heap_popmin(eheap);
1108                 BLI_assert(BM_elem_index_get(e) < tot_edge_orig);  /* handy to detect corruptions elsewhere */
1109
1110                 // printf("COST %.10f\n", value);
1111
1112                 /* under normal conditions wont be accessed again,
1113                  * but NULL just incase so we don't use freed node */
1114                 eheap_table[BM_elem_index_get(e)] = NULL;
1115
1116                 bm_decim_edge_collapse(bm, e, vquadrics, vweights, eheap, eheap_table, customdata_flag);
1117         }
1118
1119
1120 #ifdef USE_TRIANGULATE
1121         if (do_triangulate == false) {
1122                 /* its possible we only had triangles, skip this step in that case */
1123                 if (LIKELY(use_triangulate)) {
1124                         /* temp convert quads to triangles */
1125                         bm_decim_triangulate_end(bm);
1126                 }
1127         }
1128 #endif
1129
1130         /* free vars */
1131         MEM_freeN(vquadrics);
1132         MEM_freeN(eheap_table);
1133         BLI_heap_free(eheap, NULL);
1134
1135         /* testing only */
1136         // BM_mesh_validate(bm);
1137
1138         (void)tot_edge_orig;  /* quiet release build warning */
1139 }