Douglas-Peucker算法的Matlab实现

代码。

%% Douglas-Peucker 道格拉斯-普克算法

function curve = dpEdgeContour(pnts)
clc;
clear;
 
A = readmatrix('EdgeContour1.xls','Sheet','1');%输入数据
x = A(:,1)'; %输入的点横坐标
y = A(:,2)'; %输入的点纵坐标
             %也可以输入一个函数(y = f(x)还有其定义域，即x范围)

scatter(x,y,'.');
num = size(x,2); %返回x的列数

 

pnts = y;

figure; plot( x, pnts,'-.'); title('pnts');
 
head = [x(1), pnts(1)] ;
tail = [x(end), pnts(end)] ;
  
%距离阈值
dist_th = 0.1; 
 
hold on;
max_dist = 0;
max_dist_i = 0;
%把起始点和终止点加进去
curve = [head; tail];
pnt_index = [1, num];
 
while(1)         
    %目前抽稀后的曲线上的点数
    curve_num = size(curve,1);
   
    %标记是否添加了新点，若有，表明还要继续处理
    add_new_pnt = 0; 
    %  对区间头尾连线，然后计算各点到这条直线的距离 
    for nx = 1:curve_num-1
        cur_pnt = curve(nx,:);
        %下一个抽稀了的曲线上的点
        next_pnt = curve(nx+1,:);
        pnt_d = next_pnt - cur_pnt ;
         
        th = atan2(pnt_d(2), pnt_d(1));
        angle = th*180/pi;
        %直线方程
        % y = kx + b
        k = tan(th);
        b = cur_pnt(2) - k*cur_pnt(1);
        k2 = k*k;
        deno = sqrt(1 + k *k) ;
 
         max_dist = 0;
         pnt_index(nx);
         pnt_index(nx+1);
         
        %对这一区间的点计算到直线的距离
        for i= pnt_index(nx) : pnt_index(nx+1)
            dist = abs(pnts(i) - k*x(i) - b)/deno ;
 
            if(dist> max_dist)      
                max_dist = dist;
                max_dist_i = i;        
            end 
        end 
        max_dist;
        max_dist_i; 
 
        far_pnt = [x(max_dist_i), pnts(max_dist_i)];
 
        %最远的点加进去
        if(max_dist > dist_th)
             curve = [curve(1:nx,:); far_pnt; curve(nx+1:end,:)]; 
             pnt_index = [pnt_index(1:nx), max_dist_i, pnt_index(nx+1:end)];
             %标记添加了新点，可能还要继续处理
             add_new_pnt = 1;        
        end 
    end 
    close all ;
    figure; plot( x, pnts,'-.'); title('pnts');
    hold on;
    plot(curve(:,1), curve(:,2), '-g*');
    drawnow;
     
     %如果各点到直线距离都没有超过阈值则退出
     %处理完毕了，ok了，退出
     if(0 == add_new_pnt)
        break;
     end      
end 

% ans为最后取的点