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

  1 #pragma once
  2 
  3 #include "ConvexHull.h"
  4 #include "Triangle.h"
  5 #include "Circle.cpp"
  6 
  7 class TrianIndex  final
  8 {
  9 public:
 10     TrianIndex()
 11     {
 12         _isVaild = false;
 13     }
 14 
 15     ~TrianIndex()
 16     {
 17         _isVaild = false;
 18     }
 19 
 20     TrianIndex(unsigned int iA, unsigned int iB, unsigned int iC)
 21     {
 22         Init(iA,  iB,  iC);
 23     }
 24 
 25     TrianIndex(std::array<unsigned int, 3> pts) :
 26         TrianIndex(pts[0], pts[1], pts[2])
 27     {
 28     }
 29 
 30     TrianIndex(const TrianIndex& other)
 31     {
 32         this->_isVaild = other._isVaild;
 33         this->_ptIdx = other._ptIdx;
 34     }
 35 
 36     unsigned int& Get(int i)
 37     {
 38         if (i < 0 || i > 2 || !_isVaild)
 39         {
 40             ErrorThrow("Error Triangle Point Index[0, 2] Or Invaild Triangle: " + std::to_string(i));
 41             return _ptIdx[0];
 42         }
 43 
 44         return _ptIdx[i];
 45     }
 46 
 47     unsigned int& operator[](int i)
 48     {
 49         return Get(i);
 50     }
 51 
 52     const unsigned int Get(int i) const
 53     {
 54         TrianIndex* pThis = const_cast<TrianIndex*>(this);
 55         return pThis->Get(i);
 56     }
 57     
 58     const unsigned int operator[](int i) const
 59     {
 60         TrianIndex* pThis = const_cast<TrianIndex*>(this);
 61         return (*pThis)[i];
 62     }
 63 
 64     TrianIndex& operator=(const TrianIndex& other)
 65     {
 66         this->_isVaild = other._isVaild;
 67         this->_ptIdx = other._ptIdx;
 68         return *this;
 69     }
 70 
 71     bool IsVaild() const
 72     {
 73         return _isVaild;
 74     }
 75 
 76     void SetVaild(bool isVaild)
 77     {
 78         _isVaild = isVaild;
 79     }
 80 
 81     void Init(unsigned int iA, unsigned int iB, unsigned int iC)
 82     {
 83         _ptIdx[0] = iA;
 84         _ptIdx[1] = iB;
 85         _ptIdx[2] = iC;
 86         _isVaild = true;
 87     }
 88 
 89 private:
 90     bool _isVaild;
 91     std::array<unsigned int, 3> _ptIdx;
 92 };
 93 
 94 //CHECK IT: It Has Been Abandoned.
 95 //
 96 class TrianMatrix final
 97 {
 98 public:
 99     enum IsConnected : bool
100     {
101         NotConnected = false,
102         Connected = true
103     };
104 
105     TrianMatrix():
106         TrianMatrix(0)
107     {
108     }
109 
110     TrianMatrix(size_t ptCount)
111     {
112         _sign.clear();
113         _sign.resize(ptCount);
114         for (unsigned int x = 0; x < ptCount; ++x)
115         {
116             _sign[x].resize(ptCount);
117             for (unsigned int y = 0; y < ptCount; ++y)
118             {
119                 _sign[x][y] = IsConnected::NotConnected;
120             }
121         }
122     }
123 
124     ~TrianMatrix()
125     {
126         _sign.clear();
127     }
128 
129 public:
130     size_t PointCount() const
131     {
132         return _sign.size();
133     }
134 
135     bool IsEmpty() const
136     {
137         return _sign.empty();
138     }
139 
140     //
141     bool GetTrianglesByPoint(unsigned int ptIdx, 
142         std::vector<TrianIndex>& __out resTrians) const
143     {
144         bool isFound = false;
145         if (!CheckIndex(ptIdx))
146         {
147             return false;
148         }
149 
150         //Time Complexity: O(n^2).
151         for (unsigned int i = 0; i < PointCount(); ++i)
152         {
153             isFound = GetTrianglesByLine(ptIdx, i, resTrians);
154         }
155 
156         return isFound;
157     }
158 
159     bool GetTrianglesByLine(unsigned int ptIdxA, unsigned int ptIdxB, 
160         std::vector<TrianIndex>& __out resTrians) const
161     {
162         bool isFound = false;
163         short foundCount = 0;
164         if (!CheckIndex(ptIdxA, ptIdxB))
165         {
166             return false;
167         }
168 
169         //Time Complexity: O(n).
170         for (unsigned int i = 0; i < PointCount(); ++i)
171         {
172             const std::vector<IsConnected>& pts = _sign[i];
173             if (pts[ptIdxA] && pts[ptIdxB])
174             {
175                 TrianIndex trian(i, ptIdxA, ptIdxB);        //The another point at the first element of resTrians.
176                 resTrians.push_back(trian);
177                 isFound = true;
178 
179                 if (++foundCount == 2)
180                 {
181                     return true;
182                 }
183             }
184         }
185         return isFound;
186     }
187 
188     inline bool Connect(unsigned int ptIdxA, unsigned int ptIdxB)
189     {
190         return SetConnectedState(ptIdxA, ptIdxB, IsConnected::Connected);
191     }
192 
193     inline bool DisConnect(unsigned int ptIdxA, unsigned int ptIdxB)
194     {
195         return SetConnectedState(ptIdxA, ptIdxB, IsConnected::NotConnected);
196     }
197 
198     inline bool SetConnectedState(unsigned int ptIdxA, unsigned int ptIdxB, IsConnected isConn)
199     {
200         if (!CheckIndex(ptIdxA, ptIdxB))
201         {
202             return false;
203         }
204 
205         _sign[ptIdxA][ptIdxB] = isConn;
206         _sign[ptIdxB][ptIdxA] = isConn;
207         return true;
208     }
209 
210     bool IsExist(const TrianIndex& trian) const
211     {
212         if (!CheckIndex(trian[0], trian[1], trian[2]))
213         {
214             return false;
215         }
216 
217         unsigned int ptA = trian[0];
218         unsigned int ptB = trian[1];
219         unsigned int ptC = trian[2];
220 
221         return (_sign[ptA][ptB] && _sign[ptA][ptC] && _sign[ptB][ptC]);
222     }
223 
224 private:
225     bool CheckIndex(unsigned int idx, ...) const
226     {
227         std::vector<unsigned int> indeies;
228         auto CheckFunc = [&indeies, this](unsigned int newIdx) -> bool
229         {
230             if (newIdx >= this->PointCount())
231             {
232                 return false;
233             }
234 
235             for (std::vector<unsigned int>::const_iterator it = indeies.begin();
236                 it != indeies.end(); ++it)
237             {
238                 if (*it == newIdx)
239                 {
240                     return true;
241                 }
242             }
243             indeies.push_back(newIdx);
244             return false;
245         };
246         
247         va_list aps;
248         va_start(aps, idx);
249         unsigned int* pNextArg = nullptr;
250         bool isVaild = false;
251         do 
252         {
253             pNextArg = va_arg(aps, unsigned int*);
254             if (pNextArg) 
255             {
256                 unsigned int tmp = *pNextArg;
257                 if (!CheckFunc(tmp))
258                 {
259                     isVaild = true;
260                 }
261             }
262         } while (pNextArg != NULL);
263 
264         va_end(aps);
265         return isVaild;
266     }
267     
268 
269 private:
270         std::vector<
271             std::vector<IsConnected>
272     > _sign;
273 };
274 
275 typedef std::vector<TrianIndex> Triangles;
276 class Delaunay final
277 {
278 public:
279     Delaunay();
280     ~Delaunay();
281 
282 public:
283     void AddPoint(const Point2D<int>& __in newPoint);
284     void AddRandomPoints(int count, int maxX, int maxY);
285     
286     void GetAllTriangles(Triangles& __out ts) const;
287     void GetAllPoints(Points& __out pts) const;
288     void Generate();
289     
290     void Clear();
291 
292 private:
293     bool GetTrianIndexWithLine(int ptIdxA, int ptIdxB, std::array<size_t/*The Triangle Index Of m_triangles.*/, 2> & __out res, int ptIdxCur = -1) const;
294     void CutProcess(const Points& chPoints);
295     void ConvexHullCut();
296     void LOP(unsigned int trianIndex, bool* isProcessed);    //Local Optimization Procedure.
297 
298 private:
299     Points m_points;
300     Triangles m_triangles;
301 };

