[VoxelGrid] Downsampling produces wrong result

[REASY]

If my assumption that for the same input PointCloud, monotonically increasing voxel size should produce monotonically decreased downsampled PointCloud is correct then something wrong with the downsampling in VoxelGrid.

Let's take an example, the following code creates point cloud and voxelizes it with three different voxel sizes: 3.0f, 3.33f, 4.0f.

#include <pcl/test/gtest.h>
#include <pcl/point_types.h>
#include <pcl/filters/filter.h>
#include <pcl/filters/voxel_grid.h>

using namespace pcl;
using namespace pcl::io;
using namespace Eigen;

TEST (VoxelGrid, BugInFilter)
{
  const size_t size = 1000;
  float x = 0;
  float y = 360;
  float z = 720;
  PointCloud<PointXYZ> pcl {};
  for (size_t i = 0; i < size; ++i) {
    pcl.push_back(pcl::PointXYZ{x, y, z});
    x = x + 0.1f;
    y = y + 0.1f;
    z = z + 0.1f;
  }
  PointCloud<PointXYZ>::Ptr pclPtr(new PointCloud<PointXYZ>(pcl));

  VoxelGrid<PointXYZ> grid;

  std::array<float, 3> voxelSizes =  { 3.0f, 3.33f, 4.0f};
  for(const float& voxelSize: voxelSizes) {
    PointCloud<PointXYZ> output;
    grid.setLeafSize (voxelSize, voxelSize, voxelSize);
    grid.setInputCloud (pclPtr);
    grid.filter (output);

    PointXYZ& prev = output[0];
    double dx = 0;
    double dy = 0;
    double dz = 0;
    for (std::size_t i = 1; i < output.size(); ++i) {
      const PointXYZ& curr = output[i];
      dx += std::abs(curr.x - prev.x);
      dy += std::abs(curr.y - prev.y);
      dz += std::abs(curr.z - prev.z);
      prev = curr;
    }
    double avgDx = dx / output.size();
    double avgDy = dy / output.size();
    double avgDz = dz / output.size();
    std::cout << "voxelSize = " << voxelSize << ", downsampled size: " << output.size() << ", averate dx: " << avgDx
              << ", average dy: " << avgDy << ", average dz: " << avgDz<< std::endl;
  }
}

One would expect that voxelized(3.0f).size() > voxelized(3.33f).size() > voxelized(4.0f).size(), but that is not the case. (voxelized(3.0f).size() means use VoxelGrid to voxelize with leaf size 3.0 and get the size of output). It produces the following result:

voxelSize = 3, downsampled size: 67, averate dx: 1.46342, average dy: 1.46352, average dz: 1.46308
voxelSize = 3.33, downsampled size: 90, averate dx: 1.09388, average dy: 1.09396, average dz: 1.09362
voxelSize = 4, downsampled size: 49, averate dx: 1.96019, average dy: 1.96033, average dz: 1.95973

Is it expected that the average difference between x/y/z of voxelized points is around voxelSize / 2, is it a mistake to expect it to be around voxelSize?

Environment:

