10.2 Alpha Shapes

上一章介绍了凸包算法，Alpha Shapes 是凸包概念的推广，提供了一个从点集提取形状的连续谱系。

0. 算法动机：为什么需要 Alpha Shapes？

现实问题引入

想象你是一位考古学家，在沙漠中发现了一组散落的古代建筑遗迹坐标。你需要回答：

这些遗迹原本构成了什么样的聚落形状？
哪些建筑属于核心区域，哪些属于边缘？
如果考虑不同的防御半径，聚落的边界在哪里？

或者你是一位生物学家，研究蛋白质分子的表面结构：

如何定义分子的”表面”？
小分子能否进入蛋白质的某个口袋区域？
不同大小的探针分子会探测到什么样的表面？

传统方法的局限

凸包（Convex Hull） 给出了最外层的边界，但它太粗糙了——就像用橡皮筋包裹物体，会填满所有凹陷。

Delaunay 三角剖分 保留了所有可能的连接，但它太详细了——包含了内部的许多不需要的边。

Alpha Shapes 提供了一个连续的谱系，从精细到粗糙，让你可以选择合适的细节层次。

如果不使用 Alpha Shapes

形状识别错误：将分离的物体识别为一个整体
特征丢失：小的凹陷或孔洞被忽略
噪声敏感：离群点严重影响形状判断
多尺度分析困难：无法同时观察宏观和微观结构

1. 理论基础

1.1 Alpha Shapes 的定义

Alpha Shapes 是由 Edelsbrunner 和 Muecke 提出的一种几何结构，用于描述点集的”形状”。它通过参数 alpha 控制形状的精细程度，是 Delaunay 三角剖分的一个子复形。

定义 10.3 (Alpha 复形)：给定 d 维空间中的点集 S 和实数 alpha ≥ 0，alpha 复形是 Delaunay 三角剖分中所有满足以下条件的单形组成的集合：

存在一个半径为 sqrt(alpha) 的空球（不包含其他点）通过该单形的顶点
该单形的所有面也属于 alpha 复形

定义 10.4 (Alpha Shape)：alpha shape 是 alpha 复形的几何实现，即 alpha 复形中所有单形的并集。

1.2 与相关概念的关系

alpha = 0        → 点集本身
alpha 很小       → 点集的精细结构
alpha = ∞        → 凸包
Delaunay 三角剖分 → 所有 alpha 值的并集

2. 直观理解：融冰与雕刻

2.1 “融冰”类比

想象你的点集是冻结在冰块中的气泡。现在你开始加热：

Alpha = 0（冰冻状态）：

每个气泡都是孤立的点
没有任何连接
就像看着悬浮在冰中的单个气泡

    •           •
         
              •
   •
         •

Alpha 很小（开始融化）：

靠近的气泡开始通过薄冰连接
只有非常接近的点才有边相连
就像看到气泡之间的微小通道

    •           •
     \\         /
      •———•   
     /     \\  •
    •       •

Alpha 适中（半融化）：

形成了有意义的形状轮廓
小的凹陷显现，但太小的缝隙被忽略
就像看到物体的大致轮廓

    •———•
   / \\   \\
  •   •———•
   \\ /   /
    •———•

Alpha = ∞（完全融化）：

所有气泡合并成一个连通区域
这就是凸包
所有凹陷都被填满

      __________
     / •———•    \\
    / / \\   \\    \\
   |  •   •———•   |
    \\ \\ /   /   /
     \\•———•____/

2.2 “雕刻”类比

另一种理解方式是雕刻：

想象你有一块巨大的石头，里面镶嵌着许多珍珠（点）。你使用一个球形钻头（半径为 sqrt(alpha)）来雕刻：

小球钻头（小 alpha）：只能雕刻出珍珠周围很小的区域，保留很多细节
中等钻头（中等 alpha）：雕刻出珍珠之间的通道，形成轮廓
大钻头（大 alpha）：雕刻出整个外形，填满所有凹陷

2.3 数学直觉

关键洞察：Alpha Shape 不是直接定义的，而是通过”空球条件”间接定义的。

对于一个三角形（或边、顶点），问一个问题：

“能找到一个多大的空球，让它刚好接触这个单形的所有顶点？”

如果空球半径 sqrt(alpha) 足够大 → 该单形属于 alpha shape
如果太小 → 被排除在外

为什么用空球？

保证结果与 Delaunay 三角剖分兼容
提供连续的形状变化
有清晰的几何解释

3. Alpha Shapes 与相关概念的关系图

