テンプレートマッチング(template matching)

画像の中から指定した画像(テンプレート)と似ている位置を探すことを
テンプレートマッチング(template matching)と言います。

 

テンプレート 画像

 

このとき、テンプレートと画像データがどれだけ似ているか?という評価値(類似度または相違度)にはいくつかあり、次に示すような値を用います。

 

以下の式において、テンプレートの輝度値の値をT(i,j)、画像の輝度値の値をI(i,j)とします。
座標の(i,j)はテンプレートの幅をm画素、高さをn画素としたとき、左上(もしくは左下)を(0,0)、右下(もしくは右上)を(m-1、n-1)とします。

 

SSD(Sum of Squared Difference)

SSDはテンプレートをラスタスキャンし、同じ位置の画素の輝度値の差の2乗の合計が用いられます。
SSDの値が小さいほど、似ている位置となります。

 

SAD(Sum of Absolute Difference)

SADはテンプレートをラスタスキャンし、同じ位置の画素の輝度値の差の絶対値の合計が用いられます。
SADの値が小さいほど、似ている位置となります。

 

正規化相互相関【NCC:Normalized Cross-Correlation】

テンプレートと画像との類似度として、以下の正規化相互相関を用いられる場合もあります。類似度が1に近いほど、似ている位置となります。

 

ちなみに、この計算式は内積の式をCosθ = の式に変形した式と同じ事にお気づきでしょうか?上式を以下のように変形すると、M×N次元のIのベクトルとTのベクトルの内積に見えてきませんか?

 

 

正規化相互相関【ZNCC:Zero-mean Normalized Cross-Correlation】

上記NCCの相互相関係数ではテンプレートや画像の明るさが変動すると、NCCの値もふらついてしまうので、テンプレートおよび画像の輝度値の平均値をそれぞれの値から引いて計算することで、明るさの変動があっても安定的に類似度を計算することができます。

ここで、この式はテンプレートの領域内の輝度値の平均値を計算してから、さらに輝度値から平均値を引くため、標準偏差の時にも説明したように、この式でそのままプログラミングすると効率の悪いプログラムとなってしまいます。そのため、今回も同様にRZNCCの式を変形します。

 

テンプレートの平均輝度値およびテンプレートと同じ領域の画像の輝度値の平均は

 

 

で求まることから、この値をRZNCCの式に代入して整理すると、

となります。この式を用いると、プログラム的には1パスで済むでの計算効率が良くなります。

 

ΣΣ×××の部分がなんだか難しく見える方は、C言語に置き換えて

 

Sum = 0;
for(j = 0; j < N; j++){
for(i = 0; i < M; i++){
Sum += ×××;
}
}

 

と置き換えて読んで下さい。

 

ここで説明しておきながら、最近の工業用途では、この正規化相関の方法よりもパターンの輪郭に基づき、回転やスケール変動などにも対応したマッチング(輪郭サーチ)が一般的になってきています。

最近では、さらにカメラが動いてもマッチングする三次元マッチングなども登場してきています。
正規化相関では明るさの変動にも強く優れたアルゴリズムのように聞こえるかもしれませんが、パターンが重なった場合や欠けた場合に意外と弱かったりもします。

 

また、このテンプレートマッチングでは、テンプレートと画像を重ね合わせたとき、同じ画素の位置の輝度値が明暗がどれだけ似ているか?という事を評価しています。
そのため、例えば顔写真を登録しておいて、撮影した画像の中から似ている人を見つけ出す!みたいな使い方はできません。
人間の感性での似ている事とはズレがある事に注意して下さい。

 

画像処理アルゴリズムへ戻る

 

疑似カラー(Pseudo-color)

三次元データやサーモグラフィのようにデータを画像にした場合には、モノクロで表示するよりも、色を付けて表示した方が見やすい場合があります。

 

モノクロ表示 疑似カラー表示

 

このグレーの色に疑似的に色を付ける方法を疑似カラー(Pseudo-color)と言います。この色の付け方は色相を使って青~緑~赤へと変化させてもいいのですが、もう少し簡単な方法を紹介します。
青~緑~赤~青へと色を変化させるには以下のようなパターンでR,G,Bの値を変化させます。

 

 

このパターンは青~緑~赤~青の色相で言うと一周分なので、一般的に用いられるのは青~緑~赤ぐらいまでなので、上図の0~240°部分を8Bitの輝度値に割り当て

 

 

のようにすると、このようなグレースケール↓に

 

疑似カラーを割りつけると、このように↓なります。

 

このようにクレースケールにカラーを割り当てることで、より画像を見やすくしています。モノクロ8Bitの画像データの場合は、画像データはそのままに、カラーテーブルを変更するだけで、疑似カラー表示する事が出来ます。

 

画像処理アルゴリズムへ戻る

 

色相、彩度、明度の計算方法

色相Hue:色合い)、彩度Saturation:鮮やかさ)、明度Brightness,Lightness,Intensity,Value:明るさ)については、以前、変換式には色相、彩度、明度ほかのページにまとめたのですが、実は訳も分からず公式だけをまとめていました。

 

で、なんだか気持ちが悪かったので色相、彩度、明度について、よ~く調べてみました。
私なりの理解ですが、以下にまとめました。

 

R、G、Bの色空間については、下図のようにRGBをXYZのように三次元座標で表すと、一辺の長さが255で表される立方体の範囲内で全ての色を表す事が出来ます。(R,G,B各8bitの場合)

 

この立方体を白(255、255、255)に位置から黒(0、0、0)の方向へ見て、R軸を右側に取ると、

 

 

のように、正六角形となります。
この時、の方向を0°として、反時計回りにの位置が120°240°色相(0~360°(2π))を定めます。
彩度は一番外側の六角形に対して、どの割合の位置に配されているかを0~1.0で表したものが彩度となります。

 

