[surface] Unexpected results from pcl::surface::MarchingCubes

Describe the bug

The pcl::surface::MarchingCubes algorithm gives unexpected results with seemingly random patterns.

Context

I tried to reconstruct a mesh (ply or stl, but not really relevant to reproduce the issue) from a .pcd file.

Expected behavior

I expected a correct approximation of the input point cloud as mesh (at least after playing with the voxel size).

Current Behavior

The output mesh represents the input cloud to a certain degree, but has a lot of seemingly random data not present in the original cloud.

To Reproduce

Use the attached code and this example PCD file to reconstruct a mesh using MarchingCubes.

Screenshots/Code snippets

The Original PCD:

The reconstructed mesh (in MeshLab):

The code:

#include <iostream>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/io/vtk_lib_io.h>
#include <pcl/surface/marching_cubes.h>
#include <pcl/surface/marching_cubes_hoppe.h>
#include <pcl/search/kdtree.h>
#include <pcl/features/normal_3d.h>

#define VOXEL_RES 100


using namespace pcl;

int main (int argc, char** argv)
{
    if (argc < 2)
    {
        std::cerr << "Usage: " << argv[0] << " <pcdfile>" << std::endl;
        return -1;
    }
    // Load point cloud from disk
    PointCloud<PointXYZ>::Ptr cloud (new PointCloud<PointXYZ>);
    if (pcl::io::loadPCDFile<pcl::PointXYZ> (argv[1], *cloud) == -1) //* load the file
    {
        PCL_ERROR ("Couldn't read pcd file.\n");
        return (-1);
    }
    std::cout << "Loaded "
              << cloud->width * cloud->height
              << " data points from \"" << argv[1] << "\"."
              << std::endl;

    // Normal estimation for loaded point cloud
    search::KdTree<PointXYZ>::Ptr tree (new search::KdTree<PointXYZ>);
    tree->setInputCloud(cloud);
    NormalEstimation<PointXYZ, Normal> n;
    n.setInputCloud(cloud);
    n.setSearchMethod(tree);
    n.setKSearch(20);
    PointCloud<Normal>::Ptr normals (new PointCloud<Normal>);
    n.compute(*normals);

    // Concat XYZ and normal values
    PointCloud<PointNormal>::Ptr normalCloud (new PointCloud<PointNormal>);
    concatenateFields(*cloud, *normals, *normalCloud);


    // Create search tree
    search::KdTree<PointNormal>::Ptr tree1 (new search::KdTree<PointNormal>);
    tree1->setInputCloud (normalCloud);
    // Use MarchingCubes to generate a mesh
    PolygonMesh::Ptr output (new PolygonMesh);
    MarchingCubesHoppe<PointNormal> mc;
    PointCloud<PointNormal>::ConstPtr tmpCloudConstPtr (normalCloud);
    mc.setInputCloud (tmpCloudConstPtr);
    mc.setSearchMethod(tree1);
    mc.setGridResolution(VOXEL_RES, VOXEL_RES, VOXEL_RES);
    mc.reconstruct (*output);
    io::savePCDFileASCII ("output.pcd", *normalCloud);


    cout << "Created " << output->polygons.size() << " triangles." << endl;

    // Save the PolygonMesh to an .stl file
    io::savePolygonFileSTL ("output.stl", *output);


    return (0);
}

Environment

• OS: Ubuntu 18.04
• Compiler: g++ 7.5.0
• PCL Version: 1.8, but still persistent in 1.10.1

Reproduce the result with the attached C++ program and box.pcd file, starting at pcl/surface/marching_cubes_hoppe.h and the MarchingCubesHoppe entry point. Inspect how the estimated normals, search method, and grid resolution affect reconstruction, then verify that the generated STL no longer contains the reported random patterns.