C++ 基于凸包的Delaunay三角网生成算法

Delaunay三角网，写了用半天，调试BUG用了2天……醉了。

基本思路比较简单，但效率并不是很快。

1. 先生成一个凸包；

2. 只考虑凸包上的点，将凸包环切，生成一个三角网，暂时不考虑Delaunay三角网各种规则。将生成的三角形放进三角形集合 Triangles 中；

3.将其它非凸包上的点全都插入。每插入一个点 ptA，都要判断被插入的点存在于 Triangles 集合的哪个三角形（trianA）之内，并将 ptA 与此三角形的三个点进行连接，删除 trianA，并将新生成的三角形加入到集合 Triangles 中。初始三角网生成结束；

　　3.1. 若 ptA 处在三角形DEF的DE边上，那么只连接点F与ptA ；

4.遍历三角形集合 Triangles（曾考虑用邻接矩阵，但是使用矩阵的复杂度反而会更高），每遍历到一个三角形DEF，都要遍历DEF的三条边DE,EF, FD，并分别寻找另外一个三角形，此三角形与DEF存在一个公共边；

5.使用LOP（Local Optimization Procedure： 局部优化处理）处理步骤4返回的两个三角形 与DEG。

　　5.1. 做三角形DEF的外接圆cirDEF，如果点G在圆cirDEF之内，则从Triangles 集合中删除三角形DEF与DEG，并加入两个新三角形：FGD与FGE。实际上，是将两个三角形组成的四边形的对角线对调了。

6.使用新生成的两个三角形FGD与FGE执行步骤5，递归。递归出口为，对三角形进行LOP处理时，没有出现5.1的情况。

Note:

1.暂时存在的问题：

　　a. 当点数过多时，生成过程中会出现“栈溢出”情况。测试时，2000个点之内是成功的。此情况是由于三角形过多导致的递归过深，需要重新组织算法结构才能避免。

　　b. 生成时间太长。可采取的办法1是完全采用不同的生成算法，比如分块等等。其次是由于寻找三角形的复杂度过高，达到了O(n)，利用空间换时间的办法能接近O(1)，只不过多一个成员来存储三角形的拓扑结构，这个比较简单但优化效果有限。

2. 一旦遇到了一些算法上的问题，那么搬出数学知识往往非常有效（尽管我的数学水平一望见底）。使用矩阵的方式能够让算法的实现更加清晰，但复杂度有可能会升高。

*下面使用了https://www.cnblogs.com/rkexy/p/9768475.html其中的代码。*

