画像フィルタ

3. 画像フィルタ#

フィルタ という言葉は、ラテン語で「ろ過する」という意味のfiltrumに由来している。

画像を始めとするデジタル信号の処理に「フィルタ」という言葉が使われるのは、ローパスフィルタ、ハイパス・フィルタのように特定の周波数成分を取り出す (=ろ過する)という処理によるものである。

現在の画像フィルタには様々な種類があり、一概に特定の成分だけを取り出すというような単純なものばかりではないが、画像から何らかの有意な情報を取り出すための処理、という点は共通している。

参考資料

  • 『ディジタル画像処理 改定第二版』 第5章-第7章

3.1. 基本的な画像フィルタ#

画像フィルタは広い意味では、とある画素のフィルタ後の値を決定するために、周囲の画素の値を用いる計算処理であり、その処理内容によって 線形フィルタ と 非線形フィルタ の2種類に分けられる。

線形フィルタは通常、画像とフィルタを定義するカーネル関数との畳み込みによって定義される。画像を \(I(x, y)\), カーネル関数を \(K(x, y)\) とすると、線形フィルタは連続的に次のように定義される。

\[ I'(x, y) = I \ast K = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} K(\xi, \eta) I(x+ \xi, y+\eta) \mathrm{d} \xi \mathrm{d} \eta \]

畳み込み演算は定義から分かる通り、入力 \(I\) に対して線形の演算であり、次のような性質を持つ。

\[ (a I_1 + b I_2) \ast K = (a I_1) \ast K + (b I_2) \ast K \]

多くの場合は、\(K(x, y)\) は有限の定義域を持っており、かつ画像は標本化されているので、次のような画素値の重み付き和で表現できる。

\[ I'(x, y) = \sum_{i=-m}^{m} \sum_{j=-n}^{n} K(i, j) I(x+i, y+j) \]

一方、非線形フィルタは、例えばカーネル関数の定義に周囲の画素の値を用いるフィルタなどが該当する。この場合、カーネル関数は画素ごとに異なるものとなり、一般には線形性を持たない。後述のエッジ保存型の平滑化フィルタであるバイラテラル・フィルタは非線形フィルタの一種である。

3.1.1. 平滑化フィルタ#

平滑化フィルタ とは、画像の詳細度を低減させるフィルタであり、ぼかしフィルタとも呼ばれる。最も代表的な平滑化フィルタはGaussianフィルタであり、カーネル関数にGauss関数、すなわち正規分布関数を用いる。

\[ K(x, y) = \frac{1}{2\pi \sigma^2} \exp\left(-\frac{x^2 + y^2}{2\sigma^2} \right) \]

この時、Gauss関数の分散 \(\sigma\) は、画像の詳細度をどのくらい低減させるかを決定するパラメータとなり、 \(\sigma\) が大きいほど、より詳細度の低い画像が出力される。

Gaussianフィルタ (\(\sigma = 9\))

Gaussianフィルタ (\(\sigma = 19\))

../_images/0286fb82f68d1b7e9b36184a36c872847d950d1c3fd4e55986c9912106ccde7c.png

../_images/254d3a365adb7863812af1a3794a5807f1641b103ca45c153fd5c108be5e3b9e.png

上記の連続的なカーネル関数の定義を離散化するには、カーネルのサイズを画素単位で指定する必要がある。また、カーネル関数は和が定義域内で1になるという性質を持つので、次のようにカーネルを離散化する。

\[ K_{i,j} = \frac{1}{\displaystyle\sum_{i=-m}^{m} \sum_{j=-n}^{n} K(i, j)} K(i, j) \]

なお、Gaussianフィルタのカーネルは有限の窓サイズに対して定義されることも多く、一例として \(3 \times 3\) のカーネルの場合には次のような定義が用いられる。

(3.1)#\[\begin{split} K = \frac{1}{16} \begin{bmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{bmatrix} \end{split}\]

数値計算においては expの計算は負荷の高い処理であるため、カーネルのサイズが小さい場合においては、このような近似的な定義を用いる方が効率が良いだろう (もっともカーネル関数を事前計算しておけばさほど問題はない)。

Gaussianフィルタの重要な性質にカーネル関数の変数分離がある。すなわち、二次元のGauss関数は、次のように一次元のGauss関数の積で表現できる。

\[ K(x, y) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{x^2}{2\sigma^2} \right) \cdot \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{y^2}{2\sigma^2} \right) \]

カーネル関数 \(K(x, y)\) を素直にfor文を二重にして畳み込む場合、その計算量はカーネル関数を評価する窓の大きさ \((m, n)\) に比例し、全体の計算量は画素数を \(N\) として \(O(Nmn)\) となる。

しかし、カーネル関数の変数分離を利用すると、縦方向の畳み込みに対して計算量が \(O(Nm)\)、横方向の畳み込みに対して計算量が \(O(Nn)\) となるので、全体の計算量は \(O(N(m+n))\) となり、特に \(m\) と \(n\) が大きい場合には、計算量を大幅に削減できる。

この変数分離の性質を考慮すると、(3.1) が [1, 2, 1] のような一次元のカーネル関数の積によって計算できることが分かる。

