点群データのダウンサンプリング(Visual Studio 2019)

点群データの点群数が多すぎて、MeshLabなどのビューアで開けなかったり、いろいろなフィルタリング処理が重すぎたりして不便なときがあると思います。そのようなときに点群データの点群数を減らすためのダウンサンプリングとして、Point Cloud Library (PCL) で実現出来る方法をひとつご紹介します。
PCLのVoxelGridというクラスを使います。このVoxelGridの説明として下記のような説明がありました。

The VoxelGrid class creates a 3D voxel grid (think about a voxel grid as a set of tiny 3D boxes in space) over the input point cloud data. 

つまり、点群を一定サイズの三次元の箱で区切ったVocelグリッドで表現するためのクラスで、Voxel内のすべての点群はダウンサンプリングなどの方法で近似化され、Voxelの中心値で表現されます。

VoxelGridでダウンサンプリングするサンプルコードを下記に示します。このサンプルでは入力された点群全体の体積を算出し、その体積を100万個のVoxelに区切ります。そして、1つのVoxel(立方体)の1辺の長さを算出して、その値をVoxcelGridのsetLeafSizeに設定することでダウンサンプリングを実現しています。

 

#include  <iostream> 
#include  <vector>
#include <pcl/io/ply_io.h>
#include <pcl/point_cloud.h>
#include <pcl/point_types.h>
#include <pcl/common/transforms.h>
#include <pcl/common/common.h>
#include <pcl/common/io.h>
#include <pcl/filters/voxel_grid.h>
#include <pcl/features/moment_of_inertia_estimation.h>

using PointT = pcl::PointXYZRGB;

int
main()
{
    pcl::PointCloud<PointT>::Ptr cloud(new pcl::PointCloud<PointT>());
    pcl::PointCloud<PointT>::Ptr cloud_filtered(new pcl::PointCloud<PointT>());

    // Read PLY file
    pcl::PLYReader reader;
    if (reader.read<PointT>("jaguar.ply", *cloud) == -1) {
        PCL_ERROR("Couldn't read file Table.ply\n");
        return (-1);
    }

    std::cerr << "PointCloud couts before filtering: " << cloud->width * cloud->height << std::endl;

    // Get bounding box of point cloud
    pcl::MomentOfInertiaEstimation <pcl::PointXYZ> feature_extractor;
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_mono(new pcl::PointCloud<pcl::PointXYZ>());
    pcl::copyPointCloud(*cloud, *cloud_mono);
    feature_extractor.setInputCloud(cloud_mono);
    feature_extractor.compute();

    std::vector <float> moment_of_inertia;
    std::vector <float> eccentricity;
    pcl::PointXYZ min_point_AABB;
    pcl::PointXYZ max_point_AABB;
    pcl::PointXYZ min_point_OBB;
    pcl::PointXYZ max_point_OBB;
    pcl::PointXYZ position_OBB;
    Eigen::Matrix3f rotational_matrix_OBB;
    float major_value, middle_value, minor_value;
    Eigen::Vector3f major_vector, middle_vector, minor_vector;
    Eigen::Vector3f mass_center;

    feature_extractor.getMomentOfInertia(moment_of_inertia);
    feature_extractor.getEccentricity(eccentricity);
    feature_extractor.getAABB(min_point_AABB, max_point_AABB);
    feature_extractor.getOBB(min_point_OBB, max_point_OBB, position_OBB, rotational_matrix_OBB);
    feature_extractor.getEigenValues(major_value, middle_value, minor_value);
    feature_extractor.getEigenVectors(major_vector, middle_vector, minor_vector);
    feature_extractor.getMassCenter(mass_center);

    // Get bounding box volume
    const float len_x = fabs(max_point_AABB.x - min_point_AABB.x);
    const float len_y = fabs(max_point_AABB.y - min_point_AABB.y);
    const float len_z = fabs(max_point_AABB.z - min_point_AABB.z);
    const float model_volume = len_x * len_y * len_z;
    std::cerr << "len_x=" << len_x << " len_y=" << len_y << " len_z=" << len_z << " volume=" << model_volume << std::endl;

    // Compute voxel size assuming point cloud counts as 1,000,000 after downsampling.
    const float target_cout = 1000000;
    float voxel_vol = model_volume / target_cout;
    float voxel_len = pow(voxel_vol, (1.0 / 3.0));
    std::cerr << "voxel volume=" << voxel_vol << std::endl;
    std::cerr << "voxel length=" << voxel_len << std::endl;

    // Create the filtering object
    pcl::VoxelGrid<PointT> sor;
    sor.setInputCloud(cloud);
    sor.setLeafSize(voxel_len, voxel_len, voxel_len);
    sor.filter(*cloud_filtered);
    std::cerr << "PointCloud after filtering: " << cloud_filtered->width * cloud_filtered->height << std::endl;

    // Write as PLY file
    pcl::PLYWriter writer;
    writer.write("jaguar_filtered.ply", *cloud_filtered, false);

    return (0);
}

上記コードをVisual Studio 2019でコンパイルした際のCMakeLists.txtを参考用に記載します。

 

cmake_minimum_required(VERSION 2.8 FATAL_ERROR)

set(ProjectName downsample)
project(${ProjectName})
add_executable (${ProjectName} downsample.cpp)

find_package(Eigen3 REQUIRED)
include_directories(${EIGEN3_INCLUDE_DIR})
set(Eigen_LIB "Eigen3::Eigen") 
target_link_libraries (${ProjectName} ${Eigen_LIB})


find_package(PCL 1.8 REQUIRED)
include_directories(${PCL_INCLUDE_DIRS})
link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})
target_link_libraries (${ProjectName} ${PCL_LIBRARIES})