点云处理【三】(点云降采样)

点云降采样

第一章 点云数据采集
第二章 点云滤波
第二章 点云降采样

---

1. 为什么要降采样？

我们获得的数据量大，特别是几十万个以上的点云，里面有很多冗余数据，会导致处理起来比较耗时。
降采样是一种有效的减少数据、缩减计算量的方法。

2.降采样算法

2.1 随机降采样

根据设置的比例系数随机删除点云，比较接近均匀采样，但不稳定。

Open3d

import numpy as np
import open3d as o3dpcd = o3d.io.read_point_cloud("second_radius_cloud.pcd")
print(pcd)  # 输出点云点的个数
o3d.visualization.draw_geometries([pcd], window_name="原始点云",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)
downpcd = pcd.random_down_sample(sampling_ratio=0.5)
print(downpcd) #降采样后的点云数
o3d.visualization.draw_geometries([downpcd], window_name="随机降采样",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)

PCL

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/filters/random_sample.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("1697165371469.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;pcl::RandomSample<pcl::PointXYZ> random_sampling;random_sampling.setInputCloud(cloud);random_sampling.setSample(10000);  // 设置希望得到的点数random_sampling.filter(*cloud_downsampled);std::cout << "downsampled cloud size: " << cloud_downsampled->width * cloud_downsampled->height << std::endl;pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));viewer->setBackgroundColor(0, 0, 0);  // 设置背景色viewer->addPointCloud<pcl::PointXYZ>(cloud_downsampled, "sample cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud");viewer->initCameraParameters();while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}

2.2 均匀降采样

就是每隔多远采集一个点，

Open3d

import numpy as np
import open3d as o3dpcd = o3d.io.read_point_cloud("second_radius_cloud.pcd")
print(pcd)  # 输出点云点的个数
o3d.visualization.draw_geometries([pcd], window_name="原始点云",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)
downpcd = pcd.uniform_down_sample(6)
print(downpcd) #降采样后的点云数
o3d.visualization.draw_geometries([downpcd], window_name="均匀降采样",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)

PCL

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/filters/uniform_sampling.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("1697165371469.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;pcl::UniformSampling<pcl::PointXYZ> filter;		// 创建均匀采样对象filter.setInputCloud(cloud);					// 设置待采样点云filter.setRadiusSearch(10.0f);					// 设置采样半径filter.filter(*cloud_downsampled);					// 执行均匀采样，结果保存在cloud_filtered中std::cout << "downsampled cloud size: " << cloud_downsampled->width * cloud_downsampled->height << std::endl;pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));viewer->setBackgroundColor(0, 0, 0);  // 设置背景色viewer->addPointCloud<pcl::PointXYZ>(cloud_downsampled, "sample cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud");viewer->initCameraParameters();viewer->saveScreenshot("screenshot.png");while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}

2.3 体素降采样

将空间切割为均匀大小的体素网格，以非空体素的质心代替该体素内的所有点。

原点云位置使用体素降采样后会发生变化。

open3d

import numpy as np
import open3d as o3d
pcd = o3d.io.read_point_cloud("second_radius_cloud.pcd")
print(pcd)  # 输出点云点的个数
o3d.visualization.draw_geometries([pcd], window_name="原始点云",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)downpcd = pcd.voxel_down_sample(voxel_size=5)
print(downpcd)
o3d.visualization.draw_geometries([downpcd], window_name="体素降采样",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)

pcl

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/filters/voxel_grid.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("1697165371469.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;pcl::VoxelGrid<pcl::PointXYZ> sor;sor.setInputCloud(cloud);sor.setLeafSize(10.0f, 10.0f, 10.0f);sor.filter(*cloud_downsampled);std::cout << "downsampled cloud size: " << cloud_downsampled->width * cloud_downsampled->height << std::endl;pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));viewer->setBackgroundColor(0, 0, 0);  // 设置背景色viewer->addPointCloud<pcl::PointXYZ>(cloud_sampled, "sample cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud");viewer->initCameraParameters();viewer->saveScreenshot("screenshot.png");while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}

2.4 最远点降采样

首先随机选择一个点，其次，在剩下点中寻找最远的点，再去再剩下点中找到同时离这两个点最远的点，直到满足采样点个数。
Open3d

import numpy as np
import open3d as o3dpcd = o3d.io.read_point_cloud("second_radius_cloud.pcd")
print(pcd)  # 输出点云点的个数
o3d.visualization.draw_geometries([pcd], window_name="原始点云",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)
downpcd=pcd.farthest_point_down_sample(10000)
print(downpcd) #降采样后的点云数
o3d.visualization.draw_geometries([downpcd], window_name="最远点降采样",width=1024, height=768,left=50, top=50,mesh_show_back_face=True)

PCL

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/common/distances.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("1697165371469.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;size_t N = cloud->size();assert(N >= 10000);srand(time(0));size_t seed_index = rand() % N;pcl::PointXYZ p = cloud->points[seed_index];;cloud_downsampled->push_back(p);cloud->erase(cloud->begin() + seed_index);for (size_t i = 1; i < 10000; i++){float max_distance = 0;size_t max_index = 0;for (size_t j = 0; j < cloud->size(); j++){float distance = pcl::euclideanDistance(p, cloud->points[j]);if (distance > max_distance){max_distance = distance;max_index = max_index;}}p = cloud->points[max_index];cloud_downsampled->push_back(p);cloud->erase(cloud->begin() + max_index);}std::cout << "downsampled cloud size: " << cloud_downsampled->width * cloud_downsampled->height << std::endl;pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));viewer->setBackgroundColor(0, 0, 0);  // 设置背景色viewer->addPointCloud<pcl::PointXYZ>(cloud_downsampled, "sample cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud");viewer->initCameraParameters();while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}

2.5 移动最小二乘法降采样

在MLS法中，需要在一组不同位置的节点附近建立拟合曲线，每个节点都有自己的一组系数用于定义该位置附近拟合曲线的形态。因此，在计算某个节点附近的拟合曲线时，只需要计算该点的该组系数值即可。

此外，每个节点的系数取值只考虑其临近采样点，且距离节点越近的采样点贡献越大，对于未置较远的点则不予考虑。

许多文章都将移动最小二乘法作为降采样方法，我觉得这只是一种平滑，所以这里给了重建代码，不进一步实验了。

PCL

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/filters/voxel_grid.h>
#include <pcl/surface/mls.h>
#include <pcl/search/kdtree.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("second_radius_cloud.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);// 输出的PointCloud中有PointNormal类型，用来存储MLS算出的法线pcl::PointCloud<pcl::PointNormal> mls_points;// 定义MovingLeastSquares对象并设置参数pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointNormal> mls;mls.setComputeNormals(true);mls.setInputCloud(cloud);mls.setSearchMethod(tree);mls.setSearchRadius(30);// 曲面重建mls.process(mls_points);//std::cout << "downsampled cloud size: " << mls_points->width * mls_points->height << std::endl;// 使用PCLVisualizer进行可视化boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("MLS Cloud Viewer"));viewer->addPointCloud<pcl::PointNormal>(mls_points.makeShared(), "MLS Cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "MLS Cloud");viewer->addPointCloudNormals<pcl::PointNormal>(mls_points.makeShared(), 1, 0.05, "normals");  // 可选：显示法线viewer->saveScreenshot("screenshot.png");while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}

2.6 法线空间采样

通过在法向量空间内均匀随机抽样，使所选点之间的法线分布尽可能大，结果表现为地物特征变化大的地方剩余点较多，变化小的地方剩余点稀少，可有效保持地物特征。

Open3d

import open3d as o3d
import numpy as npdef normal_space_sampling(pcd, num_bins=5, num_samples=10000):# 1. 估算法线pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=10, max_nn=30))normals = np.asarray(pcd.normals)# 2. 使用法线的x、y和z分量将法线映射到一个3D直方图或“bin”空间bins = np.linspace(-1, 1, num_bins)normal_bins = np.digitize(normals, bins)unique_bins = np.unique(normal_bins, axis=0)sampled_indices = []for b in unique_bins:indices = np.all(normal_bins == b, axis=1)bin_points = np.where(indices)[0]if bin_points.size > 0:sampled_indices.append(np.random.choice(bin_points))# 如果采样点数不足，从原点云中随机选择其他点while len(sampled_indices) < num_samples:sampled_indices.append(np.random.randint(0, len(pcd.points)))# 3. 从每个bin中选择一个点进行采样sampled_points = np.asarray(pcd.points)[sampled_indices]sampled_pcd = o3d.geometry.PointCloud()sampled_pcd.points = o3d.utility.Vector3dVector(sampled_points)return sampled_pcd# 读取点云
pcd = o3d.io.read_point_cloud("second_radius_cloud.pcd")
sampled_pcd = normal_space_sampling(pcd)
o3d.visualization.draw_geometries([sampled_pcd])

PCL

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/filters/normal_space.h>
#include <pcl/features/normal_3d.h>int main(int argc, char** argv) {pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_downsampled(new pcl::PointCloud<pcl::PointXYZ>);if (pcl::io::loadPCDFile<pcl::PointXYZ> ("second_radius_cloud.pcd", *cloud) == -1){PCL_ERROR("couldn't read file");return 0;}std::cout << "Loaded " << cloud->width * cloud->height<< " data points" << std::endl;// 计算法线pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;ne.setInputCloud(cloud);pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());ne.setSearchMethod(tree);pcl::PointCloud<pcl::Normal>::Ptr cloud_normals(new pcl::PointCloud<pcl::Normal>);ne.setRadiusSearch(30);  // 设置法线估计的半径ne.compute(*cloud_normals);// 法线空间采样pcl::NormalSpaceSampling<pcl::PointXYZ, pcl::Normal> nss;nss.setInputCloud(cloud);nss.setNormals(cloud_normals);nss.setBins(5, 5, 5);  // 设置法线空间的bin数量nss.setSample(cloud->size() / 10);  // 例如，取原始点云大小的1/10nss.filter(*cloud_downsampled);std::cout << "downsampled cloud size: " << cloud_downsampled->width * cloud_downsampled->height << std::endl;pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));viewer->setBackgroundColor(0, 0, 0);  // 设置背景色viewer->addPointCloud<pcl::PointXYZ>(cloud_downsampled, "sample cloud");viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud");viewer->initCameraParameters();while (!viewer->wasStopped()){viewer->spinOnce(100);}return 0;
}