JP6173694B2 - 画像処理装置、ｘ線コンピュータ断層撮影装置および画像処理方法 - Google Patents

Description

本発明の実施形態は、画像処理装置、Ｘ線コンピュータ断層撮影装置および画像処理方法に関する。

コンピュータ断層撮影（ＣＴ）における画像の再構成のため、逐次アルゴリズムは、全変動（ＴＶ）最小化を含むことがある。たとえば、オーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）と同時代数再構成法（ＳＡＲＴ）とに基づく方法がいくつか存在している。ＴＶ最小化は、ノイズを除去するために、フィルタ補正逆投影法（ＦＢＰ）などの解析的プログラムを用いて再構成される画像に適用されることもある。

従来技術における努力にもかかわらず、いくつかの問題が未解決のまま残っており、改善が求められている。従来技術においては、エッジを除き、画像におけるノイズが極めて一様に除去されるのが一般的である。換言すると、ノイズは、全変動（ＴＶ）最小化に基づいて除去されるが、元の画像が非一様な分布を有する場合には、結果的に得られる画像は、ある領域では過渡に平滑化され、それ以外の領域では依然としてノイズを含むままである。他方で、微細なエッジを有する画像に従来技術が適用されると、それらのエッジは平滑化によって消滅してしまう場合がある。臨床的応用において日常的である非一様なノイズ分布を元の画像が有している場合、従来技術による全変動（ＴＶ）最小化では、空間的解像度を犠牲にすることなくノイズを改善することは不可能であるように思われる。

更に、ＨＵ値が平均値から２σから３σずれていて真の値から遠く離れているようなノイズを含む画像では、雑多な画素クラスタがいくつか存在する。ＴＶ最小化を実行すると、これらのクラスタは、ノイズが除去された画像において黒／白のドットとして残ってしまう。その理由は、（１）「クラスタ効果」がノイズ除去を低速化するからであり、また、（２）元のＨＵ値が真の値から遠く離れているからである。

目的は、低ノイズ領域における過剰な補償を回避し、高ノイズの画素クラスタにおけるドットを減少させるような全変動（ＴＶ）最小化アルゴリズムを実装するための実際的な解決策を提供することにある。

本実施形態に係る画像のノイズ除去装置は、
ａ）前記画像の各画素に対応付けられる方向インデックスを初期化する手段であって、前記方向インデックスは、全変動の離散勾配に関する増加と減少とのいずれかを示す、前記初期化する手段と、
ｂ）前記各画素について、前記全変動の離散勾配を決定する手段であって、前記離散勾配は前記増加又は前記減少を示す、前記決定する手段と、
ｃ）前記方向インデックスと前記離散勾配とが等しい場合には、前記各画素の画素値を所定値だけ変更して前記全変動を減少させる手段と、
ｄ）所定の条件に到達するまで、前記手段ｂ）から前記手段c）までを反復して実行させる手段と、を備える。

図１は本発明によるマルチスライスＸ線ＣＴ装置またはスキャナの一実施形態を示す図である。 図２は本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスに含まれるステップを示すフローチャートである。 図３は本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられるオーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）に含まれるステップを示すフローチャートである。 図４は本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法に含まれるステップを示すフローチャートである。 図５Ａは従来技術による全変動最小化の従来型のプロセスにおいて処理される単一のオッド画素の場合を示す図である。 図５Ｂは従来技術による全変動最小化のプロセスにおいて処理されるオッド画素クラスタの場合を示す図である。 図６は本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法によって処理される例示的な画素の場合を示す図である。 図７Ａは本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法によって処理される例を示す図である。 図７Ｂは本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法によって処理される例を示す図である。 図７Ｃは本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法によって処理される例を示す図である。 図７Ｄは本発明による全変動逐次近似法（ＴＶ−ＩＲ）の好適なプロセスにおいて用いられる離散ＴＶ最小化法によって処理される例を示す図である。 図８Ａは投影データから再構成された画像であって、その画像中に高レベルのノイズが非一様な態様で分布している画像を示す図である。 図８Ｂは全変動逐次近似法（ＴＶ−ＩＲ）を用いて再構成された画像を示す図である。 図８Ｃは本発明による一方向条件を伴う全変動逐次近似法（ＴＶ−ＩＲ）を用いて再構成された例示的な画像を示す図である。 図９Ａはエッジを横断する一連の画素とＣＴ数におけるそれらの画素値とを、元の画像と、離散全変動（ＤＴＶ）最小化の後と、全変動（ＴＶ）最小化の後とについて示すグラフである。 図９Ｂは図９Ａにおける青い矩形のボックスの中の画素の拡大されたエッジ部分を示すプロットである。

以下、本発明の一実施形態について図面を参照して説明する。同じ符号を付したものは同じ構成要素で有ることとする。図１を参照すると、本実施形態によるマルチスライスＸ線ＣＴ装置またはスキャナの構成が示されており、ガントリ１００とそれ以外のデバイスまたはユニットとが含まれる。ガントリ１００は、正面から見た様子が示されていて、Ｘ線管１０１と、環状フレーム１０２と、マルチローまたは２次元配列タイプのＸ線検出器１０３とを更に含む。Ｘ線管１０１とＸ線検出器１０３とは、被験者Ｓを横断する直径方向に、軸ＲＡを中心として回転する環状フレーム１０２の上に設置されている。回転ユニット１０７は、被験者Ｓが軸ＲＡに沿って図面の奥の方向にまたは図面から出てくる方向に移動する間、０．２７秒／回転程度の高速度でフレーム１０２を回転させる。

