几何舍入问题:对象在简单变换后不再凸

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();
                    }
...

这给出了预期的结果(前后):