TextureBake: edge flips along the height field, stage recorder, faster displace

src/libslic3r/CMakeLists.txt

@@ -470,6 +470,12 @@ set(lisbslic3r_sources
     TextureBake/TextureBakeSubdivide.hpp
     TextureBake/TextureBakeRegularize.cpp
     TextureBake/TextureBakeRegularize.hpp
+    TextureBake/TextureBakeRelocate.cpp
+    TextureBake/TextureBakeRelocate.hpp
+    TextureBake/TextureBakeFlip.cpp
+    TextureBake/TextureBakeFlip.hpp
+    TextureBake/TextureBakeDebug.cpp
+    TextureBake/TextureBakeDebug.hpp
     TextureBake/TextureBakeDisplace.cpp
     TextureBake/TextureBakeDisplace.hpp
     TextureBake/TextureBakeDecimate.cpp

src/libslic3r/TextureBake/TextureBakeDebug.cpp

@@ -0,0 +1,219 @@
+#include "TextureBakeDebug.hpp"
+
+#include <algorithm>
+#include <cinttypes>
+#include <cstdio>
+
+#include <boost/filesystem.hpp>
+#include <boost/log/trivial.hpp>
+
+namespace Slic3r {
+
+namespace {
+
+// Half-edge key, low index first so both sides of an edge form the same one.
+inline uint64_t edge_key(int a, int b)
+{
+    const uint32_t lo = uint32_t(std::min(a, b)), hi = uint32_t(std::max(a, b));
+    return (uint64_t(lo) << 32) | uint64_t(hi);
+}
+
+} // namespace
+
+void bake_stage_topology(const BakeStageMesh &mesh, size_t &open_edges, size_t &non_manifold_edges,
+                         size_t &degenerate)
+{
+    open_edges = non_manifold_edges = degenerate = 0;
+
+    // Sorted half-edges rather than a hash map: same answer, but it is one allocation and a sort
+    // instead of three million node allocations, which on a stage this size is the whole cost.
+    std::vector<uint64_t> keys;
+    keys.reserve(mesh.indices.size() * 3);
+    for (const Vec3i32 &t : mesh.indices) {
+        if (t[0] == t[1] || t[1] == t[2] || t[0] == t[2]) {
+            ++degenerate;
+            continue; // a collapsed face has no edges worth counting
+        }
+        const Vec3f &a = mesh.vertices[size_t(t[0])];
+        if ((mesh.vertices[size_t(t[1])] - a).cross(mesh.vertices[size_t(t[2])] - a).squaredNorm() <= 0.f)
+            ++degenerate; // zero area but three distinct corners: still counted as an edge carrier
+        for (int e = 0; e < 3; ++e)
+            keys.push_back(edge_key(t[e], t[(e + 1) % 3]));
+    }
+    std::sort(keys.begin(), keys.end());
+    for (size_t i = 0; i < keys.size();) {
+        size_t j = i + 1;
+        while (j < keys.size() && keys[j] == keys[i])
+            ++j;
+        const size_t incident = j - i;
+        if (incident == 1)
+            ++open_edges;
+        else if (incident > 2)
+            ++non_manifold_edges;
+        i = j;
+    }
+}
+
+void BakeStageRecorder::finish(BakeStageSnapshot &s)
+{
+    s.triangles = s.mesh.indices.size();
+    s.vertices  = s.mesh.vertices.size();
+    if (m_check_topology)
+        bake_stage_topology(s.mesh, s.open_edges, s.non_manifold_edges, s.degenerate);
+    s.topology_checked = m_check_topology;
+    if (s.triangles > m_mesh_cap) {
+        s.mesh_dropped = true;
+        s.mesh.vertices.clear();
+        s.mesh.vertices.shrink_to_fit();
+        s.mesh.indices.clear();
+        s.mesh.indices.shrink_to_fit();
+    }
+    m_stages.push_back(std::move(s));
+}
+
+void BakeStageRecorder::capture(const char *name, const TextureBake::TriSoup &soup, double ms,
+                                const std::string &detail)
+{
+    if (!m_enabled)
+        return;
+
+    BakeStageSnapshot s;
+    s.name   = name;
+    s.detail = detail;
+    s.ms     = ms;
+
+    // The pipeline works on non-indexed soup, so welding here is what turns it back into something
+    // renderable. The geometry grid, matching to_indexed_triangle_set(), so the debug view shows the
+    // same sharing the bake's own output would have.
+    const size_t                 n = soup.pos.size();
+    TextureBake::QuantizedPointMap map(TextureBake::WELD_GRID_GEOMETRY, std::min(n, size_t(1) << 22));
+    std::vector<int>             id(n);
+    s.mesh.vertices.reserve(n / 3);
+    for (size_t i = 0; i < n; ++i) {
+        id[i] = map.get_or_set(soup.pos[i], int(s.mesh.vertices.size()));
+        if (map.inserted())
+            s.mesh.vertices.push_back(soup.pos[i]);
+    }
+    s.mesh.indices.reserve(n / 3);
+    for (size_t t = 0; t + 2 < n; t += 3) {
+        // Corners that welded together carry no area; the bake's own conversion drops them too, so
+        // dropping them here keeps the stage count honest against what would be committed.
+        if (id[t] == id[t + 1] || id[t + 1] == id[t + 2] || id[t] == id[t + 2])
+            continue;
+        s.mesh.indices.emplace_back(id[t], id[t + 1], id[t + 2]);
+    }
+    finish(s);
+}
+
+void BakeStageRecorder::capture(const char *name, const std::vector<Vec3f> &vertices,
+                                const std::vector<Vec3i32> &indices, double ms,
+                                const std::string &detail)
+{
+    if (!m_enabled)
+        return;
+    BakeStageSnapshot s;
+    s.name          = name;
+    s.detail        = detail;
+    s.ms            = ms;
+    s.mesh.vertices = vertices;
+    s.mesh.indices  = indices;
+    finish(s);
+}
+
+void BakeStageRecorder::capture_note(const char *name, double ms, const std::string &detail)
+{
+    if (!m_enabled)
+        return;
+    BakeStageSnapshot s;
+    s.name             = name;
+    s.detail           = detail;
+    s.ms               = ms;
+    s.topology_checked = false;
+    m_stages.push_back(std::move(s));
+}
+
+void BakeStageRecorder::rebase(size_t from, const Transform3d *to_local, bool flip_winding)
+{
+    for (size_t i = from; i < m_stages.size(); ++i) {
+        BakeStageMesh &m = m_stages[i].mesh;
+        if (to_local != nullptr)
+            for (Vec3f &v : m.vertices)
+                v = (*to_local * v.cast<double>()).cast<float>();
+        if (flip_winding)
+            for (Vec3i32 &t : m.indices)
+                std::swap(t[1], t[2]);
+    }
+}
+
+double BakeStageRecorder::total_ms() const
+{
+    double sum = 0.0;
+    for (const BakeStageSnapshot &s : m_stages)
+        sum += s.ms;
+    return sum;
+}
+
+size_t dump_bake_stages(const std::vector<BakeStageSnapshot> &stages, const std::string &dir)
+{
+    boost::system::error_code ec;
+    boost::filesystem::create_directories(dir, ec);
+    if (ec) {
+        BOOST_LOG_TRIVIAL(error) << "BakeStageRecorder: cannot create " << dir << ": " << ec.message();
+        return 0;
+    }
+
+    size_t written = 0;
+    for (size_t i = 0; i < stages.size(); ++i) {
+        const BakeStageSnapshot &s = stages[i];
+        if (s.mesh.empty())
+            continue;
+
+        // Stage names carry spaces and punctuation; keep the filename to what every shell and viewer
+        // handles without quoting.
+        std::string safe;
+        for (const char c : s.name)
+            safe += (std::isalnum(static_cast<unsigned char>(c)) != 0) ? c : '_';
+
+        char path[1024];
+        std::snprintf(path, sizeof(path), "%s/%02zu_%s.obj", dir.c_str(), i, safe.c_str());
+        std::FILE *f = std::fopen(path, "wb");
+        if (f == nullptr) {
+            BOOST_LOG_TRIVIAL(error) << "BakeStageRecorder: cannot write " << path;
+            continue;
+        }
+        std::fprintf(f, "# texture bake stage %zu: %s\n", i, s.name.c_str());
+        if (!s.detail.empty())
+            std::fprintf(f, "# %s\n", s.detail.c_str());
+        std::fprintf(f, "# %zu triangles, %.2f ms\n", s.triangles, s.ms);
+        for (const Vec3f &v : s.mesh.vertices)
+            std::fprintf(f, "v %.6f %.6f %.6f\n", double(v.x()), double(v.y()), double(v.z()));
+        for (const Vec3i32 &t : s.mesh.indices) // OBJ indices are 1-based
+            std::fprintf(f, "f %d %d %d\n", t[0] + 1, t[1] + 1, t[2] + 1);
+        std::fclose(f);
+        ++written;
+    }
+
+    char summary[1024];
+    std::snprintf(summary, sizeof(summary), "%s/stages.txt", dir.c_str());
+    if (std::FILE *f = std::fopen(summary, "wb"); f != nullptr) {
+        std::fprintf(f, "%-3s %-24s %10s %12s %10s %8s %8s %8s  %s\n", "#", "stage", "ms", "triangles",
+                     "vertices", "open", "nonman", "degen", "detail");
+        double total = 0.0;
+        for (size_t i = 0; i < stages.size(); ++i) {
+            const BakeStageSnapshot &s = stages[i];
+            total += s.ms;
+            std::fprintf(f, "%-3zu %-24s %10.2f %12zu %10zu ", i, s.name.c_str(), s.ms, s.triangles,
+                         s.vertices);
+            if (s.topology_checked)
+                std::fprintf(f, "%8zu %8zu %8zu", s.open_edges, s.non_manifold_edges, s.degenerate);
+            else
+                std::fprintf(f, "%8s %8s %8s", "-", "-", "-");
+            std::fprintf(f, "  %s%s\n", s.detail.c_str(), s.mesh_dropped ? " [mesh over cap, not written]" : "");
+        }
+        std::fprintf(f, "\ntotal %.2f ms across %zu stages\n", total, stages.size());
+        std::fclose(f);
+    }
+    return written;
+}
+
+} // namespace Slic3r

src/libslic3r/TextureBake/TextureBakeDebug.hpp

@@ -0,0 +1,123 @@
+#pragma once
+
+// Step-by-step capture of a bake.
+//
+// A bake is a chain of stages that each rewrite the whole mesh, so when the result looks wrong the
+// only useful question is which stage made it wrong. This records the geometry, the wall time and the
+// topology after every stage, which is what the gizmo's debug view steps through and what the
+// benchmark's --dump-stages writes out.
+//
+// Deliberately independent of TriangleMesh: a stage is held as a plain vertex/index pair, which is
+// layout-compatible with indexed_triangle_set's own members (stl_vertex is Vec3f,
+// stl_triangle_vertex_indices is Vec3i32), so the GUI assigns rather than converts and the standalone
+// benchmark does not have to link admesh to use this.
+//
+// Recording is off unless enable(true) was called, and every capture site is a null-pointer check, so
+// a normal bake pays nothing for this being here.
+
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include "TextureBakeIndex.hpp"
+
+namespace Slic3r {
+
+struct BakeStageMesh
+{
+    std::vector<Vec3f>   vertices;
+    std::vector<Vec3i32> indices;
+
+    bool   empty() const { return indices.empty(); }
+    size_t triangle_count() const { return indices.size(); }
+};
+
+struct BakeStageSnapshot
+{
+    std::string   name;   // "remesh", "subdivide", ...
+    std::string   detail; // whatever the stage has to say: collapse counts, rejected moves, ...
+    BakeStageMesh mesh;   // empty when the stage was over the memory cap - see mesh_dropped
+
+    double ms        = 0.0;
+    size_t triangles = 0;
+    size_t vertices  = 0;
+
+    // Filled only when the recorder was asked to check topology: it is a sort over every half-edge,
+    // which on a multi-million triangle stage costs more than the stage being measured.
+    size_t open_edges         = 0;
+    size_t non_manifold_edges = 0;
+    size_t degenerate         = 0;
+    bool   topology_checked   = false;
+
+    // The geometry was dropped to stay inside the memory cap; every count above is still real.
+    bool mesh_dropped = false;
+};
+
+// Edge and area defects of a captured stage. Split out so a caller can run it on its own.
+void bake_stage_topology(const BakeStageMesh &mesh, size_t &open_edges, size_t &non_manifold_edges,
+                         size_t &degenerate);
+
+// Writes `<dir>/NN_name.obj` for every stage that still holds geometry, plus a `stages.txt` summary.
+// Returns how many meshes were written. Existing files with the same names are overwritten.
+//
+// A free function rather than a recorder method because by the time anyone wants the files the
+// recorder is usually gone and only the stages survive - that is how the gizmo holds them.
+size_t dump_bake_stages(const std::vector<BakeStageSnapshot> &stages, const std::string &dir);
+
+class BakeStageRecorder
+{
+public:
+    // Nothing is recorded until this is on.
+    void enable(bool on) { m_enabled = on; }
+    bool enabled() const { return m_enabled; }
+
+    // The edge scan is optional because it is O(n log n) over every half-edge, and a debug run that
+    // only wants to see the geometry should not pay for it on every stage.
+    void set_check_topology(bool on) { m_check_topology = on; }
+    bool check_topology() const { return m_check_topology; }
+
+    // Stages above this keep their counts but not their geometry. A debug run holds every stage at
+    // once, and a 4 M triangle stage is about 150 MB on its own, so without a cap stepping through a
+    // fine bake would need more memory than the bake did.
+    void   set_mesh_cap(size_t triangles) { m_mesh_cap = triangles; }
+    size_t mesh_cap() const { return m_mesh_cap; }
+
+    // `ms` is passed in rather than measured here: the caller is already timing the stage, and the
+    // capture itself (a weld, a copy, possibly an edge scan) must not land inside that measurement.
+    void capture(const char *name, const TextureBake::TriSoup &soup, double ms,
+                 const std::string &detail = {});
+    void capture(const char *name, const std::vector<Vec3f> &vertices,
+                 const std::vector<Vec3i32> &indices, double ms, const std::string &detail = {});
+    // For a stage that changed nothing a caller can still show, e.g. a skipped remesh.
+    void capture_note(const char *name, double ms, const std::string &detail);
+
+    // Index of the next stage to be recorded. Paired with rebase() to fix up a range afterwards.
+    size_t mark() const { return m_stages.size(); }
+
+    // Brings stages [from, end) into the caller's own space and winding. The bake runs in world
+    // millimetres and, for a mirrored placement, against a reversed winding; the debug view draws in
+    // the volume's local frame, so a captured range has to be brought back the same way the bake's
+    // own result is. `to_local` may be null for no transform.
+    void rebase(size_t from, const Transform3d *to_local, bool flip_winding);
+
+    const std::vector<BakeStageSnapshot> &stages() const { return m_stages; }
+    std::vector<BakeStageSnapshot>        take() { return std::move(m_stages); }
+    void                                  clear() { m_stages.clear(); }
+    bool                                  empty() const { return m_stages.empty(); }
+
+    // Total recorded wall time, which is the bake's own time minus whatever it does outside a stage.
+    double total_ms() const;
+
+    size_t dump_obj(const std::string &dir) const { return dump_bake_stages(m_stages, dir); }
+
+private:
+    void finish(BakeStageSnapshot &s);
+
+    std::vector<BakeStageSnapshot> m_stages;
+    bool                           m_enabled        = false;
+    bool                           m_check_topology = true;
+    size_t                         m_mesh_cap       = 4'000'000;
+};
+
+} // namespace Slic3r

src/libslic3r/TextureBake/TextureBakeDecimate.cpp

@@ -80,18 +80,17 @@ Vec3d face_normal_unit(const std::vector<Vec3d> &pos, int a, int b, int c)
 
 // Versions are captured at push time; a mismatch on pop means a later collapse invalidated the entry.
 // Lazy deletion, far cheaper than removing entries eagerly - but it means the heap accumulates stale
-// duplicates, so it grows to several times the edge count and its size has to be reserved up front.
-// Left to grow on its own it reallocates and copies the whole array repeatedly, which on a
-// multi-million-entry heap costs more than every collapse put together.
+// duplicates, so its size has to be reserved up front and it is compacted once the dead entries
+// dominate (see maybe_compact below).
 //
-// The collapse target is stored as float rather than double: it is a position on a mesh already held
-// in float, and halving the entry cuts the memory the sift operations drag through cache.
+// The collapse target is not stored. A matching version stamp means neither endpoint's quadric nor its
+// position has changed since the push, so the target recomputes to exactly the same value on pop - and
+// the entry drops from 40 bytes to 24. Sifting is most of this stage's time, and it is memory traffic.
 struct HeapEntry
 {
     double   cost;
     int      v1, v2;
     uint32_t ver1, ver2;
-    Vec3f    p;
     bool     operator>(const HeapEntry &o) const { return cost > o.cost; }
 };
 
@@ -289,9 +288,31 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
         heap.pop_back();
         return e;
     };
