简单的点云坡度滤波

坡度滤波"/>

简单的点云坡度滤波

8月的时候去了野外调查，一直没有时间更新博客，最近回到了学校，写了matlab中实现一种简单的坡度滤波算法，由于本人比较懒，以后的每个博客都尽量减少说明。

此外，其他的滤波：布料滤波(CSF)、简单形态学滤波、经验模态分解滤波的资源链接已经放在了

2021-8-13的博客中了。

2021-8-13-点云地面点识别-matlab示列与免费资源共享_~追风筝的猫的博客-CSDN博客

本文用到的示列数据可以在链接中下载到。

坡度滤波

 原理：坡度滤波算法中，把点云按照一定大小的格网进行划分，每个格网具有一定的坡度阈值，每个格网的点高程值减去格网中最低的高程值，然后除以它与最低点的距离得到坡度值。当坡度值超过这个阈值时，则被认定为非地面点；低于阈值则被认定为地面点。

算法步骤：

①设置格网大小（size，s），设置坡度阈值（threshold，t），计算格网的数量。

②顺序循环格网，寻找每个格网的最低点，然后格网内其他的点根据公式计算坡度值，假设最低点为min_p。

③每个格网循环每个点，用坡度值和坡度阈值t做比较，区分地面点和非地面点。

④循环结束，输出结果

MATLAB2019b-坡度滤波实现

clear % 加载点云数据，格式为txt，importdata可以打开
% load也可以打开
a = importdata("sample.txt");
% a为n*3的点云，仅仅包含了点云的xyz值% 设置格网大小，单位：米（m）
s = 5;
% 设置坡度阈值，单位：米（m）
t = 0.5;
% 设置地面点索引
gp = zeros(length(a(:,1)),1);% 计算格网数量
xcount = ceil((max(a(:,1))-min(a(:,1)))/s);
ycount = ceil((max(a(:,2))-min(a(:,2)))/s);tic
% 循环格网，按照先循环x再y的顺序进行
for j = 1:ycount% 格网的最低y值limity1 = (j-1) * s + min(a(:,2));% 格网的最高y值limity2 = j * s + min(a(:,2));% 寻找max_y和min_y之间的点b = a(a(:,2)>=limity1 & a(:,2)<limity2,:);% 索引在范围内的点idx1 = find(a(:,2)>=limity1 & a(:,2)<limity2);for i = 1:xcount% 格网的最低x值limitx1 = (i-1) * s + min(a(:,1));% 格网的最高x值limitx2 = i * s + min(a(:,1));% 寻找格网中的点c = b(b(:,1)>=limitx1 & b(:,1)<limitx2,:);% 索引格网内的点idx2 = find(b(:,1)>=limitx1 & b(:,1)<limitx2);idx3 = idx1(idx2);% 高程最小的点temp = c(c(:,3)==min(c(:,3)),:);if length(temp) > 1temp = temp(1,:);end% 计算坡度值 slopevaule = (c(:,3)-temp(:,3))./ sqrt((c(:,1)-temp(:,1)).^2+(c(:,2)-temp(:,2)).^2);% 索引满足要求的点idx4 = find(slopevaule > t);idx = idx3(idx4);gp(idx,:) = 1;end
end
tocpcshow(a,gp)

示列数据显示如下：

 滤波结果如下：其中地面点是紫色，非地面点黄色高亮显示

非地面点显示如下：

地面点显示如下：

看起来不比CSF滤波效果差，运行时间1.38秒，还不错。 

function资源已经上传，可以根据需要下载，在博主的主页。