JP6208271B2 - 医用画像処理装置 - Google Patents

Description

本発明の実施形態は、一般的にノイズ低減機能を有する医用画像処理装置に関する。

Ｘ線コンピュータ断層撮影（Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ：以下Ｘ線ＣＴと呼ぶ）装置において、より高い解像度を達成するために、ノイズは、Ｘ線ＣＴ画像が得られるようになった初期の頃からの重要な問題である。ノイズ低減技術を発展させることにおける主な関心事は、関心のある臨床医的に重要な生体組織の構造を維持することにある。この点に関して、低域通過フィルタに関する技術は、効果的なノイズ低減フィルタとしてよく知られている。しかしながら、ノイズ低減処理された画像において、臨床的な生体組織の構造は、失われる傾向にある。すなわち、低域通過フィルタの技術は、一般的に、画像（例えば、Ｘ線ＣＴ画像など）における大きな構造の検出感度を改良するために、この画像から高周波数成分を除去する技術である。しかしながら、低域通過フィルタは、画像におけるエッジ（画素値が急変する領域）を平滑化し、結果として小さい生体組織の構造の検出感度を低下させる。このため、低域通過フィルタは、画像の本質的な分解能を減少させることになる。

この低域通過フィルタに対して、いくつかの改良が実施されている。例えば、従来技術における低域通過フィルタの改良技術は、局所的にエッジが存在しないか判定し、無いとみなした近隣領域における複数のボクセル値を平均化することである。他の従来技術における低域通過フィルタの改良技術は、近隣の複数のボクセルの間の加重平均を実行するためのＸ線ＣＴのボリュームデータ内において、特別な局所的位置におけるノイズに関して繰り返し計算される標準偏差に基づいている。しかしながら、これらのノイズ低減方法は、関心のある特別な生体組織の構造を維持することを、特に目的としていない。

対照的に、近年のＸ線ＣＴにおけるノイズ低減の取り組みは、ノイズ低減フィルタとして、画像に依存する拡散係数を有する非等方拡散フィルタによって実施されている。ある画像における生体組織の構造は、非等方拡散フィルタにより指向的に平滑化される。具体例として、非等方拡散フィルタを用いたＸ線ＣＴ画像のノイズ低減方法は、長く伸びて鮮明化された管のような生体組織の構造に対し、エッジに沿う方向のみ平滑化を行う。

コンピュータによる複雑な計算にもかかわらず、非等方平滑化技術は、関心のある生体組織の構造を維持すること以上に、ノイズを低減するために注意が向けられている。いくつかの従来技術における非等方平滑化技術は、画像データにおいて関心のある小さい生体組織の構造を維持することに関して改良されている。しかしながら、これらの従来技術は、望ましくなく、フィルタ処理されないような領域を決定する必要がある。この領域は、潜在的に関連性のある維持されるべき生体組織の構造を示している。非等方拡散ノイズ低減方法とシステムは、このように、Ｘ線ＣＴ装置を含む多様なモダリティ、およびより広範囲に多次元の信号などを用いて発生された再構成画像における臨床的に重要な生体組織の構造を維持することに関して改良される問題が残っている。

目的は、医用画像における生体組織の構造を維持しつつ、ノイズを低減する処理を実行する医用画像処理装置を提供することにある。

本実施形態に係る医用画像処理装置は、再構成画像に対して非等方拡散処理を実行する非等方拡散処理部と、前記再構成画像の複数の画素において特定の画素に対応する特定の重みを発生する重み発生部とを具備し、前記非等方拡散処理部は、基準再構成画像と前記非等方拡散処理により平滑化された再構成画像との差に前記重みを乗じた項を非等方非線形拡散フィルタに加えた非同次非等方非線形拡散フィルタを用いて、前記再構成画像に対して、前記重みを有する前記非等方拡散処理を所定回数繰り返し実行すること、を特徴とする。

図１は、本実施形態に係り、Ｘ線コンピュータ断層撮影装置の構成を示す構成図である。 図２は、本実施形態に係る適応加重非等方拡散（Ａｄａｐｔｉｖｅｌｙ Ｗｅｉｇｈｔｅｄ Ａｎｉｓｏｔｒｏｐｉｃ Ｄｉｆｆｕｓｉｏｎ）方法に関する例示的なステップ、または本実施形態における処理の手順を示すフローチャートである。 図３Ａは、本実施形態に係り、あるシャープエッジに対するグラジエントとラプラシアンとに対応する強度レベルを示す図である。 図３Ｂは、本実施形態に係り、あるスローエッジに対するグラジエントとラプラシアンとに対応する強度レベルを示す図である。 図４Ａは、本実施形態に係り、適応加重非等方拡散が適用される前の最初の再構成画像の一例を示す図である。 図４Ｂは、本実施形態に係り、適応加重非等方拡散により発生された拡散画像の一例を示す図である。 図４Ｃは、本実施形態に係り、適応加重先鋭化非等方拡散により発生された拡散画像の一例を示す図である。

