PCL - MLS代碼研讀（十三）- RANDOM_UNIFORM_DENSITY上採樣方法

PCL - MLS代碼研讀（十三）- RANDOM_UNIFORM_DENSITY上採樣方法

前言

在PCL - MLS代碼研讀（十一）- computeMLSPointNormal函數中已經介紹了NONE和SAMPLE_LOCAL_PLANE上採樣方法。本篇主要介紹之前略過的RANDOM_UNIFORM_DENSITY上採樣方法。

成員變數

protected成員變數：

/** \brief Parameter that specifies the desired number of points within the search radius
  * \note Used only in the case of RANDOM_UNIFORM_DENSITY upsampling
  */
int desired_num_points_in_radius_;

desired_num_points_in_radius_指定在搜索半徑內需要達到的點的數量。

private成員變數：

/** \brief Random number generator algorithm. */
mutable std::mt19937 rng_;

/** \brief Random number generator using an uniform distribution of floats
  * \note Used only in the case of RANDOM_UNIFORM_DENSITY upsampling
  */
std::unique_ptr<std::uniform_real_distribution<>> rng_uniform_distribution_;

其中std::mt19937是一種新型的亂數產生器，可以生成品質較高的亂數。std::uniform_real_distribution則代表我們想要的亂數分布是均勻分布。他們的詳細介紹可以參考C++11 內建亂數函式庫使用教學：隨機亂數產生器與機率分布。

getter & setter

/** \brief Set the parameter that specifies the desired number of points within the search radius
  * \note Used only in the case of RANDOM_UNIFORM_DENSITY upsampling
  * \param[in] desired_num_points_in_radius the desired number of points in the output cloud in a sphere of
  * radius \ref search_radius_ around each point
  */
inline void
setPointDensity (int desired_num_points_in_radius) { desired_num_points_in_radius_ = desired_num_points_in_radius; }


/** \brief Get the parameter that specifies the desired number of points within the search radius
  * \note Used only in the case of RANDOM_UNIFORM_DENSITY upsampling
  */
inline int
getPointDensity () const { return (desired_num_points_in_radius_); }

上面這兩個函數是desired_num_points_in_radius成員變數的setter和getter。

process函數

//
template <typename PointInT, typename PointOutT> void
pcl::MovingLeastSquares<PointInT, PointOutT>::process (PointCloudOut &output)
{
  //...

  //如果使用NONE以外的上採樣方法,那麼就得對各自需要的成員變數做初始化
  switch (upsample_method_)
  {
    // Initialize random number generator if necessary
    case (RANDOM_UNIFORM_DENSITY):
    {
      std::random_device rd;
      rng_.seed (rd());
      //在(-search_radius_ / 2.0,search_radius_ / 2.0)範圍內均勻隨機生成一個數字
      const double tmp = search_radius_ / 2.0;
      rng_uniform_distribution_.reset (new std::uniform_real_distribution<> (-tmp, tmp));

      break;
    }
    case (VOXEL_GRID_DILATION):
    case (DISTINCT_CLOUD):
    {
      //...
    }
    default:
      break;
  }
  //...

  // Perform the actual surface reconstruction
  performProcessing (output);

  //...
}

初始化成員變數rng_。然後初始化rng_uniform_distribution_，讓它隨機在[-search_radius_ / 2.0, search_radius_ / 2.0]的範圍內均勻地生成浮點數。不懂這裡的取值範圍為何不是[-search_radius, search_radius]。

可以看到在process函數的最後，調用了performProcessing函數。在performProcessing函數中，進一步調用了computeMLSPointNormal函數。

