function nPntSet=Douglas(pntSet,TH)
% @author : slandarer
% pntSet  : 二维数据点
% TH      : 距离阈值

% 向量运算:计算所有点到首位两点连线距离
vertV=[pntSet(end,2)-pntSet(1,2),-pntSet(end,1)+pntSet(1,1)];
baseL=abs(sum((pntSet-pntSet(1,:)).*vertV./norm(vertV),2));

if max(baseL)<TH 
    % 若距离小于阈值则返回首尾点
    nPntSet=[pntSet(1,:);pntSet(end,:)]; 
else
    % 若距离大于阈值则左右两分支分别计算后拼接
    maxPos=find(baseL==max(baseL),1);
    L_PntSet=Douglas(pntSet(1:maxPos,:),TH);
    R_PntSet=Douglas(pntSet(maxPos:end,:),TH);
    nPntSet=[L_PntSet;R_PntSet(2:end,:)];
end
end

这是道格拉斯普克抽稀算法

% 构造一组数据
x=[50 120 200 300 350 400 450 600 660];
y=[100,150,170,200,80,120,195,215,175];
X=x'
Y=y'
pntSet=[X,Y];

% 阈值为50的dp算法
nPntSet=Douglas(pntSet,50);

% 坐标区域修饰
hold on
set(gca,'ytick',0:10:250)

% 绘制原始数据曲线
plot(pntSet(:,1),pntSet(:,2),'Color',[0 0.7 0.3],'LineWidth',2,'Marker','*');
% 绘制新数据曲线
plot(nPntSet(:,1),nPntSet(:,2),'Color',[0.3 0 0.7],'LineWidth',2,'Marker','s');

legend('原始点位置','使用道格拉斯克算法后')

这是运行文件,在设置点和阈值之后直接执行即可

GIS道格拉斯普克算法-LMLPHP

这是执行后所画出来的图片

sum1=0;
sum2=0;
a(1)=50;b(1)=100;
a(2)=120;b(2)=150;
a(3)=200;b(3)=170;
a(4)=300;b(4)=200;
a(5)=350;b(5)=80;
a(6)=400;b(6)=120;
a(7)=450;b(7)=195;
a(8)=600;b(8)=215;
a(9)=660;b(9)=175;

for i=1:8
    c(i)=sqrt((a(i+1)-a(i))^2+(b(i+1)-b(i))^2);
    sum1=sum1+c(i);
end
sum1
d(1)=50;e(1)=100;
d(2)=300;e(2)=200;
d(3)=350;e(3)=80;
d(4)=450;e(4)=195;
d(5)=660;e(5)=175;
for i=1:4
    f(i)=sqrt((d(i+1)-d(i))^2+(e(i+1)-e(i))^2);
    sum2=sum2+f(i);
end
sum2
sum2/sum1

如果要改数据的话需要一个一个改,或者直接读取文件来得更快

GIS道格拉斯普克算法-LMLPHP

这是压缩比

04-09 02:39