本実施形態は、Ｘ線コンピュータ断層撮影装置（以下、Ｘ線ＣＴと呼ぶ）に限定されない。例えば、本実施形態は、核磁気共鳴装置（Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ Ｉｍａｇｉｎｇ：以下ＭＲＩと呼ぶ）、および超音波診断装置などの他の医用画像診断装置に適用可能である。以下の実施形態では、一例として、Ｘ線ＣＴ装置に組み込まれた医用画像処理装置について説明する。なお、本実施形態は、医用画像診断装置に限らず、医用画像処理装置単体、および多次元の信号値が得られる医用画像診断装置などにも広く適用可能である。

図１は、本実施形態に係り、Ｘ線コンピュータ断層撮影装置の構成を示す構成図である。投影データ測定系は、ガントリ１を有する。ガントリ１は、Ｘ線管３と２次元アレイタイプのＸ線検出器５とを収容する。Ｘ線管３と２次元アレイタイプのＸ線検出器５とは、相対する向きに、被検体を横切って回転リング２に取り付けられる。被検体は、寝台の天板６上に載置される。回転リング２は、予め決定された回転軸の周りを回転する。このとき、寝台は、予め決定された速度で、回転軸に沿って、天板６を移動させる。架台寝台制御部９は、直線的な天板６の移動と、回転リング２の回転とを、同時に制御する。投影データ（医用画像データ）は、ヘリカル的なデータの収集中に天板６に対してＸ線管３がヘリカル軌道を描く間、または円形的なデータを収集中に天板６に対してＸ線管３が円形軌道を描く間に収集される。

Ｘ線管３は、コーンビームのような外形を有するＸ線束を放射する。コーンビームは、Ｘ線フィルタ４を通して、被検体に向けられる。Ｘ線制御部８は、トリガ信号を高電圧発生部７に供給する。高電圧発生部７は、トリガ信号を受信することにより、Ｘ線管３に高電圧を印加する。システムコントローラ１０は、Ｘ線制御部８と架台寝台制御部９との両方を制御する。この制御により、Ｘ線は、固定された角度間隔で連続的または断続的にＸ線管３から放射される。このとき、回転リング２は予め決定された回転速度で回転され、天板６は、予め決定された移動速度で移動される。

被検体を透過したＸ線は、２次元アレイタイプのＸ線検出器５によって、電気信号として検出される。Ｘ線検出器５は、複数のＸ線検出素子を有する。複数のＸ線検出素子は、１次元の列に配列される。複数の列は、２次元アレイを構成するために束ねられる。ある状況において、Ｘ線検出素子各々は、１つのチャンネルに対応する。データ収集部１１は、２次元アレイタイプのＸ線検出器５の各々のチャンネルから出力された出力信号を増幅する。データ収集部１１は、投影データを生成するために増幅した信号をディジタル信号に変換する。

処理部１２は、投影データに対して様々な処理タスクを実行する。例えば、処理部１２は、データのサンプリング、シフティング、フィルタリングおよび逆投影、再構成などの要求されたオペレーションを、投影データに対して実行する。

処理部１２は、予め決定されたアルゴリズムに従って、投影データを画像領域に逆投影する。一般的には、距離因子による重みが、逆投影されたデータに付加される。距離因子は、Ｘ線管３の位置から再構成された画素までの距離Ｌの逆数に比例する。なお、距離因子は、１／Ｌ、または１／Ｌ^２に比例してもよい。また、加重された付加データの冗長性が、画素ごとに逆投影処理が実行されている間に、任意に適用されてもよい。逆投影処理は、一般的に、Ｘ線発生源と検出器の受光素子とを結ぶレイ上に存在する画素におけるデータ値（逆投影データ）を得ることにある。なお、この処理は、様々な方法で実行されてもよい。

処理部１２は、３次元画像再構成におけるボクセル（３次元的なピクセル）各々において、Ｘ線の吸収を反映する逆投影データを決定する。コーンビームＸ線を用いた円形状のスキャン系において、画像化の領域（有効撮像視野）は、回転軸を中心とした半径Ｒの円筒形となる。処理部１２は、この画像化領域における複数のボクセルを定義する。処理部１２は、ボクセル各々に対する逆投影データを見出す。３次元画像データ、または断層画像データは、逆投影データに基づいて編集される。処理部１２、または他の部は、再構成画像のノイズ低減に関連するタスクを実行する。最終的に、表示部１４は、３次元画像、または３次元画像データに応じた再構成断層画像、後述する画像発生部１３３で発生された画像などを表示する。

本実施形態における後処理部１３は、３次元減衰マップ（ボリュームデータ）、または再構成画像などに対してノイズ低減処理を実行する。なお、後処理部１３は、投影データ、または投影データに関する信号値（例えば生データなど）に対して、ノイズ低減処理を実行してもよい。ボリュームデータ、または再構成画像は、上記で述べた測定、または投影データにより生成される。一般的に、ある生体組織の構造の分解能が少なくとも維持され、この構造またはこの構造の周囲に関するノイズが低減される場合、再構成されたＸ線ＣＴ画像の質は、改善されると予想される。Ｘ線ＣＴ画像における臨床的に重要な生体組織の構造は、大きい生体組織の構造と小さい生体組織の構造とを有する。例えば、大きい生体組織の構造（以下、大構造と呼ぶ）は、骨、大動脈、肝臓などであって、一方、小さい生体組織の構造（以下、小構造と呼ぶ）は、小さい管、石灰沈着、複数の気管支などである。

