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