File: Delaunay.cpp

  1 #include "Delaunay.h"
  2 #include <bitset>
  3 
  4 Delaunay::Delaunay()
  5 {
  6     m_points.Clear();
  7     m_triangles.clear();
  8 }
  9 
 10 
 11 Delaunay::~Delaunay()
 12 {
 13     Clear();
 14 }
 15 
 16 void Delaunay::AddPoint(const Point2D<int>& newPoint)
 17 {
 18     m_points.Add(newPoint);
 19 }
 20 
 21 void Delaunay::AddRandomPoints(int count, int maxX, int maxY)
 22 {
 23     std::vector<Point2D<int>> pointsArray;
 24     Point2D<int>::RandomPoints(count, maxX, maxY, pointsArray);
 25 
 26     m_points.Clear();
 27     for (int i = 0; i < count; ++i)
 28     {
 29         const Point2D<int>& eachpt = pointsArray.at(i);
 30         PointNode pn;    pn.Init(i, eachpt);
 31         m_points.Add(pn);
 32     }
 33 }
 34 
 35 void Delaunay::ConvexHullCut()
 36 {
 37     if (m_points.Size() < 3)
 38     {
 39         ErrorThrow("Points count too less.");
 40         return;
 41     }
 42     m_triangles.clear();
 43 
 44     ConvexHull ch(m_points);
 45     Points chPoints;
 46     ch.Generate();
 47     ch.GetConvexHullPoints(chPoints);
 48 
 49     CutProcess(chPoints);
 50 }
 51 
 52 bool Delaunay::GetTrianIndexWithLine(int ptIdxA, int ptIdxB, std::array<size_t, 2> & __out res, int ptIdxCur) const
 53 {
 54     if (ptIdxA == ptIdxB || ptIdxCur == ptIdxA || ptIdxCur == ptIdxB)
 55     {
 56         ErrorThrow("Error Parameters! \n ptIdxA: " + std::to_string(ptIdxA) + ", ptIdxB: " + std::to_string(ptIdxB) + ", ptIdxCur: " + std::to_string(ptIdxCur));
 57         return false;
 58     }
 59     res.fill(-1);
 60     uint8_t cnt = 0;
 61     bool isFound = false;
 62 
 63     for(size_t i =0; i < m_triangles.size(); ++i)
 64     {
 65         if (cnt > 2)
 66         {
 67             return isFound;
 68         }
 69 
 70         std::bitset<3> foundSign = 0;
 71         int8_t bit0Idx = -1;
 72         const TrianIndex& tIdx = m_triangles[i];
 73         for (uint8_t i = 0; i < 3; i++)
 74         {
 75             if (tIdx[i] == ptIdxA || tIdx[i] == ptIdxB)
 76             {
 77                 foundSign = (1 << i) | foundSign.to_ulong();
 78             }
 79             else
 80             {
 81                 bit0Idx = i;
 82             }
 83         }
 84 
 85         if (2 == foundSign.count() && tIdx[bit0Idx] != ptIdxCur)
 86         {
 87             res[cnt++] = i;
 88             isFound = true;
 89         }
 90     }
 91     return isFound;
 92 }
 93 
 94 void Delaunay::GetAllTriangles(Triangles & ts) const
 95 {
 96     ts = m_triangles;
 97 }
 98 
 99 void Delaunay::GetAllPoints(Points & pts) const