-    size_t pops = 0, stale_pops = 0;
+    size_t pops = 0, stale_pops = 0, compactions = 0;
 
-    const auto push_edge = [&](int v1, int v2) {
+    // An entry is stale once either endpoint has been removed or moved by a later collapse.
+    const auto is_stale = [&](const HeapEntry &e) {
+        return !active[size_t(e.v1)] || !active[size_t(e.v2)] || version[size_t(e.v1)] != e.ver1 ||
+               version[size_t(e.v2)] != e.ver2;
+    };
+    // Measured on a 2.4 M -> 750 k run, 82% of pops were stale: every collapse re-pushes the survivor's
+    // edges and orphans the old ones, so the heap grows to several times the live edge set and every
+    // sift walks that much further through memory. Dropping the dead entries and re-heapifying once
+    // they dominate costs one linear pass, amortised against the growth that triggered it.
+    //
+    // Keyed to the live face count rather than to the heap's own size: lazy popping keeps the heap from
+    // ever doubling, but the live edge set (about 1.5 per face) shrinks as decimation proceeds, so by the
+    // end the heap is several times what is still collapsible. Compact once it passes twice that.
+    const auto maybe_compact = [&]() {
+        if (heap.size() < std::max<size_t>(size_t(1) << 16, active_faces * 3))
+            return;
+        heap.erase(std::remove_if(heap.begin(), heap.end(), is_stale), heap.end());
+        std::make_heap(heap.begin(), heap.end(), std::greater<HeapEntry>());
+        ++compactions;
+    };
+
+    // Where an edge collapses to. Also re-run on pop instead of stored - see HeapEntry.
+    const auto collapse_target = [&](int v1, int v2) -> Vec3d {
         Vec3d p;
         if (!solve_q(quadrics, v1, v2, p)) {
             const Vec3d  mid = (pos[size_t(v1)] + pos[size_t(v2)]) * 0.5;
@@ -306,10 +327,16 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
             else if (e1 <= e2)          p = pos[size_t(v1)];
             else                        p = pos[size_t(v2)];
         }
+        return p;
+    };
+    const auto push_edge = [&](int v1, int v2) {
+        const Vec3d p = collapse_target(v1, v2);
+        // The cost is evaluated at the exact target and the collapse moves to its float rounding -
+        // the same split as when the rounded target was stored in the entry, so no ordering changes.
         // Where quadric costs are all near zero, shorter edges first keeps triangle quality up.
         const double len2 = (pos[size_t(v2)] - pos[size_t(v1)]).squaredNorm();
         heap_push({ eval_sum(quadrics, v1, v2, p) + len2 * 1e-8, v1, v2, version[size_t(v1)],
-                    version[size_t(v2)], p.cast<float>() });
+                    version[size_t(v2)] });
     };
 
     {
@@ -430,11 +457,7 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
             break;
 
         const int v1 = top.v1, v2 = top.v2;
-        if (!active[size_t(v1)] || !active[size_t(v2)]) {
-            ++stale_pops;
-            continue;
-        }
-        if (version[size_t(v1)] != top.ver1 || version[size_t(v2)] != top.ver2) {
+        if (is_stale(top)) {
             ++stale_pops;
             continue;
         }
@@ -443,7 +466,7 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
         lk_epoch += 2; // +2 so ep and ep+1 cannot collide with the next call
         if (has_link_violation(v1, v2, lk_epoch))
             continue;
-        const Vec3d target = top.p.cast<double>();
+        const Vec3d target = collapse_target(v1, v2).cast<float>().cast<double>();
         if (check_flipped(v1, v2, target) || check_flipped(v2, v1, target))
             continue;
 
@@ -492,6 +515,7 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
                     push_edge(v1, nb);
             }
         }