詳細は後述しますが、色相彩度はカメラやパソコンなどの性能評価(使いやすさ、価格、処理速度など)を表す時に用いるレーダーチャート(クモの巣グラフ)もどきみたいな物?!と思うと、自分の中で少し整理ができました。

 

さらに、この六角形の高さ方法に明度を割り振ると、HSV(六角錐モデル)やHLS(双六角錐モデル)となります。
それぞれの違いは明度の定義が異なり、R、G、Bの最大輝度値をImax、最小輝度値をIminとしたときに

明度V = Imax

としたものがHSV

明度L = ( Imax + Imin ) / 2

としたものがHLSとなり、明度の値は0~1.0で表されます。

これを立体で表すと

 

HSV(六角錐モデル)

HLS(双六角錐モデル)

となります。
このHSV、HLSともに、六角錐の斜面の部分が彩度が1.0となります。

 

以下、色相、彩度、明度の詳細な計算方法です。

 

HSVの計算方法

はじめにR、G、Bの輝度値の範囲を0~255から0~1.0となるように変換します。
(R、G、Bのそれぞれの値を255で割ります。)

 

【色相Hの求め方】
下図のように、0°方向にR、120°方向にG、240°方向にBだけ進み、最後の点の位置のR軸に対する角度が色相となります。

 

 

この最後の点の座標は中心を(x、y) = (0、0) とすると、R、G、Bの方向のなす角度から

 

 

となり、xとyより色相Hが求まります。

 

 

ただし、アークタンジェントの計算が出来ない場合など、この方法とは別に、近似的に求める方法もあります。(こちらの方が一般的)

 

下図を見ても分かる?ようにR、G、Bの成分の比を比べ、

Rが最大の場合、色相は-60°(300°)~60° (R方向の0°±60°)
Gが最大の場合、色相は 60°~180° (G方向の120°±60°)
Bが最大の場合、色相は 180°~300° (B方向の240°±60°)

の範囲内に色相は収まります。

 

 

以下、Rの値が最大の場合を例に取って説明したいと思います。

 

下図のように、2つの矢印の長さが分かれば、その矢印の比で角度60°を分割することで、角度(色相)を近似することが出来ます。

 

 

R、G、Bの大きさがR≧G≧Bの場合

 

色相H = 60° × (G – B) / (R – B)

R、G、Bの大きさがR≧B≧Gの場合

 

色相H = 60° × (G – B) / (R – G)

 

となります。
ただし、この場合、色相の値が負となるので、

 

色相H = 60° × (G – B) / (R – G) + 360°

 

とします。

 

と、なる理屈を理解するのに苦労しました...
図中に書いてある黄色い正三角形がポイント!
正三角形なので、三辺の長さが等しい分けで。
R以外のGやBが最大となる場合も理屈は同じです。
120°づつ回転させて考えてみると分かります。

 

この式を一般的に書くと、R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

ImaxがRのとき

ImaxがGのとき

ImaxがBのとき

となります。

 

【明度Vの求め方】

明度は、もともとHSVの定義よりR、G、Bの成分のうち、最大の成分をImaxとすると

明度V = Imax

 

とします。
明度Vの範囲は0~1.0となります。

 

【彩度Sの求め方】
R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

 

彩度S = (Imax – Imin) / Imax

 

となります。
彩度Sの範囲は0~1.0となります。

 

 

HLSの計算方法

【色相Hの求め方】
色相HはHSVの色相Hの求め方と同じです。

 

【明度Lの求め方】
明度Lは、もともとHLSの定義よりR、G、Bの成分のうち、最大の成分をImaxとすると

 

明度L = ( Imax + Imin ) / 2

 

とします。
明度Lの範囲は0~1.0となります。

 

【彩度Sの求め方】
R、G、Bの成分のうち、最大の成分をImax、最小の成分をIminとすると

 

L≦0.5のとき

彩度S = (Imax – Imin) / (Imax + Imin)

L>0.5のとき

彩度S = (Imax – Imin) / (2 – Imax – Imin)

 

となります。
彩度Sの範囲は0~1.0となります。

 

以下、補足説明です。
彩度Sは下図の外側の六角形に対して、内側の六角形の大きさの割合で求められます。

 

 

この六角形の大きさはR,G,Bの輝度値が最大となる軸上で考えると比較的分かりやすいと思います。
今回はRの値が最大となる場合とします。

 

HLSは双六角錐モデルであるため、明度Lが0.5以下の場合、外側の六角形の大きさ(上図のE’の位置)は

 

E’ = Imax + Imin

 

で求まります。

 

 

明度Lが0.5より大きい場合、外側の六角形の大きさ(上図のE’の位置)は

 

E’ = 2 – Imax – Imin

 

となります。

 

 

以上のことから、最初の彩度Sの式が求まります。

 

色相、彩度、明度を使った色判別時の注意点

色相および彩度を用いると、画像の明るさ(明度)が変動しても似た色の領域を抽出する事が可能となりますが、彩度の値が小さい場合、つまりR、G、Bの値がそれぞれ近い場合は色相の値が不安定になります。

 

例えば、
(R、G、B) = (121、120、120)の場合、 色相H = 0°
(R、G、B) = (120、121、120)の場合、 色相H = 120°
(R、G、B) = (120、120、121)の場合、 色相H = 240°

 

と、ほんの少しのR、G、Bの値の違いでも色相の値は大きく異なります。

 

また、色相Hは角度で表されるので、例えば1°も359°も値こそ離れていますが、どちらも0°±1°の範囲内で角度的には近いので、色相Hの値で単純に二値化処理することで色の領域を抽出する場合は注意して下さい。

 

画像処理アルゴリズムへ戻る

 

