end0tknr's kipple - 新web写経開発

http://d.hatena.ne.jp/end0tknr/ から移転しました

正規化相互相関を用いた Harrisコーナー検出結果の対応点探索 (python)

前回と以前のエントリの続き。というか、これまで理解不足で、さかのぼっていた感じ

Harrisコーナー検出(再び)と、pythonによる実装 - end0tknr's kipple - 新web写経開発

自己相関と相互相関の違いや、正規化相互相関 - end0tknr's kipple - 新web写経開発

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
from scipy.ndimage import filters

im1 = numpy.array(Image.open('sf_view1.jpg').convert('L'))
im2 = numpy.array(Image.open('sf_view2.jpg').convert('L'))
wid = 5

def main():
    # 画像1のHariisコーナー検出
    harrisim = compute_harris_response(im1,5)
    filtered_coords1 = get_harris_points(harrisim,wid+1)
    d1 = get_descriptors(im1,filtered_coords1,wid)
    # 画像2のHariisコーナー検出
    harrisim = compute_harris_response(im2,5)
    filtered_coords2 = get_harris_points(harrisim,wid+1)
    d2 = get_descriptors(im2,filtered_coords2,wid)

    # 画像1,2のコーナーの対応点探索
    matches = match_twosided(d1,d2)

    plb.figure()
    plb.gray()
    plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches[:100])
    plb.savefig( '2_1_1.png' )

# グレースケール画像の画素に対応する Harrisコーナー検出器の応答関数
def compute_harris_response(im,sigma=3):
    # 微分係数
    imx = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # Harris行列の成分
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    Wdet = Wxx*Wyy - Wxy**2 # 判別式   det()
    Wtr = Wxx + Wyy         # 対角成分 trace()

    return Wdet / Wtr

# Harris応答画像からコーナーを返す。
# min_distはコーナーや画像境界から分離する最小ピクセル数
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    # 閾値(threshold)を超えるコーナー候補を探す
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 候補の座標を取得. 「T」は転置method
    coords = plb.array(harrisim_t.nonzero()).T
    # 候補の値を得る
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 候補をsort
    index = plb.argsort(candidate_values)
    # 許容する点の座標を配列に格納
    allowed_locations = plb.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 最小距離を考慮しながら、最良の点を得る
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                              (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0

    return filtered_coords

# 各点について、点の周辺で幅 2*wid+1 の近傍ピクセル値を返す。
#(点の最小距離 min_distance > wid を仮定する)
def get_descriptors(image,filtered_coords,wid=5):
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1,
                      coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)

    return desc


# match()の双方向で一致を調べる
def match_twosided(desc1,desc2,threshold=0.5):

    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)

    ndx_12 = numpy.where(matches_12 >= 0)[0]

    # 非対称の場合を除去する
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1

    return matches_12

