GIS算法之DouglasPeucher算法

      在数字化过程中，需要对曲线进行采样简化，达到压缩数据的作用而又能一定程度上保持曲线形状。做法是在曲线上取有限个点，将其变为折现。

      算法过程：
      1 连接曲线首尾构造直线AB（曲线的弦）
      2 得到离该直线段距离最大的点C，并计算其与AB的距离d
      3 比较d与阈值的小小关系，如果小于阈值，结束
      4 如果大于阈值，则用C将曲线分为AC和BC两段，并分别对两端曲线作1-3步的处理；
      5 当所有的曲线都处理完毕后，一次连接各个分割点形成折现

      下面是直线类，功能比较简单，加入一个计算点到该线距离的方法：

 

代码

using  System;
using  System.Collections.Generic;
using  System.Linq;
using  System.Text;
using  System.Drawing;

namespace  DouglasPeucher
{
     class  Line
    {
         private  Point fromPoint;
         private  Point toPoint;

         public  Point FromPoint
        {
             get  {  return  fromPoint; }
             set  { fromPoint  =  value; }
        }
         public  Point ToPoint
        {
             get  {  return  toPoint; }
             set  { toPoint  =  value; }
        }


         public   void  PutCoord(Point fromPoint, Point toPoint)
        {
             this .fromPoint  =  fromPoint;
             this .toPoint  =  toPoint;
        }

         public   double  QueryDistance(Point point)
        {
             double  distance  =   0 ;
             if  (fromPoint.X  ==  toPoint.X)  // 直线垂直
            {
                distance  =  Math.Abs(toPoint.X  -  fromPoint.X);    
            }
             else   if  (fromPoint.Y  ==  toPoint.Y) // 直线水平
            {
                distance  =  Math.Abs(point.Y  -  fromPoint.Y);  
            }
             else    // 一般情况
            {
                 double  slope  =  Convert.ToDouble((fromPoint.Y  -  toPoint.Y))  /  (fromPoint.X  -  toPoint.X);
                 double  uprightSlope  =   - 1   /  slope;
                 // 得到两2直线交点
                 double  x  =  (slope  *  fromPoint.X  -  fromPoint.Y  +  point.Y  +  point.X  /  slope)
                            /  (slope  +   1.0   /  slope);
                 double  y  =  fromPoint.Y  +  slope * (x  -  fromPoint.X);
                 Point uprightPoint  = new  Point(( int )x,( int )y);
                distance  = GetDistance(point,uprightPoint); 
            }
             return  distance;
        }

         public   static   double  GetDistance(Point point1, Point point2)
        {
             double  distance  =  Math.Sqrt((point1.Y  -  point2.Y)  *  (point1.Y  -  point2.Y)  +  (point1.X  -  point2.X)  *  (point1.X  -  point2.X));
             return  distance;
        }
    }
}

 

        主要算法过程：

 

代码