色相、彩度、明度の公式

カラーの画像処理をする時には、これら色相などの知識は必須となります。
Windows標準で付いてくるペイントで、色の作成の表示をすると、雰囲気が分かると思います。

 

 

基本的に以下の色相、彩度、明度を用いて色を表すのですが、変換式にいくつかの種類があります。

 

色相(Hue)

色合いを表します。
赤や緑、青などに色を0~360°(0~2π)の角度を用いて表します。

彩度(Saturation)

鮮やかさを表します。
と一般的に言われるのですが、鮮やかさ?と言われても、いまいちピンと来ませんが、下記に示した式から見ても分かるように、R,G,Bの値にどれだけ開きがあるか?を示しています。
このことは逆にいうと、R,G,Bの値に開きが無い場合は、グレーに近い事から、彩度は如何にグレーっぽく無いか?、という事から、どれだけ純色(赤、緑、青、黄、シアン、紫など)に近いか?を表しています。

明度(Brightness,Lightness,Intensity,Value)

色の明るさを表します。

 

HSV変換

6角錐モデルとも言います。

 

【RGB⇒HSV変換】

Imax = Max(R,G,B)
Imin = Min(R,G,B)
とすると

 

R = Imaxのとき

H = 60×(G – B) / (Imax – Imin)

G = Imaxのとき

H = 60×(B – R) / (Imax – Imin) + 120

B = Imaxのとき

H = 60×(R – G) / (Imax – Imin) + 240


S = (Imax – Imin) / Imax


V = Imax


【HSV⇒RGB変換】

h = floor(H / 60)     floor()は切り捨て処理
P = V × (1 – S)
Q = V × (1 – S × (H / 60 – h))
T = V × (1 – S × (1 – H / 60 + h))
とすると

 

h= 0のとき

R = V, G = T, B = P

h= 1のとき

R = Q, G = V, B = P

h= 2のとき

R = P, G = V, B = T

h= 3のとき

R = P, G = Q, B = V

h= 4のとき

R = T, G = P, B = V

h= 5のとき

R = V, G = P, B = Q


HLS変換

双6角錐モデルとも言います。

 

【RGB⇒HLS変換】

Imax = Max(R,G,B)
Imin = Min(R,G,B)
とすると

 

R = Imaxのとき

H = 60×(G – B) / (Imax – Imin)

G = Imaxのとき

H = 60×(B – R) / (Imax – Imin) + 120

R = Imaxのとき

H = 60×(R – G) / (Imax – Imin) + 240


L = (Imax + Imin) / 2

L ≦ 0.5のとき

S = (Imax – Imin) / (Imax + Imin)

L > 0.5のとき

S = (Imax – Imin) / (2 – Imax – Imin)


【HLS⇒RGB変換】

h < 0のとき

h’ = h  + 360

h ≧ 360のとき

h’ = h  – 360

その他

h’ = h

L ≦ 0.5のとき

M2 = L × (1 + S)

L > 0.5のとき

M2 = L + S – L × S


M1 = 2 × L – M2


h’ < 60のとき

X = M1 + (M2 – M1) × h’ / 60

60 ≦ h’ < 180のとき

X = M2

180 ≦ h’ < 240のとき

X = M1 + (M2 – M1) × (240 – h’ ) / 60

240 ≦ h’ ≦ 360のとき

X = M1

 

とすると

 

R = X ただし、h = H + 120とする

G = X ただし、h = Hとする

B = X ただし、h = H – 120とする

 

カラー変換用関数

