两种曲线点抽稀算法-Python实现 附代码

何为抽稀

在处理矢量化数据时,记录中往往会有很多重复数据,对进一步数据处理带来诸多不便。多余的数据一方面浪费了较多的存储空间,另一方面造成所要表达的图形不光滑或不符合标准。因此要通过某种规则,在保证矢量曲线形状不变的情况下,
最大限度地减少数据点个数,这个过程称为抽稀。

通俗的讲就是对曲线进行采样简化,即在曲线上取有限个点,将其变为折线,并且能够在一定程度保持原有形状。比较常用的两种抽稀算法是:道格拉斯-普克(Douglas-Peuker)算法和垂距限值法。

道格拉斯-普克(Douglas-Peuker)算法

Douglas-Peuker算法(DP算法)过程如下:

  1. 连接曲线首尾两点A、B;
  2. 依次计算曲线上所有点到A、B两点所在曲线的距离;
  3. 计算最大距离D,如果D小于阈值threshold,则去掉曲线上出A、B外的所有点;如果D大于阈值threshold,则把曲线以最大距离分割成两段;
  4. 对所有曲线分段重复1-3步骤,知道所有D均小于阈值。即完成抽稀。

这种算法的抽稀精度与阈值有很大关系,阈值越大,简化程度越大,点减少的越多;反之简化程度越低,点保留的越多,形状也越趋于原曲线。

下面是Python代码实现:


  1. # -*- coding: utf-8 -*- 
  2. """ 
  3. ------------------------------------------------- 
  4.   File Name:    DouglasPeuker 
  5.   Description :  道格拉斯-普克抽稀算法 
  6.   Author :        J_hao 
  7.   date:          2017/8/16 
  8. ------------------------------------------------- 
  9.   Change Activity: 
  10.                   2017/8/16: 道格拉斯-普克抽稀算法 
  11. ------------------------------------------------- 
  12. """ 
  13. from __future__ import division 
  14.  
  15. from math import sqrt, pow 
  16.  
  17. __author__ = 'J_hao' 
  18.  
  19. THRESHOLD = 0.0001  # 阈值 
  20.  
  21.  
  22. def point2LineDistance(point_a, point_b, point_c): 
  23.     """ 
  24.     计算点a到点b c所在直线的距离 
  25.     :param point_a: 
  26.     :param point_b: 
  27.     :param point_c: 
  28.     :return: 
  29.     """ 
  30.     # 首先计算b c 所在直线的斜率和截距 
  31.     if point_b[0] == point_c[0]: 
  32.         return 9999999 
  33.     slope = (point_b[1] - point_c[1]) / (point_b[0] - point_c[0]) 
  34.     intercept = point_b[1] - slope * point_b[0] 
  35.  
  36.     # 计算点a到b c所在直线的距离 
  37.     distance = abs(slope * point_a[0] - point_a[1] + intercept) / sqrt(1 + pow(slope, 2)) 
  38.     return distance 
  39.  
  40.  
  41. class DouglasPeuker(object): 
  42.     def __init__(self): 
  43.         self.threshold = THRESHOLD 
  44.         self.qualify_list = list() 
  45.         self.disqualify_list = list() 
  46.  
  47.     def diluting(self, point_list): 
  48.         """ 
  49.         抽稀 
  50.         :param point_list:二维点列表 
  51.         :return: 
  52.         """ 
  53.         if len(point_list) < 3: 
  54.             self.qualify_list.extend(point_list[::-1]) 
  55.         else: 
  56.             # 找到与收尾两点连线距离最大的点 
  57.             max_distance_index, max_distance = 0, 0 
  58.             for index, point in enumerate(point_list): 
  59.                 if index in [0, len(point_list) - 1]: 
  60.                     continue 
  61.                 distance = point2LineDistance(point, point_list[0], point_list[-1]) 
  62.                 if distance > max_distance: 
  63.                     max_distance_index = index 
  64.                     max_distance = distance 
  65.  
  66.             # 若最大距离小于阈值,则去掉所有中间点。 反之,则将曲线按最大距离点分割 
  67.             if max_distance < self.threshold: 
  68.                 self.qualify_list.append(point_list[-1]) 
  69.                 self.qualify_list.append(point_list[0]) 
  70.             else: 
  71.                 # 将曲线按最大距离的点分割成两段 
  72.                 sequence_a = point_list[:max_distance_index] 
  73.                 sequence_b = point_list[max_distance_index:] 
  74.  
  75.                 for sequence in [sequence_a, sequence_b]: 
  76.                     if len(sequence) < 3 and sequence == sequence_b: 
  77.                         self.qualify_list.extend(sequence[::-1]) 
  78.                     else: 
  79.                         self.disqualify_list.append(sequence) 
  80.  
  81.     def main(self, point_list): 
  82.         self.diluting(point_list) 
  83.         while len(self.disqualify_list) > 0: 
  84.             self.diluting(self.disqualify_list.pop()) 
  85.         print self.qualify_list 
  86.         print len(self.qualify_list) 
  87.  
  88.  
  89. if __name__ == '__main__': 
  90.     d = DouglasPeuker() 
  91.     d.main([[104.066228, 30.644527], [104.066279, 30.643528], [104.066296, 30.642528], [104.066314, 30.641529], 
  92.             [104.066332, 30.640529], [104.066383, 30.639530], [104.066400, 30.638530], [104.066451, 30.637531], 
  93.             [104.066468, 30.636532], [104.066518, 30.635533], [104.066535, 30.634533], [104.066586, 30.633534], 
  94.             [104.066636, 30.632536], [104.066686, 30.631537], [104.066735, 30.630538], [104.066785, 30.629539], 
  95.             [104.066802, 30.628539], [104.066820, 30.627540], [104.066871, 30.626541], [104.066888, 30.625541], 
  96.             [104.066906, 30.624541], [104.066924, 30.623541], [104.066942, 30.622542], [104.066960, 30.621542], 
  97.             [104.067011, 30.620543], [104.066122, 30.620086], [104.065124, 30.620021], [104.064124, 30.620022], 
  98.             [104.063124, 30.619990], [104.062125, 30.619958], [104.061125, 30.619926], [104.060126, 30.619894], 
  99.             [104.059126, 30.619895], [104.058127, 30.619928], [104.057518, 30.620722], [104.057625, 30.621716], 
  100.             [104.057735, 30.622710], [104.057878, 30.623700], [104.057984, 30.624694], [104.058094, 30.625688], 
  101.             [104.058204, 30.626682], [104.058315, 30.627676], [104.058425, 30.628670], [104.058502, 30.629667], 
  102.             [104.058518, 30.630667], [104.058503, 30.631667], [104.058521, 30.632666], [104.057664, 30.633182], 
  103.             [104.056664, 30.633174], [104.055664, 30.633166], [104.054672, 30.633289], [104.053758, 30.633694], 
  104.             [104.052852, 30.634118], [104.052623, 30.635091], [104.053145, 30.635945], [104.053675, 30.636793], 
  105.             [104.054200, 30.637643], [104.054756, 30.638475], [104.055295, 30.639317], [104.055843, 30.640153], 
  106.             [104.056387, 30.640993], [104.056933, 30.641830], [104.057478, 30.642669], [104.058023, 30.643507], 
  107.             [104.058595, 30.644327], [104.059152, 30.645158], [104.059663, 30.646018], [104.060171, 30.646879], 
  108.             [104.061170, 30.646855], [104.062168, 30.646781], [104.063167, 30.646823], [104.064167, 30.646814], 
  109.             [104.065163, 30.646725], [104.066157, 30.646618], [104.066231, 30.645620], [104.066247, 30.644621], ]) 

