13.4 尺度空间曲面重建

尺度空间重建是曲面重建系列中处理噪声点云的有效方法，通过多尺度平滑获得高质量重建结果。

尺度空间曲面重建（Scale Space Surface Reconstruction）是一种基于多尺度分析的重建方法。它通过对点云进行迭代平滑（尺度空间），然后在不同尺度上重建表面，最后将结果传播回原始点云。

13.4.1 数学基础

尺度空间理论

尺度空间理论由 Witkin 和 Koenderink 提出，核心思想是将信号在不同尺度上表示。对于点云，尺度空间通过迭代平滑算子构建：

$P^{(k + 1)} = S (P^{(k)})$

其中 $P^{(0)} = P$ 是原始点云， $S$ 是平滑算子， $k$ 是尺度级别。

加权 PCA 平滑

CGAL 使用加权主成分分析（PCA）作为平滑算子。对于每个点 $p_{i}$ ，在其邻域 $N (p_{i})$ 内计算加权协方差矩阵：

$C = \sum_{p_{j} \in N (p_{i})} w_{ij} (p_{j} - \overset{p}{ˉ}) (p_{j} - \overset{p}{ˉ})^{T}$

其中：

$\overset{p}{ˉ} = \sum_{p_{j} \in N (p_{i})} w_{ij} p_{j}$ 是加权质心
$w_{ij} = 1/∣ N (p_{j}) ∣$ 是基于局部密度的权重

对 $C$ 进行特征分解：

$C = V Λ V^{T}$

其中 $Λ = diag (λ_{1}, λ_{2}, λ_{3})$ ， $λ_{1} \geq λ_{2} \geq λ_{3}$ 。

点 $p_{i}$ 被投影到由 $λ_{2}$ 和 $λ_{3}$ 张成的平面上：

$p_{i}^{(n e w)} = \overset{p}{ˉ} + \sum_{k = 2}^{3} ((p_{i} - \overset{p}{ˉ}) \cdot v_{k}) v_{k}$

Alpha Shape 表面提取

在特定尺度上，使用 Alpha Shape 从平滑后的点云提取表面：

$Alpha (r) = {σ \in Delaunay (P) ∣ \exists empty ball of radius \leq r through vertices of σ}$

其中 $r$ 是 alpha 半径参数。

13.4.2 算法流程

尺度空间重建包含以下阶段：

1. 尺度空间构建

输入: 点云 P, 迭代次数 K
输出: 多尺度点云序列 P^(0), P^(1), ..., P^(K)

P^(0) = P
for k = 1 to K:
    P^(k) = smooth(P^(k-1))

2. 表面重建

在最终尺度 $P^{(K)}$ 上重建表面：

构建 P^(K) 的 Delaunay 三角剖分
提取 Alpha Shape（半径 = 估计的邻域半径）
收集正则（regular）和奇异（singular）面片

3. 索引传播

将重建结果的三角形索引映射回原始点云：

对于每个三角形 (i, j, k) 在尺度 K:
    输出三角形 (index^(0)(i), index^(0)(j), index^(0)(k))

由于尺度空间保持点的顺序，索引直接对应原始点云。

13.4.3 架构分析

核心类层次

Scale_space_surface_reconstruction_3<Geom_traits>
    |
    +-- Weighted_PCA_smoother<Geom_traits> (默认平滑器)
    |       +-- Search_tree (k-d 树用于邻域查询)
    |
    +-- Alpha_shape_mesher<Geom_traits> (默认网格器)
            +-- Shape_construction_3<Geom_traits>
                    +-- Alpha_shape_3 / Fixed_alpha_shape_3

关键文件位置

文件	路径	功能
`Scale_space_surface_reconstruction_3.h`	`Scale_space_reconstruction_3/include/CGAL/`	核心重建类
`Weighted_PCA_smoother.h`	`Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/`	加权 PCA 平滑
`Alpha_shape_mesher.h`	`Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/`	Alpha Shape 网格提取
`Shape_construction_3.h`	`Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/`	形状构造工具

核心类定义

