Unexpected shape with 3D alpha shapes

[palanglois]

Hi,

Issue Details

I try to make an alpha shape approximation of a mesh. In order to do so, I used the following steps:

1. Loading the mesh (which is a connected set) in CGAL
2. Performing a dense sampling of the mesh
3. Building an alpha shape with the obtained point cloud.

Here are the representations I get at each step.

Original mesh:

Sampled point cloud:

Alpha shape:

As you can see, despite the dense sampling, I get a very raw approximation in the concave parts of my mesh. Is this result normal ? I expected the alpha shape to almost perfectly fit the original mesh.

Source Code

The mesh I'm using: model.obj

CMakeLists.txt:

cmake_minimum_required(VERSION 3.1)
project(alphashapemin)

# CGAL
find_package(CGAL REQUIRED)
include(${CGAL_USE_FILE})

# Path to the data
add_definitions(-DDATA_PATH="${CMAKE_CURRENT_SOURCE_DIR}/model.obj")

add_executable(alphashapemin main.cpp)
target_link_libraries(alphashapemin ${CGAL_LIBRARIES})

main.cpp

// STD
#include <iostream>
#include <fstream>
#include <vector>

// CGAL
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/IO/OBJ_reader.h>
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Delaunay_triangulation_3.h>

// CGAL kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

// Geometric primitives
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Triangle_3 Triangle;

// Alpha shape
typedef CGAL::Alpha_shape_vertex_base_3<Kernel>          Vb;
typedef CGAL::Alpha_shape_cell_base_3<Kernel>            Fb;
typedef CGAL::Triangulation_data_structure_3<Vb,Fb>  Tds;
typedef CGAL::Delaunay_triangulation_3<Kernel,Tds>       Triangulation_3;
typedef CGAL::Alpha_shape_3<Triangulation_3>         Alpha_shape_3;
typedef Alpha_shape_3::Alpha_iterator                Alpha_iterator;
typedef Alpha_shape_3::Facet                             Facet;

using namespace std;

vector<Point> sampleOnTriangles(vector<Triangle> triangles, int nbSamples)
{
    // Sampling
    vector<Point> sampledPoints;
    vector<double> cumulatedAreas;
    double accum = 0;
    // Build the cumulated areas histogram
    for(const auto &triangle: triangles)
    {
        accum += sqrt(CGAL::to_double(triangle.squared_area()));
        cumulatedAreas.push_back(accum);
    }
    // Normalize it
    for(size_t i = 0; i < cumulatedAreas.size(); i++)
        cumulatedAreas[i] /= accum;
    // Actual sampling
    for(int i = 0; i < nbSamples; i++)
    {
        // Select a random triangle according to the areas distribution
        double r = ((double) rand() / (RAND_MAX));
        size_t found_index = 0;
        for(size_t j=0; j<cumulatedAreas.size() && r > cumulatedAreas[j]; j++)
            found_index = j+1;

        // Draw a random point in this triangle
        double r1 = ((double) rand() / (RAND_MAX));
        double r2 = ((double) rand() / (RAND_MAX));
        Point A = triangles[found_index].vertex(0);
        Point B = triangles[found_index].vertex(1);
        Point C = triangles[found_index].vertex(2);
        Point P = CGAL::ORIGIN + (1 - sqrt(r1)) * (A - CGAL::ORIGIN) +
                  (sqrt(r1) * (1 - r2)) * (B - CGAL::ORIGIN) + (sqrt(r1) * r2) * (C - CGAL::ORIGIN);

        sampledPoints.push_back(P);
    }
    return sampledPoints;
}

