中国领先的IT技术网站
|
|

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

本文主要介绍比较常用的两种抽稀算法的Python实现:道格拉斯-普克(Douglas-Peuker)算法和垂距限值法。

作者:佚名来源:开源中国|2017-08-21 17:15

开发者大赛路演 | 12月16日,技术创新,北京不见不散


何为抽稀

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

通俗的讲就是对曲线进行采样简化,即在曲线上取有限个点,将其变为折线,并且能够在一定程度保持原有形状。比较常用的两种抽稀算法是:道格拉斯-普克(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算法是从整体上考虑一条完整的曲线,实现时较垂距限值法复杂,但垂距限值法可能会在某些情况下导致局部最优。另外在实际使用中发现采用点到另外两点所在直线距离的方法来判断偏离,在曲线弧度比较大的情况下比较准确。如果在曲线弧度比较小,弯曲程度不明显时,这种方法抽稀效果不是很理想,建议使用三点所围成的三角形面积作为判断标准。下面是抽稀效果:

曲线点抽稀算法-Python实现

曲线点抽稀算法-Python实现

【编辑推荐】

  1. 我用Python做了六百万字的歌词分析,告诉你中国Rapper都在唱些啥
  2. 学习用Python编程时要避免的3个错误
  3. Python pyspider的安装与开发
  4. Python源码理解: '+=' 和 'xx = xx + xx' 的区别
  5. 我用Python爬了12万条影评,告诉你《战狼Ⅱ》都在说些啥
【责任编辑:未丽燕 TEL:(010)68476606】

点赞 0
分享:
大家都在看
猜你喜欢

读 书 +更多

C#2005编程进阶与参考手册

本书非常详细而全面地介绍了C#程序设计语言。本书不是“5分钟学习C#”式的手册,也不是那种教您“照猫画虎”地创建一些与您的实际工作需要...

订阅51CTO邮刊

点击这里查看样刊

订阅51CTO邮刊