Shape Similarity and Self-Assembly
println("Julia ", VERSION)Julia 1.12.6Overview
Every solid in Rasmah is an SDF (signed-distance field) — a callable f(x) that is negative inside, positive outside, and zero on the boundary. Because an SDF is a differentiable function of the shape’s parameters (and of any rigid transform applied to it), we can turn “what pose puts these pieces back together?” into an ordinary optimization loop:
pick a pose θ → measure the mismatch L(θ) → gradient ∇L → θ ← θ − step·∇L → repeatThat question is not abstract: it is the one an engineer asks when reassembling a fractured part — a block broken into fragments that must slide back into their original positions. This tutorial spells out that loop by hand, end to end: the loss, the gradient, the LBFGS call, and the animation of the reassembly. Along the way it builds the two ideas the loop depends on — a similarity measure to act as the loss, and a soft union so every piece keeps gradient while still shuffled.
using Rasmah
using LinearAlgebraDownloading artifact: vtkjs Failure artifact: vtkjs Downloading artifact: vtkjs
1. Similarity measures (the loss)
The mismatch L is built from a distance between two shapes. surface_distance samples one isosurface and evaluates the other SDF there, reporting the mean / RMS / max gap plus the symmetric Chamfer and Hausdorff distances:
a = to_sdf(sphere(1.0))
b = to_sdf(translate(sphere(1.0), 0.5, 0.0, 0.0))
sd = surface_distance(a, b; n=32)
println("mean = ", round(sd.mean; digits=4))
println("rms = ", round(sd.rms; digits=4))
println("max = ", round(sd.max; digits=4))
println("chamfer = ", round(sd.chamfer; digits=4))
println("hausdorff = ", round(sd.hausdorff; digits=4))mean = 0.2498
rms = 0.2887
max = 0.5
chamfer = 0.2498
hausdorff = 0.5volume_overlap integrates a smooth Heaviside of each field over a grid and returns the intersection/union plus Dice, Jaccard, and the overlap ratio (the standard volumetric similarity scores):
lo, hi = [-2.0, -2.0, -2.0], [2.0, 2.0, 2.0]
ov = volume_overlap(a, b, lo, hi; n=32)
println("dice = ", round(ov.dice; digits=4))
println("jaccard = ", round(ov.jaccard; digits=4))
println("overlap = ", round(ov.overlap_ratio; digits=4))dice = 0.5646
jaccard = 0.3933
overlap = 0.5646l2_field_distance is the literal “sum of squared SDF differences” you might reach for first, done carefully: the integrand (f−g)² is weighted by a Gaussian band exp(−(f²+g²)/2σ²) so the integral stays near the two surfaces (where the fields are meaningful) instead of being dominated by the far field:
d = l2_field_distance(a, b, lo, hi; n=32)
println("l2 = ", round(d.l2; digits=4))
println("rms = ", round(d.rms; digits=4))l2 = 0.0803
rms = 0.2833The two shapes these measures compare — a unit sphere and a copy shifted by 0.5 along x:
vtk(surface_mesh(union(a, b); n=40))2. The assembly problem
We slice a 2×2×2 box into four slabs as SDFs (a slab is the box intersected with a thin translated box) — a block sawn into four pieces. The “shuffle” has two steps: first translate each slab so the four pieces are physically separated (they no longer touch or overlap), then give each a minor rotation about its own centroid so it sits tilted where it landed:
B = to_sdf(box(2.0, 2.0, 2.0))
slab(xc, hx) = intersection(B, translate(to_sdf(box(2hx, 2.0, 2.0)), xc, 0.0, 0.0))
base = [slab(-0.75, 0.25), slab(-0.25, 0.25), slab(0.25, 0.25), slab(0.75, 0.25)]
# step 1: translate apart (non-overlapping)
separations = [
[-0.55, 0.10, 0.05], [-0.20, 0.25, -0.10], [0.20, -0.20, 0.12], [0.55, -0.05, -0.05]
]
# step 2: a minor rotation about each piece's own centroid
rotations = [rotation_z(0.15), rotation_x(0.12), rotation_x(0.10), rotation_y(0.12)]
pieces = [
begin
moved = translate(p, t...) # separate
c = (bounds(moved)[1] .+ bounds(moved)[2]) ./ 2 # its centroid
translate(rotate(translate(moved, -c...), R), c...) # tilt in place
end for (p, t, R) in zip(base, separations, rotations)
]4-element Vector{TranslateSDF{RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}, Vector{Float64}}}:
TranslateSDF{RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}, Vector{Float64}}(RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}(TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}(TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}(IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}(BoxSDF{Float64, Float64, Float64}(1.0, 1.0, 1.0), TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}(BoxSDF{Float64, Float64, Float64}(0.25, 1.0, 1.0), [-0.75, 0.0, 0.0])), [-0.55, 0.1, 0.05]), [1.3, -0.10000000000000003, -0.050000000000000044]), [0.9887710779360422 -0.14943813247359922 0.0; 0.14943813247359922 0.9887710779360422 0.0; 0.0 0.0 1.0]), [-1.3, 0.10000000000000003, 0.050000000000000044])
TranslateSDF{RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}, Vector{Float64}}(RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}(TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}(TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}(IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}(BoxSDF{Float64, Float64, Float64}(1.0, 1.0, 1.0), TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}(BoxSDF{Float64, Float64, Float64}(0.25, 1.0, 1.0), [-0.25, 0.0, 0.0])), [-0.2, 0.25, -0.1]), [0.44999999999999996, -0.25, 0.10000000000000003]), [1.0 0.0 0.0; 0.0 0.9928086358538663 -0.11971220728891936; 0.0 0.11971220728891936 0.9928086358538663]), [-0.44999999999999996, 0.25, -0.10000000000000003])
TranslateSDF{RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}, Vector{Float64}}(RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}(TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}(TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}(IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}(BoxSDF{Float64, Float64, Float64}(1.0, 1.0, 1.0), TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}(BoxSDF{Float64, Float64, Float64}(0.25, 1.0, 1.0), [0.25, 0.0, 0.0])), [0.2, -0.2, 0.12]), [-0.44999999999999996, 0.19999999999999996, -0.12000000000000005]), [1.0 0.0 0.0; 0.0 0.9950041652780258 -0.09983341664682815; 0.0 0.09983341664682815 0.9950041652780258]), [0.44999999999999996, -0.19999999999999996, 0.12000000000000005])
TranslateSDF{RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}, Vector{Float64}}(RotateSDF{TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}, Matrix{Float64}}(TranslateSDF{TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}, Vector{Float64}}(TranslateSDF{IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}, Vector{Float64}}(IntersectionSDF{BoxSDF{Float64, Float64, Float64}, TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}}(BoxSDF{Float64, Float64, Float64}(1.0, 1.0, 1.0), TranslateSDF{BoxSDF{Float64, Float64, Float64}, Vector{Float64}}(BoxSDF{Float64, Float64, Float64}(0.25, 1.0, 1.0), [0.75, 0.0, 0.0])), [0.55, -0.05, -0.05]), [-1.3, 0.050000000000000044, 0.050000000000000044]), [0.9928086358538663 0.0 0.11971220728891936; 0.0 1.0 0.0; -0.11971220728891936 0.0 0.9928086358538663]), [1.3, -0.050000000000000044, -0.050000000000000044])The four separated, tilted pieces (the starting configuration to un-shuffle):
initial = reduce(union, pieces)
vtk(surface_mesh(initial; n=40))3. The pose parameterization
θ is always a flat vector [t₁x, t₁y, t₁z, ω₁x, ω₁y, ω₁z, t₂x, …] of length 6n. Translation is a plain 3-vector; rotation is axis-angle — the vector ω is turned into a rotation matrix by rotation_exp(skew(ω)) (Rodrigues’ formula), which is smooth at ω = 0 (no gimbal lock):
println("R(ω) at ω=[0,0,0.6]:\n", round.(rotation_exp(skew([0.0, 0.0, 0.6])); digits=3))
println("skew([1,2,3]) = \n", skew([1.0, 2.0, 3.0]))R(ω) at ω=[0,0,0.6]:
[0.825 -0.565 0.0; 0.565 0.825 0.0; 0.0 0.0 1.0]
skew([1,2,3]) =
[0.0 -3.0 2.0; 3.0 0.0 -1.0; -2.0 1.0 0.0]4. The optimization loop
The unknown θ stacks one pose [t; ω] per piece — 24 numbers for four pieces (§3). The loss L(θ) moves every piece by its pose, unions the results, and scores the match two ways: the RMS gap on the reference surface, plus a volumetric Dice term that catches internal gaps and overlaps the surface alone cannot see. The two pieces of that design — the smooth soft union and the Dice term — are the subject of §5 and §6:
# a fixed domain covering the box plus some overflow
loD, hiD = bounds(B)
pad = 0.5 .* (hiD .- loD)
loD, hiD = loD .- pad, hiD .+ pad
# reference-surface sample points (fixed; only θ is differentiated)
xs = [collect(c) for c in eachcol(surface_mesh(B; n=40).vertices)]
function loss(θ)
moved = [
rotate(
translate(pieces[i], θ[6i - 5], θ[6i - 4], θ[6i - 3]),
rotation_exp(skew(θ[(6i - 2):6i])),
) for i in eachindex(pieces)
]
A = SoftUnion(moved, 0.05)
s = zero(eltype(θ))
for x in xs
v = A(x)
s += v * v
end
surf = sqrt(s / length(xs))
ov = volume_overlap(B, A, loD, hiD; n=24, k=20.0)
surf + (1 - ov.dice)
end
grad = θ -> gradient(loss, θ)#9 (generic function with 1 method)Two things make the gradient well-defined: the sample points xs are fixed (so only the union is differentiated), and rotation_exp is smooth at ω = 0 (its small-angle branch returns I + skew(ω), not a bare identity).
Then one minimize call (LBFGS) runs the loop to convergence. store_trace=true records the sequence of poses θ the optimizer visits (the start included), so that history can be replayed as the reassembly animation in §7:
res_hh = @time minimize(
ObjectiveFunction(loss, grad),
zeros(24);
iterations=30,
store_trace=true,
)
println("minimize loss = ", round(loss(res_hh.minimizer); digits=6)) 8.807357 seconds (39.39 M allocations: 3.984 GiB, 11.24% gc time, 45.51% compilation time)
minimize loss = 0.1759The hand-written loop’s solution, assembled back into a solid:
θhh = res_hh.minimizer
assembled = reduce(
union,
[
rotate(
translate(pieces[i], θhh[6i - 5], θhh[6i - 4], θhh[6i - 3]),
rotation_exp(skew(θhh[(6i - 2):6i])),
) for i in eachindex(pieces)
],
)
println("reassembled volume = ", round(volume(surface_mesh(assembled; n=32)); digits=3))
println("reference volume = ", round(volume(surface_mesh(B; n=32)); digits=3))
vtk(surface_mesh(assembled; n=40))reassembled volume = 7.725 m^3
reference volume = 7.983 m^35. Why the soft union?
The loss unions the moved pieces and samples it on the reference. With a hard min union, each sample point evaluates only the closest piece, so a piece shuffled far from its home can receive no gradient and never move. The soft union replaces min with the smooth log-sum-exp
\[-\!k \log\sum_i \exp\!\big(-f_i(x)/k\big),\]
which pulls every piece toward every reference point it is near. As k \to 0 the soft union converges to the true min, so the converged result is the same box.
6. Why the loss must be discriminating
The RMS surface gap samples only the outer boundary of the reference and evaluates the union there. That makes it blind to the inside: two slabs that stop just short of each other leave a thin internal slit, or one slab can sit slightly rotated with a wedge gap — yet the outer surface still matches the box, so the loss is near zero while the assembly is wrong.
The volumetric Dice term fixes this. It integrates a smooth Heaviside occupancy H(−f) = 1/(1 + e^{β f}) over a grid and compares the solid of the assembled union against the solid of the reference box:
\[L_{\text{dice}} = 1 - \frac{2\,|A \cap B|}{|A| + |B|},\]
which is 0 if and only if the two solids coincide. An internal slit shrinks |A ∩ B| (and hence raises the loss) exactly as much as an overlap or a piece sticking out — so a mis-arranged assembly can no longer hide behind a small surface gap. Jaccard (|A ∩ B| / |A ∪ B|) is the stricter alternative.
7. Watching the reassembly (animated viewer)
res_hh.trace holds the pose at every LBFGS iteration; assemble_frames meshes the union at each one. Passing that frame sequence straight to vtk gives an animated viewer with a playback bar — press play (or drag the scrub slider) to watch the four slabs slide back into the box. We keep the animation light by subsampling the trace to ten evenly-spaced poses:
hist = [copy(θ) for θ in res_hh.trace]
println("iterations recorded = ", length(hist))
println("initial loss = ", round(loss(res_hh.trace[1]); digits=6))
println("final loss = ", round(loss(res_hh.trace[end]); digits=6))
idxs = round.(Int, range(1, length(hist); length=10))
frames = @time assemble_frames(pieces, hist[idxs]; n=40)
vtk(frames; times=collect(1:length(frames)), datafile="reassembly.pvd")iterations recorded = 30
initial loss = 0.689744
final loss = 0.175904
2.515778 seconds (7.05 M allocations: 387.074 MiB, 3.61% gc time, 54.28% compilation time)The first and last frames are the shuffled start and the reassembled box:
vtk(frames[1])vtk(frames[end])The datafile keyword streams the frames to disk (one .vtp per frame plus a .pvd you can open in ParaView) instead of inlining them; you can also call write_pvd(frames, path) directly.
The built-in register and assemble
The hand-written loop above spells out everything the two built-in drivers wrap. register recovers the single rigid pose that moves one shape onto another, and assemble recovers the stacked poses of n pieces — both quasi-Newton + AD over the same surface/field distances, returning structured results:
# single-object registration: move `ref` onto `src`
src = sphere(1.0)
ref = translate(sphere(1.0), 0.4, 0.0, 0.0)
reg = register(src, ref, zeros(6); store_trace=true) # a RigidAlignment
println("recovered translation: ", round.(reg.translation; digits=3))
println("final loss: ", round(reg.final_loss; digits=6))
println("pose history frames: ", length(optimization_history(reg)))recovered translation: [-0.4, -0.0, -0.0]
final loss: 0.0
pose history frames: 300register returns a RigidAlignment (translation, rotation, objective, final_loss, and the optimizer result); optimization_history reads the optimizer’s stored trace back as the sequence of pose vectors it visited (the same θ list assemble_frames animates).
# multi-body assembly: recover the 4 slab poses in one call
asm = assemble(B, pieces, zeros(24); soft=true, iterations=40) # an AssemblyResult
println("final loss: ", round(asm.final_loss; digits=4))
println("recovered poses: ", length(asm.transforms))final loss: 0.1943
recovered poses: 4assemble returns an AssemblyResult (minimizer, per-piece transforms as (t, R), final_loss, and the optimizer result), and accepts objective = :surface | :dice | :jaccard | :both for the same discriminating losses §6 motivates.
Finally, box_mesh(lo, hi) is a minimal 12-triangle axis-aligned box — the exact, low-poly counterpart of a marching-cubes box — handy as a cheap mesh target for mesh→SDF registration:
bm = box_mesh([-1.0, -1.0, -1.0], [1.0, 1.0, 1.0])
println("box_mesh: ", size(bm.vertices, 2), " vertices, ", size(bm.faces, 2), " triangles")box_mesh: 8 vertices, 12 trianglesDeclarative component assembly
register/assemble recover poses from data. When you already know the placement, build a declarative assembly instead: a Component is a named model instance with a rigid Placement (world = R·local + t), and an Assembly is a set of components whose union is one model:
bracket = assembly(
"Bracket",
component("Base", box(0.100, 0.050, 0.050)), # 100 × 50 × 50 mm block
component("Post", box(0.050, 0.050, 0.100); translate = [0.0, 0.0, 0.075]),
component("Pin", sphere(0.0125); translate = [0.0, 0.0, 0.150]), # Ø25 pivot pin
)
println(bracket)
println("assembly bounds : ", round.(bounds(bracket)[1] .* 1000; digits = 1), " .. ",
round.(bounds(bracket)[2] .* 1000; digits = 1), " mm")Assembly(Bracket: 3 components)
assembly bounds : [-50.0, -25.0, -25.0] .. [50.0, 25.0, 162.5] mmPlacements update functionally (translate/rotate return a new component), and pairwise interference (overlap volume) and clearance (minimum surface gap) reuse the §57 SDF overlap measures:
touch = assembly(
"Touch",
component("A", box(2.0, 2.0, 2.0)),
component("B", box(2.0, 2.0, 2.0); translate = [0.5, 0.0, 0.0]),
)
println("interference volume : ", round(interference(touch, "A", "B"; n = 64, k = 40.0); digits = 2))
pair = assembly(
"Pair",
component("S1", sphere(1.0)),
component("S2", sphere(1.0); translate = [3.0, 0.0, 0.0]),
)
println("minimum clearance : ", round(clearance(pair, "S1", "S2"; n = 96); digits = 3))interference volume : 5.67 m^3
minimum clearance : 1.0 mNext steps
- Assembly Optimization Strategies — robust alternatives for scattered pieces: volumetric sampling, recentering, and stochastic Adam.
- Segmentation: Thresholds, Smooth Fields, and SDFs — where the overlap measures first appear.
- Differentiable Simulation and Optimization — the gradient-based driver
minimizewraps.