後処理部１３は、非等方拡散処理部１３１、重み発生部１３２、および画像発生部１３３を有する。非等方拡散処理部１３２は、様々なサイズの関連性のある生体組織の構造、またはこの構造の周囲のエッジを維持し、かつノイズを低減するために、再構成画像またはボリュームデータに対して、適応加重非等方拡散（Ａｄａｐｔｉｖｅｌｙ Ｗｅｉｇｈｔｅｄ Ａｎｉｓｏｔｒｏｐｉｃ Ｄｉｆｆｕｓｉｏｎ：以下、ＡＷＡＤと呼ぶ）処理を実行する。この非等方拡散処理により、大構造と小構造との周囲におけるノイズは、これらの構造のエッジが維持されて、低減される。更に、ＡＷＡＤは、エッジを維持しつつノイズを低減することで、小構造に関していくつかの長所を有する。多く存在する従来の非等方拡散技術は、大構造、または相対的に緩慢に変化する生体組織の構造、およびこれらの構造の周囲におけるノイズを低減させることができるが、ＡＷＡＤは、小構造に対して特に拡散画像に画像の実質的な詳細な構造を回復することができる。

本実施形態におけるＡＷＡＤの適用方法を理解するために、まず、一般的な非等方拡散処理の概念について、次の方程式について説明する。非等方拡散処理は、以下の方程式（１ａ）で定義されたように、指向性に依存する伝導係数Ｄを加えることにより、線形拡散方程式の拡張として与えられる。以下、説明を簡単にするために、非等方拡散処理が実行される対象は、ボリュームデータのボクセル値であるとする。なお、非等方拡散処理が実行される画像の対象として、再構成画像の画素値、または生データにおける信号値であってもよい。

ここで、記号“・”は、二つのベクトル（∇とＤ（ｕ）∇（ｕ））の内積を表す。Ｄは、拡散係数である。方程式（１ａ）の一般的な解釈は、以下のとおりである。ｕは、例えばボリュームデータにおける３次元座標と時間と（ｘ、ｙ、ｚ、ｔ）に依存する関数である。すなわち、ｕは、３次元座標と時間とに依存するスカラー場（ボクセル値場）に対応する。具体的には、ｕは、位置（３次元座標）と時間とによって規定されるボクセル値を規定するための関数（以下ボクセル値関数と呼ぶ）である。なお、微分は、ボクセル値の差分に対応する。以下、ｕは、ボクセル値関数として式を説明する。

ボクセル値関数ｕのボクセル値は、時間ｔに亘って拡散される。拡散は、拡散テンソルＤ（ｕ）により支配される。拡散係数Ｄ（ｕ）は、ボクセル値関数ｕのグラジエント（∇ｕ）、または局所的な画像の幾何（ボクセル値関数の幾何）に基づいている。拡散係数Ｄは、ボクセル値関数に依存する拡散テンソルである。方程式（１ａ）の左辺は、ボクセル値関数の変化（ボクセル値の強度変化）が、拡散テンソルにボクセル値関数のグラジエントを乗じたもののダイバージェンスに等しいこと、を示している。方程式（１ａ）の等式は、ボクセルごとに成り立つ。∂ｕ／∂ｔによって示されているように、非等方拡散処理は、ボクセル値関数ｕが繰り返しの時間ｔに亘って変化することを示している。

離散化された時間形式では、非等方拡散処理は、時間に応じて繰り返しの性質を明示的に示すために方程式（１ｃ）、（１ｂ）で以下のように定義される。

ここで、ｕは、加工処理、または拡散されたボクセル値関数である。ｕ（ｔ_０）は、ｎ＝０におけるボクセル値関数であって、非等方拡散処理が実行される前のボクセル値関数である。具体的には、非等方拡散処理が実行される前のボリュームデータのボクセル値である。ｔ_ｎとｔ_ｎ＋１とは、繰り返し回数ｎとｎ＋１とにおける時間的な段階をそれぞれ示している。

上述したように、方程式（１ａ）乃至（１ｃ）で定義されるような非等方拡散処理は、ボリュームデータにおけるボクセル値が、時間と空間との関数であることを要求する。したがって、方程式（２）は、非等方拡散処理前のボリュームデータにおけるボクセル値関数ｕ_０を定義する。ボクセル値関数ｕ_０は、予め決定された時間ｔ＝０におけるボクセル値関数（以下、初期ボクセル値関数と呼ぶ）である。ｕ_０は、基準医用画像データに含まれる。

ここで、ｘ、ｙ、およびｚは、例えばボリュームデータが定義された空間上でのボクセルの座標をそれぞれ示している。同様に方程式（３）は、あらゆるボクセル値関数ｕを定義する。方程式（３）は、順次、予め決定された時間ｔ≠０におけるボクセル値関数ｕを定義する。

ここで、ｘ、ｙ、およびｚは、例えばボリュームデータが定義された空間上でのボクセルの座標をそれぞれ示している。言い換えると、方程式（２）は、初期時刻ｔ＝０におけるボリュームデータを定義する。一方、方程式（３）は、次の時刻におけるボリュームデータを定義する。

さらに方程式（４）は、ノイマン（Ｎｅｕｍａｎｎ）境界条件を定義する。ノイマン境界条件とは、境界における微分がゼロであるということを意味する。それ故に、境界において、拡散過程でボクセル値の出入り（エネルギーの喪失）を防止する。