【Win32APIの場合】 VBの表記例
‘HLS変換(Windows 2000以降、またはInternet Explorer 5.0がインストールされてある環境。(SHLWAPI.DLL Version 5.00以上)
‘h (色相)
‘赤(0)、黄(40)、緑(80)、シアン(120)、青(160)、マゼンダ(200)の順に定義0~239まで設定可
‘L (明度)
‘色の明るさをあらわす。0~240まで設定可。0が黒、240が白になる。
‘s (彩度)
‘0~240まで設定可。240が純色になる。
Public Declare Sub ColorRGBToHLS Lib “SHLWAPI.DLL” _

(ByVal clrRGB As Long, _
pwHue As Integer, _
pwLuminance As Integer, _
pwSaturation As Integer)

Public Declare Function ColorHLSToRGB Lib “SHLWAPI.DLL” _

(ByVal wHue As Integer, _
ByVal wLuminance As Integer, _
ByVal wSaturation As Integer) As Long

 

最初に紹介したペイントの色の作成では、この関数と同様の変換(設定値を含めて同じ)をしています。

 

【.NET Frameworkの場合】
System.Drawing.Color構造体にて

GetHueメソッド      HSBのH(色相)を取得
GetSaturationメソッド   HSBのS(彩度)を取得
GetBrightnessメソッド   HSBのB(明度)を取得

RGB⇒HSB変換はなし?

 

【OpenCVの場合】
cvCvtColor関数にて、以下の変換に対応
XYZ, YCrCb(YCC), HSV, HLS, L*a*b, L*u*v
その他 同じcvCvtColor関数で Bayer変換に対応

 

カラー画像処理例

そのカラー画像を色相、彩度、明度に分解し、それぞれの値(主に色相)でフィルタリング処理(バンドパスフィルタ)を行い、カラー画像に逆変換する事により、特定の色だけを抽出する事が可能になります。これにより、色の位置や個数などの検査をする事が可能となります。

処理前 処理後

ただし、彩度の値の小さな色(白やグレーに近い色)は彩度の値が不安的になりがちなので、カラー画像処理には不向きです。

 

画像処理アルゴリズムへ戻る

 

外周画像の処理

平滑化フィルタやメディアンフィルタなどの注目画素の周辺画素を用いた画像フィルタ処理では、画像の外周部分が下図のように画像の外側を参照してしまうため、処理ができなくなります。

 

 

この外周部分を処理する方法はいくつかあるのですが、代表的な方法を紹介します。
以下、5×5サイズのカーネルを用いた場合の処理を例にとって紹介します。

 

■外周部分の輝度値を画像の外側にコピーして補間する方法

5×5サイズのカーネルの場合、画像の外側に2画素分、画像の輝度値を参照
してしまうので、この2画素分の輝度値を画像の外周部分の輝度値をコピーして
輝度値を参照します。

 

おそらく?この手法が一般的だと思います。

 

■外周部分を中心にして対称の位置にある輝度値を外側にコピーして補間する方法

例えば、カーネルが座標(-1、-2)の画素を参照する場合は座標(1、2)の画素の輝度値を参照するようにします。

 

 

他にも、単純に外周部分は処理をしないで黒(輝度値=0)でマスクしたり、カーネルが画像の外側を参照する場合にカーネルを形を変えて画像の外側を参照しないようにしたりする方法などもあります。

 

画像処理アルゴリズムへ戻る

 

画像の回転

画像を回転する場合、任意点周りの回転移動でも紹介したように回転行列を使って、例えば、画像の中心周りに画像を回転させると、下図のように回転後の画像が虫食い状態になってしまいます。

 

回転前の画像 回転後の画像

 

こうならないようにするためには、回転前の画像の座標を回転行列を使って回転後の座標を計算するのではなく、回転後の座標が回転前の画像のどの座標を参照しているのかを計算し、画像を変換します。

つまり、画素を置きに行くのではなく、拾いに行くようにします。

 

【回転前の画像を回転行列を使って変換】

 

【回転後の座標が回転前のどの座標を参照しているかを計算して変換】

 

実際の変換処理は以下のように行います。

 

回転前の画像の座標を(x,y)、回転後の画像の座標を(X,Y)、画像の中心座標を(Cx,Cy)とすると単純に画像の中心周りに座標を回転すると以下のような行列で表されます。

 

 

この行列を回転前の画像の座標(x、y)に関して解けばいいので、行列の式の両辺にそれぞれ逆行列をかければいいので、以下のような手順で行列を解いていきます。

 

これで、回転前の座標(x、y)に関して解くことができます。
でも、逆行列を解くのは面倒と思う事なかれ。
(+Cx,+Cy)の平行移動の逆行列は(-Cx、-Cy)方向への移動と同じ、+θ方向への回転の逆行列は-θ方向への回転と同じなので、

 

 

とすれば、逆行列を解くことなく回転後の座標(x、y)に関して座標を解くことができます。
このようにして画像の回転処理を行うと、このようになります。

 

回転前の画像 回転後の画像

 

でも、回転後の画像はなんかギザギザしてしまっていますが、これは回転前の画像の座標を計算するとほとんどの場合、画素の画素の間の座標となってしまいますが、上図の例ではこの座標を四捨五入して輝度値を参照しているためで、bilinearやbicubicなどの補間を使うと少しは滑らかな画像となります。

今回は画像の回転について紹介していますが、画像の拡大縮小についても同様の考え方で処理することができます。

 

画像処理アルゴリズムへ戻る

 

関連記事

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

任意点周りの回転移動(アフィン変換)

画素の補間(Nearest neighbor,Bilinear,Bicubic)の計算方法

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

画像の拡大縮小、回転、平行移動などを行列を使って座標を変換する事をアフィン変換と呼びます。

X,Y座標の二次元データをアフィン変換するには、変換前の座標を(x, y)、変換後の座標を(x’,y’)とすると回転や拡大縮小用の2行2列の行列と、平行移動用に2行1列の行列を使って

$$\left(\begin{array}{c}x^{‘}\\ y^{‘}\end{array}\right)=
\left(\begin{array}{c}a & b\\ c & d\end{array}\right)
\left(\begin{array}{c}x\\ y\end{array}\right)
+ \left(\begin{array}{c}T_{x}\\ T_{y}\end{array}\right)$$

のように表現される場合もありますが、回転、拡大縮小、平行移動を1つの3x3の行列にまとめて

$$\left(\begin{array}{c}x^{‘}\\ y^{‘} \\ 1\end{array}\right)=
\left(\begin{array}{c}a & b & c\\ d & e & f\\ 0 & 0 & 1\end{array}\right)
\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)$$

と3x3の行列で表現する場合もあります。

この表現を同次座標系と呼びます。

同次座標系では一見、3行目は無駄なようにも見えるのですが、この意味の無いような1行を追加する事で、平行移動も同じ行列の積で表現でき、逆行列を使う事で、アフィン変換後の座標からアフィン変換前の座標も簡単に求める事ができるようになります。

私は3行3列の行列を用いた同次座標系のアフィン変換しかしていない、というか断然おススメなので、同次座標系で説明したいと思います。

後半でアフィン変換の実用例を示しているので、その部分で、同次座標系の恩恵を感じてもらえると嬉しいです。

 

変換前の画像を以下のようにすると、

各種変換は以下の通りとなります。

 

拡大縮小

X軸方向の拡大率をSx、Y軸方向の拡大率をSyとすると拡大縮小のアフィン変換は

 

 

と表されます。

 

例)X軸方向に2倍

 

例)Y軸方向に2倍

 

例)X軸、Y軸方向に2倍

 

例)Y軸方向に-1倍

このように、ある軸(上記の例ではX軸)に対して反転する処理の事を鏡映と呼びます。

 

平行移動

X軸方向にTx、Y軸方向にTyだけ移動するアフィン変換は

 

 

のように表されます。

 

 

回転

原点を中心に反時計回りにθ°回転する時のアフィン変換は

 

 

のように表されます。

 

 

スキュー(せん断)

