[Geometry] Alpha Shapes - 原理及我的实现

我的新博客:http://ryuzhihao.cc/


    已经快一个月没有写新文章了。其实这段时间我写了很多小程序,但是一直在忙保研的事情就没有时间写文章。一次无尽的9月(at涼宮ハルヒ)。现在终于一切事情都结束了,重新开始填坑。

    

    最近出于一些其他的用途,整理并实现了Alpha shapes,所以拿出来分享一下。首先,这篇博客以及我的程序,是根据Francois Belair写的文章而完成的,原文链接:Alpha Shape Applet

    这里贴一下在下面提供的一个测试程序的下载链接:

            http://download.csdn.net/download/mahabharata_/10001135

一、什么是Alpha Shapes?

     Alpha Shapes是从离散的空间点集(point sets)中抽象出其直观形状的一种方法,简言之,从一堆无序的点中获取大致的轮廓。

         首先,还是先贴出我的实现结果。

(1)我的某个使用Alpha Shapes形状检测的程序: 

                           

(二)专门为Alpha Shapes 设计的一个测试程序

         下图程序的下载链接在文章开头已经给出。

                          

二、Alpha Shapes算法流程

(一)单阈值Alpha Shapes算法的轮廓提取

       这是最简单的一种Alpha Shapes的实现。在一个有限离散点集S中,有n个点构成,这n个点可以组成n*(n-1)条线段,我们可以通过如下的方法判断哪条线段是在边界线上的线段:在点集S内,过任意两点P1、P2绘制半径为alpha的圆(显然,在给定半径时,过确定的两点的圆应该有两个),如果这个圆内没有其他点,则认为点P1、P2是边界点,其连线P1P2即为边界线段。

       如果粗暴的去求解,其时间复杂度为O(n^3)。但由于对称性,这种方法的复杂度可降低至O(n^2*log2n)。当然中间还可以通过一些其他的约束条件去减少计算次数,具体的计算流程如下:

       第1步:遍历每条边P1-P2。

       第2步:如果边P1-P2的长度大于直径2*Alpha,则跳过(因为无法找到过这两个点的圆)。

       第3步:根据简单的几何知识,求出两个圆的圆心C1、C2。

         (3.1)  边P1-P2的方向向量v1 =p1-p2,且设v1为(x,y),同时边p1p2的中点坐标mid=0.5*(p1+p2)。

         (3.2)  我们可以根据向量垂直时点积为0的性质,求出边P1-P2的垂直向量v2(设v2坐标为(a,b),则ax+by=0)。

         (3.3)  设边P1-P2的长度为len,求出圆心到边P1-P2的距离dist = sqrt( alpha^2 -  0.25*len^2);

         (3.4)   圆心C1 =mid+dist*v2;圆心C2 =mid-dist*v2;

       第4步:两个圆中有一个圆内部不包含点集S中的其他任何点,则边P1-P2即为边界。

       实现结果如下图:分别是轮廓的检测结果、Alpha检测圆:

                                            


                                           

       求解alpha shapes轮廓的代码如下:

struct Edge
{
    vec2 a;
    vec2 b;
}

vector<vec2> m_points;  // 点集S
vector<Edge> m_edges;   // 边


void PointSet::AlphaShapes()
{
    m_edges.clear();

    for(int i=0; i<m_points.size(); i++)
    {
        // k从i+1开始,减少重复计算
        for(int k=i+1; k<m_points.size(); k++)
        {
            // 跳过距离大于直径的情况
            if(m_points[i].distanceToPoint(m_points[k])>2*m_radius)
                continue;

            // 两个圆心
            QVector2D c1,c2;

            // 线段中点
            QVector2D center = 0.5*(m_points[i]+m_points[k]);

            // 方向向量 d = (x,y)
            QVector2D dir = m_points[i] - m_points[k];

            // 垂直向量 n = (a,b)  a*dir.x+b*dir.y = 0; a = -(b*dir.y/dir.x)
            QVector2D normal;
            normal.setY(5);  // 因为未知数有两个,随便给y附一个值5。 

            if(dir.x() != 0)  
                normal.setX(-(normal.y()*dir.y())/dir.x());
            else     // 如果方向平行于y轴
            {
                normal.setX(1);
                normal.setY(0);
            }
            normal.normalize();   // 法向量单位化

            float len = sqrt(m_radius*m_radius-(0.25*dir.length()*dir.length()));
            c1 = center + len*normal;
            c2 = center - len*normal;

            // b1、b2记录是否在圆C1、C2中找到其他点。
            bool b1 = false, b2 = false;
            for(int m=0; m<m_points.size(); m++)
            {
                if( m == i || m == k)
                    continue;

                if(b1!=true && m_points[m].distanceToPoint(c1) < m_radius)
                    b1 = true;
                if(b2!=true && m_points[m].distanceToPoint(c2) < m_radius)
                    b2 = true;

                // 如果都有内部点,不必再继续检查了
                if(b1 == true && b2 == true)
                    break;
            }

            Edge edge;
            if(b1!=true || b2!=true)
            {
                edge.a = m_points[i];
                edge.b = m_points[k];

                m_edges.push_back(edge);
            }

        }
    }

}

(二)改进的单阈值Alpha Shapes轮廓提取:

           为了减少计算量。我们可以先采取先建立三角网,再从三角网的外边界开始判断,其流程如下:

           第1步:根据点集S建立Delaunay三角网。

           第2步:若三角形中某条边的长度大于2*alpha,则删除该三角形(不符合条件)。

           第3步:对三角形的每条边进行判断:若过某条边的两点且半径为alpha的圆包含其他点,则删除该三角形。如下图所示,当前要判断的三角形为△ABC,要判断的边为边A-B,因为半径为alpha且过AB两点的圆包含点C,则从三角网中删除△ABC。

            

           第4步:在第2-3步提到的不符合Alpha shapes要求的三角形后所得到的三角网上求出三角网的边缘。该边缘即为点集S的边缘线。



©️2020 CSDN 皮肤主题: 撸撸猫 设计师:设计师小姐姐 返回首页