第19回 画像の幾何学変換処理入門（その1）で学ぶ！

今回は画像の幾何学的変換の入門的な部分であるを画像の線形変換の基礎的な事項を中心に説明します！

画像をピクセルごとの数値の行列を羅列していると考えた場合、簡単な線形代数を適用して一次変換処理が行えます。

行列演算が可能ということは、移動・回転・拡大・縮小・鏡面変換・合成変換などの線形変換による単純な一次変換から、固有値や固有ベクトルなどの特徴量計算まで、多彩な処理が可能となります。

画像変換処理と行列演算

まず、ベクトル・行列演算を復習していきます。

はじめに次の図をご参照ください。

例えば 画像間演算などで、2つの画像を足し算する場合は、単に個々のピクセル同士の足し算になります。

これは画像も行列も同じであり、上図に示した行列の和と同じ処理です。 

線形変換では主に行列やベクトルの積としての考え方が基本で、図のベクトルの積および行列の積のような計算が使用されます。

ベクトルの積は、この例ですと、

(1・4)＋(2・5)＋(3・6) = 32

行列の積は、

1・4＋4・5＋7・6 = 66

（繰り返し）

（図の引用：http://www.geisya.or.jp/~mwm48961/kou2/matrix2.html）

線形変換であるがゆえに、基本は高校の授業で習ったような一次変換です。

一次変換の基本式を次に示します。

ImageJを用いた機何変換処理

1.画像の拡大と縮小

二次元画像のx成分とy成分を独立して拡大縮小する場合の処理方法を示します。

画像をX軸方向にSx、Y軸方向にSy、だけ拡大縮小する処理は以下のとおりです。 

次の図は、オリジナル画像を二倍拡大する処理の例です。ImageJでは、Image>Scaleでスケール値を2.0(2倍拡大)や0.5(1/2縮小)などの操作ができます。Image>Zoomや、ImageJメインメニュー内の虫眼鏡マークのクリック操作でも同じような処理が可能です。

 （オリジナル画像とScale設定ウィンドウ）

（0.5倍に縮小、マトリクス数も半分に）

2.回転

原点を中心に、反時計回りに角度θだけ回転する変換式は次の通りです。

ImageJでの操作方法は、Image>Transform>Rotateで、角度θ（ここでは45）の入力をします。

3.鏡映

鏡映とは、ある任意の軸を空間軸に設定し、その軸に対象に画像を折り返す変換です。

X軸を基準にして折り返す変換は次のように表されます。

Y時の場合は次のようになります。

操作はImage>Transformで、Flip Horizontally（X軸基準鏡映：左右反転）、Flip Vertically（Y軸基準鏡映：上下反転）などの操作ができます。

（左右反転例）

4.スキュー

矩形画像を平行四辺形に変形するような変換を「スキュー」「せん断変換」といいます。

X軸方向へせん断する場合の変換式は次のようになります。

このうち、bは傾ける角度をθとすると、b=tanθとなります。

Y軸方向へのせん断は次のようになります。

ImageJでの操作は、Meijering氏が作成したTransformJを用います。

（http://imagescience.org/meijering/software/transformj/）

2015/9時点では、このプラグインはver1.5以上でないと動作しないというメッセージがでますので、Fijiという次世代のImageJを使って操作した例を示します。

Fijiには、デフォルトでこのプラグインが入っていますので、そのままPluginsからTransform>TransformJを選択して、Affineを使ってみます。

Affineを起動すると、1次変換のマトリクスをロードして、実行するウィンドウが表示されます。

はじめて起動するときには、このマトリクスファイルが作られていないので、新しく作成（create）してみます。createを起動すると、マトリクスが入力できるウィンドウが表示されます。ここでは、斜めに変換するような1次変換を行ってみます。

これをSaveして、先のウィンドウにロードして実行すると、次のようになります。

以上、今回は普段何気なく使っている機能の幾何変換処理を説明しました。

次回も続けて画像の幾何変換について説明します！

Reference

• 「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・7) 2006, p102-103」

第18回 ImageJを用いたCT画像再構成で学ぶ！

前回はフーリエ変換についてご紹介しました。実際に画像の再構成はどのように計算されるのかを見てみると、理解もより深まります。

