位相限定相関法

2020/3/30公開

はじめに

位相限定相関法（POC:Phase-Only Correlation）は高速かつ高精度にテンプレートマッチングを行う手法である。大抵は画像の横ズレ・縦ズレを検出するために用いられる。さらにLog-Polar変換をかけた画像に対して適用することで回転ズレ・スケールの違いを検出する回転不変位相限定相関法（RIPOC:Rotation Invariant Phase-Only Correlation）が導かれる。

位相限定相関法

テンプレートマッチングはテンプレート画像が与えられた画像のどこに埋め込まれているかを探す問題のことである。

位相限定相関法はその名の通り、信号の位相の情報に限定して相関をとる手法で、フーリエ変換の性質をうまく利用することで「ある位置ズレを与えたときに2つの画像がどの程度一致しているか」の評価をO(1)O(1)で行うことを可能にしている。計算量のほとんどは画像をフーリエ変換・逆フーリエ変換する部分に使われる。

画像の離散フーリエ変換

R\mathbb{R}を実数全体の集合、C\mathbb{C}を複素数全体の集合とする。またM={0,1,2,⋯ ,M−1}\mathbb{M} = \{ 0, 1, 2, \cdots, M - 1 \}、,N={0,1,2,⋯ ,N−1},\mathbb{N}= \{ 0, 1, 2, \cdots, N - 1 \}とおく。

サイズM×NM \times Nの画像は位置x,yx,yに対して画素値f(x,y)f (x,y)が定まる信号とみなすことができる（より数学的な表記が好みの人のために書いておくとf ⁣:M×N→R ⁣:(x,y)↦f(x,y)f \colon \mathbb{M} \times \mathbb{N} \to \mathbb{R} \colon (x,y) \mapsto f(x,y)である）。

この信号を離散フーリエ変換（DFT:discrete Fourier transform）した信号をF(u,v)F(u,v)とおく（こちらはF ⁣:M×N→C ⁣:(u,v)↦F(u,v)F \colon \mathbb{M} \times \mathbb{N} \to \mathbb{C} \colon (u, v) \mapsto F(u,v)となる）。すなわち

F(u,v)=1MN∑x=0M−1∑y=0N−1f(x,y)exp⁡(−2πi(uxM+vyN))(1)F(u,v) = \frac{1}{\sqrt{MN}} \sum _ {x = 0} ^ {M - 1}\sum _ {y = 0} ^ {N - 1} f(x, y) \exp \left( - 2 \pi i \left( \frac{ux}{M} + \frac{vy}{N} \right) \right) \tag{1}

である。余計な記号が多いので

j=−2πiM,k=−2πiN(2)j = - \frac{2 \pi i}{M}, k= - \frac{2 \pi i}{N} \tag{2}

を用いて整えつつ、離散フーリエ変換の作用素F\mathcal{F}を導入しておくと

F(u,v)=(Ff)(u,v)=1MN∑x∈M∑y∈Nf(x,y)exp⁡(jux+kvy)(3)F(u,v) = (\mathcal{F} f)(u,v)= \frac{1}{\sqrt{MN}} \sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} f(x, y) \exp \left( jux + kvy \right) \tag{3}

である（作用素を導入しておくとF=FfF = \mathcal{F}fと書ける）。

逆離散フーリエ変換（IDFT:inverse discrete Fourier transform）は偏角が逆向きの変換で

f(x,y)=(F−1F)(x,y)=1MN∑u∈M∑v∈NF(u,v)exp⁡(−(jux+kvy))(4)f(x,y) = (\mathcal{F} ^ {-1} F)(x,y) = \frac{1}{\sqrt{MN}} \sum _ {u \in \mathbb{M}} \sum _ {v \in \mathbb{N}} F(u, v) \exp \left( - (jux + kvy) \right) \tag{4}

で定義できる（同様にf=F−1Ff = \mathcal{F} ^ {-1}Fである）。

位相限定画像

位相限定相関法では位相限定画像の内積を扱う。DFTした後の信号F(u,v)F(u, v)は複素数信号なので、振幅A(u,v)A(u,v)と位相exp⁡(iϕ(u,v))\exp (i \phi(u,v))の情報を持っており、それらの積の形に分解できる。すなわち

F(u,v)=A(u,v)exp⁡(iϕ(u,v))(5)F(u,v) = A(u,v) \exp (i \phi(u,v)) \tag{5}

である。位相限定相関法においては、この画像の位相成分exp⁡(iϕ(u,v))\exp (i \phi(u,v))を位相限定画像（phase-only image）と呼ぶ。

位相限定画像P(u,v)P(u,v)を求めるには、DFTした画像F(u,v)F(u,v)を、それ自身の振幅A(u,v)=∣F(u,v)∣A(u, v) = | F(u,v) | で割ればよい。すなわち直感的にはP(u,v)=exp⁡(iϕ(u,v))P (u,v) = \exp(i\phi(u,v))であり、

P(u,v)=F(u,v)∣F(u,v)∣(6)P(u,v) = \frac{F(u,v)}{| F(u,v) |} \tag{6}

