Surface Reconstruction By VTK And CGAL

Point cloud data:

VTK

Use vtkSurfaceReconstructionFilter to get mesh from the point coordinates. The vtkSurfaceReconstructionFilter object generate image data, we have to use vtkContourFilter to find the surface.

CMakeLists.txt

cmake_minimum_required(VERSION 3.16)

project(VTKPractice)

find_package( VTK REQUIRED )
include( ${VTK_USE_FILE} )

add_executable(${PROJECT_NAME} "main.cpp" )

target_link_libraries( ${PROJECT_NAME} ${VTK_LIBRARIES} )

set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} /Zi")
set(CMAKE_EXE_LINKER_FLAGS_RELEASE "${CMAKE_EXE_LINKER_FLAGS_RELEASE} /DEBUG /OPT:REF /OPT:ICF")

main.cpp:

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkClipPolyData.h>
#include <vtkPlane.h>
#include <vtkPoints.h>
#include <vtkCleanPolyData.h>
#include <vtkDataSetMapper.h>
#include <vtkSTLWriter.h>
#include <vtkFillHolesFilter.h>
#include <vtkContourFilter.h>
#include <vtkSurfaceReconstructionFilter.h>
#include <vtkProperty.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkDelaunay3D.h>
#include <vtkFacetWriter.h>
#include <vtkFacetReader.h>
#include <vtkOBJExporter.h>
#include <vtkPolyData.h>
#include <vtkPolyDataWriter.h>
#include <vtkDataSetSurfaceFilter.h>

#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();

using namespace std;

int main()
{
    string filePath( "/Users/weiyang/Desktop/input.txt" );
    ifstream input;
    input.open( filePath.c_str(), ios::binary );
    if (input.is_open() == false)
    {
        cout << "open failed: " << filePath << endl;
        return -1;
    }
    vtkSmartPointer<vtkPoints> points =
            vtkSmartPointer<vtkPoints>::New();
    while( !input.eof() )
    {
        double origin[3];
        for( int i = 0; i < 3; ++i )
        {
            input >> origin[i];
        }
        points->InsertNextPoint( origin );
    }
    vtkSmartPointer<vtkPolyData> polyData =
            vtkSmartPointer<vtkPolyData>::New();
    polyData->SetPoints( points );


    // Construct the surface and create isosurface.
    vtkSmartPointer<vtkSurfaceReconstructionFilter> surf =
      vtkSmartPointer<vtkSurfaceReconstructionFilter>::New();
    surf->SetInputData( polyData );
    surf->Update();


    vtkSmartPointer<vtkContourFilter> contourFilter = vtkSmartPointer<vtkContourFilter>::New();
    contourFilter->SetInputConnection(surf->GetOutputPort());
    contourFilter->SetValue(0, 0.0);
    contourFilter->Update();


    vtkSPtrNew( mapper, vtkPolyDataMapper );
    mapper->SetInputData( contourFilter->GetOutput() );


    vtkSPtrNew( actor, vtkActor );
    actor->SetMapper( mapper );


    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor(actor);
    renderer->SetBackground( 0, 0, 0 );


    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );


    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
    renderWindowInteractor->SetRenderWindow( renderWindow );


    renderer->ResetCamera();
    renderWindow->Render();
    renderWindowInteractor->Start();


    return 0;
}

Result:

CGAL

I use the advancing front surface reconstruction way to get the mesh from point cloud information.

We have to support C++17 in the CGAL project, let’s configure it in the CMakeLists.txt.

CMakeLists.txt:

cmake_minimum_required(VERSION 3.1...3.15)
project(CGAL_Test)


find_package(CGAL REQUIRED)


# https://cmake.org/cmake/help/latest/prop_tgt/CXX_STANDARD.html
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)


# use Eigen
find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
include(CGAL_Eigen_support)