今回はCT画像の再構成パラメータをImageJのインターフェースで見てみます。

ImageJによるCT画像再構成シミュレーション

ほとんどのCT用再構成シミュレーションはほぼ共通のパラメータを持っています。

標準パラメータは次のようなものがあります。

①Number of Ray

画像の投影データ数。元画像から逆算する場合は、マトリクス数で代用することが多い。

②angular increment (投影角間隔)

X線やγ線の投影角を1°や0.5°ごとに収集する。

③Number of view

投影角と投影間隔との商で求められる。（投影角180°で投影間隔が0.5°のとき、360ビュー）

④View angle（ビュー角）

簡単なシミュレータの場合、投影角は0～180°。ファンビームでは扇状にX線やγ線が照射されるが、その扇の中心線の角度で定義。

⑤歳構成フィルタ

投影データにかけるフィルタには、Ramp、Shepp-Logann、Hann、Hamming、Cosine、Blackman。

⑥その他

上記の標準パラメータに加えて、再構成フィルタのカットオフ周波数の設定や再構成時のピクセル位置補正のために、nearest neighbor、linear、splineなどの処理を施す。

これらのパラメータに留意しながら、Damien Farrell氏がRadon Transformというプラグインを公開しているので、これで学ばせていただきましょう。（実際の操作は、第16回参照）

臨床に応用されているCT画像再構成の際に考慮されるパラメータについて

CT画像再構成を臨床に応用してみると気づくことですが、単にプロジェクションデータにフィルタで再構成しただけでは、その再構成画像のピクセル値が「CT値（Hunsfield unit）」になりません。

水と空気に合わせたCT値補正が必要になります。これを一般的に再構成時のキャリブレーションと言ったりします。

CT画像は、X線の吸収係数の分布を示していて、CT値はこの吸収係数として表されます。
CT値は水を「0」、空気を「-1000」としたときの生体組織の線減弱係数で、次の式で表されます。

CT値= (μa - μw ) / μw × K

μaは組織の減弱係数、μwは水の減弱係数、Kは定数で、通常は「1000」です。

CT値は、投影する対象の主要な組成物質の密度吸収係数と実効原子番号によって左右されます。

吸収係数の一般式は次の通りです。

記号μは線吸収係数、ρは密度（g/cm3）、Eは実効エネルギー、Zは実効原子番号を示します。

この式は一般的なX線撮影装置での投影にも当てはまります。

また、密度の高い物質をX線が通過する場合、低エネルギー部分が除去されて、周辺物質よりも泉質が硬くなるためにアーチファクトも生じてしまいます。

これはビームハードニングと呼ばれます。

このように、最終的なCT画像を生成するまでに、再構成に影響するさまざまな因子があることがわかります。

次回は、画像の幾何学変換処理について述べます！

Reference

• 「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・6) 2006, p82-83」

第17回 CT原理入門（2）-フーリエスライス理論について-で学ぶ！

前回はラドン変換における基礎的な事項について述べました。

今回は前回触れなかったバックプロジェクションの際の補足として、フーリエスライス理論について概要を説明していきます。

初歩的な、CT画像とフーリエ変換

フーリエスライス理論を一言で説明すると、「平行投影の一次元投影プロファイルの1次元フーリエ変換は、対象の二次元フーリエ空間の1方向と一致する」というものです。

筆者が一番理解しやすかった図が、前回もご紹介したこちらの参考書の図です。

（Excelによる画像再構成入門より引用）

これを踏まえて、本稿では、Kak&Slaneyの参考書（http://www.slaney.org/pct/pct-toc.html）をもとに、フーリエ変換とCT画像の関係を紐解いていきたいと思います。

まず、対象f(x,y)の二次元フーリエ変換は次の式で表されます。

同様に、投影角度をθ、そのときの回転中心からt離れたプロジェクションをPθ(t)としたフーリエ変換は、次の式で表されます。

この2つの式の構成がわかったところで、話を単純にするために、一投影単位で見ていくことにします。投影角度0°のプロジェクションを二次元フーリエ空間上の垂直周波数0のラインを選択すると、二次元フーリエ変換式が次の式に変形できます。

この式で強調している部分は、投影角度が0°のときはy方向周波数はこの計算に寄与しないということだけです。

