茂之钳
467 lines
16 KiB
C++
467 lines
16 KiB
C++
#pragma once
|
||
/**
|
||
* @file sdf_optimize.h
|
||
* @brief 基于梯度的 SDF 形状优化
|
||
*
|
||
* 提供 SDF 形状的拟合、碰撞解析、可及性分析和对称性检测等功能。
|
||
*
|
||
* ## 主要功能
|
||
*
|
||
* | 功能 | 说明 |
|
||
* |------------------------|-------------------------------------------------|
|
||
* | Shape Fitting | 将参数化形状拟合到点云或目标 SDF |
|
||
* | Collision Avoidance | 平移形状以解决穿透 |
|
||
* | Accessibility | 采样半球以评估表面上点的可及性 |
|
||
* | Symmetry Detection | 检测反射对称性并找到最佳对称平面 |
|
||
* | Volume & Mass | 通过蒙地卡罗法估计体积和质心 |
|
||
* | Parameter Optimization | 收集所有可优化的参数并计算数值梯度 |
|
||
*
|
||
* ## 使用模式
|
||
*
|
||
* ```cpp
|
||
* auto shape = SdfNode::sphere(0.5);
|
||
* auto result = fit_to_point_cloud(shape, points, 0.01, 200);
|
||
* if (result.converged) {
|
||
* // result.optimized_shape 包含优化后的参数
|
||
* }
|
||
* ```
|
||
*
|
||
* @ingroup sdf
|
||
*/
|
||
#include "vde/sdf/sdf_tree.h"
|
||
#include "vde/core/point.h"
|
||
#include <vector>
|
||
#include <functional>
|
||
|
||
namespace vde::sdf {
|
||
|
||
using core::Point3D;
|
||
using core::Vector3D;
|
||
|
||
// ────────────────────────────────────────────────
|
||
// Shape Optimization
|
||
// ────────────────────────────────────────────────
|
||
|
||
/**
|
||
* @brief 优化结果
|
||
*
|
||
* 包含优化后的形状、损失值、收敛标志和迭代历史。
|
||
*/
|
||
struct OptimizeResult {
|
||
SdfNodePtr optimized_shape; ///< 优化后的 SDF 树（参数已更新）
|
||
double final_loss; ///< 最终损失值
|
||
int iterations; ///< 实际迭代次数
|
||
bool converged; ///< 是否收敛
|
||
std::vector<double> loss_history; ///< 每次迭代的损失记录
|
||
};
|
||
|
||
/**
|
||
* @brief 损失函数类型
|
||
*
|
||
* 接受一个 SDF 求值函数（封装了参数化的形状），返回标量损失。
|
||
*
|
||
* @see fit_to_point_cloud 使用该类型的内置实现
|
||
*/
|
||
using LossFn = std::function<double(const std::function<double(const Point3D&)>&)>;
|
||
|
||
// ── Shape Fitting ───────────────────────────────
|
||
|
||
/**
|
||
* @brief 将参数化形状拟合到目标点云
|
||
*
|
||
* 最小化目标表面点到拟合形状表面的有符号距离。
|
||
*
|
||
* **算法概要:**
|
||
* 1. 遍历 target_points 计算每个点到初始形状的 SDF 值
|
||
* 2. 用有限差分计算形状参数对损失的梯度
|
||
* 3. 梯度下降更新参数
|
||
* 4. 重复直至收敛或达到最大迭代次数
|
||
*
|
||
* @param initial 初始形状（sphere, box, cylinder 等），其参数将被优化
|
||
* @param target_points 目标点云（表面采样点）
|
||
* @param learning_rate 梯度下降步长（默认 0.01）
|
||
* @param max_iterations 最大迭代次数（默认 100）
|
||
* @return 优化结果（含优化后的形状和收敛信息）
|
||
*
|
||
* @note 损失函数为 L = mean(|sdf(p_i)|)，即点到表面的平均绝对距离。
|
||
* @note initial 的共享子节点不会被复制——返回的 optimize_shape 共享未修改的子树。
|
||
*
|
||
* @code{.cpp}
|
||
* std::vector<Point3D> points = { point cloud data };
|
||
* auto sphere = SdfNode::sphere(1.0);
|
||
* auto result = fit_to_point_cloud(sphere, points, 0.01, 200);
|
||
* // result.optimized_shape 的 params.radius 已被优化
|
||
* @endcode
|
||
*
|
||
* @see fit_surface_to_points 类似但仅最小化 |sdf| 而非有符号距离
|
||
* @see fit_to_sdf 拟合到另一个 SDF
|
||
*/
|
||
[[nodiscard]] OptimizeResult fit_to_point_cloud(
|
||
const SdfNodePtr& initial,
|
||
const std::vector<Point3D>& target_points,
|
||
double learning_rate = 0.01,
|
||
int max_iterations = 100);
|
||
|
||
/**
|
||
* @brief 将形状表面移动到目标采样点
|
||
*
|
||
* 与 fit_to_point_cloud 类似，但损失函数使用 L = mean(|sdf(p_i)|)。
|
||
* 这强制采样点位于表面上（而非仅靠近表面）。
|
||
*
|
||
* @param initial 初始形状
|
||
* @param surface_points 目标表面点
|
||
* @param learning_rate 学习率
|
||
* @param max_iterations 最大迭代次数
|
||
* @return 优化结果
|
||
*
|
||
* @see fit_to_point_cloud 更通用的点云拟合（可包含内部/外部点）
|
||
*/
|
||
[[nodiscard]] OptimizeResult fit_surface_to_points(
|
||
const SdfNodePtr& initial,
|
||
const std::vector<Point3D>& surface_points,
|
||
double learning_rate = 0.01,
|
||
int max_iterations = 100);
|
||
|
||
/**
|
||
* @brief 拟合一个 SDF 形状以匹配另一个目标 SDF
|
||
*
|
||
* 在规则网格上采样两个 SDF，最小化它们之间的均方差。
|
||
*
|
||
* **损失函数:** L = mean((sdf_source(p_i) - sdf_target(p_i))²)
|
||
*
|
||
* 适合将一个参数化形状"画像"成一个更复杂的隐式形状。
|
||
*
|
||
* @param source 被优化的参数化形状
|
||
* @param target_sdf 目标 SDF 函数
|
||
* @param bmin 采样网格最小角点
|
||
* @param bmax 采样网格最大角点
|
||
* @param grid_resolution 采样网格分辨率（默认 16）
|
||
* @param learning_rate 学习率
|
||
* @param max_iterations 最大迭代次数
|
||
* @return 优化结果
|
||
*
|
||
* @warning 网格分辨率每增加一倍，计算量增加 8 倍。使用 8~32 可获得合理速度。
|
||
*
|
||
* @see fit_to_point_cloud 拟合到点云
|
||
*/
|
||
[[nodiscard]] OptimizeResult fit_to_sdf(
|
||
const SdfNodePtr& source,
|
||
const std::function<double(const Point3D&)>& target_sdf,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int grid_resolution = 16,
|
||
double learning_rate = 0.01,
|
||
int max_iterations = 100);
|
||
|
||
// ── Collision Avoidance ─────────────────────────
|
||
|
||
/**
|
||
* @brief 通过平移解决两个形状之间的穿透
|
||
*
|
||
* 在穿透区域采样点，沿梯度方向平移形状以分离它们。
|
||
* 修改传入的形状（原地修改 translate_offset 参数）。
|
||
*
|
||
* @param shape_a 形状 A（原地修改）
|
||
* @param shape_b 形状 B（原地修改）
|
||
* @param step_size 每次迭代的平移步长（默认 0.1）
|
||
* @param max_iterations 最大迭代次数（默认 50）
|
||
* @return true 如果碰撞被成功解决；false 如果在迭代次数内未能解决
|
||
*
|
||
* @code{.cpp}
|
||
* auto sphere_a = SdfNode::sphere(1.0);
|
||
* auto sphere_b = SdfNode::translate(SdfNode::sphere(1.0), Point3D(1.5, 0, 0));
|
||
* bool resolved = resolve_collision(sphere_a, sphere_b, 0.1, 100);
|
||
* @endcode
|
||
*
|
||
* @see penetration_depth 估算穿透深度
|
||
*/
|
||
[[nodiscard]] bool resolve_collision(
|
||
SdfNodePtr& shape_a, SdfNodePtr& shape_b,
|
||
double step_size = 0.1,
|
||
int max_iterations = 50);
|
||
|
||
/**
|
||
* @brief 估算两个形状之间的最小穿透距离
|
||
*
|
||
* 在包围盒内随机采样，找到两个 SDF 均 ≤ 0 的区域（重叠区），
|
||
* 计算重叠点到表面的最大距离。
|
||
*
|
||
* @param shape_a 形状 A
|
||
* @param shape_b 形状 B
|
||
* @param bmin 公共包围盒最小角点
|
||
* @param bmax 公共包围盒最大角点
|
||
* @param samples 蒙地卡罗采样数（默认 1000）
|
||
* @return 最大穿透深度（负值表示不重叠时的最小间距）
|
||
*
|
||
* @see resolve_collision 解决穿透
|
||
*/
|
||
[[nodiscard]] double penetration_depth(
|
||
const SdfNodePtr& shape_a,
|
||
const SdfNodePtr& shape_b,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int samples = 1000);
|
||
|
||
// ── Reachability / Accessibility ─────────────────
|
||
|
||
/**
|
||
* @brief 计算形状表面上点的可及性分数
|
||
*
|
||
* 从表面点 p 发射半球光线，统计未被形状遮挡的比例。
|
||
* 返回 0（完全不可及）到 1（完全可及）之间的分数。
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param p 表面点（应满足 |sdf(p)| < ε）
|
||
* @param direction 接近方向
|
||
* @param hemisphere_samples 半球采样光线数量（默认 64）
|
||
* @return 可及性分数 [0, 1]
|
||
*
|
||
* @note 采样数越多越精确，但线性增加计算量。
|
||
*
|
||
* @see find_accessible_point 找到最大可及性的表面点
|
||
*/
|
||
[[nodiscard]] double accessibility(
|
||
const SdfNodePtr& shape,
|
||
const Point3D& p, // surface point
|
||
const Vector3D& direction, // approach direction
|
||
int hemisphere_samples = 64);
|
||
|
||
/**
|
||
* @brief 找到给定接近方向下具有最大可及性的表面点
|
||
*
|
||
* 在包围盒内的规则网格上搜索，计算每个网格点处的可及性分数。
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param approach_dir 接近方向（如工具接近方向）
|
||
* @param bmin 搜索包围盒最小角点
|
||
* @param bmax 搜索包围盒最大角点
|
||
* @param grid_res 搜索网格分辨率（默认 32）
|
||
* @return 最大可及性的表面点
|
||
*
|
||
* @see accessibility 单点可及性计算
|
||
*/
|
||
[[nodiscard]] Point3D find_accessible_point(
|
||
const SdfNodePtr& shape,
|
||
const Vector3D& approach_dir,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int grid_res = 32);
|
||
|
||
// ── Symmetry Detection ──────────────────────────
|
||
|
||
/**
|
||
* @brief 检测形状在给定方向上的反射对称性
|
||
*
|
||
* 在形状的包围盒内随机采样点对 (p, p') 其中 p' 是关于给定平面的反射。
|
||
* 比较两点的 SDF 值以评估对称程度。
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param bmin 包围盒最小角点
|
||
* @param bmax 包围盒最大角点
|
||
* @param plane_normal 反射平面的法向量
|
||
* @param plane_offset 反射平面的偏移量（沿法向量方向，默认 0）
|
||
* @param samples 采样点对数（默认 1000）
|
||
* @return 对称性分数 [0, 1]（0 = 完全不对称，1 = 完美对称）
|
||
*
|
||
* @see find_symmetry_plane 自动寻找最佳对称平面
|
||
*/
|
||
[[nodiscard]] double symmetry_score(
|
||
const SdfNodePtr& shape,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
const Vector3D& plane_normal,
|
||
double plane_offset = 0.0,
|
||
int samples = 1000);
|
||
|
||
/**
|
||
* @brief 对称检测结果
|
||
*/
|
||
struct SymmetryResult {
|
||
double score; ///< 对称性分数 [0, 1]
|
||
Vector3D normal; ///< 最佳对称平面的法向量
|
||
double offset; ///< 最佳对称平面的偏移量
|
||
};
|
||
|
||
/**
|
||
* @brief 自动找到形状的最佳反射对称平面
|
||
*
|
||
* 尝试多个候选平面方向（沿主轴和若干采样方向），
|
||
* 返回对称性分数最高的平面。
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param bmin 包围盒最小角点
|
||
* @param bmax 包围盒最大角点
|
||
* @param samples 每个候选方向的采样数（默认 1000）
|
||
* @return 最佳对称平面及其分数
|
||
*
|
||
* @see symmetry_score 给定方向的对称性检测
|
||
*/
|
||
[[nodiscard]] SymmetryResult find_symmetry_plane(
|
||
const SdfNodePtr& shape,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int samples = 1000);
|
||
|
||
// ── SDF Volume Computation ──────────────────────
|
||
|
||
/**
|
||
* @brief 通过蒙地卡罗法估算隐式形状的体积
|
||
*
|
||
* 在包围盒内均匀随机采样，统计落在形状内部的点比例，
|
||
* 乘以包围盒体积得到近似体积。
|
||
*
|
||
* **公式:** V ≈ V_box · N_inside / N_total，其中 N_inside = count(sdf(p_i) < 0)
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param bmin 采样包围盒最小角点
|
||
* @param bmax 采样包围盒最大角点
|
||
* @param samples 蒙地卡罗采样数（默认 10000）
|
||
* @return 近似体积
|
||
*
|
||
* @note 精度与 √samples 成正比。10000 采样 ≈ 1% 相对误差。
|
||
* @note 包围盒应尽可能紧致以减小方差。
|
||
*
|
||
* @see center_of_mass 蒙地卡罗质心估算
|
||
*/
|
||
[[nodiscard]] double estimate_volume(
|
||
const SdfNodePtr& shape,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int samples = 10000);
|
||
|
||
/**
|
||
* @brief 通过蒙地卡罗法估算形状的质心
|
||
*
|
||
* 对形状内部的采样点取算术平均。
|
||
*
|
||
* **公式:** COM ≈ Σ p_i / N_inside，其中 sdf(p_i) < 0
|
||
*
|
||
* @param shape SDF 形状
|
||
* @param bmin 采样包围盒最小角点
|
||
* @param bmax 采样包围盒最大角点
|
||
* @param samples 蒙地卡罗采样数（默认 10000）
|
||
* @return 近似质心坐标
|
||
*
|
||
* @see estimate_volume 体积估算
|
||
*/
|
||
[[nodiscard]] Point3D center_of_mass(
|
||
const SdfNodePtr& shape,
|
||
const Point3D& bmin, const Point3D& bmax,
|
||
int samples = 10000);
|
||
|
||
// ── Parameter Collection ────────────────────────
|
||
|
||
/**
|
||
* @brief SDF 树中可变参数的引用
|
||
*
|
||
* 指向 SdfNode::params 中具体数值字段的指针，
|
||
* 配合 collect_params() 收集所有可优化参数。
|
||
*/
|
||
struct ParamRef {
|
||
double* value; ///< 指向参数值的可写指针
|
||
std::string name; ///< 参数名称（调试用）
|
||
};
|
||
|
||
/**
|
||
* @brief 从叶节点收集所有可变数值参数
|
||
*
|
||
* 遍历 SDF 树，收集所有图元节点中可优化的参数引用。
|
||
* 返回的 ParamRef 列表可传递给 numerical_gradient() 进行优化。
|
||
*
|
||
* @param root SDF 树根节点（可写引用，因为参数将被原地修改）
|
||
* @return 可变参数引用列表
|
||
*
|
||
* @code{.cpp}
|
||
* auto root = SdfNode::op_union(
|
||
* SdfNode::sphere(1.0),
|
||
* SdfNode::box(Point3D(0.5, 0.5, 0.5))
|
||
* );
|
||
* auto params = collect_params(root);
|
||
* // params[0].value 指向球体半径
|
||
* // params[1..4] 指向盒子的 extents 分量
|
||
* @endcode
|
||
*
|
||
* @see numerical_gradient 计算损失对参数的数值梯度
|
||
*/
|
||
[[nodiscard]] std::vector<ParamRef> collect_params(SdfNodePtr& root);
|
||
|
||
/**
|
||
* @brief 计算标量损失对收集到的参数的数值梯度
|
||
*
|
||
* 对每个参数逐个做中心差分，返回梯度向量。
|
||
*
|
||
* @param params collect_params() 返回的参数引用列表
|
||
* @param loss_fn 损失函数（无参数的可调用对象）
|
||
* @param eps 差分步长（默认 1e-6）
|
||
* @return 梯度向量 grad_loss，grad_loss[i] = ∂loss / ∂params[i].value
|
||
*
|
||
* @code{.cpp}
|
||
* auto params = collect_params(root);
|
||
* auto loss = [&]() {
|
||
* double s = 0;
|
||
* for (auto& pt : target_points) s += std::abs(evaluate(root, pt));
|
||
* return s / target_points.size();
|
||
* };
|
||
* auto grad = numerical_gradient(params, loss, 1e-6);
|
||
* @endcode
|
||
*
|
||
* @see collect_params 收集参数
|
||
* @see GradientDescent 梯度下降优化器
|
||
*/
|
||
[[nodiscard]] std::vector<double> numerical_gradient(
|
||
const std::vector<ParamRef>& params,
|
||
const std::function<double()>& loss_fn,
|
||
double eps = 1e-6);
|
||
|
||
// ── Gradient Descent Utility ────────────────────
|
||
|
||
/**
|
||
* @brief 固定学习率的简单梯度下降优化器
|
||
*
|
||
* 管理迭代计数和学习率，提供 step() 方法执行单步参数更新。
|
||
*
|
||
* @code{.cpp}
|
||
* auto params = collect_params(root);
|
||
* GradientDescent gd(0.01);
|
||
* for (int i = 0; i < 100; ++i) {
|
||
* auto loss = gd.step(param_vals, [&](auto& vals, auto& grad) {
|
||
* // 填充 grad，返回 loss
|
||
* });
|
||
* if (loss < tol) break;
|
||
* }
|
||
* @endcode
|
||
*/
|
||
class GradientDescent {
|
||
public:
|
||
/**
|
||
* @brief 构造优化器
|
||
* @param lr 学习率（步长因子）
|
||
*/
|
||
explicit GradientDescent(double lr) : lr_(lr) {}
|
||
|
||
/**
|
||
* @brief 执行单步梯度下降
|
||
*
|
||
* 调用 grad_fn 获取当前参数的梯度和损失值，
|
||
* 然后沿负梯度方向更新参数: param_i -= lr · grad_i。
|
||
*
|
||
* @param params 参数向量（原地修改）
|
||
* @param grad_fn 梯度计算函数，签名为
|
||
* double(const std::vector<double>& vals, std::vector<double>& grad_out)
|
||
* 接受当前参数值，填充 grad_out 并返回损失标量。
|
||
* @return 当前损失值
|
||
*/
|
||
double step(std::vector<double>& params,
|
||
const std::function<double(const std::vector<double>&,
|
||
std::vector<double>&)>& grad_fn);
|
||
|
||
/** @brief 设置学习率 */
|
||
void set_learning_rate(double lr) { lr_ = lr; }
|
||
|
||
/** @brief 获取学习率 */
|
||
[[nodiscard]] double learning_rate() const { return lr_; }
|
||
|
||
/** @brief 获取当前迭代次数 */
|
||
[[nodiscard]] int iteration() const { return iteration_; }
|
||
|
||
private:
|
||
double lr_;
|
||
int iteration_ = 0;
|
||
};
|
||
|
||
} // namespace vde::sdf
|