Point Cloud Library (PCL)  1.7.1
shot_lrf.hpp
1 /*
2  * Software License Agreement (BSD License)
3  *
4  * Point Cloud Library (PCL) - www.pointclouds.org
5  * Copyright (c) 2012, 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 #ifndef PCL_FEATURES_IMPL_SHOT_LRF_H_
41 #define PCL_FEATURES_IMPL_SHOT_LRF_H_
42 
43 #include <utility>
44 #include <pcl/features/shot_lrf.h>
45 
46 //////////////////////////////////////////////////////////////////////////////////////////////
47 // Compute a local Reference Frame for a 3D feature; the output is stored in the "rf" matrix
48 template<typename PointInT, typename PointOutT> float
49 pcl::SHOTLocalReferenceFrameEstimation<PointInT, PointOutT>::getLocalRF (const int& current_point_idx, Eigen::Matrix3f &rf)
50 {
51  const Eigen::Vector4f& central_point = (*input_)[current_point_idx].getVector4fMap ();
52  std::vector<int> n_indices;
53  std::vector<float> n_sqr_distances;
54 
55  this->searchForNeighbors (current_point_idx, search_parameter_, n_indices, n_sqr_distances);
56 
57  Eigen::Matrix<double, Eigen::Dynamic, 4> vij (n_indices.size (), 4);
58 
59  Eigen::Matrix3d cov_m = Eigen::Matrix3d::Zero ();
60 
61  double distance = 0.0;
62  double sum = 0.0;
63 
64  int valid_nn_points = 0;
65 
66  for (size_t i_idx = 0; i_idx < n_indices.size (); ++i_idx)
67  {
68  Eigen::Vector4f pt = surface_->points[n_indices[i_idx]].getVector4fMap ();
69  if (pt.head<3> () == central_point.head<3> ())
70  continue;
71 
72  // Difference between current point and origin
73  vij.row (valid_nn_points).matrix () = (pt - central_point).cast<double> ();
74  vij (valid_nn_points, 3) = 0;
75 
76  distance = search_parameter_ - sqrt (n_sqr_distances[i_idx]);
77 
78  // Multiply vij * vij'
79  cov_m += distance * (vij.row (valid_nn_points).head<3> ().transpose () * vij.row (valid_nn_points).head<3> ());
80 
81  sum += distance;
82  valid_nn_points++;
83  }
84 
85  if (valid_nn_points < 5)
86  {
87  //PCL_ERROR ("[pcl::%s::getLocalRF] Warning! Neighborhood has less than 5 vertexes. Aborting Local RF computation of feature point (%lf, %lf, %lf)\n", "SHOTLocalReferenceFrameEstimation", central_point[0], central_point[1], central_point[2]);
88  rf.setConstant (std::numeric_limits<float>::quiet_NaN ());
89 
90  return (std::numeric_limits<float>::max ());
91  }
92 
93  cov_m /= sum;
94 
95  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver (cov_m);
96 
97  const double& e1c = solver.eigenvalues ()[0];
98  const double& e2c = solver.eigenvalues ()[1];
99  const double& e3c = solver.eigenvalues ()[2];
100 
101  if (!pcl_isfinite (e1c) || !pcl_isfinite (e2c) || !pcl_isfinite (e3c))
102  {
103  //PCL_ERROR ("[pcl::%s::getLocalRF] Warning! Eigenvectors are NaN. Aborting Local RF computation of feature point (%lf, %lf, %lf)\n", "SHOTLocalReferenceFrameEstimation", central_point[0], central_point[1], central_point[2]);
104  rf.setConstant (std::numeric_limits<float>::quiet_NaN ());
105 
106  return (std::numeric_limits<float>::max ());
107  }
108 
109  // Disambiguation
110  Eigen::Vector4d v1 = Eigen::Vector4d::Zero ();
111  Eigen::Vector4d v3 = Eigen::Vector4d::Zero ();
112  v1.head<3> ().matrix () = solver.eigenvectors ().col (2);
113  v3.head<3> ().matrix () = solver.eigenvectors ().col (0);
114 
115  int plusNormal = 0, plusTangentDirection1=0;
116  for (int ne = 0; ne < valid_nn_points; ne++)
117  {
118  double dp = vij.row (ne).dot (v1);
119  if (dp >= 0)
120  plusTangentDirection1++;
121 
122  dp = vij.row (ne).dot (v3);
123  if (dp >= 0)
124  plusNormal++;
125  }
126 
127  //TANGENT
128  plusTangentDirection1 = 2*plusTangentDirection1 - valid_nn_points;
129  if (plusTangentDirection1 == 0)
130  {
131  int points = 5; //std::min(valid_nn_points*2/2+1, 11);
132  int medianIndex = valid_nn_points/2;
133 
134  for (int i = -points/2; i <= points/2; i++)
135  if ( vij.row (medianIndex - i).dot (v1) > 0)
136  plusTangentDirection1 ++;
137 
138  if (plusTangentDirection1 < points/2+1)
139  v1 *= - 1;
140  }
141  else if (plusTangentDirection1 < 0)
142  v1 *= - 1;
143 
144  //Normal
145  plusNormal = 2*plusNormal - valid_nn_points;
146  if (plusNormal == 0)
147  {
148  int points = 5; //std::min(valid_nn_points*2/2+1, 11);
149  int medianIndex = valid_nn_points/2;
150 
151  for (int i = -points/2; i <= points/2; i++)
152  if ( vij.row (medianIndex - i).dot (v3) > 0)
153  plusNormal ++;
154 
155  if (plusNormal < points/2+1)
156  v3 *= - 1;
157  } else if (plusNormal < 0)
158  v3 *= - 1;
159 
160  rf.row (0).matrix () = v1.head<3> ().cast<float> ();
161  rf.row (2).matrix () = v3.head<3> ().cast<float> ();
162  rf.row (1).matrix () = rf.row (2).cross (rf.row (0));
163 
164  return (0.0f);
165 }
166 
167 //////////////////////////////////////////////////////////////////////////////////////////////
168 template <typename PointInT, typename PointOutT> void
170 {
171  //check whether used with search radius or search k-neighbors
172  if (this->getKSearch () != 0)
173  {
174  PCL_ERROR(
175  "[pcl::%s::computeFeature] Error! Search method set to k-neighborhood. Call setKSearch(0) and setRadiusSearch( radius ) to use this class.\n",
176  getClassName().c_str ());
177  return;
178  }
179  tree_->setSortedResults (true);
180 
181  for (size_t i = 0; i < indices_->size (); ++i)
182  {
183  // point result
184  Eigen::Matrix3f rf;
185  PointOutT& output_rf = output[i];
186 
187  //output_rf.confidence = getLocalRF ((*indices_)[i], rf);
188  //if (output_rf.confidence == std::numeric_limits<float>::max ())
189  if (getLocalRF ((*indices_)[i], rf) == std::numeric_limits<float>::max ())
190  {
191  output.is_dense = false;
192  }
193 
194  for (int d = 0; d < 3; ++d)
195  {
196  output_rf.x_axis[d] = rf.row (0)[d];
197  output_rf.y_axis[d] = rf.row (1)[d];
198  output_rf.z_axis[d] = rf.row (2)[d];
199  }
200  }
201 }
202 
203 #define PCL_INSTANTIATE_SHOTLocalReferenceFrameEstimation(T,OutT) template class PCL_EXPORTS pcl::SHOTLocalReferenceFrameEstimation<T,OutT>;
204 
205 #endif // PCL_FEATURES_IMPL_SHOT_LRF_H_
206