JP2015024155A - 画像処理装置 - Google Patents

Abstract

【課題】簡易な手法にして実質的にすべての画像診断モダリティについて人間の介入なしで信頼性の高いノイズインデックスを推定すること。
【解決手段】本実施形態に係る画像処理装置は、画像内の複数の局所領域にそれぞれに対応する複数の局所ノイズインデックスを計算する計算部と、前記計算されたノイズインデックスからノイズレベルと度数との対応を表すヒストグラムを作成する作成部と、前記作成されたヒストグラムに基づいて前記画像に関するノイズインデックスとして前記局所領域の標準偏差又は擬似標準偏差を推定する推定部とを有する。
【選択図】図２

Description

本発明は、画像のノイズレベルを推定する画像処理装置に関する。

ノイズレベルの信頼性は画像品質評価に重要である。ノイズレベルの評価には、主観的視覚的評価を必要とする。主観的評価に必要な関心領域内のノイズレベルの測定には複雑な手動操作を必要とする。主観的評価は、時間がかかり、高価であり、均一な判断基準を欠く。言い替えると、主観的評価は、オペレータに依存的である。さらに、ノイズ評価の手動方法は、一般に主観的評価のための自動プロセスと互換でもない。

ある種のノイズは画像全体を十分には表さない。たとえば、標準偏差（ＳＤ）などはある種の画像では信頼性が低い。これに関して、１つの例示的技法は、画像にラプラシアンを適用してエッジ画素を除去し、エッジを除いた画像から平均値または分散のいずれかによって計算される。ラプラシアン法は、原画像から直接的ではなく、画像処理を介してノイズを計算する。

他の従来技術の自動ノイズ査定技法では、ノイズが画像診断モダリティに依存すると仮定されるが、そのノイズモデルは複雑であった。対画像減算法（ｔｗｉｎ ｉｍａｇｅ ｓｕｂｔｒａｃｔｉｏｎ ｍｅｔｈｏｄ）などの他のより単純な従来技術の自動ノイズ査定技法は、ＣＴ収集シーケンスの奇数ビューと偶数ビューとからの２つの再構成を必要としたが、これは、画像ノイズが反復ごとに変化する反復処理についてこの技法を実現性のないものにする。

目的は、簡易な手法にして実質的にすべての画像診断モダリティについて人間の介入なしで信頼性の高いノイズインデックスを推定する画像処理装置を提供することにある。

本実施形態に係る画像ノイズ推定方法は、まず、画像内の複数の局所領域にそれぞれに対応する複数の局所ノイズインデックスを計算する。計算されたノイズインデックスからノイズレベルと度数との対応を表すヒストグラムが作成される。作成されたヒストグラムに基づいて画像に関するノイズインデックスが推定される。

また、本実施形態に係る画像処理装置は、画像内の複数の局所領域にそれぞれに対応する複数の局所ノイズインデックスを計算する計算部と、前記計算されたノイズインデックスからノイズレベルと度数との対応を表すヒストグラムを作成する作成部と、前記作成されたヒストグラムに基づいて前記画像に関するノイズインデックスとして前記局所領域の標準偏差又は擬似標準偏差を推定する推定部とを有する。

本発明の一実施形態によるマルチスライスＸ線ＣＴ装置の構成を示す図。 本実施形態によるノイズ推定デバイスによって実行される処理手順の全体的な流れを示す図。 本実施形態によるノイズ推定デバイスによって実行される処理手順の詳細を示す流れ図。 本実施形態において、標準偏差を計算する際に元の画像の平均値をとることによって元信号を平滑化するいくつかのステップを示す流れ図。 本実施形態において、心筋領域内の２つのＲＯＩ（関心領域）を有する肺ＣＴ画像を示す図。 本実施形態において、図４Ａの肺ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラム。 本実施形態において、３つのＲＯＩ（関心領域）を有する頭部ＣＴ画像。 本実施形態において、図５Ａの頭部ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラム。 本実施形態において、１つのＲＯＩ（関心領域）を有する例示的な腹部ＣＴ画像。 本実施形態において、図６Ａの腹部ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラム。 本実施形態において、頭部領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフ。 本実施形態において、腹部領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフ。 本実施形態において、心臓領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフ。 本実施形態において、画像全体の等しいサイズの領域に関する複数の中間ＰＳＤ値を示す図。 本実施形態において、サイズにおいて等しくない領域に関する複数の中間ＰＳＤ値を示す図。 本実施形態において、計算されるＰＳＤベースのノイズ推定値の利用という全般的な概念を示す図。 本実施形態において、画像（原画像）の一例を示す図。 本実施形態において、図１２の原画像に対する局所フィルタリングにより発生される局所平均値画像を示す図。 本実施形態において、図１３の局所平均値画像と原画像との差画像を示す図。 本実施形態において、図１４の差画像に関する二乗画像の一例を示す図。 本実施形態において、図１５の二乗画像に対する畳み込み演算により発生される畳み込み画像を示す図。 本実施形態において、図１６の畳み込み画像に関する平方根画像を示す図。

