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