特别提示：《Matlab点云处理及可视化》系列文章旨在为初入点云处理领域的朋友提供一份较为权威、可直接用于论文写作的Matlab点云数据处理教程。所有代码均未做加密处理，部分进阶内容会有一定获取门槛，介意请直接划走。
先来看一下最终可视化结果：
1 概述
点云数据具有离散无规则性、盲目性、海量性，因此有必要建立有效的点云数据索引，一方面可提高海量点云数据的处理效率，另一方面，通过结合不同的查询方式（邻域搜索方法）和计算方法，可以为后续去噪、滤波、分类等多种数据处理提供多元化的点云特征（如法向量、曲率、坡度、类别等）。
K-D树（K-Dimension tree）是一种常用的点云数据索引方法。K-D树是一种特殊的二元空间分割树（Binary Space Partitioning trees, BSP trees），其所用分割平面需要与某一坐标轴垂直，而BSP树所用分割平面可以是任意方向的。
点云数据的邻域搜索方法主要包括固定距离（Fixed Distance Neighbors, FDN）与K最邻近（K-Nearest Neighbor, KNN）两种类型。
固定距离搜索首先设定距离阈值，通过计算当前点与待定点之间的距离，将距离小于距离阈值的所有待定点视为当前点的FDN点。常用的固定距离邻域搜索方式包括柱状邻域与球状邻域。其中，柱状邻域通过计算当前点与待定点之间的平面距离，将距离小于预设半径的所有待定点视为当前点的柱状邻域点；而球状邻域通过计算当前点与待定点之间的三维空间距离，将距离小于预设半径的所有待定点视为当前点的球状邻域点。
K最邻近搜索首先确定邻域点个数K，然后计算当前点与待定点之间的距离，并将待定点按距离大小进行排序，将距离最小的K个待定点视为当前点的KNN点。
2 代码实现
function [idx,dist] = nbselect(data,part,varargin)
% 功能：构建K-D树，选择不同邻域类型，返回邻域点索引与距离
% 输入：data - 原始数据(M*3)
% part - 待检索数据(N*3)
% varargin - 球状、柱状邻域或KNN+半径或个数
% 'Q' - 球状邻域
% 'Z' - 柱状邻域
% 'K' - KNN
% 输出：idx - 邻域点索引
% dist - 距离
% example：[sph,dist_sph] = nbselect(data,part,'Q',r_Q);
if varargin{1} == 'Q'
r_Q = varargin{2};
[idx,dist] = rangesearch(data(:,1:3),part(:,1:3),r_Q,'Distance','euclidean','NSMethod','kdtree');
elseif varargin{1} == 'Z'
r_Z = varargin{2};
[idx,dist] = rangesearch(data(:,1:2),part(:,1:2),r_Z,'Distance','euclidean','NSMethod','kdtree');
elseif varargin{1} == 'K'
k = varargin{2};
[idx,dist] = knnsearch(data(:,1:3),part(:,1:3),'Distance','euclidean','NSMethod','kdtree','K',k);
end
end
3 可视化验证
为了检测三种邻域点的搜索效果，采用ISPRS提供的LiDAR点云数据（房屋屋顶）进行验证：
%% 数据准备
% 读取数据
data = load('Vaihingen_building1.txt');
% 参数设置
k = 15; % KNN邻近点个数
r_Z = 3; % 柱状邻域半径
r_Q = 3; % 球状邻域半径
% 待搜索点确定
% 在图上用游标选点，按任意键结束，得到游标点的坐标pos
% figure_handle = figure;
% scatter3(data(:,1),data(:,2),data(:,3),15,'k','filled');
% pos = GetPOS(figure_handle,1);
% close
pos = [497336.899993900 5419701.40999997 282];
%% 邻域点搜索
% 以pos为中心，确定三种邻域内邻域点索引
[Q,dist_Q] = nbselect(data,pos,'Q',r_Q);
[Z,dist_Z] = nbselect(data,pos,'Z',r_Z);
[knn,dist_knn] = nbselect(data,pos,'K',k);
% 确定三种邻域下邻域点坐标
Pz = data(Z{1,1},1:3);
Pq = data(Q{1,1},1:3);
Pk = data(knn,1:3);
%% 可视化验证
% 窗口尺寸设置（单位：厘米）
figureUnits = 'centimeters';
figureWidth = 18;
figureHeight = 20;
figureHandle = figure;
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);
hold on
% 散点图绘制
l1 = scatter3(data(:,1),data(:,2),data(:,3),15);% 原始点云
l2 = scatter3(Pz(:,1),Pz(:,2),Pz(:,3),30);% 柱状邻域点
l3 = scatter3(Pq(:,1),Pq(:,2),Pq(:,3),20);% 球状邻域点
l4 = scatter3(Pk(:,1),Pk(:,2),Pk(:,3),15);% KNN点
l5 = scatter3(pos(:,1),pos(:,2),pos(:,3),20);% 待搜索点
hTitle = title('Neighborhood points selection');
hXLabel = xlabel('XAxis');
hYLabel = ylabel('YAxis');
hZLabel = zlabel('ZAxis');
% 细节优化
axis equal tight
set(l1,'Marker','o','MarkerFaceColor',C1,'MarkerEdgeColor',C1)
set(l2,'Marker','o','MarkerFaceColor',C2,'MarkerEdgeColor',C2)
set(l3,'Marker','^','MarkerFaceColor',C3,'MarkerEdgeColor',C3)
set(l4,'Marker','+','MarkerFaceColor',C4,'MarkerEdgeColor',C4)
set(l5,'Marker','s','MarkerFaceColor',C5,'MarkerEdgeColor',C5)
view(-27.5,46.9);% 视角
set(gca, 'Box', 'on', ... % 边框
'XGrid','on','YGrid','on','ZGrid','on', ... % 网格
'TickDir', 'out', 'TickLength', [0.01 0.01], ... % 刻度
'XMinorTick', 'off', 'YMinorTick', 'off', ... % 小刻度
'XColor', [.1 .1 .1], 'YColor', [.1 .1 .1]) % 坐标轴颜色
hLegend = legend([l1,l2,l3,l4,l5], ...
'原始点云', '柱状邻域', '球状邻域', 'KNN','当前点');
P = hLegend.Position;
hLegend.Position = P + [-0.05 -0.05 0 0];
% 字体和字号
set(gca, 'FontName', 'Arial', 'FontSize', 10)
set([hXLabel, hYLabel, hZLabel], 'FontName', 'Arial','FontSize', 11)
set(hLegend, 'FontName', '微软雅黑', 'FontSize', 11)
set(hTitle, 'FontName', 'Arial', 'FontSize', 12, 'FontWeight' , 'bold')
% 背景颜色
set(gcf,'Color',[1 1 1])
%
输出
print('test.png','-r300','-dpng')
其中，为了区分不同邻域点对象，从Matlab配色神器TheColor的XKCD颜色库中选择5种对比色：
% 颜色定义
C = TheColor('xkcd',[914 207 520 914 546]);
C1 = C(1,1:3);
C2 = C(2,1:3);
C3 = C(3,1:3);
C4 = C(4,1:3);
C5 = C(5,1:3);
(点击上图查看TheColor具体功能)
获取方式：公众号（阿昆的科研日常）后台回复 TC
最终结果如下：
大家可以比较一下三种邻域搜索方法之间的差异。
以上。
转自：“阿昆的科研日常”微信公众号
如有侵权，请联系本站删除！