template <typename Geom_traits>
class Scale_space_surface_reconstruction_3
{
public:
    // 几何类型
    typedef typename Geom_traits::FT FT;
    typedef typename Geom_traits::Point_3 Point;
    typedef std::array<std::size_t, 3> Facet;
    
    // 点云存储
    typedef std::vector<Point> Point_range;
    typedef typename Point_range::iterator Point_iterator;
    typedef typename Point_range::const_iterator Point_const_iterator;
    
    // 面片存储
    typedef std::vector<Facet> Facet_range;
    typedef typename Facet_range::iterator Facet_iterator;
    typedef typename Facet_range::const_iterator Facet_const_iterator;
    
    // 默认算法
    typedef Scale_space_reconstruction_3::Weighted_PCA_smoother<Geom_traits> 
        Weighted_PCA_smoother;
    typedef Scale_space_reconstruction_3::Alpha_shape_mesher<Geom_traits> 
        Alpha_shape_mesher;
    
    // 构造函数
    Scale_space_surface_reconstruction_3();
    
    template <typename InputIterator>
    Scale_space_surface_reconstruction_3(InputIterator begin, InputIterator end);
    
    // 点云操作
    void insert(const Point& p);
    template <typename InputIterator>
    void insert(InputIterator begin, InputIterator end);
    
    // 尺度空间操作
    template <typename Smoother>
    void increase_scale(std::size_t iterations = 1, const Smoother& smoother = Smoother());
    void increase_scale(std::size_t iterations = 1);  // 使用默认平滑器
    
    // 表面重建
    template <typename Mesher>
    void reconstruct_surface(const Mesher& mesher = Mesher());
    void reconstruct_surface();  // 使用默认网格器
    
    // 结果访问
    std::size_t number_of_points() const;
    std::size_t number_of_facets() const;
    
    const Point_range& points() const;
    Point_iterator points_begin();
    Point_iterator points_end();
    
    const Facet_range& facets() const;
    Facet_iterator facets_begin();
    Facet_iterator facets_end();
};

13.4.4 实现细节

加权 PCA 平滑器

template <typename Geom_traits, typename DiagonalizeTraits, typename ConcurrencyTag>
class Weighted_PCA_smoother
{
public:
    typedef typename Geom_traits::FT FT;
    typedef typename Geom_traits::Point_3 Point;
    typedef typename Geom_traits::Vector_3 Vector;
    
    // 构造函数
    Weighted_PCA_smoother(unsigned int neighbors = 12, unsigned int samples = 300);
    Weighted_PCA_smoother(FT squared_radius);
    
    // 平滑操作
    template <typename InputIterator>
    void operator()(InputIterator begin, InputIterator end);
    
    FT squared_radius() { return _squared_radius; }
    
private:
    // 估计邻域半径
    void estimate_neighborhood_squared_radius();
    
    // 邻域计数
    class ComputeNN {
        void operator()(const std::size_t& i) const {
            Sphere sp(_pts[i], _sq_rd);
            Inc inc(_nn[i]);
            _tree.search(boost::make_function_output_iterator(inc), sp);
        }
    };
    
    // 尺度推进
    class AdvanceSS {
        void operator()(const std::size_t& i) const {
            if (_nn[i] < 4) return;  // 邻域太小，不移动
            
            // k-近邻搜索
            Static_search search(_tree, _pts[i], _nn[i]);
            
            // 计算加权质心
            Point barycenter(0., 0., 0.);
            FT weight_sum = 0.;
            for (auto nit = search.begin(); nit \!= search.end(); ++nit) {
                weight_sum += (1.0 / _nn[boost::get<1>(nit->first)]);
            }
            for (auto nit = search.begin(); nit \!= search.end(); ++nit) {
                Vector v(CGAL::ORIGIN, boost::get<0>(nit->first));
                barycenter = barycenter + ((1.0 / _nn[boost::get<1>(nit->first)]) / weight_sum) * v;
            }
            
            // 计算加权协方差矩阵
            std::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
            for (auto nit = search.begin(); nit \!= search.end(); ++nit) {
                Vector v(barycenter, boost::get<0>(nit->first));
                FT w = (1.0 / _nn[boost::get<1>(nit->first)]);
                v = w * v;
                covariance[0] += w * v.x() * v.x();
                covariance[1] += w * v.x() * v.y();
                covariance[2] += w * v.x() * v.z();
                covariance[3] += w * v.y() * v.y();
                covariance[4] += w * v.y() * v.z();
                covariance[5] += w * v.z() * v.z();
            }
            
            // 特征分解
            std::array<FT, 9> eigenvectors = {{ 0., 0., 0., 0., 0., 0., 0., 0., 0. }};
            std::array<FT, 3> eigenvalues = {{ 0., 0., 0. }};
            DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix(
                covariance, eigenvalues, eigenvectors);
            
            // 投影到平面（去除最小特征值方向）
            Vector norm(eigenvectors[0], eigenvectors[1], eigenvectors[2]);
            Vector b2p(barycenter, _pts[i]);
            _pts[i] = barycenter + b2p - ((norm * b2p) * norm);
        }
    };
};