拡散係数Ｄは、式（５ａ）と式（６）とにより定義されるパラメータｓの関数である。

ここで、ｓは、平滑化されたボクセル値関数ｕの局所的なグラジエントの大きさである。Ｇ_σｕで示されたように、ボクセル値関数ｕは、標準偏差σを有するガウシアンフィルタＧにより平滑化される。ｋは、閾値または予め決定された値である。与えられたボクセルにおける局所的なグラジエントの大きさｓが予め決定された閾値ｋより小さい場合、ボクセル値は、ノイズとして取り扱われる。他方、与えられたボクセルにおける局所的なグラジエントの大きさｓが予め決定された閾値ｋより大きい場合、ボクセル値は、線またはエッジのような関心のある構造に対応する信号として取り扱われる。

予め決定された閾値ｋは、以下の式（５ｂ）、（５ｃ）に基づいて決定される。

一般的に、ｋ_０は、式（５ｂ）において、所定の数に設定される。所定の数とは、例えば２である。なお、ｋは１であってもよい。一方、ｄは、以下の式（５ｃ）で決定される。

式（５ａ）におけるグラジエントの局所的な大きさｓはボクセルごとに変化するため、グラジエントの局所的な大きさｓは、空間的な位置の関数である。式（５ｃ）で定義されたグラジエントの大きさの平均は、一つのボリュームデータに対して、一つの数で与えられる。

上述した非等方拡散処理は、一般的に、相対的に緩慢に変化する生体組織の構造内または相対的に大きい構造内、またはこれらの構造の周囲におけるノイズを低減するために十分にうまく機能する。

非等方拡散処理を用いるノイズ低減方法が実施されている間に、ボリュームデータに対して拡散の効果または平滑化の効果を釣り合わせるために、本実施形態における非等方拡散処理には、適切に重み付けされたソース項（ｓｏｕｒｃｅ ｔｅｒｍ：非同次項または非斉次項）が組み込まれる。適切に重み付けされたソース項は、繰り返しの各々段階に対して適用される。従って、適切に重み付されたソース項は、繰り返しの各々の段階に対して、拡散されたボクセル値を調整する。このアプローチは、非等方拡散処理の繰り返しの完了後に、拡散されたボクセル値に重み付けすることと対照的である。従って、適応加重等方拡散処理を用いた本実施形態におけるノイズ低減方法は、ボリュームデータにおける臨床的に重要な生体組織の構造を維持するために、非等方拡散処理の繰り返しにおけるいくつかの段階ごとに、非等方拡散処理前の画像において失われた詳細な生体組織の構造に関するボクセル値を、拡散されたボクセル値に加算する。

上記で定義された非等方拡散処理は、本実施形態の非等方拡散処理部１３１における非等方拡散フィルタで実行される。また、上記で定義された非等方拡散処理は、本実施形態の非等方拡散処理のステップで使用される。どちらの場合においても、非等方拡散処理は、繰り返し実行される。つまり、拡散されたボクセル値関数ｕは、繰り返しにおける最新の段階に対して求められる。最新の計算された拡散されたボクセル値関数ｕ（ｔ）は、次の繰り返しの段階で、前の段階に拡散されたボクセル値関数ｕ（ｔ−１）として用いられる。

繰り返しの各々の段階に対して、本実施形態における非等方拡散処理部１３１は、さらに、最新の段階におけるボクセル値関数と非等方拡散フィルタが適用される前のボクセル値関数との差に重みを乗じた非同次項を、非等方拡散処理が実行された結果に加える。同様に、本実施形態における非等方拡散処理の方法、または非等方拡散処理に関するステップは、最新の段階おけるボクセル値関数と非等方拡散フィルタが適用される前の初期ボクセル値関数との差に重みを乗じた非同次項を、非等方拡散処理が実行された結果に加える。

適応重み付け項（非同次項）は、以下の方程式（７ａ）、（７ｂ）および（７ｃ）のように組み込まれる。方程式（７ａ）、（７ｂ）および（７ｃ）は、適応加重非等方拡散（Ａｄａｐｔｉｖｅ Ｗｅｉｇｈｔｅｄ Ａｎｉｓｏｔｒｏｐｉｃ Ｄｉｆｆｕｓｉｏｎ：ＡＷＡＤ）を表している。ＡＷＡＤは、本実施形態により、複数のステップ各々において実行される。

ここで、Ｗは、非同次項における適応加重要素（以下重みと呼ぶ）である。離散的な時間形式では、方程式（７ａ）は、以下のように、方程式（７ｂ）および（７ｃ）で任意に定義される。

非等方拡散処理の基本的な概念は、方程式（１ａ）、（１ｂ）および（１ｃ）に関してすでに述べたため、以下、方程式（７ａ）、（７ｂ）および（７ｃ）については、非同次項に焦点をあてて解説する。後処理部１３における重み発生部１３２は、方程式（７ａ）、（７ｂ）および（７ｃ）に用いられる重みＷを発生する。重み発生部１３２は、初期ボクセル値関数（非等方拡散処理前のボリュームデータ）に基づいて、エッジマップを発生する。重み発生部１３２は、エッジマップに基づいて、ボクセルの座標にそれぞれ対応する重みを発生する。重みＷは、非等方拡散処理が実行される前の初期ボクセル値関数と、非等方拡散処理が実行された後のボクセル値関数との時間的に起因する差に掛けられる。方程式（７ａ）において、上記差は、（ｕ−ｕ_０）で表される。方程式（７ｂ）、（７ｃ）における離散的な時間形式において、同じ差は、（ｕ（ｔ_ｎ）−ｕ_０）で表される。