垂距限值法

垂距限值法其实和DP算法原理一样,但是垂距限值不是从整体角度考虑,而是依次扫描每一个点,检查是否符合要求。

算法过程如下:

  1. 以第二个点开始,计算第二个点到前一个点和后一个点所在直线的距离d;
  2. 如果d大于阈值,则保留第二个点,计算第三个点到第二个点和第四个点所在直线的距离d;若d小于阈值则舍弃第二个点,计算第三个点到第一个点和第四个点所在直线的距离d;
  3. 依次类推,直线曲线上倒数第二个点。

下面是Python代码实现:


  1. # -*- coding: utf-8 -*- 
  2. """ 
  3. ------------------------------------------------- 
  4.   File Name:    LimitVerticalDistance 
  5.   Description :  垂距限值抽稀算法 
  6.   Author :        J_hao 
  7.   date:          2017/8/17 
  8. ------------------------------------------------- 
  9.   Change Activity: 
  10.                   2017/8/17: 
  11. ------------------------------------------------- 
  12. """ 
  13. from __future__ import division 
  14.  
  15. from math import sqrt, pow 
  16.  
  17. __author__ = 'J_hao' 
  18.  
  19. THRESHOLD = 0.0001  # 阈值 
  20.  
  21.  
  22. def point2LineDistance(point_a, point_b, point_c): 
  23.     """ 
  24.     计算点a到点b c所在直线的距离 
  25.     :param point_a: 
  26.     :param point_b: 
  27.     :param point_c: 
  28.     :return: 
  29.     """ 
  30.     # 首先计算b c 所在直线的斜率和截距 
  31.     if point_b[0] == point_c[0]: 
  32.         return 9999999 
  33.     slope = (point_b[1] - point_c[1]) / (point_b[0] - point_c[0]) 
  34.     intercept = point_b[1] - slope * point_b[0] 
  35.  
  36.     # 计算点a到b c所在直线的距离 
  37.     distance = abs(slope * point_a[0] - point_a[1] + intercept) / sqrt(1 + pow(slope, 2)) 
  38.     return distance 
  39.  
  40.  
  41. class LimitVerticalDistance(object): 
  42.     def __init__(self): 
  43.         self.threshold = THRESHOLD 
  44.         self.qualify_list = list() 
  45.  
  46.     def diluting(self, point_list): 
  47.         """ 
  48.         抽稀 
  49.         :param point_list:二维点列表 
  50.         :return: 
  51.         """ 
  52.         self.qualify_list.append(point_list[0]) 
  53.         check_index = 1 
  54.         while check_index < len(point_list) - 1: 
  55.             distance = point2LineDistance(point_list[check_index], 
  56.                                           self.qualify_list[-1], 
  57.                                           point_list[check_index + 1]) 
  58.  
  59.             if distance < self.threshold: 
  60.                 check_index += 1 
  61.             else: 
  62.                 self.qualify_list.append(point_list[check_index]) 
  63.                 check_index += 1 
  64.         return self.qualify_list 
  65.  
  66.  
  67. if __name__ == '__main__': 
  68.     l = LimitVerticalDistance() 
  69.     diluting = l.diluting([[104.066228, 30.644527], [104.066279, 30.643528], [104.066296, 30.642528], [104.066314, 30.641529], 
  70.             [104.066332, 30.640529], [104.066383, 30.639530], [104.066400, 30.638530], [104.066451, 30.637531], 
  71.             [104.066468, 30.636532], [104.066518, 30.635533], [104.066535, 30.634533], [104.066586, 30.633534], 
  72.             [104.066636, 30.632536], [104.066686, 30.631537], [104.066735, 30.630538], [104.066785, 30.629539], 
  73.             [104.066802, 30.628539], [104.066820, 30.627540], [104.066871, 30.626541], [104.066888, 30.625541], 
  74.             [104.066906, 30.624541], [104.066924, 30.623541], [104.066942, 30.622542], [104.066960, 30.621542], 
  75.             [104.067011, 30.620543], [104.066122, 30.620086], [104.065124, 30.620021], [104.064124, 30.620022], 
  76.             [104.063124, 30.619990], [104.062125, 30.619958], [104.061125, 30.619926], [104.060126, 30.619894], 
  77.             [104.059126, 30.619895], [104.058127, 30.619928], [104.057518, 30.620722], [104.057625, 30.621716], 
  78.             [104.057735, 30.622710], [104.057878, 30.623700], [104.057984, 30.624694], [104.058094, 30.625688], 
  79.             [104.058204, 30.626682], [104.058315, 30.627676], [104.058425, 30.628670], [104.058502, 30.629667], 
  80.             [104.058518, 30.630667], [104.058503, 30.631667], [104.058521, 30.632666], [104.057664, 30.633182], 
  81.             [104.056664, 30.633174], [104.055664, 30.633166], [104.054672, 30.633289], [104.053758, 30.633694], 
  82.             [104.052852, 30.634118], [104.052623, 30.635091], [104.053145, 30.635945], [104.053675, 30.636793], 
  83.             [104.054200, 30.637643], [104.054756, 30.638475], [104.055295, 30.639317], [104.055843, 30.640153], 
  84.             [104.056387, 30.640993], [104.056933, 30.641830], [104.057478, 30.642669], [104.058023, 30.643507], 
  85.             [104.058595, 30.644327], [104.059152, 30.645158], [104.059663, 30.646018], [104.060171, 30.646879], 
  86.             [104.061170, 30.646855], [104.062168, 30.646781], [104.063167, 30.646823], [104.064167, 30.646814], 
  87.             [104.065163, 30.646725], [104.066157, 30.646618], [104.066231, 30.645620], [104.066247, 30.644621], ]) 
  88.     print len(diluting) 
  89.     print(diluting) 

