一种基于dem的山顶点的提取算法

0. 算法原理:

后续会持续更新, 先上效果图,

这个是局部的一个效果.

1. 计算45,145,90,180 度方向上的梯度

% File Discription:
% 45°和135° 90° 180° 山顶点提取；
% Author:loveinfall
% CreateTime:2021.5.20（What a good day!(*^__^*) ）clear
clc
close all
path = '../data/等高线/5.tif'img = imread(path);figure,imshow(img)figure(1)
[y,x] = size(img);                 % 取出图像大小
[X,Y] = meshgrid(1:x,1:y);         % 生成网格坐标
pp = double(img);                  % uint8 转换为 double 
mesh(X, Y, pp);                    % 画图
% colormap gray;                   % 选为灰度figure(2)
[C,h]= contour(X, Y, pp);
clabel(C,h)aa = imregionalmax(double(img));Threshold = 0;grayPic=im2double(img);
figure(3)
imshow(grayPic)
title('原图');
%使用构造45度角Sobel算子滤波器
a45=[-2 -1  0;-1  0  1;0  1  2];
a_SFST45=imfilter(grayPic,a45,'replicate');%功能：对任意类型数组或多维图像进行滤波。
SFST45=a_SFST45>=Threshold;
figure(4);
imshow(SFST45)
title('45度角图像边缘检测');%使用构造-45度角Sobel算子滤波器
b45=[0 -1 -2;1  0 -1;2  1  0];
b_SFST135=imfilter(grayPic,b45,'replicate');%功能：对任意类型数组或多维图像进行滤波。
SFST135=b_SFST135>=Threshold;
figure(5);
imshow(SFST135)
title('135度角图像边缘检测');%使用构造90度角Sobel算子滤波器
a90=[-1  0  1;-2  0  2;-1  0  1];
a_SFST90=imfilter(grayPic,a90,'replicate');%功能：对任意类型数组或多维图像进行滤波。
SFST90=a_SFST90>=Threshold;
figure(6);
imshow(SFST90)
title('90度角图像边缘检测');%使用构造180度角Sobel算子滤波器
b180=[1  2  1;0  0  0;-1 -2 -1];
b_SFST180=imfilter(grayPic,b180,'replicate');%功能：对任意类型数组或多维图像进行滤波。
SFST180=b_SFST180>=Threshold;
figure(7);
imshow(SFST180)
title('180度角图像边缘检测');

2. 对山顶点进行粗提取

%% 初步提取山顶点
image_height = size(img,1);
image_width = size(img,2);biase = 1;
peak_index = {};
peak_count_num_flag = 0;peak_destation.x = 0;
peak_destation.y = 0;
peak_destation.dem = 0;
peak_destation.useflag = 0;for index_x = 1+biase:1:image_height-biasefor index_y = 1+biase:1:image_width-biasenum_flag = 0;if a_SFST45(index_x+1,index_y-1) * a_SFST45(index_x-1,index_y+1)<0num_flag = num_flag + 1;endif b_SFST135(index_x-1,index_y-1) * b_SFST135(index_x+1,index_y+1)<0num_flag = num_flag + 1;endif a_SFST90(index_x,index_y-1) * a_SFST90(index_x,index_y+1)<0num_flag = num_flag + 1;endif b_SFST180(index_x-1,index_y) * b_SFST180(index_x+1,index_y)<0num_flag = num_flag + 1;endif num_flag==4 && pp(index_x,index_y)>4200 % 设置 极大值的容差peak_count_num_flag = peak_count_num_flag + 1;peak_index{peak_count_num_flag,1} = [index_x,index_y];peak_index{peak_count_num_flag,2} = pp(index_x,index_y);peak_destation(peak_count_num_flag).x = index_x;peak_destation(peak_count_num_flag).y = index_y;peak_destation(peak_count_num_flag).dem = pp(index_x,index_y);peak_destation(peak_count_num_flag).useflag = 0;endendendfigure()
mesh(X, Y, pp); for index  = 1:size((peak_index),1)hold onyy = peak_index{index,1}(1);xx = peak_index{index,1}(2);dem = peak_index{index,2};plot3(xx,yy,dem,'r+');
endfigure()
[C,h]= contour(X, Y, pp);
for index  = 1:size((peak_index),1)hold onyy = peak_index{index,1}(1);xx = peak_index{index,1}(2);dem = peak_index{index,2};plot3(xx,yy,dem,'r+');
end

3. 山顶点提取结构后处理

%% 进行后处理
% useflag 
% 0:未被使用过
% 1:添加到抽象线中
% 2:已经使用未被添加到 任何线中
other_num = 0;c_line_num = 0; % 统计抽象线的 条数
line_cell = {}; % 用于存放抽象线s_line_num = 0;
s_pot ={};
for index = 1:size((peak_destation),2)pot_num = 0; %记录 抽象线的 长度 line = {};for index_ = 1:size((peak_destation),2)if index_ == indexcontinueendif peak_destation(index_).useflag == 1continueenddistance = sqrt((peak_destation(index).x-peak_destation(index_).x)^2+ (peak_destation(index).y-peak_destation(index_).y)^2);if distance<= 5 % or    根10=3.1623  or  3*根2=4.2426if peak_destation(index_).useflag == 0peak_destation(index_).useflag =1;pot_num = pot_num + 1;line{pot_num,1} = peak_destation(index_); endendendif isempty(line) && peak_destation(index).useflag==0   % 孤立的点peak_destation(index).useflag =2;s_line_num = s_line_num + 1;s_pot{s_line_num,1} = peak_destation(index);elseif isempty(line)==0line{pot_num+1,1} = peak_destation(index);peak_destation(index).useflag =1;c_line_num = c_line_num + 1;line_cell{c_line_num,1} = line;elseif ~peak_destation(index).useflag ==1other_num =  other_num + 1;end     end%% 显示孤立 山顶figure()
mesh(X, Y, pp); for index  = 1:size((s_pot),1)hold onyy =  s_pot{index,1}.x;xx =  s_pot{index,1}.y;dem = s_pot{index,1}.dem;plot3(xx,yy,dem,'r+');
end%% 融合抽象线 山顶
merge_pot.x = 0;
merge_pot.y = 0;
merge_pot.dem=0;
merge_pot.useflag=0;
for index  = 1:size((line_cell),1)L = size((line_cell{index,1}),1);mat = line_cell{index,1};for i=1:L-1for j=1:L-iif mat{j}.dem<mat{j+1}.demtemp=mat{j};mat{j}=mat{j+1};mat{j+1}=temp;endendendmerge_pot(index).x = mat{1}.x;merge_pot(index).y = mat{1}.y;merge_pot(index).dem=mat{1}.dem;merge_pot(index).useflag=0;
endfor index  = 1:size((merge_pot),2)hold onyy =  merge_pot(index).x;xx =  merge_pot(index).y;dem = merge_pot(index).dem;plot3(xx,yy,dem,'b+');
end%% 在等高线上显示
figure()
[C,h]= contour(X, Y, pp);
for index  = 1:size((merge_pot),2)hold onyy =  merge_pot(index).x;xx =  merge_pot(index).y;dem = merge_pot(index).dem;plot3(xx,yy,dem,'b+');text(xx,yy,dem,num2str(index))
end