正交点集-外轮廓线 算法实现

landsat99

// 分布于正交栅格网的散点集合（允许有内部空洞），生成正交点集合的外轮廓线。


说明：
1.  此算法分区施加右手法则，贴合正交外轮廓。
2.  对任意散点的集合，外轮廓线不唯一。 需附加约束条件，比如面积最大、 轮廓线周长最大、最小等。
3.  Python演示算法及图形说明。
4.  此正交算法 为纯原创代码，未参考其它几何算法集。不严谨之处 多多交流指正。

1. # RecVec_grid.py
   
2. import matplotlib.pyplot as plt
   
3. import GenList as Gen
   
4. import MaxP as MP
   
5. 
   
6. 
   
7. def Direc(UpDown, A, B):
   
8.     alpha = 360 + 10
   
9.     if UpDown == 1:  # 象限号
   
10.         if B[1] == A[1] and B[0] > A[0]:
    
11.             alpha = 0  #
    
12.         elif B[1] == A[1] and B[0] < A[0]:
    
13.             alpha = 180
    
14.         elif B[0] == A[0] and B[1] > A[1]:
    
15.             alpha = 90
    
16.         elif B[0] == A[0] and B[1] < A[1]:
    
17.             alpha = 270
    
18.         else:
    
19.             pass
    
20.     if UpDown == 2:  # 象限号
    
21.         if B[0] == A[0] and B[1] > A[1]:
    
22.             alpha = 0
    
23.         elif B[0] == A[0] and B[1] < A[1]:
    
24.             alpha = 180
    
25.         elif B[1] == A[1] and B[0] > A[0]:
    
26.             alpha = 270
    
27.         elif B[1] == A[1] and B[0] < A[0]:
    
28.             alpha = 90
    
29.         else:
    
30.             pass
    
31.     if UpDown == 3:  # 象限号
    
32.         if B[1] == A[1] and B[0] < A[0]:
    
33.             alpha = 0
    
34.         elif B[1] == A[1] and B[0] > A[0]:
    
35.             alpha = 180
    
36.         elif B[0] == A[0] and B[1] < A[1]:
    
37.             alpha = 90
    
38.         elif B[0] == A[0] and B[1] > A[1]:
    
39.             alpha = 270
    
40.         else:
    
41.             pass
    
42.     if UpDown == 4:  # 象限号
    
43.         if B[0] == A[0] and B[1] < A[1]:
    
44.             alpha = 0
    
45.         elif B[0] == A[0] and B[1] > A[1]:
    
46.             alpha = 180
    
47.         elif B[1] == A[1] and B[0] < A[0]:
    
48.             alpha = 270
    
49.         elif B[1] == A[1] and B[0] > A[0]:
    
50.             alpha = 90
    
51.         else:
    
52.             pass
    
53.     return alpha
    
54. 
    
55. 
    
56. 
    
57. def Dist(A, B):
    
58.     return (B[1] - A[1])**2 + (B[0] - A[0])**2
    
59. 
    
60. 
    
61. def NextPoi(UpDown, Poi, PList):
    
62.     Dire = 360 + 10
    
63.     Unit = 1
    
64.     for id, p in enumerate(PList):
    
65.         if Direc(UpDown, Poi, p) <= Dire and Dist(Poi, p) == Unit:
    
66.             Dire = Direc(UpDown, Poi, p)
    
67.             dis = Dist(Poi, p)
    
68.             Num = id
    
69.     return PList[Num]
    
70. 
    
71. 
    
72. def plot2D(x, y, x_out, y_out):
    
73.     plt.scatter(x, y)  # 画散点图
    
74.     plt.plot(x_out, y_out)  # 画连线图
    
75.     plt.show()
    
76. 
    
77. if __name__ == '__main__':
    
78.     Points = Gen.GenList()
    
79.     XMaxPoi, Num = MP.XMaxPoint(Points)
    
80.     YMaxPoi, Num = MP.YMaxPoint(Points)
    
81.     XMinPoi, Num = MP.XMinPoint(Points)
    
82.     YMinPoi, Num = MP.YMinPoint(Points)
    
83.     PPs = Points.copy()
    
