Price a mesh cage by the control points that were dragged

benchmarks/bench_main.cpp

@@ -30,6 +30,7 @@
 #include "clay/field/relax.h"
 #include "clay/kernel/field.h"
 #include "clay/mesh/bvh.h"
+#include "clay/mesh/lattice.h"
 #include "clay/mesh/marching.h"
 #include "clay/mesh/surface_nets.h"
 #include "clay/mesh/to_field.h"
@@ -5287,6 +5288,39 @@ BENCHMARK(BM_MeshStampBoundaryAutomask)
     ->Unit(benchmark::kMicrosecond)
     ->Iterations(200);
 
+// ONE FRAME OF A CAGE DRAG, as a host previews one: the cage laid over the mesh
+// with a record, then taken back, so every iteration starts from the same form.
+// One control point dragged, at 3^3, 8^3 and the 32^3 ceiling, over ~100k
+// triangles.
+//
+// The three rows should sit close together. The cost of an evaluation is the
+// basis — O(n) an axis — plus one term per DRAGGED point, so a larger cage with
+// the same point in hand costs a little more and not n^3 more. Before the sum
+// ran over the dragged points alone it ran over all of them, and the 32^3 row
+// was ~1.7 s a frame on a 62k-vertex mesh against ~10 ms at 3^3.
+void BM_MeshLatticeDrag(benchmark::State& state) {
+    using namespace clay;
+    using namespace clay::kernel;
+    const int n = static_cast<int>(state.range(0));
+    mesh::Mesh patch = bench_surface_patch(224, 0.02f);
+    math::Aabb box;
+    for (const cfloat3& p : patch.positions) box.expand(p);
+    mesh::Lattice cage(box, n, n, n);
+    cage.set_offset(0, 0, 0, cf3(0.05f, 0.02f, 0.0f));
+    mesh::MeshSculptor sculptor(patch, 0.0f);
+
+    std::size_t moved = 0;
+    for (auto _ : state) {
+        mesh::VertexDeltas record;
+        moved = sculptor.apply_lattice(cage, &record);
+        record.revert(patch);
+    }
+    // A cage that moved nothing would time the identity early-out.
+    state.counters["moved"] = static_cast<double>(moved);
+    state.counters["dragged"] = static_cast<double>(cage.dragged_count());
+}
+BENCHMARK(BM_MeshLatticeDrag)->Arg(3)->Arg(8)->Arg(32)->Unit(benchmark::kMillisecond);
+
 // The same pair on the ADAPTIVE surface, where the topology object is a
 // different implementation of the same two methods — `one_ring` through the
 // workset's slot map, and `DynamicSurface::is_boundary_edge`. If the virtual

docs/07-brushes-and-features.md

@@ -3084,7 +3084,7 @@ difference:
 | **Alphas** | `Deformer::alpha` | `sculpt_carve_alpha` | — | On an SDF item an alpha is a **deformer**, not a primitive: an item shaped like the stamp would add material in the stamp's shape, where an alpha modulates a surface already there. The engine decodes no images on either side — a host hands over samples |
 | **Mask extrude** | `brush::mask_extrude` | `VoxelGrid::mask_extrude` | — | The mask is the region, so there is no radius. The SDF result is an ordinary item and stays an operand |
 | **Taper / twist** | `Deformer::taper`, `twist`, and the ranged forms | — | `MeshSculptor::deform` | The SDF deformers run **backwards** — they answer "where did the material at this point come from". The mesh versions are **forward** point maps, which is both easier and exact, so a tapered mesh and a tapered field are the same shape. There is deliberately no mesh `bend`: its map folds distinct points together past a gentle angle, so no forward map exists |