# 【正規化相互相関】を用い、第1の画像の各コーナー点記述子について
#  第2の画像の対応点を選択
def match(desc1,desc2,threshold=0.5):
    n = len(desc1[0])

    # 対応点ごとの距離
    d = -numpy.ones((len(desc1),len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            d1 = (desc1[i] - numpy.mean(desc1[i])) / numpy.std(desc1[i])
            d2 = (desc2[j] - numpy.mean(desc2[j])) / numpy.std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:
                d[i,j] = ncc_value

    ndx = numpy.argsort(-d)
    matchscores = ndx[:,0]
    return matchscores


# 対応点を線で結んで画像を表示する
# 入力: im1,im2(配列形式の画像)、locs1,locs2(特徴点座標)
# machescores(match()の出力)、
# show_below(対応の下に画像を表示するならTrue)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):

    im3 = appendimages(im1,im2)
    if show_below:
        im3 = numpy.vstack((im3,im3))

    plb.imshow(im3)

    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0: plb.plot([locs1[i][1],locs2[m][1]+cols1],
                         [locs1[i][0],locs2[m][0]],'c')
    plb.axis('off')

# 2つの画像を左右に並べる
def appendimages(im1,im2):

    # 行の少ない方を選び、空行を0で埋める
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]

    if rows1 < rows2:
        im1 = numpy.concatenate((im1,numpy.zeros((rows2-rows1,im1.shape[1]))),
                                axis=0)
    elif rows1 > rows2:
        im2 = numpy.concatenate((im2,numpy.zeros((rows1-rows2,im2.shape[1]))),
                                axis=0)
    # 行が同じなら、0で埋める必要なし

    return numpy.concatenate((im1,im2), axis=1)

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます。数分待ちますが... f:id:end0tknr:20171003214833p:plain

Harrisコーナー検出(再び)と、pythonによる実装

先日までのエントリで「Harrisコーナー検出は理解できた」と思い 実践コンピュータビジョンのサンプルプログラムを写経しようとしたら、 理解度が不十分だった為、やり直し。

テイラー/マクローリン展開からのHarrisコーナー検出 (その1) - end0tknr's kipple - 新web写経開発

テイラー/マクローリン展開からのHarrisコーナー検出 (その2) - end0tknr's kipple - 新web写経開発

実践コンピュータビジョン サンプルプログラム

Harrisコーナー検出の基礎

まず、ある画像上の点xについて、行列 Miを定義します。

 \Large
M_I = \nabla I \ \nabla I^{T} =
\left[ \begin{array}{c} I_x \\ I_y \end{array} \right]
\left[ \begin{array}{cc} I_x & I_y \end{array} \right] =
\left[ \begin{array}{cc} I_x^2   & I_x I_y \\\
                          I_x I_y & I_y^2  \end{array} \right]

ここで、∇については以前のエントリに記載した通りです。

ベクトル偏微分演算子 : ∇ (ナブラ) - end0tknr's kipple - 新web写経開発

また、∇ I は微分係数Ix, Iyから構成され、次のように求められます。

 \Large
I_x = I \ast G_{σx} \ 、\ I_y = I \ast G_{σy}

( 上記において、Gσx, Gσy, * はそれぞれ、標準偏差σのガウス関数と、畳込み記号 )

また、先程のMiにおいては、階数(rank)=1 で、固有値
\lambda_1 = | \nabla I |^2  と \lambda_2 = 0

次にガウス分布 Gσ 等の重み行列を用いて  \Large \overline{M_{I}} を定義します ( Harris行列 ) 。

 \Large
\overline{M_{I}} = W \ast M_I

 \overline{M_{I}} は、 成分ごとの畳み込みで周辺ピクセルの局所平均となります。

( ガウス関数や畳み込みについては、以前のエントリ参照

画像処理 - エッジ抽出フィルタ(sobel, prewitt)と、平滑化フィルタ(gaussian) - end0tknr's kipple - 新web写経開発 )

 \nabla I の値により  \overline{M_{I}} 固有値は変化し、 次のケースに分けられます。

固有値 判定
λ1, λ2 が大きな正値 コーナーあり
λ1 が大きく, λ2 ≒ 0 エッジあり
λ1 ≒0, λ2 ≒ 0 平面

ただし、コーナー判定において固有値を算出すると、 計算コスト:大 の為、次のような判定式を使用します。

 \Large
  \frac { det( \overline{M_{I}} ) }{ trace( \overline{M_{I}} ) }

( det():行列式、trace():転置 )

pythonによるHarrisコーナー検出

実践コンピュータビジョン サンプルプログラム

相変わらずの上記の写経ですが...

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
from scipy.ndimage import filters

def main():
    im = numpy.array(Image.open('empire.jpg').convert('L'))
    # Harrisコーナー検出器の応答関数
    harrisim = compute_harris_response(im)
    #
    filtered_coords = get_harris_points(harrisim,6)
    # 検出されたコーナーの表示
    plot_harris_points(im, filtered_coords)

# グレースケール画像の画素に対応する Harrisコーナー検出器の応答関数
def compute_harris_response(im,sigma=3):
    # 微分係数
    imx = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # Harris行列の成分
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    Wdet = Wxx*Wyy - Wxy**2 # 判別式   det()
    Wtr = Wxx + Wyy         # 対角成分 trace()

    return Wdet / Wtr

# Harris応答画像からコーナーを返す。
# min_distはコーナーや画像境界から分離する最小ピクセル数
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    # 閾値(threshold)を超えるコーナー候補を探す
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 候補の座標を取得. 「T」は転置method
    coords = plb.array(harrisim_t.nonzero()).T
    # 候補の値を得る
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 候補をsort
    index = plb.argsort(candidate_values)
    # 許容する点の座標を配列に格納
    allowed_locations = plb.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 最小距離を考慮しながら、最良の点を得る
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                              (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0

    return filtered_coords

def plot_harris_points(image,filtered_coords):
    plb.figure()
    plb.gray()
    plb.imshow(image)
    plb.plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    plb.axis('off')
    plb.savefig( '2_1.png' )

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます f:id:end0tknr:20171003212441p:plain

自己相関と相互相関の違いや、正規化相互相関

2画像をHarrisコーナー検知し、対応点のマッチングしようとしたら、 正規化相互相関する必要があった為。

どれも、同じと言えば同じだし、違うと言えば違う

基本のピアソン相関係数

http://end0tknr.hateblo.jp/entry/20131206/1386283015 こちらはこれまでや、以前のエントリでも扱ってきました

 \Large
\frac{ \sum_{i=1}^n(x_i- \overline{x})(y_i-\overline{y}) }
     { \sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}
       \sqrt{\sum_{i=1}^n(y_i-\overline{y})^2}}

自己相関と相互相関

この辺りは、信号処理に関わってきます

相関 説明
自己相関 時間tだけ過去の自身の波形は現在の自身にどれだけ似ているか?
相互相関 別の波形とどれだけ似ているか? 片方が正弦波ならフーリエ変換

自己相関

 \Large{
\frac{ \sum_{i=1}^{n-k}(x_i- \overline{m1_{k}})(x_{i+k}-\overline{m2_{k}}) }
     { \sqrt{\sum_{i=1}^{n-k}(x_i    -\overline{m1_{k}})^2}
       \sqrt{\sum_{i=1}^{n-k}(x_{i+k}-\overline{m2_{k}})^2}}   }
\normalsize ※
m1_k = \frac{ \sum_{i=1}^{n-k} x_{i}   }{n-k}  
m2_k = \frac{ \sum_{i=1}^{n-k} x_{i+k} }{n-k}

相互相関

 \Large{
\frac{ \sum_{i=1}^{n-k}(x_i- \overline{m1_{k}})(y_{i+k}-\overline{m2_{k}}) }
     { \sqrt{\sum_{i=1}^{n-k}(x_i    -\overline{m1_{k}})^2}
       \sqrt{\sum_{i=1}^{n-k}(y_{i+k}-\overline{m2_{k}})^2}}  }
\normalsize ※
m1_k = \frac{ \sum_{i=1}^{n-k} x_{i}   }{n-k}  
m2_k = \frac{ \sum_{i=1}^{n-k} y_{i+k} }{n-k}

参考URL http://www.ipc.tohoku-gakuin.ac.jp/nken/java/autocorr/theory.PDF

正規化相互相関

例えば、2画像(I,T)をHarrisコーナー検知した後の対応点のマッチングに使用します

 \Large
\frac
{ {\sum_{j=1}^{N-1}} {\sum_{i=1}^{M-1}}
   (I(i,j) - \overline{I}) (T(i,j) - \overline{T})
}
{
  \sqrt{
      \sum_{j=1}^{N-1} \sum_{i=1}^{M-1}
         (I(i,j) - \overline{I})^2 ×
      \sum_{j=1}^{N-1} \sum_{i=1}^{M-1}
         (T(i,j) - \overline{T})^2
  }
}

ただし


\overline{I} = \frac{1}{MN} \sum_{j=1}^{N-1} \sum_{i=1}^{M-1} I(i,j) 
\overline{T} = \frac{1}{MN} \sum_{j=1}^{N-1} \sum_{i=1}^{M-1} T(i,j)

参考URL http://www.sic.shibaura-it.ac.jp/~yaoki/mediaeng/imgret4.pdf

ベクトル偏微分演算子 : ∇ (ナブラ)

「何? この逆三角形な記号」と思った。過去、使用していたことすら忘れてた。

参考 URL

http://www.iwata-system-support.com/CAE_HomePage/vector/vectana1/vectana1.html

http://www.iwata-system-support.com/CAE_HomePage/vector/vectana11/vectana11.html

∇ (ナブラ) の前に3次元面の勾配を考える

3次元の面:f(x,y,z)があるとき、ある座標での勾配は、


\Large
\mathrm{grad} f =
\frac{\partial f}{\partial x} \vec{i} + \frac{\partial f}{\partial y} \vec{j} +
\frac{\partial f}{\partial z} \vec{k}

のように各成分を偏微分することで、得られます。 ※  \vec{i} , \vec{j} ,  \vec{k} は単位ベクトル

∇ (ナブラ) とは、ベクトルの偏微分演算子

先程のgrad式を次のように定義したものが∇です


\Large
\nabla =
\frac{\partial}{\partial x} \vec{i} + \frac{\partial}{\partial y} \vec{j} +
\frac{\partial}{\partial z} \vec{k} =
\left( \begin{array}{c}
    \frac{\partial}{\partial x} \\\
    \frac{\partial}{\partial y} \\\
    \frac{\partial}{\partial z}
 \end{array} \right)

これで、SIFT特徴量検出のDoGやLoGの理解が進むかな?

画像処理 - vlfeatによるSIFT特徴点検出

先日のHarrisコーナー検出をその後、試してみましたが、どうも検出精度がイマイチ。

どうやら Harrisコーナー検出は拡大縮小のスケール変化に弱いらしい。

テイラー/マクローリン展開からのHarrisコーナー検出 (その2) - end0tknr's kipple - 新web写経開発

そこで、今回はSIFT特徴点検出を写経(Scale-Invariant Feature Transform).

実践コンピュータビジョン サンプルプログラム

SIFT特徴点の定義等については、別の機会で

SIFT特徴点検出とは?

次のurlに丁寧が解説があります。(きちんとは理解できていません)

http://www.vision.cs.chubu.ac.jp/sift/

前準備 - install vlfeat

http://www.vlfeat.org/

SIFT特徴量検出自体は、外部ライブラリに任せます

cd ~/local
wget http://www.vlfeat.org/download/vlfeat-0.9.20.tar.gz
tar -xvf vlfeat-0.9.20.tar.gz
cd vlfeat-0.9.20
make
cd ..
ln -s vlfeat-0.9.20/bin/glnxa64 vlfeat

vlfeat を使った python script

#!/usr/local/bin/python
# -*- coding: utf-8 -*-

from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
import os

IMNAME = 'empire.jpg'
SIFT_CMD = '/home/endo/local/vlfeat/sift'
PROCESS_IMAGE_PARAMS = "--edge-thresh 10 --peak-thresh 5"

def main():
    im1 = numpy.array(Image.open(IMNAME).convert('L'))
    # SIFT特徴点検出し、保存
    process_image(IMNAME,'empire.sift')
    # SIFT特徴点を読み出し
    l1,d1 = read_features_from_file('empire.sift')

    # SIFT特徴点を描画
    plb.figure()
    plb.gray()
    plot_features(im1,l1,circle=True)
    plb.savefig( '2_3_3.png' )


# 画像の特徴点を検出しfileに保存
def process_image(imagename,resultname,params=PROCESS_IMAGE_PARAMS):
    if imagename[-3:] != 'pgm':
        # pgmファイルを作成する
        im = Image.open(imagename).convert('L')
        imagename = 'tmp.pgm'
        im.save(imagename)

    cmmd = str(SIFT_CMD +" "+imagename+" --output="+resultname+ " "+params)
    os.system(cmmd)
    print 'processed to', resultname

# SIFT特徴量を読み込んで行列形式で返す
def read_features_from_file(filename):
    f = numpy.loadtxt(filename)
    return f[:,:4],f[:,4:] # 特徴点の配置と記述子


# 画像を特徴量とともに描画.
#   入力:im(配列形式の画像)、locs(各特徴量の座標とスケール、方向)
def plot_features(im,locs,circle=False):
    plb.imshow(im)
    if circle:
        for p in locs:
            draw_circle(p[:2],p[2])
    else:
        plb.plot(locs[:,0],locs[:,1],'ob')
    plb.axis('off')

def draw_circle(c,r):
    t = numpy.arange(0,1.01,.01)*2* numpy.pi
    x = r * numpy.cos(t) + c[0]
    y = r * numpy.sin(t) + c[1]
    plb.plot(x,y,'b',linewidth=2)

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます f:id:end0tknr:20171001072748p:plain

画像処理 - エッジ抽出フィルタ(sobel, prewitt)と、平滑化フィルタ(gaussian)

先日のHarrisコーナー検出から、SIFT特徴量へ移りたかったのですが、 画像の微分ガウス関数を理解できていないようなので、寄り道。

参考URL

前処理フィルタについて | 画像処理.com | キーエンス

【画像処理】ガウシアンフィルタの原理・特徴・計算式

フィルタの種類

名称 目的 備考
ソーベル(sobel) エッジ抽出 ※1
プレヴィット(prewitt) エッジ抽出
ガウシアン(gaussian) 平滑化

※1 中央画素に2がかかる為、 コントラストの少ないエッジを強調できるがノイズも拾う

エッジ抽出における微分(離散近似)

グレースケール画像 I の明度変化は、x, y方向の微分係数 Ix,Iy で表し、 画像勾配は次のようになります。

 \Large
\nabla I = \left( \begin{array}{c} I_x \\\ I_y  \end{array} \right)
= ( I_x I_y )^T

 \nabla はナブラ、Tは転置 です

画像の微分には、次のような離散近似を行うようです。

 \Large
I_x = I * D_x  I_y = I * D_y  

上記「Dx」「Dy」の代表的なフィルタが、sobel と prewitt で、 「*」は畳込みを表します。

ソーベルとプレヴィットのフィルタ

ソーベルフィルタ(Dx, Dy)


D_x =
\left[ \begin{array}{ccc}
Dx_{11} &  Dx_{21} &  Dx_{31} \\\
Dx_{12} &  Dx_{22} &  Dx_{32} \\\
Dx_{13} &  Dx_{23} &  Dx_{33} 
\end{array} \right]
=
\left[ \begin{array}{ccc}
-1 & 0 & 1 \\\
-2 & 0 & 2 \\\
-1 & 0 & 1
\end{array} \right]

D_y =
\left[ \begin{array}{ccc}
Dy_{11} &  Dy_{21} &  Dy_{31} \\\
Dy_{12} &  Dy_{22} &  Dy_{32} \\\
Dy_{13} &  Dy_{23} &  Dy_{33} 
\end{array} \right]
=
\left[ \begin{array}{ccc}
-1 & -2 & -1 \\\
 0 &  0 &  0 \\\
 1 &  2 &  1
\end{array} \right]

プレヴィット(Dx, Dy)


D_x =
\left[ \begin{array}{ccc}
Dx_{11} &  Dx_{21} &  Dx_{31} \\\
Dx_{12} &  Dx_{22} &  Dx_{32} \\\
Dx_{13} &  Dx_{23} &  Dx_{33} 
\end{array} \right]
=
\left[ \begin{array}{ccc}
-1 & 0 & 1 \\\
-1 & 0 & 1 \\\
-1 & 0 & 1
\end{array} \right]

D_y =
\left[ \begin{array}{ccc}
Dy_{11} &  Dy_{21} &  Dy_{31} \\\
Dy_{12} &  Dy_{22} &  Dy_{32} \\\
Dy_{13} &  Dy_{23} &  Dy_{33} 
\end{array} \right]
=
\left[ \begin{array}{ccc}
-1 & -1 & -1 \\\
 0 &  0 &  0 \\\
 1 &  1 &  1
\end{array} \right]

ソーベルによる畳込み演算の例

「畳込み」とは、 9画素それぞれに上記Dx, Dyなの係数をかけ、合算 した濃度値に置き換えることです。


I =
\left[ \begin{array}{cccc}
I_{11} &  I_{21} &  I_{31} & I_{41} \\\
I_{12} &  I_{22} &  I_{32} & I_{42} \\\
I_{13} &  I_{23} &  I_{33} & I_{43} \\\
I_{14} &  I_{24} &  I_{34} & I_{44}
\end{array} \right]
=
\left[ \begin{array}{cccc}
0 &  0 &  0 & 0 \\\
0 & 10 & 10 & 0 \\\
0 &  0 & 10 & 0 \\\
0 &  0 &  0 & 0
\end{array} \right]

例えば、上記画像Iのある画素に対するIxをソーベルフィルタで算出します。


Ix_{22} \\\
= (I_{11} Dx_{11}) + (I_{21} Dx_{21}) + (I_{31} Dx_{31})\\\
+ (I_{12} Dx_{12}) + (I_{22} Dx_{22}) + (I_{32} Dx_{32})\\\
+ (I_{13} Dx_{13}) + (I_{23} Dx_{23}) + (I_{33} Dx_{33}) \\\
= (0 \cdot -1) + (0 \cdot  0) + ( 0 \cdot 1) \\\
+ (0 \cdot -2) + (10 \cdot 0) + (10 \cdot 2) \\\
+ (0 \cdot -1) + (0 \cdot  0) + (10 \cdot 1) \\\
= 30

平滑化を行うガウシアン(gaussian)フィルタ

ガウシアンフィルタでは、以下のガウス分布を使用しますが、 畳込み?を使用する面では、先程のソーベル等と同様です。

 \Large
g(x,y,\sigma) =
\frac{1}{ \sqrt{2 \pi} \sigma }
exp ( \frac{x^2+y^2}{2 \sigma^2})

例えば、 g(x,y,\sigma) 標準偏差σ=1.3の場合、


K =
\Large{ \frac{1}{16} }
\left[ \begin{array}{ccc}
 1 & 2 & 1 \\\
 2 & 4 & 2 \\\
 1 & 2 & 1
\end{array} \right]

となりますし、更に大きな標準偏差の場合


K =
\Large{ \frac{1}{256} }
\left[ \begin{array}{ccccc}
1 &  4 &  6 &  4 & 1 \\\
4 & 16 & 24 & 16 & 4 \\\
6 & 24 & 36 & 24 & 6 \\\
4 & 16 & 24 & 16 & 4 \\\
1 &  4 &  6 &  4 & 1
\end{array} \right]

となり、平滑化の効果が大きくなります。

DHTMLX Libraries (DHX) - JavaScript/HTML5 UI Components Library

jsによるgrid表示やガントチャート作成に有効かもしれないので、メモ

dhtmlx.com

テイラー/マクローリン展開からのHarrisコーナー検出 (その2)

で、Harrisコーナー検出の理論編。

以下に書いてはみたものの、十分には理解できていない印象です。

参考url

コーナー検出法 - Wikipedia

画像処理の数式を見て石になった時のための、金の針 - Qiita

Harrisコーナー検出の定義式

Harrisコーナー検出とは、 「一次微分値(差分)が一方向に大きければエッジ、多方向に大きければコーナー」 の考えに基きます。

画像上の点の強度 I(x,y)と、そこから(u,v)移動した点強度 I(x+u,y+v) の 変化量 E(u,v)を次のように表します。

 \Large
E(u,v) = \sum_{x,y} w(x,y)(I(x+u,y+v) - I(x,y))^2 (式1)

※ w(x,y)は窓関数(window function)で、ある区間外では0となるもの(例:ガウス窓)。

Iの導関数をIx , Iyとして、先程の式1にある I(x+u,y+v)部分にテイラー展開


f(x) ≒ f(x_0) + f'(x_0)(x-x_0)

を適用すると、次のように変形/近似できます

 \Large
I(x+u,y+v) ≒ I(x,y) + I_u(x,y)u + I_v(x,y)v (式2)

よって、式1は次のようになります


E(u,v) \\
\Large ≒ \sum_{x,y} w(x,y)( I(x,y) + I_u(x,y)u + I_v(x,y)v - I(x,y))^2 \\
\Large = \sum_{x,y} w(x,y)( I_u(x,y)u + I_v(x,y)v )^2 (式3)

更に、次のようにMを設定して、行列式に変形します

 \Large
M = \sum_{x,y} w(x,y)
\left[ \begin{array}{cc}
      {I_X}^{2}  & {I_X}{I_Y} \\\
      {I_X}{I_Y} & {I_Y}^{2}
\end{array} \right]
 (式4)
 \Large
E(u,v) ≒ [u,v] M 
\left[ \begin{array}{c}
      u \\\ v
\end{array} \right]
 (式5)

式4に特異点分解を行うと、エッジやコーナーの判定ができますが、 固有値算出は計算量が多い為、次のRを用いて判定を行います。

 \Large
R = detM - k(trace M)^2

ここで、 ・detM とは、行列式です。(例 ad-bc)

線形代数 : 行列式とサラスの公式、そしてクラメル式によるn元1次連立方程式の解 - end0tknr's kipple - 新web写経開発

・trace Mとは、対角和です。(例 a+d)

・kは定数 (0.04~0.06)

・R=大→コーナー、R=小さい→フラット、R<0→エッジ

テイラー/マクローリン展開からのHarrisコーナー検出 (その1)

画像処理のコーナー検出手法であるHarrisコーナー検出を行おうとしていますが、 内部でテイラー展開を使用しているので、そこから振り返り。

今回の「その1」は、「テイラー/マクローリン展開の振り返り」まで

参考url

EMANの物理学・物理数学・テイラー展開

テイラー展開 - Wikipedia

テイラー展開とは

以下の通り


f(x) = f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2!} f''(x_0)(x-x_0)^2 +
       \frac{1}{3!} f'''(x_0)(x-x_0)^3 + \cdots  (式1)

式1に  x = x_0+h を代入して、次のようにも書けます


f(x_0+h) = f(x_0) + f'(x_0)h + \frac{1}{2!} f''(x_0)h^2 +
       \frac{1}{3!} f'''(x_0)h^3 + \cdots  (式2)

また、 x_0 が十分小さい時、式1は次のように書けます


f(x) ≒ f(x_0) + f'(x_0)(x-x_0)  (式3)

マクローリン展開とは

式1において  x_0 = 0 としたもので、次のように簡略化できます。


f(x) = f(0) + f'(0)x + \frac{1}{2!} f''(0)x^2 +
       \frac{1}{3!} f'''(0)x^3 + \cdots  (式4)

テイラー展開の適用例

単純?な適用


e^x = 1 + x + \frac{1}{2!} x^2 + \frac{1}{3!}x^3 + \cdots

\sin x = x - \frac{1}{3!} x^3 + \frac{1}{5!}x^5 + \cdots

2次式の変形への適用


f(x) =  - 3 - 4x + x^2

または上記にx=1周りのテイラー展開を適用すると、以下


f(x) ≒ (1-4-3) +(2-4)(x-1) + \frac{2}{2!}(x-1)^2
     = (-3) - 2(x-1) + (x-1)^2

pythonでのk近傍法( k-nearest neighbor algorithm, k-NN )

O'Reilly Japan - 実践 コンピュータビジョン

実践コンピュータビジョン サンプルプログラム

相変わらずの上記urlの写経。

k近傍法とは?

次のurlが視覚的にもとっても分かりやすいです

qiita.com

教師データとテストデータの作成

…といっても、中身は同じである points_normal.pkl と points_normal_test.pkl が 以下のscriptで作成されます。

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import numpy
import pickle

n = 200

def main():
    #2組の点の集合が距離を置いて配置
    make_dots_pair_1()
    #1組の点の集合のまわりに、もう1組の点の集合を
    #土星の輪のように配置
    make_dots_pair_2()

def make_dots_pair_1():
    # 標準正規分布による nx2 の行列
    class_1 = 0.6 * numpy.random.randn(n,2)
    class_2 = 1.2 * numpy.random.randn(n,2) + numpy.array([5,1])
    # ones()で作成した要素が全てnの配列を hstack()で連結
    labels = numpy.hstack((numpy.ones(n),-numpy.ones(n)))

    # pickleで保存
    with open('points_normal.pkl', 'w') as f:
        pickle.dump(class_1,f)
        pickle.dump(class_2,f)
        pickle.dump(labels,f)
    with open('points_normal_test.pkl', 'w') as f:
        pickle.dump(class_1,f)
        pickle.dump(class_2,f)
        pickle.dump(labels,f)

def make_dots_pair_2():
    # 正規分布と、その周りに輪を描くように点を配置
    class_1 = 0.6 * numpy.random.randn(n,2)
    r = 0.8 * numpy.random.randn(n,1) + 5
    angle = 2*numpy.pi * numpy.random.randn(n,1)
    class_2 = numpy.hstack((r*numpy.cos(angle),r*numpy.sin(angle)))
    labels = numpy.hstack((numpy.ones(n),-numpy.ones(n)))
    # pickleで保存
    with open('points_ring.pkl', 'w') as f:
        pickle.dump(class_1,f)
        pickle.dump(class_2,f)
        pickle.dump(labels,f)

if __name__ == '__main__':
    main()

k-NN実行

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy
import pickle
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb

KNN_MODEL = None

def load_model_and_test_data():
    # pickleからmodel dataをload
    with open('points_normal.pkl', 'r') as f:
        class_1 = pickle.load(f)
        class_2 = pickle.load(f)
        labels = pickle.load(f)
        # k近傍法でクラスタリングするmodel作成
        global KNN_MODEL
        KNN_MODEL = KnnClassifier(labels, numpy.vstack((class_1,class_2)))

    # pickleからtest dataをload
    with open('points_normal_test.pkl', 'r') as f:
        class_1 = pickle.load(f)
        class_2 = pickle.load(f)
        labels = pickle.load(f)

    return class_1, class_2, labels


# 描画用関数の定義
def my_classify(x,y):
    return numpy.array([KNN_MODEL.classify([xx,yy]) for (xx,yy) in zip(x,y)])


def plot_2D_boundary(plot_range,points,decisionfcn,labels,values=[0]):
    """ plot_range:(xmin,xmax,ymin,ymax)、points:クラスの点のリスト、
    decisionfcn:評価関数、
    labels:各クラスについてdecisionfcnが返すラベルのリスト、
    values:表示する判別の輪郭のリスト """

    clist = ['b','r','g','k','m','y'] # 各classの描画色

    # 等差数列( arange() )からmeshgrid作成
    x = numpy.arange(plot_range[0],plot_range[1],.1)
    y = numpy.arange(plot_range[2],plot_range[3],.1)
    xx,yy = numpy.meshgrid(x,y)
    xxx,yyy = xx.flatten(),yy.flatten()
    zz = numpy.array(decisionfcn(xxx,yyy))
    zz = zz.reshape(xx.shape)
    
    plb.contour(xx,yy,zz)   #等高線描画

    # クラスごとに正しい点には*、間違った点には'o'を描画する
    for i in range(len(points)):
        d = decisionfcn(points[i][:,0],points[i][:,1])
        correct_ndx = labels[i]==d
        incorrect_ndx = labels[i]!=d
        plb.plot(points[i][correct_ndx,0],points[i][correct_ndx,1],
                 '*',color=clist[i])
        plb.plot(points[i][incorrect_ndx,0],points[i][incorrect_ndx,1],
                 'o',color=clist[i])

    plb.axis('equal')


class KnnClassifier(object):
  def __init__(self,labels,samples):
    """ 教師データを使って分類器を初期化する """
    self.labels = labels
    self.samples = samples

  def classify(self,point,k=3):
    """ pointを教師データ中のk個の最近傍を使って分類し、
        ラベルを返す """

    # 全教師データへの距離を計算する
    dist = numpy.array([L2dist(point,s) for s in self.samples])
    # ソートする
    ndx = dist.argsort()
    # k個の最近傍を保持するのに辞書を用いる
    votes = {}
    for i in range(k):
      label = self.labels[ndx[i]]
      votes.setdefault(label,0)
      votes[label] += 1

    return max(votes, key=lambda x: votes.get(x))

def L2dist(p1,p2):
  return numpy.sqrt( sum( (p1-p2)**2) )


def main():
    class_1, class_2, labels = load_model_and_test_data()
    # 分類の境界を描画
    plot_2D_boundary([-6,6,-6,6], [class_1,class_2], my_classify,[1,-1])
    plb.savefig( '8_1_1.png' )


if __name__ == '__main__':
    main()

↑こう書くと↓こう出力されます f:id:end0tknr:20170923112633p:plain

データの永続化(serialize)において pickle for python ≒ Storable for perl

そんな気がします

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import pickle


def main():
    test_vals = {"hoge":"hoge_val","foo":2}
    print test_vals

    # バイナリ書込みmodeでopen
    # 他modeには以下があり.
    # r:読み, w:書き, a:追記, +:読書き両方   t:text(default),b:binary
    # またwith構文により、close()不要にしています
    with open('sample.pickle', mode='wb') as f:
        pickle.dump(test_vals, f, protocol=2) #fileへの書込み(保存)

    restored_vals = {}
    with open('sample.pickle', mode='rb') as f:
        restored_vals = pickle.load(f)        #fileからの読込み

    print restored_vals


if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます

$ python test_pickle.py 
{'foo': 2, 'hoge': 'hoge_val'}
{'foo': 2, 'hoge': 'hoge_val'}

pythonでのfile open/closeやwithの勉強にもなりました

PIL / Pillow , numpy による画像処理入門 その3 - 平板化(1次元化)と主成分分析(principal component analysis; PCA)

www.oreilly.co.jp

更に↑こちらで提供されているサンプルコードのまんま

実行はできましたが、算出結果の正しさは確認していません

#!/usr/bin/python
# -*- coding: utf-8 -*-

from PIL import Image
import numpy
import os


def main():
    # 指定dir以下のjpg file名一覧を取得
    imlist = get_imlist('a_thumbs')

    im = numpy.array(Image.open(imlist[0]))
    m,n = im.shape[0:2] # jpg size

    imnbr = len(imlist) # files size


    # 平板化(1次元化)し行列へ格納
    immatrix = numpy.array([numpy.array(Image.open(im)).flatten()
                            for im in imlist],'f')

    # 主成分分析 - 返り値:写像行列(次元の重要度順), 分散, 平均
    V,S,immean = my_pca(immatrix)
    print V


def get_imlist(path):
  """ pathに指定されたフォルダのすべてのjpgファイル名のリストを返す """
  return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]


def my_pca(X):
  """ 主成分分析
   入力:X, 訓練データを平板化した配列を行として格納した行列
   出力:写像行列(次元の重要度順), 分散, 平均 """

  # 次元数を取得
  num_data,dim = X.shape

  # データをセンタリング
  mean_X = X.mean(axis=0)
  X = X - mean_X

  if dim>num_data: 
    # PCA - 高次元のときはコンパクトな裏技を用いる
    M = dot(X,X.T) # 共分散行列
    e,EV = numpy.linalg.eigh(M) # 固有値と固有ベクトル
    tmp = dot(X.T,EV).T # ここがコンパクトな裏技
    V = tmp[::-1] # 末尾の固有ベクトルほど重要なので、反転する
    S = sqrt(e)[::-1] # 固有値の並びも反転する
    for i in range(V.shape[1]):
      V[:,i] /= S
  else:
    # PCA - 低次元なら特異値分解を用いる
    U,S,V = numpy.linalg.svd(X) 
    V = V[:num_data] # 最初のnum_dataの分だけが有用

  # 写像行列と、分散、平均を返す
  return V,S,mean_X


if __name__ == '__main__':
    main()

pylab for python による画像処理入門 その2

www.oreilly.co.jp

またも↑こちらで提供されているサンプルコードのまんま

画像に対し、点や線、文字列の描画

#!/usr/bin/python
# -*- coding: utf-8 -*-

# from PIL import Image

# import numpy
# #from pylab import *
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb


# # 配列に画像読込み
# im = numpy.array(Image.open('empire.jpg'))

# # 画像表示
# # が、私のcent7環境には x11がない為、comment out
# # imshow(im)

# 点の座標
x = [100,100,400,400]
y = [200,500,200,500]

plb.plot(x,y,'r*')    # 赤い星マークを4点に描画
plb.plot(x[:2],y[:2]) #最初の2点間に線を描画

plb.title('test pylab')  # タイトル追加
plb.axis('off')
# show()

# savefig()は、画像フォーマットを拡張子で判断しますが
# 引数で指定することも可能
# https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.savefig.html

plb.savefig( '1_2_1.png' )

「import pylab」の代わりに「import matplotlib , matplotlib.pylab」

import pylab as plb

ヘッダ部を↑このように記載すると、 画面のない私のcentos on virtual boxでは実行時に

self.tk = _tkinter.create(screenName, baseName, className, interactive, wantobjects, useTk, sync, use)
_tkinter.TclError: no display name and no $DISPLAY environment variable

のようなエラーとなる為、次のようにimportしています。

import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb

等高線図とヒストグラム

#!/usr/bin/python
# -*- coding: utf-8 -*-

from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb

# 画像を配列に読み込む
im = numpy.array(Image.open('empire.jpg').convert('L'))

plb.figure()                    # 新規図
plb.gray()                      # 色を使わない
plb.contour(im, origin='image') # 左上が原点の等高線
plb.axis('equal')               # X-Y軸の増分を同じに
plb.axis('off')                 # 座標軸を非表示
plb.savefig( '1_2_2_1.png' )

plb.figure()                    # 新規図
plb.hist(im.flatten(),128)      # ヒストグラム
plb.savefig( '1_2_2_2.png' )

画像のサイズと、型の抽出

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy
from PIL import Image

im = numpy.array(Image.open('empire.jpg'))
print im.shape, im.dtype        #カラー画像のsizeと型

im = numpy.array(Image.open('empire.jpg').convert('L'),'f')
print im.shape, im.dtype        #白黒画像のsizeと型

画像のグレースケール化と、明度変換

#!/usr/bin/python
# -*- coding: utf-8 -*-

from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb


im = numpy.array(Image.open('empire.jpg').convert('L'))
im2 = 255 - im                  # 画像反転
im3 = (100.0/255) * im + 100    # 100〜200の範囲に縮小?
im4 = 255.0 * (im/255.0)**2     # 2乗

#明度の最大値/最小値 表示
for i in (im,im2,im3,im4):
  print "(MIN MAX)=(", str( int(i.min()) ), str( int(i.max())) ,")"

plb.figure()
plb.gray()
plb.axis('off')
plb.imshow(im)          #画像の登録?
plb.savefig( '1_3_2-1.png' )

plb.figure()
plb.gray()
plb.axis('off')
plb.imshow(im2)         #画像の登録?
plb.savefig( '1_3_2-2.png' )

plb.figure()
plb.gray()
plb.axis('off')
plb.imshow(im3)         #画像の登録?
plb.savefig( '1_3_2-3.png' )

plb.figure()
plb.gray()
plb.axis('off')
plb.imshow(im4)         #画像の登録?
plb.savefig( '1_3_2-4.png' )

PIL / Pillow for python による画像処理入門

www.oreilly.co.jp

と言っても、↑こちらで提供されているサンプルコードのまんま

#!/usr/bin/python
# -*- coding: utf-8 -*-

from PIL import Image

# グレースケール変換
pil_im = Image.open('empire.jpg').convert('L')
# 拡張子に応じた形式で保存
pil_im.save( 'empire_1.jpg' )

# 私のcent7環境にx11はない為、out.show() を行うと、次のようなerror
#   display: unable to open X server
#   `' @ error/display.c/DisplayImageCommand/429.
# out.show()

# サムネイル画像作成
pil_im.thumbnail((128,128))
pil_im.save( 'empire_2.jpg' )

pil_im = Image.open('empire.jpg').convert('L')
# 画像の一部を切取り→反転→貼付け
box = (100,100,400,400)
region = pil_im.crop(box)
region = region.transpose(Image.ROTATE_180)
pil_im.paste(region,box)
pil_im.save( 'empire_3.jpg' )

pil_im = Image.open('empire.jpg').convert('L')
# 画像のサイズ変更
out = pil_im.resize((128,128))
out.save( 'empire_4.jpg' )

pil_im = Image.open('empire.jpg').convert('L')
# 画像全体の回転
out = pil_im.rotate(45)
out.save( 'empire_5.jpg' )

out.show() は使用していません

画面のない私のcentos on virtual boxで、out.show() を使用すると

display: unable to open X server `' @ error/display.c/DisplayImageCommand/429.

のようなエラーとなる為、これをコメントアウトし、out.save()を使用しています。

固有値、固有ベクトルからの主成分分析 オレオレ入門 (PCA: Principal Component Analysis)

主成分分析は、算出された固有値固有ベクトルを優先度付けすることで、 多次元の行列を次元削減する為のものですが、口では上手く説明できなかったのでメモ。

固有値固有ベクトル の算出方法?

end0tknr.hateblo.jp

以前のエントリに記載している通りです。

線形変換における固有ベクトルの特徴

行列は座標の変換に用いられますが、 座標変換を行った後も固有ベクトルは、その向きが変わることはありません。

上記urlにある線形変換の図(GIFアニメーション)が分かりやすいと思います。

deeplearning4j.org

(振動においては「固有振動数」という似た用語?もありますが)

行列は、次数に応じた固有値固有ベクトルを持つ

以前の私のエントリでは、2次行列に対する固有値固有ベクトルの 算出のみ記載していますが、3次行列では3種、4次行列では4種のように、 行列は、次数に応じた固有値固有ベクトルを持ちます。

f:id:end0tknr:20170919224645p:plain

主成分分析は、固有ベクトル固有値順にソートしたもの

例えば、5次行列の固有ベクトルを算出し、 これを固有値順に表示すると次のようになると思います。

f:id:end0tknr:20170919224655p:plain

結局?、主成分分析では、固有値の大きさ(寄与度)を参照しながら、 採用する固有ベクトルを決定することで、次元削減を行います。