本文へスキップ
医用画像情報学

フーリエ変換と空間フィルタ処理(畳み込み・平滑化・エッジ検出)

医用画像情報学のなかで、出題数がいちばん多いのが画像処理です。フーリエ変換とフィルタ処理は、繰り返し問われます。

一見バラバラに見えますが、聞かれていることは実は少なく、フーリエ変換の性質と、フィルタの仕分け方の2つに集約されます。しかも仕分けは「係数の和が 1 か 0 か」を見るだけで、ほとんど片づきます。

前の記事「標本化と量子化」で出てきた「周波数成分に分解する」という考え方を、ここで正面から扱います。

フーリエ変換:実空間と周波数空間

実空間と周波数空間フーリエ変換の意味と代表的な変換ペアを示した図。実空間は位置xの関数で、画素が並んだふつうの画像である。これをフーリエ変換すると周波数空間に移り、どの細かさの縞がどれだけ入っているかが分かる。周波数空間から実空間へ戻すのが逆フーリエ変換である。代表的な変換ペアは3つある。矩形波すなわち四角い窓のフーリエ変換はsinc関数になり、幅が狭いほどsincは広がる。デルタ関数すなわち一点だけの信号のフーリエ変換はすべての周波数が等しく含まれる定数になる。余弦関数は偶関数なので実部だけの2本の線になり、虚部は出ない。フーリエ変換は線形変換であり、離散フーリエ変換も周波数領域への変換である。実空間へ戻すのは逆フーリエ変換である。フーリエ変換:実空間 ⇄ 周波数空間画像を「どの細かさの縞がどれだけ入っているか」に分解する実空間(位置 x)画素が並んだ、ふつうの画像周波数空間(空間周波数 u)どの細かさの縞が、どれだけ入っているかフーリエ変換逆フーリエ変換矩形波(四角い窓)sinc 関数デルタ関数すべて等しい定数余弦関数実部だけの2本の線※ 幅が狭いほど sinc は広がる※ 一点の信号にはあらゆる周波数が入る※ 偶関数なので虚部は出ないフーリエ変換は「実空間 → 周波数空間」への線形変換矩形波 ⇄ sinc 関数 は最頻出の変換ペア。CT の再構成でも MTF でも出てくるデルタ関数(一点)→ すべての周波数が等しく含まれる定数離散フーリエ変換も「周波数領域への」変換。実空間へ戻すのは逆フーリエ変換
フーリエ変換は実空間から周波数空間への線形変換。矩形波 ⇄ sinc 関数が最頻出の変換ペア。

フーリエ変換は、画像を「どの細かさの縞が、どれだけ入っているか」に分解する操作です。式は次のとおりです。

F(u)=∫−∞∞f(x) e−i2πux dxF(u)=\int_{-\infty}^{\infty} f(x)\,e^{-i2\pi ux}\,dx

f(x)f(x) が実空間(位置)の関数、F(u)F(u) が周波数空間の関数です。逆向き、つまり周波数空間から実空間に戻すのが逆フーリエ変換です。

離散フーリエ変換も同じで、周波数領域への変換です。「実空間領域への変換である」と書いてあれば誤りです。

式の中の e−i2πuxe^{-i2\pi ux} を核関数といいます。これはオイラーの公式

eiθ=cos⁡θ+isin⁡θe^{i\theta}=\cos\theta+i\sin\theta

に θ=−2πux\theta=-2\pi ux を代入して、cos⁡\cos が偶関数、sin⁡\sin が奇関数であることを使うと、

e−i2πux=cos⁡(2πux)−isin⁡(2πux)e^{-i2\pi ux}=\cos(2\pi ux)-i\sin(2\pi ux)

と書けます。sin⁡\sin の項に必ず虚数単位 ii が付くこと、そして符号がマイナスであること(プラスなら逆変換の核関数)の2点が問われます。

覚えておきたい変換ペアは3つです。

  • 矩形波(四角い窓)→ sinc 関数。いちばん大事なペアで、CTの再構成でも MTF でも出てきます。窓の幅が狭いほど sinc は広がります
  • デルタ関数(一点だけの信号)→ すべての周波数が等しく含まれる定数。「一瞬の信号にはあらゆる周波数が入っている」と考えると腑に落ちます
  • 余弦関数 → 実部だけ。偶関数なので虚部は出ません。虚部が出るのは奇関数である正弦関数のほうです

フーリエ変換の性質

