Skip to content

Commit

Permalink
(Minor) modifications for the registration example in the RAM article:
Browse files Browse the repository at this point in the history
* GICP now allows setting (externally computed) covariance matrices.
  If not set, GICP will compute these as usual.
* OrganizedFastMesh now accepts a complete quadratic error in the form
  of sigma_d = a*d*d + b*d + c.
* OrganizedFastMesh now allows to either use the measured depth (z) or
  the actual distance to the point (e.g., for 3D laser scans)
* from_meshes.h contains two additional inline functions for computing
  approximate local surface normals on polygons and a PolygonMesh,
  respectively, and approximate covariance matrices given points and
  normals (GICP-style).
  • Loading branch information
Dirk Holz committed Jun 16, 2015
1 parent 5ad8000 commit 3ee92d9
Show file tree
Hide file tree
Showing 6 changed files with 301 additions and 39 deletions.
1 change: 1 addition & 0 deletions features/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ if(build)
"include/pcl/${SUBSYS_NAME}/feature.h"
"include/pcl/${SUBSYS_NAME}/fpfh.h"
"include/pcl/${SUBSYS_NAME}/fpfh_omp.h"
"include/pcl/${SUBSYS_NAME}/from_meshes.h"
"include/pcl/${SUBSYS_NAME}/gfpfh.h"
"include/pcl/${SUBSYS_NAME}/integral_image2D.h"
"include/pcl/${SUBSYS_NAME}/integral_image_normal.h"
Expand Down
105 changes: 105 additions & 0 deletions features/include/pcl/features/from_meshes.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
#ifndef PCL_FEATURES_FROM_MESHES_H_
#define PCL_FEATURES_FROM_MESHES_H_

#include <pcl/features/normal_3d.h>

namespace pcl
{
namespace features
{

/** \brief Compute approximate surface normals on a mesh.
* \param[in] cloud Point cloud containing the XYZ coordinates.
* \param[in] polygons Polygons from the mesh.
* \param[out] normals Point cloud with computed surface normals
*/
template <typename PointT, typename PointNT> inline void
computeApproximateNormals(const pcl::PointCloud<PointT>& cloud, const std::vector<pcl::Vertices>& polygons, pcl::PointCloud<PointNT>& normals)
{
int nr_points = static_cast<int>(cloud.points.size());
int nr_polygons = static_cast<int>(polygons.size());

normals.header = cloud.header;
normals.width = cloud.width;
normals.height = cloud.height;
normals.points.resize(nr_points);

for ( int i = 0; i < nr_points; ++i )
normals.points[i].getNormalVector3fMap() = Eigen::Vector3f::Zero();

// NOTE: for efficiency the weight is computed implicitly by using the
// cross product, this causes inaccurate normals for meshes containing
// non-triangle polygons (quads or other types)
for ( int i = 0; i < nr_polygons; ++i )
{
const int nr_points_polygon = (int)polygons[i].vertices.size();
if (nr_points_polygon < 3) continue;

// compute normal for triangle
Eigen::Vector3f vec_a_b = cloud.points[polygons[i].vertices[0]].getVector3fMap() - cloud.points[polygons[i].vertices[1]].getVector3fMap();
Eigen::Vector3f vec_a_c = cloud.points[polygons[i].vertices[0]].getVector3fMap() - cloud.points[polygons[i].vertices[2]].getVector3fMap();
Eigen::Vector3f normal = vec_a_b.cross(vec_a_c);
pcl::flipNormalTowardsViewpoint(cloud.points[polygons[i].vertices[0]], 0.0f, 0.0f, 0.0f, normal(0), normal(1), normal(2));

// add normal to all points in polygon
for ( int j = 0; j < nr_points_polygon; ++j )
normals.points[polygons[i].vertices[j]].getNormalVector3fMap() += normal;
}

for ( int i = 0; i < nr_points; ++i )
{
normals.points[i].getNormalVector3fMap().normalize();
pcl::flipNormalTowardsViewpoint(cloud.points[i], 0.0f, 0.0f, 0.0f, normals.points[i].normal_x, normals.points[i].normal_y, normals.points[i].normal_z);
}
}


/** \brief Compute GICP-style covariance matrices given a point cloud and
* the corresponding surface normals.
* \param[in] cloud Point cloud containing the XYZ coordinates,
* \param[in] normals Point cloud containing the corresponding surface normals.
* \param[out] covariances Vector of computed covariances.
* \param[in] Optional: Epsilon for the expected noise along the surface normal (default: 0.001)
*/
template <typename PointT, typename PointNT> inline void
computeApproximateCovariances(const pcl::PointCloud<PointT>& cloud,
const pcl::PointCloud<PointNT>& normals,
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> >& covariances,
double epsilon = 0.001)
{
assert(cloud.points.size() == normals.points.size());

int nr_points = static_cast<int>(cloud.points.size());
covariances.resize(nr_points);
for (int i = 0; i < nr_points; ++i)
{
Eigen::Vector3d normal(normals.points[i].normal_x,
normals.points[i].normal_y,
normals.points[i].normal_z);

// compute rotation matrix
Eigen::Matrix3d rot;
Eigen::Vector3d y;
y << 0, 1, 0;
rot.row(2) = normal;
y = y - normal(1) * normal;
y.normalize();
rot.row(1) = y;
rot.row(0) = normal.cross(rot.row(1));

// comnpute approximate covariance
Eigen::Matrix3d cov;
cov << 1, 0, 0,
0, 1, 0,
0, 0, epsilon;
covariances[i] = rot.transpose()*cov*rot;
}
}

}
}


