Scilabを使った画像の特徴量抽出で工夫した話
目次
前提
図が黒，地が白のグレースケール画像の特徴量を抽出するプログラムを作成するときに工夫した話．
友達が考案し，自分が作成した．考案者から許諾済み．
※ 図は模様などのこと，地は背景のこと．
画像の前処理
図である黒の特徴を捉えるために，前処理として画像は画素反転後，他の画像とのスケール合わせに画素濃度の正規化処理を行う．
以降は，前処理後の画像を処理画像という．
今回は統計的な特徴量である重心について考えていく．
画像の特徴量(重心)
画像における重心とは，座標に画素値を重み付けし，その総和を取ること．つまり，画素値の平均や中心を意味する．
横方向の重心を$G_x, G_y$として数式にすると次のように表す．
ただし，$x$:横方向の座標，$y$:縦方向の座標値，$W$:横の長さ，$H$:縦の長さ，$\hat{I}$は処理画像,$\hat{I}_{x,y}$は処理画像の$(x,y)$座標の画素値を示す．
\begin{align}
G_x = \sum_{x=1}^W \sum_{y=1}^H x \hat{I}_{x,y} \\
G_y = \sum_{x=1}^W \sum_{y=1}^H y \hat{I}_{x,y}
\end{align}
画像の重心を求めるScilabプログラム
前処理のプログラムは省略する．
数式通りに画像の重心を求めると次のようなプログラムになる．
// 重心の計算
function [cgx, cgy] = calc_centerofgravity (im)
[h w] = size(im) // 画像サイズ取得
cgx = 0 // 横方向の重心
cgy = 0 // 縦方向の重心
// 重心計算
for y = 1 : h
for x = 1 : w
cgx = cgx + (x * im(y,x))
cgy = cgy + (y * im(y,x))
end
end
endfunction
しかし，これでは大量の画像を機械学習などに用いるとなると，ループのネストが深くなりすぎて，計算に時間がかかってしまう．
そこで，友達が考案したのが，画像のx座標の値を成分として持つ行列とy座標の値を成分として持つ行列を用いる方法である．
これらの行列$C_x, C_y$を次に示す．
\begin{align}
C_x &=
\begin{bmatrix}
1 & 2 & \cdots & W \\
1 & 2 & \cdots & W \\
\vdots & \vdots & & \vdots \\
1 & 2 & \cdots & W \\
\end{bmatrix} \\
C_y &=
\begin{bmatrix}
1 & 1 & \cdots & 1 \\
2 & 2 & \cdots & 2 \\
\vdots & \vdots & & \vdots \\
H & H & \cdots & H \\
\end{bmatrix}
\end{align}
この行列を用いて重心を求める方法は，
横方向についての重心について，処理画像$\hat{I}$と行列$C_x$の成分ごとの積を求めて，
得られる行列の総和を取る．
簡単に証明すると，行列$C_x$と処理画像$\hat{I}$の成分ごとの積，アダマール積を行う．
\begin{align}
C_x \otimes \hat{I} &=
\begin{bmatrix}
1 & 2 & \cdots & W \\
1 & 2 & \cdots & W \\
\vdots & \vdots & & \vdots \\
1 & 2 & \cdots & W \\
\end{bmatrix}
\otimes
\begin{bmatrix}
\hat{I}_{1,1} & \hat{I}_{1,2} & \cdots & \hat{I}_{1,W} \\
\hat{I}_{2,1} & \hat{I}_{2,2} & \cdots & \hat{I}_{2,W} \\
\vdots & \vdots & \ddots & \vdots \\
\hat{I}_{H,1} & \hat{I}_{H,2} & \cdots & \hat{I}_{H,W} \\
\end{bmatrix} \\
&=
\begin{bmatrix}
1 \times \hat{I}_{1,1} & 2 \times \hat{I}_{1,2} & \cdots & W \times \hat{I}_{1,W} \\
1 \times \hat{I}_{2,1} & 2 \times \hat{I}_{2,2} & \cdots & W \times \hat{I}_{2,W} \\
\vdots & \vdots & \ddots & \vdots \\
1 \times \hat{I}_{H,1} & 2 \times \hat{I}_{H,2} & \cdots & W \times \hat{I}_{H,W} \\
\end{bmatrix}
\end{align}
これで，得られる行列の成分は上述した横方向の重心の$x\hat{I}_{x,y}$の全成分である．
したがって，この行列の総和を求めると，$G_x$が求められる．
縦方向も同様に行う．
この方法を用いた重心を求めるScilabプログラムを次に示す．
xmatrixとymatrixが上述した行列$C_x, C_y$である．
sum()を用いて行列の成分の総和を求めることで，重心が求められる．
// 引数は画像の縦の長さと横の長さ
function [xmatrix, ymatrix] = make_xy_matrix (h,w)
// 1 〜 120 のステップ数 1 の行ベクトルを生成 (画像の各画素の x 座標の値となる)
xmatrix = linspace(1,w,w)
// xmatrix を画像の縦の長さ分コピーして行列にし、画像の x 座標の値を成分とする行列の生成完了
xmatrix = repmat(xmatrix, h, 1)
// 転置を行うことで、 y 座標の値を成分とする行列の生成完了
ymatrix = xmatrix'
endfunction
// 重心の計算
function [cgx, cgy] = calc_centerofgravity (im, xmatrix, ymatrix)
// 画像と xmatrix で成分ごとの積を求めてその総和を求めて x 軸の重心を求める
cgx = sum(im .* xmatrix)
// 画像と ymatrix で成分ごとの積を求めてその総和を求めて y 軸の重心を求める
cgy = sum(im .* ymatrix)
endfunction
計算時間の計測
100枚の図が黒，地が白の数字画像を読み込み，
for文を用いた方法と座標行列($C_x, C_y$)を用いた方法をそれぞれ実行し，
計算時間を計測する．
今回計測に用いるのは，timer()という関数．
これについては次のリンクを参考にした．
2つの方法を用いて重心を求め，それにかかる時間を計測する実行プログラムを次に示す．
前処理部と2つの方法の関数は省略する．
また，座標行列の生成にかかる時間は計測に含む．
clear // 全変数の値を消去
DATADIR = 'number/' // データディレクトリ
function measure_time ()
// for文による重心を計算する方法
timer() // 時間計測開始
for i = 0 : 9
for j = 0 : 9
filename = 'number'+string(i)+'_'+string(j)+'.pgm'
im = imread(DATADIR + filename) // 画像読み込み
rnim = rev_normalize_img(im) // 前処理
[gx gy] = calc_centerofgravity1(rnim) // 計算
end
end
t = timer() // 時間計測終了
disp("for文による重心を計算する方法: " + string(t) + "[s]")
// 座標行列による重心を計算する方法
timer() // 時間計測開始
[ xfilter yfilter ] = make_xy_filter(120,120) // 座標行列
for i = 0 : 9
for j = 0 : 9
filename = 'number'+string(i)+'_'+string(j)+'.pgm'
im = imread(DATADIR + filename) // 画像読み込み
rnim = rev_normalize_img(im) // 前処理
[gx gy] = calc_centerofgravity2(rnim, xfilter, yfilter) // 計算
end
end
t = timer() // 時間計測終了
disp("座標行列による重心を計算する方法: " + string(t) + "[s]")
endfunction
// 実行
measure_time()
計測結果は次に示す．
| 重心の計算方法 | 時間[秒] |
|---|---|
| for文による重心を計算する方法 | 5.494866 |
| 座標行列による重心を計算する方法 | 0.094056 |
だいぶ，計算時間が短縮された．
他の統計的な特徴量の分散や歪度などにも応用できる．