|
前些天做到了有关于需要简化行政边界的项目，起初用到的是ArcGIS中的概化工具与简化面工具，虽然极大地简化了行政边界，但是出现了很多重复区域与缝隙，进行拓扑修复复杂而漫长。所以只能另辟蹊径，在mapshaper.org上找到了很好的解决方案。所以，本人分三篇博客，对其用到的三种简化面算法做个简单的介绍。 首先介绍的是著名的道格拉斯-普客算法： 它的基本思路是：对每一条曲线的首末点虚连一条直线，求所有点与直线的距离，并找出最大距离值dmax，用dmax与限差D相比： 若dmax＜D，这条曲线上的中间点全部舍去； 若dmax≥D，保留dmax对应的坐标点，并以该点为界，把曲线分为两部分，对这两部分重复使用该方法。  以下是以C#为例的算法实现代码：
///
/// Uses the Douglas Peucker algorithm to reduce the number of points.
///
/// The points.
/// The tolerance.
///
public static List DouglasPeuckerReduction
(List Points, Double Tolerance)
{
if (Points == null || Points.Count < 3)
return Points;
Int32 firstPoint = 0;
Int32 lastPoint = Points.Count - 1;
List pointIndexsToKeep = new List();
//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 returnPoints = new List();
pointIndexsToKeep.Sort();
foreach (Int32 index in pointIndexsToKeep)
{
returnPoints.Add(Points[index]);
}
return returnPoints;
}
///
/// Douglases the peucker reduction.
///
/// The points.
/// The first point.
/// The last point.
/// The tolerance.
/// The point index to keep.
private static void DouglasPeuckerReduction(List
points, Int32 firstPoint, Int32 lastPoint, Double tolerance,
ref List 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);
}
}
///
/// The distance of a point from a line made from point1 and point2.
///
/// The PT1.
/// The PT2.
/// The p.
///
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)²+(x1-x2)²) *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));
}
|