Geometry rounding problems: object no longer convex after simple transformations

Question

I'm making a little app to analyze geometry. In one part of my program, I use an algorithm that has to have a convex object as input. Luckily, all my objects are initially convex, but some are just barely so (see image).

After I apply some transformations, my algorithm fails to work (it produces "infinitely" long polygons, etc), and I think this is because of rounding errors as in the image; the top vertex in the cylinder gets "pushed in" slightly because of rounding errors (very exaggerated in image) and is no longer convex.

cylinder after and before transformation

So my question is: Does anyone know of a method to "slightly convexify" an object? Here's one method I tried to implement but it didn't seem to work (or I implemented it wrong):

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)

Thanks!! I have CGAL and Boost installed so I can use any of those library functions (and I already do).

Jay Kominek · Answer 1 · 2015-05-29T14:32:50.350

3

You can certainly make the object convex by computing the convex hull of it. But that'll "convexify" anything. If you're sure your input has departed only slightly from being convex, then it shouldn't be a problem.

CGAL appears to have an implementation of 3D Quickhull in it, which would be the first thing to try. See http://doc.cgal.org/latest/Convex_hull_3/ for docs and some example programs. (I'm not sufficiently familiar with CGAL to want to reproduce any examples and claim they're correct.)

edited May 29 '15 at 14:32

answered May 29 '15 at 14:23

Jay Kominek

8,674
1
34
51

Thanks for the tip! I implemented what you suggested but the problem seems to be that even the function `convex_hull_3` doe snot produce a *strongly* convex hull (I've tested this with `is_strongly_convex_3` and the problem repeats. If only there was a way to get a strongly convex hull. – bombax May 31 '15 at 22:51
1

Looking at the definition of `is_strongly_convex_3` and considering your polyhedron, I'd look for the shared vertices of adjacent, coplanar facets. Once I found some of those vertices, I'd start with the "most shared" and try pushing them out from the shared plane of the facets by tiny (1e-14) amounts and see if the polyhedron becomes strongly convex. – Jay Kominek Jun 01 '15 at 20:03
I've just discovered that even the *strongly convex* shapes create the problem as above. Maybe the rounding error is too big for the algorithm or maybe there is another problem. I'll investigate and then mark as "correct" once I discover the problem! – bombax Jun 02 '15 at 20:26

score 1 · Accepted Answer · answered Jun 03 '15 at 00:44

In the end I discovered the root of this problem was the fact that the convex hull contained lots of triangles, whereas my input shapes were often cube-shaped, making each quadrilateral region appear as 2 triangles which had extremely similar plane equations, causing some sort of problem in the algorithm I was using.

I solved it by "detriangulating" the polyhedra, using this code. If anyone can spot any improvements or problems, let me know!

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

This gave the desired result (after and before):

after before

Geometry rounding problems: object no longer convex after simple transformations

2 Answers2