この [1, 2, 1] という一次元カーネルは、カーネルサイズが5であれば [1, 4, 6, 4, 1] のようになるが、これは紛れもなくPascalの三角形に現れる数字である。Pascalの三角形の各行は二項分布の係数を表すので \(p=0.5\)の場合を考えれば、それが一次元のGauss関数の離散的な近似となっている。

練習問題

Gaussianフィルタを実装し、変数分離を用いる場合と用いない場合とで画像のサイズに応じて計算量がどのように変化するかを調べよ。

練習問題

Gaussianフィルタのカーネル関数を二項分布の係数を用いて近似したものは、具体的にどの程度、正規分布と近くなるのだろうか。二項分布の平均と分散が試行回数 \(n\)と当たり確率\(p\)を用いて、それぞれ\(np\), \(np(1-p)\) であることを考慮して、両者の分布の差を調べよ。

3.1.2. エッジ抽出フィルタ#

画像処理において エッジ とは、画素の値に急激な変化がある部分を指す。このようなエッジは、最も単純には隣合う2つの画素の差分を取ることで計算することができる。

離散的なデータにおいて、差分を取る操作は 前進差分 (forward difference), 後退差分 (backward difference) ならびに 中心差分 (central difference) の3種類がある。画像を \(I(x, y)\) としx方向の差分を考える場合、それぞれの差分操作は次のように表せる。

\[\begin{split} \begin{align*} D_{\text{forward}}(x, y) &= \frac{1}{\Delta x} (I(x+\Delta x, y) - I(x, y)) \\ D_{\text{backward}}(x, y) &= \frac{1}{\Delta x} (I(x, y) - I(x-\Delta x, y)) \\ D_{\text{central}}(x, y) &= \frac{1}{2 \Delta x} (I(x+\Delta x, y) - I(x-\Delta x, y)) \end{align*} \end{split}\]

しかし、このような差分操作は画像のノイズの影響を受けやすく、ノイズにより発生する不要な輝度の変化がエッジとして検出されてしまうことがある。

3.1.3. Sobelフィルタ#

最も一般的なエッジ抽出フィルタに Sobelフィルタ がある。Sobelフィルタは平滑化フィルタによるノイズの低減と中心差分によるエッジ抽出を組み合わせたものである。

オリジナルのSobelフィルタは 3x3の離散的なカーネルにより定義されており、x方向、y方向の各方向のエッジを抽出するフィルタは次のように定義される。

\[\begin{split} K_x(x, y) = \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix}, \quad K_y(x, y) = \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix} \end{split}\]

このフィルタの定義に現れる [1, 2, 1] という係数はGaussianフィルタの時にも見られたものであり、上記のフィルタが平滑化フィルタと差分によるエッジ抽出を組み合わせたものであることが理解できる。

以下に、Sobelフィルタによるエッジ抽出の例を示す。本図では、正の差分値を赤色で、負の差分値を青色で表している。

x方向のエッジ

y方向のエッジ

../_images/2bd60fd91c08838dd8ff3c1a00b2794e16e8797ae2da50739f74683dccc0b3cc.png

../_images/363f707befacdb5104e1f0292285d4aa152a0df0eeac25c3d4f07bc8e8a95aed.png

Sobelフィルタを任意のカーネルサイズに拡張するために、5x5のカーネルを求めてみよう。なお、3x3以外のサイズのカーネルはいくつかの異なる定義が存在するが、ここではOpenCVに実装されている定義について説明する。

まずは、平滑化フィルタの成分として、5x5のGaussianーネル関数を用意する。

\[\begin{split} \begin{bmatrix} 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{bmatrix} \end{split}\]

次に、中心差分を取るための行列を用意するのだが、この際、中心画素の両側で差分が取られる2画素の距離に応じて、係数を調整する。この際、カーネルの端に現れる係数の絶対値は1であるものとし、中心画素に近づくに連れて等差数列的に0に近づくように係数を設定する。

\[\begin{split} D_x = \begin{bmatrix} -1 & -\frac{1}{2} & 0 & \frac{1}{2} & 1 \\ -1 & -\frac{1}{2} & 0 & \frac{1}{2} & 1 \\ -1 & -\frac{1}{2} & 0 & \frac{1}{2} & 1 \\ -1 & -\frac{1}{2} & 0 & \frac{1}{2} & 1 \\ -1 & -\frac{1}{2} & 0 & \frac{1}{2} & 1 \end{bmatrix}, \quad D_y = \begin{bmatrix} -1 & -1 & -1 & -1 & -1 \\ -\frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} \\ 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ 1 & 1 & 1 & 1 & 1 \end{bmatrix} \end{split}\]

最後に、これら二つの行列の要素積を取ると5x5のSobelフィルタとして以下の行列が得られる。

\[\begin{split} K_x = \begin{bmatrix} -1 & -2 & 0 & 2 & 1 \\ -4 & -8 & 0 & 8 & 4 \\ -6 & -12 & 0 & 12 & 6 \\ -4 & -8 & 0 & 8 & 4 \\ -1 & -2 & 0 & 2 & 1 \end{bmatrix} \quad K_y = \begin{bmatrix} -1 & -4 & -6 & -4 & -1 \\ -2 & -8 & -12 & -8 & -2 \\ 0 & 0 & 0 & 0 & 0 \\ 2 & 8 & 12 & 8 & 2 \\ 1 & 4 & 6 & 4 & 1 \end{bmatrix} \end{split}\]

