Polyline Simplification(多线段简化/离散采样)

本文参考自：https://www.codeproject.com/Articles/114797/Polyline-Simplification

前言

在计算几何中，经常我们会碰到需要简化输入Polyline的场景。多段线（Polyline）简化算法可以帮助我们减少Polyline的点数，从而降低输入规模。对多段线简化算法，通常的做法是在一定的近似精度下，删除一些点或者边。在计算机图形学中，当多段线分辨率高于显示分辨率时，多个顶点和边很可能映射到单个像素上，这意味着计算机将花费一些额外的资源去绘制不可见的顶点。在绘制之前，可以通过多段线简化算法减少点数，从而轻松地避免这种计算资源的浪费。

实践中，根据不同场景的需要，有下面几种简化算法通常比较有用：

临近点简化(Radial Distance)
垂直距离简化(Perpendicular Distance)
Douglas-Peucker算法

临近点简化

该算法是一个简单的O(n)复杂度的暴力算法。对某个顶点，接下来连续的顶点与该点（key）的距离如果都在给定的误差范围内，那么将移除这些点。过程如下图：

我们通常保留Polyline的起始点和终止点作为结果的一部分。首先将起始点标记为Key，然后沿着Polyline前进，与Key相邻的连续顶点到Key的距离如果小于给定误差就移除，遇到的第一个超过误差的点标记为下一个Key；从这个新的key开始重复上面的过程，知道到达之后一个顶点。

垂直距离简化

临近点算法将点-点距离作为误差判据，而垂直距离简化则是将点-线段的距离作为误差判据。对于每个顶点Vi，需要计算它与线段[Vi-1， Vi+1]的垂直距离，距离比给定误差小的点都将被移除。过程如下图所示：

如上图，首先，对前三个顶点进行处理，计算第二个顶点的垂直距离，将这个距离与容差进行比较，大于设定的误差，所以第二个顶点作为Key保留(简化结果的一部分)。然后算法开始处理接下来的三个顶点。这一次，计算的距离低于公差，因此中间顶点被删除。如此重复，直到处理完所有顶点。

Note：对于每个被移除的顶点Vi，下一个可能被移除的候选顶点是Vi+2。这意味着原始多段线的点数最多只能减少50%。为了获得更高的顶点减少率，需要多次遍历。

Douglas-Peucker算法

Douglas-Peucker算法使用点-边距离作为误差衡量标准。该算法从连接原始Polyline的第一个和最后一个顶点的边开始，计算所有中间顶点到边的距离，距离该边最远的顶点，如果其距离大于指定的公差，将被标记为Key并添加到简化结果中。这个过程将对当前简化中的每条边进行递归，直到原始Polyline的所有顶点与当前考察的边的距离都在允许误差范围内。该过程如下图所示:

初始时，简化结果只有一条边（起点-终点）。第一步中，将第四个顶点标记为一个Key，并相应地加入到简化结果中；第二步，处理当前结果中的第一条边，到该边的最大顶点距离低于容差，因此不添加任何点；在第三步中，当前简化的第二个边找到了一个Key（点到边的距离大于容差）。这条边在这个Key处被分割，这个新的Key添加到简化结果中。这个过程一直继续下去，直到找不到新的Key为止。注意，在每个步骤中，只处理当前简化结果的一条边。

这个算法的最坏时间复杂度是O(nm), 平均时间复杂度 O(n log m)，其中m是简化后的Polyline的点数。因此，该算法是Output-sensitive的。当m很小时，该算法会很快。

该算法在实践中有很好的近似效果，且效率很高，下面给出DP算法的实现：


