diff --git a/include/vde/sdf/sdf_operations.h b/include/vde/sdf/sdf_operations.h index 919d594..9e5522b 100644 --- a/include/vde/sdf/sdf_operations.h +++ b/include/vde/sdf/sdf_operations.h @@ -1,165 +1,138 @@ #pragma once +#include "vde/core/point.h" #include #include -#include "vde/core/point.h" +#include 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 - [[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 diff --git a/include/vde/sdf/sdf_primitives.h b/include/vde/sdf/sdf_primitives.h new file mode 100644 index 0000000..e03fc5b --- /dev/null +++ b/include/vde/sdf/sdf_primitives.h @@ -0,0 +1,163 @@ +#pragma once +#include "vde/core/point.h" +#include +#include + +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 diff --git a/include/vde/sdf/sdf_to_mesh.h b/include/vde/sdf/sdf_to_mesh.h new file mode 100644 index 0000000..2de6b51 --- /dev/null +++ b/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 + +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& f, + const Point3D& bmin, const Point3D& bmax, + int resolution = 64, double iso_level = 0.0); + +} // namespace vde::sdf diff --git a/include/vde/sdf/sdf_tree.h b/include/vde/sdf/sdf_tree.h new file mode 100644 index 0000000..0c7d7d0 --- /dev/null +++ b/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 +#include +#include + +namespace vde::sdf { +using core::Point3D; +using core::Vector3D; + +// Forward declaration +class SdfNode; +using SdfNodePtr = std::shared_ptr; + +/// 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 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 + void visit(Fn&& fn) { + fn(*this); + for (auto& c : children) c->visit(std::forward(fn)); + } + + /// Visit all nodes in the tree (const pre-order) + template + void visit(Fn&& fn) const { + fn(*this); + for (auto& c : children) c->visit(std::forward(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 diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index a410ae2..4e5f70d 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -15,6 +15,7 @@ pybind11_add_module(_vde src/bind_core.cpp src/bind_curves.cpp src/bind_mesh.cpp + src/bind_sdf.cpp ) target_include_directories(_vde diff --git a/python/src/bind_module.cpp b/python/src/bind_module.cpp index dadd7e0..36bdaf1 100644 --- a/python/src/bind_module.cpp +++ b/python/src/bind_module.cpp @@ -8,6 +8,7 @@ namespace py = pybind11; void bind_core(py::module& m); void bind_curves(py::module& m); void bind_mesh(py::module& m); +void bind_sdf(py::module& m); PYBIND11_MODULE(_vde, m) { m.doc() = "ViewDesignEngine — CAD computational geometry engine"; @@ -15,4 +16,5 @@ PYBIND11_MODULE(_vde, m) { bind_core(m); bind_curves(m); bind_mesh(m); + bind_sdf(m); } diff --git a/python/src/bind_sdf.cpp b/python/src/bind_sdf.cpp new file mode 100644 index 0000000..90b0b9e --- /dev/null +++ b/python/src/bind_sdf.cpp @@ -0,0 +1,151 @@ +#include +#include +#include +#include + +#include "vde/sdf/sdf_tree.h" +#include "vde/sdf/sdf_to_mesh.h" + +namespace py = pybind11; +using namespace vde::sdf; + +void bind_sdf(py::module& m) { + auto sdf_m = m.def_submodule("sdf", + "Signed Distance Field implicit modeling"); + + // ── SdfNode ──────────────────────────────────────── + py::class_>(sdf_m, "SdfNode", + R"pbdoc( +A node in a signed-distance-field CSG expression tree. + +SdfNode represents either a primitive shape, a boolean CSG operation, +a distance modifier, or a domain transform. Build trees by combining +nodes with the factory functions. +)pbdoc") + // Primitive factories + .def_static("sphere", &SdfNode::sphere, + py::arg("radius"), "Sphere centered at origin") + .def_static("box", &SdfNode::box, + py::arg("half_extents"), "Axis-aligned box") + .def_static("round_box", &SdfNode::round_box, + py::arg("half_extents"), py::arg("radius"), "Rounded box") + .def_static("cylinder", &SdfNode::cylinder, + py::arg("radius"), py::arg("height"), "Capped cylinder along Y") + .def_static("torus", &SdfNode::torus, + py::arg("major_radius"), py::arg("minor_radius"), "Torus in XZ plane") + .def_static("capsule", &SdfNode::capsule, + py::arg("a"), py::arg("b"), py::arg("radius"), "Capsule from a to b") + .def_static("cone", &SdfNode::cone, + py::arg("angle_rad"), py::arg("height"), "Cone along Y axis") + .def_static("plane", &SdfNode::plane, + py::arg("normal"), py::arg("offset"), "Infinite plane") + .def_static("ellipsoid", &SdfNode::ellipsoid, + py::arg("radii"), "Ellipsoid centered at origin") + .def_static("triangular_prism", &SdfNode::triangular_prism, + py::arg("height"), "Triangular prism along Y") + .def_static("hex_prism", &SdfNode::hex_prism, + py::arg("height"), "Hexagonal prism along Y") + .def_static("link", &SdfNode::link, + py::arg("radius"), py::arg("length"), py::arg("thickness"), "Link segment") + .def_static("wedge", &SdfNode::wedge, + py::arg("extents"), "Wedge shape") + + // Boolean CSG operations + .def_static("union_", &SdfNode::op_union, + py::arg("a"), py::arg("b"), "CSG union (min)") + .def_static("intersection", &SdfNode::op_intersection, + py::arg("a"), py::arg("b"), "CSG intersection (max)") + .def_static("difference", &SdfNode::op_difference, + py::arg("a"), py::arg("b"), "CSG difference (max(a, -b))") + .def_static("smooth_union", &SdfNode::smooth_union, + py::arg("a"), py::arg("b"), py::arg("k"), "Smooth blend union") + .def_static("smooth_intersection", &SdfNode::smooth_intersection, + py::arg("a"), py::arg("b"), py::arg("k"), "Smooth blend intersection") + .def_static("smooth_difference", &SdfNode::smooth_difference, + py::arg("a"), py::arg("b"), py::arg("k"), "Smooth blend difference") + + // Modifiers + .def_static("round", &SdfNode::round, + py::arg("child"), py::arg("radius"), "Offset surface outward") + .def_static("onion", &SdfNode::onion, + py::arg("child"), py::arg("thickness"), "Hollow shell") + + // Domain transforms + .def_static("repeat", &SdfNode::repeat, + py::arg("child"), py::arg("cell"), "Infinite repetition") + .def_static("mirror_x", &SdfNode::mirror_x, + py::arg("child"), "Mirror across YZ plane") + .def_static("mirror_y", &SdfNode::mirror_y, + py::arg("child"), "Mirror across XZ plane") + .def_static("mirror_z", &SdfNode::mirror_z, + py::arg("child"), "Mirror across XY plane") + .def_static("translate", &SdfNode::translate, + py::arg("child"), py::arg("offset"), "Translate in space") + .def_static("rotate", &SdfNode::rotate, + py::arg("child"), py::arg("angle_rad"), "Rotate around Y axis") + .def_static("scale", &SdfNode::scale, + py::arg("child"), py::arg("factors"), "Non-uniform scale") + .def_static("twist", &SdfNode::twist, + py::arg("child"), py::arg("amount"), "Twist around Y axis") + .def_static("bend", &SdfNode::bend, + py::arg("child"), py::arg("amount"), "Bend around Z axis") + .def_static("elongate", &SdfNode::elongate, + py::arg("child"), py::arg("h"), "Elongate (stretch space)") + .def_static("cheap_bend", &SdfNode::cheap_bend, + py::arg("child"), py::arg("amount"), "Simplified bend") + .def_static("displace", &SdfNode::displace, + py::arg("child"), py::arg("amplitude"), py::arg("frequency"), + "Sinusoidal displacement") + + // Accessors + .def_readwrite("name", &SdfNode::name, "Optional label") + .def("__repr__", [](const SdfNode& n) { + return "SdfNode(" + n.name + ")"; + }); + + // ── evaluate() ───────────────────────────────────── + sdf_m.def("evaluate", &evaluate, + py::arg("node"), py::arg("point"), + "Evaluate SDF tree at point. Returns signed distance."); + + sdf_m.def("estimate_bounds", &estimate_bounds, + py::arg("root"), py::arg("margin") = 1.0, + "Estimate half-extent bounding radius of SDF tree"); + + // ── SdfBBox ────────────────────────────────────── + py::class_(sdf_m, "SdfBBox") + .def(py::init<>()) + .def_readwrite("min", &SdfBBox::min) + .def_readwrite("max", &SdfBBox::max); + + sdf_m.def("estimate_bbox", &estimate_bbox, + py::arg("root"), py::arg("margin") = 1.0, + "Estimate full AABB of SDF tree"); + + // ── sdf_to_mesh() ────────────────────────────────── + py::class_(sdf_m, "MCMesh", + R"pbdoc( +Marching cubes mesh result. + +Attributes: + vertices: N x 3 list of Point3D + triangles: M x 3 list of index triples +)pbdoc") + .def(py::init<>()) + .def_readwrite("vertices", &vde::mesh::MCMesh::vertices) + .def_readwrite("triangles", &vde::mesh::MCMesh::triangles) + .def("__repr__", [](const vde::mesh::MCMesh& m) { + return "MCMesh(v=" + std::to_string(m.vertices.size()) + + ", t=" + std::to_string(m.triangles.size()) + ")"; + }); + + sdf_m.def("sdf_to_mesh", &sdf_to_mesh, + py::arg("root"), py::arg("resolution") = 64, + py::arg("iso_level") = 0.0, + "Convert SDF tree to triangle mesh via marching cubes"); + + sdf_m.def("sdf_to_mesh_lambda", &sdf_to_mesh_lambda, + py::arg("f"), py::arg("bmin"), py::arg("bmax"), + py::arg("resolution") = 64, py::arg("iso_level") = 0.0, + "Convert lambda SDF f(x,y,z) to triangle mesh"); +} diff --git a/python/vde/__init__.py b/python/vde/__init__.py index a1ccb91..1cb02e9 100644 --- a/python/vde/__init__.py +++ b/python/vde/__init__.py @@ -2,11 +2,12 @@ ViewDesignEngine — CAD computational geometry engine for Python. Submodules: - vde.core — Point3D, Vector3D, AABB3D, Polygon2D, convex hull, distance, transforms + vde.core — Point3D, Vector3D, AABB3D, Polygon2D, convex hull, distance, transforms vde.curves — BezierCurve, BSplineCurve, NurbsCurve - vde.mesh — delaunay_2d, delaunay_3d, HalfedgeMesh + vde.mesh — delaunay_2d, delaunay_3d, HalfedgeMesh + vde.sdf — SdfNode, CSG tree, marching cubes mesh conversion """ -from ._vde import core, curves, mesh +from ._vde import core, curves, mesh, sdf __version__ = "1.0.0" -__all__ = ["core", "curves", "mesh"] +__all__ = ["core", "curves", "mesh", "sdf"] diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c3a522e..f1df21e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -179,6 +179,20 @@ target_link_libraries(vde_sketch PUBLIC vde_core vde_compile_options ) +# ── sdf ──────────────────────────────────────────── +add_library(vde_sdf STATIC + sdf/sdf_primitives.cpp + sdf/sdf_tree.cpp + sdf/sdf_to_mesh.cpp +) +target_include_directories(vde_sdf + PUBLIC ${CMAKE_SOURCE_DIR}/include + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} +) +target_link_libraries(vde_sdf + PUBLIC vde_core vde_mesh vde_compile_options +) + # ── C API ────────────────────────────────────────── add_library(vde_capi STATIC capi/vde_capi.cpp diff --git a/src/sdf/sdf_primitives.cpp b/src/sdf/sdf_primitives.cpp new file mode 100644 index 0000000..ab0a9f5 --- /dev/null +++ b/src/sdf/sdf_primitives.cpp @@ -0,0 +1,130 @@ +#include "vde/sdf/sdf_primitives.h" +#include +#include + +namespace vde::sdf { + +// ── Helper: 2D triangle SDF (Inigo Quilez) ────── +namespace { + +struct Vec2 { double x, y; }; + +inline double dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; } +inline double cross2(Vec2 a, Vec2 b) { return a.x * b.y - a.y * b.x; } +inline Vec2 sub(Vec2 a, Vec2 b) { return {a.x - b.x, a.y - b.y}; } +inline Vec2 madd(Vec2 a, Vec2 b, double t) { return {a.x - b.x * t, a.y - b.y * t}; } +inline double len2(Vec2 a) { return a.x * a.x + a.y * a.y; } + +[[nodiscard]] double triangle_sdf_2d(Vec2 p, Vec2 a, Vec2 b, Vec2 c) { + Vec2 e0 = sub(b, a); + Vec2 e1 = sub(c, b); + Vec2 e2 = sub(a, c); + Vec2 v0 = sub(p, a); + Vec2 v1 = sub(p, b); + Vec2 v2 = sub(p, c); + + double d0 = dot(e0, e0); + double d1 = dot(e1, e1); + double d2 = dot(e2, e2); + + Vec2 pq0 = madd(v0, e0, std::clamp(dot(v0, e0) / d0, 0.0, 1.0)); + Vec2 pq1 = madd(v1, e1, std::clamp(dot(v1, e1) / d1, 0.0, 1.0)); + Vec2 pq2 = madd(v2, e2, std::clamp(dot(v2, e2) / d2, 0.0, 1.0)); + + double s = cross2(e0, e2) > 0.0 ? 1.0 : -1.0; + + double dx0 = len2(pq0); + double dx1 = len2(pq1); + double dx2 = len2(pq2); + + double sy0 = s * cross2(v0, e0); + double sy1 = s * cross2(v1, e1); + double sy2 = s * cross2(v2, e2); + + // Find minimum (d², signed) pair + double best_d2 = dx0; + double best_sd = sy0; + if (dx1 < best_d2) { best_d2 = dx1; best_sd = sy1; } + if (dx2 < best_d2) { best_d2 = dx2; best_sd = sy2; } + + if (sy1 < best_sd) { best_sd = sy1; best_d2 = dx1; } + if (sy2 < best_sd) { best_sd = sy2; best_d2 = dx2; } + + return -std::sqrt(best_d2) * (best_sd > 0.0 ? 1.0 : -1.0); +} + +} // anonymous namespace + +// ═══ cone ════════════════════════════════════════ + +double cone(const Point3D& p, double angle_rad, double height) { + double r_base = height * std::tan(angle_rad); + double h2 = height * 0.5; + + // Shift so base is at y = 0, apex at y = height + double qx = std::sqrt(p.x() * p.x() + p.z() * p.z()); + double qy = p.y() + h2; + + double b = r_base / height; // slope: (r1 - r2) / h + if (b > 1.0) b = 1.0; // clamp for wide cones (>45°) + double a = std::sqrt(1.0 - b * b); + + // Projection parameter: k = dot(q, vec2(-b, a)) + double k = a * qy - b * qx; + + if (k < 0.0) { + // Closest to base circle + return std::sqrt(qx * qx + qy * qy) - r_base; + } + if (k > a * height) { + // Closest to apex + return std::sqrt((qx - r_base) * (qx - r_base) + (qy - height) * (qy - height)); + } + // Closest to side surface + return qx * a + qy * b - r_base; +} + +// ═══ triangular_prism ═══════════════════════════ + +double triangular_prism(const Point3D& p, const Point3D& a, const Point3D& b, + const Point3D& c, double height) { + // Project triangle to XY plane + Vec2 p2{p.x(), p.y()}; + Vec2 a2{a.x(), a.y()}; + Vec2 b2{b.x(), b.y()}; + Vec2 c2{c.x(), c.y()}; + + double d_xy = triangle_sdf_2d(p2, a2, b2, c2); + return std::max(d_xy, std::abs(p.z()) - height * 0.5); +} + +// ═══ link ═══════════════════════════════════════ + +double link(const Point3D& p, double length, double major_r, double minor_r) { + // Two parallel toruses offset along X + double half_len = length * 0.5; + + // Torus at +half_len + double px1 = p.x() - half_len; + double qx1 = std::sqrt(px1 * px1 + p.z() * p.z()) - major_r; + double d1 = std::sqrt(qx1 * qx1 + p.y() * p.y()) - minor_r; + + // Torus at -half_len + double px2 = p.x() + half_len; + double qx2 = std::sqrt(px2 * px2 + p.z() * p.z()) - major_r; + double d2 = std::sqrt(qx2 * qx2 + p.y() * p.y()) - minor_r; + + return std::min(d1, d2); +} + +// ═══ infinite_cone ══════════════════════════════ + +double infinite_cone(const Point3D& p, const Point3D& apex, + const Vector3D& axis, double angle_rad) { + Vector3D d = p - apex; + double proj = d.dot(axis); // projection onto axis + double perp = (d - axis * proj).norm(); // perpendicular distance + return perp * std::cos(angle_rad) - proj * std::sin(angle_rad); +} + +} // namespace vde::sdf diff --git a/src/sdf/sdf_to_mesh.cpp b/src/sdf/sdf_to_mesh.cpp new file mode 100644 index 0000000..adca434 --- /dev/null +++ b/src/sdf/sdf_to_mesh.cpp @@ -0,0 +1,23 @@ +#include "vde/sdf/sdf_to_mesh.h" + +namespace vde::sdf { + +mesh::MCMesh sdf_to_mesh(const SdfNodePtr& root, int resolution, double iso_level) { + SdfBBox bbox = estimate_bbox(root, 1.0); + + auto sdf_fn = [&root](double x, double y, double z) -> double { + return evaluate(root, Point3D(x, y, z)); + }; + + return mesh::marching_cubes(sdf_fn, iso_level, bbox.min, bbox.max, resolution); +} + +mesh::MCMesh sdf_to_mesh_lambda( + const std::function& f, + const Point3D& bmin, const Point3D& bmax, + int resolution, double iso_level) { + + return mesh::marching_cubes(f, iso_level, bmin, bmax, resolution); +} + +} // namespace vde::sdf diff --git a/src/sdf/sdf_tree.cpp b/src/sdf/sdf_tree.cpp new file mode 100644 index 0000000..86cf27c --- /dev/null +++ b/src/sdf/sdf_tree.cpp @@ -0,0 +1,377 @@ +#include "vde/sdf/sdf_tree.h" +#include +#include + +namespace vde::sdf { + +// ── Helper: make a node ────────────────────────────────────────── +static SdfNodePtr make(SdfOp op) { return std::make_shared(op); } +static void set_one_child(SdfNodePtr& n, SdfNodePtr c) { n->children = {std::move(c)}; } +static void set_two_children(SdfNodePtr& n, SdfNodePtr a, SdfNodePtr b) { + n->children = {std::move(a), std::move(b)}; +} + +// ── Primitive factories ────────────────────────────────────────── + +SdfNodePtr SdfNode::sphere(double r) { + auto n = make(SdfOp::Sphere); + n->params.radius = r; + n->name = "sphere"; + return n; +} + +SdfNodePtr SdfNode::box(const Point3D& half_extents) { + auto n = make(SdfOp::Box); + n->params.extents = half_extents; + n->name = "box"; + return n; +} + +SdfNodePtr SdfNode::round_box(const Point3D& half_extents, double r) { + auto n = make(SdfOp::RoundBox); + n->params.extents = half_extents; + n->params.radius = r; + n->name = "round_box"; + return n; +} + +SdfNodePtr SdfNode::cylinder(double r, double h) { + auto n = make(SdfOp::Cylinder); + n->params.radius = r; + n->params.height = h; // full height + n->name = "cylinder"; + return n; +} + +SdfNodePtr SdfNode::torus(double major_r, double minor_r) { + auto n = make(SdfOp::Torus); + n->params.major_radius = major_r; + n->params.minor_radius = minor_r; + n->name = "torus"; + return n; +} + +SdfNodePtr SdfNode::capsule(const Point3D& a, const Point3D& b, double r) { + auto n = make(SdfOp::Capsule); + n->params.pt_a = a; + n->params.pt_b = b; + n->params.radius = r; + n->name = "capsule"; + return n; +} + +SdfNodePtr SdfNode::cone(double angle_rad, double h) { + auto n = make(SdfOp::Cone); + n->params.angle_rad = angle_rad; + n->params.height = h; // full height + n->name = "cone"; + return n; +} + +SdfNodePtr SdfNode::plane(const Vector3D& normal, double offset) { + auto n = make(SdfOp::Plane); + n->params.normal = normal; + n->params.offset = offset; + n->name = "plane"; + return n; +} + +SdfNodePtr SdfNode::ellipsoid(const Point3D& radii) { + auto n = make(SdfOp::Ellipsoid); + n->params.extents = radii; + n->name = "ellipsoid"; + return n; +} + +SdfNodePtr SdfNode::triangular_prism(double h) { + auto n = make(SdfOp::TriangularPrism); + n->params.height = h; // full height + n->name = "triangular_prism"; + return n; +} + +SdfNodePtr SdfNode::hex_prism(double h) { + auto n = make(SdfOp::HexPrism); + n->params.radius = h; // circumradius + n->params.height = h; // full height (= radius for equilateral) + n->name = "hex_prism"; + return n; +} + +SdfNodePtr SdfNode::link(double r, double length, double thickness) { + auto n = make(SdfOp::Link); + n->params.radius = r; // major radius + n->params.height = length; // spacing + n->params.thickness = thickness; // minor radius + n->name = "link"; + return n; +} + +SdfNodePtr SdfNode::wedge(const Point3D& extents) { + auto n = make(SdfOp::Wedge); + n->params.extents = extents; + n->name = "wedge"; + return n; +} + +// ── Boolean CSG factories ──────────────────────────────────────── + +SdfNodePtr SdfNode::op_union(SdfNodePtr a, SdfNodePtr b) { + auto n = make(SdfOp::Union); n->name="union"; + set_two_children(n, std::move(a), std::move(b)); return n; +} +SdfNodePtr SdfNode::op_intersection(SdfNodePtr a, SdfNodePtr b) { + auto n = make(SdfOp::Intersection); n->name="intersection"; + set_two_children(n, std::move(a), std::move(b)); return n; +} +SdfNodePtr SdfNode::op_difference(SdfNodePtr a, SdfNodePtr b) { + auto n = make(SdfOp::Difference); n->name="difference"; + set_two_children(n, std::move(a), std::move(b)); return n; +} +SdfNodePtr SdfNode::smooth_union(SdfNodePtr a, SdfNodePtr b, double k) { + auto n = make(SdfOp::SmoothUnion); n->params.blend_k = k; n->name="smooth_union"; + set_two_children(n, std::move(a), std::move(b)); return n; +} +SdfNodePtr SdfNode::smooth_intersection(SdfNodePtr a, SdfNodePtr b, double k) { + auto n = make(SdfOp::SmoothIntersection); n->params.blend_k = k; n->name="smooth_intersection"; + set_two_children(n, std::move(a), std::move(b)); return n; +} +SdfNodePtr SdfNode::smooth_difference(SdfNodePtr a, SdfNodePtr b, double k) { + auto n = make(SdfOp::SmoothDifference); n->params.blend_k = k; n->name="smooth_difference"; + set_two_children(n, std::move(a), std::move(b)); return n; +} + +// ── Modifier factories ─────────────────────────────────────────── + +SdfNodePtr SdfNode::round(SdfNodePtr child, double r) { + auto n = make(SdfOp::Round); n->params.radius = r; n->name="round"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::onion(SdfNodePtr child, double thickness) { + auto n = make(SdfOp::Onion); n->params.thickness = thickness; n->name="onion"; + set_one_child(n, std::move(child)); return n; +} + +// ── Domain transform factories ─────────────────────────────────── + +SdfNodePtr SdfNode::repeat(SdfNodePtr child, const Point3D& cell) { + auto n = make(SdfOp::Repeat); n->params.repeat_cell = cell; n->name="repeat"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::mirror_x(SdfNodePtr child) { + auto n = make(SdfOp::MirrorX); n->name="mirror_x"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::mirror_y(SdfNodePtr child) { + auto n = make(SdfOp::MirrorY); n->name="mirror_y"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::mirror_z(SdfNodePtr child) { + auto n = make(SdfOp::MirrorZ); n->name="mirror_z"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::translate(SdfNodePtr child, const Point3D& offset) { + auto n = make(SdfOp::Translate); n->params.translate_offset = offset; n->name="translate"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::rotate(SdfNodePtr child, double angle_rad) { + auto n = make(SdfOp::Rotate); n->params.angle_rad = angle_rad; n->name="rotate"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::scale(SdfNodePtr child, const Point3D& factors) { + auto n = make(SdfOp::Scale); n->params.scale_factors = factors; n->name="scale"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::twist(SdfNodePtr child, double amount) { + auto n = make(SdfOp::Twist); n->params.amount = amount; n->name="twist"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::bend(SdfNodePtr child, double amount) { + auto n = make(SdfOp::Bend); n->params.amount = amount; n->name="bend"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::elongate(SdfNodePtr child, const Point3D& h) { + auto n = make(SdfOp::Elongate); n->params.extents = h; n->name="elongate"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::cheap_bend(SdfNodePtr child, double amount) { + auto n = make(SdfOp::CheapBend); n->params.amount = amount; n->name="cheap_bend"; + set_one_child(n, std::move(child)); return n; +} +SdfNodePtr SdfNode::displace(SdfNodePtr child, double amplitude, double frequency) { + auto n = make(SdfOp::Displace); + n->params.amplitude = amplitude; n->params.frequency = frequency; n->name="displace"; + set_one_child(n, std::move(child)); return n; +} + +// ── Recursive evaluator ────────────────────────────────────────── + +double evaluate(const SdfNodePtr& node, const Point3D& p) { + if (!node) return 1e30; + + const auto& pr = node->params; + + switch (node->op) { + // ── Primitives (from sdf_primitives.h) ── + case SdfOp::Sphere: + return sphere(p, pr.radius); + case SdfOp::Box: + return box(p, pr.extents); + case SdfOp::RoundBox: + return round_box(p, pr.extents, pr.radius); + case SdfOp::Torus: + return torus(p, pr.major_radius, pr.minor_radius); + case SdfOp::Capsule: + return capsule(p, pr.pt_a, pr.pt_b, pr.radius); + case SdfOp::Cylinder: + return cylinder(p, pr.radius, pr.height); + case SdfOp::Cone: + return cone(p, pr.angle_rad, pr.height); + case SdfOp::Plane: + return plane(p, pr.normal, pr.offset); + case SdfOp::Ellipsoid: + return ellipsoid(p, pr.extents); + case SdfOp::TriangularPrism: { + // Equilateral triangle in XZ plane + double rh = pr.height * 0.57735; + Point3D a(-rh, 0, -rh * 0.5), b(rh, 0, -rh * 0.5), c(0, 0, rh); + return triangular_prism(p, a, b, c, pr.height); + } + case SdfOp::HexPrism: + return hex_prism(p, pr.radius, pr.height); + case SdfOp::Link: + return link(p, pr.height, pr.radius, pr.thickness); + case SdfOp::Wedge: { + const auto& e = pr.extents; + return wedge(p, e.x() * 2.0, e.y() * 2.0, e.z() * 2.0); + } + + // ── Boolean CSG ── + case SdfOp::Union: + return op_union(evaluate(node->children[0], p), + evaluate(node->children[1], p)); + case SdfOp::Intersection: + return op_intersection(evaluate(node->children[0], p), + evaluate(node->children[1], p)); + case SdfOp::Difference: + return op_difference(evaluate(node->children[0], p), + evaluate(node->children[1], p)); + case SdfOp::SmoothUnion: + return op_smooth_union(evaluate(node->children[0], p), + evaluate(node->children[1], p), pr.blend_k); + case SdfOp::SmoothIntersection: + return op_smooth_intersection(evaluate(node->children[0], p), + evaluate(node->children[1], p), pr.blend_k); + case SdfOp::SmoothDifference: + return op_smooth_difference(evaluate(node->children[0], p), + evaluate(node->children[1], p), pr.blend_k); + + // ── Modifiers (distance post-process) ── + case SdfOp::Round: + return op_round(evaluate(node->children[0], p), pr.radius); + case SdfOp::Onion: + return op_onion(evaluate(node->children[0], p), pr.thickness); + + // ── Domain transforms (point pre-process) ── + case SdfOp::Repeat: + return evaluate(node->children[0], op_repeat(p, pr.repeat_cell)); + case SdfOp::MirrorX: + return evaluate(node->children[0], op_mirror_x(p)); + case SdfOp::MirrorY: + return evaluate(node->children[0], op_mirror_y(p)); + case SdfOp::MirrorZ: + return evaluate(node->children[0], op_mirror_z(p)); + case SdfOp::Translate: + return evaluate(node->children[0], op_translate(p, pr.translate_offset)); + case SdfOp::Rotate: + return evaluate(node->children[0], op_rotate(p, pr.angle_rad)); + case SdfOp::Scale: + return evaluate(node->children[0], op_scale(p, pr.scale_factors)); + case SdfOp::Twist: + return evaluate(node->children[0], op_twist(p, pr.amount)); + case SdfOp::Bend: + return evaluate(node->children[0], op_bend(p, pr.amount)); + case SdfOp::Elongate: + return evaluate(node->children[0], op_elongate(p, pr.extents)); + case SdfOp::CheapBend: + return evaluate(node->children[0], op_cheap_bend(p, pr.amount)); + + // ── Displacement ── + case SdfOp::Displace: + return op_displace(evaluate(node->children[0], p), p, + pr.amplitude, pr.frequency); + } + + return 1e30; +} + +// ── Bounds estimation ──────────────────────────────────────────── + +static double estimate_radius(const SdfNodePtr& node) { + if (!node) return 1.0; + + switch (node->op) { + case SdfOp::Sphere: + return node->params.radius; + case SdfOp::Box: + return node->params.extents.norm(); + case SdfOp::RoundBox: + return node->params.extents.norm() + node->params.radius; + case SdfOp::Torus: + return node->params.major_radius + node->params.minor_radius; + case SdfOp::Capsule: + return (node->params.pt_b - node->params.pt_a).norm() * 0.5 + node->params.radius; + case SdfOp::Cylinder: + return std::sqrt(node->params.radius * node->params.radius + + node->params.height * node->params.height * 0.25); + case SdfOp::Cone: + return node->params.height * 0.5 + node->params.height * 0.5 * std::tan(node->params.angle_rad); + case SdfOp::Plane: + return 10.0; + case SdfOp::Ellipsoid: + return std::max({node->params.extents.x(), + node->params.extents.y(), + node->params.extents.z()}); + case SdfOp::TriangularPrism: + case SdfOp::HexPrism: + return node->params.height; + case SdfOp::Link: + return node->params.radius + node->params.height * 0.5 + node->params.thickness; + case SdfOp::Wedge: + return node->params.extents.norm(); + case SdfOp::Union: + case SdfOp::Intersection: + case SdfOp::SmoothUnion: + case SdfOp::SmoothIntersection: + return std::max(estimate_radius(node->children[0]), + estimate_radius(node->children[1])); + case SdfOp::Difference: + case SdfOp::SmoothDifference: + return estimate_radius(node->children[0]); + case SdfOp::Round: + return estimate_radius(node->children[0]) + node->params.radius; + case SdfOp::Onion: + return estimate_radius(node->children[0]) + node->params.thickness; + case SdfOp::Translate: + return estimate_radius(node->children[0]) + node->params.translate_offset.norm(); + case SdfOp::Scale: + return estimate_radius(node->children[0]) * + std::max({node->params.scale_factors.x(), + node->params.scale_factors.y(), + node->params.scale_factors.z()}); + default: + return node->children.empty() ? 1.0 : estimate_radius(node->children[0]); + } +} + +Point3D estimate_bounds(const SdfNodePtr& root, double margin) { + double r = estimate_radius(root) + margin; + return Point3D(r, r, r); +} + +SdfBBox estimate_bbox(const SdfNodePtr& root, double margin) { + Point3D half = estimate_bounds(root, margin); + return SdfBBox{-half, half}; +} + +} // namespace vde::sdf diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a96febe..3fb9807 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -15,3 +15,4 @@ add_subdirectory(brep) add_subdirectory(sdf) add_subdirectory(sketch) add_subdirectory(sdf) +add_subdirectory(sdf) diff --git a/tests/sdf/CMakeLists.txt b/tests/sdf/CMakeLists.txt index c171a34..a53f637 100644 --- a/tests/sdf/CMakeLists.txt +++ b/tests/sdf/CMakeLists.txt @@ -1,9 +1,4 @@ -function(add_vde_sdf_test name) - add_executable(${name} ${name}.cpp) - target_link_libraries(${name} PRIVATE vde_sdf GTest::gtest GTest::gtest_main) - target_include_directories(${name} PRIVATE ${CMAKE_SOURCE_DIR}/include) - gtest_discover_tests(${name}) -endfunction() - -add_vde_sdf_test(test_sdf_tree) -add_vde_sdf_test(test_sdf_to_mesh) +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) diff --git a/tests/sdf/test_sdf_primitives.cpp b/tests/sdf/test_sdf_primitives.cpp new file mode 100644 index 0000000..9503dd3 --- /dev/null +++ b/tests/sdf/test_sdf_primitives.cpp @@ -0,0 +1,511 @@ +#include +#include "vde/sdf/sdf_primitives.h" +#include + +using namespace vde::sdf; +using vde::core::Point3D; +using vde::core::Vector3D; + +constexpr double EPS = 1e-9; +constexpr double SQRT2 = 1.4142135623730951; +constexpr double SQRT3 = 1.7320508075688772; + +// ═══════════════════════════════════════════════════ +// sphere +// ═══════════════════════════════════════════════════ + +TEST(SdfSphere, Outside) { + EXPECT_NEAR(sphere(Point3D(5, 0, 0), 3.0), 2.0, EPS); + EXPECT_NEAR(sphere(Point3D(0, 0, 0), 3.0), -3.0, EPS); +} + +TEST(SdfSphere, OnSurface) { + EXPECT_NEAR(sphere(Point3D(3, 0, 0), 3.0), 0.0, EPS); + EXPECT_NEAR(sphere(Point3D(0, -3, 0), 3.0), 0.0, EPS); +} + +TEST(SdfSphere, Inside) { + // Center of a radius-5 sphere: distance should be -5 + EXPECT_NEAR(sphere(Point3D(0, 0, 0), 5.0), -5.0, EPS); + // Partway in + EXPECT_NEAR(sphere(Point3D(2, 0, 0), 5.0), -3.0, EPS); +} + +TEST(SdfSphere, Degenerate) { + // Zero radius + EXPECT_NEAR(sphere(Point3D(1, 0, 0), 0.0), 1.0, EPS); + EXPECT_NEAR(sphere(Point3D(0, 0, 0), 0.0), 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// box +// ═══════════════════════════════════════════════════ + +TEST(SdfBox, Outside) { + Point3D b(1, 1, 1); + // Outside corner + double d = box(Point3D(2, 2, 2), b); + EXPECT_GT(d, 0.0); + EXPECT_NEAR(d, std::sqrt(3.0), EPS); +} + +TEST(SdfBox, OnSurface) { + Point3D b(1, 1, 1); + EXPECT_NEAR(box(Point3D(1, 0, 0), b), 0.0, EPS); + EXPECT_NEAR(box(Point3D(0, 1, 0), b), 0.0, EPS); + EXPECT_NEAR(box(Point3D(0, 0, 1), b), 0.0, EPS); + EXPECT_NEAR(box(Point3D(-1, 0, 0), b), 0.0, EPS); +} + +TEST(SdfBox, Inside) { + Point3D b(3, 3, 3); + // Center + EXPECT_NEAR(box(Point3D(0, 0, 0), b), -3.0, EPS); + // Mid-edge + double d = box(Point3D(2, 0, 0), b); + EXPECT_NEAR(d, -1.0, EPS); +} + +TEST(SdfBox, Degenerate) { + Point3D zero(0, 0, 0); + double d = box(Point3D(1, 1, 1), zero); + EXPECT_NEAR(d, std::sqrt(3.0), EPS); +} + +// ═══════════════════════════════════════════════════ +// round_box +// ═══════════════════════════════════════════════════ + +TEST(SdfRoundBox, MatchesBoxAtZeroRadius) { + Point3D b(2, 3, 4); + for (double x = -5; x <= 5; x += 2.5) { + for (double y = -5; y <= 5; y += 2.5) { + for (double z = -5; z <= 5; z += 2.5) { + Point3D pt(x, y, z); + EXPECT_NEAR(round_box(pt, b, 0.0), box(pt, b), EPS); + } + } + } +} + +TEST(SdfRoundBox, RoundedShiftsBoundary) { + Point3D b(1, 1, 1); + // On surface of round box with r=0.5: outside moves in by 0.5 + double d_box = box(Point3D(1.5, 0, 0), b); + double d_rbox = round_box(Point3D(1.5, 0, 0), b, 0.5); + EXPECT_NEAR(d_rbox, d_box - 0.5, EPS); +} + +// ═══════════════════════════════════════════════════ +// torus +// ═══════════════════════════════════════════════════ + +TEST(SdfTorus, OnSurfaceMajor) { + // Point on major circle + double d = torus(Point3D(5, 0, 0), 5.0, 1.0); + EXPECT_NEAR(d, -1.0, EPS); // inside tube, at center of tube cross-section +} + +TEST(SdfTorus, RingCenter_Inside) { + // Center of torus (inside the hole) + double d = torus(Point3D(0, 0, 0), 5.0, 1.0); + EXPECT_NEAR(d, 4.0, EPS); // distance from origin to tube center = 5, minus tube radius = 4 +} + +TEST(SdfTorus, OutsideTube) { + // Outside the whole shape + double d = torus(Point3D(0, 3, 0), 5.0, 1.0); + // Point (0,3,0): distance from Y axis = 0, so tube center at (5,0,0) + // Distance to tube center: sqrt(25 + 9) = sqrt(34) ≈ 5.83, minus 1 = 4.83 + EXPECT_NEAR(d, std::sqrt(34.0) - 1.0, EPS); +} + +TEST(SdfTorus, Degenerate) { + double d = torus(Point3D(3, 0, 0), 3.0, 0.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// capsule +// ═══════════════════════════════════════════════════ + +TEST(SdfCapsule, OnAxis) { + Point3D a(0, 0, 0), b(0, 4, 0); + double r = 1.0; + EXPECT_NEAR(capsule(Point3D(0, 0, 0), a, b, r), -r, EPS); + EXPECT_NEAR(capsule(Point3D(0, 0, 1.0), a, b, r), 0.0, EPS); + EXPECT_NEAR(capsule(Point3D(0, 0, 0), a, b, 0.0), 0.0, EPS); +} + +TEST(SdfCapsule, OutsideCylinder) { + Point3D a(0, 0, 0), b(0, 4, 0); + double r = 1.0; + double d = capsule(Point3D(3, 2, 0), a, b, r); + EXPECT_NEAR(d, 2.0, EPS); +} + +TEST(SdfCapsule, Midpoint) { + Point3D a(0, 0, 0), b(2, 0, 0); + double r = 0.5; + EXPECT_NEAR(capsule(Point3D(1, 0, 0), a, b, r), -0.5, EPS); + EXPECT_NEAR(capsule(Point3D(1, 0.5, 0), a, b, r), 0.0, EPS); +} + +TEST(SdfCapsule, Degenerate) { + Point3D pt(0, 0, 0); + // Zero-length capsule = sphere + double d = capsule(Point3D(2, 0, 0), pt, pt, 3.0); + EXPECT_NEAR(d, -1.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// cylinder +// ═══════════════════════════════════════════════════ + +TEST(SdfCylinder, OnSurface) { + double d = cylinder(Point3D(2, 0, 0), 2.0, 4.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +TEST(SdfCylinder, Inside) { + double d = cylinder(Point3D(0, 0, 0), 2.0, 4.0); + EXPECT_NEAR(d, -2.0, EPS); +} + +TEST(SdfCylinder, OutsideCap) { + // Above top cap + double d = cylinder(Point3D(0, 3, 0), 1.0, 4.0); + EXPECT_NEAR(d, 1.0, EPS); +} + +TEST(SdfCylinder, OutsideSide) { + // Outside the side + double d = cylinder(Point3D(3, 0, 0), 1.0, 4.0); + EXPECT_NEAR(d, 2.0, EPS); +} + +TEST(SdfCylinder, Degenerate) { + // Zero radius = line segment + double d = cylinder(Point3D(0, 0, 0), 0.0, 4.0); + EXPECT_NEAR(d, 0.0, EPS); + // Zero height + double d2 = cylinder(Point3D(0, 0, 0), 2.0, 0.0); + EXPECT_NEAR(d2, 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// cone +// ═══════════════════════════════════════════════════ + +TEST(SdfCone, Apex) { + // Apex at y = height/2 = 2 + double d = cone(Point3D(0, 2, 0), std::atan(0.5), 4.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +TEST(SdfCone, OnBaseEdge) { + // Base at y = -height/2 = -2, base_radius = 4 * tan(atan(0.5)) = 2 + double d = cone(Point3D(2, -2, 0), std::atan(0.5), 4.0); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +TEST(SdfCone, Inside) { + // Inside the cone body + double d = cone(Point3D(0, 0, 0), std::atan(0.5), 4.0); + EXPECT_LT(d, 0.0); +} + +TEST(SdfCone, Outside) { + // Outside near base + double d = cone(Point3D(3, -2, 0), std::atan(0.5), 4.0); + EXPECT_GT(d, 0.0); +} + +TEST(SdfCone, BelowBase) { + // Below the base plane + double d = cone(Point3D(0, -3, 0), std::atan(0.5), 4.0); + EXPECT_GT(d, 0.0); +} + +// ═══════════════════════════════════════════════════ +// plane +// ═══════════════════════════════════════════════════ + +TEST(SdfPlane, Above) { + Vector3D n(0, 1, 0); + EXPECT_NEAR(plane(Point3D(0, 5, 0), n, 0.0), 5.0, EPS); +} + +TEST(SdfPlane, Below) { + Vector3D n(0, 1, 0); + EXPECT_NEAR(plane(Point3D(0, -3, 0), n, 0.0), -3.0, EPS); +} + +TEST(SdfPlane, WithOffset) { + Vector3D n(0, 1, 0); + EXPECT_NEAR(plane(Point3D(0, 5, 0), n, 2.0), 3.0, EPS); + EXPECT_NEAR(plane(Point3D(0, 1, 0), n, 2.0), -1.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// ellipsoid +// ═══════════════════════════════════════════════════ + +TEST(SdfEllipsoid, SphereCase) { + // Ellipsoid with equal radii = sphere + Point3D radii(3, 3, 3); + double d = ellipsoid(Point3D(3, 0, 0), radii); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +TEST(SdfEllipsoid, Inside) { + Point3D radii(3, 2, 1); + double d = ellipsoid(Point3D(0, 0, 0), radii); + EXPECT_NEAR(d, -1.0, 1e-6); +} + +TEST(SdfEllipsoid, OnSurface) { + // At (3,0,0): scaled = (1,0,0), |scaled|=1 → on surface + Point3D radii(3, 2, 1); + double d = ellipsoid(Point3D(3, 0, 0), radii); + EXPECT_NEAR(d, 0.0, 1e-6); + d = ellipsoid(Point3D(0, 2, 0), radii); + EXPECT_NEAR(d, 0.0, 1e-6); + d = ellipsoid(Point3D(0, 0, 1), radii); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +// ═══════════════════════════════════════════════════ +// triangular_prism +// ═══════════════════════════════════════════════════ + +TEST(SdfTriangularPrism, OnVertex) { + Point3D a(0, 0, 0), b(4, 0, 0), c(2, 3, 0); + double d = triangular_prism(Point3D(0, 0, 0), a, b, c, 2.0); + EXPECT_NEAR(d, -1.0, 1e-6); +} + +TEST(SdfTriangularPrism, Inside) { + Point3D a(0, 0, 0), b(4, 0, 0), c(2, 3, 0); + // Centroid of triangle, mid-Z + double d = triangular_prism(Point3D(2, 1, 0), a, b, c, 4.0); + EXPECT_LT(d, 0.0); +} + +TEST(SdfTriangularPrism, OutsideXY) { + Point3D a(0, 0, 0), b(4, 0, 0), c(2, 3, 0); + // Point far outside the triangle but within Z range + double d = triangular_prism(Point3D(10, 10, 0), a, b, c, 2.0); + EXPECT_GT(d, 0.0); +} + +TEST(SdfTriangularPrism, OutsideZ) { + Point3D a(0, 0, 0), b(4, 0, 0), c(2, 3, 0); + // Inside triangle but far above top cap + double d = triangular_prism(Point3D(2, 1, 5), a, b, c, 2.0); + EXPECT_NEAR(d, 4.0, 1e-6); +} + +// ═══════════════════════════════════════════════════ +// hex_prism +// ═══════════════════════════════════════════════════ + +TEST(SdfHexPrism, Center) { + double d = hex_prism(Point3D(0, 0, 0), 2.0, 4.0); + EXPECT_LT(d, 0.0); +} + +TEST(SdfHexPrism, OnVertex) { + // Vertex at distance r along X + double d = hex_prism(Point3D(2, 0, 0), 2.0, 2.0); + EXPECT_NEAR(d, 0.0, 2e-5); +} + +TEST(SdfHexPrism, OutsideCap) { + double d = hex_prism(Point3D(0, 3, 0), 1.0, 4.0); + EXPECT_NEAR(d, 1.0, 1e-6); +} + +TEST(SdfHexPrism, Degenerate) { + double d = hex_prism(Point3D(0, 0, 0), 0.0, 0.0); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +// ═══════════════════════════════════════════════════ +// link +// ═══════════════════════════════════════════════════ + +TEST(SdfLink, AtCenterOfFirstTorus) { + // Length 4, so torus centers at x=±2 + // Center of first torus tube: x=2+3=5? No, torus major_r=3 + // Tube center of first torus at (2+3, 0, 0) = (5, 0, 0) in XZ plane + double d = link(Point3D(5, 0, 0), 4.0, 3.0, 1.0); + EXPECT_NEAR(d, -1.0, EPS); +} + +TEST(SdfLink, BetweenToruses) { + // Point between the two toruses, inside the overlap region + double d = link(Point3D(0, 0, 0), 2.0, 2.0, 0.5); + // This is inside both toruses, so negative + EXPECT_LT(d, 0.0); +} + +TEST(SdfLink, FarOutside) { + double d = link(Point3D(100, 0, 0), 4.0, 3.0, 1.0); + EXPECT_GT(d, 90.0); +} + +// ═══════════════════════════════════════════════════ +// infinite_cylinder +// ═══════════════════════════════════════════════════ + +TEST(SdfInfiniteCylinder, OnSurface) { + Vector3D axis(0, 1, 0); + double d = infinite_cylinder(Point3D(2, 5, 0), axis, 2.0); + EXPECT_NEAR(d, 0.0, EPS); + // Any Y should work (infinite) + d = infinite_cylinder(Point3D(2, 100, 0), axis, 2.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +TEST(SdfInfiniteCylinder, Inside) { + Vector3D axis(0, 1, 0); + double d = infinite_cylinder(Point3D(0, 0, 0), axis, 3.0); + EXPECT_NEAR(d, -3.0, EPS); +} + +TEST(SdfInfiniteCylinder, Outside) { + Vector3D axis(0, 1, 0); + double d = infinite_cylinder(Point3D(5, 0, 0), axis, 1.0); + EXPECT_NEAR(d, 4.0, EPS); +} + +TEST(SdfInfiniteCylinder, SkewAxis) { + Vector3D axis(1, 0, 0); + // Along X axis: distance = sqrt(y²+z²) - r + double d = infinite_cylinder(Point3D(0, 3, 4), axis, 5.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// infinite_cone +// ═══════════════════════════════════════════════════ + +TEST(SdfInfiniteCone, AtApex) { + Point3D apex(0, 0, 0); + Vector3D axis(0, 1, 0); + double d = infinite_cone(apex, apex, axis, 0.5); + EXPECT_NEAR(d, 0.0, EPS); +} + +TEST(SdfInfiniteCone, OnSurface) { + Point3D apex(0, 0, 0); + Vector3D axis(0, 1, 0); + double angle = std::atan(1.0); // 45° cone + // At y=2, radius should be 2 (tan45° = 1) + double d = infinite_cone(Point3D(2, 2, 0), apex, axis, angle); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +TEST(SdfInfiniteCone, Inside) { + Point3D apex(0, 0, 0); + Vector3D axis(0, 1, 0); + double angle = std::atan(2.0); // wide cone + // Point close to axis, should be inside + double d = infinite_cone(Point3D(1, 2, 0), apex, axis, angle); + EXPECT_LT(d, 0.0); +} + +TEST(SdfInfiniteCone, Outside) { + Point3D apex(0, 0, 0); + Vector3D axis(0, 1, 0); + double angle = std::atan(0.5); // narrow cone + // Far from axis + double d = infinite_cone(Point3D(10, 2, 0), apex, axis, angle); + EXPECT_GT(d, 0.0); +} + +TEST(SdfInfiniteCone, BehindApex) { + Point3D apex(0, 0, 0); + Vector3D axis(0, 1, 0); + double angle = std::atan(1.0); + // Behind apex (negative Y) → outside + double d = infinite_cone(Point3D(0, -1, 0), apex, axis, angle); + EXPECT_GT(d, 0.0); +} + +// ═══════════════════════════════════════════════════ +// wedge +// ═══════════════════════════════════════════════════ + +TEST(SdfWedge, CenterOfBoxHalf) { + // Wedge: box [-w/2,w/2]×[-h/2,h/2]×[-d/2,d/2] ∩ {z ≥ x} + // Point inside the wedge region + double d = wedge(Point3D(-0.5, 0.0, 1.0), 4.0, 4.0, 4.0); + EXPECT_LT(d, 0.0); +} + +TEST(SdfWedge, OutsideDiagonal) { + // Point in box but on wrong side of diagonal (x > z) + double d = wedge(Point3D(1.5, 0.0, 0.0), 4.0, 4.0, 4.0); + EXPECT_GT(d, 0.0); +} + +TEST(SdfWedge, OnDiagonalPlane) { + // On the diagonal plane z = x + double d = wedge(Point3D(1, 0, 1), 4.0, 4.0, 4.0); + EXPECT_NEAR(d, 0.0, 1e-6); +} + +TEST(SdfWedge, OutsideBox) { + double d = wedge(Point3D(0, 10, 10), 4.0, 4.0, 4.0); + EXPECT_GT(d, 0.0); +} + +// ═══════════════════════════════════════════════════ +// extrusion +// ═══════════════════════════════════════════════════ + +TEST(SdfExtrusion, PassThrough) { + double sd2 = 3.5; + EXPECT_NEAR(extrusion(Point3D(1, 2, 100), sd2), 3.5, EPS); + EXPECT_NEAR(extrusion(Point3D(-5, 7, -3), -2.0), -2.0, EPS); +} + +TEST(SdfExtrusion, Zero) { + EXPECT_NEAR(extrusion(Point3D(0, 0, 0), 0.0), 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// extrusion_bounded +// ═══════════════════════════════════════════════════ + +TEST(SdfExtrusionBounded, InsideSlabAndShape) { + // 2D SDF says outside (positive), Z within bounds → d = max(pos, inside_Z) + double d = extrusion_bounded(Point3D(0, 0, 0), -1.0, 2.0); + EXPECT_NEAR(d, -1.0, EPS); +} + +TEST(SdfExtrusionBounded, OutsideSlab) { + // 2D SDF inside, but Z outside slab + double d = extrusion_bounded(Point3D(0, 0, 5), -1.0, 2.0); + EXPECT_NEAR(d, 3.0, EPS); +} + +TEST(SdfExtrusionBounded, Zero) { + double d = extrusion_bounded(Point3D(0, 0, 0), 0.0, 0.0); + EXPECT_NEAR(d, 0.0, EPS); +} + +// ═══════════════════════════════════════════════════ +// revolution +// ═══════════════════════════════════════════════════ + +TEST(SdfRevolution, PassThrough) { + double sd2 = -1.5; + EXPECT_NEAR(revolution(Point3D(3, 2, 4), sd2, 1.0), -1.5, EPS); +} + +TEST(SdfRevolution, Zero) { + EXPECT_NEAR(revolution(Point3D(0, 0, 0), 0.0, 0.0), 0.0, EPS); +} diff --git a/tests/sdf/test_sdf_to_mesh.cpp b/tests/sdf/test_sdf_to_mesh.cpp new file mode 100644 index 0000000..bad54a5 --- /dev/null +++ b/tests/sdf/test_sdf_to_mesh.cpp @@ -0,0 +1,86 @@ +#include "vde/sdf/sdf_to_mesh.h" +#include + +using namespace vde::sdf; +using core::Point3D; + +TEST(SdfToMesh, SphereProducesValidMesh) { + auto sphere = SdfNode::sphere(2.0); + auto mesh = sdf_to_mesh(sphere, 32, 0.0); + + EXPECT_GT(mesh.vertices.size(), 0u); + EXPECT_GT(mesh.triangles.size(), 0u); + + // All vertices should be near the sphere surface (r ≈ 2) + for (const auto& v : mesh.vertices) { + double r = v.norm(); + EXPECT_NEAR(r, 2.0, 1.0); // coarse resolution → loose tolerance + } + + // Each triangle should have 3 indices + for (const auto& tri : mesh.triangles) { + for (int i = 0; i < 3; ++i) { + EXPECT_GE(tri[i], 0); + EXPECT_LT(tri[i], static_cast(mesh.vertices.size())); + } + } +} + +TEST(SdfToMesh, UnionOfTwoSpheres) { + auto s1 = SdfNode::translate(SdfNode::sphere(1.0), Point3D(-1.5, 0, 0)); + auto s2 = SdfNode::translate(SdfNode::sphere(1.0), Point3D(1.5, 0, 0)); + auto tree = SdfNode::op_union(std::move(s1), std::move(s2)); + + auto mesh = sdf_to_mesh(tree, 32, 0.0); + + EXPECT_GT(mesh.vertices.size(), 0u); + EXPECT_GT(mesh.triangles.size(), 0u); +} + +TEST(SdfToMesh, BoxProducesValidMesh) { + auto box = SdfNode::box(Point3D(1, 1, 1)); + auto mesh = sdf_to_mesh(box, 32, 0.0); + + EXPECT_GT(mesh.vertices.size(), 0u); + EXPECT_GT(mesh.triangles.size(), 0u); +} + +TEST(SdfToMesh, HigherResolutionProducesMoreTriangles) { + auto sphere = SdfNode::sphere(1.0); + auto low = sdf_to_mesh(sphere, 16, 0.0); + auto high = sdf_to_mesh(sphere, 32, 0.0); + + EXPECT_GT(high.vertices.size(), low.vertices.size()); +} + +TEST(SdfToMesh, LambdaSphere) { + auto mesh = sdf_to_mesh_lambda( + [](double x, double y, double z) { + return std::sqrt(x*x + y*y + z*z) - 2.0; + }, + Point3D(-3, -3, -3), Point3D(3, 3, 3), + 32, 0.0 + ); + + EXPECT_GT(mesh.vertices.size(), 0u); + EXPECT_GT(mesh.triangles.size(), 0u); +} + +TEST(SdfToMesh, LambdaBox) { + auto mesh = sdf_to_mesh_lambda( + [](double x, double y, double z) { + double dx = std::abs(x) - 1.0; + double dy = std::abs(y) - 1.0; + double dz = std::abs(z) - 1.0; + return std::sqrt(std::max(dx, 0.0) * std::max(dx, 0.0) + + std::max(dy, 0.0) * std::max(dy, 0.0) + + std::max(dz, 0.0) * std::max(dz, 0.0)) + + std::min(std::max({dx, dy, dz}), 0.0); + }, + Point3D(-2, -2, -2), Point3D(2, 2, 2), + 32, 0.0 + ); + + EXPECT_GT(mesh.vertices.size(), 0u); + EXPECT_GT(mesh.triangles.size(), 0u); +} diff --git a/tests/sdf/test_sdf_tree.cpp b/tests/sdf/test_sdf_tree.cpp new file mode 100644 index 0000000..c8f422c --- /dev/null +++ b/tests/sdf/test_sdf_tree.cpp @@ -0,0 +1,285 @@ +#include "vde/sdf/sdf_tree.h" +#include +#include + +using namespace vde::sdf; +using core::Point3D; +using core::Vector3D; + +// ── Null/empty tree ────────────────────────────────────────────── + +TEST(SdfTree, NullNodeReturnsLargeValue) { + SdfNodePtr null_node = nullptr; + double d = evaluate(null_node, Point3D(0, 0, 0)); + EXPECT_GT(d, 1e10); +} + +// ── Primitive evaluation ───────────────────────────────────────── + +TEST(SdfTree, SphereInsideOutsideBoundary) { + auto s = SdfNode::sphere(1.0); + + // Origin = inside + double d_in = evaluate(s, Point3D(0, 0, 0)); + EXPECT_LT(d_in, 0.0); + EXPECT_NEAR(d_in, -1.0, 1e-6); + + // Surface (r=1 on X axis) + double d_surf = evaluate(s, Point3D(1, 0, 0)); + EXPECT_NEAR(d_surf, 0.0, 1e-6); + + // Outside + double d_out = evaluate(s, Point3D(2, 0, 0)); + EXPECT_GT(d_out, 0.0); + EXPECT_NEAR(d_out, 1.0, 1e-6); +} + +TEST(SdfTree, BoxEvaluation) { + auto b = SdfNode::box(Point3D(1, 1, 1)); + + // Center = inside + EXPECT_LT(evaluate(b, Point3D(0, 0, 0)), 0.0); + // Surface + EXPECT_NEAR(evaluate(b, Point3D(1, 0, 0)), 0.0, 1e-6); + // Outside + EXPECT_GT(evaluate(b, Point3D(2, 0, 0)), 0.0); +} + +TEST(SdfTree, TorusEvaluation) { + auto t = SdfNode::torus(1.0, 0.3); + + // Point on the ring centerline (should be -minor_radius = inside) + double d = evaluate(t, Point3D(1, 0, 0)); + EXPECT_LT(d, 0.0); + EXPECT_NEAR(d, -0.3, 1e-6); +} + +TEST(SdfTree, CylinderEvaluation) { + auto c = SdfNode::cylinder(0.5, 1.0); + + EXPECT_LT(evaluate(c, Point3D(0, 0, 0)), 0.0); // center + EXPECT_NEAR(evaluate(c, Point3D(0.5, 0, 0)), 0.0, 1e-6); // surface + EXPECT_GT(evaluate(c, Point3D(1, 0, 0)), 0.0); // outside radius + EXPECT_GT(evaluate(c, Point3D(0, 2, 0)), 0.0); // above cap +} + +TEST(SdfTree, PlaneEvaluation) { + auto p = SdfNode::plane(Vector3D(0, 1, 0), 0.0); // horizontal at y=0 + + EXPECT_GT(evaluate(p, Point3D(0, 1, 0)), 0.0); // above + EXPECT_LT(evaluate(p, Point3D(0, -1, 0)), 0.0); // below + EXPECT_NEAR(evaluate(p, Point3D(0, 0, 0)), 0.0, 1e-6); +} + +// ── CSG boolean operations ─────────────────────────────────────── + +TEST(SdfTree, UnionSphereBox) { + auto tree = SdfNode::op_union( + SdfNode::sphere(1.0), + SdfNode::box(Point3D(1, 1, 1)) + ); + + // Point inside sphere + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // Point inside box (but outside sphere) + EXPECT_LT(evaluate(tree, Point3D(1.2, 0, 0)), 0.0); + // Point outside both + EXPECT_GT(evaluate(tree, Point3D(3, 0, 0)), 0.0); +} + +TEST(SdfTree, IntersectionSphereBox) { + auto tree = SdfNode::op_intersection( + SdfNode::sphere(1.0), + SdfNode::box(Point3D(0.5, 0.5, 0.5)) + ); + + // Inside intersection region + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // Inside sphere but outside box + EXPECT_GT(evaluate(tree, Point3D(0.8, 0, 0)), 0.0); +} + +TEST(SdfTree, DifferenceSphereBox) { + auto tree = SdfNode::op_difference( + SdfNode::sphere(2.0), + SdfNode::box(Point3D(1, 1, 1)) + ); + + // Inside sphere, outside box → inside + EXPECT_LT(evaluate(tree, Point3D(1.5, 0, 0)), 0.0); + // Inside box (removed region) → outside + EXPECT_GT(evaluate(tree, Point3D(0, 0, 0)), 0.0); +} + +// ── Smooth boolean operations ──────────────────────────────────── + +TEST(SdfTree, SmoothUnion) { + auto tree = SdfNode::smooth_union( + SdfNode::sphere(1.0), + SdfNode::sphere(1.0), + 0.3 + ); + + // Both spheres at origin → should be equivalent + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); +} + +// ── Domain transforms ──────────────────────────────────────────── + +TEST(SdfTree, TranslatedSphere) { + auto tree = SdfNode::translate( + SdfNode::sphere(1.0), + Point3D(2, 0, 0) + ); + + // Original origin should be far away + EXPECT_GT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // New center should be inside + EXPECT_LT(evaluate(tree, Point3D(2, 0, 0)), 0.0); + EXPECT_NEAR(evaluate(tree, Point3D(2, 0, 0)), -1.0, 1e-6); + // Surface at new center + radius + EXPECT_NEAR(evaluate(tree, Point3D(3, 0, 0)), 0.0, 1e-6); +} + +TEST(SdfTree, ScaledSphere) { + auto tree = SdfNode::scale( + SdfNode::sphere(1.0), + Point3D(2, 1, 1) + ); + + // At (2,0,0) in world = (1,0,0) in local → on surface + EXPECT_NEAR(evaluate(tree, Point3D(2, 0, 0)), 0.0, 1e-6); + // At (0,0,0) in world = inside + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); +} + +TEST(SdfTree, RepeatedSphere) { + auto tree = SdfNode::repeat( + SdfNode::sphere(0.3), + Point3D(2, 2, 2) + ); + + // Origin → mapped to (0,0,0) inside + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // (2,0,0) → mapped to (0,0,0) inside + EXPECT_LT(evaluate(tree, Point3D(2, 0, 0)), 0.0); + // Between copies → outside + EXPECT_GT(evaluate(tree, Point3D(1, 1, 0)), 0.0); +} + +TEST(SdfTree, MirroredX) { + auto tree = SdfNode::mirror_x( + SdfNode::sphere(1.0) + ); + + // Both sides should have the sphere + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + EXPECT_LT(evaluate(tree, Point3D(-0.5, 0, 0)), 0.0); + // Surface at ±1 + EXPECT_NEAR(evaluate(tree, Point3D(1, 0, 0)), 0.0, 1e-6); + EXPECT_NEAR(evaluate(tree, Point3D(-1, 0, 0)), 0.0, 1e-6); +} + +// ── Nested operations ──────────────────────────────────────────── + +TEST(SdfTree, NestedUnionTwist) { + auto inner = SdfNode::op_union( + SdfNode::sphere(1.0), + SdfNode::box(Point3D(0.6, 0.6, 0.6)) + ); + auto tree = SdfNode::twist(std::move(inner), 0.5); + + // Should still evaluate successfully + double d = evaluate(tree, Point3D(0, 0.5, 0)); + EXPECT_LT(d, 0.0); // interior of union +} + +TEST(SdfTree, DeepNesting) { + auto a = SdfNode::sphere(1.0); + auto b = SdfNode::box(Point3D(0.5, 0.5, 0.5)); + auto uni = SdfNode::op_union(std::move(a), std::move(b)); + auto scaled = SdfNode::scale(std::move(uni), Point3D(2, 2, 2)); + auto translated = SdfNode::translate(std::move(scaled), Point3D(3, 0, 0)); + + // Original origin → should be inside after inverse transform + union + double d = evaluate(translated, Point3D(3, 0, 0)); + EXPECT_LT(d, 0.0); +} + +// ── Modifiers ──────────────────────────────────────────────────── + +TEST(SdfTree, RoundModifier) { + auto tree = SdfNode::round( + SdfNode::box(Point3D(1, 1, 1)), + 0.2 + ); + + // Center still inside + EXPECT_LT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // Original corner point (1,1,0) should now be outside (rounded away) + double d_corner = evaluate(tree, Point3D(1, 0, 0)); + // Rounded box should have slightly larger radius at corners + // The exact value depends on implementation, just verify it's farther + double d_box = evaluate(SdfNode::box(Point3D(1, 1, 1)), Point3D(1, 0, 0)); + EXPECT_NEAR(d_corner, d_box - 0.2, 1e-6); +} + +TEST(SdfTree, OnionModifier) { + auto tree = SdfNode::onion( + SdfNode::sphere(2.0), + 0.3 + ); + + // Center should be outside (hollowed out) + EXPECT_GT(evaluate(tree, Point3D(0, 0, 0)), 0.0); + // Shell surface should exist + EXPECT_NEAR(evaluate(tree, Point3D(1.85, 0, 0)), 0.0, 1e-6); +} + +// ── Bounds estimation ──────────────────────────────────────────── + +TEST(SdfTree, EstimateBoundsSphere) { + auto s = SdfNode::sphere(2.0); + Point3D b = estimate_bounds(s, 0.0); + EXPECT_NEAR(b.x(), 2.0, 1e-6); + EXPECT_NEAR(b.y(), 2.0, 1e-6); + EXPECT_NEAR(b.z(), 2.0, 1e-6); +} + +TEST(SdfTree, EstimateBoundsWithMargin) { + auto s = SdfNode::sphere(2.0); + Point3D b = estimate_bounds(s, 1.0); + EXPECT_NEAR(b.x(), 3.0, 1e-6); +} + +TEST(SdfTree, EstimateBoundsUnion) { + auto tree = SdfNode::op_union( + SdfNode::sphere(2.0), + SdfNode::box(Point3D(3, 3, 3)) + ); + Point3D b = estimate_bounds(tree, 0.0); + // Box half-extent norm = sqrt(27) ≈ 5.196 + EXPECT_GT(b.norm(), 4.0); +} + +TEST(SdfTree, EstimateBoundsTranslated) { + auto tree = SdfNode::translate( + SdfNode::sphere(1.0), + Point3D(5, 0, 0) + ); + Point3D b = estimate_bounds(tree, 0.0); + EXPECT_GT(b.norm(), 5.0); +} + +// ── Tree traversal ─────────────────────────────────────────────── + +TEST(SdfTree, VisitCountsNodes) { + auto tree = SdfNode::op_union( + SdfNode::sphere(1.0), + SdfNode::box(Point3D(1, 1, 1)) + ); + + int count = 0; + tree->visit([&count](const SdfNode&) { ++count; }); + EXPECT_EQ(count, 3); // union + sphere + box +}