方程式（７ａ）、（７ｂ）および（７ｃ）における重みＷは、時間に依存または時間に独立して決定される。なお、重みＷは、空間的に依存していてもよい。式（８）で定義されたように以下の例では、重みＷは、時間に不変な項として示される。式（８）で定義されたように重みＷは、０と１との間の規格化された値として決定される。

ここでＬは、予め決定された初期ボクセル値関数ｕ_０のエッジマップである。ある典型的なエッジマップの発生において、ガウス関数のラプラシアンが用いられる。他の典型的な実施においては、エッジマップは、ソーベル（Ｓｏｂｅｌ）エッジ検出（ガウス関数のグラジエント）、曲率エッジ検出方法、またはこれらのような知られたエッジ検出方法の組み合わせなどのような他のエッジ検出方法によって発生される。例えば、ソーベルエッジ検出方法は、行列（８ａ）で定義される以下の２種のフィルタリングマスクで与えられる。

式（９）、（１０）で定義されるようなガウス関数のラプラシアンを用いて、ある典型的な実施では、平滑化された初期ボクセル値関数ｕ_０のラプラシアンΔｕ_０に基づいたエッジマップは、単純なある種の長所を有する。このエッジマップでは、シャープなエッジとエッジの遷移帯とが検出される。エッジの遷移帯とは、エッジが現れる程度に、ボクセル値が変化する領域である。

平滑化された初期ボクセル値関数ｕ_０のラプラシアンΔｕ_０に基づいたエッジマップは、緩慢に変化する境界を減少し、ゼロまたはゼロ近傍を中央値とした領域を一様にしたマップである。

ここで、ｕ_ｓは、平滑化カーネルｈ（ｐ、ｑ、ｒ）によって平滑化された初期ボクセル値関数ｕ_０である。平滑化カーネルｈ（ｐ、ｑ、ｒ）は、任意の単純な移動平均フィルタ、またはガウシアン加重平滑化フィルタである。

単純な移動平均フィルタ、またはガウシアン加重平滑化フィルタは、以下の式（１１ａ）、（１１ｂ）により定式化される。

ここで、τは、移動平均フィルタの長さであり、σは、ガウシアンフィルタのサイズパラメータである。

上述したように、非等方拡散処理されたボクセル値関数に、非同次項が付加される。非同次項は、時間的に与えられるボクセル値関数の差と規格化された重みＷとの積である。この差は、あらゆる拡散処理前の初期ボクセル値関数と、繰り返しにおいて最新の非等方拡散処理後のボクセル値関数との差である。言い換えれば、本実施形態と本実施形態におけるステップとは、繰り返し処理の終了時に重みを加えるのではなく、繰り返しのあらゆる段階で適切にボクセル値関数に非同次項が加えられる。

図２によれば、フローチャートは、ＡＷＡＤの方法、または処理を有する典型的なステップを示している。後処理部１３と処理部１２とを有する上述した本実施形態により、下記に述べる設定またはタスクが実行される。なお、ＡＷＡＤの方法、または処理は、フローチャートに示されたこれらのステップに限定されない。

ＡＷＡＤの典型的な処理において、ステップＳ１０は、ＡＷＡＤのあらゆる繰り返しの前に、初期ボクセル値関数のエッジマップを発生する。エッジマップの発生は、式（８ａ）、（９）、（１０）、（１１ａ）、（１１ｂ）などにより、先に述べたとおりである。重み発生部１３２は、エッジマップに基づいて発生された重みを、ＡＷＡＤに関する繰り返しの段階各々で適用するために記憶する。上述した例では、エッジマップの発生は、空間に依存し、時間不変である。つまり、重みは、ボリュームデータに対して設定された座標系の一組の座標によって規定されるボクセルに対して決定される。なお、ある実施形態では、重みは、時間に亘って変化しない。また、別の実施形態では、エッジマップは、空間および時間両方に依存する重みを生成するために、繰り返しごとに発生されたボクセル値関数に基づいて発生されてもよい。

ステップＳ１０は、また、エッジ検出の前の初期ボクセル値関数に関するいくつかの前処理を有する。式（１０）、（１１ａ）、（１１ｂ）に関して上述したように、初期ボクセル値関数ｕ_０は、ｈ（ｐ、ｑ、ｒ）のような平滑化カーネルによって平滑化される。

最後にステップＳ１０は、ある特徴に関するエッジマップを発生する。発生されたエッジマップに基づいて、非等方拡散処理における非同次項で用いられる重みを非同次項に提供している間（以下、重み付けステップと呼ぶ）に、特定のエッジ検出に関する特徴は、一般的に、重み付けステップで利用される重みを決定する。このため、選択されたエッジ検出は、部分的に、ＡＷＡＤ処理の最終結果に作用する。エッジ検出に関する方法は、グラジエントマップ、ラプラシアン、曲率、構造テンソルおよびヘッシアン行列の固有値などのような様々な方法である。なお、エッジ検出は、上記検出方法の組み合わせに基づいて実行されてもよい。