四角形の画像を平行四辺形に変形する処理をスキューまたはせん断といいます。
このアフィン変換は

 

 

アフィン変換の実用方法

画像処理で使われるアフィン変換は、下図のようにピクチャボックスなどの左上が原点[座標が(0, 0)]で右方向が+X、下方向が+Y、時計周りが+θ方向となる場合が多いかと思います。

また、平行移動、拡大縮小、回転などの行列を紹介しましたが、これらの行列を1回だけで処理する事はまれで、それぞれの行列を組み合わせてアフィン変換を行います。
さらに、拡大縮小、回転、スキューのアフィン変換行列は、あくまでも原点を基点として、変換される事に注意が必要です。

例えば、画像が原点の位置に無い場合、画像をX,Y方向に2倍の大きさにしようとしたとき、単に拡大縮小のアフィン変換行列で座標を計算すると、表示位置も2倍、原点の位置から離れた位置に移動します。

これらの事を踏まえ、画像を幅方向に2倍、高さ方向に3倍し、画像を反時計方向に90°回転、さらに、画像の左上の座標が(50, 50)から(30,180)へ移動するアフィン変換を例題にとって考えたいと思います。

この変換は、一発で変換するアフィン変換行列を考えるのではなく、平行移動拡大縮小回転に分けて、アフィン変換の順番を考えます。


拡大縮小と回転が原点を基点とするため、画像を原点の位置へ移動するため、x方向に-50、y方向に-50の平行移動します。

平行移動のアフィン変換行列は

$$\begin{pmatrix} 1 & 0 & -50 \\ 0 & 1 & -50 \\ 0 & 0 & 1 \end{pmatrix}$$


拡大縮小か回転のどちらからでも構いませんが、x方向に2倍、y方向に3倍の拡大縮小を行います。

拡大縮小のアフィン変換行列は

$$\begin{pmatrix} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}$$


反時計周りに90°回転(ー90°回転)を行います。

回転のアフィン変換行列は

$$\begin{pmatrix} cos(-90) & -sin(-90) & 0 \\ sin(-90) & cos(-90) & 0 \\ 0 & 0 & 1 \end{pmatrix}$$


最後に目的の位置へ移動するのに、x方向に+30、y方向に+180の平行移動します。

平行移動のアフィン変換行列は

$$\begin{pmatrix} 1 & 0 & 30 \\ 0 & 1 & 180 \\ 0 & 0 & 1 \end{pmatrix}$$


このようにアフィン変換を平行移動、拡大縮小、回転に分解して、変換の手順を考える事が大事です。
今回の場合は

変換前 → 平行移動 → 拡大縮小 → 回転 → 平行移動 → 変換後

としています。

画像に対してアフィン変換を行う場合、考え方としては、平行移動や拡大縮小などに分解して考えますが、画像データを毎回変換するのではなく、アフィン変換行列をまとめて計算し、一発でアフィン変換処理を行います。

具体的には、今回のアフィン変換処理を行列で表すと

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 1 & 0 & 30 \\ 0 & 1 & 180 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} cos(-90) & -sin(-90) & 0 \\ sin(-90) & cos(-90) & 0 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & 0 & -50 \\ 0 & 1 & -50 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

行列部分をまとめて計算し

$$\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

となり、アフィン変換行列を一発で処理する事ができます。

さらに、変換後の座標(x’, y’)から、変換前の座標(x, y)を求める場合があるのですが、その時は、行列の左側からアフィン変換行列の逆行列を掛けて、

$$\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$$

$$\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & 3 & -120 \\ -2 & 0 & 280 \\ 0 & 0 & 1 \end{pmatrix}^{ -1 }\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}$$

$$\begin{pmatrix} x \\ y \\ 1 \end{pmatrix}=\begin{pmatrix} 0 & -0.5 & 140 \\ 0.333 & 0 & 40 \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} { x }^{ ‘ } \\ { y }^{ ‘ } \\ 1 \end{pmatrix}$$

として、変換前の座標を求める事ができます。

まとめ

  • アフィン変換は平行移動、拡大縮小、回転、スキューに分けて、変換の順番を考える
  • 実際のアフィン変換は、アフィン変換行列をまとめて計算し、一発で処理を行う
  • 変換前の座標は逆行列を使うと、求めることができる

この順番の考え方は、ほとんどの場合、画像を原点へ平行移動し、拡大縮小、回転を行ってから、目的の位置へ移動すると、アフィン変換行列が求まりますが、他にも、特定の点を基準に拡大縮小や回転を行う場合は、その点を原点へ移動すればアフィン変換行列が求まります。

(参考)

 

注意点

アフィン変換では任意の3×3(2×3)の行列で表す事ができるので、任意形状に変換できそうにも思えるのですが、四角形が平行四辺形にまでは変形できるものの、台形には変形できないのでご注意願います。
この台形に変形できる処理は射影変換(ホモグラフィ)と呼びます。

アフィン変換は今回の説明のように、画像を移動、変形させるための手法として説明されますが、もう少し汎用的に座標変換として捉えると応用範囲が広がります。

例えば、データのグラフを表示する時に、横方向、縦方向に拡大/縮小した時に、表示するデータの範囲を求める場合などに応用すると、少し便利です。

 

画像処理アルゴリズムへ戻る

 

関連記事

実際にOpenCVを使って行うアフィン変換については、こちらのページ↓にまとめました。

【OpenCV-Python】アフィン変換(同次座標系を使用)

 

マイクロソフトの.NETやDirectXで扱うアフィン変換行列は、行と列が逆になり、行列を掛ける順番も逆(右側から掛ける)になります。
この仕様については、下記ページ↓にまとめました。

【C#.NET】マイクロソフト仕様のアフィン変換

円形度