-| **Lattice / FFD** | `Deformer::lattice` on an item, `lattice_transformed` through a frame; `brush::lattice_gizmo` over a whole layer | — | `mesh::Lattice`, `MeshSculptor::apply_lattice` | **Direction.** The mesh cage runs FORWARD — `p + Bernstein(offsets, clamp(s,t,u))` — exact, once per vertex, up to 32 divisions an axis. An SDF deformer must run backwards and forward FFD has no closed-form inverse, so the SDF cage authors its offsets AS the inverse warp: not the exact inverse, differing by a term that over-travels a drag toward rising basis weight and under-travels one pointing away, measured under 1.5% of the drag. It costs `nx*ny*nz` multiply-adds **per sample inside the raymarcher**, hence 4 divisions an axis rather than 32. Both store offsets, so an untouched cage is exactly the identity, and both hold rigidly beyond the box |
+| **Lattice / FFD** | `Deformer::lattice` on an item, `lattice_transformed` through a frame; `brush::lattice_gizmo` over a whole layer | — | `mesh::Lattice`, `MeshSculptor::apply_lattice` | **Direction.** The mesh cage runs FORWARD — `p + Bernstein(offsets, clamp(s,t,u))` — exact, once per vertex, up to 32 divisions an axis, and priced by the control points actually DRAGGED rather than by the cage: a point at rest adds nothing, so the sum skips it and one corner of a 32^3 cage costs one term, not 32,768 (see `BM_MeshLatticeDrag`). An SDF deformer must run backwards and forward FFD has no closed-form inverse, so the SDF cage authors its offsets AS the inverse warp: not the exact inverse, differing by a term that over-travels a drag toward rising basis weight and under-travels one pointing away, measured under 1.5% of the drag. It costs `nx*ny*nz` multiply-adds **per sample inside the raymarcher**, hence 4 divisions an axis rather than 32. Both store offsets, so an untouched cage is exactly the identity, and both hold rigidly beyond the box |
 | **Masking** | gates any **operation**, a boolean included | freezes cells against every verb | gates every brush | The SDF gate rides the combine record rather than being a mode, which is why it protects from a boolean the same way it protects from a brush |
 | **Strokes** | `stamps_to_nodes` | `apply_to_grid` | `apply_to_mesh` | One resolved stroke, four consumers (the fourth is `apply_to_mask`). Spacing, pressure, jitter and taper mean the same thing on all of them |
 

include/clay/mesh/lattice.h

@@ -73,11 +73,24 @@ class Lattice {
     kernel::cfloat3 position(int i, int j, int k) const { return rest(i, j, k) + offset(i, j, k); }
 
     // Every offset zero: the identity, and worth asking before walking a mesh.
-    bool is_identity() const;
+    bool is_identity() const { return dragged_.empty(); }
+
+    // How many control points carry a non-zero offset — what one evaluation
+    // costs, since the points nobody dragged contribute nothing and are not
+    // visited.
+    std::size_t dragged_count() const { return dragged_.size(); }
 
     // The displacement this cage applies at a world point. Zero everywhere for
     // an untouched cage, exactly.
     //
+    // PRICED BY THE DRAGGED POINTS, NOT BY THE CAGE. The offset field is linear
+    // in the offsets, so a point left at rest adds exactly nothing; the sum
+    // runs over `dragged_` alone, and the three axis bases are built in O(n)
+    // each. Summing every control point instead made one evaluation cost
+    // nx * ny * nz multiply-adds whatever was dragged — 32,768 at the 32^3
+    // ceiling, which is 1.7 s for one pass over a 62k-vertex mesh and so
+    // seconds per frame of a host's interactive preview.
+    //
     // An axis on which the box is FLAT — which is what a cage over a plane's
     // own bounds gives — reads as the middle rather than as an end, so none of
     // the control points on it are dead. See the note in the .cpp.
@@ -89,6 +102,11 @@ class Lattice {
     math::Aabb box_;
     int nx_ = 3, ny_ = 3, nz_ = 3;
     std::vector<kernel::cfloat3> offsets_;
+    // The flat indices whose offset is non-zero, in the order they were first
+    // dragged. Kept by `set_offset`, so the set and `offsets_` cannot disagree.
+    std::vector<std::uint32_t> dragged_;
+    // Each axis's binomial coefficients C(n - 1, i), fixed by the divisions.
+    std::vector<double> binomial_x_, binomial_y_, binomial_z_;
 };
 
 }  // namespace mesh