マルチスライスＸ線ＣＴ装置は、更に、管電圧をＸ線管１０１に与えてＸ線管１０１にＸ線を発生させる高電圧発生器１０９を含む。ある実施形態では、高電圧発生器１０９は、フレーム１０２に設置される。Ｘ線は、被験者Ｓに向かって放射されるが、被験者Ｓの断面積は円によって表されている。Ｘ線検出器１０３は、被験者Ｓを横断してＸ線管１０１の反対側に配置され、被験者Ｓを通過して送信されるように放射されるＸ線を検出する。

更に図１を参照すると、Ｘ線ＣＴ装置またはスキャナは、更に、Ｘ線検出器１０３からの検出された信号を処理する他のデバイスを含む。データ取得回路またはデータ取得システム（ＤＡＳ）１０４は、それぞれのチャネルに対してＸ線検出器１０３から出力される信号を、電圧信号に変換し、増幅し、更にデジタル信号に変換する。Ｘ線検出器１０３とＤＡＳ１０４とは、１回転当たり所定の総数の投影（ＴＰＰＲ）を処理するように構成されている。

上述したデータは、非接触データ送信機１０５を介して、ガントリ１００の外部のコンソールの中に設置されている前処理デバイス１０６に送られる。前処理デバイス１０６は、生データに対して感度訂正など特定の訂正を実行する。次に、ストレージデバイス１１２が、結果として得られるデータを記憶する。このデータは、再構成処理の直前の段階においては投影データとも呼ばれる。ストレージデバイス１１２は、データ／コントロールバスを介してシステムコントローラ１１０に接続されるが、再構成デバイス１１４、離散全変動（ＴＶ）最小化デバイス１１７、表示デバイス１１６、入力デバイス１１５、そしてスキャンプランサポート装置２００も同様である。スキャンプランサポート装置２００は、画像担当技術者をサポートしてスキャンプランを作成する機能を含む。

再構成デバイス１１４と離散全変動（ＴＶ）の最小化デバイス１１７との一実施形態は、様々なソフトウェアコンポーネントとハードウェアコンポーネントとを更に含む。本発明の一態様によると、ＣＴ装置の再構成デバイス１１４は、離散ＴＶ最小化デバイス１１７が全変動を効果的に最小化する間に、画像を再構成する。一般に、本発明の一実施形態における再構成デバイス１１４と離散ＴＶ最小化デバイス１１７とは、画像を再構成する全変動逐次近似（ＴＶＩＲ）アルゴリズムを実行する。更に詳しくは、再構成デバイス１１４は、離散ＴＶ最小化デバイス１１７がＴＶ最小化プロセスを実行する間に、オーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）によるプロセスを、投影データに対して実行する。再構成デバイス１１４と離散ＴＶ最小化デバイス１１７との一実施形態は、それぞれ、第１の所定回数の反復が定められているメインループと第２の所定回数の反復が定められている内部ループとにおいて実装される２つのプロセスを実行する。別の実施形態では、離散ＴＶ最小化デバイス１１７は、予め定められた特定の条件が達成されるまで、内部ループを反復する。

ＴＶ最小化プロセスの前に、再構成デバイス１１４は、オーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）を、投影データに対して実行する。投影データは、それぞれが特定の数のビューを有する所定の数の部分集合（サブセット）Ｎにグループ分けされる。ＯＳＳＡＲＴの間に、それぞれの部分集合は、一実施形態では、シーケンシャルに処理される。別の実施形態では、これらの複数の部分集合を、複数の中央演算装置（ＣＰＵ）やグラフィクス演算装置（ＧＰＵ）などの特定のマイクロプロセッサを利用することにより、パラレルに処理することも可能である。更に別の実施形態では、再構成デバイス１１４は、緩和パラメータを用いて、オーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）を投影データに対して実行する。

オーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）の間に、再構成デバイス１１４は、２つの重要な動作を更に実行する。すなわち、それぞれの部分集合Ｎに対し、再構成デバイス１１４は、画像ボリュームを再投影して、計算された投影データを形成し、現在の画像評価に加えられると、測定された投影データと計算された投影データとの間の正規化された差を逆投影して、更新された画像ボリュームを再構成する。更に詳細には、再構成デバイス１１４の一実施形態では、システム行列の係数がどれもキャッシュされないレイトレーシング（ray tracing）法を用いることにより、画像ボリュームを再投影する。更に、再構成デバイス１１４の一実施形態では、ある部分集合の中のすべてのレイを同時に再投影するのであるが、これは、オプションであるがパラレルに実装される。他の実施形態には、システムモデルやノイズモデルが含まれる。逆投影においては、再構成デバイス１１４の一実施形態では、画素駆動法を用いて、ある部分集合の中の正規化された差投影データのすべてを逆投影して、所望の更新された画像ボリュームを形成する。再構成デバイス１１４はある部分集合の中のすべてのレイの和すなわち差投影データを逆投影して、画像ボリュームを形成するので、この動作もまた、オプションであるがパラレルに実装される。これらの動作が、すべての部分集合Ｎに対して適用されることにより、ＯＳＳＡＲＴのステップが１つ完了する。この実施形態とそれ以外の実施形態とが、オプションであるが、特許請求の範囲に記載された本発明の範囲に含まれる。

全変動（ＴＶ）最小化プロセスでは、離散ＴＶ最小化デバイス１１７は、所与の画素値を調整するため、所定のステップパラメータとステップ数とを用いて、離散ＴＶ最小化を実行する。オプションであるが、一方向条件が適用される。図６を参照して後述するように、所定の一方向条件が、特定の低ノイズ領域が過渡に補償されるまたは過剰に平滑化されることを回避しつつ、高ノイズ領域におけるノイズを減少させる。離散ＴＶ最小化では、画素がエッジ上に存在すると、その画素値はノイズ除去の前後を通じて厳密に同一のままに維持される。換言すると、離散ＴＶ最小化デバイス１１７の一実施形態では、画像の中のぼけたエッジ（blur edges）を実質的に減少させる。更に、離散ＴＶ最小化デバイス１１７の一実施形態では、特に、強力なノイズ除去の後でのノイズ除去された不自然なノイズテクスチャが実質的に改善される。離散ＴＶ最小化デバイス１１７の一実施形態では、所定の条件に到達すると、離散ＴＶ最小化を終了する。例示的な終了条件として、一方向条件に起因するゼロ勾配があり、他の例としては、離散ＴＶ最小化が所定の回数だけ反復されることがある。

