c道格拉斯-普克算法 [1]  (Douglas–Peucker algorithm，亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点，并减少点的数量的一种算法。它的优点是具有平移和旋转不变性，给定曲线与阈值后，抽样结果一定。

算法的基本思路是:对每一条曲线的首末点虚连一条直线,求所有点与直线的距离,并找出最大距离值dmax ,用dmax与限差D相比:若dmax <D,这条曲线上的中间点全部舍去;若dmax ≥D,保留dmax 对应的坐标点,并以该点为界,把曲线分为两部分,对这两部分重复使用该方法。

详细步骤:

(1) 在曲线首尾两点间虚连一条直线,求出其余各点到该直线的距离，如右图（1）。

（2）选其最大者与阈值相比较,若大于阈值,则离该直线距离最大的点保留,否则将直线两端点间各点全部舍去，如右图（2），第4点保留。

（3）依据所保留的点,将已知曲线分成两部分处理,重复第1、2步操作，迭代操作，即仍选距离最大者与阈值比较,依次取舍,直到无点可舍去，最后得到满足给定精度限差的曲线点坐标，如图（3）、（4）依次保留第6点、第7点，舍去其他点,即完成线的化简。 

matlab 代码实现： 

function curve = dp(pnts)
%dp点抽稀算法
clc;
close all;
clear all;
 
x =  1:0.01:  2;
num = size(x,2);
 
pnts = 2*sin(8*pi*x) ; % + rand(1,num)*0.8;
pnts = x .* sin(8*pi*x) + rand(1,num)*0.5;
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 
 
 

c++ 代码实现： 

#include
#include
#include "DouglasPeucker.h"
using namespace std;
void readin(vector &Points,const char * filename)
{
MyPointStruct SinglePoint;
FILE *fp = fopen(filename,"r");
while(fscanf(fp,"%lf%lf",&SinglePoint.X,&SinglePoint.Y)!=EOF)
{
   Points.push_back(SinglePoint);
}
}

void DouglasPeuckerAlgorithm(vector &Points,int

tolerance,const char*filename)
{
DouglasPeucker Instance(Points,tolerance);
Instance.WriteData(filename);
}

void DumpOut1()
{
printf("done!\n");
}

void DumpOut2()
{
printf("need 3 command line parameter:\n[0]executable file name;\n[1]

file name of the input data;\n[2]file name of the output data;\n[3]

threshold.\n");
}

int main(int argc, const char *argv[])
{
if(argc==4)
{
   vector Points;
   readin(Points,argv[1]);

   int threshold = atoi(argv[3]); 
   DouglasPeuckerAlgorithm(Points,threshold,argv[2]);
   DumpOut1();
}
else
{
   DumpOut2();
}
return 0;
}


///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////

////////////

 

#ifndef DOUGLAGPEUCKER

#include
#include
using namespace std;

struct MyPointStruct // 点的结构
{
public: 
double X;
double Y;
double Z;
MyPointStruct()
{
   this->X = 0;
   this->Y = 0;
   this->Z = 0;
};

MyPointStruct(double x, double y, double z) // 点的构造函数
{
   this->X = x;
   this->Y = y;
   this->Z = z;
};
~MyPointStruct(){};
};

class DouglasPeucker
{
public:
vector PointStruct;
vector myTag; // 标记特征点的一个bool数组
vector PointNum;//离散化得到的点号
DouglasPeucker(void){};
DouglasPeucker(vector &Points,int tolerance);
~DouglasPeucker(){};

void WriteData(const char *filename);
private:
void DouglasPeuckerReduction(int firstPoint, int lastPoint, double

tolerance);
double PerpendicularDistance(MyPointStruct &point1, MyPointStruct

&point2, MyPointStruct &point3);
MyPointStruct myConvert(int index);
};

#define DOUGLAGPEUCKER
#endif


///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////

//////////////////

 

#include "DouglasPeucker.h"

double DouglasPeucker::PerpendicularDistance(MyPointStruct &point1,

MyPointStruct &point2, MyPointStruct &point3)
{
// 点到直线的距离公式法
double A, B, C, maxDist = 0;
A = point2.Y - point1.Y;
B = point1.X - point2.X;
C = point2.X * point1.Y - point1.X * point2.Y;
maxDist = fabs((A * point3.X + B * point3.Y + C) / sqrt(A * A + B *

B));
return maxDist;
}

MyPointStruct DouglasPeucker::myConvert(int index)
{
return PointStruct[index];
}

void DouglasPeucker::DouglasPeuckerReduction(int firstPoint, int

lastPoint, double tolerance)
{
double maxDistance = 0;
int indexFarthest = 0; // 记录最大值时点元素在数组中的下标

for (int index = firstPoint; index < lastPoint; index++)
{
   double distance = PerpendicularDistance(myConvert(firstPoint),

myConvert(lastPoint), myConvert(index));

   if (distance > maxDistance)
   {
    maxDistance = distance;
    indexFarthest = index;
   }
}
if (maxDistance > tolerance && indexFarthest != 0) 
{
   myTag[indexFarthest] = true; // 记录特征点的索引信息

   DouglasPeuckerReduction(firstPoint, indexFarthest, tolerance);
   DouglasPeuckerReduction(indexFarthest, lastPoint, tolerance);
}
}

DouglasPeucker::DouglasPeucker(vector &Points,int

tolerance)
{
PointStruct = Points;
int totalPointNum =Points.size();

myTag.resize(totalPointNum,0);

DouglasPeuckerReduction(0, totalPointNum-1, tolerance);

for (int index = 0; index
{
   if(myTag[index])PointNum.push_back(index);
}
}
void DouglasPeucker::WriteData(const char *filename)
{
FILE *fp = fopen(filename,"w");
int pSize = PointNum.size();
for(int index=0;index
{
   fprintf(fp,"%lf\t%lf\n",PointStruct[PointNum[index]].X,PointStruct

[PointNum[index]].Y);
}
}