两种曲线点抽稀算法-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