3.1 概念层次图

                    凸包 (Convex Hull)
                         |
                         | alpha = ∞
                         v
              +---------------------+
              |                     |
         Alpha Shape            Alpha Shape
         (大 alpha)             (小 alpha)
              |                     |
              |                     |
              v                     v
       +--------------+      +--------------+
       |   粗粒度形状  |      |   细粒度形状  |
       |   填充凹陷   |      |   保留凹陷   |
       +--------------+      +--------------+
              |                     |
              |                     |
              v                     v
              +---------------------+
                         |
                         v
              Delaunay 三角剖分
                  (所有边)
                         |
                         v
              点集本身 (alpha = 0)

3.2 不同 Alpha 值效果对比

考虑一个简单点集：U形分布的点

Alpha = 0:                    Alpha = 小:
•       •                     •———•
•       •                     |   |
•       •                     |   |
•       •                     |   |
 •     •                       \\ /
  •   •                         •
   • •                         / \\
    •                         •   •

Alpha = 中:                   Alpha = ∞ (凸包):
•———•                         •———•
|   |                        /     \\
|   |                       /       \\
|   |                      |         |
|   |                      |         |
 \\ /                        \\       /
  •                          \\     /
 / \\                          •———•
•   •

3.3 与 Gabriel 图的关系

Gabriel 图是 Alpha Shape 在 alpha → 0 极限时的特例：

Gabriel 边：以该边为直径的圆内不包含其他点
对应于 alpha_min 有定义的边
是最”紧致”的连接方式

4. CGAL Alpha Shapes 架构

4.1 模块组织

Alpha_shapes_2/         # 二维 Alpha Shapes
  include/CGAL/
    Alpha_shape_2.h           # 主类
    Alpha_shape_vertex_base_2.h  # 顶点基类
    Alpha_shape_face_base_2.h    # 面基类
    internal/Lazy_alpha_nt_2.h   # 惰性数值类型

Alpha_shapes_3/         # 三维 Alpha Shapes
  include/CGAL/
    Alpha_shape_3.h           # 主类
    Alpha_shape_vertex_base_3.h  # 顶点基类
    Alpha_shape_cell_base_3.h    # 单元基类
    internal/Lazy_alpha_nt_3.h   # 惰性数值类型

4.2 类层次结构

Delaunay_triangulation_2/3
        ↑
    Alpha_shape_2/3
        
关键模板参数：
- Dt: 底层三角剖分类型 (Delaunay_triangulation_2/3)
- ExactAlphaComparisonTag: 是否使用精确 alpha 比较

5. Alpha Shape 2D 实现详解

5.1 核心类定义

// Alpha_shape_2.h
template <class Dt, class ExactAlphaComparisonTag = Tag_false>
class Alpha_shape_2 : public Dt {
public:
    typedef Dt Triangulation;
    typedef typename Dt::Geom_traits Gt;
    typedef typename Dt::Triangulation_data_structure Tds;
    
    // Alpha 值类型
    typedef typename internal::Alpha_nt_selector_2<
        Gt, ExactAlphaComparisonTag, typename Dt::Weighted_tag>::Type_of_alpha Type_of_alpha;
    
    // 分类类型
    enum Classification_type {EXTERIOR, SINGULAR, REGULAR, INTERIOR};
    
    // 模式
    enum Mode {GENERAL, REGULARIZED};
    
private:
    // 区间映射表
    Interval_face_map _interval_face_map;    // 面 → alpha 区间
    Interval_edge_map _interval_edge_map;    // 边 → alpha 区间
    Interval_vertex_map _interval_vertex_map; // 顶点 → alpha 区间
    
    Alpha_spectrum _alpha_spectrum;  // 排序后的所有 alpha 值
    Type_of_alpha _alpha;            // 当前 alpha 值
    Mode _mode;                      // 当前模式
};

5.2 区间计算

// 面的 alpha 区间初始化
template <class Dt, class EACT>
void Alpha_shape_2<Dt, EACT>::initialize_interval_face_map() {
    Type_of_alpha alpha_f;
    
    // 遍历所有有限面
    for (Finite_faces_iterator face_it = faces_begin(); 
         face_it \!= faces_end(); ++face_it) {
        // 计算外接圆半径的平方
        alpha_f = squared_radius(face_it);
        _interval_face_map.insert(Interval_face(alpha_f, face_it));
        
        // 设置交叉引用
        face_it->set_alpha(alpha_f);
    }
}
 
