pythonとopencvを使って画像処理を勉強していきます。
前回
python+opencvで画像処理の勉強4 周波数領域におけるフィルタリング
周波数フィルタリングについて学びました。
python+opencvで画像処理の勉強5 幾何学的変換
今回は2値画像処理を勉強していきます。
1つ目の使用画像は旭川動物園のペンギンの散歩の写真です。
2値化
2値化の意味
例えばグレースケール画像などで、中間値をなくし、白または黒の2値の画像に変換することを2値化とよびます。
2値画像は、例えばある値（しきい値）以上の画素値を1、それ未満の画像を0に変換することで得られます。
cv2.thresholdで2値画像を作成します。
fig, ax = plt.subplots(2, 4, figsize=(18, 6))
ax[0][0].imshow(img_gray, 'gray')
ax[0][0].set_xticks([])
ax[0][0].set_yticks([])
ax[1][0].hist(img_gray.flatten(), 256, [0,256]);
threshholds = [160, 100, 40]
for i in range((len(threshholds))):
ret, dst = cv2.threshold(img_gray, threshholds[i], 255, cv2.THRESH_BINARY)
ax[0][i+1].imshow(dst, 'gray')
ax[0][i+1].set_title('thresh={}'.format(threshholds[i]))
ax[0][i+1].set_xticks([])
ax[0][i+1].set_yticks([])
ax[1][i+1].hist(dst.flatten(), 256, [0,256]);
p-タイル法
p-タイル法とは、対象物の占める領域の画素数が、あらかじめ予測できる場合に、予測された画素数に応じてしきい値を決める手法です。
画素値の低いところから頻度値を積算し、予測された画素数を超えたときの画素値をしきい値とします。
対象画像に対する事前の知識が必要であり、未知の画像を処理する場合には適していません。
モード法
グレースケール画像は2つの山を持つことが多い。モード法は山と山の谷の底をしきい値とします。
画像に十分な画素数がないと山がはっきりせず、安定に底を検出することができません。
判別分析法
黒白どちらかのクラスに属する画素は、あるしきい値tの左右に分かれて分布します。
その分布の分離度が大きくなるようにしきい値を決める方法を**判別分析法（大津の方法）**と呼びます。
分離度は、クラス間分散$\sigma_{b}^2$とクラス内分散$\sigma_{w}^2$の日$\sigma_{b}^2/\sigma_{w}^2$で定義します。
全体の画素値の平均と分散を$m_{t}$と$\sigma_{t}^2$、黒画素クラスを$m_{1}$と$\sigma_{1}^2$、黒画素クラスを$m_{2}$と$\sigma_{2}^2$とし、
$\omega_{1}$と$\omega_{2}$をそれぞれ黒画素クラスと白画素クラスに属する画素の数とします。
クラス内分散は、
$$
\sigma_{w}^2=\frac{\omega_{1}\sigma{1}^2+\omega_{2}\sigma{2}^2}{\omega_{1}+\omega_{2}}
$$
クラス間分散は、以下の式で表されます。
$$
\sigma_{w}^2=\frac{\omega_{1}(m_{1}-m_{t})^2+\omega_{2}(m_{2}-m_{t})^2}{\omega_{1}+\omega_{2}}=
\frac{\omega_{1}\omega_{2}(m_{1}-m_{2})^2}{(\omega_{1}+\omega_{2})^2}
$$
全分散$\sigma_{t}^2$とクラス間分散$\sigma_{b}^2$、クラス内分散$\sigma_{w}^2$には、以下のような関係があります。
$$
\sigma_{t}^2=\sigma_{b}^2+\sigma_{w}^2
$$
したがって、以下の式が成り立ち、しきい値tに関係なく全分散$\sigma_{t}^2$は一定のため、クラス間分散$\sigma_{b}^2$が最大になるようにtを決めればよいことがわかります。
$$
\frac{\sigma_{b}^2}{\sigma_{w}^2}=\frac{\sigma_{b}^2}{\sigma_{t}^2-\sigma_{b}^2}
$$
oprncvに2値化処理の手法がいくつか用意されているので見ていきます。
methods = [cv2.THRESH_BINARY, cv2.THRESH_OTSU, cv2.THRESH_TRIANGLE,
cv2.THRESH_BINARY_INV, cv2.THRESH_TOZERO, cv2.THRESH_TRUNC,
cv2.THRESH_TOZERO_INV]
method_name = ['THRESH_BINARY', 'THRESH_OTSU', 'THRESH_TRIANGLE',
'THRESH_BINARY_INV', 'THRESH_TOZERO', 'THRESH_TRUNC',
'THRESH_TOZERO_INV']
fig, ax = plt.subplots(2, 4, figsize=(18, 5), subplot_kw=({"xticks":(), "yticks":()}))
ax[0][0].imshow(img_gray, 'gray')
for i in range(1,(len(methods)+1)):
ret, dst = cv2.threshold(img_gray, 125, 200, methods[i-1])
ax[i//4][i%4].imshow(dst, 'gray')
ax[i//4][i%4].set_title('type={}'.format(method_name[i-1]))
適応的2値化処理
閾値を固定せず、画素ごとに変化させる二値化処理。画素毎に異なる閾値を算出する手法です。
fig, ax = plt.subplots(2, 2, figsize=(11, 6), subplot_kw=({"xticks":(), "yticks":()}))
methods = [cv2.ADAPTIVE_THRESH_MEAN_C, cv2.ADAPTIVE_THRESH_GAUSSIAN_C]
types = [cv2.THRESH_BINARY, cv2.THRESH_BINARY_INV]
m_names = ['ADAPTIVE_THRESH_MEAN_C', 'ADAPTIVE_THRESH_GAUSSIAN_C']
t_names = ['THRESH_BINARY', 'THRESH_BINARY_INV']
for i in range(2):
for j in range(2):
dst = cv2.adaptiveThreshold(img_gray, 200, methods[i], types[j], 11, 2)
ax[i][j].imshow(dst, 'gray')
ax[i][j].set_title('{} : {}'.format(m_names[i], t_names[j]))
2値画像の基本処理と計測
連結性
注目画素に対して上下左右の画素を4近傍と呼び、その近傍に対して注目画素の連結を定義したものを4連結と呼びます。
4近傍に斜め方向の近傍を加えたものを8近傍と呼び、連結を定義したものを8連結と呼びます。
連結している画素の集合は連結成分と呼びます。
対象の連結成分のなかにあり、背景に連結していない白画素の集合を、穴と呼びます。
例を以下に示します。
kinbo4 = np.array([[[255,255,255], [0,0,255], [255,255,255]],
[[0,0,255], [255,0,0], [0,0,255]],
[[255,255,255], [0,0,255], [255,255,255]]])
kinbo8 = np.array([[[0,0,255], [0,0,255], [0,0,255]],
[[0,0,255], [255,0,0], [0,0,255]],
[[0,0,255], [0,0,255], [0,0,255]]])
ana = np.array([[0,0,0,0,0,0,0],
[0,0,1,1,1,1,0],
[0,1,1,0,1,1,0],
[0,1,0,1,1,1,0],
[0,1,1,1,0,1,0],
[0,1,1,1,1,0,0],
[0,0,0,0,0,0,0]])
fig, ax = plt.subplots(1, 3, figsize=(10, 3), subplot_kw=({"xticks":(), "yticks":()}))
ax[0].imshow(kinbo4);
ax[1].imshow(kinbo8);
ax[2].imshow(-ana, 'gray');
1つ目と2つ目の画像は注目画像(赤)と近傍点(青)の画素を表しています。
3つ目の画像に対してそれぞれの近傍点を見ていくと、4連結(1つ目の画像)では1つの穴、8連結では3つの穴となります。
4連結では斜めのつながりを連結とみなしません。
輪郭追跡
連結成分の境界を求めることを、輪郭追跡と呼びます。
-
ラスタスキャンによって、白画素から黒画素に変わる画素を探索します。
ラスタスキャンとは、画像の左上を起点に、左端から右に画素を調べ、右端についたら行を1つ下がって左端から右に画素を調べる操作です。
探索した方向を進入方向とします。 - 進入方向を基点に番号順に右回りに黒画素を探索します。
- 見つかった黒画素に移動する。黒画素が開始点で、かつ、次の移動点が追跡済みの場合は、処理を終了し、追跡結果を登録します。
そうでない場合は、2の処理を繰り返します。
連結成分の外輪郭の場合、追跡は右回りになり、内輪郭すなわち穴の輪郭の場合は、追跡は左回りとなります。
輪郭追跡時に、開始点座標と方向コード（注目点右の画素位置を0とし時計回りに7まで）を並べたものをチェインコードと呼びます。
収縮・膨張処理
収縮
背景または穴に接する対象の画素をひとまわりはぎとる処理を、収縮と呼びます。
この処理を繰り返すと、対象の領域は小さくなります。
この処理で、1つの連結成分が複数に分割されることがある。すなわち、連結性を保持しない処理です。
cv2.dilateを使って収縮の様子を確認します。
kernels = [2, 3, 4]
iters = [1, 2, 3]
fig, ax = plt.subplots(3, 3, figsize=(16, 9), subplot_kw=({"xticks":(), "yticks":()}))
for i in range(3):
for j in range(3):
kernel = np.ones((kernels[i], kernels[i]), np.uint8)
dst = cv2.dilate(img_gray, kernel, iterations = iters[j])
ax[i][j].imshow(dst, 'gray')
ax[i][j].set_title('kernel={} : iter={}'.format(kernels[i], iters[j]))
膨張
背景または穴に接する対象の画素をひとまわり加える処理を、膨張と呼びます。
この処理を繰り返すと、穴は小さくなります。
この処理で、複数の連結成分が1つの連結成分となることがある。すなわち、膨張も連結性を保持しない処理です。
膨張の場合はcv2.erodeを使用します。
kernels = [2, 3, 4]
iters = [1, 2, 3]
fig, ax = plt.subplots(3, 3, figsize=(16, 10), subplot_kw=({"xticks":(), "yticks":()}))
for i in range(3):
for j in range(3):
kernel = np.ones((kernels[i], kernels[i]), np.uint8)
dst = cv2.erode(img_gray, kernel, iterations = iters[j])
ax[i][j].imshow(dst, 'gray')
ax[i][j].set_title('kernel={} : iter={}'.format(kernels[i], iters[j]))
クロージング・オープニング
同じ回数膨張して収縮する処理をクロージングと呼び、画像中の小さな穴を除くことができます。
同じ回数収縮して膨張する処理をオープニングと呼び、画像中の小さな連悦部分を除くことができます。
cv2.erodeとcv2.dilateを使って実装します。
n = [1, 2, 3]
fig, ax = plt.subplots(2, 3, figsize=(16, 6), subplot_kw=({"xticks":(), "yticks":()}))
for i in range(2):
for j in range(3):
kernel = np.ones((2, 2), np.uint8)
if i == 0:
dst = cv2.erode(img_gray, kernel, iterations = n[j])
dst = cv2.dilate(dst, kernel, iterations = n[j])
typ = 'closing'
else:
dst = cv2.dilate(dst, kernel, iterations = n[j])
dst = cv2.erode(img_gray, kernel, iterations = n[j])
typ = 'opening'
ax[i][j].imshow(dst, 'gray')
ax[i][j].set_title('type={} iter={}'.format(typ, n[j]))
また、cv2.morphologyExを使っても同様の処理ができます。
methods = [cv2.MORPH_CLOSE, cv2.MORPH_OPEN, cv2.MORPH_GRADIENT]
method_names = ['MORPH_CLOSE', 'MORPH_OPEN', 'MORPH_GRADIENT']
kernel = np.ones((3,3), np.uint8)
fig, ax = plt.subplots(1, 3, figsize=(16, 6), subplot_kw=({"xticks":(), "yticks":()}))
for i in range(3):
dst = cv2.morphologyEx(img_gray, methods[i], kernel)
ax[i].imshow(dst, 'gray')
ax[i].set_title('method={}'.format(method_names[i]))
ラベリング
同じ連結成分を構成する画素に同じ番号を付け、異なる連結部分を構成する画素に異なる番号を付ける処理をラベリングと呼びます。
ラスタスキャンによるラベリング
- ラスタスキャンを行い、ラベルの付いていない画素を調べ、見つかったら注目画素とします。
- 注目画素の上の画素(または左上の画素)がラベルを持つとき、上の画素(または左上の画素)の(左上の画素)ラベルを注目画素に付けます。
左の画素(及び上、右上の画素)がラベルを持ち、注目画素のラベルと異なるときルックアップテーブルにそれらのラベルが同一連結成分に属することを記録します。
(右上の画素がラベルを持ち、注目画素のラベルが異なり、かつ右の画素が白画素のとき、ルックアップテーブルに2つのラベルが同一連結成分に属することを記録します。) - 注目画素の上の画素(および左上の画素)が白画素で、左の画素がラベルを持つとき、そのラベルを注目画素に付けます。
- 注目画素の上も左も(左上の画素)も白画素のとき、新しいラベルを注目画素に付けます。
- 走査する画素がなくなるまで1から繰り返します。
つぎに、再度ラスタスキャンを行い、ルックアップテーブルを参照しながら、同一の連結成分に属するラベル群から、最も小さいラベルを付け直します。
輪郭追跡によるラベリング
輪郭追跡した外輪郭の結果で、1というラベルが付いているとき、
ラベル1を除いた連結成分に対して、左上の画素から再度輪郭追跡結果にラベル2を付けます。
さらに、ラベル2を除いた連結成分に対して、同様の処理を行い、ラベルを付けるべき画素がなくなるまで処理を繰り返します。
同じ連結成分に付いたラベルを統合して、同じラベル番号を振り直すことによって、ラベリングを行うことができます。
ここでは、クリオネの写真で例を示していきます。
cv2.thresholdで2値化を行い、cv2.connectedComponentsで連結した領域を抽出します。
import random
mono_src = cv2.threshold(img_gray2, 86, 255, cv2.THRESH_BINARY)[1].astype('uint8')
ret, markers = cv2.connectedComponents(mono_src)
color_src = img_rgb2.copy()
height, width = mono_src.shape[:2]
colors = []
# 色を領域の数だけ用意
for i in range(1, ret + 1):
colors.append(np.array([random.randint(0, 255), random.randint(0, 255), random.randint(0, 255)]))
# markersの値に応じて色を決めていく
for y in range(0, height):
for x in range(0, width):
if markers[y, x] > 0:
color_src[y, x] = colors[markers[y, x]]
else:
color_src[y, x] = [0, 0, 0]
fig, ax = plt.subplots(1, 2, figsize=(12, 6), subplot_kw=({"xticks":(), "yticks":()}))
ax[0].imshow(img_rgb2);
ax[1].imshow(color_src);
ラベリング処理の結果では、ノイズの部分にもラベルが割り振られています。
詳しい情報も得られるcv2.connectedComponentsWithStatsがあるので、
これを利用して面積でフィルタをかけ、ノイズを除きます。
mono_src = cv2.threshold(img_gray2, 86, 255, cv2.THRESH_BINARY)[1].astype('uint8')
nlabels, labeling, contours, CoGs = cv2.connectedComponentsWithStats(mono_src)
color_src = img_rgb2.copy()
height, width = mono_src.shape[:2]
colors = []
# 色を領域の数だけ用意
for i in range(1, nlabels + 1):
colors.append(np.array([random.randint(0, 255), random.randint(0, 255), random.randint(0, 255)]))
for y in range(0, height):
for x in range(0, width):
if labeling[y, x] > 0:
_,_,_,_,size = contours[labeling[y, x]]
# 面積が小さいものは対象外
if size >= 20:
color_src[y, x] = colors[labeling[y, x]]
else:
color_src[y, x] = [0, 0, 0]
else:
color_src[y, x] = [0, 0, 0]
fig, ax = plt.subplots(1, 2, figsize=(12, 6), subplot_kw=({"xticks":(), "yticks":()}))
ax[0].imshow(img_rgb2);
ax[1].imshow(color_src);
ディープラーニングじゃなくてもこんな画像がだせました。
形状特徴パラメータ
2値化画像中の連結成分は、形状の特徴によって分類できます。この特徴を数値化したものを形状特徴パラメータと呼ばれます。
クリオネ(1匹)の画像で見ていきます。
img_g = img_tmp_gray.copy()
img_c = img_tmp_rgb.copy()
# 2値化
ret, thresh = cv2.threshold(img_g,110,255,cv2.THRESH_BINARY)
# 膨張・収縮処理
kernel = np.ones((3, 3), np.uint8)
thresh = cv2.erode(thresh, kernel, iterations = 2)
thresh = cv2.dilate(thresh, kernel, iterations = 3)
thresh = cv2.erode(thresh, kernel, iterations = 1)
# 輪郭の抽出
contours, hierarchy = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
# 輪郭の特徴(モーメント)の算出
cnt = contours[0]
M = cv2.moments(cnt)
fig, ax = plt.subplots(1, 2, figsize=(10, 4), subplot_kw=({"xticks":(), "yticks":()}))
ax[0].imshow(img_tmp_rgb);
ax[1].imshow(thresh, 'gray');
重心
連結成分の画素に等しい重さがあると仮定したときの、図形全体の重さの中心
cx = int(M['m10']/M['m00'])
cy = int(M['m01']/M['m00'])
print(cx, cy)
54 55
面積
連結成分を構成する画素の数
area = cv2.contourArea(cnt)
area
1975.0
周囲長
輪郭追跡し一周する移動量
perimeter = cv2.arcLength(cnt,True)
perimeter
233.82
円形度
図形がどれだけ円に近いかを表す尺度。円に近いほど1に近づく。
$$4\pi S/L^2$$
roundness = 4*np.pi*area / perimeter**2
roundness
0.4539
オイラー数
連結成分の数から穴の数を引いたもの
主軸
画素の位置に重み付けをして合計した数値、モーメント特徴を利用して主軸方向を求めることができます。
$$
\theta=\frac{1}{2}\tan^{-1}(\frac{2M(1,1)}{M(2,0)-M(0,2)})
$$
theta=-0.5*np.arctan((2*M['m11'])/(M['m20']-M['m02'])) * 180/np.pi
theta
44.12
外輪郭・外接長方形・外接円・主軸
fig, ax = plt.subplots(2, 4, figsize=(14, 6), subplot_kw=({"xticks":(), "yticks":()}))
# 主軸
img_c1 = img_c.copy()
rows,cols = img_c1.shape[:2]
[vx, vy, x, y] = cv2.fitLine(cnt, cv2.DIST_L2, 0, 0.01, 0.01)
lefty = int((-x*vy/vx) + y)
righty = int(((cols-x)*vy/vx)+y)
dst = cv2.line(img_c1, (cols-1, righty), (0, lefty), (0, 255, 0), 1)
ax[0][0].imshow(dst);
# 外輪郭1
img_c1 = img_c.copy()
epsilon = 0.01*cv2.arcLength(cnt,True)
approx = cv2.approxPolyDP(cnt,epsilon,True)
dst = cv2.polylines(img_c1, [approx], True, (255, 255, 0))
ax[0][1].imshow(dst);
# 外輪郭2(凸包(Convex Hull))
img_c1 = img_c.copy()
hull = cv2.convexHull(cnt)
dst = cv2.polylines(img_c1, [hull], True, (255, 255, 0))
ax[0][2].imshow(dst);
# 外接長方形
img_c1 = img_c.copy()
x,y,w,h = cv2.boundingRect(cnt)
img_rect = cv2.rectangle(img_c1, (x, y), (x+w, y+h), (255, 0, 0), 1)
ax[0][3].imshow(img_rect);
# 回転を考慮した外接長方形
img_c1 = img_c.copy()
rect = cv2.minAreaRect(cnt)
box = cv2.boxPoints(rect)
box = np.int0(box)
img_rot_rect = cv2.drawContours(img_c1, [box], 0, (0, 255, 0), 1)
ax[1][0].imshow(img_rot_rect);
# 外接円
img_c1 = img_c.copy()
(x,y),radius = cv2.minEnclosingCircle(cnt)
center = (int(x),int(y))
radius = int(radius)
img_circle = cv2.circle(img_c1, center, radius, (255, 0, 0), 1)
ax[1][1].imshow(img_circle);
# 外接楕円
img_c1 = img_c.copy()
ellipse = cv2.fitEllipse(cnt)
img_ellipse = cv2.ellipse(img_c1, ellipse, (0, 255, 0), 1)
ax[1][2].imshow(img_ellipse);
# 重心
img_c1 = img_c.copy()
img_center = cv2.circle(img_c1, (cx,cy), 3, (255, 0, 0), -1)
ax[1][3].imshow(img_center);
距離と距離変換画像
距離の定義
2つの画素$A(x_{a},y_{a})$と$B(x_{b},y_{b})$の間の距離は、経路の制約によって以下のような定義があります。
ユークリッド距離：$\sqrt{(x_{a}-x_{b})^2+(y_{a}-y_{b})^2}$
市街地距離：$|x_{a}-x_{b}|+|y_{a}-y_{b}|$
チェス盤距離：$max(|x_{a}-x_{b}|,|y_{a}-y_{b}|)$
ユークリッド距離は、AとBを直線で結んだ距離、市街地距離は、4近傍距離とも呼び上下左右の4近傍の移動経路による距離です。
チェス盤距離は、8近傍距離とも呼び、8近傍の移動経路による距離です。
距離変換画像
連結成分の各画素の背景からの最短距離を与える変換を距離変換と呼び、変換後の画像を距離変換画像と呼びます。
2値化した画像にcv2.distanceTransformを使用します。
img_g1 = img_g.copy()
ret, thresh = cv2.threshold(img_g1,110,255,cv2.THRESH_BINARY)
kernel = np.ones((2, 2), np.uint8)
thresh = cv2.erode(thresh, kernel, iterations = 2)
thresh = cv2.dilate(thresh, kernel, iterations = 3)
dist_transform = cv2.distanceTransform(thresh, cv2.DIST_L2, 5)
plt.imshow(dist_transform)
plt.xticks([]);
plt.yticks([]);
線画像のベクトル化
ベクトル化の流れ
連結成分の連結性を保存したまま画素を削る処理を細線化と呼びます。
細線化された2値画像の画素は、特徴点である端点、分岐点、通過点の3種類に分類されます。
図に例を示す。赤い点が端点、緑色が分岐点、青の点が通過点を表しています。
blue=np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,255,0,0,0,0,0,0],
[0,0,0,255,0,255,0,255,0,255,0,0,0],
[0,0,255,0,255,0,255,0,255,0,255,0,0],
[0,0,0,0,255,0,255,0,255,0,0,0,0],
[0,0,0,0,0,0,255,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0]])
green=np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,255,0,255,0,255,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0]])
red =np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,255,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,255,0,0,0,0,0,0,0,0,0,255,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,255,0,0,0,255,0,0,0,0],
[0,0,0,0,0,0,255,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0]])
tmp_re = cv2.merge([red, green, blue])
plt.imshow(tmp_re);
plt.xticks([]);
plt.yticks([]);
細線化手法
対象を黒画素、背景を白画素とします。
注目画素を中心に3×3の画素値のパターンを観測し、注目画素が以下の3つの条件を満たすとき、その値を白画素にします。
- 境界上にある黒画素であること。
- 白画素に変更しても連結性が保持されること。
- 端点でないこと。
画素値を更新する方法には、逐次法と並列法があります。
逐次法は、更新した結果を用いて、すぐに次の処理を行います。
並列法は、消去可能かどうかを検証する画像と更新を行う画像を別々に持ち、更新用画像に処理結果を反映していきます。
高速道路の図を使って説明していきます。
細線化は2値化した画像に対してcv2.ximgproc.thinningで処理します。
ret, thresh = cv2.threshold(img_gray3, 120, 255, cv2.THRESH_BINARY_INV)
skeleton1 = cv2.ximgproc.thinning(thresh, thinningType=cv2.ximgproc.THINNING_ZHANGSUEN)
skeleton2 = cv2.ximgproc.thinning(thresh, thinningType=cv2.ximgproc.THINNING_GUOHALL)
fig, ax = plt.subplots(2, 2, figsize=(10, 7), subplot_kw=({"xticks":(), "yticks":()}))
ax[0][0].imshow(img_rgb3);
ax[0][1].imshow(thresh);
ax[1][0].imshow(skeleton1);
ax[1][1].imshow(skeleton2);
細線の特徴点の抽出
細線化処理の結果は、端点、分岐点、孤立点、通過点の特徴量に分類されます。
特徴点と特徴点によって区切られた画素列を、セグメントと呼びます。
テンプレートマッチングにより特徴点の抽出を行います。
各特徴点のテンプレートを用意します。
# 端点
tanten1 = np.array([[[0,0,0], [0,255,0], [0,255,0]],
[[0,0,0], [0,255,0], [0,255,255]],
[[0,0,0], [0,255,0], [255,255,0]],
[[0,0,0], [0,255,0], [255,255,255]]])
tanten2 = np.array([[0,0,0], [0,255,0], [255,0,0]])
tmp1 = np.array([[np.rot90(tanten1[j], i) for i in range(4) for j in range(4)]]).reshape([16,3,3]).astype('uint8')
tmp2 = np.array([[np.rot90(tanten2, i) for i in range(4)]]).reshape([4,3,3]).astype('uint8')
tanten = np.concatenate([tmp1,tmp2],axis=0)
# 3分岐
bunki3_1 = np.array([[[0,255,0], [255,255,255], [0,0,0]],
[[255,255,0], [255,255,255], [0,0,0]],
[[0,255,255], [255,255,255], [0,0,0]],
[[255,255,255], [255,255,255], [0,0,0]],
[[0,255,0], [0,255,0], [255,0,255]]])
#bunki3_1 = np.array([[[0,255,0], [255,255,255], [0,0,0]]])
bunki3_2 = np.array([[255,0,0], [0,255,0], [255,0,255]])
tmp1 = np.array([[np.rot90(bunki3_1[j], i) for i in range(4) for j in range(5)]]).reshape([20,3,3]).astype('uint8')
tmp2 = np.array([[np.rot90(bunki3_2, i) for i in range(4)]]).reshape([4,3,3]).astype('uint8')
bunki3 = np.concatenate([tmp1,tmp2],axis=0)
# 4分岐
bunki4_1 = np.array([[[255,255,0],[255,255,255],[0,255,0]],
[[255,255,255],[255,255,255],[0,255,0]],
[[255,255,255],[255,255,255],[255,255,0]]])
bunki4_2 = np.array([[[255,255,255],[255,255,255],[255,255,255]],
[[0,255,0],[255,255,255],[0,255,0]],
[[255,0,255],[0,255,0],[255,0,255]],
[[255,0,255],[0,255,0],[255,255,255]]])
bunki4_1 = np.array([[np.rot90(bunki4_1[j], i) for i in range(4) for j in range(3)]]).reshape([12,3,3]).astype('uint8')
bunki4=np.concatenate([bunki4_1,bunki4_2],axis=0).astype('uint8')
例えば端点であれば
のような画像です。
dst = img_rgb3.copy()
test = skeleton2.copy()
for i in range(20):
templ = tanten[i]
result = cv2.matchTemplate(test, templ, cv2.TM_CCOEFF_NORMED)
loc = np.where(result >= 0.9)
for pt in zip(*loc[::-1]):
cv2.rectangle(dst, (pt[0]-templ.shape[1],pt[1]-templ.shape[0]), (pt[0] + templ.shape[1], pt[1] + templ.shape[0]), (0,0,255), 2)
for i in range(20):
templ = bunki3[i]
result = cv2.matchTemplate(test, templ, cv2.TM_CCOEFF_NORMED)
loc = np.where(result >= 0.8)
for pt in zip(*loc[::-1]):
cv2.rectangle(dst, (pt[0]-templ.shape[1],pt[1]-templ.shape[0]), (pt[0] + templ.shape[1], pt[1] + templ.shape[0]), (0,255,0), 2)
for i in range(16):
templ = bunki4[i]
result = cv2.matchTemplate(test, templ, cv2.TM_SQDIFF_NORMED)
loc = np.where(result <= 0.5)
for pt in zip(*loc[::-1]):
cv2.rectangle(dst, (pt[0]-templ.shape[1],pt[1]-templ.shape[0]), (pt[0] + templ.shape[1], pt[1] + templ.shape[0]), (255,0,0), 1)
plt.imshow(dst)
plt.xticks([])
plt.yticks([])
3分岐と4分岐がうまく分けられていませんが、ある程度検出できました。
これで2値画像処理についての勉強は終わります。
次回
領域処理
参考
ディジタル画像処理[改訂第二版] | ディジタル画像処理編集委員会 |本 | 通販 | Amazon