描述

使用阿尔法形状算法的情况是这样的：当我们拥有一些二维或三维的点云时，点云执行过聚类操作后会呈现物体本身的形状特征。我们希望通过一些简化的方式描述出这个形状。这时阿尔法形状算法会是一个很好的选择。

1.基本概念

点云是指在三维空间中的一组点，通常用于表示三维物体的形状。而在计算几何和计算机图形学中，点云的凸包和凹包是非常重要的概念，它们被广泛应用于三维建模、虚拟现实、机器人技术等领域。

凸包是指包围点云的最小凸多面体，即最小的凸多面体，使得点云中的所有点都在多面体的内部或表面上。凸包可以用于计算点云的表面积和体积，以及进行碰撞检测和物体识别等应用。

凹包（也称为 alpha shape）是指包围点云的最小凹多面体，即最小的多面体，使得点云中的所有点都在多面体的内部或表面上。与凸包不同，凹包可以包含凸部分，因此可以更好地表示一些复杂的形状。凹包的形状取决于一个参数 alpha，该参数控制着多面体的平滑程度。较小的 alpha 值将导致多面体较为精细地描述点云的形状，而较大的 alpha 值将导致多面体趋向于凸包。

2.原文内容

阿尔法形状算法最初是由 Herbert Edelsbrunner 和 Ernst P. Mücke 在1994年的一篇名为“Three-Dimensional Alpha Shapes”的论文中介绍的。

文章链接：https://dl.acm.org/doi/10.1145/174462.156635

文章的两个截图如下

意思很明显。每组图都是由6个图构成的，最后一个图就是点云的真实形状。从第一个图到第六个图，依次是α值从正无穷降低到0的表现。我们从图中很容易发现，α值较大时，多面体就是点云的凸包形状。随着α值的减少，多面体越来越趋近于点云形状，当α值较为合适时，多面体可以较完美的表征点云的形状（图4）。当α值很小时，多面体表征形状的能力会重新下降。

3.2D算法

2D阿尔法形状算法的操作对象，是二维点云。

算法的基本思想是这样的：用一个半径为α的圆，绕点集S滚动。当α值取值合适时，这个圆是无法从点集的边界点滚进点集内部的。能够组织圆滚进内部的点，就是这个点集的边界点，这些边界点相邻两点构成的就是边缘轮廓线。在点集S内任取两个点p1和p2，过这两点绘制半径为α的圆，若圆内没有其他点，则可以认为p1p2是边界线。

无法实现的算法步骤：显然一个简单的算法就是，对于一个大小为n的点集，我们可以遍历n(n-1)条两两点构成的边，判断这些边的长度是否大于圆直径2α。大于直径时，这两点构成的边不是边界线；如果小于直径时，再去检查其余点是否存在于过该边的两个圆中。

https://blog.csdn.net/qq_32867925/article/details/113461263
这篇文章写的就是这样的实现过程。这样做，算法复杂度粗略就是一个o(n^3)。这在实际工程中是无法承受的开销，唯一的好处就是理解阿尔法算法容易一些。

在计算点距的过程中实际有很多可以优化的部分，一个新的知识——Delaunay三角网（德劳内三角网），我们利用它就可以大大简化计算量了。

3.1Delaunay三角网

关于Delaunay三角网，这篇博客写的非常不错，也有很清晰的图片，先贴在这里可以参考https://blog.csdn.net/qq_49838656/article/details/119641392
Delaunay三角网是一系列相连的但不重叠的三角形的集合, 而且这些三角形的外接圆不包含这个面域的其他任何点。Delaunay三角网还有一个对偶形式，称为 Voronoi 图（维诺图）。维诺图上包含了一系列多边形和一系列点，其中维诺图上的点就是Delaunay三角网中全部三角形的外接圆圆心，维诺图上的多边形的边是任意两个点连线的垂直平分线。Delaunay 三角剖分和 Voronoi 图是一一对应的，它们之间的对偶关系可以用于计算几何和计算机图形学中的许多问题。

看一下下面这张别人给出的图，就很好理解了。

上图就是一个Delaunay三角网和它的维诺图。

