1.. _pfh_estimation: 2 3Point Feature Histograms (PFH) descriptors 4------------------------------------------ 5 6As point feature representations go, surface normals and curvature estimates 7are somewhat basic in their representations of the geometry around a specific 8point. Though extremely fast and easy to compute, they cannot capture too much 9detail, as they approximate the geometry of a point's k-neighborhood with only 10a few values. As a direct consequence, most scenes will contain many points 11with the same or very similar feature values, thus reducing their informative 12characteristics. 13 14This tutorial introduces a family of 3D feature descriptors coined PFH (Point 15Feature Histograms) for simplicity, presents their theoretical advantages and 16discusses implementation details from PCL's perspective. As a prerequisite, 17please go ahead and read the :ref:`normal_estimation` tutorial first, as PFH 18signatures rely on both xyz 3D data as well as surface normals. 19 20Theoretical primer 21------------------ 22 23The goal of the PFH formulation is to encode a point's k-neighborhood 24geometrical properties by generalizing the mean curvature around the point 25using a multi-dimensional histogram of values. This highly dimensional 26hyperspace provides an informative signature for the feature representation, is 27invariant to the 6D pose of the underlying surface, and copes very well with 28different sampling densities or noise levels present in the neighborhood. 29 30A Point Feature Histogram representation is based on the relationships between 31the points in the k-neighborhood and their estimated surface normals. Simply 32put, it attempts to capture as best as possible the sampled surface variations 33by taking into account all the interactions between the directions of the 34estimated normals. **The resultant hyperspace is thus dependent on the quality of 35the surface normal estimations at each point.** 36 37The figure below presents an influence region diagram of the PFH computation 38for a query point (:math:`p_q`), marked with red and placed in the middle of a 39circle (sphere in 3D) with radius **r**, and all its **k** neighbors (points 40with distances smaller than the radius **r**) are fully interconnected in a 41mesh. The final PFH descriptor is computed as a histogram of relationships 42between all pairs of points in the neighborhood, and thus has a computational 43complexity of :math:`O(k^2)`. 44 45.. image:: images/pfh_estimation/pfh_diagram.png 46 :align: center 47 48 49To compute the relative difference between two points :math:`p_i` and 50:math:`p_j` and their associated normals :math:`n_i` and :math:`n_j`, we 51define a fixed coordinate frame at one of the points (see the figure below). 52 53.. math:: 54 55 {\mathsf u} =& \boldsymbol{n}_s \\ 56 {\mathsf v} =& {\mathsf u} \times \frac{(\boldsymbol{p}_t-\boldsymbol{p}_s)}{{\|\boldsymbol{p}_t-\boldsymbol{p}_s\|}_{2}} \\ 57 {\mathsf w} =& {\mathsf u} \times {\mathsf v} 58 59.. image:: images/pfh_estimation/pfh_frame.png 60 :align: center 61 62Using the above **uvw** frame, the difference between the two normals 63:math:`n_s` and :math:`n_t` can be expressed as a set of angular features as 64follows: 65 66.. math:: 67 68 \alpha &= {\mathsf v} \cdot \boldsymbol{n}_t \\ 69 \phi &= {\mathsf u} \cdot \frac{(\boldsymbol{p}_t - \boldsymbol{p}_s)}{d}\\ 70 \theta &= \arctan ({\mathsf w} \cdot \boldsymbol{n}_t, {\mathsf u} \cdot \boldsymbol{n}_t) \\ 71 72where **d** is the Euclidean distance between the two points 73:math:`\boldsymbol{p}_s` and :math:`\boldsymbol{p}_t`, 74:math:`d={\|\boldsymbol{p}_t-\boldsymbol{p}_s\|}_2`. The quadruplet 75:math:`\langle\alpha, \phi, \theta, d\rangle` is computed for each pair of two 76points in k-neighborhood, therefore reducing the 12 values (xyz and normal 77information) of the two points and their normals to 4. 78 79To estimate a PFH quadruplet for a pair of points, use: 80 81.. code-block:: cpp 82 83 computePairFeatures (const Eigen::Vector4f &p1, const Eigen::Vector4f &n1, 84 const Eigen::Vector4f &p2, const Eigen::Vector4f &n2, 85 float &f1, float &f2, float &f3, float &f4); 86 87See the API documentation for additional details. 88 89To create the final PFH representation for the query point, the set of all 90quadruplets is binned into a histogram. The binning process divides each 91feature’s value range into **b** subdivisions, and counts the number of 92occurrences in each subinterval. Since three out of the four features presented 93above are measures of the angles between normals, their values can easily be 94normalized to the same interval on the trigonometric circle. A binning example 95is to divide each feature interval into the same number of equal parts, and 96therefore create a histogram with :math:`b^4` bins in a fully correlated space. 97In this space, a histogram bin increment corresponds to a point having certain 98values for all its 4 features. The figure below presents examples of Point 99Feature Histograms representations for different points in a cloud. 100 101In some cases, the fourth feature, **d**, does not present an extreme 102significance for 2.5D datasets, usually acquired in robotics, as the distance 103between neighboring points increases from the viewpoint. Therefore, omitting 104**d** for scans where the local point density influences this feature dimension 105has proved to be beneficial. 106 107.. image:: images/pfh_estimation/example_pfhs.jpg 108 :align: center 109 110.. note:: 111 112 For more information and mathematical derivations, including an analysis of PFH signatures for different surface geometries please see [RusuDissertation]_. 113 114 115Estimating PFH features 116----------------------- 117 118Point Feature Histograms are implemented in PCL as part of the `pcl_features 119<http://docs.pointclouds.org/trunk/a02944.html>`_ library. 120 121The default PFH implementation uses 5 binning subdivisions (e.g., each of the 122four feature values will use this many bins from its value interval), and does 123not include the distances (as explained above -- although the 124**computePairFeatures** method can be called by the user to obtain the 125distances too, if desired) which results in a 125-byte array (:math:`5^3`) of 126float values. These are stored in a **pcl::PFHSignature125** point type. 127 128The following code snippet will estimate a set of PFH features for all the 129points in the input dataset. 130 131.. code-block:: cpp 132 :linenos: 133 134 #include <pcl/point_types.h> 135 #include <pcl/features/pfh.h> 136 137 { 138 pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>); 139 pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal> ()); 140 141 ... read, pass in or create a point cloud with normals ... 142 ... (note: you can create a single PointCloud<PointNormal> if you want) ... 143 144 // Create the PFH estimation class, and pass the input dataset+normals to it 145 pcl::PFHEstimation<pcl::PointXYZ, pcl::Normal, pcl::PFHSignature125> pfh; 146 pfh.setInputCloud (cloud); 147 pfh.setInputNormals (normals); 148 // alternatively, if cloud is of tpe PointNormal, do pfh.setInputNormals (cloud); 149 150 // Create an empty kdtree representation, and pass it to the PFH estimation object. 151 // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given). 152 pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ()); 153 //pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr tree (new pcl::KdTreeFLANN<pcl::PointXYZ> ()); -- older call for PCL 1.5- 154 pfh.setSearchMethod (tree); 155 156 // Output datasets 157 pcl::PointCloud<pcl::PFHSignature125>::Ptr pfhs (new pcl::PointCloud<pcl::PFHSignature125> ()); 158 159 // Use all neighbors in a sphere of radius 5cm 160 // IMPORTANT: the radius used here has to be larger than the radius used to estimate the surface normals!!! 161 pfh.setRadiusSearch (0.05); 162 163 // Compute the features 164 pfh.compute (*pfhs); 165 166 // pfhs->size () should have the same size as the input cloud->size ()* 167 } 168 169The actual **compute** call from the **PFHEstimation** class does nothing internally but:: 170 171 for each point p in cloud P 172 173 1. get the nearest neighbors of p 174 175 2. for each pair of neighbors, compute the three angular values 176 177 3. bin all the results in an output histogram 178 179 180To compute a single PFH representation from a k-neighborhood, use: 181 182.. code-block:: cpp 183 184 computePointPFHSignature (const pcl::PointCloud<PointInT> &cloud, 185 const pcl::PointCloud<PointNT> &normals, 186 const std::vector<int> &indices, 187 int nr_split, 188 Eigen::VectorXf &pfh_histogram); 189 190Where *cloud* is the input point cloud that contains the points, *normals* is 191the input point cloud that contains the normals (could be equal to cloud if 192*PointInT=PointNT=PointNormal*), *indices* represents the set of k-nearest 193neighbors from *cloud*, *nr_split* is the number of subdivisions to use for the 194binning process for each feature interval, and *pfh_histogram* is the output 195resultant histogram as an array of float values. 196 197.. note:: 198 199 For efficiency reasons, the **compute** method in **PFHEstimation** does not check if the normals contain NaN or infinite values. 200 Passing such values to **compute()** will result in undefined output. 201 It is advisable to check the normals, at least during the design of the processing chain or when setting the parameters. 202 This can be done by inserting the following code before the call to **compute()**: 203 204 .. code-block:: cpp 205 206 for (int i = 0; i < normals->size(); i++) 207 { 208 if (!pcl::isFinite<pcl::Normal>((*normals)[i])) 209 { 210 PCL_WARN("normals[%d] is not finite\n", i); 211 } 212 } 213 214 In production code, preprocessing steps and parameters should be set so that normals are finite or raise an error. 215 216