+        maybe_compact();
 
         if (on_progress) {
             const double p = std::min(1.0, double(init_faces - active_faces) / double(to_remove));
@@ -504,7 +528,8 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
     }
 
     BOOST_LOG_TRIVIAL(info) << "TextureBake decimate: pops=" << pops << " stale=" << stale_pops
-                            << " heap_peak=" << heap.capacity() << " faces=" << active_faces;
+                            << " compactions=" << compactions << " heap_peak=" << heap.capacity()
+                            << " faces=" << active_faces;
 
     // Rebuild from the surviving faces, with per-face normals.
     TriSoup &out = result.geometry;

src/libslic3r/TextureBake/TextureBakeDisplace.cpp

@@ -120,20 +120,23 @@ TriSoup apply_displacement(const TriSoup &geometry, const HeightSampleFn &sample
 
         std::vector<Vec3d> cur = smooth_nrm, nxt(unique_count, Vec3d::Zero());
         for (int iter = 0; iter < settings.blend_normal_smoothing; ++iter) {
-            for (size_t id = 0; id < unique_count; ++id) {
-                const uint32_t s = csr_start[id], e = csr_start[id + 1];
-                if (e == s) {
-                    nxt[id] = cur[id];
-                    continue;
+            // Jacobi, so every vertex reads the previous iteration and the rows are independent.
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, unique_count, 4096), [&](const tbb::blocked_range<size_t> &r) {
+                for (size_t id = r.begin(); id < r.end(); ++id) {
+                    const uint32_t s = csr_start[id], e = csr_start[id + 1];
+                    if (e == s) {
+                        nxt[id] = cur[id];
+                        continue;
+                    }
+                    Vec3d sum = Vec3d::Zero();
+                    for (uint32_t k = s; k < e; ++k)
+                        sum += cur[neighbors[k]];
+                    sum /= double(e - s);
+                    const double len = sum.norm();
+                    // Cancelling neighbours mean a knife edge; keep what we had.
+                    nxt[id] = (len > 1e-12) ? Vec3d(sum / len) : cur[id];
                 }
-                Vec3d sum = Vec3d::Zero();
-                for (uint32_t k = s; k < e; ++k)
-                    sum += cur[neighbors[k]];
-                sum /= double(e - s);
-                const double len = sum.norm();
-                // Cancelling neighbours mean a knife edge; keep what we had.
-                nxt[id] = (len > 1e-12) ? Vec3d(sum / len) : cur[id];
-            }
+            });
             cur.swap(nxt);
         }
         blend_nrm = std::move(cur);

src/libslic3r/TextureBake/TextureBakeFlip.cpp

