轨迹数据压缩算法-经典Douglas-Peucker算法 c++实现

时间原因，仅贴代码记录。
数据为一个司机的轨迹数据。

#include <bits/stdc++.h>
using namespace std;

#define pi 3.1415926535897932384626433832795
const double EARTH_RADIUS = 6378137; //地球半径 KM
const int MAXN_V = 1e5 + 5;          //最大轨迹数据点个数
const double MAXN_DIST = 8.0;        //最大距离阈值
int trajec_num = 0,trajec_num_compress = 0;
struct Point
{
    double lon; //经度
    double lat; //纬度
    double time; //时间
} point[MAXN_V], point_compress[MAXN_V];
bool cmp_time(Point A,Point B)
{
    return A.time<B.time;
}
int atoi(string x)
{
    int ret = 0;
    for (int i = 0; i < x.length(); ++i)
        ret = ret * 10 + (x[i] - '0');
    return ret;
}

double atod(string x)
{
    double ret = 0;
    int i;
    for (i = 0; i < x.length(); ++i)
    {
        if (x[i] == '.')
            break;
        ret = ret * 10 + (x[i] - '0');
    }
    double xs = 0.0, ws = 1;
    for (i++; i < x.length(); ++i)
    {
        xs = xs * 10 + x[i] - '0';
        ws *= 10;
    }
    ret += xs / ws;
    return ret;
}

void Read()
{
    freopen("C:\\Users\\14551\\Desktop\\毕设\\数据集\\Dguiji.txt", "r", stdin);
    string trajec; //单条轨迹
    int cnt = 0;
    while (cin >> trajec)
    {
        string ele = "";
        int cnt = 0, pos = 0;
        for (int i = 0; i < trajec.length(); ++i)
        {
            if (trajec[i] == ',')
            {
                if (cnt == 0)
                {
                    pos = atoi(ele);
                    point[pos].time=pos;
                }
                else if (cnt == 1)
                    point[pos].lon = atod(ele);
                ++cnt;
                ele = "";
            }
            else
            {
                ele += trajec[i];
            }
        }
        if (cnt == 2)
            point[pos].lat = atod(ele);
        ++trajec_num;
    }
    freopen("CON","r",stdin);
}

void Print()
{
    freopen("C:\\Users\\14551\\Desktop\\毕设\\数据集\\Dguiji_compress.txt", "w", stdout);
    cout.unsetf(std::ios_base::floatfield);
    for(int i=0;i<trajec_num_compress;++i)
    {
        cout<<point_compress[i].time<<","<<point_compress[i].lon<<","<<point_compress[i].lat<<"\n";
    }
    freopen("CON","w",stdout);
}


inline double rad(double d)
{
    return d * pi / 180.0;
}
double getdist(Point A, Point B) //计算经纬度球面距离
{
    double a = rad(A.lat) - rad(B.lat);
    double b = rad(A.lon) - rad(B.lon);
    double s = 2 * asin(sqrt(pow(sin(a / 2), 2) + cos(rad(A.lat)) * cos(rad(B.lat)) * pow(sin(b / 2), 2)));
    s = s * EARTH_RADIUS;
    return fabs(s);
}

double get_vdist(Point A, Point B, Point X) //海伦公式求垂直距离
{
    double a = getdist(A, B);
    if (a == 0)
        return getdist(A, X);
    double b = getdist(A, X);
    double c = getdist(B, X);
    double p = (a + b + c) / 2;
    double s = sqrt(fabs(p * (p - a) * (p - b) * (p - c)));
    return s * 2 / a;
}

void DP_Compress(int st, int ed, double Dmax)
{
    if(st<ed)
    {
        int index = st+1,keypoint_index = st;
        double max_vdist = 0;
        while(index<ed)
        {
            double vdist = get_vdist(point[st],point[ed],point[index]);
            if(vdist>max_vdist)
            {
                max_vdist=vdist;
                keypoint_index = index ;
            }
            ++index;
        }
        if(max_vdist>=Dmax)
        {
            point_compress[trajec_num_compress++]=point[keypoint_index];
            DP_Compress(st,keypoint_index,Dmax);
            DP_Compress(keypoint_index,ed,Dmax);
        }
    }
}

double get_averErr()
{
    sort(point_compress,point_compress+trajec_num_compress,cmp_time);
    double sumErr = 0.0;
    for(int i=0;i<trajec_num_compress-1;++i)
    {
        int Aid = point_compress[i].time;
        int Bid = point_compress[i+1].time;
        int id = Aid+1;
        while(id<Bid)
        {
            sumErr+=get_vdist(point_compress[i],point_compress[i+1],point[id++]);
        }
    }
    return sumErr/trajec_num;
}

int main()
{
    Read();
    point_compress[trajec_num_compress++]=point[0];
    point_compress[trajec_num_compress++]=point[trajec_num-1];
    DP_Compress(0,trajec_num,MAXN_DIST);
    cout<<"compression rate = "<<trajec_num_compress<<"/"<<trajec_num<<" = "<<(double)1.0*trajec_num_compress/trajec_num<<"\n";
    cout<<"average error =  "<<fixed<<setprecision(10)<<get_averErr()<<"\n";
    Print();
    return 0;
}

运行结果：
compression rate = 189/3160 = 0.0598101
average error =  1.5649999698

效果还是不错的，算法的思想可以算是分治。
最差时间复杂度O(N²)，平均可以达到O(Nlog2N)，最好可以达到O(N)（即只需首尾两点）。这取决于每次keypoint的选取以及自身轨迹特点。
欢迎指正和评论！