Fast Triangulation of Unordered Point Clouds (OpenMP Version)#

This tutorial explains how to run a greedy surface triangulation algorithm (OpenMP version) on a PointCloud with normals, to obtain a triangle mesh based on projections of the local neighborhoods. The OpenMP version divides the point cloud into multiple segments and invokes the original version of the greedy surface triangulation based on the segment count. By running the segments independently and concurrently, there will be a disjoint result between each segment. To overcome the disjoint issue, additional overlap segments are created to go through greedy surface triangulation and combine with other processed segments.

Note

The output of OpenMP version may differ from serial version of greedy projection triangulation.

Background: Algorithm and Parameters#

Refer to original greedy_triangulation for more detail original greedy triangulation parameters. All greedy triangulation parameters are supported except getPointStates(), getPartIDs(), getSFN() and getFFN().

Additional parameters are defined to provide control segments and overlap.

setNumberofThreads(unsigned) controls how many segments are created from the input point cloud. Each segment size is determined using the formula (maximum x - minimum x) / number of threads, and all segment sizes are equal. Each segment will be assigned to a CPU thread/core. It is recommended to obtain the number of threads/cores using omp_get_max_threads().
setBlockOverlapPercentage(double) controls how wide to create the overlap region among two segments. The overlap region is this parameter percentage multiply the size of segment. If there is still gap between segments, can overcome by increase the overlap percentage.
setRemoveDuplicateMesh(bool) controls whether to clean up the duplicate point cloud and polygons after combining all the segments output from GP3.

Note

This tutorial is applicable for execution for both within inside and outside a Docker image. It assumes that the pcl-oneapi-tutorials Deb package is installed, and the user has copied the tutorial directory from /opt/intel/pcl/oneapi/tutorials/ to a user-writable directory.

Prepare the environment:

cd <path-to-oneapi-tutorials>/greedy_projection

greedy_projectoin.cpp should be in the directory with following content:

// SPDX-License-Identifier: Apache-2.0
// Copyright (C) 2025 Intel Corporation
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/io/vtk_io.h>
#include <pcl/search/kdtree.h> // for KdTree
#include <pcl/features/normal_3d.h>
#include <pcl/oneapi/surface_omp/gp3.h>

int
main ()
{
  // 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::io::loadPCDFile ("bun0.pcd", 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::GreedyProjectionTriangulationOMP<pcl::PointNormal> gp3;
  pcl::PolygonMesh triangles;

  // Set the maximum distance between connected points (maximum edge length)
  gp3.setSearchRadius (0.025);

  // 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);
  gp3.setNumberOfThreads(3);
  gp3.setBlockOverlapPercentage(0.04);
  gp3.setRemoveDuplicateMesh(true);

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

  pcl::io::saveVTKFile("mesh.vtk", triangles);

  // Finish
  return (0);
}

Source the Intel® oneAPI Base Toolkit environment:
```
source /opt/intel/oneapi/setvars.sh
```

(Optional) Set up proxy setting to download test data:

export http_proxy="http://<http_proxy>:port"
export https_proxy="http://<https_proxy>:port"

Build the code:

mkdir build && cd build
cmake ../
make -j

Run the binary:
```
./greedy_projection
```
Output:

The program will save the output as “mesh.vtk”. View the VTK file by:
```
pcl_viewer mesh.vtk
```

Code Explanation#

As the example PCD has only XYZ coordinates, load it into a PointCloud<PointXYZ>.

// 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::io::loadPCDFile ("bun0.pcd", cloud_blob);
pcl::fromPCLPointCloud2 (cloud_blob, *cloud);

The method requires normals, so normals are estimated using the standard method from PCL.

// 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);

Since coordinates and normals are required to be in the same PointCloud, a PointNormal type point cloud is created.

// 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);

The lines below deals with the initialization of the required objects. In the OpenMP version, include <pcl/oneapi/surface_omp/gp3.h> and declare the class by appending OMP to the original class.”

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

// Initialize objects
pcl::GreedyProjectionTriangulationOMP<pcl::PointNormal> gp3;

The lines below set the parameters, as explained above, to configure the algorithm.

// Set the maximum distance between connected points (maximum edge length)
gp3.setSearchRadius (0.025);

// 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);
gp3.setNumberOfThreads(3);
gp3.setBlockOverlapPercentage(0.04);
gp3.setRemoveDuplicateMesh(true);

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

gp3.setInputCloud (cloud_with_normals);
gp3.setSearchMethod (tree2);
gp3.reconstruct (triangles);