点云数据分割 - 盆和冠层 - 常规

点云数据分割 - 盆和冠层

本代码使用区域生长算法将点云数据分割为两类：盆和冠层。

1. 加载点云数据

clear
close all;
clc

%% **********获取点云数据************
[fileName,pathName]=uigetfile('*.txt','文本文件','E:\我的所有研究文件\2022-03\番茄图像+实测值\扫描仪数据\20220620\CK1');   %选择要进行计算的三维点云数据文件路径

if isempty(fileName) || length(fileName) == 1
    fprintf('未选择txt文件！\n');
    return;
end
tic
data=importdata([pathName,fileName]);   %加载txt文件
pc = pointCloud(data(:, 1:3));
pc2 = pointCloud([pc.Location(:, 1), pc.Location(:, 2), pc.Location(:, 3)]);
figure;
% figure;       % 创建新的图形窗口
% 将点云z轴反转
set(gcf,'Color',[1 1 1]);  % 设置图形窗口背景颜色为白色
pcshow(pc2,'MarkerSize',20);   % 显示点云
% xlim([-200 300])
% ylim([200 700])
%  zlim([-1000 -500])
% pcshow(pc,BackgroundColor == [1,1,1],MarkerSize == 20)      %MATLAB2021设置背景颜色，其他版本设置颜色可参阅MATLAB帮助文档
title('原始点云')

%% **********设置相关参数************
min_pts_per_cluster = 150;%较小的会导致更多、更小的簇

max_pts_per_cluster = 25000;%较大的会导致较少、较大的簇

theta_threshold = 3.0/ 180 * pi;%表示相邻点之间的法线夹角阈值，用于区分另一个簇。
%当相邻点之间的法线夹角超过 theta_threshold 时，则认为它们属于不同的簇。
%一般来说，较小的 theta_threshold 值会导致更多的簇形状变化，而较大值会将相邻聚类融合为一个

residual_threshold = 6.5;%表示分类器的残差（即法线方向上的误差）最大值。
%当某个点到属于同一簇的平面法线的距离小于该阈值，则认为该点属于该簇。
%如果阈值较小，则可能会将互不相干的点合并为一簇，而大的则可能会导致簇间的断裂

curvature_threshold = 1.5;%表示当前点的曲率与周围点平均曲率的比值，用于区分平面和曲面。
%当该比值低于曲率阈值时，认为该点是属于平面的。
%较小的曲率阈值会导致曲面的分割数减少，而较大的会增加曲面数量。

neighbour_number =78;%表示每个点的邻居数，即在找到相邻点时的最大数量。
%如果该值较小，则相邻点形成的簇较小；如果该值较大，则相邻点形成的簇较大。

%% **********区域生长算法************
points = double(pc.Location);
% normals = pcnormals(pc,neighbour_number);

% 计算点云法线
normals = pcnormals(pc,neighbour_number);

% 翻转法线方向
%  normals = -normals;

rg = RG();
rg.points = points;
rg.normals = normals;
rg.min_pts_per_cluster = min_pts_per_cluster;
rg.max_pts_per_cluster = max_pts_per_cluster;
rg.theta_threshold = theta_threshold;
rg.residual_threshold = residual_threshold;
rg.curvature_threshold = curvature_threshold;
rg.neighbour_number = neighbour_number;
res = rg.extract();

labels = res.point_labels;
labels(labels==0) = -1;
classes = unique(labels);
num = length(classes);
for i=1:num
    fprintf('第%d类别数量：%d\n',i,length(labels(labels==classes(i))))
%     labels(labels==classes(i)) = i;
end
fprintf('分割出的类别数目为：%d\n',num)

figure;


% set(gcf,'Color',[1 1 1]);  % 设置图形窗口背景颜色为白色

 pc2 = pointCloud([pc.Location(:, 1), pc.Location(:, 2), pc.Location(:, 3)]);


 pcshow(pc2.Location,labels,'MarkerSize',20);
%    view(0,180);
%  camup([0 0 -1]);
% daspect([1 1 1]);
% view(0,90);

% %   xlim([-200 300])
% %  ylim([200 700])
% %  zlim([-1000 -500])
 % pcshow(pc.Location,labels,'MarkerSize',20);   % 显示点云
title('分割点云')
% pcshow(pc.Location,labels,BackgroundColor==[1,1,1],MarkerSize==20)
colormap(hsv(num))
toc

2. 分割点云为盆和冠层

% 设置高度阈值和曲率阈值
height_threshold = 2.5; % 根据您的点云数据进行适当调整
curvature_threshold = 0.1; % 根据您的点云数据进行适当调整

% 提取点云数据中的Z坐标
Z = pc.Location(:, 3);

% 计算高度差
height_diff = Z - min(Z);

% 计算滤波阈值
filter_threshold = height_threshold;

% 计算点云数据的曲率
curvatures = computeCurvature(pc);

% 创建标记矩阵，用于记录每个点的分类结果
labels = zeros(size(pc.Location, 1), 1);

% 标记冠层和盆点
labels(height_diff > filter_threshold & curvatures > curvature_threshold) = 1; % 标记为冠层
labels(height_diff <= filter_threshold | curvatures <= curvature_threshold) = 2; % 标记为盆点

% 可选：可视化结果
figure;
scatter3(pc.Location(labels == 1, 1), pc.Location(labels == 1, 2), pc.Location(labels == 1, 3), 'r', 'filled'); % 以红色显示冠层点云
hold on;
scatter3(pc.Location(labels == 2, 1), pc.Location(labels == 2, 2), pc.Location(labels == 2, 3), 'b', 'filled'); % 以蓝色显示盆点云
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Curved Point Cloud Segmentation');

% 相关函数和区域生长算法的定义与之前的代码相同，请参考之前的代码。

3. 解释

height_threshold 和 curvature_threshold 是关键参数，需要根据您的点云数据进行调整以获得最佳分割结果。
computeCurvature(pc) 函数用于计算点云的曲率，您需要根据实际情况定义该函数。
代码使用 scatter3 函数可视化分割后的点云，您可以根据需要使用其他可视化工具。

4. 注意事项

该示例代码只是一个简单的分割方法，可能无法满足所有情况下的需求。
您可以根据实际情况调整代码，例如使用其他分割方法或修改参数。
建议您使用其他方法，例如基于特征的方法、机器学习方法等，来提高分割精度。

5. 扩展

您可以根据实际需求对分割结果进行调整。例如，您可以将分割后的数据进一步处理，例如计算每个类别点的数量、体积等。