特别提示：《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

最终结果如下：

大家可以比较一下三种邻域搜索方法之间的差异。

以上。

转自：“阿昆的科研日常”微信公众号

如有侵权，请联系本站删除！