// 边的 alpha 区间初始化（核心算法）
template <class Dt, class EACT>
void Alpha_shape_2<Dt, EACT>::initialize_interval_edge_map() {
    for (edge_it = edges_begin(); edge_it \!= edges_end(); ++edge_it) {
        Interval3 interval;
        Edge edge = *edge_it;
        
        Face_handle pFace = edge.first;
        int i = edge.second;
        Face_handle pNeighbor = pFace->neighbor(i);
        int Neigh_i = pNeighbor->index(pFace);
        
        if (\!is_infinite(pFace) && \!is_infinite(pNeighbor)) {
            // 内部边：两个相邻面都存在
            Type_of_alpha squared_radius_Face = find_interval(pFace);
            Type_of_alpha squared_radius_Neighbor = find_interval(pNeighbor);
            
            // 确保正确的顺序
            if (squared_radius_Neighbor < squared_radius_Face) {
                std::swap(edge, squared_radius_Face, squared_radius_Neighbor);
            }
            
            // 判断是否为 Gabriel 边
            interval = (is_attached(pFace, i) || is_attached(pNeighbor, Neigh_i)) ?
                make_triple(UNDEFINED, squared_radius_Face, squared_radius_Neighbor) :
                make_triple(squared_radius(pFace, i), squared_radius_Face, 
                           squared_radius_Neighbor);
        } else {
            // 凸包边：只有一个相邻面
            if (is_infinite(pFace)) {
                if (\!is_infinite(pNeighbor)) {
                    interval = is_attached(pNeighbor, Neigh_i) ?
                        make_triple(UNDEFINED, find_interval(pNeighbor), Infinity) :
                        make_triple(squared_radius(pNeighbor, Neigh_i), 
                                   find_interval(pNeighbor), Infinity);
                } else {
                    // 两个面都是无限的（退化情况）
                    interval = make_triple(squared_radius(pFace, i), 
                                          Infinity, Infinity);
                }
            }
            // ...
        }
        
        _interval_edge_map.insert(Interval_edge(interval, edge));
        (edge.first)->set_ranges(edge.second, interval);
    }
}

5.3 区间表示详解

对于每个单形，使用三元组表示其 alpha 区间：

// Interval3 = (alpha_min, alpha_mid, alpha_max)
typedef Triple<Type_of_alpha, Type_of_alpha, Type_of_alpha> Interval_3;
 
// 面的区间：只有一个值（外接圆半径）
// face_alpha = 外接圆半径的平方
 
// 边的区间：
// - alpha_min: Gabriel 半径（如果是 Gabriel 边）或 UNDEFINED
// - alpha_mid: 较小相邻面的外接圆半径
// - alpha_max: 较大相邻面的外接圆半径或 Infinity
 
// 顶点的区间：
// - alpha_mid: 最小 incident 面的 alpha
// - alpha_max: 最大 incident 面的 alpha 或 Infinity

5.4 分类算法

// 面的分类
template <class Dt, class EACT>
typename Alpha_shape_2<Dt, EACT>::Classification_type
Alpha_shape_2<Dt, EACT>::classify(const Face_handle& f, 
                                    const Type_of_alpha& alpha) const {
    if (is_infinite(f)) return EXTERIOR;
    
    // 面的区间就是其外接圆半径
    return (find_interval(f) <= alpha) ? INTERIOR : EXTERIOR;
}
 
// 边的分类（核心逻辑）
template <class Dt, class EACT>
typename Alpha_shape_2<Dt, EACT>::Classification_type
Alpha_shape_2<Dt, EACT>::classify(const Face_handle& f, int i,
                                    const Type_of_alpha& alpha) const {
    if (is_infinite(f, i)) return EXTERIOR;
    
    Interval3 interval = find_interval(std::make_pair(f, i));
    
    if (alpha < interval.second) {
        // alpha < alpha_mid
        if (get_mode() == REGULARIZED ||
            interval.first == UNDEFINED ||
            alpha < interval.first) {
            return EXTERIOR;
        } else {
            return SINGULAR;  // Gabriel 边且 alpha_min <= alpha < alpha_mid
        }
    } else {
        // alpha >= alpha_mid
        if (interval.third == Infinity || alpha < interval.third) {
            return REGULAR;  // 在边界上
        } else {
            return INTERIOR;  // 在内部
        }
    }
}
 
// 顶点的分类
template <class Dt, class EACT>
typename Alpha_shape_2<Dt, EACT>::Classification_type
Alpha_shape_2<Dt, EACT>::classify(const Vertex_handle& v,
                                    const Type_of_alpha& alpha) const {
    Interval2 interval = v->get_range();
    
    if (alpha < interval.first) {
        return (get_mode() == REGULARIZED) ? EXTERIOR : SINGULAR;
    } else if (interval.second == Infinity || alpha < interval.second) {
        return REGULAR;
    } else {
        return INTERIOR;
    }
}

