几何舍入问题:对象在简单变换后不再凸
Geometry rounding problems: object no longer convex after simple transformations
我正在制作一个分析几何的小应用程序。在我的程序的一部分中,我使用 has 的算法将凸对象作为输入。幸运的是,我所有的对象最初都是凸的,但有些只是 勉强 所以(见图)。
在我应用一些转换后,我的算法无法工作(它产生 "infinitely" 长多边形等),我认为这是因为图像中的舍入误差;由于舍入误差(在图像中非常夸张),圆柱体的顶部顶点略有 "pushed in",不再凸起。
所以我的问题是:有人知道 "slightly convexify" 对象的方法吗?这是我尝试实现的一种方法,但它似乎不起作用(或者我实现错误):
1. Average all vertices together to create a vertex C inside the convex shape.
2. Let d[v] be the distance from C to vertex v.
3. Scale each vertex v from the center C with the scale factor 1 / (1+d[v] * CONVEXIFICATION_FACTOR)
谢谢!!我安装了 CGAL 和 Boost,所以我可以使用这些库函数中的任何一个(我已经这样做了)。
您当然可以通过计算对象的 convex hull 使对象成为凸面。但这将 "convexify" 任何东西。如果您确定您的输入只是稍微偏离凸面,那么这应该不是问题。
CGAL 似乎在其中实现了 3D Quickhull,这将是第一个尝试的事情。有关文档和一些示例程序,请参阅 http://doc.cgal.org/latest/Convex_hull_3/。 (我对 CGAL 不够熟悉,无法重现任何示例并声称它们是正确的。)
最后我发现这个问题的根源是凸包包含很多三角形,而我的输入形状通常是立方体,使得每个四边形区域显示为 2 个具有极其相似平面的三角形方程式,导致我使用的算法出现某种问题。
我用 "detriangulating" 多面体解决了它,使用这段代码。如果有人能发现任何改进或问题,请告诉我!
#include <algorithm>
#include <cmath>
#include <vector>
#include <CGAL/convex_hull_traits_3.h>
#include <CGAL/convex_hull_3.h>
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Aff_transformation_3 Transformation;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
struct Plane_from_facet {
Polyhedron::Plane_3 operator()(Polyhedron::Facet& f) {
Polyhedron::Halfedge_handle h = f.halfedge();
return Polyhedron::Plane_3(h->vertex()->point(),
h->next()->vertex()->point(),
h->opposite()->vertex()->point());
}
};
inline static double planeDistance(Plane &p, Plane &q) {
double sc1 = max(abs(p.a()),
max(abs(p.b()),
max(abs(p.c()),
abs(p.d()))));
double sc2 = max(abs(q.a()),
max(abs(q.b()),
max(abs(q.c()),
abs(q.d()))));
Plane r(p.a() * sc2,
p.b() * sc2,
p.c() * sc2,
p.d() * sc2);
Plane s(q.a() * sc1,
q.b() * sc1,
q.c() * sc1,
q.d() * sc1);
return ((r.a() - s.a()) * (r.a() - s.a()) +
(r.b() - s.b()) * (r.b() - s.b()) +
(r.c() - s.c()) * (r.c() - s.c()) +
(r.d() - s.d()) * (r.d() - s.d())) / (sc1 * sc2);
}
static void detriangulatePolyhedron(Polyhedron &poly) {
vector<Polyhedron::Halfedge_handle> toJoin;
for (auto edge = poly.edges_begin(); edge != poly.edges_end(); edge++) {
auto f1 = edge->facet();
auto f2 = edge->opposite()->facet();
if (planeDistance(f1->plane(), f2->plane()) < 1E-5) {
toJoin.push_back(edge);
}
}
for (auto edge = toJoin.begin(); edge != toJoin.end(); edge++) {
poly.join_facet(*edge);
}
}
...
Polyhedron convexHull;
CGAL::convex_hull_3(shape.begin(),
shape.end(),
convexHull);
transform(convexHull.facets_begin(),
convexHull.facets_end(),
convexHull.planes_begin(),
Plane_from_facet());
detriangulatePolyhedron(convexHull);
Plane bounds[convexHull.size_of_facets()];
int boundCount = 0;
for (auto facet = convexHull.facets_begin(); facet != convexHull.facets_end(); facet++) {
bounds[boundCount++] = facet->plane();
}
...
这给出了预期的结果(前后):
我正在制作一个分析几何的小应用程序。在我的程序的一部分中,我使用 has 的算法将凸对象作为输入。幸运的是,我所有的对象最初都是凸的,但有些只是 勉强 所以(见图)。
在我应用一些转换后,我的算法无法工作(它产生 "infinitely" 长多边形等),我认为这是因为图像中的舍入误差;由于舍入误差(在图像中非常夸张),圆柱体的顶部顶点略有 "pushed in",不再凸起。
所以我的问题是:有人知道 "slightly convexify" 对象的方法吗?这是我尝试实现的一种方法,但它似乎不起作用(或者我实现错误):
1. Average all vertices together to create a vertex C inside the convex shape.
2. Let d[v] be the distance from C to vertex v.
3. Scale each vertex v from the center C with the scale factor 1 / (1+d[v] * CONVEXIFICATION_FACTOR)
谢谢!!我安装了 CGAL 和 Boost,所以我可以使用这些库函数中的任何一个(我已经这样做了)。
您当然可以通过计算对象的 convex hull 使对象成为凸面。但这将 "convexify" 任何东西。如果您确定您的输入只是稍微偏离凸面,那么这应该不是问题。
CGAL 似乎在其中实现了 3D Quickhull,这将是第一个尝试的事情。有关文档和一些示例程序,请参阅 http://doc.cgal.org/latest/Convex_hull_3/。 (我对 CGAL 不够熟悉,无法重现任何示例并声称它们是正确的。)
最后我发现这个问题的根源是凸包包含很多三角形,而我的输入形状通常是立方体,使得每个四边形区域显示为 2 个具有极其相似平面的三角形方程式,导致我使用的算法出现某种问题。
我用 "detriangulating" 多面体解决了它,使用这段代码。如果有人能发现任何改进或问题,请告诉我!
#include <algorithm>
#include <cmath>
#include <vector>
#include <CGAL/convex_hull_traits_3.h>
#include <CGAL/convex_hull_3.h>
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Aff_transformation_3 Transformation;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
struct Plane_from_facet {
Polyhedron::Plane_3 operator()(Polyhedron::Facet& f) {
Polyhedron::Halfedge_handle h = f.halfedge();
return Polyhedron::Plane_3(h->vertex()->point(),
h->next()->vertex()->point(),
h->opposite()->vertex()->point());
}
};
inline static double planeDistance(Plane &p, Plane &q) {
double sc1 = max(abs(p.a()),
max(abs(p.b()),
max(abs(p.c()),
abs(p.d()))));
double sc2 = max(abs(q.a()),
max(abs(q.b()),
max(abs(q.c()),
abs(q.d()))));
Plane r(p.a() * sc2,
p.b() * sc2,
p.c() * sc2,
p.d() * sc2);
Plane s(q.a() * sc1,
q.b() * sc1,
q.c() * sc1,
q.d() * sc1);
return ((r.a() - s.a()) * (r.a() - s.a()) +
(r.b() - s.b()) * (r.b() - s.b()) +
(r.c() - s.c()) * (r.c() - s.c()) +
(r.d() - s.d()) * (r.d() - s.d())) / (sc1 * sc2);
}
static void detriangulatePolyhedron(Polyhedron &poly) {
vector<Polyhedron::Halfedge_handle> toJoin;
for (auto edge = poly.edges_begin(); edge != poly.edges_end(); edge++) {
auto f1 = edge->facet();
auto f2 = edge->opposite()->facet();
if (planeDistance(f1->plane(), f2->plane()) < 1E-5) {
toJoin.push_back(edge);
}
}
for (auto edge = toJoin.begin(); edge != toJoin.end(); edge++) {
poly.join_facet(*edge);
}
}
...
Polyhedron convexHull;
CGAL::convex_hull_3(shape.begin(),
shape.end(),
convexHull);
transform(convexHull.facets_begin(),
convexHull.facets_end(),
convexHull.planes_begin(),
Plane_from_facet());
detriangulatePolyhedron(convexHull);
Plane bounds[convexHull.size_of_facets()];
int boundCount = 0;
for (auto facet = convexHull.facets_begin(); facet != convexHull.facets_end(); facet++) {
bounds[boundCount++] = facet->plane();
}
...
这给出了预期的结果(前后):