本発明による上述したマルチスライスＸ線ＣＴ装置またはスキャナでは、再構成デバイス１１４の一実施形態が、逐次近似法に基づいて画像を再構成する。再構成デバイス１１４の別の実施形態では、フィルタ補正逆投影（ＦＢＰ）などの解析的な再構成法に基づいて、画像を再構成する。更に、本発明によるマルチスライスＸ線ＣＴ装置またはスキャナの上述した実施形態には、再構成デバイス１１４と離散ＴＶ最小化デバイス１１７とのために別々のデバイスが含まれている。他の実施形態では、再構成デバイス１１４と離散ＴＶ最小化デバイス１１７とは、オプションであるが、ソフトウェアとハードウェアとの組み合わせを含む単一のデバイスである。本発明は、どのような特定のハードウェアとソフトウェアとの組み合わせにも限定されることなく、一方向条件を伴う離散ＴＶ最小化に基づいてノイズ除去の後のノイズテクスチャを改善するように実現される。

次に図２を参照すると、フローチャートが、本発明による全変動逐次近似法（ＴＶ−ＩＲ）の例示的なプロセスに含まれるステップを示している。初期化ステップＳ１０では、画像Ｘ０が０に初期化される。同じステップにおいて、反復回数（number of iteration）を示す値が第１の所定の数Ｋに初期化され、現在の反復回数が反復カウンタまたはメインループカウンタｋに維持される。ステップＳ２０では、緩和パラメータを用いたオーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）が、本発明によるＴＶ−ＩＲの例示的なプロセスにおける測定された投影データに対して実行される。ステップＳ２０は、中間的な画像Ｘ１を出力する。ＯＳＳＡＲＴの詳細なステップに関しては、図３を参照して後述する。

更に図２を参照すると、示されている例示的なプロセスには、特定のステップが反復されるおよび／または逐次的に実行される２つのループが含まれる。第１のループは、ステップＳ２０とステップＳ３０とステップＳ４０とを含むメインループまたは外部ループである。メインループは、この例示的なプロセスによるステップＳ１０においてメインループカウンタｋに割り当てられる第１の所定の数Ｋに従って反復される。メインループのそれぞれの反復が終了すると、メインループカウンタｋは、初期値Ｋから１ずつデクリメントされる。ある例示的なプロセスでは、第１の所定の数Ｋは、オプションであるが、投影データの順序のつけられた部分集合（オーダードサブセット）Ｎ_setsの数に基づいて決定される。これに関しては、図３を参照して後述する。メインループの中で、緩和パラメータのラムダを用いたオーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）がステップＳ２０において投影データに適用される。ウェイト値が、ＯＳＳＡＲＴのための緩和パラメータλに記憶される。

第２のループは、ステップＳ３０とステップＳ４０とを含むサブループまたは内側ループである。サブループは、全変動ステップ（ＴＶステップ）ループとも呼ばれる。サブループは、一実施形態では全変動カウンタまたはサブループカウンタｓに割り当てられた第２の所定の数Ｓに従って反復される。ステップＳ３０とステップＳ４０とは、この実施形態では、割り当てられた第２の回数Ｓだけサブループにおいて反復される。サブループのそれぞれの反復が終了すると、サブループカウンタｓは、初期値Ｓから１ずつデクリメントされる。別の実施形態では、ステップＳ３０とＳ４０とが、所定の条件に到達するまでサブループにおいて反復される。

ステップＳ３０では、離散全変動最小化法が、所定の技術に基づいて全変動を最小化するように、逐次的に実行される。ステップＳ３０が終了すると、画像Ｘ２が、ステップＳ２０からの中間画像Ｘ１から生成される。ステップＳ３０において生成される画像は、ステップＳ４０において中間画像ホルダまたは変数Ｘ１に割り当てられる。

ステップＳ４２では、サブループが終了したかどうかが、サブループカウンタｓまたは所定の終了条件に基づいて判断される。ステップＳ４２においてサブループが終了したと判断される場合には、Ｓ４４において、ステップＳ２０からステップＳ４０がメインループにおいて反復されるべきかどうかが更に判断される。他方で、ステップＳ４２においてサブループがまだ終了していないと判断される場合には、ステップＳ３０とステップＳ４０とがサブループにおいて反復される。

ステップＳ４４では、メインループカウンタｋに基づいて、メインループが終了したかどうかが判断される。ステップＳ４４においてメインループが終了したと判断される場合には、画像Ｘ２がステップＳ５０において出力される。他方で、ステップＳ４４においてメインループがまだ終了していないと判断される場合には、ステップＳ２０とステップＳ３０とステップＳ４０とが、メインループにおいて反復される。

次に図３を参照すると、フローチャートが図２のステップＳ２０における更なる行為またはサブステップを示しており、これらのステップは、本発明による全変動逐次近似法（ＴＶ−ＩＲ）の例示的なプロセスにおいて用いられる緩和パラメータを用いたオーダードサブセット同時代数再構成法（ＯＳＳＡＲＴ）に含まれる。ＯＳＳＡＲＴの示されているプロセスは、特定のステップが反復され逐次的に実行される１つのループを含む。このループにはステップＳ２３とステップＳ２４とステップＳ２５とが含まれ、特定の条件が満たされるようにではなく、所定の数Ｎ_setsに従って反復される。ループが終了すると、ＯＳＳＡＲＴは、ステップＳ２７において、変数Ｘ１に中間画像を出力する。

