JP2004535872A

JP2004535872A - Ｃｔ画像形成におけるストリークアーチファクトおよび照明不均一性を低減または除去するシステムおよび方法

Info

Publication number: JP2004535872A
Application number: JP2003514485A
Authority: JP
Inventors: タメツ−ペナジョセ; マージャッタソフィアトッターマンサーラ; パーカーケビン
Original assignee: ヴァーチャルスコピックスリミテッドライアビリティカンパニー
Priority date: 2001-07-19
Filing date: 2002-07-18
Publication date: 2004-12-02
Also published as: US6801646B1; US7406211B2; CA2452528A1; WO2003009215A1; EP1423819A1; US20050094857A1

Abstract

コンピュータ断層撮影画像（１０２）において、金属インプラントの存在によって生じるストリークアーチファクトおよび組織の密度が不均一に推定されることを低減または解消する。このアルゴリズムは、１）照明補正（１０４）と、２）適応３次元フィルタリング（１０６）の２つの基本ステップを有する。このアルゴリズムは、形態的なグレースケール拡張によって、ストリークの方向および不均一密度の程度を推定することから始める。引き続き、これらの推定値に基づいて正しい密度を推定し、ストリークの方向および局所的な画像コントラストによって決まるパラメータをもつ適応３次元フィルタリング（１０６）によってストリークを低減する。

Description

【技術分野】
【０００１】
本発明は、ＣＴ（コンピュータ断層撮影法）などによって得られた画像を処理するシステムおよび方法を対象とし、より詳細には、金属インプラントの存在によって生じるアーチファクトを低減または除去するこのようなシステムおよび方法を対象とする。
【背景技術】
【０００２】
Ｘ線ＣＴ（コンピュータ断層撮影法）は、特に医学および歯科の応用分野では重要な画像形成技術である。多数の異なる角度からの一連のＸ線ビームを使用して、患者の体の断面像を生成する。コンピュータ内でこれら複数のスライス画像を合成して、器官、骨格および組織を極めて詳細に表示することができる３次元像にする。ＣＴにより、高空間分解能が得られ、３次元的に重ね合わせることができ、動きによって生じるぼけが最小限に抑えられる。
【０００３】
しかし、金属インプラントや充填物など減衰が強く生じる対象物が存在すると、画像内に、放射状アーチファクトとも呼ばれるストリークアーチファクト（ｓｔｒｅａｋａｒｔｉｆａｃｔ）が生じる。ＣＴで生じる別の問題は、不均一な照明によって組織の密度が不均一に推定されることである。これらの問題を低減または解消するために多くの研究がなされてきたが、未だ満足な技術は見つかっていない。
【０００４】
金属インプラントによって生じるアーチファクトの問題を解決しようとする１つの試みがある（例えば、Ｈｕの特許文献１参照）。この特許は、ヘリカル再構築からのデータをバックグラウンド成分と鮮明構造成分に分離するＣＴ画質改善方法を教示している。鮮明な構造と画像のバックグラウンドは、通常、ＣＴ数が大きく異なるので、グレースケール閾値処理を用いて分離することができる。画像のバックグラウンドをフィルタリング処理して高周波アーチファクトを除去し、これらの画像を再び合成する。しかし、Ｈｕの技法では、望まれるものよりも多くの計算量が必要であり、Ｈｕの手法は、放射線医学分野では広く採用されていないようである。
【０００５】
【特許文献１】
米国特許第５，５６１，６９５号明細書
【発明の開示】
【発明が解決しようとする課題】
【０００６】
当技術分野では、ストリークアーチファクトおよび照明の不均一性を低減または除去する計算上効率的な技術が求められていることは上記から容易に明らかであろう。
【課題を解決するための手段】
【０００７】
したがって、本発明の目的は、金属インプラントなどによって生じるストリークアーチファクトを低減または除去することである。
【０００８】
本発明の別の目的は、照明の不均一性を補正することである。
【０００９】
本発明の別の目的は、上記目的を計算上効率的な方法で実現することである。
【００１０】
上記その他の目的を実現するために、本発明は、再構築データを直接取り扱うＣＴその他の画像形成技術を対象とする。この技術は、１）照明補正と、２）適応３Ｄフィルタリング（ａｄａｐｔｉｖｅ３Ｄｆｉｌｔｅｒｉｎｇ）の２つの基本ステップを有する。このアルゴリズムは、形態的なグレースケール拡張（ｇｒａｙｓｃａｌｅｍｏｒｐｈｏｌｏｇｙｄｉｌａｔｉｏｎ）によって、ストリークの方向および不均一密度の程度を推定することから始める。引き続き、この推定値に基づいて正しい密度を推定し、ストリークの方向および局所的な画像コントラストによって決まるパラメータをもつ適応３Ｄフィルタリングによってストリークを低減する。
【発明を実施するための最良の形態】
【００１１】
次に、図面を参照して、本発明の好ましい実施形態を詳細に述べる。
【００１２】
図１に、好ましい実施形態で実施するプロセスの概要を示す。このプロセスは、ステップ１０２で生データを取得した後で、照明補正ステップ１０４および適応３Ｄフィルタリングステップ１０６の２つのステップを含む。次に、これら２つのステップを詳細に説明する。
【００１３】
このアルゴリズムでは、金属インプラントを含む１組の再構築ＣＴ画像ｆ（ｘ，ｙ）から、それに続いて、３Ｄ（３次元）ボリューム画像ｆ（ｘ，ｙ，ｚ）内でこれらの画像をスタックする。次いで、このアルゴリズムでは、すべての単一スライスをグレースケール拡張して、ストリークアーチファクトの伝播ポテンシャル場（ｐｒｏｐａｇａｔｉｏｎｐｏｔｅｎｔｉａｌｆｉｅｌｄ）ｐ（ｘ，ｙ，ｚ）を推定し、かつ不均一照明Ｉ（ｘ，ｙ，ｚ）を推定する。グレースケール拡張によるポテンシャル場ｐ（ｘ，ｙ，ｚ）の推定は、図２に示す以下のステップによって実施する。
【００１４】
ステップ２０２：生データの複製を生成する。ｒ（ｘ，ｙ，ｚ）＝ｆ（ｘ，ｙ，ｚ）
ステップ２０４：ｐ（ｘ，ｙ，ｚ）＝０となるように初期化を行い、カウンタｃを１にセットする。
【００１５】
ステップ２０６：すべての再構築ＣＴスライス画像ｆ（ｘ，ｙ）を次のようにグレースケール拡張する。
【００１６】
【数１】