二項分布の係数 \({}_n C_k\)の性質から、以上の方法で計算されるフィルタカーネルに現れる係数は任意のカーネルサイズ \(N\)について必ず整数になる。

3.1.4. Cannyフィルタ#

Sobelフィルタと並んで、よく用いられるエッジ検出フィルタに Cannyフィルタがある。Sobelフィルタが、「勾配がある場所」をエッジだと考える一方で、エッジを物体の輪郭線であるとみなすと、その線幅は1ピクセルであることが好ましい場合もある。

CannyフィルタはSobelフィルタの結果から、「物体の輪郭線としてのエッジ」をおよそ1ピクセル幅で検出するアルゴリズムである。Canny フィルタは主につの処理からなる。

  1. Sobel フィルタ (3x3 のフィルタカーネル)による勾配強度の計算

  2. 勾配が極大となる画素の検出

  3. ヒステリシス閾値処理

以下では、これら三つの処理について順に見ていこう。

Sobel フィルタによる勾配強度の計算#

まず、Sobel フィルタを用いて、x方向の勾配 \(\Delta x\) とy方向の勾配 \(\Delta y\) を計算する。具体的にはx方向とy方向の各方向に対してSobelフィルタを用いてエッジ検出を行い、その勾配のノルムを各画素のエッジ強度 \(E\) であるとする。

\[ E = \sqrt{(\Delta x)^2 + (\Delta y)^2} \]

なお、Cannyフィルタは、Sobelフィルタをかける前に、あらかじめ入力画像にGaussianフィルタをかけてノイズを低減する処理が入る。このような処理で得られたエッジ強度画像を以下に示す。

入力画像

エッジ強度

../_images/472bf3643be473b3c0cf8d8652d2ee735eeae6efc117e47b04ec0db75b839748.png

../_images/7a808494c37cadf63d95e154534aa372509d8b87ef8f5eddfff9c96f96274ffc.png

勾配が極大の画素の検出#

次に、各画素について、勾配方向を計算し、その勾配方向に沿った方向に中心画素の輝度値が極大になっている画素だけを取り出す。

まずは、先ほどの x方向勾配 \(\Delta x\) と y方向勾配 \(\Delta y\) から、余接関数の逆関数を用いて勾配方向 \(\theta\) を求める。数学的には、このような勾配方向 \(\theta\) は

\[ \theta = \tan^{-1}\left(\frac{\Delta y}{\Delta x}\right) \]

のように求められる。

しかし、計算機的には \(\Delta x\) が0の場合にゼロ除算によって正しく計算ができない場合があるので、単純な余接関数の逆関数を計算する math.atanやnumpy.arctanの代わりにxとyの値を別々に引数として取る math.atan2 や numpy.arctan2 などの関数を用いるほうが好ましい。

なお、math.atan などは \(0\)から\(\pi\)の範囲で値を返すのに対して、math.atan2 などは \(-\pi\)から\(\pi\)の範囲で値を返すので注意が必要である。今回は「向きを無視した方向」が分かれば良いので、arctanの戻り値が負の場合には\(\pi\)を足し算しておく。

最後に得られた角度 \(\theta\)を量子化して、0°、45°、90°、135° いずれかに割り当てる。

角度 \(\theta\) が0°, 45°, 90°, 135° のいずれかに割り当てられたら、その方向と画素の周囲3x3の画素 (このような隣接画素を8隣接画素という)を用いて、注目画素がエッジの中心にあるかどうかを調べる。

具体的には、0°であれば左右の画素に比べて中心の画素のエッジ強度が大きいかどうかを調べ、45°であれば左上と右下の画素に比べて中心の画素のエッジ強度が大きいかどうかを調べるといった具合になる (画像の画素の座標は左上が原点であることに注意)。

このようにして、エッジ強度が極大の画素だけを取り出した結果を以下に示す。この結果から、エッジ強度極大の画素だけを取り出すと、かなりエッジの線幅が細くなったことが確認できる。

エッジ強度

エッジ強度極大画素

../_images/7a808494c37cadf63d95e154534aa372509d8b87ef8f5eddfff9c96f96274ffc.png

../_images/78ba5c638b0690bdb5c6f5b577a7964ae55fdbc4637df539535b9abdfdddb810.png

ヒステリシス閾値処理#

ヒステリシス というのは、元々、とある状態から変化を加えていったときに、変化の過程と、元に戻る過程で異なる状態経路を取るような現象のことを指す。Canny フィルタで用いられる「ヒステリシス閾値処理」は、これに着想を得て、二段階の閾値処理をかける手法である。

ヒステリシス閾値処理には二種類の閾値 \(\tau_1\)と\(\tau_2\)を使用する (ただし\(\tau_1 < \tau_2\))。もし上記の処理で取り出したエッジ強度極大画素の勾配の値が\(\tau_2\)を超えていたら、その画素は確定で輪郭線であると見なす。また、反対に極大勾配の値が\(\tau_1\)を下回っていたら、その画素は確定で輪郭線ではないと見なす。

残ったエッジ強度極大画素の勾配 \([ \tau_1, \tau_2 ]\) の区間に入っているが、そのような画素については8隣接の画素をみて、少なくとも1つの画素が\(\tau_2\)以上の勾配値を持つエッジ強度極大画素であるなら、その画素も輪郭線と見なす。

