# Fast triangulation of unordered point clouds

This tutorial explains how to run a greedy surface triangulation algorithm on aPointCloud with normals, to obtain a triangle mesh based on projections of thelocal neighborhoods. An example of the method’s output can be seen here:

# Background: algorithm and parameters

The method works by maintaining a list of points from which the mesh can begrown (“fringe” points) and extending it until all possible points areconnected. It can deal with unorganized points, coming from one or multiplescans, and having multiple connected parts. It works best if the surface islocally smooth and there are smooth transitions between areas with differentpoint densities.

Triangulation is performed locally, by projecting the local neighborhood of apoint along the point’s normal, and connecting unconnected points. Thus, thefollowing parameters can be set:

• setMaximumNearestNeighbors(unsigned) and setMu(double) control the size ofthe neighborhood. The former defines how many neighbors are searched for,while the latter specifies the maximum acceptable distance for a point to beconsidered, relative to the distance of the nearest point (in order to adjustto changing densities). Typical values are 50-100 and 2.5-3 (or 1.5 forgrids).
• setSearchRadius(double) is practically the maximum edge length for everytriangle. This has to be set by the user such that to allow for the biggesttriangles that should be possible.
• setMinimumAngle(double) and setMaximumAngle(double) are the minimum andmaximum angles in each triangle. While the first is not guaranteed, thesecond is. Typical values are 10 and 120 degrees (in radians).
• setMaximumSurfaceAgle(double) and setNormalConsistency(bool) are meant todeal with the cases where there are sharp edges or corners and where twosides of a surface run very close to each other. To achieve this, points arenot connected to the current point if their normals deviate more than thespecified angle (note that most surface normal estimation methods producesmooth transitions between normal angles even at sharp edges). This angle iscomputed as the angle between the lines defined by the normals (disregardingthe normal’s direction) if the normal-consistency-flag is not set, as not allnormal estimation methods can guarantee consistently oriented normals.Typically, 45 degrees (in radians) and false works on most datasets.

Please see the example below, and you can consult the following paper and itsreferences for more details:

@InProceedings{Marton09ICRA,
author    = {Zoltan Csaba Marton and Radu Bogdan Rusu and Michael Beetz},
title     = {{On Fast Surface Reconstruction Methods for Large and Noisy Datasets}},
booktitle = {Proceedings of the IEEE International Conference on Robotics and Automation (ICRA)},
month     = {May 12-17},
year      = {2009},
}

# The code

First, create a file, let’s say, greedy_projection.cpp in your favoriteeditor, and place the following code inside it:

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/normal_3d.h>
#include <pcl/surface/gp3.h>
#include <pcl/io/vtk_io.h>

int
main (int argc, char** argv)
{
// Load input file into a PointCloud<T> with an appropriate type
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PCLPointCloud2 cloud_blob;
pcl::fromPCLPointCloud2 (cloud_blob, *cloud);
//* the data should be available in cloud

// Normal estimation*
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> n;
pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal>);
pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
tree->setInputCloud (cloud);
n.setInputCloud (cloud);
n.setSearchMethod (tree);
n.setKSearch (20);
n.compute (*normals);
//* normals should not contain the point normals + surface curvatures

// Concatenate the XYZ and normal fields*
pcl::PointCloud<pcl::PointNormal>::Ptr cloud_with_normals (new pcl::PointCloud<pcl::PointNormal>);
pcl::concatenateFields (*cloud, *normals, *cloud_with_normals);
//* cloud_with_normals = cloud + normals

// Create search tree*
pcl::search::KdTree<pcl::PointNormal>::Ptr tree2 (new pcl::search::KdTree<pcl::PointNormal>);
tree2->setInputCloud (cloud_with_normals);

// Initialize objects
pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;
pcl::PolygonMesh triangles;

// Set the maximum distance between connected points (maximum edge length)

