Surface normals are important properties of a geometric surface, and are
heavily used in many areas such as computer graphics applications, to apply the
correct light sources that generate shadings and other visual effects.
Given a geometric surface, it’s usually trivial to infer the direction of the
normal at a certain point on the surface as the vector perpendicular to the
surface at that point. However, since the point cloud datasets that we acquire
represent a set of point samples on the real surface, there are two
possibilities:
obtain the underlying surface from the acquired point cloud dataset, using
surface meshing techniques, and then compute the surface normals from the
mesh;
use approximations to infer the surface normals from the point cloud dataset
directly.
This tutorial will address the latter, that is, given a point cloud dataset,
directly compute the surface normals at each point in the cloud.
Though many different normal estimation methods exist, the one that we will
concentrate on this tutorial is one of the simplest, and is formulated as
follows. The problem of determining the normal to a point on the surface is
approximated by the problem of estimating the normal of a plane tangent to the
surface, which in turn becomes a least-square plane fitting estimation problem.
Note
For more information, including the mathematical equations of the least-squares problem, see [RusuDissertation].
The solution for estimating the surface normal is therefore reduced to an
analysis of the eigenvectors and eigenvalues (or PCA – Principal Component
Analysis) of a covariance matrix created from the nearest neighbors of the
query point. More specifically, for each point
, we
assemble the covariance matrix
as follows:
Where
is the number of point neighbors considered in the
neighborhood of
,
represents the 3D centroid of the nearest neighbors,
is the
-th eigenvalue of the covariance matrix, and
the
-th eigenvector.
To estimate a covariance matrix from a set of points in PCL, you can use:
1// Placeholder for the 3x3 covariance matrix at each surface patch 2Eigen::Matrix3fcovariance_matrix; 3// 16-bytes aligned placeholder for the XYZ centroid of a surface patch 4Eigen::Vector4fxyz_centroid; 5 6// Estimate the XYZ centroid 7compute3DCentroid(cloud,xyz_centroid); 8 9// Compute the 3x3 covariance matrix10computeCovarianceMatrix(cloud,xyz_centroid,covariance_matrix);
In general, because there is no mathematical way to solve for the sign of the
normal, its orientation computed via Principal Component Analysis (PCA) as
shown above is ambiguous, and not consistently oriented over an entire point
cloud dataset. The figure below presents these effects on two sections of a
larger dataset representing a part of a kitchen environment. The right part of
the figure presents the Extended Gaussian Image (EGI), also known as the normal
sphere, which describes the orientation of all normals from the point cloud.
Since the datasets are 2.5D and have thus been acquired from a single
viewpoint, normals should be present only on half of the sphere in the EGI.
However, due to the orientation inconsistency, they are spread across the
entire sphere.
The solution to this problem is trivial if the viewpoint
is in fact known. To orient all normals
consistently towards the viewpoint, they need to satisfy the equation:
The figure below presents the results after all normals in the datasets from
the above figure have been consistently oriented towards the viewpoint.
To re-orient a given point normal manually in PCL, you can use:
If the dataset has multiple acquisition viewpoints, then the above normal re-orientation method does not hold, and more complex algorithms need to be implemented. Please see [RusuDissertation] for more information.
As previously explained, a surface normal at a point needs to be estimated from
the surrounding point neighborhood support of the point (also called
k-neighborhood).
The specifics of the nearest-neighbor estimation problem raise the question of
the right scale factor: given a sampled point cloud dataset, what are the
correct k (given via pcl::Feature::setKSearch) or r (given via
pcl::Feature::setRadiusSearch) values that should be used in determining
the set of nearest neighbors of a point?
This issue is of extreme importance and constitutes a limiting factor in the
automatic estimation (i.e., without user given thresholds) of a point feature
representation. To better illustrate this issue, the figure below presents the
effects of selecting a smaller scale (i.e., small r or k) versus a
larger scale (i.e., large r or k). The left part of the figures depicts
a reasonable well-chosen scale factor, with estimated surface normals
approximately perpendicular for the two planar surfaces and small edges
visible all across the table. If the scale factor however is too big (right
part), and thus the set of neighbors is larger covering points from adjacent
surfaces, the estimated point feature representations get distorted, with
rotated surface normals at the edges of the two planar surfaces, and smeared
edges and suppressed fine details.
Without going into too many details, it suffices to assume that for now, the
scale for the determination of a point’s neighborhood has to be selected based
on the level of detail required by the application. Simply put, if the
curvature at the edge between the handle of a mug and the cylindrical part is
important, the scale factor needs to be small enough to capture those details,
and large otherwise.
Though an example of normal estimation has already been given in
Features, we will revise one of them here for the purpose of
better explaining what goes on behind the scenes.
The following code snippet will estimate a set of surface normals for all the
points in the input dataset.
1#include<pcl/point_types.h> 2#include<pcl/features/normal_3d.h> 3 4{ 5pcl::PointCloud<pcl::PointXYZ>::Ptrcloud(newpcl::PointCloud<pcl::PointXYZ>); 6 7...read,passinorcreateapointcloud... 8 9// Create the normal estimation class, and pass the input dataset to it10pcl::NormalEstimation<pcl::PointXYZ,pcl::Normal>ne;11ne.setInputCloud(cloud);1213// Create an empty kdtree representation, and pass it to the normal estimation object.14// Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).15pcl::search::KdTree<pcl::PointXYZ>::Ptrtree(newpcl::search::KdTree<pcl::PointXYZ>());16ne.setSearchMethod(tree);1718// Output datasets19pcl::PointCloud<pcl::Normal>::Ptrcloud_normals(newpcl::PointCloud<pcl::Normal>);2021// Use all neighbors in a sphere of radius 3cm22ne.setRadiusSearch(0.03);2324// Compute the features25ne.compute(*cloud_normals);2627// cloud_normals->size () should have the same size as the input cloud->size ()*28}
The actual compute call from the NormalEstimation class does nothing internally but:
Where cloud is the input point cloud that contains the points, indices
represents the set of k-nearest neighbors from cloud, and plane_parameters
and curvature represent the output of the normal estimation, with
plane_parameters holding the normal (nx, ny, nz) on the first 3 coordinates,
and the fourth coordinate is D = nc . p_plane (centroid here) + p. The output surface curvature is estimated as a relationship between the eigenvalues of the covariance matrix (as presented above), as:
For a detailed description of this property, see M. Pauly, M. Gross and L. P. Kobbelt, “Efficient simplification of point-sampled surfaces”.
For the speed-savvy users, PCL provides an additional implementation of surface
normal estimation which uses multi-core/multi-threaded paradigms using OpenMP
to speed the computation. The name of the class is
pcl::NormalEstimationOMP, and its API is 100% compatible to the
single-threaded pcl::NormalEstimation, which makes it suitable as a drop-in
replacement. On a system with 8 cores, you should get anything between 6-8
times faster computation times.
Note
If your dataset is organized (e.g., acquired using a TOF camera, stereo camera, etc – that is, it has a width and a height), for even faster results see the Normal Estimation Using Integral Images.