更に図３を参照しながら、ステップＳ２１からステップＳ２６までのそれぞれについて、以下で一般的に説明する。ステップＳ２１では、画像Ｘ０と測定された投影データｐとが、図２のステップＳ２０を通じて得られる。あるオプショナルな実装形態では、第１の所定の記録構造の中に記憶されているウェイト値の中の１つが、メインループのメインループカウンタｋの値に対応するインデックスｋに基づいて選択され、緩和パラメータλに割り当てられる。

ステップＳ２２では、測定された投影データｐが、部分集合と称されるＮ_sets個のグループに分割される。部分集合Ｎ_setsの数は、ビューの総数がサブビューの所定の数によって除算されるという特定の規則に基づく。たとえば、投影データが全部で１２００個のビューを有し、サブビューの所定の数が１５である場合には、部分集合の数は１２００を１５で除算して得られる８０である。あるいは、それぞれの部分集合の中のビューの数は、オプションであるが、変化する。ビューの総数が固定されている場合に、サブビューが多数あることに起因して部分集合の数が比較的少ない場合には、画像の収束は遅く、メインループにおいてより多くの反復が必要となる。他方で、サブビューが少数であることに起因して部分集合の数が比較的多い場合には、画像は高速で収束し、より少数の反復しか必要でない。同時に、画像はより多くのノイズを含む。上述した関係に基づくと、第１の所定の数Ｋの値は、オプションであるが、図２のステップＳ１０において修正される。

部分集合Ｎ_setsはそれぞれが多数のサブレイＮ_subraysを含み、この数は、アレイＮ_raysの総数を部分集合Ｎ_setsによって除算することによって決定される。部分集合Ｎ_setsのそれぞれに対して、ステップＳ２３は、画像Ｘ０を再投影し計算された投影データを形成する。それに対して、ステップＳ２４は、測定された投影データと計算された投影データとの間の正規化された差を逆投影する。更に詳細には、ステップＳ２３のある例示的なプロセスでは、システム行列の係数がどれもキャッシュされないレイトレーシング法を用いることにより、画像ボリュームを再投影する。更に、ステップＳ２３の一実施形態では、ある部分集合の中のすべてのレイを同時に再投影するのであるが、これは、オプションであるがパラレルに実装される。逆投影においては、ステップＳ２４の一実施形態では、画素駆動法を用いて、ある部分集合の中の正規化された差投影データのすべてを逆投影して、所望の更新された画像ボリュームを形成する。ステップＳ２４は画像ボリュームを形成するある部分集合の中のすべてのレイの和すなわち正規化された差投影データを逆投影するので、この動作もまた、オプションであるが、パラレルに実装される。これらの動作がすべての部分集合Ｎ_setsに対して適用されることにより、ＯＳＳＡＲＴのステップが１つ終了する。部分集合カウンタは、Ｎ_setsに初期化された後で、それぞれの動作が終了するたびに、１ずつデクリメントされる。

更に図３を参照すると、ステップＳ２５は、ステップＳ２４からの逆投影され重み付けされた画像を画像Ｘ０に加算することによって、画像Ｘ０を更新する。逆投影された画像に重み付けするには、図６を参照して後述するように緩和パラメータλが用いられる。結果的に、ステップＳ２６において、Ｎ_setsが既に終了したかどうかが判断される。Ｎ_setsがまだ終了していないとステップＳ２６で判断された場合には、ループは、ステップＳ２３に進むことによって、反復される。他方で、Ｎ_setsが既に終了したとステップＳ２６で判断された場合には、ステップＳ２５で生成された画像が、画像Ｘ１がステップＳ２７において出力される前に、中間画像ホルダまたは変数Ｘ１に割り当てられる。

次に図４を参照すると、フローチャートが図２のステップＳ３０の更なるステップを示しており、これらのステップは、本発明による全変動逐次近似法（ＴＶ−ＩＲ）のプロセスにおいて用いられる離散全変動（ＴＶ）最小化法を含む。離散ＴＶ最小化法の示されているプロセスは、特定のステップが反復され逐次的に実行されるサブループを含む。このサブループにはステップＳ１１０からステップＳ１６０までが含まれ、所定の回数および／または特定の条件が満たされるまで反復される。サブループが終了すると、離散ＴＶ最小化法は、画像を変数Ｘ２に戻す。

更に図４を参照しながら、ステップＳ１００からステップＳ１７０までのそれぞれについて、以下で一般的に説明する。ステップＳ１００では、方向インデックスＤＩが、それぞれの画素に対して、０などの所定の値に初期化される。この場合の所定の値０は、所与の画素に対する全変動に対する方向変化における”増加傾向と減少傾向との混在（未知数ともいう）”を意味する。同様にして、所定の値１は、所与の画素に対する全変動に対する方向変化における増加を意味し、他方で、所定の値−１は、所与の画素に対する全変動に対する方向変化における減少を意味する。ステップＳ１１０では、全変動（ＴＶ）の勾配ＧＲが所与の画素ｉに対して決定され、対応する値がＧＲ［ｉ］に記憶される。全変動の勾配ＧＲの符号は、所与の画素ｉに対し、＋１によって増加、−１によって減少、または、０によって混在を示す。対応する画素ｉに対し、方向インデックスＤＩ［ｉ］が初期化されたものとして０すなわち混在であるかどうかが、ステップＳ１２０において判断される。ステップＳ１２０において方向インデックスＤＩ［ｉ］が０すなわち混在であると判断されると、方向インデックスＤＩ［ｉ］は、ステップＳ１３０において、勾配ＧＲ［ｉ］の対応する符号に更新される。すなわち、全変動の勾配ＧＲがゼロよりも大きい場合には、増加値＋１が方向インデックスＤＩ［ｉ］に割り当てられる。同様に、全変動の勾配ＧＲがゼロよりも小さい場合には、減少値−１が方向インデックスＤＩ［ｉ］に割り当てられる。最後に、全変動の勾配ＧＲがゼロである場合には、混在０が方向インデックスＤＩ［ｉ］に割り当てられる。ステップＳ１３０の後では、離散全変動最小化の例示的なプロセスは、ステップＳ１４０に進む。他方、方向インデックスＤＩ［ｉ］が０すなわち混在ではないとステップＳ１２０で判断される場合には、方向インデックスＤＩ［ｉ］はステップＳ１３０で更新されず、離散全変動最小化の例示的なプロセスは、ステップＳ１４０に進む。