輪郭追跡を行うと周囲長を求めることができますが、この周囲長を用いた代表的な特徴量の円形度を紹介します。

 

円形度とは円らしさを表す値で値が1となる時、もっとも円に近くなります。
定義は面積(画素数)をS周囲長Lとすると、

 

円形度 = 4πS/L2

 

となります。

 

【円形度の例】

面積S = 96
周囲長L = 34.627円形度 = 1.006
面積S = 87
周囲長L = 33.456円形度 = 0.977

なぜそうなるのか?理屈は簡単で、
円の半径とすると、円の面積Sπr周囲長L2πrより、

 

円形度 = 4πS / L2 = 4π(πr2) / (2πr)2 = 1

 

となり、円形度4πS/L2に近いときもっとも円らしいということになります。

 

画像処理アルゴリズムへ戻る

関連記事

【OpenCV-Python】円形度

輪郭追跡処理アルゴリズム

二値化された画像において、各連結部分の境界部分を求める事を輪郭追跡といいます。
この輪郭追跡を行う事で周囲長や連結領域の高さなどを求めることが出来ます。

 

以下、この輪郭追跡の処理アルゴリズムを紹介します。

 

まず、左上からラスタスキャンを行い、白の部分を見つけます。

 

 

最初の白の部分(上図の赤丸の部分)をスタート地点として、反時計回りに輪郭の部分を検出します。
最初の白の部分は左上からラスタスキャンしたため、スタート地点の右上、上、左上、左には白の部分が無い事が確定しているので、最初に左下の部分から反時計回りに白の部分を検索し、白の部分が見つかったら次の輪郭へ検索を続けます。
最後にスタート地点に戻ったら処理は終了です。

 

この時、輪郭の向きを下図のように定義しておくと、

 

最初の図の輪郭は

 

0→0→1→0→2→3→4→3→2→5→5→6→4→0→6→5

 

の順で変化しています。(この番号はチェインコードといいます。)
ここにちょっとした規則があります。
例えば、の次に5、1の次に6、・・・などは絶対に来ません。(なぜかは、ちょっと考えると分かります。)
一般的に表現すると、
一つ前の輪郭の向きから時計回りに3つ分の向きには輪郭が存在しない。
という事になります。
この性質を利用して、輪郭を検索するときは、
一つ前の輪郭の向きの値をVold
輪郭の検索を開始する向きの値を Vnew
とすると、
Vnew = (Vold + 6) % 8; (%は8で割った時の余りの計算)
となるVnew の向きから輪郭を検索すれば効率的に輪郭を検索することが可能となります。

 

また、輪郭の周囲長を求める時、縦、横方向の輪郭は長さ、斜め方向の輪郭は√2 として
長さを足していくと、周囲長が求まります。

 

そして、ここからがあまり書かれていない注意点!
スタート地点の白の画素のパターンが以下のパターンとなる場合、スタート地点を2回通る事になるので、ご注意下さい。

※グレーの部分は白でも黒でもよい部分

 

画像処理アルゴリズムへ戻る

 

ラベリング

二値化画像処理された画像において、白の部分(または黒の部分)が連続した画素に同じ番号を割り振る処理を
ラベリングと言います。

 

通常、同じ番号ごとの面積(画素数)や幅、高さなどの特徴量を求めて欠陥検査や分類処理などに用いられます。

 

ラベリングには、二値化された画像の縦、横方向に連続している部分を同じラベルにする4連結と、
縦、横、斜め方向に連続している部分を同じラベルにする8連結との2種類の処理があります。
(4近傍、8近傍という場合もあります。)

【二値化画像】
【4連結】 【8連結】

以下、8連結の場合において、ラベリング処理アルゴリズムを紹介します。

 

まず、画像全ての画素のラベル番号を0(ゼロ)で初期化しておき、ラベリングで番号を割り付けるための
ラベリング番号のルックアップテーブルを用意しておきます。(テーブルの使い方の詳細は後ほど)

 

そして、画像の左上からラスタスキャンを行い、画素の色がの位置を検索します。

 

白の画素の左上、上、右上、左の画素のラベル番号を参照し、全て0(ゼロ)の場合は、最後に割り振った番号+1のラベル番号を割り振ります。
もし、参照した画素のラベル番号が複数存在した場合は、最小の番号を割り振ります。

 

最初に番号を割り振った様子↓

 

この処理を左上からラスタスキャンして続けていきます。

 

新しい番号が割り振られるとき↓

 

新しい番号が割り振られるとき↓ その2

 

参照した画素のラベル番号が複数存在した場合、最小の番号を割り振ります。
このとき、用いなかったラベル番号(下図の例では3)のルックアップテーブルの番号を最小の番号に書き換えます。

 

同様に、処理を続けていきます↓

 

全ての白の画素に番号が割り振れれた様子↓

 

ただし、よ~く見ると、ラベル番号の2や3など、連続していない部分が存在している事に気がつきます。
そこで、ルックアップテーブルの出番です。

 

処理の途中でルックアップテーブルの番号を変更した番号を
2→1、3→1、6→5
に変更すると、連続した領域は同じ番号になります。

 

このまま番号を変更しても良いのですが、ラベル番号が虫食いになるのも気持ちが悪いので、以下のようにルックアップテーブルの番号を詰め直してから、ルックアップテーブルを参照し、全画素のラベル番号を修正すると、以下のようになります。

 

基本的なアルゴリズムは上記、紹介したアルゴリズムで良いのですが、ラベリング処理を行う画像のパターンによってはうまくラベル番号が割り振られない場合があります。
そこは、何とかして(ちょっと説明が難しい部分)うまく割り振られるようにチャレンジしてみて下さい。

 

参考までに、ラベリング処理をすると失敗しやすい特徴的な画像を添付しておきます。↓