である。ただし、ある(u,v)(u,v)において∣F(u,v)∣=0| F(u,v) | = 0のとき(6)(6)式の左辺は定義できない。より正確には(5)(5)式を満たしさえすればよいので振幅が 0 の周波数(u,v)(u,v)に対応する位相P(u,v)P(u,v)は任意である。実際の計算では分母に微小な正数ε\varepsilonを足すことで 0 割りによるエラーを防ぎつつ

P(u,v)=F(u,v)∣F(u,v)∣+ε=00+ε=0(7)P(u,v) = \frac{F(u,v)}{| F(u,v) | + \varepsilon} = \frac{0}{0 + \varepsilon} = 0 \tag{7}

としておく。

位相限定画像は、直感的には画像の形状に関する情報を保持している。実際に位相限定画像を逆フーリエ変換で空間領域の関数p=F−1Pp = \mathcal{F} ^ {-1} Pに戻してからその実部を表示すると、元の画像が浮かび上がる。一応丁寧に書いておくと、

p(x,y)=1MN∑x∈M∑y∈NP(u,v)exp⁡(−(jux+kvy))(8)p(x,y) = \frac{1}{\sqrt{MN}} \sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} P(u, v) \exp \left( - (jux + kvy) \right) \tag{8}

である。

一般に、画像は低周波成分の振幅が大きく、高周波成分の振幅が小さい。そして高周波成分には画像の輪郭の情報が含まれているから、周波数領域での画像F(u,v)F(u,v)をその振幅で割った位相限定画像への変換は元の画像の輪郭を強調する変換であるとも解釈できる。

位相限定相関法

位相限定相関法は理論的には空間領域における位相限定画像を用いたテンプレートマッチングである。与えられたサイズM′×N′M ^ \prime \times N ^ \primeの画像f(x′,y′)f(x ^ \prime,y ^ \prime)の中から、サイズM×N  (M≤M′,N≤N′)M \times N \,\,(M \leq M ^ \prime, N \leq N ^ \prime)のテンプレート画像g(x,y)g(x, y)が一致する場所を探す。

テンプレートが完全に一致する点(Δx,Δy)(\Delta x, \Delta y)では画像どうしの内積が最大になる。すなわちテンプレートマッチングの問題は

arg⁡max⁡(Δx,Δy)∑x∈M∑y∈Nf(x+Δx,y+Δy)g(x,y)(9)\underset{(\Delta x, \Delta y)}{\operatorname{arg}\operatorname{max}} \sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} f(x + \Delta x, y + \Delta y) g(x, y) \tag{9}

と書ける。位相限定相関法では代わりに位相限定画像どうしの内積が最大になる点を探す。つまり与えられた画像、テンプレート画像のそれぞれの位相限定画像をp(x′,y′),q(x,y)p(x ^ \prime, y ^ \prime), q(x,y)とおけば

arg⁡max⁡(Δx,Δy)∑x∈M∑y∈Np(x+Δx,y+Δy)q(x,y)(10)\underset{(\Delta x, \Delta y)}{\operatorname{arg}\operatorname{max}} \sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} p(x + \Delta x, y + \Delta y) q(x, y) \tag{10}

である。位相限定画像は元の画像の輪郭の情報が強調されているから、画像の輪郭がぴったり合う点を探せという問題に置き換わったとも解釈できる。

ところで(9)(9)式や(10)(10)式の

∑x∈M∑y∈Nf(x+Δx,y+Δy)g(x,y)(11)\sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} f(x + \Delta x, y + \Delta y) g(x, y) \tag{11}

の部分は、畳み込み積分によく似ている。実際これは2つの信号の間の相互相関関数（cross-correlation function）であり、片方の信号を反転して畳み込み積分をかけたものに一致する。信号ffとggの相互相関関数は

(f⋆g)(s,t)=∑x∈M∑y∈Nf(x+s,y+t)g(x,y)(12)(f \star g)(s, t) = \sum _ {x \in \mathbb{M}} \sum _ {y \in \mathbb{N}} f(x + s, y + t) g(x, y) \tag{12}

で定義され、s=Δx,t=Δys= \Delta x, t = \Delta yとおけば(11)(11)式に一致する。さて、相互相関関数には畳み込み積分とほぼ同様の性質が成り立ち、

F(f⋆g)(u,v)=(Ff)(u,v)(Fˉg)(u,v)=F(u,v)G(u,v)‾(13)\mathcal{F} (f \star g)(u,v) = (\mathcal{F} f)(u,v) (\bar{\mathcal{F}} g)(u,v) = F(u,v)\overline{G(u,v)} \tag{13}

である。両辺を逆フーリエ変換すれば

(f⋆g)(s,t)=(F−1(F(u,v)G(u,v)‾))(s,t)(14)(f \star g)(s, t) = \left( \mathcal{F} ^ {-1} (F (u,v) \overline{G(u,v)}) \right)(s,t) \tag{14}

