/* * ***** BEGIN GPL LICENSE BLOCK ***** * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software Foundation, * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * * ***** END GPL LICENSE BLOCK ***** */ /** \file blender/blenlib/intern/polyfill2d_beautify.c * \ingroup bli * * This function is to improve the tessellation resulting from polyfill2d, * creating optimal topology. * * The functionality here matches #BM_mesh_beautify_fill, * but its far simpler to perform this operation in 2d, * on a simple polygon representation where we _know_: * * - The polygon is primitive with no holes with a continuous boundary. * - Tris have consistent winding. * - 2d (saves some hassles projecting face pairs on an axis for every edge-rotation) * also saves us having to store all previous edge-states (see #EdRotState in bmesh_beautify.c) * * \note * * No globals - keep threadsafe. */ #include "BLI_utildefines.h" #include "BLI_math.h" #include "BLI_memarena.h" #include "BLI_heap.h" #include "BLI_polyfill2d_beautify.h" /* own include */ #include "BLI_strict_flags.h" /* Used to find matching edges. */ struct OrderEdge { uint verts[2]; uint e_half; }; /* Half edge used for rotating in-place. */ struct HalfEdge { uint v; uint e_next; uint e_radial; uint base_index; }; static int oedge_cmp(const void *a1, const void *a2) { const struct OrderEdge *x1 = a1, *x2 = a2; if (x1->verts[0] > x2->verts[0]) { return 1; } else if (x1->verts[0] < x2->verts[0]) { return -1; } if (x1->verts[1] > x2->verts[1]) { return 1; } else if (x1->verts[1] < x2->verts[1]) { return -1; } /* only for pradictability */ if (x1->e_half > x2->e_half) { return 1; } else if (x1->e_half < x2->e_half) { return -1; } /* Should never get here, no two edges should be the same. */ BLI_assert(false); return 0; } BLI_INLINE bool is_boundary_edge(uint i_a, uint i_b, const uint coord_last) { BLI_assert(i_a < i_b); return ((i_a + 1 == i_b) || UNLIKELY((i_a == 0) && (i_b == coord_last))); } /** * Assuming we have 2 triangles sharing an edge (2 - 4), * check if the edge running from (1 - 3) gives better results. * * \param lock_degenerate: Use to avoid rotating out of a degenerate state. * - When true, an existing zero area face on either side of the (2 - 4) split will return a positive value. * - When false, the check must be non-biased towards either split direction. * * \return (negative number means the edge can be rotated, lager == better). */ float BLI_polyfill_beautify_quad_rotate_calc_ex( const float v1[2], const float v2[2], const float v3[2], const float v4[2], const bool lock_degenerate) { /* not a loop (only to be able to break out) */ do { /* Allow very small faces to be considered non-zero. */ const float eps_zero_area = 1e-12f; const float area_2x_234 = cross_tri_v2(v2, v3, v4); const float area_2x_241 = cross_tri_v2(v2, v4, v1); const float area_2x_123 = cross_tri_v2(v1, v2, v3); const float area_2x_134 = cross_tri_v2(v1, v3, v4); BLI_assert((ELEM(v1, v2, v3, v4) == false) && (ELEM(v2, v1, v3, v4) == false) && (ELEM(v3, v1, v2, v4) == false) && (ELEM(v4, v1, v2, v3) == false)); /* * Test for unusable (1-3) state. * - Area sign flipping to check faces aren't going to point in opposite directions. * - Area epsilon check that the one of the faces won't be zero area. */ if ((area_2x_123 >= 0.0f) != (area_2x_134 >= 0.0f)) { break; } else if ((fabsf(area_2x_123) <= eps_zero_area) || (fabsf(area_2x_134) <= eps_zero_area)) { break; } /* Test for unusable (2-4) state (same as above). */ if ((area_2x_234 >= 0.0f) != (area_2x_241 >= 0.0f)) { if (lock_degenerate) { break; } else { return -FLT_MAX; /* always rotate */ } } else if ((fabsf(area_2x_234) <= eps_zero_area) || (fabsf(area_2x_241) <= eps_zero_area)) { return -FLT_MAX; /* always rotate */ } { /* testing rule: the area divided by the perimeter, * check if (1-3) beats the existing (2-4) edge rotation */ float area_a, area_b; float prim_a, prim_b; float fac_24, fac_13; float len_12, len_23, len_34, len_41, len_24, len_13; /* edges around the quad */ len_12 = len_v2v2(v1, v2); len_23 = len_v2v2(v2, v3); len_34 = len_v2v2(v3, v4); len_41 = len_v2v2(v4, v1); /* edges crossing the quad interior */ len_13 = len_v2v2(v1, v3); len_24 = len_v2v2(v2, v4); /* note, area is in fact (area * 2), * but in this case its OK, since we're comparing ratios */ /* edge (2-4), current state */ area_a = fabsf(area_2x_234); area_b = fabsf(area_2x_241); prim_a = len_23 + len_34 + len_24; prim_b = len_41 + len_12 + len_24; fac_24 = (area_a / prim_a) + (area_b / prim_b); /* edge (1-3), new state */ area_a = fabsf(area_2x_123); area_b = fabsf(area_2x_134); prim_a = len_12 + len_23 + len_13; prim_b = len_34 + len_41 + len_13; fac_13 = (area_a / prim_a) + (area_b / prim_b); /* negative number if (1-3) is an improved state */ return fac_24 - fac_13; } } while (false); return FLT_MAX; } static float polyedge_rotate_beauty_calc( const float (*coords)[2], const struct HalfEdge *edges, const struct HalfEdge *e_a) { const struct HalfEdge *e_b = &edges[e_a->e_radial]; const struct HalfEdge *e_a_other = &edges[edges[e_a->e_next].e_next]; const struct HalfEdge *e_b_other = &edges[edges[e_b->e_next].e_next]; const float *v1, *v2, *v3, *v4; v1 = coords[e_a_other->v]; v2 = coords[e_a->v]; v3 = coords[e_b_other->v]; v4 = coords[e_b->v]; return BLI_polyfill_beautify_quad_rotate_calc(v1, v2, v3, v4); } static void polyedge_beauty_cost_update_single( const float (*coords)[2], const struct HalfEdge *edges, struct HalfEdge *e, Heap *eheap, HeapNode **eheap_table) { const uint i = e->base_index; if (eheap_table[i]) { BLI_heap_remove(eheap, eheap_table[i]); eheap_table[i] = NULL; } { /* recalculate edge */ const float cost = polyedge_rotate_beauty_calc(coords, edges, e); /* We can get cases where both choices generate very small negative costs, which leads to infinite loop. * Anyway, costs above that are not worth recomputing, maybe we could even optimize it to a smaller limit? * Actually, FLT_EPSILON is too small in some cases, 1e-6f seems to work OK hopefully? * See T43578, T49478. */ if (cost < -1e-6f) { eheap_table[i] = BLI_heap_insert(eheap, cost, e); } else { eheap_table[i] = NULL; } } } static void polyedge_beauty_cost_update( const float (*coords)[2], struct HalfEdge *edges, struct HalfEdge *e, Heap *eheap, HeapNode **eheap_table) { struct HalfEdge *e_arr[4]; e_arr[0] = &edges[e->e_next]; e_arr[1] = &edges[e_arr[0]->e_next]; e = &edges[e->e_radial]; e_arr[2] = &edges[e->e_next]; e_arr[3] = &edges[e_arr[2]->e_next]; for (uint i = 0; i < 4; i++) { if (e_arr[i] && e_arr[i]->base_index != UINT_MAX) { polyedge_beauty_cost_update_single( coords, edges, e_arr[i], eheap, eheap_table); } } } static void polyedge_rotate( struct HalfEdge *edges, struct HalfEdge *e) { /** CCW winding, rotate internal edge to new vertical state. *
* Before After * X X * / \ /|\ * e4/ \e5 e4/ | \e5 * / e3 \ / | \ * X ------- X -> X e0|e3 X * \ e0 / \ | / * e2\ /e1 e2\ | /e1 * \ / \|/ * X X **/ struct HalfEdge *ed[6]; uint ed_index[6]; ed_index[0] = (uint)(e - edges); ed[0] = &edges[ed_index[0]]; ed_index[1] = ed[0]->e_next; ed[1] = &edges[ed_index[1]]; ed_index[2] = ed[1]->e_next; ed[2] = &edges[ed_index[2]]; ed_index[3] = e->e_radial; ed[3] = &edges[ed_index[3]]; ed_index[4] = ed[3]->e_next; ed[4] = &edges[ed_index[4]]; ed_index[5] = ed[4]->e_next; ed[5] = &edges[ed_index[5]]; ed[0]->e_next = ed_index[2]; ed[1]->e_next = ed_index[3]; ed[2]->e_next = ed_index[4]; ed[3]->e_next = ed_index[5]; ed[4]->e_next = ed_index[0]; ed[5]->e_next = ed_index[1]; ed[0]->v = ed[5]->v; ed[3]->v = ed[2]->v; } /** * The intention is that this calculates the output of #BLI_polyfill_calc * * * \note assumes the \a coords form a boundary, * so any edges running along contiguous (wrapped) indices, * are ignored since the edges wont share 2 faces. */ void BLI_polyfill_beautify( const float (*coords)[2], const uint coords_tot, uint (*tris)[3], /* structs for reuse */ MemArena *arena, Heap *eheap) { const uint coord_last = coords_tot - 1; const uint tris_len = coords_tot - 2; /* internal edges only (between 2 tris) */ const uint edges_len = tris_len - 1; HeapNode **eheap_table; const uint half_edges_len = 3 * tris_len; struct HalfEdge *half_edges = BLI_memarena_alloc(arena, sizeof(*half_edges) * half_edges_len); struct OrderEdge *order_edges = BLI_memarena_alloc(arena, sizeof(struct OrderEdge) * 2 * edges_len); uint order_edges_len = 0; /* first build edges */ for (uint i = 0; i < tris_len; i++) { for (uint j_curr = 0, j_prev = 2; j_curr < 3; j_prev = j_curr++) { const uint e_index_prev = (i * 3) + j_prev; const uint e_index_curr = (i * 3) + j_curr; half_edges[e_index_prev].v = tris[i][j_prev]; half_edges[e_index_prev].e_next = e_index_curr; half_edges[e_index_prev].e_radial = UINT_MAX; half_edges[e_index_prev].base_index = UINT_MAX; uint e_pair[2] = {tris[i][j_prev], tris[i][j_curr]}; if (e_pair[0] > e_pair[1]) { SWAP(uint, e_pair[0], e_pair[1]); } /* ensure internal edges. */ if (!is_boundary_edge(e_pair[0], e_pair[1], coord_last)) { order_edges[order_edges_len].verts[0] = e_pair[0]; order_edges[order_edges_len].verts[1] = e_pair[1]; order_edges[order_edges_len].e_half = e_index_prev; order_edges_len += 1; } } } BLI_assert(edges_len * 2 == order_edges_len); qsort(order_edges, order_edges_len, sizeof(struct OrderEdge), oedge_cmp); for (uint i = 0, base_index = 0; i < order_edges_len; base_index++) { const struct OrderEdge *oe_a = &order_edges[i++]; const struct OrderEdge *oe_b = &order_edges[i++]; BLI_assert(oe_a->verts[0] == oe_a->verts[0] && oe_a->verts[1] == oe_a->verts[1]); half_edges[oe_a->e_half].e_radial = oe_b->e_half; half_edges[oe_b->e_half].e_radial = oe_a->e_half; half_edges[oe_a->e_half].base_index = base_index; half_edges[oe_b->e_half].base_index = base_index; } /* order_edges could be freed now. */ /* Now perform iterative rotations. */ #if 0 eheap_table = BLI_memarena_alloc(arena, sizeof(HeapNode *) * (size_t)edges_len); #else /* We can re-use this since its big enough. */ eheap_table = (void *)order_edges; order_edges = NULL; #endif /* Build heap. */ { struct HalfEdge *e = half_edges; for (uint i = 0; i < half_edges_len; i++, e++) { /* Accounts for boundary edged too (UINT_MAX). */ if (e->e_radial < i) { const float cost = polyedge_rotate_beauty_calc(coords, half_edges, e); if (cost < 0.0f) { eheap_table[e->base_index] = BLI_heap_insert(eheap, cost, e); } else { eheap_table[e->base_index] = NULL; } } } } while (BLI_heap_is_empty(eheap) == false) { struct HalfEdge *e = BLI_heap_popmin(eheap); eheap_table[e->base_index] = NULL; polyedge_rotate(half_edges, e); /* recalculate faces connected on the heap */ polyedge_beauty_cost_update( coords, half_edges, e, eheap, eheap_table); } BLI_heap_clear(eheap, NULL); /* MEM_freeN(eheap_table); */ /* arena */ /* get tris from half edge. */ uint tri_index = 0; for (uint i = 0; i < half_edges_len; i++) { struct HalfEdge *e = &half_edges[i]; if (e->v != UINT_MAX) { uint *tri = tris[tri_index++]; tri[0] = e->v; e->v = UINT_MAX; e = &half_edges[e->e_next]; tri[1] = e->v; e->v = UINT_MAX; e = &half_edges[e->e_next]; tri[2] = e->v; e->v = UINT_MAX; } } }