6. 画像特徴量

画像に何が写っているのかを認識するためには、画像特徴量が良く用いられる。画像特徴量、特に「局所」画像特徴量とは、画像のある画素を中心とした周辺領域の持つ画像的な特徴をベクトルとして表したものである。

画像特徴量には非常に多くの種類があり、その全てを本節で紹介することは難しい。そこで、本項では画像認識に用いられる代表的な特徴量であるSIFT (scale-invariant feature transform)を中心に、その技術的背景から応用までを解説する。

参考資料

• 『ディジタル画像処理』 P.232 第11章 パターン・図形・特徴の検出とマッチング

6.1. 画像特徴の考え方

画像の特徴を考えるうえで、「特徴」というのが何によって表されるのかを考える必要がある。

まずは思考実験から始めてみよう。真っ白な画像と市松模様の画像を考えた時に、どちらがより「特徴的」な画像と言えるだろうか？

おそらく多くの読者は市松模様の画像を選ぶだろう。なぜなら、真っ白な画像はどの画素を見ても同じ値であって、特に画像上に変化がない一方で、市松模様の画像は画素の値が一定の間隔で変化しているからである。

この思考実験から分かることは、画像の特徴の一つとして「輝度が変化すること」を使えそうだということである。輝度が変化する箇所というのは、いわゆる エッジ が存在する場所である。

実際、エッジというのはものの構造を理解するうえで大変有用であり、その抽出法として、すでにSobelフィルタやCanny法などの手法について紹介した。また、Histogram of Oriented Gradients [18] (HOG)という特徴量では、ある範囲に存在するエッジの向きをヒストグラム化することで画像の特徴を表す。

ところで、エッジが存在する場所というのは、エッジを横切る方向に画像の輝度が変化する箇所である。しかしエッジに沿う方向には輝度は変化していない。

では、2本のエッジが交差する箇所ではどうだろうか？エッジの交差する箇所のことを コーナー と呼ぶが、コーナーの存在する箇所は2つのエッジが交差するため、より多くの方向に沿った輝度変化が生じていると考えられる。したがって、一つの見方としてはコーナーはエッジよりも、より「特徴的」な場所であると言える。

このような考えに基づき、画像の特徴抽出においては、より特徴的な点としてコーナーの箇所を 特徴点 (keypoint) として抽出し、その周囲の輝度変化等を 特徴量 (feature descriptor) として記述する。

6.2. コーナー検出

特徴点としてのコーナーを検出するにあたり、前節で紹介した Hough変換 などを用いて、直線同士が硬化する箇所を検出しても良いのだが、そもそも各画素がコーナーなのか、エッジなのかというのは、やや曖昧で判断が難しい面もある。

従って、各画素が平坦なのか、エッジなのか、コーナーなのかを、画素ごとに定量的に表す方法があれば有用だろう。このような画素のコーナーらしさのことを コーナネス (cornerness) と呼ぶ。

6.2.1. Moravecのコーナー検出

このようなコーナー点の検出方法についていち早く提案したのがHans Moravecで、彼の博士論文の中で、とある画素がコーナーであるかを検出するために画像パッチの類似度を用いる方法を提案した [19]。

彼の手法は、とある画素 \((x, y)\) を中心とする \(N \times N\) の画像パッチをを考え、そのパッチをとある方向 \((u, v)\) だけずらした位置にある画像パッチとの類似度 \(E(u, v)\) を計算する。この\(E(u, v)\) をいくつかの \((u, v)\) (通常は8近傍の画素)に対して計算し、その最小値をコーナーらしさ \(C_{\text{Moravec}}\) として用いた。

これを数式で表すと次のようになる。

\[ C_{\text{Moravec}}(x, y) = \min_{u, v} \sum_{\delta x, \delta y} ( I(x + \delta x, y + \delta y) - I(x + u + \delta x, y + v + \delta y))^2 \]

Moravecの手法について、もう少し深く見てみよう。画像が平坦で輝度の変化が少ない場所では、画像パッチをどの方向にずらしても輝度の変化は小さい。したがって、 \(C_{\text{Moravec}}\) の値は小さくなる。また、画素がエッジの上にある場合には、パッチをエッジに沿ってずらした場合に輝度の変化が小さくなる。故に、\(C_{\text{Moravec}}\) の値も小さくなる。

