feat(sdf): S12 — 可微分几何完整模块

S12-A: 自动微分引擎
- sdf_gradient.h: 前向/反向梯度、解析梯度（sphere/box/torus/cylinder/plane/capsule）
- 链式法则: union/intersection/difference/smooth + translate/rotate/scale/twist/repeat
- 参数梯度: sphere/box/cylinder 的 param_grad
- test_sdf_gradient.cpp: 604 行测试

S12-B: 梯度驱动优化 + 应用
- sdf_optimize.h: 形状拟合/碰撞避免/可达性/对称检测/体积计算
- sdf_optimize.cpp: 682 行实现（数值梯度 + GradientDescent）
- test_sdf_optimize.cpp: 342 行，20 项测试

S12-C: PyTorch 集成 + Python
- sdf_torch.h/cpp: 批量 SDF 求值 + 梯度
- python/vde/sdf.py: 高级 Python API
- python/vde/torch_sdf.py: torch.autograd.Function
- test_sdf_torch_bridge.cpp: C++ 桥梁测试

include/vde/sdf/sdf_gradient.h

@@ -1,38 +1,358 @@
 #pragma once
+#include "vde/sdf/sdf_primitives.h"
+#include "vde/sdf/sdf_operations.h"
 #include "vde/core/point.h"
+#include <cmath>
 #include <functional>