ただし、ｈ（ｘ，ｙ）はカーネルであり、ｃが奇数の場合は以下のように与えられ、
【００１７】
【数２】

ｃが偶数の場合は以下のように与えられる。
【００１８】
【数３】

この拡張処理により、強度が高い領域から拡大する前線（ｇｒｏｗｉｎｇｆｒｏｎｔ）が確立される。
【００１９】
ステップ２０８：ｐ（ｘ，ｙ，ｚ）＝ｃとし、ｇ（ｘ，ｙ，ｚ）＞ｒ（ｘ，ｙ，ｚ）を満たすすべての点について拡大経路（ｇｒｏｗｉｎｇｐａｔｈ）に沿った照明関数（ｉｌｌｕｍｉｎａｔｉｏｎｄｅｒｉｖａｔｉｖｅ）Ｉ（ｘ，ｙ，ｚ）の導関数を推定する。この方向導関数（ｄｉｒｅｃｔｉｏｎａｌｄｅｒｉｖａｔｉｖｅ）から照明関数を再構築する。このステップは、後で図３を参照して説明する。
【００２０】
ステップ２１０：ｒ（ｘ，ｙ，ｚ）＝ｇ（ｘ，ｙ，ｚ）とする。
【００２１】
ステップ２１２：ｃを次のように増やす。ｃ＝ｃ＋１
ステップ２１４：ｐ（ｘ，ｙ，ｚ）内の点の値が変化したかどうかを判定する。変化があった場合には、ステップ２０６〜２１２を繰り返す。そうでない場合には、このプロセスをステップ２１６で終了する。
【００２２】
次に、ステップ２０８に関して上記で述べた照明関数の再構築を説明する。体内の金属により、画像の照明プロフィールになんらかの変化が生じることになる。画像は、それぞれ均一な密度の小片状要素からなると考える。したがって、任意の組織内部の画像密度に大きな変動がないはずであり、照明関数は、推定照明導関数（ｅｓｔｉｍａｔｉｏｎｏｆｔｈｅｉｌｌｕｍｉｎａｔｉｏｎｄｅｒｉｖａｔｉｖｅ）
【００２３】
【数４】