最后

其实DP算法和垂距限值法原理一样,DP算法是从整体上考虑一条完整的曲线,实现时较垂距限值法复杂,但垂距限值法可能会在某些情况下导致局部最优。另外在实际使用中发现采用点到另外两点所在直线距离的方法来判断偏离,在曲线弧度比较大的情况下比较准确。如果在曲线弧度比较小,弯曲程度不明显时,这种方法抽稀效果不是很理想,建议使用三点所围成的三角形面积作为判断标准。下面是抽稀效果:

作者:佚名

来源:51CTO

时间: 2024-08-28 12:42:50

两种曲线点抽稀算法-Python实现 附代码的相关文章

两种php给图片加水印的实现代码_php实例

PHP最简单的加水印方法 <?php $img = imagecreatefromjpeg($filename); $logo = imagecreatefromjpeg($filename); /*imagecraetefromjpeg-由文件或URL创建一个新图像 imagecreatefromjpeg(string $filename) 如果启用了fopen包装器,URL可以作为文件名*/ imagecopy($img,$logo,15,15,0,0,$width,$height); /*

两种简单实现菜单高亮显示的JS类代码_javascript技巧

记得当年写静态页时,为了实现高亮都是在每个页面加不同的样式,呵.高亮显示我觉得对于web前端来说,是比较常用到的效果,正好此次又要用到,特地整理出我所写的两种高亮类. 其实思路很简单,第一种方法是通过遍历链接组的href值,通过indexOf判断href值是否被包含在浏览器当前url值中.此方法有一定局限,比如对于iframe内的菜单是不能这样判断的; 第二种方法适用范围更广一样,实现思路也比较简单,即通过判断点击,给点击项加载高亮样式. 第一种判断当前URL值高亮类代码: 复制代码 代码如下:

JavaScript判断两种格式的输入日期的正确性的代码_javascript技巧

最简单的 复制代码 代码如下: function isValidDate(dateStr) {             var matchArray = dateStr.match(/^[0-9]+-[0-1][0-9]-[0-3][0-9]$/)             if (matchArray == null) {               alert("Invalid date: " + dateStr);               return false;      

python实现中文输出的两种方法

  这篇文章主要介绍了python实现中文输出的两种方法,实例分析了Python操作中文输出的技巧,需要的朋友可以参考下 方法一: 用encode和decode 如: ? 1 2 3 4 5 6 7 8 9 10 11 import os.path import xlrd,sys Filename='/home/tom/Desktop/1234.xls' if not os.path.isfile(Filename): raise NameError,"%s is not a valid fil

python 列表转化为字符串的两种方式

Python 列表转化为字符串的两种方式 (1)方式一 Python代码   >>> juice=['orange','a','b']   >>> ''.join(juice)   'orangeab'       (2)方式二: Python代码   >>> juice=['orange','a','b']   >>> content='%s'*len(juice) % tuple(juice)   >>> pri

gridview实现服务器端和客户端全选的两种方法

 这篇文章主要介绍了gridview实现服务器端和客户端全选的两种方法,需要的朋友可以参考下  代码如下: <%@ Page Language="C#" AutoEventWireup="true"%>   <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-tran

js判断字符是否是汉字的两种方法小结

 本篇文章主要是对js判断字符是否是汉字的两种方法进行了详细的总结介绍,需要的朋友可以过来参考下,希望对大家有所帮助 有时需要判断一个字符是不是汉字,比如在用户输入含有中英文的内容时,需要判断是否超过规定长度就要用到.用 Javascript 判断通常有两种方法.    1.用正则表达式判断    代码如下: <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org

ReentrantLock和synchronized两种锁定机制

ReentrantLock和synchronized两种锁定机制 应用synchronized同步锁 把代码块声明为 synchronized,使得该代码具有 原子性(atomicity)和 可见性(visibility). 原子性意味着一个线程一次只能执行由一个指定监控对象(lock)保护的代码,从而防止多个线程在更新共享状态时相互冲突. 可见性类似volatile关键字. 应用ReentrantLock显示锁 ReentrantLock 类实现了 Lock ,它拥有与 synchronize

Android通过AsyncTask与ThreadPool(线程池)两种方式异步加载大量数据的分析与对比

   在加载大量数据的时候,经常会用到异步加载,所谓异步加载,就是把耗时的工作放到子线程里执行,当数据加载完毕的时候再到主线程进行UI刷新.在数据量非常大的情况下,我们通常会使用两种技术来进行异步加载,一是通过AsyncTask来实现,另一种方式则是通过ThreadPool来实现,今天我们就通过一个例子来讲解和对比这两种实现方式.     本文原创,如需转载,请注明转载地址http://blog.csdn.net/carrey1989/article/details/12002033     项