approximate_triangle_mesh() causes some vertices to shoot off beyond original box

[dmitrykudinov]

Please use the following template to help us solving your issue.

Issue Details

When applying VSA::approximate_triangle_mesh() on a triangular mesh, I'm getting some of the resulting vertices shooting off the original geometry. I tried using PMP::orient_polygon_soup() as a preprocessing step, but it made no difference.

Before VSA:

After (look at the far-left wall/corner) :

min_error_drop = 0.03
number_of_iterations = 40
subdivision_ratio = 3.

Source OBJ can be found here

Source Code

// requires CGAL 5.3+

#include <iostream>
#include <fstream>
#include <string.h>
#include <stdio.h>
#include <dirent.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_approximation/approximate_triangle_mesh.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>

#include <CGAL/IO/OBJ_reader.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;

namespace VSA = CGAL::Surface_mesh_approximation;
namespace PMP = CGAL::Polygon_mesh_processing;

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using namespace std;

Mesh read_mesh(string fname)
{
    Mesh mesh;
    // read input surface triangle mesh
    vector<K::Point_3> verts;
    vector<vector<size_t> > faces;
    std::ifstream in(fname);
    if(!in || !CGAL::read_OBJ(in, verts, faces))
    {
        cout << "!!! error reading " + fname << endl;
        return mesh;
    }
//    PMP::orient_polygon_soup(verts, faces);
    PMP::polygon_soup_to_polygon_mesh(verts, faces, mesh);
    return mesh;
}

void write_mesh(string fname, Mesh &mesh)
{
    // saving OBJ. Requires CGAL 5.3+
    vector<K::Point_3> verts;
    vector<vector<size_t> > faces;
    PMP::polygon_mesh_to_polygon_soup(mesh, verts, faces);
    CGAL::write_OBJ(fname, verts, faces);
}

Mesh approximate_bld3d_mesh(
    Mesh &mesh_in,
    float min_error_drop,
    int number_of_iterations,
    float subdivision_ratio
)
{
    Mesh mesh_out;
    vector<Kernel::Point_3> verts;
    vector<array<size_t, 3> > faces; // triplets of indices (always triangles)

    VSA::approximate_triangle_mesh(mesh_in,
        CGAL::parameters::min_error_drop(min_error_drop). // seeding with minimum error drop
        number_of_iterations(number_of_iterations). // set number of clustering iterations after seeding
        subdivision_ratio(subdivision_ratio). // set chord subdivision ratio threshold when meshing
        anchors(back_inserter(verts)). // output vertices
        triangles(back_inserter(faces))// output indexed faces
    );

    // convert from soup to surface mesh
    PMP::orient_polygon_soup(verts, faces);
    PMP::polygon_soup_to_polygon_mesh(verts, faces, mesh_out);
    cout << "v: " << verts.size() << "\tf: " << faces.size() << "\t";
    return mesh_out;
}

int main(int argc, char **argv)
{
    if (argc != 5 && argc != 2)
    {
        cout << "Usage:\n\tvsa <folder to OBJs> [<min_error_drop(0.03)> <number_of_iterations(10)> <subdivision_ratio(1.)>]" << endl;
        return EXIT_SUCCESS;
    }
    const char* folder_in = argv[1];
    float min_error_drop = argc==5? atof(argv[2]): 0.03;
    int number_of_iterations = argc==5? atoi(argv[3]): 10;
    float subdivision_ratio = argc==5? atof(argv[4]): 1.;

    struct dirent *entry = nullptr;
    DIR *dp = nullptr;
    dp = opendir(folder_in);
    string fname_in, fname_out, fname_ext;
    Mesh mesh_in, mesh_out;

    if (dp != nullptr)
         while ((entry = readdir(dp)))
        {
            fname_in = entry->d_name;
            fname_ext =  fname_in.length()>4? fname_in.substr(fname_in.length() - 4) : "";
            if (fname_ext == ".obj" || fname_ext == ".OBJ")
                try
                {
                    fname_in = string(folder_in) + (folder_in[strlen(folder_in)-1]=='/'? "": "/") + fname_in;
                    fname_out = fname_in.substr(0, fname_in.length() - 4) + "-vsa.obj";
                    //reading mesh from file
                    mesh_in = read_mesh(fname_in);

                    //mesh approximation
                    mesh_out = approximate_bld3d_mesh(
                        mesh_in,
                        min_error_drop,
                        number_of_iterations,
                        subdivision_ratio
                    );
                    //saving resulting mesh
                    write_mesh(fname_out, mesh_out);
                    cout << fname_out << endl;
                }
                catch(const exception& e)
                {
                    cout << "!!! error processing " << fname_in << ": " << e.what() << endl;
                }
        }

    closedir(dp);

    return EXIT_SUCCESS;
}

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 16.04 LTS
• Compiler: gcc 9.3.0
• Release or debug mode: release
• Specific flags used (if any): -
• CGAL version: seen on 5.2.1, 5.3-i-92
• Boost version: 1.69.0
• Other libraries versions if used (Eigen, TBB, etc.):

[sloriot]

I think this is the "expected" behavior. In the original publication you see: Since the proxies can be seen as approximate faces of the final mesh, we must put vertices at the intersection of the proxies. Thus, we create an anchor vertex at every original vertex where three or more regions meet. In order to account for every region, we then check whether each region boundary has at least three anchor vertices attached to it; if not, we simply add anchor vertices accordingly as it will guarantee the presence of at least one polygon per region. The spatial position of these anchor vertices is determined as follows: foreach neighboring proxy of an anchor, we compute the projection of the associated vertex from which the anchor was created onto the proxy(i.e.,its ideal position for this proxy); we then simply average these projections.

You have such an issue because you have several sheets at the corner with many self-intersections that the package try to approximate:

[dmitrykudinov]

Thank you for reviewing the case. Do you think CGAL::Polygon_mesh_processing::smooth_mesh(), as a preprocessing step to VSA, may help here? Sorry, cannot check it myself easily - for some reason, having a hard time making Ceres/Eigen being properly picked up on my ubuntu 16.04.

A bit more context here: we are evaluating CGAL's VSA implementation as an alternative to our in-house solution, potentially acquiring a commercial CGAL license.

In general, the self-overlapping in our input meshes cannot be completely avoided as they are generated by a neural network. We need a robust and efficient solution allowing to reduce the number of vertices through fitting approximating planar proxies, yet without the shoot-offs like in the above case.

[dmitrykudinov]

Actually, figured out the Ceres issue and tried smooth_mesh() - no luck here either. I guess one way to get rid of self-intersection could be through resampling and mesh reconstruction from point cloud as a preprocessing before VSA. Or is there a better CGAL-based way to get rid of self-intersections?