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