Improving Physical Object Reconstruction with Bending Models

# Illustrative sketch (not from the paper)
import torch

# 1. dense surface (placeholder tensor) from 3D Gaussian splatting
surface = torch.randn(1, 3, 1000, requires_grad=False)  # [B,3,N]

# 2. discretize to graph G = {V,E,T}
V = torch.randn(100, 3, requires_grad=True)               # mass nodes
E = [(i, i+1) for i in range(99)]                         # axial springs
T = [(i, i+1, i+2) for i in range(98)]                    # bending triplets (O(N))

# 3. mechanical parameters
k_axial = torch.full((len(E),), 1.0, requires_grad=True)
γ_axial = torch.full((len(E),), 0.1, requires_grad=True)
k_bend = torch.full((len(T),), 0.5, requires_grad=True)
c_bend = torch.full((len(T),), 0.05, requires_grad=True)

# 4. energy functions
def axial_energy(v, e, k):
    i, j = e
    return 0.5 * k * (v[i] - v[j]).norm()**2

def bend_energy(v, t, k):
    i, j, k_idx = t
    # angle between (i-j) and (j-k)
    a = v[i] - v[j]
    b = v[k_idx] - v[j]
    cos_theta = (a @ b) / (a.norm() * b.norm() + 1e-8)
    theta = torch.acos(torch.clamp(cos_theta, -1.0, 1.0))
    return 0.5 * k * theta**2

# 5. compute forces
f = torch.zeros_like(V)
for idx, e in enumerate(E):
    f[e[0]] += -torch.autograd.grad(axial_energy(V, e, k_axial[idx]), V, retain_graph=True)[0][e[0]]
    f[e[1]] += -torch.autograd.grad(axial_energy(V, e, k_axial[idx]), V, retain_graph=True)[0][e[1]]
for idx, t in enumerate(T):
    f[t[0]] += -torch.autograd.grad(bend_energy(V, t, k_bend[idx]), V, retain_graph=True)[0][t[0]]
    f[t[1]] += -torch.autograd.grad(bend_energy(V, t, k_bend[idx]), V, retain_graph=True)[0][t[1]]
    f[t[2]] += -torch.autograd.grad(bend_energy(V, t, k_bend[idx]), V, retain_graph=True)[0][t[2]]

# 6. explicit Euler step (differentiable)
Δt = 0.01
V_next = V + Δt * f

# 7. loss (Chamfer placeholder + tracking placeholder)
L_geo = torch.tensor(0.0)  # replace with actual Chamfer
L_trk = torch.nn.functional.mse_loss(V_next, V)  # dummy tracking
loss = L_geo + L_trk

# 8. optimize parameters
opt = torch.optim.Adam([k_axial, γ_axial, k_bend, c_bend], lr=1e-3)
opt.zero_grad()
loss.backward()
opt.step()

// Illustrative sketch (not from the paper)
const tf = require('@tensorflow/tfjs-node');

// 1. dense surface placeholder (tensor of 3D points)
const surface = tf.randomNormal([1, 3, 1000]); // [B,3,N]

// 2. discretize to graph G = {V,E,T}
let V = tf.variable(tf.randomNormal([100, 3])); // mass nodes
const E = Array.from({length: 99}, (_,i) => [i, i+1]); // axial springs
const T = Array.from({length: 98}, (_,i) => [i, i+1, i+2]); // bending triplets (O(N))

// 3. mechanical parameters (trainable scalars)
let kAxial = tf.variable(tf.fill([E.length], 1.0));
let gammaAxial = tf.variable(tf.fill([E.length], 0.1));
let kBend = tf.variable(tf.fill([T.length], 0.5));
let cBend = tf.variable(tf.fill([T.length], 0.05));

// 4. energy helpers
function axialEnergy(v, e, k) {
  const vi = v.gather([e[0]]);
  const vj = v.gather([e[1]]);
  return tf.mul(0.5 * k, tf.norm(tf.sub(vi, vj)).square());
}
function bendEnergy(v, t, k) {
  const vi = v.gather([t[0]]);
  const vj = v.gather([t[1]]);
  const vk = v.gather([t[2]]);
  const a = tf.sub(vi, vj);
  const b = tf.sub(vk, vj);
  const cosTheta = tf.div(tf.sum(tf.mul(a, b)), tf.mul(tf.norm(a), tf.norm(b)).add(1e-8));
  const theta = tf.acos(tf.clipByValue(cosTheta, -1, 1));
  return tf.mul(0.5 * k, tf.square(theta));
}

// 5. compute forces (gradient of energies)
function computeForces(V) {
  return tf.tidy(() => {
    let f = tf.zerosLike(V);
    E.forEach((e, idx) => {
      const energy = axialEnergy(V, e, kAxial.arraySync()[idx]);
      const grads = tf.grad(v => axialEnergy(v, e, kAxial.arraySync()[idx]))(V);
      f = tf.add(f, grads);
    });
    T.forEach((t, idx) => {
      const energy = bendEnergy(V, t, kBend.arraySync()[idx]);
      const grads = tf.grad(v => bendEnergy(v, t, kBend.arraySync()[idx]))(V);
      f = tf.add(f, grads);
    });
    return f.neg(); // f = -∇E
  });
}

// 6. explicit Euler integration (differentiable)
const dt = 0.01;
const forces = computeForces(V);
const Vnext = tf.add(V, tf.mul(dt, forces));

// 7. loss (placeholder Chamfer + tracking)
const Lgeo = tf.scalar(0.0); // replace with actual Chamfer implementation
const Ltrk = tf.losses.meanSquaredError(Vnext, V); // dummy tracking loss
const loss = tf.add(Lgeo, Ltrk);

// 8. optimizer step
const optimizer = tf.train.adam(1e-3);
optimizer.minimize(() => loss, true, [kAxial, gammaAxial, kBend, cBend]);