18 #ifndef POINT_CLOUD_COMMON_UTILITIES_H_ 
   19 #define POINT_CLOUD_COMMON_UTILITIES_H_ 
   25 #include <gtsam/geometry/Pose3.h> 
   26 #include <gtsam/geometry/Point3.h> 
   28 #include <pcl/features/fpfh.h> 
   29 #include <pcl/features/integral_image_normal.h> 
   30 #include <pcl/features/normal_3d.h> 
   31 #include <pcl/features/impl/normal_3d.hpp> 
   32 #include <pcl/filters/fast_bilateral.h> 
   33 #include <pcl/filters/impl/fast_bilateral.hpp> 
   34 #include <pcl/filters/normal_space.h> 
   35 #include <pcl/filters/voxel_grid.h> 
   36 #include <pcl/kdtree/impl/kdtree_flann.hpp> 
   37 #include <pcl/point_cloud.h> 
   38 #include <pcl/point_types.h> 
   39 #include <pcl/search/impl/kdtree.hpp> 
   40 #include <pcl/search/impl/search.hpp> 
   46 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
   47 void EstimateNormals(
const typename pcl::PointCloud<PointType>::Ptr cloud, 
const double search_radius,
 
   48                      pcl::PointCloud<PointWithNormalType>& cloud_with_normals);
 
   50 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
   52   const typename pcl::PointCloud<PointType>::Ptr cloud,
 
   53   const typename pcl::IntegralImageNormalEstimation<PointType, pcl::Normal>::NormalEstimationMethod method,
 
   54   const bool depth_dependent_smoothing, 
const double max_depth_change_factor, 
const double normal_smoothing_size,
 
   55   pcl::PointCloud<PointWithNormalType>& cloud_with_normals);
 
   57 template <
typename Po
intType>
 
   59                             const int num_samples);
 
   61 Eigen::Matrix4f 
RansacIA(
const pcl::PointCloud<pcl::PointXYZINormal>::Ptr source_cloud,
 
   62                          const pcl::PointCloud<pcl::PointXYZINormal>::Ptr target_cloud);
 
   65   const pcl::PointCloud<pcl::PointXYZINormal>::Ptr cloud_with_normals);
 
   67 Eigen::Matrix<double, 1, 6> 
PointToPlaneJacobian(
const gtsam::Point3& source_point, 
const gtsam::Vector3& normal,
 
   68                                                  const gtsam::Pose3& target_T_source);
 
   71                                                  const gtsam::Pose3& target_T_source);
 
   74                                            const std::vector<Eigen::Vector3d>& target_points);
 
   76 template <
typename Po
intType>
 
   78                                            const pcl::search::KdTree<PointType>& kdtree,
 
   79                                            const double search_radius = 0.03);
 
   81 template <
typename Po
intType>
 
   82 bool computePointNormal(
const pcl::PointCloud<PointType>& cloud, 
const std::vector<int>& indices, 
float& normal_x,
 
   83                         float& normal_y, 
float& normal_z, 
float& curvature);
 
   85 boost::optional<Eigen::Matrix<double, 6, 6>> 
PointToPointCovariance(
const std::vector<Eigen::Vector3d>& source_points,
 
   88 boost::optional<Eigen::Matrix<double, 6, 6>> 
PointToPlaneCovariance(
const std::vector<Eigen::Vector3d>& source_points,
 
   89                                                                     const std::vector<Eigen::Vector3d>& target_normals,
 
   92 template <
typename Po
intType>
 
   95 template <
typename Po
intType>
 
   98 pcl::PointXYZI 
Interpolate(
const pcl::PointXYZI& point_a, 
const pcl::PointXYZI& point_b, 
const double alpha);
 
  100 template <
typename Po
intType>
 
  101 bool ValidPoint(
const PointType& point) = 
delete;
 
  104 bool ValidPoint<pcl::PointXYZ>(
const pcl::PointXYZ& point);
 
  107 bool ValidPoint<pcl::PointXYZI>(
const pcl::PointXYZI& point);
 
  110 bool ValidPoint<pcl::PointNormal>(
const pcl::PointNormal& point);
 
  113 bool ValidPoint<pcl::PointXYZINormal>(
const pcl::PointXYZINormal& point);
 
  115 template <
typename Type>
 
  116 bool ApproxZero(
const Type& point, 
const double epsilon = 1e-5) {
 
  117   return std::abs(point) <= epsilon;
 
  120 template <
typename Po
intXYZType>
 
  123 template <
typename Po
intXYZType>
 
  126 template <
typename Po
intNormalType>
 
  129 template <
typename Po
intIntensityType>
 
  132 template <
typename Po
intType>
 
  137 template <
typename Po
intType>
 
  140 template <
typename Po
intType>
 
  141 typename pcl::PointCloud<PointType>::Ptr 
DownsamplePointCloud(
const typename pcl::PointCloud<PointType>::Ptr cloud,
 
  142                                                               const double leaf_size);
 
  144 template <
typename Po
intType>
 
  146   const typename pcl::PointCloud<PointType>::Ptr cloud, 
const double sigma_s, 
const double sigma_r);
 
  147 template <
typename Po
intType>
 
  149                            const typename pcl::PointCloud<PointType>& target_cloud,
 
  150                            pcl::Correspondences& correspondences);
 
  152 template <
typename Po
intType>
 
  154   const typename pcl::PointCloud<PointType>::Ptr unfiltered_cloud);
 
  156 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
  158   const typename pcl::PointCloud<PointType>::Ptr unfiltered_cloud, 
