485 lines
15 KiB
C++
485 lines
15 KiB
C++
/* SPDX-FileCopyrightText: 2023 Blender Authors
|
|
*
|
|
* SPDX-License-Identifier: GPL-2.0-or-later */
|
|
|
|
/** \file
|
|
* \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.
|
|
* - Triangles 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.cc`).
|
|
*
|
|
* \note
|
|
*
|
|
* No globals - keep threadsafe.
|
|
*/
|
|
|
|
#include "BLI_utildefines.h"
|
|
|
|
#include "BLI_heap.h"
|
|
#include "BLI_math_geom.h"
|
|
#include "BLI_math_matrix.h"
|
|
#include "BLI_memarena.h"
|
|
|
|
#include "BLI_polyfill_2d_beautify.h" /* own include */
|
|
|
|
#include "BLI_strict_flags.h" /* IWYU pragma: keep. Keep last. */
|
|
|
|
/* 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 OrderEdge *x1 = static_cast<const OrderEdge *>(a1);
|
|
const OrderEdge *x2 = static_cast<const OrderEdge *>(a2);
|
|
if (x1->verts[0] > x2->verts[0]) {
|
|
return 1;
|
|
}
|
|
if (x1->verts[0] < x2->verts[0]) {
|
|
return -1;
|
|
}
|
|
|
|
if (x1->verts[1] > x2->verts[1]) {
|
|
return 1;
|
|
}
|
|
if (x1->verts[1] < x2->verts[1]) {
|
|
return -1;
|
|
}
|
|
|
|
/* Only for predictability. */
|
|
if (x1->e_half > x2->e_half) {
|
|
return 1;
|
|
}
|
|
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)));
|
|
}
|
|
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,
|
|
float *r_area)
|
|
{
|
|
/* 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));
|
|
|
|
if (r_area) {
|
|
*r_area = (fabsf(area_2x_234) + fabsf(area_2x_241) +
|
|
/* Include both pairs for predictable results. */
|
|
fabsf(area_2x_123) + fabsf(area_2x_134)) /
|
|
8.0f;
|
|
}
|
|
|
|
/*
|
|
* 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;
|
|
}
|
|
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;
|
|
}
|
|
|
|
return -FLT_MAX; /* always rotate */
|
|
}
|
|
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;
|
|
}
|
|
|
|
float BLI_polyfill_edge_calc_rotate_beauty__area(const float v1[3],
|
|
const float v2[3],
|
|
const float v3[3],
|
|
const float v4[3],
|
|
const bool lock_degenerate)
|
|
{
|
|
/* not a loop (only to be able to break out) */
|
|
do {
|
|
float v1_xy[2], v2_xy[2], v3_xy[2], v4_xy[2];
|
|
|
|
/* first get the 2d values */
|
|
{
|
|
const float eps = 1e-5f;
|
|
float no_a[3], no_b[3];
|
|
float no[3];
|
|
float axis_mat[3][3];
|
|
float no_scale;
|
|
cross_tri_v3(no_a, v2, v3, v4);
|
|
cross_tri_v3(no_b, v2, v4, v1);
|
|
|
|
// printf("%p %p %p %p - %p %p\n", v1, v2, v3, v4, e->l->f, e->l->radial_next->f);
|
|
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));
|
|
|
|
add_v3_v3v3(no, no_a, no_b);
|
|
if (UNLIKELY((no_scale = normalize_v3(no)) == 0.0f)) {
|
|
break;
|
|
}
|
|
|
|
axis_dominant_v3_to_m3(axis_mat, no);
|
|
mul_v2_m3v3(v1_xy, axis_mat, v1);
|
|
mul_v2_m3v3(v2_xy, axis_mat, v2);
|
|
mul_v2_m3v3(v3_xy, axis_mat, v3);
|
|
mul_v2_m3v3(v4_xy, axis_mat, v4);
|
|
|
|
/**
|
|
* Check if input faces are already flipped.
|
|
* Logic for 'signum_i' addition is:
|
|
*
|
|
* Accept:
|
|
* - (1, 1) or (-1, -1): same side (common case).
|
|
* - (-1/1, 0): one degenerate, OK since we may rotate into a valid state.
|
|
*
|
|
* Ignore:
|
|
* - (-1, 1): opposite winding, ignore.
|
|
* - ( 0, 0): both degenerate, ignore.
|
|
*
|
|
* \note The cross product is divided by 'no_scale'
|
|
* so the rotation calculation is scale independent.
|
|
*/
|
|
if (!(signum_i_ex(cross_tri_v2(v2_xy, v3_xy, v4_xy) / no_scale, eps) +
|
|
signum_i_ex(cross_tri_v2(v2_xy, v4_xy, v1_xy) / no_scale, eps)))
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
/**
|
|
* Important to lock degenerate here,
|
|
* since the triangle pars will be projected into different 2D spaces.
|
|
* Allowing to rotate out of a degenerate state can flip the faces
|
|
* (when performed iteratively).
|
|
*/
|
|
return BLI_polyfill_beautify_quad_rotate_calc_ex(
|
|
v1_xy, v2_xy, v3_xy, v4_xy, lock_degenerate, nullptr);
|
|
} while (false);
|
|
|
|
return FLT_MAX;
|
|
}
|
|
|
|
static float polyedge_rotate_beauty_calc(const float (*coords)[2],
|
|
const HalfEdge *edges,
|
|
const HalfEdge *e_a,
|
|
float *r_area)
|
|
{
|
|
const HalfEdge *e_b = &edges[e_a->e_radial];
|
|
|
|
const HalfEdge *e_a_other = &edges[edges[e_a->e_next].e_next];
|
|
const 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_ex(v1, v2, v3, v4, false, r_area);
|
|
}
|
|
|
|
static void polyedge_beauty_cost_update_single(const float (*coords)[2],
|
|
const HalfEdge *edges,
|
|
HalfEdge *e,
|
|
Heap *eheap,
|
|
HeapNode **eheap_table)
|
|
{
|
|
const uint i = e->base_index;
|
|
/* recalculate edge */
|
|
float area;
|
|
const float cost = polyedge_rotate_beauty_calc(coords, edges, e, &area);
|
|
/* 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 #43578, #49478.
|
|
*
|
|
* In fact a larger epsilon can still fail when the area of the face is very large,
|
|
* now the epsilon is scaled by the face area.
|
|
* See #56532. */
|
|
if (cost < -1e-6f * max_ff(area, 1.0f)) {
|
|
BLI_heap_insert_or_update(eheap, &eheap_table[i], cost, e);
|
|
}
|
|
else {
|
|
if (eheap_table[i]) {
|
|
BLI_heap_remove(eheap, eheap_table[i]);
|
|
eheap_table[i] = nullptr;
|
|
}
|
|
}
|
|
}
|
|
|
|
static void polyedge_beauty_cost_update(
|
|
const float (*coords)[2], HalfEdge *edges, HalfEdge *e, Heap *eheap, HeapNode **eheap_table)
|
|
{
|
|
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(HalfEdge *edges, const HalfEdge *e)
|
|
{
|
|
/** CCW winding, rotate internal edge to new vertical state.
|
|
*
|
|
* \code{.unparsed}
|
|
* Before After
|
|
* X X
|
|
* / \ /|\
|
|
* e4/ \e5 e4/ | \e5
|
|
* / e3 \ / | \
|
|
* X ------- X -> X e0|e3 X
|
|
* \ e0 / \ | /
|
|
* e2\ /e1 e2\ | /e1
|
|
* \ / \|/
|
|
* X X
|
|
* \endcode
|
|
*/
|
|
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;
|
|
}
|
|
|
|
void BLI_polyfill_beautify(const float (*coords)[2],
|
|
const uint coords_num,
|
|
uint (*tris)[3],
|
|
|
|
/* structs for reuse */
|
|
MemArena *arena,
|
|
Heap *eheap)
|
|
{
|
|
const uint coord_last = coords_num - 1;
|
|
const uint tris_len = coords_num - 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;
|
|
HalfEdge *half_edges = static_cast<HalfEdge *>(
|
|
BLI_memarena_alloc(arena, sizeof(*half_edges) * half_edges_len));
|
|
OrderEdge *order_edges = static_cast<OrderEdge *>(
|
|
BLI_memarena_alloc(arena, sizeof(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(OrderEdge), oedge_cmp);
|
|
|
|
for (uint i = 0, base_index = 0; i < order_edges_len; base_index++) {
|
|
const OrderEdge *oe_a = &order_edges[i++];
|
|
const OrderEdge *oe_b = &order_edges[i++];
|
|
BLI_assert(oe_a->verts[0] == oe_b->verts[0] && oe_a->verts[1] == oe_b->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 = (HeapNode **)order_edges;
|
|
order_edges = nullptr;
|
|
#endif
|
|
|
|
/* Build heap. */
|
|
{
|
|
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, nullptr);
|
|
if (cost < 0.0f) {
|
|
eheap_table[e->base_index] = BLI_heap_insert(eheap, cost, e);
|
|
}
|
|
else {
|
|
eheap_table[e->base_index] = nullptr;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
while (BLI_heap_is_empty(eheap) == false) {
|
|
HalfEdge *e = static_cast<HalfEdge *>(BLI_heap_pop_min(eheap));
|
|
eheap_table[e->base_index] = nullptr;
|
|
|
|
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, nullptr);
|
|
|
|
// MEM_freeN(eheap_table); /* arena */
|
|
|
|
/* Get triangles from half edge. */
|
|
uint tri_index = 0;
|
|
for (uint i = 0; i < half_edges_len; i++) {
|
|
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;
|
|
}
|
|
}
|
|
}
|