以下に、画素値が\([0, 255]\)の256階調で表されている場合に \(\tau_1 = 50\), \(\tau_2 = 100\) としてヒステリシス閾値処理を実行した結果を示す。

エッジ強度極大画素

Cannyフィルタ

../_images/78ba5c638b0690bdb5c6f5b577a7964ae55fdbc4637df539535b9abdfdddb810.png

../_images/b89220b287495e2d01be0a81995b404f436db4094482dd8599213384064ecb51.png

以上がCannyフィルタの処理の概要で、結果の画像に示す通り、線幅がおよそ1ピクセルの輪郭線が抽出できていることが分かる。

練習問題

SobelフィルタとCannyフィルタの各フィルタは、それぞれどのような場合に用いるのが有用だろうか。得られるエッジの特徴を踏まえて考察せよ。

3.1.5. エッジ保存型の平滑化フィルタ#

前述したGaussianフィルタなどの平滑化フィルタは画像全体を同じ強さでぼかすため、出力される画像からは元の画像に見られるような物体の詳細は一様に失われてしまう。しかし、単に画像からノイズを除去したい場合などは、画像としては平滑化しつつもエッジのような画像中の物体の特徴を表すような情報は保持したい。

このような場合に用いられるのが エッジ保存型の平滑化フィルタ である。代表的なものには メディアンフィルタ や バイラテラルフィルタ (bilateral filter) [10] がある。

メディアンフィルタ#

メディアンフィルタは画像のRチャネル, Gチャネル, Bチャネルそれぞれに対して個別に実行される。

各チャネルに対するフィルタ処理では、各画素の周囲を含む正方形領域内の画素値をソートし、その中央値をフィルタ後の画素値として採用する。こうすることでノイズ含む画素のような周囲と比べて極端に明るい画素や極端に暗い画素の影響を低減することができる。

以下の図に示す通り、メディアンフィルタはノイズがかなり強い画像に対しても元の画像を上手く再現できていることが分かる。

入力画像

メディアンフィルタ

../_images/f2f085b630c0e68e08a421822564818313580518d51ebcf40fa7cf966dcc2f99.png

../_images/59dcd5032cedd630bee3fa83e4508b99d6228a299377456cb65e247adff45f09.png

一方、メディアンフィルタが上手く動作するのは、ノイズが空間的に離れた箇所に生じている場合に限られる。このようなノイズをごま塩ノイズ (salt-and-pepper noise) という。

そのため、画像全体に一様に発生するノイズ (このようなノイズをGaussノイズあるいはホワイトノイズと呼ぶ) 画像の輝度に比例するような強さのノイズ (このようなノイズをPoissonノイズと呼ぶ)に対してはメディアンフィルタはそれほど効果的ではない。

バイラテラルフィルタ#

バイラテラルフィルタ (bilateral filter)は Tomasi と Manduchi によって 1998 年に提案されたフィルタ[10]で、Gaussian フィルタが元画像の信号までぼかしてしまう問題を緩和したものである。

具体的には、色が近い画素同士の間でだけぼかしがかかるように Gauss 関数に色を考慮する項を追加している。バイラテラルフィルタのフィルタカーネルは以下のように定義される。

\[ f(I, x, y, \xi, \eta) = \frac{1}{K} \exp\left(- \frac{(\xi^2 + \eta^2)}{2 \sigma_s^2} \right) \exp\left( -\frac{\| I(x, y) - I(x + \xi, y + \eta)\|^2}{2 \sigma_r^2} \right) \]

ただし、\(\sigma_s\)は画素の位置情報に関するぼかしの強さを調整するパラメータ (s は spatial の s)、\(\sigma_r\)は画像の輝度に関するぼかしの強さを調整するパラメータ (r は range の r)である。このとき、フィルタカーネルは画像の輝度値を参照するため、原点からのずれ \((\xi, \eta)\)以外にフィルタ対象の画像 \(I\)と画素の座標 \((x, y)\)も取っている。

したがって、バイラテラルフィルタは畳み込み積分に見られたような線形性を持たない非線形フィルタの一種である。以下に、Poissonノイズを付加したノイズ画像に対するバイラテラルフィルタの適用結果を示す。

入力画像

バイラテラルフィルタ

../_images/95886453d841e09b5f1e122c649c569451615dc23d00ba99d3a7343bd64e8a9b.png

../_images/a679486f70732fd092ebd22519d78e9a87e63dca61df3aa8e9d900e5397aeed8.png

この結果から、バイラテラルフィルタは画像全体に一様に生じるPoissonノイズ (ならびにホワイトノイズ)に対して効果的であることが分かる。

通常、カメラを通して撮影される写真はISO感度を高くするとPoissonノイズ (ショットノイズとも呼ばれる)が、画素ごとの電気信号の読み出し時にホワイトノイズが生じる。そのため、バイラテラルフィルタはガメラ画像のノイズ除去によって広く用いられている。

しかし、バイラテラルフィルタはその非線形性から比較的計算量の大きなフィルタである。特にカーネルサイズを大きく取る場合にはその計算時間が増大する。補足: Pythonでの画像フィルタの実装 にPythonによるバイラテラルフィルタの素直な実装例を示したが、 450x300のようなかなり小さな画像に対して、(Pythonの計算速度が遅いことを加味しても) かなりの時間がかかることが分かる。