@@ -0,0 +1,212 @@
+#include "TextureBakeFlip.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+#include <tbb/parallel_sort.h>
+
+namespace Slic3r {
+namespace TextureBake {
+
+FlipResult flip_edges_to_height(const TriSoup &geometry, const std::vector<int> &face_parent_id,
+                                const HeightSampleFn &sample, const FlipSettings &settings,
+                                const std::vector<uint8_t> &locked)
+{
+    FlipResult result;
+    result.geometry       = geometry;
+    result.face_parent_id = face_parent_id;
+    const size_t count    = geometry.pos.size();
+    const size_t tri_ct   = count / 3;
+    if (count == 0 || !sample || settings.passes <= 0)
+        return result;
+
+    // Weld, so a flip rewrites triangles in terms of shared vertices rather than positions.
+    QuantizedPointMap  weld(WELD_GRID_GEOMETRY, std::min(count, size_t(1) << 22));
+    std::vector<int>   vid(count);
+    std::vector<Vec3f> pos;
+    std::vector<float> weight; // exclude weight per unique vertex, the max over its copies
+    const bool         has_weight = !geometry.exclude_weight.empty();
+    for (size_t i = 0; i < count; ++i) {
+        vid[i] = weld.get_or_set(geometry.pos[i], int(pos.size()));
+        if (weld.inserted()) {
+            pos.push_back(geometry.pos[i]);
+            weight.push_back(has_weight ? geometry.exclude_weight[i] : 0.f);
+        } else if (has_weight) {
+            weight[size_t(vid[i])] = std::max(weight[size_t(vid[i])], geometry.exclude_weight[i]);
+        }
+    }
+    const size_t nv = pos.size();
+
+    // Triangles as vertex ids; a locked or excluded triangle never takes part.
+    std::vector<std::array<int, 3>> tri(tri_ct);
+    std::vector<uint8_t>            fixed(tri_ct, 0);
+    for (size_t t = 0; t < tri_ct; ++t) {
+        tri[t] = { vid[t * 3], vid[t * 3 + 1], vid[t * 3 + 2] };
+        if (!locked.empty() && t < locked.size() && locked[t])
+            fixed[t] = 1;
+        if (has_weight && geometry.exclude_weight[t * 3] > 0.99f)
+            fixed[t] = 1;
+    }
+
+    // Height per unique vertex along its area-weighted normal, the direction displacement will use.
+    std::vector<Vec3f> nrm(nv, Vec3f::Zero());
+    std::vector<Vec3f> face_n(tri_ct);
+    const auto rebuild_normals = [&]() {
+        std::fill(nrm.begin(), nrm.end(), Vec3f::Zero());
+        for (size_t t = 0; t < tri_ct; ++t) {
+            const Vec3f fn = (pos[size_t(tri[t][1])] - pos[size_t(tri[t][0])]).cross(pos[size_t(tri[t][2])] - pos[size_t(tri[t][0])]);
+            face_n[t]      = fn;
+            for (int k = 0; k < 3; ++k)
+                nrm[size_t(tri[t][size_t(k)])] += fn;
+        }
+        for (Vec3f &n : nrm) {
+            const float l = n.norm();
+            n = (l > 0.f) ? Vec3f(n / l) : Vec3f(0.f, 0.f, 1.f);
+        }
+    };
+    rebuild_normals();
+    std::vector<float> h(nv, 0.f);
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, nv), [&](const tbb::blocked_range<size_t> &r) {
+        for (size_t v = r.begin(); v < r.end(); ++v)
+            h[v] = sample(pos[v], nrm[v], nrm[v]);
+    });
+    float h_lo = std::numeric_limits<float>::max(), h_hi = -h_lo;
+    for (const float x : h) { h_lo = std::min(h_lo, x); h_hi = std::max(h_hi, x); }
+    const float range = h_hi - h_lo;
+    if (!(range > 0.f))
+        return result; // flat: every diagonal is as good as the other
+    const float min_gain = float(settings.min_gain_fraction) * range;
+    const float planar   = float(settings.min_planar_cos);
+
+    struct Candidate
+    {
+        uint32_t t1, t2;   // the two triangles
+        uint8_t  k1, k2;   // corner index in each where the shared edge starts (t1: a->c, t2: c->a)
+        float    gain;
+    };
+
+    for (int pass = 0; pass < settings.passes; ++pass) {
+        // Half-edges keyed by their undirected edge, sorted so the two halves of an interior edge land
+        // next to each other; a run of exactly two with opposite directions is a manifold interior
+        // edge. Sorting beats hashing here by an order of magnitude on a few million triangles.
+        struct Half { uint64_t key; uint32_t corner; };
+        std::vector<Half> half;
+        half.reserve(count);
+        for (size_t t = 0; t < tri_ct; ++t)
+            for (int k = 0; k < 3; ++k) {
+                const int from = tri[t][size_t(k)], to = tri[t][size_t((k + 1) % 3)];
+                if (from == to) continue;
+                const uint32_t lo = uint32_t(std::min(from, to)), hi = uint32_t(std::max(from, to));
+                half.push_back({ (uint64_t(lo) << 32) | hi, uint32_t(t * 3 + size_t(k)) });
+            }
+        tbb::parallel_sort(half.begin(), half.end(), [](const Half &x, const Half &y) {
+            return x.key != y.key ? x.key < y.key : x.corner < y.corner;
+        });
+        std::vector<std::pair<uint32_t, uint32_t>> edges; // (corner in t1, corner in t2), t1 < t2
+        edges.reserve(half.size() / 2);
+        for (size_t i = 0; i < half.size();) {
+            size_t j = i + 1;
+            while (j < half.size() && half[j].key == half[i].key) ++j;
+            if (j - i == 2) {
+                // Opposite directions: the lower vertex id is `from` in exactly one of the two.
+                const uint32_t c1 = half[i].corner, c2 = half[i + 1].corner;
+                const int f1 = tri[c1 / 3][size_t(c1 % 3)], f2 = tri[c2 / 3][size_t(c2 % 3)];
+                if (f1 != f2)
+                    edges.emplace_back(c1, c2);
+            }
+            i = j;
+        }
+        std::vector<Candidate> cands(edges.size());
+        std::vector<uint8_t>   valid(edges.size(), 0);
+        tbb::parallel_for(tbb::blocked_range<size_t>(0, edges.size()), [&](const tbb::blocked_range<size_t> &r) {
+            for (size_t i = r.begin(); i < r.end(); ++i) {
+                const uint32_t c1 = uint32_t(edges[i].first), c2 = uint32_t(edges[i].second);
+                const uint32_t t1 = c1 / 3, t2 = c2 / 3;
+                const int      k1 = int(c1 % 3), k2 = int(c2 % 3);
+                if (fixed[t1] || fixed[t2]) continue;
+                if (!face_parent_id.empty() && face_parent_id[t1] != face_parent_id[t2]) continue;
+                // a->c is the shared edge in t1, with b opposite; t2 runs c->a with d opposite.
+                const int a = tri[t1][size_t(k1)], c = tri[t1][size_t((k1 + 1) % 3)], b = tri[t1][size_t((k1 + 2) % 3)];
+                const int d = tri[t2][size_t((k2 + 2) % 3)];
+                if (b == d) continue;
+                // Level quads have nothing to gain; the test is on the corners, before any sampling.
+                const float hmin = std::min({ h[size_t(a)], h[size_t(b)], h[size_t(c)], h[size_t(d)] });
+                const float hmax = std::max({ h[size_t(a)], h[size_t(b)], h[size_t(c)], h[size_t(d)] });
+                if (hmax - hmin < min_gain) continue;
+                // Coplanar enough to have a real alternative, and convex so the alternative is valid:
+                // the new triangles (a, d, b) and (d, c, b) must both face the way the quad does, with
+                // a decent share of its area.
+                const Vec3f n1 = face_n[t1], n2 = face_n[t2];
+                const float l1 = n1.norm(), l2 = n2.norm();
+                if (l1 <= 0.f || l2 <= 0.f || n1.dot(n2) < planar * l1 * l2) continue;
+                const Vec3f quad_n = (n1 + n2).normalized();
+                const Vec3f &pa = pos[size_t(a)], &pb = pos[size_t(b)], &pc = pos[size_t(c)], &pd = pos[size_t(d)];
+                const Vec3f m1 = (pd - pa).cross(pb - pa), m2 = (pc - pd).cross(pb - pd);
+                const float area_old = l1 + l2, area_new = m1.dot(quad_n) + m2.dot(quad_n);
+                const float min_part = 0.05f * area_old;
+                if (m1.dot(quad_n) < min_part || m2.dot(quad_n) < min_part) continue;
+                if (std::abs(area_new - area_old) > 0.02f * area_old) continue; // not the same quad: folded
+                // The diagonals' midpoint errors.
+                const auto mid_err = [&](int u, int w) {
+                    const Vec3f p = 0.5f * (pos[size_t(u)] + pos[size_t(w)]);
+                    Vec3f       n = nrm[size_t(u)] + nrm[size_t(w)];
+                    const float l = n.norm();
+                    n = (l > 0.f) ? Vec3f(n / l) : quad_n;
+                    return std::abs(sample(p, n, n) - 0.5f * (h[size_t(u)] + h[size_t(w)]));
+                };
+                const float gain = mid_err(a, c) - mid_err(b, d);
+                if (gain < min_gain) continue;
+                cands[i] = { t1, t2, uint8_t(k1), uint8_t(k2), gain };
+                valid[i] = 1;
+            }
+        });
+        std::vector<Candidate> chosen;
+        for (size_t i = 0; i < cands.size(); ++i)
+            if (valid[i]) chosen.push_back(cands[i]);
+        std::sort(chosen.begin(), chosen.end(), [](const Candidate &x, const Candidate &y) {
+            return x.gain != y.gain ? x.gain > y.gain : (x.t1 != y.t1 ? x.t1 < y.t1 : x.t2 < y.t2);
+        });
+        // Best first, and a triangle changes at most once per pass.
+        std::vector<uint8_t> touched(tri_ct, 0);
+        size_t               applied = 0;
+        for (const Candidate &cd : chosen) {
+            if (touched[cd.t1] || touched[cd.t2]) continue;
+            const int a = tri[cd.t1][cd.k1], c = tri[cd.t1][size_t((cd.k1 + 1) % 3)], b = tri[cd.t1][size_t((cd.k1 + 2) % 3)];
+            const int d = tri[cd.t2][size_t((cd.k2 + 2) % 3)];
+            tri[cd.t1] = { a, d, b };
+            tri[cd.t2] = { d, c, b };
+            touched[cd.t1] = touched[cd.t2] = 1;
+            ++applied;
+        }
+        result.flipped += applied;
+        if (applied == 0)
+            break;
+        rebuild_normals(); // face normals feed the planarity test of the next pass
+    }
+
+    if (result.flipped == 0)
+        return result;
+
+    // Back to the soup: positions, weights and per-face normals from the (possibly rewritten) triangles.
+    for (size_t t = 0; t < tri_ct; ++t) {
+        for (int k = 0; k < 3; ++k) {
+            const size_t i = t * 3 + size_t(k);
+            const int    v = tri[t][size_t(k)];
+            result.geometry.pos[i] = pos[size_t(v)];
+            if (has_weight)
+                result.geometry.exclude_weight[i] = weight[size_t(v)];
+        }
+        Vec3f       n   = (result.geometry.pos[t * 3 + 1] - result.geometry.pos[t * 3]).cross(result.geometry.pos[t * 3 + 2] - result.geometry.pos[t * 3]);
+        const float len = n.norm();
+        n = (len > 0.f) ? Vec3f(n / len) : Vec3f(0.f, 0.f, 1.f);
+        result.geometry.nrm[t * 3] = result.geometry.nrm[t * 3 + 1] = result.geometry.nrm[t * 3 + 2] = n;
+    }
+    return result;
+}
+
+} // namespace TextureBake
+} // namespace Slic3r