本実施形態は、画像に関する高信頼性のノイズインデックスを発生するために、ノイズインデックスの信頼性を低下させる領域内の画素の影響を抑制し、又はノイズインデックスの信頼性を高める領域（特徴領域）を対象としてノイズインデックスを決定する。特徴領域は、画像内に必ずしも明瞭に存在するわけではない。特徴領域は、局所の平均値またはその時間変化により探索される。例えば平均値は、画像内又は局所の画素値中央である。したがって、本実施形態は、ほとんどの種類の画像で実質的に信頼できるノイズインデックスを計算するために、最もノイズを知覚できる画素の中でノイズレベルを評価する方法とシステムとを対象とする。

本実施形態では、ノイズインデックスの計算は、システムノイズの低い組織領域を対象に実行されることが好ましい。組織領域は、一般に、ＣＴ画像又は超音波画像では軟部組織として比較的容易に特定される。ＣＴ画像では、軟部組織領域は、ＨＵ範囲（（−１０００ＨＵ〜１０００ＨＵ）の中央内（−１００ＨＵ〜１００ＨＵ）に対応する。その一方で、肺画像内の極小構造は、ノイズインデックスの計算上やっかいである。腫瘍および外傷は、一般に周囲の健全な組織に対して低いコントラストを有するので、診断の観点からはより狭い表示範囲または表示ウィンドウを要求する。ＣＴ画像は、画像内の組織領域を特定するのに都合がよい。超音波画像の場合には、統計的に安定したスペックルノイズが明確に展開される軟部組織領域内の画素を対象にノイズインデックスが計算される。

信頼性のあるノイズインデックスの計算は、ノイズフィルタの有効性評価に重要である。ノイズフィルタは、一般に、画像内のある種のノイズまたはアーチファクトを減らすために適用される。しかし、ノイズフィルタの有効性評価は、主観に左右されることが多い。他の状況では、ある種のノイズを減らすフィルタが、ターゲットノイズを実質的に減らすが、同一のノイズフィルタが、画像に別の種類のノイズまたはアーチファクトを生じさせることがある。これらの理由から、信頼できるノイズインデックスの計算は、ノイズを減らすフィルタの有効性を評価するために重要である。

信頼できるノイズインデックスは、反復ノイズ低減（ｉｔｅｒａｔｉｖｅ ｎｏｉｓｅ ｒｅｄｕｃｔｉｏｎ）プロセスの終了基準にも使われる。ノイズインデックスは、中央レベル輝度でゆっくり変化する平均値を有する画素（有効画素）を使って計算される。これらの有効画素は、一般的には、画像内の画素の大多数である。

本実施形態は、信頼性が高く実用的なノイズインデックスを計算（推定）するために、擬似標準偏差（ＰＳＤ：ｐｓｅｕｄｏ−ｓｔａｎｄａｒｄ ｄｅｖｉａｔｉｏｎ）とヒストグラム手法とを導入する。

ヒストグラム手法は、上記有効画素を識別する。エッジ、線、または小さいパッチなどのすばやく変化する領域と比較して、ゆっくり変化する領域は、通常は大きく、より多くの画素母集団を形成する。これらの要因は、ゆっくり変化する領域を分離するためのヒストグラム手法を正当化し、分離された領域がノイズレベル分析をうけることを可能にする。ヒストグラム手法において、頻度計数の区分（ビンサイズ）は、ノイズ画像に関する最終的なインデックス計算結果の頑健性および正確さに影響する。一般に、大きいビンは、ヒストグラム手法を安定させるが、より不正確にする。その一方で、小さいビンは、計算を正確にするが、より不安定にする。本実施形態では、当初は比較的小さいビンサイズが利用され、その後、ビン位置の精度を維持しながら移動平均を利用することによって、ビンサイズが拡大される。

一般に、不規則ノイズの影響を大きく受けたデータは、本実施形態においては、再構成画像データである。本実施形態では、再構成画像データを使って、ノイズインデックスとしてＰＳＤが計算される。画像データ内のノイズインデックス（ＰＳＤ）を計算する方法は、例えば測定データ内のノイズインデックス（ＰＳＤ）を計算する方法に実質的に適用可能である。

画像データの使用は、本実施形態のいくつかの態様によれば、測定されたデータに対するいくつかの実用的利益を有する。画像データは、検出器などの入力デバイスに関わらず、画素値の集合体という統一的な性質を有する。従って画像データは、測定データより一般的であり、多用途である。少なくともこれらの理由から、画像領域内の画像データの使用は、完全に自動的で機能的な単純なノイズ画像に関する最終的なインデックス計算方法を実施する際に有利である。

図１に示すように、本実施形態のマルチスライスＸ線ＣＴ装置は、ガントリ１００と他のデバイスまたはユニットとを有する。ガントリ１００は、Ｘ線管１０１と、環状フレーム１０２と、マルチロウタイプ（多列型）または２次元アレイタイプのＸ線検出器１０３とを備える。Ｘ線管１０１およびＸ線検出器１０３は、環状フレーム１０２上に被検体Ｓを挟んで対向に取り付けられる。環状フレーム１０２は、軸ＲＡの回りで回転する。回転ユニット１０７は、被検体Ｓが回転軸ＲＡに沿って移動される間に０．４秒／回転などの高速でフレーム１０２を連続的に回転させる。

マルチスライスＸ線ＣＴ装置は、電流調整器１１３と、Ｘ線管１０１からＸ線を発生させるためにＸ線管１０１にランプ電圧を印加する高電圧発生器１０９とを有する。例えば、高電圧発生器１０９は、フレーム１０２に取り付けられる。Ｘ線は、被検体Ｓに向かって放射される。被検体Ｓの断面は、円に近似的である。Ｘ線検出器１０３は、被検体Ｓを透過したＸ線を検出するために、被検体Ｓを介してＸ線管１０１に対向する。