そのため、Tomasiらのバイラテラルフィルタを高速化する手法は非常に広く研究されている。典型的なアイディアとしては、画素の空間的な位置によるフィルタ部分と画素の輝度に対するフィルタ部分を何らかの方法で、Guassianフィルタのような線形フィルタで近似するというものである。

例えば、OpenCVの cv2.bilateralFilter では、Bilateral Gridと呼ばれる画像を輝度方向に拡張したような三次元ボリューム画像を用意し、Gaussianフィルタの分離可能性を利用して、各方向沿って順に平滑化処理を行うことで大幅に処理時間を短縮している。処理の詳細については原著論文を参照されたい [11]。

3.2. Fourier変換の基本#

画像のフィルタ処理を理解するうえで重要な数学的な道具に Fourier変換 がある。Fourier変換は、数学者のJoseph Fourierが熱方程式の解法を考えている中で発見したものである。Fourier変換そのものの説明に入る前に、その思考の過程を簡単に説明しておこう。

3.2.1. 熱拡散方程式の解法#

熱拡散方程式は時刻 \(t\) における位置 \(x\)での熱量を表す関数 \(u(x, t)\) に関する以下の偏微分方程式を指す。

(3.2)#\[ \frac{\partial u(x, t)}{\partial t} = \alpha \frac{\partial^2 u(x, t)}{\partial x^2} \]

この偏微分方程式を解く際には、\(u(x, t)\)が次のような形式でかけると仮定する。

\[ u(x, t) = X(x) T(t) \]

この表現を用いると、\(x\), \(t\) に関する偏微分は \(X(x)\) と \(T(t)\) に対する全微分で書き換えられるので、 (3.2) は次のように書き換えられる。

\[\begin{split} \begin{align*} & X(x) T'(t) = \alpha X''(x) T(t) \\ \therefore & \frac{X(x)}{X''(x)} = \frac{T'(t)}{\alpha T(t)} = - \lambda \end{align*} \end{split}\]

本式は第1項, 第2項が \(x\), \(t\)の変化に依らないことから定数であり、これを \(-\lambda\) とする。すると、次の二式が得られる。

\[\begin{split} \begin{cases} X(x) = -\lambda X''(x) \\ T(t) = -\alpha \lambda T'(t) \end{cases} \end{split}\]

この式から \(X(x)\), \(T(t)\) は次のように書ける。

\[\begin{split} \begin{align*} X(x) &= A \sin(- \sqrt{\lambda} x) + B \cos(-\sqrt{\lambda} x) \\ T(t) &= C \exp(-\alpha \lambda t) \end{align*} \end{split}\]

ここで、熱分布\(u(x, t)\)が\([0, L]\)の範囲で定義されていて、境界条件として \(u(0, t) = u(L, t) = 0\) が与えられているとする。さらに、\(u(x, 0)\)が単純な正弦波として

\[ u(x, 0) = \sin\left(\frac{n \pi x}{L}\right) \]

と与えられている (\(n\) は適当な自然数)とする。すると、\(X(x)\), \(T(t)\) は次のように書ける。

\[\begin{split} \begin{align*} X(x) &= \sin\left(\frac{n \pi x}{L}\right) \\ T(t) &= C \exp\left(-\frac{\alpha n^2 \pi^2}{L^2} t\right) \end{align*} \end{split}\]

この式から \(u(x, t)\)が初期条件として、正弦波で表される場合には、熱拡散方程式の解が得られることが分かる。しかし、これをより一般的な熱分布を初期条件とする場合、すなわち \(u(x, 0) = f(x)\) のような場合にどのように解けばよいか、というのが当時の関心事であった。

3.2.2. Fourier級数展開#

前述の熱方程式の例では、 \(f(x)\) は \([0, L]\) の範囲で定義されており、境界条件 \(u(0, t) = u(L, t) = 0\) から周期関数であるとみなせる。このような任意の周期関数が異なる正弦波、余弦波の重ね合わせで表現できる、というのがFourierのアイディアである。

以下、説明を簡単にするために\(f(x)\)が、\([-L, L]\) の範囲で定義された周期 \(2L\) の周期関数であるとする。この時、\(f(x)\)は次のFourier級数によって表せる。

\[ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos\left(\frac{n \pi x}{L}\right) + b_n \sin\left(\frac{n \pi x}{L}\right) \right) \]

ここで、三角関数を1周期分の定義域で積分した際の性質として以下の等式を用いる。

\[\begin{split} \begin{align*} \int_{-L}^{L} \sin\left(\frac{n \pi x}{L}\right) \mathrm{d}x &= 0 \\ \int_{-L}^{L} \cos\left(\frac{n \pi x}{L}\right) \mathrm{d}x &= 0 \end{align*} \end{split}\]

この性質を用いると、両辺を\([0, L]\)の範囲で積分することで、係数 \(a_0\) は次のように書ける。

\[ a_0 = \int_{-L}^L f(x) \mathrm{d}x \]

さらに、三角関数の積を1周期分の定義域で積分した際の次の等式を用いる。