100 {
101     pts = m_points;
102 }
103 
104 void Delaunay::CutProcess(const Points & chPoints)
105 {
106     int chSize = static_cast<int>(chPoints.Size());
107 
108     //Next Convex Hull Points.
109     Points nextPts;
110 
111     std::array<PointNode, 3> aTriangle;
112     int i = 0;
113     int atSize = 0;    /*atSize: aTriangle size.*/
114 
115     while (i <= chSize) /*Including the first point.*/
116     {
117         const PointNode*const cutpt = chPoints.GetBySequence(i % chSize);
118         aTriangle[atSize++] = *cutpt;
119 
120         if (atSize == 3)
121         {
122             if (Line<int>::IsCollinear(3, aTriangle[0]._node, aTriangle[1]._node, aTriangle[2]._node))
123             {
124                 //Delete the second point.
125                 aTriangle[1] = aTriangle[2];
126                 --atSize;
127                 ++i;
128                 continue;
129             }
130 
131             //Form a triangle.
132             TrianIndex ttmp(aTriangle[0]._index, aTriangle[1]._index, aTriangle[2]._index);
133             m_triangles.push_back(ttmp);
134 
135             //The last serves as the starting point for the next triangle.
136             aTriangle[0] = aTriangle[2];
137             atSize = 1;
138             
139             nextPts.Add(aTriangle[2]);
140         }
141         ++i;
142     }
143 
144     //Add the remaining points.
145     //The first point in aTriangle has been added.
146     for (int i = 1; i < atSize; ++i)
147     {
148         nextPts.Add(aTriangle[i]);
149     }
150     
151 
152     if (nextPts.Size() >= 3)
153     {
154         CutProcess(nextPts);
155     }
156 }
157 
158 //Generate the delaunay.
159 void Delaunay::Generate()
160 {
161     ConvexHullCut();
162 
163     for (int i = 0 ; i < m_points.Size(); ++i)
164     {
165         const Point2D<int>& eachPt = m_points(i);
166         for(std::vector<TrianIndex>::iterator it = m_triangles.begin(); it != m_triangles.end(); )
167         {
168             const TrianIndex eachTri = *it;
169             const Point2D<int>& ptA = m_points[eachTri[0]];
170             const Point2D<int>& ptB = m_points[eachTri[1]];
171             const Point2D<int>& ptC = m_points[eachTri[2]];
172             Triangle tmpTri(ptA, ptB, ptC);
173 
174             bool isBreak = false;
175             Triangle::EnumPTRelation eptr = tmpTri.RelationPT(eachPt);
176             switch (eptr)
177             {
178             case Triangle::EnumPTRelation::POINT_ON_VERTICE:
179             case Triangle::EnumPTRelation::POINT_OUTSIDE:
180                 ++it;
181                 isBreak = false;
182                 break;
183 
184             case Triangle::EnumPTRelation::POINT_ON_LINE:
185                 for (unsigned int j = 0; j < 3; ++j)
186                 {
187                     unsigned int ptIdxA = j % 3;
188                     unsigned int ptIdxB = (j + 1) % 3;
189                     unsigned int ptIdxC = (j + 2) % 3;
190                     if (Line<int>::IsCollinear(3, m_points[eachTri[ptIdxA]], m_points[eachTri[ptIdxB]], eachPt))
191                     {
192                         TrianIndex tmpAPC(eachTri[ptIdxC], i, eachTri[ptIdxA]);
193                         TrianIndex tmpAPB(eachTri[ptIdxC], i, eachTri[ptIdxB]);
194                         it = m_triangles.erase(it);    //Delete ABC.
195                         m_triangles.push_back(tmpAPC);
196                         m_triangles.push_back(tmpAPB);
197                         break;
198                     }
199                 }
200                 isBreak = false;
201                 break;
202 
203             case Triangle::EnumPTRelation::POINT_INSIDE:
204                 it = m_triangles.erase(it);
205                 for (unsigned int j = 0; j < 3; ++j)
206                 {
207                     unsigned int ptIdx = j % 3;
208                     unsigned int ptIdxAnother = (j + 1) % 3;
209 
210                     TrianIndex tmp(eachTri[ptIdx], i, eachTri[ptIdxAnother]);
211                     m_triangles.push_back(tmp);
212                 }
213                 isBreak = true;
214                 break;
215             }
216 
217             if(isBreak)
218             {
219                 break;
220             }
221         }
222     }
223 
224     //LOP
225     const size_t TrianglesSize = m_triangles.size();
226     bool* isProcessed = nullptr;
227     if (isProcessed == nullptr)
228     {
229         isProcessed = new bool[TrianglesSize];
230         memset(isProcessed, false, TrianglesSize);
231     }
232     
233     for (unsigned int i = 0; i < TrianglesSize; ++i)
234     {
235         LOP(i, isProcessed);
236     }
237     
238     if (isProcessed != nullptr)
239     {
240         delete[] isProcessed;
241     }
242 }
243 
244 void Delaunay::LOP(unsigned int trianIndex, bool* isProcessed = nullptr)
245 {
246     auto SetProcessedFunc = [&isProcessed](unsigned int idx, bool val)
247     {
248         if (nullptr != isProcessed)
249         {
250             isProcessed[idx] = val;
251         }
252     };
253 
254     if (trianIndex >= m_triangles.size())
255     {
256         ErrorThrow("The Index Of Triangles Out Of Range. trianIndex: " + std::to_string(trianIndex));
257         return;
258     }
259 
260     if (isProcessed != nullptr && isProcessed[trianIndex])
261     {
262         return;
263     }
264 
265     TrianIndex* const curIdxT = &m_triangles[trianIndex];
266     for (unsigned int i = 0; i < 3; ++i)
267     {
268         unsigned int ptIdxA = curIdxT->Get(i % 3);
269         unsigned int ptIdxB = curIdxT->Get((i + 1) % 3);
270         unsigned int ptIdxC = curIdxT->Get((i + 2) % 3);
271 
272         //Looking For Triangles Containing Line AB. 
273         std::array<size_t, 2> trianFound;
274         bool isFound = GetTrianIndexWithLine(ptIdxA, ptIdxB, trianFound, ptIdxC);
275         if (!isFound)
276         {
277             //TODO: Not Found The Triangles.
278             //
279             continue;
280         }
281 
282         for (int eachIdx = 0; eachIdx < trianFound.size()/*eachIdx < 2*/; ++eachIdx)
283         {
284             int foundIdx = static_cast<int>(trianFound[eachIdx]);
285             if (foundIdx < 0/* == -1*/)
286             {
287                 continue;
288             }
289 
290             const TrianIndex& trianABD = m_triangles[foundIdx];
291 
292             //Find The Point In eachT, And Not Containing A & B.
293             unsigned int ptIdxD = 0;
294             for (int epi = 0; epi < 3; ++epi)
295             {
296                 unsigned int ep = trianABD[epi];
297                 if (ep != ptIdxA && ep != ptIdxB)
298                 {
299                     ptIdxD = ep;
300                 }
301             }
302 
303             const Triangle curT(m_points[trianABD[0]], m_points[trianABD[1]], m_points[trianABD[2]]);
304             Circle<double, double> cir;
305             curT.Circumcircle(cir);
306             if (MyMathTools::LessThan(cir.radius/*Length OD*/, Point2D<double>::Distance(cir.center, m_points[ptIdxC])/*Length OC*/))
307             {
308                 continue;
309             }
310 
311             //ABC(curIdxT) & ABD(foundIdx) -> ACD & BCD
312             //A -> ptIdxA; B -> ptIdxB; C -> ptIdxC; D -> ptIdxD
313             curIdxT->Init(ptIdxA, ptIdxC, ptIdxD);
314             m_triangles[foundIdx].Init(ptIdxB, ptIdxC, ptIdxD);
315 
316             SetProcessedFunc(trianIndex, false);
317             SetProcessedFunc(foundIdx, false);
318             LOP(trianIndex, isProcessed);
319             LOP(foundIdx, isProcessed);
320         }
321     }
322     SetProcessedFunc(trianIndex, true);
323 }
324 
325 void Delaunay::Clear()
326 {
327     m_points.Clear();
328     m_triangles.clear();
329 }

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

 1 case WM_PAINT:
 2         hdc = BeginPaint(hWnd, &ps);
 3         {
 4             HBRUSH hBrush = CreateSolidBrush(RGB(255, 0, 45));
 5             HBRUSH hOldBrush = (HBRUSH)SelectObject(hdc, hBrush);
 6 
 7             const int count = 277;
 8             Delaunay dela;
 9             
10             dela.AddRandomPoints(count, 1200, 720);
11             dela.Generate();
12 
13             Points points;
14             Triangles ts;
15             dela.GetAllTriangles(ts);
16             dela.GetAllPoints(points);
17 
18             const int pointWidth = 3;
19             int i = 0;
20             for (int i = 0; i < points.Size(); ++i)
21             {
22                 const Point2D<int>& pt = points[i];
23                 Rectangle(hdc, pt.x - pointWidth, pt.y - pointWidth, pt.x + pointWidth, pt.y + pointWidth);
24                 i++;
25             }
26     
27             for (const TrianIndex& et : ts)
28             {
29                 Triangle tmpTri(points[et[0]], points[et[1]], points[et[2]]);
30 
31                 std::array<Point2D<int>, 3> pts;
32                 tmpTri.GetVertices(pts);
33 
34                 std::array<Line<int>, 3> ls;
35                 tmpTri.GetLines(ls);
36 
37                 for (int idx = 0; idx < 3; idx++)
38                 {
39                     const Point2D<int>& pt = pts[idx];
40                     Rectangle(hdc, pt.x - pointWidth, pt.y - pointWidth, pt.x + pointWidth, pt.y + pointWidth);
41 
42                     const Line<int>& l = ls[idx];
43                     const Point2D<int>& ptA = l.GetPointA();
44                     const Point2D<int>& ptB = l.GetPointB();
45 
46                     MoveToEx(hdc, ptA.x, ptA.y, NULL);
47                     LineTo(hdc, ptB.x, ptB.y);
48                 }
49             }
50 
51             SelectObject(hdc, hOldBrush);
52             DeleteObject(hBrush);
53         }
54 
55         EndPaint(hWnd, &ps);
56         break;

生成结果：

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

秒客网

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

相关文章