道格拉斯-普克抽稀(Douglas-Peucker)算法是一个常用的线抽稀算法，将一条线在不改变大概形状的情况下，使用最少数量的点描述。基本原理如下图：
通过Google我找到了两段代码，一段是c#的，一段是python的。现将其分享出来
C#版的:
/// <summary>
/// Uses the Douglas Peucker algorithm to reduce the number of points.
/// </summary>
/// <param name="Points">The points.</param>
/// <param name="Tolerance">The tolerance.</param>
/// <returns></returns>
public static List<Point> DouglasPeuckerReduction
(List<Point> Points, Double Tolerance)
{
if (Points == null || Points.Count < 3)
return Points;
Int32 firstPoint = 0;
Int32 lastPoint = Points.Count - 1;
List<Int32> pointIndexsToKeep = new List<Int32>();
//Add the first and last index to the keepers
pointIndexsToKeep.Add(firstPoint);
pointIndexsToKeep.Add(lastPoint);
//The first and the last point cannot be the same
while (Points[firstPoint].Equals(Points[lastPoint]))
{
lastPoint--;
}
DouglasPeuckerReduction(Points, firstPoint, lastPoint,
Tolerance, ref pointIndexsToKeep);
List<Point> returnPoints = new List<Point>();
pointIndexsToKeep.Sort();
foreach (Int32 index in pointIndexsToKeep)
{
returnPoints.Add(Points[index]);
}
return returnPoints;
}
/// <summary>
/// Douglases the peucker reduction.
/// </summary>
/// <param name="points">The points.</param>
/// <param name="firstPoint">The first point.</param>
/// <param name="lastPoint">The last point.</param>
/// <param name="tolerance">The tolerance.</param>
/// <param name="pointIndexsToKeep">The point index to keep.</param>
private static void DouglasPeuckerReduction(List<Point>
points, Int32 firstPoint, Int32 lastPoint, Double tolerance,
ref List<Int32> pointIndexsToKeep)
{
Double maxDistance = 0;
Int32 indexFarthest = 0;
for (Int32 index = firstPoint; index < lastPoint; index++)
{
Double distance = PerpendicularDistance
(points[firstPoint], points[lastPoint], points[index]);
if (distance > maxDistance)
{
maxDistance = distance;
indexFarthest = index;
}
}
if (maxDistance > tolerance && indexFarthest != 0)
{
//Add the largest point that exceeds the tolerance
pointIndexsToKeep.Add(indexFarthest);
DouglasPeuckerReduction(points, firstPoint,
indexFarthest, tolerance, ref pointIndexsToKeep);
DouglasPeuckerReduction(points, indexFarthest,
lastPoint, tolerance, ref pointIndexsToKeep);
}
}
/// <summary>
/// The distance of a point from a line made from point1 and point2.
/// </summary>
/// <param name="pt1">The PT1.</param>
/// <param name="pt2">The PT2.</param>
/// <param name="p">The p.</param>
/// <returns></returns>
public static Double PerpendicularDistance
(Point Point1, Point Point2, Point Point)
{
//Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)| *Area of triangle
//Base = v((x1-x2)2+(x1-x2)2) *Base of Triangle*
//Area = .5*Base*H *Solve for height
//Height = Area/.5/Base
Double area = Math.Abs(.5 * (Point1.X * Point2.Y + Point2.X *
Point.Y + Point.X * Point1.Y - Point2.X * Point1.Y - Point.X *
Point2.Y - Point1.X * Point.Y));
Double bottom = Math.Sqrt(Math.Pow(Point1.X - Point2.X, 2) +
Math.Pow(Point1.Y - Point2.Y, 2));
Double height = area / bottom * 2;
return height;
//Another option
//Double A = Point.X - Point1.X;
//Double B = Point.Y - Point1.Y;
//Double C = Point2.X - Point1.X;
//Double D = Point2.Y - Point1.Y;
//Double dot = A * C + B * D;
//Double len_sq = C * C + D * D;
//Double param = dot / len_sq;
//Double xx, yy;
//if (param < 0)
//{
// xx = Point1.X;
// yy = Point1.Y;
//}
//else if (param > 1)
//{
// xx = Point2.X;
// yy = Point2.Y;
//}
//else
//{
// xx = Point1.X + param * C;
// yy = Point1.Y + param * D;
//}
//Double d = DistanceBetweenOn2DPlane(Point, new Point(xx, yy));
}
python 版的:
# pure-Python Douglas-Peucker line simplification/generalization
#
# this code was written by Schuyler Erle <schuyler@nocat.net> and is
# made available in the public domain.
#
# the code was ported from a freely-licensed example at
# http://www.3dsoftware.com/Cartography/Programming/PolyLineReduction/
#
# the original page is no longer available, but is mirrored at
# http://www.mappinghacks.com/code/PolyLineReduction/
"""
>>> line = [(0,0),(1,0),(2,0),(2,1),(2,2),(1,2),(0,2),(0,1),(0,0)]
>>> simplify_points(line, 1.0)
[(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)]
>>> line = [(0,0),(0.5,0.5),(1,0),(1.25,-0.25),(1.5,.5)]
>>> simplify_points(line, 0.25)
[(0, 0), (0.5, 0.5), (1.25, -0.25), (1.5, 0.5)]
"""
import math
def simplify_points (pts, tolerance):
anchor = 0
floater = len(pts) - 1
stack = []
keep = set()
stack.append((anchor, floater))
while stack:
anchor, floater = stack.pop()
# initialize line segment
if pts[floater] != pts[anchor]:
anchorX = float(pts[floater][0] - pts[anchor][0])
anchorY = float(pts[floater][1] - pts[anchor][1])
seg_len = math.sqrt(anchorX ** 2 + anchorY ** 2)
# get the unit vector
anchorX /= seg_len
anchorY /= seg_len
else:
anchorX = anchorY = seg_len = 0.0
# inner loop:
max_dist = 0.0
farthest = anchor + 1
for i in range(anchor + 1, floater):
dist_to_seg = 0.0
# compare to anchor
vecX = float(pts[i][0] - pts[anchor][0])
vecY = float(pts[i][1] - pts[anchor][1])
seg_len = math.sqrt( vecX ** 2 + vecY ** 2 )
# dot product:
proj = vecX * anchorX + vecY * anchorY
if proj < 0.0:
dist_to_seg = seg_len
else:
# compare to floater
vecX = float(pts[i][0] - pts[floater][0])
vecY = float(pts[i][1] - pts[floater][1])
seg_len = math.sqrt( vecX ** 2 + vecY ** 2 )
# dot product:
proj = vecX * (-anchorX) + vecY * (-anchorY)
if proj < 0.0:
dist_to_seg = seg_len
else: # calculate perpendicular distance to line (pythagorean theorem):
dist_to_seg = math.sqrt(abs(seg_len ** 2 - proj ** 2))
if max_dist < dist_to_seg:
max_dist = dist_to_seg
farthest = i
if max_dist <= tolerance: # use line segment
keep.add(anchor)
keep.add(floater)
else:
stack.append((anchor, farthest))
stack.append((farthest, floater))
keep = list(keep)
keep.sort()
return [pts[i] for i in keep]
if __name__ == "__main__":
import doctest
doctest.testmod()