【００２４】
を伝播する前線の方向に積分することによって再構築することができる。
【００２５】
【数５】

【００２６】
の推定は、近傍差分の適応重みつき差分（ａｄａｐｔｉｖｅｗｅｉｇｈｔｅｄｄｉｆｆｅｒｅｎｃｅ）と、近傍の前推定
【００２７】
【数６】

【００２８】
によって行われる。
【００２９】
伝播する前線を推定するステップ２０２で、
【００３０】
【数７】

【００３１】
を初期化してゼロにする場合、ステップ２０８では、以下の式を用いて方向導関数を推定する。
【００３２】
【数８】

ただし、
【００３３】
【数９】

であり、ｕは
【００３４】
【数１０】

【００３５】
となる正規化定数であり、
【００３６】
【数１１】

である。上式により、ｐ（ｘ，ｙ，ｚ）＝ｃであるすべての点で、推定照明方向導関数が得られる。
【００３７】
数量σは不連続性の閾値であり、画像化する体によって決まる。人体の研究に適した閾値はσ＝０．１４である。一方、定数αは画像分解能によって決まる。０．５×０．５ｍｍ^２のボクセル（ｖｏｘｅｌ）分解能に適した値はα＝０．０３である。上記から、重み関数ω（ｉ，ｊ，ｃ）は拡大前線によって決まり、γ（ｉ，ｊ）は画像の不連続性によって決まることが理解されよう。
【００３８】
照明関数の導関数を推定した後で、それを積分して照明関数を求める。この積分は、Ｉ（ｘ，ｙ，ｚ）＝ｋと設定することによって開始する。ただし、ｋは、平均画像強度に等しい定数である。Ｉ（ｘ，ｙ，ｚ）の初期化後、この積分は、画像の外側の境界から始めて、ポテンシャル画像がゼロになる点まで行う。このプロセスは、図３に示す以下のステップによって実施する。
【００３９】
ステップ３０２：Ｉ（ｘ，ｙ，ｚ）＝ｋとなるように初期化を行い、ｃ＝ｍａｘ（ｐ（ｘ，ｙ，ｚ））となるようにカウンタｃを最大ポテンシャルにセットする。
【００４０】
ステップ３０４：ｐ（ｘ，ｙ，ｚ）＝ｃであるすべての点の近傍に注目し、それらの点で、次式を行う。
【００４１】
【数１２】

ただし、
【００４２】
【数１３】

であり、ｕは、ｕΣΣβ（ｉ，ｊ）＝１となる正規化定数である。
【００４３】
ステップ３０６：ｃを次のように減らす。ｃ＝ｃ−１
ステップ３０８：ｃ＜１かどうかを判定する。ｃ＜１になるまでステップ３０４と３０６を繰り返し、ｃ＜１になったとき、ステップ３１０でこのプロセスを終了する。
【００４４】
照明場が得られた後で、ｆ（ｘ，ｙ，ｚ）における照明の不均一性を次のように補正する。
【００４５】
【数１４】

