Phase Congruency、位相一致は比較的に新しい概念です。一つの数学的評価法として、まだ30年ちょっとの年月しか経ってないようです。数千年の歴史を持つ有名な数学概念と比べれば、完全に新しいですね。こんなものですが、デジタル信号からの特徴抽出に有望かもしれません。とりあえず勉強してみました。
一本のデジタル信号は仮に\( x(n) \)とすれば、\( f(n) \)を\( x(n) \)のフーリエ数列の近似とします。つまり：
\(x(n) \doteq f(n) = \frac{\alpha_0}{2} + \sum_{k=1}^{N-1}\alpha_k cos(\frac{2\pi n k}{N}+\phi_k)\)
\( \alpha_k \ge 0\)
<式1>
この場合は、特定点の位相一致\( PC \)は下記の式で求められます。
\(
PC(n)=\max_{\phi \in [0,2\pi]}\frac{\sum_{k=1}^{N-1}{\alpha_k}cos(\frac{2\pi n k}{N}+\phi_k-\phi)}{\sum_{k=1}^{N-1}{\alpha_k}}
\)
<式2>
この式の計算はすごく大変ですが、加速する方法はあります。しかし、その前に\( f(n)\)を算出する必要があります。運良く、フーリエ変換からフーリエ数列を導出することはできます。
欲しがるフーリエ数列は<式1>のようなものですが、一気にそこまで行くというのはちょっと難しいので、第一歩としては下記の<式3>を得るようにします。
\(
f(n)=a_0+\sum_{k=1}^{N-1}(a_k cos(\frac{2\pi k n}{N}) + b_k sin(\frac{2\pi k n}{N}))
\)
<式3>
\( x(n) \)のフーリエ変換を\( X(n) \)とすれば、スケールを無視して、下記の各式で<式3>の各パラメータを求められます。
\(a_k=Re(X_k) , k=1..N-1\)
\(b_k=-Im(X_k) , k=1..N-1\)
\(a_0=Re(X_0)\)
次のステップは<式3>を<式2>に変換する。ここでも問題点が出てきます。\( a cos(t)+ b sin(t) \) は \( \gamma cos(t + \delta) \) になれば、\( \gamma \)と\( \delta \)はそれぞれなんでしょうか？方法はパラメータの比較です。
\( a cos(t)+ b sin(t) = \gamma cos(t + \delta) \)
\( a cos(t)+ b sin(t) = \gamma cos(t)cos(\delta) – \gamma sin(t)sin(\delta) \)
式の左右を比較して、下記のを得られます。
\( a = \gamma cos(\delta) \)
\( b = – \gamma sin(\delta) \)
両式を平方して足せば、\( \gamma = \pm \sqrt{a^2+b^2} \)になります。
両式を割り算すれば、\( \delta = \arctan(-\frac{b}{a}) = – \arctan(\frac{b}{a}) \)になります。
この結果で、式3から式1に戻すように組んでみますと、
\( \alpha_0 = 2 a_0 = 2 |{X_0}|\)
\( \alpha_k = \sqrt{a_k^2 + b_k^2} = \sqrt{(Re(X_k))^2 + (-Im(X_k))^2} = |X_k| \)
\( \phi_k = – \arctan( \frac{-Im(X_k)}{Re(X_k)} ) \) \(= \) \(\arctan( \frac{Im(X_k)}{Re(X_k)} ) \) \( = \angle {X_k} \)
これで、式1は\( X \)から完全に推定できます。
\( f(n) = {|X_0|} + \sum_{k=1}^{N-1}{|X_k|} cos(\frac{2\pi n k}{N}+\angle {X_k}) \)
<式4>
ひゅ～、できましたね。まだ道中場ですが、注意点があります。式4は\( X \)のスケールにって\( x \)とスケールが違うかもしれません。ただ、位相一致の計算には関係ない話です。
次の問題としては式2の\( \phi \)を知りたいです。
これは当該点の位相で分かります。
一点の位相を知るには、この点の同位相(in phase)と直角位相(quadrature phase)を知る必要があります。同位相は既に\( f(n) \)にあります。直角位相は\( f(n) \)のヒルベルト変換で得られます。しかし、\( f(n) \)は既に周波数領域にあるため、\( f(n) \)の各項に\( \frac{\pi}{2} \)を足して\( h(n) \)を得ました。そして、\( PC(n) \)に対して\( \phi = \arctan(\frac{h(n)}{f(n)}) \)になります。これで、\( PC(n) \)は計算できます。つまり：
\(
PC(n)=\frac{\sum_{k=1}^{N-1}{\alpha_k}cos(\frac{2\pi n k}{N}+\phi_k-\arctan(\frac{h(n)}{f(n)}))}{\sum_{k=1}^{N-1}{\alpha_k}}
\)
<式5>
式5は理論的に行けるが、実際に実行しますと、まだ結構時間が掛かります。実行時間を短縮するためでは、\( N \)を項を全部足さなくても一定の結果は得られます。結果を見てみましょう。
図1: 位相一致の精度検証、単位ステップ関数
図1から、低周波数の項目だけでもかなり正確な結果を獲得できます。なので、高周波を省くことで、かなりスピードアップできます。
次は、ノイズをいれた場合を見てましょう。同様に単位ステップですが、下記の結果を選られました。
図2: 位相一致の精度検証、単位ステップ関数、ノイズ入り
図2からは、安心できる結果を見れます。低周波数項では、ほぼノイズを無視してステップを検出できます。
位相一致の関数\( PC \)を画像に掛けてみました。横と縦両方向のみに適用しました。どっちも64個の周波数分量を計算に入れました。結果は下記の図です。
図3: 位相一致性の映像
図3からは、人物の輪郭が出ていますが、横と縦の線は多く見れます。もっと多くの方向から計算しないとうまく消せないと考えられています。この方法は理論上で行けますが、現実上では、スピードの問題で暫くはなかなか使えないかもしれません。とりあえず、いろいろともっと勉強してから、位相一致法をもう一度考える日が来ると思います。
参考：
http://www.myphysicslab.com/trig_identity1.html
https://en.wikipedia.org/wiki/Fourier_series
http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT7/node2.html
Peter Kovesi, Image Features From Phase Congruency, Robotics and Vision Reasearch Group, 1995
Peter Kovesi, Image Features From Phase Congruency, Videre: Journal of Computer Vision Research, Quarterly Journal, Summer 1999, Volume 1, Number 3, The MIT Press