在之前的文章中，分享了Matlab基于KD树的邻域点搜索方法：
在此基础上，进一步分享一下基于KD树的离散点密度特征提取方法。
先来看一下成品效果：
特别提示：《Matlab点云处理及可视化》系列文章整理自作者博士期间的部分成果，旨在为初入点云处理领域的朋友提供一份较为权威、可直接用于论文写作的Matlab点云数据处理教程。如果觉得有用可以分享给更多的人。
1 概述
点云密度特征一般用单位面积/体积内的离散点数量表示。
其中，二维(2-Dimension, 2D)平面密度可以通过柱状邻域提取；三维(3-Dimension, 3D)体密度可以通过球邻域提取。
2 代码实现
function density_2D = density2D_KD(data,radius)
% 功能：利用KD树提取离散点2D密度特征
% 输入：data - 原始数据(m*3)
% 输出：planes - 拟合所得平面参数
M = size(data,1);
density_2D = zeros(M,1);
idx = rangesearch(data(:,1:2),data(:,1:2),radius,'Distance','euclidean','NSMethod','kdtree');
for i = 1:M
density_2D(i,1) = length(idx{i})/(pi*radius^2);
end
end
function density_3D = density3D_KD(data,radius)
% 功能：利用KD树提取离散点3D密度特征
% 输入：data - 原始数据(m*3)
% 输出：planes - 拟合所得平面参数
M = size(data,1);
density_3D = zeros(M,1);
idx = rangesearch(data(:,1:3),data(:,1:3),radius,'Distance','euclidean','NSMethod','kdtree');
for i = 1:M
density_3D(i,1) = length(idx{i})/(4/3*pi*radius^3);
end
end
3 可视化验证
为了检测密度特征提取的效果，采用某地机载LiDAR实测点云数据进行验证：
%% 数据准备
% 读取数据
data = load('data1.txt');
%% 点云密度特征提取
% 定义半径
radius = 1;
% 2D平面密度
density_2D = density2D_KD(data(:,1:2),radius);
% 3D体密度
density_3D = density3D_KD(data(:,1:3),radius);
%% 原始点云可视化
% 窗口尺寸设置（单位：厘米）
figureUnits = 'centimeters';
figureWidth = 2*10;
figureHeight = 1*10;
figureHandle1 = figure;
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);
t1 = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
nexttile
scatter(data(:,1),data(:,2),7,data(:,3),'filled')
hTitle1 = title('原始点云(俯视)');
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca,'ztick',[])
axis equal off tight
colormap(map)
colorbar
nexttile
scatter3(data(:,1),data(:,2),data(:,3),5,data(:,3),'filled')
hTitle2 = title('原始点云(侧视)');
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca,'ztick',[])
axis equal off tight
view(-1.157471311476228e+02,53.408837245676686)
colorbar
% 细节优化
set([hTitle1,hTitle2], 'FontName', '微软雅黑', 'FontSize', 12, 'FontWeight' , 'bold')
% 背景颜色
set(gcf,'Color',[1 1 1])
%
输出
figW = figureWidth;
figH = figureHeight;
set(figureHandle1,'PaperUnits',figureUnits);
set(figureHandle1,'PaperPosition',[0 0 figW figH]);
fileout = 'test0';
print(figureHandle1,[fileout,'.png'],'-r300','-dpng');
%% 平面、体密度提取结果可视化验证
% 窗口尺寸设置（单位：厘米）
figureUnits = 'centimeters';
figureWidth = 2.35*10;
figureHeight = 1*10;
figureHandle2 = figure;
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);
t2 = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
nexttile
scatter(data(:,1),data(:,2),15,density_2D(:,1),'filled')
hTitle1 = title('基于KD树的2D平面密度特征可视化');
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca,'ztick',[])
axis equal off tight
colormap(map)
colorbar
nexttile
scatter(data(:,1),data(:,2),15,density_3D(:,1),'filled')
hTitle2 = title('基于KD树的3D体密度特征可视化');
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca,'ztick',[])
axis equal off tight
colorbar
% 细节优化
set([hTitle1,hTitle2], 'FontName', '微软雅黑', 'FontSize', 12, 'FontWeight' , 'bold')
% 背景颜色
set(gcf,'Color',[1 1 1])
%
输出
figW = figureWidth;
figH = figureHeight;
set(figureHandle2,'PaperUnits',figureUnits);
set(figureHandle2,'PaperPosition',[0 0 figW figH]);
fileout = 'test';
print(figureHandle2,[fileout,'.png'],'-r300','-dpng');
其中，为了区分不同对象，从Matlab配色神器TheColor的SCI颜色库中选择渐变色：
%% 颜色定义
map = TheColor('sci',2068);
% map = flipud(map);
(点击上图查看TheColor具体功能)
获取方式：公众号（阿昆的科研日常）后台回复 TC
最终结果如下：
以上。
转自：“阿昆的科研日常”微信公众号
如有侵权，请联系本站删除！