private   void  btnSimplify_Click( object  sender, EventArgs e)
        {
            lemen  =   int .Parse(tbLimen.Text);
            CurveDisperse(points,lemen); 
             // 对点集合进行顺序纠正
           

            resultPoints.Insert( 0 , points[ 0 ]);
            resultPoints.Insert(resultPoints.Count  -   1 , points[points.Count  -   1 ]);
             // 删除重复的点
             for  ( int  k1  =   0 ; k1  <  resultPoints.Count; k1 ++ )
            {
                Point p1  =  points[k1];
                 for  ( int  k2  =  k1  +   1 ; k2  <  resultPoints.Count; k2 ++ )
                {
                    Point p2  =  resultPoints[k2];
                     if  (p1.X  ==  p2.X  &&  p1.Y  ==  p2.Y)
                        resultPoints.RemoveAt(k2);
                }
            }
             // 纠正结果点集的排序
            CorrectOrder();
             // 绘制到屏幕上
            DrawToScreen(resultPoints);
        }

         private   void  CorrectOrder()
        {
             int [] indexs  =   new   int [resultPoints.Count];
             for  ( int  i  =   0 ; i  <  resultPoints.Count; i ++ )
            {
                Point rst  =  resultPoints[i];
                 for  ( int  j  =   0 ; j  <  points.Count; j ++ )
                {
                    Point pt  =  points[j];
                     if  (pt.X  ==  rst.X  &&  pt.Y  ==  rst.Y)
                    {
                        indexs[i]  =  j;
                         break ;
                    }
                }
            }

             int  temp;
            Point tempPoint;
             for  ( int  i  =   0 ; i  <  indexs.Length; i ++ )
            {
                 for  ( int  j  =  i  +   1 ; j  <  indexs.Length; j ++ )
                {
                     if  (indexs[i]  >=  indexs[j])
                    {
                        temp  =  indexs[i];
                        indexs[i]  =  indexs[j];
                        indexs[j]  =  temp;
                        tempPoint  =  resultPoints[i];
                        resultPoints[i]  =  resultPoints[j];
                        resultPoints[j]  =  tempPoint;
                    }
                }
            }
        }

         private   void  CurveDisperse(List < Point >  pointColl, double  lemen)
        {
             /* 算法过程
             1 连接曲线首尾构造直线AB（曲线的弦）
             2 得到离该直线段距离最大的点C，并计算其与AB的距离d
             3 比较d与阈值的小小关系，如果小于阈值，结束
             4 如果大于阈值，则用C将曲线分为AC和BC两段，并分别对两端曲线作1-3步的处理；
             5 当所有的曲线都处理完毕后，一次连接各个分割点形成折现
             *  */
            Line line  =   new  Line();
            line.PutCoord(pointColl[ 0 ], pointColl[pointColl.Count  -   1 ]);
             double  maxDistance;
             int  index;
            SetMaxDistance(pointColl, line,  out  maxDistance,  out  index);
             if  (maxDistance  >=  lemen)
            {
                resultPoints.Add(pointColl[index]);
                 // 构造2个新的曲线
                List < Point >  curveLeft  =  GetSubList( 0 , index, pointColl);
                List < Point >  curveRight  =  GetSubList(index, pointColl.Count  -   1 , pointColl);
                CurveDisperse(curveLeft, lemen);
                CurveDisperse(curveRight, lemen);
            }

        }

         private  List < Point >  GetSubList( int  sIndex, int  eIndex,List < Point >  points)
        {
            List < Point >  subList  = new  List < Point > ();
             if  (sIndex  <   0   ||  eIndex  <   0   ||  sIndex  >  points.Count  -   1   ||  eIndex  >  points.Count  -   1 )
                 throw   new  Exception( " 索引出错！ " );

             for  ( int  i  =  sIndex; i  <=  eIndex; i ++ )
            {
                subList.Add(points[i]);
            }
             return  subList;
        }

         // 设置最大距离和索引值
         private   void  SetMaxDistance(List < Point >  pointColl,Line line, out   double  maxDistance, out   int  index)
        {
             double  distance;
            maxDistance  =   double .MinValue;
            index  =   0 ;
             for  ( int  i  =   0 ; i  <  pointColl.Count;i ++ )
            {
                Point pt  =  pointColl[i];
                distance  =  line.QueryDistance(pt);
                 if  (distance  >  maxDistance)
                {
                    index  =  i;
                    maxDistance  =  distance;
                }
            }
        }

private   void  DrawToScreen(List < Point >  pointColl)
        {
             for  ( int  i  =   0 ; i  <  pointColl.Count - 1 ; i ++ )
            {
                g.DrawLine(Pens.Blue, pointColl[i], pointColl[i + 1 ]);  
            }
            pictureBox1.Refresh();
        }

 

      下面是测试结果图示，折现即为曲线离散化的结果：

转载于:https://www.cnblogs.com/gis_sky/archive/2010/05/11/DouglasPeucher.html