• OS: Ubuntu 20.04
• Compiler: c++ (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
• PCL Version: 7d52b10

[mvieth]

The data you are using seems to be the problem. I extended your test:

#include <pcl/point_types.h>
#include <pcl/filters/filter.h>
#include <pcl/filters/voxel_grid.h>

using namespace pcl;
using namespace pcl::io;
using namespace Eigen;

int main()
{
  const size_t size = 1000;
#if 0
  PointCloud<PointXYZ> pcl {};
  for (size_t i = 0; i < size; ++i) {
    pcl.push_back(pcl::PointXYZ{10.0f * static_cast<float> (rand ()) / RAND_MAX,
                                10.0f * static_cast<float> (rand ()) / RAND_MAX,
                                10.0f * static_cast<float> (rand ()) / RAND_MAX});
  }
#else
  float x = 0;
  float y = 360;
  float z = 720;
  PointCloud<PointXYZ> pcl {};
  for (size_t i = 0; i < size; ++i) {
    pcl.push_back(pcl::PointXYZ{x, y, z});
    x = x + 0.1f;
    y = y + 0.1f;
    z = z + 0.1f;
  }
#endif
  PointCloud<PointXYZ>::Ptr pclPtr(new PointCloud<PointXYZ>(pcl));

  std::vector<float> voxelSizes =  {2.0123f, 2.3423f, 2.6723f, 3.0123f, 3.3423f, 3.6723f, 4.0123f, 4.3423f, 4.6723f};
  for(const float& voxelSize: voxelSizes) {
    PointCloud<PointXYZ> output;
    VoxelGrid<PointXYZ> grid;
    grid.setLeafSize (voxelSize, voxelSize, voxelSize);
    grid.setInputCloud (pclPtr);
    grid.filter (output);

    PointXYZ& prev = output[0];
    double dx = 0;
    double dy = 0;
    double dz = 0;
    for (std::size_t i = 1; i < output.size(); ++i) {
      const PointXYZ& curr = output[i];
      dx += std::abs(curr.x - prev.x);
      dy += std::abs(curr.y - prev.y);
      dz += std::abs(curr.z - prev.z);
      prev = curr;
    }
    double avgDx = dx / output.size();
    double avgDy = dy / output.size();
    double avgDz = dz / output.size();
    std::cout << "voxelSize = " << voxelSize << ", downsampled size: " << output.size() << ", averate dx: " << avgDx
              << ", average dy: " << avgDy << ", average dz: " << avgDz<< std::endl;
  }
}

You are using a straight line with regular spacing between points, and voxel sizes that are a multiple of this spacing. That means that some points are exactly on the border of two voxels, and there may be voxels with only one point, leading to strange results. This is related to Aliasing. If you use slightly different voxel sizes (no exact multiples of the point-to-point spacing), the results are what you would expect. Alternatively, you can use a random cloud, which also gives you a monotonically decreasing downsampling size.
So I don't see any reason to assume a bug

[REASY]

@mvieth thank you for the response!

Hm, but the algorithm should be resilient to the data issues, no? And it should be able to properly handle edge cases?

If I change the way how the index is computed to compute the delta between X/Y/Z and minX/Y/Z and then multiply it by inverse_leaf_size, then the result is monotonically decreasing. Here is the patch:

Index: filters/include/pcl/filters/impl/voxel_grid.hpp
IDEA additional info:
Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP
<+>UTF-8
===================================================================
diff --git a/filters/include/pcl/filters/impl/voxel_grid.hpp b/filters/include/pcl/filters/impl/voxel_grid.hpp
--- a/filters/include/pcl/filters/impl/voxel_grid.hpp	(revision 7d52b107223ef8adec990608e8685392b48d79c6)
+++ b/filters/include/pcl/filters/impl/voxel_grid.hpp	(date 1658394056813)
@@ -325,9 +325,9 @@
         if (!isXYZFinite ((*input_)[index]))
           continue;
 
-      int ijk0 = static_cast<int> (std::floor ((*input_)[index].x * inverse_leaf_size_[0]) - static_cast<float> (min_b_[0]));
-      int ijk1 = static_cast<int> (std::floor ((*input_)[index].y * inverse_leaf_size_[1]) - static_cast<float> (min_b_[1]));
-      int ijk2 = static_cast<int> (std::floor ((*input_)[index].z * inverse_leaf_size_[2]) - static_cast<float> (min_b_[2]));
+      int ijk0 = static_cast<int> (std::floor ((*input_)[index].x - min_p[0]) * inverse_leaf_size_[0]);
+      int ijk1 = static_cast<int> (std::floor ((*input_)[index].y - min_p[1]) * inverse_leaf_size_[1]);
+      int ijk2 = static_cast<int> (std::floor ((*input_)[index].z - min_p[2]) * inverse_leaf_size_[2]);
 
       // Compute the centroid leaf index
       int idx = ijk0 * divb_mul_[0] + ijk1 * divb_mul_[1] + ijk2 * divb_mul_[2];

It produces:

voxelSize = 3, downsampled size: 67, averate dx: 1.46342, average dy: 1.46352, average dz: 1.46308
voxelSize = 3.33, downsampled size: 59, averate dx: 1.63643, average dy: 1.63654, average dz: 1.63604
voxelSize = 4, downsampled size: 49, averate dx: 1.96019, average dy: 1.96033, average dz: 1.95973