画像フィルタ

3. 画像フィルタ#

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

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

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

画像フィルタは、画像を分析する前に、分析に悪影響を及ぼすノイズを除去するなどの前処理に用いられることもあれば、フィルタ結果そのものを分析対象とすることもある。また、現在、画像に対する機械学習で広く用いられる畳み込みニューラルネットは、画像フィルタの一種とも考えられるなど、画像フィルタは画像処理の非常に重要な処理の一つである。

クイズ1

身近にある「フィルタ」にはどのようなものがあるだろうか？必ずしもデジタルなものばかりではなく、アナログなものも含めて思いつくものを挙げてみよう。その上で、それらのフィルタがどのような目的で用いられているかを考えてみよう。

例: コーヒーフィルタ → コーヒーを淹れる時にコーヒーの液と粉を分ける役割

クイズ2

プライバシーの保護のために、画像の中の人物の顔をぼかすことがある。しかし、このようなぼかし処理は本当にプライバシーを保護できているのだろうか？仮に保護できているとすれば、それはどのような理由によるものだろうか？逆に保護できていないとすれば、どのような理由によるものだろうか？

参考資料

  • 『ディジタル画像処理 改定第二版』 第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\) が大きいほど、より詳細度の低い画像が出力される。

表 3.1 異なる \(\sigma\) によるGaussianフィルタの結果#

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

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

../_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フィルタにおいて、変数分離を用いる場合と用いない場合とで、フィルタ処理に必要な計算量は、画像のサイズ \(N\) に応じて計算量がどのように変化するかを考察せよ。また、変数分離が可能なカーネル関数がどういう性質を持つかを考察せよ。

練習問題

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フィルタによるエッジ抽出の例を示す。本図では、正の差分値を赤色で、負の差分値を青色で表している。

表 3.2 Sobelフィルタによるx, y各方向のエッジ抽出結果#

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

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

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

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

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

エッジ強度

エッジ強度極大画素

../_images/7a808494c37cadf63d95e154534aa372509d8b87ef8f5eddfff9c96f96274ffc.png

../_images/de8e820e4f325468c3ec99c05858b8edd3fbdffdc6e388dc0bfa73d8dc56e3c0.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/de8e820e4f325468c3ec99c05858b8edd3fbdffdc6e388dc0bfa73d8dc56e3c0.png

../_images/0714f16db6a24438e35b8bdd4a0a4701b1fa51dc443cea2dda5167e14dad88ac.png

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

練習問題

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

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

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

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

メディアンフィルタ#

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

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

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

表 3.3 メディアンフィルタの例#

入力画像

メディアンフィルタ

../_images/4f2587e313b2d334bdaf40cdd42414142cec49192a5f97cc519f6726a04770f6.png

../_images/85aba38fdafd4bc87a40af07d5e9dda14ce5c0639e9c0b6b3a52861a780e1f4d.png

一方、メディアンフィルタが上手く動作するのは、ノイズが空間的に離れた箇所に生じている場合である。例えば、ごま塩ノイズ (salt-and-pepper noise) と呼ばれるノイズは、画像の中のいくつかの画素が、周囲の画素と比べて極端に明るい値を取るか、極端に暗い値を取るようなノイズであり、メディアンフィルタが有効に働く。

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

バイラテラルフィルタ#

バイラテラルフィルタ (bilateral filter)は Tomasi と Manduchi によって 1998 年に提案されたフィルタ[9]で、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ノイズを付加したノイズ画像に対するバイラテラルフィルタの適用結果を示す。

表 3.4 バイラテラルフィルタの例#

入力画像

バイラテラルフィルタ

../_images/d47f118c9968159b03cfa694bb73b087bff58cc9f0a0b6212a5c7e9447fe52ae.png

../_images/b00744632d6412fb4a2539e51213fed82cf86799dd0d3e90aef3be40a6e0935d.png

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

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

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

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

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

ノイズの種類

画像に限らず、一般のデジタルデータに混入するノイズ (雑音) には様々な種類がある。最も一般的なものはホワイトノイズや白色雑音と呼ばれるもので、ノイズ量が正規分布に従う。また、画像撮影の文脈では、より輝度の高い箇所で分散が大きくなるPoissonノイズもよく見られる。

また、近年は一般的なカメラから得られる画像では見かけることは少なくなっているが、ごま塩ノイズも一部の画像センサの不良などが原因で発生することがある。

これら、多様なノイズのそれぞれに対して、ごま塩ノイズであればメディアンフィルタ、という具合に適切なフィルタを選択することが応用上は重要だろう。

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}\]

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