84.     PPs.remove(XMaxPoi)
    
85.     Poi_A = XMaxPoi
    
86.     RList = [XMaxPoi]
    
87.     Poi_B = []
    
88.     while Poi_B != YMaxPoi:
    
89.         Poi_B = NextPoi(1, Poi_A, PPs)
    
90.         RList.append(Poi_B)
    
91.         PPs.remove(Poi_B)
    
92.         Poi_A = Poi_B
    
93.     Poi_A = YMaxPoi
    
94.     Poi_B = []
    
95.     while Poi_B != XMinPoi:
    
96.         Poi_B = NextPoi(2, Poi_A, PPs)
    
97.         RList.append(Poi_B)
    
98.         PPs.remove(Poi_B)
    
99.         Poi_A = Poi_B
    
100.     Poi_A = XMinPoi
     
101.     Poi_B = []
     
102.     while Poi_B != YMinPoi:
     
103.         Poi_B = NextPoi(3, Poi_A, PPs)
     
104.         RList.append(Poi_B)
     
105.         PPs.remove(Poi_B)
     
106.         Poi_A = Poi_B
     
107.     Poi_A = YMinPoi
     
108.     Poi_B = []
     
109.     PPs = [XMaxPoi] + PPs  # 加回第一点
     
110.     while Poi_B != XMaxPoi:
     
111.         Poi_B = NextPoi(4, Poi_A, PPs)
     
112.         RList.append(Poi_B)
     
113.         PPs.remove(Poi_B)
     
114.         Poi_A = Poi_B
     
115.     print(RList)
     
116. 
     
117.     XList = [i[0] for i in Points]
     
118.     YList = [i[1] for i in Points]
     
119.     X_OutList = [i[0] for i in RList]
     
120.     Y_OutList = [i[1] for i in RList]
     
121.     plot2D(XList, YList, X_OutList, Y_OutList)
     