更に図４を参照すると、離散全変動最小化の例示的なプロセスは、方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とが所与の画素ｉに対して同じ符号を有するかどうか、ステップＳ１４０で判断する。ステップＳ１４０において、方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とが同じ符号を有すると判断される場合には、ステップＳ１５０で、この所与の画素ｉに対する画素値は所定の定数を用いて更新される。この所定の定数の例示的な値は、元の画像の評価された標準偏差ＳｔｄＤをＮで除算した値である。ただし、ここでＮは１０から１００までの範囲にある。ステップＳ１５０の後では、離散全変動最小化の例示的なプロセスは、ステップＳ１６０に進む。他方で、方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とが同じ符号を有していないとステップＳ１４０で判断される場合には、画素値ｉはステップＳ１５０で更新されずに、離散全変動最小化の例示的なプロセスはステップＳ１６０に進む。

離散全変動最小化は、所定の規則に従って終了する。ステップＳ１６０において、画像の中のすべての画素が処理されたかどうかが判断される。画像の中のすべての画素がまだ処理されてはいないと判断される場合には、離散全変動最小化の例示的なプロセスは、残っている未処理の画素のためにステップＳ１１０に戻る。他方で、既に画像の中の画素がすべて処理されたと判断される場合には、離散全変動最小化の例示的なプロセスはステップＳ１９０に進み、ステップＳ１９０で追加的な終了条件が調べられる。離散全変動最小化プロセスのある実施形態では、所定の終了条件は、ステップＳ１７０に示されているように、方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とがすべての画素において異なる符号を有すること、である。ステップＳ１７０で、画像の中の画素のいずれかが画素方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とにおいて異なる符号を有していないと判断される場合には、離散全変動最小化の例示的なプロセスは、追加的な処理のためにステップＳ１１０に戻る。しかし、ステップＳ１７０で、方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とにおいてすべての画素が異なる符号を有していると判断される場合には、離散全変動最小化の例示的なプロセスは終了する。

離散全変動最小化プロセスの代替の実施形態では、ステップＳ１７０における終了条件は、すべての画素の方向インデックスＤＩ［ｉ］と勾配ＧＲ［ｉ］とにおける符号が異なることの代わりに、ステップＳ１１０からステップＳ１６０までを所定の回数反復することによって達成される。ステップＳ１７０において、所定の回数の反復がまだ達成されていないと判断される場合には、離散全変動最小化の例示的なプロセスは、追加的な処理のためにステップＳ１１０に戻る。例示的な所定の反復回数Ｎは、元の画像の評価された標準偏差ＳｔｄＤを定数Ｃで除算した値である。ただし、ここでＣはノイズが除去された画像における評価された標準偏差である。他方で、所定の回数の反復が達成されたとステップＳ１７０で判断される場合には、離散全変動最小化の例示的なプロセスは終了する。上述した例示的および別の実施形態は、本発明による離散全変動最小化プロセスを限定しない。

図５Ａは、従来技術による全変動最小化の従来型のプロセスにおいて処理される単一の奇妙な画素の場合を示している。一番上の図では、単一の画素ｂが、隣接する画素ａおよびｂとの関係で示されており、単一の画素ｂの値は、初期状態では、その真の値よりも相当に大きいと想定されている。離散全変動最小化における最初の更新の後では、単一の画素ｂの値は、更新Ｉにおいて実線で示されているように、破線で表されている元の値と比較して減少している。離散全変動最小化における２回目の更新の後では、単一の画素ｂの値は、更新ＩＩにおいて実線で示されているように、破線で表されている元の値や二重の破線で表されている最初に更新された値と比較して更に減少している。

図５Ｂは、本発明による離散全変動最小化のプロセスにおいて処理される奇妙な画素クラスタの場合を示している。一番上の図では、画素２および３が、隣接する画素１および４との関係で示されており、画素２および３の値は、初期状態では、その真の値よりも相当に大きいと想定されている。離散全変動最小化における最初の更新の後では、単一の画素２の値は、更新Ｉにおいて実線で示されているように、破線で表されている元の値と比較して減少している。離散全変動最小化における２回目の更新の後では、今度は、単一の画素３の値が、更新ＩＩにおいて実線で示されているように、破線で表されている元の値と比較して減少している。図５Ａおよび５Ｂにおいて示されているように、画素のクラスタが調整される場合には、ただ１つの画素値だけが更新されるため、クラスタを処理するためにはより長い時間を要する。すなわち、雑多なクラスタをカバーする方が、単一の画素をカバーするよりも時間がかかるのが一般的である。

図６は、本発明による全変動逐次近似法（ＴＶ−ＩＲ）の例示的なプロセスにおいて用いられる一方向条件を伴う離散ＴＶ最小化法によって処理される例示的な画素の場合を示している。一般的に、画素値は、それぞれのステップでの勾配によって決定されるように、増加または減少する。一方向条件とは、それぞれの画素値は、第１のステップで決定される非ゼロである勾配の値に従って単に増加するまたは減少するというものである。換言すると、本発明による一方向条件では、離散全変動最小化の間を通じて方向の変化は同一のまま維持され、途中で変更されることはあり得ない。