# Link with Boost.ProgramOptions (optional)
find_package(Boost QUIET COMPONENTS program_options)
if(Boost_PROGRAM_OPTIONS_FOUND)
  if(TARGET Boost::program_options)
    set(Boost_PROGRAM_OPTIONS_LIBRARY Boost::program_options)
  endif()
  if(CGAL_AUTO_LINK_ENABLED)
    message(STATUS "Boost.ProgramOptions library: found")
  else()
    message(
      STATUS "Boost.ProgramOptions library: ${Boost_PROGRAM_OPTIONS_LIBRARY}")
  endif()
  add_definitions("-DCGAL_USE_BOOST_PROGRAM_OPTIONS")
  list(APPEND CGAL_3RD_PARTY_LIBRARIES ${Boost_PROGRAM_OPTIONS_LIBRARY})
endif()

create_single_source_cgal_program( "main.cpp" )  #defined in CGAL_CreateSingleSourceCGALProgram.cmake

target_link_libraries(main PUBLIC CGAL::Eigen_support)

main.cpp:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/remove_outliers.h>
#include <CGAL/grid_simplify_point_set.h>
#include <CGAL/jet_smooth_point_set.h>
#include <CGAL/jet_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/poisson_surface_reconstruction.h>
#include <CGAL/Advancing_front_surface_reconstruction.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <cstdlib>
#include <vector>
#include <fstream>
#include <iostream>


