X線写真のMTFをPythonで計算その１
はじめに
このページはX線写真におけるMTF(modulation transfer function:変調伝達関数)について、まとめたものです。
その１にてMTFの概要を説明をしています。
その２にてPythonを使って実際に計算をしています。
参考にしてください。
MTFについては、下記の書籍を参考にしました。
X線の医用画像に関する本はいくつか読みましたが、MTFの記述はこれが詳しいと思います。
１．空間周波数とはなにか
MTFを説明する前に空間周波数について説明します。
一般的に、周波数といえば、音や電波、機械的な振動の世界ではおなじみのHzが存在します。Hzという単位はcycles/secを意味します。
画像の世界では、時間的な変化を意味するHzではなく、空間的に変化する信号として、空間周波数が扱われます。
単位は、cycles/mmやline-pairs/mm(LP/mm)が使われます。
上記の単位をみると空間周波数が何を意味するか、理解しやすいです。
ここで、画像はさまざまな正弦波と余弦波の合成で表せます。（フーリエの定理）
そのため、１枚の画像にも高空間周波数が支配的な領域と、低空間周波数が支配的な領域が複雑に存在します。
空間周波数成分の違いと解像度の違いは対応しているため、高い解像度が欲しい場合は、高い空間周波数成分で十分に応答する必要があります。
２．MTFとはなにか
MTFは、その空間周波数毎の応答性を示したものです。
実用的には、書籍にある以下の説明がわかりやすいです。
MTFとは、各空間周波数における応答値の関数である。
具体的な方法として、インパルス応答をフーリエ変換した後、その絶対値をゼロ周波数の値で正規化することで得られる
MTFは空間周波数の増加とともに減少し、高い周波数成分のレスポンスが悪くなるのが普通です。
それは、細かい信号ほど画像に表しにくいことを意味しており、経験的にも一致すると思います。
ざっくり言ってしまえば、撮影する装置のMTFを得ることで、その装置がどれだけ細かいものをきちんと撮影できるかが分かることになります。
３．最高空間周波数とナイキスト周波数
今後のために、ナイキスト周波数について説明をします。
ナイキスト周波数は、サンプリング間隔d（画素ピッチにあたる）を定めた場合に、忠実に再現することができる最高の周波数のことで、1/(2d)と定義されます。
ナイキスト周波数は、装置が忠実に再現できる最高の周波数ということになります。
しかし、被写体には、様々な周波数成分が存在します。
そのため、撮影した画像には、ナイキスト周波数よりも高い周波数成分が存在することが多いです。
このナイキスト周波数より高い成分は、エリアシング誤差、または折返し誤差と呼ばれる誤差を発生させ、画像にシャギーのような模様やモアレを表示させる原因となります。
対策として、実用上十分なサンプリング間隔（画素ピッチを細かくする）にしたり、アンチエイリアス処理が施されます。
４．DRにおけるMTF
X線写真は一昔前はフイルムで撮影されていましたが、今ではデジタル的に撮影できるDRが一般的になりました。
ここでは、DRでのMTFについて説明します。
DRになるとデジタル化されたがゆえに、サンプリングの標本化による誤差を含みます。
上で説明したナイキスト周波数以上の空間周波数のMTFが測定することができず、ナイキスト周波数による折返し誤差も生じ、正しい測定が困難となります。
その結果、X線フイルムのMTFとの比較が困難となるという問題も発生しました。
そこで、DRシステムでもX線フィルムのMTFと比較できるようにpresampled MTFが提案されました。
測定方法は非常に簡単で、金属エッジを碁盤目のように並んでいる画素ピッチに対してやや傾ける（1.5〜3.0°が推奨）だけです。
その後、得られた画像をもとに、計算をすることで、MTFの結果がでます。
傾けることにより、サブピクセルレベルで金属エッジ部の画像を推定することができるため、サンプリングによる誤差を無視できるようになるというイメージです。
具体的なところは、下の計算を見ながら説明していきます。
（エッジだけでなく、スリットやチャートを用いる計算方法もあります。）
５．MTFの計算方法
実際の計算は、その２（近日公開）に譲りますが、MTFの計算の流れの概要を説明したいと思います。
ここでは、エッジ画像を使ってpresampled MTFを計算する場合を想定します。
金属エッジを使って得られたエッジ画像から、合成プロファイル（ESF）を作成します。そのデータを微分し、針のように凸になったエッジ拡がり関数（LSF）を作成します。次に、それをフーリエ変換と正規化し、求めたいMTFが出力されます。イメージとしては下記のような流れです。
実際には、画像データの線形化とX線量への変換処理が必要になります。また、エッジから微分によってESFを求める場合は、ノイズを除去する処理を行うことが書籍では推奨されています。
ちなみに、上の図でLSFの右側が乱れているのは、明るい方が画像のノイズが多いためです。