File: Delaunay.h

 #pragma once

 #include "ConvexHull.h"
 #include "Triangle.h"
 #include "Circle.cpp"

 class TrianIndex  final
 {
 public:
     TrianIndex()
     {
         _isVaild = false;
     }

     ~TrianIndex()
     {
         _isVaild = false;
     }

     TrianIndex(unsigned int iA, unsigned int iB, unsigned int iC)
     {
         Init(iA,  iB,  iC);
     }

     TrianIndex(std::array<unsigned int, > pts) :
         TrianIndex(pts[], pts[], pts[])
     {
     }

     TrianIndex(const TrianIndex& other)
     {
         this->_isVaild = other._isVaild;
         this->_ptIdx = other._ptIdx;
     }

     unsigned int& Get(int i)
     {
         if (i <  || i >  || !_isVaild)
         {
             ErrorThrow("Error Triangle Point Index[0, 2] Or Invaild Triangle: " + std::to_string(i));
             return _ptIdx[];
         }

         return _ptIdx[i];
     }

     unsigned int& operator[](int i)
     {
         return Get(i);
     }

     const unsigned int Get(int i) const
     {
         TrianIndex* pThis = const_cast<TrianIndex*>(this);
         return pThis->Get(i);
     }

     const unsigned int operator[](int i) const
     {
         TrianIndex* pThis = const_cast<TrianIndex*>(this);
         return (*pThis)[i];
     }

     TrianIndex& operator=(const TrianIndex& other)
     {
         this->_isVaild = other._isVaild;
         this->_ptIdx = other._ptIdx;
         return *this;
     }

     bool IsVaild() const
     {
         return _isVaild;
     }

     void SetVaild(bool isVaild)
     {
         _isVaild = isVaild;
     }

     void Init(unsigned int iA, unsigned int iB, unsigned int iC)
     {
         _ptIdx[] = iA;
         _ptIdx[] = iB;
         _ptIdx[] = iC;
         _isVaild = true;
     }

 private:
     bool _isVaild;
     std::array<unsigned int, > _ptIdx;
 };

 //CHECK IT: It Has Been Abandoned.
 //
 class TrianMatrix final
 {
 public:
     enum IsConnected : bool
     {
         NotConnected = false,
         Connected = true
     };

     TrianMatrix():
         TrianMatrix()
     {
     }

     TrianMatrix(size_t ptCount)
     {
         _sign.clear();
         _sign.resize(ptCount);
         for (unsigned int x = ; x < ptCount; ++x)
         {
             _sign[x].resize(ptCount);
             for (unsigned int y = ; y < ptCount; ++y)
             {
                 _sign[x][y] = IsConnected::NotConnected;
             }
         }
     }

     ~TrianMatrix()
     {
         _sign.clear();
     }

 public:
     size_t PointCount() const
     {
         return _sign.size();
     }

     bool IsEmpty() const
     {
         return _sign.empty();
     }

     //
     bool GetTrianglesByPoint(unsigned int ptIdx,
         std::vector<TrianIndex>& __out resTrians) const
     {
         bool isFound = false;
         if (!CheckIndex(ptIdx))
         {
             return false;
         }

         //Time Complexity: O(n^2).
         for (unsigned int i = ; i < PointCount(); ++i)
         {
             isFound = GetTrianglesByLine(ptIdx, i, resTrians);
         }

         return isFound;
     }

     bool GetTrianglesByLine(unsigned int ptIdxA, unsigned int ptIdxB,
         std::vector<TrianIndex>& __out resTrians) const
     {
         bool isFound = false;
         short foundCount = ;
         if (!CheckIndex(ptIdxA, ptIdxB))
         {
             return false;
         }

         //Time Complexity: O(n).
         for (unsigned int i = ; i < PointCount(); ++i)
         {
             const std::vector<IsConnected>& pts = _sign[i];
             if (pts[ptIdxA] && pts[ptIdxB])
             {
                 TrianIndex trian(i, ptIdxA, ptIdxB);        //The another point at the first element of resTrians.
                 resTrians.push_back(trian);
                 isFound = true;

                 if (++foundCount == )
                 {
                     return true;
                 }
             }
         }
         return isFound;
     }

     inline bool Connect(unsigned int ptIdxA, unsigned int ptIdxB)
     {
         return SetConnectedState(ptIdxA, ptIdxB, IsConnected::Connected);
     }

     inline bool DisConnect(unsigned int ptIdxA, unsigned int ptIdxB)
     {
         return SetConnectedState(ptIdxA, ptIdxB, IsConnected::NotConnected);
     }

     inline bool SetConnectedState(unsigned int ptIdxA, unsigned int ptIdxB, IsConnected isConn)
     {
         if (!CheckIndex(ptIdxA, ptIdxB))
         {
             return false;
         }

         _sign[ptIdxA][ptIdxB] = isConn;
         _sign[ptIdxB][ptIdxA] = isConn;
         return true;
     }

     bool IsExist(const TrianIndex& trian) const
     {
         if (!CheckIndex(trian[], trian[], trian[]))
         {
             return false;
         }

         unsigned int ptA = trian[];
         unsigned int ptB = trian[];
         unsigned int ptC = trian[];

         return (_sign[ptA][ptB] && _sign[ptA][ptC] && _sign[ptB][ptC]);
     }

 private:
     bool CheckIndex(unsigned int idx, ...) const
     {
         std::vector<unsigned int> indeies;
         auto CheckFunc = [&indeies, this](unsigned int newIdx) -> bool
         {
             if (newIdx >= this->PointCount())
             {
                 return false;
             }

             for (std::vector<unsigned int>::const_iterator it = indeies.begin();
                 it != indeies.end(); ++it)
             {
                 if (*it == newIdx)
                 {
                     return true;
                 }
             }
             indeies.push_back(newIdx);
             return false;
         };

         va_list aps;
         va_start(aps, idx);
         unsigned int* pNextArg = nullptr;
         bool isVaild = false;
         do
         {
             pNextArg = va_arg(aps, unsigned int*);
             if (pNextArg)
             {
                 unsigned int tmp = *pNextArg;
                 if (!CheckFunc(tmp))
                 {
                     isVaild = true;
                 }
             }
         } while (pNextArg != NULL);

         va_end(aps);
         return isVaild;
     }

 private:
         std::vector<
             std::vector<IsConnected>
     > _sign;
 };

 typedef std::vector<TrianIndex> Triangles;
 class Delaunay final
 {
 public:
     Delaunay();
     ~Delaunay();

 public:
     void AddPoint(const Point2D<int>& __in newPoint);
     void AddRandomPoints(int count, int maxX, int maxY);

     void GetAllTriangles(Triangles& __out ts) const;
     void GetAllPoints(Points& __out pts) const;
     void Generate();

     void Clear();

 private:
     bool GetTrianIndexWithLine(int ptIdxA, int ptIdxB, std::array<size_t/*The Triangle Index Of m_triangles.*/, > & __out res, int ptIdxCur = -) const;
     void CutProcess(const Points& chPoints);
     void ConvexHullCut();
     void LOP(unsigned int trianIndex, bool* isProcessed);    //Local Optimization Procedure.

 private:
     Points m_points;
     Triangles m_triangles;
 };