フーリエ変換の性質フーリエ変換の性質を4つ示した図。第1に線形変換であり、a倍のf(x)とb倍のg(x)の和はaF(u)とbG(u)の和になる。定数倍も足し算もそのまま通る。第2に平行移動しても絶対値は変わらない。f(x−a)の絶対値はF(u)と同じで、変わるのは位相だけであり振幅スペクトルは不変である。第3に実空間の畳み込みは周波数空間の掛け算になる。和のフーリエ変換が積というのは誤りで、積になるのは畳み込みのほうである。第4に横に縮めると周波数側は広がる。f(2x)のフーリエ変換は2分の1のF(uを2で割ったもの)であり、F(2u)ではない。面積が等しいという記述は誤りで、正しくはf(x)の面積がF(0)、F(u)の面積がf(0)である。畳み込み定理が使えるから周波数空間での処理は速い。実空間では重畳積分、周波数空間では単純な掛け算である。フーリエ変換の性質:試験に出る4つとくに「畳み込み ⇄ 掛け算」の対応が土台になる線形変換であるa f(x) + b g(x) → a F(u) + b G(u)定数倍も足し算も、そのまま通る平行移動しても、絶対値は変わらないf(x − a) → 絶対値は F(u) と同じ変わるのは位相だけ。振幅スペクトルは不変実空間の畳み込み = 周波数空間の掛け算f(x) * g(x) → F(u) × G(u)「和のフーリエ変換が積」は誤り。積になるのは畳み込みのほう横に縮めると、周波数側は広がるf(2x) → ½ F(u ÷ 2)F(2u) ではない。縮めたぶん周波数は高いほうへ伸びる「面積が等しい」は誤り。正しくは f(x) の面積 = F(0)、F(u) の面積 = f(0)畳み込み定理が使えるから、周波数空間での処理は速い(掛け算だけで済む)実空間では重畳積分(畳み込み)、周波数空間では単純な掛け算。この対応が全体の土台
とくに実空間の畳み込み = 周波数空間の掛け算が土台。平行移動しても絶対値は変わらない。

性質は4つ押さえれば足ります。

1. 線形変換である

af(x)+bg(x) ⟶ aF(u)+bG(u)a f(x)+b g(x)\ \longrightarrow\ a F(u)+b G(u)

定数倍も足し算も、そのまま通ります。2f(x)2f(x) の変換は 2F(u)2F(u) であって、12F(u)\frac{1}{2}F(u) ではありません。

2. 平行移動しても絶対値は変わらない

f(x)f(x) と f(x−a)f(x-a) のフーリエ変換は、絶対値(振幅スペクトル)が同じです。変わるのは位相だけです。

3. 実空間の畳み込み = 周波数空間の掛け算

f(x)∗g(x) ⟶ F(u)×G(u)f(x)*g(x)\ \longrightarrow\ F(u)\times G(u)

これが畳み込み定理で、画像処理の土台です。実空間でフィルタをかける操作(重畳積分)は、周波数空間では単純な掛け算になります。だから周波数空間での処理は速いのです。

なお、和のフーリエ変換は和(F+GF+G)です。「f(x)+g(x)f(x)+g(x) のフーリエ変換は F(u)⋅G(u)F(u)\cdot G(u)」という選択肢は誤りで、積になるのは畳み込みのほうです。

4. 横に縮めると、周波数側は広がる

f(2x)f(2x) の変換は 12F(u/2)\frac{1}{2}F(u/2) です。F(2u)F(2u) ではありません。実空間で縮めたぶん、周波数は高いほうへ伸びます。

また、「f(x)f(x) の面積と F(u)F(u) の面積は等しい」も誤りです。正しくは f(x)f(x) の面積 = F(0)F(0)、F(u)F(u) の面積 = f(0)f(0) という関係です。

空間フィルタは「重ねて、掛けて、全部足す」