◆Level.1

 

◆Level.2

 

◆Level.Max!?

 

ラベリング処理のアルゴリズムには、この手法の他にも輪郭線を追跡し、閉じた領域を1つのラベルとする手法や、ランレングス的に求める方法もあります。

 

画像処理アルゴリズムへ戻る

 

細線化

二値化された画像において、線の中心1画素分だけを残すように線を細くする処理を細線化といいます。細線化のアルゴリズムにはHilditch、田村の方法、Zhang Suenなど、いろいろあるのですが、ここでは田村の方法について説明します。
細線化は黒の線を細くする場合と白の線を細くする場合が考えられますが、以下では白の線を細くするのを前提として説明しますので、黒の線を細くする場合は白と黒を読み替えて下さい。

 

細線化処理前 細線化処理後

 

ざっくり言うと、収縮処理をすると線幅が細くなりますが、線の長さが短くなったり、細い線は消えてしまったりするので、線の端点や交点、線幅が1となった画素などを除去しないように条件を付けながら収縮処理を行います。

 

処理をフローチャートを以下に示します。

 

 

画像全体にラスタスキャンを行い、注目画素の3×3画素の並びがパターン1に該当する場合は中心の画素を除去(白から黒に変える)します。ただし、除去しないパターンに該当する場合は除去を行いません。

パターン1に該当する画素が1画素も無ければその時点で終了です。

次にパターン1の時と同様に3×3画素の並びがパターン2に該当する場合は中心の画素を除去(白から黒に変える)します。ただし、除去しないパターンに該当する場合は除去を行いません。

パターン2に該当する画素が1画素も無ければその時点で終了です。そうでない場合はパターン1の処理に戻り、除去する画素がなくなるまでパターン1、パターン2の除去処理を繰り返します。

 

 

 

画像処理アルゴリズムへ戻る

 

膨張・収縮・オープニング・クロージング

膨張・収縮処理では一般的に二値化された白黒の画像に対して処理が行われ、

注目画素の周辺に1画素でも白い画素があれば白に置き換える処理を膨張(Dilation)
逆に周辺に1画素でも黒い画素があれば黒に置き換える処理を収縮(Erosion)といいます。

 

【元画像】

 

【膨張処理】

 

【収縮処理】

二値化された画像ではなく、グレー画像に対して処理を行う場合は、膨張の場合、注目画素の近傍の最大輝度値を注目画素の輝度値に置き換えます。
収縮の場合は最小輝度値に置き換えることでグレー画像に対して処理を行います。
この処理を二値化画像と区別するために、それぞれ最大値フィルタ最小値フィルタと呼ぶ場合もあります。

 

【処理例】

元画像 膨張 膨張
膨張 収縮 収縮
収縮 収縮 収縮
収縮 膨張 膨張
膨張 連続処理

膨張・収縮処理は単独で処理を行う事はまれで、膨張・収縮を繰り返し処理を行う場合が多くあります。
とくに、同じ回数分だけ膨張して収縮する処理をクロージング(Closing)

同じ回数分だけ収縮して膨張する処理をオープニング(Opening)とよびます。

 

元画像 クロージング
元画像 オープニング

 

この処理例を見てもわかるように、オープニング、クロージング処理では小さいパターン細いパターンの除去を行います。

 

さらに、クロージングした画像から元画像を差し引いた処理をブラックハット(Black-Hat)
元画像からオープニングした画像を差し引いた処理をトップハット(Top-Hat)といいます。

 

元画像 ブラックハット
元画像 トップハット

 

これらの処理を見ても分かるように、膨張・収縮処理を用いることで、欠け、断線、ホコリやゴミといった欠陥の検査や、細い線状のパターンの文字や線などの抽出への応用が考えられます。

 

画像処理アルゴリズムへ戻る

 

アンシャープマスキング(鮮鋭化フィルタ)

アンシャープマスキング【unsharp masking】とは、画像のボヤけた輪郭を強調するフィルタ処理です。

アンシャープマスキング処理前 アンシャープマスキング処理後

 

【処理アルゴリズム】

元画像データ
元画像データの平滑化処理を行い、元画像から平滑化データを差し引きます。
差し引いた分を元画像に上乗せします

 

この処理により、平滑化処理でボヤけてしまう分だけ、逆に輪郭をくっきりさせる事ができます。

 

この処理をカーネルで表現すると

となります。
kは任意レートで、値が大きいほど、輪郭が鮮明になる効果が高くなります。
(後半の元画像から平滑化画像を引いている部分の処理はラプラシアンフィルタと呼ばれます。)

 

今回は3×3のマスクサイズを例にとって紹介しましたが、マスクサイズを大きくするとマスクサイズよりも小さなパターンを強調する効果を得る事もできます。

 

 

画像処理アルゴリズムへ戻る

 

エッジ抽出(Sobel,Prewittなど)

下図のように画像の輪郭(エッジ)を抽出する方法を紹介します。

エッジ抽出前 エッジ抽出後

上図はソーベルフィルタの例です。

 

微分フィルタ

画像の輝度値に対して、隣り合う画素の輝度差が大きいほど、画像のエッジだと考えることができ、
画像の輝度値データ 画像輝度値データに対して

 

横方向の差分(偏微分)以下のカーネルで求めることができる。

0 0 0
0 -1 1
0 0 0

 

縦方向の差分(偏微分)下のは以下のカーネルで求めることができる。

0 1 0
0 -1 0
0 0 0

または

0 0 0
0 -1 0
0 1 0

※画像データの並びがボトムアップかボトムダウンかにより使い分けます。

 

このとき、エッジの強さは

 

エッジの傾きは

として求めることができる。

 

プリューウィットフィルタ(Prewitt filter)