#endif // PCL_FEATURES_FROM_MESHES_H_


42 changes: 32 additions & 10 deletions registration/include/pcl/registration/gicp.h
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,10 @@ namespace pcl
typedef PointIndices::Ptr PointIndicesPtr;
typedef PointIndices::ConstPtr PointIndicesConstPtr;

typedef std::vector< Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> > MatricesVector;
typedef boost::shared_ptr< MatricesVector > MatricesVectorPtr;
typedef boost::shared_ptr< const MatricesVector > MatricesVectorConstPtr;

typedef typename Registration<PointSource, PointTarget>::KdTree InputKdTree;
typedef typename Registration<PointSource, PointTarget>::KdTreePtr InputKdTreePtr;

Expand All @@ -104,8 +108,6 @@ namespace pcl
: k_correspondences_(20)
, gicp_epsilon_(0.001)
, rotation_epsilon_(2e-3)
, input_covariances_(0)
, target_covariances_(0)
, mahalanobis_(0)
, max_inner_iterations_(20)
{
Expand Down Expand Up @@ -144,21 +146,41 @@ namespace pcl
input[i].data[3] = 1.0;

pcl::IterativeClosestPoint<PointSource, PointTarget>::setInputSource (cloud);
input_covariances_.clear ();
input_covariances_.reserve (input_->size ());
input_covariances_.reset ();
}

/** \brief Provide a pointer to the covariances of the input source (if computed externally!).
* If not set, GeneralizedIterativeClosestPoint will compute the covariances itself.
* Make sure to set the covariances AFTER setting the input source point cloud (setting the input source point cloud will reset the covariances).
* \param[in] target the input point cloud target
*/
inline void
setSourceCovariances (const MatricesVectorPtr& covariances)
{
input_covariances_ = covariances;
}

/** \brief Provide a pointer to the input target (e.g., the point cloud that we want to align the input source to)
* \param[in] target the input point cloud target
*/
inline void
setInputTarget (const PointCloudTargetConstPtr &target)
{
pcl::IterativeClosestPoint<PointSource, PointTarget>::setInputTarget(target);
target_covariances_.clear ();
target_covariances_.reserve (target_->size ());
target_covariances_.reset ();
}

/** \brief Provide a pointer to the covariances of the input target (if computed externally!).
* If not set, GeneralizedIterativeClosestPoint will compute the covariances itself.
* Make sure to set the covariances AFTER setting the input source point cloud (setting the input source point cloud will reset the covariances).
* \param[in] target the input point cloud target
*/
inline void
setTargetCovariances (const MatricesVectorPtr& covariances)
{
target_covariances_ = covariances;
}

/** \brief Estimate a rigid rotation transformation between a source and a target point cloud using an iterative
* non-linear Levenberg-Marquardt approach.
* \param[in] cloud_src the source point cloud dataset
Expand Down Expand Up @@ -266,13 +288,13 @@ namespace pcl


/** \brief Input cloud points covariances. */
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> > input_covariances_;
MatricesVectorPtr input_covariances_;

/** \brief Target cloud points covariances. */
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> > target_covariances_;
MatricesVectorPtr target_covariances_;

/** \brief Mahalanobis matrices holder. */
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> > mahalanobis_;
std::vector<Eigen::Matrix3d> mahalanobis_;

/** \brief maximum number of optimizations */
int max_inner_iterations_;
Expand All @@ -286,7 +308,7 @@ namespace pcl
template<typename PointT>
void computeCovariances(typename pcl::PointCloud<PointT>::ConstPtr cloud,
const typename pcl::search::KdTree<PointT>::Ptr tree,
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> >& cloud_covariances);
MatricesVector& cloud_covariances);

