From 757d62dee2efd96adf64698364d67554c7a337d9 Mon Sep 17 00:00:00 2001 From: Leonardo Araujo dos Santos Date: Sat, 26 Sep 2026 12:38:53 -0300 Subject: [PATCH] Price a mesh cage by the control points that were dragged Lattice::displacement summed every control point at every vertex, so one evaluation cost nx*ny*nz multiply-adds whatever had been dragged: 32,768 terms per vertex at the 32^3 ceiling, about 1.7 s per application over a 62k-vertex mesh with a single corner in hand. A host previewing a cage drag pays that on every pointer move. The offset field is linear in the offsets, so a point at rest adds nothing. The cage now keeps the set of dragged points (maintained by set_offset) and sums over it alone; is_identity is that set being empty. Each axis basis is built in O(n) from powers of t and 1-t with precomputed binomials, in double. Results agree with the full Bernstein sum to float resolution; no ABI, Python or format change. BM_MeshLatticeDrag times one preview frame at 3^3, 8^3 and 32^3. --- benchmarks/bench_main.cpp | 34 ++++++ docs/07-brushes-and-features.md | 2 +- include/clay/mesh/lattice.h | 20 +++- .../proposal.md | 35 ++++++ .../specs/meshing/spec.md | 30 +++++ .../tasks.md | 6 + src/mesh/lattice.cpp | 99 ++++++++++------- tests/unit/test_mesh_lattice.cpp | 105 ++++++++++++++++++ 8 files changed, 289 insertions(+), 42 deletions(-) create mode 100644 openspec/changes/price-a-mesh-cage-by-its-dragged-points/proposal.md create mode 100644 openspec/changes/price-a-mesh-cage-by-its-dragged-points/specs/meshing/spec.md create mode 100644 openspec/changes/price-a-mesh-cage-by-its-dragged-points/tasks.md diff --git a/benchmarks/bench_main.cpp b/benchmarks/bench_main.cpp index ce57d359..13ffe5c6 100644 --- a/benchmarks/bench_main.cpp +++ b/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(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(moved); + state.counters["dragged"] = static_cast(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 diff --git a/docs/07-brushes-and-features.md b/docs/07-brushes-and-features.md index 4956d392..9a78f1b1 100644 --- a/docs/07-brushes-and-features.md +++ b/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 | diff --git a/include/clay/mesh/lattice.h b/include/clay/mesh/lattice.h index f0b38518..2a779b46 100644 --- a/include/clay/mesh/lattice.h +++ b/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 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 dragged_; + // Each axis's binomial coefficients C(n - 1, i), fixed by the divisions. + std::vector binomial_x_, binomial_y_, binomial_z_; }; } // namespace mesh diff --git a/openspec/changes/price-a-mesh-cage-by-its-dragged-points/proposal.md b/openspec/changes/price-a-mesh-cage-by-its-dragged-points/proposal.md new file mode 100644 index 00000000..c9fd8952 --- /dev/null +++ b/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. diff --git a/openspec/changes/price-a-mesh-cage-by-its-dragged-points/specs/meshing/spec.md b/openspec/changes/price-a-mesh-cage-by-its-dragged-points/specs/meshing/spec.md new file mode 100644 index 00000000..7c1c298f --- /dev/null +++ b/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 diff --git a/openspec/changes/price-a-mesh-cage-by-its-dragged-points/tasks.md b/openspec/changes/price-a-mesh-cage-by-its-dragged-points/tasks.md new file mode 100644 index 00000000..a852df75 --- /dev/null +++ b/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). diff --git a/src/mesh/lattice.cpp b/src/mesh/lattice.cpp index 044f2fe3..cdaaa0cd 100644 --- a/src/mesh/lattice.cpp +++ b/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 binomials(int count) { + std::vector c(static_cast(count), 1.0); + for (int i = 1; i < count; ++i) c[i] = c[i - 1] * static_cast(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& binomial, double t, double* out) { + const int count = static_cast(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(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(at)); + } else if (was && !now) { + dragged_.erase(std::find(dragged_.begin(), dragged_.end(), static_cast(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(nx_); + const std::size_t plane = row * static_cast(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(sx), static_cast(sy), static_cast(sz)); } } // namespace mesh diff --git a/tests/unit/test_mesh_lattice.cpp b/tests/unit/test_mesh_lattice.cpp index 42342863..32f23c6e 100644 --- a/tests/unit/test_mesh_lattice.cpp +++ b/tests/unit/test_mesh_lattice.cpp @@ -302,3 +302,108 @@ TEST_CASE("a flat cage still works on the axes that are not flat") { for (const cfloat3& p : m.positions) highest = std::max(highest, p.y); CHECK(highest > 0.1f); } + +namespace { + +// The offset field summed over EVERY control point, with the basis built by de +// Casteljau's recurrence in double — the formula as the spec states it, written +// the slow and obvious way, to hold the dragged-point evaluation against. +cfloat3 reference_displacement(const Lattice& cage, cfloat3 p) { + const math::Aabb& box = cage.box(); + const auto parameter = [](float v, float lo, float hi) { + const float span = hi - lo; + if (!(span > 1e-9f)) return 0.5; + return static_cast(std::min(1.0f, std::max(0.0f, (v - lo) / span))); + }; + const auto basis = [](int count, double t) { + std::vector b(static_cast(count), 0.0); + b[0] = 1.0; + for (int d = 1; d < count; ++d) { + b[d] = t * b[d - 1]; + for (int i = d - 1; i > 0; --i) b[i] = (1.0 - t) * b[i] + t * b[i - 1]; + b[0] = (1.0 - t) * b[0]; + } + return b; + }; + const std::vector bx = basis(cage.nx(), parameter(p.x, box.min.x, box.max.x)); + const std::vector by = basis(cage.ny(), parameter(p.y, box.min.y, box.max.y)); + const std::vector bz = basis(cage.nz(), parameter(p.z, box.min.z, box.max.z)); + double sx = 0.0, sy = 0.0, sz = 0.0; + for (int k = 0; k < cage.nz(); ++k) + for (int j = 0; j < cage.ny(); ++j) + for (int i = 0; i < cage.nx(); ++i) { + const double w = bx[i] * by[j] * bz[k]; + const cfloat3 o = cage.offset(i, j, k); + sx += w * o.x; + sy += w * o.y; + sz += w * o.z; + } + return cf3(static_cast(sx), static_cast(sy), static_cast(sz)); +} + +} // namespace + +TEST_CASE("a cage counts the points that were dragged, and forgets one put back") { + // What one evaluation costs is the number of dragged points, so the count + // has to follow every write — including one that returns a point to rest. + Lattice cage(math::Aabb{cf3(-1, -1, -1), cf3(1, 1, 1)}, 32, 32, 32); + CHECK(cage.dragged_count() == 0); + cage.set_offset(0, 0, 0, cf3(0.05f, 0, 0)); + cage.set_offset(31, 31, 31, cf3(0, 0.1f, 0)); + CHECK(cage.dragged_count() == 2); + // Dragged again is still one point. + cage.set_offset(0, 0, 0, cf3(0.07f, 0, 0)); + CHECK(cage.dragged_count() == 2); + // Back to rest leaves the sum, and the cage is the identity again. + cage.set_offset(0, 0, 0, cf3(0, 0, 0)); + cage.set_offset(31, 31, 31, cf3(0, 0, 0)); + CHECK(cage.dragged_count() == 0); + CHECK(cage.is_identity()); + const cfloat3 d = cage.displacement(cf3(-1, -1, -1)); + CHECK(d.x == 0.0f); + CHECK(d.y == 0.0f); + CHECK(d.z == 0.0f); +} + +TEST_CASE("summing only the dragged points is the whole Bernstein sum") { + // The optimisation is only that a zero offset adds nothing, so it must + // agree with the full trivariate sum everywhere — inside the box, on its + // faces, and past it where the parameters clamp — up to the 32^3 ceiling + // and on an uneven cage. + const math::Aabb box{cf3(-0.6f, -0.5f, -0.4f), cf3(0.7f, 0.5f, 0.45f)}; + for (const int n : {2, 3, 8, 32}) { + CAPTURE(n); + Lattice cage(box, n, n == 2 ? 5 : n, n); + cage.set_offset(0, 0, 0, cf3(0.05f, -0.02f, 0.01f)); + cage.set_offset(n - 1, cage.ny() / 2, n / 3, cf3(-0.1f, 0.2f, 0.03f)); + cage.set_offset(n / 2, cage.ny() - 1, n - 1, cf3(0.0f, 0.0f, -0.3f)); + for (int step = 0; step <= 20; ++step) { + const float a = -0.9f + 1.8f * static_cast(step) / 20.0f; + const cfloat3 p = cf3(a, 0.8f * a - 0.1f, 0.3f - 0.5f * a); + CAPTURE(step); + const cfloat3 got = cage.displacement(p); + const cfloat3 want = reference_displacement(cage, p); + CHECK(got.x == doctest::Approx(want.x).epsilon(1e-6).scale(1e-6)); + CHECK(got.y == doctest::Approx(want.y).epsilon(1e-6).scale(1e-6)); + CHECK(got.z == doctest::Approx(want.z).epsilon(1e-6).scale(1e-6)); + } + } +} + +TEST_CASE("a corner drag on the largest cage still moves its corner exactly") { + // At degree 31 the corner basis is t^31 — the term most exposed to how the + // basis is built — and Bernstein interpolates the corners, so the corner of + // the box moves by exactly the corner's offset and the far one not at all. + const math::Aabb box{cf3(-1, -1, -1), cf3(1, 1, 1)}; + Lattice cage(box, 32, 32, 32); + const cfloat3 pull = cf3(0.05f, -0.04f, 0.03f); + cage.set_offset(31, 31, 31, pull); + const cfloat3 at_corner = cage.displacement(cf3(1, 1, 1)); + CHECK(at_corner.x == pull.x); + CHECK(at_corner.y == pull.y); + CHECK(at_corner.z == pull.z); + const cfloat3 far = cage.displacement(cf3(-1, -1, -1)); + CHECK(far.x == 0.0f); + CHECK(far.y == 0.0f); + CHECK(far.z == 0.0f); +}