b563e286a48c55b40079e9f3c10c216a78758001
[blender.git] / source / blender / blenlib / intern / polyfill2d_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  * ***** END GPL LICENSE BLOCK *****
19  */
20
21 /** \file blender/blenlib/intern/polyfill2d_beautify.c
22  *  \ingroup bli
23  *
24  * This function is to improve the tessellation resulting from polyfill2d,
25  * creating optimal topology.
26  *
27  * The functionality here matches #BM_mesh_beautify_fill,
28  * but its far simpler to perform this operation in 2d,
29  * on a simple polygon representation where we _know_:
30  *
31  * - The polygon is primitive with no holes with a continuous boundary.
32  * - Tris have consistent winding.
33  * - 2d (saves some hassles projecting face pairs on an axis for every edge-rotation)
34  *   also saves us having to store all previous edge-states (see #EdRotState in bmesh_beautify.c)
35  *
36  * \note
37  *
38  * No globals - keep threadsafe.
39  */
40
41 #include "BLI_utildefines.h"
42 #include "BLI_math.h"
43
44 #include "BLI_memarena.h"
45 #include "BLI_edgehash.h"
46 #include "BLI_heap.h"
47
48 #include "BLI_polyfill2d_beautify.h"  /* own include */
49
50 #include "BLI_strict_flags.h"
51
52 struct PolyEdge {
53         /** ordered vert indices (smaller first) */
54         unsigned int verts[2];
55         /** ordered face indices (depends on winding compared to the edge verts)
56          * - (verts[0], verts[1])  == faces[0]
57          * - (verts[1], verts[0])  == faces[1]
58          */
59         unsigned int faces[2];
60         /**
61          * The face-index which isn't used by either of the edges verts [0 - 2].
62          * could be calculated each time, but cleaner to store for reuse.
63          */
64         unsigned int faces_other_v[2];
65 };
66
67
68 #ifndef NDEBUG
69 /**
70  * Only to check for error-cases.
71  */
72 static void polyfill_validate_tri(unsigned int (*tris)[3], unsigned int tri_index, EdgeHash *ehash)
73 {
74         const unsigned int *tri = tris[tri_index];
75         int j_curr;
76
77         BLI_assert(!ELEM(tri[0], tri[1], tri[2]) &&
78                    !ELEM(tri[1], tri[0], tri[2]) &&
79                    !ELEM(tri[2], tri[0], tri[1]));
80
81         for (j_curr = 0; j_curr < 3; j_curr++) {
82                 struct PolyEdge *e;
83                 unsigned int e_v1 = tri[(j_curr    )    ];
84                 unsigned int e_v2 = tri[(j_curr + 1) % 3];
85                 e = BLI_edgehash_lookup(ehash, e_v1, e_v2);
86                 if (e) {
87                         if (e->faces[0] == tri_index) {
88                                 BLI_assert(e->verts[0] == e_v1);
89                                 BLI_assert(e->verts[1] == e_v2);
90                         }
91                         else if (e->faces[1] == tri_index) {
92                                 BLI_assert(e->verts[0] == e_v2);
93                                 BLI_assert(e->verts[1] == e_v1);
94                         }
95                         else {
96                                 BLI_assert(0);
97                         }
98
99                         BLI_assert(e->faces[0] != e->faces[1]);
100                         BLI_assert(ELEM(e_v1, UNPACK3(tri)));
101                         BLI_assert(ELEM(e_v2, UNPACK3(tri)));
102                         BLI_assert(ELEM(e_v1, UNPACK2(e->verts)));
103                         BLI_assert(ELEM(e_v2, UNPACK2(e->verts)));
104                         BLI_assert(e_v1 != tris[e->faces[0]][e->faces_other_v[0]]);
105                         BLI_assert(e_v1 != tris[e->faces[1]][e->faces_other_v[1]]);
106                         BLI_assert(e_v2 != tris[e->faces[0]][e->faces_other_v[0]]);
107                         BLI_assert(e_v2 != tris[e->faces[1]][e->faces_other_v[1]]);
108
109                         BLI_assert(ELEM(tri_index, UNPACK2(e->faces)));
110                 }
111         }
112 }
113 #endif
114
115 BLI_INLINE bool is_boundary_edge(unsigned int i_a, unsigned int i_b, const unsigned int coord_last)
116 {
117         BLI_assert(i_a < i_b);
118         return ((i_a + 1 == i_b) || UNLIKELY((i_a == 0) && (i_b == coord_last)));
119 }
120 /**
121  * Assuming we have 2 triangles sharing an edge (2 - 4),
122  * check if the edge running from (1 - 3) gives better results.
123  *
124  * \param lock_degenerate: Use to avoid rotating out of a degenerate state.
125  * - When true, an existing zero area face on either side of the (2 - 4) split will return a positive value.
126  * - When false, the check must be non-biased towards either split direction.
127  *
128  * \return (negative number means the edge can be rotated, lager == better).
129  */
130 float BLI_polyfill_beautify_quad_rotate_calc_ex(
131         const float v1[2], const float v2[2], const float v3[2], const float v4[2],
132         const bool lock_degenerate)
133 {
134         /* not a loop (only to be able to break out) */
135         do {
136                 const float area_2x_234 = cross_tri_v2(v2, v3, v4);
137                 const float area_2x_241 = cross_tri_v2(v2, v4, v1);
138
139                 const float area_2x_123 = cross_tri_v2(v1, v2, v3);
140                 const float area_2x_134 = cross_tri_v2(v1, v3, v4);
141
142                 BLI_assert((ELEM(v1, v2, v3, v4) == false) &&
143                            (ELEM(v2, v1, v3, v4) == false) &&
144                            (ELEM(v3, v1, v2, v4) == false) &&
145                            (ELEM(v4, v1, v2, v3) == false));
146
147                 /*
148                  * Test for unusable (1-3) state.
149                  * - Area sign flipping to check faces aren't going to point in opposite directions.
150                  * - Area epsilon check that the one of the faces won't be zero area.
151                  */
152                 if ((area_2x_123 >= 0.0f) != (area_2x_134 >= 0.0f)) {
153                         break;
154                 }
155                 else if ((fabsf(area_2x_123) <= FLT_EPSILON) || (fabsf(area_2x_134) <= FLT_EPSILON)) {
156                         break;
157                 }
158
159                 /* Test for unusable (2-4) state (same as above). */
160                 if ((area_2x_234 >= 0.0f) != (area_2x_241 >= 0.0f)) {
161                         if (lock_degenerate) {
162                                 break;
163                         }
164                         else {
165                                 return -FLT_MAX;  /* always rotate */
166                         }
167                 }
168                 else if ((fabsf(area_2x_234) <= FLT_EPSILON) || (fabsf(area_2x_241) <= FLT_EPSILON)) {
169                         return -FLT_MAX;  /* always rotate */
170                 }
171
172                 {
173                         /* testing rule: the area divided by the perimeter,
174                          * check if (1-3) beats the existing (2-4) edge rotation */
175                         float area_a, area_b;
176                         float prim_a, prim_b;
177                         float fac_24, fac_13;
178
179                         float len_12, len_23, len_34, len_41, len_24, len_13;
180
181                         /* edges around the quad */
182                         len_12 = len_v2v2(v1, v2);
183                         len_23 = len_v2v2(v2, v3);
184                         len_34 = len_v2v2(v3, v4);
185                         len_41 = len_v2v2(v4, v1);
186                         /* edges crossing the quad interior */
187                         len_13 = len_v2v2(v1, v3);
188                         len_24 = len_v2v2(v2, v4);
189
190                         /* note, area is in fact (area * 2),
191                          * but in this case its OK, since we're comparing ratios */
192
193                         /* edge (2-4), current state */
194                         area_a = fabsf(area_2x_234);
195                         area_b = fabsf(area_2x_241);
196                         prim_a = len_23 + len_34 + len_24;
197                         prim_b = len_41 + len_12 + len_24;
198                         fac_24 = (area_a / prim_a) + (area_b / prim_b);
199
200                         /* edge (1-3), new state */
201                         area_a = fabsf(area_2x_123);
202                         area_b = fabsf(area_2x_134);
203                         prim_a = len_12 + len_23 + len_13;
204                         prim_b = len_34 + len_41 + len_13;
205                         fac_13 = (area_a / prim_a) + (area_b / prim_b);
206
207                         /* negative number if (1-3) is an improved state */
208                         return fac_24 - fac_13;
209                 }
210         } while (false);
211
212         return FLT_MAX;
213 }
214
215 static float polyedge_rotate_beauty_calc(
216         const float (*coords)[2],
217         const unsigned int (*tris)[3],
218         const struct PolyEdge *e)
219 {
220         const float *v1, *v2, *v3, *v4;
221
222         v1 = coords[tris[e->faces[0]][e->faces_other_v[0]]];
223         v3 = coords[tris[e->faces[1]][e->faces_other_v[1]]];
224         v2 = coords[e->verts[0]];
225         v4 = coords[e->verts[1]];
226
227         return BLI_polyfill_beautify_quad_rotate_calc(v1, v2, v3, v4);
228 }
229
230 static void polyedge_beauty_cost_update_single(
231         const float (*coords)[2],
232         const unsigned int (*tris)[3],
233         const struct PolyEdge *edges,
234         struct PolyEdge *e,
235         Heap *eheap, HeapNode **eheap_table)
236 {
237         const unsigned int i = (unsigned int)(e - edges);
238
239         if (eheap_table[i]) {
240                 BLI_heap_remove(eheap, eheap_table[i]);
241                 eheap_table[i] = NULL;
242         }
243
244         {
245                 /* recalculate edge */
246                 const float cost = polyedge_rotate_beauty_calc(coords, tris, e);
247                 /* We can get cases where both choices generate very small negative costs, which leads to infinite loop.
248                  * Anyway, costs above that are not worth recomputing, maybe we could even optimize it to a smaller limit?
249                  * Actually, FLT_EPSILON is too small in some cases, 1e-6f seems to work OK hopefully?
250                  * See T43578, T49478. */
251                 if (cost < -1e-6f) {
252                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
253                 }
254                 else {
255                         eheap_table[i] = NULL;
256                 }
257         }
258 }
259
260 static void polyedge_beauty_cost_update(
261         const float (*coords)[2],
262         const unsigned int (*tris)[3],
263         const struct PolyEdge *edges,
264         struct PolyEdge *e,
265         Heap *eheap, HeapNode **eheap_table,
266         EdgeHash *ehash)
267 {
268         const unsigned int *tri_0 = tris[e->faces[0]];
269         const unsigned int *tri_1 = tris[e->faces[1]];
270         unsigned int i;
271
272         struct PolyEdge *e_arr[4] = {
273                 BLI_edgehash_lookup(ehash,
274                         tri_0[(e->faces_other_v[0]    ) % 3],
275                         tri_0[(e->faces_other_v[0] + 1) % 3]),
276                 BLI_edgehash_lookup(ehash,
277                         tri_0[(e->faces_other_v[0] + 2) % 3],
278                         tri_0[(e->faces_other_v[0]    ) % 3]),
279                 BLI_edgehash_lookup(ehash,
280                         tri_1[(e->faces_other_v[1]    ) % 3],
281                         tri_1[(e->faces_other_v[1] + 1) % 3]),
282                 BLI_edgehash_lookup(ehash,
283                         tri_1[(e->faces_other_v[1] + 2) % 3],
284                         tri_1[(e->faces_other_v[1]    ) % 3]),
285         };
286
287
288         for (i = 0; i < 4; i++) {
289                 if (e_arr[i]) {
290                         BLI_assert(!(ELEM(e_arr[i]->faces[0], UNPACK2(e->faces)) &&
291                                      ELEM(e_arr[i]->faces[1], UNPACK2(e->faces))));
292
293                         polyedge_beauty_cost_update_single(
294                                 coords, tris, edges,
295                                 e_arr[i],
296                                 eheap, eheap_table);
297                 }
298         }
299 }
300
301 static void polyedge_rotate(
302         unsigned int (*tris)[3],
303         struct PolyEdge *e,
304         EdgeHash *ehash)
305 {
306         unsigned int e_v1_new = tris[e->faces[0]][e->faces_other_v[0]];
307         unsigned int e_v2_new = tris[e->faces[1]][e->faces_other_v[1]];
308
309 #ifndef NDEBUG
310         polyfill_validate_tri(tris, e->faces[0], ehash);
311         polyfill_validate_tri(tris, e->faces[1], ehash);
312 #endif
313
314         BLI_assert(e_v1_new != e_v2_new);
315         BLI_assert(!ELEM(e_v2_new, UNPACK3(tris[e->faces[0]])));
316         BLI_assert(!ELEM(e_v1_new, UNPACK3(tris[e->faces[1]])));
317
318         tris[e->faces[0]][(e->faces_other_v[0] + 1) % 3] = e_v2_new;
319         tris[e->faces[1]][(e->faces_other_v[1] + 1) % 3] = e_v1_new;
320
321         e->faces_other_v[0] = (e->faces_other_v[0] + 2) % 3;
322         e->faces_other_v[1] = (e->faces_other_v[1] + 2) % 3;
323
324         BLI_assert((tris[e->faces[0]][e->faces_other_v[0]] != e_v1_new) &&
325                    (tris[e->faces[0]][e->faces_other_v[0]] != e_v2_new));
326         BLI_assert((tris[e->faces[1]][e->faces_other_v[1]] != e_v1_new) &&
327                    (tris[e->faces[1]][e->faces_other_v[1]] != e_v2_new));
328
329         BLI_edgehash_remove(ehash, e->verts[0], e->verts[1], NULL);
330         BLI_edgehash_insert(ehash, e_v1_new, e_v2_new, e);
331
332         if (e_v1_new < e_v2_new) {
333                 e->verts[0] = e_v1_new;
334                 e->verts[1] = e_v2_new;
335         }
336         else {
337                 /* maintain winding info */
338                 e->verts[0] = e_v2_new;
339                 e->verts[1] = e_v1_new;
340
341                 SWAP(unsigned int, e->faces[0], e->faces[1]);
342                 SWAP(unsigned int, e->faces_other_v[0], e->faces_other_v[1]);
343         }
344
345         /* update adjacent data */
346         {
347                 unsigned int e_side = 0;
348
349                 for (e_side = 0; e_side < 2; e_side++) {
350                         /* 't_other' which we need to swap out is always the same edge-order */
351                         const unsigned int t_other = (((e->faces_other_v[e_side]) + 2)) % 3;
352                         unsigned int t_index = e->faces[e_side];
353                         unsigned int t_index_other = e->faces[!e_side];
354                         unsigned int *tri = tris[t_index];
355
356                         struct PolyEdge *e_other;
357                         unsigned int e_v1 = tri[(t_other    )    ];
358                         unsigned int e_v2 = tri[(t_other + 1) % 3];
359
360                         e_other = BLI_edgehash_lookup(ehash, e_v1, e_v2);
361                         if (e_other) {
362                                 BLI_assert(t_index != e_other->faces[0] && t_index != e_other->faces[1]);
363                                 if (t_index_other == e_other->faces[0]) {
364                                         e_other->faces[0] = t_index;
365                                         e_other->faces_other_v[0] = (t_other + 2) % 3;
366                                         BLI_assert(!ELEM(tri[e_other->faces_other_v[0]], e_v1, e_v2));
367                                 }
368                                 else if (t_index_other == e_other->faces[1]) {
369                                         e_other->faces[1] = t_index;
370                                         e_other->faces_other_v[1] = (t_other + 2) % 3;
371                                         BLI_assert(!ELEM(tri[e_other->faces_other_v[1]], e_v1, e_v2));
372                                 }
373                                 else {
374                                         BLI_assert(0);
375                                 }
376                         }
377                 }
378         }
379
380 #ifndef NDEBUG
381         polyfill_validate_tri(tris, e->faces[0], ehash);
382         polyfill_validate_tri(tris, e->faces[1], ehash);
383 #endif
384
385         BLI_assert(!ELEM(tris[e->faces[0]][e->faces_other_v[0]], UNPACK2(e->verts)));
386         BLI_assert(!ELEM(tris[e->faces[1]][e->faces_other_v[1]], UNPACK2(e->verts)));
387 }
388
389 /**
390  * The intention is that this calculates the output of #BLI_polyfill_calc
391  *
392  *
393  * \note assumes the \a coords form a boundary,
394  * so any edges running along contiguous (wrapped) indices,
395  * are ignored since the edges wont share 2 faces.
396  */
397 void BLI_polyfill_beautify(
398         const float (*coords)[2],
399         const unsigned int coords_tot,
400         unsigned int (*tris)[3],
401
402         /* structs for reuse */
403         MemArena *arena, Heap *eheap, EdgeHash *ehash)
404 {
405         const unsigned int coord_last = coords_tot - 1;
406         const unsigned int tris_tot = coords_tot - 2;
407         /* internal edges only (between 2 tris) */
408         const unsigned int edges_tot = tris_tot - 1;
409         unsigned int edges_tot_used = 0;
410         unsigned int i;
411
412         HeapNode **eheap_table;
413
414         struct PolyEdge *edges = BLI_memarena_alloc(arena, edges_tot * sizeof(*edges));
415
416         BLI_assert(BLI_heap_size(eheap) == 0);
417         BLI_assert(BLI_edgehash_size(ehash) == 0);
418
419         /* first build edges */
420         for (i = 0; i < tris_tot; i++) {
421                 unsigned int j_prev, j_curr, j_next;
422                 j_prev = 2;
423                 j_next = 1;
424                 for (j_curr = 0; j_curr < 3; j_next = j_prev, j_prev = j_curr++) {
425                         int e_index;
426
427                         unsigned int e_pair[2] = {
428                                 tris[i][j_prev],
429                                 tris[i][j_curr],
430                         };
431
432                         if (e_pair[0] > e_pair[1]) {
433                                 SWAP(unsigned int, e_pair[0], e_pair[1]);
434                                 e_index = 1;
435                         }
436                         else {
437                                 e_index = 0;
438                         }
439
440                         if (!is_boundary_edge(e_pair[0], e_pair[1], coord_last)) {
441                                 struct PolyEdge *e;
442                                 void **val_p;
443
444                                 if (!BLI_edgehash_ensure_p(ehash, e_pair[0], e_pair[1], &val_p)) {
445                                         e = &edges[edges_tot_used++];
446                                         *val_p = e;
447                                         memcpy(e->verts, e_pair, sizeof(e->verts));
448 #ifndef NDEBUG
449                                         e->faces[!e_index] = (unsigned int)-1;
450 #endif
451                                 }
452                                 else {
453                                         e = *val_p;
454                                         /* ensure each edge only ever has 2x users */
455 #ifndef NDEBUG
456                                         BLI_assert(e->faces[e_index] == (unsigned int)-1);
457                                         BLI_assert((e->verts[0] == e_pair[0]) &&
458                                                    (e->verts[1] == e_pair[1]));
459 #endif
460                                 }
461
462                                 e->faces[e_index] = i;
463                                 e->faces_other_v[e_index] = j_next;
464                         }
465                 }
466         }
467
468         /* now perform iterative rotations */
469         eheap_table = BLI_memarena_alloc(arena, sizeof(HeapNode *) * (size_t)edges_tot);
470
471         // for (i = 0; i < tris_tot; i++) { polyfill_validate_tri(tris, i, eh); }
472
473         /* build heap */
474         for (i = 0; i < edges_tot; i++) {
475                 struct PolyEdge *e = &edges[i];
476                 const float cost = polyedge_rotate_beauty_calc(coords, (const unsigned int (*)[3])tris, e);
477                 if (cost < 0.0f) {
478                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
479                 }
480                 else {
481                         eheap_table[i] = NULL;
482                 }
483         }
484
485         while (BLI_heap_is_empty(eheap) == false) {
486                 struct PolyEdge *e = BLI_heap_popmin(eheap);
487                 i = (unsigned int)(e - edges);
488                 eheap_table[i] = NULL;
489
490                 polyedge_rotate(tris, e, ehash);
491
492                 /* recalculate faces connected on the heap */
493                 polyedge_beauty_cost_update(
494                         coords, (const unsigned int (*)[3])tris, edges,
495                         e,
496                         eheap, eheap_table, ehash);
497         }
498
499         BLI_heap_clear(eheap, NULL);
500         BLI_edgehash_clear_ex(ehash, NULL, BLI_POLYFILL_ALLOC_NGON_RESERVE);
501
502         /* MEM_freeN(eheap_table); */  /* arena */
503 }