78e3e66b70aa39451259b4255ff5575c3d88c856
[blender.git] / source / blender / bmesh / tools / bmesh_beautify.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): Joseph Eagar.
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/tools/bmesh_beautify.c
24  *  \ingroup bmesh
25  *
26  * Beautify the mesh by rotating edges between triangles
27  * to more attractive positions until no more rotations can be made.
28  *
29  * In principle this is very simple however there is the possibility of
30  * going into an eternal loop where edges keep rotating.
31  * To avoid this - each edge stores a set of it previous
32  * states so as not to rotate back.
33  *
34  * TODO
35  * - Take face normals into account.
36  */
37
38 #include "BLI_math.h"
39 #include "BLI_heap.h"
40 #include "BLI_polyfill2d_beautify.h"
41
42 #include "MEM_guardedalloc.h"
43
44 #include "bmesh.h"
45 #include "bmesh_beautify.h"  /* own include */
46
47
48 // #define DEBUG_TIME
49
50 #ifdef DEBUG_TIME
51 #  include "PIL_time.h"
52 #  include "PIL_time_utildefines.h"
53 #endif
54
55 /* -------------------------------------------------------------------- */
56 /* GSet for edge rotation */
57
58 typedef struct EdRotState {
59         int v1, v2; /*  edge vert, small -> large */
60         int f1, f2; /*  face vert, small -> large */
61 } EdRotState;
62
63 #if 0
64 /* use BLI_ghashutil_inthash_v4 direct */
65 static uint erot_gsetutil_hash(const void *ptr)
66 {
67         const EdRotState *e_state = (const EdRotState *)ptr;
68         return BLI_ghashutil_inthash_v4(&e_state->v1);
69 }
70 #endif
71 #if 0
72 static int erot_gsetutil_cmp(const void *a, const void *b)
73 {
74         const EdRotState *e_state_a = (const EdRotState *)a;
75         const EdRotState *e_state_b = (const EdRotState *)b;
76         if      (e_state_a->v1 < e_state_b->v1) return -1;
77         else if (e_state_a->v1 > e_state_b->v1) return  1;
78         else if (e_state_a->v2 < e_state_b->v2) return -1;
79         else if (e_state_a->v2 > e_state_b->v2) return  1;
80         else if (e_state_a->f1 < e_state_b->f1) return -1;
81         else if (e_state_a->f1 > e_state_b->f1) return  1;
82         else if (e_state_a->f2 < e_state_b->f2) return -1;
83         else if (e_state_a->f2 > e_state_b->f2) return  1;
84         else                                    return  0;
85 }
86 #endif
87 static GSet *erot_gset_new(void)
88 {
89         return BLI_gset_new(BLI_ghashutil_inthash_v4_p, BLI_ghashutil_inthash_v4_cmp, __func__);
90 }
91
92 /* ensure v0 is smaller */
93 #define EDGE_ORD(v0, v1)   \
94         if (v0 > v1) {         \
95                 SWAP(int, v0, v1); \
96         } (void)0
97
98 static void erot_state_ex(const BMEdge *e, int v_index[2], int f_index[2])
99 {
100         BLI_assert(BM_edge_is_manifold(e));
101         BLI_assert(BM_vert_in_edge(e, e->l->prev->v)              == false);
102         BLI_assert(BM_vert_in_edge(e, e->l->radial_next->prev->v) == false);
103
104         /* verts of the edge */
105         v_index[0] = BM_elem_index_get(e->v1);
106         v_index[1] = BM_elem_index_get(e->v2);
107         EDGE_ORD(v_index[0], v_index[1]);
108
109         /* verts of each of the 2 faces attached to this edge
110          * (that are not apart of this edge) */
111         f_index[0] = BM_elem_index_get(e->l->prev->v);
112         f_index[1] = BM_elem_index_get(e->l->radial_next->prev->v);
113         EDGE_ORD(f_index[0], f_index[1]);
114 }
115
116 static void erot_state_current(const BMEdge *e, EdRotState *e_state)
117 {
118         erot_state_ex(e, &e_state->v1, &e_state->f1);
119 }
120
121 static void erot_state_alternate(const BMEdge *e, EdRotState *e_state)
122 {
123         erot_state_ex(e, &e_state->f1, &e_state->v1);
124 }
125
126 /* -------------------------------------------------------------------- */
127 /* Calculate the improvement of rotating the edge */
128
129 static float bm_edge_calc_rotate_beauty__area(
130         const float v1[3], const float v2[3], const float v3[3], const float v4[3])
131 {
132         /* not a loop (only to be able to break out) */
133         do {
134                 float v1_xy[2], v2_xy[2], v3_xy[2], v4_xy[2];
135
136                 /* first get the 2d values */
137                 {
138                         const float eps = 1e-5;
139                         float no_a[3], no_b[3];
140                         float no[3];
141                         float axis_mat[3][3];
142                         float no_scale;
143                         cross_tri_v3(no_a, v2, v3, v4);
144                         cross_tri_v3(no_b, v2, v4, v1);
145
146                         // printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
147                         BLI_assert((ELEM(v1, v2, v3, v4) == false) &&
148                                    (ELEM(v2, v1, v3, v4) == false) &&
149                                    (ELEM(v3, v1, v2, v4) == false) &&
150                                    (ELEM(v4, v1, v2, v3) == false));
151
152                         add_v3_v3v3(no, no_a, no_b);
153                         if (UNLIKELY((no_scale = normalize_v3(no)) <= FLT_EPSILON)) {
154                                 break;
155                         }
156
157                         axis_dominant_v3_to_m3(axis_mat, no);
158                         mul_v2_m3v3(v1_xy, axis_mat, v1);
159                         mul_v2_m3v3(v2_xy, axis_mat, v2);
160                         mul_v2_m3v3(v3_xy, axis_mat, v3);
161                         mul_v2_m3v3(v4_xy, axis_mat, v4);
162
163                         /**
164                          * Check if input faces are already flipped.
165                          * Logic for 'signum_i' addition is:
166                          *
167                          * Accept:
168                          * - (1, 1) or (-1, -1): same side (common case).
169                          * - (-1/1, 0): one degenerate, OK since we may rotate into a valid state.
170                          *
171                          * Ignore:
172                          * - (-1, 1): opposite winding, ignore.
173                          * - ( 0, 0): both degenerate, ignore.
174                          *
175                          * \note The cross product is divided by 'no_scale'
176                          * so the rotation calculation is scale independent.
177                          */
178                         if (!(signum_i_ex(cross_tri_v2(v2_xy, v3_xy, v4_xy) / no_scale, eps) +
179                               signum_i_ex(cross_tri_v2(v2_xy, v4_xy, v1_xy) / no_scale, eps)))
180                         {
181                                 break;
182                         }
183                 }
184
185                 /**
186                  * Important to lock degenerate here,
187                  * since the triangle pars will be projected into different 2D spaces.
188                  * Allowing to rotate out of a degenerate state can flip the faces (when performed iteratively).
189                  */
190                 return BLI_polyfill_beautify_quad_rotate_calc_ex(v1_xy, v2_xy, v3_xy, v4_xy, true);
191         } while (false);
192
193         return FLT_MAX;
194 }
195
196 static float bm_edge_calc_rotate_beauty__angle(
197         const float v1[3], const float v2[3], const float v3[3], const float v4[3])
198 {
199         /* not a loop (only to be able to break out) */
200         do {
201                 float no_a[3], no_b[3];
202                 float angle_24, angle_13;
203
204                 /* edge (2-4), current state */
205                 normal_tri_v3(no_a, v2, v3, v4);
206                 normal_tri_v3(no_b, v2, v4, v1);
207                 angle_24 = angle_normalized_v3v3(no_a, no_b);
208
209                 /* edge (1-3), new state */
210                 /* only check new state for degenerate outcome */
211                 if ((normal_tri_v3(no_a, v1, v2, v3) == 0.0f) ||
212                     (normal_tri_v3(no_b, v1, v3, v4) == 0.0f))
213                 {
214                         break;
215                 }
216                 angle_13 = angle_normalized_v3v3(no_a, no_b);
217
218                 return angle_13 - angle_24;
219         } while (false);
220
221         return FLT_MAX;
222 }
223
224 /**
225  * Assuming we have 2 triangles sharing an edge (2 - 4),
226  * check if the edge running from (1 - 3) gives better results.
227  *
228  * \return (negative number means the edge can be rotated, lager == better).
229  */
230 float BM_verts_calc_rotate_beauty(
231         const BMVert *v1, const BMVert *v2, const BMVert *v3, const BMVert *v4,
232         const short flag, const short method)
233 {
234         /* not a loop (only to be able to break out) */
235         do {
236                 if (flag & VERT_RESTRICT_TAG) {
237                         const BMVert *v_a = v1, *v_b = v3;
238                         if (BM_elem_flag_test(v_a, BM_ELEM_TAG) ==  BM_elem_flag_test(v_b, BM_ELEM_TAG)) {
239                                 break;
240                         }
241                 }
242
243                 if (UNLIKELY(v1 == v3)) {
244                         // printf("This should never happen, but does sometimes!\n");
245                         break;
246                 }
247
248                 switch (method) {
249                         case 0:
250                                 return bm_edge_calc_rotate_beauty__area(v1->co, v2->co, v3->co, v4->co);
251                         default:
252                                 return bm_edge_calc_rotate_beauty__angle(v1->co, v2->co, v3->co, v4->co);
253                 }
254         } while (false);
255
256         return FLT_MAX;
257 }
258
259 static float bm_edge_calc_rotate_beauty(const BMEdge *e, const short flag, const short method)
260 {
261         const BMVert *v1, *v2, *v3, *v4;
262         v1 = e->l->prev->v;               /* first vert co */
263         v2 = e->l->v;                     /* e->v1 or e->v2*/
264         v3 = e->l->radial_next->prev->v;  /* second vert co */
265         v4 = e->l->next->v;               /* e->v1 or e->v2*/
266
267         return BM_verts_calc_rotate_beauty(v1, v2, v3, v4, flag, method);
268 }
269
270 /* -------------------------------------------------------------------- */
271 /* Update the edge cost of rotation in the heap */
272
273 BLI_INLINE bool edge_in_array(const BMEdge *e, const BMEdge **edge_array, const int edge_array_len)
274 {
275         const int index = BM_elem_index_get(e);
276         return ((index >= 0) &&
277                 (index < edge_array_len) &&
278                 (e == edge_array[index]));
279 }
280
281 /* recalc an edge in the heap (surrounding geometry has changed) */
282 static void bm_edge_update_beauty_cost_single(
283         BMEdge *e, Heap *eheap, HeapNode **eheap_table, GSet **edge_state_arr,
284         /* only for testing the edge is in the array */
285         const BMEdge **edge_array, const int edge_array_len,
286
287         const short flag, const short method)
288 {
289         if (edge_in_array(e, edge_array, edge_array_len)) {
290                 const int i = BM_elem_index_get(e);
291                 GSet *e_state_set = edge_state_arr[i];
292
293                 if (eheap_table[i]) {
294                         BLI_heap_remove(eheap, eheap_table[i]);
295                         eheap_table[i] = NULL;
296                 }
297
298                 /* check if we can add it back */
299                 BLI_assert(BM_edge_is_manifold(e) == true);
300
301                 /* check we're not moving back into a state we have been in before */
302                 if (e_state_set != NULL) {
303                         EdRotState e_state_alt;
304                         erot_state_alternate(e, &e_state_alt);
305                         if (BLI_gset_haskey(e_state_set, (void *)&e_state_alt)) {
306                                 // printf("  skipping, we already have this state\n");
307                                 return;
308                         }
309                 }
310
311                 {
312                         /* recalculate edge */
313                         const float cost = bm_edge_calc_rotate_beauty(e, flag, method);
314                         if (cost < 0.0f) {
315                                 eheap_table[i] = BLI_heap_insert(eheap, cost, e);
316                         }
317                         else {
318                                 eheap_table[i] = NULL;
319                         }
320                 }
321         }
322 }
323
324 /* we have rotated an edge, tag other edges and clear this one */
325 static void bm_edge_update_beauty_cost(
326         BMEdge *e, Heap *eheap, HeapNode **eheap_table, GSet **edge_state_arr,
327         const BMEdge **edge_array, const int edge_array_len,
328         /* only for testing the edge is in the array */
329         const short flag, const short method)
330 {
331         int i;
332
333         BMEdge *e_arr[4] = {
334             e->l->next->e,
335             e->l->prev->e,
336             e->l->radial_next->next->e,
337             e->l->radial_next->prev->e,
338         };
339
340         BLI_assert(e->l->f->len == 3 &&
341                    e->l->radial_next->f->len == 3);
342
343         BLI_assert(BM_edge_face_count_is_equal(e, 2));
344
345         for (i = 0; i < 4; i++) {
346                 bm_edge_update_beauty_cost_single(
347                         e_arr[i],
348                         eheap, eheap_table, edge_state_arr,
349                         edge_array, edge_array_len,
350                         flag, method);
351         }
352 }
353
354 /* -------------------------------------------------------------------- */
355 /* Beautify Fill */
356
357 /**
358  * \note This function sets the edge indices to invalid values.
359  */
360 void BM_mesh_beautify_fill(
361         BMesh *bm, BMEdge **edge_array, const int edge_array_len,
362         const short flag, const short method,
363         const short oflag_edge, const short oflag_face)
364 {
365         Heap *eheap;             /* edge heap */
366         HeapNode **eheap_table;  /* edge index aligned table pointing to the eheap */
367
368         GSet       **edge_state_arr  = MEM_callocN((size_t)edge_array_len * sizeof(GSet *), __func__);
369         BLI_mempool *edge_state_pool = BLI_mempool_create(sizeof(EdRotState), 0, 512, BLI_MEMPOOL_NOP);
370         int i;
371
372 #ifdef DEBUG_TIME
373         TIMEIT_START(beautify_fill);
374 #endif
375
376         eheap = BLI_heap_new_ex((uint)edge_array_len);
377         eheap_table = MEM_mallocN(sizeof(HeapNode *) * (size_t)edge_array_len, __func__);
378
379         /* build heap */
380         for (i = 0; i < edge_array_len; i++) {
381                 BMEdge *e = edge_array[i];
382                 const float cost = bm_edge_calc_rotate_beauty(e, flag, method);
383                 if (cost < 0.0f) {
384                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
385                 }
386                 else {
387                         eheap_table[i] = NULL;
388                 }
389
390                 BM_elem_index_set(e, i);  /* set_dirty */
391         }
392         bm->elem_index_dirty |= BM_EDGE;
393
394         while (BLI_heap_is_empty(eheap) == false) {
395                 BMEdge *e = BLI_heap_popmin(eheap);
396                 i = BM_elem_index_get(e);
397                 eheap_table[i] = NULL;
398
399                 BLI_assert(BM_edge_face_count_is_equal(e, 2));
400
401                 e = BM_edge_rotate(bm, e, false, BM_EDGEROT_CHECK_EXISTS);
402
403                 BLI_assert(e == NULL || BM_edge_face_count_is_equal(e, 2));
404
405                 if (LIKELY(e)) {
406                         GSet *e_state_set = edge_state_arr[i];
407
408                         /* add the new state into the set so we don't move into this state again
409                          * note: we could add the previous state too but this isn't essential)
410                          *       for avoiding eternal loops */
411                         EdRotState *e_state = BLI_mempool_alloc(edge_state_pool);
412                         erot_state_current(e, e_state);
413                         if (UNLIKELY(e_state_set == NULL)) {
414                                 edge_state_arr[i] = e_state_set = erot_gset_new();  /* store previous state */
415                         }
416                         BLI_assert(BLI_gset_haskey(e_state_set, (void *)e_state) == false);
417                         BLI_gset_insert(e_state_set, e_state);
418
419
420                         // printf("  %d -> %d, %d\n", i, BM_elem_index_get(e->v1), BM_elem_index_get(e->v2));
421
422                         /* maintain the index array */
423                         edge_array[i] = e;
424                         BM_elem_index_set(e, i);
425
426                         /* recalculate faces connected on the heap */
427                         bm_edge_update_beauty_cost(e, eheap, eheap_table, edge_state_arr,
428                                                    (const BMEdge **)edge_array, edge_array_len,
429                                                    flag, method);
430
431                         /* update flags */
432                         if (oflag_edge) {
433                                 BMO_edge_flag_enable(bm, e, oflag_edge);
434                         }
435
436                         if (oflag_face) {
437                                 BMO_face_flag_enable(bm, e->l->f, oflag_face);
438                                 BMO_face_flag_enable(bm, e->l->radial_next->f, oflag_face);
439                         }
440                 }
441         }
442
443         BLI_heap_free(eheap, NULL);
444         MEM_freeN(eheap_table);
445
446         for (i = 0; i < edge_array_len; i++) {
447                 if (edge_state_arr[i]) {
448                         BLI_gset_free(edge_state_arr[i], NULL);
449                 }
450         }
451
452         MEM_freeN(edge_state_arr);
453         BLI_mempool_destroy(edge_state_pool);
454
455 #ifdef DEBUG_TIME
456         TIMEIT_END(beautify_fill);
457 #endif
458 }
459