// Set typical values for the parameters
gp3.setMu (2.5);
gp3.setMaximumNearestNeighbors (100);
gp3.setMaximumSurfaceAngle(M_PI/4); // 45 degrees
gp3.setMinimumAngle(M_PI/18); // 10 degrees
gp3.setMaximumAngle(2*M_PI/3); // 120 degrees
gp3.setNormalConsistency(false);

// Get result
gp3.setInputCloud (cloud_with_normals);
gp3.setSearchMethod (tree2);
gp3.reconstruct (triangles);

std::vector<int> parts = gp3.getPartIDs();
std::vector<int> states = gp3.getPointStates();

pcl::io::saveVTKFile("../data/testmesh.vtk", triangles);
// Finish
return (0);
}

The input file you can find at pcl/test/bun0.pcd

# The explanation

Now, let’s break down the code piece by piece.

// Load input file into a PointCloud<T> with an appropriate type
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PCLPointCloud2 cloud_blob;
pcl::fromPCLPointCloud2 (cloud_blob, *cloud);
//* the data should be available in cloud

as the example PCD has only XYZ coordinates, we load it into aPointCloud<PointXYZ>.

// Normal estimation*
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> n;
pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal>);
pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
tree->setInputCloud (cloud);
n.setInputCloud (cloud);
n.setSearchMethod (tree);
n.setKSearch (20);
n.compute (*normals);
//* normals should not contain the point normals + surface curvatures

the method requires normals, so they are estimated using the standard methodfrom PCL.

 // Concatenate the XYZ and normal fields*
pcl::PointCloud<pcl::PointNormal>::Ptr cloud_with_normals (new pcl::PointCloud<pcl::PointNormal>);
pcl::concatenateFields (*cloud, *normals, *cloud_with_normals);
//* cloud_with_normals = cloud + normals

Since coordinates and normals need to be in the same PointCloud, we create a PointNormal type point cloud.

  // Create search tree*
pcl::search::KdTree<pcl::PointNormal>::Ptr tree2 (new pcl::search::KdTree<pcl::PointNormal>);
tree2->setInputCloud (cloud_with_normals);

// Initialize objects
pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;
pcl::PolygonMesh triangles;

The above lines deal with the initialization of the required objects.

// Set the maximum distance between connected points (maximum edge length)

// Set typical values for the parameters
gp3.setMu (2.5);
gp3.setMaximumNearestNeighbors (100);
gp3.setMaximumSurfaceAngle(M_PI/4); // 45 degrees
gp3.setMinimumAngle(M_PI/18); // 10 degrees
gp3.setMaximumAngle(2*M_PI/3); // 120 degrees
gp3.setNormalConsistency(false);

The above lines set the parameters, as explained above.

 // Get result
gp3.setInputCloud (cloud_with_normals);
gp3.setSearchMethod (tree2);
gp3.reconstruct (triangles);

The lines above set the input objects and perform the actual triangulation.

// Additional vertex information
std::vector<int> parts = gp3.getPartIDs();
std::vector<int> states = gp3.getPointStates();


for each point, the ID of the containing connected component and its “state”(i.e. gp3.FREE, gp3.BOUNDARY or gp3.COMPLETED) can be retrieved.

# Saving and viewing the result

You can view the smoothed cloud for example by saving into a VTK file by:

#include <pcl/io/vtk_io.h>
...
saveVTKFile ("mesh.vtk", triangles);

# Compiling and running the program

cmake_minimum_required(VERSION 2.8 FATAL_ERROR)

project(greedy_projection)

find_package(PCL 1.2 REQUIRED)

include_directories(${PCL_INCLUDE_DIRS}) link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS}) add_executable (greedy_projection greedy_projection.cpp) target_link_libraries (greedy_projection${PCL_LIBRARIES})


After you have made the executable, you can run it. Simply do:

\$ ./greedy_projection

09-06 421

09-09 101

#### segmentation using delaunay <em>Triangulation</em>

2017年06月02日 00:00

07-18 4363

11-06 1127

12-05 3960

01-15 306

10-15 688

10-13 458