#はじめに
長野県の地図を見ていると気になるものがありました。
この町、めっちゃ細いところがありますね。
これは長さを測ってみるしかないですね
どれくらい細いのか、調べてみました。
手順
立科町の画像を用意し、グレースケールにします。
img = cv2.imread("tateshina.png",3)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
thresh,bin_img = cv2.threshold(gray,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
plt.imshow(bin_img,cmap='gray')
よくある、モノクロ立科町が出来ました。
2つ以上に分割する
この状態から任意の二点を取得して、最短距離を求めても望みのものは出ないので
二分割しましょう。
itter=1
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))
while True :
erosion = cv2.erode(bin_img,kernel,iterations = itter)
contours_erosion, _ = cv2.findContours(erosion, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
if len(contours_erosion) > 1:
break
itter += 1
plt.imshow(erosion)
plt.show()
kernelを用意して、少しずつerodeをかけていき、二つに分割されたところで切ります。
2つのブロブの最短距離を探す
分かりにくいですが、切科町になりました。
この時の２つのブロブの最短距離を求めましょう。
def get_distance(x1, y1, x2, y2):
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
distance_list = []
ch_list = []
for i in (contours_erosion[0]):
for j in (contours_erosion[1]):
distance_list.append(get_distance(*i[0],*j[0]))
ch_list.append([i,j])
ch_list[np.argmin(distance_list)]
x1 = ch_list[np.argmin(distance_list)][0][0][0]
y1 = ch_list[np.argmin(distance_list)][0][0][1]
x2 = ch_list[np.argmin(distance_list)][1][0][0]
y2 = ch_list[np.argmin(distance_list)][1][0][1]
fig, ax = plt.subplots(figsize=(8, 8))
fig = 0
cv2.line(erosion, (x1, y1), (x2, y2), (128, 128, 255), thickness=1)
plt.imshow(erosion,cmap="gray")
最短となる線分が引けました。
何がしたかったと申しますと、この線を通る任意の二点を結ぶ線分がもっとも細い箇所となるので
それらを求めたいと思います。
##エッジを取得する
エッジを取得しましょう。
edge = cv2.morphologyEx(bin_img, cv2.MORPH_GRADIENT, kernel)
plt.imshow(edge,cmap="gray")
サイバー立科町が取れました。これで輪郭の座標が取れますね。
##分割した線に近い点を集める
任意の二点で計算してもいいですが、計算量を減らす工夫をしましょう。
今回は先ほど引いた２点の重心からの距離が近い点を集めましょう。
row_list = np.where(edge == 255)[0]
col_list = np.where(edge == 255)[1]
near_list = []
for i in range(len(col_list)):
u = np.array([(x2 + x1)/2, (y2 + y1)/2])
v = np.array([col_list[i], row_list[i]])
L = np.linalg.norm(u-v)
if L > itter*2:
continue
near_list.append(v)
これで、多少点は減らせたはずです。
##最短の線を通る2点を探す
2つの線分が、交わっているか否かは
def intersect(p1, p2, p3, p4):
tc1 = (p1[0] - p2[0]) * (p3[1] - p1[1]) + (p1[1] - p2[1]) * (p1[0] - p3[0])
tc2 = (p1[0] - p2[0]) * (p4[1] - p1[1]) + (p1[1] - p2[1]) * (p1[0] - p4[0])
td1 = (p3[0] - p4[0]) * (p1[1] - p3[1]) + (p3[1] - p4[1]) * (p3[0] - p1[0])
td2 = (p3[0] - p4[0]) * (p2[1] - p3[1]) + (p3[1] - p4[1]) * (p3[0] - p2[0])
return tc1*tc2<0 and td1*td2<0
shortest_list = []
from itertools import combinations
comb = list(combinations(near_list,2))
for tmp in comb:
if (intersect([x1,y1],[x2,y2],tmp[0],tmp[1])) == 1:
L = np.linalg.norm(np.array(tmp[0])-np.array(tmp[1]))
shortest_list.append(L)
else:
shortest_list.append(50000)
もし交わっているなら長さを保存し、そうでないなら、あり得ない距離を入れています。
short = np.argmin(shortest_list)
f1x,f1y,f2x,f2y=*comb[short][0],*comb[short][1]
cv2.line(erosion,(f1x,f1y) ,(f2x,f2y), (192,0, 192), thickness=1)
res = cv2.addWeighted(gray, 0.5, erosion, 0.5, 0)
fig, ax = plt.subplots(figsize=(32, 32))
plt.imshow(res)
小さくてわかりにくいですが、最短の場所が分かりました。
これで、実際に立科町で測定する必要がなくなりました。
総括
画像処理は単純な処理が多い分、複雑なことが出来て面白いですね。
縮小や拡大の方法に不安があるので、ほかに良い方法がないかは考えてみたいと思います。
## 参考