6. Alpha Shape 3D 实现详解

6.1 核心类定义

// Alpha_shape_3.h
template <class Dt, class ExactAlphaComparisonTag = Tag_false>
class Alpha_shape_3 : public Dt {
public:
    typedef CGAL::Alpha_status<NT> Alpha_status;
    typedef Compact_container<Alpha_status> Alpha_status_container;
    
    // 分类类型和模式
    enum Classification_type {EXTERIOR, SINGULAR, REGULAR, INTERIOR};
    enum Mode {GENERAL, REGULARIZED};
    
private:
    // 使用 Alpha_status 存储每个单形的状态
    Alpha_status_container alpha_status_container;
    
    // 映射表
    Alpha_cell_map alpha_cell_map;      // 单元 → alpha
    Alpha_facet_map alpha_min_facet_map; // 面 → alpha_min
    Alpha_edge_map alpha_min_edge_map;   // 边 → alpha_min
    Alpha_vertex_map alpha_min_vertex_map; // 顶点 → alpha_min
    
    Edge_alpha_map edge_alpha_map;  // 边 → Alpha_status 迭代器
};
 
// Alpha_status 定义（Alpha_shape_cell_base_3.h）
template <class NT>
class Alpha_status {
    NT _alpha_min;   // 最小 alpha（Gabriel 半径）
    NT _alpha_mid;   // 中间 alpha（进入 alpha 复形的阈值）
    NT _alpha_max;   // 最大 alpha（离开 alpha 复形的阈值）
    bool _is_Gabriel;    // 是否是 Gabriel 单形
    bool _is_on_chull;   // 是否在凸包上
};

6.2 初始化流程

// 完整的 alpha shape 初始化
template <class Dt, class EACT>
void Alpha_shape_3<Dt, EACT>::initialize_alpha(bool reinitialize) {
    if (\!reinitialize) initialize_alpha_cell_map();
    initialize_alpha_facet_maps(reinitialize);
    initialize_alpha_edge_maps(reinitialize);
    initialize_alpha_vertex_maps(reinitialize);
    initialize_alpha_spectrum();
}
 
// 单元映射初始化
template <class Dt, class EACT>
void Alpha_shape_3<Dt, EACT>::initialize_alpha_cell_map() {
    for (cell_it = finite_cells_begin(); cell_it \!= done; ++cell_it) {
        alpha = squared_radius(cell_it);  // 外接球半径
        alpha_cell_map.insert(typename Alpha_cell_map::value_type(alpha, cell_it));
        cell_it->set_alpha(alpha);
    }
}
 
// 面映射初始化
template <class Dt, class EACT>
void Alpha_shape_3<Dt, EACT>::initialize_alpha_facet_maps(bool reinitialize) {
    for (fit = finite_facets_begin(); fit \!= finite_facets_end(); ++fit) {
        as = alpha_status_container.insert(Alpha_status());
        
        pCell = fit->first;
        i = fit->second;
        pNeighbor = pCell->neighbor(i);
        iNeigh = pNeighbor->index(pCell);
        
        if (\!is_infinite(pCell) && \!is_infinite(pNeighbor)) {
            // 内部面
            NT alpha_Cell = pCell->get_alpha();
            NT alpha_Neighbor = pNeighbor->get_alpha();
            
            alpha_mid = std::min(alpha_Cell, alpha_Neighbor);
            alpha_max = std::max(alpha_Cell, alpha_Neighbor);
            
            as->set_is_on_chull(false);
            as->set_alpha_mid(alpha_mid);
            as->set_alpha_max(alpha_max);
        } else {
            // 凸包面
            alpha_mid = \!is_infinite(pCell) ? pCell->get_alpha() 
                                            : pNeighbor->get_alpha();
            as->set_alpha_mid(alpha_mid);
            as->set_is_on_chull(true);
        }
        
        // 设置交叉引用
        pCell->set_facet_status(i, as);
        pNeighbor->set_facet_status(iNeigh, as);
    }
    
    // GENERAL 模式：计算 Gabriel 半径
    if (get_mode() == GENERAL && alpha_min_facet_map.empty()) {
        for (fit = finite_facets_begin(); fit \!= finite_facets_end(); ++fit) {
            as = fit->first->get_facet_status(fit->second);
            if (is_Gabriel(*fit)) {
                as->set_is_Gabriel(true);
                alpha_min = squared_radius(*fit);
                as->set_alpha_min(alpha_min);
                alpha_min_facet_map.insert(...);
            }
        }
    }
}

