茂之钳
692 lines
24 KiB
C++
692 lines
24 KiB
C++
#pragma once
|
||
/**
|
||
* @file sdf_gradient.h
|
||
* @brief SDF 梯度计算（数值微分、解析梯度、链式法则）
|
||
*
|
||
* 提供两类梯度计算：
|
||
*
|
||
* 1. **空间梯度** ∇f(p) — SDF 在场中的方向导数（法向量方向）
|
||
* - 有限差分近似
|
||
* - 各图元的解析梯度
|
||
* - CSG 操作的链式法则
|
||
*
|
||
* 2. **参数梯度** ∂f/∂θ — SDF 对形状参数的敏感度
|
||
* - 有限差分计算 ∂f/∂radius, ∂f/∂height 等
|
||
* - 用于基于梯度的形状优化
|
||
*
|
||
* ## 梯度约定
|
||
*
|
||
* 对于 SDF，梯度指向距离增加最快的方向。
|
||
* 在表面上 (f=0)，∇f 等于外法向量（单位长度）。
|
||
* 在外部 (f>0)，∇f 指向远离表面的方向。
|
||
* 在内部 (f<0)，∇f 指向最近的表面。
|
||
*
|
||
* @ingroup sdf
|
||
*/
|
||
#include "vde/sdf/sdf_primitives.h"
|
||
#include "vde/sdf/sdf_operations.h"
|
||
#include "vde/core/point.h"
|
||
#include <cmath>
|
||
#include <functional>
|
||
|
||
namespace vde::sdf {
|
||
|
||
using core::Point3D;
|
||
using core::Vector3D;
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Finite-Difference Gradient
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief SDF 的空间梯度（中心差分近似）
|
||
*
|
||
* 使用三轴中心有限差分:
|
||
* ∂f/∂x ≈ (f(x+h, y, z) - f(x-h, y, z)) / (2h)
|
||
*
|
||
* 二阶精度，适用于任意黑盒 SDF 函数。
|
||
*
|
||
* @param f SDF 函数（接受 Point3D 参数）
|
||
* @param p 求值点
|
||
* @param h 差分步长（默认 1e-6）
|
||
* @return 梯度向量 ∇f(p)
|
||
*
|
||
* @note 步长太小会导致舍入误差主导，太大导致截断误差主导。
|
||
* 默认值 1e-6 在单精度场景下是较好的平衡。
|
||
* @note 对于含不连续性的 SDF（如 sharp CSG），有限差分可能在边界产生错误梯度。
|
||
*
|
||
* @code{.cpp}
|
||
* auto f = [](const Point3D& p) { return sphere(p, 1.0); };
|
||
* Vector3D n = gradient(f, Point3D(1, 0, 0)); // ≈ (1, 0, 0)
|
||
* @endcode
|
||
*
|
||
* @see evaluate_with_gradient 一次性获取值和梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D gradient(
|
||
const std::function<double(const Point3D&)>& f,
|
||
const Point3D& p, double h = 1e-6)
|
||
{
|
||
Vector3D g;
|
||
double inv_2h = 0.5 / h;
|
||
g.x() = (f(Point3D(p.x() + h, p.y(), p.z())) - f(Point3D(p.x() - h, p.y(), p.z()))) * inv_2h;
|
||
g.y() = (f(Point3D(p.x(), p.y() + h, p.z())) - f(Point3D(p.x(), p.y() - h, p.z()))) * inv_2h;
|
||
g.z() = (f(Point3D(p.x(), p.y(), p.z() + h)) - f(Point3D(p.x(), p.y(), p.z() - h))) * inv_2h;
|
||
return g;
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Gradient Result Struct
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 前向模梯度结果: SDF 值 + 梯度
|
||
*
|
||
* 在一次计算中返回 SDF 值和空间梯度，避免重复求值。
|
||
* 梯度向量在表面上即为单位外法向量的近似。
|
||
*/
|
||
struct GradResult {
|
||
double value; ///< SDF 值 f(p)
|
||
Vector3D grad; ///< 空间梯度 ∇f(p) = (df/dx, df/dy, df/dz)
|
||
};
|
||
|
||
/**
|
||
* @brief 同时计算 SDF 值和梯度的便捷函数
|
||
*
|
||
* @param f SDF 函数
|
||
* @param p 求值点
|
||
* @param h 差分步长
|
||
* @return value + grad 结构
|
||
*
|
||
* @see gradient 仅计算梯度
|
||
*/
|
||
[[nodiscard]] inline GradResult evaluate_with_gradient(
|
||
const std::function<double(const Point3D&)>& f,
|
||
const Point3D& p, double h = 1e-6)
|
||
{
|
||
return GradResult{f(p), gradient(f, p, h)};
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Parameter Gradients (Finite Difference)
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 球体 SDF 对半径的参数梯度 ∂f/∂r
|
||
*
|
||
* 用中心差分计算: (sphere(p, r+h) - sphere(p, r-h)) / (2h)
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 球体半径
|
||
* @param h 差分步长（默认 1e-6）
|
||
* @return ∂f/∂radius
|
||
*
|
||
* @note 对于球体，解析梯度为 -1（球面上），但此函数提供与实现无关的数值验证。
|
||
*
|
||
* @see sphere_param_grad 封装版本
|
||
*/
|
||
[[nodiscard]] inline double sphere_radius_gradient(const Point3D& p, double radius, double h = 1e-6)
|
||
{
|
||
return (sphere(p, radius + h) - sphere(p, radius - h)) / (2.0 * h);
|
||
}
|
||
|
||
/**
|
||
* @brief 长方体 SDF 对各轴半尺寸的参数梯度 ∂f/∂extent[i]
|
||
*
|
||
* @param p 采样点
|
||
* @param extents 半边长
|
||
* @param axis 轴索引: 0=x, 1=y, 2=z
|
||
* @param h 差分步长
|
||
* @return ∂f/∂extent[axis]
|
||
*
|
||
* @see box_param_grad 封装版本
|
||
*/
|
||
[[nodiscard]] inline double box_extent_gradient(const Point3D& p, const Point3D& extents,
|
||
int axis, double h = 1e-6)
|
||
{
|
||
Point3D ext_plus = extents;
|
||
Point3D ext_minus = extents;
|
||
ext_plus[axis] += h;
|
||
ext_minus[axis] -= h;
|
||
return (box(p, ext_plus) - box(p, ext_minus)) / (2.0 * h);
|
||
}
|
||
|
||
/**
|
||
* @brief 圆环 SDF 对主半径的参数梯度 ∂f/∂major_r
|
||
*
|
||
* @param p 采样点
|
||
* @param major_r 主半径
|
||
* @param minor_r 副半径
|
||
* @param h 差分步长
|
||
* @return 参数梯度
|
||
*/
|
||
[[nodiscard]] inline double torus_major_gradient(const Point3D& p, double major_r, double minor_r,
|
||
double h = 1e-6)
|
||
{
|
||
return (torus(p, major_r + h, minor_r) - torus(p, major_r - h, minor_r)) / (2.0 * h);
|
||
}
|
||
|
||
/**
|
||
* @brief 圆环 SDF 对副半径的参数梯度 ∂f/∂minor_r
|
||
*
|
||
* @return 参数梯度
|
||
*/
|
||
[[nodiscard]] inline double torus_minor_gradient(const Point3D& p, double major_r, double minor_r,
|
||
double h = 1e-6)
|
||
{
|
||
return (torus(p, major_r, minor_r + h) - torus(p, major_r, minor_r - h)) / (2.0 * h);
|
||
}
|
||
|
||
/**
|
||
* @brief 圆柱 SDF 对高度的参数梯度 ∂f/∂height
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 圆柱半径
|
||
* @param height 圆柱高度
|
||
* @param h 差分步长
|
||
* @return 参数梯度
|
||
*/
|
||
[[nodiscard]] inline double cylinder_height_gradient(const Point3D& p, double radius, double height,
|
||
double h = 1e-6)
|
||
{
|
||
return (cylinder(p, radius, height + h) - cylinder(p, radius, height - h)) / (2.0 * h);
|
||
}
|
||
|
||
/**
|
||
* @brief 圆柱 SDF 对半径的参数梯度 ∂f/∂radius
|
||
*
|
||
* @return 参数梯度
|
||
*/
|
||
[[nodiscard]] inline double cylinder_radius_gradient(const Point3D& p, double radius, double height,
|
||
double h = 1e-6)
|
||
{
|
||
return (cylinder(p, radius + h, height) - cylinder(p, radius - h, height)) / (2.0 * h);
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Analytical Gradients (Exact)
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 球体 SDF 的解析梯度
|
||
*
|
||
* **公式:** ∇f = p / |p|
|
||
*
|
||
* 即从球心指向采样点的单位径向向量。
|
||
* 注意：球心处 p = (0,0,0)，梯度退化，返回零向量。
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 球半径（当前未使用，保留以保持接口一致）
|
||
* @return 解析梯度（单位向量）
|
||
*
|
||
* @note 在表面上 (|p| = radius)，梯度即为外法向量。
|
||
*/
|
||
[[nodiscard]] inline Vector3D sphere_gradient_analytic(const Point3D& p, double /*radius*/)
|
||
{
|
||
double n = p.norm();
|
||
if (n < 1e-30) return Vector3D::Zero();
|
||
return p / n;
|
||
}
|
||
|
||
/**
|
||
* @brief 长方体 SDF 的解析梯度
|
||
*
|
||
* **算法:**
|
||
* - 外部: 梯度分量非零当且仅当对应坐标超出 semi-extent，方向向外
|
||
* - 内部: 梯度指向最近的表面（单分量）
|
||
*
|
||
* 注意：长方体 SDF 使用平滑近似，边缘处的梯度在内部/外部是分段常数。
|
||
*
|
||
* @param p 采样点
|
||
* @param half_extents 半边长
|
||
* @return 解析梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D box_gradient_analytic(const Point3D& p, const Point3D& half_extents)
|
||
{
|
||
Vector3D g = Vector3D::Zero();
|
||
if (std::abs(p.x()) > half_extents.x()) g.x() = (p.x() > 0.0 ? 1.0 : -1.0);
|
||
if (std::abs(p.y()) > half_extents.y()) g.y() = (p.y() > 0.0 ? 1.0 : -1.0);
|
||
if (std::abs(p.z()) > half_extents.z()) g.z() = (p.z() > 0.0 ? 1.0 : -1.0);
|
||
|
||
// For interior: gradient points toward nearest face
|
||
if (g.norm() < 1e-30 && half_extents.norm() > 1e-30) {
|
||
// Find the dimension where the point is closest to a face
|
||
Point3D q = Point3D(
|
||
half_extents.x() - std::abs(p.x()),
|
||
half_extents.y() - std::abs(p.y()),
|
||
half_extents.z() - std::abs(p.z())
|
||
);
|
||
// Smallest q component = closest face
|
||
if (q.x() <= q.y() && q.x() <= q.z()) g.x() = (p.x() >= 0.0 ? 1.0 : -1.0);
|
||
else if (q.y() <= q.x() && q.y() <= q.z()) g.y() = (p.y() >= 0.0 ? 1.0 : -1.0);
|
||
else g.z() = (p.z() >= 0.0 ? 1.0 : -1.0);
|
||
}
|
||
return g;
|
||
}
|
||
|
||
/**
|
||
* @brief 圆环 SDF 的解析梯度
|
||
*
|
||
* **SDF:** f = √((√(pₓ²+p𝓏²) - R_major)² + p_y²) - R_minor
|
||
*
|
||
* 梯度通过链式法则推导，将环向距离和 Y 方向分别求导。
|
||
*
|
||
* **特殊情况:** 在 Y 轴上 (pₓ=p𝓏=0)，x 和 z 梯度分量退化为 0。
|
||
*
|
||
* @param p 采样点
|
||
* @param major_r 主半径
|
||
* @param minor_r 副半径
|
||
* @return 解析梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D torus_gradient_analytic(const Point3D& p, double major_r, double minor_r)
|
||
{
|
||
double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
|
||
if (rho < 1e-30) {
|
||
// On the Y axis: x and z gradients are 0 by symmetry;
|
||
// y gradient = py / sqrt(major² + py²)
|
||
double inner = std::sqrt(major_r * major_r + p.y() * p.y());
|
||
if (inner < 1e-30) return Vector3D::Zero();
|
||
return Vector3D(0.0, p.y() / inner, 0.0);
|
||
}
|
||
double inner = std::sqrt((rho - major_r) * (rho - major_r) + p.y() * p.y());
|
||
if (inner < 1e-30) return Vector3D::Zero();
|
||
double d_rho = (rho - major_r) / inner;
|
||
return Vector3D(
|
||
d_rho * p.x() / rho,
|
||
p.y() / inner,
|
||
d_rho * p.z() / rho
|
||
);
|
||
}
|
||
|
||
/**
|
||
* @brief 无限圆柱（沿 Y 轴）的解析梯度
|
||
*
|
||
* **公式:** ∇f = (pₓ, 0, p𝓏) / √(pₓ² + p𝓏²)
|
||
*
|
||
* 梯度仅含水平分量，方向为从轴线径向向外。
|
||
* 在轴线上退化返回 (0,1,0)。
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 圆柱半径（未使用，保留接口一致）
|
||
* @return 解析梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D cylinder_gradient_analytic(const Point3D& p, double /*radius*/)
|
||
{
|
||
double rho = std::sqrt(p.x() * p.x() + p.z() * p.z());
|
||
if (rho < 1e-30) return Vector3D::UnitY(); // degenerate case: along the axis
|
||
return Vector3D(p.x() / rho, 0.0, p.z() / rho);
|
||
}
|
||
|
||
/**
|
||
* @brief 平面的解析梯度
|
||
*
|
||
* **公式:** ∇f = normal
|
||
*
|
||
* 平面的梯度为常向量（法向量），与位置无关。
|
||
*
|
||
* @param normal 单位法向量
|
||
* @return 梯度（等于 normal）
|
||
*
|
||
* @note 这是最简单的解析梯度，因为平面 SDF 是线性的。
|
||
*/
|
||
[[nodiscard]] inline Vector3D plane_gradient_analytic(const Vector3D& normal)
|
||
{
|
||
return normal;
|
||
}
|
||
|
||
/**
|
||
* @brief 胶囊体 SDF 的解析梯度
|
||
*
|
||
* 梯度指向点 p 到最近线段上的投影方向。
|
||
*
|
||
* **算法:**
|
||
* 1. 找到 p 在线段 ab 上的投影点 closest
|
||
* 2. 若 |p-closest| ≈ 0（点在轴上），返回零向量
|
||
* 3. 否则梯度 = (p - closest) / |p - closest|
|
||
*
|
||
* 退化情况: 当 a ≈ b（线段退化为点）时，退化为球体的径向梯度。
|
||
*
|
||
* @param p 采样点
|
||
* @param a 线段起点
|
||
* @param b 线段终点
|
||
* @param radius 扫掠球半径
|
||
* @return 解析梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D capsule_gradient_analytic(const Point3D& p,
|
||
const Point3D& a, const Point3D& b,
|
||
double radius)
|
||
{
|
||
Point3D pa = p - a;
|
||
Point3D ba = b - a;
|
||
double ba_sq = ba.squaredNorm();
|
||
if (ba_sq < 1e-30) {
|
||
// Degenerate to sphere
|
||
double n = pa.norm();
|
||
if (n < 1e-30) return Vector3D::Zero();
|
||
return pa / n;
|
||
}
|
||
double h = std::clamp(pa.dot(ba) / ba_sq, 0.0, 1.0);
|
||
Point3D closest = a + ba * h;
|
||
Point3D delta = p - closest;
|
||
double d = delta.norm();
|
||
if (d < 1e-30) return Vector3D::Zero();
|
||
return delta / d;
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Chain Rule for CSG Operations
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 并集操作后的梯度传递
|
||
*
|
||
* 取两个子对象中较近者的梯度（因为距离值来自较近者）。
|
||
*
|
||
* @param grad_a 子对象 A 的梯度
|
||
* @param grad_b 子对象 B 的梯度
|
||
* @param d1 A 的 SDF 值
|
||
* @param d2 B 的 SDF 值
|
||
* @return 合并后的梯度
|
||
*
|
||
* @warning 当 d1 ≈ d2 时，梯度可能在两个子对象间不连续跳变。
|
||
*
|
||
* @see chain_smooth_union 避免不连续的平滑版本
|
||
*/
|
||
[[nodiscard]] inline Vector3D chain_union(const Vector3D& grad_a, const Vector3D& grad_b,
|
||
double d1, double d2)
|
||
{
|
||
return (d1 < d2) ? grad_a : grad_b;
|
||
}
|
||
|
||
/**
|
||
* @brief 交集操作后的梯度传递
|
||
*
|
||
* 取两个子对象中较远者的梯度（因为距离值来自较远者）。
|
||
*
|
||
* @return 合并后的梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D chain_intersection(const Vector3D& grad_a, const Vector3D& grad_b,
|
||
double d1, double d2)
|
||
{
|
||
return (d1 > d2) ? grad_a : grad_b;
|
||
}
|
||
|
||
/**
|
||
* @brief 差集操作后的梯度传递
|
||
*
|
||
* 若 A 决定结果 (d1 > -d2)，使用 ∇A；否则使用 -∇B（B 的内部变成外部）。
|
||
*
|
||
* @param grad_a 子对象 A 的梯度
|
||
* @param grad_b 子对象 B 的梯度
|
||
* @param d1 f_A
|
||
* @param d2 f_B
|
||
* @return 合并后的梯度
|
||
*/
|
||
[[nodiscard]] inline Vector3D chain_difference(const Vector3D& grad_a, const Vector3D& grad_b,
|
||
double d1, double d2)
|
||
{
|
||
return (d1 > -d2) ? grad_a : Vector3D(-grad_b.x(), -grad_b.y(), -grad_b.z());
|
||
}
|
||
|
||
/**
|
||
* @brief 平滑并集操作的梯度传递（含完整链式法则）
|
||
*
|
||
* 与普通并集的硬切换不同，本函数使用多项式混合函数的导数来平滑插值两个子对象的梯度。
|
||
*
|
||
* **推导:**
|
||
* f = d₁(1-h) + d₂h - k·h(1-h)，h = clamp(0.5+0.5(d₂-d₁)/k, 0, 1)
|
||
* ∂f/∂d₁ = 2 - 3h
|
||
* ∂f/∂d₂ = 3h - 1
|
||
*
|
||
* 在混合区域 (0 < h < 1) 内，权重平滑地从 A 过渡到 B。
|
||
*
|
||
* @param grad_a 子对象 A 的梯度
|
||
* @param grad_b 子对象 B 的梯度
|
||
* @param d1 f_A
|
||
* @param d2 f_B
|
||
* @param k 混合参数
|
||
* @return 平滑合并后的梯度
|
||
*
|
||
* @see op_smooth_union 对应的平滑并集 SDF 操作
|
||
* @see chain_union 硬切换版本
|
||
*/
|
||
[[nodiscard]] inline Vector3D chain_smooth_union(const Vector3D& grad_a, const Vector3D& grad_b,
|
||
double d1, double d2, double k)
|
||
{
|
||
if (k <= 0.0) return chain_union(grad_a, grad_b, d1, d2);
|
||
double h = std::clamp(0.5 + 0.5 * (d1 - d2) / k, 0.0, 1.0);
|
||
if (h <= 0.0) return grad_a;
|
||
if (h >= 1.0) return grad_b;
|
||
|
||
// Full derivatives w.r.t. d1, d2:
|
||
// ∂f/∂d1 = 2 - 3h, ∂f/∂d2 = 3h - 1
|
||
double w1 = 2.0 - 3.0 * h;
|
||
double w2 = 3.0 * h - 1.0;
|
||
return grad_a * w1 + grad_b * w2;
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Chain Rule for Domain Transforms
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 平移变换后的梯度传递（梯度不变）
|
||
*
|
||
* 平移是等距变换（isometry），函数值不变，
|
||
* 因此子对象的梯度直接穿透。
|
||
*
|
||
* @param child_result 变换空间中的子对象求值结果
|
||
* @param offset 平移量（仅接口保留）
|
||
* @return 相同梯度
|
||
*/
|
||
[[nodiscard]] inline GradResult chain_translate(const GradResult& child_result, const Point3D& /*offset*/)
|
||
{
|
||
return child_result;
|
||
}
|
||
|
||
/**
|
||
* @brief 绕 Y 轴旋转后的梯度传递
|
||
*
|
||
* 由于 op_rotate 对采样点施加 R(-θ)，链式法则要求对梯度施加 R(θ)。
|
||
* 即：∇f_world = R(θ) · ∇f_local
|
||
*
|
||
* @param child_result 旋转后坐标系中的子对象求值
|
||
* @param angle_rad 旋转角度
|
||
* @return 世界坐标系的梯度
|
||
*
|
||
* @see op_rotate 对应的域变形操作
|
||
*/
|
||
[[nodiscard]] inline GradResult chain_rotate(const GradResult& child_result, double angle_rad)
|
||
{
|
||
GradResult result = child_result;
|
||
double c = std::cos(angle_rad);
|
||
double s = std::sin(angle_rad);
|
||
// R(θ) = [[c, 0, -s], [0, 1, 0], [s, 0, c]]
|
||
result.grad = Vector3D(
|
||
c * child_result.grad.x() - s * child_result.grad.z(),
|
||
child_result.grad.y(),
|
||
s * child_result.grad.x() + c * child_result.grad.z()
|
||
);
|
||
return result;
|
||
}
|
||
|
||
/**
|
||
* @brief 非均匀缩放后的梯度传递
|
||
*
|
||
* **公式:** ∇f_world = (∇f_child.x / sx, ∇f_child.y / sy, ∇f_child.z / sz)
|
||
*
|
||
* 因为 f_world(p) = f_child(p / s)，由链式法则得 ∂f_world/∂p = ∇f_child · diag(1/s)。
|
||
*
|
||
* @note 缩放后梯度不再是单位长度，除非是均匀缩放。
|
||
*
|
||
* @param child_result 缩放后坐标系中的子对象求值
|
||
* @param factors 缩放因子
|
||
* @return 世界坐标系的梯度
|
||
*/
|
||
[[nodiscard]] inline GradResult chain_scale(const GradResult& child_result, const Point3D& factors)
|
||
{
|
||
GradResult result = child_result;
|
||
result.grad = Vector3D(
|
||
child_result.grad.x() / factors.x(),
|
||
child_result.grad.y() / factors.y(),
|
||
child_result.grad.z() / factors.z()
|
||
);
|
||
return result;
|
||
}
|
||
|
||
/**
|
||
* @brief 扭曲变换后的梯度传递（含完整雅可比）
|
||
*
|
||
* 扭曲 T(p) 绕 Y 轴旋转角度 amount·p_y。
|
||
* 梯度传递需要对雅可比矩阵的转置做链式法则：
|
||
*
|
||
* **雅可比 J:**
|
||
* ```
|
||
* ∂T_x/∂p_x = c, ∂T_x/∂p_y = -a·(p_x·s + p_z·c), ∂T_x/∂p_z = -s
|
||
* ∂T_y/∂p_x = 0, ∂T_y/∂p_y = 1, ∂T_y/∂p_z = 0
|
||
* ∂T_z/∂p_x = s, ∂T_z/∂p_y = a·(p_x·c - p_z·s), ∂T_z/∂p_z = c
|
||
* ```
|
||
* 其中 c = cos(a·p_y), s = sin(a·p_y), a = amount。
|
||
*
|
||
* @param child_result 扭曲空间中的子对象求值
|
||
* @param amount 单位高度的旋转量
|
||
* @param p 原始采样点
|
||
* @return 世界坐标系的梯度
|
||
*
|
||
* @see op_twist 对应的域变形操作
|
||
* @see chain_rotate 不含 Y 依赖的简化旋转梯度传递
|
||
*/
|
||
[[nodiscard]] inline GradResult chain_twist(const GradResult& child_result,
|
||
double amount, const Point3D& p)
|
||
{
|
||
double angle = amount * p.y();
|
||
double c = std::cos(angle);
|
||
double s = std::sin(angle);
|
||
|
||
// Forward-twisted position
|
||
double tp_x = p.x() * c - p.z() * s;
|
||
double tp_z = p.x() * s + p.z() * c;
|
||
|
||
GradResult result = child_result;
|
||
|
||
// Spatial part: inverse rotation R(-angle) applied to child gradient
|
||
double gx = c * child_result.grad.x() + s * child_result.grad.z();
|
||
double gz = -s * child_result.grad.x() + c * child_result.grad.z();
|
||
|
||
// Y-component: extra term from angle's dependence on p.y
|
||
// ∂/∂y correction: amount * (tp_x * gz_child - tp_z * gx_child)
|
||
double gy = child_result.grad.y() + amount * (tp_x * child_result.grad.z()
|
||
- tp_z * child_result.grad.x());
|
||
|
||
result.grad = Vector3D(gx, gy, gz);
|
||
return result;
|
||
}
|
||
|
||
/**
|
||
* @brief 周期重复后的梯度传递（梯度不变）
|
||
*
|
||
* 重复映射是平移变换的叠加，每个晶格内梯度相同（周期域）。
|
||
* 在晶格边界处存在不连续性。
|
||
*
|
||
* @param child_result 中心晶格内的子对象求值
|
||
* @param cell 晶格尺寸（仅接口保留）
|
||
* @param p 原始采样点（仅接口保留）
|
||
* @return 相同梯度
|
||
*
|
||
* @warning 在晶格边界处梯度不连续，可能影响基于梯度的优化。
|
||
*/
|
||
[[nodiscard]] inline GradResult chain_repeat(const GradResult& child_result,
|
||
const Point3D& /*cell*/, const Point3D& /*p*/)
|
||
{
|
||
return child_result;
|
||
}
|
||
|
||
// ═══════════════════════════════════════════════
|
||
// Parameter Gradient Vector (All Shape Params)
|
||
// ═══════════════════════════════════════════════
|
||
|
||
/**
|
||
* @brief 参数梯度结构体
|
||
*
|
||
* 存储 SDF 对所有形状参数的偏导数：
|
||
* ∂f/∂radius, ∂f/∂height, ∂f/∂extent[3], ∂f/∂major, ∂f/∂minor,
|
||
* ∂f/∂angle, ∂f/∂translate[3], ∂f/∂rotate, ∂f/∂scale[3]
|
||
*
|
||
* 用于基于梯度的形状优化和参数拟合。
|
||
*
|
||
* @see sphere_param_grad, box_param_grad, cylinder_param_grad
|
||
*/
|
||
struct ParamGrad {
|
||
double d_radius = 0.0; ///< ∂f/∂radius
|
||
double d_height = 0.0; ///< ∂f/∂height
|
||
Point3D d_extents = Point3D::Zero(); ///< ∂f/∂(extent.x, extent.y, extent.z)
|
||
double d_major = 0.0; ///< ∂f/∂major_radius
|
||
double d_minor = 0.0; ///< ∂f/∂minor_radius
|
||
double d_angle = 0.0; ///< ∂f/∂angle_rad
|
||
Vector3D d_normal = Vector3D::Zero(); ///< ∂f/∂normal (无效分量清零)
|
||
double d_offset = 0.0; ///< ∂f/∂offset
|
||
Point3D d_translate = Point3D::Zero(); ///< ∂f/∂translate_offset
|
||
double d_rotate = 0.0; ///< ∂f/∂rotation_angle
|
||
Point3D d_scale = Point3D::Zero(); ///< ∂f/∂scale_factors
|
||
};
|
||
|
||
/**
|
||
* @brief 球体所有参数的梯度（当前仅 radius 有意义）
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 球体半径
|
||
* @param h 差分步长
|
||
* @return 参数梯度（仅 d_radius 非零）
|
||
*
|
||
* @see sphere_radius_gradient 底层数值差分
|
||
*/
|
||
[[nodiscard]] inline ParamGrad sphere_param_grad(const Point3D& p, double radius, double h = 1e-6)
|
||
{
|
||
ParamGrad pg;
|
||
pg.d_radius = sphere_radius_gradient(p, radius, h);
|
||
return pg;
|
||
}
|
||
|
||
/**
|
||
* @brief 长方体所有参数的梯度
|
||
*
|
||
* 分别计算三个轴的半尺寸对 SDF 值的影响。
|
||
*
|
||
* @param p 采样点
|
||
* @param half_extents 半边长
|
||
* @param h 差分步长
|
||
* @return 参数梯度（d_extents 分量填充）
|
||
*
|
||
* @see box_extent_gradient 底层数值差分
|
||
*/
|
||
[[nodiscard]] inline ParamGrad box_param_grad(const Point3D& p, const Point3D& half_extents,
|
||
double h = 1e-6)
|
||
{
|
||
ParamGrad pg;
|
||
pg.d_extents.x() = box_extent_gradient(p, half_extents, 0, h);
|
||
pg.d_extents.y() = box_extent_gradient(p, half_extents, 1, h);
|
||
pg.d_extents.z() = box_extent_gradient(p, half_extents, 2, h);
|
||
return pg;
|
||
}
|
||
|
||
/**
|
||
* @brief 圆柱所有参数的梯度
|
||
*
|
||
* 分别计算半径和高度对 SDF 值的影响。
|
||
*
|
||
* @param p 采样点
|
||
* @param radius 圆柱半径
|
||
* @param height 圆柱高度
|
||
* @param h 差分步长
|
||
* @return 参数梯度（d_radius, d_height 填充）
|
||
*/
|
||
[[nodiscard]] inline ParamGrad cylinder_param_grad(const Point3D& p, double radius, double height,
|
||
double h = 1e-6)
|
||
{
|
||
ParamGrad pg;
|
||
pg.d_radius = cylinder_radius_gradient(p, radius, height, h);
|
||
pg.d_height = cylinder_height_gradient(p, radius, height, h);
|
||
return pg;
|
||
}
|
||
|
||
} // namespace vde::sdf
|