/** \return trace of mat1^t . mat2
* \param mat1 matrix of dimension nxm
Expand Down
22 changes: 14 additions & 8 deletions registration/include/pcl/registration/impl/gicp.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ template <typename PointSource, typename PointTarget>
template<typename PointT> void
pcl::GeneralizedIterativeClosestPoint<PointSource, PointTarget>::computeCovariances(typename pcl::PointCloud<PointT>::ConstPtr cloud,
const typename pcl::search::KdTree<PointT>::Ptr kdtree,
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> >& cloud_covariances)
MatricesVector& cloud_covariances)
{
if (k_correspondences_ > int (cloud->size ()))
{
Expand All @@ -73,7 +73,7 @@ pcl::GeneralizedIterativeClosestPoint<PointSource, PointTarget>::computeCovarian
cloud_covariances.resize (cloud->size ());

typename pcl::PointCloud<PointT>::const_iterator points_iterator = cloud->begin ();
std::vector<Eigen::Matrix3d, Eigen::aligned_allocator<Eigen::Matrix3d> >::iterator matrices_iterator = cloud_covariances.begin ();
MatricesVector::iterator matrices_iterator = cloud_covariances.begin ();
for(;
points_iterator != cloud->end ();
++points_iterator, ++matrices_iterator)
Expand Down Expand Up @@ -361,11 +361,17 @@ pcl::GeneralizedIterativeClosestPoint<PointSource, PointTarget>::computeTransfor
// Set the mahalanobis matrices to identity
mahalanobis_.resize (N, Eigen::Matrix3d::Identity ());
// Compute target cloud covariance matrices
if (target_covariances_.empty ())
computeCovariances<PointTarget> (target_, tree_, target_covariances_);
if ((!target_covariances_) || (target_covariances_->empty ()))
{
target_covariances_.reset (new MatricesVector);
computeCovariances<PointTarget> (target_, tree_, *target_covariances_);
}
// Compute input cloud covariance matrices
if (input_covariances_.empty ())
computeCovariances<PointSource> (input_, tree_reciprocal_, input_covariances_);
if ((!input_covariances_) || (input_covariances_->empty ()))
{
input_covariances_.reset (new MatricesVector);
computeCovariances<PointSource> (input_, tree_reciprocal_, *input_covariances_);
}

base_transformation_ = guess;
nr_iterations_ = 0;
Expand Down Expand Up @@ -404,8 +410,8 @@ pcl::GeneralizedIterativeClosestPoint<PointSource, PointTarget>::computeTransfor
// Check if the distance to the nearest neighbor is smaller than the user imposed threshold
if (nn_dists[0] < dist_threshold)
{
Eigen::Matrix3d &C1 = input_covariances_[i];
Eigen::Matrix3d &C2 = target_covariances_[nn_indices[0]];
Eigen::Matrix3d &C1 = (*input_covariances_)[i];
Eigen::Matrix3d &C2 = (*target_covariances_)[nn_indices[0]];
Eigen::Matrix3d &M = mahalanobis_[i];
// M = R*C1
M = R * C1;
Expand Down
8 changes: 4 additions & 4 deletions registration/src/gicp6d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -191,10 +191,10 @@ namespace pcl
mahalanobis_.resize (N, Eigen::Matrix3d::Identity ());

// Compute target cloud covariance matrices
computeCovariances<PointTarget> (target_, tree_, target_covariances_);
computeCovariances<PointTarget> (target_, tree_, *target_covariances_);
// Compute input cloud covariance matrices
computeCovariances<PointSource> (input_, tree_reciprocal_,
input_covariances_);
*input_covariances_);

base_transformation_ = guess;
nr_iterations_ = 0;
Expand Down Expand Up @@ -237,8 +237,8 @@ namespace pcl
// Check if the distance to the nearest neighbor is smaller than the user imposed threshold
if (nn_dists[0] < dist_threshold)
{
Eigen::Matrix3d &C1 = input_covariances_[i];
Eigen::Matrix3d &C2 = target_covariances_[nn_indices[0]];
Eigen::Matrix3d &C1 = (*input_covariances_)[i];
Eigen::Matrix3d &C2 = (*target_covariances_)[nn_indices[0]];
Eigen::Matrix3d &M = mahalanobis_[i];
// M = R*C1
M = R * C1;
Expand Down
Loading

0 comments on commit 3ee92d9

Please sign in to comment.