さらに、図２を参照すると、ステップＳ２０は、ボリュームデータの各々のボクセル値に対して、予め決定された非等方拡散処理を実行する。非等方拡散処理の典型的なステップは、方程式（１ａ）、（１ｂ）、（１ｃ）、（２）、（３）、（４）、（５ａ）、（５ｂ）、（５ｃ）と（６）に関して、先に述べたとおりである。

ステップＳ３０において、非等方拡散処理は、ステップＳ２０で拡散された結果に、非同次項を加える。ステップＳ２０で拡散された結果に非同次項を加えることは、式（７ａ）、（７ｂ）、（７ｃ）と（８）で関連づけられたタスクによって完成される。式（７ａ）、（７ｂ）、（７ｃ）は、非等方拡散処理自身に関連する項と、非同次項とを有する。ステップＳ２０とステップＳ３０とは、典型的なフローチャートにおいて二つの一連のステップとして示されているが、ＡＷＡＤの処理は、様々な方法で実行される。一つの代替的な実施形態または処理において、ステップＳ２０とＳ３０とは、並行して同時に実施されてもよい。どんな場合でも、ステップＳ２０における拡散された結果に、繰り返しの各々の段階に対して、ステップＳ３０において、非同次項が付加される。式（８）は、ステップＳ１０で発生されたエッジマップに基づいて規格化された重みＷを定義する。

ステップＳ３０に関して述べたように、ＡＷＡＤの処理は、ある好ましい処理において繰り返し実行される。これに関して、ステップＳ４０は、ステップＳ２０とＳ３０とにおけるＡＷＡＤの処理を終わらせるか否かについて決定する。ステップＳ４０が、ステップＳ２０とＳ３０とにおけるＡＷＡＤの処理を終わらせると決定した場合、実例となるステップは、その処理を終える。一方、ステップＳ４０が、ステップＳ２０とＳ３０とにおけるＡＷＡＤの処理を終わらせないと決定した場合、実例となるステップは、ステップＳ２０とＳ３０とを繰り返し実行することによりその処理を続行する。式（７ａ）、（７ｂ）、（７ｃ）に関して述べたように、最新のボクセル値関数ｕ（ｔ＋１）は、繰り返しの次の段階が実行される前に、繰り返しの現在の段階における非等方拡散処理された結果ｕ（ｔ）に代入される。

代わりの実施形態において、ステップＳ３０は、予め決定された先鋭化フィルタを適用する任意のステップを有していてもよい。式（１２）と（１３）は、適応加重先鋭化非等方拡散（Ａｄａｐｔｉｖｅ Ｗｅｉｇｈｔｅｄ Ｓｈａｒｐ Ｓｏｕｒｃｅ Ｄｉｆｆｕｓｉｏｎ：以下、ＡＷＳＳＡＤと呼ぶ）に関する処理を定義する。

ここで、Ｓ（ｕ_０）は、式（１３）で定義される予め決定されたアンシャープフィルタを用いて、初期ボクセル値関数ｕ_０を先鋭化した先鋭化ボクセル値関数である。

ここでｗは、アンシャープフィルタのパラメータである。ｗ＝０．８で、適度に良い結果が得られる。非等方拡散処理部１３１は、非等方拡散処理が所定回数に亘って実行されたボクセル値関数のうち、最新のボクセル値関数と先鋭化したボクセル値関数との差に重みを乗じた非同次項を、最新の非等方拡散処理の結果に加える。これにより、非等方拡散処理後におけるエッジが先鋭化される。

ステップＳ４０は、予め決定された終了条件に基づいて、ＡＷＡＤまたはＡＷＳＳＡＤの処理がさらに繰り返し実行されるべきか否かを、決定する。一般に、終了条件または繰り返しの回数を決定する方法は、二つある。一つの終了条件は、最新の繰り返しの結果による平均２乗誤差と、前の繰り返しの結果による平均２乗誤差とを比較することに基づくものである。これら２つの平均二乗誤差が予め設定された終了値より小さい場合、処理は停止される。他の終了条件は、繰り返しの各々の段階でノイズレベルを決定するためのノイズ評価ツールに基づくものである。最新の繰り返しの結果によるノイズレベルが、初期のノイズレベルから、５０％のノイズレベルのような予め設定されたノイズ低減目標へ、減少した場合、繰り返しは停止される。

非等方拡散処理部１３１は、予め決定された時間スケールｔを記憶する。非等方拡散処理部１３１は、ＡＷＡＤまたはＡＷＳＳＡＤの処理を、前の繰り返しに関するボクセル値関数ｕ（ｔ_ｎ）に基づいてボクセル値関数ｕ（ｔ_ｎ＋１）を計算するために、特定の時間間隔Δｔで実行する。ＡＷＡＤまたはＡＷＳＳＡＤの処理における繰り返しの回数は、ｔ／Δｔで与えられる。特定の時間間隔Δｔは、処理の繰り返しの安定性のために以下の式（１４）により決定される。

