Fix for triangulate and beauty-fill
[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
41 #include "MEM_guardedalloc.h"
42
43 #include "bmesh.h"
44 #include "bmesh_beautify.h"  /* own include */
45
46
47 // #define DEBUG_TIME
48
49 #ifdef DEBUG_TIME
50 #  include "PIL_time.h"
51 #  include "PIL_time_utildefines.h"
52 #endif
53
54 /* -------------------------------------------------------------------- */
55 /* GSet for edge rotation */
56
57 typedef struct EdRotState {
58         int v1, v2; /*  edge vert, small -> large */
59         int f1, f2; /*  face vert, small -> large */
60 } EdRotState;
61
62 static unsigned int erot_gsetutil_hash(const void *ptr)
63 {
64         const EdRotState *e_state = (const EdRotState *)ptr;
65         unsigned int
66         hash  = BLI_ghashutil_inthash(SET_INT_IN_POINTER(e_state->v1));
67         hash ^= BLI_ghashutil_inthash(SET_INT_IN_POINTER(e_state->v2));
68         hash ^= BLI_ghashutil_inthash(SET_INT_IN_POINTER(e_state->f1));
69         hash ^= BLI_ghashutil_inthash(SET_INT_IN_POINTER(e_state->f2));
70         return hash;
71 }
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
87 static GSet *erot_gset_new(void)
88 {
89         return BLI_gset_new(erot_gsetutil_hash, erot_gsetutil_cmp, __func__);
90 }
91
92 /* ensure v0 is smaller */
93 #define EDGE_ORD(v0, v1) \
94         if (v0 > v1) {       \
95                 v0 ^= v1;        \
96                 v1 ^= v0;        \
97                 v0 ^= v1;        \
98         } (void)0
99
100 static void erot_state_ex(const BMEdge *e, int v_index[2], int f_index[2])
101 {
102         BLI_assert(BM_edge_is_manifold((BMEdge *)e));
103         BLI_assert(BM_vert_in_edge(e, e->l->prev->v)              == false);
104         BLI_assert(BM_vert_in_edge(e, e->l->radial_next->prev->v) == false);
105
106         /* verts of the edge */
107         v_index[0] = BM_elem_index_get(e->v1);
108         v_index[1] = BM_elem_index_get(e->v2);
109         EDGE_ORD(v_index[0], v_index[1]);
110
111         /* verts of each of the 2 faces attached to this edge
112          * (that are not apart of this edge) */
113         f_index[0] = BM_elem_index_get(e->l->prev->v);
114         f_index[1] = BM_elem_index_get(e->l->radial_next->prev->v);
115         EDGE_ORD(f_index[0], f_index[1]);
116 }
117
118 static void erot_state_current(const BMEdge *e, EdRotState *e_state)
119 {
120         erot_state_ex(e, &e_state->v1, &e_state->f1);
121 }
122
123 static void erot_state_alternate(const BMEdge *e, EdRotState *e_state)
124 {
125         erot_state_ex(e, &e_state->f1, &e_state->v1);
126 }
127
128 /* -------------------------------------------------------------------- */
129 /* Calculate the improvement of rotating the edge */
130
131 /**
132  * \return a negative value means the edge can be rotated.
133  */
134 static float bm_edge_calc_rotate_beauty__area(
135         const float v1[3], const float v2[3], const float v3[3], const float v4[3])
136 {
137         /* not a loop (only to be able to break out) */
138         do {
139                 float v1_xy[2], v2_xy[2], v3_xy[2], v4_xy[2];
140
141                 /* first get the 2d values */
142                 {
143                         bool is_zero_a, is_zero_b;
144                         float no[3];
145                         float axis_mat[3][3];
146
147                         // printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
148                         BLI_assert((ELEM3(v1, v2, v3, v4) == false) &&
149                                    (ELEM3(v2, v1, v3, v4) == false) &&
150                                    (ELEM3(v3, v1, v2, v4) == false) &&
151                                    (ELEM3(v4, v1, v2, v3) == false));
152
153                         is_zero_a = area_tri_v3(v2, v3, v4) <= FLT_EPSILON;
154                         is_zero_b = area_tri_v3(v2, v4, v1) <= FLT_EPSILON;
155
156                         if (LIKELY(is_zero_a == false && is_zero_b == false)) {
157                                 float no_a[3], no_b[3];
158                                 normal_tri_v3(no_a, v2, v3, v4);  /* a */
159                                 normal_tri_v3(no_b, v2, v4, v1);  /* b */
160                                 add_v3_v3v3(no, no_a, no_b);
161                                 if (UNLIKELY(normalize_v3(no) <= FLT_EPSILON)) {
162                                         break;
163                                 }
164                         }
165                         else if (is_zero_a == false) {
166                                 normal_tri_v3(no, v2, v3, v4);  /* a */
167                         }
168                         else if (is_zero_b == false) {
169                                 normal_tri_v3(no, v2, v4, v1);  /* b */
170                         }
171                         else {
172                                 /* both zero area, no useful normal can be calculated */
173                                 break;
174                         }
175
176                         // { float a = angle_normalized_v3v3(no_a, no_b); printf("~ %.7f\n", a); fflush(stdout);}
177
178                         axis_dominant_v3_to_m3(axis_mat, no);
179                         mul_v2_m3v3(v1_xy, axis_mat, v1);
180                         mul_v2_m3v3(v2_xy, axis_mat, v2);
181                         mul_v2_m3v3(v3_xy, axis_mat, v3);
182                         mul_v2_m3v3(v4_xy, axis_mat, v4);
183                 }
184
185                 // printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
186
187                 if (is_quad_convex_v2(v1_xy, v2_xy, v3_xy, v4_xy)) {
188                         /* testing rule: the area divided by the perimeter,
189                          * check if (1-3) beats the existing (2-4) edge rotation */
190                         float area_a, area_b;
191                         float prim_a, prim_b;
192                         float fac_24, fac_13;
193
194                         float len_12, len_23, len_34, len_41, len_24, len_13;
195
196                         /* edges around the quad */
197                         len_12 = len_v2v2(v1_xy, v2_xy);
198                         len_23 = len_v2v2(v2_xy, v3_xy);
199                         len_34 = len_v2v2(v3_xy, v4_xy);
200                         len_41 = len_v2v2(v4_xy, v1_xy);
201                         /* edges crossing the quad interior */
202                         len_13 = len_v2v2(v1_xy, v3_xy);
203                         len_24 = len_v2v2(v2_xy, v4_xy);
204
205                         /* edge (2-4), current state */
206                         area_a = area_tri_v2(v2_xy, v3_xy, v4_xy);
207                         area_b = area_tri_v2(v2_xy, v4_xy, v1_xy);
208                         prim_a = len_23 + len_34 + len_24;
209                         prim_b = len_24 + len_41 + len_12;
210                         fac_24 = (area_a / prim_a) + (area_b / prim_b);
211
212                         /* edge (1-3), new state */
213                         area_a = area_tri_v2(v1_xy, v2_xy, v3_xy);
214                         area_b = area_tri_v2(v1_xy, v3_xy, v4_xy);
215                         prim_a = len_12 + len_23 + len_13;
216                         prim_b = len_34 + len_41 + len_13;
217                         fac_13 = (area_a / prim_a) + (area_b / prim_b);
218
219                         /* negative number if (1-3) is an improved state */
220                         return fac_24 - fac_13;
221                 }
222         } while (false);
223
224         return FLT_MAX;
225 }
226
227 static float bm_edge_calc_rotate_beauty__angle(
228         const float v1[3], const float v2[3], const float v3[3], const float v4[3])
229 {
230         /* not a loop (only to be able to break out) */
231         do {
232                 float no_a[3], no_b[3];
233                 float angle_24, angle_13;
234
235                 /* edge (2-4), current state */
236                 normal_tri_v3(no_a, v2, v3, v4);
237                 normal_tri_v3(no_b, v2, v4, v1);
238                 angle_24 = angle_normalized_v3v3(no_a, no_b);
239
240                 /* edge (1-3), new state */
241                 /* only check new state for degenerate outcome */
242                 if ((normal_tri_v3(no_a, v1, v2, v3) == 0.0f) ||
243                     (normal_tri_v3(no_b, v1, v3, v4) == 0.0f))
244                 {
245                         break;
246                 }
247                 angle_13 = angle_normalized_v3v3(no_a, no_b);
248
249                 return angle_13 - angle_24;
250         } while (false);
251
252         return FLT_MAX;
253 }
254
255 float BM_verts_calc_rotate_beauty(
256 const BMVert *v1, const BMVert *v2, const BMVert *v3, const BMVert *v4, const short flag, const short method)
257 {
258         /* not a loop (only to be able to break out) */
259         do {
260                 if (flag & VERT_RESTRICT_TAG) {
261                         const BMVert *v_a = v1, *v_b = v3;
262                         if (BM_elem_flag_test(v_a, BM_ELEM_TAG) ==  BM_elem_flag_test(v_b, BM_ELEM_TAG)) {
263                                 break;
264                         }
265                 }
266
267                 if (UNLIKELY(v1 == v3)) {
268                         // printf("This should never happen, but does sometimes!\n");
269                         break;
270                 }
271
272                 switch (method) {
273                         case 0:
274                                 return bm_edge_calc_rotate_beauty__area(v1->co, v2->co, v3->co, v4->co);
275                         default:
276                                 return bm_edge_calc_rotate_beauty__angle(v1->co, v2->co, v3->co, v4->co);
277                 }
278         } while (false);
279
280         return FLT_MAX;
281 }
282
283 static float bm_edge_calc_rotate_beauty(const BMEdge *e, const short flag, const short method)
284 {
285         const BMVert *v1, *v2, *v3, *v4;
286         v1 = e->l->prev->v;               /* first vert co */
287         v2 = e->l->v;                     /* e->v1 or e->v2*/
288         v3 = e->l->radial_next->prev->v;  /* second vert co */
289         v4 = e->l->next->v;               /* e->v1 or e->v2*/
290
291         return BM_verts_calc_rotate_beauty(v1, v2, v3, v4, flag, method);
292 }
293
294 /* -------------------------------------------------------------------- */
295 /* Update the edge cost of rotation in the heap */
296
297 BLI_INLINE bool edge_in_array(const BMEdge *e, const BMEdge **edge_array, const int edge_array_len)
298 {
299         const int index = BM_elem_index_get(e);
300         return ((index >= 0) &&
301                 (index < edge_array_len) &&
302                 (e == edge_array[index]));
303 }
304
305 /* recalc an edge in the heap (surrounding geometry has changed) */
306 static void bm_edge_update_beauty_cost_single(BMEdge *e, Heap *eheap, HeapNode **eheap_table, GSet **edge_state_arr,
307                                               /* only for testing the edge is in the array */
308                                               const BMEdge **edge_array, const int edge_array_len,
309
310                                               const short flag, const short method)
311 {
312         if (edge_in_array(e, edge_array, edge_array_len)) {
313                 const int i = BM_elem_index_get(e);
314                 GSet *e_state_set = edge_state_arr[i];
315
316                 if (eheap_table[i]) {
317                         BLI_heap_remove(eheap, eheap_table[i]);
318                         eheap_table[i] = NULL;
319                 }
320
321                 /* check if we can add it back */
322                 BLI_assert(BM_edge_is_manifold(e) == true);
323
324                 /* check we're not moving back into a state we have been in before */
325                 if (e_state_set != NULL) {
326                         EdRotState e_state_alt;
327                         erot_state_alternate(e, &e_state_alt);
328                         if (BLI_gset_haskey(e_state_set, (void *)&e_state_alt)) {
329                                 // printf("  skipping, we already have this state\n");
330                                 return;
331                         }
332                 }
333
334                 {
335                         /* recalculate edge */
336                         const float cost = bm_edge_calc_rotate_beauty(e, flag, method);
337                         if (cost < 0.0f) {
338                                 eheap_table[i] = BLI_heap_insert(eheap, cost, e);
339                         }
340                         else {
341                                 eheap_table[i] = NULL;
342                         }
343                 }
344         }
345 }
346
347 /* we have rotated an edge, tag other edges and clear this one */
348 static void bm_edge_update_beauty_cost(BMEdge *e, Heap *eheap, HeapNode **eheap_table, GSet **edge_state_arr,
349                                        const BMEdge **edge_array, const int edge_array_len,
350                                        /* only for testing the edge is in the array */
351                                        const short flag, const short method)
352 {
353         int i;
354
355         BMEdge *e_arr[4] = {
356             e->l->next->e,
357             e->l->prev->e,
358             e->l->radial_next->next->e,
359             e->l->radial_next->prev->e,
360         };
361
362         BLI_assert(e->l->f->len == 3 &&
363                    e->l->radial_next->f->len == 3);
364
365         BLI_assert(BM_edge_face_count(e) == 2);
366
367         for (i = 0; i < 4; i++) {
368                 bm_edge_update_beauty_cost_single(
369                         e_arr[i],
370                         eheap, eheap_table, edge_state_arr,
371                         edge_array, edge_array_len,
372                         flag, method);
373         }
374 }
375
376 /* -------------------------------------------------------------------- */
377 /* Beautify Fill */
378
379 /**
380  * \note This function sets the edge indicies to invalid values.
381  */
382 void BM_mesh_beautify_fill(BMesh *bm, BMEdge **edge_array, const int edge_array_len,
383                                   const short flag, const short method,
384                                   const short oflag_edge, const short oflag_face)
385 {
386         Heap *eheap;             /* edge heap */
387         HeapNode **eheap_table;  /* edge index aligned table pointing to the eheap */
388
389         GSet       **edge_state_arr  = MEM_callocN((size_t)edge_array_len * sizeof(GSet *), __func__);
390         BLI_mempool *edge_state_pool = BLI_mempool_create(sizeof(EdRotState), 512, 512, BLI_MEMPOOL_SYSMALLOC);
391         int i;
392
393 #ifdef DEBUG_TIME
394         TIMEIT_START(beautify_fill);
395 #endif
396
397         eheap = BLI_heap_new_ex((unsigned int)edge_array_len);
398         eheap_table = MEM_mallocN(sizeof(HeapNode *) * (size_t)edge_array_len, __func__);
399
400         /* build heap */
401         for (i = 0; i < edge_array_len; i++) {
402                 BMEdge *e = edge_array[i];
403                 const float cost = bm_edge_calc_rotate_beauty(e, flag, method);
404                 if (cost < 0.0f) {
405                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
406                 }
407                 else {
408                         eheap_table[i] = NULL;
409                 }
410
411                 BM_elem_index_set(e, i);  /* set_dirty */
412         }
413         bm->elem_index_dirty |= BM_EDGE;
414
415         while (BLI_heap_is_empty(eheap) == false) {
416                 BMEdge *e = BLI_heap_popmin(eheap);
417                 i = BM_elem_index_get(e);
418                 eheap_table[i] = NULL;
419
420                 BLI_assert(BM_edge_face_count(e) == 2);
421
422                 e = BM_edge_rotate(bm, e, false, BM_EDGEROT_CHECK_EXISTS);
423
424                 BLI_assert(e == NULL || BM_edge_face_count(e) == 2);
425
426                 if (LIKELY(e)) {
427                         GSet *e_state_set = edge_state_arr[i];
428
429                         /* add the new state into the set so we don't move into this state again
430                          * note: we could add the previous state too but this isn't essential)
431                          *       for avoiding eternal loops */
432                         EdRotState *e_state = BLI_mempool_alloc(edge_state_pool);
433                         erot_state_current(e, e_state);
434                         if (UNLIKELY(e_state_set == NULL)) {
435                                 edge_state_arr[i] = e_state_set = erot_gset_new();  /* store previous state */
436                         }
437                         BLI_assert(BLI_gset_haskey(e_state_set, (void *)e_state) == false);
438                         BLI_gset_insert(e_state_set, e_state);
439
440
441                         // printf("  %d -> %d, %d\n", i, BM_elem_index_get(e->v1), BM_elem_index_get(e->v2));
442
443                         /* maintain the index array */
444                         edge_array[i] = e;
445                         BM_elem_index_set(e, i);
446
447                         /* recalculate faces connected on the heap */
448                         bm_edge_update_beauty_cost(e, eheap, eheap_table, edge_state_arr,
449                                                    (const BMEdge **)edge_array, edge_array_len,
450                                                    flag, method);
451
452                         /* update flags */
453                         if (oflag_edge)
454                                 BMO_elem_flag_enable(bm, e, oflag_edge);
455                         if (oflag_face) {
456                                 BMO_elem_flag_enable(bm, e->l->f, oflag_face);
457                                 BMO_elem_flag_enable(bm, e->l->radial_next->f, oflag_face);
458                         }
459                 }
460         }
461
462         BLI_heap_free(eheap, NULL);
463         MEM_freeN(eheap_table);
464
465         for (i = 0; i < edge_array_len; i++) {
466                 if (edge_state_arr[i]) {
467                         BLI_gset_free(edge_state_arr[i], NULL);
468                 }
469         }
470
471         MEM_freeN(edge_state_arr);
472         BLI_mempool_destroy(edge_state_pool);
473
474 #ifdef DEBUG_TIME
475         TIMEIT_END(beautify_fill);
476 #endif
477 }
478