微分フィルタではノイズに弱いため、横方向(または縦方向)のエッジを計算してから
縦方向(または横方向)に平均化処理を行う手法です。カーネルには以下とおりです。

 

横方向の差分

-1 0 1
-1 0 1
-1 0 1

 

縦方向の差分

1 1 1
0 0 0
-1 -1 -1

または

-1 -1 -1
0 0 0
1 1 1

 

ソーベルフィルタ(Sobel filter)

プリューウィットフィルタではノイズを除去するのに平滑化処理を行っていましたが、ガウシアン平滑化処理を行ったのがソーベルフィルタとなります。カーネルには以下とおりです。

 

横方向の差分

-1 0 1
-2 0 2
-1 0 1

 

縦方向の差分

1 2 1
0 0 0
-1 -2 -1

または

-1 -2 -1
0 0 0
1 2 1

 

画像処理アルゴリズムへ戻る

 

メディアンフィルタ

平滑化フィルタでは注目画素の周辺画素の輝度値を平均し、ノイズを除去していましたが、画像の輪郭もボケてしまう欠点がありました。それに対し、メディアンフィルタでは周辺輝度値の大きさを順に並べ、メディアン(中央値)を注目画素に置き換えることでノイズを除去します。
特に周辺画素の輝度値よりも大きく異なるノイズ(ゴマ塩ノイズとかスパイクノイズという)を除去するのに効果を発揮します。

 

メディアンフィルタ処理前 メディアンフィルタ処理後

輝度値の3D表示

輝度値の3D表示

 

処理の詳細


 

注目画素(画像中央の輝度値165の部分)周辺の輝度値を取得します。
61、96、41、57、165、34、24、30、31

 

この輝度値を順番に並べます。

 

24、30、31、34、41、57、61、96、165

 

並べた輝度値のメディアン(中央値)の41で輝度値165を置き換えます。

この処理を全画素について行うと、ノイズを除去することができます。

 

 

ただし、平滑化処理に比べ、処理が重い...

 

もう少し具体例でいうとこんな感じ

【メディアン処理前】

 

【メディアン処理後】


画像処理アルゴリズムへ戻る

 

ガウシアンフィルタの処理アルゴリズムとその効果

移動平均フィルタでは注目画素周辺の輝度値を単に平均していましが、一般的な画像では 注目画素に近い画素の輝度値は注目画素の輝度値と近い場合が多いですが、注目画素から遠くなればなるほど、注目画素の輝度値とは差が大きくなる場合が多くなります。
この事を考慮し、注目画素に近いほど、平均値を計算するときの重みを大きくし、遠くなるほど重みを小さくなるようにガウス分布の関数

を用いてレートを計算しているのがガウシアンフィルタです。
σの値が小さいほど平滑化の効果は小さくなり、大きいほど効果が大きくなりますが、
よく以下のカーネルが用いられます。

 

3×3の場合

 

5×5の場合

 

ガウシアンフィルタにはローパスフィルタと同様の効果があるそうです。

ということで、フィルタ処理した画像をフーリエ変換し、確かめてみました。

 

オリジナル(処理前)の画像はこちら↓

オリジナル画像 二次元フーリエ変換画像

 

上記の画像にノイズ除去系のガウシアンフィルタ、移動平均フィルタ、メディアンフィルタの
処理を行い、二次元フーリエ変換を行い確認してみました。

 

上記画像に5×5のガウシアンフィルタ処理を行うと

5×5のガウシアンフィルタ処理 二次元フーリエ変換画像

 

5×5の移動平均処理を行うと

5×5の移動平均処理 二次元フーリエ変換画像

 

5×5のメディアンフィルタ処理を行うと

5×5のメディアンフィルタ処理 二次元フーリエ変換画像

 

確かに処理結果を見てみると、ガウシアンフィルタが最も高周波成分を除去できているように
思います。

 

なぜ、そうなるのか?
本には数式においても、この効果が証明できるような事が書いてありましたが、ちょっと難しいので、
カーネルの値について見てみたいと思います。

 

そもそも、ある特定の周期の成分を消すためには、どうすれば良いか?というと、
周期の半分の離れた2点のデータを平均していけば、その周波数の成分を消す事が出来ます。

 

 

ここで、画像で表すことのできる最も高い周波数は

 

 

となる2画素周期のパターンで、このパターンを消すためには、隣り合う2画素の輝度値を
平均すると、画像データから最も高周波の成分を除去する事が出来ます。
しかし、隣り合う2画素の平均の結果は、画素間の位置の輝度値を示してしまうので、
この平均のさらに平均値を取ります。

 

連続する5つの輝度値の値をI-2~I2とすると、平均した値の平均値は

 

となり、なんとこの結果が3×3のガウシアンフィルタの横方向(縦方向)の係数と等しくなります。
さらに、平均の平均を取ると、5×5のガウシアンフィルタの係数と同じ。
7×7のガウシアンフィルタまで確認してみましたが、この平均の平均で表すことが出来るようです。

 

という訳で、ガウシアンフィルタの式

 

 

で表される係数なら何でも良いと思いつつも、

 

 

の係数が良く使われるのかも?しれませんね。

 

さらにこの係数、半分の半分で求めるのは、ちょっと面倒。
これがまた都合よく、(a+b)nの式を展開したときの係数に使われる
パスカルの三角形
の値

 

 

と一次元のガウシアンフィルタの係数と等しくなります。(少なくとも7×7ぐらいまでは)
これを以下のように二次元的に掛け合わせて、係数の合計で割ればガウシアンフィルタの係数となります。
例えば7×7のガウシアンフィルタの係数では下図のようになります。

 

 

この値を係数の合計(=4096)で割ればガウシアンフィルタの係数の出来上がり!!!
(このときσの値は約1.3となります。)

 

画像処理アルゴリズムへ戻る