feat(sdf): S11 — SDF 隐式建模完整模块

S11-A: SDF 图元库
- 14 种基础形状（inline header）
- sphere/box/round_box/torus/capsule/cylinder/cone/plane/ellipsoid
- triangular_prism/hex_prism/link/wedge
- 2D 挤出（extrusion/extrusion_bounded/revolution）

S11-B: SDF 操作 + 域变形
- 锐利布尔: union/intersection/difference
- 平滑布尔: smooth_union/smooth_intersection/smooth_difference
- 修饰器: round/onion
- 域变形: repeat/mirror/rotate/translate/scale/twist/bend/elongate/displace/cheap_bend
- 24 项测试

S11-C: CSG 表达式树 + 网格转换 + Python 绑定
- SdfNode 树结构 — 工厂构造函数，递归求值
- sdf_to_mesh — 基于 marching_cubes 的 SDF→网格
- bind_sdf — pybind11 Python 绑定

文件: 14 文件，2,545 行（测试 1,296 行）

include/vde/sdf/sdf_operations.h

@@ -1,165 +1,138 @@
 #pragma once
+#include "vde/core/point.h"
 #include <cmath>
 #include <algorithm>
-#include "vde/core/point.h"
+#include <numbers>
 
 namespace vde::sdf {
 
-// ═══════════════════════════════════════════════════
-//  Boolean Operations (sharp)
-// ═══════════════════════════════════════════════════
+using core::Point3D;
+using core::Vector3D;
+
+// ── Boolean operations on distance values ──
 
-/// Union (min) — combine two shapes
 [[nodiscard]] inline double op_union(double d1, double d2) {
     return std::min(d1, d2);
 }
 
-/// Intersection (max) — region common to both shapes
 [[nodiscard]] inline double op_intersection(double d1, double d2) {
     return std::max(d1, d2);
 }
 
-/// Difference — subtract d2 from d1 (d1 \ d2)
 [[nodiscard]] inline double op_difference(double d1, double d2) {
     return std::max(-d1, d2);
 }
 
-// ═══════════════════════════════════════════════════
-//  Smooth Boolean Operations (Inigo Quilez)
-// ═══════════════════════════════════════════════════
-
-namespace detail {
-    template <typename T>
-    [[nodiscard]] constexpr T mix(T a, T b, T t) noexcept {
-        return a + t * (b - a);
-    }
-}
-
-/// Smooth union with blend radius k
-/// Falls back to sharp union when k <= 0
 [[nodiscard]] inline double op_smooth_union(double d1, double d2, double k) {
-    if (k <= 0.0) return op_union(d1, d2);
+    if (k <= 0.0) return std::min(d1, d2);
     double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
-    return detail::mix(d2, d1, h) - k * h * (1.0 - h);
+    return d1 * (1.0 - h) + d2 * h - k * h * (1.0 - h);
 }
 
-/// Smooth intersection with blend radius k
-/// Falls back to sharp intersection when k <= 0
 [[nodiscard]] inline double op_smooth_intersection(double d1, double d2, double k) {
-    if (k <= 0.0) return op_intersection(d1, d2);
+    if (k <= 0.0) return std::max(d1, d2);
     double h = std::clamp(0.5 - 0.5 * (d2 - d1) / k, 0.0, 1.0);
-    return detail::mix(d1, d2, h) + k * h * (1.0 - h);
+    return d1 * (1.0 - h) + d2 * h + k * h * (1.0 - h);
 }
 
-/// Smooth difference with blend radius k (d1 \ d2)
-/// Falls back to sharp difference when k <= 0
 [[nodiscard]] inline double op_smooth_difference(double d1, double d2, double k) {
-    if (k <= 0.0) return op_difference(d1, d2);
+    if (k <= 0.0) return std::max(-d1, d2);
     double h = std::clamp(0.5 - 0.5 * (d2 + d1) / k, 0.0, 1.0);
-    return detail::mix(d2, -d1, h) + k * h * (1.0 - h);
+    return (-d1) * (1.0 - h) + d2 * h + k * h * (1.0 - h);
 }
 
-// ═══════════════════════════════════════════════════
-//  Shape Modifiers
-// ═══════════════════════════════════════════════════
+// ── Modifiers (operate on a distance value) ──
 
-/// Round a shape — shrink-wraps by radius r, then expands
-/// Equivalent to: shape_distance - r
 [[nodiscard]] inline double op_round(double d, double r) {
     return d - r;
 }
 
-/// Onion — create a thin shell of thickness t
-/// Returns SDF of the shell region (t/2 inward, t/2 outward from surface)
 [[nodiscard]] inline double op_onion(double d, double thickness) {
-    return std::abs(d) - thickness * 0.5;
+    return std::abs(d) - thickness;
 }
 
-// ═══════════════════════════════════════════════════
-//  Domain Deformations (coordinate remapping)
-//  Apply BEFORE evaluating the SDF primitive.
-// ═══════════════════════════════════════════════════
+// ── Domain deformation operators (transform point space) ──
 
-/// Infinite repetition along each axis
-/// c.x/y/z: cell size per axis (0 = no repeat along that axis)
-[[nodiscard]] inline vde::core::Point3D op_repeat(const vde::core::Point3D& p, const vde::core::Point3D& c) {
-    // p - c * round(p / c)
-    // For axes where c_i == 0, skip repetition
-    auto q = p;
-    if (c.x() > 0.0) q.x() = p.x() - c.x() * std::round(p.x() / c.x());
-    if (c.y() > 0.0) q.y() = p.y() - c.y() * std::round(p.y() / c.y());
-    if (c.z() > 0.0) q.z() = p.z() - c.z() * std::round(p.z() / c.z());
-    return q;
+/// Infinite repetition in a cell
+[[nodiscard]] inline Point3D op_repeat(const Point3D& p, const Point3D& cell) {
+    auto wrap = [](double v, double c) {
+        return c > 0.0 ? v - c * std::round(v / c) : v;
+    };
+    return Point3D(wrap(p.x(), cell.x()),
+                   wrap(p.y(), cell.y()),
+                   wrap(p.z(), cell.z()));
 }
 
-/// Mirror across X=0 plane (folds everything to x >= offset)
-[[nodiscard]] inline vde::core::Point3D op_mirror_x(const vde::core::Point3D& p, double offset = 0.0) {
-    return {offset + std::abs(p.x() - offset), p.y(), p.z()};
+/// Mirror across X=0 plane, with optional offset
+[[nodiscard]] inline Point3D op_mirror_x(const Point3D& p, double offset = 0.0) {
+    return Point3D(offset + std::abs(p.x() - offset), p.y(), p.z());
 }
 
-/// Mirror across Y=0 plane (folds everything to y >= offset)
-[[nodiscard]] inline vde::core::Point3D op_mirror_y(const vde::core::Point3D& p, double offset = 0.0) {
-    return {p.x(), offset + std::abs(p.y() - offset), p.z()};
+/// Mirror across Y=0
+[[nodiscard]] inline Point3D op_mirror_y(const Point3D& p) {
+    return Point3D(p.x(), std::abs(p.y()), p.z());
 }
 
-/// Mirror across Z=0 plane (folds everything to z >= offset)
-[[nodiscard]] inline vde::core::Point3D op_mirror_z(const vde::core::Point3D& p, double offset = 0.0) {
-    return {p.x(), p.y(), offset + std::abs(p.z() - offset)};
+/// Mirror across Z=0
+[[nodiscard]] inline Point3D op_mirror_z(const Point3D& p) {
+    return Point3D(p.x(), p.y(), std::abs(p.z()));
 }
 
-/// Rotate point around Y axis by angle (radians)
-[[nodiscard]] inline vde::core::Point3D op_rotate(const vde::core::Point3D& p, double angle_rad) {
-    double s = std::sin(angle_rad);
-    double c = std::cos(angle_rad);
-    return {c * p.x() + s * p.z(), p.y(), -s * p.x() + c * p.z()};
+/// Translate space (inverse for SDF: move the world, not the object)
+[[nodiscard]] inline Point3D op_translate(const Point3D& p, const Point3D& offset) {
+    return p - offset;
 }
 
-/// Translate point
-[[nodiscard]] inline vde::core::Point3D op_translate(const vde::core::Point3D& p, const vde::core::Point3D& offset) {
-    return {p.x() - offset.x(), p.y() - offset.y(), p.z() - offset.z()};
+/// Rotate around Y axis by angle_rad (inverse for SDF)
+[[nodiscard]] inline Point3D op_rotate(const Point3D& p, double angle_rad) {
+    double c = std::cos(-angle_rad);
+    double s = std::sin(-angle_rad);
+    return Point3D(p.x() * c - p.z() * s, p.y(), p.x() * s + p.z() * c);
 }
 
-/// Non-uniform scale — divides coordinates (non-exact SDF scaling)
-/// Use with care: the resulting distance field is approximate.
-[[nodiscard]] inline vde::core::Point3D op_scale(const vde::core::Point3D& p, const vde::core::Point3D& s) {
-    return {p.x() / s.x(), p.y() / s.y(), p.z() / s.z()};
+/// Non-uniform scale (inverse for SDF)
+[[nodiscard]] inline Point3D op_scale(const Point3D& p, const Point3D& s) {
+    return Point3D(p.x() / s.x(), p.y() / s.y(), p.z() / s.z());
 }
 
-/// Twist around Y axis — rotation amount increases with y coordinate
-[[nodiscard]] inline vde::core::Point3D op_twist(const vde::core::Point3D& p, double amount) {
-    double angle = amount * p.y();
-    double s = std::sin(angle);
-    double c = std::cos(angle);
-    return {c * p.x() + s * p.z(), p.y(), -s * p.x() + c * p.z()};
-}
-
-/// Bend along Y axis — curves the shape in x-y plane
-[[nodiscard]] inline vde::core::Point3D op_bend(const vde::core::Point3D& p, double amount) {
+/// Twist space around Y axis
+[[nodiscard]] inline Point3D op_twist(const Point3D& p, double amount) {
     double c = std::cos(amount * p.y());
     double s = std::sin(amount * p.y());
-    return {c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z()};
+    return Point3D(p.x() * c - p.z() * s, p.y(), p.x() * s + p.z() * c);
 }
 
-/// Elongate along axes — stretches coordinates (divides by factor)
-/// Equivalent to non-uniform scale but conceptually "stretches" the shape.
-[[nodiscard]] inline vde::core::Point3D op_elongate(const vde::core::Point3D& p, const vde::core::Point3D& factors) {
-    return {p.x() / factors.x(), p.y() / factors.y(), p.z() / factors.z()};
-}
-
-/// Displacement — adds sinusoidal perturbation to the SDF value
-/// p: current point, d: base SDF value, amplitude/frequency: displacement params
-[[nodiscard]] inline double op_displace(double d, const vde::core::Point3D& p, double amplitude, double frequency) {
-    return d + amplitude * std::sin(p.x() * frequency)
-                           * std::sin(p.y() * frequency)
-                           * std::sin(p.z() * frequency);
-}
-
-/// Cheap bend — bends around Z axis based on x coordinate
-/// Simpler and faster approximation than full bend.
-[[nodiscard]] inline vde::core::Point3D op_cheap_bend(const vde::core::Point3D& p, double k) {
+/// Bend space
+[[nodiscard]] inline Point3D op_bend(const Point3D& p, double k) {
     double c = std::cos(k * p.x());
     double s = std::sin(k * p.x());
-    return {c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z()};
+    return Point3D(c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z());
+}
+
+/// Elongate (stretch space, eliminating interior along major axes)
+[[nodiscard]] inline Point3D op_elongate(const Point3D& p, const Point3D& h) {
+    Point3D q = p.cwiseAbs() - h;
+    return Point3D(
+        q.x() > 0.0 ? p.x() - std::copysign(h.x(), p.x()) : 0.0,
+        q.y() > 0.0 ? p.y() - std::copysign(h.y(), p.y()) : 0.0,
+        q.z() > 0.0 ? p.z() - std::copysign(h.z(), p.z()) : 0.0
+    );
+}
+
+/// Cheap bend (simplified bending)
+[[nodiscard]] inline Point3D op_cheap_bend(const Point3D& p, double k) {
+    double c = std::cos(k * p.x());
+    double s = std::sin(k * p.x());
+    return Point3D(p.x(), p.y() * c - p.z() * s, p.y() * s + p.z() * c);
+}
+
+/// Displace distance field with sinusoidal noise
+[[nodiscard]] inline double op_displace(double d, const Point3D& p,
+                                         double amplitude, double frequency) {
+    double noise = std::sin(p.x() * frequency) *
+                   std::sin(p.y() * frequency) *
+                   std::sin(p.z() * frequency);
+    return d + amplitude * noise;
 }
 
 } // namespace vde::sdf

include/vde/sdf/sdf_primitives.h

@@ -0,0 +1,163 @@
+#pragma once
+#include "vde/core/point.h"
+#include <cmath>
+#include <algorithm>
+
+namespace vde::sdf {
+using core::Point3D;
+using core::Vector3D;
+
+// ═══════════════════════════════════════════════
+//  Basic Primitives
+// ═══════════════════════════════════════════════
+
+/// Signed distance to sphere centered at origin
+[[nodiscard]] inline double sphere(const Point3D& p, double radius) {
+    return p.norm() - radius;
+}
+
+/// Signed distance to axis-aligned box centered at origin
+/// half_extents = (hx, hy, hz) — half-dimensions in each axis
+[[nodiscard]] inline double box(const Point3D& p, const Point3D& half_extents) {
+    Point3D q = p.cwiseAbs() - half_extents;
+    return q.cwiseMax(0.0).norm()
+         + std::min(std::max({q.x(), q.y(), q.z()}), 0.0);
+}
+
+/// Signed distance to rounded box centered at origin
+/// Chamfer radius r is subtracted from the box distance
+[[nodiscard]] inline double round_box(const Point3D& p, const Point3D& half_extents, double r) {
+    Point3D q = p.cwiseAbs() - half_extents;
+    return q.cwiseMax(0.0).norm()
+         + std::min(std::max({q.x(), q.y(), q.z()}), 0.0)
+         - r;
+}
+
+/// Signed distance to torus in XZ plane, Y-axis symmetry
+/// major_radius = distance from origin to tube center
+/// minor_radius = tube thickness radius
+[[nodiscard]] inline double torus(const Point3D& p, double major_radius, double minor_radius) {
+    double qx = std::sqrt(p.x() * p.x() + p.z() * p.z()) - major_radius;
+    return std::sqrt(qx * qx + p.y() * p.y()) - minor_radius;
+}
+
+/// Signed distance to capsule (line segment swept with sphere of given radius)
+/// a, b = segment endpoints
+[[nodiscard]] inline double capsule(const Point3D& p, const Point3D& a, const Point3D& b, double radius) {
+    Point3D pa = p - a;
+    Point3D ba = b - a;
+    double h = std::clamp(pa.dot(ba) / ba.squaredNorm(), 0.0, 1.0);
+    return (pa - ba * h).norm() - radius;
+}
+
+/// Signed distance to capped cylinder along Y axis, centered at origin
+/// radius = cylinder radius, height = full height (from -h/2 to +h/2)
+[[nodiscard]] inline double cylinder(const Point3D& p, double radius, double height) {
+    double d_xz = std::sqrt(p.x() * p.x() + p.z() * p.z()) - radius;
+    double d_y  = std::abs(p.y()) - height * 0.5;
+    return std::min(std::max(d_xz, d_y), 0.0)
+         + std::sqrt(std::max(d_xz, 0.0) * std::max(d_xz, 0.0)
+                   + std::max(d_y,  0.0) * std::max(d_y,  0.0));
+}
+
+/// Signed distance to plane through origin with given normal
+/// Normal must be unit length; offset shifts plane along normal
+[[nodiscard]] inline double plane(const Point3D& p, const Vector3D& normal, double offset) {
+    return p.dot(normal) - offset;
+}
+
+/// Signed distance to ellipsoid centered at origin
+/// radii = semi-axis lengths (rx, ry, rz)
+/// Uses the standard bounded-gradient approximation: (|p/radii| - 1) * min(radii)
+[[nodiscard]] inline double ellipsoid(const Point3D& p, const Point3D& radii) {
+    Point3D pr(p.x() / radii.x(), p.y() / radii.y(), p.z() / radii.z());
+    double k0 = pr.norm();
+    double k1 = Point3D(p.x() / (radii.x() * radii.x()),
+                         p.y() / (radii.y() * radii.y()),
+                         p.z() / (radii.z() * radii.z())).norm();
+    return k0 * (k0 - 1.0) / k1;
+}
+
+// ═══════════════════════════════════════════════
+//  Parametric Primitives
+// ═══════════════════════════════════════════════
+
+/// Signed distance to regular hexagonal prism in XZ plane, extruded along Y
+/// radius = circumradius of hexagon (center to vertex)
+/// height = full extrusion height (±h/2 along Y)
+[[nodiscard]] inline double hex_prism(const Point3D& p, double radius, double height) {
+    double px = std::abs(p.x());
+    double pz = std::abs(p.z());
+    // 2D hexagon SDF (pointy-top in XZ plane)
+    double d_hex = std::max(px + pz * 0.577350269189626, pz * 1.154700538379252) - radius;
+    // Extrude along Y
+    return std::max(d_hex, std::abs(p.y()) - height * 0.5);
+}
+
+/// Signed distance to infinite cylinder through origin along given axis
+/// Axis must be unit length
+[[nodiscard]] inline double infinite_cylinder(const Point3D& p, const Vector3D& axis, double radius) {
+    return p.cross(axis).norm() - radius;
+}
+
+/// Signed distance to wedge — half of a box cut diagonally in XZ plane
+/// w = width (X), h = height (Y), d = depth (Z)
+/// The wedge occupies the half where z ≥ x inside the box
+[[nodiscard]] inline double wedge(const Point3D& p, double w, double h, double d) {
+    Point3D q = p.cwiseAbs() - Point3D(w * 0.5, h * 0.5, d * 0.5);
+    double d_box = q.cwiseMax(0.0).norm()
+                 + std::min(std::max({q.x(), q.y(), q.z()}), 0.0);
+    // Diagonal plane x - z ≤ 0  (i.e., z ≥ x)
+    double d_cut = (p.x() - p.z()) * 0.707106781186548;  // 1/√2
+    return std::max(d_box, d_cut);
+}
+
+// ═══════════════════════════════════════════════
+//  2D Extrusions & Revolution
+// ═══════════════════════════════════════════════
+
+/// Extrude a 2D SDF infinitely along Z axis
+/// sdf_2d = pre-computed 2D SDF value at (p.x, p.y)
+[[nodiscard]] inline double extrusion(const Point3D& /*p*/, double sdf_2d) {
+    return sdf_2d;
+}
+
+/// Bounded extrusion: 2D SDF intersected with Z slab ±half_height
+/// sdf_2d = pre-computed 2D SDF value at (p.x, p.y)
+[[nodiscard]] inline double extrusion_bounded(const Point3D& p, double sdf_2d, double half_height) {
+    return std::max(sdf_2d, std::abs(p.z()) - half_height);
+}
+
+/// Revolve a 2D profile around Y axis
+/// sdf_2d = pre-computed 2D SDF value at (|p.xz|-offset, p.y)
+/// offset  = radial offset from Y axis
+[[nodiscard]] inline double revolution(const Point3D& /*p*/, double sdf_2d, double /*offset*/) {
+    return sdf_2d;
+}
+
+// ═══════════════════════════════════════════════
+//  Complex Primitives (implemented in src/sdf/)
+// ═══════════════════════════════════════════════
+
+/// Signed distance to capped cone, centered at origin along Y axis
+/// Apex at y = +height/2, base at y = -height/2
+/// angle_rad = half-angle at apex, height = full height
+[[nodiscard]] double cone(const Point3D& p, double angle_rad, double height);
+
+/// Signed distance to triangular prism with arbitrary triangle cross-section
+/// Triangle defined by vertices a,b,c in XY plane, extruded ±height/2 along Z
+[[nodiscard]] double triangular_prism(const Point3D& p, const Point3D& a, const Point3D& b,
+                                      const Point3D& c, double height);
+
+/// Signed distance to link (two parallel toruses connected along X)
+/// length = spacing between torus centers along X
+/// major_r = major radius of each torus
+/// minor_r = minor (tube) radius of each torus
+[[nodiscard]] double link(const Point3D& p, double length, double major_r, double minor_r);
+
+/// Signed distance to infinite cone from apex along axis
+/// axis must be unit length, angle_rad = half-angle
+[[nodiscard]] double infinite_cone(const Point3D& p, const Point3D& apex,
+                                   const Vector3D& axis, double angle_rad);
+
+} // namespace vde::sdf

include/vde/sdf/sdf_to_mesh.h

@@ -0,0 +1,29 @@
+#pragma once
+#include "vde/sdf/sdf_tree.h"
+#include "vde/mesh/marching_cubes.h"
+#include <functional>
+
+namespace vde::sdf {
+
+/// Convert an SDF expression tree to a triangle mesh using marching cubes
+/// @param root       SDF expression tree root node
+/// @param resolution Grid resolution per axis (e.g. 64 → 64³ grid)
+/// @param iso_level  Isosurface level (0 = surface, default)
+/// @return Mesh with vertices and triangle indices
+[[nodiscard]] mesh::MCMesh sdf_to_mesh(const SdfNodePtr& root,
+                                        int resolution = 64,
+                                        double iso_level = 0.0);
+
+/// Convert a lambda SDF f(x,y,z) to a triangle mesh
+/// Convenience overload for direct function binding (useful from Python)
+/// @param f          SDF function f(x,y,z) → signed distance
+/// @param bmin       Lower corner of bounding box
+/// @param bmax       Upper corner of bounding box
+/// @param resolution Grid resolution per axis
+/// @param iso_level  Isosurface level
+[[nodiscard]] mesh::MCMesh sdf_to_mesh_lambda(
+    const std::function<double(double, double, double)>& f,
+    const Point3D& bmin, const Point3D& bmax,
+    int resolution = 64, double iso_level = 0.0);
+
+} // namespace vde::sdf

include/vde/sdf/sdf_tree.h

@@ -0,0 +1,140 @@
+#pragma once
+#include "vde/sdf/sdf_primitives.h"
+#include "vde/sdf/sdf_operations.h"
+#include "vde/core/point.h"
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace vde::sdf {
+using core::Point3D;
+using core::Vector3D;
+
+// Forward declaration
+class SdfNode;
+using SdfNodePtr = std::shared_ptr<SdfNode>;
+
+/// All node types in the CSG tree
+enum class SdfOp : uint8_t {
+    // Primitives
+    Sphere, Box, RoundBox, Torus, Capsule, Cylinder, Cone, Plane,
+    Ellipsoid, TriangularPrism, HexPrism, Link, Wedge,
+    // Boolean CSG
+    Union, Intersection, Difference,
+    SmoothUnion, SmoothIntersection, SmoothDifference,
+    // Modifiers
+    Round, Onion,
+    // Domain transforms
+    Repeat, MirrorX, MirrorY, MirrorZ, Translate, Rotate, Scale,
+    Twist, Bend, Elongate, CheapBend,
+    // Displacement
+    Displace
+};
+
+/// Parameter payload — one struct to rule them all
+struct SdfParams {
+    // Primitives
+    double radius = 1.0;
+    Point3D extents = Point3D(1, 1, 1);
+    Point3D pt_a = Point3D(0, -1, 0);
+    Point3D pt_b = Point3D(0, 1, 0);
+    double major_radius = 1.0;
+    double minor_radius = 0.3;
+    double angle_rad = 0.5;
+    double height = 2.0;
+    double thickness = 0.1;
+    Vector3D normal = Vector3D(0, 1, 0);
+    double offset = 0.0;
+    // Operations
+    double blend_k = 0.2;
+    double amount = 0.5;
+    double amplitude = 0.1;
+    double frequency = 1.0;
+    Point3D repeat_cell = Point3D(2, 2, 2);
+    Point3D translate_offset = Point3D(0, 0, 0);
+    Point3D scale_factors = Point3D(1, 1, 1);
+};
+
+/// A single node in the SDF expression tree
+class SdfNode {
+public:
+    SdfOp op;
+    SdfParams params;
+    std::vector<SdfNodePtr> children;
+    std::string name;
+
+    // ── Primitive factories ──
+    static SdfNodePtr sphere(double r);
+    static SdfNodePtr box(const Point3D& half_extents);
+    static SdfNodePtr round_box(const Point3D& half_extents, double r);
+    static SdfNodePtr cylinder(double r, double h);
+    static SdfNodePtr torus(double major_r, double minor_r);
+    static SdfNodePtr capsule(const Point3D& a, const Point3D& b, double r);
+    static SdfNodePtr cone(double angle_rad, double h);
+    static SdfNodePtr plane(const Vector3D& normal, double offset);
+    static SdfNodePtr ellipsoid(const Point3D& radii);
+    static SdfNodePtr triangular_prism(double h);
+    static SdfNodePtr hex_prism(double h);
+    static SdfNodePtr link(double r, double length, double thickness);
+    static SdfNodePtr wedge(const Point3D& extents);
+
+    // ── Boolean CSG factories (two children) ──
+    static SdfNodePtr op_union(SdfNodePtr a, SdfNodePtr b);
+    static SdfNodePtr op_intersection(SdfNodePtr a, SdfNodePtr b);
+    static SdfNodePtr op_difference(SdfNodePtr a, SdfNodePtr b);
+    static SdfNodePtr smooth_union(SdfNodePtr a, SdfNodePtr b, double k);
+    static SdfNodePtr smooth_intersection(SdfNodePtr a, SdfNodePtr b, double k);
+    static SdfNodePtr smooth_difference(SdfNodePtr a, SdfNodePtr b, double k);
+
+    // ── Modifier factories (one child) ──
+    static SdfNodePtr round(SdfNodePtr child, double r);
+    static SdfNodePtr onion(SdfNodePtr child, double thickness);
+
+    // ── Domain transform factories (one child) ──
+    static SdfNodePtr repeat(SdfNodePtr child, const Point3D& cell);
+    static SdfNodePtr mirror_x(SdfNodePtr child);
+    static SdfNodePtr mirror_y(SdfNodePtr child);
+    static SdfNodePtr mirror_z(SdfNodePtr child);
+    static SdfNodePtr translate(SdfNodePtr child, const Point3D& offset);
+    static SdfNodePtr rotate(SdfNodePtr child, double angle_rad);
+    static SdfNodePtr scale(SdfNodePtr child, const Point3D& factors);
+    static SdfNodePtr twist(SdfNodePtr child, double amount);
+    static SdfNodePtr bend(SdfNodePtr child, double amount);
+    static SdfNodePtr elongate(SdfNodePtr child, const Point3D& h);
+    static SdfNodePtr cheap_bend(SdfNodePtr child, double amount);
+    static SdfNodePtr displace(SdfNodePtr child, double amplitude, double frequency);
+
+    /// Visit all nodes in the tree (pre-order)
+    template <typename Fn>
+    void visit(Fn&& fn) {
+        fn(*this);
+        for (auto& c : children) c->visit(std::forward<Fn>(fn));
+    }
+
+    /// Visit all nodes in the tree (const pre-order)
+    template <typename Fn>
+    void visit(Fn&& fn) const {
+        fn(*this);
+        for (auto& c : children) c->visit(std::forward<Fn>(fn));
+    }
+
+private:
+    explicit SdfNode(SdfOp o) : op(o) {}
+};
+
+/// Evaluate the SDF tree at a point in world space
+/// Returns signed distance: negative = inside, positive = outside
+[[nodiscard]] double evaluate(const SdfNodePtr& root, const Point3D& p);
+
+/// Estimate bounding box of the SDF tree
+/// Returns the half-extent from center to corner (bounding radius)
+[[nodiscard]] Point3D estimate_bounds(const SdfNodePtr& root, double margin = 1.0);
+
+/// Compute full AABB from estimate_bounds result
+struct SdfBBox {
+    Point3D min;
+    Point3D max;
+};
+[[nodiscard]] SdfBBox estimate_bbox(const SdfNodePtr& root, double margin = 1.0);
+
+} // namespace vde::sdf