离散化曲线.doc

Douglas-Peucker算法曲线离散化
在数字化过程中，需要对曲线进行采样简化，即在曲线上取有限个点，将其变为折线，并且能够在一定程度上保持原有的形状。

经典的Douglas-Peucker算法描述如下：
（1）在曲线首尾两点A，B之间连接一条直线AB，该直线为曲线的弦；
（2）得到曲线上离该直线段距离最大的点C，计算其与AB的距离d；
（3）比较该距离与预先给定的阈值threshold的大小，如果小于threshold，则该直线段作为曲线的近似，该段曲线处理完毕。

（4）如果距离大于阈值，则用C将曲线分为两段AC和BC，并分别对两段取信进行1~3的处理。

（5）当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即可以作为曲线的近似。

Douglas-Peucker算法源代码下载
/*sample.cpp*/
#include <iostream>
#include <vector>
#include "DouglasPeucker.h"
using namespace std;
void readin(vector<MyPointStruct> &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<MyPointStruct> &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<MyPointStruct> Points;
readin(Points,argv[1]);
int threshold = atoi(argv[3]);
DouglasPeuckerAlgorithm(Points,threshold,argv[2]);
DumpOut1();
}
else
{
DumpOut2();
}
return 0;
}
///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ////////////////
/*DouglasPeucker.h*/
#ifndef DOUGLAGPEUCKER
#include <math.h>
#include <vector>
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<MyPointStruct> PointStruct;
vector<bool> myTag; // 标记特征点的一个bool数组
vector<int> PointNum;//离散化得到的点号
DouglasPeucker(void){};
DouglasPeucker(vector<MyPointStruct> &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
///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// //////////////////////
/*DouglasPeucker.cpp*/
#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<MyPointStruct> &Points,int tolerance)
{
PointStruct = Points;
int totalPointNum =Points.size();
myTag.resize(totalPointNum,0);
DouglasPeuckerReduction(0, totalPointNum-1, tolerance);
for (int index = 0; index<totalPointNum; 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<pSize;index++)
{
fprintf(fp,"%lf\t%lf\n",PointStruct[PointNum[index]].X,PointStruct [PointNum[index]].Y);
}
}
///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// /////////////////////。