// Assume that classes are already given for the objects:
//    Point and Vector with
//        coordinates {float x, y, z;}     // as many as are needed
//        operators for:
//            == to test equality
//            != to test inequality
//            (Vector)0 = (0,0,0)         (null vector)
//            Point  = Point ± Vector
//            Vector = Point - Point
//            Vector = Vector ± Vector
//            Vector = Scalar * Vector    (scalar product)
//            Vector = Vector * Vector    (cross product)
//    Segment with defining endpoints {Point P0, P1;}
//===================================================================// dot product (3D) which allows vector operations in arguments
#define dot(u,v)  ((u).x * (v).x + (u).y * (v).y + (u).z * (v).z)
#define norm2(v)  dot(v,v)         // norm2 = squared length of vector
#define norm(v)   sqrt(norm2(v))   // norm = length of vector
#define d2(u,v)   norm2(u-v)       // distance squared = norm2 of difference
#define d(u,v)    norm(u-v)        // distance = norm of difference// poly_decimate(): - remove vertices to get a smaller approximate polygon
//    Input:  tol = approximation tolerance
//            V[] = polyline array of vertex points
//            n   = the number of points in V[]
//    Output: sV[]= reduced polyline vertexes (max is n)
//    Return: m   = the number of points in sV[]
int
poly_decimate( float tol, Point* V, int n, Point* sV )
{int    i, k, m, pv;             // misc countersfloat  tol2 = tol * tol;        // tolerance squaredPoint* vt = new Point[n];       // vertex bufferint*   mk = new int[n] = {0};   // marker  buffer// STAGE 1.  Vertex Reduction within tolerance of  prior vertex clustervt[0] = V[0];               // start at the beginningfor (i=k=1, pv=0; i<n; i++) {if (d2(V[i], V[pv]) < tol2)continue;vt[k++] = V[i];pv = i;}if (pv < n-1)vt[k++] = V[n-1];       // finish at the end// STAGE 2.  Douglas-Peucker polyline reductionmk[0] = mk[k-1] = 1;       //  mark the first and last vertexespoly_decimateDP( tol, vt, 0, k-1, mk  );// copy marked vertices to the reduced polylinefor (i=m=0; i<k; i++) {if (mk[i])sV[m++] =  vt[i];}delete vt;delete mk;return m;         //  m vertices in reduced polyline
}// poly_decimateDP():
//  This is the Douglas-Peucker recursive reduction routine
//  It marks vertexes that are part of the reduced polyline
//  for approximating the polyline subchain v[j] to v[k].
//    Input:  tol  = approximation tolerance
//            v[]  = polyline array of vertex points
//            j,k  = indices for the subchain v[j] to v[k]
//    Output: mk[] = array of markers matching vertex array v[]
void
poly_decimateDP( float tol, Point* v, int j, int k, int* mk )
{if (k <= j+1) // there is nothing to decimatereturn;// check for adequate approximation by segment S from v[j] to v[k]int     maxi = j;           // index of vertex farthest from Sfloat   maxd2 = 0;          // distance squared of farthest vertexfloat   tol2 = tol * tol;   // tolerance squaredSegment S = {v[j], v[k]};   // segment from v[j] to v[k]Vector  u = S.P1 - S.P0;    // segment direction vectordouble  cu = dot(u,u);      // segment length squared// test each vertex v[i] for max distance from S// compute using the Algorithm dist_Point_to_Segment()// Note: this works in any dimension (2D, 3D, ...)Vector  w;Point   Pb;                 // base of perpendicular from v[i] to Sdouble  b, cw, dv2;         // dv2 = distance v[i] to S squaredfor (int i=j+1; i<k; i++){// compute distance squaredw = v[i] - S.P0;cw = dot(w,u);if ( cw <= 0 )dv2 =d2(v[i], S.P0);else if ( cu <= cw )dv2 =d2(v[i], S.P1);else {b = cw / cu;Pb = S.P0 + b  * u;dv2 =d2(v[i], Pb);}// test with current max distance  squaredif (dv2 <= maxd2)continue;// v[i] is a new max vertexmaxi = i;maxd2 = dv2;}if (maxd2 > tol2)         // error is worse than the tolerance{// split the polyline at the farthest  vertex from Smk[maxi] = 1;       // mark v[maxi] for the reduced polyline// recursively decimate the two subpolylines at v[maxi]poly_decimateDP(  tol, v, j, maxi, mk );  // polyline v[j] to v[maxi]poly_decimateDP(  tol, v, maxi, k, mk );  // polyline v[maxi] to v[k]}// else the approximation is OK, so ignore intermediate vertexesreturn;
}