Ｘ線検出器１０３からの検出信号を処理するための他のデバイスが設けられる。データ収集回路またはデータ収集システム（ＤＡＳ）１０４は、チャネルごとにＸ線検出器１０３からの出力信号を電圧信号に変換し、これを増幅し、さらに、これをディジタル信号に変換する。Ｘ線検出器１０３およびＤＡＳ１０４は、１回転あたりの膨大な投影データセットを処理するように構成される。

ＤＡＳ１０４の出力データは、非接触データ送信器１０５を介して、ガントリ１００の外部のコンソール内に収容された前処理デバイス１０６に送られる。前処理デバイス１０６は、生データに対して、感度補正などのある種の補正を実行する。その後、ストレージデバイス１１２は、再構成処理の直前段階にある投影データとも呼ばれるデータを格納する。ストレージデバイス１１２は、再構成デバイス１１４、表示デバイス１１６、入力デバイス１１５、およびスキャンプランサポート装置２００と一緒に、データ／制御バスを介してシステムコントローラ１１０に接続される。スキャンプランサポート装置２００は、スキャンプランを開発するために画像診断技師をサポートする機能を含む。

ノイズ推定デバイス１１７は、さまざまなソフトウェアモジュールとハードウェアコンポーネントとの組合せからなる。ノイズ推定デバイス１１７は、ノイズインデックスと擬似標準偏差（ＰＳＤ）とヒストグラムの関連機能を備える。ノイズ推定デバイス１１７は、データ／制御バスを介して再構成デバイス１１４とストレージデバイス１１２とに接続される。再構成デバイス１１４は、ストレージデバイス１１２に格納される投影データに基づいて、画像を再構成する。投影データは、検出器１０３の出力に基づいてデータ収集システム１０４と前処理デバイス１０６とにより発生される。画像データと測定データは画像信号と測定信号にそれぞれ同義である。単にデータ、信号というときは、それぞれ画像データと測定データを表している。

ノイズ推定デバイス１１７は、再構成デバイス１１４および／またはストレージデバイス１１２から再構成画像データを受け取る。再構成された画像のデータは、すでに画像領域にあるので、ノイズ推定デバイス１１７は、画像データを用いたノイズ画像に関する最終的なインデックス計算処理以外の追加処理を必要としない。

ノイズ推定デバイス１１７が再構成デバイス１１４および／または１１２から測定データを受け取るとき、測定データは再構成処理前のデータであるので、ノイズ推定デバイス１１７はノイズ画像に関する最終的なインデックス計算処理以外の追加処理を必要とするかもしれない。

ノイズ推定デバイス１１７は、実質的に人間の介入なしで、ノイズインデックスとしてのＰＳＤを計算する際に、受け取られた信号に対してタスクの所定のセットを実行する。自動プロセスまたは自動方法は、ヒストグラム分析のビンサイズおよび望まれない背景データを除去するためのしきい値などのパラメータの所定の条件データセットを要する。さらに、パラメータは、ノイズインデックスとしてＰＳＤを計算する際に信号またはデータをフィルタリングの処理単位（局所）を規定するカーネルのサイズまたは特性も定義する。

図２には、本実施形態によるノイズ推定デバイス１１７の処理手順が示されている。ステップＳ１０では、元の信号Ｄ１は、所定の一様カーネルＫＩを畳み込まれる。それにより元の信号Ｄ１は、平滑化される。原画像Ｄ１を構成する画素各々を中心とした局所にカーネルＫＩが適用される。局所内の画素の画素値にカーネルＫＩの構成要素が乗算され、合計される。この計算結果が、当該局所の中心画素の元の画素値から置き換えられる。全ての画素について同様の局所処理が実行される。元の信号Ｄ１は、測定データまたは画像データのいずれかである。カーネルＫＩは、典型的にはＮ×Ｎなどの正方形カーネルである。カーネルＫＩの他の例は、Ｎ×Ｍなどの長方形カーネルである。ＮとＭはともに整数である。ステップＳ１０での畳込みにより画像は平滑化される。ステップＳ２０では、ステップＳ１０による平滑化信号（一次平滑化画像）と元の信号（原画像）Ｄ１との間の差が計算される。平滑化画像を構成する複数の画素各々に対して原画像Ｄ１上の同位置にある画素との間で、画素値の差が計算される。この計算結果は、必要に応じて「差画像」という。

次にノイズ推定デバイス１１７は、ステップＳ３０を実行する。ステップＳ３０では、ステップＳ２０で計算された差信号が、二乗される。画像であれば、差画像を構成する画素の画素値がそれぞれ二乗される。この計算結果は、必要に応じて「二乗画像」という。

ステップＳ４０は、２回目の平滑化処理である。ステップＳ４０では、ステップＳ３０の二乗結果（二乗画像）に対して、所定の一様カーネルＫＩＩを用いて畳み込むことによってもう一度平滑化される。二乗画像を構成する画素各々を中心とした局所にカーネルＫＩＩが適用される。局所内の画素の画素値にカーネルＫＩＩの構成要素が乗算され、合計される。この計算結果が、当該局所の中心画素の元の画素値から置き換えられる。全ての画素について同様の局所処理が実行される。２回目の平滑化処理により、二次平滑化画像が発生される。カーネルＫＩＩは、典型的にはＰ×Ｐなどの正方形カーネルである。カーネルＫＩＩの他の例は、Ｐ×Ｑなどの長方形カーネルである。ＰとＱはともに整数である。カーネルＫＩとＫＩＩは、典型的には相違するが、同一であることは否定されない。