Delaunay三角网有很多优秀的数学性质，这篇博客介绍了很多特性, https://blog.csdn.net/qixun7099/article/details/100739039。这些性质使得 Delaunay 三角剖分在计算几何、计算机图形学、地理信息系统等领域得到广泛应用。

其中以空外接圆特性最为重要：三角网中任意一个三角形的外接圆，圆内不会再包含三角网中其他三角形的顶点。2D阿尔法形状算法要用的也是这个特性。

3.2利用Delaunay三角网的2D算法

上面说过Delaunay三角网具备空外接圆特性，我们要利用这个特性，来降低计算点云中点距的计算复杂度。针对点数为n的点云，我们希望构建一个属于它的Delaunay三角网，这样我们在计算阿尔法形状算法的边界线时，只用计算对应三角形的边长是否满足要求即可。

接下来要分两个部分来介绍算法。一是：为什么要利用Delaunay三角网来对点云预处理，它的好处是什么？二是：利用Delaunay三角网的阿尔法形状算法如何实现？

3.2.1. Delaunay三角网预处理点云

首先增量法是构建Delaunay三角网最常用的算法之一，在原有的三角网中插入新的点，是不需要每次都遍历整个三角网查找不符合三角剖分的三角形，并进行优化的。因为除了新插入的点外的三角网都是符合三角剖分的。除此之外，还有分割合并算法、递归生长算法等。

基于Bowyer-Watson算法就是增量化构建Delaunay三角网的经典算法，是一种逐点插入的算法，它从一个空三角形开始，逐步将点插入到三角形中，同时保持 Delaunay 三角剖分的性质。与普通的增量法不同，Bowyer-Watson 算法使用了一种更高效的数据结构，称为增量三角形网格（incremental triangulation），可以在 O(log n) 的时间内完成点的插入和删除操作。

sorry，这里埋个坑，构建Delaunay三角网的算法有很多，无法在这篇里详细介绍。以后找个时间专门学习一下。这里只能先略过了……

构建了Delaunay三角网，就相当于我们对给定的无序点云集S，生成了一个有序的数据结构。利用Delaunay三角网的空外接圆特性和这个有序结构，我们可以在计算中，减去那些不相关点的计算过程。针对阿尔法形状算法，利用Delaunay三角网来对点云预处理的好处，就是当我们针对点云集S构建了属于它的三角网后。对于半径长为α的圆，能否滚进点云中就变成了简单的两个判断逻辑：一是：三角形的边长不可以超过直径2α，二是过三角形两点半径为α的圆，不可以包含三角网的其他顶点。因此，预处理构建三角网后，不符合条件的三角形会随着计算逐步删除，阿尔法形状算法的计算量将大大降低。

3.2.2. 2D阿尔法形状算法步骤

（1）根据点集S建立它的Delaunay三角网

（2）遍历三角网，如果该三角形中任意一条边的长度大于2α，则删除该三角形。意味着半径为α的圆是可以滚进这个三角形内的。

（3）对剩余三角网的每条边进行判断：若过某条边的两点且半径为α的圆包含其他点，则删除该三角形。意味着半径为α的圆虽然没通过该边完全滚进点云形状内，但滚进去的部分仍然侵略了形状的部分点。

（4）求出剩余三角网的边缘边，这些边构成了点集S的边缘线。

3.2.3. C++代码实现

针对opencv中的Delaunay三角剖分计算，这篇文章写的很好：https://blog.csdn.net/GoodLi199309/article/details/80191677
我写了一个2D阿尔法形状算法的程序，贴在下面

#include <iostream>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <algorithm>

#include <opencv2/opencv.hpp>
#include <opencv2/core.hpp>

#define random(a,b) (rand()%(b-a)+a)
#define alpha 100

