8月的时候去了野外调查，一直没有时间更新博客，最近回到了学校，写了matlab中实现一种简单的坡度滤波算法，由于本人比较懒，以后的每个博客都尽量减少说明。
此外，其他的滤波：布料滤波(CSF)、简单形态学滤波、经验模态分解滤波的资源链接已经放在了
2021-8-13的博客中了。
https://blog.csdn.net/qq_39632121/article/details/119671068?spm=1001.2014.3001.5501
本文用到的示列数据可以在链接中下载到。
坡度滤波
?原理：坡度滤波算法中，把点云按照一定大小的格网进行划分，每个格网具有一定的坡度阈值，每个格网的点高程值减去格网中最低的高程值，然后除以它与最低点的距离得到坡度值。当坡度值超过这个阈值时，则被认定为非地面点；低于阈值则被认定为地面点。
算法步骤：
①设置格网大小（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) > 1
temp = 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
toc
pcshow(a,gp)示列数据显示如下：
?滤波结果如下：其中地面点是紫色，非地面点黄色高亮显示
?
非地面点显示如下：
地面点显示如下：
看起来不比CSF滤波效果差，运行时间1.38秒，还不错。?
function资源已经上传，可以根据需要下载，在博主的主页。