\[\begin{split} \begin{align*} \int_{-L}^{L} \sin\left(\frac{n \pi x}{L}\right) \sin\left(\frac{m \pi x}{L}\right) \mathrm{d}x &= L \delta_{nm} \\ \int_{-L}^{L} \cos\left(\frac{n \pi x}{L}\right) \cos\left(\frac{m \pi x}{L}\right) \mathrm{d}x &= L \delta_{nm} \\ \int_{-L}^{L} \sin\left(\frac{n \pi x}{L}\right) \cos\left(\frac{m \pi x}{L}\right) \mathrm{d}x &= 0 \end{align*} \end{split}\]

なお \(\delta_{nm}\) は \(n = m\) のときに1を、それ以外のときに0を取るものとする。これらの式を用いると各 \(a_n\), \(b_n\) は次のように書ける。

\[\begin{split} \begin{align*} a_n &= \frac{1}{L} \int_{-L}^{L} f(x) \cos\left(\frac{n \pi x}{L}\right) \mathrm{d}x \\ b_n &= \frac{1}{L} \int_{-L}^{L} f(x) \sin\left(\frac{n \pi x}{L}\right) \mathrm{d}x \end{align*} \end{split}\]

また、 \(e^{i \theta} = \cos(\theta) + i \sin(\theta)\) というオイラーの公式を用いると、Fourier級数展開は次のように書き直すこともできる。

(3.3)#\[ f(x) = \sum_{n=-\infty}^{\infty} c_n e^{i \frac{n \pi x}{L}} \]

この際、係数 \(c_n\) は次のように定義される。

(3.4)#\[ c_n = \frac{1}{2L} \int_{-L}^{L} f(x) e^{-i \frac{n \pi x}{L}} \mathrm{d}x \]

このような複素数を用いたFourier級数展開を 複素Fourier級数展開 と呼ぶ。

以上の説明は、Fourier級数展開の式の必要条件を調べたに過ぎず、実際に右辺の無限級数が\(f(x)\)に収束するかどうかを示してはいない。Fourier級数展開の収束性に関する議論は少々複雑な数学的な議論を必要とするため、より専門的な書籍 (「関数解析」に関する書籍を見るとよいだろう)に譲ることとする。

3.2.3. Fourier変換#

ここでFourier級数展開に関して一つの疑問が生じる。Fourier級数展開できるのは一定の周期を持つ周期関数のみであるが、これを周期を持たないような一般の関数にも拡張できるのだろうか。

周期を持たない、ということは、言い換えれば周期が \((-\infty, \infty)\) である、とみなすことができる。そこで、周期 \(L\) を無限大に近づけるとFourier級数展開がどのようになるかを考えてみる。

まず複数Fourier級数展開の係数 \(c_n\) の式 (3.4) を複素Fourier級数展開の式 (3.3) に代入してみよう。すると、次のような式が得られる。

\[ f(x) = \sum_{n=-\infty}^{\infty} \left( \frac{1}{2L} \int_{-L}^{L} f(\xi) e^{-i \frac{n \pi \xi}{L}} \mathrm{d} \xi \right) e^{i \frac{n \pi x}{L}} \]

ここで \(L \rightarrow \infty\) の極限を考えるために、 \(\omega_n = \frac{n \pi}{L}\), \(\Delta \omega = \frac{\pi}{L}\) と置き換えると、次の式が得られる。

\[ f(x) = \frac{1}{2\pi} \sum_{n=-\infty}^{\infty} \Delta \omega \int_{-L}^{L} f(\xi) e^{-i \omega_n \xi} e^{i \omega_n x} \mathrm{d} \xi \]

この式はRiemann積分の定義式そのものなので、\(L \rightarrow \infty\) の極限を取ると和の部分を積分によって置き換えることができる。

\[ f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\xi) e^{-i \omega \xi} e^{i \omega x} \mathrm{d} \omega \mathrm{d} \xi \]

このような操作により関数の周期が \((-\infty, \infty)\) である、すなわち特定の周期を持たない一般的な関数に対しては、Fourier級数展開において数列として表されていた係数 \(c_n\) が \(\omega\) に対する関数として、

\[ \mathcal{F} [f] = F(\omega) = \frac{1}{2\pi} \int_{-\infty}^{\infty} f(x) e^{-i \omega x} \mathrm{d} x \]

という形で表せることが分かる。この操作を Fourier変換 と呼ぶ。反対にFourier変換によって得られる周波数 \(\omega\) の関数を空間的な位置 \(x\) に関する関数に戻す操作

(3.5)#\[ \mathcal{F}^{-1} [F] = f(x) = \int_{-\infty}^{\infty} F(\omega) e^{i \omega x} \mathrm{d} \omega \]

を 逆Fourier変換 と呼ぶ。

3.2.4. Fourier変換と畳み込み積分#

Fourier変換はいくつもの興味深い性質を持つが、本項では画像信号処理において重要な畳み込み積分とFourier変換の関係について説明する。

一般に、二つの連続関数 \(f(x)\), \(g(x)\) の畳み込み積分は次のように定義される。

\[ (f \ast g)(x) = \int_{-\infty}^{\infty} f(\xi) g(x - \xi) \mathrm{d} \xi \]

したがって、このFourier変換は次のように書ける。

\[ \mathcal{F}[f \ast g] = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\xi) g(x - \xi) e^{-i \omega x} \mathrm{d}\xi \mathrm{d} x \]