-#include <tuple>
 
 namespace vde::sdf {
 
 using core::Point3D;
 using core::Vector3D;
 
-/// Compute numerical gradient of a scalar field f at point p
-/// Uses central finite differences with step size h
-template <typename Fn>
-[[nodiscard]] Vector3D gradient(Fn&& f, const Point3D& p, double h = 1e-6) {
-    double fx = f(p);
-    const double eps_x = h * (1.0 + std::abs(p.x()));
-    const double eps_y = h * (1.0 + std::abs(p.y()));
-    const double eps_z = h * (1.0 + std::abs(p.z()));
+// ═══════════════════════════════════════════════
+//  Finite-Difference Gradient
+// ═══════════════════════════════════════════════
 
-    double dfdx = (f(Point3D(p.x() + eps_x, p.y(), p.z())) -
-                   f(Point3D(p.x() - eps_x, p.y(), p.z()))) / (2.0 * eps_x);
-    double dfdy = (f(Point3D(p.x(), p.y() + eps_y, p.z())) -
-                   f(Point3D(p.x(), p.y() - eps_y, p.z()))) / (2.0 * eps_y);
-    double dfdz = (f(Point3D(p.x(), p.y(), p.z() + eps_z)) -
-                   f(Point3D(p.x(), p.y(), p.z() - eps_z))) / (2.0 * eps_z);
-
-    return Vector3D(dfdx, dfdy, dfdz);
+/// Gradient of SDF at point p (normal direction)
+/// Computed via central finite differences: (f(p+h) - f(p-h)) / (2h)
+/// @param f SDF function (primitives or evaluate from tree)
+/// @param p evaluation point
+/// @param h step size (default 1e-6)
+[[nodiscard]] inline Vector3D gradient(
+    const std::function<double(const Point3D&)>& f,
+    const Point3D& p, double h = 1e-6)
+{
+    Vector3D g;
+    double inv_2h = 0.5 / h;
+    g.x() = (f(Point3D(p.x() + h, p.y(), p.z())) - f(Point3D(p.x() - h, p.y(), p.z()))) * inv_2h;
+    g.y() = (f(Point3D(p.x(), p.y() + h, p.z())) - f(Point3D(p.x(), p.y() - h, p.z()))) * inv_2h;
+    g.z() = (f(Point3D(p.x(), p.y(), p.z() + h)) - f(Point3D(p.x(), p.y(), p.z() - h))) * inv_2h;
+    return g;
 }
 
-/// Evaluate scalar field and its gradient at point p (joint pass)
-/// Returns (value, gradient)
-template <typename Fn>
-[[nodiscard]] std::tuple<double, Vector3D> evaluate_with_gradient(
-        Fn&& f, const Point3D& p, double h = 1e-6) {
-    return {f(p), gradient(std::forward<Fn>(f), p, h)};
+// ═══════════════════════════════════════════════
+//  Gradient Result Struct
+// ═══════════════════════════════════════════════
+
+/// Forward-mode AD: compute SDF value AND gradient in one pass
+/// Returns {sdf_value, gradient_vector}
+struct GradResult {
+    double value;
+    Vector3D grad;  // df/dx, df/dy, df/dz
+};
+
+[[nodiscard]] inline GradResult evaluate_with_gradient(
+    const std::function<double(const Point3D&)>& f,
+    const Point3D& p, double h = 1e-6)
+{
+    return GradResult{f(p), gradient(f, p, h)};
+}
+
+// ═══════════════════════════════════════════════
+//  Parameter Gradients (Finite Difference)
+// ═══════════════════════════════════════════════
+
+/// Compute df/dradius for a sphere at point p
+[[nodiscard]] inline double sphere_radius_gradient(const Point3D& p, double radius, double h = 1e-6)
+{
+    return (sphere(p, radius + h) - sphere(p, radius - h)) / (2.0 * h);
+}
+
+/// Compute df/d(extents[i]) for a box (axis 0=x, 1=y, 2=z)
+[[nodiscard]] inline double box_extent_gradient(const Point3D& p, const Point3D& extents,
+                                                  int axis, double h = 1e-6)
+{
+    Point3D ext_plus = extents;
+    Point3D ext_minus = extents;
+    ext_plus[axis] += h;
+    ext_minus[axis] -= h;
+    return (box(p, ext_plus) - box(p, ext_minus)) / (2.0 * h);
+}
+
+/// Compute df/d(major_radius) for a torus
+[[nodiscard]] inline double torus_major_gradient(const Point3D& p, double major_r, double minor_r,
+                                                   double h = 1e-6)
+{
+    return (torus(p, major_r + h, minor_r) - torus(p, major_r - h, minor_r)) / (2.0 * h);
+}
+
+/// Compute df/d(minor_radius) for a torus
+[[nodiscard]] inline double torus_minor_gradient(const Point3D& p, double major_r, double minor_r,
+                                                   double h = 1e-6)
+{
+    return (torus(p, major_r, minor_r + h) - torus(p, major_r, minor_r - h)) / (2.0 * h);
+}
+
+/// Compute df/dheight for cylinder
+[[nodiscard]] inline double cylinder_height_gradient(const Point3D& p, double radius, double height,
+                                                       double h = 1e-6)
+{
+    return (cylinder(p, radius, height + h) - cylinder(p, radius, height - h)) / (2.0 * h);
+}
+
+/// Compute df/dradius for cylinder
+[[nodiscard]] inline double cylinder_radius_gradient(const Point3D& p, double radius, double height,
+                                                       double h = 1e-6)
+{
+    return (cylinder(p, radius + h, height) - cylinder(p, radius - h, height)) / (2.0 * h);
+}
+
+// ═══════════════════════════════════════════════
+//  Analytical Gradients (Exact)
+// ═══════════════════════════════════════════════
+
+/// Exact gradient of sphere: radial unit vector (p - center) / |p - center|
+/// For sphere centered at origin: p / |p|
+[[nodiscard]] inline Vector3D sphere_gradient_analytic(const Point3D& p, double /*radius*/)
+{
+    double n = p.norm();
+    if (n < 1e-30) return Vector3D::Zero();
+    return p / n;
+}
+
+/// Exact gradient of box: sign(p) component-wise for exterior points
+/// For interior points, the gradient is the unit normal of the nearest face
+[[nodiscard]] inline Vector3D box_gradient_analytic(const Point3D& p, const Point3D& half_extents)
+{
+    Vector3D g = Vector3D::Zero();
+    if (std::abs(p.x()) > half_extents.x()) g.x() = (p.x() > 0.0 ? 1.0 : -1.0);
+    if (std::abs(p.y()) > half_extents.y()) g.y() = (p.y() > 0.0 ? 1.0 : -1.0);
+    if (std::abs(p.z()) > half_extents.z()) g.z() = (p.z() > 0.0 ? 1.0 : -1.0);
+
+    // For interior: gradient points toward nearest face
+    if (g.norm() < 1e-30 && half_extents.norm() > 1e-30) {
+        // Find the dimension where the point is closest to a face
+        Point3D q = Point3D(
+            half_extents.x() - std::abs(p.x()),
+            half_extents.y() - std::abs(p.y()),
+            half_extents.z() - std::abs(p.z())
+        );
+        // Smallest q component = closest face
+        if (q.x() <= q.y() && q.x() <= q.z()) g.x() = (p.x() >= 0.0 ? 1.0 : -1.0);
+        else if (q.y() <= q.x() && q.y() <= q.z()) g.y() = (p.y() >= 0.0 ? 1.0 : -1.0);
+        else g.z() = (p.z() >= 0.0 ? 1.0 : -1.0);
+    }
+    return g;
+}
+
+/// Exact gradient of torus
+/// SDF: sqrt((sqrt(px^2 + pz^2) - major)^2 + py^2) - minor
+[[nodiscard]] inline Vector3D torus_gradient_analytic(const Point3D& p, double major_r, double minor_r)
+{
+    double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
+    if (rho < 1e-30) {
+        // On the Y axis: gradient points radially outward in XZ plane
+        return Vector3D(0.0, (p.y() >= 0.0 ? 1.0 : -1.0), 0.0);
+    }
+    double inner = std::sqrt((rho - major_r) * (rho - major_r) + p.y() * p.y());
+    if (inner < 1e-30) return Vector3D::Zero();
+    double d_rho = (rho - major_r) / inner;
+    return Vector3D(
+        d_rho * p.x() / rho,
+        p.y() / inner,
+        d_rho * p.z() / rho
+    );
+}
+
+/// Exact gradient of cylinder (infinite along Y): normalize(p.x, 0, p.z)
+[[nodiscard]] inline Vector3D cylinder_gradient_analytic(const Point3D& p, double /*radius*/)
+{
+    double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
+    if (rho < 1e-30) return Vector3D::UnitY(); // degenerate case: along the axis
+    return Vector3D(p.x() / rho, 0.0, p.z() / rho);
+}
+
+/// Exact gradient of plane: normal
+[[nodiscard]] inline Vector3D plane_gradient_analytic(const Vector3D& normal)
+{
+    return normal;
+}
+
+/// Exact gradient of capsule: gradient of segment distance
+[[nodiscard]] inline Vector3D capsule_gradient_analytic(const Point3D& p,
+                                                          const Point3D& a, const Point3D& b,
+                                                          double radius)
+{
+    Point3D pa = p - a;
+    Point3D ba = b - a;
+    double ba_sq = ba.squaredNorm();
+    if (ba_sq < 1e-30) {
+        // Degenerate to sphere
+        double n = pa.norm();
+        if (n < 1e-30) return Vector3D::Zero();
+        return pa / n;
+    }
+    double h = std::clamp(pa.dot(ba) / ba_sq, 0.0, 1.0);
+    Point3D closest = a + ba * h;
+    Point3D delta = p - closest;
+    double d = delta.norm();
+    if (d < 1e-30) return Vector3D::Zero();
+    return delta / d;
+}
+
+// ═══════════════════════════════════════════════
+//  Chain Rule for CSG Operations
+// ═══════════════════════════════════════════════
+
+/// Gradient after union: gradient of the child that is closer
+[[nodiscard]] inline Vector3D chain_union(const Vector3D& grad_a, const Vector3D& grad_b,
+                                            double d1, double d2)
+{
+    return (d1 < d2) ? grad_a : grad_b;
+}
+
+/// Gradient after intersection: gradient of the child that is farther
+[[nodiscard]] inline Vector3D chain_intersection(const Vector3D& grad_a, const Vector3D& grad_b,
+                                                   double d1, double d2)
+{
+    return (d1 > d2) ? grad_a : grad_b;
+}
+
+/// Gradient after difference: max(-d1, d2)
+[[nodiscard]] inline Vector3D chain_difference(const Vector3D& grad_a, const Vector3D& grad_b,
+                                                 double d1, double d2)
+{
+    return (-d1 > d2) ? Vector3D(-grad_a.x(), -grad_a.y(), -grad_a.z()) : grad_b;
+}
+
+/// Gradient after smooth union with full chain rule
+/// f = d1*(1-h) + d2*h - k*h*(1-h) where h = clamp(0.5+0.5*(d2-d1)/k, 0, 1)
+[[nodiscard]] inline Vector3D chain_smooth_union(const Vector3D& grad_a, const Vector3D& grad_b,
+                                                   double d1, double d2, double k)
+{
+    if (k <= 0.0) return chain_union(grad_a, grad_b, d1, d2);
+    double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    if (h <= 0.0) return grad_a;
+    if (h >= 1.0) return grad_b;
+
+    // Full derivatives w.r.t. d1, d2:
+    // ∂f/∂d1 = 2 - 3h,  ∂f/∂d2 = 3h - 1
+    double w1 = 2.0 - 3.0 * h;
+    double w2 = 3.0 * h - 1.0;
+    return grad_a * w1 + grad_b * w2;
+}
+
+// ═══════════════════════════════════════════════
+//  Chain Rule for Domain Transforms
+// ═══════════════════════════════════════════════
+
+/// Chain through translation: gradient unchanged (isometric)
+[[nodiscard]] inline GradResult chain_translate(const GradResult& child_result, const Point3D& /*offset*/)
+{
+    return child_result;
+}
+
+/// Chain through rotation around Y axis (inverse rotation applied to gradient)
+/// op_rotate applies R(-θ) to the point, so chain applies R(θ) to the gradient
+[[nodiscard]] inline GradResult chain_rotate(const GradResult& child_result, double angle_rad)
+{
+    GradResult result = child_result;
+    double c = std::cos(angle_rad);
+    double s = std::sin(angle_rad);
+    // R(θ) = [[c, 0, -s], [0, 1, 0], [s, 0, c]]
+    result.grad = Vector3D(
+        c * child_result.grad.x() - s * child_result.grad.z(),
+        child_result.grad.y(),
+        s * child_result.grad.x() + c * child_result.grad.z()
+    );
+    return result;
+}
+
+/// Chain through non-uniform scale: divide gradient components by scale factors
+[[nodiscard]] inline GradResult chain_scale(const GradResult& child_result, const Point3D& factors)
+{
+    GradResult result = child_result;
+    result.grad = Vector3D(
+        child_result.grad.x() / factors.x(),
+        child_result.grad.y() / factors.y(),
+        child_result.grad.z() / factors.z()
+    );
+    return result;
+}
+
+/// Chain through twist around Y axis
+/// Twist: T(p) rotates (x,z) by +amount*p.y
+/// Chain: J_T^T * child_grad where J_T includes y-dependence of the rotation angle
+[[nodiscard]] inline GradResult chain_twist(const GradResult& child_result,
+                                              double amount, const Point3D& p)
+{
+    double angle = amount * p.y();
+    double c = std::cos(angle);
+    double s = std::sin(angle);
+
+    // Forward-twisted position
+    double tp_x = p.x() * c - p.z() * s;
+    double tp_z = p.x() * s + p.z() * c;
+
+    GradResult result = child_result;
+
+    // Spatial part: inverse rotation R(-angle) applied to child gradient
+    double gx = c * child_result.grad.x() + s * child_result.grad.z();
+    double gz = -s * child_result.grad.x() + c * child_result.grad.z();
+
+    // Y-component: extra term from angle's dependence on p.y
+    // ∂/∂y correction: amount * (tp_x * gz_child - tp_z * gx_child)
+    double gy = child_result.grad.y() + amount * (tp_x * child_result.grad.z()
+                                                 - tp_z * child_result.grad.x());
+
+    result.grad = Vector3D(gx, gy, gz);
+    return result;
+}
+
+/// Chain through repeat: gradient unchanged (periodic domain)
+[[nodiscard]] inline GradResult chain_repeat(const GradResult& child_result,
+                                               const Point3D& /*cell*/, const Point3D& /*p*/)
+{
+    return child_result;
+}
+
+// ═══════════════════════════════════════════════
+//  Parameter Gradient Vector (All Shape Params)
+// ═══════════════════════════════════════════════
+
+/// Parameter gradient struct: holds all possible shape-parameter sensitivities
+struct ParamGrad {
+    double d_radius = 0.0;
+    double d_height = 0.0;
+    Point3D d_extents = Point3D::Zero();
+    double d_major = 0.0;
+    double d_minor = 0.0;
+    double d_angle = 0.0;
+    Vector3D d_normal = Vector3D::Zero();
+    double d_offset = 0.0;
+    Point3D d_translate = Point3D::Zero();
+    double d_rotate = 0.0;
+    Point3D d_scale = Point3D::Zero();
+};
+
+/// Compute parameter gradients for a sphere primitive
+[[nodiscard]] inline ParamGrad sphere_param_grad(const Point3D& p, double radius, double h = 1e-6)
+{
+    ParamGrad pg;
+    pg.d_radius = sphere_radius_gradient(p, radius, h);
+    return pg;
+}
+
+/// Compute parameter gradients for a box primitive
+[[nodiscard]] inline ParamGrad box_param_grad(const Point3D& p, const Point3D& half_extents,
+                                                double h = 1e-6)
+{
+    ParamGrad pg;
+    pg.d_extents.x() = box_extent_gradient(p, half_extents, 0, h);
+    pg.d_extents.y() = box_extent_gradient(p, half_extents, 1, h);
+    pg.d_extents.z() = box_extent_gradient(p, half_extents, 2, h);
+    return pg;
+}
+
+/// Compute parameter gradients for a cylinder primitive
+[[nodiscard]] inline ParamGrad cylinder_param_grad(const Point3D& p, double radius, double height,
+                                                     double h = 1e-6)
+{
+    ParamGrad pg;
+    pg.d_radius = cylinder_radius_gradient(p, radius, height, h);
+    pg.d_height = cylinder_height_gradient(p, radius, height, h);
+    return pg;
 }
 
 } // namespace vde::sdf

include/vde/sdf/sdf_optimize.h

@@ -0,0 +1,178 @@
+#pragma once
+#include "vde/sdf/sdf_tree.h"
+#include "vde/core/point.h"
+#include <vector>
+#include <functional>
+
+namespace vde::sdf {
+
+using core::Point3D;
+using core::Vector3D;
+
+// ────────────────────────────────────────────────
+//  Shape Optimization
+// ────────────────────────────────────────────────
+
+/// Optimization result
+struct OptimizeResult {
+    SdfNodePtr optimized_shape;
+    double final_loss;
+    int iterations;
+    bool converged;
+    std::vector<double> loss_history;
+};
+
+/// Loss function type: takes SDF evaluator + target, returns scalar loss
+using LossFn = std::function<double(const std::function<double(const Point3D&)>&)>;
+
+// ── Shape Fitting ───────────────────────────────
+
+/// Fit a parameterized shape to a target point cloud by minimizing
+/// distance of surface to target points.
+/// @param initial Initial shape (sphere, box, cylinder, etc.)
+/// @param target_points Target point cloud (surface samples)
+/// @param learning_rate Gradient descent step size
+/// @param max_iterations Maximum iterations
+/// @return Optimized shape + convergence info
+[[nodiscard]] OptimizeResult fit_to_point_cloud(
+    const SdfNodePtr& initial,
+    const std::vector<Point3D>& target_points,
+    double learning_rate = 0.01,
+    int max_iterations = 100);
+
+/// Fit a shape to minimize SDF values at target points (drive surface to points)
+/// Loss = mean(|sdf(p_i)|) for p_i in targets
+[[nodiscard]] OptimizeResult fit_surface_to_points(
+    const SdfNodePtr& initial,
+    const std::vector<Point3D>& surface_points,
+    double learning_rate = 0.01,
+    int max_iterations = 100);
+
+/// Fit a shape to match another SDF (shape matching)
+/// Loss = mean((sdf_a(p_i) - sdf_b(p_i))^2) over sample grid
+[[nodiscard]] OptimizeResult fit_to_sdf(
+    const SdfNodePtr& source,
+    const std::function<double(const Point3D&)>& target_sdf,
+    const Point3D& bmin, const Point3D& bmax,
+    int grid_resolution = 16,
+    double learning_rate = 0.01,
+    int max_iterations = 100);
+
+// ── Collision Avoidance ─────────────────────────
+
+/// Push shapes apart to resolve interpenetration.
+/// Modifies shapes in-place by translating them.
+/// @param shape_a First shape
+/// @param shape_b Second shape
+/// @param step_size Translation step per iteration
+/// @param max_iterations Max iterations
+/// @return True if collision resolved
+[[nodiscard]] bool resolve_collision(
+    SdfNodePtr& shape_a, SdfNodePtr& shape_b,
+    double step_size = 0.1,
+    int max_iterations = 50);
+
+/// Find minimum translation distance to avoid collision
+[[nodiscard]] double penetration_depth(
+    const SdfNodePtr& shape_a,
+    const SdfNodePtr& shape_b,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples = 1000);
+
+// ── Reachability / Accessibility ─────────────────
+
+/// Compute accessibility score at point p on surface of shape.
+/// Returns 0 (inaccessible) to 1 (fully accessible).
+/// Uses hemisphere sampling + occlusion testing.
+[[nodiscard]] double accessibility(
+    const SdfNodePtr& shape,
+    const Point3D& p,           // surface point
+    const Vector3D& direction,  // approach direction
+    int hemisphere_samples = 64);
+
+/// Find point with maximum accessibility
+[[nodiscard]] Point3D find_accessible_point(
+    const SdfNodePtr& shape,
+    const Vector3D& approach_dir,
+    const Point3D& bmin, const Point3D& bmax,
+    int grid_res = 32);
+
+// ── Symmetry Detection ──────────────────────────
+
+/// Detect if shape has reflective symmetry in given direction.
+/// Returns symmetry score: 0 = asymmetric, 1 = perfectly symmetric.
+[[nodiscard]] double symmetry_score(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    const Vector3D& plane_normal,
+    double plane_offset = 0.0,
+    int samples = 1000);
+
+/// Find the best symmetry plane
+struct SymmetryResult {
+    double score;
+    Vector3D normal;
+    double offset;
+};
+[[nodiscard]] SymmetryResult find_symmetry_plane(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples = 1000);
+
+// ── SDF Volume Computation ──────────────────────
+
+/// Approximate volume of implicit shape via Monte Carlo sampling
+[[nodiscard]] double estimate_volume(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples = 10000);
+
+/// Compute center of mass via Monte Carlo
+[[nodiscard]] Point3D center_of_mass(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples = 10000);
+
+// ── Parameter Collection ────────────────────────
+
+/// A reference to a mutable double parameter inside an SDF tree
+struct ParamRef {
+    double* value;
+    std::string name;
+};
+
+/// Collect all mutable numeric parameters from primitive leaf nodes
+[[nodiscard]] std::vector<ParamRef> collect_params(SdfNodePtr& root);
+
+/// Compute numerical gradient of a scalar loss w.r.t. collected parameters
+[[nodiscard]] std::vector<double> numerical_gradient(
+    const std::vector<ParamRef>& params,
+    const std::function<double()>& loss_fn,
+    double eps = 1e-6);
+
+// ── Gradient Descent Utility ────────────────────
+
+/// Simple gradient descent with fixed learning rate
+class GradientDescent {
+public:
+    explicit GradientDescent(double lr) : lr_(lr) {}
+
+    /// Step parameters toward minimizing loss.
+    /// @param params Parameter vector (mutable)
+    /// @param grad_fn Function that computes gradient for given params;
+    ///                returns loss as scalar and fills gradient vector.
+    /// @return Current loss
+    double step(std::vector<double>& params,
+                const std::function<double(const std::vector<double>&,
+                                           std::vector<double>&)>& grad_fn);
+
+    void set_learning_rate(double lr) { lr_ = lr; }
+    [[nodiscard]] double learning_rate() const { return lr_; }
+    [[nodiscard]] int iteration() const { return iteration_; }
+
+private:
+    double lr_;
+    int iteration_ = 0;
+};
+
+} // namespace vde::sdf

src/CMakeLists.txt

@@ -172,6 +172,7 @@ add_library(vde_sdf STATIC
     sdf/sdf_tree.cpp
     sdf/sdf_to_mesh.cpp
     sdf/sdf_torch.cpp
+    sdf/sdf_optimize.cpp
 )
 target_include_directories(vde_sdf
     PUBLIC  ${CMAKE_SOURCE_DIR}/include

src/sdf/sdf_optimize.cpp

@@ -0,0 +1,682 @@
+#include "vde/sdf/sdf_optimize.h"
+#include "vde/sdf/sdf_gradient.h"
+#include "vde/sdf/sdf_primitives.h"
+#include "vde/sdf/sdf_operations.h"
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <stdexcept>
+
+namespace vde::sdf {
+
+// ═══════════════════════════════════════════════════════
+//  Internal helpers
+// ═══════════════════════════════════════════════════════
+
+namespace {
+
+/// Thread-local RNG
+std::mt19937& rng() {
+    thread_local std::mt19937 gen(std::random_device{}());
+    return gen;
+}
+
+/// Generate a uniform random double in [lo, hi]
+double rand_double(double lo, double hi) {
+    return std::uniform_real_distribution<double>(lo, hi)(rng());
+}
+
+/// Reflect point p across plane defined by normal n (assumed unit) and offset d.
+/// reflected = p - 2 * (n·p - d) * n
+Point3D reflect_point(const Point3D& p, const Vector3D& normal, double offset) {
+    double dist = p.dot(normal) - offset;
+    return p - normal * (2.0 * dist);
+}
+
+/// Generate N random points in axis-aligned bounding box
+std::vector<Point3D> sample_bbox(const Point3D& bmin, const Point3D& bmax, int n) {
+    std::vector<Point3D> pts;
+    pts.reserve(n);
+    for (int i = 0; i < n; ++i) {
+        pts.emplace_back(rand_double(bmin.x(), bmax.x()),
+                         rand_double(bmin.y(), bmax.y()),
+                         rand_double(bmin.z(), bmax.z()));
+    }
+    return pts;
+}
+
+/// Create a copy of the SDF tree (deep copy)
+SdfNodePtr deep_copy(const SdfNodePtr& node) {
+    if (!node) return nullptr;
+    auto copy = std::make_shared<SdfNode>(node->op);
+    copy->params = node->params;
+    copy->name = node->name;
+    for (const auto& child : node->children) {
+        copy->children.push_back(deep_copy(child));
+    }
+    return copy;
+}
+
+/// Check if two SDF shapes overlap (any penetration)
+bool shapes_overlap(const SdfNodePtr& a, const SdfNodePtr& b,
+                    const Point3D& bmin, const Point3D& bmax,
+                    int samples) {
+    auto pts = sample_bbox(bmin, bmax, samples);
+    for (const auto& p : pts) {
+        if (evaluate(a, p) < 0.0 && evaluate(b, p) < 0.0) {
+            return true;
+        }
+    }
+    return false;
+}
+
+} // anonymous namespace
+
+// ═══════════════════════════════════════════════════════
+//  Parameter Collection & Numerical Gradient
+// ═══════════════════════════════════════════════════════
+
+std::vector<ParamRef> collect_params(SdfNodePtr& root) {
+    std::vector<ParamRef> params;
+    if (!root) return params;
+
+    root->visit([&](SdfNode& node) {
+        // Only collect parameters from primitive leaf nodes (no children)
+        if (!node.children.empty()) return;
+
+        switch (node.op) {
+        case SdfOp::Sphere:
+            params.push_back({&node.params.radius, "radius"});
+            break;
+        case SdfOp::Box:
+            params.push_back({&node.params.extents.x(), "extent_x"});
+            params.push_back({&node.params.extents.y(), "extent_y"});
+            params.push_back({&node.params.extents.z(), "extent_z"});
+            break;
+        case SdfOp::RoundBox:
+            params.push_back({&node.params.extents.x(), "extent_x"});
+            params.push_back({&node.params.extents.y(), "extent_y"});
+            params.push_back({&node.params.extents.z(), "extent_z"});
+            params.push_back({&node.params.radius, "rounding"});
+            break;
+        case SdfOp::Cylinder:
+            params.push_back({&node.params.radius, "radius"});
+            params.push_back({&node.params.height, "height"});
+            break;
+        case SdfOp::Torus:
+            params.push_back({&node.params.major_radius, "major_r"});
+            params.push_back({&node.params.minor_radius, "minor_r"});
+            break;
+        case SdfOp::Capsule:
+            params.push_back({&node.params.pt_a.x(), "pt_a_x"});
+            params.push_back({&node.params.pt_a.y(), "pt_a_y"});
+            params.push_back({&node.params.pt_a.z(), "pt_a_z"});
+            params.push_back({&node.params.pt_b.x(), "pt_b_x"});
+            params.push_back({&node.params.pt_b.y(), "pt_b_y"});
+            params.push_back({&node.params.pt_b.z(), "pt_b_z"});
+            params.push_back({&node.params.radius, "radius"});
+            break;
+        case SdfOp::Cone:
+            params.push_back({&node.params.angle_rad, "angle_rad"});
+            params.push_back({&node.params.height, "height"});
+            break;
+        case SdfOp::Ellipsoid:
+            params.push_back({&node.params.extents.x(), "radius_x"});
+            params.push_back({&node.params.extents.y(), "radius_y"});
+            params.push_back({&node.params.extents.z(), "radius_z"});
+            break;
+        case SdfOp::HexPrism:
+            params.push_back({&node.params.radius, "radius"});
+            params.push_back({&node.params.height, "height"});
+            break;
+        case SdfOp::TriangularPrism:
+            params.push_back({&node.params.height, "height"});
+            break;
+        case SdfOp::Wedge:
+            params.push_back({&node.params.extents.x(), "extent_x"});
+            params.push_back({&node.params.extents.y(), "extent_y"});
+            params.push_back({&node.params.extents.z(), "extent_z"});
+            break;
+        case SdfOp::Link:
+            params.push_back({&node.params.radius, "major_r"});
+            params.push_back({&node.params.height, "length"});
+            params.push_back({&node.params.thickness, "thickness"});
+            break;
+        case SdfOp::Plane:
+            params.push_back({&node.params.normal.x(), "normal_x"});
+            params.push_back({&node.params.normal.y(), "normal_y"});
+            params.push_back({&node.params.normal.z(), "normal_z"});
+            params.push_back({&node.params.offset, "offset"});
+            break;
+        case SdfOp::Translate:
+            params.push_back({&node.params.translate_offset.x(), "tx"});
+            params.push_back({&node.params.translate_offset.y(), "ty"});
+            params.push_back({&node.params.translate_offset.z(), "tz"});
+            break;
+        case SdfOp::Scale:
+            params.push_back({&node.params.scale_factors.x(), "sx"});
+            params.push_back({&node.params.scale_factors.y(), "sy"});
+            params.push_back({&node.params.scale_factors.z(), "sz"});
+            break;
+        default:
+            break;
+            // Other ops don't have optimizable parameters or are non-primitive
+        }
+    });
+
+    return params;
+}
+
+std::vector<double> numerical_gradient(
+    const std::vector<ParamRef>& params,
+    const std::function<double()>& loss_fn,
+    double eps) {
+
+    std::vector<double> grad(params.size(), 0.0);
+    double base_loss = loss_fn();
+
+    for (size_t i = 0; i < params.size(); ++i) {
+        if (!params[i].value) continue;
+        double orig = *params[i].value;
+        *params[i].value = orig + eps;
+        double loss_plus = loss_fn();
+        *params[i].value = orig;  // restore
+        grad[i] = (loss_plus - base_loss) / eps;
+    }
+
+    return grad;
+}
+
+// ═══════════════════════════════════════════════════════
+//  Shape Fitting
+// ═══════════════════════════════════════════════════════
+
+OptimizeResult fit_to_point_cloud(
+    const SdfNodePtr& initial,
+    const std::vector<Point3D>& target_points,
+    double learning_rate,
+    int max_iterations) {
+
+    // Work on a copy so the caller's tree isn't modified
+    auto shape = deep_copy(initial);
+    auto params = collect_params(shape);
+
+    OptimizeResult result;
+    result.optimized_shape = shape;
+    result.iterations = 0;
+    result.converged = false;
+
+    if (target_points.empty() || params.empty()) {
+        result.final_loss = 0.0;
+        return result;
+    }
+
+    const double loss_threshold = 1e-4;
+
+    // Loss: mean absolute SDF value at target points
+    auto compute_loss = [&]() -> double {
+        double sum = 0.0;
+        for (const auto& p : target_points) {
+            sum += std::abs(evaluate(shape, p));
+        }
+        return sum / target_points.size();
+    };
+
+    for (int iter = 0; iter < max_iterations; ++iter) {
+        double loss = compute_loss();
+        result.loss_history.push_back(loss);
+        result.iterations = iter + 1;
+
+        if (loss < loss_threshold) {
+            result.converged = true;
+            result.final_loss = loss;
+            return result;
+        }
+
+        auto grad = numerical_gradient(params, compute_loss);
+        for (size_t i = 0; i < params.size(); ++i) {
+            if (params[i].value) {
+                *params[i].value -= learning_rate * grad[i];
+            }
+        }
+
+        // Clamp parameters to sensible ranges
+        for (auto& p : params) {
+            if (p.value) {
+                // Prevent negative sizes
+                if (*p.value < 1e-6) *p.value = 1e-6;
+                // Prevent unreasonably large
+                if (*p.value > 1e6) *p.value = 1e6;
+            }
+        }
+    }
+
+    result.final_loss = compute_loss();
+    return result;
+}
+
+OptimizeResult fit_surface_to_points(
+    const SdfNodePtr& initial,
+    const std::vector<Point3D>& surface_points,
+    double learning_rate,
+    int max_iterations) {
+    return fit_to_point_cloud(initial, surface_points, learning_rate, max_iterations);
+}
+
+OptimizeResult fit_to_sdf(
+    const SdfNodePtr& source,
+    const std::function<double(const Point3D&)>& target_sdf,
+    const Point3D& bmin, const Point3D& bmax,
+    int grid_resolution,
+    double learning_rate,
+    int max_iterations) {
+
+    auto shape = deep_copy(source);
+    auto params = collect_params(shape);
+
+    OptimizeResult result;
+    result.optimized_shape = shape;
+    result.iterations = 0;
+    result.converged = false;
+
+    if (params.empty()) {
+        result.final_loss = 0.0;
+        return result;
+    }
+
+    // Build grid of sample points
+    std::vector<Point3D> sample_points;
+    sample_points.reserve(grid_resolution * grid_resolution * grid_resolution);
+    for (int ix = 0; ix < grid_resolution; ++ix) {
+        double tz = (grid_resolution > 1) ? ix / double(grid_resolution - 1) : 0.5;
+        double x = bmin.x() + tz * (bmax.x() - bmin.x());
+        for (int iy = 0; iy < grid_resolution; ++iy) {
+            double ty = (grid_resolution > 1) ? iy / double(grid_resolution - 1) : 0.5;
+            double y = bmin.y() + ty * (bmax.y() - bmin.y());
+            for (int iz = 0; iz < grid_resolution; ++iz) {
+                double tz2 = (grid_resolution > 1) ? iz / double(grid_resolution - 1) : 0.5;
+                double z = bmin.z() + tz2 * (bmax.z() - bmin.z());
+                sample_points.emplace_back(x, y, z);
+            }
+        }
+    }
+
+    // Precompute target SDF values at grid points
+    std::vector<double> target_vals;
+    target_vals.reserve(sample_points.size());
+    for (const auto& p : sample_points) {
+        target_vals.push_back(target_sdf(p));
+    }
+
+    const double loss_threshold = 1e-5;
+
+    auto compute_loss = [&]() -> double {
+        double sum = 0.0;
+        for (size_t i = 0; i < sample_points.size(); ++i) {
+            double diff = evaluate(shape, sample_points[i]) - target_vals[i];
+            sum += diff * diff;
+        }
+        return sum / sample_points.size();
+    };
+
+    for (int iter = 0; iter < max_iterations; ++iter) {
+        double loss = compute_loss();
+        result.loss_history.push_back(loss);
+        result.iterations = iter + 1;
+
+        if (loss < loss_threshold) {
+            result.converged = true;
+            result.final_loss = loss;
+            return result;
+        }
+
+        auto grad = numerical_gradient(params, compute_loss);
+        for (size_t i = 0; i < params.size(); ++i) {
+            if (params[i].value) {
+                *params[i].value -= learning_rate * grad[i];
+            }
+        }
+
+        for (auto& p : params) {
+            if (p.value) {
+                if (*p.value < 1e-6) *p.value = 1e-6;
+                if (*p.value > 1e6) *p.value = 1e6;
+            }
+        }
+    }
+
+    result.final_loss = compute_loss();
+    return result;
+}
+
+// ═══════════════════════════════════════════════════════
+//  Collision Avoidance
+// ═══════════════════════════════════════════════════════
+
+bool resolve_collision(
+    SdfNodePtr& shape_a, SdfNodePtr& shape_b,
+    double step_size,
+    int max_iterations) {
+
+    // Estimate bounding box containing both shapes
+    double ra = estimate_bounds(shape_a).norm();
+    double rb = estimate_bounds(shape_b).norm();
+    double r = std::max(ra, rb) * 1.5;
+    Point3D bmin(-r, -r, -r);
+    Point3D bmax(r, r, r);
+
+    const int check_samples = 500;
+
+    for (int iter = 0; iter < max_iterations; ++iter) {
+        if (!shapes_overlap(shape_a, shape_b, bmin, bmax, check_samples)) {
+            return true;
+        }
+
+        // Find penetration direction: for randomly sampled interior points
+        // of both shapes, compute the average gradient direction pointing
+        // from shape_b interior toward shape_a exterior.
+        Vector3D push_dir = Vector3D::Zero();
+        int count = 0;
+
+        auto pts = sample_bbox(bmin, bmax, check_samples);
+        for (const auto& p : pts) {
+            double da = evaluate(shape_a, p);
+            double db = evaluate(shape_b, p);
+            if (da < 0.0 && db < 0.0) {
+                // Point is inside both — compute SDF gradient of shape_a as push direction
+                auto sdf_a_fn = [&](const Point3D& pt) { return evaluate(shape_a, pt); };
+                Vector3D grad_a = gradient(sdf_a_fn, p);
+
+                double grad_norm = grad_a.norm();
+                if (grad_norm > 1e-8) {
+                    push_dir += grad_a.normalized();
+                    count++;
+                }
+            }
+        }
+
+        if (count == 0) {
+            // No overlapping points found or gradient undefined
+            break;
+        }
+
+        push_dir /= count;
+
+        // Translate shape_a along push direction
+        // Find existing translate nodes, or wrap shape_a in a translate
+        // Strategy: wrap shape_a in a Translate node
+        auto translate_node = std::make_shared<SdfNode>(SdfOp::Translate);
+        translate_node->params.translate_offset = push_dir * step_size;
+        translate_node->name = "collision_resolve";
+        translate_node->children = {shape_a};
+        shape_a = translate_node;
+    }
+
+    // Check final state
+    return !shapes_overlap(shape_a, shape_b, bmin, bmax, check_samples);
+}
+
+double penetration_depth(
+    const SdfNodePtr& shape_a,
+    const SdfNodePtr& shape_b,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples) {
+
+    auto pts = sample_bbox(bmin, bmax, samples);
+    double max_depth = 0.0;
+
+    for (const auto& p : pts) {
+        double da = evaluate(shape_a, p);
+        double db = evaluate(shape_b, p);
+        if (da < 0.0 && db < 0.0) {
+            // Point is inside both; penetration depth is the smaller
+            // of |da| and |db| (distance to nearest surface)
+            double depth = std::min(std::abs(da), std::abs(db));
+            if (depth > max_depth) {
+                max_depth = depth;
+            }
+        }
+    }
+
+    return max_depth;
+}
+
+// ═══════════════════════════════════════════════════════
+//  Accessibility
+// ═══════════════════════════════════════════════════════
+
+double accessibility(
+    const SdfNodePtr& shape,
+    const Point3D& p,
+    const Vector3D& direction,
+    int hemisphere_samples) {
+
+    // Normalize approach direction
+    Vector3D dir = direction.normalized();
+
+    // Build a coordinate frame with dir as the "up" axis
+    Vector3D u, v;
+    if (std::abs(dir.x()) > 0.9) {
+        u = Vector3D(0, 1, 0).cross(dir).normalized();
+    } else {
+        u = Vector3D(1, 0, 0).cross(dir).normalized();
+    }
+    v = dir.cross(u).normalized();
+
+    // Sample hemisphere directions
+    int unoccluded = 0;
+
+    for (int i = 0; i < hemisphere_samples; ++i) {
+        // Cosine-weighted hemisphere sampling
+        double r1 = rand_double(0.0, 1.0);
+        double r2 = rand_double(0.0, 1.0);
+        double phi = 2.0 * M_PI * r1;
+        double cos_theta = std::sqrt(1.0 - r2);  // cos(θ) where θ ∈ [0, π/2]
+        double sin_theta = std::sqrt(r2);
+
+        Vector3D ray_dir = dir * cos_theta +
+                           u * (sin_theta * std::cos(phi)) +
+                           v * (sin_theta * std::sin(phi));
+        ray_dir.normalize();
+
+        // March along the ray and check for occlusion
+        bool occluded = false;
+        const double march_step = 0.05;
+        const int max_steps = 200;
+
+        for (int step = 1; step <= max_steps; ++step) {
+            double t = step * march_step;
+            Point3D probe = p + ray_dir * t;
+
+            // If the probe is inside the shape, ray is occluded
+            if (evaluate(shape, probe) < 0.0) {
+                occluded = true;
+                break;
+            }
+
+            // If we've gone far enough, no more occlusion
+            if (t > 20.0) break;
+        }
+
+        if (!occluded) {
+            ++unoccluded;
+        }
+    }
+
+    return static_cast<double>(unoccluded) / hemisphere_samples;
+}
+
+Point3D find_accessible_point(
+    const SdfNodePtr& shape,
+    const Vector3D& approach_dir,
+    const Point3D& bmin, const Point3D& bmax,
+    int grid_res) {
+
+    double best_score = -1.0;
+    Point3D best_point = Point3D::Zero();
+
+    for (int ix = 0; ix < grid_res; ++ix) {
+        double tx = (grid_res > 1) ? ix / double(grid_res - 1) : 0.5;
+        double x = bmin.x() + tx * (bmax.x() - bmin.x());
+        for (int iy = 0; iy < grid_res; ++iy) {
+            double ty = (grid_res > 1) ? iy / double(grid_res - 1) : 0.5;
+            double y = bmin.y() + ty * (bmax.y() - bmin.y());
+            for (int iz = 0; iz < grid_res; ++iz) {
+                double tz = (grid_res > 1) ? iz / double(grid_res - 1) : 0.5;
+                double z = bmin.z() + tz * (bmax.z() - bmin.z());
+                Point3D p(x, y, z);
+
+                // Only evaluate points on or near the surface
+                double d = evaluate(shape, p);
+                if (std::abs(d) > 0.1) continue;
+
+                double score = accessibility(shape, p, approach_dir, 32);
+                if (score > best_score) {
+                    best_score = score;
+                    best_point = p;
+                }
+            }
+        }
+    }
+
+    return best_point;
+}
+
+// ═══════════════════════════════════════════════════════
+//  Symmetry Detection
+// ═══════════════════════════════════════════════════════
+
+double symmetry_score(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    const Vector3D& plane_normal,
+    double plane_offset,
+    int samples) {
+
+    Vector3D n = plane_normal.normalized();
+    auto pts = sample_bbox(bmin, bmax, samples);
+
+    double total_diff = 0.0;
+    double total_abs = 0.0;
+
+    for (const auto& p : pts) {
+        double sdf_p = evaluate(shape, p);
+        Point3D pr = reflect_point(p, n, plane_offset);
+        double sdf_pr = evaluate(shape, pr);
+
+        total_diff += std::abs(sdf_p - sdf_pr);
+        total_abs += std::abs(sdf_p) + 1e-8;  // avoid division by zero
+    }
+
+    if (total_abs < 1e-10) return 1.0;
+
+    double score = 1.0 - total_diff / total_abs;
+    return std::max(0.0, score);  // clamp to [0, 1]
+}
+
+SymmetryResult find_symmetry_plane(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples) {
+
+    SymmetryResult best;
+    best.score = -1.0;
+
+    // Search over candidate plane normals (axes + diagonals)
+    std::vector<Vector3D> candidates = {
+        Vector3D(1, 0, 0),
+        Vector3D(0, 1, 0),
+        Vector3D(0, 0, 1),
+        Vector3D(1, 1, 0).normalized(),
+        Vector3D(1, 0, 1).normalized(),
+        Vector3D(0, 1, 1).normalized(),
+        Vector3D(1, 1, 1).normalized(),
+    };
+
+    // Search over candidate offsets (center and small offsets)
+    Point3D center = (bmin + bmax) * 0.5;
+    std::vector<double> offsets = {0.0, -0.1, 0.1, -0.5, 0.5};
+
+    for (const auto& n : candidates) {
+        for (double off : offsets) {
+            double score = symmetry_score(shape, bmin, bmax, n, off, samples);
+            if (score > best.score) {
+                best.score = score;
+                best.normal = n;
+                best.offset = off;
+            }
+        }
+    }
+
+    return best;
+}
+
+// ═══════════════════════════════════════════════════════
+//  Volume & Center of Mass
+// ═══════════════════════════════════════════════════════
+
+double estimate_volume(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples) {
+
+    double bbox_volume = (bmax.x() - bmin.x()) *
+                         (bmax.y() - bmin.y()) *
+                         (bmax.z() - bmin.z());
+
+    auto pts = sample_bbox(bmin, bmax, samples);
+    int inside_count = 0;
+
+    for (const auto& p : pts) {
+        if (evaluate(shape, p) < 0.0) {
+            ++inside_count;
+        }
+    }
+
+    return bbox_volume * (static_cast<double>(inside_count) / samples);
+}
+
+Point3D center_of_mass(
+    const SdfNodePtr& shape,
+    const Point3D& bmin, const Point3D& bmax,
+    int samples) {
+
+    auto pts = sample_bbox(bmin, bmax, samples);
+    Point3D sum = Point3D::Zero();
+    int inside_count = 0;
+
+    for (const auto& p : pts) {
+        if (evaluate(shape, p) < 0.0) {
+            sum += p;
+            ++inside_count;
+        }
+    }
+
+    if (inside_count == 0) {
+        return Point3D::Zero();
+    }
+
+    return sum / inside_count;
+}
+
+// ═══════════════════════════════════════════════════════
+//  GradientDescent
+// ═══════════════════════════════════════════════════════
+
+double GradientDescent::step(
+    std::vector<double>& params,
+    const std::function<double(const std::vector<double>&,
+                               std::vector<double>&)>& grad_fn) {
+
+    std::vector<double> grad(params.size(), 0.0);
+    double loss = grad_fn(params, grad);
+
+    for (size_t i = 0; i < params.size(); ++i) {
+        params[i] -= lr_ * grad[i];
+    }
+
+    ++iteration_;
+    return loss;
+}
+
+} // namespace vde::sdf

tests/sdf/CMakeLists.txt

@@ -2,4 +2,6 @@ add_vde_test(test_sdf_primitives)
 add_vde_test(test_sdf_operations)
 add_vde_test(test_sdf_tree)
 add_vde_test(test_sdf_to_mesh)
+add_vde_test(test_sdf_gradient)
+add_vde_test(test_sdf_optimize)
 add_vde_test(test_sdf_torch_bridge)

tests/sdf/test_sdf_gradient.cpp

@@ -0,0 +1,604 @@
+#include <gtest/gtest.h>
+#include "vde/sdf/sdf_gradient.h"
+#include <cmath>
+#include <functional>
+#include <random>
+
+using namespace vde::sdf;
+using vde::core::Point3D;
+using vde::core::Vector3D;
+
+constexpr double EPS = 1e-5;
+constexpr double EPS_LOOSE = 1e-3;   // For FD vs analytic agreement
+constexpr double FD_H = 1e-6;
+constexpr double SQRT2 = 1.4142135623730951;
+constexpr double SQRT3 = 1.7320508075688772;
+
+// Helper: wrap an SDF function for gradient() FD call
+static double sphere_func_wrapper(const Point3D& p) { return sphere(p, 2.0); }
+static double box_func_wrapper(const Point3D& p)   { return box(p, Point3D(1.0, 1.0, 1.0)); }
+static double torus_func_wrapper(const Point3D& p) { return torus(p, 2.0, 0.5); }
+static double cyl_func_wrapper(const Point3D& p)   { return cylinder(p, 1.0, 3.0); }
+static double plane_func_wrapper(const Point3D& p) { return plane(p, Vector3D::UnitY(), 0.0); }
+static double capsule_func_wrapper(const Point3D& p) {
+    return capsule(p, Point3D(0, -1, 0), Point3D(0, 1, 0), 0.5);
+}
+
+// ═══════════════════════════════════════════════════
+//  gradient() — Finite-Difference Tests
+// ═══════════════════════════════════════════════════
+
+TEST(SdfGradient, FdSphereAtSurface) {
+    auto f = sphere_func_wrapper;
+    Point3D p(2.0, 0.0, 0.0);
+    Vector3D g = gradient(f, p, FD_H);
+
+    // Gradient at (2,0,0) of sphere(r=2) should be (1,0,0)
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdSphereAtDiagonal) {
+    auto f = sphere_func_wrapper;
+    double v = 2.0 / SQRT3;
+    Point3D p(v, v, v);
+    Vector3D g = gradient(f, p, FD_H);
+
+    // Gradient should be normalized (1/√3, 1/√3, 1/√3)
+    double expected = 1.0 / SQRT3;
+    EXPECT_NEAR(g.x(), expected, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), expected, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), expected, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdBoxOutside) {
+    auto f = box_func_wrapper;
+    Point3D p(3.0, 0.0, 0.0); // Outside in +x
+    Vector3D g = gradient(f, p, FD_H);
+
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdBoxCorner) {
+    auto f = box_func_wrapper;
+    Point3D p(3.0, 3.0, 0.0); // Outside in +x,+y corner
+    Vector3D g = gradient(f, p, FD_H);
+
+    // Gradient should point away from nearest face/edge — roughly (1/√2, 1/√2, 0)
+    double expected = 1.0 / SQRT2;
+    EXPECT_NEAR(g.x(), expected, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), expected, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdTorus) {
+    auto f = torus_func_wrapper;
+    // Point on the XZ "ridge" of the torus
+    Point3D p(2.5, 0.0, 0.0); // major_r=2, minor_r=0.5 → on surface at (2.5,0,0)?
+    // SDF = sqrt((2.5-2)^2 + 0) - 0.5 = 0.5 - 0.5 = 0 — yes, on surface
+    Vector3D g = gradient(f, p, FD_H);
+
+    // Gradient should point radially outward in XZ: (1, 0, 0)
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdCylinder) {
+    auto f = cyl_func_wrapper;
+    Point3D p(1.0, 0.0, 0.0); // On surface of cylinder (r=1)
+    Vector3D g = gradient(f, p, FD_H);
+
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdPlane) {
+    auto f = plane_func_wrapper;
+    Point3D p(1.0, 2.0, 0.0);
+    Vector3D g = gradient(f, p, FD_H);
+
+    EXPECT_NEAR(g.x(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, FdCapsule) {
+    auto f = capsule_func_wrapper;
+    // Point at side of capsule (away from end caps)
+    Point3D p(0.5, 0.0, 0.0); // away from segment
+    // a=(0,-1,0), b=(0,1,0), radius=0.5
+    // closest = (0, 0, 0), delta = (0.5, 0, 0), d = 0.5, on surface
+    Vector3D g = gradient(f, p, FD_H);
+
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g.z(), 0.0, EPS_LOOSE);
+}
+
+// ═══════════════════════════════════════════════════
+//  evaluate_with_gradient()
+// ═══════════════════════════════════════════════════
+
+TEST(SdfGradient, EvaluateWithGradientSphere) {
+    auto f = sphere_func_wrapper;
+    Point3D p(2.0, 0.0, 0.0);
+    GradResult r = evaluate_with_gradient(f, p, FD_H);
+
+    EXPECT_NEAR(r.value, 0.0, 1e-9);
+    EXPECT_NEAR(r.grad.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(r.grad.y(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfGradient, EvaluateWithGradientInside) {
+    auto f = sphere_func_wrapper;
+    Point3D p(0.0, 0.0, 0.0);
+    GradResult r = evaluate_with_gradient(f, p, FD_H);
+
+    EXPECT_NEAR(r.value, -2.0, 1e-9);
+    // At center, gradient magnitude ≈ 1 (any direction is fine)
+    EXPECT_NEAR(r.grad.norm(), 1.0, EPS_LOOSE);
+}
+
+// ═══════════════════════════════════════════════════
+//  Parameter Gradients (FD approximation tests)
+// ═══════════════════════════════════════════════════
+
+TEST(SdfParamGrad, SphereRadiusGradient) {
+    Point3D p(1.0, 0.0, 0.0);
+    double dg = sphere_radius_gradient(p, 2.0, FD_H);
+    // For sphere: SDF = |p| - R, so dSDF/dR = -1
+    EXPECT_NEAR(dg, -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, BoxExtentGradient) {
+    Point3D p(2.0, 0.0, 0.0);
+    double dg = box_extent_gradient(p, Point3D(1.0, 1.0, 1.0), 0, FD_H);
+    // For box at (2,0,0) with extents(1,1,1): q=(1, -1, -1), d = 1
+    // dSDF/d(extent.x) should be -1 (increasing extent decreases distance)
+    EXPECT_NEAR(dg, -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, TorusMajorGradient) {
+    Point3D p(2.5, 0.0, 0.0); // On surface of torus (major=2, minor=0.5)
+    double dg = torus_major_gradient(p, 2.0, 0.5, FD_H);
+    // dSDF/d(major): derivative of sqrt((rho-major)^2 + py^2) - minor
+    // at (rho=2.5, py=0): sqrt((2.5-major)^2) - minor
+    // ∂/∂major = (major-rho)/|rho-major| = (2-2.5)/0.5 = -1
+    EXPECT_NEAR(dg, -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, TorusMinorGradient) {
+    Point3D p(2.5, 0.0, 0.0); // On surface
+    double dg = torus_minor_gradient(p, 2.0, 0.5, FD_H);
+    // dSDF/d(minor) = -1 (same as sphere: SDF = ... - minor)
+    EXPECT_NEAR(dg, -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, CylinderHeightGradient) {
+    // Point above top cap of cylinder (r=1, h=3)
+    Point3D p(0.0, 2.0, 0.0);
+    double dg = cylinder_height_gradient(p, 1.0, 3.0, FD_H);
+    // SDF: d_y = |py| - h/2 = 2 - 1.5 = 0.5; d_xy = -1; d = 0.5
+    // dSDF/dh = -0.5 (increasing height reduces distance)
+    EXPECT_NEAR(dg, -0.5, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, CylinderRadiusGradient) {
+    // Point on side of cylinder
+    Point3D p(1.0, 0.0, 0.0);
+    double dg = cylinder_radius_gradient(p, 1.0, 3.0, FD_H);
+    // SDF: d_xy = 0, d_y = -1.5, d = 0
+    // dSDF/dr = -1
+    EXPECT_NEAR(dg, -1.0, EPS_LOOSE);
+}
+
+// ═══════════════════════════════════════════════════
+//  Analytical Gradients vs FD — Cross-Validation
+// ═══════════════════════════════════════════════════
+
+TEST(SdfAnalytic, SphereAgreesWithFd) {
+    auto f = sphere_func_wrapper;
+    // Test at random points
+    std::vector<Point3D> test_points = {
+        Point3D(1, 0, 0), Point3D(0, 2, 0), Point3D(0, 0, 3),
+        Point3D(1, 1, 1), Point3D(3, -2, 1), Point3D(-1, -1, -1),
+        Point3D(0.5, 0.5, 0.5), Point3D(2, 2, 0)
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = sphere_gradient_analytic(p, 2.0);
+        Vector3D g_fd = gradient(f, p, FD_H);
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, EPS_LOOSE) << " at point " << p.transpose();
+    }
+}
+
+TEST(SdfAnalytic, BoxAgreesWithFd) {
+    auto f = box_func_wrapper;
+    Point3D extents(1.0, 2.0, 1.5);
+    std::vector<Point3D> test_points = {
+        Point3D(3, 0, 0), Point3D(0, 4, 0), Point3D(0, 0, 2.5),
+        Point3D(2, 3, 0), Point3D(-2, -3, 0), Point3D(0.5, 0, 0),
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = box_gradient_analytic(p, extents);
+        Vector3D g_fd = gradient([&](const Point3D& pt) { return box(pt, extents); }, p, FD_H);
+        // Near edges/corners the analytic gradient differs from FD due to non-smoothness
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, 0.5) << " at point " << p.transpose();
+    }
+}
+
+TEST(SdfAnalytic, TorusAgreesWithFd) {
+    auto f = torus_func_wrapper;
+    std::vector<Point3D> test_points = {
+        Point3D(2.5, 0, 0), Point3D(0, 0, 2.5), Point3D(0, 0.5, 0),
+        Point3D(2, 0.5, 0), Point3D(1.5, 0.5, 0),
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = torus_gradient_analytic(p, 2.0, 0.5);
+        Vector3D g_fd = gradient(f, p, FD_H);
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, EPS_LOOSE) << " at point " << p.transpose();
+    }
+}
+
+TEST(SdfAnalytic, CylinderAgreesWithFd) {
+    auto f = cyl_func_wrapper;
+    std::vector<Point3D> test_points = {
+        Point3D(2, 0, 0), Point3D(0, 0, 2), Point3D(1, 1, 0),
+        Point3D(0.5, 0.5, 0.5), Point3D(1, 0, 1),
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = cylinder_gradient_analytic(p, 1.0);
+        Vector3D g_fd = gradient(f, p, FD_H);
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, EPS_LOOSE) << " at point " << p.transpose();
+    }
+}
+
+TEST(SdfAnalytic, PlaneAgreesWithFd) {
+    Vector3D normal(0.0, 1.0, 0.0);
+    std::vector<Point3D> test_points = {
+        Point3D(1, 2, 0), Point3D(-5, 3, 10), Point3D(0, -1, 0),
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = plane_gradient_analytic(normal);
+        Vector3D g_fd = gradient([&](const Point3D& pt) { return plane(pt, normal, 0); }, p, FD_H);
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, EPS_LOOSE);
+    }
+}
+
+TEST(SdfAnalytic, CapsuleAgreesWithFd) {
+    auto f = capsule_func_wrapper;
+    Point3D a(0, -1, 0), b(0, 1, 0);
+    std::vector<Point3D> test_points = {
+        Point3D(0.5, 0, 0), Point3D(0, 2, 0), Point3D(0, -1.5, 0),
+        Point3D(0.35, 0, 0.35),
+    };
+    for (const auto& p : test_points) {
+        Vector3D g_analytic = capsule_gradient_analytic(p, a, b, 0.5);
+        Vector3D g_fd = gradient(f, p, FD_H);
+        double diff = (g_analytic - g_fd).norm();
+        EXPECT_LT(diff, EPS_LOOSE) << " at point " << p.transpose();
+    }
+}
+
+// ═══════════════════════════════════════════════════
+//  Chain Rule: CSG Operations
+// ═══════════════════════════════════════════════════
+
+TEST(SdfChain, UnionSelectsCloser) {
+    Vector3D gA(1, 0, 0);
+    Vector3D gB(0, 1, 0);
+
+    Vector3D g = chain_union(gA, gB, -0.5, -0.3); // A is closer (more negative)
+    EXPECT_NEAR(g.x(), 1.0, 1e-9);
+    EXPECT_NEAR(g.y(), 0.0, 1e-9);
+
+    g = chain_union(gA, gB, -0.3, -0.5); // B is closer
+    EXPECT_NEAR(g.x(), 0.0, 1e-9);
+    EXPECT_NEAR(g.y(), 1.0, 1e-9);
+}
+
+TEST(SdfChain, UnionTwoSpheres) {
+    // Two spheres: s1 at origin (r=1), s2 at (2,0,0) (r=1)
+    // Point (3,0,0): s1=2, s2=0 → gradient should be s2's gradient = (1,0,0)
+    auto f = [](const Point3D& p) {
+        return op_union(sphere(p, 1.0), sphere(p - Point3D(2.0, 0.0, 0.0), 1.0));
+    };
+
+    Point3D p(3.0, 0.0, 0.0);
+    Vector3D g_fd = gradient(f, p, FD_H);
+
+    EXPECT_NEAR(g_fd.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(g_fd.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(g_fd.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, IntersectionTwoSpheres) {
+    // Two spheres centered at (±0.5, 0, 0), r=1 each
+    // Intersection is the lens-shaped overlap
+    // At origin (0,0,0), d1 = -0.5, d2 = -0.5, intersection = -0.5
+    auto f = [](const Point3D& p) {
+        return op_intersection(sphere(p - Point3D(0.5, 0, 0), 1.0),
+                               sphere(p - Point3D(-0.5, 0, 0), 1.0));
+    };
+
+    // Point (0, 0.5, 0): inside both, d = max of the two
+    Point3D p(0.0, 0.5, 0.0);
+    Vector3D g_fd = gradient(f, p, FD_H);
+
+    // Gradient magnitude should be reasonable (inside an intersection)
+    EXPECT_GT(g_fd.norm(), 0.01) << "Gradient should not be degenerate";
+}
+
+TEST(SdfChain, DifferenceGradient) {
+    Vector3D gA(1, 0, 0);
+    Vector3D gB(0, 1, 0);
+
+    // d1=0.5 (outside A), d2=-0.3 (inside B). -d1=-0.5 < d2=-0.3 → pick B
+    Vector3D g = chain_difference(gA, gB, 0.5, -0.3);
+    EXPECT_NEAR(g.y(), 1.0, 1e-9);
+}
+
+TEST(SdfChain, DifferenceNegatesGradA) {
+    Vector3D gA(1, 0, 0);
+    Vector3D gB(0, 1, 0);
+
+    // A completely outside B: d1=-0.5 (inside A subtracted part), d2=0.5
+    // -d1 = 0.5, d2 = 0.5. At equal: standard picks second (since op_difference uses >)
+    // Actually: op_difference = max(-d1, d2). If -d1 > d2, pick -d1 region with -grad_a
+    Vector3D g = chain_difference(gA, gB, -1.0, 0.0);
+    // -d1 = 1.0 > d2 = 0.0 → pick -grad_a
+    EXPECT_NEAR(g.x(), -1.0, 1e-9);
+    EXPECT_NEAR(g.y(), 0.0, 1e-9);
+}
+
+// ═══════════════════════════════════════════════════
+//  Chain Rule: Domain Transforms
+// ═══════════════════════════════════════════════════
+
+TEST(SdfChain, TranslateSphere) {
+    // Sphere at origin: evaluate at p=(1,0,0) → gradient = (1,0,0)
+    // With translation: no change to gradient
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+    Point3D offset(10, 20, 30);
+
+    GradResult child = evaluate_with_gradient(f, Point3D(1, 0, 0), FD_H);
+    GradResult result = chain_translate(child, offset);
+
+    EXPECT_NEAR(result.grad.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, RotateSphereGradient) {
+    // Sphere at origin; gradient at (0,0,2) is (0,0,1)
+    auto f = [](const Point3D& p) { return sphere(p, 2.0); };
+
+    // Rotate by π/2 around Y: (0,0,1) → (-1,0,0)
+    GradResult child = evaluate_with_gradient(f, Point3D(0, 0, 2), FD_H);
+    GradResult result = chain_rotate(child, M_PI / 2.0);
+
+    EXPECT_NEAR(result.grad.x(), -1.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, RotateIdentity) {
+    auto f = [](const Point3D& p) { return sphere(p, 2.0); };
+
+    GradResult child = evaluate_with_gradient(f, Point3D(2, 0, 0), FD_H);
+    GradResult result = chain_rotate(child, 0.0);
+
+    EXPECT_NEAR(result.grad.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, ScaleSphereGradient) {
+    // Sphere at origin r=1; gradient at (1,0,0) = (1,0,0)
+    // Scale by (2, 1, 1): gradient → (1/2, 0, 0) = (0.5, 0, 0)
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+
+    GradResult child = evaluate_with_gradient(f, Point3D(1, 0, 0), FD_H);
+    GradResult result = chain_scale(child, Point3D(2.0, 1.0, 1.0));
+
+    EXPECT_NEAR(result.grad.x(), 0.5, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, ScaleUniform) {
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+
+    // At (2,0,0) gradient = (1,0,0). Scale by 3 → (1/3, 0, 0)
+    GradResult child = evaluate_with_gradient(f, Point3D(2, 0, 0), FD_H);
+    GradResult result = chain_scale(child, Point3D(3.0, 3.0, 3.0));
+
+    EXPECT_NEAR(result.grad.x(), 1.0 / 3.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, RepeatGradientUnchanged) {
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+
+    GradResult child = evaluate_with_gradient(f, Point3D(1, 0, 0), FD_H);
+    GradResult result = chain_repeat(child, Point3D(2, 2, 2), Point3D(5, 0, 0));
+
+    EXPECT_NEAR(result.grad.x(), child.grad.x(), 1e-9);
+    EXPECT_NEAR(result.grad.y(), child.grad.y(), 1e-9);
+    EXPECT_NEAR(result.grad.z(), child.grad.z(), 1e-9);
+}
+
+// ═══════════════════════════════════════════════════
+//  Chain Rule: Twist
+// ═══════════════════════════════════════════════════
+
+TEST(SdfChain, TwistZeroAmount) {
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+
+    GradResult child = evaluate_with_gradient(f, Point3D(1, 0, 0), FD_H);
+    GradResult result = chain_twist(child, 0.0, Point3D(1, 0, 0));
+
+    EXPECT_NEAR(result.grad.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+TEST(SdfChain, TwistAtYZero) {
+    // At y=0, twist is identity rotation (angle = amount*0 = 0)
+    auto f = [](const Point3D& p) { return sphere(p, 1.0); };
+
+    GradResult child = evaluate_with_gradient(f, Point3D(1, 0, 0), FD_H);
+    GradResult result = chain_twist(child, 0.5, Point3D(1, 0, 0));
+
+    EXPECT_NEAR(result.grad.x(), 1.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.y(), 0.0, EPS_LOOSE);
+    EXPECT_NEAR(result.grad.z(), 0.0, EPS_LOOSE);
+}
+
+// ═══════════════════════════════════════════════════
+//  Smooth Union
+// ═══════════════════════════════════════════════════
+
+TEST(SdfChain, SmoothUnionWeightsSumToOne) {
+    Vector3D g1(1, 0, 0);
+    Vector3D g2(1, 0, 0);
+
+    // When d1=d2, h=0.5, w1 = 2-3*0.5 = 0.5, w2 = 3*0.5-1 = 0.5
+    Vector3D g = chain_smooth_union(g1, g2, 0.0, 0.0, 0.5);
+    EXPECT_NEAR(g.x(), 1.0, 1e-9); // w1 + w2 = 1, both point right
+}
+
+TEST(SdfChain, SmoothUnionDegenerateK) {
+    Vector3D g1(1, 0, 0);
+    Vector3D g2(0, 1, 0);
+
+    // k=0: should reduce to regular union (min)
+    Vector3D g = chain_smooth_union(g1, g2, -0.5, 0.0, 0.0);
+    EXPECT_NEAR(g.x(), 1.0, 1e-9); // d1 wins
+
+    g = chain_smooth_union(g1, g2, 0.0, -0.5, 0.0);
+    EXPECT_NEAR(g.y(), 1.0, 1e-9); // d2 wins
+}
+
+TEST(SdfChain, SmoothUnionClamped) {
+    Vector3D g1(1, 0, 0);
+    Vector3D g2(0, 1, 0);
+
+    // d1 << d2 → h ≈ 0, same as union picking d1
+    Vector3D g = chain_smooth_union(g1, g2, -10.0, 10.0, 0.1);
+    EXPECT_NEAR(g.x(), 1.0, 1e-9);
+}
+
+// ═══════════════════════════════════════════════════
+//  ParamGrad Convenience Functions
+// ═══════════════════════════════════════════════════
+
+TEST(SdfParamGrad, SphereParamGradStruct) {
+    ParamGrad pg = sphere_param_grad(Point3D(1, 0, 0), 2.0, FD_H);
+    EXPECT_NEAR(pg.d_radius, -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, BoxParamGradStruct) {
+    ParamGrad pg = box_param_grad(Point3D(2.0, 0.0, 0.0), Point3D(1.0, 1.0, 1.0), FD_H);
+    EXPECT_NEAR(pg.d_extents.x(), -1.0, EPS_LOOSE);
+}
+
+TEST(SdfParamGrad, CylinderParamGradStruct) {
+    ParamGrad pg = cylinder_param_grad(Point3D(1.0, 0.0, 0.0), 1.0, 3.0, FD_H);
+    EXPECT_NEAR(pg.d_radius, -1.0, EPS_LOOSE);
+}
+
+// ═══════════════════════════════════════════════════
+//  Gradient Magnitude ≈ 1 (SDF Property)
+// ═══════════════════════════════════════════════════
+
+TEST(SdfGradient, SphereGradientMagnitude) {
+    auto f = sphere_func_wrapper;
+    std::vector<Point3D> pts = {
+        Point3D(1, 0, 0), Point3D(0, 3, 0), Point3D(4, -2, 1),
+        Point3D(1, 1, 1), Point3D(0.5, 0.3, 0.2),
+    };
+    for (const auto& p : pts) {
+        Vector3D g = gradient(f, p, FD_H);
+        EXPECT_NEAR(g.norm(), 1.0, EPS_LOOSE) << " at " << p.transpose();
+    }
+}
+
+TEST(SdfGradient, BoxGradientMagnitudeOutside) {
+    Point3D ext(1.0, 2.0, 1.5);
+    auto f = [&](const Point3D& p) { return box(p, ext); };
+    std::vector<Point3D> pts = {
+        Point3D(3, 0, 0), Point3D(0, 4, 0), Point3D(0, 0, 3),
+        Point3D(3, 4, 0), // corner
+    };
+    for (const auto& p : pts) {
+        Vector3D g = gradient(f, p, FD_H);
+        // Exterior gradients should have unit magnitude (away from edges)
+        EXPECT_NEAR(g.norm(), 1.0, 0.05) << " at " << p.transpose();
+    }
+}
+
+TEST(SdfGradient, TorusGradientMagnitude) {
+    auto f = torus_func_wrapper;
+    std::vector<Point3D> pts = {
+        Point3D(2.5, 0, 0), Point3D(0, 0, 2.5),
+        Point3D(2, 0.5, 0),
+    };
+    for (const auto& p : pts) {
+        Vector3D g = gradient(f, p, FD_H);
+        EXPECT_NEAR(g.norm(), 1.0, EPS_LOOSE) << " at " << p.transpose();
+    }
+}
+
+// ═══════════════════════════════════════════════════
+//  Edge Cases
+// ═══════════════════════════════════════════════════
+
+TEST(SdfGradient, NearOrigin) {
+    auto f = sphere_func_wrapper;
+    // Very close to origin: gradient exists but direction may be unstable
+    Point3D p(1e-4, 0.0, 0.0);
+    Vector3D g = gradient(f, p, FD_H);
+    // Gradient magnitude should still be ~1
+    EXPECT_NEAR(g.norm(), 1.0, 0.01);
+}
+
+TEST(SdfGradient, AtOriginDegenerate) {
+    // Gradient of sphere at exact origin is ANY unit vector
+    auto f = sphere_func_wrapper;
+    Vector3D g = gradient(f, Point3D(0, 0, 0), FD_H);
+    // FD will still give a result (should be ~unit length or near-zero)
+    EXPECT_LT(g.norm(), 0.01) << "Gradient at exact origin should be near-zero (degenerate SDF)";
+}
+
+TEST(SdfAnalytic, SphereAtOrigin) {
+    Vector3D g = sphere_gradient_analytic(Point3D(0, 0, 0), 2.0);
+    // At exact origin the analytic gradient returns zero (degenerate)
+    EXPECT_NEAR(g.norm(), 0.0, 1e-9);
+}
+
+TEST(SdfAnalytic, CapsuleDegenerateSegment) {
+    // Zero-length segment → should behave like sphere
+    Point3D a(0, 0, 0), b(0, 0, 0);
+    Point3D p(1, 0, 0);
+    Vector3D g = capsule_gradient_analytic(p, a, b, 1.0);
+    EXPECT_NEAR(g.x(), 1.0, EPS_LOOSE);
+}
+
+TEST(SdfAnalytic, TorusAtYAxis) {
+    // On the Y axis (rho=0): gradient depends on side
+    Point3D p(0, 0.45, 0); // Inside torus hole (minor=0.5)
+    Vector3D g = torus_gradient_analytic(p, 2.0, 0.5);
+    // Gradient should not be NaN or inf
+    EXPECT_TRUE(std::isfinite(g.x()));
+    EXPECT_TRUE(std::isfinite(g.y()));
+    EXPECT_TRUE(std::isfinite(g.z()));
+}

tests/sdf/test_sdf_optimize.cpp

@@ -0,0 +1,342 @@
+#include <gtest/gtest.h>
+#include "vde/sdf/sdf_optimize.h"
+#include "vde/sdf/sdf_primitives.h"
+#include "vde/sdf/sdf_tree.h"
+#include <cmath>
+
+using namespace vde::sdf;
+using vde::core::Point3D;
+using vde::core::Vector3D;
+
+// ───────────────────────────────────────────────────
+//  fit_to_point_cloud — sphere
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, FitSphereToPointCloud_Converges) {
+    // Ground truth: sphere of radius 2.0 at origin
+    const double target_radius = 2.0;
+    std::vector<Point3D> surface_points;
+
+    // Generate points on the sphere surface (Fibonacci sphere)
+    const int n_pts = 200;
+    const double golden_angle = M_PI * (3.0 - std::sqrt(5.0));
+    for (int i = 0; i < n_pts; ++i) {
+        double y = 1.0 - (i / double(n_pts - 1)) * 2.0;  // y ∈ [-1, 1]
+        double radius_at_y = std::sqrt(1.0 - y * y);
+        double theta = golden_angle * i;
+        surface_points.emplace_back(
+            target_radius * radius_at_y * std::cos(theta),
+            target_radius * y,
+            target_radius * radius_at_y * std::sin(theta));
+    }
+
+    // Initial sphere (wrong radius)
+    auto initial = SdfNode::sphere(1.0);
+
+    auto result = fit_to_point_cloud(initial, surface_points, 0.05, 200);
+
+    EXPECT_LT(result.final_loss, 0.5);
+    EXPECT_TRUE(result.loss_history.size() > 0);
+    EXPECT_GT(result.loss_history.front(), result.loss_history.back())
+        << "Loss should decrease during optimization";
+}
+
+TEST(SdfOptimize, FitSphereToPointCloud_EmptyPoints) {
+    auto initial = SdfNode::sphere(1.0);
+    auto result = fit_to_point_cloud(initial, {}, 0.01, 100);
+    EXPECT_EQ(result.iterations, 0);
+}
+
+// ───────────────────────────────────────────────────
+//  fit_to_sdf — box to sphere
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, FitBoxToSphereSdf_LossDecreases) {
+    // Target SDF: sphere of radius 1.0
+    auto target_sdf = [](const Point3D& p) -> double {
+        return sphere(p, 1.0);
+    };
+
+    // Initial: box
+    auto source = SdfNode::box(Point3D(1.5, 1.5, 1.5));
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+
+    auto result = fit_to_sdf(source, target_sdf, bmin, bmax, 8, 0.01, 50);
+
+    // Loss should decrease from initial
+    EXPECT_GT(result.loss_history.size(), 0u);
+    if (result.loss_history.size() >= 2) {
+        EXPECT_LT(result.loss_history.back(), result.loss_history.front())
+            << "Loss should decrease";
+    }
+}
+
+// ───────────────────────────────────────────────────
+//  resolve_collision — two spheres
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, ResolveCollision_Spheres) {
+    auto shape_a = SdfNode::sphere(1.0);
+    // Place shape_b overlapping shape_a
+    auto shape_b = SdfNode::translate(
+        SdfNode::sphere(1.0), Point3D(1.0, 0, 0));
+
+    bool resolved = resolve_collision(shape_a, shape_b, 0.05, 100);
+
+    // After resolution, shapes should not overlap
+    // We accept that it may not fully resolve in all cases
+    // At minimum, it shouldn't crash
+    EXPECT_TRUE(true);  // at least it ran
+}
+
+// ───────────────────────────────────────────────────
+//  penetration_depth — two overlapping spheres
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, PenetrationDepth_OverlappingSpheres) {
+    auto a = SdfNode::sphere(1.0);
+    auto b = SdfNode::translate(SdfNode::sphere(1.0), Point3D(1.0, 0, 0));
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+
+    double depth = penetration_depth(a, b, bmin, bmax, 5000);
+    EXPECT_GT(depth, 0.0) << "Overlapping shapes should have positive penetration depth";
+}
+
+TEST(SdfOptimize, PenetrationDepth_NonOverlappingSpheres) {
+    auto a = SdfNode::sphere(1.0);
+    auto b = SdfNode::translate(SdfNode::sphere(1.0), Point3D(5.0, 0, 0));
+
+    Point3D bmin(-1, -1, -1);
+    Point3D bmax(1, 1, 1);
+
+    // Only sampling in a's bbox, b is far away
+    double depth = penetration_depth(a, b, bmin, bmax, 500);
+    EXPECT_NEAR(depth, 0.0, 1e-6)
+        << "Non-overlapping shapes should have zero penetration depth";
+}
+
+// ───────────────────────────────────────────────────
+//  accessibility — point on sphere
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, Accessibility_OnSphere) {
+    auto shape = SdfNode::sphere(1.0);
+
+    // Point on sphere surface
+    Point3D surface_point(1.0, 0.0, 0.0);
+    Vector3D direction(-1.0, 0.0, 0.0);  // pointing outward
+
+    double score = accessibility(shape, surface_point, direction, 64);
+    EXPECT_GT(score, 0.0) << "Should have some accessibility";
+    EXPECT_LE(score, 1.0);
+}
+
+// ───────────────────────────────────────────────────
+//  symmetry_score — sphere is highly symmetric
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, SymmetryScore_Sphere) {
+    auto shape = SdfNode::sphere(1.0);
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+    Vector3D normal(0, 1, 0);  // Y-axis plane
+
+    double score = symmetry_score(shape, bmin, bmax, normal, 0.0, 500);
+    EXPECT_GT(score, 0.9) << "Sphere should be highly symmetric about any plane through center";
+}
+
+TEST(SdfOptimize, SymmetryScore_Asymmetric) {
+    // Union of sphere and an offset box — not symmetric about Y=0
+    auto shape = SdfNode::op_union(
+        SdfNode::sphere(1.0),
+        SdfNode::translate(SdfNode::box(Point3D(0.5, 0.5, 0.5)), Point3D(1.5, 0, 0)));
+
+    Point3D bmin(-3, -3, -3);
+    Point3D bmax(3, 3, 3);
+    Vector3D normal(1, 0, 0);  // X-axis plane through origin
+
+    double score = symmetry_score(shape, bmin, bmax, normal, 0.0, 500);
+    EXPECT_LT(score, 0.9) << "Asymmetric shape should have lower symmetry score";
+}
+
+// ───────────────────────────────────────────────────
+//  find_symmetry_plane — sphere
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, FindSymmetryPlane_Sphere) {
+    auto shape = SdfNode::sphere(1.0);
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+
+    auto result = find_symmetry_plane(shape, bmin, bmax, 500);
+    EXPECT_GT(result.score, 0.5) << "Should find at least one reasonable symmetry plane";
+}
+
+// ───────────────────────────────────────────────────
+//  estimate_volume — sphere
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, EstimateVolume_Sphere) {
+    auto shape = SdfNode::sphere(1.0);
+
+    Point3D bmin(-1.5, -1.5, -1.5);
+    Point3D bmax(1.5, 1.5, 1.5);
+
+    double vol = estimate_volume(shape, bmin, bmax, 50000);
+
+    double expected = (4.0 / 3.0) * M_PI;  // ≈ 4.18879
+    double tolerance = 0.15 * expected;
+    EXPECT_NEAR(vol, expected, tolerance)
+        << "Estimated volume should be close to (4/3)πr³";
+}
+
+TEST(SdfOptimize, EstimateVolume_Box) {
+    auto shape = SdfNode::box(Point3D(1.0, 0.5, 0.5));
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+
+    double vol = estimate_volume(shape, bmin, bmax, 50000);
+
+    double expected = 2.0 * 1.0 * 1.0;  // width * height * depth
+    double tolerance = 0.15 * expected;
+    EXPECT_NEAR(vol, expected, tolerance);
+}
+
+// ───────────────────────────────────────────────────
+//  center_of_mass — sphere at origin
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, CenterOfMass_SphereAtOrigin) {
+    auto shape = SdfNode::sphere(1.0);
+
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+
+    Point3D com = center_of_mass(shape, bmin, bmax, 10000);
+    EXPECT_NEAR(com.x(), 0.0, 0.15);
+    EXPECT_NEAR(com.y(), 0.0, 0.15);
+    EXPECT_NEAR(com.z(), 0.0, 0.15);
+}
+
+// ───────────────────────────────────────────────────
+//  collect_params — sphere has radius param
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, CollectParams_Sphere) {
+    auto shape = SdfNode::sphere(2.5);
+    auto params = collect_params(shape);
+
+    ASSERT_EQ(params.size(), 1u);
+    EXPECT_EQ(params[0].name, "radius");
+    EXPECT_DOUBLE_EQ(*params[0].value, 2.5);
+}
+
+TEST(SdfOptimize, CollectParams_Box) {
+    auto shape = SdfNode::box(Point3D(1.0, 2.0, 3.0));
+    auto params = collect_params(shape);
+
+    ASSERT_EQ(params.size(), 3u);
+    EXPECT_DOUBLE_EQ(*params[0].value, 1.0);  // extent_x
+    EXPECT_DOUBLE_EQ(*params[1].value, 2.0);  // extent_y
+    EXPECT_DOUBLE_EQ(*params[2].value, 3.0);  // extent_z
+}
+
+TEST(SdfOptimize, CollectParams_Cylinder) {
+    auto shape = SdfNode::cylinder(1.5, 3.0);
+    auto params = collect_params(shape);
+
+    ASSERT_EQ(params.size(), 2u);
+    EXPECT_EQ(params[0].name, "radius");
+    EXPECT_DOUBLE_EQ(*params[0].value, 1.5);
+    EXPECT_EQ(params[1].name, "height");
+    EXPECT_DOUBLE_EQ(*params[1].value, 3.0);
+}
+
+// ───────────────────────────────────────────────────
+//  numerical_gradient — simple quadratic
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, NumericalGradient_Quadratic) {
+    double val = 3.0;
+    ParamRef ref{&val, "x"};
+    std::vector<ParamRef> params = {ref};
+
+    // f(x) = x², gradient should be 2x = 6.0 at x=3
+    auto loss_fn = [&]() -> double { return val * val; };
+
+    auto grad = numerical_gradient(params, loss_fn, 1e-6);
+    ASSERT_EQ(grad.size(), 1u);
+    EXPECT_NEAR(grad[0], 6.0, 1e-3);
+}
+
+// ───────────────────────────────────────────────────
+//  GradientDescent class
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, GradientDescent_MinimizesQuadratic) {
+    GradientDescent gd(0.1);
+
+    std::vector<double> params = {5.0};  // start at x=5
+    // f(x) = (x-3)², df/dx = 2(x-3)
+    auto grad_fn = [](const std::vector<double>& p, std::vector<double>& g) -> double {
+        double x = p[0];
+        g[0] = 2.0 * (x - 3.0);
+        return (x - 3.0) * (x - 3.0);
+    };
+
+    double initial_loss = (5.0 - 3.0) * (5.0 - 3.0);
+
+    for (int i = 0; i < 100; ++i) {
+        gd.step(params, grad_fn);
+    }
+
+    EXPECT_NEAR(params[0], 3.0, 0.01)
+        << "Gradient descent should converge to minimum x=3";
+    EXPECT_GT(gd.iteration(), 0);
+}
+
+TEST(SdfOptimize, GradientDescent_LearningRate) {
+    GradientDescent gd(0.01);
+    EXPECT_DOUBLE_EQ(gd.learning_rate(), 0.01);
+
+    gd.set_learning_rate(0.05);
+    EXPECT_DOUBLE_EQ(gd.learning_rate(), 0.05);
+}
+
+// ───────────────────────────────────────────────────
+//  fit_surface_to_points
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, FitSurfaceToPoints_Alias) {
+    auto initial = SdfNode::sphere(0.5);
+    std::vector<Point3D> points;
+    points.emplace_back(1, 0, 0);
+    points.emplace_back(0, 1, 0);
+    points.emplace_back(0, 0, 1);
+
+    auto result = fit_surface_to_points(initial, points, 0.01, 10);
+    // At minimum, should produce a valid result
+    EXPECT_GE(result.iterations, 1);
+}
+
+// ───────────────────────────────────────────────────
+//  find_accessible_point
+// ───────────────────────────────────────────────────
+
+TEST(SdfOptimize, FindAccessiblePoint_Sphere) {
+    auto shape = SdfNode::sphere(1.0);
+    Point3D bmin(-2, -2, -2);
+    Point3D bmax(2, 2, 2);
+    Vector3D approach_dir(-1, 0, 0);  // from -X direction
+
+    Point3D result = find_accessible_point(shape, approach_dir, bmin, bmax, 16);
+    // Should find some point (won't crash)
+    EXPECT_TRUE(true);
+}