1f4b598be4ce207e88bf6205f7b8ffb932a4293a
[blender-staging.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  * \return (negative number means the edge can be rotated, lager == better).
125  */
126 float BLI_polyfill_beautify_quad_rotate_calc(
127         const float v1[2], const float v2[2], const float v3[2], const float v4[2])
128 {
129         /* not a loop (only to be able to break out) */
130         do {
131                 bool is_zero_a, is_zero_b;
132
133                 const float area_2x_234 = cross_tri_v2(v2, v3, v4);
134                 const float area_2x_241 = cross_tri_v2(v2, v4, v1);
135
136                 const float area_2x_123 = cross_tri_v2(v1, v2, v3);
137                 const float area_2x_134 = cross_tri_v2(v1, v3, v4);
138
139                 {
140                         BLI_assert((ELEM(v1, v2, v3, v4) == false) &&
141                                    (ELEM(v2, v1, v3, v4) == false) &&
142                                    (ELEM(v3, v1, v2, v4) == false) &&
143                                    (ELEM(v4, v1, v2, v3) == false));
144
145                         is_zero_a = (fabsf(area_2x_234) <= FLT_EPSILON);
146                         is_zero_b = (fabsf(area_2x_241) <= FLT_EPSILON);
147
148                         if (is_zero_a && is_zero_b) {
149                                 break;
150                         }
151                 }
152
153                 /* one of the tri's was degenerate, check we're not rotating
154                  * into a different degenerate shape or flipping the face */
155                 if ((fabsf(area_2x_123) <= FLT_EPSILON) || (fabsf(area_2x_134) <= FLT_EPSILON)) {
156                         /* one of the new rotations is degenerate */
157                         break;
158                 }
159
160                 if ((area_2x_123 >= 0.0f) != (area_2x_134 >= 0.0f)) {
161                         /* rotation would cause flipping */
162                         break;
163                 }
164
165                 {
166                         /* testing rule: the area divided by the perimeter,
167                          * check if (1-3) beats the existing (2-4) edge rotation */
168                         float area_a, area_b;
169                         float prim_a, prim_b;
170                         float fac_24, fac_13;
171
172                         float len_12, len_23, len_34, len_41, len_24, len_13;
173
174                         /* edges around the quad */
175                         len_12 = len_v2v2(v1, v2);
176                         len_23 = len_v2v2(v2, v3);
177                         len_34 = len_v2v2(v3, v4);
178                         len_41 = len_v2v2(v4, v1);
179                         /* edges crossing the quad interior */
180                         len_13 = len_v2v2(v1, v3);
181                         len_24 = len_v2v2(v2, v4);
182
183                         /* note, area is in fact (area * 2),
184                          * but in this case its OK, since we're comparing ratios */
185
186                         /* edge (2-4), current state */
187                         area_a = fabsf(area_2x_234);
188                         area_b = fabsf(area_2x_241);
189                         prim_a = len_23 + len_34 + len_24;
190                         prim_b = len_41 + len_12 + len_24;
191                         fac_24 = (area_a / prim_a) + (area_b / prim_b);
192
193                         /* edge (1-3), new state */
194                         area_a = fabsf(area_2x_123);
195                         area_b = fabsf(area_2x_134);
196                         prim_a = len_12 + len_23 + len_13;
197                         prim_b = len_34 + len_41 + len_13;
198                         fac_13 = (area_a / prim_a) + (area_b / prim_b);
199
200                         /* negative number if (1-3) is an improved state */
201                         return fac_24 - fac_13;
202                 }
203         } while (false);
204
205         return FLT_MAX;
206 }
207
208 static float polyedge_rotate_beauty_calc(
209         const float (*coords)[2],
210         const unsigned int (*tris)[3],
211         const struct PolyEdge *e)
212 {
213         const float *v1, *v2, *v3, *v4;
214
215         v1 = coords[tris[e->faces[0]][e->faces_other_v[0]]];
216         v3 = coords[tris[e->faces[1]][e->faces_other_v[1]]];
217         v2 = coords[e->verts[0]];
218         v4 = coords[e->verts[1]];
219
220         return BLI_polyfill_beautify_quad_rotate_calc(v1, v2, v3, v4);
221 }
222
223 static void polyedge_beauty_cost_update_single(
224         const float (*coords)[2],
225         const unsigned int (*tris)[3],
226         const struct PolyEdge *edges,
227         struct PolyEdge *e,
228         Heap *eheap, HeapNode **eheap_table)
229 {
230         const unsigned int i = (unsigned int)(e - edges);
231
232         if (eheap_table[i]) {
233                 BLI_heap_remove(eheap, eheap_table[i]);
234                 eheap_table[i] = NULL;
235         }
236
237         {
238                 /* recalculate edge */
239                 const float cost = polyedge_rotate_beauty_calc(coords, tris, e);
240                 /* We can get cases where both choices generate very small negative costs, which leads to infinite loop.
241                  * Anyway, costs above that are not worth recomputing, maybe we could even optimze it to a smaller limit?
242                  * See T43578. */
243                 if (cost < -FLT_EPSILON) {
244                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
245                 }
246                 else {
247                         eheap_table[i] = NULL;
248                 }
249         }
250 }
251
252 static void polyedge_beauty_cost_update(
253         const float (*coords)[2],
254         const unsigned int (*tris)[3],
255         const struct PolyEdge *edges,
256         struct PolyEdge *e,
257         Heap *eheap, HeapNode **eheap_table,
258         EdgeHash *ehash)
259 {
260         const unsigned int *tri_0 = tris[e->faces[0]];
261         const unsigned int *tri_1 = tris[e->faces[1]];
262         unsigned int i;
263
264         struct PolyEdge *e_arr[4] = {
265                 BLI_edgehash_lookup(ehash,
266                         tri_0[(e->faces_other_v[0]    ) % 3],
267                         tri_0[(e->faces_other_v[0] + 1) % 3]),
268                 BLI_edgehash_lookup(ehash,
269                         tri_0[(e->faces_other_v[0] + 2) % 3],
270                         tri_0[(e->faces_other_v[0]    ) % 3]),
271                 BLI_edgehash_lookup(ehash,
272                         tri_1[(e->faces_other_v[1]    ) % 3],
273                         tri_1[(e->faces_other_v[1] + 1) % 3]),
274                 BLI_edgehash_lookup(ehash,
275                         tri_1[(e->faces_other_v[1] + 2) % 3],
276                         tri_1[(e->faces_other_v[1]    ) % 3]),
277         };
278
279
280         for (i = 0; i < 4; i++) {
281                 if (e_arr[i]) {
282                         BLI_assert(!(ELEM(e_arr[i]->faces[0], UNPACK2(e->faces)) &&
283                                      ELEM(e_arr[i]->faces[1], UNPACK2(e->faces))));
284
285                         polyedge_beauty_cost_update_single(
286                                 coords, tris, edges,
287                                 e_arr[i],
288                                 eheap, eheap_table);
289                 }
290         }
291 }
292
293 static void polyedge_rotate(
294         unsigned int (*tris)[3],
295         struct PolyEdge *e,
296         EdgeHash *ehash)
297 {
298         unsigned int e_v1_new = tris[e->faces[0]][e->faces_other_v[0]];
299         unsigned int e_v2_new = tris[e->faces[1]][e->faces_other_v[1]];
300
301 #ifndef NDEBUG
302         polyfill_validate_tri(tris, e->faces[0], ehash);
303         polyfill_validate_tri(tris, e->faces[1], ehash);
304 #endif
305
306         BLI_assert(e_v1_new != e_v2_new);
307         BLI_assert(!ELEM(e_v2_new, UNPACK3(tris[e->faces[0]])));
308         BLI_assert(!ELEM(e_v1_new, UNPACK3(tris[e->faces[1]])));
309
310         tris[e->faces[0]][(e->faces_other_v[0] + 1) % 3] = e_v2_new;
311         tris[e->faces[1]][(e->faces_other_v[1] + 1) % 3] = e_v1_new;
312
313         e->faces_other_v[0] = (e->faces_other_v[0] + 2) % 3;
314         e->faces_other_v[1] = (e->faces_other_v[1] + 2) % 3;
315
316         BLI_assert((tris[e->faces[0]][e->faces_other_v[0]] != e_v1_new) &&
317                    (tris[e->faces[0]][e->faces_other_v[0]] != e_v2_new));
318         BLI_assert((tris[e->faces[1]][e->faces_other_v[1]] != e_v1_new) &&
319                    (tris[e->faces[1]][e->faces_other_v[1]] != e_v2_new));
320
321         BLI_edgehash_remove(ehash, e->verts[0], e->verts[1], NULL);
322         BLI_edgehash_insert(ehash, e_v1_new, e_v2_new, e);
323
324         if (e_v1_new < e_v2_new) {
325                 e->verts[0] = e_v1_new;
326                 e->verts[1] = e_v2_new;
327         }
328         else {
329                 /* maintain winding info */
330                 e->verts[0] = e_v2_new;
331                 e->verts[1] = e_v1_new;
332
333                 SWAP(unsigned int, e->faces[0], e->faces[1]);
334                 SWAP(unsigned int, e->faces_other_v[0], e->faces_other_v[1]);
335         }
336
337         /* update adjacent data */
338         {
339                 unsigned int e_side = 0;
340
341                 for (e_side = 0; e_side < 2; e_side++) {
342                         /* 't_other' which we need to swap out is always the same edge-order */
343                         const unsigned int t_other = (((e->faces_other_v[e_side]) + 2)) % 3;
344                         unsigned int t_index = e->faces[e_side];
345                         unsigned int t_index_other = e->faces[!e_side];
346                         unsigned int *tri = tris[t_index];
347
348                         struct PolyEdge *e_other;
349                         unsigned int e_v1 = tri[(t_other    )    ];
350                         unsigned int e_v2 = tri[(t_other + 1) % 3];
351
352                         e_other = BLI_edgehash_lookup(ehash, e_v1, e_v2);
353                         if (e_other) {
354                                 BLI_assert(t_index != e_other->faces[0] && t_index != e_other->faces[1]);
355                                 if (t_index_other == e_other->faces[0]) {
356                                         e_other->faces[0] = t_index;
357                                         e_other->faces_other_v[0] = (t_other + 2) % 3;
358                                         BLI_assert(!ELEM(tri[e_other->faces_other_v[0]], e_v1, e_v2));
359                                 }
360                                 else if (t_index_other == e_other->faces[1]) {
361                                         e_other->faces[1] = t_index;
362                                         e_other->faces_other_v[1] = (t_other + 2) % 3;
363                                         BLI_assert(!ELEM(tri[e_other->faces_other_v[1]], e_v1, e_v2));
364                                 }
365                                 else {
366                                         BLI_assert(0);
367                                 }
368                         }
369                 }
370         }
371
372 #ifndef NDEBUG
373         polyfill_validate_tri(tris, e->faces[0], ehash);
374         polyfill_validate_tri(tris, e->faces[1], ehash);
375 #endif
376
377         BLI_assert(!ELEM(tris[e->faces[0]][e->faces_other_v[0]], UNPACK2(e->verts)));
378         BLI_assert(!ELEM(tris[e->faces[1]][e->faces_other_v[1]], UNPACK2(e->verts)));
379 }
380
381 /**
382  * The intention is that this calculates the output of #BLI_polyfill_calc
383  *
384  *
385  * \note assumes the \a coords form a boundary,
386  * so any edges running along contiguous (wrapped) indices,
387  * are ignored since the edges wont share 2 faces.
388  */
389 void BLI_polyfill_beautify(
390         const float (*coords)[2],
391         const unsigned int coords_tot,
392         unsigned int (*tris)[3],
393
394         /* structs for reuse */
395         MemArena *arena, Heap *eheap, EdgeHash *ehash)
396 {
397         const unsigned int coord_last = coords_tot - 1;
398         const unsigned int tris_tot = coords_tot - 2;
399         /* internal edges only (between 2 tris) */
400         const unsigned int edges_tot = tris_tot - 1;
401         unsigned int edges_tot_used = 0;
402         unsigned int i;
403
404         HeapNode **eheap_table;
405
406         struct PolyEdge *edges = BLI_memarena_alloc(arena, edges_tot * sizeof(*edges));
407
408         BLI_assert(BLI_heap_size(eheap) == 0);
409         BLI_assert(BLI_edgehash_size(ehash) == 0);
410
411         /* first build edges */
412         for (i = 0; i < tris_tot; i++) {
413                 unsigned int j_prev, j_curr, j_next;
414                 j_prev = 2;
415                 j_next = 1;
416                 for (j_curr = 0; j_curr < 3; j_next = j_prev, j_prev = j_curr++) {
417                         int e_index;
418
419                         unsigned int e_pair[2] = {
420                                 tris[i][j_prev],
421                                 tris[i][j_curr],
422                         };
423
424                         if (e_pair[0] > e_pair[1]) {
425                                 SWAP(unsigned int, e_pair[0], e_pair[1]);
426                                 e_index = 1;
427                         }
428                         else {
429                                 e_index = 0;
430                         }
431
432                         if (!is_boundary_edge(e_pair[0], e_pair[1], coord_last)) {
433                                 struct PolyEdge *e = BLI_edgehash_lookup(ehash, e_pair[0], e_pair[1]);
434                                 if (e == NULL) {
435                                         e = &edges[edges_tot_used++];
436                                         BLI_edgehash_insert(ehash, e_pair[0], e_pair[1], e);
437                                         memcpy(e->verts, e_pair, sizeof(e->verts));
438 #ifndef NDEBUG
439                                         e->faces[!e_index] = (unsigned int)-1;
440 #endif
441                                 }
442                                 else {
443
444                                         /* ensure each edge only ever has 2x users */
445 #ifndef NDEBUG
446                                         BLI_assert(e->faces[e_index] == (unsigned int)-1);
447                                         BLI_assert((e->verts[0] == e_pair[0]) &&
448                                                    (e->verts[1] == e_pair[1]));
449 #endif
450                                 }
451
452                                 e->faces[e_index] = i;
453                                 e->faces_other_v[e_index] = j_next;
454                         }
455                 }
456         }
457
458         /* now perform iterative rotations */
459         eheap_table = BLI_memarena_alloc(arena, sizeof(HeapNode *) * (size_t)edges_tot);
460
461         // for (i = 0; i < tris_tot; i++) { polyfill_validate_tri(tris, i, eh); }
462
463         /* build heap */
464         for (i = 0; i < edges_tot; i++) {
465                 struct PolyEdge *e = &edges[i];
466                 const float cost = polyedge_rotate_beauty_calc(coords, (const unsigned int (*)[3])tris, e);
467                 if (cost < 0.0f) {
468                         eheap_table[i] = BLI_heap_insert(eheap, cost, e);
469                 }
470                 else {
471                         eheap_table[i] = NULL;
472                 }
473         }
474
475         while (BLI_heap_is_empty(eheap) == false) {
476                 struct PolyEdge *e = BLI_heap_popmin(eheap);
477                 i = (unsigned int)(e - edges);
478                 eheap_table[i] = NULL;
479
480                 polyedge_rotate(tris, e, ehash);
481
482                 /* recalculate faces connected on the heap */
483                 polyedge_beauty_cost_update(
484                         coords, (const unsigned int (*)[3])tris, edges,
485                         e,
486                         eheap, eheap_table, ehash);
487         }
488
489         BLI_heap_clear(eheap, NULL);
490         BLI_edgehash_clear_ex(ehash, NULL, BLI_POLYFILL_ALLOC_NGON_RESERVE);
491
492         /* MEM_freeN(eheap_table); */  /* arena */
493 }