c++ - 为什么我的 Bowyer-Watson 有时会得到不与超级三角形共享顶点的三角形？
问题描述
我尝试了几种不同的方法来计算一个点是否在外接圆内，但我不断得到，偶尔，没有三角形与原始超三角形不共享一个点。
我的代码：
#include <iostream>
#include <vector>
#include <random>
#include <algorithm>
#include <chrono>
bool is_point_inside_circumcircle(const int32_t pX, const int32_t pY,
int32_t aX, int32_t aY, int32_t bX, int32_t bY, int32_t cX, int32_t cY)
{
const uint32_t ab = aX * aX + aY * aY;
const uint32_t bc = bX * bX + bY * bY;
const uint32_t ca = cX * cX + cY * cY;
uint32_t circumX = (ab * (cY - bY) + bc * (aY - cY) + ca * (bY - aY)) / (aX * (cY - bY) + bX * (aY - cY) + cX * (bY - aY));
uint32_t circumY = (ab * (cX - bX) + bc * (aX - cX) + ca * (bX - aX)) / (aY * (cX - bX) + bY * (aX - cX) + cY * (bX - aX));
circumX /= 2;
circumY /= 2;
uint32_t dX = aX - circumX;
uint32_t dY = aY - circumY;
uint32_t circumR = dX * dX + dY * dY;
dX = pX - circumX;
dY = pY - circumY;
uint32_t d = dX * dX + dY * dY;
return d <= circumR;
}
struct Edge
{
int32_t x1;
int32_t y1;
int32_t x2;
int32_t y2;
};
bool operator==(const Edge& lhs, const Edge& rhs)
{
return (lhs.x1 == rhs.x1 && lhs.y1 == rhs.y1 && lhs.x2 == rhs.x2 && lhs.y2 == rhs.y2)
|| (lhs.x1 == rhs.x2 && lhs.y1 == rhs.y2 && lhs.x2 == rhs.x1 && lhs.y2 == rhs.y1);
}
struct Triangle
{
Edge ab;
Edge bc;
Edge ca;
};
bool triangle_contains(const Triangle& tri, int32_t x, int32_t y)
{
return (tri.ab.x1 == x && tri.ab.y1 == y) ||
(tri.bc.x1 == x && tri.bc.y1 == y) ||
(tri.ca.x1 == x && tri.ca.y1 == y);
}
bool operator==(const Triangle& lhs, const Triangle& rhs)
{
return
((lhs.ab.x1 == rhs.ab.x1 && lhs.ab.y1 == rhs.ab.y1) ||
(lhs.ab.x1 == rhs.bc.x1 && lhs.ab.y1 == rhs.bc.y1) ||
(lhs.ab.x1 == rhs.ca.x1 && lhs.ca.y1 == rhs.ca.y1)) &&
((lhs.bc.x1 == rhs.ab.x1 && lhs.bc.y1 == rhs.ab.y1) ||
(lhs.bc.x1 == rhs.bc.x1 && lhs.bc.y1 == rhs.bc.y1) ||
(lhs.bc.x1 == rhs.ca.x1 && lhs.bc.y1 == rhs.ca.y1)) &&
((lhs.ca.x1 == rhs.ab.x1 && lhs.ca.y1 == rhs.ab.y1) ||
(lhs.ca.x1 == rhs.bc.x1 && lhs.ca.y1 == rhs.bc.y1) ||
(lhs.ca.x1 == rhs.ca.x1 && lhs.ca.y1 == rhs.ca.y1));
}
bool operator!=(const Triangle& lhs, const Triangle& rhs)
{
return !(lhs==rhs);
}
std::vector<Triangle> generate_galaxy(const std::size_t numSystems)
{
std::vector<std::pair<int32_t,int32_t>> galaxy;
std::mt19937 gen(std::chrono::system_clock::now().time_since_epoch().count());
std::uniform_int_distribution<> dimension(0, numSystems*3-1);
for(std::size_t i = 0; i < numSystems; ++i)
{
galaxy.emplace_back(dimension(gen),dimension(gen));
}
std::vector<Triangle> triangulation;
Triangle super;
int32_t minX = std::numeric_limits<int32_t>::max();
int32_t minY = std::numeric_limits<int32_t>::max();
int32_t maxX = 0;
int32_t maxY = 0;
for(const std::pair<int32_t,int32_t>& n : galaxy)
{
if(n.first < minX) minX = n.first;
if(n.second < minY) minY = n.second;
if(n.first > maxX) maxX = n.first;
if(n.second > maxY) maxY = n.second;
}
int32_t dX = maxX - minX;
int32_t dY = maxY - minY;
int32_t dM = std::max(dX,dY);
int32_t midX = (minX + maxX) / 2;
int32_t midY = (minY + maxY) / 2;
super.ab.x1 = super.ca.x2 = midX - 20 * dM;
super.ab.y1 = super.ca.y2 = midY - dM;
super.ab.x2 = super.bc.x1 = midX;
super.ab.y2 = super.bc.y1 = midY + 20 * dM;
super.bc.x2 = super.ca.x1 = midX + 20 * dM;
super.bc.y2 = super.ca.y1 = midY - dM;
triangulation.push_back(super);
for(const std::pair<int32_t,int32_t>& n : galaxy)
{
std::vector<Triangle> badTriangles;
for(const Triangle& triangle : triangulation)
{
if(is_point_inside_circumcircle(n.first, n.second, triangle.ab.x1, triangle.ab.x2, triangle.bc.x1, triangle.bc.y1,
triangle.ca.x1, triangle.ca.y1))
{
badTriangles.push_back(triangle);
}
}
std::cout << "\n--------------------------------------------" << std::endl;
std::cout << "\n\nAFTER NEW NODE CREATING BAD TRIANGLES" << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nTriangulation (" << triangulation.size() << ")";
for(const Triangle& tri : triangulation)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nBad (" << badTriangles.size() << ")";
for(const Triangle& tri : badTriangles)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::vector<Edge> polygon;
std::vector<Triangle>::iterator triangleItr = badTriangles.begin();
while(triangleItr != badTriangles.end())
{
if(std::find_if(badTriangles.begin(),badTriangles.end(),[&](const Triangle& tri)
{
return
*triangleItr != tri && (
triangleItr->ab == tri.ab ||
triangleItr->ab == tri.bc ||
triangleItr->ab == tri.ca);
}) == badTriangles.end())
{
polygon.push_back(triangleItr->ab);
}
if(std::find_if(badTriangles.begin(),badTriangles.end(),[&](const Triangle& tri)
{
return
*triangleItr != tri && (
triangleItr->bc == tri.ab ||
triangleItr->bc == tri.bc ||
triangleItr->bc == tri.ca);
}) == badTriangles.end())
{
polygon.push_back(triangleItr->bc);
}
if(std::find_if(badTriangles.begin(),badTriangles.end(),[&](const Triangle& tri)
{
return
*triangleItr != tri && (
triangleItr->ca == tri.ab ||
triangleItr->ca == tri.bc ||
triangleItr->ca == tri.ca);
}) == badTriangles.end())
{
polygon.push_back(triangleItr->ca);
}
++triangleItr;
}
triangulation.erase(std::remove_if(triangulation.begin(),triangulation.end(),[&](const Triangle& tri)
{
for(const Triangle& triangle : badTriangles)
{
if(tri == triangle)
return true;
}
return false;
}),triangulation.end());
std::cout << "\n\nAFTER ADDING EDGES TO POLYGON BECAUSE NOT SHARED IN BAD TRIANGLES AND REMOVE BAD TRIANGLES FROM TRIANGULATION" << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nTriangulation (" << triangulation.size() << ")";
for(const Triangle& tri : triangulation)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nBad (" << badTriangles.size() << ")";
for(const Triangle& tri : badTriangles)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nPoly (" << polygon.size() << ")";
for(const Edge& e : polygon)
std::cout << "\n\t" << e.x1 << "," << e.y1 << " > " << e.x2 << "," << e.y2;
std::cout << std::endl;
for(const Edge& edge : polygon)
{
Triangle newTri;
newTri.ab = edge;
Edge bc;
bc.x1 = edge.x2;
bc.y1 = edge.y2;
bc.x2 = n.first;
bc.y2 = n.second;
Edge ca;
ca.x1 = n.first;
ca.y1 = n.second;
ca.x2 = edge.x1;
ca.y2 = edge.y1;
newTri.bc = bc;
newTri.ca = ca;
triangulation.push_back(newTri);
}
std::cout << "\n\nADDING NEW TRIANGLES FROM EDGES TO TRIANGULATION" << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nTriangulation (" << triangulation.size() << ")";
for(const Triangle& tri : triangulation)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nBad (" << badTriangles.size() << ")";
for(const Triangle& tri : badTriangles)
std::cout << "\n\t" << tri.ab.x1 << "," << tri.ab.y1 << " "
<< tri.bc.x1 << "," << tri.bc.y1 << " "
<< tri.ca.x1 << "," << tri.ca.y1;
std::cout << std::endl;
std::cout << "NODE X,Y: " << n.first << "," << n.second
<< "\nPoly (" << polygon.size() << ")";
for(const Edge& e : polygon)
std::cout << "\n\t" << e.x1 << "," << e.y1 << " > " << e.x2 << "," << e.y2;
std::cout << std::endl;
}
triangulation.erase(std::remove_if(triangulation.begin(),triangulation.end(),[&](const Triangle& tri)
{
return triangle_contains(tri,super.ab.x1,super.ab.y1) ||
triangle_contains(tri,super.bc.x1,super.bc.y1) ||
triangle_contains(tri,super.ca.x1,super.ca.y1);
}),triangulation.end());
// every edge in triangulation is now the connection of vertices
return triangulation;
}
int main()
{
std::vector<Triangle> galaxy = generate_galaxy(4);
std::cout << "\n\nOUTCOMES:\n" << std::endl;
for(const Triangle& tri : galaxy)
{
std::cout << tri.ab.x1 << "," << tri.ab.y1 << " > " << tri.ab.x2 << "," << tri.ab.y2 << std::endl;
std::cout << tri.bc.x1 << "," << tri.bc.y1 << " > " << tri.bc.x2 << "," << tri.bc.y2 << std::endl;
std::cout << tri.ca.x1 << "," << tri.ca.y1 << " > " << tri.ca.x2 << "," << tri.ca.y2 << std::endl;
}
return 0;
}
我也尝试过从https://github.com/Bl4ckb0ne/delaunay-triangulation/blob/master/dt/delaunay.cpp采取一些方法，希望它会有所帮助，但似乎没有有很大的不同。
当没有生成三角形时，上面的输出如下所示：
--------------------------------------------
AFTER NEW NODE CREATING BAD TRIANGLES
NODE X,Y: 2,10
Triangulation (1)
-154,-1 6,167 166,-1
NODE X,Y: 2,10
Bad (1)
-154,-1 6,167 166,-1
AFTER ADDING EDGES TO POLYGON BECAUSE NOT SHARED IN BAD TRIANGLES AND REMOVE BAD TRIANGLES FROM TRIANGULATION
NODE X,Y: 2,10
Triangulation (0)
NODE X,Y: 2,10
Bad (1)
-154,-1 6,167 166,-1
NODE X,Y: 2,10
Poly (3)
-154,-1 > 6,167
6,167 > 166,-1
166,-1 > -154,-1
ADDING NEW TRIANGLES FROM EDGES TO TRIANGULATION
NODE X,Y: 2,10
Triangulation (3)
-154,-1 6,167 2,10
6,167 166,-1 2,10
166,-1 -154,-1 2,10
NODE X,Y: 2,10
Bad (1)
-154,-1 6,167 166,-1
NODE X,Y: 2,10
Poly (3)
-154,-1 > 6,167
6,167 > 166,-1
166,-1 > -154,-1
--------------------------------------------
AFTER NEW NODE CREATING BAD TRIANGLES
NODE X,Y: 3,6
Triangulation (3)
-154,-1 6,167 2,10
6,167 166,-1 2,10
166,-1 -154,-1 2,10
NODE X,Y: 3,6
Bad (2)
-154,-1 6,167 2,10
6,167 166,-1 2,10
AFTER ADDING EDGES TO POLYGON BECAUSE NOT SHARED IN BAD TRIANGLES AND REMOVE BAD TRIANGLES FROM TRIANGULATION
NODE X,Y: 3,6
Triangulation (1)
166,-1 -154,-1 2,10
NODE X,Y: 3,6
Bad (2)
-154,-1 6,167 2,10
6,167 166,-1 2,10
NODE X,Y: 3,6
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
ADDING NEW TRIANGLES FROM EDGES TO TRIANGULATION
NODE X,Y: 3,6
Triangulation (5)
166,-1 -154,-1 2,10
-154,-1 6,167 3,6
2,10 -154,-1 3,6
6,167 166,-1 3,6
166,-1 2,10 3,6
NODE X,Y: 3,6
Bad (2)
-154,-1 6,167 2,10
6,167 166,-1 2,10
NODE X,Y: 3,6
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
--------------------------------------------
AFTER NEW NODE CREATING BAD TRIANGLES
NODE X,Y: 4,5
Triangulation (5)
166,-1 -154,-1 2,10
-154,-1 6,167 3,6
2,10 -154,-1 3,6
6,167 166,-1 3,6
166,-1 2,10 3,6
NODE X,Y: 4,5
Bad (4)
-154,-1 6,167 3,6
2,10 -154,-1 3,6
6,167 166,-1 3,6
166,-1 2,10 3,6
AFTER ADDING EDGES TO POLYGON BECAUSE NOT SHARED IN BAD TRIANGLES AND REMOVE BAD TRIANGLES FROM TRIANGULATION
NODE X,Y: 4,5
Triangulation (1)
166,-1 -154,-1 2,10
NODE X,Y: 4,5
Bad (4)
-154,-1 6,167 3,6
2,10 -154,-1 3,6
6,167 166,-1 3,6
166,-1 2,10 3,6
NODE X,Y: 4,5
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
ADDING NEW TRIANGLES FROM EDGES TO TRIANGULATION
NODE X,Y: 4,5
Triangulation (5)
166,-1 -154,-1 2,10
-154,-1 6,167 4,5
2,10 -154,-1 4,5
6,167 166,-1 4,5
166,-1 2,10 4,5
NODE X,Y: 4,5
Bad (4)
-154,-1 6,167 3,6
2,10 -154,-1 3,6
6,167 166,-1 3,6
166,-1 2,10 3,6
NODE X,Y: 4,5
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
--------------------------------------------
AFTER NEW NODE CREATING BAD TRIANGLES
NODE X,Y: 10,8
Triangulation (5)
166,-1 -154,-1 2,10
-154,-1 6,167 4,5
2,10 -154,-1 4,5
6,167 166,-1 4,5
166,-1 2,10 4,5
NODE X,Y: 10,8
Bad (4)
-154,-1 6,167 4,5
2,10 -154,-1 4,5
6,167 166,-1 4,5
166,-1 2,10 4,5
AFTER ADDING EDGES TO POLYGON BECAUSE NOT SHARED IN BAD TRIANGLES AND REMOVE BAD TRIANGLES FROM TRIANGULATION
NODE X,Y: 10,8
Triangulation (1)
166,-1 -154,-1 2,10
NODE X,Y: 10,8
Bad (4)
-154,-1 6,167 4,5
2,10 -154,-1 4,5
6,167 166,-1 4,5
166,-1 2,10 4,5
NODE X,Y: 10,8
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
ADDING NEW TRIANGLES FROM EDGES TO TRIANGULATION
NODE X,Y: 10,8
Triangulation (5)
166,-1 -154,-1 2,10
-154,-1 6,167 10,8
2,10 -154,-1 10,8
6,167 166,-1 10,8
166,-1 2,10 10,8
NODE X,Y: 10,8
Bad (4)
-154,-1 6,167 4,5
2,10 -154,-1 4,5
6,167 166,-1 4,5
166,-1 2,10 4,5
NODE X,Y: 10,8
Poly (4)
-154,-1 > 6,167
2,10 > -154,-1
6,167 > 166,-1
166,-1 > 2,10
OUTCOMES:
Process returned 0 (0x0) execution time : 0.620 s
Press any key to continue.
我可以看到三角剖分只包含与计算的原始超三角形共享一个顶点的三角形。所以我知道为什么它们都被删除了，但我在网上找不到任何关于为什么会发生这种情况的信息。我逐步完成了该过程，并将其与Bowyer–Watson 算法的描述和图像相匹配，我看不到任何问题。
解决方案
推荐阅读
- python - 使用 Python 读取 Exif：MakerNote IFD 中可能损坏的字段标记 0x0000
- twitter-bootstrap - Twitter 按钮仅显示推文文本，无按钮
- python - 使用 Python 脚本运行完全独立的 powershell 脚本
- python - Python二叉搜索树删除函数
- android - Android - 确定显示尺寸
- swift - 向下转换“任何”时出现 Swift 错误
- hive - 代理键是如何在 hive 中处理的
- c# - RSS阅读器部分读取节点
- dart - Dart/Flutter：iOS 应用因 Firebase ML Vision 崩溃
- c - 如何在 C 中使用可变数量的变量？