Point Cloud Library (PCL)  1.8.1-dev
fpfh.hpp
1 /*
2  * Software License Agreement (BSD License)
3  *
4  * Point Cloud Library (PCL) - www.pointclouds.org
5  * Copyright (c) 2010-2011, Willow Garage, Inc.
6  * Copyright (c) 2012-, Open Perception, Inc.
7  *
8  * All rights reserved.
9  *
10  * Redistribution and use in source and binary forms, with or without
11  * modification, are permitted provided that the following conditions
12  * are met:
13  *
14  * * Redistributions of source code must retain the above copyright
15  * notice, this list of conditions and the following disclaimer.
16  * * Redistributions in binary form must reproduce the above
17  * copyright notice, this list of conditions and the following
18  * disclaimer in the documentation and/or other materials provided
19  * with the distribution.
20  * * Neither the name of the copyright holder(s) nor the names of its
21  * contributors may be used to endorse or promote products derived
22  * from this software without specific prior written permission.
23  *
24  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
27  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
28  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
29  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
30  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
31  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32  * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
33  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
34  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
35  * POSSIBILITY OF SUCH DAMAGE.
36  *
37  * $Id$
38  *
39  */
40 
41 #ifndef PCL_FEATURES_IMPL_FPFH_H_
42 #define PCL_FEATURES_IMPL_FPFH_H_
43 
44 #include <pcl/features/fpfh.h>
45 #include <pcl/features/pfh_tools.h>
46 
47 //////////////////////////////////////////////////////////////////////////////////////////////
48 template <typename PointInT, typename PointNT, typename PointOutT> bool
50  const pcl::PointCloud<PointInT> &cloud, const pcl::PointCloud<PointNT> &normals,
51  int p_idx, int q_idx, float &f1, float &f2, float &f3, float &f4)
52 {
53  pcl::computePairFeatures (cloud.points[p_idx].getVector4fMap (), normals.points[p_idx].getNormalVector4fMap (),
54  cloud.points[q_idx].getVector4fMap (), normals.points[q_idx].getNormalVector4fMap (),
55  f1, f2, f3, f4);
56  return (true);
57 }
58 
59 //////////////////////////////////////////////////////////////////////////////////////////////
60 template <typename PointInT, typename PointNT, typename PointOutT> void
62  const pcl::PointCloud<PointInT> &cloud, const pcl::PointCloud<PointNT> &normals,
63  int p_idx, int row, const std::vector<int> &indices,
64  Eigen::MatrixXf &hist_f1, Eigen::MatrixXf &hist_f2, Eigen::MatrixXf &hist_f3)
65 {
66  Eigen::Vector4f pfh_tuple;
67  // Get the number of bins from the histograms size
68  int nr_bins_f1 = static_cast<int> (hist_f1.cols ());
69  int nr_bins_f2 = static_cast<int> (hist_f2.cols ());
70  int nr_bins_f3 = static_cast<int> (hist_f3.cols ());
71 
72  // Factorization constant
73  float hist_incr = 100.0f / static_cast<float>(indices.size () - 1);
74 
75  // Iterate over all the points in the neighborhood
76  for (size_t idx = 0; idx < indices.size (); ++idx)
77  {
78  // Avoid unnecessary returns
79  if (p_idx == indices[idx])
80  continue;
81 
82  // Compute the pair P to NNi
83  if (!computePairFeatures (cloud, normals, p_idx, indices[idx], pfh_tuple[0], pfh_tuple[1], pfh_tuple[2], pfh_tuple[3]))
84  continue;
85 
86  // Normalize the f1, f2, f3 features and push them in the histogram
87  int h_index = static_cast<int> (floor (nr_bins_f1 * ((pfh_tuple[0] + M_PI) * d_pi_)));
88  if (h_index < 0) h_index = 0;
89  if (h_index >= nr_bins_f1) h_index = nr_bins_f1 - 1;
90  hist_f1 (row, h_index) += hist_incr;
91 
92  h_index = static_cast<int> (floor (nr_bins_f2 * ((pfh_tuple[1] + 1.0) * 0.5)));
93  if (h_index < 0) h_index = 0;
94  if (h_index >= nr_bins_f2) h_index = nr_bins_f2 - 1;
95  hist_f2 (row, h_index) += hist_incr;
96 
97  h_index = static_cast<int> (floor (nr_bins_f3 * ((pfh_tuple[2] + 1.0) * 0.5)));
98  if (h_index < 0) h_index = 0;
99  if (h_index >= nr_bins_f3) h_index = nr_bins_f3 - 1;
100  hist_f3 (row, h_index) += hist_incr;
101  }
102 }
103 
104 //////////////////////////////////////////////////////////////////////////////////////////////
105 template <typename PointInT, typename PointNT, typename PointOutT> void
107  const Eigen::MatrixXf &hist_f1, const Eigen::MatrixXf &hist_f2, const Eigen::MatrixXf &hist_f3,
108  const std::vector<int> &indices, const std::vector<float> &dists, Eigen::VectorXf &fpfh_histogram)
109 {
110  assert (indices.size () == dists.size ());
111  double sum_f1 = 0.0, sum_f2 = 0.0, sum_f3 = 0.0;
112  float weight = 0.0, val_f1, val_f2, val_f3;
113 
114  // Get the number of bins from the histograms size
115  int nr_bins_f1 = static_cast<int> (hist_f1.cols ());
116  int nr_bins_f2 = static_cast<int> (hist_f2.cols ());
117  int nr_bins_f3 = static_cast<int> (hist_f3.cols ());
118  int nr_bins_f12 = nr_bins_f1 + nr_bins_f2;
119 
120  // Clear the histogram
121  fpfh_histogram.setZero (nr_bins_f1 + nr_bins_f2 + nr_bins_f3);
122 
123  // Use the entire patch
124  for (size_t idx = 0, data_size = indices.size (); idx < data_size; ++idx)
125  {
126  // Minus the query point itself
127  if (dists[idx] == 0)
128  continue;
129 
130  // Standard weighting function used
131  weight = 1.0f / dists[idx];
132 
133  // Weight the SPFH of the query point with the SPFH of its neighbors
134  for (int f1_i = 0; f1_i < nr_bins_f1; ++f1_i)
135  {
136  val_f1 = hist_f1 (indices[idx], f1_i) * weight;
137  sum_f1 += val_f1;
138  fpfh_histogram[f1_i] += val_f1;
139  }
140 
141  for (int f2_i = 0; f2_i < nr_bins_f2; ++f2_i)
142  {
143  val_f2 = hist_f2 (indices[idx], f2_i) * weight;
144  sum_f2 += val_f2;
145  fpfh_histogram[f2_i + nr_bins_f1] += val_f2;
146  }
147 
148  for (int f3_i = 0; f3_i < nr_bins_f3; ++f3_i)
149  {
150  val_f3 = hist_f3 (indices[idx], f3_i) * weight;
151  sum_f3 += val_f3;
152  fpfh_histogram[f3_i + nr_bins_f12] += val_f3;
153  }
154  }
155 
156  if (sum_f1 != 0)
157  sum_f1 = 100.0 / sum_f1; // histogram values sum up to 100
158  if (sum_f2 != 0)
159  sum_f2 = 100.0 / sum_f2; // histogram values sum up to 100
160  if (sum_f3 != 0)
161  sum_f3 = 100.0 / sum_f3; // histogram values sum up to 100
162 
163  // Adjust final FPFH values
164  for (int f1_i = 0; f1_i < nr_bins_f1; ++f1_i)
165  fpfh_histogram[f1_i] *= static_cast<float> (sum_f1);
166  for (int f2_i = 0; f2_i < nr_bins_f2; ++f2_i)
167  fpfh_histogram[f2_i + nr_bins_f1] *= static_cast<float> (sum_f2);
168  for (int f3_i = 0; f3_i < nr_bins_f3; ++f3_i)
169  fpfh_histogram[f3_i + nr_bins_f12] *= static_cast<float> (sum_f3);
170 }
171 
172 //////////////////////////////////////////////////////////////////////////////////////////////
173 template <typename PointInT, typename PointNT, typename PointOutT> void
175  Eigen::MatrixXf &hist_f1, Eigen::MatrixXf &hist_f2, Eigen::MatrixXf &hist_f3)
176 {
177  // Allocate enough space to hold the NN search results
178  // \note This resize is irrelevant for a radiusSearch ().
179  std::vector<int> nn_indices (k_);
180  std::vector<float> nn_dists (k_);
181 
182  std::set<int> spfh_indices;
183  spfh_hist_lookup.resize (surface_->points.size ());
184 
185  // Build a list of (unique) indices for which we will need to compute SPFH signatures
186  // (We need an SPFH signature for every point that is a neighbor of any point in input_[indices_])
187  if (surface_ != input_ ||
188  indices_->size () != surface_->points.size ())
189  {
190  for (size_t idx = 0; idx < indices_->size (); ++idx)
191  {
192  int p_idx = (*indices_)[idx];
193  if (this->searchForNeighbors (p_idx, search_parameter_, nn_indices, nn_dists) == 0)
194  continue;
195 
196  spfh_indices.insert (nn_indices.begin (), nn_indices.end ());
197  }
198  }
199  else
200  {
201  // Special case: When a feature must be computed at every point, there is no need for a neighborhood search
202  for (size_t idx = 0; idx < indices_->size (); ++idx)
203  spfh_indices.insert (static_cast<int> (idx));
204  }
205 
206  // Initialize the arrays that will store the SPFH signatures
207  size_t data_size = spfh_indices.size ();
208  hist_f1.setZero (data_size, nr_bins_f1_);
209  hist_f2.setZero (data_size, nr_bins_f2_);
210  hist_f3.setZero (data_size, nr_bins_f3_);
211 
212  // Compute SPFH signatures for every point that needs them
213  std::set<int>::iterator spfh_indices_itr = spfh_indices.begin ();
214  for (int i = 0; i < static_cast<int> (spfh_indices.size ()); ++i)
215  {
216  // Get the next point index
217  int p_idx = *spfh_indices_itr;
218  ++spfh_indices_itr;
219 
220  // Find the neighborhood around p_idx
221  if (this->searchForNeighbors (*surface_, p_idx, search_parameter_, nn_indices, nn_dists) == 0)
222  continue;
223 
224  // Estimate the SPFH signature around p_idx
225  computePointSPFHSignature (*surface_, *normals_, p_idx, i, nn_indices, hist_f1, hist_f2, hist_f3);
226 
227  // Populate a lookup table for converting a point index to its corresponding row in the spfh_hist_* matrices
228  spfh_hist_lookup[p_idx] = i;
229  }
230 }
231 
232 //////////////////////////////////////////////////////////////////////////////////////////////
233 template <typename PointInT, typename PointNT, typename PointOutT> void
235 {
236  // Allocate enough space to hold the NN search results
237  // \note This resize is irrelevant for a radiusSearch ().
238  std::vector<int> nn_indices (k_);
239  std::vector<float> nn_dists (k_);
240 
241  std::vector<int> spfh_hist_lookup;
242  computeSPFHSignatures (spfh_hist_lookup, hist_f1_, hist_f2_, hist_f3_);
243 
244  output.is_dense = true;
245  // Save a few cycles by not checking every point for NaN/Inf values if the cloud is set to dense
246  if (input_->is_dense)
247  {
248  // Iterate over the entire index vector
249  for (size_t idx = 0; idx < indices_->size (); ++idx)
250  {
251  if (this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0)
252  {
253  for (int d = 0; d < fpfh_histogram_.size (); ++d)
254  output.points[idx].histogram[d] = std::numeric_limits<float>::quiet_NaN ();
255 
256  output.is_dense = false;
257  continue;
258  }
259 
260  // ... and remap the nn_indices values so that they represent row indices in the spfh_hist_* matrices
261  // instead of indices into surface_->points
262  for (size_t i = 0; i < nn_indices.size (); ++i)
263  nn_indices[i] = spfh_hist_lookup[nn_indices[i]];
264 
265  // Compute the FPFH signature (i.e. compute a weighted combination of local SPFH signatures) ...
266  weightPointSPFHSignature (hist_f1_, hist_f2_, hist_f3_, nn_indices, nn_dists, fpfh_histogram_);
267 
268  // ...and copy it into the output cloud
269  for (int d = 0; d < fpfh_histogram_.size (); ++d)
270  output.points[idx].histogram[d] = fpfh_histogram_[d];
271  }
272  }
273  else
274  {
275  // Iterate over the entire index vector
276  for (size_t idx = 0; idx < indices_->size (); ++idx)
277  {
278  if (!isFinite ((*input_)[(*indices_)[idx]]) ||
279  this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0)
280  {
281  for (int d = 0; d < fpfh_histogram_.size (); ++d)
282  output.points[idx].histogram[d] = std::numeric_limits<float>::quiet_NaN ();
283 
284  output.is_dense = false;
285  continue;
286  }
287 
288  // ... and remap the nn_indices values so that they represent row indices in the spfh_hist_* matrices
289  // instead of indices into surface_->points
290  for (size_t i = 0; i < nn_indices.size (); ++i)
291  nn_indices[i] = spfh_hist_lookup[nn_indices[i]];
292 
293  // Compute the FPFH signature (i.e. compute a weighted combination of local SPFH signatures) ...
294  weightPointSPFHSignature (hist_f1_, hist_f2_, hist_f3_, nn_indices, nn_dists, fpfh_histogram_);
295 
296  // ...and copy it into the output cloud
297  for (int d = 0; d < fpfh_histogram_.size (); ++d)
298  output.points[idx].histogram[d] = fpfh_histogram_[d];
299  }
300  }
301 }
302 
303 #define PCL_INSTANTIATE_FPFHEstimation(T,NT,OutT) template class PCL_EXPORTS pcl::FPFHEstimation<T,NT,OutT>;
304 
305 #endif // PCL_FEATURES_IMPL_FPFH_H_
306 
bool isFinite(const PointT &pt)
Tests if the 3D components of a point are all finite param[in] pt point to be tested return true if f...
Definition: point_tests.h:54
void computeSPFHSignatures(std::vector< int > &spf_hist_lookup, Eigen::MatrixXf &hist_f1, Eigen::MatrixXf &hist_f2, Eigen::MatrixXf &hist_f3)
Estimate the set of all SPFH (Simple Point Feature Histograms) signatures for the input cloud...
Definition: fpfh.hpp:174
void computeFeature(PointCloudOut &output)
Estimate the Fast Point Feature Histograms (FPFH) descriptors at a set of points given by <setInputCl...
Definition: fpfh.hpp:234
std::vector< PointT, Eigen::aligned_allocator< PointT > > points
The point data.
Definition: point_cloud.h:410
void weightPointSPFHSignature(const Eigen::MatrixXf &hist_f1, const Eigen::MatrixXf &hist_f2, const Eigen::MatrixXf &hist_f3, const std::vector< int > &indices, const std::vector< float > &dists, Eigen::VectorXf &fpfh_histogram)
Weight the SPFH (Simple Point Feature Histograms) individual histograms to create the final FPFH (Fas...
Definition: fpfh.hpp:106
bool computePairFeatures(const pcl::PointCloud< PointInT > &cloud, const pcl::PointCloud< PointNT > &normals, int p_idx, int q_idx, float &f1, float &f2, float &f3, float &f4)
Compute the 4-tuple representation containing the three angles and one distance between two points re...
Definition: fpfh.hpp:49
void computePointSPFHSignature(const pcl::PointCloud< PointInT > &cloud, const pcl::PointCloud< PointNT > &normals, int p_idx, int row, const std::vector< int > &indices, Eigen::MatrixXf &hist_f1, Eigen::MatrixXf &hist_f2, Eigen::MatrixXf &hist_f3)
Estimate the SPFH (Simple Point Feature Histograms) individual signatures of the three angular (f1...
Definition: fpfh.hpp:61
PCL_EXPORTS bool computePairFeatures(const Eigen::Vector4f &p1, const Eigen::Vector4f &n1, const Eigen::Vector4f &p2, const Eigen::Vector4f &n2, float &f1, float &f2, float &f3, float &f4)
Compute the 4-tuple representation containing the three angles and one distance between two points re...
bool is_dense
True if no points are invalid (e.g., have NaN or Inf values).
Definition: point_cloud.h:418