File: Delaunay.cpp

 #include "Delaunay.h"
 #include <bitset>

 Delaunay::Delaunay()
 {
     m_points.Clear();
     m_triangles.clear();
 }

 Delaunay::~Delaunay()
 {
     Clear();
 }

 void Delaunay::AddPoint(const Point2D<int>& newPoint)
 {
     m_points.Add(newPoint);
 }

 void Delaunay::AddRandomPoints(int count, int maxX, int maxY)
 {
     std::vector<Point2D<int>> pointsArray;
     Point2D<int>::RandomPoints(count, maxX, maxY, pointsArray);

     m_points.Clear();
     for (int i = ; i < count; ++i)
     {
         const Point2D<int>& eachpt = pointsArray.at(i);
         PointNode pn;    pn.Init(i, eachpt);
         m_points.Add(pn);
     }
 }

 void Delaunay::ConvexHullCut()
 {
     if (m_points.Size() < )
     {
         ErrorThrow("Points count too less.");
         return;
     }
     m_triangles.clear();

     ConvexHull ch(m_points);
     Points chPoints;
     ch.Generate();
     ch.GetConvexHullPoints(chPoints);

     CutProcess(chPoints);
 }

 bool Delaunay::GetTrianIndexWithLine(int ptIdxA, int ptIdxB, std::array<size_t, > & __out res, int ptIdxCur) const
 {
     if (ptIdxA == ptIdxB || ptIdxCur == ptIdxA || ptIdxCur == ptIdxB)
     {
         ErrorThrow("Error Parameters! \n ptIdxA: " + std::to_string(ptIdxA) + ", ptIdxB: " + std::to_string(ptIdxB) + ", ptIdxCur: " + std::to_string(ptIdxCur));
         return false;
     }
     res.fill(-);
     uint8_t cnt = ;
     bool isFound = false;

     for(size_t i =; i < m_triangles.size(); ++i)
     {
         if (cnt > )
         {
             return isFound;
         }

         std::bitset<> foundSign = ;
         int8_t bit0Idx = -;
         const TrianIndex& tIdx = m_triangles[i];
         for (uint8_t i = ; i < ; i++)
         {
             if (tIdx[i] == ptIdxA || tIdx[i] == ptIdxB)
             {
                 foundSign = ( << i) | foundSign.to_ulong();
             }
             else
             {
                 bit0Idx = i;
             }
         }

         if ( == foundSign.count() && tIdx[bit0Idx] != ptIdxCur)
         {
             res[cnt++] = i;
             isFound = true;
         }
     }
     return isFound;
 }

 void Delaunay::GetAllTriangles(Triangles & ts) const
 {
     ts = m_triangles;
 }

 void Delaunay::GetAllPoints(Points & pts) const
 {
     pts = m_points;
 }

 void Delaunay::CutProcess(const Points & chPoints)
 {
     int chSize = static_cast<int>(chPoints.Size());

     //Next Convex Hull Points.
     Points nextPts;

     std::array<PointNode, > aTriangle;
     int i = ;
     int atSize = ;    /*atSize: aTriangle size.*/

     while (i <= chSize) /*Including the first point.*/
     {
         const PointNode*const cutpt = chPoints.GetBySequence(i % chSize);
         aTriangle[atSize++] = *cutpt;

         if (atSize == )
         {
             if (Line<int>::IsCollinear(, aTriangle[]._node, aTriangle[]._node, aTriangle[]._node))
             {
                 //Delete the second point.
                 aTriangle[] = aTriangle[];
                 --atSize;
                 ++i;
                 continue;
             }

             //Form a triangle.
             TrianIndex ttmp(aTriangle[]._index, aTriangle[]._index, aTriangle[]._index);
             m_triangles.push_back(ttmp);

             //The last serves as the starting point for the next triangle.
             aTriangle[] = aTriangle[];
             atSize = ;

             nextPts.Add(aTriangle[]);
         }
         ++i;
     }

     //Add the remaining points.
     //The first point in aTriangle has been added.
     for (int i = ; i < atSize; ++i)
     {
         nextPts.Add(aTriangle[i]);
     }

     if (nextPts.Size() >= )
     {
         CutProcess(nextPts);
     }
 }

 //Generate the delaunay.
 void Delaunay::Generate()
 {
     ConvexHullCut();

     for (int i =  ; i < m_points.Size(); ++i)
     {
         const Point2D<int>& eachPt = m_points(i);
         for(std::vector<TrianIndex>::iterator it = m_triangles.begin(); it != m_triangles.end(); )
         {
             const TrianIndex eachTri = *it;
             const Point2D<int>& ptA = m_points[eachTri[]];
             const Point2D<int>& ptB = m_points[eachTri[]];
             const Point2D<int>& ptC = m_points[eachTri[]];
             Triangle tmpTri(ptA, ptB, ptC);

             bool isBreak = false;
             Triangle::EnumPTRelation eptr = tmpTri.RelationPT(eachPt);
             switch (eptr)
             {
             case Triangle::EnumPTRelation::POINT_ON_VERTICE:
             case Triangle::EnumPTRelation::POINT_OUTSIDE:
                 ++it;
                 isBreak = false;
                 break;

             case Triangle::EnumPTRelation::POINT_ON_LINE:
                 for (unsigned int j = ; j < ; ++j)
                 {
                     unsigned int ptIdxA = j % ;
                     unsigned int ptIdxB = (j + ) % ;
                     unsigned int ptIdxC = (j + ) % ;
                     if (Line<int>::IsCollinear(, m_points[eachTri[ptIdxA]], m_points[eachTri[ptIdxB]], eachPt))
                     {
                         TrianIndex tmpAPC(eachTri[ptIdxC], i, eachTri[ptIdxA]);
                         TrianIndex tmpAPB(eachTri[ptIdxC], i, eachTri[ptIdxB]);
                         it = m_triangles.erase(it);    //Delete ABC.
                         m_triangles.push_back(tmpAPC);
                         m_triangles.push_back(tmpAPB);
                         break;
                     }
                 }
                 isBreak = false;
                 break;

             case Triangle::EnumPTRelation::POINT_INSIDE:
                 it = m_triangles.erase(it);
                 for (unsigned int j = ; j < ; ++j)
                 {
                     unsigned int ptIdx = j % ;
                     unsigned int ptIdxAnother = (j + ) % ;

                     TrianIndex tmp(eachTri[ptIdx], i, eachTri[ptIdxAnother]);
                     m_triangles.push_back(tmp);
                 }
                 isBreak = true;
                 break;
             }

             if(isBreak)
             {
                 break;
             }
         }
     }

     //LOP
     const size_t TrianglesSize = m_triangles.size();
     bool* isProcessed = nullptr;
     if (isProcessed == nullptr)
     {
         isProcessed = new bool[TrianglesSize];
         memset(isProcessed, false, TrianglesSize);
     }

     for (unsigned int i = ; i < TrianglesSize; ++i)
     {
         LOP(i, isProcessed);
     }

     if (isProcessed != nullptr)
     {
         delete[] isProcessed;
     }
 }

 void Delaunay::LOP(unsigned int trianIndex, bool* isProcessed = nullptr)
 {
     auto SetProcessedFunc = [&isProcessed](unsigned int idx, bool val)
     {
         if (nullptr != isProcessed)
         {
             isProcessed[idx] = val;
         }
     };

     if (trianIndex >= m_triangles.size())
     {
         ErrorThrow("The Index Of Triangles Out Of Range. trianIndex: " + std::to_string(trianIndex));
         return;
     }

     if (isProcessed != nullptr && isProcessed[trianIndex])
     {
         return;
     }

     TrianIndex* const curIdxT = &m_triangles[trianIndex];
     for (unsigned int i = ; i < ; ++i)
     {
         unsigned int ptIdxA = curIdxT->Get(i % );
         unsigned int ptIdxB = curIdxT->Get((i + ) % );
         unsigned int ptIdxC = curIdxT->Get((i + ) % );

         //Looking For Triangles Containing Line AB.
         std::array<size_t, > trianFound;
         bool isFound = GetTrianIndexWithLine(ptIdxA, ptIdxB, trianFound, ptIdxC);
         if (!isFound)
         {
             //TODO: Not Found The Triangles.
             //
             continue;
         }

         for (int eachIdx = ; eachIdx < trianFound.size()/*eachIdx < 2*/; ++eachIdx)
         {
             int foundIdx = static_cast<int>(trianFound[eachIdx]);
             if (foundIdx < /* == -1*/)
             {
                 continue;
             }

             const TrianIndex& trianABD = m_triangles[foundIdx];

             //Find The Point In eachT, And Not Containing A & B.
             unsigned int ptIdxD = ;
             for (int epi = ; epi < ; ++epi)
             {
                 unsigned int ep = trianABD[epi];
                 if (ep != ptIdxA && ep != ptIdxB)
                 {
                     ptIdxD = ep;
                 }
             }

             const Triangle curT(m_points[trianABD[]], m_points[trianABD[]], m_points[trianABD[]]);
             Circle<double, double> cir;
             curT.Circumcircle(cir);
             if (MyMathTools::LessThan(cir.radius/*Length OD*/, Point2D<double>::Distance(cir.center, m_points[ptIdxC])/*Length OC*/))
             {
                 continue;
             }

             //ABC(curIdxT) & ABD(foundIdx) -> ACD & BCD
             //A -> ptIdxA; B -> ptIdxB; C -> ptIdxC; D -> ptIdxD
             curIdxT->Init(ptIdxA, ptIdxC, ptIdxD);
             m_triangles[foundIdx].Init(ptIdxB, ptIdxC, ptIdxD);

             SetProcessedFunc(trianIndex, false);
             SetProcessedFunc(foundIdx, false);
             LOP(trianIndex, isProcessed);
             LOP(foundIdx, isProcessed);
         }
     }
     SetProcessedFunc(trianIndex, true);
 }

 void Delaunay::Clear()
 {
     m_points.Clear();
     m_triangles.clear();
 }