更に図６を参照すると、一番上の図では、単一の画素ｉが、隣接する画素ｉ−１およびｉ＋１との関係で示されており、単一の画素ｉの値は、初期状態では、その真の値よりも相当に大きいと想定されている。すなわち、画素ｉは減少するはずであり、方向インデックスは最初の全変動最小化ステップにおいて減少を示すように設定される。最初の全変動最小化ステップの後では、単一の画素ｉの値は、更新Ｉにおいて実線で示されているように、破線で表されている元の値と比較して減少している。第２の全変動最小化ステップの後では、単一の画素ｉの値は、更新ＩＩにおいて実線で示されているように、破線で表されている元の値や二重破線で表されている最初に更新された値と比較して更に減少している。第３の全変動最小化ステップの後では、一方向条件のために、単一の画素ｉの値が、更新ＩＩＩにおいて破線で示されているように、二重破線で表されている最初に更新された値と比較して増加することが回避されている。換言すると、第３の更新ＩＩＩの後では、単一の画素ｉの値は、更新ＩＩにおいて実線で示されているように、破線で表されている元の値や二重破線で表されている最初に更新された値と比較して、変化せずに維持されている。

図７Ａから７Ｄには、本発明による全変動逐次近似（ＴＶ−ＩＲ）法の一実施形態において、離散ＴＶ最小化技術が様々な例示的な場合を処理することが示されている。一般的に、１次元離散分布における離散全変動は、次の数式（１）によって定義される。

全変動最小化は、次の数式（２）によって定義される。

ここで、ｎは、１、２、３、・・・であり、ａは定数である。全変動の離散勾配は離散全変動の勾配と実質的に同義的であり、本出願では数式（３）および（４）によって定義される。

１次元の場合には、数式（４）の値は、−２、０または＋２の中のいずれかである。

同様に、３次元画像ｆ（ｉ，ｊ，ｋ）に対する勾配方向は、数式（５）に見られるように所定の方向に沿った関数ｇｔｖによって決定される。

２次元の場合には、インデックスｋが画像スライスインデックスであるから、数式（５）の値は、−４、−２、０、＋２、＋４の中のいずれか１つである。同様に、所定の点における３次元画像ｆ（ｉ，ｊ，ｋ）に対する３次元ＴＶ勾配は、−６、−４、−２、０、＋２、＋４、＋６の中のいずれか１つである。

更に、整数ｍは、（ｉ，ｊ，ｋ）におけるドットを削除するように決定されるが、ｍは次の数式（６）によって定義される。

ここで、ａは所定の値であり、ｍは所定の閾値ｍ₀よりも大きく、点（ｉ，ｊ，ｋ）は孤立したノイズドットである。孤立したノイズドットは、ｆ（ｉ，ｊ，ｋ）＝ｆ（ｉ，ｊ，ｋ）＋ｍａ／２ｇｔｖ（ｉ，ｊ，ｋ）によって削除される。同様に、ドット値ｍは、（ｉ，ｊ，ｋ）におけるドットを削除するように決定され、ｍは数式（７）によって定義される。

ここで、ａは所定の値であり、ｍは所定の閾値ｍ₀よりも大きく、点（ｉ，ｊ，ｋ）は孤立したノイズドットである。孤立したノイズドットは、ｆ（ｉ，ｊ，ｋ）＝ｆ（ｉ，ｊ，ｋ）＋ｍａ／４ｇｔｖ（ｉ，ｊ，ｋ）によって削除される。同様に、数式（５）および（７）は、３次元画像ｆ（ｉ，ｊ，ｋ）に対する所定の体積に一般化することができる。

ここで図７Ａを参照すると、第１の信号保全パターンが、離散全変動最小化における例示的な１組の画素値に対して示されている。画素値の初期状態では、図７Ａの上側の図に示されているように、３つの画素値は、第１の画素ｉ−１から第２の画素ｉに増加し、第２の画素ｉから第３の画素ｉ＋１に増加する。すなわち、数式（８）は、勾配がゼロである間は、画素値は増加することを示している。

結果的に、離散全変動最小化の後では、第２の画素ｉでの更新された画素値は、図７Ａの更新後の下側の図に示されているように、同じままである。従来型のＴＶ最小化では、図７Ａのパターンに対するＴＶ勾配値は、連続的な画像モデルのためにゼロではないのが通常である。画素値はゼロでない値を用いて更新され、パターンが変更されるので、結果的にエッジのぶれが生じる。

次に図７Ｂを参照すると、第２の信号保全パターンが、全変動最小化における例示的な１組の画素値に対して示されている。画素値の初期状態では、図７Ｂの上側の図に示されているように、画素値は、第１の画素ｉ−１から第２の画素ｉに増加するが、第２の画素ｉから第３の画素ｉ＋１まで減少する。すなわち、数式（９）は、数式（４）により勾配が２である間は、画素値は増加した後で減少することを示している。

結果的に、全変動最小化の後では、図７Ｂの更新後の下側の図に示されているように、第２の画素ｉはその値が減少する。

次に図７Ｃを参照すると、第３の信号保全パターンが、全変動最小化における例示的な１組の画素値に対して示されている。画素値の初期状態では、図７Ｃの上側の図に示されているように、３つの画素値は、第１の画素ｉ−１から第２の画素ｉまで減少し、第２の画素ｉから第３の画素ｉ＋１まで減少する。すなわち、数式（１０）は、勾配がゼロである間は、画素値は減少することを示している。

結果的に、離散全変動最小化の後では、第２の画素ｉでの更新された画素値は、図７Ｃの更新後の下側の図に示されているように、同じままである。

次に図７Ｄを参照すると、第４の信号保全パターンが、全変動最小化における例示的な１組の画素値に対して示されている。画素値の初期状態では、図７Ｄの上側の図に示されているように、画素値は、第１の画素ｉ−１から第２の画素ｉまで減少するが、第２の画素ｉから第３の画素ｉ＋１に増加する。すなわち、数式（１１）は、勾配が−２である間は、画素値は減少した後で増加することを示す。