ところが、画素がコーナーの上にある場合には、パッチをどの方向にずらしても、輝度の変化が生じる。従って、\(C_{\text{Moravec}}\) の値は大きくなる。この考えに従い、Moravecの手法では、 \(C_{\text{Moravec}}\) の値が大きい画素をコーナーとして検出する。

Moravecは、2枚の画像のペアから各画像の深度を推定する ステレオ・マッチング にこの手法を応用し、その有用性を実証した (ステレオについては ステレオビジョン で扱う)。

しかし、Moravecの手法は、画像パッチの類似度を計算するために計算コストが大きく、また、パッチの大きさや、画像をずらす方向の取り方に結果が強く依存するため、コーナー検出の頑健性が低いという課題があった。

6.2.2. Harrisのコーナー検出

Moravecのコーナー検出の課題を解決したアルゴリズムに Harrisのコーナー検出 [20] がある。

Harrisのコーナー検出では画像パッチの類似度を計算する代わりに、局所的な輝度の変化から計算される 構造化テンソル (structural tensor)という2×2の行列によって、各画素の上のエッジらしさ (edgeness)やコーナーらしさ (cornerness)を表す。

構造化テンソル \(\mathbf{M}\) は各画素 \(I(x, y)\) における勾配 \(I_x(x, y)\) と \(I_y(x, y)\) を用いて、次のように定義される。