# RecVec_grid.py<br /> import matplotlib.pyplot as plt<br /> import GenList as Gen<br /> import MaxP as MP<br /> <br /> <br /> def Direc(UpDown, A, B):<br />     alpha = 360 + 10<br />     if UpDown == 1:  # 象限号<br />         if B[1] == A[1] and B[0] > A[0]:<br />             alpha = 0  #<br />         elif B[1] == A[1] and B[0] < A[0]:<br />             alpha = 180<br />         elif B[0] == A[0] and B[1] > A[1]:<br />             alpha = 90<br />         elif B[0] == A[0] and B[1] < A[1]:<br />             alpha = 270<br />         else:<br />             pass<br />     if UpDown == 2:  # 象限号<br />         if B[0] == A[0] and B[1] > A[1]:<br />             alpha = 0<br />         elif B[0] == A[0] and B[1] < A[1]:<br />             alpha = 180<br />         elif B[1] == A[1] and B[0] > A[0]:<br />             alpha = 270<br />         elif B[1] == A[1] and B[0] < A[0]:<br />             alpha = 90<br />         else:<br />             pass<br />     if UpDown == 3:  # 象限号<br />         if B[1] == A[1] and B[0] < A[0]:<br />             alpha = 0<br />         elif B[1] == A[1] and B[0] > A[0]:<br />             alpha = 180<br />         elif B[0] == A[0] and B[1] < A[1]:<br />             alpha = 90<br />         elif B[0] == A[0] and B[1] > A[1]:<br />             alpha = 270<br />         else:<br />             pass<br />     if UpDown == 4:  # 象限号<br />         if B[0] == A[0] and B[1] < A[1]:<br />             alpha = 0<br />         elif B[0] == A[0] and B[1] > A[1]:<br />             alpha = 180<br />         elif B[1] == A[1] and B[0] < A[0]:<br />             alpha = 270<br />         elif B[1] == A[1] and B[0] > A[0]:<br />             alpha = 90<br />         else:<br />             pass<br />     return alpha<br /> <br /> <br /> <br /> def Dist(A, B):<br />     return (B[1] - A[1])**2 + (B[0] - A[0])**2<br /> <br /> <br /> def NextPoi(UpDown, Poi, PList):<br />     Dire = 360 + 10<br />     Unit = 1<br />     for id, p in enumerate(PList):<br />         if Direc(UpDown, Poi, p) <= Dire and Dist(Poi, p) == Unit:<br />             Dire = Direc(UpDown, Poi, p)<br />             dis = Dist(Poi, p)<br />             Num = id<br />     return PList[Num]<br /> <br /> <br /> def plot2D(x, y, x_out, y_out):<br />     plt.scatter(x, y)  # 画散点图<br />     plt.plot(x_out, y_out)  # 画连线图<br />     plt.show()<br /> <br /> if __name__ == '__main__':<br />     Points = Gen.GenList()<br />     XMaxPoi, Num = MP.XMaxPoint(Points)<br />     YMaxPoi, Num = MP.YMaxPoint(Points)<br />     XMinPoi, Num = MP.XMinPoint(Points)<br />     YMinPoi, Num = MP.YMinPoint(Points)<br />     PPs = Points.copy()<br />     PPs.remove(XMaxPoi)<br />     Poi_A = XMaxPoi<br />     RList = [XMaxPoi]<br />     Poi_B = []<br />     while Poi_B != YMaxPoi:<br />         Poi_B = NextPoi(1, Poi_A, PPs)<br />         RList.append(Poi_B)<br />         PPs.remove(Poi_B)<br />         Poi_A = Poi_B<br />     Poi_A = YMaxPoi<br />     Poi_B = []<br />     while Poi_B != XMinPoi:<br />         Poi_B = NextPoi(2, Poi_A, PPs)<br />         RList.append(Poi_B)<br />         PPs.remove(Poi_B)<br />         Poi_A = Poi_B<br />     Poi_A = XMinPoi<br />     Poi_B = []<br />     while Poi_B != YMinPoi:<br />         Poi_B = NextPoi(3, Poi_A, PPs)<br />         RList.append(Poi_B)<br />         PPs.remove(Poi_B)<br />         Poi_A = Poi_B<br />     Poi_A = YMinPoi<br />     Poi_B = []<br />     PPs = [XMaxPoi] + PPs  # 加回第一点<br />     while Poi_B != XMaxPoi:<br />         Poi_B = NextPoi(4, Poi_A, PPs)<br />         RList.append(Poi_B)<br />         PPs.remove(Poi_B)<br />         Poi_A = Poi_B<br />     print(RList)<br /> <br />     XList = [i[0] for i in Points]<br />     YList = [i[1] for i in Points]<br />     X_OutList = [i[0] for i in RList]<br />     Y_OutList = [i[1] for i in RList]<br />     plot2D(XList, YList, X_OutList, Y_OutList)<br />

1. # MaxP.py
   
2. def XMaxPoint(PList):
   
3.     Xmax = PList[0][0]
   
4.     Ymax = PList[0][1]
   
5.     for id, poi in enumerate(PList):
   
6.         if poi[0] > Xmax:
   
7.             Xmax = poi[0]
   
8.             Ymax = poi[1]
   
9.             Num = id
   
10.         if poi[0] >= Xmax and poi[1] > Ymax:
    
11.             Xmax = poi[0]
    
12.             Num = id
    
13.     return PList[Num], Num
    
14. 
    
15. 
    
16. def YMaxPoint(PList):
    
17.     Xmax = PList[0][0]
    
18.     Ymax = PList[0][1]
    
19.     for id, poi in enumerate(PList):
    
20.         if poi[1] > Ymax:
    
21.             Ymax = poi[1]
    
22.             Xmax = poi[0]
    
23.             Num = id
    
24.         if poi[1] >= Ymax and poi[0] < Xmax:
    
25.             Ymax = poi[1]
    
26.             Num = id
    
27.     return PList[Num], Num
    
28. 
    
29. 
    
30. def XMinPoint(PList):
    
31.     Xmin = PList[0][0]
    
32.     Ymin = PList[0][1]
    
33.     for id, poi in enumerate(PList):
    
34.         if poi[0] < Xmin:
    
35.             Xmin = poi[0]
    
36.             Ymin = poi[1]
    
37.             Num = id
    
38.         if poi[0] <= Xmin and poi[1] < Ymin:
    
39.             Xmin = poi[0]
    