この新たに補正された画像ｇ（ｘ，ｙ，ｚ）では、特に密度が一定であるべき点で、元の画像よりも画像強度の変動が少なくなる。
【００４６】
現実には、元の画像は、不連続性の検出が極めて困難となり得るほど劣化していた。この問題を回避するために、ストリークアーチファクトの一部の周りのノイズを減らす非線形な構造保存フィルタ（ｎｏｎ−ｌｉｎｅａｒｓｔｒｕｃｔｕｒｅｐｒｅｓｅｒｖｉｎｇｆｉｌｔｅｒ）によって、画像ｆ（ｘ，ｙ，ｚ）をあらかじめフィルタリングする。照明プロフィールの推定には、このフィルタリングした画像を使用した。
【００４７】
上記で説明したように画像を補正した後で、ストリークアーチファクトを低減することができる。３×３×３の適応フィルタＬ（ｘ，ｙ，ｚ）を用いてストリークアーチファクトを低減する。このフィルタの係数は、伝播ポテンシャルおよび局所的な画像コントラストの関数であり、次式で与えられる。
【００４８】
【数１５】

【００４９】
ただし、ｓは、点（ｘ_０，ｙ_０，ｚ_０）におけるノイズの推定変動値であり、ηは正規化定数である。
【００５０】
異なるポテンシャルをもつ隣接点ならびに密度が極めて異なる隣接点における係数は、非常に小さくなる。一方、ノイズの推定変動値とポテンシャルが類似しており、かつ密度が極めて類似している隣接点における係数は、より大きくなる。このように定式化すると、アーチファクトが少なくなり、組織間で強いコントラストをもつ点が保存される。
【００５１】
点（ｘ_０，ｙ_０，ｚ_０）における最終処理画像ｈ（ｘ_０，ｙ_０，ｚ_０）は、次式で与えられる。
【００５２】
【数１６】

この適応フィルタリング処理を数回行って、画像中に存在する大部分のストリークを効果的に除去することができる。
【００５３】
面内分解能がスライスの厚さよりもはるかに小さい場合、画像を劣化させないようにこの適応フィルタを改変することができる。スペース／ストリークを対象としたノイズ推定値（ｓｐａｃｅａｎｄｓｔｒｅａｋｏｒｉｅｎｔｅｄｎｏｉｓｅｅｓｔｉｍａｔｉｏｎ）を与えることによって、画質を改善することもできる。このようにして、ストリークアーチファクトがより大きな領域にはフィルタリング処理を行い、ストリークがそれほど強くない領域ではフィルタリング処理を避けることができる。
【００５４】
腰部人工器官（ｈｉｐｐｒｏｓｔｈｅｓｉｓ）を含むいくつかのＣＴ画像に対してこのアルゴリズムをテストした。図４Ａに、これらの画像の１枚の単一スライスを示す。この画像スライスは、大腿骨基部および人工カップ（ｐｒｏｓｔｈｅｓｉｓｃｕｐ）を含む。図４Ｂに、推定ポテンシャル場を示し、図４Ｃに、推定照明場を示す。図４Ｄは、照明アーチファクト除去後のＣＴ画像である。図４Ｅに、適応フィルタによってアーチファクトを平滑化した後のフィルタ処理画像を示す。図５Ａおよび５Ｂにそれぞれ、生データおよび処理データからの腰部ボリュームレンダリング（ｖｏｌｕｍｅｔｒｉｃｒｅｎｄｅｒｉｎｇ）を示す。処理データからの腰部ボリュームレンダリングでは、金属人工器官ならびにそれを取り囲む骨格のより細部が見えることがわかる。
【００５５】
上記で開示した実施形態および他の実施形態を、図６のブロック図で示すようなシステムで実施することができる。システム６００は、画像データなどを入力する入力装置６０２を含む。入力装置６０２から入力された情報は、ハードドライブなどの記憶装置６０６と、上記で開示した処理を実施する処理ユニット６０８と、例えばサーフェスレンダリング（ｓｕｒｆａｃｅｒｅｎｄｅｒｉｎｇ）によって観察用の最終処理画像を用意（ｐｒｅｐａｒｅ：作成）するグラフィックレンダリングエンジン６１０とを有するワークステーション６０４で受け取られる。出力装置６１２は、レンダリングエンジン６１０によってレンダリングされた画像を観察するためのモニタ、または画像を記録するビデオレコーダなど別の記憶装置、あるいはそれら両方を含み得る。ワークステーション６０４およびグラフィックレンダリングエンジン６１０の例は、ＳｉｌｉｃｏｎＧｒａｐｈｉｃｓ社のＩｎｄｉｇｏワークステーションおよびＩｒｉｘＥｘｐｌｏｒｅｒ３Ｄグラフィックエンジンである。
【００５６】
上記で好ましい実施形態を詳細に述べてきたが、本開示を検討した当業者には、本発明の範囲内で他の実施形態を実現できることが容易に理解されよう。例えば、αおよびσについて与えられた数値は、限定的なものではなく例である。また、画像化する対象物は、人間または動物の体の任意の部分その他の任意の対象物とすることができる。さらに、本発明は、ＣＴ以外の画像形成技術にも適用可能である。したがって、本発明は、添付の特許請求の範囲によってのみ限定されると解釈すべきである。
【図面の簡単な説明】
【００５７】
【図１】好ましい実施形態によるプロセスを示す流れ図である。
【図２】好ましい実施形態によるプロセスを示す流れ図である。
【図３】好ましい実施形態によるプロセスを示す流れ図である。
【図４Ａ】画像スライスの処理におけるステップを示す写真を表す図である。
【図４Ｂ】画像スライスの処理におけるステップを示す写真を表す図である。
【図４Ｃ】画像スライスの処理におけるステップを示す写真を表す図である。
【図４Ｄ】画像スライスの処理におけるステップを示す写真を表す図である。
【図４Ｅ】画像スライスの処理におけるステップを示す写真を表す図である。
【図５Ａ】生データからの画像のボリュームレンダリングを示す写真を表す図である。
【図５Ｂ】処理データからの画像のボリュームレンダリングを示す写真を表す図である。
【図６】好ましい実施形態を実施することができるシステムを示すブロック図である。