ステップＳ５０では、ステップＳ４０からの二次平滑化画像を構成する画素各々の画素値の平方根が計算される。これら一連Ｓ１０−Ｓ５０の計算処理により、ノイズインデックスとしての疑似標準偏差（ＰＳＤ）が計算される。計算結果は、全体としてＰＳＤ分布を示す。

ヒストグラムを生成するために、ＰＳＤ分布を構成する複数のＰＳＤは、例えば昇順で並び替えられる。ＰＳＤデータは、１次元（１Ｄ）ベクトルデータとして記録される（Ｄ３）。次に、ＰＳＤ全範囲を均等に分割された複数のヒストグラムビン（複数のブラケット、複数の区分ともいう）各々に含まれるＰＳＤの個数が計数される。垂直軸は、ＰＳＤの個数または出現度数であり、水平軸は、ヒストグラムビン列、つまり各ヒストグラムビンのＰＳＤ範囲又は中心値の配列軸である。

ステップＳ６０では、ヒストグラムの対象から画像の背景領域を示す画素が除去される。背景画素は、他の有用な組織画素と比較して、ノイズインデックスの変動に対する影響は大きい。言い替えると、これらの望ましくない背景画素の選択的除去の失敗は、ヒストグラム内のいくつかの誤ったピークをもたらす可能性があり、つまり画像に関する最終的なノイズインデックスの信頼性を、臨床的に関連しない背景を反映するこれらの誤ったピークにより低下させるおそれがある。最後に、ステップＳ７０では、画像に関する最終的なノイズインデックスとしてのＰＳＤベースのノイズレベルが、ヒストグラムに基づく分布モードとして決定される。

ノイズインデックス推定方法の上述の方法は、ゆっくり変化する領域内の画素母集団が、エッジなどのすばやく変化する領域内の画素母集団よりノイズインデックスの信頼性向上に有用であると仮定するものである。分布モードを選択することによって、最も一般的なまたは頻繁に発生する擬似標準偏差（ＰＳＤ）の値が、画像全体または信号全体の実質的に真のノイズレベルを表す単一のノイズインデックスとして選択される。

図３Ａには、本実施形態によるノイズ推定デバイス１１７によって実行されるノイズ推定処理手順の詳細が示されている。一般に、ノイズは、光子計数での量子ノイズと、データ収集システム（ＤＡＳ）１０４での電子ノイズと、アナログ−ディジタル変換（Ａ／Ｄ）での量子化ノイズと、再構成処理の近似誤差と、その他の発生源に起因して生じる。Ｉ（ｘ，ｙ）が、ある加法性ノイズによって壊された画像を表すと仮定すると、画像データＩ（ｘ，ｙ）は、次の式（１）によって表される。

ここで、Ｉ_０は、ノイズのない真値画像であり、ｎは様々なノイズ発生源に起因して真値に加算されたノイズ成分である。ｘおよびｙは、画素の２次元座標である。図３Ａの流れ図に含まれるステップは、所定のノイズインデックスＰＳＤを計算するために、式（２）で要約される。

ここで、

は、畳込み演算子を表し、ｗ（ｘ，ｙ）は、正規化された一様な移動平均カーネルである。ブラケットの内部の移動平均は、各画素近傍の平均値を計算するが、ブラケットの外部の移動平均は、平均二乗誤差（ＭＳＥ）の近傍平均値を計算する。ＭＳＥが、標準偏差の公式が普通に行うようにサンプルの信号平均値を減算するのではなく、フィルタリングされたサンプルを減算することによって計算されることに留意されたい。ＰＳＤは、（１）の仮定が有効である時に「真の」標準偏差に接近する。ＰＳＤは、標準偏差（ＳＤ）の直接計算の計算効率を改善するために使用される。

ヒストグラム手法は、領域を分離し、ノイズレベルを分析するのに利用されてきた。この手法は、ゆっくり変化する領域が、通常は大きく、エッジ、線、または小さいパッチなどのすばやく変化する領域より大幅に大きい画素母集団を有するという、ほとんどの画像または信号での一般的な発見に基づく。ヒストグラム手法を用いると、ビンサイズは、ノイズ画像に関する最終的なインデックス計算結果の頑健性と正確さとにとって重要である。上述したように、ビンサイズは、ノイズ画像に関する最終的なインデックス計算の安定性と正確さとのトレードオフを有するので、比較的小さいビンサイズが、当初にヒストグラムを入手するのに使用され、移動平均が、その後、ビン位置の精度を保ちながらビンサイズを拡大するために適用される。これに関して、ｈ（ｕ）がヒストグラムであり、ｍ（ｕ）が、１Ｄ正規化された移動平均カーネルであると仮定すると、フィルタリングされたヒストグラムＨ（ｕ）は、下の式（３）のように畳込みによって定義される。

ここで、ｕは、ＰＤＳ値内のビン位置である。移動平均ウィンドウの長さ（Ｌ）は、ヒストグラム内の最大値の位置（Ｕmax）の単調増加関数である。