openspec/changes/price-a-mesh-cage-by-its-dragged-points/proposal.md

@@ -0,0 +1,35 @@
+# Price a mesh cage by its dragged points
+
+## Why
+
+`mesh::Lattice::displacement` summed every control point of the cage at every
+vertex: `nx * ny * nz` multiply-adds whatever had been dragged. At the 32^3
+ceiling that is 32,768 terms per vertex, and a host that previews a cage drag
+by laying the cage over the mesh on every pointer move paid it on every frame.
+ClaySpace measured one frame at ~1.7 s on a 62k-vertex mesh (4.7-5.2 s on an
+80k-triangle fixture), against ~10 ms at 3^3, with a single corner in hand
+(CyberdyneCorp/ClaySpaceDesktop#176).
+
+The offset field is linear in the offsets, so a control point left at rest adds
+exactly nothing to the sum. Its cost is work with no result.
+
+## What changes
+
+- The cage keeps the set of control points carrying a non-zero offset,
+  maintained by `set_offset`, and `displacement` sums over that set alone.
+  `is_identity` becomes that set being empty, in constant time.
+- Each axis's Bernstein basis is built in O(n) from powers of `t` and `1 - t`
+  with precomputed binomials, in double, instead of de Casteljau's O(n^2)
+  recurrence in float — cheap beside the old n^3 sum, and not beside the new
+  one.
+- `dragged_count()` reports the size of that set, which is what one
+  evaluation costs.
+- `BM_MeshLatticeDrag` times one preview frame (apply with a record, then
+  revert) at 3^3, 8^3 and 32^3 with one point dragged over ~100k triangles.
+
+## Scope
+
+`mesh::Lattice` only. No C ABI, Python or file-format change: every entry point
+keeps its signature and meaning, and the results agree with the full sum to
+float resolution. The SDF lattice deformer (`clattice_point`) is a different
+evaluator, capped at 4 divisions an axis, and is untouched.

openspec/changes/price-a-mesh-cage-by-its-dragged-points/specs/meshing/spec.md

@@ -0,0 +1,30 @@
+## ADDED Requirements
+
+### Requirement: A mesh cage is priced by the control points that were dragged
+Evaluating a mesh lattice cage SHALL cost in proportion to the control points
+that carry a non-zero offset, not to the number of control points the cage
+holds. A control point at rest contributes exactly nothing to the offset field,
+so the evaluation SHALL NOT visit it.
+
+The per-axis basis SHALL cost O(n) in the divisions on that axis, so that a
+cage of the largest size with one point dragged costs a small multiple of the
+smallest cage with the same point dragged, rather than n^3 times as much.
+
+The result SHALL agree with the full trivariate Bernstein sum over every
+control point to float resolution, inside the box and past it where the
+parameters clamp; Bernstein SHALL still interpolate the corner control points
+exactly; and a cage whose every point has been returned to rest SHALL again be
+exactly the identity.
+
+#### Scenario: One corner of the largest cage
+- **WHEN** one control point of a 32x32x32 cage is dragged and the cage is applied to a mesh
+- **THEN** the evaluation sums one term per vertex, not 32,768
+- **AND** the time per application is within a small multiple of the same drag on a 3x3x3 cage
+
+#### Scenario: The sparse sum is the whole sum
+- **WHEN** a cage of 2, 3, 8 or 32 divisions an axis has a few points dragged
+- **THEN** the displacement at points inside, on and outside its box matches the full Bernstein sum over every control point to float resolution
+
+#### Scenario: A point put back is forgotten
+- **WHEN** a dragged control point is set back to a zero offset
+- **THEN** it no longer counts as dragged, and a cage with none left is the identity

openspec/changes/price-a-mesh-cage-by-its-dragged-points/tasks.md

@@ -0,0 +1,6 @@
+- [x] Keep the dragged control points in `mesh::Lattice`, maintained by `set_offset`.
+- [x] Sum `displacement` over the dragged points only, with an O(n) basis per axis.
+- [x] Hold the new evaluation against the full trivariate Bernstein sum at 2, 3, 8 and 32 divisions, inside and outside the box.
+- [x] Pin the corner interpolation at degree 31 and the dragged count through a point returned to rest.
+- [x] Add `BM_MeshLatticeDrag` at 3^3, 8^3 and 32^3.
+- [ ] Ship in a release, and move ClaySpace's engine pin to it (CyberdyneCorp/ClaySpaceDesktop#176).

src/mesh/lattice.cpp

@@ -14,21 +14,36 @@ int clamp_div(int n) {
     return std::max(kMinLatticeDivisions, std::min(kMaxLatticeDivisions, n));
 }
 
-// Bernstein basis of degree n at t, written into `out` (n + 1 values).
+// Each C(count - 1, i), the coefficients of the degree-(count - 1) Bernstein
+// basis. Exact in a double for every count up to kMaxLatticeDivisions: the
+// largest, C(31, 15), is about 3e8.
+std::vector<double> binomials(int count) {
+    std::vector<double> c(static_cast<std::size_t>(count), 1.0);
+    for (int i = 1; i < count; ++i) c[i] = c[i - 1] * static_cast<double>(count - i) / i;
+    return c;
+}
+
+// Bernstein basis of degree count - 1 at t, written into `out` (count values).
 //
-// De Casteljau's recurrence rather than binomial coefficients times powers:
-// it needs no factorials, cannot overflow at the degrees a cage reaches, and
-// is numerically the better-behaved of the two. Each row is built from the one
-// above by b[i] = (1 - t) * b[i] + t * b[i - 1].
-void bernstein(int degree, float t, float* out) {
-    out[0] = 1.0f;
-    for (int d = 1; d <= degree; ++d) {
-        out[d] = t * out[d - 1];
-        for (int i = d - 1; i > 0; --i) out[i] = (1.0f - t) * out[i] + t * out[i - 1];
-        out[0] = (1.0f - t) * out[0];
+// Built from the powers of t and 1 - t, so an axis costs O(n) rather than the
+// O(n^2) of de Casteljau's recurrence — which was cheap beside the n^3 sum it
+// fed, and is not beside the sum over dragged points that replaced it. In
+// double, so the partition of unity holds to far below a float's resolution
+// and a uniformly dragged cage still translates exactly once rounded back.
+void bernstein(const std::vector<double>& binomial, double t, double* out) {
+    const int count = static_cast<int>(binomial.size());
+    double rest[kMaxLatticeDivisions];
+    rest[0] = 1.0;
+    for (int m = 1; m < count; ++m) rest[m] = rest[m - 1] * (1.0 - t);
+    double power = 1.0;
+    for (int i = 0; i < count; ++i) {
+        out[i] = binomial[i] * power * rest[count - 1 - i];
+        power *= t;
     }
 }
 
+bool is_zero(cfloat3 v) { return v.x == 0.0f && v.y == 0.0f && v.z == 0.0f; }
+
 // Where p sits in the box, on one axis, clamped to [0, 1].
 //
 // Clamped because the offset field is only DEFINED over the cage; a vertex
@@ -54,6 +69,9 @@ float axis_parameter(float v, float lo, float hi) {
 Lattice::Lattice(const math::Aabb& box, int nx, int ny, int nz)
     : box_(box), nx_(clamp_div(nx)), ny_(clamp_div(ny)), nz_(clamp_div(nz)) {
     offsets_.assign(static_cast<std::size_t>(nx_) * ny_ * nz_, cf3(0, 0, 0));
+    binomial_x_ = binomials(nx_);
+    binomial_y_ = binomials(ny_);
+    binomial_z_ = binomials(nz_);
 }
 
 std::size_t Lattice::index(int i, int j, int k) const {
@@ -67,7 +85,15 @@ cfloat3 Lattice::offset(int i, int j, int k) const {
 
 void Lattice::set_offset(int i, int j, int k, cfloat3 v) {
     if (i < 0 || j < 0 || k < 0 || i >= nx_ || j >= ny_ || k >= nz_) return;
-    offsets_[index(i, j, k)] = v;
+    const std::size_t at = index(i, j, k);
+    const bool was = !is_zero(offsets_[at]);
+    const bool now = !is_zero(v);
+    offsets_[at] = v;
+    if (now && !was) {
+        dragged_.push_back(static_cast<std::uint32_t>(at));
+    } else if (was && !now) {
+        dragged_.erase(std::find(dragged_.begin(), dragged_.end(), static_cast<std::uint32_t>(at)));
+    }
 }
 
 cfloat3 Lattice::rest(int i, int j, int k) const {
@@ -81,39 +107,32 @@ cfloat3 Lattice::rest(int i, int j, int k) const {
                along(box_.min.z, box_.max.z, k, nz_));
 }
 
-bool Lattice::is_identity() const {
-    for (const cfloat3& o : offsets_)
-        if (o.x != 0.0f || o.y != 0.0f || o.z != 0.0f) return false;
-    return true;
-}
-
 cfloat3 Lattice::displacement(cfloat3 p) const {
-    if (box_.empty()) return cf3(0, 0, 0);
-    const float s = axis_parameter(p.x, box_.min.x, box_.max.x);
-    const float t = axis_parameter(p.y, box_.min.y, box_.max.y);
-    const float u = axis_parameter(p.z, box_.min.z, box_.max.z);
+    if (box_.empty() || dragged_.empty()) return cf3(0, 0, 0);
+    const double s = axis_parameter(p.x, box_.min.x, box_.max.x);
+    const double t = axis_parameter(p.y, box_.min.y, box_.max.y);
+    const double u = axis_parameter(p.z, box_.min.z, box_.max.z);
 
     // Degree is one less than the point count, so a two-point axis is degree
     // one — exactly linear — and no separate trilinear path is needed.
-    float bx[kMaxLatticeDivisions], by[kMaxLatticeDivisions], bz[kMaxLatticeDivisions];
-    bernstein(nx_ - 1, s, bx);
-    bernstein(ny_ - 1, t, by);
-    bernstein(nz_ - 1, u, bz);
-
-    cfloat3 sum = cf3(0, 0, 0);
-    for (int k = 0; k < nz_; ++k) {
-        if (bz[k] == 0.0f) continue;
-        for (int j = 0; j < ny_; ++j) {
-            const float wjk = by[j] * bz[k];
-            if (wjk == 0.0f) continue;
-            for (int i = 0; i < nx_; ++i) {
-                const float w = bx[i] * wjk;
-                if (w == 0.0f) continue;
-                sum = sum + offsets_[index(i, j, k)] * w;
-            }
-        }
+    double bx[kMaxLatticeDivisions], by[kMaxLatticeDivisions], bz[kMaxLatticeDivisions];
+    bernstein(binomial_x_, s, bx);
+    bernstein(binomial_y_, t, by);
+    bernstein(binomial_z_, u, bz);
+
+    // Only the points somebody dragged: every other offset is zero and would
+    // add exactly nothing. See the header for what the full sum cost.
+    const std::size_t row = static_cast<std::size_t>(nx_);
+    const std::size_t plane = row * static_cast<std::size_t>(ny_);
+    double sx = 0.0, sy = 0.0, sz = 0.0;
+    for (const std::uint32_t at : dragged_) {
+        const double w = bx[at % row] * by[(at / row) % ny_] * bz[at / plane];
+        const cfloat3& o = offsets_[at];
+        sx += w * o.x;
+        sy += w * o.y;
+        sz += w * o.z;
     }
-    return sum;
+    return cf3(static_cast<float>(sx), static_cast<float>(sy), static_cast<float>(sz));
 }
 
 }  // namespace mesh