Alpha Shape 网格提取

template <typename Geom_traits, typename FixedSurface>
class Alpha_shape_mesher
{
public:
    typedef typename Geom_traits::FT FT;
    typedef typename Geom_traits::Point_3 Point;
    typedef std::array<unsigned int, 3> Facet;
    
    // 构造函数
    Alpha_shape_mesher(FT squared_radius,
                       bool separate_shells = false,
                       bool force_manifold = false,
                       FT border_angle = 45.);
    
    // 网格提取
    template <typename InputIterator, typename OutputIterator>
    void operator()(InputIterator begin, InputIterator end, OutputIterator output);
    
    // 结果访问
    std::size_t number_of_triangles() const;
    std::size_t number_of_shells() const;
    
    Facet_const_iterator surface_begin() const;
    Facet_const_iterator surface_end() const;
    Facet_const_iterator shell_begin(std::size_t shell) const;
    Facet_const_iterator shell_end(std::size_t shell) const;
    
private:
    // 快速收集面片（无序）
    void collect_facets_quick();
    
    // 按壳收集面片
    void collect_facets();
    void collect_shell(const SFacet& f);
    
    // 气泡检测（用于流形约束）
    void detect_bubbles();
    
    // 修复非流形边/顶点
    void fix_nonmanifold_edges();
    void fix_nonmanifold_vertices();
    
    // 标记处理
    inline bool is_handled(Cell_handle c, unsigned int li) const;
    inline void mark_handled(Cell_handle c, unsigned int li);
    inline void mark_handled(SFacet f);
};

面片收集算法

void collect_shell(Cell_handle c, unsigned int li)
{
    std::stack<SFacet> stack;
    stack.push(SFacet(c, li));
    
    SFacet f;
    bool processed = false;
    
    while (\!stack.empty()) {
        f = stack.top();
        stack.pop();
        
        if (is_handled(f))
            continue;
        
        // 添加面片到表面
        mark_handled(f);
        _surface.push_back(ordered_facet_indices(f));
        
        if (\!processed) {
            if (_separate_shells || _shells.size() == 0) {
                _shells.push_back(--_surface.end());
            }
            processed = true;
        }
        
        // 保存壳索引
        _map_f2s[f] = _index - 1;
        
        // 标记对侧（用于流形约束）
        if (_force_manifold)
            mark_opposite_handled(f);
        
        // 沿边旋转，找到相邻的正则/奇异面片
        for (int i = 0; i < 4; ++i) {
            if (i == f.second || is_handled(f.first, i))
                continue;
            
            Cell_handle n = f.first;
            int ni = i;
            Vertex_handle a = f.first->vertex(f.second);
            Classification_type cl = _shape->classify(SFacet(n, ni));
            
            // 旋转直到找到正则或奇异面片
            while (cl \!= Shape::REGULAR && cl \!= Shape::SINGULAR) {
                Cell_handle p = n;
                n = n->neighbor(ni);
                ni = n->index(a);
                int pi = n->index(p);
                a = n->vertex(pi);
                cl = _shape->classify(SFacet(n, ni));
            }
            
            stack.push(SFacet(n, ni));
        }
    }
}

