Cleanup: doxygen comments
[blender.git] / source / blender / blenlib / intern / polyfill_2d_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/polyfill_2d_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_heap.h"
46
47 #include "BLI_polyfill_2d_beautify.h"  /* own include */
48
49 #include "BLI_strict_flags.h"
50
51 /* Used to find matching edges. */
52 struct OrderEdge {
53         uint verts[2];
54         uint e_half;
55 };
56
57 /* Half edge used for rotating in-place. */
58 struct HalfEdge {
59         uint v;
60         uint e_next;
61         uint e_radial;
62         uint base_index;
63 };
64
65 static int oedge_cmp(const void *a1, const void *a2)
66 {
67         const struct OrderEdge *x1 = a1, *x2 = a2;
68         if (x1->verts[0] > x2->verts[0]) {
69                 return 1;
70         }
71         else if (x1->verts[0] < x2->verts[0]) {
72                 return -1;
73         }
74
75         if (x1->verts[1] > x2->verts[1]) {
76                 return 1;
77         }
78         else if (x1->verts[1] < x2->verts[1]) {
79                 return -1;
80         }
81
82         /* only for pradictability */
83         if (x1->e_half > x2->e_half) {
84                 return 1;
85         }
86         else if (x1->e_half < x2->e_half) {
87                 return -1;
88         }
89         /* Should never get here, no two edges should be the same. */
90         BLI_assert(false);
91         return 0;
92 }
93
94 BLI_INLINE bool is_boundary_edge(uint i_a, uint i_b, const uint coord_last)
95 {
96         BLI_assert(i_a < i_b);
97         return ((i_a + 1 == i_b) || UNLIKELY((i_a == 0) && (i_b == coord_last)));
98 }
99 /**
100  * Assuming we have 2 triangles sharing an edge (2 - 4),
101  * check if the edge running from (1 - 3) gives better results.
102  *
103  * \param lock_degenerate: Use to avoid rotating out of a degenerate state.
104  * - When true, an existing zero area face on either side of the (2 - 4) split will return a positive value.
105  * - When false, the check must be non-biased towards either split direction.
106  *
107  * \return (negative number means the edge can be rotated, lager == better).
108  */
109 float BLI_polyfill_beautify_quad_rotate_calc_ex(
110         const float v1[2], const float v2[2], const float v3[2], const float v4[2],
111         const bool lock_degenerate)
112 {
113         /* not a loop (only to be able to break out) */
114         do {
115                 /* Allow very small faces to be considered non-zero. */
116                 const float eps_zero_area = 1e-12f;
117                 const float area_2x_234 = cross_tri_v2(v2, v3, v4);
118                 const float area_2x_241 = cross_tri_v2(v2, v4, v1);
119
120                 const float area_2x_123 = cross_tri_v2(v1, v2, v3);
121                 const float area_2x_134 = cross_tri_v2(v1, v3, v4);
122
123                 BLI_assert((ELEM(v1, v2, v3, v4) == false) &&
124                            (ELEM(v2, v1, v3, v4) == false) &&
125                            (ELEM(v3, v1, v2, v4) == false) &&
126                            (ELEM(v4, v1, v2, v3) == false));
127                 /*
128                  * Test for unusable (1-3) state.
129                  * - Area sign flipping to check faces aren't going to point in opposite directions.
130                  * - Area epsilon check that the one of the faces won't be zero area.
131                  */
132                 if ((area_2x_123 >= 0.0f) != (area_2x_134 >= 0.0f)) {
133                         break;
134                 }
135                 else if ((fabsf(area_2x_123) <= eps_zero_area) || (fabsf(area_2x_134) <= eps_zero_area)) {
136                         break;
137                 }
138
139                 /* Test for unusable (2-4) state (same as above). */
140                 if ((area_2x_234 >= 0.0f) != (area_2x_241 >= 0.0f)) {
141                         if (lock_degenerate) {
142                                 break;
143                         }
144                         else {
145                                 return -FLT_MAX;  /* always rotate */
146                         }
147                 }
148                 else if ((fabsf(area_2x_234) <= eps_zero_area) || (fabsf(area_2x_241) <= eps_zero_area)) {
149                         return -FLT_MAX;  /* always rotate */
150                 }
151
152                 {
153                         /* testing rule: the area divided by the perimeter,
154                          * check if (1-3) beats the existing (2-4) edge rotation */
155                         float area_a, area_b;
156                         float prim_a, prim_b;
157                         float fac_24, fac_13;
158
159                         float len_12, len_23, len_34, len_41, len_24, len_13;
160
161                         /* edges around the quad */
162                         len_12 = len_v2v2(v1, v2);
163                         len_23 = len_v2v2(v2, v3);
164                         len_34 = len_v2v2(v3, v4);
165                         len_41 = len_v2v2(v4, v1);
166                         /* edges crossing the quad interior */
167                         len_13 = len_v2v2(v1, v3);
168                         len_24 = len_v2v2(v2, v4);
169
170                         /* note, area is in fact (area * 2),
171                          * but in this case its OK, since we're comparing ratios */
172
173                         /* edge (2-4), current state */
174                         area_a = fabsf(area_2x_234);
175                         area_b = fabsf(area_2x_241);
176                         prim_a = len_23 + len_34 + len_24;
177                         prim_b = len_41 + len_12 + len_24;
178                         fac_24 = (area_a / prim_a) + (area_b / prim_b);
179
180                         /* edge (1-3), new state */
181                         area_a = fabsf(area_2x_123);
182                         area_b = fabsf(area_2x_134);
183                         prim_a = len_12 + len_23 + len_13;
184                         prim_b = len_34 + len_41 + len_13;
185                         fac_13 = (area_a / prim_a) + (area_b / prim_b);
186
187                         /* negative number if (1-3) is an improved state */
188                         return fac_24 - fac_13;
189                 }
190         } while (false);
191
192         return FLT_MAX;
193 }
194
195 static float polyedge_rotate_beauty_calc(
196         const float (*coords)[2],
197         const struct HalfEdge *edges,
198         const struct HalfEdge *e_a)
199 {
200         const struct HalfEdge *e_b = &edges[e_a->e_radial];
201
202         const struct HalfEdge *e_a_other = &edges[edges[e_a->e_next].e_next];
203         const struct HalfEdge *e_b_other = &edges[edges[e_b->e_next].e_next];
204
205         const float *v1, *v2, *v3, *v4;
206
207         v1 = coords[e_a_other->v];
208         v2 = coords[e_a->v];
209         v3 = coords[e_b_other->v];
210         v4 = coords[e_b->v];
211
212         return BLI_polyfill_beautify_quad_rotate_calc(v1, v2, v3, v4);
213 }
214
215 static void polyedge_beauty_cost_update_single(
216         const float (*coords)[2],
217         const struct HalfEdge *edges,
218         struct HalfEdge *e,
219         Heap *eheap, HeapNode **eheap_table)
220 {
221         const uint i = e->base_index;
222         /* recalculate edge */
223         const float cost = polyedge_rotate_beauty_calc(coords, edges, e);
224         /* We can get cases where both choices generate very small negative costs, which leads to infinite loop.
225          * Anyway, costs above that are not worth recomputing, maybe we could even optimize it to a smaller limit?
226          * Actually, FLT_EPSILON is too small in some cases, 1e-6f seems to work OK hopefully?
227          * See T43578, T49478. */
228         if (cost < -1e-6f) {
229                 BLI_heap_insert_or_update(eheap, &eheap_table[i], cost, e);
230         }
231         else {
232                 if (eheap_table[i]) {
233                         BLI_heap_remove(eheap, eheap_table[i]);
234                         eheap_table[i] = NULL;
235                 }
236         }
237 }
238
239 static void polyedge_beauty_cost_update(
240         const float (*coords)[2],
241         struct HalfEdge *edges,
242         struct HalfEdge *e,
243         Heap *eheap, HeapNode **eheap_table)
244 {
245         struct HalfEdge *e_arr[4];
246         e_arr[0] = &edges[e->e_next];
247         e_arr[1] = &edges[e_arr[0]->e_next];
248
249         e = &edges[e->e_radial];
250         e_arr[2] = &edges[e->e_next];
251         e_arr[3] = &edges[e_arr[2]->e_next];
252
253         for (uint i = 0; i < 4; i++) {
254                 if (e_arr[i] && e_arr[i]->base_index != UINT_MAX) {
255                         polyedge_beauty_cost_update_single(
256                                 coords, edges,
257                                 e_arr[i],
258                                 eheap, eheap_table);
259                 }
260         }
261 }
262
263 static void polyedge_rotate(
264         struct HalfEdge *edges,
265         struct HalfEdge *e)
266 {
267         /** CCW winding, rotate internal edge to new vertical state.
268          *
269          * \code{.unparsed}
270          *   Before         After
271          *      X             X
272          *     / \           /|\
273          *  e4/   \e5     e4/ | \e5
274          *   / e3  \       /  |  \
275          * X ------- X -> X e0|e3 X
276          *   \ e0  /       \  |  /
277          *  e2\   /e1     e2\ | /e1
278          *     \ /           \|/
279          *      X             X
280          * \endcode
281          */
282         struct HalfEdge *ed[6];
283         uint ed_index[6];
284
285         ed_index[0] = (uint)(e - edges);
286         ed[0] = &edges[ed_index[0]];
287         ed_index[1] = ed[0]->e_next;
288         ed[1] = &edges[ed_index[1]];
289         ed_index[2] = ed[1]->e_next;
290         ed[2] = &edges[ed_index[2]];
291
292         ed_index[3] = e->e_radial;
293         ed[3] = &edges[ed_index[3]];
294         ed_index[4] = ed[3]->e_next;
295         ed[4] = &edges[ed_index[4]];
296         ed_index[5] = ed[4]->e_next;
297         ed[5] = &edges[ed_index[5]];
298
299         ed[0]->e_next = ed_index[2];
300         ed[1]->e_next = ed_index[3];
301         ed[2]->e_next = ed_index[4];
302         ed[3]->e_next = ed_index[5];
303         ed[4]->e_next = ed_index[0];
304         ed[5]->e_next = ed_index[1];
305
306         ed[0]->v = ed[5]->v;
307         ed[3]->v = ed[2]->v;
308 }
309
310 /**
311  * The intention is that this calculates the output of #BLI_polyfill_calc
312  *
313  *
314  * \note assumes the \a coords form a boundary,
315  * so any edges running along contiguous (wrapped) indices,
316  * are ignored since the edges wont share 2 faces.
317  */
318 void BLI_polyfill_beautify(
319         const float (*coords)[2],
320         const uint coords_tot,
321         uint (*tris)[3],
322
323         /* structs for reuse */
324         MemArena *arena, Heap *eheap)
325 {
326         const uint coord_last = coords_tot - 1;
327         const uint tris_len = coords_tot - 2;
328         /* internal edges only (between 2 tris) */
329         const uint edges_len = tris_len - 1;
330
331         HeapNode **eheap_table;
332
333         const uint half_edges_len = 3 * tris_len;
334         struct HalfEdge *half_edges = BLI_memarena_alloc(arena, sizeof(*half_edges) * half_edges_len);
335         struct OrderEdge *order_edges = BLI_memarena_alloc(arena, sizeof(struct OrderEdge) * 2 * edges_len);
336         uint order_edges_len = 0;
337
338         /* first build edges */
339         for (uint i = 0; i < tris_len; i++) {
340                 for (uint j_curr = 0, j_prev = 2; j_curr < 3; j_prev = j_curr++) {
341                         const uint e_index_prev = (i * 3) + j_prev;
342                         const uint e_index_curr = (i * 3) + j_curr;
343
344                         half_edges[e_index_prev].v = tris[i][j_prev];
345                         half_edges[e_index_prev].e_next = e_index_curr;
346                         half_edges[e_index_prev].e_radial = UINT_MAX;
347                         half_edges[e_index_prev].base_index = UINT_MAX;
348
349                         uint e_pair[2] = {tris[i][j_prev], tris[i][j_curr]};
350                         if (e_pair[0] > e_pair[1]) {
351                                 SWAP(uint, e_pair[0], e_pair[1]);
352                         }
353
354                         /* ensure internal edges. */
355                         if (!is_boundary_edge(e_pair[0], e_pair[1], coord_last)) {
356                                 order_edges[order_edges_len].verts[0] = e_pair[0];
357                                 order_edges[order_edges_len].verts[1] = e_pair[1];
358                                 order_edges[order_edges_len].e_half = e_index_prev;
359                                 order_edges_len += 1;
360                         }
361                 }
362         }
363         BLI_assert(edges_len * 2 == order_edges_len);
364
365         qsort(order_edges, order_edges_len, sizeof(struct OrderEdge), oedge_cmp);
366
367         for (uint i = 0, base_index = 0; i < order_edges_len; base_index++) {
368                 const struct OrderEdge *oe_a = &order_edges[i++];
369                 const struct OrderEdge *oe_b = &order_edges[i++];
370                 BLI_assert(oe_a->verts[0] == oe_a->verts[0] && oe_a->verts[1] == oe_a->verts[1]);
371                 half_edges[oe_a->e_half].e_radial = oe_b->e_half;
372                 half_edges[oe_b->e_half].e_radial = oe_a->e_half;
373                 half_edges[oe_a->e_half].base_index = base_index;
374                 half_edges[oe_b->e_half].base_index = base_index;
375         }
376         /* order_edges could be freed now. */
377
378         /* Now perform iterative rotations. */
379 #if 0
380         eheap_table = BLI_memarena_alloc(arena, sizeof(HeapNode *) * (size_t)edges_len);
381 #else
382         /* We can re-use this since its big enough. */
383         eheap_table = (void *)order_edges;
384         order_edges = NULL;
385 #endif
386
387         /* Build heap. */
388         {
389                 struct HalfEdge *e = half_edges;
390                 for (uint i = 0; i < half_edges_len; i++, e++) {
391                         /* Accounts for boundary edged too (UINT_MAX). */
392                         if (e->e_radial < i) {
393                                 const float cost = polyedge_rotate_beauty_calc(coords, half_edges, e);
394                                 if (cost < 0.0f) {
395                                         eheap_table[e->base_index] = BLI_heap_insert(eheap, cost, e);
396                                 }
397                                 else {
398                                         eheap_table[e->base_index] = NULL;
399                                 }
400                         }
401                 }
402         }
403
404         while (BLI_heap_is_empty(eheap) == false) {
405                 struct HalfEdge *e = BLI_heap_pop_min(eheap);
406                 eheap_table[e->base_index] = NULL;
407
408                 polyedge_rotate(half_edges, e);
409
410                 /* recalculate faces connected on the heap */
411                 polyedge_beauty_cost_update(
412                         coords, half_edges,
413                         e,
414                         eheap, eheap_table);
415         }
416
417         BLI_heap_clear(eheap, NULL);
418
419         /* MEM_freeN(eheap_table); */  /* arena */
420
421         /* get tris from half edge. */
422         uint tri_index = 0;
423         for (uint i = 0; i < half_edges_len; i++) {
424                 struct HalfEdge *e = &half_edges[i];
425                 if (e->v != UINT_MAX) {
426                         uint *tri = tris[tri_index++];
427
428                         tri[0] = e->v;
429                         e->v = UINT_MAX;
430
431                         e = &half_edges[e->e_next];
432                         tri[1] = e->v;
433                         e->v = UINT_MAX;
434
435                         e = &half_edges[e->e_next];
436                         tri[2] = e->v;
437                         e->v = UINT_MAX;
438                 }
439         }
440 }