diff --git a/src/libslic3r/MeshBoolean.cpp b/src/libslic3r/MeshBoolean.cpp
index 5cd1968e22..4d9dda0c2e 100644
--- a/src/libslic3r/MeshBoolean.cpp
+++ b/src/libslic3r/MeshBoolean.cpp
@@ -261,7 +261,8 @@ indexed_triangle_set cgal_to_indexed_triangle_set(const CGALMesh &cgalmesh)
 // /////////////////////////////////////////////////////////////////////////////
 
 indexed_triangle_set remesh_isotropic(const indexed_triangle_set &mesh, double target_edge_length,
-                                      unsigned n_iterations, double sharp_angle_deg)
+                                      unsigned n_iterations, double sharp_angle_deg,
+                                      unsigned n_relaxation_steps)
 {
     if (mesh.indices.empty() || target_edge_length <= 0.0)
         return mesh;
@@ -305,6 +306,7 @@ indexed_triangle_set remesh_isotropic(const indexed_triangle_set &mesh, double t
 
         CGALProc::isotropic_remeshing(faces(cgal_mesh), target_edge_length, cgal_mesh,
                                       CGALParams::number_of_iterations(n_iterations)
+                                          .number_of_relaxation_steps(n_relaxation_steps)
                                           .edge_is_constrained_map(ecm)
                                           .protect_constraints(true));
     } catch (const std::exception &) {
diff --git a/src/libslic3r/MeshBoolean.hpp b/src/libslic3r/MeshBoolean.hpp
index efaf553420..c307d1528b 100644
--- a/src/libslic3r/MeshBoolean.hpp
+++ b/src/libslic3r/MeshBoolean.hpp
@@ -89,8 +89,13 @@ std::optional<std::vector<Vec2f>> parameterize_lscm(const indexed_triangle_set &
 // Edges whose dihedral angle exceeds `sharp_angle_deg`, and any open border, are held fixed so hard
 // features survive instead of being eroded by the relaxation pass; pass 0 to remesh everything.
 // Returns the input unchanged if remeshing fails (e.g. a non-manifold or self-intersecting input).
+// `n_relaxation_steps` is the number of tangential relaxation passes run inside each iteration. That
+// relaxation is what actually evens out the triangle distribution - splitting and collapsing alone
+// only bring edge *lengths* near the target, leaving the vertices wherever they happened to land. CGAL
+// defaults it to 1, which on a few iterations is not enough to look uniform.
 indexed_triangle_set remesh_isotropic(const indexed_triangle_set &mesh, double target_edge_length,
-                                      unsigned n_iterations = 3, double sharp_angle_deg = 40.0);
+                                      unsigned n_iterations = 3, double sharp_angle_deg = 40.0,
+                                      unsigned n_relaxation_steps = 1);
 }
 
 namespace mcut {
diff --git a/src/libslic3r/TextureBake/TextureBakeDecimate.cpp b/src/libslic3r/TextureBake/TextureBakeDecimate.cpp
index d4744f40a7..e591f7f191 100644
--- a/src/libslic3r/TextureBake/TextureBakeDecimate.cpp
+++ b/src/libslic3r/TextureBake/TextureBakeDecimate.cpp
@@ -6,6 +6,11 @@
 #include <limits>
 #include <queue>
 
+#include <boost/log/trivial.hpp>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+
 namespace Slic3r {
 namespace TextureBake {
 
@@ -74,13 +79,19 @@ 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.
+// 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.
+//
+// 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.
 struct HeapEntry
 {
     double   cost;
     int      v1, v2;
     uint32_t ver1, ver2;
-    Vec3d    p;
+    Vec3f    p;
     bool     operator>(const HeapEntry &o) const { return cost > o.cost; }
 };
 
@@ -143,16 +154,29 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
     }
 
     std::vector<Quadric> quadrics(vert_count);
-    for (size_t f = 0; f < face_count; ++f) {
-        const int a = faces[f * 3], b = faces[f * 3 + 1], c = faces[f * 3 + 2];
-        if (a < 0)
-            continue;
-        const Vec3d nrm = face_normal_unit(pos, a, b, c);
-        if (nrm.isZero())
-            continue;
-        const double d = -nrm.dot(pos[size_t(a)]);
-        for (const int v : { a, b, c })
-            quadrics[size_t(v)].add_plane(nrm.x(), nrm.y(), nrm.z(), d);
+    {
+        // The plane per face is independent; accumulating it into the three incident vertices is not,
+        // so only the first half is parallel.
+        std::vector<Vec4d> planes(face_count, Vec4d::Zero());
+        tbb::parallel_for(tbb::blocked_range<size_t>(0, face_count),
+                          [&](const tbb::blocked_range<size_t> &range) {
+                              for (size_t f = range.begin(); f < range.end(); ++f) {
+                                  const int a = faces[f * 3], b = faces[f * 3 + 1], c = faces[f * 3 + 2];
+                                  if (a < 0)
+                                      continue;
+                                  const Vec3d nrm = face_normal_unit(pos, a, b, c);
+                                  if (nrm.isZero())
+                                      continue;
+                                  planes[f] = Vec4d(nrm.x(), nrm.y(), nrm.z(), -nrm.dot(pos[size_t(a)]));
+                              }
+                          });
+        for (size_t f = 0; f < face_count; ++f) {
+            const Vec4d &pl = planes[f];
+            if (pl.head<3>().isZero())
+                continue;
+            for (int k = 0; k < 3; ++k)
+                quadrics[size_t(faces[f * 3 + size_t(k)])].add_plane(pl.x(), pl.y(), pl.z(), pl.w());
+        }
     }
 
     // Two penalty planes per endpoint on a sharp interior edge, each perpendicular to one adjacent
@@ -250,7 +274,22 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
     uint32_t              epoch = 1, lk_epoch = 1;
     size_t                active_faces = face_count;
 
-    std::priority_queue<HeapEntry, std::vector<HeapEntry>, std::greater<HeapEntry>> heap;
+    // A plain vector driven by the heap algorithms, so the capacity can be reserved. Lazy deletion
+    // means roughly one entry per edge plus one per re-push after each collapse; the reserve below is
+    // sized from the edge count and simply grows if a mesh needs more.
+    std::vector<HeapEntry> heap;
+    heap.reserve(std::min<size_t>(face_count * 3, size_t(1) << 24));
+    const auto heap_push = [&](HeapEntry e) {
+        heap.push_back(e);
+        std::push_heap(heap.begin(), heap.end(), std::greater<HeapEntry>());
+    };
+    const auto heap_pop = [&]() {
+        std::pop_heap(heap.begin(), heap.end(), std::greater<HeapEntry>());
+        const HeapEntry e = heap.back();
+        heap.pop_back();
+        return e;
+    };
+    size_t pops = 0, stale_pops = 0;
 
     const auto push_edge = [&](int v1, int v2) {
         Vec3d p;
@@ -269,8 +308,8 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
         }
         // 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 });
+        heap_push({ eval_sum(quadrics, v1, v2, p) + len2 * 1e-8, v1, v2, version[size_t(v1)],
+                    version[size_t(v2)], p.cast<float>() });
     };
 
     {
@@ -384,27 +423,32 @@ DecimateResult decimate(const TriSoup &geometry, size_t target_triangles, bool h
             reached_target = true;
         }
 
-        const HeapEntry top = heap.top();
-        heap.pop();
+        const HeapEntry top = heap_pop();
+        ++pops;
         // The popped entry is the cheapest left, so exceeding the tolerance ends the run.
         if (reached_target && top.cost > harvest_ceil)
             break;
 
         const int v1 = top.v1, v2 = top.v2;
-        if (!active[size_t(v1)] || !active[size_t(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 (version[size_t(v1)] != top.ver1 || version[size_t(v2)] != top.ver2) {
+            ++stale_pops;
             continue;
+        }
         if (shared_face_count(v1, v2) < 2)
             continue;
         lk_epoch += 2; // +2 so ep and ep+1 cannot collide with the next call
         if (has_link_violation(v1, v2, lk_epoch))
             continue;
-        if (check_flipped(v1, v2, top.p) || check_flipped(v2, v1, top.p))
+        const Vec3d target = top.p.cast<double>();
+        if (check_flipped(v1, v2, target) || check_flipped(v2, v1, target))
             continue;
 
         // v1 survives at the new position, v2 goes.
-        pos[size_t(v1)] = top.p;
+        pos[size_t(v1)] = target;
         quadrics[size_t(v1)] += quadrics[size_t(v2)];
         ++version[size_t(v1)];
 
@@ -459,6 +503,9 @@ 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;
+
     // Rebuild from the surviving faces, with per-face normals.
     TriSoup &out = result.geometry;
     for (size_t f = 0; f < face_count; ++f) {
diff --git a/src/libslic3r/TextureBake/TextureBakeDisplace.cpp b/src/libslic3r/TextureBake/TextureBakeDisplace.cpp
index ef417a7bef..540cb7dda8 100644
--- a/src/libslic3r/TextureBake/TextureBakeDisplace.cpp
+++ b/src/libslic3r/TextureBake/TextureBakeDisplace.cpp
@@ -4,6 +4,9 @@
 #include <cmath>
 #include <limits>
 
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+
 namespace Slic3r {
 namespace TextureBake {
 
@@ -204,20 +207,30 @@ TriSoup apply_displacement(const TriSoup &geometry, const HeightSampleFn &sample
         }
     }
 
-    // Pass 2: one sample per unique position.
-    std::vector<double>  grey(unique_count, 0.0);
-    std::vector<uint8_t> grey_set(unique_count, 0);
-    for (size_t i = 0; i < count; ++i) {
-        const size_t vid = size_t(vertex_id[i]);
-        if (grey_set[vid])
-            continue;
-        grey_set[vid] = 1;
-        grey[vid] = double(sample(geometry.pos[i], smooth_nrm[vid].cast<float>(),
-                                  blend_nrm[vid].cast<float>()));
-    }
+    // Pass 2: one sample per unique position. A representative corner is picked first so the sampling
+    // itself is a flat parallel loop - it is a texture fetch plus projection maths per layer, and by
+    // far the most expensive thing in this stage.
+    std::vector<double> grey(unique_count, 0.0);
+    std::vector<int>    representative(unique_count, -1);
+    for (size_t i = 0; i < count; ++i)
+        if (representative[size_t(vertex_id[i])] < 0)
+            representative[size_t(vertex_id[i])] = int(i);
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, unique_count),
+                      [&](const tbb::blocked_range<size_t> &range) {
+                          for (size_t vid = range.begin(); vid < range.end(); ++vid) {
+                              const int rep = representative[vid];
+                              if (rep < 0)
+                                  continue;
+                              grey[vid] = double(sample(geometry.pos[size_t(rep)],
+                                                        smooth_nrm[vid].cast<float>(),
+                                                        blend_nrm[vid].cast<float>()));
+                          }
+                      });
 
-    // Pass 3: move every copy of a position by the identical vector.
-    for (size_t i = 0; i < count; ++i) {
+    // Pass 3: move every copy of a position by the identical vector. Each iteration writes only its
+    // own output slot, so the loop is independent per corner.
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, count), [&](const tbb::blocked_range<size_t> &range) {
+    for (size_t i = range.begin(); i < range.end(); ++i) {
         const Vec3f &p   = geometry.pos[i];
         const size_t vid = size_t(vertex_id[i]);
 
@@ -249,17 +262,18 @@ TriSoup apply_displacement(const TriSoup &geometry, const HeightSampleFn &sample
             moved.z() = double(p.z());
 
         out.pos[i] = moved.cast<float>();
-
-        if (on_progress && (i % 5000) == 0 && !on_progress(double(i) / double(count)))
-            return geometry; // cancelled: hand back the input untouched
     }
+    });
 
     // Per-face, not averaged across shared positions: averaging can flip an excluded face's normal
     // when its neighbours moved outward.
-    for (size_t t = 0; t + 2 < count; t += 3) {
-        const Vec3f n = (out.pos[t + 1] - out.pos[t]).cross(out.pos[t + 2] - out.pos[t]).normalized();
-        out.nrm[t] = out.nrm[t + 1] = out.nrm[t + 2] = n;
-    }
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, count / 3), [&](const tbb::blocked_range<size_t> &r) {
+        for (size_t f = r.begin(); f < r.end(); ++f) {
+            const size_t t = f * 3;
+            const Vec3f  n = (out.pos[t + 1] - out.pos[t]).cross(out.pos[t + 2] - out.pos[t]).normalized();
+            out.nrm[t] = out.nrm[t + 1] = out.nrm[t + 2] = n;
+        }
+    });
     out.exclude_weight = geometry.exclude_weight;
     return out;
 }
diff --git a/src/libslic3r/TextureBake/TextureBakePipeline.cpp b/src/libslic3r/TextureBake/TextureBakePipeline.cpp
index a86fda0e0d..510c44e844 100644
--- a/src/libslic3r/TextureBake/TextureBakePipeline.cpp
+++ b/src/libslic3r/TextureBake/TextureBakePipeline.cpp
@@ -1,8 +1,11 @@
 #include "TextureBakePipeline.hpp"
 
 #include <algorithm>
+#include <chrono>
 #include <cmath>
 
+#include <boost/log/trivial.hpp>
+
 namespace Slic3r {
 namespace TextureBake {
 
@@ -111,6 +114,15 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
     const auto     report = [&](const char *stage, double f) {
         return !on_progress || on_progress(stage, f);
     };
+    // Per-stage wall time. The stages differ in cost by orders of magnitude depending on the model, so
+    // 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) {
+        const double ms = std::chrono::duration<double, std::milli>(clock_now() - t_stage).count();
+        BOOST_LOG_TRIVIAL(info) << "TextureBake " << stage << ": " << ms << " ms, " << tris << " tris";
+        t_stage = clock_now();
+    };
 
     if (input.empty() || !sample) {
         result.geometry = input;
@@ -122,6 +134,7 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
         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());
     if (!report("subdivide", 1.0)) {
         result.canceled = true;
         return result;
@@ -134,6 +147,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());
         if (!report("regularize", 1.0)) {
             result.canceled = true;
             return result;
@@ -159,15 +173,31 @@ 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());
         } else {
             sub.geometry       = std::move(reg.geometry);
             sub.face_parent_id = std::move(reg.face_parent_id);
         }
     }
 
-    // 3. Displace.
+    // 3. Align the mesh to the height field's edges, then displace.
+    if (settings.relocate) {
+        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;
+        }
+        RelocateResult rel = relocate_to_contours(sub.geometry, sample, settings.relocate_opts, locked);
+        BOOST_LOG_TRIVIAL(info) << "TextureBake relocate: moved=" << rel.moved
+                                << " rejected=" << rel.rejected;
+        sub.geometry = std::move(rel.geometry);
+        lap("relocate", sub.geometry.triangle_count());
+    }
+
     TriSoup displaced = apply_displacement(sub.geometry, sample, settings.displace, bounds,
                                            [&](double f) { return report("displace", f); });
+    lap("displace", displaced.triangle_count());
     if (!report("displace", 1.0)) {
         result.canceled = true;
         return result;
@@ -189,6 +219,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());
             parent.clear(); // no longer meaningful
         }
         if (!report("decimate", 1.0)) {
@@ -198,14 +229,16 @@ PipelineResult run_pipeline(const TriSoup &input, const HeightSampleFn &sample,
     }
 
     // 5. Flatten the bed-contact surface.
-    if (settings.displace.bottom_angle_limit > 0.f)
+    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);
 
     // 6. Close the T-junctions decimation left behind. Only meaningful when it ran.
-    if (mode == PipelineMode::Export && parent.empty())
+    if (mode == PipelineMode::Export && parent.empty()) {
         displaced = resolve_t_junctions(displaced);
+        lap("repair", displaced.triangle_count());
+    }
 
     result.geometry       = std::move(displaced);
     result.face_parent_id = std::move(parent);
diff --git a/src/libslic3r/TextureBake/TextureBakePipeline.hpp b/src/libslic3r/TextureBake/TextureBakePipeline.hpp
index a9a7db1957..a1cee46c1c 100644
--- a/src/libslic3r/TextureBake/TextureBakePipeline.hpp
+++ b/src/libslic3r/TextureBake/TextureBakePipeline.hpp
@@ -2,7 +2,7 @@
 
 // The bake pipeline:
 //
-//   subdivide -> [regularize -> re-subdivide] -> displace -> [decimate]
+//   subdivide -> [regularize -> re-subdivide] -> [relocate] -> displace -> [decimate]
 //             -> bottom clamp -> bottom snap -> [resolve T-junctions]
 //
 // Regularization sits between two subdivisions on purpose: it dissolves the slivers refinement
@@ -21,6 +21,7 @@
 #include "TextureBakeDisplace.hpp"
 #include "TextureBakeIndex.hpp"
 #include "TextureBakeRegularize.hpp"
+#include "TextureBakeRelocate.hpp"
 #include "TextureBakeRepair.hpp"
 #include "TextureBakeSubdivide.hpp"
 
@@ -47,6 +48,11 @@ struct PipelineSettings
     // re-refining what it just merged.
     double regularize_second_pass_mul = 1.1;
 
+    // Slide vertices onto the height map's own edges before displacing, so a step lands on a mesh
+    // edge instead of being quantised to wherever the grid fell.
+    bool             relocate = false;
+    RelocateSettings relocate_opts;
+
     DisplaceSettings displace;
 
     // Export mode only.
@@ -56,6 +62,10 @@ struct PipelineSettings
     // Lock the untextured region against both regularization and decimation.
     bool preserve_untextured = true;
 
+    // Push anything that displaced below the plate back up to it. Downward movement is otherwise left
+    // alone, so relief on the underside is kept - only what would sink through the plate is stopped.
+    bool clamp_below_plate = false;
+
     // Snap vertices within this of the bottom plane onto it. 0 disables.
     double bottom_snap_tol = 0.1;
 
diff --git a/src/libslic3r/TextureBake/TextureBakeRelocate.cpp b/src/libslic3r/TextureBake/TextureBakeRelocate.cpp
new file mode 100644
index 0000000000..5ea33b353b
--- /dev/null
+++ b/src/libslic3r/TextureBake/TextureBakeRelocate.cpp
@@ -0,0 +1,205 @@
+#include "TextureBakeRelocate.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+
+namespace Slic3r {
+namespace TextureBake {
+
+namespace {
+
+// Any two unit vectors orthogonal to n. Which two does not matter - the gradient is expressed in this
+// basis and converted straight back, so the result is basis independent.
+void tangent_basis(const Vec3f &n, Vec3f &t1, Vec3f &t2)
+{
+    const Vec3f a = (std::abs(n.x()) < 0.9f) ? Vec3f(1.f, 0.f, 0.f) : Vec3f(0.f, 1.f, 0.f);
+    t1 = n.cross(a).normalized();
+    t2 = n.cross(t1).normalized();
+}
+
+} // namespace
+
+RelocateResult relocate_to_contours(const TriSoup &geometry, const HeightSampleFn &sample,
+                                    const RelocateSettings &settings, const std::vector<uint8_t> &locked)
+{
+    RelocateResult result;
+    result.geometry     = geometry;
+    const size_t count  = geometry.pos.size();
+    const size_t tri_ct = count / 3;
+    if (count == 0 || !sample || settings.iterations <= 0)
+        return result;
+
+    // Weld, so every copy of a position moves together and the mesh cannot come apart.
+    QuantizedPointMap  weld(WELD_GRID_GEOMETRY, std::min(count, size_t(1) << 22));
+    std::vector<int>   vid(count);
+    std::vector<Vec3f> pos;
+    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]);
+    }
+    const size_t nv = pos.size();
+
+    // Incident corners per position, CSR style, plus the mean incident edge length that sets the scale
+    // for both the finite difference and the move limit.
+    std::vector<uint32_t> start(nv + 1, 0);
+    for (size_t i = 0; i < count; ++i)
+        ++start[size_t(vid[i]) + 1];
+    for (size_t v = 0; v < nv; ++v)
+        start[v + 1] += start[v];
+    std::vector<uint32_t> inc(count), cursor(nv, 0);
+    for (size_t i = 0; i < count; ++i)
+        inc[start[size_t(vid[i])] + cursor[size_t(vid[i])]++] = uint32_t(i);
+
+    std::vector<uint8_t> frozen(nv, 0);
+    if (!locked.empty())
+        for (size_t t = 0; t < tri_ct && t < locked.size(); ++t)
+            if (locked[t])
+                for (int k = 0; k < 3; ++k)
+                    frozen[size_t(vid[t * 3 + size_t(k)])] = 1;
+
+    std::vector<float> edge_len(nv, 0.f), normal_len(nv, 0.f);
+    std::vector<Vec3f> nrm(nv, Vec3f::Zero());
+    const auto rebuild_frames = [&]() {
+        std::fill(nrm.begin(), nrm.end(), Vec3f::Zero());
+        std::fill(edge_len.begin(), edge_len.end(), 0.f);
+        std::vector<uint32_t> deg(nv, 0);
+        for (size_t t = 0; t < tri_ct; ++t) {
+            const int a = vid[t * 3], b = vid[t * 3 + 1], c = vid[t * 3 + 2];
+            const Vec3f fn = (pos[size_t(b)] - pos[size_t(a)]).cross(pos[size_t(c)] - pos[size_t(a)]);
+            for (int k = 0; k < 3; ++k) {
+                const int u = vid[t * 3 + size_t(k)], w = vid[t * 3 + size_t((k + 1) % 3)];
+                nrm[size_t(u)] += fn;
+                edge_len[size_t(u)] += (pos[size_t(w)] - pos[size_t(u)]).norm();
+                ++deg[size_t(u)];
+            }
+        }
+        for (size_t v = 0; v < nv; ++v) {
+            const float l = nrm[v].norm();
+            nrm[v]      = (l > 0.f) ? Vec3f(nrm[v] / l) : Vec3f(0.f, 0.f, 1.f);
+            edge_len[v] = deg[v] > 0 ? edge_len[v] / float(deg[v]) : 0.f;
+        }
+    };
+    rebuild_frames();
+
+    // The level to snap onto, taken from the height actually present on this patch rather than assumed.
+    // A texture that never reaches full black or white would otherwise be measured against a range it
+    // does not occupy.
+    double h_lo = std::numeric_limits<double>::max(), h_hi = -h_lo;
+    {
+        std::vector<float> h0(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)
+                h0[v] = sample(pos[v], nrm[v], nrm[v]);
+        });
+        for (const float h : h0) {
+            h_lo = std::min(h_lo, double(h));
+            h_hi = std::max(h_hi, double(h));
+        }
+    }
+    const double h_range = h_hi - h_lo;
+    if (!(h_range > 0.0))
+        return result; // a flat height field has no contour to snap to
+    const double target   = h_lo + h_range * settings.contour_level;
+    // A gradient is worth acting on when the height changes by this much across one edge length.
+    const double min_grad = h_range * settings.min_gradient_fraction;
+
+    std::vector<uint8_t> ever_moved(nv, 0);
+
+    for (int iter = 0; iter < settings.iterations; ++iter) {
+        std::vector<Vec3f> proposal(nv);
+        std::vector<uint8_t> want(nv, 0);
+
+        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) {
+                if (frozen[v] || edge_len[v] <= 0.f)
+                    continue;
+                const Vec3f n = nrm[v];
+                Vec3f       t1, t2;
+                tangent_basis(n, t1, t2);
+                const float eps = edge_len[v] * float(settings.gradient_step_fraction);
+                if (eps <= 0.f)
+                    continue;
+
+                // Central differences in the tangent plane. Sampling the field itself, not the mesh,
+                // so the gradient is the image's, at whatever resolution the mesh happens to have.
+                const double h  = double(sample(pos[v], n, n));
+                const double gx = (double(sample(pos[v] + t1 * eps, n, n)) -
+                                   double(sample(pos[v] - t1 * eps, n, n))) / (2.0 * double(eps));
+                const double gy = (double(sample(pos[v] + t2 * eps, n, n)) -
+                                   double(sample(pos[v] - t2 * eps, n, n))) / (2.0 * double(eps));
+                const double g2 = gx * gx + gy * gy;
+                if (g2 <= 0.0)
+                    continue;
+                // Scale-free test: how much the height changes across one edge, versus the patch range.
+                if (std::sqrt(g2) * double(edge_len[v]) < min_grad)
+                    continue;
+
+                // Newton step onto the level set h = target, expressed back in 3D.
+                const double s   = -(h - target) / g2;
+                Vec3f        d   = t1 * float(s * gx) + t2 * float(s * gy);
+                const float  cap = edge_len[v] * float(settings.max_move_fraction);
+                const float  len = d.norm();
+                if (len <= 0.f)
+                    continue;
+                if (len > cap)
+                    d *= cap / len;
+                proposal[v] = pos[v] + d;
+                want[v]     = 1;
+            }
+        });
+
+        // Apply one at a time: a move is only valid against the neighbourhood as it stands, and two
+        // adjacent vertices moving together can invert a triangle neither would have on its own.
+        size_t applied = 0;
+        for (size_t v = 0; v < nv; ++v) {
+            if (!want[v])
+                continue;
+            const Vec3f old = pos[v];
+            pos[v]          = proposal[v];
+            bool ok = true;
+            for (uint32_t k = start[v]; k < start[v + 1] && ok; ++k) {
+                const size_t t = size_t(inc[k]) / 3;
+                const Vec3f &a = pos[size_t(vid[t * 3])];
+                const Vec3f  n2 = (pos[size_t(vid[t * 3 + 1])] - a).cross(pos[size_t(vid[t * 3 + 2])] - a);
+                // Compared against the frame this vertex carried before the move: a triangle that
+                // flips or collapses means the move crossed a neighbour.
+                if (n2.squaredNorm() <= 0.f || n2.normalized().dot(nrm[v]) < 0.f)
+                    ok = false;
+            }
+            if (ok) {
+                ++applied;
+                ever_moved[v] = 1;
+            } else {
+                pos[v] = old;
+                ++result.rejected;
+            }
+        }
+        if (applied == 0)
+            break;
+        rebuild_frames();
+    }
+
+    for (size_t v = 0; v < nv; ++v)
+        if (ever_moved[v])
+            ++result.moved;
+
+    // Write the relocated positions back to every copy, and rebuild the per-face normals.
+    for (size_t i = 0; i < count; ++i)
+        result.geometry.pos[i] = pos[size_t(vid[i])];
+    for (size_t t = 0; t < tri_ct; ++t) {
+        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
diff --git a/src/libslic3r/TextureBake/TextureBakeRelocate.hpp b/src/libslic3r/TextureBake/TextureBakeRelocate.hpp
new file mode 100644
index 0000000000..cd674d03ae
--- /dev/null
+++ b/src/libslic3r/TextureBake/TextureBakeRelocate.hpp
@@ -0,0 +1,72 @@
+#pragma once
+
+// Tangential relocation: slide vertices along the surface so triangle edges land on the height map's
+// own edges, before any displacement happens.
+//
+// Displacement moves vertices along the normal only, so a step in the height map - the wall of a
+// mortar groove, the rim of an embossed shape - is reproduced wherever the triangle grid happens to
+// fall, as a staircase quantised to triangle boundaries. Refining further only makes the steps
+// smaller; it never straightens them, because the edge in the image still does not coincide with any
+// edge in the mesh.
+//
+// This pass fixes the cause rather than the symptom. For a vertex sitting near a step it takes a
+// Newton step onto the contour: with h the sampled height and g its tangential gradient, the move
+//
+//     d = -(h - target) * g / |g|^2
+//
+// lands on the level set h = target to first order. The ring of vertices nearest each step therefore
+// snaps onto it, the triangle edges between them follow the step, and the displaced result has a
+// straight wall instead of a sawtooth.
+//
+// Vertices in flat regions have no gradient to speak of and are left alone, so the pass costs nothing
+// where there is nothing to align.
+
+#include <cstdint>
+#include <functional>
+#include <vector>
+
+#include "TextureBakeDisplace.hpp"
+#include "TextureBakeIndex.hpp"
+
+namespace Slic3r {
+namespace TextureBake {
+
+struct RelocateSettings
+{
+    // Passes. Each is a Newton step, so a couple converge for vertices that start reasonably close;
+    // more mainly helps ones that begin further away.
+    int iterations = 3;
+
+    // How far a vertex may move in one pass, as a fraction of the mean length of its incident edges.
+    // Below a half it cannot pass a neighbour, which is what keeps the triangulation valid without
+    // needing a full topological check.
+    double max_move_fraction = 0.35;
+
+    // A vertex is only pulled when the height varies enough across its own footprint to mean
+    // something - as a fraction of the height range over the whole patch. Below this the gradient is
+    // noise, and chasing it would scramble flat regions.
+    double min_gradient_fraction = 0.05;
+
+    // The level to snap onto, as a position in the sampled height range: 0.5 is midway between the
+    // lowest and highest point of the relief, which is where the wall of a step is steepest.
+    double contour_level = 0.5;
+
+    // Sampling offset for the finite-difference gradient, as a fraction of the local edge length.
+    double gradient_step_fraction = 0.25;
+};
+
+struct RelocateResult
+{
+    TriSoup geometry;
+    size_t  moved = 0;      // positions that were relocated at least once
+    size_t  rejected = 0;   // moves refused because a triangle would have inverted
+};
+
+// `locked`, when non-empty, has one entry per input triangle; a vertex touching a locked triangle is
+// never moved, so an excluded region keeps its exact vertex positions.
+RelocateResult relocate_to_contours(const TriSoup &geometry, const HeightSampleFn &sample,
+                                    const RelocateSettings &settings = {},
+                                    const std::vector<uint8_t> &locked = {});
+
+} // namespace TextureBake
+} // namespace Slic3r
diff --git a/src/libslic3r/TextureBake/TextureBakeSubdivide.cpp b/src/libslic3r/TextureBake/TextureBakeSubdivide.cpp
index 7a0082e257..0a5932e962 100644
--- a/src/libslic3r/TextureBake/TextureBakeSubdivide.cpp
+++ b/src/libslic3r/TextureBake/TextureBakeSubdivide.cpp
@@ -4,6 +4,10 @@
 #include <cmath>
 #include <unordered_map>
 
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_for.h>
+#include <tbb/parallel_reduce.h>
+
 namespace Slic3r {
 namespace TextureBake {
 
@@ -104,13 +108,18 @@ PassResult subdivide_pass(VertStore &verts, const std::vector<int> &indices, dou
         return out; // changed stays false: nothing left to refine
     }
 
-    // Step 1.5.
-    size_t predicted = 0;
-    for (size_t t = 0; t < tri_count; ++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));
-        predicted += (n == 0) ? 1 : size_t(n + 1);
-    }
+    // Step 1.5. Read-only against the finished mark set, so it reduces in parallel.
+    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));
+                acc += (n == 0) ? 1 : size_t(n + 1);
+            }
+            return acc;
+        },
+        std::plus<size_t>());
     if (predicted > size_t(safety_cap)) {
         out.indices        = indices;
         out.face_excluded  = face_excluded;
@@ -255,18 +264,20 @@ IndexedMesh to_indexed(const TriSoup &geometry)
 
     // Per-face normals: unit for the angle test, raw for the area-weighted accumulation.
     std::vector<Vec3d> face_unit(n), face_raw(n);
-    for (size_t t = 0; t + 2 < n; t += 3) {
-        const Vec3d  a   = geometry.pos[t].cast<double>();
-        const Vec3d  b   = geometry.pos[t + 1].cast<double>();
-        const Vec3d  c   = geometry.pos[t + 2].cast<double>();
-        const Vec3d  r   = (b - a).cross(c - a);
-        const double len = r.norm();
-        const Vec3d  u   = (len > 0.0) ? Vec3d(r / len) : Vec3d(0.0, 0.0, 1.0);
-        for (int v = 0; v < 3; ++v) {
-            face_unit[t + size_t(v)] = u;
-            face_raw[t + size_t(v)]  = r;
+    tbb::parallel_for(tbb::blocked_range<size_t>(0, n / 3), [&](const tbb::blocked_range<size_t> &range) {
+        for (size_t f = range.begin(); f < range.end(); ++f) {
+            const size_t t   = f * 3;
+            const Vec3d  a   = geometry.pos[t].cast<double>();
+            const Vec3d  r   = (geometry.pos[t + 1].cast<double>() - a).cross(
+                                 geometry.pos[t + 2].cast<double>() - a);
+            const double len = r.norm();
+            const Vec3d  u   = (len > 0.0) ? Vec3d(r / len) : Vec3d(0.0, 0.0, 1.0);
+            for (int v = 0; v < 3; ++v) {
+                face_unit[t + size_t(v)] = u;
+                face_raw[t + size_t(v)]  = r;
+            }
         }
-    }
+    });
 
     out.indices.resize(n);
     out.pos_canon_map = QuantizedPointMap(WELD_GRID_GEOMETRY, std::min(n, size_t(1) << 22));
diff --git a/src/libslic3r/TextureDisplacement.cpp b/src/libslic3r/TextureDisplacement.cpp
index 4fcf00b462..a129b503ff 100644
--- a/src/libslic3r/TextureDisplacement.cpp
+++ b/src/libslic3r/TextureDisplacement.cpp
@@ -1405,12 +1405,17 @@ indexed_triangle_set build_texture_displacement_v2(const indexed_triangle_set
     settings.regularize    = options.v2_regularize;
     settings.max_triangles = size_t(std::max(0, options.v2_max_triangles_k)) * 1000;
     settings.preserve_untextured = true;
+    settings.clamp_below_plate   = options.v2_clamp_below_plate;
+    settings.relocate            = options.v2_relocate;
     // 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;
 
     TextureBake::DisplaceBounds bounds;
     bounds.min = bounds.max = mesh.vertices.empty() ? Vec3f::Zero() : mesh.vertices.front();
diff --git a/src/libslic3r/TextureDisplacement.hpp b/src/libslic3r/TextureDisplacement.hpp
index c187c3be48..63cb443420 100644
--- a/src/libslic3r/TextureDisplacement.hpp
+++ b/src/libslic3r/TextureDisplacement.hpp
@@ -373,6 +373,13 @@ struct TextureDisplacementOptions
     float v2_refine_mm       = 0.3f;
     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;
+    // 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;
 
     // 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.
@@ -390,8 +397,8 @@ 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, color_mix_enabled, mix_mode,
-           color_despeckle);
+           pipeline_v2, v2_refine_mm, v2_regularize, v2_max_triangles_k, v2_clamp_below_plate,
+           v2_relocate, color_mix_enabled, mix_mode, color_despeckle);
         color_mix_mode = ColorMixMode(mix_mode);
     }
 };
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoTextureDisplacement.cpp b/src/slic3r/GUI/Gizmos/GLGizmoTextureDisplacement.cpp
index 8de71180f0..70e4b83e7f 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoTextureDisplacement.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoTextureDisplacement.cpp
@@ -3828,8 +3828,13 @@ bool GLGizmoTextureDisplacement::plan_remesh(const indexed_triangle_set &src, fl
         target_edge_mm = float(budget_edge);
     }
 