src/libslic3r/TextureBake/TextureBakeFlip.hpp

@@ -0,0 +1,57 @@
+#pragma once
+
+// Data-dependent edge flipping: choose each quad's diagonal to follow the height field, before any
+// displacement happens.
+//
+// Refinement produces a regular grid, and a step in the height map that crosses that grid at an
+// angle lands on alternating corners: one triangle of a quad gets a raised corner, the next does not,
+// and the displaced wall comes out as a sawtooth the size of the grid. Finer triangles make the teeth
+// smaller, never straight. The cause is the diagonal, not the density: with the diagonal running
+// along the step both triangles of the quad sit cleanly on one side or the other, and the wall is a
+// straight line between them.
+//
+// So, for every interior edge, compare the two diagonals of the quad it spans by how well each
+// interpolates the field at its own midpoint - the sampled height there against the mean of its two
+// endpoints - and keep the better one. A quad whose four corners are level is skipped outright, so
+// flat regions cost nothing; a non-planar quad (a model crease) is never touched, nor is one whose
+// triangles belong to different source faces or straddle the painted boundary.
+
+#include <cstdint>
+#include <vector>
+
+#include "TextureBakeDisplace.hpp"
+#include "TextureBakeIndex.hpp"
+
+namespace Slic3r {
+namespace TextureBake {
+
+struct FlipSettings
+{
+    // Passes over all edges. Flips interact through shared triangles, so a pass applies non-conflicting
+    // ones and the next pass picks up the rest; two or three settle a grid.
+    int passes = 3;
+
+    // A flip has to reduce the midpoint error by at least this fraction of the height range, so noise
+    // on a rough surface does not toggle diagonals for nothing.
+    double min_gain_fraction = 0.02;
+
+    // The two triangles have to be this coplanar (cosine of their normals' angle) for the quad to have
+    // a meaningful alternative diagonal at all: across a real crease there is none.
+    double min_planar_cos = 0.985; // ~10 degrees
+};
+
+struct FlipResult
+{
+    TriSoup          geometry;
+    std::vector<int> face_parent_id;
+    size_t           flipped = 0;
+};
+
+// `face_parent_id` may be empty; when given it is carried through unchanged (a flip never crosses a
+// parent boundary). `locked` flags triangles that must not change (per triangle, may be empty).
+FlipResult flip_edges_to_height(const TriSoup &geometry, const std::vector<int> &face_parent_id,
+                                const HeightSampleFn &sample, const FlipSettings &settings,
+                                const std::vector<uint8_t> &locked);
+
+} // namespace TextureBake
+} // namespace Slic3r

src/libslic3r/TextureBake/TextureBakeIndex.hpp

@@ -9,6 +9,7 @@
 // vertices distinct (they merge at 100 um, giving needle artifacts after displacement) while still
 // absorbing float noise; 1 um is what collapse positioning needs.
 
+#include <algorithm>
 #include <cmath>
 #include <cstdint>
 #include <vector>
@@ -26,8 +27,12 @@ static constexpr double WELD_GRID_DECIMATION = 1e6; // 1 um
 // tie rule matters.
 inline int64_t grid_round(double v) { return int64_t(std::floor(v + 0.5)); }
 
-// Open-addressing table over flat arrays: no allocation per lookup, exact integer key comparison.
+// Open-addressing table, linear probing: no allocation per lookup, exact integer key comparison.
 // Values must be non-negative; -1 is the empty sentinel and what get() returns on a miss.
