From 2355bb998d5bb39dcebe1d83daf70e11c4c8e158 Mon Sep 17 00:00:00 2001 From: ExPikaPaka Date: Tue, 15 Sep 2026 08:58:26 +0200 Subject: [PATCH] Texture displacement: step cutter, v2 colours, auto resolution, anchored bake frame --- src/libslic3r/TextureDisplacement.cpp | 2422 +++++++++++++++++++++++-- src/libslic3r/TextureDisplacement.hpp | 186 +- 2 files changed, 2398 insertions(+), 210 deletions(-) diff --git a/src/libslic3r/TextureDisplacement.cpp b/src/libslic3r/TextureDisplacement.cpp index a129b503ff..cfc75c286f 100644 --- a/src/libslic3r/TextureDisplacement.cpp +++ b/src/libslic3r/TextureDisplacement.cpp @@ -1,7 +1,11 @@ #include "TextureDisplacement.hpp" #include +#include +#include +#include #include +#include #include #include #include @@ -10,17 +14,25 @@ #include #include #include +#include +#include #include #include +#include + #include #include #include +#include + +#include "AABBTreeIndirect.hpp" #include "MeshBoolean.hpp" #include "Model.hpp" #include "PNGReadWrite.hpp" #include "TriangleSelector.hpp" +#include "TextureBake/TextureBakeDebug.hpp" #include "TextureBake/TextureBakeMesh.hpp" #include "TextureBake/TextureBakePipeline.hpp" @@ -120,48 +132,106 @@ struct DecodedTextureCache std::unordered_map>, DecodedHeightTexture>> entries; }; DecodedTextureCache g_decoded_texture_cache; + +// The smoothed copy is cached too, one per image: the thumbnail, the bump-preview height texture, the +// colour texture, the projector texture, the UV editor's background and the preview job all ask for +// the same (image, smoothing) pair in the same frame while the Smoothing slider moves, and each of them +// blurring its own copy is what froze the UI. Keyed like the raw cache; a different smoothing value +// simply replaces the entry. +struct SmoothedTextureCache +{ + struct Entry + { + std::weak_ptr> source; + float smoothing = 0.f; + DecodedHeightTexture texture; + }; + std::mutex mutex; + std::unordered_map entries; +}; +SmoothedTextureCache g_smoothed_texture_cache; } // namespace namespace { -// A few passes of a separable box blur approximate a Gaussian, cheaply. `radius` is in pixels; 0 is -// a no-op. Wraps at the edges so a tiling height map stays seamless after smoothing. Operates on the -// grayscale byte buffer in place. -void smooth_height_pixels(std::vector &pixels, int width, int height, int radius) +// A few passes of a separable box blur approximate a Gaussian, cheaply. `radius` is in whole texels; +// 0 is a no-op. Wraps at the edges so a tiling height map stays seamless after smoothing. Operates on +// the grayscale byte buffer in place. +// +// Each pass is a sliding window - one add and one subtract per texel - so the cost is the image size, +// not the image size times the radius. The Smoothing slider drives this on every frame it moves, for +// every consumer of the layer, and the radius goes up to 48 texels: the per-window loop this replaces +// stalled the UI for seconds on a large map. Rows (and column blocks) run in parallel. The rounding is +// the old code's exactly, so a blur gives the same bytes as before. +void smooth_height_pixels_box(std::vector &pixels, int width, int height, int radius) { if (radius <= 0 || width <= 0 || height <= 0 || pixels.size() != size_t(width) * size_t(height)) return; const int window = 2 * radius + 1; - const float inv = 1.f / float(window); + const float inv = 1.f / float(window); + const auto wrap = [](int i, int n) { return (i % n + n) % n; }; std::vector tmp(pixels.size()); for (int pass = 0; pass < 2; ++pass) { // two passes -> smoother than a single box - // Horizontal. - for (int y = 0; y < height; ++y) { - const size_t row = size_t(y) * width; - for (int x = 0; x < width; ++x) { - float sum = 0.f; - for (int k = -radius; k <= radius; ++k) { - int sx = x + k; - sx = (sx % width + width) % width; // wrap - sum += float(pixels[row + size_t(sx)]); + // Horizontal: a running sum per row. + tbb::parallel_for(tbb::blocked_range(0, height), [&](const tbb::blocked_range &r) { + for (int y = r.begin(); y < r.end(); ++y) { + const uint8_t *src = pixels.data() + size_t(y) * size_t(width); + uint8_t *dst = tmp.data() + size_t(y) * size_t(width); + uint32_t sum = 0; + for (int k = -radius; k <= radius; ++k) + sum += src[wrap(k, width)]; + for (int x = 0; x < width; ++x) { + dst[x] = uint8_t(std::lround(float(sum) * inv)); + sum += src[wrap(x + radius + 1, width)]; + sum -= src[wrap(x - radius, width)]; } - tmp[row + size_t(x)] = uint8_t(std::lround(sum * inv)); } - } - // Vertical. - for (int x = 0; x < width; ++x) + }); + // Vertical: one running sum per column, advanced row by row so the reads stay row-major. + tbb::parallel_for(tbb::blocked_range(0, width, 256), [&](const tbb::blocked_range &r) { + std::vector sum(size_t(r.size()), 0); + for (int k = -radius; k <= radius; ++k) { + const uint8_t *row = tmp.data() + size_t(wrap(k, height)) * size_t(width); + for (int x = r.begin(); x < r.end(); ++x) + sum[size_t(x - r.begin())] += row[x]; + } for (int y = 0; y < height; ++y) { - float sum = 0.f; - for (int k = -radius; k <= radius; ++k) { - int sy = y + k; - sy = (sy % height + height) % height; // wrap - sum += float(tmp[size_t(sy) * width + size_t(x)]); + uint8_t *dst = pixels.data() + size_t(y) * size_t(width); + const uint8_t *add = tmp.data() + size_t(wrap(y + radius + 1, height)) * size_t(width); + const uint8_t *sub = tmp.data() + size_t(wrap(y - radius, height)) * size_t(width); + for (int x = r.begin(); x < r.end(); ++x) { + uint32_t &sx = sum[size_t(x - r.begin())]; + dst[x] = uint8_t(std::lround(float(sx) * inv)); + sx += add[x]; + sx -= sub[x]; } - pixels[size_t(y) * width + size_t(x)] = uint8_t(std::lround(sum * inv)); } + }); } } + +// The same blur with a *continuous* radius, which is what the Smoothing slider drives. +// +// A box blur can only work in whole texels, so mapping the slider straight onto a rounded radius +// made it move in visible jumps - and its very first step off zero was a full one-texel blur rather +// than a hint of one, which is what made the control feel like it switched on rather than ramped up. +// Blur at the next whole texel up and cross-fade the raw image back in by the fraction left over: +// below one texel that fade *is* the sub-texel kernel, and above it it turns each integer step into +// a continuous ramp. +void smooth_height_pixels(std::vector &pixels, int width, int height, float radius) +{ + if (radius <= 0.f || width <= 0 || height <= 0 || pixels.size() != size_t(width) * size_t(height)) + return; + + const int whole = std::max(1, int(std::ceil(radius))); + const float mix = std::clamp(radius / float(whole), 0.f, 1.f); + const std::vector raw = (mix < 0.999f) ? pixels : std::vector{}; + smooth_height_pixels_box(pixels, width, height, whole); + if (!raw.empty()) + for (size_t i = 0; i < pixels.size(); ++i) + pixels[i] = uint8_t(std::lround(float(raw[i]) + (float(pixels[i]) - float(raw[i])) * mix)); +} } // namespace DecodedHeightTexture decode_height_texture(const TextureDisplacementLayer &layer) @@ -231,11 +301,21 @@ DecodedHeightTexture decode_height_texture(const TextureDisplacementLayer &layer entries[key] = { std::weak_ptr>(layer.image_data), result }; } - // Smoothing radius scales with the texture so the effect is resolution-independent; capped so the - // blur stays affordable on large maps. + // Smoothing radius scales with the texture so the same slider value blurs the same *fraction* of + // the image whatever resolution it came in at, and stays continuous in the slider - see + // smooth_height_pixels(). The cap is a cost limit, not part of the mapping: the blur is + // O(width * height * radius) per pass, so a large map with the slider at the top would otherwise + // stall every preview rebuild. if (layer.smoothing > 0.f) { - const int max_radius = std::clamp(int(std::lround(0.02f * std::min(result.width, result.height))), 1, 32); - const int radius = std::max(1, int(std::lround(layer.smoothing * float(max_radius)))); + { + std::lock_guard lock(g_smoothed_texture_cache.mutex); + auto it = g_smoothed_texture_cache.entries.find(key); + if (it != g_smoothed_texture_cache.entries.end() && it->second.smoothing == layer.smoothing && + it->second.source.lock() == layer.image_data) + return it->second.texture; + } + const float span = 0.05f * float(std::min(result.width, result.height)); + const float radius = std::clamp(layer.smoothing, 0.f, 1.f) * std::min(span, 48.f); smooth_height_pixels(result.pixels, result.width, result.height, radius); // Colour gets the same blur, per channel. It is the same knob for the same reason: detail in // the image finer than the mesh can carry is noise either way, and low-passing it here is the @@ -251,10 +331,131 @@ DecodedHeightTexture decode_height_texture(const TextureDisplacementLayer &layer result.rgb[i * 3 + size_t(c)] = channel[i]; } } + std::lock_guard lock(g_smoothed_texture_cache.mutex); + auto &entries = g_smoothed_texture_cache.entries; + for (auto it = entries.begin(); it != entries.end();) + it = it->second.source.expired() ? entries.erase(it) : std::next(it); + entries[key] = { std::weak_ptr>(layer.image_data), layer.smoothing, result }; } return result; } +namespace { +struct TextureDetailCache +{ + std::mutex mutex; + std::unordered_map>, TextureDetail>> entries; +}; +TextureDetailCache g_texture_detail_cache; +} // namespace + +TextureDetail analyze_texture_detail(const TextureDisplacementLayer &layer) +{ + TextureDetail out; + if (layer.empty()) + return out; + const void *key = layer.image_data.get(); + { + std::lock_guard lock(g_texture_detail_cache.mutex); + auto it = g_texture_detail_cache.entries.find(key); + if (it != g_texture_detail_cache.entries.end() && it->second.first.lock() == layer.image_data) + return it->second.second; + } + // On the image as imported: the Smoothing slider must not move the recommendation around. + TextureDisplacementLayer raw = layer; + raw.smoothing = 0.f; + const DecodedHeightTexture tex = decode_height_texture(raw); + const int w = tex.width, h = tex.height; + if (w >= 3 && h >= 3) { + double sum = 0.0; + size_t sharp = 0, n = 0; + for (int y = 1; y < h - 1; ++y) { + const uint8_t *row = tex.pixels.data() + size_t(y) * size_t(w); + for (int x = 1; x < w - 1; ++x) { + const float dx = 0.5f * (float(row[x + 1]) - float(row[x - 1])); + const float dy = 0.5f * (float(row[x + w]) - float(row[x - w])); + const float mag = std::sqrt(dx * dx + dy * dy); + sum += mag; + sharp += mag > 30.f; + ++n; + } + } + out.mean_gradient = float(sum / double(n)); + out.sharp_fraction = float(sharp) / float(n); + if (out.sharp_fraction > 0.15f || out.mean_gradient > 50.f) out.pixels_per_edge = 1.f; + else if (out.sharp_fraction > 0.05f || out.mean_gradient > 20.f) out.pixels_per_edge = 1.5f; + else if (out.mean_gradient > 8.f) out.pixels_per_edge = 2.5f; + else out.pixels_per_edge = 4.f; + } + std::lock_guard lock(g_texture_detail_cache.mutex); + auto &entries = g_texture_detail_cache.entries; + for (auto it = entries.begin(); it != entries.end();) + it = it->second.first.expired() ? entries.erase(it) : std::next(it); + entries[key] = { std::weak_ptr>(layer.image_data), out }; + return out; +} + +V2Resolution recommend_v2_resolution(const indexed_triangle_set &mesh, + const std::vector &layers, + const Transform3d &volume_to_world) +{ + // BumpMesh's smart resolution, the numbers included: equilateral-cover triangle density, a 16 M + // triangle refinement cap taken at 75 %, a 0.5 mm reference relief for the budget. + constexpr double TRIS_PER_AREA = 2.309, CAP_TRIANGLES = 16e6 * 0.75; + constexpr double EDGE_MIN = 0.05, EDGE_MAX = 5.0; + constexpr double BUDGET_MIN = 10e3, BUDGET_MAX = 2000e3, REF_DEPTH = 0.5, MIN_DEPTH = 0.1; + + V2Resolution out; + // The finest layer decides: the smallest detail edge (texel x pixels per edge) across the layers. + double detail_edge = std::numeric_limits::max(), depth = 0.0; + for (const TextureDisplacementLayer &layer : layers) { + if (layer.empty() || layer.tiling_scale <= 0.f) + continue; + const DecodedHeightTexture &tex = decode_height_texture(layer); + if (tex.width <= 0) + continue; + const TextureDetail detail = analyze_texture_detail(layer); + const double texel = double(layer.tiling_scale) / double(tex.width); + const double edge = texel * double(detail.pixels_per_edge); + if (edge < detail_edge) { + detail_edge = edge; + out.texel_mm = float(texel); + out.pixels_per_edge = detail.pixels_per_edge; + depth = std::abs(double(layer.depth_mm)); + } + } + if (out.texel_mm <= 0.f || mesh.vertices.empty()) + return out; + + // Surface area and diagonal in world mm: the tile is in world mm and the pipeline refines there. + double area = 0.0; + Vec3d bmin = Vec3d::Constant(std::numeric_limits::max()), bmax = -bmin; + std::vector world(mesh.vertices.size()); + for (size_t i = 0; i < world.size(); ++i) { + world[i] = volume_to_world * mesh.vertices[i].cast(); + bmin = bmin.cwiseMin(world[i]); + bmax = bmax.cwiseMax(world[i]); + } + for (const stl_triangle_vertex_indices &t : mesh.indices) + area += 0.5 * (world[size_t(t[1])] - world[size_t(t[0])]).cross(world[size_t(t[2])] - world[size_t(t[0])]).norm(); + const double diag = (bmax - bmin).norm(); + + const double budget_edge = std::sqrt(TRIS_PER_AREA * area / CAP_TRIANGLES); + double edge = std::max(detail_edge, budget_edge); + out.budget_bound = budget_edge > detail_edge; + const double hi = std::max(EDGE_MIN, std::min(EDGE_MAX, diag / 50.0)); + edge = std::clamp(edge, EDGE_MIN, hi); + edge = std::max(EDGE_MIN, std::ceil(edge * 100.0) / 100.0); // up, so the cap holds + out.edge_mm = float(edge); + + const double depth_scale = std::sqrt(REF_DEPTH / std::max(depth, MIN_DEPTH)); + const double target_edge = double(out.pixels_per_edge) * double(out.texel_mm) * depth_scale; + const double raw = TRIS_PER_AREA * area / (target_edge * target_edge); + const double stepped = std::round(raw / 10e3) * 10e3; + out.budget_k = int(std::clamp(stepped, BUDGET_MIN, BUDGET_MAX) / 1000.0); + return out; +} + Vec2f project_planar(const Vec3f &position, const Vec3f &normal) { // Planar-project onto the two axes orthogonal to the dominant component of `normal`. Called @@ -375,15 +576,114 @@ Vec3f face_normal(const indexed_triangle_set &mesh, const stl_triangle_vertex_in return (len > 1e-12f) ? Vec3f(n / len) : Vec3f::UnitZ(); } -// Groups triangles into charts: two triangles sharing an edge join the same chart only if the angle -// between their normals is below `seam_angle_deg`. Everything sharper than that is a seam, and the -// unwrap gets cut there instead of being forced to flatten across it. +// Whether a chart can be laid flat as one island: a topological disk (V - E + F = 1 with a single boundary loop, +// which LSCM requires) whose faces all point within ~100 degrees of their average - past that even a disk folds +// over itself when flattened. Degenerate faces carry no direction and are left out of the normal test. +bool chart_is_flattenable(const indexed_triangle_set &mesh, const std::vector &normals, const std::vector &faces) +{ + constexpr float MIN_NORMAL_COS = -0.17f; // cos(100 deg) + + std::unordered_map edge_use; + std::unordered_map local_vertex; + edge_use.reserve(faces.size() * 2); + local_vertex.reserve(faces.size()); + Vec3f normal_sum = Vec3f::Zero(); + for (const int f : faces) { + const stl_triangle_vertex_indices &tri = mesh.indices[size_t(f)]; + for (int i = 0; i < 3; ++i) { + ++edge_use[undirected_edge_key(tri[i], tri[(i + 1) % 3])]; + local_vertex.emplace(tri[i], int(local_vertex.size())); + } + normal_sum += normals[size_t(f)]; + } + if (int(local_vertex.size()) - int(edge_use.size()) + int(faces.size()) != 1) + return false; + + UnionFind loops(local_vertex.size()); + int boundary_seed = -1; + for (const auto &[key, uses] : edge_use) { + if (uses > 2) + return false; // non-manifold + if (uses == 1) { + const int a = local_vertex[int(key >> 32)], b = local_vertex[int(uint32_t(key))]; + loops.unite(a, b); + boundary_seed = a; + } + } + if (boundary_seed < 0) + return false; // closed + const int loop = loops.find(boundary_seed); + for (const auto &[key, uses] : edge_use) + if (uses == 1 && loops.find(local_vertex[int(key >> 32)]) != loop) + return false; // a second boundary loop: a ring, or a disk with a hole + + const float len = normal_sum.norm(); + if (len < 1e-6f) + return false; + const Vec3f mean = normal_sum / len; + for (const int f : faces) { + const stl_triangle_vertex_indices &tri = mesh.indices[size_t(f)]; + const bool degenerate = (mesh.vertices[tri[1]] - mesh.vertices[tri[0]]).cross(mesh.vertices[tri[2]] - mesh.vertices[tri[0]]) + .squaredNorm() < 1e-20f; + if (!degenerate && normals[size_t(f)].dot(mean) < MIN_NORMAL_COS) + return false; + } + return true; +} + +// Which side of a cut each face of an unflattenable chart goes to (true / false, parallel to `faces`). The cut +// separates the faces' normals along their widest spread, so a tube splits lengthwise into two half-tubes and a +// closed sphere into two hemispheres. Where the normals barely vary (a flat ring) it splits the face centroids +// along their longest axis instead. +std::vector split_chart_sides(const indexed_triangle_set &mesh, const std::vector &normals, const std::vector &faces) +{ + const auto principal_axis = [](const std::vector &samples, Vec3f &mean, Vec3f &axis) { + mean = Vec3f::Zero(); + for (const Vec3f &s : samples) + mean += s; + mean /= float(samples.size()); + Eigen::Matrix3f covariance = Eigen::Matrix3f::Zero(); + for (const Vec3f &s : samples) { + const Vec3f d = s - mean; + covariance += d * d.transpose(); + } + const Eigen::SelfAdjointEigenSolver solver(covariance); + axis = solver.eigenvectors().col(2); // eigenvalues come sorted ascending + return solver.eigenvalues()(2) / float(samples.size()); + }; + + std::vector samples(faces.size()); + for (size_t i = 0; i < faces.size(); ++i) + samples[i] = normals[size_t(faces[i])]; + Vec3f mean, axis; + if (principal_axis(samples, mean, axis) < 1e-3f) { + for (size_t i = 0; i < faces.size(); ++i) { + const stl_triangle_vertex_indices &tri = mesh.indices[size_t(faces[i])]; + samples[i] = (mesh.vertices[tri[0]] + mesh.vertices[tri[1]] + mesh.vertices[tri[2]]) / 3.f; + } + principal_axis(samples, mean, axis); + } + std::vector side(faces.size()); + for (size_t i = 0; i < faces.size(); ++i) + side[i] = (samples[i] - mean).dot(axis) > 0.f; + return side; +} + +// Groups triangles into charts. First two triangles sharing an edge join the same chart only if the angle between +// their normals is below `seam_angle_deg` (and the edge is not a marked seam): everything sharper is a seam. The +// test uses each face's own normal - averaging in its neighbours, as this once did, cuts low-poly flat faces apart, +// because the two triangles of one cube face average in different neighbouring faces. +// +// Then every chart that cannot be laid flat as one island (see chart_is_flattenable()) is cut in two and each +// connected piece checked again: a cylinder's smooth side is one chart by angle but a ring, which no flattening +// can open without a cut, and a sphere is closed. std::vector segment_into_charts(const indexed_triangle_set &mesh, const std::vector &normals, float seam_angle_deg, const std::unordered_set &seam_keys, int &chart_count) { - const int n_faces = int(mesh.indices.size()); - const float cos_threshold = std::cos(std::clamp(seam_angle_deg, 0.f, 180.f) * float(M_PI) / 180.f); + constexpr int MAX_SPLIT_DEPTH = 8; + const int n_faces = int(mesh.indices.size()); + const float cos_threshold = std::cos(std::clamp(seam_angle_deg, 0.f, 180.f) * float(M_PI) / 180.f); // Each shared edge, with the (up to two) faces on it. std::unordered_map> edge_faces; @@ -398,47 +698,102 @@ std::vector segment_into_charts(const indexed_triangle_set &mesh, const std } } - // The chart-join test compares a *neighbourhood-averaged* normal per face - the face plus its - // edge-adjacent neighbours (~5 samples) - rather than the single face normal. On a finely - // tessellated curved surface this stops one noisy triangle from spuriously cutting (or a lone - // near-flat sliver from wrongly merging) a chart, while a genuine sharp crease, where the whole - // neighbourhood on each side agrees, still cuts. This is the "use 5 points, not one" refinement. - std::vector smoothed(mesh.indices.size()); - for (int f = 0; f < n_faces; ++f) { - Vec3f acc = normals[size_t(f)]; - const stl_triangle_vertex_indices &tri = mesh.indices[size_t(f)]; - for (int i = 0; i < 3; ++i) { - const auto it = edge_faces.find(undirected_edge_key(tri[i], tri[(i + 1) % 3])); - if (it == edge_faces.end()) - continue; - const int nb = (it->second.first == f) ? it->second.second : it->second.first; - if (nb >= 0) - acc += normals[size_t(nb)]; - } - smoothed[size_t(f)] = (acc.norm() > 1e-8f) ? Vec3f(acc.normalized()) : normals[size_t(f)]; - } - - UnionFind uf(mesh.indices.size()); + // Faces joined across each edge that does not cut. A face has three edges, so at most three such neighbours. + std::vector> adjacent(static_cast(n_faces), { -1, -1, -1 }); + const auto link = [&adjacent](int f, int nb) { + for (int &slot : adjacent[size_t(f)]) + if (slot < 0) { + slot = nb; + return; + } + }; for (const auto &[key, fp] : edge_faces) { - if (fp.second < 0) + if (fp.second < 0 || fp.first == fp.second) continue; // a boundary edge of the patch, nothing on the far side to join // A manually/auto marked seam always cuts, whatever the dihedral angle - that is exactly // what lets "mark seam" / "cut island" split a chart that is otherwise flat enough to merge. if (!seam_keys.empty() && seam_keys.count(key)) continue; - if (smoothed[size_t(fp.first)].dot(smoothed[size_t(fp.second)]) >= cos_threshold) - uf.unite(fp.first, fp.second); + if (normals[size_t(fp.first)].dot(normals[size_t(fp.second)]) >= cos_threshold) { + link(fp.first, fp.second); + link(fp.second, fp.first); + } } - std::vector chart_of(mesh.indices.size(), -1); - std::unordered_map root_to_chart; - chart_count = 0; - for (int f = 0; f < n_faces; ++f) { - const auto [it, inserted] = root_to_chart.emplace(uf.find(f), chart_count); - if (inserted) - ++chart_count; - chart_of[size_t(f)] = it->second; + // Connected pieces of `faces` over `adjacent`, staying within each face's current `group`. + std::vector group(static_cast(n_faces), 0), visited(static_cast(n_faces), -1); + int visit_pass = 0; + const auto connected_pieces = [&](const std::vector &faces, std::vector> &out) { + const int pass = visit_pass++; + std::vector stack; + for (const int seed : faces) { + if (visited[size_t(seed)] == pass) + continue; + visited[size_t(seed)] = pass; + std::vector piece{ seed }; + stack.assign(1, seed); + while (!stack.empty()) { + const int f = stack.back(); + stack.pop_back(); + for (const int nb : adjacent[size_t(f)]) + if (nb >= 0 && visited[size_t(nb)] != pass && group[size_t(nb)] == group[size_t(f)]) { + visited[size_t(nb)] = pass; + piece.push_back(nb); + stack.push_back(nb); + } + } + out.push_back(std::move(piece)); + } + }; + + std::vector all_faces(static_cast(n_faces)); + std::iota(all_faces.begin(), all_faces.end(), 0); + std::vector> pieces; + connected_pieces(all_faces, pieces); + std::vector, int>> pending; // a piece, and how many cuts produced it + for (std::vector &piece : pieces) + pending.emplace_back(std::move(piece), 0); + + std::vector> charts; + int next_group = 1; + while (!pending.empty()) { + std::vector faces = std::move(pending.back().first); + const int depth = pending.back().second; + pending.pop_back(); + if (faces.size() < 2 || depth >= MAX_SPLIT_DEPTH || chart_is_flattenable(mesh, normals, faces)) { + charts.push_back(std::move(faces)); + continue; + } + const std::vector side = split_chart_sides(mesh, normals, faces); + std::vector halves[2]; + for (size_t i = 0; i < faces.size(); ++i) + halves[side[i] ? 1 : 0].push_back(faces[i]); + if (halves[0].empty() || halves[1].empty()) { + charts.push_back(std::move(faces)); + continue; + } + for (std::vector &half : halves) { + for (const int f : half) + group[size_t(f)] = next_group; + ++next_group; + pieces.clear(); + connected_pieces(half, pieces); + for (std::vector &piece : pieces) + pending.emplace_back(std::move(piece), depth + 1); + } } + + // Chart ids in first-encountered-triangle order, which TextureIsland indexing documents. + std::vector> first_face(charts.size()); + for (size_t c = 0; c < charts.size(); ++c) + first_face[c] = { *std::min_element(charts[c].begin(), charts[c].end()), c }; + std::sort(first_face.begin(), first_face.end()); + + std::vector chart_of(static_cast(n_faces), -1); + for (size_t id = 0; id < first_face.size(); ++id) + for (const int f : charts[first_face[id].second]) + chart_of[size_t(f)] = int(id); + chart_count = int(charts.size()); return chart_of; } @@ -833,22 +1188,107 @@ bool join_chart_placement(const PatchUnwrap &unwrap, const std::vector; + +// Whether two triangles overlap by more than `eps` (separating axis test). Triangles that merely share an edge or +// a corner, as neighbours in a net do, do not. +bool triangles_overlap(const Tri2 &a, const Tri2 &b, float eps) +{ + for (const Tri2 *t : { &a, &b }) + for (int i = 0; i < 3; ++i) { + const Vec2f edge = (*t)[(i + 1) % 3] - (*t)[i]; + const float len = edge.norm(); + if (len < 1e-12f) + continue; + const Vec2f axis(-edge.y() / len, edge.x() / len); + float a_min = std::numeric_limits::max(), a_max = std::numeric_limits::lowest(); + float b_min = a_min, b_max = a_max; + for (int k = 0; k < 3; ++k) { + const float pa = axis.dot(a[k]), pb = axis.dot(b[k]); + a_min = std::min(a_min, pa); + a_max = std::max(a_max, pa); + b_min = std::min(b_min, pb); + b_max = std::max(b_max, pb); + } + if (a_max <= b_min + eps || b_max <= a_min + eps) + return false; + } + return true; +} + +// The triangles already placed in one net, bucketed in a uniform grid so a candidate chart is only tested against +// its neighbourhood. Triangles spanning many cells go into a list that is tested against everything instead. +struct NetGrid +{ + static constexpr int BIG_SPAN = 16; + float cell; + float eps; + std::unordered_map> cells; + std::vector big; + + static uint64_t key(int x, int y) { return (uint64_t(uint32_t(x)) << 32) | uint32_t(y); } + bool range(const Tri2 &t, int &x0, int &y0, int &x1, int &y1) const + { + const Vec2f lo = t[0].cwiseMin(t[1]).cwiseMin(t[2]), hi = t[0].cwiseMax(t[1]).cwiseMax(t[2]); + x0 = int(std::floor(lo.x() / cell)); + y0 = int(std::floor(lo.y() / cell)); + x1 = int(std::floor(hi.x() / cell)); + y1 = int(std::floor(hi.y() / cell)); + return x1 - x0 <= BIG_SPAN && y1 - y0 <= BIG_SPAN; + } + bool overlaps(const Tri2 &t) const + { + for (const Tri2 &b : big) + if (triangles_overlap(t, b, eps)) + return true; + int x0, y0, x1, y1; + if (!range(t, x0, y0, x1, y1)) { + for (const auto &[k, tris] : cells) + for (const Tri2 &b : tris) + if (triangles_overlap(t, b, eps)) + return true; + return false; + } + for (int x = x0; x <= x1; ++x) + for (int y = y0; y <= y1; ++y) + if (const auto it = cells.find(key(x, y)); it != cells.end()) + for (const Tri2 &b : it->second) + if (triangles_overlap(t, b, eps)) + return true; + return false; + } + void insert(const Tri2 &t) + { + int x0, y0, x1, y1; + if (!range(t, x0, y0, x1, y1)) { + big.push_back(t); + return; + } + for (int x = x0; x <= x1; ++x) + for (int y = y0; y <= y1; ++y) + cells[key(x, y)].push_back(t); + } +}; +} // namespace + std::vector compute_connected_net(const PatchUnwrap &unwrap) { - std::vector islands(size_t(std::max(unwrap.chart_count, 0))); - if (unwrap.chart_count <= 1) + const int n = std::max(unwrap.chart_count, 0); + std::vector islands(static_cast(n)); + if (n <= 1) return islands; // Chart adjacency, with one representative shared edge per adjacent pair. const auto edges = build_shared_edges(unwrap); struct PairEdge { ChartEdge a, b; }; std::map, PairEdge> pair_edge; - std::vector> adj(size_t(unwrap.chart_count)); + std::vector> adj(static_cast(n)); for (const auto &[base_edge, list] : edges) { for (size_t i = 0; i < list.size(); ++i) for (size_t j = i + 1; j < list.size(); ++j) { const int c1 = list[i].chart, c2 = list[j].chart; - if (c1 == c2 || c1 < 0 || c2 < 0 || c1 >= unwrap.chart_count || c2 >= unwrap.chart_count) + if (c1 == c2 || c1 < 0 || c2 < 0 || c1 >= n || c2 >= n) continue; const std::pair pk{ std::min(c1, c2), std::max(c1, c2) }; if (pair_edge.count(pk)) @@ -859,63 +1299,146 @@ std::vector compute_connected_net(const PatchUnwrap &unwrap) } } - std::vector placed(size_t(unwrap.chart_count), false); - std::vector pmin(size_t(unwrap.chart_count)), pmax(size_t(unwrap.chart_count)); - const auto chart_bbox = [&](int c, Vec2f &lo, Vec2f &hi) { - lo = Vec2f(std::numeric_limits::max(), std::numeric_limits::max()); - hi = Vec2f(std::numeric_limits::lowest(), std::numeric_limits::lowest()); - for (size_t i = 0; i < unwrap.uvs.size(); ++i) - if (unwrap.vertex_chart[i] == c) { - const Vec2f p = apply_island_transform(unwrap.uvs[i], c, unwrap, islands); - lo = lo.cwiseMin(p); - hi = hi.cwiseMax(p); - } - }; - const auto overlaps_placed = [&](int c, const Vec2f &lo, const Vec2f &hi) { - // A small inset, so charts that merely share an edge (touching bboxes) aren't judged to overlap. - const Vec2f eps = 0.02f * (hi - lo).cwiseAbs(); - for (int o = 0; o < unwrap.chart_count; ++o) - if (placed[size_t(o)] && o != c && - lo.x() + eps.x() < pmax[size_t(o)].x() && hi.x() - eps.x() > pmin[size_t(o)].x() && - lo.y() + eps.y() < pmax[size_t(o)].y() && hi.y() - eps.y() > pmin[size_t(o)].y()) - return true; - return false; + // Per chart: its vertices, its triangles and its flattened area. + std::vector> chart_verts(static_cast(n)), chart_tris(static_cast(n)); + std::vector chart_area(static_cast(n), 0.f); + for (size_t i = 0; i < unwrap.uvs.size(); ++i) + if (const int c = unwrap.vertex_chart[i]; c >= 0 && c < n) + chart_verts[size_t(c)].push_back(int(i)); + float extent_sum = 0.f, extent = 0.f; + for (size_t t = 0; t < unwrap.indices.size(); ++t) { + const stl_triangle_vertex_indices &tri = unwrap.indices[t]; + const int c = unwrap.vertex_chart[size_t(tri[0])]; + if (c < 0 || c >= n) + continue; + chart_tris[size_t(c)].push_back(int(t)); + const Vec2f &p0 = unwrap.uvs[size_t(tri[0])], &p1 = unwrap.uvs[size_t(tri[1])], &p2 = unwrap.uvs[size_t(tri[2])]; + const Vec2f e0 = p1 - p0, e1 = p2 - p0; + chart_area[size_t(c)] += 0.5f * std::abs(e0.x() * e1.y() - e0.y() * e1.x()); + const Vec2f size = p0.cwiseMax(p1).cwiseMax(p2) - p0.cwiseMin(p1).cwiseMin(p2); + extent_sum += std::max(size.x(), size.y()); + extent = std::max({ extent, p0.cwiseAbs().maxCoeff(), p1.cwiseAbs().maxCoeff(), p2.cwiseAbs().maxCoeff() }); + } + const float cell = std::max(extent_sum / float(std::max(unwrap.indices.size(), 1)), 1e-6f); + // Touching neighbours may overlap by rounding: a fraction of a typical triangle, but at least what float + // rounding of a placement at this distance from the origin can produce. + const float eps = std::max(0.02f * cell, 4e-6f * extent); + + const auto placed = [&](const Eigen::Matrix &m, int t) { + const stl_triangle_vertex_indices &tri = unwrap.indices[size_t(t)]; + Tri2 out; + for (int k = 0; k < 3; ++k) + out[size_t(k)] = m.block<2, 2>(0, 0) * unwrap.uvs[size_t(tri[k])] + m.col(2); + return out; }; - // BFS from chart 0 (kept at its packed position), unfolding each newly reached chart onto its parent. - std::queue q; - placed[0] = true; - chart_bbox(0, pmin[0], pmax[0]); - q.push(0); - while (!q.empty()) { - const int p = q.front(); - q.pop(); - for (const int c : adj[size_t(p)]) { - if (placed[size_t(c)]) - continue; - const std::pair pk{ std::min(p, c), std::max(p, c) }; - const auto it = pair_edge.find(pk); - if (it == pair_edge.end()) - continue; - const ChartEdge &pe = (p < c) ? it->second.a : it->second.b; - const ChartEdge &ce = (p < c) ? it->second.b : it->second.a; - const Vec2f pp0 = apply_island_transform(unwrap.uvs[size_t(pe.uv_lo)], p, unwrap, islands); - const Vec2f pp1 = apply_island_transform(unwrap.uvs[size_t(pe.uv_hi)], p, unwrap, islands); - islands[size_t(c)] = solve_edge_alignment(unwrap.uvs[size_t(ce.uv_lo)], unwrap.uvs[size_t(ce.uv_hi)], pp0, - pp1, unwrap.chart_centroid[size_t(c)]); + // Grow a net from the largest chart not yet in one, unfolding each neighbour onto the chart it was reached from + // (bigger neighbours first, so slivers don't claim the good spots) unless its triangles would overlap the net. + // A chart that doesn't fit stays out and roots a net of its own later, so nothing is left in a random spot. + std::vector by_area(static_cast(n)); + std::iota(by_area.begin(), by_area.end(), 0); + std::stable_sort(by_area.begin(), by_area.end(), [&chart_area](int a, int b) { return chart_area[size_t(a)] > chart_area[size_t(b)]; }); - Vec2f lo, hi; - chart_bbox(c, lo, hi); - if (overlaps_placed(c, lo, hi)) { - islands[size_t(c)] = TextureIsland{}; // would collide: leave it where the packing put it - continue; // and don't unfold its subtree off a packed chart + std::vector net_of(static_cast(n), -1); + int net_count = 0; + for (const int root : by_area) { + if (net_of[size_t(root)] >= 0 || chart_tris[size_t(root)].empty()) + continue; + const int net = net_count++; + NetGrid grid{ cell, eps, {}, {} }; + net_of[size_t(root)] = net; + { + const Eigen::Matrix m = island_transform_matrix(root, unwrap, islands); + for (const int t : chart_tris[size_t(root)]) + grid.insert(placed(m, t)); + } + std::queue q; + q.push(root); + while (!q.empty()) { + const int p = q.front(); + q.pop(); + std::vector neighbours = adj[size_t(p)]; + std::stable_sort(neighbours.begin(), neighbours.end(), + [&chart_area](int a, int b) { return chart_area[size_t(a)] > chart_area[size_t(b)]; }); + for (const int c : neighbours) { + if (net_of[size_t(c)] >= 0 || chart_tris[size_t(c)].empty()) + continue; + const auto it = pair_edge.find({ std::min(p, c), std::max(p, c) }); + if (it == pair_edge.end()) + continue; + const ChartEdge &pe = (p < c) ? it->second.a : it->second.b; + const ChartEdge &ce = (p < c) ? it->second.b : it->second.a; + const Vec2f pp0 = apply_island_transform(unwrap.uvs[size_t(pe.uv_lo)], p, unwrap, islands); + const Vec2f pp1 = apply_island_transform(unwrap.uvs[size_t(pe.uv_hi)], p, unwrap, islands); + islands[size_t(c)] = solve_edge_alignment(unwrap.uvs[size_t(ce.uv_lo)], unwrap.uvs[size_t(ce.uv_hi)], pp0, + pp1, unwrap.chart_centroid[size_t(c)]); + + const Eigen::Matrix m = island_transform_matrix(c, unwrap, islands); + std::vector tris; + tris.reserve(chart_tris[size_t(c)].size()); + bool fits = true; + for (const int t : chart_tris[size_t(c)]) { + tris.push_back(placed(m, t)); + if (grid.overlaps(tris.back())) { + fits = false; + break; + } + } + if (!fits) { + islands[size_t(c)] = TextureIsland{}; + continue; + } + for (const Tri2 &t : tris) + grid.insert(t); + net_of[size_t(c)] = net; + q.push(c); } - placed[size_t(c)] = true; - pmin[size_t(c)] = lo; - pmax[size_t(c)] = hi; - q.push(c); } } + + // Shelf-pack the nets side by side, tallest first, the way compute_patch_unwrap() packs charts. + std::vector lo(static_cast(net_count), Vec2f::Constant(std::numeric_limits::max())); + std::vector hi(static_cast(net_count), Vec2f::Constant(std::numeric_limits::lowest())); + for (int c = 0; c < n; ++c) { + const int net = net_of[size_t(c)]; + if (net < 0) + continue; + const Eigen::Matrix m = island_transform_matrix(c, unwrap, islands); + for (const int v : chart_verts[size_t(c)]) { + const Vec2f p = m.block<2, 2>(0, 0) * unwrap.uvs[size_t(v)] + m.col(2); + lo[size_t(net)] = lo[size_t(net)].cwiseMin(p); + hi[size_t(net)] = hi[size_t(net)].cwiseMax(p); + } + } + float total_area = 0.f, widest = 0.f; + for (int k = 0; k < net_count; ++k) { + const Vec2f size = hi[size_t(k)] - lo[size_t(k)]; + total_area += size.x() * size.y(); + widest = std::max(widest, size.x()); + } + const float shelf_width = std::max(widest, std::sqrt(std::max(total_area, 0.f)) * 1.4f); + const float margin = std::max(shelf_width * 0.02f, 1e-4f); + std::vector order(static_cast(net_count)); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), [&](int a, int b) { + return hi[size_t(a)].y() - lo[size_t(a)].y() > hi[size_t(b)].y() - lo[size_t(b)].y(); + }); + std::vector shift(static_cast(net_count), Vec2f::Zero()); + float cursor_x = 0.f, cursor_y = 0.f, row_height = 0.f; + for (const int k : order) { + const Vec2f size = hi[size_t(k)] - lo[size_t(k)]; + if (cursor_x > 0.f && cursor_x + size.x() > shelf_width) { + cursor_x = 0.f; + cursor_y += row_height + margin; + row_height = 0.f; + } + shift[size_t(k)] = Vec2f(cursor_x, cursor_y) - lo[size_t(k)]; + cursor_x += size.x() + margin; + row_height = std::max(row_height, size.y()); + } + for (int c = 0; c < n; ++c) + if (net_of[size_t(c)] >= 0) + islands[size_t(c)].offset += shift[size_t(net_of[size_t(c)])]; return islands; } @@ -931,6 +1454,30 @@ void average_island_scales(std::vector &islands) island.scale = mean; } +std::vector apply_lscm_uv_overrides(PatchUnwrap &unwrap, const std::vector> &overrides) +{ + std::vector edited(unwrap.uvs.size(), false); + if (overrides.empty()) + return edited; + // Mesh-vertex keys first, so an edit on one copy of the same vertex takes precedence over them. + std::map by_mesh_vertex; + for (const auto &[key, uv] : overrides) + if (key >= 0) + by_mesh_vertex[key] = uv; + if (!by_mesh_vertex.empty()) + for (size_t i = 0; i < unwrap.uvs.size(); ++i) + if (const auto it = by_mesh_vertex.find(unwrap.source_vertex[i]); it != by_mesh_vertex.end()) { + unwrap.uvs[i] = it->second; + edited[i] = true; + } + for (const auto &[key, uv] : overrides) + if (const int i = -key - 1; key < 0 && size_t(i) < unwrap.uvs.size()) { + unwrap.uvs[size_t(i)] = uv; + edited[size_t(i)] = true; + } + return edited; +} + std::vector compute_lscm_uvs(const indexed_triangle_set &patch, const TextureDisplacementLayer &layer) { // Padding disabled (0), matching the UV editor: the packed islands the editor shows and the ones the @@ -940,33 +1487,68 @@ std::vector compute_lscm_uvs(const indexed_triangle_set &patch, const Tex if (unwrap.empty()) return {}; - // Manual per-vertex UV edits (UV editor Vertex/Edge modes) override the automatic raw unwrap - // coordinate for a mesh vertex, before the island transform - so the edit rides along with any - // island move/rotate exactly like the rest of the island. See TextureDisplacementLayer:: - // lscm_uv_overrides. Small (hand edits), so a plain map is ample. - std::map overrides; - for (const auto &[mv, uv] : layer.lscm_uv_overrides) - overrides[mv] = uv; + // Manual UV edits (UV editor Vertex/Edge modes) replace the automatic raw unwrap coordinate, before the island + // transform - so the edit rides along with any island move/rotate exactly like the rest of the island. + PatchUnwrap edited_unwrap = unwrap; + const std::vector edited = apply_lscm_uv_overrides(edited_unwrap, layer.lscm_uv_overrides); // One UV per patch vertex: a seam vertex has several (one per chart it touches) and has to // settle on one, since it can only be displaced to a single position. See compute_lscm_uvs()'s - // header comment - the surface stays watertight regardless. + // header comment - the surface stays watertight regardless. A copy that was edited by hand wins, + // so the edit is what bakes; otherwise the first copy does. std::vector per_vertex(patch.vertices.size(), Vec2f::Zero()); std::vector assigned(patch.vertices.size(), false); - for (size_t i = 0; i < unwrap.uvs.size(); ++i) { - const int pv = unwrap.source_vertex[i]; - if (pv >= 0 && !assigned[size_t(pv)]) { - Vec2f raw = unwrap.uvs[i]; - const auto it = overrides.find(pv); - if (it != overrides.end()) - raw = it->second; - per_vertex[size_t(pv)] = apply_island_transform(raw, unwrap.vertex_chart[i], unwrap, layer.islands); + for (const bool edited_pass : { true, false }) + for (size_t i = 0; i < edited_unwrap.uvs.size(); ++i) { + const int pv = edited_unwrap.source_vertex[i]; + if (pv < 0 || assigned[size_t(pv)] || edited[i] != edited_pass) + continue; + per_vertex[size_t(pv)] = apply_island_transform(edited_unwrap.uvs[i], edited_unwrap.vertex_chart[i], unwrap, layer.islands); assigned[size_t(pv)] = true; } - } return per_vertex; } +namespace { +// apply_uv_transform()'s per-layer constants, worked out once. Triplanar sampling runs the transform +// three times per point, and recomputing the rotation's cos/sin and the tiling reciprocal on every one +// of them was most of the per-sample arithmetic. +struct UVTransform +{ + float scale, cs, sn, aspect; + Vec2f offset; + + UVTransform(const TextureDisplacementLayer &layer, float aspect_) + { + scale = (layer.tiling_scale > 1e-6f) ? (1.f / layer.tiling_scale) : 1.f; + const float rad = layer.rotation_deg * float(M_PI) / 180.f; + cs = std::cos(rad); + sn = std::sin(rad); + aspect = aspect_; + offset = layer.offset; + } + Vec2f operator()(const Vec2f &planar) const + { + const Vec2f scaled = planar * scale; + Vec2f rotated(scaled.x() * cs - scaled.y() * sn, scaled.x() * sn + scaled.y() * cs); + // See apply_uv_transform() for why the aspect correction follows the rotation. + if (aspect > 0.f && aspect != 1.f) + rotated.y() *= aspect; + return rotated + offset; + } +}; + +// |n|^TRIPLANAR_BLEND_SHARPNESS per component. The exponent is a compile-time 4, so two squarings do +// what three std::pow calls per sample did. +static_assert(TRIPLANAR_BLEND_SHARPNESS == 4.f, "triplanar_weights() hard-codes the fourth power"); +inline Vec3f triplanar_weights(const Vec3f &normal) +{ + Vec3f w = normal.cwiseAbs(); + w = w.cwiseProduct(w); + return w.cwiseProduct(w); +} +} // namespace + Vec2f apply_uv_transform(const Vec2f &planar, const TextureDisplacementLayer &layer, float aspect) { const float scale = (layer.tiling_scale > 1e-6f) ? (1.f / layer.tiling_scale) : 1.f; @@ -1040,9 +1622,10 @@ float sample_layer_height(const DecodedHeightTexture &texture, const TextureDisp // width / height of the height map, so a non-square image keeps its proportions (see // apply_uv_transform()). Every projection except the projective "from view" one funnels through // here, so this one line is what makes them all aspect-correct. - const float aspect = (texture.height > 0) ? float(texture.width) / float(texture.height) : 1.f; + const float aspect = (texture.height > 0) ? float(texture.width) / float(texture.height) : 1.f; + const UVTransform xf(layer, aspect); auto sample_at = [&](const Vec2f &planar) { - return texture.sample(apply_uv_transform(planar, layer, aspect), layer.tile_enabled, layer.tile_method); + return texture.sample(xf(planar), layer.tile_enabled, layer.tile_method); }; // Precomputed per-patch LSCM solve wins over the layer's own method (see the header): the @@ -1080,9 +1663,7 @@ float sample_layer_height(const DecodedHeightTexture &texture, const TextureDisp // hard switch is what produced a visible seam wherever the dominant axis flips (see // TextureProjectionMethod::Triplanar); a weighted blend is continuous across that transition // by construction, since the weight of the axis being left behind falls smoothly to zero. - Vec3f w = normal.cwiseAbs(); - w = Vec3f(std::pow(w.x(), TRIPLANAR_BLEND_SHARPNESS), std::pow(w.y(), TRIPLANAR_BLEND_SHARPNESS), - std::pow(w.z(), TRIPLANAR_BLEND_SHARPNESS)); + Vec3f w = triplanar_weights(normal); const float w_sum = w.x() + w.y() + w.z(); if (w_sum < 1e-8f) // Degenerate normal: no axis is meaningfully dominant, so no blend is meaningful either. @@ -1107,17 +1688,17 @@ bool sample_layer_color(const DecodedHeightTexture &texture, const TextureDispla // the two differ in what "nothing here" means. Height returns 0, which is a perfectly good height // (no displacement); colour has no such neutral value - black is a colour - so every path that // returns 0 there has to report false here instead, and the caller leaves the triangle uncoloured. - const float aspect = (texture.height > 0) ? float(texture.width) / float(texture.height) : 1.f; + const float aspect = (texture.height > 0) ? float(texture.width) / float(texture.height) : 1.f; + const UVTransform xf(layer, aspect); auto sample_at = [&](const Vec2f &planar) { - return texture.sample_color(apply_uv_transform(planar, layer, aspect), layer.tile_enabled, - layer.tile_method); + return texture.sample_color(xf(planar), layer.tile_enabled, layer.tile_method); }; // Outside a non-tiled placement there is no texture at all - the same hard edge sample() gives the // height. Checked explicitly because sample_color() reports it as black, which is a real colour. auto covered = [&](const Vec2f &planar) { if (layer.tile_enabled) return true; - const Vec2f uv = apply_uv_transform(planar, layer, aspect); + const Vec2f uv = xf(planar); return uv.x() >= 0.f && uv.x() < 1.f && uv.y() >= 0.f && uv.y() < 1.f; }; @@ -1168,9 +1749,7 @@ bool sample_layer_color(const DecodedHeightTexture &texture, const TextureDispla // Blended tri-planar, weighted exactly as the height is, so colour and relief stay registered // across the cross-fade band at a 90-degree edge. - Vec3f w = normal.cwiseAbs(); - w = Vec3f(std::pow(w.x(), TRIPLANAR_BLEND_SHARPNESS), std::pow(w.y(), TRIPLANAR_BLEND_SHARPNESS), - std::pow(w.z(), TRIPLANAR_BLEND_SHARPNESS)); + Vec3f w = triplanar_weights(normal); const float w_sum = w.x() + w.y() + w.z(); if (w_sum < 1e-8f) { const Vec2f p(position.x(), position.y()); @@ -1372,7 +1951,10 @@ indexed_triangle_set build_texture_displacement_v2(const indexed_triangle_set const std::vector &layers, const TextureDisplacementFacetsData &facets_data, const TextureDisplacementOptions &options, - const DisplacementProgressFn &progress) + const DisplacementProgressFn &progress, + const TextureColorRequest *color, + bool flip_normals, + BakeStageRecorder *debug) { HeightFieldSampler combined = make_combined_displacement_sampler(mesh, layers, facets_data); if (!combined) @@ -1400,22 +1982,36 @@ indexed_triangle_set build_texture_displacement_v2(const indexed_triangle_set if (std::all_of(excluded.begin(), excluded.end(), [](uint8_t e) { return e != 0; })) return mesh; // nothing painted + // "Auto" resolution and budget (0 and -1) resolve here, from the texture and the model - the mesh + // is already in world mm at this point, so no transform is needed. + const bool auto_edge = options.v2_refine_mm <= 0.f, auto_budget = options.v2_max_triangles_k < 0; + const V2Resolution rec = (auto_edge || auto_budget) ? recommend_v2_resolution(mesh, layers) : V2Resolution{}; TextureBake::PipelineSettings settings; - settings.refine_length = std::max(0.01f, options.v2_refine_mm); + settings.refine_length = auto_edge ? std::max(0.05, double(rec.edge_mm)) : std::max(0.01f, options.v2_refine_mm); settings.regularize = options.v2_regularize; - settings.max_triangles = size_t(std::max(0, options.v2_max_triangles_k)) * 1000; + settings.max_triangles = size_t(auto_budget ? std::max(0, rec.budget_k) : options.v2_max_triangles_k) * 1000; + BOOST_LOG_TRIVIAL(info) << "TextureBake resolution: " << settings.refine_length << " mm" + << (auto_edge ? " (auto)" : "") << ", budget " << settings.max_triangles / 1000 << " k" + << (auto_budget ? " (auto)" : ""); settings.preserve_untextured = true; - settings.clamp_below_plate = options.v2_clamp_below_plate; + // Always on: relief driven under the plate is unprintable whichever pipeline produced it, so this + // is no longer a choice the user has to make. Only geometry that ends up below the model's own + // bottom is moved - downward relief that stays clear of the plate is untouched. + settings.clamp_below_plate = true; settings.relocate = options.v2_relocate; + settings.flip_edges = options.v2_flip_edges; // The sampler already returns millimetres, so the displacement stage must not scale it again. settings.displace.amplitude = 1.f; settings.displace.symmetric = false; // The paint decides what moves here, so the angle limits stay off. settings.displace.bottom_angle_limit = 0.f; settings.displace.top_angle_limit = 0.f; - // The sampler above takes the smooth normal, never the smoothed blend normal, so computing the - // latter would build a full adjacency graph and run its iterations for a value nothing reads. - settings.displace.blend_normal_smoothing = 0; + // The sampler's normal only picks the projection blend (triplanar weights); the displacement + // direction is the pipeline's own smooth normal. Handing it the Laplacian-smoothed blend normal + // spreads a crease's 50/50 blend over a band of rows instead of one: on a cube edge the single row + // of vertices that samples both faces' patterns half and half otherwise comes out as a row of + // notches, since it matches neither face. + settings.displace.blend_normal_smoothing = 32; TextureBake::DisplaceBounds bounds; bounds.min = bounds.max = mesh.vertices.empty() ? Vec3f::Zero() : mesh.vertices.front(); @@ -1424,39 +2020,130 @@ indexed_triangle_set build_texture_displacement_v2(const indexed_triangle_set bounds.max = bounds.max.cwiseMax(v); } - const auto sample = [&combined](const Vec3f &pos, const Vec3f &smooth_normal, const Vec3f &) { - // The smooth normal: the direction the vertex actually moves along. - return combined(pos, smooth_normal); + const auto sample = [&combined](const Vec3f &pos, const Vec3f &, const Vec3f &blend_normal) { + // The blend normal chooses the projection; the move itself is along the smooth normal, which + // the pipeline applies on its own. + return combined(pos, blend_normal); }; + // The pipeline takes its displacement direction from the soup's winding, so a mirrored placement + // would drive the whole relief inwards. The paint masks were read off `mesh` above, against its + // own vertex order, so the winding can only be turned round after that - here, on the copy that + // becomes the soup - and has to be turned back on the way out, since the caller undoes the same + // mirror when it maps the result back into the volume's coordinates. + indexed_triangle_set oriented = mesh; + if (flip_normals) + for (stl_triangle_vertex_indices &t : oriented.indices) + std::swap(t[1], t[2]); + // 0 means no simplification, i.e. Bake mode. const TextureBake::PipelineMode mode = settings.max_triangles > 0 ? TextureBake::PipelineMode::Export : TextureBake::PipelineMode::Bake; - TextureBake::PipelineResult result = TextureBake::run_pipeline( - TextureBake::to_soup(mesh, excluded), sample, settings, bounds, mode, excluded, + // The pipeline works on `oriented`, whose winding was reversed above for a mirrored placement, so + // the stages it records are wound the same way. Note where they start and turn the whole range + // back afterwards, exactly as the result itself is turned back below. + const size_t debug_mark = (debug != nullptr) ? debug->mark() : 0; + TextureBake::PipelineResult result = TextureBake::run_pipeline( + TextureBake::to_soup(oriented, excluded), sample, settings, bounds, mode, excluded, [&progress](const char *, double f) { return !progress || progress(std::clamp(int(f * 100.0), 0, 99)); - }); + }, + debug); + if (debug != nullptr && flip_normals) + debug->rebase(debug_mark, nullptr, /* flip_winding */ true); if (result.canceled || result.geometry.empty()) return {}; indexed_triangle_set out = TextureBake::to_indexed_triangle_set(result.geometry); - return out.indices.empty() ? mesh : out; + if (out.indices.empty()) + return mesh; + if (flip_normals) + for (stl_triangle_vertex_indices &t : out.indices) + std::swap(t[1], t[2]); + + // Colour, per output triangle. The topology is new, so unlike the classic path there is no base + // triangle to inherit a colour from: each output triangle samples the colour stack at its own + // centroid, and takes colour only where the base surface under it is painted - found by the nearest + // base triangle, which is never more than the relief depth away. Then the same despeckle and + // filament resolution as the classic path. + if (color != nullptr && color->out_triangle != nullptr && bool(color->quantize)) { + std::vector out_color(out.indices.size(), 0); + const ColorFieldSampler sampler = make_combined_color_sampler(mesh, layers, facets_data, color->quantize); + if (sampler) { + const bool all_painted = std::none_of(excluded.begin(), excluded.end(), [](uint8_t e) { return e != 0; }); + AABBTreeIndirect::Tree3f tree; + if (!all_painted) + tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(mesh.vertices, mesh.indices); + std::vector palette(out.indices.size(), -1); + tbb::parallel_for(tbb::blocked_range(0, out.indices.size()), [&](const tbb::blocked_range &r) { + for (size_t i = r.begin(); i < r.end(); ++i) { + const stl_triangle_vertex_indices &t = out.indices[i]; + const Vec3f &a = out.vertices[size_t(t[0])], &b = out.vertices[size_t(t[1])], &c = out.vertices[size_t(t[2])]; + const Vec3f centroid = (a + b + c) / 3.f; + if (!all_painted) { + size_t hit = 0; + Vec3f hit_point; + AABBTreeIndirect::squared_distance_to_indexed_triangle_set(mesh.vertices, mesh.indices, tree, + centroid, hit, hit_point); + if (hit >= excluded.size() || excluded[hit] != 0) + continue; + } + Vec3f n = (b - a).cross(c - a); + const float l = n.norm(); + n = (l > 0.f) ? Vec3f(n / l) : Vec3f::UnitZ(); + palette[i] = sampler(centroid, n); + } + }); + despeckle_triangle_colors(out, palette, color->despeckle_passes); + for (size_t i = 0; i < out.indices.size(); ++i) { + if (palette[i] < 0) + continue; + const stl_triangle_vertex_indices &t = out.indices[i]; + const Vec3f centroid = (out.vertices[size_t(t[0])] + out.vertices[size_t(t[1])] + out.vertices[size_t(t[2])]) / 3.f; + const int filament = color->resolve ? color->resolve(palette[i], centroid) : palette[i]; + if (filament >= 0) + out_color[i] = uint8_t(std::min(filament + 1, 255)); + } + } + *color->out_triangle = std::move(out_color); + } + return out; } } // namespace -indexed_triangle_set build_texture_displacement(const indexed_triangle_set &base_mesh, - const std::vector &layers, - const TextureDisplacementFacetsData &facets_data, - const TextureDisplacementOptions &options, - const DisplacementProgressFn &progress, - const TextureColorRequest *color) +// The bake proper. Runs entirely in whatever space `base_mesh` is given in; the public entry point +// below is what puts it in world space and brings the result back. +// `flip_normals` says the mesh is wound the opposite way round from its outward direction, which is +// what a mirroring world transform leaves behind: the positions are right, but every normal derived +// from the winding points into the model. See build_texture_displacement(). +static indexed_triangle_set build_texture_displacement_in_place( + const indexed_triangle_set &base_mesh, + const std::vector &layers, + const TextureDisplacementFacetsData &facets_data, + const TextureDisplacementOptions &options, + const DisplacementProgressFn &progress, + const TextureColorRequest *color, + bool flip_normals, + BakeStageRecorder *debug) { // Returns true to keep going. An aborted run returns {} (see the header): an empty mesh is the // one result no caller can mistake for a finished bake and commit onto the volume. const auto report = [&progress](int percent) { return !progress || progress(percent); }; + // Stage capture for the debug view. A stage is timed from the end of the previous one, so the + // capture itself - a copy plus an edge scan - sits outside every measurement it reports. + auto stage_clock = std::chrono::steady_clock::now(); + const auto capture = [&](const char *name, const indexed_triangle_set &m, + const std::string &detail = {}) { + if (debug == nullptr) + return; + const double ms = + std::chrono::duration(std::chrono::steady_clock::now() - stage_clock).count(); + debug->capture(name, m.vertices, m.indices, ms, detail); + stage_clock = std::chrono::steady_clock::now(); + }; + indexed_triangle_set mesh = base_mesh; // TriangleSelector's vertex array starts with the mesh's own vertices (any extra ones, created // where a brush stroke split a triangle, are appended after them), and get_facets_strict() @@ -1469,8 +2156,11 @@ indexed_triangle_set build_texture_displacement(const indexed_triangle_set if (mesh.vertices.empty() || mesh.indices.empty()) return mesh; + capture("input", mesh, "as handed to the bake"); + if (options.pipeline_v2) - return build_texture_displacement_v2(mesh, layers, facets_data, options, progress); + return build_texture_displacement_v2(mesh, layers, facets_data, options, progress, color, flip_normals, + debug); // Layers are combined in slot order, like stacked layers in an image editor: each one folds its // own displacement into the running total via its blend mode (see TextureBlendMode). @@ -1523,6 +2213,14 @@ indexed_triangle_set build_texture_displacement(const indexed_triangle_set } } + // Both normal passes above read their direction out of the triangle winding, so a mirrored + // placement leaves every one of them pointing into the model - the relief would be carved rather + // than raised. Correct them once, here, where every later stage (displacement direction, the + // planar/triplanar projection axes, the cylinder axis) picks them up already right. + if (flip_normals) + for (Vec3f &n : vertex_normals) + n = -n; + if (!report(5)) return {}; @@ -1787,9 +2485,19 @@ indexed_triangle_set build_texture_displacement(const indexed_triangle_set if (!any_displacement) return mesh; + // The model's own resting plane, taken before anything moves - see the clamp after smoothing. + float resting_z = std::numeric_limits::max(); + for (const Vec3f &v : mesh.vertices) + resting_z = std::min(resting_z, v.z()); + + size_t moved_count = 0; for (size_t vi = 0; vi < mesh.vertices.size(); ++vi) - if (displaced[vi]) + if (displaced[vi]) { mesh.vertices[vi] += vertex_normals[vi] * displacement[vi]; + ++moved_count; + } + capture("displace", mesh, std::to_string(moved_count) + " of " + + std::to_string(mesh.vertices.size()) + " vertices moved"); // Post-process relaxation of what the height maps left behind, restricted to the vertices that // actually moved - the untouched part of the model keeps its exact geometry, and the ring of @@ -1810,8 +2518,28 @@ indexed_triangle_set build_texture_displacement(const indexed_triangle_set progress ? DisplacementProgressFn([&report, it = options.smooth_iterations](int pass) { return report(70 + (29 * (pass + 1)) / std::max(it, 1)); }) : DisplacementProgressFn{}); + capture("smooth", mesh, + std::to_string(options.smooth_iterations) + " iterations at strength " + + std::to_string(options.smooth_strength)); } + // Nothing driven below the model's own resting plane can be printed: it is either through the + // build plate or, once the slicer drops the part back down onto it, holding the whole model up in + // the air. Push it back up to the plane. Only vertices the displacement actually moved are + // eligible - untouched geometry is already exactly where it started - and only those that ended + // up below it, so downward relief that stays clear of the plate is left alone. Runs in world + // space (see build_texture_displacement()), so this really is the plate and not some scaled + // stand-in for it. + size_t clamped_count = 0; + if (resting_z < std::numeric_limits::max()) + for (size_t vi = 0; vi < mesh.vertices.size(); ++vi) + if (displaced[vi] && mesh.vertices[vi].z() < resting_z) { + mesh.vertices[vi].z() = resting_z; + ++clamped_count; + } + capture("clamp to plate", mesh, std::to_string(clamped_count) + " vertices raised to z = " + + std::to_string(resting_z)); + if (!report(99)) return {}; if (want_color) { @@ -1840,6 +2568,70 @@ indexed_triangle_set build_texture_displacement(const indexed_triangle_set return mesh; } +indexed_triangle_set build_texture_displacement(const indexed_triangle_set &base_mesh, + const std::vector &layers, + const TextureDisplacementFacetsData &facets_data, + const TextureDisplacementOptions &options, + const DisplacementProgressFn &progress, + const TextureColorRequest *color, + const Transform3d &volume_to_world, + BakeStageRecorder *debug) +{ + // An untransformed volume on an untransformed instance is by far the common case, and the round + // trip costs two matrix multiplies per vertex on a mesh that can carry millions of them - so take + // the identity out of the way rather than paying for it. + // The texture is projected in the bake frame - world orientation and scale, anchored at the + // volume's origin (see texture_displacement_bake_frame()). + const Transform3d frame = texture_displacement_bake_frame(volume_to_world); + if (frame.matrix().isApprox(Transform3d::Identity().matrix())) + return build_texture_displacement_in_place(base_mesh, layers, facets_data, options, progress, color, + false, debug); + + const Transform3d to_local = frame.inverse(); + // A mirroring placement leaves the positions correct but every winding-derived normal pointing + // the wrong way. The winding itself is deliberately *not* touched here: the paint masks encode + // each split triangle against its own vertex order, so reordering a triangle's vertices would + // mirror the paint inside it. The bake is told instead, and negates the normals it computes. + const bool mirrored = frame.linear().determinant() < 0.0; + + indexed_triangle_set world = base_mesh; + for (Vec3f &v : world.vertices) + v = (frame * v.cast()).cast(); + + const size_t debug_mark = (debug != nullptr) ? debug->mark() : 0; + indexed_triangle_set out = build_texture_displacement_in_place(world, layers, facets_data, options, + progress, color, mirrored, debug); + // Everything the bake recorded is in world millimetres, like `out` itself. The debug view draws in + // the volume's local frame, so the stages are brought back the same way the result is. + if (debug != nullptr) + debug->rebase(debug_mark, &to_local, /* flip_winding */ false); + // A cancelled run returns {} and must stay {} - an empty mesh is the signal the caller checks + // before committing anything onto the volume. + if (out.vertices.empty()) + return out; + + for (Vec3f &v : out.vertices) + v = (to_local * v.cast()).cast(); + return out; +} + +Transform3d texture_displacement_bake_frame(const Transform3d &volume_to_world) +{ + // World orientation and scale, but the origin moved to where the volume's own origin sits: the + // texture then rides with the model when it is moved about the plate (and a single, untiled stamp + // starts on the model rather than at the plate's corner), while a tile is still `tiling_scale` + // printed millimetres whatever the instance's scale. + return Eigen::Translation3d(-volume_to_world.translation()) * volume_to_world; +} + +Transform3d texture_displacement_volume_to_world(const ModelVolume &volume) +{ + const ModelObject *object = volume.get_object(); + if (object == nullptr || object->instances.empty() || object->instances.front() == nullptr) + return volume.get_matrix(); + return object->instances.front()->get_matrix() * volume.get_matrix(); +} + indexed_triangle_set build_texture_displacement(const ModelVolume &volume) { TextureDisplacementFacetsData facets_data; @@ -1847,7 +2639,8 @@ indexed_triangle_set build_texture_displacement(const ModelVolume &volume) facets_data[size_t(i)] = volume.texture_displacement_facet(i).get_data(); return build_texture_displacement(volume.mesh().its, volume.texture_displacement_layers, facets_data, - volume.texture_displacement_options); + volume.texture_displacement_options, {}, nullptr, + texture_displacement_volume_to_world(volume)); } void smooth_mesh_vertices(indexed_triangle_set &mesh, const std::vector &movable, float strength, @@ -2099,7 +2892,8 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, float chord_tolerance_mm, float min_edge_length_mm, float border_edge_length_mm, const DisplacementProgressFn &progress, - const ColorFieldSampler &color, float color_edge_length_mm) + const ColorFieldSampler &color, float color_edge_length_mm, + bool split_multi_crossings) { // Neighbour slots that are not a triangle index. constexpr int NB_BOUNDARY = -1; // open edge: terminal on its own, bisected from this side alone @@ -2136,6 +2930,9 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, // driven by the length baseline alone. const bool feature_mode = bool(sampler) && chord_tolerance_mm > 0.f; const bool color_mode = bool(color) && color_edge_length_mm > 0.f; + const bool multi_mode = feature_mode && split_multi_crossings; + float step_iso_mm = 0.f; // the mid-level a step crosses; set once the region is sampled + float step_dev_mm = std::numeric_limits::max(); // deviations above this are steps const float color_sq = color_edge_length_mm > 0.f ? color_edge_length_mm * color_edge_length_mm : 0.f; const float min_floor_sq = min_edge_length_mm > 0.f ? min_edge_length_mm * min_edge_length_mm : 0.f; const float target_sq = target_edge_length_mm > 0.f ? target_edge_length_mm * target_edge_length_mm : 0.f; @@ -2284,6 +3081,14 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, } } }); + if (multi_mode) { + // The level a step crosses: midway through the relief actually present on the region. + float lo = std::numeric_limits::max(), hi = -lo; + for (size_t v = 0; v < verts.size(); ++v) + if (wanted[v]) { lo = std::min(lo, vheight[v]); hi = std::max(hi, vheight[v]); } + step_iso_mm = (lo < hi) ? 0.5f * (lo + hi) : 0.f; + step_dev_mm = (lo < hi) ? 0.4f * (hi - lo) : std::numeric_limits::max(); + } } auto elen_sq = [&](int a, int b) -> float { return (verts[a] - verts[b]).squaredNorm(); }; @@ -2349,6 +3154,12 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, return split; }; + // The spacing at which the chord test samples a triangle. Features narrower than the resolution + // floor cannot be resolved by refinement anyway, so sampling any finer than that only costs time; + // sampling any coarser is what let thin features fall between the samples. Below the floor's own + // scale it is held at 0.05 mm, since finer is invisible on an FDM part. + const float sample_spacing = std::max(min_edge_length_mm, 0.05f); + auto detail_error = [&](int ti) -> float { if (tri_err[ti] >= 0.f) return tri_err[ti]; @@ -2356,16 +3167,73 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, const Vec3f pa = verts[t.v[0]], pb = verts[t.v[1]], pc = verts[t.v[2]]; const Vec3f na = vnormal[t.v[0]], nb = vnormal[t.v[1]], nc = vnormal[t.v[2]]; const float ha = height_of(t.v[0]), hb = height_of(t.v[1]), hc = height_of(t.v[2]); - static const float BARY[4][3] = { { 0.5f, 0.5f, 0.f }, { 0.f, 0.5f, 0.5f }, - { 0.5f, 0.f, 0.5f }, { 1.f / 3, 1.f / 3, 1.f / 3 } }; - float maxerr = 0.f; - for (const auto &w : BARY) { - Vec3f n = w[0] * na + w[1] * nb + w[2] * nc; - const float nl = n.norm(); - n = (nl > 1e-12f) ? Vec3f(n / nl) : na; - const float actual = sampler(w[0] * pa + w[1] * pb + w[2] * pc, n); - maxerr = std::max(maxerr, std::abs(actual - (w[0] * ha + w[1] * hb + w[2] * hc))); + + // A barycentric lattice over the whole triangle, dense enough that no feature wider than the + // sample spacing can hide between two samples. Four fixed samples (the three edge midpoints + // and the centroid) were the previous test, and on a 1 mm starting triangle they left 0.5 mm + // gaps - wider than a grid line or a knurl ridge - so whole features were never seen, the + // triangle scored zero, and it was never refined. The lattice is bounded at 8 subdivisions + // (45 points) so a large triangle does not cost thousands of samples; anything it misses at + // that size is caught once its children are small enough for the lattice to reach it. + const float longest = std::sqrt(std::max({ (pb - pa).squaredNorm(), (pc - pb).squaredNorm(), (pa - pc).squaredNorm() })); + const int n = std::clamp(int(std::ceil(longest / sample_spacing)), 2, 8); + float maxerr = 0.f; + // Whether the mid-level contour crosses an edge with a jump that is a step's worth, walked at + // the sample spacing (the lattice below is too coarse on a long edge: a 1 mm edge gets 8 + // samples, and a 0.15 mm grid line slips between them). Such a triangle is the cutter's. + bool sharp_cross = false; + if (multi_mode) { + const Vec3f *P[3] = { &pa, &pb, &pc }; + const Vec3f *Nn[3] = { &na, &nb, &nc }; + const float H[3] = { ha, hb, hc }; + for (int e = 0; e < 3; ++e) { + const Vec3f &p0 = *P[e], &p1 = *P[(e + 1) % 3], &n0 = *Nn[e], &n1 = *Nn[(e + 1) % 3]; + const int m = std::clamp(int(std::ceil((p1 - p0).norm() / sample_spacing)), 2, 32); + float last = H[e]; + for (int i = 1; i <= m; ++i) { + float h; + if (i == m) h = H[(e + 1) % 3]; + else { + const float t = float(i) / float(m); + Vec3f nn = n0 + (n1 - n0) * t; + const float nl = nn.norm(); + nn = (nl > 1e-12f) ? Vec3f(nn / nl) : n0; + h = sampler(p0 + (p1 - p0) * t, nn); + } + if ((h > step_iso_mm) != (last > step_iso_mm) && std::abs(h - last) > step_dev_mm) { + sharp_cross = true; + break; + } + last = h; + } + if (sharp_cross) + break; + } } + for (int i = 0; i <= n; ++i) + for (int j = 0; j <= n - i; ++j) { + const int k = n - i - j; + const float wa = float(i) / float(n), wb = float(j) / float(n), wc = float(k) / float(n); + float actual; + if (i == n) actual = ha; + else if (j == n) actual = hb; + else if (k == n) actual = hc; + else { + Vec3f nn = wa * na + wb * nb + wc * nc; + const float nl = nn.norm(); + nn = (nl > 1e-12f) ? Vec3f(nn / nl) : na; + actual = sampler(wa * pa + wb * pb + wc * pc, nn); + // Corners lie on the plane by construction; only the rest carry chord error. + maxerr = std::max(maxerr, std::abs(actual - (wa * ha + wb * hb + wc * hc))); + } + } + // A sharp step crossing the triangle is the cutter's job. Refinement cannot bring such a + // triangle under tolerance - a discontinuity has no chord - and would only carpet the step + // down to the floor, so its chord error is not counted. A feature lying entirely inside the + // triangle crosses no edge and still counts, which is what makes it surface until it does + // cross one. + if (multi_mode && sharp_cross) + maxerr = 0.f; tri_err[ti] = maxerr; return maxerr; }; @@ -2387,8 +3255,9 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, float p = 0.f; if (flags & REFINE_PAINTED) { p = (target_sq > 0.f) ? ll / target_sq : 0.f; - if (feature_mode) + if (feature_mode) { p = std::max(p, detail_error(ti) / chord_tolerance_mm); + } // Colour is per facet, so a colour boundary can only be drawn where there are edges along // it. Length target, floored by min_edge_length_mm above, exactly like the border band. if (color_mode && straddles_color(ti)) @@ -2404,6 +3273,18 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, return p; }; + // What ordering the heap by. priority() alone saturates: a triangle straddling a hard step in the + // texture scores relief / tolerance however small it gets, so its children always come back to the + // top and the whole budget is spent carpeting step edges down to the floor while a large triangle + // with a moderate error a few millimetres away is never reached. Weighting by area makes the key + // the surface error a split actually removes, which is what a fixed budget should be spent on - + // and once the budget is gone, what remains is spread evenly instead of piled on one feature. + auto tri_area = [&](int ti) -> float { + const Tri &t = tris[ti]; + return 0.5f * (verts[t.v[1]] - verts[t.v[0]]).cross(verts[t.v[2]] - verts[t.v[0]]).norm(); + }; + auto heap_key = [&](int ti, float p) -> float { return p * tri_area(ti); }; + auto set_nb = [&](int ti, int u, int v, int val) { if (ti < 0) return; @@ -2536,7 +3417,7 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, }); for (int ti = 0; ti < int(tris.size()); ++ti) if (initial[size_t(ti)] > 1.f) - queue.emplace(initial[size_t(ti)], ti); + queue.emplace(heap_key(ti, initial[size_t(ti)]), ti); } // Every iteration either drops one satisfied triangle from the queue or performs exactly one @@ -2584,15 +3465,1184 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, for (const int t : touched) if (const float p = priority(t); p > 1.f) - queue.emplace(p, t); + queue.emplace(heap_key(t, p), t); // The bisection may have been a step on the way to ti rather than ti itself, in which case ti // is not in `touched` and has to go back on the heap to be walked again. if (cur != ti) if (const float p = priority(ti); p > 1.f) - queue.emplace(p, ti); + queue.emplace(heap_key(ti, p), ti); } return emit(); } +namespace { + +// Ear clipping of a simple polygon given counter-clockwise in 2D; `pts` are indexed by `poly`. Small +// polygons only (a triangle clipped by a few contour polylines), so the plain O(n^2) form is enough. +// A vertex that lies on a candidate ear's boundary does not block it, and when no ear can be found +// (numerically degenerate input) the flattest vertex is dropped rather than fanning blindly, which +// on a concave region would put triangles outside it; the zero-area triangle keeps the edges paired. +void ear_clip(const std::vector &pts, std::vector poly, const std::function &emit) +{ + auto cross2 = [](const Vec2f &o, const Vec2f &a, const Vec2f &b) { + return (a.x() - o.x()) * (b.y() - o.y()) - (a.y() - o.y()) * (b.x() - o.x()); + }; + // A corner flatter than this is collinear: several vertices along one edge of the triangle. Such + // a corner is never an ear - the zero-area triangle would lie on the edge, where the neighbour, + // seeing the same points from its side, may emit the same one. + float flat = 0.f; + { + Vec2f lo = pts[size_t(poly[0])], hi = lo; + for (int id : poly) { lo = lo.cwiseMin(pts[size_t(id)]); hi = hi.cwiseMax(pts[size_t(id)]); } + flat = 1e-5f * (hi - lo).squaredNorm(); + } + while (poly.size() > 3) { + const size_t n = poly.size(); + bool found = false; + for (size_t i = 0; i < n && !found; ++i) { + const int ia = poly[(i + n - 1) % n], ib = poly[i], ic = poly[(i + 1) % n]; + const Vec2f &a = pts[size_t(ia)], &b = pts[size_t(ib)], &c = pts[size_t(ic)]; + const float area2 = cross2(a, b, c); + if (area2 <= flat) + continue; // reflex or flat corner + const float eps = 1e-6f * area2; + bool clean = true; + size_t on_base = 0; + for (size_t j = 0; j < n && clean; ++j) { + const int ij = poly[j]; + if (ij == ia || ij == ib || ij == ic) + continue; + const Vec2f &p = pts[size_t(ij)]; + clean = !(cross2(a, b, p) > eps && cross2(b, c, p) > eps && cross2(c, a, p) > eps); + // On the ear's base, between its ends: a run of vertices along the base, as when + // several sit on one edge of the triangle. + if (std::abs(cross2(c, a, p)) <= flat && (p - c).dot(a - c) > 0.f && (p - a).dot(c - a) > 0.f) + ++on_base; + } + if (!clean) + continue; + if (on_base > 0) { + // The base carries other vertices. When those are all that is left of the polygon, + // in order from c round to a, fan them from b and finish; otherwise the ear is not + // clear of the rest of the polygon. + if (on_base + 3 != n) + continue; + for (size_t j = 0; j + 1 < n - 1; ++j) { + const int u = poly[(i + 1 + j) % n], v = poly[(i + 2 + j) % n]; + emit(ib, u, v); + } + poly.clear(); + return; + } + emit(ia, ib, ic); + poly.erase(poly.begin() + long(i)); + found = true; + } + if (!found) { + // Clip the smallest convex corner anyway, or failing that the flattest: a thin triangle + // keeps every edge accounted for, and a convex one at least stays on the polygon's side. + size_t pick = 0; + float best = std::numeric_limits::max(); + bool convex_found = false; + for (size_t i = 0; i < n; ++i) { + const Vec2f &a = pts[size_t(poly[(i + n - 1) % n])], &b = pts[size_t(poly[i])], &c = pts[size_t(poly[(i + 1) % n])]; + const float m = cross2(a, b, c); + if (m > flat && (!convex_found || m < best)) { best = m; pick = i; convex_found = true; } + else if (!convex_found && std::abs(m) < best) { best = std::abs(m); pick = i; } + } + emit(poly[(pick + n - 1) % n], poly[pick], poly[(pick + 1) % n]); + poly.erase(poly.begin() + long(pick)); + } + } + if (poly.size() == 3) + emit(poly[0], poly[1], poly[2]); +} + +// Douglas-Peucker on an open polyline; endpoints always survive. +void simplify_polyline(const std::vector &in, float tol, std::vector &out) +{ + out.clear(); + if (in.size() < 3) { out = in; return; } + std::vector keep(in.size(), 0); + keep.front() = keep.back() = 1; + std::vector> stack = { { 0, in.size() - 1 } }; + while (!stack.empty()) { + const auto [i0, i1] = stack.back(); + stack.pop_back(); + if (i1 <= i0 + 1) + continue; + const Vec2f d = in[i1] - in[i0]; + const float len = d.norm(); + float best = -1.f; + size_t bi = i0; + for (size_t i = i0 + 1; i < i1; ++i) { + const Vec2f r = in[i] - in[i0]; + const float dist = (len > 1e-12f) ? std::abs(d.x() * r.y() - d.y() * r.x()) / len : r.norm(); + if (dist > best) { best = dist; bi = i; } + } + if (best > tol) { + keep[bi] = 1; + stack.push_back({ i0, bi }); + stack.push_back({ bi, i1 }); + } + } + for (size_t i = 0; i < in.size(); ++i) + if (keep[i]) out.push_back(in[i]); +} + +} // namespace + +namespace { + +// The cutter's view of one crossing of the mid-level contour with a mesh edge. +struct Crossing +{ + Vec3f p, n; + float t; // parameter from the edge's lower vertex id toward the higher + int a = -1, b = -1; // the edge, lower vertex id first + float sin_min = 1.f; // sine of the smallest angle a contour makes with the edge + float dt_max = std::numeric_limits::max(); // room along the edge before a copy could + // pass the copies of the contour's next vertex + bool used = false; // some cut triangle's contour ends here + bool no_seam = false; // a contour joins it straight to its neighbour on the edge: nothing between + bool sharp = false; + float thick = 0.f; // narrower side of the feature at the crossing, along the gradient + int single = -1, lo = -1, hi = -1; + int near_a = -1, near_b = -1; // the copy on a's side and on b's, when doubled +}; +struct Edge +{ + std::vector xs; // sorted by t +}; + +inline uint64_t step_edge_key(int a, int b) +{ + if (a > b) std::swap(a, b); + return (uint64_t(uint32_t(a)) << 32) | uint32_t(b); +} + +// What cut_mesh_at_steps() learns about the mesh before it cuts anything: the vertex heights (and the +// nudged vertex positions), the crossings on every painted edge, and the verdict `ok` - whether this +// is a step texture worth cutting at all. The verdict is what the prepare path probes on the coarse +// mesh before it chooses how to refine, so refinement and cut agree on which textures are stepped. +struct StepScan +{ + bool ok = false; + std::vector pos, vnormal; + std::vector vh; + float iso = 0.f, range = 0.f; + std::unordered_map edges; + + std::pair sample_between(int a, int b, float t) const + { + const Vec3f p = pos[size_t(a)] + (pos[size_t(b)] - pos[size_t(a)]) * t; + Vec3f n = vnormal[size_t(a)] + (vnormal[size_t(b)] - vnormal[size_t(a)]) * t; + const float l = n.norm(); + n = (l > 1e-12f) ? Vec3f(n / l) : vnormal[size_t(a)]; + return std::make_pair(p, n); + } + bool side_of(float h) const { return h > iso; } +}; + +StepScan scan_steps(const indexed_triangle_set &mesh, const std::vector ®ion, + const HeightFieldSampler &sampler, float step_width_mm, float seam_gap_mm, + float min_feature_mm, bool nudge) +{ + StepScan scan; + const size_t nv = mesh.vertices.size(), nt = mesh.indices.size(); + if (!sampler || region.size() != nt || step_width_mm <= 0.f) + return scan; + + // 1. Heights at the vertices of the cuttable triangles, along the same normals the bake will use. + // Vertex positions are worked on in `pos`: a vertex that happens to sit within a texel of a step + // is nudged off it below, and everything from the edge march to the output uses the nudged place. + std::vector &pos = scan.pos; + std::vector &vnormal = scan.vnormal; + pos = mesh.vertices; + vnormal = texture_displacement_vertex_normals(mesh); + std::vector wanted(nv, 0); + for (size_t t = 0; t < nt; ++t) + if (region[t] != 0) + for (int k = 0; k < 3; ++k) + wanted[size_t(mesh.indices[t][k])] = 1; + std::vector &vh = scan.vh; + vh.assign(nv, 0.f); + tbb::parallel_for(tbb::blocked_range(0, nv), [&](const tbb::blocked_range &r) { + for (size_t v = r.begin(); v < r.end(); ++v) + if (wanted[v]) + vh[v] = sampler(pos[v], vnormal[v]); + }); + float hmin = std::numeric_limits::max(), hmax = -hmin; + for (size_t v = 0; v < nv; ++v) + if (wanted[v]) { hmin = std::min(hmin, vh[v]); hmax = std::max(hmax, vh[v]); } + const float range = scan.range = hmax - hmin; + if (!(range > 1e-6f)) + return scan; + const float iso = scan.iso = 0.5f * (hmin + hmax); + // Only a texture that is mostly steps is cut: one whose heights sit at two levels with little in + // between (a grid, a knurl, a logo). On a smooth or noisy relief the mid-level contour runs through + // every triangle without being a step anywhere, and cutting along it would only multiply the + // triangles; the chord-based refinement is the right tool there. Judged from the vertex heights, + // which sample the whole painted surface. + { + size_t near_level = 0, total = 0; + for (size_t v = 0; v < nv; ++v) + if (wanted[v]) { ++total; near_level += std::abs(vh[v] - iso) > 0.35f * range; } + if (total == 0 || double(near_level) < 0.6 * double(total)) + return scan; + } + + // 1b. A vertex inside the blend of a sharp step - the texture's bilinear ramp is a texel wide - would + // be displaced to a height between the two sides, and every triangle at it would ramp. Such a vertex + // is nudged along the surface away from the step, down its own side's slope, until it samples a + // pure value; a step's width or so. Only where the surface is flat around it: a vertex on a crease + // of the model (a box edge) stays where it is. Skipped by a probe, which only wants the verdict. + if (nudge) { + std::vector creased(nv, 0); + std::vector first_normal(nv, Vec3f::Zero()); + std::vector shortest(nv, std::numeric_limits::max()); // shortest edge at the vertex + for (size_t t = 0; t < nt; ++t) { + const auto &tri = mesh.indices[t]; + const Vec3f &a = pos[size_t(tri[0])]; + Vec3f fn = (pos[size_t(tri[1])] - a).cross(pos[size_t(tri[2])] - a); + const float fl = fn.norm(); + if (fl < 1e-12f) continue; + fn /= fl; + for (int k = 0; k < 3; ++k) { + Vec3f &f0 = first_normal[size_t(tri[k])]; + if (f0.isZero()) f0 = fn; + else if (f0.dot(fn) < 0.985f) creased[size_t(tri[k])] = 1; // ~10 degrees + const float el = (pos[size_t(tri[(k + 1) % 3])] - pos[size_t(tri[k])]).norm(); + shortest[size_t(tri[k])] = std::min(shortest[size_t(tri[k])], el); + shortest[size_t(tri[(k + 1) % 3])] = std::min(shortest[size_t(tri[(k + 1) % 3])], el); + } + } + tbb::parallel_for(tbb::blocked_range(0, nv), [&](const tbb::blocked_range &r) { + for (size_t v = r.begin(); v < r.end(); ++v) { + if (!wanted[v] || creased[v] || std::abs(vh[v] - iso) > 0.47f * range) + continue; + const Vec3f &n = vnormal[v]; + const Vec3f ax = (std::abs(n.x()) < 0.9f) ? Vec3f::UnitX() : Vec3f::UnitY(); + const Vec3f t1 = n.cross(ax).normalized(), t2 = n.cross(t1).normalized(); + const float d = step_width_mm; + const float gx = sampler(pos[v] + t1 * d, n) - sampler(pos[v] - t1 * d, n); + const float gy = sampler(pos[v] + t2 * d, n) - sampler(pos[v] - t2 * d, n); + Vec3f g = t1 * gx + t2 * gy; + const float gl = g.norm(); + if (gl < 0.25f * range) + continue; // not at a step: a slope, which the chord test handles + g /= gl; + // Down the slope for a low vertex, up it for a high one, in steps of half a texel, and + // never by more than a fraction of the shortest edge at the vertex: on a finely refined + // mesh a texel-sized move would fold the triangles around it. + const Vec3f dir = (vh[v] > iso) ? g : Vec3f(-g); + Vec3f best = pos[v]; + float bh = vh[v]; + const float limit = std::min(2.f * step_width_mm, 0.3f * shortest[v]); + for (int i = 1; i <= 8; ++i) { + const float dist = 0.25f * step_width_mm * float(i); + if (dist > limit) + break; + const Vec3f q = pos[v] + dir * dist; + const float h = sampler(q, n); + if ((h > iso) != (vh[v] > iso)) + break; // crossed the step: the nudge would change the vertex's side + best = q; bh = h; + if (std::abs(h - iso) > 0.48f * range) + break; + } + if (std::abs(bh - iso) > std::abs(vh[v] - iso)) { pos[v] = best; vh[v] = bh; } + } + }); + } + + const auto sample_between = [&scan](int a, int b, float t) { return scan.sample_between(a, b, t); }; + const auto side_of = [&scan](float h) { return scan.side_of(h); }; + + // 2. Crossings on the edges, found on the field itself - on a 1 mm edge over a texture with 30 um + // texels, interpolating the endpoint heights could put a step anywhere along the edge. Both + // triangles on an edge see the same crossings, which is what keeps the cut conformal. + auto &edges = scan.edges; + edges.reserve(nt); + + const auto find_crossings = [&](int a, int b, Edge &edge) { + if (a > b) std::swap(a, b); + const float len = (pos[size_t(b)] - pos[size_t(a)]).norm(); + const int steps = std::clamp(int(std::ceil(len / (0.5f * step_width_mm))), 4, 128); + float t0 = 0.f, h0 = vh[size_t(a)] - iso; + for (int i = 1; i <= steps; ++i) { + const float ti = float(i) / float(steps); + const auto [pp, nn] = sample_between(a, b, ti); + const float hi_ = sampler(pp, nn) - iso; + if (h0 * hi_ < 0.f || (hi_ == 0.f && h0 != 0.f)) { + float lo_t = t0, hi_t = ti, lo_h = h0; + for (int k = 0; k < 10; ++k) { + const float tm = 0.5f * (lo_t + hi_t); + const auto [pm, nm] = sample_between(a, b, tm); + if (lo_h * (sampler(pm, nm) - iso) <= 0.f) hi_t = tm; else { lo_t = tm; } + } + Crossing c; + c.a = a; c.b = b; + // Never exactly on a vertex: that would make a zero-area triangle of the split. + c.t = std::clamp(0.5f * (lo_t + hi_t), 1e-3f, 1.f - 1e-3f); + std::tie(c.p, c.n) = sample_between(a, b, c.t); + // Sharp means the field changes by at least half its range across `step_width_mm` in + // the steepest tangent direction - a step, not a slope that merely passes mid-level. + const Vec3f ax = (std::abs(c.n.x()) < 0.9f) ? Vec3f::UnitX() : Vec3f::UnitY(); + const Vec3f t1v = c.n.cross(ax).normalized(), t2v = c.n.cross(t1v).normalized(); + const float d = 0.5f * step_width_mm; + const float gx = sampler(c.p + t1v * d, c.n) - sampler(c.p - t1v * d, c.n); + const float gy = sampler(c.p + t2v * d, c.n) - sampler(c.p - t2v * d, c.n); + const float gl = std::sqrt(gx * gx + gy * gy); + c.sharp = gl >= 0.5f * range; + // How wide the feature is on either side of the step, along the gradient: the + // distance until the field returns to mid-level, in steps of half the step width, up + // to `thick_max`. A texture whose features are only a few texels wide has no pure + // interior for a seam copy to land in and is left to refinement (see the gate below). + { + const float thick_max = 4.f * step_width_mm; + c.thick = thick_max; + if (gl > 1e-12f) { + const Vec3f g = (t1v * gx + t2v * gy) / gl; + for (int sgn = -1; sgn <= 1; sgn += 2) { + const bool high = sgn > 0; // uphill along the gradient + for (int i = 1; i <= 8; ++i) { + const float dist = 0.5f * step_width_mm * float(i); + if (dist >= c.thick) + break; + if (side_of(sampler(c.p + g * (float(sgn) * dist), c.n)) != high) { + c.thick = std::min(c.thick, dist); + break; + } + } + } + } + } + edge.xs.push_back(c); + } + t0 = ti; h0 = hi_; + } + // A pair of crossings closer than the smallest feature the mesh is meant to carry is a feature + // too thin to print, and one closer than the seam gap has no room for two seams: dropping both + // leaves the surface flat there, at the side around it. + const float thinnest = std::max(min_feature_mm, seam_gap_mm); + for (size_t i = 0; i + 1 < edge.xs.size();) { + if ((edge.xs[i + 1].t - edge.xs[i].t) * len < thinnest) + edge.xs.erase(edge.xs.begin() + long(i), edge.xs.begin() + long(i) + 2); + else + ++i; + } + // Parity: the endpoints' sides say whether the count must be odd or even; the march starts and + // ends on the vertex heights themselves, so a mismatch is a numerical accident. Drop the lot + // then - a triangle split on a wrong count is worse than a ramp there. + const bool odd_expected = side_of(vh[size_t(a)]) != side_of(vh[size_t(b)]); + if ((edge.xs.size() % 2 == 1) != odd_expected) + edge.xs.clear(); + }; + + // Edges to examine, in a fixed order so the crossings can be found in parallel. + std::vector edge_keys; + for (size_t t = 0; t < nt; ++t) { + if (region[t] == 0) + continue; + const auto &tri = mesh.indices[t]; + for (int e = 0; e < 3; ++e) + edge_keys.push_back(step_edge_key(tri[e], tri[(e + 1) % 3])); + } + std::sort(edge_keys.begin(), edge_keys.end()); + edge_keys.erase(std::unique(edge_keys.begin(), edge_keys.end()), edge_keys.end()); + std::vector found(edge_keys.size()); + tbb::parallel_for(tbb::blocked_range(0, edge_keys.size()), [&](const tbb::blocked_range &r) { + for (size_t i = r.begin(); i < r.end(); ++i) + find_crossings(int(edge_keys[i] >> 32), int(uint32_t(edge_keys[i])), found[i]); + }); + size_t crossings = 0, sharp = 0; + double thick_sum = 0.; + for (size_t i = 0; i < edge_keys.size(); ++i) { + for (const Crossing &c : found[i].xs) { ++crossings; sharp += c.sharp; thick_sum += c.thick; } + edges.emplace(edge_keys[i], std::move(found[i])); + } + // The mid-level contour has to be a step nearly everywhere it is crossed, or this is not a step + // texture: a noisy relief that happens to sit at two levels is crossed all over, mostly gently. + if (crossings == 0 || double(sharp) < 0.8 * double(crossings)) + return scan; + // And the features have to be wider than the step itself. A texture whose blobs are only a few + // texels across (noise, a fine grain) is bimodal and sharp at every crossing, yet has no pure + // interior for the seam copies to land in, and its contour runs through every triangle: cutting it + // doubles the triangle count for walls the same size as the interpolation blur. Judged from the + // mean feature thickness at the crossings, measured in step widths (2.25 = four to five texels). + if (thick_sum < 2.25 * double(step_width_mm) * double(crossings)) + return scan; + + scan.ok = true; + return scan; + +} + +} // namespace + +bool texture_has_steps_to_cut(const indexed_triangle_set &mesh, const std::vector ®ion, + const HeightFieldSampler &sampler, float step_width_mm, float seam_gap_mm, + float min_feature_mm) +{ + return scan_steps(mesh, region, sampler, step_width_mm, seam_gap_mm, min_feature_mm, /*nudge*/ false).ok; +} + +indexed_triangle_set cut_mesh_at_steps(const indexed_triangle_set &mesh, const std::vector ®ion, + const HeightFieldSampler &sampler, float step_width_mm, + float seam_gap_mm, float min_feature_mm, std::vector *out_source, + size_t *out_cut_count) +{ + const size_t nv = mesh.vertices.size(), nt = mesh.indices.size(); + if (out_cut_count) + *out_cut_count = 0; + const auto passthrough = [&]() { + if (out_source) { + out_source->resize(nt); + std::iota(out_source->begin(), out_source->end(), 0); + } + return mesh; + }; + if (!sampler || region.size() != nt || step_width_mm <= 0.f) + return passthrough(); + + StepScan scan = scan_steps(mesh, region, sampler, step_width_mm, seam_gap_mm, min_feature_mm, /*nudge*/ true); + if (!scan.ok) + return passthrough(); + std::vector &pos = scan.pos; + const std::vector &vnormal = scan.vnormal; + std::vector &vh = scan.vh; + const float iso = scan.iso, range = scan.range; + auto &edges = scan.edges; + const auto sample_between = [&scan](int a, int b, float t) { return scan.sample_between(a, b, t); }; + const auto side_of = [&scan](float h) { return scan.side_of(h); }; + const auto edge_key = [](int a, int b) { return step_edge_key(a, b); }; + (void) vnormal; (void) range; (void) iso; + + // 3. Per triangle: the perimeter as a loop of corners and crossings, and the contour inside the + // triangle traced from a local raster of the field (marching squares at half the step width), as + // polylines from one crossing to another. The trace is what makes a corner of the pattern come out + // as a corner instead of a chord clipping it, and what pairs the crossings up - no guessing. A + // polyline is kept only when its two ends land on distinct crossings the edge march found; + // otherwise the triangle is split at its crossings without a seam (the neighbours' cuts still meet + // no T-junction) and the surface ramps there. Unpainted neighbours of a crossed edge are split the + // same way. + struct Chain + { + int i = -1, j = -1; // crossing indices in the loop, from i to j + std::vector pts; // interior vertices, from i toward j + std::vector lo, hi; // output copies of the interior vertices (same when single) + float clearance = std::numeric_limits::max(); // distance to the nearest other contour + }; + struct Loop + { + std::vector ids; // perimeter entries: corner vertex ids (>= 0) or -(1 + crossing index) + std::vector xs; // crossings in perimeter order + std::vector xpos; // their positions in ids + std::vector chains; + std::vector partner; // per crossing: the crossing its chain leads to + std::vector chain_of; + std::vector no_seam; // crossings joined to their neighbour on the edge with nothing between + Vec3f N = Vec3f::Zero(); + bool ok = false; + }; + const auto walk_edge = [&](int a, int b, std::vector &out) { + auto it = edges.find(edge_key(a, b)); + if (it == edges.end()) return; + Edge &edge = it->second; + if (a < b) for (auto &c : edge.xs) out.push_back(&c); + else for (auto ci = edge.xs.rbegin(); ci != edge.xs.rend(); ++ci) out.push_back(&*ci); + }; + std::vector loops(nt); + for (size_t t = 0; t < nt; ++t) { + const auto &tri = mesh.indices[t]; + Loop &L = loops[t]; + for (int e = 0; e < 3; ++e) { + L.ids.push_back(tri[e]); + std::vector on_edge; + walk_edge(tri[e], tri[(e + 1) % 3], on_edge); + for (Crossing *c : on_edge) { + L.xpos.push_back(int(L.ids.size())); + L.ids.push_back(-(1 + int(L.xs.size()))); + L.xs.push_back(c); + } + } + const Vec3f &a0 = pos[size_t(tri[0])]; + L.N = (pos[size_t(tri[1])] - a0).cross(pos[size_t(tri[2])] - a0); + const float nl = L.N.norm(); + if (nl > 1e-12f) L.N /= nl; + } + + const float cell = 0.5f * step_width_mm; + const auto trace = [&](size_t t) { + Loop &L = loops[t]; + const int n = int(L.xs.size()); + const auto &tri = mesh.indices[t]; + if (n == 0 || region[t] == 0) + return; + if (n % 2 != 0 || L.N.squaredNorm() < 0.5f) + return; + // Local frame on the triangle's plane; the loop runs counter-clockwise in it. + const Vec3f &A = pos[size_t(tri[0])], &B = pos[size_t(tri[1])], &C = pos[size_t(tri[2])]; + const Vec3f U = (B - A).normalized(), V = L.N.cross(U); + const auto to2 = [&](const Vec3f &p) { return Vec2f((p - A).dot(U), (p - A).dot(V)); }; + const auto to3 = [&](const Vec2f &q) { return Vec3f(A + U * q.x() + V * q.y()); }; + const Vec2f a2 = to2(A), b2 = to2(B), c2 = to2(C); + const float det = (b2.x() - a2.x()) * (c2.y() - a2.y()) - (c2.x() - a2.x()) * (b2.y() - a2.y()); + if (std::abs(det) < 1e-12f) + return; + + // A barycentric lattice over the triangle, dense enough that its boundary rows are at least as + // fine as the edge march. Its boundary points take their side from the march - the lower + // corner's side, flipped at every crossing passed - so the contour traced through the lattice + // ends exactly at the crossings the neighbours share; no clipping, no matching. A conflict + // between that and the far corner's own side (an edge whose crossings the march dropped for + // parity) leaves the triangle uncut. + const float longest = std::sqrt(std::max({ (b2 - a2).squaredNorm(), (c2 - b2).squaredNorm(), (a2 - c2).squaredNorm() })); + const int N = std::clamp(int(std::ceil(longest / cell)), 4, 96); + const float h = longest / float(N); + const int W = N + 1; + // Lattice point (j, k): A + (B - A) j/N + (C - A) k/N, j + k <= N. Edge 0 (A->B) is k == 0, + // edge 1 (B->C) is j + k == N, edge 2 (C->A) is j == 0. + const auto at2 = [&](int j, int k) { return Vec2f(a2 + (b2 - a2) * (float(j) / float(N)) + (c2 - a2) * (float(k) / float(N))); }; + std::vector f(size_t(W) * size_t(W), 0.f); + std::vector sign(size_t(W) * size_t(W), 0); + const auto id = [&](int j, int k) { return size_t(k) * size_t(W) + size_t(j); }; + for (int k = 0; k <= N; ++k) + for (int j = 0; j + k <= N; ++j) { + const float wb = float(j) / float(N), wc = float(k) / float(N), wa = 1.f - wb - wc; + Vec3f nn = wa * vnormal[size_t(tri[0])] + wb * vnormal[size_t(tri[1])] + wc * vnormal[size_t(tri[2])]; + const float nl = nn.norm(); + nn = (nl > 1e-12f) ? Vec3f(nn / nl) : L.N; + const float v = sampler(A + (B - A) * wb + (C - A) * wc, nn) - iso; + f[id(j, k)] = v; + sign[id(j, k)] = v > 0.f ? 1 : -1; + } + // Boundary sides from the march. Each edge walks from its first corner in loop order; the + // crossings on it are in that order in L.xs. Position along edge e of lattice point m/N. + const auto edge_point = [&](int e, int m) -> std::pair { // (j, k) + return e == 0 ? std::make_pair(m, 0) : e == 1 ? std::make_pair(N - m, m) : std::make_pair(0, N - m); + }; + // The crossings of edge e, with their parameter from the edge's first corner in loop order. + std::vector>> edge_x(3); + { + int e = -1, k = 0; + for (int idv : L.ids) { + if (idv >= 0) { ++e; continue; } + const Crossing &c = *L.xs[size_t(k)]; + // c.t runs from the lower vertex id; the loop walks tri[e] -> tri[(e + 1) % 3]. + const float u = (tri[e] < tri[(e + 1) % 3]) ? c.t : 1.f - c.t; + edge_x[size_t(e)].push_back({ u, k }); + ++k; + } + } + for (int e = 0; e < 3; ++e) { + int8_t s = side_of(vh[size_t(tri[e])]) ? 1 : -1; + size_t next = 0; + for (int m = 0; m <= N; ++m) { + const float u = float(m) / float(N); + while (next < edge_x[size_t(e)].size() && edge_x[size_t(e)][next].first < u) { s = -s; ++next; } + const auto [j, k] = edge_point(e, m); + sign[id(j, k)] = s; + // Keep the value on the march's side of zero, so interpolation toward the interior + // does not put a crossing where the march has none - and clear of zero, so a contour + // that does pass between this point and the next row is not put on the boundary. + if ((f[id(j, k)] > 0.f) != (s > 0)) + f[id(j, k)] = float(s) * std::max(std::abs(f[id(j, k)]), 0.05f * range); + } + if ((s > 0) != side_of(vh[size_t(tri[(e + 1) % 3])])) + return; // parity conflict on this edge + } + + // Marching triangles. A contour point lives on a lattice edge, keyed by its lower point and + // direction: 0 (j,k)-(j+1,k), 1 (j,k)-(j,k+1), 2 (j+1,k)-(j,k+1). A boundary lattice edge with + // a sign change carries exactly the march crossing in its interval, and is that crossing. + std::vector pts; + std::vector pt_x; // crossing index, or -1 for an interior point + std::unordered_map point_of; + std::vector> link; + bool bad = false; + const auto point_on = [&](int j, int k, int dir) -> int { + const uint32_t key = (uint32_t(j) << 16) | (uint32_t(k) << 2) | uint32_t(dir); + auto it = point_of.find(key); + if (it != point_of.end()) + return it->second; + // Endpoints (j0, k0) and (jb, kb). + const int j0 = dir == 2 ? j + 1 : j, k0 = k; + const int jb = dir == 0 ? j + 1 : j, kb = dir == 0 ? k : k + 1; + int xk = -1; + Vec2f p; + int e = -1; + float u0 = 0.f, u1 = 0.f; + if (k0 == 0 && kb == 0) { e = 0; u0 = float(j0) / float(N); u1 = float(jb) / float(N); } + else if (j0 + k0 == N && jb + kb == N) { e = 1; u0 = float(k0) / float(N); u1 = float(kb) / float(N); } + else if (j0 == 0 && jb == 0) { e = 2; u0 = float(N - k0) / float(N); u1 = float(N - kb) / float(N); } + if (e >= 0) { + const float lo = std::min(u0, u1), hi = std::max(u0, u1); + for (const auto &[u, k_] : edge_x[size_t(e)]) + if (u >= lo && u < hi) { if (xk >= 0) bad = true; xk = k_; } + if (xk < 0) { bad = true; p = 0.5f * (at2(j0, k0) + at2(jb, kb)); } + else p = to2(L.xs[size_t(xk)]->p); + } else { + const float fa = f[id(j0, k0)], fb = f[id(jb, kb)]; + const float s = (fa != fb) ? std::clamp(fa / (fa - fb), 0.f, 1.f) : 0.5f; + p = at2(j0, k0) + (at2(jb, kb) - at2(j0, k0)) * s; + } + const int pid = int(pts.size()); + pts.push_back(p); + pt_x.push_back(xk); + link.push_back({ -1, -1 }); + point_of.emplace(key, pid); + return pid; + }; + const auto connect = [&](int p, int q) { + for (int *slot : { &link[size_t(p)][0], &link[size_t(p)][1] }) + if (*slot < 0) { *slot = q; break; } + for (int *slot : { &link[size_t(q)][0], &link[size_t(q)][1] }) + if (*slot < 0) { *slot = p; break; } + }; + // Up triangle (j,k),(j+1,k),(j,k+1): edges dir 0 at (j,k), dir 1 at (j,k), dir 2 at (j,k). + // Down triangle (j+1,k),(j+1,k+1),(j,k+1): dir 1 at (j+1,k), dir 0 at (j,k+1), dir 2 at (j,k). + for (int k = 0; k < N; ++k) + for (int j = 0; j + k < N; ++j) { + { + const int8_t s0_ = sign[id(j, k)], s1_ = sign[id(j + 1, k)], s2_ = sign[id(j, k + 1)]; + if (!(s0_ == s1_ && s1_ == s2_)) { + std::vector ends; + if (s0_ != s1_) ends.push_back(point_on(j, k, 0)); + if (s0_ != s2_) ends.push_back(point_on(j, k, 1)); + if (s1_ != s2_) ends.push_back(point_on(j, k, 2)); + if (ends.size() == 2) connect(ends[0], ends[1]); + } + } + if (j + k + 1 < N) { + const int8_t s0_ = sign[id(j + 1, k)], s1_ = sign[id(j + 1, k + 1)], s2_ = sign[id(j, k + 1)]; + if (!(s0_ == s1_ && s1_ == s2_)) { + std::vector ends; + if (s0_ != s1_) ends.push_back(point_on(j + 1, k, 1)); + if (s1_ != s2_) ends.push_back(point_on(j, k + 1, 0)); + if (s0_ != s2_) ends.push_back(point_on(j, k, 2)); + if (ends.size() == 2) connect(ends[0], ends[1]); + } + } + } + if (bad) + return; + + // Chains from every crossing to the crossing it leads to. Closed loops never touch the + // boundary and are left alone - refinement resolves a feature that sits entirely inside a + // triangle. + std::vector matched(static_cast(n), -1); + std::vector chains; + std::vector> raws; // the traced polyline of each chain, crossing to crossing + std::vector seen(pts.size(), 0); + for (size_t s = 0; s < pts.size(); ++s) { + if (seen[s] || pt_x[s] < 0) + continue; + std::vector full; + int prev = -1, cur = int(s); + while (cur >= 0 && !seen[size_t(cur)]) { + seen[size_t(cur)] = 1; + full.push_back(pts[size_t(cur)]); + if (pt_x[size_t(cur)] >= 0 && cur != int(s)) + break; + const int nx0 = link[size_t(cur)][0], nx1 = link[size_t(cur)][1]; + const int next = (nx0 != prev) ? nx0 : nx1; + prev = cur; cur = next; + } + if (cur < 0 || pt_x[size_t(cur)] < 0 || cur == int(s)) + return; // a contour that starts at a crossing and ends nowhere: the lattice is inconsistent + Chain ch; + ch.i = pt_x[s]; + ch.j = pt_x[size_t(cur)]; + if (matched[size_t(ch.i)] >= 0 || matched[size_t(ch.j)] >= 0 || ch.i == ch.j) + return; + matched[size_t(ch.i)] = ch.j; + matched[size_t(ch.j)] = ch.i; + chains.push_back(std::move(ch)); + raws.push_back(std::move(full)); + } + // How close each contour comes to another: two sides of a thin feature. Simplifying either by + // the usual tolerance, or offsetting its copies by the usual gap, could then push it across + // the other, so both are scaled to the clearance. Segment-to-segment over every pair. + { + const auto seg_dist = [](const Vec2f &a, const Vec2f &b, const Vec2f &c, const Vec2f &d) { + const auto pt_seg = [](const Vec2f &q, const Vec2f &u, const Vec2f &v) { + const Vec2f uv = v - u; + const float t = std::clamp((q - u).dot(uv) / std::max(uv.squaredNorm(), 1e-12f), 0.f, 1.f); + return (q - (u + uv * t)).norm(); + }; + return std::min({ pt_seg(a, c, d), pt_seg(b, c, d), pt_seg(c, a, b), pt_seg(d, a, b) }); + }; + for (size_t x = 0; x < raws.size(); ++x) + for (size_t y = x + 1; y < raws.size(); ++y) { + float best = std::min(chains[x].clearance, chains[y].clearance); + for (size_t i = 0; i + 1 < raws[x].size(); ++i) + for (size_t j = 0; j + 1 < raws[y].size(); ++j) + best = std::min(best, seg_dist(raws[x][i], raws[x][i + 1], raws[y][j], raws[y][j + 1])); + chains[x].clearance = std::min(chains[x].clearance, best); + chains[y].clearance = std::min(chains[y].clearance, best); + } + } + for (size_t c = 0; c < chains.size(); ++c) { + Chain &ch = chains[c]; + const std::vector &full = raws[c]; + const float tol = std::min(h, 0.25f * ch.clearance); + const float gap = std::min(seam_gap_mm, 0.5f * ch.clearance); + std::vector simple; + simplify_polyline(full, tol, simple); + // A vertex within the seam gap of the previous one, or of the far end, leaves its copies no + // room to clear the blend of the step, and a corner that close to an edge is within the + // seam's own width anyway: it is merged into its neighbour. + { + std::vector spaced; + spaced.push_back(simple.front()); + for (size_t m = 1; m + 1 < simple.size(); ++m) + if ((simple[m] - spaced.back()).norm() >= gap && (simple[m] - simple.back()).norm() >= gap) + spaced.push_back(simple[m]); + spaced.push_back(simple.back()); + simple = std::move(spaced); + } + const int m_loop = int(L.ids.size()); + const bool adjacent = (L.xpos[size_t(ch.j)] - L.xpos[size_t(ch.i)] + m_loop) % m_loop == 1 || + (L.xpos[size_t(ch.i)] - L.xpos[size_t(ch.j)] + m_loop) % m_loop == 1; + if (simple.size() == 2 && full.size() > 2 && adjacent) { + // Never down to a bare chord between two crossings next to each other on one edge: the + // shallow pocket between them would be joined along the edge itself, and the regions on + // either side would overlap. The raw point farthest from the chord stays. + const Vec2f d = full.back() - full.front(); + const float dl = std::max(d.norm(), 1e-9f); + size_t bi = full.size() / 2; + float bd = 0.1f * tol; // below this the contour is straight: take its middle point + for (size_t m = 1; m + 1 < full.size(); ++m) { + const Vec2f r = full[m] - full.front(); + const float dist = std::abs(d.x() * r.y() - d.y() * r.x()) / dl; + if (dist > bd) { bd = dist; bi = m; } + } + simple.insert(simple.begin() + 1, full[bi]); + } + // Hairpins - the contour doubling back on itself with less than the seam gap between the + // two sides, the tip of a grain line thinner than the gap - would give a region polygon + // whose two sides cross. The spike vertex goes, and the feature is flattened there, as a + // thin feature on an edge is. + for (bool again = true; again;) { + again = false; + for (size_t m = 1; m + 1 < simple.size(); ++m) { + const Vec2f d1 = simple[m] - simple[m - 1], d2 = simple[m + 1] - simple[m]; + const float l1 = d1.norm(), l2 = d2.norm(); + if (l1 < 1e-9f || l2 < 1e-9f) { simple.erase(simple.begin() + long(m)); again = true; break; } + if (d1.dot(d2) / (l1 * l2) > -0.5f) + continue; // turns less than 120 degrees: a corner, not a hairpin + // Distance between the two sides: the shorter side's far end to the other side. + const Vec2f &tip = simple[m]; + const auto dist_to = [&](const Vec2f &q, const Vec2f &a, const Vec2f &b) { + const Vec2f ab = b - a; + const float t = std::clamp((q - a).dot(ab) / std::max(ab.squaredNorm(), 1e-12f), 0.f, 1.f); + return (q - (a + ab * t)).norm(); + }; + const float between = (l1 < l2) ? dist_to(simple[m - 1], tip, simple[m + 1]) : dist_to(simple[m + 1], simple[m - 1], tip); + if (between < seam_gap_mm) { + simple.erase(simple.begin() + long(m)); + again = true; + break; + } + } + } + for (size_t m = 1; m + 1 < simple.size(); ++m) + ch.pts.push_back(to3(simple[m])); + if (ch.pts.empty() && adjacent) { + // Two crossings next to each other on one edge and the contour, once simplified, going + // straight from one to the other: an empty pocket. There is nothing to wall off - a wall + // along the edge would be a zero-area strip, and the neighbour would lay its own on top + // - so neither crossing gets a seam. Recorded here and applied once the tracing is done. + L.no_seam.push_back(ch.i); + L.no_seam.push_back(ch.j); + } + } + for (int k = 0; k < n; ++k) + if (matched[size_t(k)] < 0) + return; + L.chains = std::move(chains); + L.partner = std::move(matched); + L.chain_of.assign(static_cast(n), -1); + for (size_t c = 0; c < L.chains.size(); ++c) { + L.chain_of[size_t(L.chains[c].i)] = int(c); + L.chain_of[size_t(L.chains[c].j)] = int(c); + } + L.ok = true; + }; + tbb::parallel_for(tbb::blocked_range(0, nt), [&](const tbb::blocked_range &r) { + for (size_t t = r.begin(); t < r.end(); ++t) + trace(t); + }); + size_t cut_count = 0; + for (size_t t = 0; t < nt; ++t) { + cut_count += loops[t].ok ? 1 : 0; + if (loops[t].ok) + for (int k : loops[t].no_seam) + loops[t].xs[size_t(k)]->no_seam = true; + } + + // 4. Which crossings a contour reaches, and how obliquely: the copies of a doubled crossing are moved + // apart along the edge, so the shallower the contour meets the edge the further apart they need to + // be to sit the same distance clear of it. + for (size_t t = 0; t < nt; ++t) { + Loop &L = loops[t]; + if (!L.ok) + continue; + for (auto &ch : L.chains) { + Crossing &p = *L.xs[size_t(ch.i)], &q = *L.xs[size_t(ch.j)]; + p.used = q.used = true; + const Vec3f p_next = ch.pts.empty() ? q.p : ch.pts.front(); + const Vec3f q_prev = ch.pts.empty() ? p.p : ch.pts.back(); + for (auto [c, toward] : { std::make_pair(&p, p_next), std::make_pair(&q, q_prev) }) { + Vec3f d = toward - c->p; + const float dl = d.norm(); + if (dl < 1e-12f) continue; + const Vec3f ed = (pos[size_t(c->b)] - pos[size_t(c->a)]).normalized(); + c->sin_min = std::min(c->sin_min, ed.cross(d / dl).norm()); + // Only the part of a move along the edge that runs along the contour brings the copy + // toward the next vertex's copies. With no interior vertex the chord's two ends share + // it; otherwise the interior copies, placed first, say how much room there is. + if (ch.pts.empty()) + c->dt_max = std::min(c->dt_max, 0.45f * dl / std::max(std::abs(ed.dot(d / dl)), 1e-3f)); + } + } + } + + // 5. Output vertices: the originals, then one or two per crossing and per interior contour vertex. + // A sharp crossing that a contour reaches is doubled, into two vertices on the edge either side of + // it, each on its own side of the contour; a triangle on that edge that is not cut itself simply + // carries both on its perimeter, one after the other, so the cut neighbour's wall still meets a + // closed surface. Keeping the copies on the edge is what keeps every triangle's split planar and + // inside the triangle. A crossing that is not sharp, or that no contour reaches, stays one vertex + // and the surface ramps there. + indexed_triangle_set out; + out.vertices = pos; + const float half_gap = 0.5f * seam_gap_mm; + // Interior contour vertices: doubled when the contour's ends are, offset across the contour along + // the bisector of the two segments meeting there, toward the high side; mitre limited, and held + // inside the triangle - a copy pushed across an edge would fold the region over the neighbour. + for (size_t t = 0; t < nt; ++t) { + Loop &L = loops[t]; + if (!L.ok) + continue; + const auto &tri = mesh.indices[t]; + const Vec3f &A = pos[size_t(tri[0])], &B = pos[size_t(tri[1])], &C = pos[size_t(tri[2])]; + const Vec3f U = (B - A).normalized(), V = L.N.cross(U); + const Vec2f a2(0.f, 0.f), b2((B - A).dot(U), (B - A).dot(V)), c2((C - A).dot(U), (C - A).dot(V)); + const float det = (b2.x() - a2.x()) * (c2.y() - a2.y()) - (c2.x() - a2.x()) * (b2.y() - a2.y()); + const auto bary2 = [&](const Vec3f &p) { + const Vec2f q((p - A).dot(U), (p - A).dot(V)); + const float wb = ((q.x() - a2.x()) * (c2.y() - a2.y()) - (c2.x() - a2.x()) * (q.y() - a2.y())) / det; + const float wc = ((b2.x() - a2.x()) * (q.y() - a2.y()) - (q.x() - a2.x()) * (b2.y() - a2.y())) / det; + return Vec3f(1.f - wb - wc, wb, wc); + }; + // The largest fraction of the offset that keeps the copy inside, with a little to spare. + const auto inside_fraction = [&](const Vec3f &p, const Vec3f &off) { + const Vec3f w0 = bary2(p), w1 = bary2(p + off); + float lambda = 1.f; + for (int i = 0; i < 3; ++i) + if (w1[i] < 0.f && w0[i] > w1[i]) + lambda = std::min(lambda, w0[i] / (w0[i] - w1[i])); + return 0.9f * std::max(lambda, 0.f); + }; + const bool s0 = side_of(vh[size_t(mesh.indices[t][0])]); + for (auto &ch : L.chains) { + const Crossing &p = *L.xs[size_t(ch.i)], &q = *L.xs[size_t(ch.j)]; + const bool doubled = p.sharp && p.used && !p.no_seam && q.sharp && q.used && !q.no_seam; + // The arc right after crossing ch.i lies on one side; which way that is from the contour + // says where "high" is. The arc's side: the loop's first corner's, flipped once per + // crossing passed, so after crossing index k it is !s0 for even k and s0 for odd. + const bool arc_high = (ch.i % 2 == 0) ? !s0 : s0; + const int after_pos = L.xpos[size_t(ch.i)] + 1; + const int entry = L.ids[size_t(after_pos) % L.ids.size()]; + const Vec3f ref = entry >= 0 ? pos[size_t(entry)] : L.xs[size_t(-entry - 1)]->p; + const size_t m = ch.pts.size(); + ch.lo.resize(m); ch.hi.resize(m); + // The left of the walk from i to j is one side of the contour all along it; whether that + // side is high is settled once, at the first segment, where the arc is a known reference. + // The direction the contour leaves the crossing in, taken over a usable length: the first + // interior vertex may sit right next to it. + Vec3f d_first = q.p - p.p; + for (size_t k = 0; k < m; ++k) + if ((ch.pts[k] - p.p).norm() > 0.25f * seam_gap_mm) { d_first = ch.pts[k] - p.p; break; } + d_first.normalize(); + const bool left_high = (L.N.cross(d_first).dot(ref - p.p) > 0.f) == arc_high; + for (size_t k = 0; k < m; ++k) { + if (!doubled) { + ch.lo[k] = ch.hi[k] = int(out.vertices.size()); + out.vertices.push_back(ch.pts[k]); + continue; + } + const Vec3f prev = (k == 0) ? p.p : ch.pts[k - 1]; + const Vec3f next = (k + 1 == m) ? q.p : ch.pts[k + 1]; + Vec3f d1 = (ch.pts[k] - prev).normalized(), d2 = (next - ch.pts[k]).normalized(); + Vec3f n1 = L.N.cross(d1), n2 = L.N.cross(d2); + Vec3f bis = n1 + n2; + float bl = bis.norm(); + Vec3f off; + if (bl < 1e-6f) off = n1; + else { + bis /= bl; + // Mitre: the offset polyline stays `half_gap` from both segments, up to twice that. + const float cos_half = std::max(bis.dot(n1), 0.5f); + off = bis / cos_half; + } + if (!left_high) off = -off; + // No further along either segment than half its length, so copies never cross their + // neighbours'; across the segments the offset is free. + const float ol = off.norm(); + const float reach = std::min({ half_gap, 0.45f * ch.clearance, + 0.45f * (ch.pts[k] - prev).norm() / std::max(std::abs(off.dot(d1)) / ol, 1e-3f), + 0.45f * (next - ch.pts[k]).norm() / std::max(std::abs(off.dot(d2)) / ol, 1e-3f) }); + const Vec3f off_lo = -off * (reach * inside_fraction(ch.pts[k], -off * reach)); + const Vec3f off_hi = off * (reach * inside_fraction(ch.pts[k], off * reach)); + // As for the crossings: each copy at the first of a few positions outward that samples + // pure, or the vertex stays single. + const auto place = [&](const Vec3f &o, bool high) -> int { + for (int step = 1; step <= 4; ++step) { + const Vec3f q = ch.pts[k] + o * (0.5f * float(step)); + if (step > 2 && inside_fraction(ch.pts[k], q - ch.pts[k]) < 0.999f) break; + const float h = sampler(q, L.N) - iso; + if ((h > 0.f) == high && std::abs(h) >= 0.3f * range) { + out.vertices.push_back(q); + return int(out.vertices.size()) - 1; + } + } + return -1; + }; + const int il = place(off_lo, false), ih = place(off_hi, true); + if (il < 0 || ih < 0) { + out.vertices.resize(out.vertices.size() - (il >= 0 ? 1 : 0) - (ih >= 0 ? 1 : 0)); + ch.lo[k] = ch.hi[k] = int(out.vertices.size()); + out.vertices.push_back(ch.pts[k]); + continue; + } + ch.lo[k] = il; ch.hi[k] = ih; + } + // Room for the end crossings' copies along the edge: the nearer of the first interior + // vertex's copies, measured along the contour, less a margin, over the edge's share of + // that direction. + if (m > 0) + for (auto [c, k, other] : { std::make_tuple(&p, size_t(0), q.p), std::make_tuple(&q, m - 1, p.p) }) { + (void) other; + Crossing &cr = const_cast(*c); + const Vec3f d = ch.pts[k] - cr.p; + const float dl = d.norm(); + if (dl < 1e-9f) continue; + const Vec3f dn = d / dl; + const float u_lo = (out.vertices[size_t(ch.lo[k])] - cr.p).dot(dn); + const float u_hi = (out.vertices[size_t(ch.hi[k])] - cr.p).dot(dn); + const Vec3f ed = (pos[size_t(cr.b)] - pos[size_t(cr.a)]).normalized(); + cr.dt_max = std::min(cr.dt_max, 0.9f * std::max(std::min(u_lo, u_hi), 0.f) / std::max(std::abs(ed.dot(dn)), 1e-3f)); + } + } + } + + for (auto &kv : edges) { + Edge &edge = kv.second; + for (size_t i = 0; i < edge.xs.size(); ++i) { + Crossing &c = edge.xs[i]; + if (!(c.sharp && c.used) || c.no_seam) { + c.single = int(out.vertices.size()); + out.vertices.push_back(c.p); + continue; + } + const float len = (pos[size_t(c.b)] - pos[size_t(c.a)]).norm(); + // Room along the edge: a third of the way to the previous and the next crossing or vertex, + // and less than half the shortest contour segment leaving the crossing, so the copies of + // the contour's next vertex cannot cross these. + const float t_prev = (i == 0) ? 0.f : edge.xs[i - 1].t; + const float t_next = (i + 1 == edge.xs.size()) ? 1.f : edge.xs[i + 1].t; + float dt = half_gap / std::max(c.sin_min, 0.15f) / std::max(len, 1e-6f); + dt = std::min({ dt, (c.t - t_prev) / 3.f, (t_next - c.t) / 3.f, c.dt_max / std::max(len, 1e-6f) }); + if (!(dt > 1e-6f)) { + c.single = int(out.vertices.size()); + out.vertices.push_back(c.p); + continue; + } + // The side of the edge before this crossing: the lower vertex's, flipped once per crossing + // passed on the way. + const bool before_high = side_of(vh[size_t(c.a)]) != (i % 2 == 1); + // Each copy has to sample a pure value of its own side - inside the step's blend it would be + // displaced to a height between the sides, and one such vertex tilts every triangle at it. + // The copy is placed at the first of a few positions along the edge, from dt outward, that + // samples pure; when none does the crossing stays one vertex and the surface ramps across + // that edge only. + // The search never goes past what the neighbours allow: less than half way to the next + // crossing along the edge (so two crossings' copies never meet), most of the way to the + // edge's vertex, and not past dt_max toward the contour's next vertex. + const float lo_share = (i == 0) ? 0.9f : 0.49f, hi_share = (i + 1 == edge.xs.size()) ? 0.9f : 0.49f; + const float t_lo_lim = c.t - std::min(lo_share * (c.t - t_prev), c.dt_max / std::max(len, 1e-6f)); + const float t_hi_lim = c.t + std::min(hi_share * (t_next - c.t), c.dt_max / std::max(len, 1e-6f)); + const auto place = [&](float sgn, bool high) -> int { + for (int step = 1; step <= 6; ++step) { + const float tt = c.t + sgn * dt * (0.5f * float(step)); + if (tt < t_lo_lim || tt > t_hi_lim) break; + const auto [q, qn] = sample_between(c.a, c.b, tt); + const float h = sampler(q, qn) - iso; + if ((h > 0.f) == high && std::abs(h) >= 0.3f * range) { + out.vertices.push_back(q); + return int(out.vertices.size()) - 1; + } + } + return -1; + }; + const int na = place(-1.f, before_high), nb = place(+1.f, !before_high); + if (na < 0 || nb < 0) { + out.vertices.resize(out.vertices.size() - (na >= 0 ? 1 : 0) - (nb >= 0 ? 1 : 0)); + c.single = int(out.vertices.size()); + out.vertices.push_back(c.p); + continue; + } + c.near_a = na; c.near_b = nb; + c.lo = before_high ? c.near_b : c.near_a; + c.hi = before_high ? c.near_a : c.near_b; + } + } + // 6. Output triangles: each region of a cut triangle ear-clipped on its own side's copies, a wall + // strip per contour wound with the surface, everything else passed through or plainly split. + out.indices.reserve(nt * 2); + std::vector source; + source.reserve(nt * 2); + const auto emit = [&](int a, int b, int c, int src) { out.indices.emplace_back(a, b, c); source.push_back(src); }; + const auto copy_for = [&](const Crossing &c, bool high) { return c.single >= 0 ? c.single : (high ? c.hi : c.lo); }; + // Triangulates a triangle that carries extra vertices on its edges without a zero-area sliver: the + // first extra vertex found is joined to the opposite corner, which splits the triangle in two that + // carry fewer extras each, and so on. A fan from a corner cannot do this - the extras on that + // corner's own edges are collinear with it. Each side is the corner followed by the extras on the + // way to the next corner. + const std::function, 3>, int)> split_tri = + [&](std::array, 3> side, int src) { + for (int e = 0; e < 3; ++e) + if (side[size_t(e)].size() > 1) { + const auto &S = side[size_t(e)]; + const int a = S[0], x = S[1], c = side[size_t((e + 2) % 3)][0]; + std::vector rest(S.begin() + 1, S.end()); // x and the extras after it + split_tri({ std::vector{ a }, std::vector{ x }, side[size_t((e + 2) % 3)] }, src); + split_tri({ std::move(rest), side[size_t((e + 1) % 3)], std::vector{ c } }, src); + return; + } + emit(side[0][0], side[1][0], side[2][0], src); + }; + // A wall strip between the low and high copies of a contour, from crossing p to crossing q. + const auto emit_wall = [&](const Crossing &p, const Chain &ch, const Crossing &q, const Vec3f &N, int src) { + std::vector lo, hi; + lo.push_back(copy_for(p, false)); hi.push_back(copy_for(p, true)); + for (size_t k = 0; k < ch.pts.size(); ++k) { lo.push_back(ch.lo[k]); hi.push_back(ch.hi[k]); } + lo.push_back(copy_for(q, false)); hi.push_back(copy_for(q, true)); + auto tri_up = [&](int i, int j, int k) { + if (i == j || j == k || k == i) return; + const Vec3f &A = out.vertices[size_t(i)], &B = out.vertices[size_t(j)], &C = out.vertices[size_t(k)]; + if ((B - A).cross(C - A).dot(N) >= 0.f) emit(i, j, k, src); else emit(i, k, j, src); + }; + for (size_t k = 0; k + 1 < lo.size(); ++k) { + tri_up(lo[k], lo[k + 1], hi[k + 1]); + tri_up(lo[k], hi[k + 1], hi[k]); + } + }; + + for (size_t t = 0; t < nt; ++t) { + const auto &tri = mesh.indices[t]; + const Loop &L = loops[t]; + const int n = int(L.xs.size()); + if (n == 0) { emit(tri[0], tri[1], tri[2], int(t)); continue; } + const size_t m = L.ids.size(); + const bool s0 = side_of(vh[size_t(tri[0])]); + + if (L.ok) { + // Every crossing starts one region: the run of the loop after it, on the side the field has + // there, following any other contour it meets back to the loop, until it returns to the + // start. Each region is tracked per (crossing, side) so the one between two contours is + // emitted once, not once per contour. The region's vertices are collected in the triangle's + // plane for ear clipping, since a traced contour can make it concave. + const Vec3f &A = pos[size_t(tri[0])], &B = pos[size_t(tri[1])]; + const Vec3f U = (B - A).normalized(), V = L.N.cross(U); + std::vector pts2; + std::vector ids2; + const auto push = [&](int id) { + const Vec3f &p = out.vertices[size_t(id)]; + pts2.emplace_back((p - A).dot(U), (p - A).dot(V)); + ids2.push_back(id); + }; + // Interior vertices of chain c, walking from crossing `from`, on side S. + const auto push_chain = [&](int c, int from, bool S) { + const Chain &ch = L.chains[size_t(c)]; + const auto &cp = S ? ch.hi : ch.lo; + if (from == ch.i) for (size_t k = 0; k < cp.size(); ++k) push(cp[k]); + else for (size_t k = cp.size(); k-- > 0;) push(cp[k]); + }; + const auto side_after = [&](int k) { return (k % 2 == 0) ? !s0 : s0; }; + std::vector> consumed(static_cast(n), { false, false }); + for (int k = 0; k < n; ++k) { + const bool S = side_after(k); + if (consumed[size_t(k)][S]) + continue; + pts2.clear(); ids2.clear(); + push(copy_for(*L.xs[size_t(k)], S)); + consumed[size_t(k)][S] = true; + size_t i = size_t(L.xpos[size_t(k)] + 1) % m; + for (size_t guard = 0; guard < 2 * m; ++guard) { + const int id = L.ids[i]; + if (id >= 0) { push(id); i = (i + 1) % m; continue; } + const int c = -id - 1; + consumed[size_t(c)][S] = true; + push(copy_for(*L.xs[size_t(c)], S)); + if (c == L.partner[size_t(k)]) + break; + // Follow the contour from c to its partner, then continue the loop after it. + const int pc = L.partner[size_t(c)]; + push_chain(L.chain_of[size_t(c)], c, S); + consumed[size_t(pc)][S] = true; + push(copy_for(*L.xs[size_t(pc)], S)); + if (pc == L.partner[size_t(k)]) + break; + i = size_t(L.xpos[size_t(pc)] + 1) % m; + } + // Close the polygon along the region's own contour back to k, in reverse. + push_chain(L.chain_of[size_t(k)], L.partner[size_t(k)], S); + std::vector poly(ids2.size()); + std::iota(poly.begin(), poly.end(), 0); + ear_clip(pts2, poly, [&](int a, int b, int c) { emit(ids2[size_t(a)], ids2[size_t(b)], ids2[size_t(c)], int(t)); }); + } + for (const auto &ch : L.chains) + emit_wall(*L.xs[size_t(ch.i)], ch, *L.xs[size_t(ch.j)], L.N, int(t)); + continue; + } + // Not cut: split at the crossings so the neighbours' cuts meet no T-junction. A doubled crossing + // contributes both copies, the one nearer the vertex the walk comes from first. + std::array, 3> side; + int e = -1; + for (size_t i = 0; i < m; ++i) { + if (L.ids[i] >= 0) { ++e; side[size_t(e)].push_back(L.ids[i]); continue; } + const Crossing &c = *L.xs[size_t(-L.ids[i] - 1)]; + if (c.single >= 0) { side[size_t(e)].push_back(c.single); continue; } + const bool from_a = tri[e] == c.a; + side[size_t(e)].push_back(from_a ? c.near_a : c.near_b); + side[size_t(e)].push_back(from_a ? c.near_b : c.near_a); + } + split_tri(std::move(side), int(t)); + } + // Safety net: a texture whose features are of the seam's own scale can fold copies over each + // other in spite of everything above. When more than a trace of the output faces the wrong way, + // the cut is not trusted and the mesh goes out as it came in. + { + size_t inverted = 0; + for (size_t i = 0; i < out.indices.size(); ++i) { + const auto &f = out.indices[i]; + const auto &g = mesh.indices[size_t(source[i])]; + const Vec3f n = (out.vertices[size_t(f[1])] - out.vertices[size_t(f[0])]).cross(out.vertices[size_t(f[2])] - out.vertices[size_t(f[0])]); + const Vec3f ns = (mesh.vertices[size_t(g[1])] - mesh.vertices[size_t(g[0])]).cross(mesh.vertices[size_t(g[2])] - mesh.vertices[size_t(g[0])]); + inverted += n.dot(ns) < 0.f && n.norm() > 1e-7f; + } + if (double(inverted) > 0.002 * double(out.indices.size())) + return passthrough(); + } + if (out_source) *out_source = std::move(source); + if (out_cut_count) *out_cut_count = cut_count; + return out; +} + } // namespace Slic3r diff --git a/src/libslic3r/TextureDisplacement.hpp b/src/libslic3r/TextureDisplacement.hpp index 63cb443420..ea208fbcb3 100644 --- a/src/libslic3r/TextureDisplacement.hpp +++ b/src/libslic3r/TextureDisplacement.hpp @@ -20,6 +20,11 @@ namespace Slic3r { +// Optional step-by-step capture of a bake; see TextureBake/TextureBakeDebug.hpp. Forward declared and +// taken by pointer so this header, which most of the texture feature includes, does not grow a +// dependency for something only the debug view and the benchmarks use. +class BakeStageRecorder; + class ModelVolume; // Bits of subdivide_mesh_adaptive()'s per-triangle `refine_region` mask. See that function. @@ -259,14 +264,14 @@ struct TextureDisplacementLayer // chart ids, so this is meaningful only against the unwrap it was made on. std::vector island_groups; - // Only used by TextureProjectionMethod::LSCM: manual per-vertex UV edits made in the UV editor's - // Vertex/Edge select modes. Each pair is (mesh vertex index, its overriding raw-unwrap coordinate in - // mm) - the *raw* unwrap position, i.e. before the island transform, so the edited vertex still - // moves and rotates with its island. In compute_lscm_uvs() this replaces the automatic unwrap - // coordinate for that vertex; in the editor it edits the displayed geometry directly. Keyed in mesh- - // vertex space like lscm_seam_edges (dropped on a topology change). The raw coordinate is only - // meaningful against the current unwrap, so a re-unwrap clears these. A mesh vertex shared by several - // charts (a seam vertex) settles on one, matching compute_lscm_uvs()'s single-UV-per-vertex rule. + // Only used by TextureProjectionMethod::LSCM: manual UV edits made in the UV editor's Vertex/Edge + // select modes. Each pair is (key, its overriding raw-unwrap coordinate in mm) - the *raw* unwrap + // position, i.e. before the island transform, so the edited vertex still moves and rotates with its + // island. The key names what is edited (see apply_lscm_uv_overrides()): + // - negative: one unwrapped vertex of the current unwrap, as lscm_uv_override_key(index). This is what + // the editor stores, so a seam vertex dragged in one island leaves its copies in the neighbouring + // islands where they are. Like `islands`, meaningful only against the unwrap it was made on. + // - zero or positive: a mesh vertex, every unwrapped copy of it (how edits were stored before). std::vector> lscm_uv_overrides; // How this layer folds into the displacement accumulated by the layers below it. Ignored for @@ -365,21 +370,27 @@ struct TextureDisplacementOptions // deliberately (which is a blunter version of the per-layer edge-smoothing falloff). bool smooth_skip_border = true; - // Alternative bake pipeline, for side-by-side comparison. The path above is topology-preserving - // and needs the mesh prepared first; this one refines, removes slivers, displaces and optionally - // simplifies in one run. Off by default, and it produces no colour - it rebuilds the topology, so - // the per-facet assignment has nothing stable to attach to. - bool pipeline_v2 = false; - float v2_refine_mm = 0.3f; + // Which bake pipeline. On (the default): refine, remove slivers, displace and simplify in one run, + // nothing to prepare first. Off, the classic path: the mesh is prepared first (remesh, adaptive + // subdivision, step cut) and then displaced vertex by vertex, keeping the topology - which is what + // colour needs, since the per-facet assignment has nothing stable to attach to once the topology is + // rebuilt. Projects saved with the classic path keep it: the flag is stored per volume. + bool pipeline_v2 = true; + // Refinement edge length, mm. 0 (the default) means automatic: chosen from the texture's texel + // size and sharpness and the model's size, see recommend_v2_resolution(). A saved project with an + // explicit value keeps it. + float v2_refine_mm = 0.f; bool v2_regularize = false; - int v2_max_triangles_k = 750; // 0 skips simplification, which is worth comparing on its own - // Stop displacement pushing geometry through the build plate. Only what would end up below the - // model's own bottom is moved; downward relief above that is untouched. - bool v2_clamp_below_plate = false; + // Simplification target in thousands of triangles. -1 (the default) means automatic, from the + // same recommendation; 0 skips simplification, which is worth comparing on its own. + int v2_max_triangles_k = -1; // Slide vertices onto the texture's own edges before displacing. Displacement moves vertices along // the normal only, so without this a step in the image is reproduced wherever the triangle grid // happens to fall, as a staircase rather than a straight wall. bool v2_relocate = false; + // Data-dependent edge flips before displacement (see TextureBakeFlip.hpp). On by default and not a + // user setting; deliberately left out of serialize() so project files are unaffected. + bool v2_flip_edges = true; // Colour, all of which belongs to the stack rather than to any one layer: it is about how the // printer will realise the colours, not about which image they came from. @@ -397,12 +408,44 @@ struct TextureDisplacementOptions { int mix_mode = int(color_mix_mode); ar(displace_border, smooth_enabled, smooth_strength, smooth_iterations, smooth_skip_border, - pipeline_v2, v2_refine_mm, v2_regularize, v2_max_triangles_k, v2_clamp_below_plate, + pipeline_v2, v2_refine_mm, v2_regularize, v2_max_triangles_k, v2_relocate, color_mix_enabled, mix_mode, color_despeckle); color_mix_mode = ColorMixMode(mix_mode); } }; +// How much detail a height texture carries, as BumpMesh's smart resolution measures it: central +// differences of the grey image, the mean gradient and the share of texels steeper than 30 grey +// levels, mapped to how many texels one mesh edge may span (1 for a hard-edged image, 4 for a smooth +// one). Cached per image, like the decode. +struct TextureDetail +{ + float mean_gradient = 0.f; + float sharp_fraction = 0.f; + float pixels_per_edge = 4.f; +}; +TextureDetail analyze_texture_detail(const TextureDisplacementLayer &layer); + +// The default pipeline's automatic resolution: the refinement edge and the simplification budget the +// texture and the model call for, when the options leave them at "auto". +// - edge = texel size (tile / image width, in world mm, over the finest layer) x pixels per edge, +// but no finer than keeps the refinement under a 12 M triangle cap for this surface area, clamped +// to [0.05 mm, min(5 mm, diagonal / 50)] and rounded up to 0.01 mm; +// - budget = the triangle count an edge of that texel size needs over the surface, scaled by the +// relief depth (a gentle relief needs fewer), stepped to 10 k and clamped to [10 k, 2000 k]. +// `edge_mm` is 0 when no layer has a usable texture. +struct V2Resolution +{ + float edge_mm = 0.f; + int budget_k = 0; + float texel_mm = 0.f; + float pixels_per_edge = 0.f; + bool budget_bound = false; // the edge came from the triangle cap, not from the texture +}; +V2Resolution recommend_v2_resolution(const indexed_triangle_set &mesh, + const std::vector &layers, + const Transform3d &volume_to_world = Transform3d::Identity()); + // Decoded height (and, for a colour source image, colour) samples, independent of any GUI/OpenGL // texture object so they can be evaluated from a background bake Job as well as from GUI-side // preview code. @@ -601,11 +644,12 @@ Vec2f apply_island_transform(const Vec2f &uv, int chart, const PatchUnwrap &unwr // its vertex buffer at all. Eigen::Matrix island_transform_matrix(int chart, const PatchUnwrap &unwrap, const std::vector &islands); -// Lays the unwrap's charts out as a connected net: charts that share a mesh edge are unfolded so +// Lays the unwrap's charts out as connected nets: charts that share a mesh edge are unfolded so // their shared edge coincides (a cube -> its six faces joined along a spanning tree of edges, the rest // left as free borders). Charts stay separate islands, so their borders still show and any of them can -// still be moved by hand afterwards. A chart whose unfold would overlap one already placed is left -// where the packing put it. Returns one placement per chart. See the gizmo's auto-connect option. +// still be moved by hand afterwards. Each net grows from the largest chart not yet placed; a chart whose +// triangles would overlap the net stays out of it and starts a net of its own. The nets are then packed +// side by side. Returns one placement per chart. See the gizmo's auto-connect option. std::vector compute_connected_net(const PatchUnwrap &unwrap); // The placement that unfolds `child` onto `parent` along their shared mesh edge, honouring `parent`'s @@ -642,6 +686,13 @@ PatchUnwrap compute_patch_unwrap(const indexed_triangle_set &patch, float seam_a // both the layer's seam angle and its hand-placed islands. std::vector compute_lscm_uvs(const indexed_triangle_set &patch, const TextureDisplacementLayer &layer); +// The TextureDisplacementLayer::lscm_uv_overrides key for one unwrapped vertex (an index into PatchUnwrap::uvs). +inline int lscm_uv_override_key(int unwrapped_vertex) { return -(unwrapped_vertex + 1); } + +// Writes the overrides into `unwrap.uvs`: mesh-vertex keys onto every copy of their vertex, then unwrapped-vertex +// keys onto their one copy. Returns, per unwrapped vertex, whether an override set it. +std::vector apply_lscm_uv_overrides(PatchUnwrap &unwrap, const std::vector> &overrides); + // One paint mask (as stored by ModelVolume::texture_displacement_facets) per possible layer slot. using TextureDisplacementFacetsData = std::array; @@ -716,12 +767,34 @@ struct TextureColorRequest // straight to a TriangleSelector without a second mapping table. std::vector *out_triangle = nullptr; }; +// Where the volume sits on the plate: its instance transform times its own volume transform, i.e. +// mesh coordinates -> world millimetres. +// +// Every number the user sets is in real millimetres on the printed part - "Depth (mm)", "Tile size +// (mm)" - and the build plate is a world plane, so the bake runs in world space and transforms the +// result back. Doing it in the volume's own coordinates instead made a scaled instance stretch both +// the relief depth and the tiling by the scale factor, and under a non-uniform scale it also +// displaced along the wrong direction: a mesh normal maps to the world normal through the inverse +// transpose, not through the transform itself, so the relief leaned. Identity - the default - is +// exactly the old behaviour and is what an untransformed volume gives. indexed_triangle_set build_texture_displacement(const indexed_triangle_set &base_mesh, const std::vector &layers, const TextureDisplacementFacetsData &facets_data, const TextureDisplacementOptions &options = {}, const DisplacementProgressFn &progress = {}, - const TextureColorRequest *color = nullptr); + const TextureColorRequest *color = nullptr, + const Transform3d &volume_to_world = Transform3d::Identity(), + // When given and enabled, receives the mesh after each + // stage, already brought back into `base_mesh`'s frame. + BakeStageRecorder *debug = nullptr); + +// `volume`'s mesh coordinates -> world millimetres: its first instance's transform times its own. +// The mesh is shared by every instance, so a multi-instance object can only be baked for one of +// them; the first is what the gizmo edits against. Identity when the volume has no object yet. +Transform3d texture_displacement_volume_to_world(const ModelVolume &volume); +// The frame the bake and the previews project the texture in: `volume_to_world` with its translation +// removed, i.e. world orientation and scale about the volume's own origin. See build_texture_displacement(). +Transform3d texture_displacement_bake_frame(const Transform3d &volume_to_world); // Convenience overload for main-thread callers: extracts the mesh/layers/paint data/options from // `volume` and forwards to the overload above. @@ -881,7 +954,69 @@ indexed_triangle_set subdivide_mesh_adaptive(const indexed_triangle_set &mesh, float border_edge_length_mm = 0.f, const DisplacementProgressFn &progress = nullptr, const ColorFieldSampler &color = nullptr, - float color_edge_length_mm = 0.f); + float color_edge_length_mm = 0.f, + // Step mode, for a mesh that cut_mesh_at_steps() will cut + // next: a triangle one of whose edges crosses a *sharp* + // step of the height field (a jump of 40 % of the relief + // between two samples at the sample spacing, across the + // mid-level) gets a chord error of zero, since the step is + // the cutter's to reproduce and refining it only carpets + // the edge of the pattern. Features that no edge crosses + // yet still refine until one does. No effect outside + // feature mode. + bool split_multi_crossings = false); + +// Cuts `mesh` along the height field's mid-level contour wherever the field steps sharply across it, +// and doubles the seam, so that displacing the result produces a vertical wall at the step instead of +// a ramp across whichever triangle the step happened to fall in. +// +// This is what a binary height map - a grid, stripes, a knurl, wood grain as black-and-white bands - +// needs, and what refinement alone cannot give it. Refining a triangle that straddles a step never +// brings the chord error under any tolerance: the surface has a discontinuity, and a finer triangle only +// makes the ramp narrower. Every refinement level then leaves a band of its own size along every step, +// and the triangle budget ends up spent on carpeting the edges of the pattern, one to two orders of +// magnitude more triangles than the pattern needs, while the ramps are still visible. +// +// What it does, in order: +// - Samples the field at the vertices along the bake's normals and takes the mid-level as the contour. +// A vertex that sits inside a step's blend is nudged along the surface, down its own side's slope, +// until it samples a pure value, so no vertex bakes to a half height. +// - Marches every painted edge for crossings of the mid-level (sampled at half `step_width_mm`, then +// bisected). Both triangles at an edge see the same crossings, which keeps the result conformal. A +// crossing is *sharp* when the field changes by half its range within `step_width_mm`. +// - Traces the contour inside each triangle on a local lattice (marching triangles) and keeps the +// polylines that join one crossing to another; a closed loop inside a triangle is left to +// refinement. Each polyline is simplified and kept clear of its neighbours. +// - Doubles the seam: every polyline vertex and every crossing gets a copy on each side, moved +// `seam_gap_mm` apart across the contour to where each samples a pure value of its own side. The +// regions between the seams are ear-clipped; the strip between the two copies is triangulated flat +// and becomes the wall once the high side is displaced. A triangle whose contour could not be traced +// is split at its crossings without a seam, so the neighbours still meet it without a T-junction. +// +// Nothing is cut unless the texture is a step texture: most vertex heights sit at one of two levels, +// most crossings are sharp, and the features are wider than a couple of step widths (a noisy or a +// smooth relief is passed through unchanged, as is a mesh with nothing painted). If the cut would leave +// more than a trace of inverted triangles the input is returned unchanged as well. +// +// `region` flags the triangles that may be cut (the painted ones, same encoding as +// subdivide_mesh_adaptive(), any non-zero value). An edge may carry any number of crossings; two closer +// than `max(min_feature_mm, seam_gap_mm)` are a feature too thin to carry a seam and are dropped as a +// pair, leaving the surface flat there. `out_source` receives the input triangle each output triangle +// descends from, so paint carries over; `out_cut_count` the number of input triangles that were cut +// (0 means the mesh came back unchanged). +indexed_triangle_set cut_mesh_at_steps(const indexed_triangle_set &mesh, const std::vector ®ion, + const HeightFieldSampler &sampler, float step_width_mm, + float seam_gap_mm, float min_feature_mm = 0.f, + std::vector *out_source = nullptr, size_t *out_cut_count = nullptr); + +// The cutter's verdict alone: whether cut_mesh_at_steps() with the same arguments would cut anything, +// judged the same way (two levels, sharp crossings, features wider than the step) but without the +// vertex nudge or any tracing, so it is cheap enough to ask on the coarse mesh. The prepare path asks +// it *before* refining: a subdivision run in step mode leaves the steps alone for the cutter, which is +// only right if the cutter is then going to cut them. +bool texture_has_steps_to_cut(const indexed_triangle_set &mesh, const std::vector ®ion, + const HeightFieldSampler &sampler, float step_width_mm, float seam_gap_mm, + float min_feature_mm = 0.f); // The recipe for getting a mesh ready to receive displacement: even out the triangle density, then // refine it where the texture bends. Either stage is skipped when its target is <= 0. Pure data, and @@ -906,6 +1041,9 @@ struct TextureDisplacementPrepareParams // Separate from the height criteria because colour lands per facet: a flat surface carrying a // sharp colour edge needs triangles along that edge even though its height is perfectly smooth. float subdiv_color_edge_mm = 0.f; + // Cut the mesh along sharp steps in the texture after refining, so they bake as walls rather than + // ramps. See cut_mesh_at_steps(). Off leaves the old ramp behaviour. + bool cut_steps = false; }; // What a preparation run produced. An empty `mesh` means there was nothing to do and the caller must