void saveAlphaShape(Alpha_shape_3 &as, string path)
{
    vector<Facet> facets;
    as.get_alpha_shape_facets(std::back_inserter(facets),
                              Alpha_shape_3::REGULAR);
    std::size_t nbf=facets.size();
    vector<Point> outPoints;
    vector<vector<size_t>> outFaces;
    for (std::size_t i=0;i<nbf;++i)
    {
        //To have a consistent orientation of the facet, always consider an exterior cell
        if ( as.classify( facets[i].first )!=Alpha_shape_3::EXTERIOR )
            facets[i]=as.mirror_facet( facets[i] );
        CGAL_assertion(  as.classify( facets[i].first )==Alpha_shape_3::EXTERIOR  );

        int indices[3]={
                (facets[i].second+1)%4,
                (facets[i].second+2)%4,
                (facets[i].second+3)%4,
        };

        // according to the encoding of vertex indices, this is needed to get
        // a consistent orienation
        if ( facets[i].second%2==0 ) std::swap(indices[0], indices[1]);

        outPoints.push_back(facets[i].first->vertex(indices[0])->point());
        outPoints.push_back(facets[i].first->vertex(indices[1])->point());
        outPoints.push_back(facets[i].first->vertex(indices[2])->point());

        vector<size_t> curFace = {3*i, 3*i+1, 3*i+2};
        outFaces.push_back(curFace);
    }

    ofstream outFile(path);
    for(const auto &point: outPoints)
        outFile << "v " << point.x() << " " << point.y() << " " << point.z() << endl;

    for(const auto &face: outFaces)
        outFile << "f " << face[0]+1 << " " << face[1]+1 << " " << face[2]+1 << endl;
    outFile.close();

}

int main() {

    //Read obj input file
    vector<Point> points;
    vector<vector<size_t> > faces;

    ifstream inStream(DATA_PATH);
    read_OBJ(inStream, points, faces);
    vector<Triangle> meshTriangles;
    for(const auto face: faces)
        meshTriangles.emplace_back(points[face[0]], points[face[1]], points[face[2]]);
    vector<Point> sampledPoints = sampleOnTriangles(meshTriangles, 400000);

    // Output sampled point cloud
    ofstream outFile("sample.obj");
    for(const auto &point: sampledPoints)
        outFile << "v " << point.x() << " " << point.y() << " " << point.z() << endl;
    outFile.close();

    // Compute alpha shape
    Alpha_shape_3 target(sampledPoints.begin(),sampledPoints.end());
    Alpha_iterator opts = target.find_optimal_alpha(1);
    cout << "alpha: " << *opts << endl;
    target.set_alpha(*opts);

    // Save alpha shape
    saveAlphaShape(target, "target.obj");
    return 0;
}

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 16.04 64 bits
• Compiler: gcc 4.9
• Release or debug mode: Release
• Specific flags used (if any):
• CGAL version: 5.0.2
• Boost version: 1.58
• Other libraries versions if used (Eigen, TBB, etc.):

[MaelRL]

The find_optimal_alpha() is based on the regularized alpha shape. In the regularized alpha shapes, k-simplicies (vertices, edges, faces) which belong to the alpha shape but do not have an incident (k+1)-simplex (respectively edges, faces, cells) in the alpha shape are ignored. Therefore the optimal alpha happens later than what you expected here.

What you are interested in is more like a general alpha shapes, such that faces in the alpha complex are not removed even if there is no incident cells. You will then get something like:

This is generated using the 3D alpha shapes CGAL demo, which has a useful slider for alpha value. I modified it to draw faces that are singular in addition to the regular ones (see https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html#ae887e9ecdaca28790f6ccdd73b84e40c) and hand picked a sensible alpha value.

I am not exactly sure how you could deduce which alpha value best fits your need for reconstruction purposes. The CGAL package Scale space reconstruction seems to be doing what you would like to do. Maybe @sgiraudot can expand on that as I am not very familiar with surface reconstruction.

[sgiraudot]

The alpha shape mesher of Scale Space relies on a scale parameter, similarly to the Alpha Shape themselves in CGAL, so it won't help.

My advice would be to use something like CGAL::compute_average_spacing() on the point cloud to have an estimation on the spacing between points, and then use a small factor (x2 or x3) or this value for alpha.