6.3 边状态计算

// 计算边的 alpha 状态（动态计算）
template <class Dt, class EACT>
void Alpha_shape_3<Dt, EACT>::compute_edge_status(
    const Cell_handle& c, int i, int j, 
    Alpha_status& as) const {
    
    as.set_is_on_chull(false);
    
    // 遍历所有 incident 单元
    Cell_circulator ccirc = incident_cells(c, i, j);
    last = ccirc;
    while (is_infinite(ccirc)) ++ccirc;
    
    alpha = (*ccirc).get_alpha();
    as.set_alpha_mid(alpha);
    as.set_alpha_max(alpha);
    
    while (++ccirc \!= last) {
        if (\!is_infinite(ccirc)) {
            alpha = (*ccirc).get_alpha();
            if (alpha < as.alpha_mid())
                as.set_alpha_mid(alpha);
            if (as.alpha_max() < alpha)
                as.set_alpha_max(alpha);
        }
    }
    
    // 检查 incident 面
    Facet_circulator fcirc = incident_facets(c, i, j);
    do {
        if (\!is_infinite(*fcirc)) {
            asf = (*fcirc).first->get_facet_status((*fcirc).second);
            if (asf->is_on_chull())
                as.set_is_on_chull(true);
        }
    } while (++fcirc \!= done);
    
    // GENERAL 模式：计算 Gabriel 半径
    if (get_mode() == GENERAL) {
        if (is_Gabriel(c, i, j)) {
            as.set_is_Gabriel(true);
            as.set_alpha_min(squared_radius(c, i, j));
        }
    }
}

7. 使用示例

7.1 基本 2D Alpha Shape

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_2.h>
#include <CGAL/Alpha_shape_vertex_base_2.h>
#include <CGAL/Alpha_shape_face_base_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <fstream>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_2 Point;
typedef K::Segment_2 Segment;
 
// 定义三角剖分数据结构
typedef CGAL::Alpha_shape_vertex_base_2<K> Vb;
typedef CGAL::Alpha_shape_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Triangulation_2;
typedef CGAL::Alpha_shape_2<Triangulation_2> Alpha_shape_2;
 
// 提取 alpha shape 边
template <class OutputIterator>
void alpha_edges(const Alpha_shape_2& A, OutputIterator out) {
    for (auto it = A.alpha_shape_edges_begin(); 
         it \!= A.alpha_shape_edges_end(); ++it) {
        *out++ = A.segment(*it);
    }
}
 
int main() {
    // 读取点集
    std::list<Point> points;
    std::ifstream is("data/fin");
    int n;
    is >> n;
    std::copy_n(std::istream_iterator<Point>(is), n, 
                std::back_inserter(points));
    
    // 构造 alpha shape
    Alpha_shape_2 A(points.begin(), points.end(),
                    FT(10000),           // alpha 值
                    Alpha_shape_2::GENERAL);  // 模式
    
    // 提取边
    std::vector<Segment> segments;
    alpha_edges(A, std::back_inserter(segments));
    
    std::cout << segments.size() << " alpha shape edges" << std::endl;
    
    // 查找最优 alpha（使所有点都在边界或内部）
    auto optimal = A.find_optimal_alpha(1);
    std::cout << "Optimal alpha: " << *optimal << std::endl;
    
    return 0;
}

7.2 3D Alpha Shape

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
 
// 定义 3D 三角剖分
typedef CGAL::Alpha_shape_vertex_base_3<K> Vb;
typedef CGAL::Alpha_shape_cell_base_3<K> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Dt;
typedef CGAL::Alpha_shape_3<Dt> Alpha_shape_3;
 
int main() {
    std::vector<Point> points;
    // ... 填充点集 ...
    
    // 构造 alpha shape（REGULARIZED 模式）
    Alpha_shape_3 A(points.begin(), points.end(),
                    0.0,  // 初始 alpha
                    Alpha_shape_3::REGULARIZED);
    
    // 设置 alpha 值
    A.set_alpha(10.0);
    
    // 遍历 alpha shape 面
    for (auto fit = A.alpha_shape_facets_begin();
         fit \!= A.alpha_shape_facets_end(); ++fit) {
        // 处理每个面
        auto classification = A.classify(*fit);
        if (classification == Alpha_shape_3::REGULAR) {
            // 在边界上
        }
    }
    
    // 获取连通分量数
    std::size_t n_components = A.number_of_solid_components();
    std::cout << "Solid components: " << n_components << std::endl;
    
    return 0;
}

7.3 使用 Weighted Alpha Shapes