13.4.5 使用示例

完整示例代码

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Timer.h>
 
#include <fstream>
#include <iostream>
 
// 类型定义
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Scale_space_surface_reconstruction_3<Kernel> Reconstruction;
 
typedef Reconstruction::Facet_const_iterator Facet_iterator;
 
int main(int argc, char** argv)
{
    // 1. 读取点云
    std::string fname = (argc == 1) ? 
        CGAL::data_file_path("points_3/kitten.off") : argv[1];
    
    std::cout << "Reading point cloud..." << std::flush;
    std::vector<Point> points;
    
    if (\!CGAL::IO::read_points(fname, std::back_inserter(points))) {
        std::cerr << "Error: cannot read file" << std::endl;
        return EXIT_FAILURE;
    }
    std::cout << " done: " << points.size() << " points." << std::endl;
    
    // 2. 创建重建对象
    CGAL::Timer t;
    t.start();
    
    Reconstruction reconstruct(points.begin(), points.end());
    
    // 3. 增加尺度（平滑）
    std::cout << "Increasing scale..." << std::flush;
    reconstruct.increase_scale(4);  // 4 次迭代
    std::cout << " done in " << t.time() << " sec." << std::endl;
    t.reset();
    
    // 4. 重建表面
    std::cout << "Reconstructing surface..." << std::flush;
    reconstruct.reconstruct_surface();
    std::cout << " done in " << t.time() << " sec." << std::endl;
    
    // 5. 输出统计
    std::cout << "Result:" << std::endl;
    std::cout << "  Points: " << reconstruct.number_of_points() << std::endl;
    std::cout << "  Facets: " << reconstruct.number_of_facets() << std::endl;
    
    // 6. 保存结果
    std::ofstream out("output.off");
    out << reconstruct;  // 使用重载的 << 操作符
    out.close();
    
    std::cout << "Output saved to output.off" << std::endl;
    
    return EXIT_SUCCESS;
}

使用自定义平滑器

#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
 
// 使用 Jet 拟合平滑器
typedef CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel> Jet_smoother;
 
Reconstruction reconstruct(points.begin(), points.end());
 
// 使用 Jet 平滑器，4 次迭代
Jet_smoother jet_smoother;
reconstruct.increase_scale(4, jet_smoother);

使用自定义网格器

#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
 
// 使用推进前沿网格器
typedef CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel> 
    Advancing_front_mesher;
 
Reconstruction reconstruct(points.begin(), points.end());
reconstruct.increase_scale(4);
 
// 使用推进前沿重建表面
Advancing_front_mesher af_mesher;
reconstruct.reconstruct_surface(af_mesher);

增量式重建

// 增量式处理，逐步增加尺度并观察结果
Reconstruction reconstruct;
 
// 插入点
reconstruct.insert(points.begin(), points.end());
 
// 逐步增加尺度
for (std::size_t i = 1; i <= 4; ++i) {
    reconstruct.increase_scale(1);  // 每次增加 1 级
    
    // 在当前尺度重建
    reconstruct.reconstruct_surface();
    
    // 保存中间结果
    std::ofstream out("scale_" + std::to_string(i) + ".off");
    out << reconstruct;
    out.close();
}

13.4.6 复杂度分析

时间复杂度

步骤	复杂度	说明
邻域半径估计	$O (n lo g n)$	基于 k-d 树
每次平滑迭代	$O (n \cdot k lo g n)$	$k$ 为平均邻域点数
Delaunay 三角剖分	$O (n lo g n)$	最终尺度
Alpha Shape 提取	$O (n lo g n)$	面片收集

总体复杂度为 $O (K \cdot n \cdot k lo g n)$ ，其中 $K$ 是迭代次数。

空间复杂度

点云存储： $O (n)$ （每尺度）
k-d 树： $O (n)$
Delaunay 三角剖分： $O (n)$
表面存储： $O (m)$ ， $m$ 为输出面片数

总体空间复杂度为 $O (n + m)$ 。

参数影响

