Проблемы округления геометрии: объект больше не выпуклый после простых преобразований
Я делаю небольшое приложение для анализа геометрии. В одной части моей программы я использую алгоритм, который имеет иметь выпуклый объект в качестве входных данных. К счастью, все мои объекты изначально выпуклые, но некоторые просто едва так (см. изображение).
После того, как я применил некоторые преобразования, мой алгоритм не работает (он производит «бесконечно» длинные полигоны и т. Д.), И я думаю, что это из-за ошибок округления, как на изображении; верхняя вершина в цилиндре слегка «вдавливается» из-за ошибок округления (очень преувеличена на изображении) и больше не является выпуклой.
Итак, мой вопрос: кто-нибудь знает метод «слегка выпуклости» объекта? Вот один метод, который я пытался реализовать, но он, похоже, не работал (или я реализовал его неправильно):
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, поэтому я могу использовать любую из этих библиотечных функций (и я уже это делаю).
Решение
Вы, конечно, можете сделать объект выпуклым, вычислив выпуклый корпус этого Но это «выпукло» все. Если вы уверены, что ваш вклад немного отклонился от выпуклости, то это не должно быть проблемой.
Кажется, в CGAL есть реализация 3D Quickhull, которая будет первой попыткой. Увидеть http://doc.cgal.org/latest/Convex_hull_3/ для документов и некоторых примеров программ. (Я недостаточно знаком с CGAL, чтобы захотеть воспроизвести какие-либо примеры и утверждать, что они верны.)
Другие решения
В конце концов я обнаружил, что корень этой проблемы заключался в том, что выпуклая оболочка содержала множество треугольников, тогда как мои входные формы часто имели форму куба, в результате чего каждая четырехугольная область представлялась в виде 2 треугольников, которые имели чрезвычайно похожие плоские уравнения, вызывая некоторый вид проблемы в алгоритме, который я использовал.
Я решил это, «детриангулируя» многогранники, используя этот код. Если кто-то может заметить какие-либо улучшения или проблемы, дайте мне знать!
#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();
}
...
Это дало желаемый результат (после и до):