// 带权重的 alpha shapes（使用 Regular 三角剖分）
#include <CGAL/Regular_triangulation_3.h>
 
typedef K::Weighted_point_3 Weighted_point;
typedef CGAL::Regular_triangulation_3<K, Tds> Rt;
typedef CGAL::Alpha_shape_3<Rt> Weighted_alpha_shape_3;
 
int main() {
    std::vector<Weighted_point> weighted_points;
    // ... 填充带权重点 ...
    
    Weighted_alpha_shape_3 A(weighted_points.begin(), 
                             weighted_points.end());
    
    // 使用方式与普通 alpha shape 相同
    return 0;
}

7.4 遍历 Alpha 谱

// 遍历所有不同的 alpha 值
void explore_alpha_spectrum(const Alpha_shape_2& A) {
    std::cout << "Alpha spectrum has " << A.number_of_alphas() 
              << " values" << std::endl;
    
    for (auto it = A.alpha_begin(); it \!= A.alpha_end(); ++it) {
        FT alpha = *it;
        std::cout << "Alpha = " << alpha << std::endl;
        
        // 临时设置 alpha 值进行分类
        FT old_alpha = A.set_alpha(alpha);
        
        // 统计各类单形数量
        int regular_edges = 0, singular_edges = 0;
        for (auto eit = A.alpha_shape_edges_begin();
             eit \!= A.alpha_shape_edges_end(); ++eit) {
            auto cls = A.classify(*eit);
            if (cls == Alpha_shape_2::REGULAR) ++regular_edges;
            else if (cls == Alpha_shape_2::SINGULAR) ++singular_edges;
        }
        
        std::cout << "  Regular edges: " << regular_edges << std::endl;
        std::cout << "  Singular edges: " << singular_edges << std::endl;
    }
}

7.5 实际应用：分子表面积计算

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <iostream>
#include <vector>
#include <fstream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef K::FT FT;
 
typedef CGAL::Alpha_shape_vertex_base_3<K> Vb;
typedef CGAL::Alpha_shape_cell_base_3<K> Cb;
typedef CGAL::Triangulation_data_structure_3<Vb, Cb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Dt;
typedef CGAL::Alpha_shape_3<Dt> Alpha_shape_3;
 
// 计算分子表面积（使用不同探针半径）
void compute_molecular_surface(const std::vector<Point>& atoms, 
                               double probe_radius) {
    // 构建 Alpha Shape
    Alpha_shape_3 alpha_shape(atoms.begin(), atoms.end(),
                              0.0,  // 初始 alpha
                              Alpha_shape_3::REGULARIZED);
    
    // 将探针半径转换为 alpha 值
    // alpha = probe_radius^2
    FT alpha = probe_radius * probe_radius;
    alpha_shape.set_alpha(alpha);
    
    // 计算表面积
    double surface_area = 0.0;
    for (auto fit = alpha_shape.alpha_shape_facets_begin();
         fit \!= alpha_shape.alpha_shape_facets_end(); ++fit) {
        if (alpha_shape.classify(*fit) == Alpha_shape_3::REGULAR) {
            // 获取面的三个顶点
            auto f = *fit;
            Point p1 = f.first->vertex((f.second+1)%4)->point();
            Point p2 = f.first->vertex((f.second+2)%4)->point();
            Point p3 = f.first->vertex((f.second+3)%4)->point();
            
            // 计算三角形面积
            K::Triangle_3 tri(p1, p2, p3);
            surface_area += std::sqrt(tri.squared_area());
        }
    }
    
    std::cout << "Probe radius: " << probe_radius << " Å" << std::endl;
    std::cout << "Accessible surface area: " << surface_area << " Å²" << std::endl;
    
    // 计算连通分量（分子口袋）
    std::size_t n_components = alpha_shape.number_of_solid_components();
    std::cout << "Number of cavities: " << n_components - 1 << std::endl;
}
 
int main(int argc, char* argv[]) {
    std::vector<Point> atoms;
    
    // 从文件读取原子坐标（PDB格式简化）
    std::ifstream in((argc > 1) ? argv[1] : "molecule.xyz");
    double x, y, z;
    while (in >> x >> y >> z) {
        atoms.emplace_back(x, y, z);
    }
    
    std::cout << "Loaded " << atoms.size() << " atoms" << std::endl;
    
    // 使用不同探针半径计算表面积
    std::vector<double> probe_radii = {1.0, 1.4, 2.0, 3.0};  // 单位：埃
    
    for (double radius : probe_radii) {
        compute_molecular_surface(atoms, radius);
        std::cout << "------------------------" << std::endl;
    }
    
    return 0;
}