40.             Num = id
    
41.     return PList[Num], Num
    
42. 
    
43. 
    
44. def YMinPoint(PList):
    
45.     Xmin = PList[0][0]
    
46.     Ymin = PList[0][1]
    
47.     for id, poi in enumerate(PList):
    
48.         if poi[1] < Ymin:
    
49.             Ymin = poi[1]
    
50.             Xmin = poi[0]
    
51.             Num = id
    
52.         if poi[1] <= Ymin and poi[0] > Xmin:
    
53.             Ymin = poi[1]
    
54.             Num = id
    
55.     return PList[Num], Num
    

# MaxP.py<br /> def XMaxPoint(PList):<br />     Xmax = PList[0][0]<br />     Ymax = PList[0][1]<br />     for id, poi in enumerate(PList):<br />         if poi[0] > Xmax:<br />             Xmax = poi[0]<br />             Ymax = poi[1]<br />             Num = id<br />         if poi[0] >= Xmax and poi[1] > Ymax:<br />             Xmax = poi[0]<br />             Num = id<br />     return PList[Num], Num<br /> <br /> <br /> def YMaxPoint(PList):<br />     Xmax = PList[0][0]<br />     Ymax = PList[0][1]<br />     for id, poi in enumerate(PList):<br />         if poi[1] > Ymax:<br />             Ymax = poi[1]<br />             Xmax = poi[0]<br />             Num = id<br />         if poi[1] >= Ymax and poi[0] < Xmax:<br />             Ymax = poi[1]<br />             Num = id<br />     return PList[Num], Num<br /> <br /> <br /> def XMinPoint(PList):<br />     Xmin = PList[0][0]<br />     Ymin = PList[0][1]<br />     for id, poi in enumerate(PList):<br />         if poi[0] < Xmin:<br />             Xmin = poi[0]<br />             Ymin = poi[1]<br />             Num = id<br />         if poi[0] <= Xmin and poi[1] < Ymin:<br />             Xmin = poi[0]<br />             Num = id<br />     return PList[Num], Num<br /> <br /> <br /> def YMinPoint(PList):<br />     Xmin = PList[0][0]<br />     Ymin = PList[0][1]<br />     for id, poi in enumerate(PList):<br />         if poi[1] < Ymin:<br />             Ymin = poi[1]<br />             Xmin = poi[0]<br />             Num = id<br />         if poi[1] <= Ymin and poi[0] > Xmin:<br />             Ymin = poi[1]<br />             Num = id<br />     return PList[Num], Num<br />

1. # GenList.py
   
2. def GenList():
   
3.     Points = []
   
4.     for i in range(6, 9):
   
5.         for j in range(14, 17):
   
6.             Points.append([i, j])
   
7.     for i in range(4, 12):
   
8.         for j in range(12, 14):
   
9.             Points.append([i, j])
   
10.     for i in range(5, 10):
    
11.         for j in range(10, 12):
    
12.             Points.append([i, j])
    
13.     for i in range(2, 13):
    
14.         for j in range(6, 10):
    
15.             Points.append([i, j])
    
16.     for i in range(4, 13):
    
17.         for j in range(4, 6):
    
18.             Points.append([i, j])
    
19.     for i in range(3, 8):
    
20.         for j in range(0, 4):
    
21.             Points.append([i, j])
    
22.     for i in range(9, 12):
    
23.         for j in range(2, 5):
    
24.             Points.append([i, j])
    
25.     return Points
    

# GenList.py<br /> def GenList():<br />     Points = []<br />     for i in range(6, 9):<br />         for j in range(14, 17):<br />             Points.append([i, j])<br />     for i in range(4, 12):<br />         for j in range(12, 14):<br />             Points.append([i, j])<br />     for i in range(5, 10):<br />         for j in range(10, 12):<br />             Points.append([i, j])<br />     for i in range(2, 13):<br />         for j in range(6, 10):<br />             Points.append([i, j])<br />     for i in range(4, 13):<br />         for j in range(4, 6):<br />             Points.append([i, j])<br />     for i in range(3, 8):<br />         for j in range(0, 4):<br />             Points.append([i, j])<br />     for i in range(9, 12):<br />         for j in range(2, 5):<br />             Points.append([i, j])<br />     return Points<br />