長さＬは、次の式（４）を用いて計算される。

なお、Ｌは｛Ｕ_max／Ｋ｝を超える最も近い奇数の整数である。ここで、係数ｋは、データノイズレベルの関数でもある。Ｕmaxが１０未満である時に、ｋ＝５は、Ｕmaxが９０を超える間に安定した解を保証し、ｋ＝３は、安定した解を保証する。その一方で、Ｕmaxが１０と９０との間である時には、ｋは、次の単調減少する式（５）を使用して選択されなければならない。

ｋの選択は、比較的静かなヒストグラムを軽いカーネルによってまたは平滑化カーネルなしで処理しながら、ノイズの多いヒストグラムをより重い平滑化カーネルによって処理することを可能にする。

さらに、正確なノイズ画像に関する最終的なインデックス計算を保証するために、上の手順の複数のパスが、オプションで、本実施形態のある種の例示的な方法とシステムとで実行される。一般に、２パスが十分である。第１のパスでは、第１のヒストグラムピークでのカウントが、平滑化マスクｗ（ｕ，ｖ）の最良のサイズを判断するのに使用される。より大きいマスクが、大きいカウントについて選択され、逆も同様である。マスクサイズに関する経験的に導き出されたセットの一例は次の通りである。
カウント＜５０００、マスクサイズ＝７×７
カウント＞１００００、マスクサイズ＝１１×１１
さらに図３Ａを参照すると、ステップＳ８０では、元の画像Ｉが、その画像の移動平均をとることによって平滑化される。すなわち、所定のサイズの正規化された一様移動平均カーネルが、適用されている。破線は、近傍の外の画像の広がりを示し、実線は、中心画素及びその近傍の画素とからなる局所画素群を示す。これに関して、この例でのＤ１の局所は、９画素を含む３×３近傍である。画像Ｉにわたって移動平均をとることによって、新しい画像Ｉ_ＭＡが作成される。この画像Ｉ_ＭＡを構成する画素各々の値は、その３×３近傍の平均値である。特定の近傍領域Ｄ１について、対応する移動平均結果は、Ｄ４である。図２に関して既に述べたように、ステップＳ８０での正規化された一様カーネルを用いる畳込み処理の結果は、ここでは、各画素での画像の近傍平均値である。ステップＳ２０では、マイナス記号によって示されるように、ステップＳ８０からのＤ４内の値のうちの対応するものと元の信号Ｄ１との間の差が、計算される。さらに、ステップＳ３０では、差を二乗する。ステップＳ８０Ａでは、二乗された差の近傍平均が、正規化された一様移動平均カーネルを適用することによって計算される。最後に、ステップＳ５０で、ステップＳ８０Ａからの結果の平方根が計算される。したがって、画像内の各画素の擬似標準偏差（ＰＳＤ）が計算された。ステップＳ８０およびＳ８０Ａでの正規化された一様移動平均カーネルは、オプションで、１つの方法では同一であり、もう１つの方法では異なる。

図３Ｂを参照すると、画像Ｉは、複数のブロックに分割される。ブロックサイズはＤ１に等価である。ステップＳ９０では、各ブロック内の画素平均値が計算される。平均値は、Ｄ５内のすべての３×３近傍画素に割り当てられる。ステップ１００では、ステップＳ９０からの平均値と元の信号Ｄ１との間の差（誤差）が計算される。さらに、ステップＳ１１０では、誤差が二乗される。ステップＳ９０Ａでは、二乗された誤差のブロック平均が計算される。最後に、ステップＳ１２０で、ステップＳ９０Ａからの結果の平方根が、計算される。つまり局所領域ごとに通常の標準偏差（ＳＤ）が計算されている。

通常のＳＤに対してＰＳＤは利点がある。ＰＳＤは標準偏差の計算処理対象からエッジ変動要因を除外することによって、ＳＤが本質的に備えるエッジ変動による信頼性の低下を抑制できる。ＰＳＤは、画像が一定の平均値またはゆっくり変化する平均値を有する領域内で、ＳＤに近づく。

ファントムデータおよび患者データを含む多数のデータセットに上述のノイズ推定方法を適用して推定されたノイズインデックスは、注意深い手動操作を介入させて計算された標準偏差に等価な結果を示す。手動評価では、手動で選択された最良の関心領域（ＲＯＩ）が、究極の判断基準として使用される。上で説明したように、次の例示的画像は、再構成されたＣＴ画像であるが、本実施形態によるノイズ推定の方法および実施形態は、再構成されたＣＴ画像に限定されず、ＣＴからの他の測定されたデータおよび他のモダリティに適用可能である。

ここで図４Ａおよび４Ｂを参照すると、本実施形態によるＰＳＤベースの汎用ノイズ推定値を計算する上で説明された例示的方法が、肺ＣＴ画像に適用されている。図４Ａは、心筋領域内の白線によって示される２つの注意深く選択されたＲＯＩ（関心領域）を有する例示的なズームされた肺ＣＴ画像である。手動ノイズ画像に関する最終的なインデックス計算が、ＲＯＩ１１とＲＯＩ１２とを含む２つのＲＯＩ内で行われた。測定されたノイズの標準偏差（ＳＤ）は、ＲＯＩ１１内では３３．２５ＨＵであり、ＲＯＩ１２内では３３．７９ＨＵであった。これらのＳＤ値は、本実施形態による上で説明された例示的方法に基づく同一の肺画像のＰＳＤのノイズインデックスと比較される。

