目次
はじめに
境界検出は画像処理の中で特に重要な手法です.画像の中に含まれるエッヂを強調する方法は多くの画像処理の研究者によって研究されました.このサイトで紹介した Total Variationに基づく画像再構成や,画像とフィルタの畳み込み演算なども,こういった目的でなされることが多いです.また,古典的な方法であれば,Hough変換などが考えられるでしょう,
スパースモデリングに基づく画像の再構成 Part2. Total Variation最小化(Split Bregman)に基づく画像再構成 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-| エルピクセル株式会社
問題設定
例えば,上のような画像を見てみましょう.画像の中心付近に左側と右側とを区別する線が見えると思いますが,このような境界線を検出するにはどうすればよいのでしょうか.
今回はJ.Brian Burnsらによって1986年に提案された,「Extracting Straight Lines」と,それの発展系の「LDS: a Line Segment Detector」を紹介します.
Extracting Straight Lines
Extracting Straight Lines のアルゴリズムでは,Gradient orientation と,Line Supportという概念が必要になります.それぞれに関して,具体的に説明していきます.
Gradient Orientation
Gradient とは勾配のことです.もし画像の輝度が大きく変化している部分があるなら,それは境界である可能性が高くなるでしょう.この輝度値変化のベクトルをGradient Orientation と論文では呼んでいます.
画像の輝度値を$I$として, \begin{align} \nabla u = \begin{pmatrix} \nabla_x u, \nabla_y u \end{pmatrix}^\top \end{align}
を各ピクセルで計算します. 先ほどの画像でこれを計算すると,以下のような結果が得られます.
Gradient Orientation
取り出す領域
Gradient Orientation
ただし,ノイズによってGradient は歪められてしまうので,$|\nabla u|$がある閾値よりも小さいものはカットします.
Line Support
次の段階として,このGradient の向きが同じ方向を向いている領域を探します. 境界線付近のGradient は,境界線に垂直な方向に揃っているはずです. 向きの揃っている領域,というのがLine Support です.
LSD: a Line Segment Detector [2] からの引用.
LSDのアルゴリズムでは,このLine Support 領域がある閾値よりも大きい場合,それを境界として認識します.
本記事では,Line Support 領域を求め,Line Support 領域を連結していくことにより,エッジを検出します.
実装
import numpy as np import matplotlib.pyplot as plt import cv2 import os import scipy as sp from scipy import signal import skimage import sys pi = np.pi
画像の作成
Is = np.zeros([50,80]) shp = np.shape(Is) for i in range(shp[0]): for j in range(int(shp[1]/2)): Is[i,j] = 125 def addGaussianNoise(src): [row,col] = np.shape(src) mean = 0 var = 5.5 sigma = 135 gauss = np.random.normal(mean,sigma,(row,col)) gauss = gauss.reshape(row,col) noisy = src + gauss return noisy Is = Is + addGaussianNoise(Is) ## Normalize Is = np.round(255*(Is - np.min(Is))/(np.max(Is) - np.min(Is)))
Gradient の計算
kernel_y = np.array([[-1, -1], [1, 1]]) kernel_x = np.array([[-1, 1], [-1, 1]]) Is_grady = cv2.filter2D(src=Is, kernel=kernel_y, ddepth=-1) Is_gradx = cv2.filter2D(src=Is, kernel=kernel_x, ddepth=-1) import matplotlib.patches as pat plt.figure() ax = plt.subplot(1,1,1) plt.quiver(Is_gradx,Is_grady,angles="xy",headwidth=3,color="#444444") rec1 = pat.Rectangle(xy = (35, 15), width = 10, height = 10, angle = 0, color = "blue", alpha = 0.5) ax.add_patch(rec1) ax.set_aspect('equal') plt.figure() ax2 = plt.subplot(1,1,1) plt.quiver(Is_gradx[15:25,35:45],Is_grady[15:25,35:45],angles="xy",headwidth=3,color="#444444") ax2.set_aspect('equal')
Gradient Orientation を16方向へ分類
Is_grad_norm = np.sqrt(Is_gradx**2 + Is_grady**2) ## Parameter1 TOL_GRAD = 100 Is_grad_theta = np.ones(shp)*pi/2 Is_grad_theta[Is_gradx != 0] = np.arctan(Is_grady[Is_gradx != 0]/Is_gradx[Is_gradx != 0]) Is_grad_theta = Is_grad_theta * (Is_grad_norm > Thre_grad) ## theta = pi/8 theta_all = np.arange(-pi,pi,theta) Is_grad_cls = -pi*np.ones(shp) Is_grad_theta_k = np.ones(shp) * sys.maxsize ones = np.ones(shp) for theta_i in theta_all: abs_grad_thetai = np.abs(Is_grad_theta - theta_i) grad_cls = abs_grad_thetai < Is_grad_theta_k Is_grad_cls[grad_cls] = theta_i*ones[grad_cls] Is_grad_theta_k[grad_cls] = abs_grad_thetai[grad_cls] tg = Is_grad_norm > Thre_grad Is_grad_cls = Is_grad_cls * tg
Line Support & 境界検出
TOL_ANG = pi/6 theta_ch = [-pi,-pi/2,0,pi/2] TOL_PIX = np.min(shp)*0.3 for theta_i in theta_ch: Is_grad_region = np.logical_and(Is_grad_cls > (theta_i - TOL_ANG),Is_grad_cls < (theta_i + TOL_ANG)) Is_grad_region = (Is_grad_region * tg + 0).astype("uint8") plt.figure() plt.imshow(Is,cmap="gray") plt.title(str(theta_i)) if theta_i in [-pi,0]: Is_grad_region_TOL = np.sum(Is_grad_region,axis=0) > TOL_PIX if any(Is_grad_region_TOL) == True: Is_grad_region_TOL_ind = np.where(Is_grad_region_TOL)[0] for Is_grad_region_TOL_i in Is_grad_region_TOL_ind: x1 = Is_grad_region_TOL_i x2 = Is_grad_region_TOL_i y1 = np.min(np.where([Is_grad_region[:,Is_grad_region_TOL_i]][0])) y2 = np.max(np.where([Is_grad_region[:,Is_grad_region_TOL_i]][0])) plt.plot([x1,x2],[y1,y2],'--r') if theta_i in [-pi/2,pi/2]: Is_grad_region_TOL = np.sum(Is_grad_region,axis=1) > TOL_PIX if any(Is_grad_region_TOL) == True: Is_grad_region_TOL_ind = np.where(Is_grad_region_TOL)[0] for Is_grad_region_TOL_i in Is_grad_region_TOL_ind: x1 = Is_grad_region_TOL_i x2 = Is_grad_region_TOL_i y1 = np.min(np.where([Is_grad_region[Is_grad_region_TOL_i,:]][0])) y2 = np.max(np.where([Is_grad_region[Is_grad_region_TOL_i,:]][0])) plt.plot([x1,x2],[y1,y2],'--k')
Output
上のコードは縦線,横線の境界線を検出するアルゴリズムですが,もう少し拡張すれば,斜めを向いている場合などでも出来そうですね.
参考文献
[1] Extracting Straight Lines
https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4767808
J. BRIAN BURNS, ALLEN R. HANSON, MEMBER, IEEE, AND EDWARD M. RISEMAN, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-8, NO. 4, JULY 1986 E
[2] LSD: a Line Segment Detector
http://www.ipol.im/pub/art/2012/gjmr-lsd/?utm_source=doi
Rafael Grompone von Gioi, J´er´emie Jakubowicz, Jean-Michel Morel, Gregory Randall Published in Image Processing On Line on 2012–03–24. Submitted