空間フィルタの計算空間フィルタの計算手順を示した図。もとの画素値の3掛ける3は上から3と9と4、7と6と2、5と8と1である。フィルタすなわちカーネルは上からマイナス1と0とマイナス1、0と5と0、マイナス1と0とマイナス1で、係数の合計は5引く4で1になる。計算はフィルタの同じ位置どうしを掛けて9個すべてを足すだけで、フィルタが0の場所は計算しなくてよい。この例ではマイナス1掛ける3でマイナス3、マイナス1掛ける4でマイナス4、5掛ける6でプラス30、マイナス1掛ける5でマイナス5、マイナス1掛ける1でマイナス1となり、合計は17である。効くのは四隅と中心の5か所だけである。係数の合計が1なら画像全体の明るさは変わらない。この例は中心がプラス5、四隅がマイナス1で、周囲との差を強調する鮮鋭化フィルタである。空間フィルタ:重ねて、掛けて、全部足す注目画素のまわり3×3と、フィルタの同じ位置どうしを掛けて合計するもとの画素値(3×3)394762581フィルタ(カーネル)-10-1050-10-1係数の合計 = 5 − 4 = 1掛けて、全部足す(−1)×3 = −3(−1)×4 = −4( 5 )×6 = +30(−1)×5 = −5(−1)×1 = −1合計 = 17色の付いた5か所だけが効く空間フィルタは「重ねて、掛けて、全部足す」だけ注目画素を中心に3×3を切り出し、フィルタの同じ位置どうしを掛けて、9個すべてを足すフィルタが 0 の場所は計算しなくてよい。この例で効くのは四隅と中心の5か所だけ係数の合計が 1 なら、画像全体の明るさは変わらないこの例は「中心+5、四隅−1」。周囲との差を強調する鮮鋭化フィルタ
注目画素のまわり3×3と、フィルタの同じ位置どうしを掛けて、9個すべてを足す。フィルタが 0 の場所は計算しなくてよい。

空間フィルタ(重畳積分・畳み込み)の計算は、手順としては単純です。

  1. 注目画素を中心に 3×3 を切り出す
  2. フィルタの同じ位置どうしを掛ける
  3. 9個すべてを足す

図の例では、フィルタが 0 の場所が4つあるので、実際に効くのは四隅と中心の5か所だけです。

(−1)(3)+(−1)(4)+(5)(6)+(−1)(5)+(−1)(1)=30−13=17(-1)(3)+(-1)(4)+(5)(6)+(-1)(5)+(-1)(1)=30-13=17

計算問題では、注目画素が端にあるかどうかを最初に確認してください。5×5 の中央(i=3, j=3)であれば、はみ出しの処理を考えずに済みます。

もうひとつ、メディアンフィルタだけは掛け算をしません。範囲内の画素値を小さい順に並べて、ちょうど真ん中の値を採るフィルタです。平均ではありません。9個なら5番目の値です。

係数の和で仕分ける

係数の和によるフィルタの仕分けフィルタの係数の和による仕分けを示した図。係数の和が1のものは平滑化フィルタで、3掛ける3がすべて1の移動平均フィルタが例である。あとで9で割る。係数の和が0のものは微分すなわちエッジ検出フィルタで、上から0と1と0、1とマイナス4と1、0と1と0というラプラシアンフィルタが例である。これは2次微分にあたる。1次元フィルタのマイナス1、1、0を台形の信号にかけると、平らな部分はすべて0になり、立ち上がりでは負、立ち下がりでは正の値が出る。係数の和を見るだけで平滑化か微分かが分かる。和が1なら平滑化で移動平均や加重平均やガウシアン、和が0なら微分やエッジ検出でソーベルやプレヴィットやラプラシアンである。微分は変化しているところだけを取り出し、向きを符号で表す。係数の和が 1 なら平滑化、0 なら微分これだけで、ほとんどのフィルタ問題は仕分けできる係数の和 = 1 → 平滑化111111111移動平均フィルタ(全部 1、あとで 9 で割る)係数の和 = 0 → 微分(エッジ)0101-41010ラプラシアンフィルタ(2次微分)1次元フィルタ(−1, 1, 0)を台形にかけるともとの信号(台形)フィルタ後(微分)立ち上がり → 負平らな部分 → 0立ち下がり → 正係数の和を見るだけで、平滑化か微分かが分かる和 = 1 → 平滑化(移動平均・加重平均・ガウシアン)和 = 0 → 微分・エッジ検出(ソーベル・プレヴィット・ラプラシアン)微分は「変化しているところだけを取り出し、向きを符号で表す」
係数の和が 1 なら平滑化、0 なら微分(エッジ検出)。微分フィルタは、平らな部分を 0 にして、変化の向きを符号で表す。

フィルタの種類を見分ける、いちばん速い方法がこれです。

  • 係数の和 = 1 → 平滑化。全体の明るさが変わらない。移動平均、加重平均、ガウシアン
  • 係数の和 = 0 → 微分(エッジ検出)。平らな部分の出力が 0 になる。ソーベル、プレヴィット、ラプラシアン

なぜこうなるのか。一様な面(どの画素も同じ値 vv)にフィルタをかけると、出力は v×(係数の和)v \times (\text{係数の和}) になります。和が 1 なら vv のまま、和が 0 なら 0 です。つまり和が 0 のフィルタは、変化のないところを消すわけで、これがまさに微分のふるまいです。

