Surface Reconstruction By VTK And CGAL

Published by StephenWei on

Point cloud data:

input.txt

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:

Categories: CGALVTK
Tags:
0 0 votes
Article Rating
Subscribe
Notify of
guest

0 Comments
Inline Feedbacks
View all comments

Related Posts

VTK

VTK – How to make top renderer background color visible?

I added two renderers in renderWindow, and set a blue background for the top renderer. But the color is not visible, here is source code. We can find the description for the usage of the Read more…

VTK

Change ColorLevel Of vtkImageMapper To Adjust The Exposure

If you have used vtkImageActor or vtkImageMapper + vtkActor2D to display an image in VTK renderer window, you can change the value of ColorLevel of vtkImageMapper to adjust the exposure to get a better effect. Read more…

VTK

Make vtkImageActor Object Always Face Camera And Be Under vtkFollower Object

Here is a simple example about make vtkImageActor always face current active camera and keep it under a vtkFollower object. It’s impossible to send image data to vtkFollower directly, so we have to add the Read more…

Tex To PDF
 : convert the Latex file which suffix is tex to a PDF file

Got it!
X
0
Would love your thoughts, please comment.x
()
x