// 重心の計算
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

だいぶ，計算時間が短縮された．
他の統計的な特徴量の分散や歪度などにも応用できる．

Scilabによる画像の特徴量抽出で工夫した話

Scilabによる画像の特徴量抽出で工夫した話

目次

前提

画像の前処理

画像の特徴量(重心)

画像の重心を求めるScilabプログラム

計算時間の計測

Discussion