であり、再びs=Δx,t=Δys= \Delta x, t = \Delta yとおいて、位相限定画像での議論に立ち戻れば、

(p⋆q)(Δx,Δy)=(F−1(P(u,v)Q(u,v)‾))(Δx,Δy)(15)(p \star q)(\Delta x, \Delta y) = \left( \mathcal{F} ^ {-1} (P (u,v) \overline{Q(u,v)}) \right)(\Delta x, \Delta y) \tag{15}

となって、テンプレート画像を(Δx,Δy)(\Delta x, \Delta y)だけずらしたときの相互相関関数が求まることになる。あとは相互相関が最大になる(Δx,Δy)(\Delta x, \Delta y)を見つければよい。

窓かけ

テンプレートマッチングでは画像の輪郭を合わせることが重要になるから、位相限定画像のうちでも重要なのは高周波成分である。したがってハイパスフィルタをかけてから相互相関をとることで精度の向上を狙うことがある。

回転不変位相限定相関法

回転不変位相限定相関法（RIPOC）は、位相限定相関法（POC）を用いて回転ズレとスケールの違いを検出する手法である。位相限定相関法は縦ズレと横ズレを検出できる手法なので、事前に Log-Polar 変換によって回転ズレを縦ズレに、スケールの違いを横ズレに変換してから POC を適用すれば回転ズレとスケールの違いを検出できることになる。

Log-Polar 変換

Log-Polar 変換は以下の式で定義される座標変換である。

{ρ=log⁡x2+y2θ=arctan⁡(y/x)\begin{cases} \rho = \log \sqrt{x ^ 2+ y ^ 2} \\ \theta = \arctan(y / x) \end{cases}

Log-Polar 変換の逆変換は以下の式で定義できる。

{x=eρcos⁡θy=eρsin⁡θ\begin{cases} x = e ^ \rho \cos \theta \\ y = e ^ \rho \sin \theta \end{cases}

解説に合わせて、画像としてプロットするときはρ\rhoを横軸、θ\thetaを縦軸としておこう。

元の画像の縦横がともにaa倍に引き伸ばされたとき、すなわちx′=1ax,y′=1ayx' = \frac{1}{a} x, y' = \frac{1}{a} yという変換を行ったとき、

ρ′=log⁡(x′)2+(y′)2=log⁡(1/a)2(x2+y2)=log⁡x2+y2−log⁡a2=ρ−log⁡a2\begin{aligned} \rho ' &= \log \sqrt{(x') ^ 2 + (y ') ^ 2} \\ &= \log \sqrt{ (1/a) ^ 2 (x ^ 2+ y ^ 2)} \\ &=\log \sqrt{x ^ 2+ y ^ 2} - \log a ^ 2 \\ &= \rho - \log a ^ 2 \end{aligned}

である。つまりもとの画像をaa倍に拡大すると、Log-Polar 変換後の画像は左側にlog⁡a2\log a ^ 2だけシフトする。

元の画像を反時計回りにθ0\theta _ 0だけ回転させる座標変換は

{x′=xcos⁡θ0−ysin⁡θ0y′=xsin⁡θ0+ycos⁡θ0\begin{cases} x' = x \cos \theta _ 0 - y \sin \theta _ 0 \\ y' = x \sin \theta _ 0 + y \cos \theta _ 0 \end{cases}

と表せる。このとき

ρ′=log⁡(x′)2+(y′)2=log⁡x2+y2=ρ\rho ' = \log \sqrt{(x') ^ 2 + (y') ^ 2} = \log \sqrt{x^ 2 + y ^ 2} = \rho

で変化はないが、

θ′=arctan⁡(y′/x′)=arctan⁡(xsin⁡θ0+ycos⁡θ0xcos⁡θ0−ysin⁡θ0)=arctan⁡(tan⁡θ0+yx1−yxtan⁡θ0)=arctan⁡(y/x)+arctan⁡(tan⁡θ0)=θ+θ0\begin{aligned} \theta ' &= \arctan(y' / x') \\ &= \arctan \left( \frac{x \sin \theta _ 0 + y \cos \theta _ 0}{x \cos \theta _ 0 - y \sin \theta _ 0} \right) \\ &= \arctan \left( \frac{\tan \theta _ 0 + \frac{y}{x}}{1 - \frac{y}{x} \tan \theta _ 0} \right) \\ &= \arctan(y/x) + \arctan(\tan \theta _ 0) \\ &= \theta + \theta _ 0 \end{aligned}

となる。したがって Log-Polar 変換後の画像は上にθ0\theta _ 0だけシフトする。

あとは POC によって得られた

{Δx=−log⁡a2Δy=θ0\begin{cases} \Delta x = - \log a ^ 2 \\ \Delta y = \theta _ 0 \end{cases}

の関係から

{a=exp⁡(−12Δx)θ0=Δy\begin{cases} a = \exp ( - \frac{1}{2} \Delta x) \\ \theta _ 0 = \Delta y \end{cases}

が求まる。

参考文献

Previous特異値分解Nextカーネル法（概要）

Last updated