File: TestMain.cpp (win32 默认程序)

 case WM_PAINT:
         hdc = BeginPaint(hWnd, &ps);
         {
             HBRUSH hBrush = CreateSolidBrush(RGB(, , ));
             HBRUSH hOldBrush = (HBRUSH)SelectObject(hdc, hBrush);

             const int count = ;
             Delaunay dela;

             dela.AddRandomPoints(count, , );
             dela.Generate();

             Points points;
             Triangles ts;
             dela.GetAllTriangles(ts);
             dela.GetAllPoints(points);

             const int pointWidth = ;
             int i = ;
             for (int i = ; i < points.Size(); ++i)
             {
                 const Point2D<int>& pt = points[i];
                 Rectangle(hdc, pt.x - pointWidth, pt.y - pointWidth, pt.x + pointWidth, pt.y + pointWidth);
                 i++;
             }

             for (const TrianIndex& et : ts)
             {
                 Triangle tmpTri(points[et[]], points[et[]], points[et[]]);

                 std::array<Point2D<int>, > pts;
                 tmpTri.GetVertices(pts);

                 std::array<Line<int>, > ls;
                 tmpTri.GetLines(ls);

                 for (int idx = ; idx < ; idx++)
                 {
                     const Point2D<int>& pt = pts[idx];
                     Rectangle(hdc, pt.x - pointWidth, pt.y - pointWidth, pt.x + pointWidth, pt.y + pointWidth);

                     const Line<int>& l = ls[idx];
                     const Point2D<int>& ptA = l.GetPointA();
                     const Point2D<int>& ptB = l.GetPointB();

                     MoveToEx(hdc, ptA.x, ptA.y, NULL);
                     LineTo(hdc, ptB.x, ptB.y);
                 }
             }

             SelectObject(hdc, hOldBrush);
             DeleteObject(hBrush);
         }

         EndPaint(hWnd, &ps);
         break;

生成结果：