翻譯自:https://www.codeproject.com/Articles/114797/Polyline-Simplification#headingDPN
參考：https://www.codeproject.com/Articles/114797/Polyline-Simplification
整個記錄關于演算法的部分是翻譯來的，原作者實作的語言是C++，但是我看不懂這類代碼，于是自己用Python實作了一遍，其中可能會有錯誤的地方，歡迎指出來讓我改正，
In [1]:
from shapefile import Reader, Writer import numpy as np import math import matplotlib.pyplot as plt
In [26]:
# 該shapefile檔案內含一條線段，坐標系為：3857
shp = Reader('./vector/lyr_line.shp', encoding='utf-8')
In [3]:
feat = shp.shape(0) line = np.array(feat.points)
In [4]:
# 保存為shapefile檔案
def shpGen(fileName, pts):
w = Writer(fileName, encoding='utf-8', shapeType=3)
w.fields = shp.fields
w.line([pts])
w.record("測驗", 3)
w.close()
In [5]:
# 點到線距離
def point2lineDistance(x0,y0,x1,y1,x2,y2):
a =x2 -x1
b = y2 - y1
c = b*x1- a*y1
return abs(a*y0 - b*x0 + c) / math.sqrt(a*a + b*b)
# return abs((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1)) / math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
# 點到點距離
def point2pointDistance(x1,y1,x2,y2):
dx = x2-x1
dy = y2-y1
return math.sqrt(dx*dx+dy*dy)
# 繪制
def plotLine(pts, simplifyLst, newlabel:str):
plt.figure(figsize=(15, 6))
plt.plot(pts[:,0], pts[:,1], 'r-*', linewidth=0.5, mfc='red', mec='red', label='origin')
plt.plot(simplifyLst[:,0], simplifyLst[:,1], 'b-s', linewidth=0.5, mfc='blue', mec='blue', label=newlabel)
plt.legend(loc='upper left')
plt.show()
Radial Distance?
該演算法是一個時間復雜度為${O(n)}$的暴力演算法，在連續的頂點中，與當前某一個key（簡化程序中標記保留的頂點）的距離太近時會被去掉，
起止頂點一般會保留作為簡化線的一部分，故標記為key，從第一個key（第一個頂點）開始，演算法將遍歷整條線段，所有連續頂點中，移除那些距離不超過容忍值的頂點，而第一個超過容忍值的頂點將標記為下一個key，從新key開始，演算法重新遍歷余下的點并重復上述程序，直到最后一個點，
In [6]:
def NearestSimplify(tolorence:float, pts:list) -> list:
start = 0
validFlags = np.zeros((len(pts), ), dtype="bool")
for i in range(1, len(pts)):
dist = math.sqrt(math.pow(pts[i][0]- pts[start][0], 2) +
math.pow(pts[i][1] - pts[start][1], 2))
if dist > tolorence:
validFlags[i] = True
start = i
validFlags[0] = True
return validFlags
In [7]:
flags = NearestSimplify(20, feat.points) plotLine(line, line[flags], "Radial Distance")
In [8]:
flags = NearestSimplify(10, feat.points) plotLine(line, line[flags], "Radial Distance")
Perpendicular Distance?
臨近點演算法將點-點距離作為誤差判據，而垂直距離簡化則是將點-線段的距離作為誤差判據，對于每個頂點Vi，需要計算它與線段[Vi-1， Vi+1]的垂直距離，距離小于給定容忍值的點將被移除，程序如下圖所示：
最開始處理前三個頂點，計算第二個頂點的垂直距離，與容忍值作比較，大于則標記為key，然后演算法往前移動開始處理下一個前三個頂點，第二次距離小于容忍值，移除中間頂點，重復直到處理所有頂點，
注意：對于每個被移除的頂點Vi，下一個可能被移除的候選頂點是Vi+2，這意味著原始多段線的點數最多只能減少50%，為了獲得更高的頂點減少率，需要多次遍歷，
In [27]:
def PerpendicularSimplify(tolorence:float, pts:list) -> list:
sp = 0
ep = 2
validFlags = np.zeros((len(pts), ), dtype="bool")
for i in range(1, len(pts) - 1):
verDist = point2lineDistance(
pts[i][0], pts[i][1],
pts[sp][0], pts[sp][1],
pts[ep][0], pts[ep][1])
if verDist > tolorence:
validFlags[i] = True
sp = i
ep = i + 2
validFlags[0] = validFlags[-1] = True
return validFlags
In [10]:
perLine = PerpendicularSimplify(20, feat.points) plotLine(line, line[perLine], "Perpendicular Distance")
In [11]:
perLine = PerpendicularSimplify(10, feat.points) plotLine(line, line[perLine], "Perpendicular Distance")
Reumann-Witkam?
該演算法使用點-線距離作距離判斷，首先定義一條線段，起止點為原始線段的前兩個頂點，對于每一個連續的頂點vi，計算它與前述線段的垂直距離，當距離超過容忍值時，將頂點vi - 1標記為key，頂點vi和vi + 1用來定義新的線段，演算法重復直到最后一個點，下圖為演算法演示：
紅色的條帶由指定的容忍值生成，而線段則由原始線的前兩個頂點生成，第三個頂點不在紅色條帶內，意味著該頂點是一個key，通過第二，三頂點定義新的條帶，位于條帶內的最后一個頂點也標記為key，而其他點則移除，演算法重復直到新吊帶包含原始線段的最后一個頂點，
In [12]:
def ReuWitkamSimplify(tolorence:float, pts:list) -> list:
sp = 0
ep = 1
validFlags = np.zeros((len(pts), ), dtype="bool")
for i in range(2, len(pts)):
verDist = point2lineDistance(
pts[i][0], pts[i][1],
pts[sp][0], pts[sp][1],
pts[ep][0], pts[ep][1])
if verDist > tolorence:
validFlags[i-1] = True
sp = ep
ep = i
validFlags[0] = validFlags[-1] = True
return validFlags
In [13]:
rLine = ReuWitkamSimplify(20, feat.points) plotLine(line, line[rLine], "Reumann-Witkam")
In [14]:
rLine = ReuWitkamSimplify(10, feat.points) plotLine(line, line[rLine], "Reumann-Witkam")
Opheim?
Ophein 使用 最小和最大距離容忍值來約束搜索區域，對于連續的頂點vi，它與當前key（通常是v0的徑向距離可以被計算出來，在容忍值范圍內的最后一個頂點被用來定義一條射線R(vkey, vi)，如果這樣的頂點不存在，那么射線將被定義為R(vkey, vkey + 1)，計算每一個連續的頂點vj（在vi后邊）與射線R間的垂直距離， 當該距離超過最小容忍值或者當vj與vkey間的輻射距離大于最大容忍值時，在vj-1位置上的頂點，將被認為是一個新key，
值得注意的是，關于如何定義射線${R(v_{key},v_i)}$ 去控制搜索區域的方法存在多個，可能的方法有如下幾種：
- Reumann-Witkam 的方法 : ${i}$ = ${key}$ + 1.
- 外部第一個點：${key}$ < ${i}$ 且${v_i}$ 是第一個掉落在最小輻射距離外部的點.
- 內部最后一個點：${key}$ < ${i}$ 且 ${v_i}$ 是掉落在最小輻射距離內的最后一個點；若該點不存在，則使用方法1.
In [31]:
def opheimSimplify(minTolorence:int, maxTolorence:int, pts) -> list:
vKey = 0
validPtsLst = np.zeros((len(pts),), dtype='bool')
for i in range(vKey, len(pts)):
if point2pointDistance(pts[vKey][0], pts[vKey][1], pts[i][0], pts[i][1]) > minTolorence:
vi = i
break
for j in range(vi + 1, len(pts)):
vj2rayDist = point2lineDistance(
pts[j][0], pts[j][1],
pts[vKey][0], pts[vKey][1],
pts[vi][0], pts[vi][1])
radialDist = point2pointDistance(pts[vKey][0], pts[vKey][1], pts[j][0], pts[j][1])
if vj2rayDist <= minTolorence or radialDist > maxTolorence:
vKey = j - 1
validPtsLst[j-1] = True
else:
continue
validPtsLst[0] = validPtsLst[-1]= True
return np.array(validPtsLst)
In [16]:
opheim = opheimSimplify(2, 20, np.array(feat.points)) plotLine(line, line[opheim], "opheim")
In [17]:
opheim = opheimSimplify(2, 10, np.array(feat.points)) plotLine(line, line[opheim], "opheim")
Lang?
Lang 簡化演算法定義一個固定大小的搜索區域，它的第一個和最后一個點來自線段的同一個片段——用來計算處于中間的頂點的垂直距離，若任何計算出的距離大于指定的容忍值，那么搜索區域將通過排除其最后一個頂點進行收縮，這個程序會一直重復直到所有垂直距離小于特定的容忍值，或者中間頂點全部排除，所有中間節點移除后，新的搜索區域將以舊搜索區域的最后一個點作為起點進行構建，
搜索區域使用前四個頂點來生成，計算每一個中間頂點到片段S(v0, v4)的垂直距離，只要至少存在一個頂點的垂直距離超過容忍值，搜索區域便簡化一個頂點，再重復計算中間頂點到片段S(v0, v3)，所有的中間頂點到片段S的垂直距離都位于容忍值內，搜索區域的最后一個頂點v3作為一個新key同時它也用來定義一個新的區域（片段）S(v3, v7)，
- 第一個點作為vKey,創建片段S(v0, v4)
- 求中間頂點到片段的垂直距離
- 片段依次排除最后一個頂點（當距離超過容忍值）
- 片段內無頂點的垂直距離超過容忍值，此時最后一個點作為新key，同時構建新片段
- 重復1-4
In [18]:
def langSimplify(tolorence, pts, steps=0) -> list:
vKey = 0
end = 2
validLst = np.zeros((len(pts),), dtype='bool')
if(steps == 0):
maxSteps = len(pts) - 1
else:
maxSteps = steps
while vKey < maxSteps and end < maxSteps:
isBigger = False
# print("segment:(%s，%s)" % (vKey, end))
for i in range(vKey + 1, end):
perDist = point2lineDistance(
pts[i][0], pts[i][1],
pts[vKey][0], pts[vKey][1],
pts[end][0], pts[end][1])
if perDist > tolorence:
isBigger = True
end -= 1
break
else:
continue
if not isBigger:
vKey = end
end += 2
validLst[vKey] = True
validLst[0] = validLst[-1] = True
return validLst
In [19]:
langLine = langSimplify(20, feat.points, steps=0) plotLine(line, line[langLine], "Lang")
In [20]:
langLine = langSimplify(10, feat.points, steps=0) plotLine(line, line[langLine], "Lang")
Douglas-Peucker?
Douglas-Peucker演算法使用點-邊距離作為誤差衡量標準，該演算法從連接原始Polyline的第一個和最后一個頂點的邊開始，計算所有中間頂點到邊的距離，距離該邊最遠的頂點，如果其距離大于指定的公差，將被標記為Key并添加到簡化結果中，這個程序將對當前簡化中的每條邊進行遞回，直到原始Polyline的所有頂點與當前考察的邊的距離都在允許誤差范圍內，該程序如下圖所示:
初始時，簡化結果只有一條邊（起點-終點），第一步中，將第四個頂點標記為一個Key，并相應地加入到簡化結果中；第二步，處理當前結果中的第一條邊，到該邊的最大頂點距離低于容差，因此不添加任何點；在第三步中，當前簡化的第二個邊找到了一個Key（點到邊的距離大于容差），這條邊在這個Key處被分割，這個新的Key添加到簡化結果中，這個程序一直繼續下去，直到找不到新的Key為止，注意，在每個步驟中，只處理當前簡化結果的一條邊，
In [22]:
def douglasPeuckerSimplify(tolorence, pts) -> list:
start = 0
end = len(pts) - 1
def findPt(start, end):
curPt = 0
maxDist = 0
for i in range(start + 1, end):
perpenDist = point2lineDistance(
pts[i][0], pts[i][1],
pts[start][0], pts[start][1],
pts[end][0], pts[end][1])
if perpenDist > maxDist:
maxDist = perpenDist
curPt = i
if maxDist > tolorence:
preKey = findPt(start, curPt)
postKey = findPt(curPt, end)
return [] + preKey + postKey
else:
return [pts[end]]
return [pts[0]] + findPt(start, end)
In [23]:
res = douglasPeuckerSimplify(20, feat.points) plotLine(line, np.array(res), "Douglas-Peucker")
In [24]:
res = douglasPeuckerSimplify(10, feat.points) plotLine(line, np.array(res), "Douglas-Peucker")
In [25]:
res = douglasPeuckerSimplify(5, feat.points) plotLine(line, np.array(res), "Douglas-Peucker")
In [ ]:
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/539738.html
標籤：其他