8. 复杂度分析

8.1 时间复杂度

操作	2D	3D
构造	O(n log n)	O(n log n)
初始化区间映射	O(n log n)	O(n log n)
设置 alpha 值	O(1)	O(1)
分类查询	O(1)	O(1)
遍历 alpha shape 边/面	O(k)	O(k)
查找最优 alpha	O(n log n)	O(n log n)
计算连通分量	O(n)	O(n)

8.2 空间复杂度

数据结构	2D	3D
Delaunay 三角剖分	O(n)	O(n)
区间映射表	O(n)	O(n)
Alpha 谱	O(n)	O(n)
总计	O(n)	O(n)

9. 实际应用案例

9.1 案例一：点云形状分析

问题：从激光扫描获取的建筑点云，需要提取不同细节层次的轮廓。

解决方案：

使用小 alpha 提取精细结构（窗户、装饰）
使用中等 alpha 提取主要结构（墙体、屋顶）
使用大 alpha 提取整体外形

代码示例：

void analyze_building_point_cloud(const std::vector<Point_3>& points) {
    Alpha_shape_3 alpha(points.begin(), points.end());
    
    // 精细结构
    alpha.set_alpha(0.5);
    auto fine_components = alpha.number_of_solid_components();
    std::cout << "Fine structure components: " << fine_components << std::endl;
    
    // 中等结构
    alpha.set_alpha(5.0);
    auto medium_components = alpha.number_of_solid_components();
    std::cout << "Medium structure components: " << medium_components << std::endl;
    
    // 整体结构
    alpha.set_alpha(50.0);
    auto coarse_components = alpha.number_of_solid_components();
    std::cout << "Coarse structure components: " << coarse_components << std::endl;
}

9.2 案例二：生物分子 pocket 检测

问题：寻找蛋白质表面的小分子可结合口袋。

方法：

使用小探针（1.0 Å）检测所有表面凹陷
使用大探针（3.0 Å）过滤掉太小的口袋
口袋 = 大探针无法进入但小探针可以进入的区域

9.3 案例三：传感器网络覆盖分析

问题：分析无线传感器网络的覆盖范围和空洞。

方法：

每个传感器位置作为一个点
Alpha shape 显示覆盖区域的连通性
空洞对应于未被覆盖的区域

10. 参数调优指导

10.1 如何选择 Alpha 值

自动选择方法：

// 方法1：使用最优 alpha
auto optimal = alpha_shape.find_optimal_alpha(1);
alpha_shape.set_alpha(*optimal);
 
// 方法2：基于点云密度
FT average_spacing = CGAL::compute_average_spacing(points, 6);
FT alpha = 4.0 * average_spacing * average_spacing;  // 2倍间距的平方
alpha_shape.set_alpha(alpha);
 
// 方法3：基于百分位数
std::vector<FT> alpha_values;
for (auto it = alpha_shape.alpha_begin(); 
     it \!= alpha_shape.alpha_end(); ++it) {
    alpha_values.push_back(*it);
}
std::sort(alpha_values.begin(), alpha_values.end());
FT median_alpha = alpha_values[alpha_values.size() / 2];
alpha_shape.set_alpha(median_alpha);

10.2 模式选择指南

应用场景	推荐模式	原因
形状重建	REGULARIZED	避免孤立的边和点
拓扑分析	GENERAL	保留所有拓扑信息
分子建模	REGULARIZED	更自然的表面
网络分析	GENERAL	保留所有连接信息

10.3 常见问题与解决

问题1：结果有太多孔洞

原因：alpha 值太小
解决：逐步增大 alpha，或使用 find_optimal_alpha

问题2：重要特征丢失

原因：alpha 值太大
解决：减小 alpha，或使用多尺度分析

问题3：数值不稳定

原因：点集有极端分布
解决：使用 ExactAlphaComparisonTag = Tag_true

11. 关键文件位置