ここで図４Ｂを参照すると、本実施形態の上で説明された例示的方法による、図４Ａの肺ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラムが、生成される。このヒストグラムは、肺ＣＴ画像内の画素ごとに計算されたＰＳＤから生成される。３３ＨＵの弱いピークが、最も頻繁に発生するノイズ値を表すので、このピークが単一の推定されたノイズインデックスとして自動的に選択される。自動的に計算されるノイズレベルは、図４ＡのＲＯＩ１１とＲＯＩ１２との手動測定された結果と実質的に一致する３３ＨＵである。

ここで図５Ａおよび５Ｂを参照すると、本実施形態によるＰＳＤベースの汎用ノイズ推定値を計算する上で説明された例示的方法が、肺ＣＴ画像に適用されている。図５Ａは、白線によって示される３つの注意深く選択されたＲＯＩ（関心領域）を有する例示的な頭部ＣＴ画像である。手動ノイズ画像に関する最終的なインデックス計算が、ＲＯＩ１１とＲＯＩ１２とＲＯＩ１３とを含む３つのＲＯＩ内で行われた。測定されたノイズの標準偏差（ＳＤ）は、ＲＯＩ１１内では２１．４９ＨＵであり、ＲＯＩ１２内では２１．７６ＨＵであり、ＲＯＩ１３内では２３．７０ＨＵであった。これらのＳＤ値は、本実施形態による上で説明された例示的方法に基づく同一の頭部画像のＰＳＤのノイズインデックスと比較される。

ここで図５Ｂを参照すると、本実施形態の上で説明された例示的方法による、図５Ａの頭部ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラムが、生成される。このヒストグラムは、頭部ＣＴ画像内の画素ごとに計算されたＰＳＤから生成される。２２ＨＵの弱いピークが、最も頻繁に発生するノイズ値を表すので、このピークが単一の推定されたノイズインデックスとして自動的に選択される。自動的に計算されるノイズレベルは、図５ＡのＲＯＩ１１とＲＯＩ１１とＲＯＩ１３との手動測定された結果と実質的に一致する２２ＨＵである。

ここで図６Ａおよび６Ｂを参照すると、本実施形態によるＰＳＤベースの汎用ノイズ推定値を計算する上で説明された例示的方法が、腹部ＣＴ画像に適用されている。図６Ａは、白線によって示される１つの注意深く選択されたＲＯＩ（関心領域）を有する例示的な腹部ＣＴ画像である。手動ノイズ画像に関する最終的なインデックス計算が、ＲＯＩ内で行われた。測定されたノイズの標準偏差（ＳＤ）は、１１．７７ＨＵであった。このＳＤ値は、本実施形態による上で説明された例示的方法に基づく同一の腹部画像のＰＳＤのノイズインデックスと比較される。

ここで図６Ｂを参照すると、本実施形態の上で説明された例示的方法による、図６Ａの腹部ＣＴ画像を利用するＰＳＤに基づいて単一のノイズインデックスを計算するためのヒストグラムが、生成される。このヒストグラムは、腹部ＣＴ画像内の画素ごとに計算されたＰＳＤから生成される。１３ＨＵの弱いピークが、最も頻繁に発生するノイズ値を表すので、このピークが単一の推定されたノイズインデックスとして自動的に選択される。自動的に計算されるノイズレベルは、図６ＡのＲＯＩの手動測定された結果と実質的に一致する１３ＨＵである。

上で説明したＲＯＩ手法に加えて、ＰＳＤ値は、ＣＴ画像内の異なる身体領域内の体積全体を通じてスライス上の対応する標準偏差（ＳＤ）と比較される。図７は、頭部領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフである。ｘ軸は、スライス番号を示し、ｙ軸は、ＰＳＤ値またはＳＤ値のいずれかを示す。実線は、対応するスライス上の自動計算されたＰＳＤ値を示し、正方形を伴う実線は、手動計算されたＳＤ値を示す。体積全体の頭部内のスライス全体を通じて、自動計算されたＰＳＤ値は、手動計算されたＳＤ値と実質的に一致する。

図８は、腹部領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフである。ｘ軸は、スライス番号を示し、ｙ軸は、ＰＳＤ値またはＳＤ値のいずれかを示す。実線は、対応するスライス上の自動計算されたＰＳＤ値を示し、正方形を伴う実線は、手動計算されたＳＤ値を示す。体積全体の腹部内のスライス全体を通じて、自動計算されたＰＳＤ値は、手動計算されたＳＤ値と実質的に一致する。

図９は、心臓領域内のＣＴスライスにまたがるＰＳＤ値とＳＤ値との比較を示すグラフである。ｘ軸は、スライス番号を示し、ｙ軸は、ＰＳＤ値またはＳＤ値のいずれかを示す。実線は、対応するスライス上の自動計算されたＰＳＤ値を示し、正方形を伴う実線は、手動計算されたＳＤ値を示す。体積全体の心臓内のスライス全体を通じて、自動計算されたＰＳＤ値は、手動計算されたＳＤ値と実質的に一致する。