landsat99

此算法的右手航向角 取为正交判断（针对正交栅格散点集合）。

对任意散点集合，需设定另外的约束条件。

mahuan1279

删掉的点都是以该点为中心八方向都有点的点。

landsat99

mahuan1279 发表于 2023-2-7 08:44
删掉的点都是以该点为中心八方向都有点的点。

您这个思路很棒！确实可以解决很多阴角歧点、回头线等问题。
目前我用的四分区的办法（集合4个极值分区）分别制定航向规则，来解决阴角歧点、回头线问题。

您的方案， 估计算法上会更佳。赞一个

qjchen

感觉题目定义可能会导致一些无解的情况
比如矩形格点只有3个点的情况，或者有单独几点伸出的情况，也可能如楼主说的，需要限制条件
https://karobben.github.io/2022/03/07/Python/point_outline/
此处也有一个类似的问题，用了α来限制

landsat99

qjchen 发表于 2023-2-7 09:54
感觉题目定义可能会导致一些无解的情况
比如矩形格点只有3个点的情况，或者有单独几点伸出的情况，也可能 ...

同意您的观点。

对有单点的分支凸起，此算法会中断。单点凸起的轮廓线，是需要单独定义。

如果轮廓线定义为： 散点集合的最大面积轮廓线，则算法将会大大简化。

mahuan1279

landsat99 发表于 2023-2-7 11:27
同意您的观点。

对有单点的分支凸起，此算法会中断。单点凸起的轮廓线，是需要单独定义。

那不就是求凸包么？

landsat99

是的。凸包算法将大为简化算法。即便手写算法，也无需分区制定规则。尤其Convex的成熟包很多，稳定 效率优秀。


Scipy的ConvexHull实例

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. from scipy.spatial import ConvexHull
   
4. 
   
5. x = np.random.normal(0, 1, 100)
   
6. y = np.random.normal(0, 1, 100)
   
7. xy = np.hstack((x[:, np.newaxis], y[:, np.newaxis]))
   
8. hull = ConvexHull(xy)
   
9. plt.scatter(x, y)
   
10. plt.plot(x[hull.vertices], y[hull.vertices])
    
11. plt.plot(np.append(x[hull.vertices], x[hull.vertices][0]), np.append(y[hull.vertices], y[hull.vertices][0]))
    
12. plt.show()
    

import numpy as np<br /> import matplotlib.pyplot as plt<br /> from scipy.spatial import ConvexHull<br /> <br /> x = np.random.normal(0, 1, 100)<br /> y = np.random.normal(0, 1, 100)<br /> xy = np.hstack((x[:, np.newaxis], y[:, np.newaxis]))<br /> hull = ConvexHull(xy)<br /> plt.scatter(x, y)<br /> plt.plot(x[hull.vertices], y[hull.vertices])<br /> plt.plot(np.append(x[hull.vertices], x[hull.vertices][0]), np.append(y[hull.vertices], y[hull.vertices][0]))<br /> plt.show()<br />

landsat99

// 2023.02.14
  
改为正则矢量法，可得到较优方案。普式性、简洁性均较大改善。
点集/线集序列应用，在较多情况下 正则矢量法会优于航向角法。
改进过程 不定期更新中.  https://landsat99.github.io/GridContour.html

同时欢迎您提供 数据测试

1. #
   
2. # 不支持 “单线分支”的图形集合
   
3. #
   
4. import matplotlib.pyplot as plt
   
5. import Data1 as Data
   
6. 
   
7. def XYMaxPoint(PList):
   
8.     Xmax = max(i[0] for i in PList)
   
9.     Xlist = []
   
10.     for poi in PList:
    
11.         if poi[0] == Xmax:
    
12.             Xlist.append(poi)
    
13.     Ymax = max(i[1] for i in Xlist)
    
14.     return [Xmax, Ymax]
    
15. 
    
16. def NextPoi(P, PList):
    
17.     Unit = 1.0
    
18.     Poi_B = []
    
19.     TList = []
    
