MATLAB实现点云地面点非地面点滤波

MATLAB实现点云地面点非地面点滤波

本次滤波采用了自主设计的均值-方差滤波法，通过方差均值计算高程异常值，然后将异常值选取为非地面点。

Cloudcompare打开数据文件

Matlab打开数据文件，并显示地面点和非地面点

整体数据显示，明显可以看出存在很多非地面点

找随机点及最邻近的50个点

通过直方图的结果，将高程大于300的点的高程用中值代替

使用区域均值替代异常高程值

使用区域最小值替代异常高程值

真实非地面点和均值-方差滤波法提取的非地面点

真实的非地面点

均值-方差滤波法得到的非地面点

通过上述的结果可以看出，大量的非地面点已经被提取到了并且滤掉了，但是还有一些零散的点没有被滤波掉，这个问题的话还得通过参数得调整来解决，后面可能会有比较好的办法解答。
附录：部分代码
%% 地面点滤波
clc,clear,close all
%% 加载路径
addpath(genpath(‘F:\专业课程\三维激光雷达\实验四\test_4’));
pcdata=importdata(’.\samp22_选做题数据.txt’);
%% 显示
ptCloud=[pcdata(:,1)-mean(pcdata(:,1)),pcdata(:,2)-mean(pcdata(:,2)),pcdata(:,3)];
figure
pcshow(ptCloud);
%% 找其中的部分点
index=logical(pcdata(:,4)0);
ptCloud_1=ptCloud(index1,:);%可以随机选
ptCloud_2=ptCloud(index0,:);
figure
subplot(1,2,1)
pcshow(ptCloud_1);
colormap(jet);
title(‘ground point cloud’);
subplot(1,2,2)
pcshow(ptCloud_2);
colormap(jet);
title(‘no ground cloud filtering’);
%% KD检索ptCloud_1中距离最近的10个点
kdtreeobj = KDTreeSearcher(ptCloud,‘distance’,‘euclidean’);
[idx,dist]=knnsearch(ptCloud,ptCloud_1,‘dist’,‘euclidean’,‘k’,50);
% 显示编号为2000的50个近邻点
points=[];
for i=1:50
id=idx(2000,i);
points=[points;ptCloud(id,1),ptCloud(id,2),ptCloud(id,3)];
end
plot3(ptCloud(:,1),ptCloud(:,2),ptCloud(:,3),’.g’);
hold on;
plot3(points(:,1),points(:,2),points(:,3),’.r’);
hold on;
plot3(ptCloud_1(2000,1),ptCloud_1(2000,2),ptCloud_1(2000,3),‘ob’);
%% 按最近的五十个点进行方差均值滤波
new_ptCloud=ptCloud;
new_ptCloud_median=median(new_ptCloud);
threshold=300;% 地面点高程阈值
index=logical(new_ptCloud(:,3)>threshold);
% new_ptCloud(index1,3)=new_ptCloud_median(1,3);
for i=1:size(ptCloud_1,1)
index=idx(i,:);
Elevation=ptCloud(index,3);
Elevation_mean=mean(Elevation);%求均值
Elevation_var=var(Elevation)/20;%求方差
Elevation_min=min(Elevation);%最小值
temp=new_ptCloud(index,3);
for j=1:size(temp,1)
if abs(temp(j,1)-Elevation_mean)>Elevation_var
temp(j,1)=nan;
end
end
new_ptCloud(index,3)=temp;
end
figure
subplot(1,2,1)
pcshow(ptCloud);
title(‘Original point cloud’);
subplot(1,2,2)
pcshow(new_ptCloud);
colormap(jet);
title(‘point cloud filtering’);
%
figure
subplot(1,2,1)
pcshow(ptCloud_2);
colormap(jet);
title(‘really non ground point’);
subplot(1,2,2)
pcshow(new_ptCloud);
colormap(jet);
title(‘filtering non ground point’);