そして、投影角度0°で、投影関数はy軸に沿った線積分で表されるので、次のようになります。

これを先ほどの式に代入すると、次のように表すことができます。

これで、ある角度のプロジェクションのフーリエ変換が、対象の二次元フーリエ変換の1ラインとなることが理解できます。

そして、他の角度から投影されたプロジェクションについては、対応する角度θ分だけ二次元フーリエ空間で放射状に配置されます。

この図を見ているとわかりますが、二次元フーリエ空間での中心から離れた部分（外周）は、データの密度が粗いですね。この外周部分は、高周波成分でしたね。

よく、CT画像の画質を評価するときに、画像の中心から離れた外周にノイズが多かったりするのは、CT画像が主に低周波成分で構成されていて、それ以外の高周波成分は、二次元フーリエ空間上でデータが疎になっているので、情報量が低くなり、実画像でのノイズが増えるんです。（実際には、この部分に然るべきフィルタがかけられ、補正されます。）

逆フーリエ変換

このように、二次元フーリエ空間が放射状に埋められていることがわかりました。あとはこれを逆フーリエ変換すると、もとの画像が求められます。

逆フーリエ変換の一般式は、次の式で表されます。

参考までに、実際コンピュータ計算を行う際には、対象f(x,y)は離散値で扱われます。この場合、もし対象f(x,y)の定義位置が下図に示す条件であるならば、逆フーリエ変換は次のように表されます

離散逆フーリエ変換の一般式は次の式です。

この計算、、たいへんですね、、。

でも、すでにImageJプラグインに実装されていますので、ソースコードから学ばせてもらうのが一番効率が良さそうです。

研究成果や教材をオープンにしてくれる世界中の研究者に感謝。

次回もCTの原理について触れていきます！

Reference

• 「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・4) 2006, p78-79」

第16回 投影による画像再構成法の原理についてで学ぶ！

今回は、CTの断層像（コンベンショナル）を得るための画像再構成の原理について、初歩的な内容から説明していきます！

ラドン変換

次の図のように、ある対象物（ここでは、CT装置のガントリー中心に置かれている被写体と思ってください）が存在する空間上にx-y座標空間を設定し、その空間上で任意角度θの方向（図中のn軸の方向）からX線が照射されたと仮定すると、その投影面SにX線強度値の分布が投影（projection）されます。

このような任意の方向のprojectionデータを集めると、投影した対象物の内部構造まで可視化できそうな気がしてきます。そんな知的好奇心を持ちながら、実際は、どうなるのかみていきましょう。

実際のラドン変換の数式は、次のようになります。

このラドン変換は、下図左上の関係を数式で表したものです。g(s,θ)は下図右上のサイノグラムの要素になります。

（Excelによる画像再構成入門より引用）

ラドン変換によって変換された座標は、1つの投影方向を表しています。これがすべての投影方向に繰り返されることになります。そしてそれぞれの投影分布を正規化し、回転角度（θ）ごとに順に並べたものが、サイノグラムと呼ばれるものです。

この元画像（本当は、何かの物体）を二次元（x-y平面）の任意の角度から複数投影して、サイノグラムを得る変換を一般的にラドン変換と言ってもいいと思います。

ではこのラドン変換で、変な画像（サイノグラム）を作成しましたが、これが何の役に立つのかをみていきましょう。

バックプロジェクション

ラドン変換を用いれば、元画像に対する投影分布を信号強度として正規化し、回転角度ごとに順に並べ、サイノグラムを作ることができます。

このサイノグラムをx-y平面に戻すと、対象物の内部の情報を含んだ断層面が得られます。このx-y平面に戻す方法を、逆投影と呼びます。

逆投影結果画像b(x,y)は、投影方向ごとにマトリクス上に配置した投影変換結果（サイノグラム）をp(k,m)とすると、回転角度ごとのラドン変換の結果から、次の式で表されます。

この処理をバックプロジェクション（逆投影）といいます。

しかし、このb(x,y)は単にすべて回転角度（投影角度）の投影データを加算して二次元分布に変換したのみで、元の画像f(x,y)とはまだ一致しません。元の画像と比べるとぼけた画像になります。このボケの原因は、投影データ間の重なりです。

単に矩形のビームを重ねていくと、中心周辺が点状になり、ボケが広がっていきます。