20.     Grid_4 = [[P[0] + Unit, P[1]], [P[0] - Unit, P[1]],
    
21.               [P[0], P[1] + Unit], [P[0], P[1] - Unit]]
    
22.     for i in Grid_4:
    
23.         if i in PList:
    
24.             TList.append(i)
    
25.     for A in TList:
    
26.         k = 1
    
27.         for B in TList:
    
28.             Wise = (A[0] - P[0]) * (B[1] - P[1]) - \
    
29.                 (A[1] - P[1]) * (B[0] - P[0])
    
30.             if Wise < 0:
    
31.                 k = 0
    
32.         if k == 1:
    
33.             Poi_B = A
    
34.     return Poi_B
    
35. 
    
36. def plot2D(x, y, x_out, y_out):
    
37.     plt.scatter(x, y)
    
38.     plt.plot(x_out, y_out)
    
39.     plt.show()
    
40. 
    
41. if __name__ == '__main__':
    
42.     PPs = Data.Points.copy()
    
43.     PoiStart = XYMaxPoint(Data.Points)
    
44.     Poi_A = PoiStart
    
45.     RList = []
    
46.     Poi_B = []
    
47.     n = 0
    
48.     while Poi_B != PoiStart:
    
49.         Poi_B = NextPoi(Poi_A, PPs)
    
50.         RList.append(Poi_A)
    
51.         PPs.remove(Poi_A)
    
52.         Poi_A = Poi_B
    
53.         n = n + 1
    
54.         if n == 2:
    
55.             PPs.append(PoiStart)
    
56.     RList.append(PoiStart)
    
57.     XList = [i[0] for i in Data.Points]
    
58.     YList = [i[1] for i in Data.Points]
    
59.     X_OutList = [i[0] for i in RList]
    
60.     Y_OutList = [i[1] for i in RList]
    
61.     plot2D(XList, YList, X_OutList, Y_OutList)
    

#<br /> # 不支持 “单线分支”的图形集合<br /> #<br /> import matplotlib.pyplot as plt<br /> import Data1 as Data<br /> <br /> def XYMaxPoint(PList):<br />     Xmax = max(i[0] for i in PList)<br />     Xlist = []<br />     for poi in PList:<br />         if poi[0] == Xmax:<br />             Xlist.append(poi)<br />     Ymax = max(i[1] for i in Xlist)<br />     return [Xmax, Ymax]<br /> <br /> def NextPoi(P, PList):<br />     Unit = 1.0<br />     Poi_B = []<br />     TList = []<br />     Grid_4 = [[P[0] + Unit, P[1]], [P[0] - Unit, P[1]],<br />               [P[0], P[1] + Unit], [P[0], P[1] - Unit]]<br />     for i in Grid_4:<br />         if i in PList:<br />             TList.append(i)<br />     for A in TList:<br />         k = 1<br />         for B in TList:<br />             Wise = (A[0] - P[0]) * (B[1] - P[1]) - \<br />                 (A[1] - P[1]) * (B[0] - P[0])<br />             if Wise < 0:<br />                 k = 0<br />         if k == 1:<br />             Poi_B = A<br />     return Poi_B<br /> <br /> def plot2D(x, y, x_out, y_out):<br />     plt.scatter(x, y)<br />     plt.plot(x_out, y_out)<br />     plt.show()<br /> <br /> if __name__ == '__main__':<br />     PPs = Data.Points.copy()<br />     PoiStart = XYMaxPoint(Data.Points)<br />     Poi_A = PoiStart<br />     RList = []<br />     Poi_B = []<br />     n = 0<br />     while Poi_B != PoiStart:<br />         Poi_B = NextPoi(Poi_A, PPs)<br />         RList.append(Poi_A)<br />         PPs.remove(Poi_A)<br />         Poi_A = Poi_B<br />         n = n + 1<br />         if n == 2:<br />             PPs.append(PoiStart)<br />     RList.append(PoiStart)<br />     XList = [i[0] for i in Data.Points]<br />     YList = [i[1] for i in Data.Points]<br />     X_OutList = [i[0] for i in RList]<br />     Y_OutList = [i[1] for i in RList]<br />     plot2D(XList, YList, X_OutList, Y_OutList)<br />