図の下段は、台形の信号に 1次元フィルタ (−1,1,0)(-1, 1, 0) をかけた結果です。

  • 平らな部分はすべて 0
  • 立ち上がりでは負、立ち下がりでは正と、符号が逆になる

「変化しているところだけを取り出し、変化の向きを符号で表す」という、微分そのものの性質が見てとれます。

微分の次数も問われます。ソーベル・プレヴィットは1次微分で、傾きの大きさと向きが得られます。ラプラシアンは2次微分で、エッジの位置がゼロ交差として現れます。

フィルタ早見表:実空間と周波数空間

フィルタの早見表フィルタの早見表。実空間で使うフィルタは重畳積分で処理する。平滑化は係数の和が1で、移動平均や加重平均は単純にならし、ガウシアンフィルタは自然にぼかす。非線形のメディアンフィルタはスパイク除去とエッジ保存に使う。微分すなわちエッジ検出は係数の和が0で、ソーベルとプレヴィットは1次微分、ラプラシアンは2次微分である。鮮鋭化にはアンシャープマスキングすなわちボケマスク処理がある。周波数空間で使うフィルタは掛け算で処理する。ローパスすなわち平滑化にはバターワースがあり遮断周波数と次数で調整する。ガウシアンやハニングは核医学で多用される。高周波を強調するrampフィルタはFBPのボケ取りに使うが、高周波を強調するのでノイズも増えるためローパスと併用する。復元にはウィーナーフィルタがある。実空間の重畳積分は周波数空間の掛け算にあたる。周波数空間でのフィルタ処理が重畳積分で行われるという記述は誤りである。バターワースはローパスでありハイパスではない。ガウシアンは平滑化で輪郭は強調されない。繰り返し数とサブセット数はOSEMのパラメータでウィーナーとは無関係である。フィルタ早見表:実空間と周波数空間どこで使うフィルタか、何をするフィルタかの2軸で整理する実空間で使う(重畳積分)平滑化(係数の和 = 1)移動平均・加重平均単純にならすガウシアンフィルタ自然にぼかす非線形メディアンフィルタスパイク除去・エッジ保存微分=エッジ検出(和 = 0)ソーベル・プレヴィット1次微分ラプラシアン2次微分鮮鋭化アンシャープマスキングボケマスク処理周波数空間で使う(掛け算)ローパス(平滑化)バターワース遮断周波数と次数ガウシアン・ハニング核医学で多用高周波を強調ramp フィルタFBP のボケ取り復元ウィーナーフィルタノイズを考えた復元ramp は高周波を強調するのでノイズも増える → ローパスと併用実空間の重畳積分 = 周波数空間の掛け算「周波数空間でのフィルタ処理は重畳積分で行う」は誤り。周波数空間では単純な掛け算バターワースはローパス。ハイパスではないガウシアンは平滑化。輪郭は強調されず、むしろ不明瞭になる繰り返し数とサブセット数は OSEM(逐次近似)のパラメータ。ウィーナーとは無関係
実空間の重畳積分 = 周波数空間の掛け算。バターワースはローパスで、ハイパスではない。

フィルタは「どこで使うか」と「何をするか」の2軸で整理すると混乱しません。

実空間で使うフィルタは重畳積分(畳み込み)で処理します。平滑化(移動平均・加重平均・ガウシアン)、非線形のメディアン、微分系(ソーベル・プレヴィット・ラプラシアン)、鮮鋭化のアンシャープマスキング(ボケマスク処理)です。

周波数空間で使うフィルタは掛け算で処理します。核医学でよく出るのがここです。

  • バターワースフィルタ … ローパス(平滑化)。遮断周波数と次数の2つのパラメータで、なだらかさを調整します。ハイパスではありません
  • ガウシアン・ハニング … 同じくローパス
  • ramp フィルタ … 高周波ほど強調する特性。単純逆投影で生じる 1/r1/r のボケを打ち消すために使い、これがフィルタ補正逆投影法(FBP)です。ただし高周波を強調するのでノイズも増えるため、実際にはバターワースなどのローパスと組み合わせます
  • ウィーナーフィルタ … ノイズを考慮した復元フィルタ

まぎらわしいのが、繰り返し数(iteration)とサブセット数です。これは OSEM などの逐次近似法のパラメータで、ウィーナーフィルタとは無関係です。

ヒストグラムから分かること・分からないこと