結果的に、離散全変動最小化の後では、図７Ｄの更新後の下側の図に示されているように、第２の画素ｉはその値が増加する。

図８Ａは、投影データから再構成された画像を示しているが、この画像は、非一様な態様でその画像の中に分布する高レベルのノイズを有している。

図８Ｂは、離散全変動最小化を用いて図８Ａからノイズを除去した画像を示している。全変動最小化は、上述した本発明による一方向条件を伴うことなく用いられている。結果的に、全変動最小化がノイズをある程度まで減少させたことをこの画像は示しているが、画像の中のいくつかの領域は、ノイズ除去の後で過渡な補償がなされている、すなわち、過剰に平滑化されているように見える。これらの過剰な平滑化がなされた領域は、当初に、比較的低レベルのノイズを有していたのである。

図８Ｃは、本発明による一方向条件を伴う離散全変動最小化を用いて図８Ａからノイズを除去した例示的な画像を示している。離散全変動最小化は、上述した本発明による一方向条件を伴って用いられている。結果的に、再構成された画像は、当初に比較的低レベルのノイズを含んでいた領域を、過渡な補償すなわち過剰な平滑化を行うことなく離散全変動最小化によってノイズが実質的に減少されたことを示している。

次に図９Ａを参照すると、グラフが、あるエッジを横断する一連の画素について、それらの画素値とＣＴ数との関係を、元の画像の場合は実線によって、離散全変動（ＤＴＶ）最小化の後の場合はドット状の線によって、全変動（ＴＶ）最小化の後の場合はポイントドット状の線によって、示している。これらの画素の中でエッジの上にある部分は、図９Ａにおける青色の矩形の箱によって示されている。図９Ｂは、図９Ａの青色の矩形の箱の中の画素のエッジ部分を拡大して示している。画素値の拡大されている部分は、元の画像とドット状の線によって示されている離散全変動（ＤＴＶ）の後との間で、実質的に等しい。他方で、画素値の拡大された部分は、ポイントドット線によって示されている全変動最小化の後のものとは若干異なっている。

特許請求の範囲に記載されるように、本発明の発明的な概念は、全変動最小化の離散的勾配（ＤＴＶ）と特定の再構成アルゴリズムに関する一方向条件とを含む特徴の組み合わせを、必ずしも要求するものではない。たとえば、全変動の離散的勾配（ＤＴＶ）の最小化を一方向条件を伴わずに適用する技術が、所定の逐次近似アルゴリズムまたは所定の解析的再構成アルゴリズムのいずれかと共に用いられる。同様にして、全変動の離散勾配（ＲＴＶ）の最小化を伴うことなく、一方向条件だけを適用する技術が、所定の逐次近似アルゴリズムまたは所定の解析的再構成アルゴリズムのいずれかと共に用いられる。更に、一方向条件だけを適用する技術が、離散全変動（ＤＴＶ）最小化技術または従来型の全変動（ＴＶ）最小化技術のいずれかと組み合わされる。これらの例は、単に例示のためのものであり、本発明による画像のノイズ除去に関する発明概念を限定することはない。

一般的に、上述した特徴は、本発明による画像に対していくつかの区別することができる効果を有する。たとえば、一方向条件は、元の画像が非均一的なノイズ分布を有していると仮定した場合に、ノイズ除去後の画像におけるノイズの一様性を向上させ、画像における低ノイズ領域における過剰な平滑化を回避する。他方で、ＤＴＶ最小化は、従来型のＴＶ最小化と比較すると、エッジ保全効果に関して優れた結果を生じる。ＤＴＶを用いると、エッジの上の画素値は、特定のノイズ除去処理の後で、実質的に影響を受けることがない。ＤＴＶ最小化と従来型のＴＶ最小化との差は、強いエッジに関しては実質的に無視できるが、他方で、ＤＴＶ最小化は、弱いエッジの保全を実質的に向上させる。もう一つの区別として、ＤＴＶは元の値を維持するが、ＴＶは、エッジの底部および頂部に近い画素について、値を若干平滑化するという点がある。

しかし、本発明の多くの特徴や効果について、本発明の構成および機能の詳細と共に以上で説明してきたが、本開示は説明目的であるに過ぎず、特に構成要素の形状や大きさや配列に関して、また、ソフトウェアやハードウェアや両者の組み合わせによる実装形態に関して詳細に変更がなされたとしても、それらの変化は、特許請求の範囲の記載において用いられている用語の広く一般的な意味によって示される全範囲に関し、本発明の原理の範囲に含まれることを理解されたい。

Claims (17)

1. 画像のノイズを除去する画像処理装置であって、
   ａ）前記画像の各画素に対応付けられる方向インデックスを初期化する手段であって、前記方向インデックスは、全変動の離散勾配に関する増加と減少とのいずれかを示す、前記初期化する手段と、
   ｂ）前記各画素について、前記全変動の離散勾配を決定する手段であって、前記離散勾配は前記増加又は前記減少を示す、前記決定する手段と、
   ｃ）前記方向インデックスと前記離散勾配とが等しい場合には、前記各画素の画素値を所定値だけ変更して前記全変動を減少させる手段と、
   ｄ）所定の条件に到達するまで、前記手段ｂ）から前記手段c）までを反復して実行させる手段と、を備える画像処理装置。
2. 前記所定の条件とは、前記画素のそれぞれに対する勾配方向において前記方向インデックスと前記離散勾配とが異なることである、請求項１記載の画像処理装置。
3. 前記所定の条件とは、前記手段ｂ）から前記手段ｃ）までが所定の回数反復して実行されることである、請求項１記載の画像処理装置。
4. 所定の方向における２次元画像ｆ（ｉ，ｊ，ｋ）に対する勾配方向は次の関数ｇｔｖによって決定される、
   請求項１記載の画像処理装置。