const double search_radius);
 
  161 template <
typename Po
intXYZType>
 
  164   return nonzero_point;
 
  167 template <
typename Po
intXYZType>
 
  169   const bool finite_point = pcl_isfinite(point.x) && pcl_isfinite(point.y) && pcl_isfinite(point.z);
 
  171   return finite_point && nonzero_point;
 
  174 template <
typename Po
intNormalType>
 
  176   const bool finite_normal =
 
  177     pcl_isfinite(point.normal_x) && pcl_isfinite(point.normal_y) && pcl_isfinite(point.normal_z);
 
  179   return finite_normal && nonzero_normal;
 
  182 template <
typename Po
intIntensityType>
 
  184   return pcl_isfinite(point.intensity);
 
  187 template <
typename Po
intType>
 
  189   size_t new_index = 0;
 
  190   for (
const auto& point : cloud.points) {
 
  192     if (!valid_point) 
continue;
 
  193     cloud.points[new_index++] = point;
 
  195   if (new_index != cloud.points.size()) {
 
  196     cloud.points.resize(new_index);
 
  200   cloud.width = 
static_cast<uint32_t
>(new_index);
 
  201   cloud.is_dense = 
true;
 
  204 template <
typename Po
intType>
 
  206   for (
auto& point : cloud.points) {
 
  211       point.x = std::numeric_limits<double>::quiet_NaN();
 
  212       point.y = std::numeric_limits<double>::quiet_NaN();
 
  213       point.z = std::numeric_limits<double>::quiet_NaN();
 
  218 template <
typename Po
intType>
 
  219 typename pcl::PointCloud<PointType>::Ptr 
DownsamplePointCloud(
const typename pcl::PointCloud<PointType>::Ptr cloud,
 
  220                                                               const double leaf_size) {
 
  221   pcl::VoxelGrid<PointType> voxel_grid;
 
  222   voxel_grid.setInputCloud(cloud);
 
  223   voxel_grid.setLeafSize(leaf_size, leaf_size, leaf_size);
 
  224   typename pcl::PointCloud<PointType>::Ptr downsampled_cloud(
new pcl::PointCloud<PointType>());
 
  225   voxel_grid.filter(*downsampled_cloud);
 
  226   return downsampled_cloud;
 
  229 template <
typename Po
intType>
 
  231   const typename pcl::PointCloud<PointType>::Ptr cloud, 
const double sigma_s, 
const double sigma_r) {
 
  232   typename pcl::PointCloud<PointType>::Ptr downsampled_cloud(
new pcl::PointCloud<PointType>());
 
  233   typename pcl::FastBilateralFilter<PointType> bilateral_filter;
 
  234   bilateral_filter.setInputCloud(cloud);
 
  235   bilateral_filter.setSigmaS(sigma_s);
 
  236   bilateral_filter.setSigmaR(sigma_r);
 
  237   bilateral_filter.filter(*downsampled_cloud);
 
  238   return downsampled_cloud;
 
  241 template <
typename Po
intType>
 
  243                            const typename pcl::PointCloud<PointType>& target_cloud,
 
  244                            pcl::Correspondences& correspondences) {
 
  245   for (
auto correspondence_it = correspondences.begin(); correspondence_it != correspondences.end();) {
 
  246     const auto& input_point = (input_cloud)[correspondence_it->index_query];
 
  247     const auto& target_point = (target_cloud)[correspondence_it->index_match];
 
  248     const bool invalid_correspondence = !std::isfinite(input_point.x) || !std::isfinite(input_point.y) ||
 
  249                                         !std::isfinite(input_point.z) || !std::isfinite(target_point.x) ||
 
  250                                         !std::isfinite(target_point.y) || !std::isfinite(target_point.z) ||
 
  251                                         !std::isfinite(target_point.normal_x) ||
 
  252                                         !std::isfinite(target_point.normal_y) || !std::isfinite(target_point.normal_z);
 
  253     if (invalid_correspondence) {
 
  254       correspondence_it = correspondences.erase(correspondence_it);
 
  261 template <
typename Po
intType>
 
  263   const typename pcl::PointCloud<PointType>::Ptr unfiltered_cloud) {
 
  264   typename pcl::PointCloud<PointType>::Ptr filtered_cloud(
new pcl::PointCloud<PointType>());
 
  265   pcl::copyPointCloud(*unfiltered_cloud, *filtered_cloud);
 
  267   return filtered_cloud;
 
  270 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
  271 void EstimateNormals(
const typename pcl::PointCloud<PointType>::Ptr cloud, 
const double search_radius,
 
  272                      typename pcl::PointCloud<PointWithNormalType>& cloud_with_normals) {
 
  273   typename pcl::NormalEstimation<PointType, pcl::Normal> ne;
 
  274   ne.setInputCloud(cloud);
 
  275   typename pcl::search::KdTree<PointType>::Ptr tree(
new pcl::search::KdTree<PointType>());
 
  276   ne.setSearchMethod(tree);
 
  277   ne.setRadiusSearch(search_radius);
 
  278   pcl::PointCloud<pcl::Normal> cloud_normals;
 
  279   ne.compute(cloud_normals);
 
  280   pcl::concatenateFields(*cloud, cloud_normals, cloud_with_normals);
 
  283 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
  285   const typename pcl::PointCloud<PointType>::Ptr cloud,
 
  286   const typename pcl::IntegralImageNormalEstimation<PointType, pcl::Normal>::NormalEstimationMethod method,
 
  287   const bool depth_dependent_smoothing, 
const double max_depth_change_factor, 
const double normal_smoothing_size,
 
  288   typename pcl::PointCloud<PointWithNormalType>& cloud_with_normals) {
 
  289   typename pcl::IntegralImageNormalEstimation<PointType, pcl::Normal> ne;
 
  290   ne.setNormalEstimationMethod(method);
 
  291   ne.setMaxDepthChangeFactor(max_depth_change_factor);
 
  292   ne.setNormalSmoothingSize(normal_smoothing_size);
 
  293   ne.setDepthDependentSmoothing(depth_dependent_smoothing);
 
  294   ne.setInputCloud(cloud);
 
  299   ne.setSearchMethod(organized_neighbor);
 
  300   pcl::PointCloud<pcl::Normal> cloud_normals;
 
  301   ne.compute(cloud_normals);
 
  302   pcl::concatenateFields(*cloud, cloud_normals, cloud_with_normals);
 
  305 template <
typename Po
intType>
 
  307                             const int num_samples) {
 
  308   typename pcl::NormalSpaceSampling<PointType, PointType> normal_space_sampling;
 
  309   normal_space_sampling.setInputCloud(cloud);
 
  310   normal_space_sampling.setNormals(cloud);
 
  311   normal_space_sampling.setBins(bins_per_axis, bins_per_axis, bins_per_axis);
 
  312   normal_space_sampling.setSample(num_samples);
 
  313   normal_space_sampling.filter(*cloud);
 
  316 template <
typename Po
intType, 
typename Po
intWithNormalType>
 
  318   const typename pcl::PointCloud<PointType>::Ptr unfiltered_cloud, 
const double search_radius) {
 
  319   typename pcl::PointCloud<PointType>::Ptr filtered_cloud = FilteredPointCloud<PointType>(unfiltered_cloud);
 
  320   typename pcl::PointCloud<PointWithNormalType>::Ptr filtered_cloud_with_normals(
 
  321     new pcl::PointCloud<PointWithNormalType>());
 
  322   EstimateNormals<PointType, PointWithNormalType>(filtered_cloud, search_radius, *filtered_cloud_with_normals);
 
  324   return filtered_cloud_with_normals;
 
  327 template <
typename Po
intType>
 
  332 template <
typename Po
intType>
 
  334   return Eigen::Vector3d(point_with_normal.normal[0], point_with_normal.normal[1], point_with_normal.normal[2]);
 
  337 template <
typename Po
intType>
 
  339                                            const pcl::search::KdTree<PointType>& kdtree, 
const double search_radius) {
 
  342   pcl_point.x = point.x();
 
  343   pcl_point.y = point.y();
 
  344   pcl_point.z = point.z();
 
  346   std::vector<int> nn_indices;
 
  347   std::vector<float> distances;
 
  348   if (kdtree.radiusSearch(pcl_point, search_radius, nn_indices, distances, 0) < 3) {
 
  349     LogDebug(
"GetNormal: Failed to get enough neighboring points for query point.");
 
  357   if (!
computePointNormal(cloud, nn_indices, normal_x, normal_y, normal_z, curvature)) {
 
  358     LogDebug(
"GetNormal: Failed to compute point normal.");
 
  362   const double vpx = cloud.sensor_origin_.coeff(0);
 
  363   const double vpy = cloud.sensor_origin_.coeff(1);
 
  364   const double vpz = cloud.sensor_origin_.coeff(2);
 
  365   flipNormalTowardsViewpoint(pcl_point, vpx, vpy, vpz, normal_x, normal_y, normal_z);
 
  369 template <
typename Po
intType>
 
  370 bool computePointNormal(
const pcl::PointCloud<PointType>& cloud, 
const std::vector<int>& indices, 
float& normal_x,
 
  371                         float& normal_y, 
float& normal_z, 
float& curvature) {
 
  374   static EIGEN_ALIGN16 Eigen::Matrix3f covariance_matrix_;
 
  376   static Eigen::Vector4f xyz_centroid_;
 
  377   if (indices.size() < 3 ||
 
  378       pcl::computeMeanAndCovarianceMatrix(cloud, indices, covariance_matrix_, xyz_centroid_) == 0) {
 
  383   pcl::solvePlaneParameters(covariance_matrix_, normal_x, normal_y, normal_z, curvature);
 
  387 #endif  // POINT_CLOUD_COMMON_UTILITIES_H_