この式は\(y = x - \xi\) と変数変換すると、次のように書き直せる。

\[ \mathcal{F}[f \ast g] = \int_{-\infty}^{\infty} f(\xi) \left( \int_{-\infty}^{\infty} g(y) e^{-i \omega (y + \xi)} \mathrm{d} y \right) \mathrm{d} \xi = \int_{-\infty}^{\infty} f(\xi) e^{-i \omega \xi} \mathrm{d} \xi \int_{-\infty}^{\infty} g(y) e^{-i \omega y} \mathrm{d} y \]

以上より、畳み込み積分のFourier変換は畳み込まれる二つの関数のFourier変換の積として、

\[ \mathcal{F}[f \ast g] = \mathcal{F}[f] \cdot \mathcal{F}[g] \]

のように書けることが分かる。

この性質は計算機上で画像フィルタを行う上で非常に重要である。これは、Fourier変換を離散的な信号に対して計算する離散Fourier変換 (discrete Fourier transform, DFT)の計算量が、信号の長さ \(N\) に対して \(O(N \log N)\) で実現できることに起因する。このような高速な離散Fourier変換を実現するアルゴリズムを 高速Fourier変換 (fast Fourier transform)と呼び、しばしばFFTと略される。

仮に画像と離散的なカーネル関数が、どちらも \(N \times N\) の大きさで表されているとすると、これらをFourier変換するのに必要な計算量は \(O(N^2 \log N)\) である。周波数空間でこの結果の積を取る計算量は \(O(N^2)\), 逆Fourier変換にかかる計算量は \(O(N^2 \log N)\) なので、全体の計算量は \(O(N^2 \log N)\) に抑えられる。

同様の計算を位置空間での畳み込みで実現する場合、各画素におけるフィルタ後の画素値の計算に \(O(N^2)\)の計算量が必要で、それを \(N^2\) 個の画素に行うので全体の計算量は \(O(N^4)\) となる。以上より、Fourier変換を用いたフィルタ処理は、位置空間での畳み込みに比べて計算量を大幅に削減できることが分かる。

ただし、Fourier変換により高速化できるのは、あくまで畳み込み積分で表せる線形フィルタに限られる点には注意が必要である。

練習問題

異なるサイズの画像に二次元フーリエ変換 (NumPyの場合は np.fft.fft2)を適用して、サイズによって計算量がどのように変化するかを調べなさい。

3.3. 周波数空間における画像フィルタ#

画像を二次元Fourier変換によって周波数画像に変換すると、周波数領域の処理によって様々なフィルタ処理を実現することができる。ここでは、その代表例として、一部の周波数だけを取り出すことで実現される ローパス・フィルタ ならびに ハイパス・フィルタ, 加えて、畳み込まれたカーネル関数が既知の場合にフィルタ前の画像を復元できる Wienerフィルタ について紹介する。

3.3.1. ローパス・フィルタとハイパス・フィルタ#

まずは画像を二次元Fourier変換すると、周波数領域でどのような情報が得られるのかを確認しておこう。次の画像は左に示す宇宙飛行士の画像を二次元Fourier変換したときに得られる周波数成分の強度 (パワースペクトル)を示した画像である。

元画像

周波数画像

../_images/95886453d841e09b5f1e122c649c569451615dc23d00ba99d3a7343bd64e8a9b.png

../_images/735ec7d4912c70085b526b928caa0120ad42b6f0289216190390a9a1ffe7f01f.png

右側の周波数画像は画像の中心が原点となっていて、各画素の画素値がx, y各方向に沿った周波数成分 \((\omega_x, \omega_y)\) の強度を示している。一般的に、自然画像を周波数成分に直すと、原点付近の周波数の小さな成分が強く、周波数の大きな領域は弱くなっていることが多い。

ローパス・フィルタでは、低周波数の領域、すなわち原点に近い部分の周波数成分だけを残し、それ以外の成分を0とするようにマスクする。この低周波数成分だけを残した周波数画像を逆Fourier変換すると、次のように全体的に見た目がボケた画像が得られる。

マスク済み周波数画像

ローパスフィルタの結果

../_images/a54d0671594dc55f6350a63ff2668acc02c016bca0d368799ef6172fb47e868a.png

../_images/57faab658c1175298cc9842f87d546e1060e702214a06c7620e15e5213d2b01f.png

反対に、高周波数成分だけを残すようにマスク (=ハイパス)して得られた周波数画像を逆Fourier変換すると、次のようにエッジの部分だけが残ったハイパス・フィルタ結果が得られる。

マスク済み周波数画像

ハイパスフィルタの結果

../_images/397e4388e67092a9aa0aea6ce21bde4e3295da8c19dd1dc09fb19e0db8d08893.png

../_images/3e15e9ab6f6483a6986879af09df2d18e508ed322aed6332f4affe8f2520cf6f.png

3.3.2. Wienerフィルタ#

畳み込み積分により定義される線形フィルタで、なおかつフィルタカーネルが既知の場合には、フィルタ前の画像を復元することもできる。その最も基本的なフィルタがWienerフィルタである。

とある画像 \(f\) がカーネル関数 \(h\) によってフィルタ処理された結果 \(g\) は次のような畳み込み積分で表されるのだった。

\[ g = f \ast h \]