単一のノイズ推定インデックスを計算する上で説明された方法では、最も代表的なＰＳＤ値が、本実施形態による例示的な方法またはシステムで、画像全体、体積全体、またはデータの完全なセットについて自動的に選択される。単一選択手法を保護するために、本実施形態によるもう１つの例示的な方法またはシステムでは、画像全体、体積全体、またはデータの完全なセットではなく、より特定の関心領域または関心区域内で、複数のノイズ推定インデックス値がＰＳＤ値から計算される。

ここで図１０Ａおよび１０Ｂを参照すると、図に、本実施形態によるもう１つの例示的な方法またはシステムでより特定の関心領域または関心区域内でＰＳＤから選択される複数のノイズ推定インデックス値が示されている。図１０Ａに、画像全体の等しいサイズの領域の複数の中間ＰＳＤ値を示す。画像画素区域Ｉ１全体が、９つの等しいサイズの部分Ａ１からＡ９に分割される。ＰＳＤ値は、最初に、本実施形態による上で説明された方法とシステムとに従って画像Ｉ１内の画素ごとに計算される。その後、単一の最も代表的なＰＳＤ値が、画像Ｉ１全体ではなく所定の等しいサイズの部分Ａ１からＡ９のそれぞれで自動的に選択される。複数のＰＳＤ値は、部分Ａ１からＡ９のそれぞれの実質的に局所的なノイズ特性を保持するので、選択されたＰＳＤノイズ推定値のそれぞれが、オプションで、ある種の状況で計算される。複数の中間ノイズ推定値は、背景ノイズを手動で除去する必要を有利になくすことができる。有利なことに、複数の中間ノイズ推定値は、エッジなどのすばやく変化する画素母集団を有力に含む画像内にある可能性がある。

ここで図１０Ｂを参照すると、図に、サイズにおいて等しくない領域の複数の中間ＰＳＤ値が示されている。画像画素区域Ｉ２全体が、所定の個数の等しくないサイズの部分Ｂ１とＢ２とに分割される。たとえば、等しくないサイズの部分Ｂ１とＢ２とは、それぞれ、１つの関心領域（ＲＯＩ）であり、それぞれが医療的に関心を持たれる特定の構造を含む。ＰＳＤ値は、最初に、本実施形態による上で説明された方法とシステムとに従って画像Ｉ２内の画素ごとに計算される。その後、単一の最も代表的なＰＳＤ値が、画像Ｉ２全体ではなく所定の等しくないサイズの部分Ｂ１とＢ２とのそれぞれで自動的に選択される。複数のＰＳＤ値は、医療的に重要な区域Ｂ１とＢ２とのそれぞれの実質的に局所的なノイズ特性を保持するので、選択されたＰＳＤノイズ推定値のそれぞれが、オプションで、ある種の状況で計算される。複数の中間ノイズ推定値は、背景ノイズを手動で除去する必要を有利になくすことができる。有利なことに、複数の中間ノイズ推定値は、エッジなどのすばやく変化する画素母集団を有力に含む画像内にある可能性がある。

ここで図１１を参照すると、図に、本実施形態の上で説明された例示的な方法とシステムとによって計算されるＰＳＤベースのノイズ推定値の利用という全般的な概念が示されている。この概念は、ＰＳＤ値が本実施形態に従って画像全体の画素のそれぞれについて計算された後に、それらのＰＳＤ値を処理する流れで示されてもいる。ステップＳ１００では、ＰＳＤ値が、所定の画素近傍サイズに従って画素のそれぞれについて計算される。これらのＰＳＤ値のそれぞれは、実質的に局所データに敏感なので、これらのＰＳＤ値を、ステップＳ１００に示されているように局所ＰＳＤ値と称する。

後続のステップＳ１１０では、中間ＰＳＤ値が、オプションで、ステップＳ１００の局所ＰＳＤ値に基づいて計算される。図１０Ａと１０Ｂとに関して説明したように、中間ＰＳＤ値は、ルールの所定のセットに従って計算される。例で示したように、中間ＰＳＤ値は、オプションで、画像の等しいサイズまたは等しくないサイズの部分について選択される。どの場合でも、中間ＰＳＤ値は、画像全体、体積全体、スライス全体、またはデータ全体の単一の代表的なノイズ推定値ではない。中間ＰＳＤ値の計算を実施する際には、複数のヒストグラムが、本実施形態に従って、画像全体、体積全体、スライス全体、またはデータ全体の諸部分または諸区域について生成される。複数のＰＳＤベースの値は、上の判断基準に限定されず、オプションで、測定された信号の角度依存性などの他のデータ特性に基づいて計算される。

さらに図１１を参照すると、単一の代表的なＰＳＤ値が、オプションで、ステップＳ１００の局所ＰＳＤ値および／またはステップＳ１１０の中間ＰＳＤ値に基づいて計算される。図１０Ａと１０Ｂとに関して説明したように、中間ＰＳＤ値は、ルールの所定のセットに従って計算される。例で示したように、複数の中間ＰＳＤ値が、オプションで、画像の等しいサイズまたは等しくないサイズの部分について選択される。対照的に、単一の代表的なＰＳＤ値は、本実施形態のもう１つの実施形態に従って、画像全体、体積全体、スライス全体、またはデータ全体に広域的なＰＳＤベースのノイズインデックスとして割り当てられる。広域的なＰＳＤベースのノイズインデックスは、オプションで、中間ＰＳＤ値と組み合わせて計算される。代替案では、広域的なＰＳＤベースのノイズインデックスは、オプションで、単一のノイズ推定インデックスとして計算される。