neighbors：邻域大小。较大的值产生更平滑的结果，但计算量更大。
iterations：迭代次数。更多的迭代产生更平滑的表面，但可能丢失细节。
squared_radius：Alpha 半径。控制表面细节程度。

13.4.7 应用建议

适用场景

尺度空间重建最适合以下场景：

含有噪声的点云
需要插值原始点的应用
多尺度分析需求
无向点云（不需要法向）

预处理建议

离群点移除：虽然尺度空间对噪声鲁棒，但严重离群点仍会影响结果
采样密度：确保采样足够密集以捕捉表面细节
迭代次数选择：根据噪声水平和细节保留需求调整

常见问题

问题 1：重建结果过于平滑，丢失细节

原因：迭代次数过多或邻域过大
解决：减少迭代次数，或减小 neighbors 参数

问题 2：重建结果有孔洞

原因：Alpha 半径过大或点云采样不足
解决：调整 Alpha 半径，或使用更小的邻域

问题 3：非流形几何

原因：点云自相交或过于复杂
解决：启用 force_manifold 选项，或预处理分割点云

参考文献

Digne, J., Morel, J. M., Souzani, C. M., & Lartigue, C. (2011). Scale space meshing of raw data point sets. Computer Graphics Forum, 30(6), 1630-1642.
Witkin, A. P. (1983). Scale-space filtering. In Readings in Computer Vision (pp. 329-332).
Edelsbrunner, H., & Mücke, E. P. (1994). Three-dimensional alpha shapes. ACM Transactions on Graphics, 13(1), 43-72.
CGAL Documentation: Scale Space Surface Reconstruction. https://doc.cgal.org/latest/Scale_space_reconstruction_3/

Chunibyo

Explorer

13.4 尺度空间曲面重建

13.4 尺度空间曲面重建

13.4.1 数学基础

尺度空间理论

加权 PCA 平滑

Alpha Shape 表面提取

13.4.2 算法流程

1. 尺度空间构建

2. 表面重建

3. 索引传播

13.4.3 架构分析

核心类层次

关键文件位置

核心类定义

13.4.4 实现细节

加权 PCA 平滑器

Alpha Shape 网格提取

面片收集算法

13.4.5 使用示例

完整示例代码

使用自定义平滑器

使用自定义网格器

增量式重建

13.4.6 复杂度分析

时间复杂度

空间复杂度

参数影响

13.4.7 应用建议

适用场景

预处理建议

常见问题

参考文献

Graph View

Table of Contents

Backlinks

Chunibyo

Explorer

13.4 尺度空间曲面重建

13.4 尺度空间曲面重建 §

13.4.1 数学基础 §

尺度空间理论 §

加权 PCA 平滑 §

Alpha Shape 表面提取 §

13.4.2 算法流程 §

1. 尺度空间构建 §

2. 表面重建 §

3. 索引传播 §

13.4.3 架构分析 §

核心类层次 §

关键文件位置 §

核心类定义 §

13.4.4 实现细节 §

加权 PCA 平滑器 §

Alpha Shape 网格提取 §

面片收集算法 §

13.4.5 使用示例 §

完整示例代码 §

使用自定义平滑器 §

使用自定义网格器 §

增量式重建 §

13.4.6 复杂度分析 §

时间复杂度 §

空间复杂度 §

参数影响 §

13.4.7 应用建议 §

适用场景 §

预处理建议 §

常见问题 §

参考文献 §

Graph View

Table of Contents

Backlinks

13.4 尺度空间曲面重建

13.4.1 数学基础

尺度空间理论

加权 PCA 平滑

Alpha Shape 表面提取

13.4.2 算法流程

1. 尺度空间构建

2. 表面重建

3. 索引传播

13.4.3 架构分析

核心类层次

关键文件位置

核心类定义

13.4.4 实现细节

加权 PCA 平滑器

Alpha Shape 网格提取

面片收集算法

13.4.5 使用示例

完整示例代码

使用自定义平滑器

使用自定义网格器

增量式重建

13.4.6 复杂度分析

时间复杂度

空间复杂度

参数影响

13.4.7 应用建议

适用场景

预处理建议

常见问题

参考文献