-    indexed_triangle_set remeshed =
-        MeshBoolean::cgal::remesh_isotropic(src, double(target_edge_mm), 3, double(sharp_angle_deg));
+    // Five iterations with three relaxation passes each, rather than CGAL's three-and-one. Splitting
+    // and collapsing bring edge lengths near the target but leave the vertices where they fell, so it
+    // is the relaxation count that decides how even the result looks - and three passes in total was
+    // nowhere near enough. Fifteen costs proportionally more, but only on an explicit Remesh.
+    indexed_triangle_set remeshed = MeshBoolean::cgal::remesh_isotropic(
+        src, double(target_edge_mm), /* iterations */ 5, double(sharp_angle_deg),
+        /* relaxation steps */ 3);
     // remesh_isotropic() signals failure by handing the input straight back, so compare against it
     // structurally. Vertex count alone is not enough: a remesh that only redistributes triangles at
     // roughly the current density legitimately lands on the same count, and treating that as failure
@@ -5405,6 +5410,24 @@ void GLGizmoTextureDisplacement::on_render_input_window(float x, float y, float
                                   "simplification off, which is worth comparing on its own."),
                               m_imgui->scaled(20.f));
         ImGui::PopItemWidth();
+        m_preview_params_dirty |= ImGui::Checkbox(_u8L("Keep above build plate").c_str(),
+                                                  &opts.v2_clamp_below_plate);
+        if (ImGui::IsItemHovered())
+            m_imgui->tooltip(_u8L("Push anything the displacement drove below the build plate back up to it. "
+                                  "Only geometry that would end up under the model's own bottom is moved - "
+                                  "relief that goes downward but stays above the plate is left as it is."),
+                              m_imgui->scaled(20.f));
+
+        m_preview_params_dirty |= ImGui::Checkbox(_u8L("Align mesh to texture edges").c_str(),
+                                                  &opts.v2_relocate);
+        if (ImGui::IsItemHovered())
+            m_imgui->tooltip(_u8L("Slide vertices sideways onto the edges in the texture before displacing them. "
+                                  "Displacement can only move vertices up and down, so without this a sharp step "
+                                  "in the image lands wherever the triangles happen to be and comes out as a "
+                                  "staircase. Moving the vertices onto the step first gives a straight wall at the "
+                                  "same triangle count."),
+                              m_imgui->scaled(20.f));
+
         m_preview_params_dirty |= ImGui::Checkbox(_u8L("Clean up slivers").c_str(), &opts.v2_regularize);
         if (ImGui::IsItemHovered())
             m_imgui->tooltip(_u8L("Collapse the thin triangles refinement inherits from the model's own "