std::vector<cv::Point2f> getCircleCenter(cv::Point2f p1, cv::Point2f p2, float radius) {
    std::vector<cv::Point2f> result;

    if (p1.x == p2.x) {
        float y = (p1.y + p2.y)/2;
        float y_ = (p1.y - p2.y)/2;
        float x1 = p1.x - sqrt(radius*radius - y_*y_);
        float x2 = p1.x + sqrt(radius*radius - y_*y_);
        result.push_back(cv::Point2f(x1, y));
        result.push_back(cv::Point2f(x2, y));
        return result;
    }

    float c1 = (p2.x*p2.x - p1.x*p1.x + p2.y*p2.y - p1.y*p1.y) / (2 *(p2.x - p1.x));  
    float c2 = (p2.y - p1.y) / (p2.x - p1.x);  //斜率
    float A = (c2*c2 + 1);  
    float B = (2 * p1.x*c2 - 2 * c1*c2 - 2 * p1.y);  
    float C = p1.x*p1.x - 2 * p1.x*c1 + c1*c1 + p1.y*p1.y - radius*radius;  
    float y1 = (-B + sqrt(B*B - 4 * A*C)) / (2 * A);
    float y2 = (-B - sqrt(B*B - 4 * A*C)) / (2 * A);
    float x1 = c1 - c2 * y1;
    float x2 = c1 - c2 * y2;


    result.push_back(cv::Point2f(x1, y1));
    result.push_back(cv::Point2f(x2, y2));

    return result;
}

std::vector<cv::Point2f> alpha_shape(const std::vector<cv::Point2f>& points) {

    float minX = points[0].x, maxX = points[0].x;
    float minY = points[0].y, maxY = points[0].y;
    for (const auto& point : points) {
        if (point.x < minX) {
            minX = point.x;
        }
        if (point.x > maxX) {
            maxX = point.x;
        }
        if (point.y < minY) {
            minY = point.y;
        }
        if (point.y > maxY) {
            maxY = point.y;
        }
    }
    std::cout<<"点云 x_min:"<<minX<<" x_max:"<<maxX<<" y_min:"<<minY<<" y_max:"<<maxY<<std::endl;

    int subdiv_size_u = maxX - minX + 30;
    int subdiv_size_v = maxY - minY + 30;
    // Create a Delaunay triangulation of the points
    cv::Subdiv2D subdiv(cv::Rect(0, 0, subdiv_size_u, subdiv_size_v));
    for (int i = 0; i < points.size(); ++i) {
        subdiv.insert(cv::Point2f(points[i].x-minX+10, points[i].y-minY+10));
    }

    // Get the edges of the Delaunay triangulation
    std::vector<cv::Vec6f> triangleList;
    subdiv.getTriangleList(triangleList);

    std::vector<cv::Vec6f> triangleList_result;
    // Iterate over the edges of the Delaunay triangulation
    for (const auto& triangle : triangleList) {
        // Get the vertices of the triangle
        cv::Point2f p1(triangle[0], triangle[1]);
        cv::Point2f p2(triangle[2], triangle[3]);
        cv::Point2f p3(triangle[4], triangle[5]);

        if (cv::norm(p1 - p2) > 2 * alpha 
         || cv::norm(p1 - p3) > 2 * alpha 
         || cv::norm(p2 - p3) > 2 * alpha) {
            continue;
        }
        triangleList_result.push_back(triangle);
    }

    // Create a vector to store the edges of the alpha shape
    std::vector<cv::Point2f> alphaEdges;

    // Iterate over the edges of the Delaunay triangulation
    for (const auto& triangle : triangleList_result) {
        cv::Point2f p1(triangle[0], triangle[1]);
        cv::Point2f p2(triangle[2], triangle[3]);
        cv::Point2f p3(triangle[4], triangle[5]);   

        std::vector<cv::Point2f> c1, c2, c3;

        std::vector<cv::Point2f> centers;
        std::vector<int> valids;

        int label = 1;
        c1 = getCircleCenter(p1, p2, alpha);
        label = 1;
        for (int i = 0; i < c1.size(); i++) {
            centers.push_back(c1[i]);
            valids.push_back(label);
            label++;
        }

        c2 = getCircleCenter(p1, p3, alpha);
        label = 3;
        for (int i = 0; i < c2.size(); i++) {
            centers.push_back(c2[i]);
            valids.push_back(label);
            label++;
        }

        c3 = getCircleCenter(p2, p3, alpha);
        label = 5;
        for (int i = 0; i < c3.size(); i++) {
            centers.push_back(c3[i]);
            valids.push_back(label);
            label++;
        }

        if (centers.size() == 0) {
            continue;
        }

        for (const auto& otherTriangle : triangleList_result) {
            if (triangle == otherTriangle) {
                continue;
            }

            // Get the vertices of the other triangle
            cv::Point2f op1(otherTriangle[0], otherTriangle[1]);
            cv::Point2f op2(otherTriangle[2], otherTriangle[3]);
            cv::Point2f op3(otherTriangle[4], otherTriangle[5]);

            for (int i = 0; i < centers.size(); i++) {
                if (cv::norm(centers[i] - op1) < alpha
                 || cv::norm(centers[i] - op2) < alpha 
                 || cv::norm(centers[i] - op3) < alpha) {
                    valids[i] = 0;
                }
            }
        }

        for (int i = 0; i < centers.size(); i++) {
            // std::cout<<p1<<" "<<p2<<" "<<p3<<" "<<i<<std::endl;
            if (valids[i] == 1 || valids[i] == 2) {
                alphaEdges.push_back(p1);
                alphaEdges.push_back(p2);
                break;
            } else if (valids[i] == 3 || valids[i] == 4) {
                alphaEdges.push_back(p1);
                alphaEdges.push_back(p3);
                break;
            } else if (valids[i] == 5 || valids[i] == 6) {
                alphaEdges.push_back(p2);
                alphaEdges.push_back(p3);
                break;
            } else {
                continue;
            }
        }
    }

    // Remove duplicate edges
    std::sort(alphaEdges.begin(), alphaEdges.end(), [&](const cv::Point2f& p1, const cv::Point2f& p2) {
        return p1.x < p2.x;
    });
    alphaEdges.erase(std::unique(alphaEdges.begin(), alphaEdges.end()), alphaEdges.end());

    for (int i = 0; i < alphaEdges.size(); ++i) {
        alphaEdges[i].x = alphaEdges[i].x + minX - 10;
        alphaEdges[i].y = alphaEdges[i].y + minY - 10;
    }

    return alphaEdges;
}