5. ドット値ｍは（ｉ，ｊ，ｋ）におけるドットを削除するために決定され、前記ドット値ｍは、ａを前記所定の値として、下記の数式によって定義される、
   請求項１記載の画像処理装置。
6. 孤立したノイズドットは、ｆ（ｉ，ｊ，ｋ）＝ｆ（ｉ，ｊ，ｋ）＋（ｍａ／２）ｇｔｖ（ｉ，ｊ，ｋ）によって削除される、請求項５記載の画像処理装置。
7. ドット値ｍは（ｉ，ｊ，ｋ）におけるドットを削除するために決定され、前記ドット値ｍは、ａを前記所定の値として、下記の数式によって定義される、
   請求項１記載の画像処理装置。
8. ドット値ｍが所定の閾値ｍ₀よりも大きいときには、ｆ（ｉ，ｊ，ｋ）の点は孤立したノイズドットである、請求項５乃至７のうちいずれか一項記載の画像処理装置。
9. 前記孤立したノイズドットは、ｆ（ｉ，ｊ，ｋ）＝ｆ（ｉ，ｊ，ｋ）＋（ｍａ／４）ｇｔｖ（ｉ，ｊ，ｋ）によって削除される、請求項８記載の画像処理装置。
10. 前記手段ａ）から前記手段ｄ）までは逐次近似法によって実行される、請求項１記載の画像処理装置。
11. 前記手段ａ）から前記手段ｄ）までは解析的再構成法によって実行される、請求項１記載の画像処理装置。
12. ｅ）測定されたデータから画像データを取得する手段と、
    ｆ）少なくとも単一のレイを用いて前記画像データを前方へ投影し、再投影されたデータを生成する手段と、
    ｇ）前記再投影されたデータと前記測定されたデータとの間の差を決定する手段と、
    ｈ）前記差に基づいて前記画像データを更新し、更新された画像を生成する手段であって、前記手段ｅ）から前記手段ｈ）までは前記手段ａ）の前に実行される、前記画像を生成する手段と、
    ｉ）前記更新された画像を前記手段ｆ）と前記手段ａ）から前記手段ｄ）までにおける前記画像データとして用い、前記手段ｆ）から前記手段ｈ）までを反復して実行させる手段と、を更に備える請求項１記載の画像処理装置。
13. ｅ）測定されたデータから画像データを取得する手段と、
    ｆ）少なくとも単一のレイを用いて前記画像データを前方へ投影し、再投影されたデータを生成する手段と、
    ｇ）前記再投影されたデータと前記測定されたデータとの間のデータ差を決定する手段と、
    ｈ）前記データ差を逆投影し、差画像を生成する手段と、
    ｉ）前記差画像に対して、前記手段ａ）から前記手段ｄ）までを実行する手段と、
    ｊ）前記差画像に基づいて前記画像データを更新し、更新された画像を生成する手段であって、前記手段ｅ）から前記手段ｉ）までは前記手段ａ）の前に実行される、前記画像を生成する手段と、
    ｋ）前記更新された画像を前記手段ｆ）と前記手段ａ）から前記手段ｄ）までにおける前記画像データとして用い、前記手段ｆ）から前記手段ｉ）までを反復して実行させる手段と、を更に備える請求項１記載の画像処理装置。
14. 画像のノイズを除去するためのＸ線コンピュータ断層撮影装置であって、
    画素を有する画像を投影データから再構成するための再構成デバイスと、
    所与の画素における全変動の勾配方向における増加又は減少が示す１組の所定の値を設定するための離散全変動（ＴＶ）最小化ユニットであって、前記画素のそれぞれに対して前記全変動の離散勾配を決定し、前記離散勾配は前記増加又は前記減少を示し、前記画素のそれぞれの画素値を、前記画素の中の前記１つの画素に対する前記勾配方向において方向インデックスと前記離散勾配とが等しい場合には、所定の値だけ変更して前記全変動を減少させ、所定の条件に到達するまで反復する、離散全変動（ＴＶ）最小化ユニットと、を備えるＸ線コンピュータ断層撮影装置。
15. 画像のノイズを除去する方法であって、
    ａ）前記画像の各画素に対応付けられる方向インデックスを初期化するステップと、前記方向インデックスは、全変動の離散勾配に関する増加又は減少を示す、
    ｂ）前記各画素について、前記全変動の離散勾配を決定するステップと、前記離散勾配は前記増加又は前記減少を示す、
    ｃ）前記方向インデックスと前記離散勾配とが等しい場合には、前記各画素の画素値を所定値だけ変更して前記全変動を減少させるステップと、
    ｄ）所定の条件に到達するまで、前記ステップｂ）から前記ステップｃ）までを反復するステップと、を備える画像処理方法。
16. 画像のノイズを除去する画像処理方法であって、
    ｉ）画素のそれぞれに対して全変動の離散勾配を決定するステップであって、前記離散勾配は増加又は減少を示す、ステップと、
    ｉｉ）前記画素のそれぞれの画素値を、前記離散勾配に従って、所定の値だけ変更するステップと、
    ｉｉｉ）所定の条件に到達するまで前記ステップｉ）と前記ステップｉｉ）とを反復するステップと、を備える画像処理方法。
17. 画像のノイズを除去する画像処理装置であって、
    ｉ）画素のそれぞれに対して全変動の離散勾配を決定する手段であって、前記離散勾配は増加又は減少を示す、前記決定する手段と、
    ｉｉ）前記画素のそれぞれの画素値を、前記離散勾配に従って、所定の値だけ変更する手段と、
    ｉｉｉ）所定の条件に到達するまで前記手段ｉ）と前記手段ｉｉ）とを反復して実行させる手段と、を備える画像処理装置。