文件路径	说明
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_2/include/CGAL/Alpha_shape_2.h`	2D Alpha Shape 主类
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_2/include/CGAL/Alpha_shape_face_base_2.h`	2D 面基类
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_2/include/CGAL/Alpha_shape_vertex_base_2.h`	2D 顶点基类
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_3/include/CGAL/Alpha_shape_3.h`	3D Alpha Shape 主类
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_3/include/CGAL/Alpha_shape_cell_base_3.h`	3D 单元基类
`/home/chunibyo/workspace/osc/cgal/Alpha_shapes_3/include/CGAL/Alpha_shape_vertex_base_3.h`	3D 顶点基类

12. 注意事项

模式选择：
- GENERAL 模式保留所有维度的单形，包括孤立的顶点和边
- REGULARIZED 模式只保留与最高维单元相邻的单形
数值稳定性：
- 使用 Exact_predicates_inexact_constructions_kernel (Epick)
- 对于需要精确比较的场景，使用 Tag_true 作为第二个模板参数
Gabriel 单形：
- 只在 GENERAL 模式下计算
- Gabriel 边/面的 alpha_min 是其 Gabriel 球半径
性能优化：
- 顶点/边/面列表有缓存机制
- 修改 alpha 值后需要重新计算列表

参考文献

Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of a set of points in the plane. IEEE Transactions on Information Theory, 29(4), 551-559.
Edelsbrunner, H., & Muecke, E. P. (1994). Three-dimensional alpha shapes. ACM Transactions on Graphics, 13(1), 43-72.
CGAL Documentation: 2D Alpha Shapes. https://doc.cgal.org/latest/Alpha_shapes_2/
CGAL Documentation: 3D Alpha Shapes. https://doc.cgal.org/latest/Alpha_shapes_3/

Chunibyo

Explorer

10.2 Alpha Shapes

10.2 Alpha Shapes §

0. 算法动机：为什么需要 Alpha Shapes？ §

现实问题引入 §

传统方法的局限 §

如果不使用 Alpha Shapes §

1. 理论基础 §

1.1 Alpha Shapes 的定义 §

1.2 与相关概念的关系 §

2. 直观理解：融冰与雕刻 §

2.1 “融冰”类比 §

2.2 “雕刻”类比 §

2.3 数学直觉 §

3. Alpha Shapes 与相关概念的关系图 §

3.1 概念层次图 §

3.2 不同 Alpha 值效果对比 §

3.3 与 Gabriel 图的关系 §

4. CGAL Alpha Shapes 架构 §

4.1 模块组织 §

4.2 类层次结构 §

5. Alpha Shape 2D 实现详解 §

5.1 核心类定义 §

5.2 区间计算 §

5.3 区间表示详解 §

5.4 分类算法 §

6. Alpha Shape 3D 实现详解 §

6.1 核心类定义 §

6.2 初始化流程 §

6.3 边状态计算 §

7. 使用示例 §

7.1 基本 2D Alpha Shape §

7.2 3D Alpha Shape §

7.3 使用 Weighted Alpha Shapes §

7.4 遍历 Alpha 谱 §

7.5 实际应用：分子表面积计算 §

8. 复杂度分析 §

8.1 时间复杂度 §

8.2 空间复杂度 §

9. 实际应用案例 §

9.1 案例一：点云形状分析 §

9.2 案例二：生物分子 pocket 检测 §

9.3 案例三：传感器网络覆盖分析 §

10. 参数调优指导 §

10.1 如何选择 Alpha 值 §

10.2 模式选择指南 §

10.3 常见问题与解决 §

11. 关键文件位置 §

12. 注意事项 §

参考文献 §

Graph View

Table of Contents

Backlinks

10.2 Alpha Shapes

0. 算法动机：为什么需要 Alpha Shapes？

现实问题引入

传统方法的局限

如果不使用 Alpha Shapes

1. 理论基础

1.1 Alpha Shapes 的定义

1.2 与相关概念的关系

2. 直观理解：融冰与雕刻

2.1 “融冰”类比

2.2 “雕刻”类比

2.3 数学直觉

3. Alpha Shapes 与相关概念的关系图

3.1 概念层次图

3.2 不同 Alpha 值效果对比

3.3 与 Gabriel 图的关系

4. CGAL Alpha Shapes 架构

4.1 模块组织

4.2 类层次结构

5. Alpha Shape 2D 实现详解

5.1 核心类定义

5.2 区间计算

5.3 区间表示详解

5.4 分类算法

6. Alpha Shape 3D 实现详解

6.1 核心类定义

6.2 初始化流程

6.3 边状态计算

7. 使用示例

7.1 基本 2D Alpha Shape

7.2 3D Alpha Shape

7.3 使用 Weighted Alpha Shapes

7.4 遍历 Alpha 谱

7.5 实际应用：分子表面积计算

8. 复杂度分析

8.1 时间复杂度

8.2 空间复杂度

9. 实际应用案例

9.1 案例一：点云形状分析

9.2 案例二：生物分子 pocket 检测

9.3 案例三：传感器网络覆盖分析

10. 参数调优指导

10.1 如何选择 Alpha 值

10.2 模式选择指南

10.3 常见问题与解决

11. 关键文件位置

12. 注意事项

参考文献