関数 \(f\), \(h\), \(g\)のFourier変換を \(F\), \(H\), \(G\)のように書くことにすると、周波数領域では単純な関数同士の積として

\[ G = F H \]

という関係式が得られる。このことからカーネル関数 \(H\) が既知であればフィルタ前の画像 \(F\) はフィルタ後の画像 \(G\) を用いて次のように書ける。

\[ F = \mathcal{F}^{-1} \left[ \frac{G}{H} \right] \]

このような原理に基づくフィルタの打ち消し処理を 逆フィルタ と呼ぶ。

デジタル画像の上でボケフィルタをかけている場合にはこの処理で大きな問題はないが、通常、画像のボケを除去したい場面では、画像に含まれるノイズの影響を考える必要がある。

例えば、カメラで撮影した写真が多少ピンボケしていて、そのボケを撮影後に除去したい場合を考えてみよう。この場合、ボケは光学的に生じるので、カメラのイメージセンサに入ってくる光が表す画像はその時点でボケている。そのボケ画像に対して、イメージセンサによる読み出し時にホワイトノイズ等のノイズが加わる。従って、ノイズ成分を \(n\), そのFourier変換を \(N\) で表すと、カメラによって得られるピンボケ画像は次のような式で表せる。

\[ G = F \ast H + N \]

このような場合に単純に逆フィルタを適用すると復元画像 \(\hat{F}\) は次のような式が得られる。

\[ \hat{F} = \mathcal{F}^{-1} \left[ \frac{G}{H} \right] = \mathcal{F}^{-1} \left[ \frac{F H + N}{H} \right] = F + \mathcal{F}^{-1} \left[ \frac{N}{H} \right] \]

一般にカーネル関数のFourier変換 \(H\) は、特に高周波数の領域で小さな値を取っていることが多い。一方で、ノイズの成分は全周波数領域の成分が混合しているものとみなせるため、 \(N\) を \(H\) で除してしまうと、ノイズの成分が意図せず強調されてしまう。

この問題を防ぐため Wienerフィルタではノイズが確率的に発生するとして、誤差の期待値を最小にするような復元用のフィルタ \(K\) を求める。理想的な復元フィルタは次の式を満たすはずである。

\[ F = K (FH + N) \]

ここで \(N\) は確率的に発生すると考えて、両辺の差を取り、その二乗誤差の期待値 \(E\) を考える。

\[\begin{split} \begin{align*} E &= \mathbb{E} \left[ | F - K (FH + N) |^2 \right] \\ &= \mathbb{E} \left[ | F |^2 \right] - 2 \mathbb{E} \left[ (FH + N)^{*} FK \right] + \mathbb{E} \left[ | K (FH + N) |^2 \right] \end{align*} \end{split}\]

なお、 \(F^{*}\) は複素数 \(F\) の複素共役を表す。この誤差の期待値が最小になるときには \(K\) での微分が0になるはずなので、次のような等式が得られる。

\[ \frac{\partial E}{\partial K} = -2 \mathbb{E} \left[ (FH + N)^{*} F \right] + 2 \mathbb{E} \left[ K |FH + N |^2 \right] = 0 \]

ここで \(F\) や \(H\) は定数であり、 \(N\) の期待値が0である、すなわち \(\mathbb{E}[N] = 0\) であることを考慮すると、上式は次のように書き直せる。

\[ K | F |^2 | H |^2 + K \mathbb{E} \left[ |N|^2 \right] = | F |^2 H^{*} \]

この式を \(K\) について整理すると、次の式が得られる。

\[ K = \frac{|F|^2 H^{*}}{|F|^2 |H|^2 + \mathbb{E} \left[ |N|^2 \right]} = \frac{1}{H} \frac{|H|^2}{|H|^2 + \frac{\mathbb{E} \left[ |N|^2 \right]}{|F|^2}} \]

通常、真の画像 \(F\) とノイズ成分 \(N\) は未知であるので、 \(\mathbb{E}[|N|^2] / |F|^2\) の値は分からない。そこで、Wienerフィルタではこの比をパラメータ \(\Gamma\) として定数で近似する。

\[ K \approx \frac{1}{H} \frac{|H|^2}{|H|^2 + \Gamma} \]

この式で周波数表示が表されるフィルタをWienerフィルタと呼ぶ。Wienerフィルタは \(\Gamma\) の値を上手く調整することで、ノイズの影響を抑えつつ画像を復元することができる。以下に、Wienerフィルタを用いた画像復元の例を示す。

元画像

ボケ + ノイズ

../_images/1788d55c276b3c888e3cc5b63d0d54460b9ae4bc38007539c13d4952a9ce6b59.png

../_images/4c13fb6030d2de20b51f274d90f95e74cd3128d30411ba0108c8ce9987470269.png

逆フィルタ

Wienerフィルタ (\(\gamma = 0.001\))

../_images/8c9e6a52e1e5c3a10ca51f1f1452d3104de491c23fe7429a3cdf39b9f88e33b7.png

../_images/fad861b1c1cfde84fdc637606da087de72f0479238aa838a28dfe95c4394c80e.png

練習問題

Wienerフィルタを実装し、 \(\gamma\) の値を変化させたときに、どのような結果の違いが得られるかを確認し、そのような違いが生じる原因を考察せよ。

目次