上記式（２）について具体的に説明する。ここでは図１２に示す画像（原画像）を例示する。原画像Ｉ(x,y)の局所にフィルタｗ(x,y)が畳み込まれる。フィルタｗ(x,y)は、３×３画素のサイズであり、各要素は１／９である。フィルタされた画像の各画素は、３×３画素の局所内の平均値を示す。図１２に破線で示す局所の平均値は、７．７であり、図１３に示すように当該平均値が当該局所の中心画素（ハッチングされた画素）の元値から置き換えられる。図１４に示すように、局所平均値と、当該中心画素の元値との差が計算される。差は、式（２）で[]内の値を示す。

当該差の二乗による二乗画像を図１５に示す。二乗画像をフィルタｗ(x、y)で畳み込み処理をすると、図１６に示す画像が得られる。その平方根としてのＰＳＤ分布は図１７に示している。

図１２−図１７の例では、組織境界がＰＳＤ計算単位（破線で示す局所）を横断している。仮に通常の手順でＳＤを計算すると２．３５を示す。これは境界もノイズと同様に誤差として貢献するからである。しかしＰＳＤは画素毎に周りの平均値を求めているので、境界の影響をより受け難いと言う特徴がある。

ＰＳＤ以外のノイズフィッティングを求めても良い。例えばＰＳＤの代わりにＳＤを使用してもよい。ＳＤを使用することにより、組織境界でもＳＤ値が大きくなるが、ヒストグラムを計算すると、組織境界は画像全体の中で一部分しかないので、ＰＳＤによる結果とＳＤによる結果はほぼ同じになる。あるいは画像に対しノイズ低減フィルターをかける、高次の関数でフィッティングするなどしてノイズ分布を求めても良い。ノイズ低減フィルターは、拡散フィルター(Diffusion Filter)、コヒーレントフィルター(Coherent Filter)、中央値フィルター(Median Filter)、ガウスフィルター(Gauss Filter)、平均値フィルターなどから任意に選択される。フィッティング方法としては、例えばＢ−ｓｐｌｉｎｅでフィッティングする、３次元関数でフィッティングするなどし任意の種類のフィッティングが用いられる。このようにしてノイズ分布を元にＳＤを計算することで局所のノイズインデックスを求めても良い。

しかし、本発明の構造と機能との詳細と一緒に、本発明の多数の特性と特徴とが、前述の説明で示されたが、本開示が例示にすぎないことと、詳細において、特に部分の形状、サイズ、および配置の問題とソフトウェア、ハードウェア、または両方の組合せでの実施において、変更を行うことができるが、その変更が、添付の特許請求の範囲がその中で表される用語の広義の一般的な意味によって示される完全な範囲までの本発明の原理に含まれることを理解されたい。

Claims (14)

1. 画像内の複数の局所領域にそれぞれに対応する複数の局所ノイズインデックスを計算する計算部と、
   前記計算されたノイズインデックスからノイズレベルと度数との対応を表すヒストグラムを作成する作成部と、
   前記作成されたヒストグラムに基づいて前記画像に関するノイズインデックスとして前記局所領域の標準偏差又は擬似標準偏差を推定する推定部とを有することを特徴とする画像処理装置。
2. 前記計算部は、前記画像から平滑化画像を発生することを特徴とする請求項１記載の画像処理装置。
3. 前記擬似標準偏差は、前記画像と前記平滑化画像との差画像に対して画素毎に二乗処理、平滑化処理、平方根処理を順にかけることにより発生されることを特徴とする請求項２記載の画像処理装置。
4. 前記作成されたヒストグラムに基づいて前記局所領域のサイズの良否を決定する決定部をさらに有することを特徴とする請求項１記載の画像処理装置。
5. 前記複数の局所領域は同サイズを有する請求項１記載の画像処理装置。
6. 前記複数の局所領域は異なるサイズを有する請求項１記載の画像処理装置。
7. 前記複数の局所領域は臨床的に重要な構造を含む請求項１記載の画像処理装置。
8. 前記計算部は、前記画像をノイズ低減画像に変換し、前記ノイズ低減画像と前記画像とを引き算してノイズ画像を発生し、前記ノイズ画像の標準偏差を計算することを特徴とする請求項１記載の画像処理装置。
9. 前記ノイズ低減画像は、前記画像にノイズ低減フィルタをかけることにより発生される請求項１記載の画像処理装置。
10. 前記ノイズ低減画像は、前記画像を特定の関数にフィッティングすることにより発生される請求項１記載の画像処理装置。
11. 前記推定部は、前記作成されたヒストグラムから前記画像の背景部分に対応する度数を除外する請求項１記載の画像処理装置。
12. 前記背景部分は、視野外の部分と骨に対応する部分を含む請求項１１記載の画像処理装置。
13. 前記作成部は、前記ヒストグラムの最高度数に基づいて前記局所領域のサイズが調整される請求項１記載の画像処理装置。
14. 前記画像は３次元画像であり、
    前記複数の局所ノイズインデックスは、前記３次元画像を構成する複数の２次元画像に対して計算され、
    前記ヒストグラムは、前記複数の２次元画像に対して計算された全ての前記複数の局所ノイズインデックスを対象として作成され、
    前記３次元画像に関するノイズインデックスが、前記作成されたヒストグラムに基づいて推定される請求項１記載の画像処理装置。