template <typename PointInT, typename PointOutT> void
pcl::MovingLeastSquares<PointInT, PointOutT>::performProcessing (PointCloudOut &output)
{
  // ...
  for (int cp = 0; cp < static_cast<int> (indices_->size ()); ++cp)
  {
    // ...

    // Get the initial estimates of point positions and their neighborhoods
    if (searchForNeighbors ((*indices_)[cp], nn_indices, nn_sqr_dists))
    {
      // Check the number of nearest neighbors for normal estimation (and later for polynomial fit as well)
      if (nn_indices.size () >= 3)
      {
        // ...
        computeMLSPointNormal (index, nn_indices, projected_points, projected_points_normals, *corresponding_input_indices_, mls_results_[mls_result_index]);
      }
    }
  }

  //...
}

computeMLSPointNormal函數

在computeMLSPointNormal函數中，會先計算上採樣的目標點數為num_points_to_add，然後由rng_和rng_uniform_distribution_在範圍內生成點，對點做投影後，再調用addProjectedPointNormal函數將投影點加入輸出點雲中。

computeMLSPointNormal函數中跟RANDOM_UNIFORM_DENSITY上採樣方法相關的部分如下：

case (RANDOM_UNIFORM_DENSITY):
{

首先計算需要填入的點數：

// Compute the local point density and add more samples if necessary
//半徑為何還要除以2?
//為什麼不是desired_num_points_in_radius_-nn_indices.size ()
//desired_num_points_in_radius與nn_indices共用半徑
const int num_points_to_add = static_cast<int> (std::floor (desired_num_points_in_radius_ / 2.0 / static_cast<double> (nn_indices.size ())));

==為何是用除的不是用減的？==以下為筆者的猜測：假設每一個點在上採樣後會變成x個點，所以半徑內nn_indices.size()個點會變成nn_indices.size()乘x個點，如果取x為desired_num_points_in_radius_/nn_indices.size()，那麼上採樣後搜索半徑內的點數就會約略等於desired_num_points_in_radius_。代碼中的/2.0是做什麼用的？

如果不需要做上採樣，則對query_point投影後將它加入projected_points數據結構中：

// Just add the query point, because the density is good
if (num_points_to_add <= 0)
{
  // Just add the current point
  const MLSResult::MLSProjectionResults proj = mls_result.projectQueryPoint (projection_method_, nr_coeff_);
  addProjectedPointNormal (index, proj.point, proj.normal, mls_result.curvature, projected_points, projected_points_normals, corresponding_input_indices);
}

如果需要做上採樣，則每次在query_point附近[-search_radius_ / 2.0, search_radius_ / 2.0]的範圍內隨機生成一個點，將它投影到曲面或平面上，最後再把投影點加入projected_points，projected_points_normals等數據結構中。

else
{
  // Sample the local plane
  for (int num_added = 0; num_added < num_points_to_add;)
  {
    //在(-search_radius_ / 2.0,search_radius_ / 2.0)範圍內均勻隨機生成一個數字
    const double u = (*rng_uniform_distribution_) (rng_);
    const double v = (*rng_uniform_distribution_) (rng_);

    // Check if inside circle; if not, try another coin flip
    // 如果把"search_radius/2"當成半徑,那麼這句就是合理的
    // 這邊continue後,就會跳過最底下的"num_added++"
    if (u * u + v * v > search_radius_ * search_radius_ / 4)
      continue;

    MLSResult::MLSProjectionResults proj;
    //為何不做跟projectPointSimpleToPolynomialSurface一樣的檢查:
    //if (order > 1 && c_vec.size () >= (order + 1) * (order + 2) / 2 && std::isfinite (c_vec[0]))
    if (order_ > 1 && mls_result.num_neighbors >= 5 * nr_coeff_)
      proj = mls_result.projectPointSimpleToPolynomialSurface (u, v);
    else
      proj = mls_result.projectPointToMLSPlane (u, v);

    addProjectedPointNormal (index, proj.point, proj.normal, mls_result.curvature, projected_points, projected_points_normals, corresponding_input_indices);

    num_added++;
  }
}
break;

上採樣結果

圖一是原始點雲，圖二是上採樣後的點雲。