Claims

画像中のアーチファクトを減らし、照明の均一性を改善しながら対象物を画像化する方法であって、
（ａ）前記対象物から生画像データを受け取るステップと、
（ｂ）前記生画像データから照明関数を推定するステップと、
（ｃ）前記推定照明関数に従って前記生画像データを補正して照明補正画像を得るステップと、
（ｄ）前記照明補正画像から適応フィルタを作成するステップと、
（ｅ）前記適応フィルタおよび前記照明補正画像から最終処理画像を作成するステップとを含むことを特徴とする方法。
前記生画像データを受け取る前記ステップを、コンピュータ断層撮影法によって実施することを特徴とする請求項１に記載の方法。
前記対象物は体の部分であることを特徴とする請求項１に記載の方法。
前記照明関数を推定する前記ステップは、
（ｉ）前記生画像データの複製を生成するステップと、
（ｉｉ）ポテンシャル場をゼロに初期化するステップと、
（ｉｉｉ）第１カウンタを１に初期化するステップと、
（ｉｖ）前記カウンタ値に従って前記生画像データのグレースケール拡張を行って、グレースケール拡張データを得るステップと、
（ｖ）前記ポテンシャル場を前記第１カウンタに等しく設定するステップと、
（ｖｉ）前記照明関数の方向導関数を推定するステップと、
（ｖｉｉ）前記方向導関数から前記照明関数を再構築するステップと、
（ｖｉｉｉ）前記生画像データを前記グレースケール拡張データに等しく設定するステップと、
（ｉｘ）前記第１カウンタを増やすステップと、
（ｘ）前記ポテンシャル場内の点の値が変化しなくなるまで、ステップ（ｉｖ）から（ｉｘ）を繰り返すステップとを含むことを特徴とする請求項１に記載の方法。
各点における前記方向導関数を、その点の生画像データおよび隣接する点の生画像データから推定することを特徴とする請求項４に記載の方法。
前記照明関数は、
（Ａ）平均画像強度に等しい定数に前記照明関数を初期化するステップと、
（Ｂ）前記ポテンシャル場の最大値に第２カウンタを初期化するステップと、
（Ｃ）前記強度関数および前記強度関数の方向導関数から前記強度関数を計算し直すステップと、
（Ｄ）前記第２カウンタを減らすステップと、
（Ｅ）前記第２カウンタがカウンタ閾値未満になるまで、ステップ（Ｃ）および（Ｄ）を繰り返すステップとによって再構築されることを特徴とする請求項５に記載の方法。
前記カウンタ閾値は１であることを特徴とする請求項６に記載の方法。
各点における前記適応フィルタは、前記ポテンシャル場および前記照明補正画像の関数であることを特徴とする請求項４に記載の方法。
各点における前記最終処理画像は、その点および隣接する点の前記適応フィルタおよび前記照明補正画像の関数であることを特徴とする請求項８に記載の方法。
各点における前記最終処理画像は、その点および隣接する点の前記適応フィルタおよび前記照明補正画像の関数であることを特徴とする請求項１に記載の方法。
ステップ（ｅ）は、前記照明補正画像に前記適応フィルタを複数回適用するステップを含むことを特徴とする請求項１に記載の方法。
画像中のアーチファクトを減らし、照明の均一性を改善しながら対象物を画像化するシステムであって、
生画像データを受け取る入力部と、
前記生画像データから照明関数を推定し、前記推定照明関数に従って前記生画像データを補正して照明補正画像を得、前記照明補正画像から適応フィルタを作成し、前記適応フィルタおよび前記照明補正画像から最終処理画像を作成するプロセッサと、
前記最終処理画像を出力する出力部とを備えることを特徴とするシステム。
前記プロセッサは、
（ｉ）前記生画像データの複製を生成するステップと、
（ｉｉ）ポテンシャル場をゼロに初期化するステップと、
（ｉｉｉ）第１カウンタを１に初期化するステップと、
（ｉｖ）前記カウンタ値に従って前記生画像データのグレースケール拡張を行って、グレースケール拡張データを得るステップと、
（ｖ）前記ポテンシャル場を前記第１カウンタに等しく設定するステップと、
（ｖｉ）前記照明関数の方向導関数を推定するステップと、
（ｖｉｉ）前記方向導関数から前記照明関数を再構築するステップと、
（ｖｉｉｉ）前記生画像データを前記グレースケール拡張データに等しく設定するステップと、
（ｉｘ）前記第１カウンタを増やすステップと、
（ｘ）前記ポテンシャル場内の点の値が変化しなくなるまで、ステップ（ｉｖ）から（ｉｘ）を繰り返すステップとによって、前記照明関数を推定することを特徴とする請求項１２に記載のシステム。
各点における前記方向導関数を、その点の生画像データおよび隣接する点の生画像データから推定することを特徴とする請求項１３に記載のシステム。
前記照明関数は、
（Ａ）平均画像強度に等しい定数に前記照明関数を初期化するステップと、
（Ｂ）前記ポテンシャル場の最大値に第２カウンタを初期化するステップと、
（Ｃ）前記強度関数および前記強度関数の方向導関数から前記強度関数を計算し直すステップと、
（Ｄ）前記第２カウンタを減らすステップと、
（Ｅ）前記第２カウンタがカウンタ閾値未満になるまで、ステップ（Ｃ）および（Ｄ）を繰り返すステップとによって再構築されることを特徴とする請求項１４に記載のシステム。
前記カウンタ閾値は１であることを特徴とする請求項１５に記載のシステム。
各点における前記適応フィルタは、前記ポテンシャル場および前記照明補正画像の関数であることを特徴とする請求項１３に記載のシステム。
各点における前記最終処理画像は、その点および隣接する点の前記適応フィルタおよび前記照明補正画像の関数であることを特徴とする請求項１７に記載のシステム。
各点における前記最終処理画像は、その点および隣接する点の前記適応フィルタおよび前記照明補正画像の関数であることを特徴とする請求項１２に記載のシステム。
前記プロセッサは、前記照明補正画像に前記適応フィルタを複数回適用することによって、前記最終処理画像を作成することを特徴とする請求項１２に記載のシステム。