\[ a_0 = \frac{1}{L} \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変換により高速化できるのは、あくまで畳み込み積分で表せる線形フィルタに限られる点には注意が必要である。

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

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

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

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

元画像

周波数画像

../_images/d47f118c9968159b03cfa694bb73b087bff58cc9f0a0b6212a5c7e9447fe52ae.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

練習問題

Gaussianフィルタのカーネル関数は平滑化の強度を調整するパラメータとして \(\sigma\) を持つ。このカーネル関数のFourier変換を求め、周波数領域での \(\sigma\) の意味を考察せよ。

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 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] \]

従って、この式により得られる \(\hat{F}\) を逆Fourier変換することで、復元画像 \(\hat{f}\) が得られる。

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

この問題を防ぐため Wienerフィルタではノイズが確率的に発生するとして、誤差の期待値を最小にするような復元用のフィルタ \(K\) を求める。理想的な復元フィルタ \(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\) として定数で近似する。

(3.6)#\[ K \approx \frac{1}{H} \frac{|H|^2}{|H|^2 + \Gamma} \]

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

元画像

ボケ + ノイズ

../_images/1788d55c276b3c888e3cc5b63d0d54460b9ae4bc38007539c13d4952a9ce6b59.png

../_images/37c76709c80aa0f61c9dd75633bae7c5e6e20e53b60983acdd120936a08c83b4.png

逆フィルタ

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

../_images/734a2dc583529e965bd3e9b82b0542f30796238e48706643640e31cf4b9c1908.png

../_images/d7e7ce0cf5f9774a5e257308385e36cdf2fa684b372d0ef5f19a54349b21cd8b.png

3.4. 画像品質の評価#

ノイズ除去やボケ除去などのフィルタによって元画像を復元した時に、その画像の品質がどの程度であるのかを定量的に評価することは、これらの手法の性能を評価するうえで非常に重要である。

例えば、ノイズのない参照用画像に対して、人工的にノイズを付与し、性能を評価したいノイズ除去フィルタを適用した結果を参照画像と比較すれば、ノイズ除去の性能を評価することができる。

このように評価用画像 \(I(x, y)\) と参照用画像 \(J(x, y)\) が与えられた時に、評価用画像が参照用画像に対して、どの程度の品質にあるかを評価する指標には様々なものがある。近年ではAlexNetやVGG等の深層ニューラルネットワークの出力する中間層の特徴から計算される LPIPS [11] (learned perceptual image patch similarity) のような指標も、画像の類似度評価に広く使われている。

以下では、画像の品質評価の標準的な手法である ピーク信号対雑音比 (peak signal-to-noise ratio, PSNR) と 構造類似度指数 (structural similarity index measure, SSIM) について説明する。

3.4.1. ピーク信号対雑音比 (PSNR)#

PSNRは、量子化されているデジタル画像に対する品質の評価尺度である。

PSNRは、特定の量子化レベルにおいて、取りうるエネルギーの最大値と、画像劣化によって評価用画像に生じしたノイズの比率を表す。

ここで、画像が取りうる最大のエネルギー \(I_\text{MAX}\) とは、量子化レベルの二乗で表され、画像が8bit符号なし整数で量子化されている場合には \(I_\text{MAX} = (2^8 - 1)^2\) で表される。

一方で、画像の劣化により各生じたノイズは、参照用画像と評価用画像の差を取ることで抽出できるので、この画素ごとの差の二乗平均を取ることで、ノイズ量を見積もる。

\[ \text{MSE}(I, J) = \sum_{x = 1}^W \sum_{y = 1}^H \| I(x, y) - J(x, y) \|^2 \]

この量を 平均二乗誤差 (mean square error, MSE) と呼ぶ。また、平均二乗誤差の平方根を取ったものを 平方根平均二乗誤差 (root mean square error, RMSE)と呼ぶ。

PSNRは、その名前にある通り、信号に対する雑音の比率なので、この比の対数を取って次のように定義される。

\[ \text{PSNR}(I, J) = 10 \log_{10} \frac{I_\text{MAX}^2}{\text{MSE}(I, J)} = 20 \log_{10} \frac{I_\text{MAX}}{\text{RMSE}(I, J)} \]

この定義を、シンプルに言い換えると、画素値の最大値に対する平均のノイズ量の対数をとって20倍したもの、であると考えられる。

例えば、量子化の階調数に対して10%の誤差が存在する場合、PSNRは20となり、1%しか誤差がなければ、PSNRは40となる。この考察からも分かる通り、PSNRが、おおむね30以上の値をとっていれば、評価用画像は十分に高いと言って良い。

3.4.2. 構造類似度指標 (SSIM)#

PSNRは画素ごとの誤差を用いた指標であるため、例えば、画像が右に1ピクセルだけずれた場合など、人の目にはほとんど同じにみえるような場合でも、大きな誤差が生じていると判断されてしまう場合がある (その場合、PSNRは小さな値を取る)。