+//
+// Key and value live together in one 32-byte cell. They used to be four parallel arrays, which made a
+// single probe touch four cache lines - and probing this table was 13% of a whole bake, because every
+// stage welds the full soup through it.
 class QuantizedPointMap
 {
 public:
@@ -46,8 +51,8 @@ public:
 
     int get(float x, float y, float z)
     {
-        return m_val[slot(grid_round(double(x) * m_quant), grid_round(double(y) * m_quant),
-                          grid_round(double(z) * m_quant))];
+        return m_cells[slot(grid_round(double(x) * m_quant), grid_round(double(y) * m_quant),
+                            grid_round(double(z) * m_quant))].val;
     }
     int get(const Vec3f &p) { return get(p.x(), p.y(), p.z()); }
 
@@ -55,35 +60,24 @@ public:
     // return it. inserted() then says which of the two happened.
     int get_or_set(float x, float y, float z, int value)
     {
-        const int64_t qx = grid_round(double(x) * m_quant);
-        const int64_t qy = grid_round(double(y) * m_quant);
-        const int64_t qz = grid_round(double(z) * m_quant);
-        const size_t  i  = slot(qx, qy, qz);
-        if (m_val[i] != -1) {
-            m_inserted = false;
-            return m_val[i];
-        }
-        m_qx[i] = qx; m_qy[i] = qy; m_qz[i] = qz;
-        m_val[i]   = value;
-        m_inserted = true;
-        if (++m_size > size_t(double(m_cap) * 0.7))
-            grow();
-        return value;
+        return get_or_set_key(grid_round(double(x) * m_quant), grid_round(double(y) * m_quant),
+                              grid_round(double(z) * m_quant), value);
     }
     int get_or_set(const Vec3f &p, int value) { return get_or_set(p.x(), p.y(), p.z(), value); }
 
     // The same table as a set of integer tuples (edge marking, midpoint cache). Quantisation is
     // bypassed: routing ids through the float overloads loses precision above 2^24.
-    int get_key(int64_t a, int64_t b, int64_t c) { return m_val[slot(a, b, c)]; }
+    int get_key(int64_t a, int64_t b, int64_t c) { return m_cells[slot(a, b, c)].val; }
     int get_or_set_key(int64_t a, int64_t b, int64_t c, int value)
     {
         const size_t i = slot(a, b, c);
-        if (m_val[i] != -1) {
+        Cell        &cell = m_cells[i];
+        if (cell.val != -1) {
             m_inserted = false;
-            return m_val[i];
+            return cell.val;
         }
-        m_qx[i] = a; m_qy[i] = b; m_qz[i] = c;
-        m_val[i]   = value;
+        cell.qx = a; cell.qy = b; cell.qz = c;
+        cell.val   = value;
         m_inserted = true;
         if (++m_size > size_t(double(m_cap) * 0.7))
             grow();
@@ -91,27 +85,34 @@ public:
     }
 
 private:
+    struct Cell
+    {
+        int64_t qx = 0, qy = 0, qz = 0;
+        int32_t val = -1;
+    };
+
     void alloc(size_t cap)
     {
         m_cap  = cap;
         m_mask = cap - 1;
-        m_qx.assign(cap, 0);
-        m_qy.assign(cap, 0);
-        m_qz.assign(cap, 0);
-        m_val.assign(cap, -1);
+        m_cells.assign(cap, Cell{});
     }
 
     size_t slot(int64_t qx, int64_t qy, int64_t qz) const
     {
-        uint32_t h = uint32_t(int32_t(qx) * int32_t(0x9E3779B1)) ^
-                     uint32_t(int32_t(qy) * int32_t(0x85EBCA77)) ^
-                     uint32_t(int32_t(qz) * int32_t(0xC2B2AE3D));
+        // Unsigned multiplies: the signed versions overflowed on nearly every key, which is undefined
+        // behaviour. The resulting bits are identical on every target OrcaSlicer builds for.
+        //
+        // A stronger 64-bit finalizer was tried and measured no faster - the probing that shows up in a
+        // profile is subdivide's parallel mark count, spread over every core, not long probe chains.
+        uint32_t h = (uint32_t(qx) * 0x9E3779B1u) ^ (uint32_t(qy) * 0x85EBCA77u) ^ (uint32_t(qz) * 0xC2B2AE3Du);
         h ^= h >> 15;
         size_t i = size_t(h) & m_mask;
         // Equality is checked against the stored 64-bit keys, so truncating to 32 bits for the hash
         // costs collisions at worst, never a wrong answer.
-        while (m_val[i] != -1) {
-            if (m_qx[i] == qx && m_qy[i] == qy && m_qz[i] == qz)
+        while (m_cells[i].val != -1) {
+            const Cell &c = m_cells[i];
+            if (c.qx == qx && c.qy == qy && c.qz == qz)
                 return i;
             i = (i + 1) & m_mask;
         }
@@ -120,24 +121,17 @@ private:
 
     void grow()
     {
-        std::vector<int64_t> oqx = std::move(m_qx), oqy = std::move(m_qy), oqz = std::move(m_qz);
-        std::vector<int>     oval = std::move(m_val);
-        const size_t         ocap = m_cap;
-        alloc(ocap * 2);
-        for (size_t i = 0; i < ocap; ++i) {
-            if (oval[i] == -1)
-                continue;
-            const size_t s = slot(oqx[i], oqy[i], oqz[i]);
-            m_qx[s] = oqx[i]; m_qy[s] = oqy[i]; m_qz[s] = oqz[i];
-            m_val[s] = oval[i];
-        }
+        std::vector<Cell> old = std::move(m_cells);
+        alloc(m_cap * 2);
+        for (const Cell &c : old)
+            if (c.val != -1)
+                m_cells[slot(c.qx, c.qy, c.qz)] = c;
     }
 
-    double               m_quant;
-    size_t               m_cap = 0, m_mask = 0, m_size = 0;
-    bool                 m_inserted = false;
-    std::vector<int64_t> m_qx, m_qy, m_qz;
-    std::vector<int>     m_val;
+    double            m_quant;
+    size_t            m_cap = 0, m_mask = 0, m_size = 0;
+    bool              m_inserted = false;
+    std::vector<Cell> m_cells;
 };
 
 // Three consecutive entries per triangle. The indexers turn this into shared vertices where a stage

src/libslic3r/TextureBake/TextureBakePipeline.cpp

@@ -1,8 +1,11 @@
 #include "TextureBakePipeline.hpp"
 
+#include "TextureBakeDebug.hpp"
+
 #include <algorithm>
 #include <chrono>
 #include <cmath>
+#include <string>
 
 #include <boost/log/trivial.hpp>
 
@@ -108,7 +111,7 @@ size_t snap_bottom_to_flat(TriSoup &geometry, float bottom_z, double tol)
 PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
                             const PipelineSettings &settings, const DisplaceBounds &bounds,
                             PipelineMode mode, const std::vector<uint8_t> &face_excluded,
-                            const PipelineProgressFn &on_progress)
+                            const PipelineProgressFn &on_progress, BakeStageRecorder *debug)
 {
     PipelineResult result;
     const auto     report = [&](const char *stage, double f) {
@@ -118,9 +121,16 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
     // without this it is guesswork which one to attack.
     auto       clock_now = [] { return std::chrono::steady_clock::now(); };
     auto       t_stage   = clock_now();
-    const auto lap       = [&](const char *stage, size_t tris) {
+    // One call site for both the log line and the debug capture, so a stage cannot appear in one and
+    // be missing from the other. The capture happens after the elapsed time is read: welding the soup
+    // and scanning its edges costs more than some of the stages do, and must not land inside the
+    // measurement it is reporting.
+    const auto lap = [&](const char *stage, const TriSoup &geometry, const std::string &detail = {}) {
         const double ms = std::chrono::duration<double, std::milli>(clock_now() - t_stage).count();
-        BOOST_LOG_TRIVIAL(info) << "TextureBake " << stage << ": " << ms << " ms, " << tris << " tris";
+        BOOST_LOG_TRIVIAL(info) << "TextureBake " << stage << ": " << ms << " ms, "
+                                << geometry.triangle_count() << " tris";
+        if (debug != nullptr)
+            debug->capture(stage, geometry, ms, detail);
         t_stage = clock_now();
     };
 
@@ -128,13 +138,16 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
         result.geometry = input;
         return result;
     }
+    if (debug != nullptr)
+        debug->capture("input", input, 0.0, "as handed to the pipeline");
+    t_stage = clock_now(); // the capture above is not part of the first stage
 
     // 1. Refine to the target edge length.
     SubdivideResult sub = subdivide(
         input, settings.refine_length, face_excluded, /* fast */ false, settings.safety_cap,
         [&](double f, size_t, double) { return report("subdivide", f); });
     result.safety_cap_hit = sub.safety_cap_hit;
-    lap("subdivide", sub.geometry.triangle_count());
+    lap("subdivide", sub.geometry);
     if (!report("subdivide", 1.0)) {
         result.canceled = true;
         return result;
@@ -147,7 +160,7 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
         RegularizeResult reg = regularize_mesh(sub.geometry, sub.face_parent_id,
                                                settings.refine_length, ropts);
         result.collapse_count = reg.collapse_count;
-        lap("regularize", reg.geometry.triangle_count());
+        lap("regularize", reg.geometry, std::to_string(reg.collapse_count) + " collapses");
         if (!report("regularize", 1.0)) {
             result.canceled = true;
             return result;
@@ -173,7 +186,7 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
                                   ? reg.face_parent_id[size_t(mid)] : -1;
             }
             sub.face_parent_id = std::move(composed);
-            lap("re-subdivide", sub.geometry.triangle_count());
+            lap("re-subdivide", sub.geometry);
         } else {
             sub.geometry       = std::move(reg.geometry);
             sub.face_parent_id = std::move(reg.face_parent_id);
@@ -192,12 +205,31 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
         BOOST_LOG_TRIVIAL(info) << "TextureBake relocate: moved=" << rel.moved
                                 << " rejected=" << rel.rejected;
         sub.geometry = std::move(rel.geometry);
-        lap("relocate", sub.geometry.triangle_count());
+        lap("relocate", sub.geometry,
+            "moved " + std::to_string(rel.moved) + ", rejected " + std::to_string(rel.rejected));
+    }
+
+    // 3b. Diagonals along the height field's steps, so they displace into straight walls.
+    if (settings.flip_edges) {
+        std::vector<uint8_t> locked;
+        if (settings.preserve_untextured && !sub.geometry.exclude_weight.empty()) {
+            locked.assign(sub.geometry.triangle_count(), 0);
+            for (size_t t = 0; t < locked.size(); ++t)
+                locked[t] = sub.geometry.exclude_weight[t * 3] > 0.99f ? 1 : 0;
+        }
+        FlipResult fl = flip_edges_to_height(sub.geometry, sub.face_parent_id, sample, settings.flip_opts, locked);
+        sub.geometry       = std::move(fl.geometry);
+        sub.face_parent_id = std::move(fl.face_parent_id);
+        lap("align edges", sub.geometry, std::to_string(fl.flipped) + " flips");
+        if (!report("align edges", 1.0)) {
+            result.canceled = true;
+            return result;
+        }
     }
 
     TriSoup displaced = apply_displacement(sub.geometry, sample, settings.displace, bounds,
                                            [&](double f) { return report("displace", f); });
-    lap("displace", displaced.triangle_count());
+    lap("displace", displaced);
     if (!report("displace", 1.0)) {
         result.canceled = true;
         return result;
@@ -206,8 +238,12 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
     // 4. Decimate - export only. A bake needs the face-parent map, which a collapse destroys.
     std::vector<int> parent = std::move(sub.face_parent_id);
     if (mode == PipelineMode::Export) {
+        // Only when the mesh is actually over budget. Harvesting flat faces on a mesh that already fits
+        // cost several times the decimation itself and degraded the relief it was handed; it now only
+        // runs as part of a decimation that has to happen anyway. The repair pass below keys off the
+        // same decision (it runs only when decimation did), so an under-budget bake skips both.
         const bool needs_decimation = displaced.triangle_count() > settings.max_triangles;
-        if (needs_decimation || settings.harvest_flat) {
+        if (needs_decimation) {
             std::vector<uint8_t> locked;
             if (settings.preserve_untextured && !displaced.exclude_weight.empty()) {
                 locked.assign(displaced.triangle_count(), 0);
@@ -219,7 +255,7 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
                                           [&](double f) { return report("decimate", f); });
             result.locked_over_budget = dec.locked_over_budget;
             displaced                 = std::move(dec.geometry);
-            lap("decimate", displaced.triangle_count());
+            lap("decimate", displaced, "over budget, simplified");
             parent.clear(); // no longer meaningful
         }
         if (!report("decimate", 1.0)) {
@@ -229,15 +265,21 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
     }
 
     // 5. Flatten the bed-contact surface.
-    if (settings.clamp_below_plate || settings.displace.bottom_angle_limit > 0.f)
-        clamp_below_bottom(displaced, bounds.min.z());
-    if (settings.bottom_snap_tol > 0.0)
-        snap_bottom_to_flat(displaced, bounds.min.z(), settings.bottom_snap_tol);
+    {
+        const bool clamped = settings.clamp_below_plate || settings.displace.bottom_angle_limit > 0.f;
+        if (clamped)
+            clamp_below_bottom(displaced, bounds.min.z());
+        size_t snapped = 0;
+        if (settings.bottom_snap_tol > 0.0)
+            snapped = snap_bottom_to_flat(displaced, bounds.min.z(), settings.bottom_snap_tol);
+        if (clamped || settings.bottom_snap_tol > 0.0)
+            lap("bottom clamp + snap", displaced, std::to_string(snapped) + " triangles snapped flat");
+    }
 
     // 6. Close the T-junctions decimation left behind. Only meaningful when it ran.
     if (mode == PipelineMode::Export && parent.empty()) {
         displaced = resolve_t_junctions(displaced);
-        lap("repair", displaced.triangle_count());
+        lap("repair", displaced);
     }
 
     result.geometry       = std::move(displaced);

src/libslic3r/TextureBake/TextureBakePipeline.hpp

@@ -21,11 +21,17 @@
 #include "TextureBakeDisplace.hpp"
 #include "TextureBakeIndex.hpp"
 #include "TextureBakeRegularize.hpp"
+#include "TextureBakeFlip.hpp"
 #include "TextureBakeRelocate.hpp"
 #include "TextureBakeRepair.hpp"
 #include "TextureBakeSubdivide.hpp"
 
 namespace Slic3r {
+
+// Optional step-by-step capture; see TextureBakeDebug.hpp. A pointer, and forward declared, so the
+// pipeline header stays free of the mesh types the recorder converts into.
+class BakeStageRecorder;
+
 namespace TextureBake {
 
 enum class PipelineMode
@@ -53,10 +59,17 @@ struct PipelineSettings
     bool             relocate = false;
     RelocateSettings relocate_opts;
 
+    // Choose each quad's diagonal to follow the height field before displacing, so a step that crosses
+    // the grid at an angle comes out as a straight wall instead of a sawtooth. See TextureBakeFlip.hpp.
+    bool         flip_edges = true;
+    FlipSettings flip_opts;
+
     DisplaceSettings displace;
 
     // Export mode only.
     size_t max_triangles = 750'000;
+    // Keep removing zero-cost flat faces past the target. Only applies when decimation runs, i.e. when
+    // the displaced mesh is over max_triangles - an under-budget mesh is never decimated.
     bool   harvest_flat  = true;
     double harvest_tol   = DECIMATE_DEFAULT_HARVEST_TOL;
     // Lock the untextured region against both regularization and decimation.
@@ -86,10 +99,13 @@ struct PipelineResult
     bool             canceled           = false;
 };
 
+// `debug`, when given and enabled, receives the mesh after every stage that ran - which is the only
+// way to tell which stage a bad result came from, since each one rewrites the whole mesh.
 PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
                             const PipelineSettings &settings, const DisplaceBounds &bounds,
                             PipelineMode mode, const std::vector<uint8_t> &face_excluded = {},
-                            const PipelineProgressFn &on_progress = {});
+                            const PipelineProgressFn &on_progress = {},
+                            BakeStageRecorder *debug = nullptr);
 
 // Snap anything that ended below the model's original bottom back up to it.
 void clamp_below_bottom(TriSoup &geometry, float bottom_z);

src/libslic3r/TextureBake/TextureBakeRegularize.cpp

@@ -3,7 +3,9 @@
 #include <algorithm>
 #include <cmath>
 #include <limits>
-#include <numeric>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
 
 namespace Slic3r {
 namespace TextureBake {
@@ -286,31 +288,38 @@ RegularizeResult regularize_mesh(const TriSoup &geometry, const std::vector<int>
         return true;
     };
 
+    // Per-triangle thinness for a round's candidate scan; -1 marks "not a candidate". Scored in
+    // parallel, since the scan only reads the mesh, then gathered serially in index order.
+    std::vector<double> round_aspect(tri_count);
     for (int round = 0; round < opts.maxrounds; ++round) {
         // Rebuilt each round so earlier collapses inform the priorities.
-        std::vector<int>    cand;
-        std::vector<double> cand_aspect;
-        for (size_t t = 0; t < tri_count; ++t) {
-            if (tri_deleted[t])
-                continue;
-            const int a = corners[t * 3], b = corners[t * 3 + 1], c = corners[t * 3 + 2];
-            if (std::min({ sq_dist(a, b), sq_dist(b, c), sq_dist(c, a) }) <= 0.0)
-                continue;
-            const double aspect2 = tri_aspect_sq(t);
-            if (aspect2 < aspect_thr2)
-                continue;
-            cand.push_back(int(t));
-            cand_aspect.push_back(aspect2);
-        }
-        // Worst first; ties keep ascending order so the pass is deterministic.
-        std::vector<int> order(cand.size());
-        std::iota(order.begin(), order.end(), 0);
-        std::stable_sort(order.begin(), order.end(),
-                         [&](int x, int y) { return cand_aspect[size_t(x)] > cand_aspect[size_t(y)]; });
+        tbb::parallel_for(tbb::blocked_range<size_t>(0, tri_count), [&](const tbb::blocked_range<size_t> &r) {
+            for (size_t t = r.begin(); t < r.end(); ++t) {
+                round_aspect[t] = -1.0;
+                if (tri_deleted[t])
+                    continue;
+                const int a = corners[t * 3], b = corners[t * 3 + 1], c = corners[t * 3 + 2];
+                if (std::min({ sq_dist(a, b), sq_dist(b, c), sq_dist(c, a) }) <= 0.0)
+                    continue;
+                const double aspect2 = tri_aspect_sq(t);
+                if (aspect2 >= aspect_thr2)
+                    round_aspect[t] = aspect2;
+            }
+        });
+        // Worst first; ties keep ascending triangle order so the pass is deterministic. Sorting the
+        // (thinness, triangle) pairs themselves gives exactly the order the index sort with its
+        // indirect comparator did, without that comparator's extra lookup on every comparison.
+        std::vector<std::pair<double, int>> cand;
+        for (size_t t = 0; t < tri_count; ++t)
+            if (round_aspect[t] >= 0.0)
+                cand.emplace_back(round_aspect[t], int(t));
+        std::sort(cand.begin(), cand.end(), [](const std::pair<double, int> &x, const std::pair<double, int> &y) {
+            return x.first != y.first ? x.first > y.first : x.second < y.second;
+        });
 
         size_t round_collapses = 0;
-        for (const int oi : order) {
-            const size_t t = size_t(cand[size_t(oi)]);
+        for (const std::pair<double, int> &entry : cand) {
+            const size_t t = size_t(entry.second);
             if (tri_deleted[t])
                 continue;
             const int a = corners[t * 3], b = corners[t * 3 + 1], c = corners[t * 3 + 2];
@@ -332,8 +341,13 @@ RegularizeResult regularize_mesh(const TriSoup &geometry, const std::vector<int>
 
     // Drop deleted triangles and rebuild the soup.
     const bool have_weights = !geometry.exclude_weight.empty();
+    const size_t survivors = tri_count - size_t(std::count(tri_deleted.begin(), tri_deleted.end(), uint8_t(1)));
     std::vector<int> out_parent;
     TriSoup         &out = result.geometry;
+    out_parent.reserve(survivors);
+    out.pos.reserve(survivors * 3);
+    if (have_weights)
+        out.exclude_weight.reserve(survivors * 3);
     for (size_t t = 0; t < tri_count; ++t) {
         if (tri_deleted[t])
             continue;
@@ -351,7 +365,6 @@ RegularizeResult regularize_mesh(const TriSoup &geometry, const std::vector<int>
     // Rebuilt from the compacted geometry - the collapses moved vertices.
     out.nrm.assign(out.pos.size(), Vec3f::Zero());
     {
-        std::vector<Vec3d> accum(out.pos.size(), Vec3d::Zero());
         QuantizedPointMap  weld(WELD_GRID_GEOMETRY, out.pos.size());
         std::vector<int>   vid(out.pos.size());
         int                next = 0;

src/libslic3r/TextureBake/TextureBakeSubdivide.cpp

@@ -75,8 +75,11 @@ PassResult subdivide_pass(VertStore &verts, const std::vector<int> &indices, dou
     const size_t tri_count  = indices.size() / 3;
     const bool   have_canon = !verts.canon.empty();
 
-    QuantizedPointMap mid_cache(1.0, 1 << 16);
-    QuantizedPointMap split_edges(1.0, 1 << 16);
+    // Sized from the pass rather than grown from 64 k: a fine pass marks on the order of one edge per
+    // triangle, and growing to that by doubling rehashes the whole table a dozen times.
+    const size_t      expect = std::max<size_t>(size_t(1) << 16, tri_count);
+    QuantizedPointMap mid_cache(1.0, expect);
+    QuantizedPointMap split_edges(1.0, expect);
 
     // With canonical ids the key is the canonical *position* id, so split copies either side of a
     // sharp edge see one another's decision; without them the vertex index serves.
@@ -108,13 +111,19 @@ PassResult subdivide_pass(VertStore &verts, const std::vector<int> &indices, dou
         return out; // changed stays false: nothing left to refine
     }
 
-    // Step 1.5. Read-only against the finished mark set, so it reduces in parallel.
-    const size_t predicted = tbb::parallel_reduce(
+    // Step 1.5. Read-only against the finished mark set, so it runs in parallel - and it keeps each
+    // triangle's three marks (bit 0 = ab, 1 = bc, 2 = ca), so the serial rebuild below reads a byte
+    // instead of probing the hash map three more times per triangle.
+    std::vector<uint8_t> tri_marks(tri_count);
+    const size_t         predicted = tbb::parallel_reduce(
         tbb::blocked_range<size_t>(0, tri_count), size_t(0),
         [&](const tbb::blocked_range<size_t> &range, size_t acc) {
             for (size_t t = range.begin(); t < range.end(); ++t) {
-                const int a = indices[t * 3], b = indices[t * 3 + 1], c = indices[t * 3 + 2];
-                const int n = int(is_marked(a, b)) + int(is_marked(b, c)) + int(is_marked(c, a));
+                const int     a = indices[t * 3], b = indices[t * 3 + 1], c = indices[t * 3 + 2];
+                const uint8_t m = uint8_t(is_marked(a, b)) | uint8_t(is_marked(b, c) << 1) |
+                                  uint8_t(is_marked(c, a) << 2);
+                tri_marks[t]    = m;
+                const int n     = (m & 1) + ((m >> 1) & 1) + ((m >> 2) & 1);
                 acc += (n == 0) ? 1 : size_t(n + 1);
             }
             return acc;
@@ -150,7 +159,7 @@ PassResult subdivide_pass(VertStore &verts, const std::vector<int> &indices, dou
         const int     a    = indices[t * 3], b = indices[t * 3 + 1], c = indices[t * 3 + 2];
         const uint8_t excl = face_excluded.empty() ? uint8_t(0) : face_excluded[t];
         const int     pid  = face_parent_id.empty() ? 0 : face_parent_id[t];
-        const bool    s_ab = is_marked(a, b), s_bc = is_marked(b, c), s_ca = is_marked(c, a);
+        const bool    s_ab = (tri_marks[t] & 1) != 0, s_bc = (tri_marks[t] & 2) != 0, s_ca = (tri_marks[t] & 4) != 0;
         const int     n    = int(s_ab) + int(s_bc) + int(s_ca);
 
         if (n == 0) {
@@ -345,20 +354,23 @@ TriSoup to_non_indexed(const VertStore &verts, const std::vector<int> &indices,
     if (want_weights)
         out.exclude_weight.resize(tri_count * 3);
 
-    for (size_t t = 0; t < tri_count; ++t) {
-        // The per-face flag, not the interpolated weight: merging by maximum can push an *included*
-        // face's corners to 1 when it borders two excluded neighbours, wrongly excluding it.
-        const bool  have_face_flag = !face_excluded.empty();
-        const float face_w         = have_face_flag ? (face_excluded[t] ? 1.f : 0.f) : 0.f;
-        for (int v = 0; v < 3; ++v) {
-            const size_t vidx = size_t(indices[t * 3 + size_t(v)]);
-            out.pos[t * 3 + size_t(v)] = verts.pos[vidx].cast<float>();
-            out.nrm[t * 3 + size_t(v)] = verts.nrm[vidx].cast<float>();
-            if (want_weights)
-                out.exclude_weight[t * 3 + size_t(v)] =
-                    have_face_flag ? face_w : float(verts.wgt[vidx]);
+    // Each triangle writes only its own three slots, so this is a plain parallel gather.
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, tri_count), [&](const tbb::blocked_range<size_t> &r) {
+        for (size_t t = r.begin(); t < r.end(); ++t) {
+            // The per-face flag, not the interpolated weight: merging by maximum can push an *included*
+            // face's corners to 1 when it borders two excluded neighbours, wrongly excluding it.
+            const bool  have_face_flag = !face_excluded.empty();
+            const float face_w         = have_face_flag ? (face_excluded[t] ? 1.f : 0.f) : 0.f;
+            for (int v = 0; v < 3; ++v) {
+                const size_t vidx = size_t(indices[t * 3 + size_t(v)]);
+                out.pos[t * 3 + size_t(v)] = verts.pos[vidx].cast<float>();
+                out.nrm[t * 3 + size_t(v)] = verts.nrm[vidx].cast<float>();
+                if (want_weights)
+                    out.exclude_weight[t * 3 + size_t(v)] =
+                        have_face_flag ? face_w : float(verts.wgt[vidx]);
+            }
         }
-    }
+    });
     return out;
 }
 
@@ -375,9 +387,9 @@ SubdivideResult subdivide(const TriSoup &geometry, double max_edge_length,
     IndexedMesh        indexed   = fast ? to_indexed_fast(geometry) : to_indexed(geometry);
     QuantizedPointMap *canon_map = indexed.has_canon ? &indexed.pos_canon_map : nullptr;
 
-    std::vector<int>     current_indices  = indexed.indices;
-    std::vector<uint8_t> current_excluded = face_excluded;
     const size_t         initial_tris     = indexed.indices.size() / 3;
+    std::vector<int>     current_indices  = std::move(indexed.indices); // nothing reads it again
+    std::vector<uint8_t> current_excluded = face_excluded;
     std::vector<int>     current_parent(initial_tris);
     for (size_t i = 0; i < initial_tris; ++i)
         current_parent[i] = int(i);
@@ -400,13 +412,19 @@ SubdivideResult subdivide(const TriSoup &geometry, double max_edge_length,
 
         if (on_progress) {
             // Reported after the pass, so the value falls each iteration instead of lagging a step.
-            double max_edge_sq = 0.0;
-            for (size_t t = 0; t + 2 < current_indices.size(); t += 3) {
-                const int a = current_indices[t], b = current_indices[t + 1], c = current_indices[t + 2];
-                max_edge_sq = std::max({ max_edge_sq, edge_len_sq(indexed.verts, a, b),
-                                         edge_len_sq(indexed.verts, b, c),
+            // A max is order-independent, so the scan reduces in parallel to the same value.
+            const double max_edge_sq = tbb::parallel_reduce(
+                tbb::blocked_range<size_t>(0, current_indices.size() / 3), 0.0,
+                [&](const tbb::blocked_range<size_t> &r, double acc) {
+                    for (size_t f = r.begin(); f < r.end(); ++f) {
+                        const int a = current_indices[f * 3], b = current_indices[f * 3 + 1],
+                                  c = current_indices[f * 3 + 2];
+                        acc = std::max({ acc, edge_len_sq(indexed.verts, a, b), edge_len_sq(indexed.verts, b, c),
                                          edge_len_sq(indexed.verts, c, a) });
-            }
+                    }
+                    return acc;
+                },
+                [](double x, double y) { return std::max(x, y); });
             if (!on_progress(std::min(0.95, double(iter + 1) / SUBDIVIDE_MAX_ITERATIONS),
                              current_indices.size() / 3, std::sqrt(max_edge_sq)))
                 break; // whole passes only, so what we have is still crack-free