（参考記事から引用）

そのため、数式的に表現すると、単にバックプロジェクションした画像はボケを含む次の式で表されます。

このh(x,y)がボケの関数（点広がり関数）です。

では、このボケの影響を抑えて、元画像を作る方法がないか考えるわけですが、その方法が、フィルタ補正逆投影法（フィルタ補正バックプロジェクション）といいます。

フィルタ補正逆投影法

画像の中心領域がボケてしまうということがわかっているので、先にフィルタをかけて、ボケを抑える方法がフィルタ補正逆投影法です。

まず、サイノグラムをなす個々の投影データを1つずつフーリエ変換した後、高い周波数領域になるほど線形に増加するフィルタを乗じて、バックプロジェクションをします。（サイノグラム自体をフーリエ変換してからこの処理をかけてしまうとサンプリング誤差が生じます。）

この線形フィルタをLakフィルタ（またはrampフィルタ）といいます。Lakフィルタは開発者の名前の先頭の文字からこの名前が付けられています。

このフィルタは全く画像ノイズがない理想の状態の画像にかけることを想定されているため、実際には、Shepp-LoganフィルタやHammingフィルタなどがCT画像再構成の標準フィルタとして利用される場合が多いようです。

（Excelによる画像再構成入門より引用）

※こちらの図は逐次近似法の例ですが、フィルタリングをイメージしやすいので載せています。

それでは、既存のCT画像を使って、この操作を辿ってみます。

まず、ImageJのプラグイン「Radon Transform.jar」をダウンロードして、プラグインフォルダに入れて、ImageJを起動します。

（http://rsb.info.nih.gov/ij/plugins/radon-transform.html）

今回はサンプルのCT胸部画像を使っていきます。

Plugins>Radon Transformを起動して、そのままデフォルト設定で、「Calculate」してみます。

すると、投影角度が1度ごとの投影データを行として表したサイノグラムが計算されます。

サイノグラムは先述の投影データを縦軸が角度、横軸をサンプリング位置として上から順に正規化した投影データの濃淡値を並べたデータです。

そして、このサイノグラムを分解して、個々の投影データにフィルタをかけて、元に戻すフィルタ補正逆投影をやってみます。

まずはramp(Lak)フィルタで逆投影で戻してみます。「Use Filtering」のチェックがチェックされていることを確認し、フィルタを選択して、「Reconstruct」を押すだけです。

次は、Shepp-Loganフィルタです。

次はHammingフィルタです。

（それぞれのフィルタ）

最後に、なにもフィルタをかけない、単純なバックプロジェクションの結果です。「Use Filtering」のチェックを外して計算します。

フィルタ補正をして逆投影するおかげで、こんなに画質が改善されているということがわかりますね。

（知っ得）

よく、画像再構成の理論を学ぶ時に目にするデジタルファントムがありますよね。
これです。

Shepp-Loganファントムと呼ばれるものです。
この画像を自分で作りたい方、下記URLが参考になるかもしれません。

http://bigwww.epfl.ch/thevenaz/shepplogan/

次回はCTの再構成法と画質評価法について述べていきます。

Reference

• 「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・3) 2006, p118-119」

第15回 画像処理フィルタ入門（その4）で学ぶ！

今回は、フィルタリングの中でも周波数フィルタリングについて説明していきます。

医療においてもよく使用される基本的なフィルタリング技法ですので、ImageJを使って実体験していきましょう！

周波数フィルタリング

「画像の周波数表現とは何か」「画像を周波数として表現することの利点は何か」など、本ブログアーティクルではImageJを使いながら紹介してきました。

周波数世界でのフィルタリング（周波数フィルタリングまたは空間周波数フィルタリング）は、第11回でその一連の処理過程について数学的な説明も含めて紹介しましたが、ウィナーフィルタを中心とした話でした。

本稿では、周波数フィルタリングの基本に戻り、特定の周波数成分を単純に遮断したり、通過したりするフィルタリングの説明を中心に、ImageJで実践しながら説明を行っていきます。

ローパスフィルタ

画像の周波数成分の低周波数成分だけを残し、高周波成分を除去するようなフィルタをローパスフィルタと称します。