このような問題を解決するために、画像の画素値の統計的な差を使って画像の品質を評価する指標があり、その代表例がSSIMである。

SSIMは、明度、コントラスト、構造の3つの要素を用いて、評価用画像 \(I\) と参照用画像 \(J\) の類似度を評価する。

明度に関しては、各画像の輝度の平均 \(\mu_I\), \(\mu_J\) を用いて、次のように定義される。

\[ l(I, J) = \frac{2 \mu_I \mu_J + c_1}{\mu_I^2 + \mu_J^2 + c_1} \]

コントラストに関しては、各画像の輝度の標準偏差 \(\sigma_I\), \(\sigma_J\) を用いて定義される。画像のコントラストには、いくつかの定義があるが、一般的に言って、画像内の輝度の範囲が広いほどコントラストが大きいとされる。SSIMにおいては、この輝度の標準偏差を用いて、次のようにコントラストの類似度を評価する。

\[ c(I, J) = \frac{2 \sigma_I \sigma_J + c_2}{\sigma_I^2 + \sigma_J^2 + c_2} \]

最後に、構造の類似度は、評価用画像と参照用画像の輝度の相関を用いて定義される。具体的には、2つの画像の輝度の共分散を \(\sigma_{IJ}\) として、次のように定義される。

\[ s(I, J) = \frac{\sigma_{IJ} + c_3}{\sigma_I \sigma_J + c_3} \]

ここで、 \(c_1\), \(c_2\), \(c_3\) は、分母が0になることを防ぐための定数である。一般的には、画像の輝度の量子化レベル \(L\) を用いて、 \(c_1 = (k_1 L)^2\), \(c_2 = (k_2 L)^2\), \(c_3 = c_2 / 2\) として定義されることが多い。この際、 \(k_1 = 0.01\), \(k_2 = 0.03\) がよく用いられる。

一般的なSSIMは、これら3つの類似度の適当な指数を取ったものの積を用いて、次のように表される。

\[ \text{SSIM}(I, J) = [l(I, J)]^\alpha \cdot [c(I, J)]^\beta \cdot [s(I, J)]^\gamma \]

ここで、 \(\alpha = \beta = \gamma = 1\) とすること、広く用いられているSSIMの定義式である次式が得られる。

\[ \text{SSIM}(I, J) = \frac{(2 \mu_I \mu_J + c_1)(2 \sigma_{IJ} + c_2)}{(\mu_I^2 + \mu_J^2 + c_1)(\sigma_I^2 + \sigma_J^2 + c_2)} \]

SSIMの3要素はその定義から0から1の範囲の値を取ることが分かる (相加平均と相乗平均の大小関係を用いて示せる)。従って、SSIMも0から1の範囲の値を取り、1に近いほど、評価用画像と参照用画像の類似度が高いことを示す。

なお、SSIMは輝度の平均や標準偏差を用いて定義されるため、厳密にはグレイスケール画像に対して定義される指標である。従って、カラー画像に対してSSIMを用いる場合には、RGBの各チャンネルに対してSSIMを計算し、その平均を取るなどの方法が用いられることが多い。

3.5. プログラミング演習#

レポートは テンプレートファイル を使用して作成してください。また、ファイル名は 「(7桁の学籍番号)_第x回_画像処理レポート.docx」 (xの部分は何回目の課題なのかを記入)に変更してください。

課題作成上の注意

課題を作成する際には、プログラムは別に .py ファイルで作成して、本レポートと一緒に圧縮したうえで提出してください。また、Jupyter Notebook形式のファイル (拡張子が.ipynb)のものは受け付けません。

加えて、プログラムを添付したのみで内容に関する説明や結果に関する考察のないもの、単なる結果の羅列になっているもの(またはそのように見えるもの)は採点しませんのでご注意ください。

手作り特徴量

scikit-learnには sklearn.datasets.load_digits() でロード可能な手書き数字の練習用データセットが用意されている。このデータセットを使って、単純な手作り特徴量に関する実験を行ってみよう。

今回作成する手作り特徴量は、次のような手順で作成されるものとする。

  1. 画像をグレースケール化する (上記のデータセットは既にグレースケール)

  2. 画像に対してエッジ抽出フィルタを適用する

  3. エッジ抽出フィルタの結果を二値化する

  4. 二値化された画像からエッジとして抽出された画素の数を数え、これを画像の特徴とする

この特徴を0-9の各数字に対して計算し、数字ごとのヒストグラムを作成してみよう。ここから、数字の識別可能性について考察しなさい。

Wienerフィルタのパラメータ

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

この際、異なる \(\Gamma\) で得られた復元画像の品質をピーク信号対雑音比 (PSNR) と構造画像類似度 (SSIM) の二つの指標を用いて定量的に評価すること。

目次