ここで、Ｎは、非等方拡散処理の対象となるデータの次元の数である。例えば、非等方拡散処理の対象となるデータがボリュームデータである場合、ボリュームデータの次元は３なため、Ｎ＝３となる。なお、非等方拡散処理の対象となるデータが再構成画像などの２次元画像データである場合、２次元画像データの次元は２なため、Ｎ＝２となる。式（１４）に基づいて、Δｔの値は、繰り返し間隔の上限を制限する。繰り返し間隔Δｔがあまりにも小さい場合、ＡＷＡＤまたはＡＷＳＳＡＤの処理は、全般的なノイズ低減処理に係る時間を実質的に低下させる。ある実施形態において、Δｔは、およそ１／４から上限を１として設定される。例えば、Δｔは、３次元再構成画像のノイズ低減で、Ｎ＝３に対して１／８、１／１２または１／１６のような値として選択される。

画像発生部１３３は、非等方拡散処理部１３１で処理された画像データ（ボリュームデータなど）に基づいて、エッジが強調または先鋭化され、ノイズが低減された医用画像を発生する。

図３Ａ、図３Ｂは、シャープエッジまたは緩慢なエッジに対するグラジエント、ラプラシアンでの強度レベルを示している。強度レベルとは、例えばボクセル値、画素値、信号値などである。ステップＳ１０に関して上述したように、エッジマップは、特定のエッジ検出によりある特徴とともに発生される。この点に関して、選択されたエッジ検出は、ＡＷＡＤまたはＡＷＳＳＡＤの処理に関する最終的な結果に、影響を部分的に与える。エッジマップは、厚いまたは緩慢に変化するエッジまたは境界よりも小さいまたはシャープなエッジを強調するために発生される。ＡＷＡＤまたはＡＷＳＳＡＤの一つの利点は、画像がノイズ低減されると同時に小さいエッジを保持することにある。この目的のために、重み発生部１３２は、ラプラシアンを用いてエッジマップを発生する。

図３Ａは、シャープなエッジのグラジエントとラプラシアンとに対応する強度レベルを示す。一番上のグラフは、一つの具体例としての鮮明なエッジに関する強度レベルを示す。強度レベルは、ゼロ強度レベルを横切って速く変化する。中間のグラフは、グラジエントの強度レベルを示す。グラジエントの強度レベルの凸部は、一番上の典型的なシャープなエッジに対応する。エッジにおける強度レベルが負の値から正の値へゼロ点で変化するとき、グラジエントの強度レベルは、最高レベルにある。一番下のグラフは、ラプラシアンの強度レベルを示す。ラプラシアンの強度レベルは、一次元の場合において２階微分となる。典型的なシャープエッジに関するエッジを増大または減少させるために、１階微分は、凸部または凹部を与える。凸部または凹部の頂点は、１階部分の絶対値が最大となる場所にある。１階微分の頂点は、二階微分においてゼロとなる。２階微分がゼロとなる両側で、２階微分またはラプラシアンにおいて凸部または凹部が現れる。１階微分が凸部の場合、２階微分の凸部は正であり、１階微分が凹部の場合、２階微分の凹部は負となる。ある実施形態において、２階微分の絶対値は、エッジにおける重みとして選択される。その結果、重みは、エッジ情報を維持するためにエッジに沿った幅に含まれるボクセル値関数に対応付けられる。エッジにおけるラプラシアンに対応する重みはゼロだが、エッジは、非等方拡散処理において実質的に維持される。シャープエッジだけは、より大きいラプラシアン値を有する。これらのエッジは、ラプラシアンに関する重みを用いることで実質的に維持される。

図３Ｂは、ある緩慢なエッジのグラジエントとラプラシアンとに対応する強度レベルを示す。一番上のグラフは、一つの具体例としての緩慢なエッジに関する強度レベルを示す。強度は、図３Ａにおける一番上の図と比較して、ゼロ強度レベルを横切ってゆっくり変化する。中間のグラフは、グラジエントの強度レベルを示す。グラジエントの強度レベルは、一番上の典型的な緩慢なエッジに対応する。エッジにおける強度レベルが負の値から正の値へゼロ点で変化するとき、グラジエントの強度レベルは、より小さい強度レベルでゆっくりと変化する。一番下のグラフは、ラプラシアンの強度レベルを示す。一番下のグラフは、ラプラシアンの強度レベルがほとんど変化しないことを示している。シャープエッジと比較して、緩慢に変化するエッジ、または厚い境界を有する大きいエッジは、ノイズに比べて一般的に小さいラプラシアン値を有する。この理由のため、ある実施形態においては、ラプラシアンは、実質的に無視される。

図４Ａ、図４Ｂ、図４Ｃは、画像発生部１３３により発生された画像の一例を示す図である。図４Ｂ、図４Ｃは、適応加重非等方拡散の処理に関するいくつかの効果を示す図である。図４Ａは、適応加重非等方拡散が適用される前の初期の再構成画像を示す図である。図４Ｂは、繰り返し回数は４であって、繰り返しの間隔は１／８であるＡＷＡＤの処理により発生された拡散画像を示す図である。なお、閾値として、２を用いてもよい。図４Ｃは、繰り返し回数は１３、繰り返しの間隔は１／８であって、ｗ＝０．８であるＡＷＳＳＡＤの処理により発生された拡散画像を示す図である。なお、閾値として、２を用いてもよい。

