Point Cloud Library (PCL)  1.11.1
lum.hpp
1 /*
2  * Software License Agreement (BSD License)
3  *
4  * Point Cloud Library (PCL) - www.pointclouds.org
5  * Copyright (c) 2010-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: lum.hpp 5663 2012-05-02 13:49:39Z florentinus $
38  *
39  */
40 
41 #ifndef PCL_REGISTRATION_IMPL_LUM_HPP_
42 #define PCL_REGISTRATION_IMPL_LUM_HPP_
43 
44 #include <tuple>
45 
46 
47 namespace pcl
48 {
49 
50 namespace registration
51 {
52 
53 template<typename PointT> inline void
55 {
56  slam_graph_ = slam_graph;
57 }
58 
59 
60 template<typename PointT> inline typename LUM<PointT>::SLAMGraphPtr
62 {
63  return (slam_graph_);
64 }
65 
66 
67 template<typename PointT> typename LUM<PointT>::SLAMGraph::vertices_size_type
69 {
70  return (num_vertices (*slam_graph_));
71 }
72 
73 
74 template<typename PointT> void
75 LUM<PointT>::setMaxIterations (int max_iterations)
76 {
77  max_iterations_ = max_iterations;
78 }
79 
80 
81 template<typename PointT> inline int
83 {
84  return (max_iterations_);
85 }
86 
87 
88 template<typename PointT> void
89 LUM<PointT>::setConvergenceThreshold (float convergence_threshold)
90 {
91  convergence_threshold_ = convergence_threshold;
92 }
93 
94 
95 template<typename PointT> inline float
97 {
98  return (convergence_threshold_);
99 }
100 
101 
102 template<typename PointT> typename LUM<PointT>::Vertex
104 {
105  Vertex v = add_vertex (*slam_graph_);
106  (*slam_graph_)[v].cloud_ = cloud;
107  if (v == 0 && pose != Eigen::Vector6f::Zero ())
108  {
109  PCL_WARN("[pcl::registration::LUM::addPointCloud] The pose estimate is ignored for the first cloud in the graph since that will become the reference pose.\n");
110  (*slam_graph_)[v].pose_ = Eigen::Vector6f::Zero ();
111  return (v);
112  }
113  (*slam_graph_)[v].pose_ = pose;
114  return (v);
115 }
116 
117 
118 template<typename PointT> inline void
119 LUM<PointT>::setPointCloud (const Vertex &vertex, const PointCloudPtr &cloud)
120 {
121  if (vertex >= getNumVertices ())
122  {
123  PCL_ERROR("[pcl::registration::LUM::setPointCloud] You are attempting to set a point cloud to a non-existing graph vertex.\n");
124  return;
125  }
126  (*slam_graph_)[vertex].cloud_ = cloud;
127 }
128 
129 
130 template<typename PointT> inline typename LUM<PointT>::PointCloudPtr
131 LUM<PointT>::getPointCloud (const Vertex &vertex) const
132 {
133  if (vertex >= getNumVertices ())
134  {
135  PCL_ERROR("[pcl::registration::LUM::getPointCloud] You are attempting to get a point cloud from a non-existing graph vertex.\n");
136  return (PointCloudPtr ());
137  }
138  return ((*slam_graph_)[vertex].cloud_);
139 }
140 
141 
142 template<typename PointT> inline void
143 LUM<PointT>::setPose (const Vertex &vertex, const Eigen::Vector6f &pose)
144 {
145  if (vertex >= getNumVertices ())
146  {
147  PCL_ERROR("[pcl::registration::LUM::setPose] You are attempting to set a pose estimate to a non-existing graph vertex.\n");
148  return;
149  }
150  if (vertex == 0)
151  {
152  PCL_ERROR("[pcl::registration::LUM::setPose] The pose estimate is ignored for the first cloud in the graph since that will become the reference pose.\n");
153  return;
154  }
155  (*slam_graph_)[vertex].pose_ = pose;
156 }
157 
158 
159 template<typename PointT> inline Eigen::Vector6f
160 LUM<PointT>::getPose (const Vertex &vertex) const
161 {
162  if (vertex >= getNumVertices ())
163  {
164  PCL_ERROR("[pcl::registration::LUM::getPose] You are attempting to get a pose estimate from a non-existing graph vertex.\n");
165  return (Eigen::Vector6f::Zero ());
166  }
167  return ((*slam_graph_)[vertex].pose_);
168 }
169 
170 
171 template<typename PointT> inline Eigen::Affine3f
173 {
174  Eigen::Vector6f pose = getPose (vertex);
175  return (pcl::getTransformation (pose (0), pose (1), pose (2), pose (3), pose (4), pose (5)));
176 }
177 
178 
179 template<typename PointT> void
180 LUM<PointT>::setCorrespondences (const Vertex &source_vertex, const Vertex &target_vertex, const pcl::CorrespondencesPtr &corrs)
181 {
182  if (source_vertex >= getNumVertices () || target_vertex >= getNumVertices () || source_vertex == target_vertex)
183  {
184  PCL_ERROR("[pcl::registration::LUM::setCorrespondences] You are attempting to set a set of correspondences between non-existing or identical graph vertices.\n");
185  return;
186  }
187  Edge e;
188  bool present;
189  std::tie (e, present) = edge (source_vertex, target_vertex, *slam_graph_);
190  if (!present)
191  std::tie (e, present) = add_edge (source_vertex, target_vertex, *slam_graph_);
192  (*slam_graph_)[e].corrs_ = corrs;
193 }
194 
195 
196 template<typename PointT> inline pcl::CorrespondencesPtr
197 LUM<PointT>::getCorrespondences (const Vertex &source_vertex, const Vertex &target_vertex) const
198 {
199  if (source_vertex >= getNumVertices () || target_vertex >= getNumVertices ())
200  {
201  PCL_ERROR("[pcl::registration::LUM::getCorrespondences] You are attempting to get a set of correspondences between non-existing graph vertices.\n");
202  return (pcl::CorrespondencesPtr ());
203  }
204  Edge e;
205  bool present;
206  std::tie (e, present) = edge (source_vertex, target_vertex, *slam_graph_);
207  if (!present)
208  {
209  PCL_ERROR("[pcl::registration::LUM::getCorrespondences] You are attempting to get a set of correspondences from a non-existing graph edge.\n");
210  return (pcl::CorrespondencesPtr ());
211  }
212  return ((*slam_graph_)[e].corrs_);
213 }
214 
215 
216 template<typename PointT> void
218 {
219  int n = static_cast<int> (getNumVertices ());
220  if (n < 2)
221  {
222  PCL_ERROR("[pcl::registration::LUM::compute] The slam graph needs at least 2 vertices.\n");
223  return;
224  }
225  for (int i = 0; i < max_iterations_; ++i)
226  {
227  // Linearized computation of C^-1 and C^-1*D and convergence checking for all edges in the graph (results stored in slam_graph_)
228  typename SLAMGraph::edge_iterator e, e_end;
229  for (std::tie (e, e_end) = edges (*slam_graph_); e != e_end; ++e)
230  computeEdge (*e);
231 
232  // Declare matrices G and B
233  Eigen::MatrixXf G = Eigen::MatrixXf::Zero (6 * (n - 1), 6 * (n - 1));
234  Eigen::VectorXf B = Eigen::VectorXf::Zero (6 * (n - 1));
235 
236  // Start at 1 because 0 is the reference pose
237  for (int vi = 1; vi != n; ++vi)
238  {
239  for (int vj = 0; vj != n; ++vj)
240  {
241  // Attempt to use the forward edge, otherwise use backward edge, otherwise there was no edge
242  Edge e;
243  bool present1;
244  std::tie (e, present1) = edge (vi, vj, *slam_graph_);
245  if (!present1)
246  {
247  bool present2;
248  std::tie (e, present2) = edge (vj, vi, *slam_graph_);
249  if (!present2)
250  continue;
251  }
252 
253  // Fill in elements of G and B
254  if (vj > 0)
255  G.block (6 * (vi - 1), 6 * (vj - 1), 6, 6) = -(*slam_graph_)[e].cinv_;
256  G.block (6 * (vi - 1), 6 * (vi - 1), 6, 6) += (*slam_graph_)[e].cinv_;
257  B.segment (6 * (vi - 1), 6) += (present1 ? 1 : -1) * (*slam_graph_)[e].cinvd_;
258  }
259  }
260 
261  // Computation of the linear equation system: GX = B
262  // TODO investigate accuracy vs. speed tradeoff and find the best solving method for our type of linear equation (sparse)
263  Eigen::VectorXf X = G.colPivHouseholderQr ().solve (B);
264 
265  // Update the poses
266  float sum = 0.0;
267  for (int vi = 1; vi != n; ++vi)
268  {
269  Eigen::Vector6f difference_pose = static_cast<Eigen::Vector6f> (-incidenceCorrection (getPose (vi)).inverse () * X.segment (6 * (vi - 1), 6));
270  sum += difference_pose.norm ();
271  setPose (vi, getPose (vi) + difference_pose);
272  }
273 
274  // Convergence check
275  if (sum <= convergence_threshold_ * static_cast<float> (n - 1))
276  return;
277  }
278 }
279 
280 
281 template<typename PointT> typename LUM<PointT>::PointCloudPtr
283 {
284  PointCloudPtr out (new PointCloud);
285  pcl::transformPointCloud (*getPointCloud (vertex), *out, getTransformation (vertex));
286  return (out);
287 }
288 
289 
290 template<typename PointT> typename LUM<PointT>::PointCloudPtr
292 {
293  PointCloudPtr out (new PointCloud);
294  typename SLAMGraph::vertex_iterator v, v_end;
295  for (std::tie (v, v_end) = vertices (*slam_graph_); v != v_end; ++v)
296  {
297  PointCloud temp;
298  pcl::transformPointCloud (*getPointCloud (*v), temp, getTransformation (*v));
299  *out += temp;
300  }
301  return (out);
302 }
303 
304 
305 template<typename PointT> void
307 {
308  // Get necessary local data from graph
309  PointCloudPtr source_cloud = (*slam_graph_)[source (e, *slam_graph_)].cloud_;
310  PointCloudPtr target_cloud = (*slam_graph_)[target (e, *slam_graph_)].cloud_;
311  Eigen::Vector6f source_pose = (*slam_graph_)[source (e, *slam_graph_)].pose_;
312  Eigen::Vector6f target_pose = (*slam_graph_)[target (e, *slam_graph_)].pose_;
313  pcl::CorrespondencesPtr corrs = (*slam_graph_)[e].corrs_;
314 
315  // Build the average and difference vectors for all correspondences
316  std::vector<Eigen::Vector3f, Eigen::aligned_allocator<Eigen::Vector3f> > corrs_aver (corrs->size ());
317  std::vector<Eigen::Vector3f, Eigen::aligned_allocator<Eigen::Vector3f> > corrs_diff (corrs->size ());
318  int oci = 0; // oci = output correspondence iterator
319  for (int ici = 0; ici != static_cast<int> (corrs->size ()); ++ici) // ici = input correspondence iterator
320  {
321  // Compound the point pair onto the current pose
322  Eigen::Vector3f source_compounded = pcl::getTransformation (source_pose (0), source_pose (1), source_pose (2), source_pose (3), source_pose (4), source_pose (5)) * (*source_cloud)[(*corrs)[ici].index_query].getVector3fMap ();
323  Eigen::Vector3f target_compounded = pcl::getTransformation (target_pose (0), target_pose (1), target_pose (2), target_pose (3), target_pose (4), target_pose (5)) * (*target_cloud)[(*corrs)[ici].index_match].getVector3fMap ();
324 
325  // NaN points can not be passed to the remaining computational pipeline
326  if (!std::isfinite (source_compounded (0)) || !std::isfinite (source_compounded (1)) || !std::isfinite (source_compounded (2)) || !std::isfinite (target_compounded (0)) || !std::isfinite (target_compounded (1)) || !std::isfinite (target_compounded (2)))
327  continue;
328 
329  // Compute the point pair average and difference and store for later use
330  corrs_aver[oci] = 0.5 * (source_compounded + target_compounded);
331  corrs_diff[oci] = source_compounded - target_compounded;
332  oci++;
333  }
334  corrs_aver.resize (oci);
335  corrs_diff.resize (oci);
336 
337  // Need enough valid correspondences to get a good triangulation
338  if (oci < 3)
339  {
340  PCL_WARN("[pcl::registration::LUM::compute] The correspondences between vertex %d and %d do not contain enough valid correspondences to be considered for computation.\n", source (e, *slam_graph_), target (e, *slam_graph_));
341  (*slam_graph_)[e].cinv_ = Eigen::Matrix6f::Zero ();
342  (*slam_graph_)[e].cinvd_ = Eigen::Vector6f::Zero ();
343  return;
344  }
345 
346  // Build the matrices M'M and M'Z
347  Eigen::Matrix6f MM = Eigen::Matrix6f::Zero ();
348  Eigen::Vector6f MZ = Eigen::Vector6f::Zero ();
349  for (int ci = 0; ci != oci; ++ci) // ci = correspondence iterator
350  {
351  // Fast computation of summation elements of M'M
352  MM (0, 4) -= corrs_aver[ci] (1);
353  MM (0, 5) += corrs_aver[ci] (2);
354  MM (1, 3) -= corrs_aver[ci] (2);
355  MM (1, 4) += corrs_aver[ci] (0);
356  MM (2, 3) += corrs_aver[ci] (1);
357  MM (2, 5) -= corrs_aver[ci] (0);
358  MM (3, 4) -= corrs_aver[ci] (0) * corrs_aver[ci] (2);
359  MM (3, 5) -= corrs_aver[ci] (0) * corrs_aver[ci] (1);
360  MM (4, 5) -= corrs_aver[ci] (1) * corrs_aver[ci] (2);
361  MM (3, 3) += corrs_aver[ci] (1) * corrs_aver[ci] (1) + corrs_aver[ci] (2) * corrs_aver[ci] (2);
362  MM (4, 4) += corrs_aver[ci] (0) * corrs_aver[ci] (0) + corrs_aver[ci] (1) * corrs_aver[ci] (1);
363  MM (5, 5) += corrs_aver[ci] (0) * corrs_aver[ci] (0) + corrs_aver[ci] (2) * corrs_aver[ci] (2);
364 
365  // Fast computation of M'Z
366  MZ (0) += corrs_diff[ci] (0);
367  MZ (1) += corrs_diff[ci] (1);
368  MZ (2) += corrs_diff[ci] (2);
369  MZ (3) += corrs_aver[ci] (1) * corrs_diff[ci] (2) - corrs_aver[ci] (2) * corrs_diff[ci] (1);
370  MZ (4) += corrs_aver[ci] (0) * corrs_diff[ci] (1) - corrs_aver[ci] (1) * corrs_diff[ci] (0);
371  MZ (5) += corrs_aver[ci] (2) * corrs_diff[ci] (0) - corrs_aver[ci] (0) * corrs_diff[ci] (2);
372  }
373  // Remaining elements of M'M
374  MM (0, 0) = MM (1, 1) = MM (2, 2) = static_cast<float> (oci);
375  MM (4, 0) = MM (0, 4);
376  MM (5, 0) = MM (0, 5);
377  MM (3, 1) = MM (1, 3);
378  MM (4, 1) = MM (1, 4);
379  MM (3, 2) = MM (2, 3);
380  MM (5, 2) = MM (2, 5);
381  MM (4, 3) = MM (3, 4);
382  MM (5, 3) = MM (3, 5);
383  MM (5, 4) = MM (4, 5);
384 
385  // Compute pose difference estimation
386  Eigen::Vector6f D = static_cast<Eigen::Vector6f> (MM.inverse () * MZ);
387 
388  // Compute s^2
389  float ss = 0.0f;
390  for (int ci = 0; ci != oci; ++ci) // ci = correspondence iterator
391  ss += static_cast<float> (std::pow (corrs_diff[ci] (0) - (D (0) + corrs_aver[ci] (2) * D (5) - corrs_aver[ci] (1) * D (4)), 2.0f)
392  + std::pow (corrs_diff[ci] (1) - (D (1) + corrs_aver[ci] (0) * D (4) - corrs_aver[ci] (2) * D (3)), 2.0f)
393  + std::pow (corrs_diff[ci] (2) - (D (2) + corrs_aver[ci] (1) * D (3) - corrs_aver[ci] (0) * D (5)), 2.0f));
394 
395  // When reaching the limitations of computation due to linearization
396  if (ss < 0.0000000000001 || !std::isfinite (ss))
397  {
398  (*slam_graph_)[e].cinv_ = Eigen::Matrix6f::Zero ();
399  (*slam_graph_)[e].cinvd_ = Eigen::Vector6f::Zero ();
400  return;
401  }
402 
403  // Store the results in the slam graph
404  (*slam_graph_)[e].cinv_ = MM * (1.0f / ss);
405  (*slam_graph_)[e].cinvd_ = MZ * (1.0f / ss);
406 }
407 
408 
409 template<typename PointT> inline Eigen::Matrix6f
411 {
412  Eigen::Matrix6f out = Eigen::Matrix6f::Identity ();
413  float cx = std::cos (pose (3)), sx = sinf (pose (3)), cy = std::cos (pose (4)), sy = sinf (pose (4));
414  out (0, 4) = pose (1) * sx - pose (2) * cx;
415  out (0, 5) = pose (1) * cx * cy + pose (2) * sx * cy;
416  out (1, 3) = pose (2);
417  out (1, 4) = -pose (0) * sx;
418  out (1, 5) = -pose (0) * cx * cy + pose (2) * sy;
419  out (2, 3) = -pose (1);
420  out (2, 4) = pose (0) * cx;
421  out (2, 5) = -pose (0) * sx * cy - pose (1) * sy;
422  out (3, 5) = sy;
423  out (4, 4) = sx;
424  out (4, 5) = cx * cy;
425  out (5, 4) = cx;
426  out (5, 5) = -sx * cy;
427  return (out);
428 }
429 
430 } // namespace registration
431 } // namespace pcl
432 
433 #define PCL_INSTANTIATE_LUM(T) template class PCL_EXPORTS pcl::registration::LUM<T>;
434 
435 #endif // PCL_REGISTRATION_IMPL_LUM_HPP_
436 
PointCloudPtr getConcatenatedCloud() const
Return a concatenated point cloud of all the SLAM graph&#39;s point clouds compounded onto their current ...
Definition: lum.hpp:291
void setCorrespondences(const Vertex &source_vertex, const Vertex &target_vertex, const pcl::CorrespondencesPtr &corrs)
Add/change a set of correspondences for one of the SLAM graph&#39;s edges.
Definition: lum.hpp:180
typename PointCloud::Ptr PointCloudPtr
Definition: lum.h:118
PointCloudPtr getPointCloud(const Vertex &vertex) const
Return a point cloud from one of the SLAM graph&#39;s vertices.
Definition: lum.hpp:131
Eigen::Matrix< float, 6, 6 > Matrix6f
Definition: lum.h:55
void compute()
Perform LUM&#39;s globally consistent scan matching.
Definition: lum.hpp:217
Globally Consistent Scan Matching based on an algorithm by Lu and Milios.
Definition: lum.h:111
typename SLAMGraph::edge_descriptor Edge
Definition: lum.h:138
void setLoopGraph(const SLAMGraphPtr &slam_graph)
Set the internal SLAM graph structure.
Definition: lum.hpp:54
typename SLAMGraph::vertex_descriptor Vertex
Definition: lum.h:137
Eigen::Vector6f getPose(const Vertex &vertex) const
Return a pose estimate from one of the SLAM graph&#39;s vertices.
Definition: lum.hpp:160
Eigen::Matrix6f incidenceCorrection(const Eigen::Vector6f &pose)
Returns a pose corrected 6DoF incidence matrix.
Definition: lum.hpp:410
Eigen::Affine3f getTransformation(const Vertex &vertex) const
Return a pose estimate from one of the SLAM graph&#39;s vertices as an affine transformation matrix...
Definition: lum.hpp:172
Eigen::Matrix< float, 6, 1 > Vector6f
Definition: lum.h:54
void transformPointCloud(const pcl::PointCloud< PointT > &cloud_in, pcl::PointCloud< PointT > &cloud_out, const Eigen::Transform< Scalar, 3, Eigen::Affine > &transform, bool copy_all_fields)
Apply an affine transform defined by an Eigen Transform.
Definition: transforms.hpp:221
Vertex addPointCloud(const PointCloudPtr &cloud, const Eigen::Vector6f &pose=Eigen::Vector6f::Zero())
Add a new point cloud to the SLAM graph.
Definition: lum.hpp:103
void setPointCloud(const Vertex &vertex, const PointCloudPtr &cloud)
Change a point cloud on one of the SLAM graph&#39;s vertices.
Definition: lum.hpp:119
PointCloud represents the base class in PCL for storing collections of 3D points. ...
Definition: distances.h:55
int getMaxIterations() const
Get the maximum number of iterations for the compute() method.
Definition: lum.hpp:82
shared_ptr< SLAMGraph > SLAMGraphPtr
Definition: lum.h:136
SLAMGraphPtr getLoopGraph() const
Get the internal SLAM graph structure.
Definition: lum.hpp:61
void getTransformation(Scalar x, Scalar y, Scalar z, Scalar roll, Scalar pitch, Scalar yaw, Eigen::Transform< Scalar, 3, Eigen::Affine > &t)
Create a transformation from the given translation and Euler angles (XYZ-convention) ...
Definition: eigen.hpp:608
void setMaxIterations(int max_iterations)
Set the maximum number of iterations for the compute() method.
Definition: lum.hpp:75
void computeEdge(const Edge &e)
Linearized computation of C^-1 and C^-1*D (results stored in slam_graph_).
Definition: lum.hpp:306
void setConvergenceThreshold(float convergence_threshold)
Set the convergence threshold for the compute() method.
Definition: lum.hpp:89
float getConvergenceThreshold() const
Get the convergence threshold for the compute() method.
Definition: lum.hpp:96
Definition: norms.h:54
shared_ptr< Correspondences > CorrespondencesPtr
SLAMGraph::vertices_size_type getNumVertices() const
Get the number of vertices in the SLAM graph.
Definition: lum.hpp:68
PointCloudPtr getTransformedCloud(const Vertex &vertex) const
Return a point cloud from one of the SLAM graph&#39;s vertices compounded onto its current pose estimate...
Definition: lum.hpp:282
pcl::CorrespondencesPtr getCorrespondences(const Vertex &source_vertex, const Vertex &target_vertex) const
Return a set of correspondences from one of the SLAM graph&#39;s edges.
Definition: lum.hpp:197
void setPose(const Vertex &vertex, const Eigen::Vector6f &pose)
Change a pose estimate on one of the SLAM graph&#39;s vertices.
Definition: lum.hpp:143