\[\begin{split} \mathbf{M} = \begin{pmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{pmatrix} \end{split}\]

勾配 \(I_x\) や \(I_y\) は、Sobelフィルタによって計算できるので、構造化テンソルもまた容易に計算できる。ただし、実用的には、離散的な勾配はノイズの影響を多く受けるので、構造化テンソルに対して Gaussianフィルタ等を適用して平滑化しておくことも多い。

構造化テンソルは テンソル という名前が示す通り、とある方向ベクトルに対して行列の二次形式を取ると、その方向に沿ったエッジの強度を与える。このような行列 (テンソル)の特徴として、行列の固有ベクトルが特徴的なエッジが現れる方向を、固有値がそのエッジの強さを表す。

実際、構造化テンソルの固有値を \(\lambda_1\), \(\lambda_2\) (\(\lambda_1 \geq \lambda_2\)) とすると、\(\lambda_1\) に対応する固有ベクトル \(\mathbf{u}_1\) の方向に沿って特徴的なエッジが現れると言える。また、言うまでもなく \(\lambda_2\) に対応する固有ベクトル \(\mathbf{u}_2\) は \(\mathbf{u}_1\) と直交していて、\(\mathbf{u}_2\)の方向は最もエッジ強度が弱い方向である。なお、構造化テンソルは正定値行列なので、 \(\lambda_1\) と \(\lambda_2\) はともに正の値である。

これを踏まえると、平坦、エッジ、コーナーの3つの状態は、構造化テンソルの固有値の値によって以下のように特徴づけられる。

• 平坦: 2つの固有値がともに小さい

• エッジ: 1つの固有値が大きく、もう1つの固有値が小さい

• コーナー: 2つの固有値がともに大きい

そこで、Harrisのコーナー検出においては、次のHarris応答と呼ばれる値 \(R\) を用いて、これらの状態を定量化する。

\[ R = \lambda_1 \lambda_2 - k (\lambda_1 + \lambda_2)^2 = \det(\mathbf{M}) - k \cdot \text{trace}(\mathbf{M})^2 \]

ここで、 \(k\) は応答の強さを調整するパラメータで、一般的には \(k = 0.04\) が用いられる。

さて、相加平均と相乗平均の関係を考えてみると、

\[ \lambda_1 \lambda_2 - \frac{1}{4} (\lambda_1 + \lambda_2)^2 \leq 0 \]

であると言える。この式を元に係数の \(1/4\) の部分を小さな値 \(k\) に置き換えたものがHarris応答で、Harris応答は正負の両方の値を取る。

Harris応答は両方の固有値が同程度の大きさであるときには正の値を取るが、その大きさは固有値の大きさに比例する。また、片方の固有値が大きく、もう片方の固有値が小さいときには、Harris応答は負の値を取る可能性がある。このことから、Harris応答が正の大きな値を取る場合に、その画素がコーナーであると考えられる。

しかし、Harris応答は、その計算過程において、構造化テンソルの要素に平滑化フィルタを適用しているため、ある画素でHarris応答が大きい時には、その周囲の画素でもHarris応答が大きくなる傾向がある。

そのため、Harrisのコーナー検出では、

• Harris応答が一定のしきい値以上である

• Harris応答が極大になっている

という2つの条件を満たす画素をコーナーとして検出する。

以下にHarrisのコーナー検出結果と各画素のHarris応答の大きさを示す。パラメータの設定の問題で、多少の取り逃がしがあるものの、概ねコーナーらしい箇所を正しく検出できていることが確認できる。

まとめ: コーナー検出

• 画像の中で特徴的な箇所を表す指標としてコーナーらしさ (cornerness) がよく用いられる

• コーナー検出アルゴリズムの元祖はMoravecの手法で、画像パッチの類似度を用いてコーナーらしさを定量化した

• Harrisのコーナー検出は、局所的な輝度変化を表す構造化テンソルを基にコーナーらしさを定量化した

6.3. Scale-Invariant Feature Transform (SIFT)

MoravecやHarrisの手法で検出されたコーナーは画像と画像の特徴的な点のマッチングにも応用できたが、画像が回転したり、スケールが変わったりすると、コーナーらしさも変化するという問題があった。

SIFTは、2001年にDavid Loweによって提案された手法で、このようなスケールや回転に対して不変な画像特徴量を目指して開発された [21] (2004年にはジャーナル論文としてまとめ直された [22])。

SIFTは画像のスケールに依存しない特徴点の検出のために Difference of Gaussian (DoG) を用いる。また、特徴量の記述においても、DoGを用いて得られた特徴点の「スケール」を考慮する。

6.3.1. DoGとLoG

SIFT以前の研究の中で、Gauss関数を二階微分して得られるLaplacian of Gaussian (LoG)を用いると、スケールを考慮した「blob」と呼ばれる特徴が得られることが分かっていた [23]。

Laplacian of Gaussianは数式としては以下のような関数である。

\[ \text{LoG}(x, y, \sigma) = \Delta G(x, y, \sigma) = -\frac{1}{\pi \sigma^4} \left( 1 - \frac{x^2 + y^2}{2 \sigma^2} \right) e^{-\frac{x^2 + y^2}{2\sigma^2}} \]

しかし、Gauss関数の画像の二階微分をフィルタカーネルとする場合、効率的なフィルタ処理が難しい。また、同様の処理はGaussianフィルタ後にLaplacianフィルタをかけることでも実現できるが、この場合はLaplacianフィルタにおける差分計算がノイズに敏感であるという問題がある。

このような背景から、SIFTではLoGの近似として、平滑化の強さ \(\sigma\) の異なる2つのGuassianフィルタ結果の差分を取る処理であるDoGを用いた。

Gaussianフィルタは計算コストが小さいので、異なる \(\sigma\) で画像をフィルタしておくこともそれほど問題にはならない。また、一度、フィルタ結果が計算できたなら、あとはフィルタ後の画像の単純な差分を取るだけでDoG画像が計算できる。

実際、Gauss関数を \(\sigma\) で微分してみると、

\[ \frac{\partial}{\partial \sigma} G(x, y, \sigma) = -\frac{1}{\pi\sigma^3} \left( 1 - \frac{x^2 + y^2}{2\sigma^2} \right) e^{-\frac{x^2 + y^2}{2\sigma^2}} \]

という式が得られるので、Gauss関数は次のような関係式を満たす。

\[ \sigma \Delta G(x, y, \sigma) = \frac{\partial}{\partial \sigma} G(x, y, \sigma) \]

このことから、差分を取る2つのGaussianフィルタのスケール \(\sigma_1\) と \(\sigma_2\) の差が比較的小さいという仮定において、DoGはLoGの近似として機能することが分かる。

6.3.2. DoGによる特徴点検出

DoGを用いた特徴量検出においては、 位置とスケールを含む座標 \((x, y, \sigma)\) におけるDoGの値を考える。このスケールを含むDoGの値を \(D(x, y, \sigma)\) とする時、 \(D(x, y, \sigma)\) が極大値を取るような \((x, y, \sigma)\) を特徴点として検出する。

この際、画像には異なる \(\sigma\) であらかじめGaussianフィルタをかけておき、隣あう \(\sigma\) の画像の差分を後で計算すれば、より効率的にDoG画像を取得できる。

今、入力画像を \(I(x, y)\) とし、最も弱いGaussianフィルタのスケールを \(\sigma_0\) としよう。また、Gaussianフィルタのスケールは平滑化のレベルを1つ上げるごとに \(k\) 倍になるとする。

この場合、\(n\) 番目のGaussianフィルタのスケールは次のように表される。

\[ L(x, y, k^n \sigma_0) = G(x, y, k^n \sigma_0) \ast I(x, y) \]

実際の画像に \(\sigma_0 = 1.6\), \(k = \sqrt{2}\) として、8レベルのGuassianフィルタを掛けると、次のような画像が得られる。

この3次元画像上で、極大になっている点を抽出し、そのときの \(\sigma\) の値を基に、特徴点を表示すると、次のような画像が得られる。

画像のスケールが変化する、ということは、画像の空間解像度が変化することと等価である。

従って、画像が小さくなるケースは、上記のDoG画像のレベルを1つ上げることとと同じ処理で、反対に画像が大きくなるケースは、DoG画像のレベルを1つ下げることと同じ処理である。

故に、DoGを利用して得られる特徴点は、スケール不変な特徴点になっている。

練習問題

実際にDoG特徴を計算するプログラムを作成し、画像のサイズが変化した場合にも、概ね同じ箇所に特徴点が得られることを確認せよ。

6.3.3. 特徴量の記述

SIFTでは、特徴量を計算するにあたり、特徴点周囲の輝度から、特徴の向きをはじめに計算する。

DoG画像上の点 \(D(x, y, \sigma)\) における特徴の向きは、対応するGaussianフィルタ画像 \(L(x, y, \sigma)\) を用いて計算される。具体的には、点 \((x, y)\) の周囲の画素 (例えば8×8の領域) において、各画素の勾配方向と、勾配強度を以下の差分式を用いて計算する。

\[\begin{split} \begin{align*} m(x, y) &= \sqrt{I_x^2(x, y) + I_y^2(x, y)} \\ \theta(x, y) &= \tan^{-1} \left( \frac{I_y(x, y)}{I_x(x, y)} \right) \end{align*} \end{split}\]

この結果を標本化された36方向に対してヒストグラム化する。各角度に投票される値は、勾配強度 \(m(x + \delta x, y + \delta y)\) に、中心点 \((x, y)\) からの距離を考慮したGauss関数による重みを掛けた値を用いる。

このようにして得られたヒストグラムに対して「最大値 × 0.8」以上の得票を得ているビンを特徴の向きとして保存する。この方法では、特徴の向きが複数検出される可能性があるが、その場合は、各向きに対して、以後の特徴量計算の処理を行う。

各特徴点について、その特徴の向きを可視化すると、次の画像のようになる。

このようにして各特徴点における特徴の向きが求まったら、特徴量を回転不変にするために、特徴点付近の画像をその向きに回転させる。その上で、スケールに沿った画像の周辺領域を16×16の領域に分割し、各領域における勾配強度と勾配の向きを計算する。

その後、各4×4領域ごとに勾配方向を8方向に量子化したヒストグラムを作成する。すると、8個のビンを持つヒストグラムが16個手に入るので、それらを連結して128次元のベクトルとする。

最後に、このベクトルを要素の最大値で割って正規化すると、輝度の変化に対しても頑健な特徴量としてSIFTの特徴記述子が得られる。

最後に得られた特徴を用いて、異なる画像同士で特徴点をマッチングしてみる。

以下の結果では、2つの微妙に向きの異なる画像に対してSIFT特徴量をマッチングしている。対応する特徴点同士は特徴量の距離が近い順に32個を直線で結んである。

この結果を見てみると、屋根や飾りの部分を中心に概ね正しい場所同士がマッチングできていることが分かる。

6.3.4. SIFTのその後

SIFTは画像の理解のための特徴量として広く普及したが、この技術は特許が取られていたために、商用の利用は難しく、それを回避するために類似の技術がいくつか開発された (SIFTの特許は2020年3月に失効した)。

また、SIFTはLoGの代わりにDoGを用いるなど、高速化のための工夫をいくつか取り入れていたものの、それでも計算量が小さいとは言えず、より高速に抽出可能な特徴量の開発が進められた。

このような流れの中で開発された手法の代表例の一つに SURF (Speeded Up Robust Features) [24] がある。SURFは、SIFTにおける特徴点の抽出で用いられていたDoGを、さらにBoxフィルタで近似して高速化した。さらに、特徴量の記述においても、特徴の向きの計算等にHaar-likeウェーブレット変換を用いることで、さらなる高速化を図った。

SURFはSIFTよりも高速な抽出が可能ではあったが、物体検出の精度はSIFTに劣ると言われている。また、SURFの技術にも特許が取られていたため、こちらも商用の利用には制限があった (SURFの特許は2029年4月に切れる)。

よりリアルタイム性を重視した画像特徴量の代表格には ORB (Oriented FAST and Rotated BRIEF) [25] がある。ORBは、FASTコーナー検出器を用いて特徴点の検出を高速かつ正確にし、さらにBRIEFと呼ばれる方法を応用して回転に頑健な特徴量を記述する。ORBはSIFTやSURFよりも計算量が小さく、リアルタイム処理が可能である。

また、ORBはSIFTやSURFと異なり、特許が取られていないため、商用利用にも制限がない。

しかし、今日においては、このような画像特徴量の利用は、深層学習によって概ね置き換えられたと言って良い。深層学習以前の画像特徴量は hand-crafted 特徴などと呼ばれている通り、どのような特徴が画像認識に有意なのかを人間が考えるものであったが、深層学習においては、画像の特徴を目的に応じて自動的に学習することが可能であることがその理由である。

まとめ: SIFT

• SIFTは、それ以前の特徴点検出やスケール不変特徴量の考え方を発展させた手法である

• 特徴点の検出には、Laplacian of Gaussian (LoG) の近似として、Difference of Gaussian (DoG) を用いることで高速化を図っている

• 各特徴点は「特徴の向き」の情報を持ち、特徴量の記述は、この向きに回転させた画像の輝度を用いて行われる

6.4. SIFTを用いた画像認識

ここまでで、SIFTのような特徴量を用いると、異なる画像間の対応する特徴点をマッチングできることが分かった。このような特徴点のマッチングは、物体追跡やステレオによる深度推定などに利用でき、それ自体が大変有用である。

さらに、SIFTのような画像特徴量は画像認識のために利用することもできる。そのためには、特徴点ごとに計算される特徴量から、画像全体の特徴を計算する必要がある。そこで、画像が「特徴点ごとの特徴量を単語とするような文章」であると考え、文書分類の分野で用いられる Bag of Words (BoW) の考え方を応用する。この手法のことを Bag of Visual Words (BoVW) や Bag of Features (BoF) と呼ぶ。

自然言語分野のBag of Wordsは辞書に登録された各単語が文書の中に何回登場するかをヒストグラムとして表わしたものであるが、SIFTのような画像特徴量は、自然言語の単語のようにピッタリと一致することはないので、そのままで単語のように扱うことはできない。

そこで、訓練画像全てから得られる特徴量をあらかじめクラスタ分割しておき、各特徴が属するクラスタを単語、クラスタの集まりを辞書と見なすことでBag of Wordsと同様の処理を可能にする。このような辞書のことを コードブック (codebook) と呼ぶ。

6.4.1. コードブックを利用した画像認識

コードブックを作成するための、最も基本的な方法は、訓練画像から得られる全ての特徴ベクトルをあらかじめクラスタリングする方法である。

クラスタリング とは、複数のデータをとある基準に従って、グループ分けする手法で、代表的な手法にはk平均法 (k-means clustering)などがある。

SIFTを例にとると、各特徴ベクトルは128次元であるので、128次元空間における点の集合に対して、距離が近いもの同士が同じクラスタに属するようにグループ分けをする。このようにしてグループ分けができたら、各クラスタの重心をそのクラスタ中心として、コードブック上に保存する (つまりクラスタ中心がコードになる)。

このようにしてコードブックが得られたら、各特徴ベクトルは、最も近いクラスタ中心に対応すコード、すなわち単語に置き換えられる。すると、とある画像1枚から得られた特徴量から、コードブックに登録されている単語の出現頻度をヒストグラム化することができる。これが、その画像全体の特徴として用いられる。

あとは、訓練画像内の全ての画像をコードブックを用いてヒストグラム化し、そのヒストグラムを用いて、画像分類等の機械学習器を学習すれば良い。

6.4.2. Gauss混合モデルを用いた改良

上記の例では、離散的に特徴量をクラスタ分割して、最も近いクラスタに対して、各特量を割り当てることでヒストグラムを作成した。ただし、複数のクラスタの境界に特徴量が存在している可能性もあるため、その情報を捨てて、一番近いクラスタに特徴量を割り当てるのは、あまり適切とは言えない。

そこで、クラスタリングを用いる代わりに、特徴量の分布をGauss混合モデルによって近似しておき、各ガウス分布への寄与 (小数で表わされる)の和を画像の特徴ベクトルとして与えることを考える。

Gauss混合モデル (Gaussian mixture model)とは、複数のGauss分布を和が1となるような正の重みによって重み付けたものであり、確率密度関数 \(P(\mathbf{x} | \boldsymbol\Theta)\), \(\mathbf{x} \in \mathbb{R}^D\) が以下の形で書ける。

(6.1)\[\begin{split} \begin{align*} P(\mathbf{x} | \boldsymbol\Theta) &= \sum_{k=1}^K \alpha_k \mathcal{N}(\mathbf{x} | \boldsymbol\mu_k, \boldsymbol\Sigma_k) \\ &= \sum_{k=1}^K \frac{\alpha_k}{(2 \pi)^{D / 2} \left| \boldsymbol\Sigma_k \right|^{1/2}} \exp \left( - \frac{1}{2}(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol\Sigma_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) \right) \end{align*} \end{split}\]

のように書ける。この式で \(\boldsymbol\Theta\) はパラメータの集合 \(\{ \alpha_k, \boldsymbol\mu_k, \boldsymbol\Sigma_k : k = 1, \ldots, K \}\) を表わし、 \(\alpha_k \in [0, 1)\), \(\boldsymbol\mu_k \in \mathbb{R}^D\), \(\boldsymbol\Sigma_k \in \mathbb{R}^{D \times D}\) は、それぞれ\(k\)番目のGauss分布に対する混合率、分布中心、共分散行列を表わす。

Gauss混合分布を用いると、とある特徴量 \(\mathbf{x}\) が \(k\)番目のGauss分布に対してどの程度の寄与率 \(\gamma_k(\mathbf{x})\) を持つかを以下のように計算できる。

(6.2)\[ \gamma_k(\mathbf{x}) = \frac{\alpha_k}{(2\pi)^{D/2} \left| \boldsymbol\Sigma_k \right|^{1/2}}\frac{\exp \left( - \frac{1}{2}(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol\Sigma_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) \right) }{P(\mathbf{x} | \boldsymbol\Theta)} \]

これにより、特徴量\(\mathbf{x}\)に対する、特徴表現として、

\[ \boldsymbol\gamma = (\gamma_1(\mathbf{x}), \ldots, \gamma_K(\mathbf{x}))^\top \]

が得られる (\(P(\mathbf{x} | \boldsymbol\Theta)\) での除算により、全要素の和が1になることに注意)。これを画像に含まれる全ての特徴量 \(\mathcal{X} = \{ \mathbf{x}_i : i = 1, \ldots, N \}\) に対して計算し、その平均値を画像の特徴ベクトルとする。

\[ \mathbf{f}(\mathcal{X}) = \frac{1}{N} \sum_{i=1}^N \boldsymbol\gamma(\mathbf{x}_i) \]

このように画像全体を特徴化すれば、各特徴量の寄与を全てのクラスタに対して考慮できるため、単純に最近傍のクラスタを選択する方法よりも良い特徴表現であると考えられる。

6.4.3. Fisherベクトルの利用

BoVWの改良としてFisherベクトルを利用する手法は2007年にPerronninらによって提案された[26]。Fisherベクトルは、統計学等でも用いられるFisher情報量に基づく特徴表現である。とある確率密度分布が\(\mathbf{x} \in \mathcal{X}\)のパラメータ (母数)\(\boldsymbol\Theta\)に関する事後分布として \(P(\mathbf{x} | \boldsymbol\Theta)\)のように与えられる場合を考える。

Fisherベクトルの導出

Fisher情報量を求めるに当たり、Fisherベクトル\(V(\mathbf{x}; \boldsymbol\Theta)\)を

\[ V(\mathbf{x}; \boldsymbol\Theta) = \nabla_{\boldsymbol\Theta} \log P(\mathbf{x} | \boldsymbol\Theta) \]

のように表わす。ここで、とあるパラメータ \(\theta\) に注目すると、次の等式が満たされる。

\[ \frac{\partial}{\partial \theta} \log P(\mathbf{x} | \boldsymbol\Theta) = \frac{1}{P(\mathbf{x} | \boldsymbol\Theta)} \frac{\partial}{\partial \theta} P(\mathbf{x} | \boldsymbol\Theta) \]

従って、Fisherベクトルのパラメータ \(\theta\) に関する要素の平均を考えると、次のようになる。

\[\begin{split} \begin{align*} \mathbb{E}_{\mathbf{x} \sim P(\mathbf{x} | \boldsymbol\Theta)} \left[ \frac{\partial}{\partial \theta} \log P(\mathbf{x} | \boldsymbol\Theta) \right] &= \int \frac{1}{P(\mathbf{x} | \boldsymbol\Theta)} \frac{\partial P(\mathbf{x} | \boldsymbol\Theta)}{\partial \theta} P(\mathbf{x} | \boldsymbol\Theta) \mathrm{d}\mathbf{x} \\ &= \frac{\partial}{\partial \theta} \int P(\mathbf{x} | \boldsymbol\Theta) \mathrm{d}\mathbf{x} \\ &= \frac{\partial}{\partial \theta} 1 = 0 \end{align*} \end{split}\]

以上から、Fisherベクトルの期待値はゼロベクトルに一致することが分かる。

Fisher情報量 \(\mathcal{I}(\mathcal{X})\) は上記のFisherベクトルの分散を用いて定義される。この際、Fisherベクトルの平均がゼロベクトルであることを用いると次の式で与えられる。

\[ \mathcal{I}(\mathcal{X}) = \mathbb{E}_{\mathbf{x} \sim P(\mathbf{x} | \boldsymbol\Theta)} \left[ (\nabla_{\boldsymbol\Theta} \log P(\mathbf{x} | \boldsymbol\Theta))^2 \right] \]

さて、ここで、\(P(\mathbf{x} | \boldsymbol\Theta)\) が(6.1) のGauss混合分布で与えられる場合を考える。以後、計算を簡単にするために、共分散行列は対角成分を \(\boldsymbol\sigma_k\) とする対角行列であるとする。すると、 \(P(\mathbf{x} | \boldsymbol\Theta)\)の\(\pi_k\), \(\mu_{k,d}\), \(\sigma_{k,d}\) に関する偏微分はそれぞれ以下のように書ける。

\[\begin{split} \begin{align} \frac{\partial P(\mathbf{x} | \boldsymbol\Theta)}{\partial \alpha_k} &= \frac{\exp \left( -\frac{1}{2}(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol\Sigma_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) \right)}{(2\pi)^{D/2} | \boldsymbol\Sigma_k |^{1/2}} = \mathcal{N}(\mathbf{x} | \boldsymbol\mu_{k}, \boldsymbol\sigma_k) \\ % \frac{\partial P(\mathbf{x} | \boldsymbol\Theta)}{\partial \mu_{k,d}} &= \alpha_k \left[ \frac{x_d - \mu_{k,d}}{\sigma_{k,d}^2} \right] \mathcal{N}(\mathbf{x} | \boldsymbol\mu_{k}, \boldsymbol\sigma_k) \\ % \frac{\partial P(\mathbf{x} | \boldsymbol\Theta)}{\partial \sigma_{k,d}} &= \alpha_k \left[ \frac{(x_d - \mu_{k,d})^2}{\sigma_{k,d}^3} - \frac{1}{\sigma_{k,d}} \right] \mathcal{N}(\mathbf{x} | \boldsymbol\mu_{k}, \boldsymbol\sigma_k) \end{align} \end{split}\]

したがって、(6.2)の寄与率\(\gamma_k\)を用いると、フィッシャーベクトルの各次元は以下のように書き直せる。

\[\begin{split} \begin{align} \frac{\partial \log P(\mathbf{x}| \boldsymbol\Theta)}{\partial \alpha_k} &= \frac{\gamma_k(\mathbf{x})}{\alpha_k} \\ % \frac{\partial \log P(\mathbf{x}| \boldsymbol\Theta)}{\partial \mu_{k,d}} &= \gamma_k(\mathbf{x}) \left[ \frac{x_d - \mu_{k,d}}{\sigma_{k,d}^2} \right] \\ % \frac{\partial \log P(\mathbf{x}| \boldsymbol\Theta)}{\partial \sigma_{k,d}} &= \gamma_k(\mathbf{x}) \left[ \frac{(x_d - \mu_{k,d})^2}{\sigma_{k,d}^3} - \frac{1}{\sigma_{k,d}} \right] \end{align} \end{split}\]

ここで、\(\mathcal{X} = \{ \mathbf{x}_1, \ldots, \mathbf{x}_N \}\)として、これらの同時分布\(P(\mathcal{X} | \boldsymbol\Theta)\)を考える。\(\mathbf{x}_i\)は互いに独立なサンプルであるので、

\[ P(\mathcal{X} | \boldsymbol\Theta) = \prod_{i=1}^N P(\mathbf{x}_i | \boldsymbol\Theta) \]

故に、この対数尤度 \(L(\mathcal{X} | \boldsymbol\Theta)\) は、

\[ L(\mathcal{X} | \boldsymbol\Theta) = \sum_{i=1}^N \log P(\mathbf{x}_i | \boldsymbol\Theta) \]

となることが分かる。この尤度関数のパラメータに関する偏微分からFisherベクトルを求めると、以下のようになる。

(6.3)\[ \frac{\partial L(\mathcal{X} | \boldsymbol\Theta)}{\partial \alpha_k} = \sum_{i=1}^N \left[ \frac{\gamma_k(\mathbf{x}_i)}{\alpha_k} - \frac{\gamma_1(\mathbf{x}_i)}{\alpha_1} \right] \quad \text{for} ~ k \geq 2 \]

(6.4)\[\begin{split} \frac{\partial L(\mathcal{X} | \boldsymbol\Theta)}{\partial \mu_{k,d}} = \sum_{i=1}^N \gamma_k(\mathbf{x}_i) \left[ \frac{x_d - \mu_{k,d}}{\sigma_{k,d}^2} \right] \\ \end{split}\]

(6.5)\[ \frac{\partial L(\mathcal{X} | \boldsymbol\Theta)}{\partial \sigma_{k,d}} = \sum_{i=1}^N \gamma_k(\mathbf{x}_i) \left[ \frac{(x_d - \mu_{k,d})^2}{\sigma_{k,d}^3} - \frac{1}{\sigma_{k,d}} \right] \]

ただし、(6.3)では、\(\sum_{k} \alpha_k = 1\)であることを用いて、自由度を一つ減らしてある。

ここまでの式から、Fisherベクトルの次元は\(\alpha_k\)に関する成分が\(K - 1\)次元、\(\mu_{k,d}\)に関する成分が\(DK\)次元、\(\sigma_{k,d}\)に関する成分が\(DK\)次元あるので、合計で\((2D + 1) K - 1\)次元となる。

Fisherベクトルの意味

上記の通り、Fisherベクトルはとあるデータ集合 \(\mathcal{X}\) を近似的に表わすGauss混合分布の対数尤度 \(\log P(\mathbf{x} | \boldsymbol\Theta)\) について、そのパラメータ \(\boldsymbol\Theta\) に関する勾配を要素に持つ。この意味を考えてみよう。

通常、確率密度関数をサンプル \(\mathcal{X}\) に沿うようにフィッティングする場合、その確率密度関数の対数尤度が最大になるようにパラメータを最適化する。この場合、パラメータ \(\theta\) に関して、対数尤度の勾配が大きい、ということは、そのパラメータの変化が尤度関数の変化に大きな影響を与えることを意味する。

言い換えれば、Fisherベクトルの各要素は、そのパラメータの変化が、どの程度尤度に影響を与えるかを示している。

この影響の分散を表わしているのがFisher情報量であり、パラメータ \(\theta\) に関しての分散が小さければ、そのパラメータに関しては、全体として尤度関数はほとんど影響を受けないことを意味する。一方で、分散が大きいときには、サンプル \(\mathbf{x} \in \mathcal{X}\) ごとに、パラメータ \(\theta\) の変化が尤度関数に与える影響にばらつきがあることを意味する。

従って、とある画像に対して計算されるFisherベクトルは、その画像(から得られた特徴量の集合)が、画像全体(のGauss混合分布)から見て、どのようなの意外性を持つか、と言い換えても良い。その画像が全体のGauss混合分布にとって十分尤もらしいものであれば、尤度の勾配から定まるFisherベクトルの要素は小さな値(特に意外ではない)を取り、そうでなければ、大きな値(意外である)を取る、という訳である。

そのため、Fisherベクトルを用いると、特定の画像にしか含まれないような特徴量を強調することができ、より分類に対して有利な画像の特徴表現が得られると期待できる。

定性的な議論にはなるが、以上が、Fisherベクトルが画像分類のタスクにおいて有用な理由である。

まとめ: SIFTを用いた画像認識

• SIFTをはじめとする特徴点ごとの特徴量は、コードブックを用いることで画像全体の特徴量を表現できる

• コードブックは、単純な特徴量のクラスタリングの他、Gauss混合分布によって記述される

• より良い画像の特徴表現を得るために、Gauss混合分布から計算されるFisherベクトルを利用する方法がある