要約すると、本実施形態は、ＡＷＡＤまたはＡＷＳＳＡＤの処理が実行される前のデータ（ボリュームデータ、再構成画像などの画像データ）に基づいて発生されたエッジマップを用いて、画素またはボクセルごとに重みを発生することができる。続いて、所定回数に亘って繰り返し実行される非等方拡散処理ごとに、非等方拡散処理の前の画像データと最新の非等方拡散処理の結果との差にボクセルごとに重みを乗じた非同次項を、最新の非等方拡散処理の結果に加えることができる。すなわち、本実施形態は、初期の画像データにおける情報に基づいてエッジおよび構造の特徴を調節するために、負のフィードバックメカニズムを適用する。これにより、最終的な結果は、ノイズが低減された画像であって、バランスがとれ、収斂され、特徴づけられた画像が得られる。すなわち、ノイズが低減され、かつ生体組織の細かな構造が維持された画像（エッジが強調された画像、エッジが先鋭化された画像）を発生することができる。

また、本実施形態によれば、エッジ検出を、非等方拡散処理の前の画像データに合わせて選択することができる。これにより、非等方拡散処理の前の画像データに適したエッジ検出が可能となる。加えて、本実施形態によれば、非等方拡散処理の前の画像データの次元に追うおじた繰り返し回数を設定することができる。これにより、非等方拡散処理の処理時間を適したものにすることができる。なお、本実施形態の技術的思想を医用画像診断装置で実現する場合には、医用画像診断装置は、非等方拡散処理部１３１、重み発生部１３２、画像発生部１３３を有する後処理部１３を具備することを特徴とする。

本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

１…ガントリ、２…回転リング、３…Ｘ線管、４…Ｘ線フィルタ、５…Ｘ線検出器、６…天板、７…高電圧発生部、８…Ｘ線制御部、９…架台寝台制御部、１０…システムコントローラ、１１…データ収集部、１２…処理部、１３…後処理部、１４…表示部、１３１…非等方拡散処理部、１３２…重み発生部、１３３…画像発生部。

Claims (11)

1. 再構成画像に対して非等方拡散処理を実行する非等方拡散処理部と、
   前記再構成画像の複数の画素において特定の画素に対応する特定の重みを発生する重み発生部とを具備し、
   前記非等方拡散処理部は、
   基準再構成画像と前記非等方拡散処理により平滑化された再構成画像との差に前記重みを乗じた項を非等方非線形拡散フィルタに加えた非同次非等方非線形拡散フィルタを用いて、前記再構成画像に対して、前記重みを有する前記非等方拡散処理を所定回数繰り返し実行すること、
   を特徴とする医用画像処理装置。
2. 前記重み発生部は、前記複数の重みを一度に発生する
   請求項１に記載の医用画像処理装置。
3. 前記重み発生部は、前記複数の重みを、前記非等方拡散処理の繰り返し毎に発生する
   請求項１に記載の医用画像処理装置。
4. 前記非同次非等方非線形拡散フィルタは、
   により与えられ、ｕ（ｔ_ｎ＋１）は、前記非等方拡散処理の繰り返し回数がｎ＋１回目における画素値であり、ｕ（ｔ_ｎ）は、前記繰り返し回数がｎ回目における画素値であり、ｕ０は前記非同次非等方非線形拡散フィルタにかけられる前の初期画素値であり、Ｗは、前記初期画素値に対応する前記重みであること、
   を特徴とする請求項１に記載の医用画像処理装置。
5. 前記基準再構成画像は、前記非等方拡散処理が実行される前の再構成画像であること、
   を特徴とする請求項１または４に記載の医用画像処理装置。
6. 前記重み発生部は、ラプラシアンエッジ検出を用いて、前記重みを発生すること、
   を特徴とする請求項１乃至５のうちいずれか一項に記載の医用画像処理装置。
7. 前記重み発生部は、ソーベルエッジ検出を用いて、前記重みを発生すること、
   を特徴とする請求項１乃至５のうちいずれか一項に記載の医用画像処理装置。
8. 前記重み発生部は、曲率エッジ検出を用いて、前記重みを発生すること、
   を特徴とする請求項１乃至５のうちいずれか一項に記載の医用画像処理装置。
9. 前記再構成画像は、前記非等方拡散処理の前に予め平滑化された再構成画像であること、
   を特徴とする請求項１乃至８のうちいずれか一項に記載の医用画像処理装置。
10. 前記基準再構成画像は、前記非等方拡散処理される前の再構成画像とアンシャープフィルタのパラメータとに基づいて、先鋭化された再構成画像であること、
    を特徴とする請求項１と請求項４と請求項５とのうちいずれか一項に記載の医用画像処理装置。
11. 前記先鋭化された再構成画像は、
    により与えられ、Ｓ（ｕ_０）は、前記先鋭化された再構成画像の画素値であり、ｕ_０は、前記非同次非等方非線形拡散フィルタにかけられる前の初期画素値であり、ｕは、前記非等方拡散処理の繰り返し回数がｎ回目における画素値であり、ｗは、前記パラメータであること、
    を特徴とする請求項１０に記載の医用画像処理装置。