int main(int argc, char** argv) {
    std::vector<cv::Point2f> points;

    // 输入数据1：随机生成一组float类型的点云
    srand((int)time(0));
    int number = 100;
    for (int i = 0; i < number; i++) {
        float x = float(random(-30, 30));
        float x_ = float(random(1, 99));
        float y = float(random(-30, 30));
        float y_ = float(random(1, 99));
        cv::Point2f p(x+x_*0.01, y+y_*0.01);
        // std::cout<<p.x<<" "<<p.y<<std::endl;
        points.push_back(p);
    }

    // 输入数据2：一组矩形规则点云，可用来检查效果
    // for (int i = 0; i < 20; i++) {
    //     for (int j = 0; j < 20; j++) {
    //         cv::Point2f p(-10+10*i, -10+10*j);
    //         points.push_back(p);
    //     }
    // }

    std::vector<cv::Point2f> points_edge = alpha_shape(points);

    cv::Mat image = cv::Mat::zeros(600, 600, CV_8UC3);
    for (int i = 0; i < points.size(); i++){
        int u = 300 + points[i].x;
        int v = 300 + points[i].y;
        cv::circle(image, cv::Point(u, v), 1, cv::Scalar(255, 0, 0), -1);
    }
    std::cout<<"阿尔法边缘点个数: "<<points_edge.size()<<std::endl;
    for (int i = 0; i < points_edge.size(); i++){
        int u = 300 + points_edge[i].x;
        int v = 300 + points_edge[i].y;
        cv::circle(image, cv::Point(u, v), 3, cv::Scalar(0, 255, 0), -1);
    }


    cv::imshow("alpha", image);
    cv::waitKey(0);

    return 0;
}

代码的效果展示：
当点云是规则形状时

当点云随机时

代码中的alpha值被设置为了100，所以会出现边界点。可以自行调整它的值，值越小，边界点越少。代码的效果我只是简单验证了，原理无误。运算时间待优化，有问题随时联系。

总结

3D阿尔法形状算法的操作对象，是三维点云。与2D形状算法类似，下一篇文章来介绍。