ヒストグラムから分かることヒストグラムから分かることと分からないことを示した図。左には画素の並びだけが違う2枚の画像がある。ひとつはなめらかな画像、もうひとつは同じ画素値をばらばらに並べ替えた画像で、見た目はまったく違う。しかし右に示すとおりヒストグラムは完全に同じである。ヒストグラムから求められるのは期待値すなわち平均、分散や標準偏差、歪度や尖度、中央値である。ヒストグラムから求められないのは鮮鋭度、MTF、ノイズの空間周波数特性で、どれも位置の情報が要る。ヒストグラムは位置の情報を完全に捨てており、どの画素値が何個あったかだけの集計である。だから画素値の統計は出せても、鮮鋭度のような空間の情報は出せない。ヒストグラムは「位置の情報」を捨てている同じヒストグラムでも、並べ替えれば見た目はまったく変わる画素の「並び」だけが違う2枚なめらかな画像ばらばらな画像見た目はまったく違うでもヒストグラムは同じ25665462画素値(暗い ← → 明るい)どちらの画像も、この1つのヒストグラムヒストグラムから求められる期待値(平均)分散・標準偏差歪度・尖度中央値ヒストグラムからは求められない鮮鋭度MTFノイズの空間周波数特性(どれも位置の情報が要る)ヒストグラムは「位置の情報」を完全に捨てているどの画素値が何個あったかだけの集計。どこにあったかは残らないだから画素値の統計は出せても、鮮鋭度のような空間の情報は出せない
ヒストグラムは位置の情報を完全に捨てている。だから鮮鋭度は求められない。

ヒストグラムは「どの画素値が何個あったか」の集計です。ここで決定的なのは、画素がどこにあったかという情報が、完全に捨てられていることです。

図の2枚の画像は、まったく同じ画素値の集まりを、並べ替えただけのものです。見た目はまるで違いますが、ヒストグラムは完全に一致します。

したがって、

  • 求められる … 期待値(平均)、分散・標準偏差、歪度・尖度、中央値。どれも画素値の統計です
  • 求められない … 鮮鋭度、MTF、ノイズの空間周波数特性。どれも隣の画素との関係、つまり位置の情報が要ります

試験で問われやすいポイント

論点覚える形
フーリエ変換の向き実空間 → 周波数領域。戻すのが逆変換
線形性線形変換。2f(x)→2F(u)2f(x) \to 2F(u)
核関数cos⁡(2πux)−isin⁡(2πux)\cos(2\pi ux) - i\sin(2\pi ux)。ii が付き、符号はマイナス
変換ペア矩形波 → sinc、デルタ → 定数、余弦 → 実部だけ
平行移動絶対値は不変。変わるのは位相だけ
畳み込み定理実空間の畳み込み = 周波数空間の掛け算
和の変換和は和。積になるのは畳み込み
縮小f(2x)→12F(u/2)f(2x) \to \frac{1}{2}F(u/2)。F(2u)F(2u) ではない
面積f(x)f(x) の面積 = F(0)F(0)、F(u)F(u) の面積 = f(0)f(0)
空間フィルタの計算同じ位置どうしを掛けて9個すべて足す
メディアン並べ替えて真ん中。平均ではない
係数の和 = 1平滑化
係数の和 = 0微分・エッジ検出
1次微分ソーベル・プレヴィット
2次微分ラプラシアン
バターワース周波数空間のローパス
ramp高周波強調。FBP のボケ取り。ノイズも増える
繰り返し数・サブセット数OSEM のパラメータ
ヒストグラム位置の情報が残らない。鮮鋭度は出せない

まとめ

  • フーリエ変換は実空間から周波数空間への線形変換。離散フーリエ変換も同じ向き
  • 核関数は cos⁡(2πux)−isin⁡(2πux)\cos(2\pi ux) - i\sin(2\pi ux)。sin⁡\sin に ii が付き、符号はマイナス
  • 変換ペアは矩形波 ⇄ sinc、デルタ ⇄ 定数、余弦 ⇄ 実部だけ
  • 平行移動しても絶対値は不変。変わるのは位相だけ
  • 実空間の畳み込み = 周波数空間の掛け算。これが画像処理の土台
  • 空間フィルタは同じ位置どうしを掛けて全部足す。メディアンだけは並べ替えて真ん中
  • 係数の和が 1 なら平滑化、0 なら微分。これでほとんど仕分けできる
  • バターワースはローパス。ramp は高周波強調でノイズも増えるためローパスと併用する
  • ヒストグラムは位置の情報を捨てているので、鮮鋭度や MTF は求められない

次は、ここまでとは毛色が変わります。人が画像を見て診断できたかどうかを測る指標、ROC 解析です。

この分野のほかの記事

← トップに戻る | 医用画像情報学の記事一覧