Point Cloud Library (PCL)  1.7.0
susan.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  *
7  * All rights reserved.
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions
11  * are met:
12  *
13  * * Redistributions of source code must retain the above copyright
14  * notice, this list of conditions and the following disclaimer.
15  * * Redistributions in binary form must reproduce the above
16  * copyright notice, this list of conditions and the following
17  * disclaimer in the documentation and/or other materials provided
18  * with the distribution.
19  * * Neither the name of Willow Garage, Inc. nor the names of its
20  * contributors may be used to endorse or promote products derived
21  * from this software without specific prior written permission.
22  *
23  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
29  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31  * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
33  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34  * POSSIBILITY OF SUCH DAMAGE.
35  *
36  */
37 
38 #ifndef PCL_SUSAN_IMPL_HPP_
39 #define PCL_SUSAN_IMPL_HPP_
40 
41 #include <pcl/keypoints/susan.h>
42 #include <pcl/features/normal_3d.h>
43 #include <pcl/features/integral_image_normal.h>
44 
45 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
46 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
48 {
49  nonmax_ = nonmax;
50 }
51 
52 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
53 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
55 {
56  geometric_validation_ = validate;
57 }
58 
59 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
60 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
62 {
63  search_radius_ = radius;
64 }
65 
66 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
67 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
69 {
70  distance_threshold_ = distance_threshold;
71 }
72 
73 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
74 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
76 {
77  angular_threshold_ = angular_threshold;
78 }
79 
80 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
81 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
83 {
84  intensity_threshold_ = intensity_threshold;
85 }
86 
87 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
88 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
90 {
91  normals_ = normals;
92 }
93 
94 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
95 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
97 {
98  surface_ = cloud;
99  normals_.reset (new pcl::PointCloud<NormalT>);
100 }
101 
102 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
103 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
105 {
106  threads_ = nr_threads;
107 }
108 
109 
110 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
111 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
112 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::USAN (const PointInT& nucleus,
113 // const NormalT& nucleus_normal,
114 // const std::vector<int>& neighbors,
115 // const float& t,
116 // float& response,
117 // Eigen::Vector3f& centroid) const
118 // {
119 // float area = 0;
120 // response = 0;
121 // float x0 = nucleus.x;
122 // float y0 = nucleus.y;
123 // float z0 = nucleus.z;
124 // //xx xy xz yy yz zz
125 // std::vector<float> coefficients(6);
126 // memset (&coefficients[0], 0, sizeof (float) * 6);
127 // for (std::vector<int>::const_iterator index = neighbors.begin(); index != neighbors.end(); ++index)
128 // {
129 // if (pcl_isfinite (normals_->points[*index].normal_x))
130 // {
131 // Eigen::Vector3f diff = normals_->points[*index].getNormal3fMap () - nucleus_normal.getNormal3fMap ();
132 // float c = diff.norm () / t;
133 // c = -1 * pow (c, 6.0);
134 // c = exp (c);
135 // Eigen::Vector3f xyz = surface_->points[*index].getVector3fMap ();
136 // centroid += c * xyz;
137 // area += c;
138 // coefficients[0] += c * (x0 - xyz.x) * (x0 - xyz.x);
139 // coefficients[1] += c * (x0 - xyz.x) * (y0 - xyz.y);
140 // coefficients[2] += c * (x0 - xyz.x) * (z0 - xyz.z);
141 // coefficients[3] += c * (y0 - xyz.y) * (y0 - xyz.y);
142 // coefficients[4] += c * (y0 - xyz.y) * (z0 - xyz.z);
143 // coefficients[5] += c * (z0 - xyz.z) * (z0 - xyz.z);
144 // }
145 // }
146 
147 // if (area > 0)
148 // {
149 // centroid /= area;
150 // if (area < geometric_threshold)
151 // response = geometric_threshold - area;
152 // // Look for edge direction
153 // // X direction
154 // if ((coefficients[3]/coefficients[0]) < 1 && (coefficients[5]/coefficients[0]) < 1)
155 // direction = Eigen::Vector3f (1, 0, 0);
156 // else
157 // {
158 // // Y direction
159 // if ((coefficients[0]/coefficients[3]) < 1 && (coefficients[5]/coefficients[3]) < 1)
160 // direction = Eigen::Vector3f (0, 1, 0);
161 // else
162 // {
163 // // Z direction
164 // if ((coefficients[0]/coefficients[5]) < 1 && (coefficients[3]/coefficients[5]) < 1)
165 // direction = Eigen::Vector3f (0, 1, 0);
166 // // Diagonal edge
167 // else
168 // {
169 // //XY direction
170 // if ((coefficients[2]/coeffcients[1]) < 1 && (coeffcients[4]/coeffcients[1]) < 1)
171 // {
172 // if (coeffcients[1] > 0)
173 // direction = Eigen::Vector3f (1,1,0);
174 // else
175 // direction = Eigen::Vector3f (-1,1,0);
176 // }
177 // else
178 // {
179 // //XZ direction
180 // if ((coefficients[1]/coeffcients[2]) > 1 && (coeffcients[4]/coeffcients[2]) < 1)
181 // {
182 // if (coeffcients[2] > 0)
183 // direction = Eigen::Vector3f (1,0,1);
184 // else
185 // direction = Eigen::Vector3f (-1,0,1);
186 // }
187 // //YZ direction
188 // else
189 // {
190 // if (coeffcients[4] > 0)
191 // direction = Eigen::Vector3f (0,1,1);
192 // else
193 // direction = Eigen::Vector3f (0,-1,1);
194 // }
195 // }
196 // }
197 // }
198 // }
199 
200 // // std::size_t max_index = std::distance (coefficients.begin (), max);
201 // // switch (max_index)
202 // // {
203 // // case 0 : direction = Eigen::Vector3f (1, 0, 0); break;
204 // // case 1 : direction = Eigen::Vector3f (1, 1, 0); break;
205 // // case 2 : direction = Eigen::Vector3f (1, 0, 1); break;
206 // // case 3 : direction = Eigen::Vector3f (0, 1, 0); break;
207 // // case 4 : direction = Eigen::Vector3f (0, 1, 1); break;
208 // // case 5 : direction = Eigen::Vector3f (0, 0, 1); break;
209 // // }
210 // }
211 // }
212 
213 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
214 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
216 {
218  {
219  PCL_ERROR ("[pcl::%s::initCompute] init failed!\n", name_.c_str ());
220  return (false);
221  }
222 
223  if (normals_->empty ())
224  {
225  PointCloudNPtr normals (new PointCloudN ());
226  normals->reserve (normals->size ());
227  if (!surface_->isOrganized ())
228  {
230  normal_estimation.setInputCloud (surface_);
231  normal_estimation.setRadiusSearch (search_radius_);
232  normal_estimation.compute (*normals);
233  }
234  else
235  {
238  normal_estimation.setInputCloud (surface_);
239  normal_estimation.setNormalSmoothingSize (5.0);
240  normal_estimation.compute (*normals);
241  }
242  normals_ = normals;
243  }
244  if (normals_->size () != surface_->size ())
245  {
246  PCL_ERROR ("[pcl::%s::initCompute] normals given, but the number of normals does not match the number of input points!\n", name_.c_str ());
247  return (false);
248  }
249  return (true);
250 }
251 
252 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
253 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
255  const Eigen::Vector3f& centroid,
256  const Eigen::Vector3f& nc,
257  const PointInT& point) const
258 {
259  Eigen::Vector3f pc = centroid - point.getVector3fMap ();
260  Eigen::Vector3f pn = nucleus - point.getVector3fMap ();
261  Eigen::Vector3f pc_cross_nc = pc.cross (nc);
262  return ((pc_cross_nc.norm () <= tolerance_) && (pc.dot (nc) >= 0) && (pn.dot (nc) <= 0));
263 }
264 
265 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
266 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::isValidQueryPoint3D (int point_index) const
267 // {
268 // return (isFinite (surface_->points [point_index]) &&
269 // isFinite (normals_->points [point_index]));
270 // }
271 
272 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
273 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::isValidQueryPoint2D (int point_index) const
274 // {
275 // return (isFinite (surface_->points [point_index]));
276 // }
277 
278 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
279 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::isWithinSusan2D (int nucleus, int neighbor) const
280 // {
281 // return (fabs (intensity_ (surface_->points[nucleus]) -
282 // intensity_ (surface_->points[neighbor])) <= intensity_threshold_);
283 // }
284 
285 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
286 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::isWithinSusan3D (int nucleus, int neighbor) const
287 // {
288 // Eigen::Vector3f nucleus_normal = normals_->point[nucleus].getVector3fMap ();
289 // return (1 - nucleus_normal.dot (normals_->points[*index].getNormalVector3fMap ()) <= angular_threshold_);
290 // }
291 
292 // template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> bool
293 // pcl::SUSANKeypoint<PointInT, PointOutT, NormalT, IntensityT>::isWithinSusanH (int nucleus, int neighbor) const
294 // {
295 // Eigen::Vector3f nucleus_normal = normals_->point[nucleus].getVector3fMap ();
296 // return ((1 - nucleus_normal.dot (normals_->points[*index].getNormalVector3fMap ()) <= angular_threshold_) ||
297 // (fabs (intensity_ (surface_->points[nucleus]) - intensity_ (surface_->points[neighbor])) <= intensity_threshold_));
298 // }
299 
300 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
301 template <typename PointInT, typename PointOutT, typename NormalT, typename IntensityT> void
303 {
304  boost::shared_ptr<pcl::PointCloud<PointOutT> > response (new pcl::PointCloud<PointOutT> ());
305  response->reserve (surface_->size ());
306 
307  // Check if the output has a "label" field
308  label_idx_ = pcl::getFieldIndex<PointOutT> (output, "label", out_fields_);
309 
310  const int input_size = static_cast<int> (input_->size ());
311 //#ifdef _OPENMP
312 //#pragma omp parallel for shared (response) num_threads(threads_)
313 //#endif
314  for (int point_index = 0; point_index < input_size; ++point_index)
315  {
316  const PointInT& point_in = input_->points [point_index];
317  const NormalT& normal_in = normals_->points [point_index];
318  if (!isFinite (point_in) || !isFinite (normal_in))
319  continue;
320  else
321  {
322  Eigen::Vector3f nucleus = point_in.getVector3fMap ();
323  Eigen::Vector3f nucleus_normal = normals_->points [point_index].getNormalVector3fMap ();
324  float nucleus_intensity = intensity_ (point_in);
325  std::vector<int> nn_indices;
326  std::vector<float> nn_dists;
327  tree_->radiusSearch (point_in, search_radius_, nn_indices, nn_dists);
328  float area = 0;
329  Eigen::Vector3f centroid = Eigen::Vector3f::Zero ();
330  // Exclude nucleus from the usan
331  std::vector<int> usan; usan.reserve (nn_indices.size () - 1);
332  for (std::vector<int>::const_iterator index = nn_indices.begin(); index != nn_indices.end(); ++index)
333  {
334  if ((*index != point_index) && pcl_isfinite (normals_->points[*index].normal_x))
335  {
336  // if the point fulfill condition
337  if ((fabs (nucleus_intensity - intensity_ (input_->points[*index])) <= intensity_threshold_) ||
338  (1 - nucleus_normal.dot (normals_->points[*index].getNormalVector3fMap ()) <= angular_threshold_))
339  {
340  ++area;
341  centroid += input_->points[*index].getVector3fMap ();
342  usan.push_back (*index);
343  }
344  }
345  }
346 
347  float geometric_threshold = 0.5f * (static_cast<float> (nn_indices.size () - 1));
348  if ((area > 0) && (area < geometric_threshold))
349  {
350  // if no geometric validation required add the point to the response
351  if (!geometric_validation_)
352  {
353  PointOutT point_out;
354  point_out.getVector3fMap () = point_in.getVector3fMap ();
355  //point_out.intensity = geometric_threshold - area;
356  intensity_out_.set (point_out, geometric_threshold - area);
357  // if a label field is found use it to save the index
358  if (label_idx_ != -1)
359  {
360  // save the index in the cloud
361  uint32_t label = static_cast<uint32_t> (point_index);
362  memcpy (reinterpret_cast<char*> (&point_out) + out_fields_[label_idx_].offset,
363  &label, sizeof (uint32_t));
364  }
365 //#ifdef _OPENMP
366 //#pragma omp critical
367 //#endif
368  response->push_back (point_out);
369  }
370  else
371  {
372  centroid /= area;
373  Eigen::Vector3f dist = nucleus - centroid;
374  // Check the distance <= distance_threshold_
375  if (dist.norm () >= distance_threshold_)
376  {
377  // point is valid from distance point of view
378  Eigen::Vector3f nc = centroid - nucleus;
379  // Check the contiguity
380  std::vector<int>::const_iterator usan_it = usan.begin ();
381  for (; usan_it != usan.end (); ++usan_it)
382  {
383  if (!isWithinNucleusCentroid (nucleus, centroid, nc, input_->points[*usan_it]))
384  break;
385  }
386  // All points within usan lies on the segment [nucleus centroid]
387  if (usan_it == usan.end ())
388  {
389  PointOutT point_out;
390  point_out.getVector3fMap () = point_in.getVector3fMap ();
391  // point_out.intensity = geometric_threshold - area;
392  intensity_out_.set (point_out, geometric_threshold - area);
393  // if a label field is found use it to save the index
394  if (label_idx_ != -1)
395  {
396  // save the index in the cloud
397  uint32_t label = static_cast<uint32_t> (point_index);
398  memcpy (reinterpret_cast<char*> (&point_out) + out_fields_[label_idx_].offset,
399  &label, sizeof (uint32_t));
400  }
401 //#ifdef _OPENMP
402 //#pragma omp critical
403 //#endif
404  response->push_back (point_out);
405  }
406  }
407  }
408  }
409  }
410  }
411 
412  response->height = 1;
413  response->width = static_cast<uint32_t> (response->size ());
414 
415  if (!nonmax_)
416  output = *response;
417  else
418  {
419  output.points.clear ();
420  output.points.reserve (response->points.size());
421 
422 //#ifdef _OPENMP
423 //#pragma omp parallel for shared (output) num_threads(threads_)
424 //#endif
425  for (int idx = 0; idx < static_cast<int> (response->points.size ()); ++idx)
426  {
427  const PointOutT& point_in = response->points [idx];
428  const NormalT& normal_in = normals_->points [idx];
429  //const float intensity = response->points[idx].intensity;
430  const float intensity = intensity_out_ (response->points[idx]);
431  if (!isFinite (point_in) || !isFinite (normal_in) || (intensity == 0))
432  continue;
433  std::vector<int> nn_indices;
434  std::vector<float> nn_dists;
435  tree_->radiusSearch (idx, search_radius_, nn_indices, nn_dists);
436  bool is_minima = true;
437  for (std::vector<int>::const_iterator nn_it = nn_indices.begin(); nn_it != nn_indices.end(); ++nn_it)
438  {
439 // if (intensity > response->points[*nn_it].intensity)
440  if (intensity > intensity_out_ (response->points[*nn_it]))
441  {
442  is_minima = false;
443  break;
444  }
445  }
446  if (is_minima)
447 //#ifdef _OPENMP
448 //#pragma omp critical
449 //#endif
450  output.points.push_back (response->points[idx]);
451  }
452 
453  output.height = 1;
454  output.width = static_cast<uint32_t> (output.points.size());
455  }
456  // we don not change the denseness
457  output.is_dense = input_->is_dense;
458 }
459 
460 #define PCL_INSTANTIATE_SUSAN(T,U,N) template class PCL_EXPORTS pcl::SUSANKeypoint<T,U,N>;
461 #endif // #ifndef PCL_HARRIS_KEYPOINT_3D_IMPL_H_
462