一般的な周波数フィルタリングの定義からすると、周波数領域で表現した元画像をf(u,v)、周波数フィルタをH(u,v)、フィルタリングした出力画像をG(u,v)とすると、周波数フィルタリングの式は、次の式で表されます。

この定義通りに操作を行うには、空間周波数領域での低周波付近（中央部分）が1（つまり、真っ白）、高周波付近が0（真っ黒）であるような周波数フィルタを作れば、上述の式より、元画像とこのローパスフィルタを乗じることで、フィルタリングされた出力画像を得ます。

しかし、もっと簡単な方法として、元画像の画像中央部分である低周波領域だけを抽出して、逆フーリエ変換を行うことで、この出力画像を得ることができます。

では実際にこの低周波領域だけを残す方法を詳細に述べていきます。

まず、任意の画像を用意（今回は腰椎MRI画像）します。

この元画像をProcess>FFT>FFTでフーリエ変換します。

FFTで得られた空間周波数画像の中央に円形ROIを置きます。このとき、ImageJにカーソルの位置が表示されるので、真に中心に持って来たい場合は、注意しながら操作を行います。

ImageJに表示されている情報は、カーソルの位置がゼロ周波数からどれくらいの周波数分離れているのかを示す（mm/c）や、周波数平面の方向成分である角度情報およびその強度が表示されます。この（mm/c）は、例えばFOVが345mmで512×512マトリックスの場合、1ピクセルの大きさは、345/512=0.673mmで、1周期あたりの最小の長さがピクセル2つ分なので、約1.35mm/cycleと表示されるはずです。

次に、低周波成分を残すための円形のROIの外側を削除し、内側だけを残すには、ROIを作成した後に、Edit>Clear outsideで操作します。

そうすると、ROI以外の領域が白（濃度値255）に置き換わるので、これを反転します。

反転をするには、まず、ROI以外の領域を選択する必要があるので、Edit>Selection>Make inverseで、ROI以外の領域を選択してから、Edit>Invertで黒く（濃度値0）反転します。

これで、次の図のような空間周波数画像が得られます。

そして、この画像をProcess>FFT>Inverse FFTで変換しましょう。

次のようなローパスフィルタ処理のかかった画像が得られます。

（補足）

まれに、周波数座標系の画像を空間座標系の画像として表現したフィルタを利用する場合がありますが、こういった場合には、Process>FFT>Custom Filter>FFT filterでこのようなフィルタ画像を選択することで、元画像にフィルタリングできます。

ハイパスフィルタ

ローパスフィルタと逆の周波数成分濾過フィルタをハイパスフィルタと言います。

ハイパスフィルタの例を示します。

（ローパスフィルタ作成同様に、中央にROIをおき、Edit>Cut）

（上記画像を逆フーリエ変換）

ここで、ローパス、ハイパスを見比べると、ローパスフィルタの方がハイパスフィルタよりも画像の情報が多いことがわかります。

これは、画像を観察する場合の周波数成分のほとんどが低周波成分であることを意味しています。

この性質は巧く応用され、画像圧縮などに利用されています。

特定周波数領域フィルタ

これらのローパスフィルタやハイパスフィルタは、ある特定の領域の周波数だけを通すあるいは遮断するフィルタです。このような特定周波数領域を通すようなフィルタはバンドパスフィルタと呼ばれています。特定周波数なので、例えば、あるハイパスフィルタAとハイパスフィルタBを差し引いたドーナツ型の特定周波数領域をフィルタリングすることもできますね。

このような処理を可能にするImageJの機能に、Process>FFT>Bandpass Filterがあります。

このBandpass Filter機能は、低周波成分のLarge Structureと高周波成分のSmall Structureのピクセルサイズを決めて、ガウシアンフィルタを施して、それぞれをサブトラクションしてフィルタを作成し、元画像に適用します。縞用のアーチファクトを抑制するSuppress Stripes機能も使えます。

実際に、このバンドパスフィルタをデフォルト設定値で元画像に適用し、元画像のFFTとフィルタリング画像のFFTとの差分をとると、どのくらいの周波数情報が遮断されたのかがわかります。

（バンドパスフィルタで遮断された周波数）

次回は、少し内容を変えて、断層画像の投影再構成を説明します！

Reference

• 「山本修司：ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・2) 2006, p123-125」