// types
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Sphere_3 Sphere_3;
typedef CGAL::Point_set_3<Point_3, Vector_3> Point_set;
int main(int argc, char*argv[])
{
  Point_set points;
  if (argc < 2)
  {
    std::cerr << "Usage: " << argv[0] << " [input.xyz/off/ply/las]" << std::endl;
    return EXIT_FAILURE;
  }
  const char* input_file = argv[1];
  std::ifstream stream (input_file, std::ios_base::binary);
  if (!stream)
  {
    std::cerr << "Error: cannot read file " << input_file << std::endl;
    return EXIT_FAILURE;
  }
  stream >> points;
  std::cout << "Read " << points.size () << " point(s)" << std::endl;


  typedef CGAL::Point_set_3<Point_3, Vector_3>::const_iterator Iterator;
  for( Iterator it1 = points.begin(); it1!= points.end(); ++it1 )
  {
      const Point_3 &p1 = points.point( *it1 );
      std::cout << *it1 << ": " << p1 << std::endl;
  }


  if (points.empty())
    return EXIT_FAILURE;
  CGAL::remove_outliers<CGAL::Sequential_tag>
    (points,
     24, // Number of neighbors considered for evaluation
     points.parameters().threshold_percent (5.0)); // Percentage of points to remove
  std::cout << points.number_of_removed_points()
            << " point(s) are outliers." << std::endl;

  points.collect_garbage();

  // Compute average spacing using neighborhood of 6 points
  double spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag> (points, 6);

  // Simplify using a grid of size 2 * average spacing
  CGAL::grid_simplify_point_set (points, 2. * spacing);
  std::cout << points.number_of_removed_points()
            << " point(s) removed after simplification." << std::endl;
  points.collect_garbage();
  CGAL::jet_smooth_point_set<CGAL::Sequential_tag> (points, 24);
  unsigned int reconstruction_choice
    = (argc < 3 ? 0 : atoi(argv[2]));
  if (reconstruction_choice == 1) // Advancing front
  {
    typedef std::array<std::size_t, 3> Facet; // Triple of indices
    std::vector<Facet> facets;
    // The function is called using directly the points raw iterators
    CGAL::advancing_front_surface_reconstruction(points.points().begin(),
                                                 points.points().end(),
                                                 std::back_inserter(facets));
    std::cout << facets.size ()
              << " facet(s) generated by reconstruction." << std::endl;
    // copy points for random access
    std::vector<Point_3> vertices;
    vertices.reserve (points.size());
    std::copy (points.points().begin(), points.points().end(), std::back_inserter (vertices));
    CGAL::Surface_mesh<Point_3> output_mesh;
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh (vertices, facets, output_mesh);
    std::ofstream f ("out.off");
    f << output_mesh;
    f.close ();
  }
  else // Handle error
  {
    std::cerr << "Error: invalid reconstruction id: " << reconstruction_choice << std::endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}

The code comes from the tutorial, https://doc.cgal.org/latest/Manual/tuto_reconstruction.html#TutorialsReconstruction_algorithms

Build:

build > find /home/stephen/Downloads/CGAL-5.5/lib/cmake/CGAL/ -name CGALConfig.cmake
/home/stephen/Downloads/CGAL-5.5/lib/cmake/CGAL/CGALConfig.cmake
build > cmake  -DCGAL_DIR=/home/stephen/Downloads/CGAL-5.5/lib/cmake/CGAL/ ../      
-- Using header-only CGAL
-- Targetting Unix Makefiles
-- Using /usr/bin/c++ compiler.

Run it：

./main /Users/weiyang/Desktop/input.txt 1

We got the file out.off and convert it to out2.off to display the mesh easily.

out.off   
out2.off

Show the mesh after reading the point cloud from the file out2.off.

New project for display the result.

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkActor.h>
#include <vtkConeSource.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkClipPolyData.h>
#include <vtkPlane.h>
#include <vtkPoints.h>
#include <vtkPolyDataMapper.h>
#include <vtkDataSetMapper.h>
#include <vtkSTLWriter.h>
#include <vtkFillHolesFilter.h>
#include <vtkProperty.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkDelaunay3D.h>
#include <vtkFacetWriter.h>
#include <vtkFacetReader.h>
#include <vtkOBJExporter.h>
#include <vtkPolyData.h>
#include <vtkPolyDataWriter.h>


using namespace std;


#define vtkSPtr vtkSmartPointer
#define vtkSPtrNew(Var, Type) vtkSPtr<Type> Var = vtkSPtr<Type>::New();


int main()
{
    string filePath( "/Users/weiyang/Desktop/out2.off" );
    ifstream input;
    input.open( filePath.c_str(), ios::binary );
    if (input.is_open() == false)
    {
        cout << "open failed: " << filePath << endl;
        return -1;
    }
    int pointCount, cellCount, normalCount;
    input >> pointCount;
    input >> cellCount;
    input >> normalCount;


    vtkSmartPointer<vtkPoints> points =
            vtkSmartPointer<vtkPoints>::New();
    for( int i = 0; i < pointCount; ++i )
    {
        double value[3];
        for( int j = 0; j < 3; ++j )
        {
            input >> value[j];
        }
        points->InsertNextPoint( value );
    }
    vtkSmartPointer<vtkPolyData> polyData =
            vtkSmartPointer<vtkPolyData>::New();
    polyData->SetPoints( points );


    int ptCountInCell, tmp;
    vtkSmartPointer<vtkCellArray> cellArray = vtkSmartPointer<vtkCellArray>::New();
    for( int i = 0; i < cellCount; ++i )
    {
        input >> ptCountInCell;


        vtkIdType pts[3];
        for( int j = 0; j < ptCountInCell; ++j )
        {
            input >> tmp;
            pts[j] = tmp;
        }
        cellArray->InsertNextCell( 3, pts );
    }
    polyData->SetPolys( cellArray );
    polyData->Modified();


    vtkSPtrNew( mapper, vtkPolyDataMapper );
    mapper->SetInputData( polyData );


    vtkSPtrNew( actor, vtkActor );
    actor->SetMapper( mapper );


    vtkSPtrNew( renderer, vtkRenderer );
    renderer->AddActor(actor);
    renderer->SetBackground( 0, 0, 0 );


    vtkSPtrNew( renderWindow, vtkRenderWindow );
    renderWindow->AddRenderer( renderer );


    vtkSPtrNew( renderWindowInteractor, vtkRenderWindowInteractor );
    renderWindowInteractor->SetRenderWindow( renderWindow );


    renderer->ResetCamera();
    renderWindow->Render();
    renderWindowInteractor->Start();


    return 0;
}

Result: