道格拉斯-普克算法

道格拉斯-普克算法 [1]  (Douglas–Peucker algorithm，亦稱為拉默-道格拉斯-普克算法、迭代適應點算法、分裂與合並算法)是將曲線近似表示為一系列點，並減少點的數量的一種算法。它的優點是具有平移和旋轉不變性，給定曲線與閾值后，抽樣結果一定。

算法步驟

1. 連接曲線首尾兩點A、B形成一條直線AB；
2. 計算曲線上離該直線段距離最大的點C，計算其與AB的距離d；
3. 比較該距離與預先給定的閾值threshold的大小，如果小於threshold，則以該直線作為曲線的近似，該段曲線處理完畢。
4. 如果距離大於閾值，則用點C將曲線分為兩段AC和BC，並分別對兩段曲線進行步驟[1~3]的處理。
5. 當所有曲線都處理完畢后，依次連接各個分割點形成折線，作為原曲線的近似。

//DouglasPeucker.h
#pragma once
#include <vector>
#include <stdio.h>
#include <iostream>

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 MyPointStruct
{
public:
    //MyPointStruct pointXYZ;
    vector<MyPointStruct> PointStruct;
    vector<bool> myTag; // 標記特征點的一個bool數組
    vector<int> PointNum;//離散化得到的點號

public:
    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);
};

//DouglasPeucker.cpp

#include "DouglasPeucker.h"
#include <stdio.h>

DouglasPeucker::DouglasPeucker()
{

}

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);
        }
    }
}


DouglasPeucker::~DouglasPeucker()
{

}

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);
    }
}

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);
    }
}

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];
}