JP2016198469A - 画像処理装置、画像処理システム、画像処理方法、及びプログラム - Google Patents

Abstract

【課題】散乱線の分布が変化する境界も含めた散乱線像を精度良く推定する。【解決手段】画像処理装置１０５は、放射線源から被写体に放射線を照射して撮影された放射線画像を取得し、放射線画像に含まれる散乱線像を推定し放射線画像から散乱線像が低減された補正画像を出力する。散乱線像の推定は、推定モジュール２２０において、第一の関数の出力についての特定の周波数成分と、第一の関数とは異なる第二の関数の出力についての特定の周波数成分と、に基づき散乱線成分を得る散乱線モデルを放射線画像に適用することにより行う。【選択図】図２

Description

本明細書の開示は、放射線画像に画像処理を施す画像処理装置、及び画像処理プログラムに関する。

被写体に放射線を照射することにより得られる放射線画像には、放射線源から直進する一次放射線による成分のほかに、被写体内で散乱された放射線である散乱線による成分が含まれる。かかる散乱線の成分は放射線画像のコントラストを低下させるおそれがある。
Ｘ線画像の散乱線成分を低減する技術として、特許文献１には、一次Ｘ線で表される１つの式に対して２つのガウス関数をコンボリューションした式を用いて散乱線成分をモデル化し、１次Ｘ線画像を得ることが開示されている。

特開２０１０−１８８１１３号公報

ところで、図９（ａ）は胸部ファントムを用いて測定し、実験的に放射線画像に含まれる散乱線成分を取得した例を示す図である。図９（ａ）によると、散乱線成分には、被写体を透過せず被写体の周辺を通過したＸ線が入射した領域（素抜け領域）９０１と被写体領域９０２との境界に、いわゆるスキンライン９０３が存在する。これは、被写体の厚みによって、被写体に入射した放射線が散乱する態様が異なり、素抜け領域では散乱しないため、散乱Ｘ線の強度分布が不連続となるためである。

本発明の発明者らは、特許文献１に開示の、上述した一次Ｘ線で表される１つの式は被写体内で一回散乱した一次散乱Ｘ線の強度を表すことを見出し、上述した一次Ｘ線で表される１つの式だけでは、図９（ａ）に例示されるような散乱Ｘ線の強度分布を正確に表現するには不十分であることを見出した。本発明の発明者らは、上記の課題に鑑み鋭意検討した結果、以下の構成により散乱線成分をより正確に推定できることを見出した。

本発明に係る画像処理装置は、被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、前記被写体内で前記放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、前記被写体内で前記放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づいて、前記画像取得手段により取得された放射線画像に含まれる散乱線成分を推定する推定手段と、前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力手段と、を有することを特徴とする。

本発明によれば、被写体内で放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、被写体内で放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づくことにより、散乱線の強度分布をより正確に推定できる。

本発明の実施形態に係る画像処理装置を含む情報システムの構成を示す図である。 本発明の実施形態に係る画像処理装置の構成を示す図である。 本発明の実施形態に係る処理を示すフローチャートである。 本発明の実施形態に係る第一の関数で近似する散乱Ｘ線を例示する図である。 本発明の実施形態に係る第二の関数で近似する散乱Ｘ線を例示する図である。 本発明の別の実施形態に係る処理を示すフローチャートである。 本発明の別の実施形態に係るラプラシアンピラミッドによる周波数分解および画像再構成を示す図である。 本発明の別の実施形態に係るラプラシアンピラミッドによる周波数分解の周波数特性を示す図である。 ファントムを用いて測定した放射線画像に含まれる散乱線成分の一例を示す図である。

図１に基づいて本発明の実施形態に係る画像処理装置と、放射線撮影システム１００を含む情報システム１２０について説明する。一実施形態に係る画像処理装置は、放射線撮影システム１００に含まれる制御部１０５に対応する。ここで言う情報システム１２０とは、たとえばＨＩＳ（Ｈｏｓｐｉｔａｌ Ｉｎｆｏｒｍａｔｉｏｎ Ｓｙｓｔｅｍ）１１１と、ＲＩＳ（Ｒａｄｉｏｇｒａｐｈｙ Ｉｎｆｏｒｍａｔｉｏｎ Ｓｙｓｔｅｍ）１１２と、ＷＳ（Ｗｏｒｋ Ｓｔａｔｉｏｎ）１１３と、ＰＡＣＳ（Ｐｉｃｔｕｒｅ Ａｒｃｈｉｖｉｎｇ ａｎｄ Ｃｏｍｍｕｎｉｃａｔｉｏｎ Ｓｙｓｔｅｍ）１１４と、Ｖｉｅｗｅｒ１１５と、Ｐｒｉｎｔｅｒ１１６とを有する。ＨＩＳ１１１は患者情報や放射線撮影による検査等を含む診療情報を総合的に管理するシステムである。ＲＩＳ１１２は放射線撮影のオーダを管理するシステムである。ＷＳ１１３は画像処理端末であり、放射線撮影システム１００で撮影された放射線画像に画像処理を施す。ＷＳは同様の機能を有するソフトウェアをインストールした一又は複数のコンピュータで代替してもよい。ＰＡＣＳ１１４は当該情報システム１２０内の放射線撮影やその他の医用画像撮影装置で得られた画像を保持するデータベースシステムである。ＰＡＣＳ１１４は医用画像及びかかる医用画像の撮影条件や患者情報等の付帯情報を記憶する記憶部（不図示）と、当該記憶部に記憶される情報を管理するコントローラ（不図示）とを有する。Ｖｉｅｗｅｒ１１５は、画像診断用の端末であり、ＰＡＣＳ１１４等に記憶された画像を読み出し、診断のために表示する。Ｐｒｉｎｔｅｒ１１６はたとえばフィルムプリンタであり、ＰＡＣＳ１１４等に記憶された画像をフィルムに出力する。

実施形態における放射線撮影システム１００では、放射線としてＸ線を用いる。放射線撮影システム１００は、放射線発生装置の例であるＸ線源１０１と、ＦＰＤ（Ｆｌａｔ Ｐａｎｅｌ Ｄｅｔｅｃｔｏｒ）１０２と、制御部１０５とを有する。これらはケーブルもしくは通信システムを介して接続されている。制御部１０５は撮影された放射線画像に、当該撮影の撮影条件や患者情報等を付帯させる。たとえばＤＩＣＯＭ（Ｄｉｇｉｔａｌ Ｉｍａｇｉｎｇ ａｎｄ Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ ｉｎ Ｍｅｄｉｃｉｎｅ）規格に則って情報を付帯させ、放射線画像のデータ、患者情報、及び撮影条件等の情報を含むＤＩＣＯＭ画像ファイルを生成する。制御部１０５は当該画像をＷＳ１１３やＰＡＣＳ１１４に送信する。当該撮影のオーダはたとえばＲＩＳ１１２から制御部１０５に送信される。制御部１０５はＲＩＳ１１２からの入力情報に応じて撮影条件を記憶部（不図示）から読み出す。

Ｘ線源１０１は、Ｘ線管、又は医用画像もしくは他の画像の取得に適した他の任意の放射線源であってもよい。高圧発生部１０４は、操作者が曝射スイッチを押下するとＸ線源１０１に高圧パルスを与え、Ｘ線源１０１から被写体１０３が配置されている領域にＸ線を曝射させる。被写体１０３を透過もしくは周囲を通過したＸ線は、Ｘ線検出器であるＦＰＤ１０２に入射する。ＦＰＤ１０２は制御部１０５に制御され、入射したＸ線を電気信号に変換した後、デジタル画像として制御部１０５に送信する。たとえば、ＦＰＤ１０２は入射したＸ線を蛍光体（不図示）が可視光に変換し、可視光をフォトダイオード（不図示）が検出し、Ａ／Ｄ変換器（不図示）により電気信号に変換する。あるいは、ＦＰＤ１０２はＸ線をアモルファスセレン（不図示）により電気信号に変換する。放射線画像の画素値はＦＰＤ１０２を構成する放射線検出素子１０２ａからの出力により得られる。放射線検出素子１０２ａは、例えば蛍光体（不図示）とフォトダイオード（不図示）で構成される。別の例では、アモルファスセレン（不図示）で構成される。

デジタル画像は制御部１０５やＷＳ１１３で画像処理が施され、ＰＡＣＳ１１４等に保存される。情報システム１２０に含まれる各部はバスやその他通信システムにより相互に接続されていればよく、それぞれを遠隔に設置することもできる。

続いて、図２に基づいて本発明の実施形態に係る画像処理装置の構成を詳細に説明する。本発明の実施形態に係る画像処理装置は、情報システム１２０放射線撮影システム１００と接続されている制御部１０５であり、一又は複数のコンピュータで構成される。制御部１０５を構成するコンピュータは、主制御部であるＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）２０１、記憶部であるＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）２０２、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）２０５、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）２０６、グラフィック制御部であるＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）２０８、通信部であるＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）２０３及び２０４、接続部であるＵＳＢ（Ｕｎｉｖｅｒｓａｌ Ｓｅｒｉａｌ Ｂｕｓ）２０７、ＨＤＭＩ（登録商標）（Ｈｉｇｈ Ｄｅｆｉｎｉｔｉｏｎ Ｍｕｌｔｉｍｅｄｉａ Ｉｎｔｅｒｆａｃｅ）２０９を有し、これらは内部バスにより通信可能に接続されている。ＣＰＵ２０１は制御部１０５及びこれに接続する各部を統合的に制御する制御回路である。ＲＡＭ２０２は制御部１０５を及びこれに接続する各部における処理を実行するためのプログラムや、画像処理で用いる各種パラメータを記憶するためのメモリである。ＲＡＭ２０２に展開されたプログラムに含まれる命令がＣＰＵ２０１で逐次実行されることにより、後述する画像処理が実現される。たとえば、通信部である第一のＮＩＣ２０３は放射線撮影を行う施設のアクセスポイントに接続し、第二のＮＩＣ２０４は情報システム１２０内の通信を中継するアクセスポイントに接続する。ＳＳＤ２０６は上述したようなプログラムや、撮影により得られる放射線画像、付帯情報、その他各種パラメータが記憶される。ＵＳＢ２０７は操作部１０７と接続する。ＧＰＵ２０８は画像処理ユニットであり、ＣＰＵ２０１からの制御に応じて画像処理を実行する。画像処理の結果得られる画像はＨＤＭＩ（登録商標）２０９を介してモニタ１０６に出力され、表示される。モニタ１０６及び操作部１０７はタッチパネルモニタに統合されていても良い。

ＳＳＤ２０６に記憶されるプログラムは、たとえば撮影制御モジュール２１１、通信制御モジュール２１２、画像取得モジュール２１３、出力モジュール２１４、表示制御モジュール２１５、推定モジュール２２０である。

撮影制御モジュール２１１は、ＣＰＵ２０１に放射線撮影の実行を統合的に制御させるためのプログラムである。撮影制御モジュール２１１により、たとえば操作入力に応じた撮影条件の指定や、ＦＰＤ１０２の状態を要求する信号の送信などが行われる。

通信制御モジュール２１２は、第一のＮＩＣ２０３及び第二のＮＩＣ２０４による通信を制御する。通信制御モジュール２１２により、たとえば操作部１０７からの入力に応じてＦＰＤ１０２を撮影可能な状態に遷移させるための信号を送信させる。

画像取得モジュール２１３は、ＣＰＵ２０１により実行されることにより、本発明の実施形態に係る画像処理に供する画像を取得する工程の制御を行う。たとえば、ＮＩＣ２０３に、ＦＰＤ１０２により撮影された放射線画像を受信させる。放射線画像の受信において、データ量の小さい、放射線画像の縮小画像を先行的に受信させた後、放射線画像のうち当該縮小画像以外のデータを受信させて放射線画像の受信を完了させることとしてもよい。縮小画像は、たとえばＦＰＤ１０２を構成し、放射線画像の画素値を与える複数の放射線検出素子のうちの偶数列を読み出すなど、選択的に一部の素子から読み出した出力信号だけを用いて得られる。または、いくつかの素子をまとめて読みだしてもよい。読みだした画像を複数の小領域に分割して、小領域の代表知を用いて縮小画像としてもよい。あるいは、ＮＩＣ２０３に、ＰＡＣＳ１１４やその他ネットワーク上の記憶部に記憶された放射線画像を受信させる。あるいは画像処理装置１０５のＳＳＤ２０６やその他の記憶部に記憶されている放射線画像を読み出す。また、本発明の実施形態に係る画像処理に供する前に、周知の画像処理を施すように制御してもよい。

推定モジュール２２０は、ＣＰＵ２０１により実行され、放射線画像が一次放射線成分と散乱線成分との合計で表現されるという条件と、散乱線モデルにより散乱線成分を一次放射線成分から求めることができるという条件と、に基づいて、放射線画像の一次放射線成分あるいは散乱線成分を推定する。推定モジュール２２０は例えば散乱線強度近似モジュール２２１、周波数処理モジュール２２２、及び最適化モジュール２２３を有する。これらのモジュールがＣＰＵ２０１により実行される。これにより、ＣＰＵ２０１はＧＰＵ２０８を制御し、ＧＰＵ２０８により散乱線像の推定処理が実行される。

散乱線強度近似モジュール２２１は、ＣＰＵ２０１により実行され、一次放射線成分Ｐを入力することにより散乱線成分の強度を求める。例えば、第一の関数について、一次放射線成分のｎ次近似解Ｐ_ｎを入力し、第一の出力データを得る。また第二の関数について、一次放射線成分のｎ次近似解Ｐ_ｎを入力し、第二の出力データを得る。

周波数処理モジュール２２２は、ＣＰＵ２０１により実行されることにより、関数の出力についての特定の周波数成分を得る。例えば第一の出力データにガウス関数Ｇ_１をコンボリューションし、第一の出力データについて特定の周波数成分を得る。また第二の出力データにガウス関数Ｇ_２をコンボリューションし、第二の出力データについて特定の周波数成分を得る。これらを所定の比率で重みづけして合成することにより、散乱線成分のｎ次近似解Ｓ_ｎを得る。

最適化モジュール２２３は、ＣＰＵ２０１により実行され、逐次近似法に基づく漸化式により、Ｐ_ｎとＳ_ｎに基づいてＰ_ｎ＋１が求められる。求められた値は、散乱線強度近似モジュール２２１に再び入力される。最適化モジュール２２３は収束判定し、収束したと判定されるまで、同様の順序で処理が繰り返される。収束したと判定された時点の近似解、例えばＰ_ｎ、Ｓ_ｎを、放射線画像の一次放射線成分、散乱線成分の推定結果とする。かかる処理の詳細については図３を参照しながら後述する。

出力モジュール２１４は、ＣＰＵ２０１により実行されることにより、本発明の実施形態に係る画像処理を施して散乱線成分を低減した補正画像の出力を制御する。たとえば、出力モジュール２１４はモニタ１０６に補正画像を出力することにより、モニタ１０６に補正画像を表示させる。またたとえば、出力モジュール２１４はＮＩＣ２０４を介して補正画像をＰＡＣＳ１１４やＰｒｉｎｔｅｒ１１６に出力する。これにより、ＰＡＣＳ１１４に補正画像が保存され、Ｐｒｉｎｔｅｒ１１６によって補正画像がフィルム等に出力される。また、出力モジュール２１４は、制御部１０５内外のその他の記憶部（不図示）に補正画像を出力し、補正画像を記憶させてもよい。さらに、ＤＩＣＯＭ規格に則って種々の情報を付帯させて補正画像を出力することが好ましい。モダリティとは、患者を撮影し、医用画像を生成する画像生成装置である。本発明の実施形態に係る情報システム１２０においては、たとえばＸ線源１０１とＦＰＤ１０２とを有する放射線撮影システム１００がモダリティに該当する。このとき、Ｍｏｄａｌｉｔｙタグ（０００８，００６０）としてＤｉｇｉｔａｌ Ｒａｄｉｏｇｒａｐｈｙを示すＤＸを付帯させる。動画撮影の場合にはＲａｄｉｏ Ｆｌｕｏｒｏｓｃｏｐｙを示すＲＦを付帯させる。さらにＰＡＣＳ１１４に保存する場合にはＳｅｒｖｉｃｅとＯｂｊｅｃｔのＰａｉｒを指定するタグであるＳＯＰ Ｃｌａｓｓ ＵＩＤ（０００８，００１６）タグとして、ＯｂｊｅｃｔのＤｉｇｉｔａｌ Ｘ−ｒａｙ ＩｍａｇｅとＳｅｒｖｉｃｅのＳｔｏｒａｇｅの組合せを示す１．２．８４０．１０００８．５．１．４．１．１．１．１を付帯させる。

表示制御モジュール２１５は、モニタ１０６に表示される内容を制御する。例えば、患者情報、撮影条件の情報、ＦＰＤ１０２の状態を示す情報等をモニタ１０６に表示させる制御を行う。これらは、上述の補正画像とともにモニタ１０６に表示される。

別の実施形態では、出力制御モジュール２１４による補正画像をモニタ１０６に表示させる表示制御を、表示制御モジュール２１５で実行することとしてもよい。この場合、表示制御モジュール２１５で、表示画面に撮影された放射線画像を表示させたり、補正画像を表示させたりする。

本発明の実施形態に係る画像処理を実行する主体は制御部１０５である。以降、本発明の実施形態に係る画像処理を説明する観点から制御部１０５を画像処理装置１０５と呼ぶことがある。

なお、制御部１０５が有する構成要素の一部又は全部は、制御部１０５に固定されるものではなく、情報システム１２０に含まれる画像処理システムとして実現してもよい。たとえば、画像取得モジュール２１３、出力制御モジュール２１４、及び推定モジュール２２０を有する画像処理プログラムを実行する画像処理装置を、撮影制御モジュール２１１を実行する制御部１０５とは別に設けることとしてもよい。またたとえば、ＷＳ１１３がこれらのモジュールの一部または全部を有していてもよい。ＰＡＣＳ１１４がこれらのモジュールの一部または全部を有していてもよい。ＦＰＤ１０２が、たとえば推定モジュール２２０を有するＦＰＧＡ（ｆｉｅｌｄ−ｐｒｏｇｒａｍｍａｂｌｅ ｇａｔｅ ａｒｒａｙ）を備えていてもよい。制御部１０５が有する構成要素が異なる装置に重複して含まれていてもよく、操作者の指示に従って処理を行う装置を選択できるようにしてもよい。さらに、ネットワークを介して接続されたワークステーション、サーバ、記憶装置によって構成されていてもよく、必要に応じてこれらの装置と通信して本発明の実施形態に係る画像処理を行うようにしてもよい。

次に、図３に基づいて制御部１０５が実行する画像処理を詳細に説明する。下記の処理において、特に断りがない場合、処理の主体は、ＣＰＵ２０１またはＧＰＵ２０８であり、推定モジュール２２０をＣＰＵ２０１で実行することにより、下記の散乱線成分を推定する処理が実行される。

放射線撮影により得られた画像は、Ｘ線源１０１からＦＰＤ１０２の各素子に直線的に到達した一次Ｘ線による一次Ｘ線像に、被写体１０３内で散乱した散乱Ｘ線による散乱Ｘ線像が重畳している。

ここで放射線撮影により得られた画像（入力画像）をＭ^´（ｘ，ｙ）、入力画像に含まれる一次Ｘ線成分をＰ^´（ｘ，ｙ）、散乱Ｘ線成分をＳ^´（ｘ，ｙ）とする。（ｘ，ｙ）は画像の画素位置もしくは、その画素位置の画素値を与えるＦＰＤ１０２の素子位置を表す。これらの画像の関係は式１で表すことができる。

散乱Ｘ線は、Ｘ線源から照射されたＸ線が被写体１０３を透過する過程で発生し、散乱を受けなかった成分が一次Ｘ線であるので、Ｓ^´（ｘ，ｙ）はＰ^´（ｘ，ｙ）と相関がある。しかしながら、この相関は後述する式で表されるように、非線形である。したがって、画像処理装置１０５においては反復的な処理によりＳ^´（ｘ，ｙ）とＰ^´（ｘ，ｙ）を最適化することとし、最尤推定法を用いる。他にも非線形最適化には様々な手法があり、最小二乗法、ニュートン法、凸解析法等の手法を用いても良い。

ステップＳ３００で、ＣＰＵ２０１が、画像取得モジュール２１３を実行することにより、被写体１０３に放射線を照射して撮影された放射線画像を取得する。かかる放射線画像は、散乱線推定処理の対象である入力画像となる。かかる入力画像として、データ量の小さい縮小画像を取得し、これを推定処理に用いることとすれば、ＦＰＤ１０２からのデータ送信とそれに続く画像処理をより高速に行うことができる。また散乱線成分は低周波成分が主であるので、縮小画像から推定しても散乱線成分の推定精度に与える影響が小さい。

ステップＳ３０１ではＣＰＵ２０１が、散乱線強度近似モジュール２２１を実行することにより、入力画像を正規化する。たとえば、ＦＰＤ１０２に到達している放射線量である一次Ｘ線量と散乱Ｘ線量の合計値を被写体１０３に対する入射線量で正規化する。入射線量は、被写体１０３が存在しないと仮定したときにＦＰＤ１０２で検出される放射線量としてよい。入射線量は管電圧、管電流、照射時間、撮影部位といった撮影条件から推測してもよいし、被写体１０３を透過しなかったＸ線が到達した領域（素抜け領域）の画素値から求めてもよい。入射線量Ｑで正規化するならば、それぞれ次のように表せる。

以下、Ｍ（ｘ，ｙ）、Ｐ（ｘ，ｙ）、Ｓ（ｘ，ｙ）は正規化されたものとして扱う。

ステップＳ３０２ではＣＰＵ２０１が、散乱線強度近似モジュール２２１を実行することにより、最尤推定法におけるＰ（ｘ，ｙ）の初期値Ｐ^０（ｘ，ｙ）を設定する。Ｐ^０（ｘ，ｙ）は正値であればよいが、一次Ｘ線成分Ｐ（ｘ，ｙ）は入力画像Ｍ（ｘ，ｙ）に近いと考えられるため、本実施形態ではＭ（ｘ，ｙ）を初期値とする。

Ｓ３０３では、散乱線強度近似モジュール２２１がＣＰＵ２０１に、第一の関数に例えば入力画像Ｍを代入して得られる第一の出力データとして、第一の散乱Ｘ線強度を計算させる。ここで第一の関数とは、Ｘ線が散乱して生じた散乱Ｘ線であって、Ｘ線が散乱された後に減衰されないと仮定した散乱Ｘ線による散乱線強度を近似した関数（近似式）である。第一の出力データは、第一の関数の出力である。図４は第一の関数で近似する散乱Ｘ線を例示する図である。Ｘ線源１０１の焦点４０６から被写体１０３内のある位置４０２を通過するパスＬ４０３上で生じた散乱Ｘ線が、散乱後に減衰されずに，放射線検出素子１０２ａのうちのある素子位置４０４（ｘ，ｙ）に到達すると仮定する。ここでは、焦点４０６から画素位置４０４までのパス上にｚ軸を設定し、ｚ軸正方向は一次Ｘ線の進行方向とする。原点はＸ線が被写体１０３に入射する位置とする。減弱係数は被写体１０３内で一定の値μとする。このとき被写体１０３内部のパスＬ４０３上の位置ｚ４０２で発生する散乱Ｘ線は、

に比例するとする。また、正規化された一次Ｘ線Ｐについて

が成り立つ。したがって、パスＬ４０３上で発生し、散乱後に減衰されずに素子位置４０４（ｘ，ｙ）に到達する第一の散乱Ｘ線の全強度は

と表すことができ、一次Ｘ線成分Ｐ（ｘ，ｙ）の関数として与えられる。すなわち、たとえば１−Ｐが第一の関数である。

第一の関数で近似される散乱Ｘ線は、被写体に入射した放射線が位置４０２で散乱され、例えば、散乱後にほとんど減衰されることなく素抜け領域の画素位置４０５に到達する散乱線を、Ｘ線源から画素位置４０４までの直線経路に近似したものと考えてもよい。

画像処理装置１０５においては、第一の関数とは異なる複数の関数に基づいて散乱線像を推定することが好ましい。ステップＳ３０４では、散乱線強度近似モジュール２２１がＣＰＵ２０１に、第二の関数に例えば入力画像Ｍを代入して得られる第二の出力データとして、第二の散乱Ｘ線強度を計算させる。ここで第二の関数とは、Ｘ線が散乱して生じた散乱Ｘ線であって、散乱の後に被写体１０３内でさらに減衰すると仮定した散乱Ｘ線による散乱線強度を近似した関数（近似式）である。第二の出力データは、第二の関数の出力である。図５は第二の関数で近似する散乱Ｘ線を例示する図である。ここでは、主として被写体１０３内を透過してきた散乱Ｘ線を近似している。Ｘ線源１０１の焦点５０６から、放射線検出素子１０２ａのうちのある素子位置５０４（ｘ，ｙ）を通過するパスＬ５０３上で発生した散乱Ｘ線が、パスＬ５０３を通過、減衰してＦＰＤ１０２に入射すると仮定する。ここでは、焦点５０６から素子位置５０４までのパス上にｚ軸を設定し、ｚ軸正方向は一次放射線の進行方向とする。原点はＸ線が被写体１０３に入射する位置とする。減弱係数はμである。このとき被写体１０３内部のパスＬ５０３上の位置ｚ５０２で発生する散乱Ｘ線は式２に比例するとする。さらに被写体１０３中を距離Ｌ−ｚ透過して減衰されると仮定し、パスＬで発生する第二の散乱Ｘ線の全強度は

と表すことができ、一次Ｘ線成分Ｐ（ｘ，ｙ）の関数として与えられる。すなわち、たとえば−ＰｌｎＰが第二の関数である。

散乱Ｘ線は被写体１０３内で広い角度に散乱されることを考慮すると、それぞれの関数の出力データについて特定の周波数成分を合成して散乱線成分を推定することが好ましい。ステップＳ３０５では、周波数処理モジュール２２２がＣＰＵ２０１に周波数処理を行わせることにより、ステップＳ３０４までの処理で得られる出力データを複数の周波数成分に分解する。まず、散乱Ｘ線成分は低周波の成分が主であるので、ガウス関数Ｇ_１をコンボリューションする。第一の出力データについての低周波成分は

と表せる。さらに、散乱Ｘ線の不連続性を考慮すると、低周波以外の成分を加味することが好ましい。第一の散乱Ｘ線モデルによる散乱線像の低周波以外の成分は、

と表せる。ここで、Ｆはフーリエ変換、Ｆ^−１はフーリエ逆変換を表す。

ＣＰＵ２０１は、周波数処理モジュール２２２を実行することにより、上述の処理で得られた複数の周波数成分を、パラメータＡ_１、Ｃ_１によって重み付けして合成する。これにより、第一の出力データについての周波数成分Ｓ_１を得る。
ここまでの処理結果によれば、Ｓ_１は

と表せる。ここでＡ_１は第一の散乱Ｘ線の低周波成分の重み、Ｂ_１は第一の散乱Ｘ線の広がり、Ｃ_１は第一の散乱Ｘ線の低周波以外の成分の重みを表すパラメータである。

ステップＳ３０６では、第二の出力データについてもＳ３０５と同様にして特定の周波数成分を得る。低周波成分はガウス関数Ｇ_２をコンボリューションして、

と表せる。低周波以外の成分は、

と表せる。

ＣＰＵ２０１は、周波数処理モジュール２２２を実行することにより、上述の処理で得られた複数の周波数成分を、パラメータＡ_２、Ｃ_２によって重み付けして合成する。これにより、第二の出力データについての周波数成分Ｓ_２を得る。
ここまでの処理結果によれば、Ｓ_２は

と表せる。ここでＡ_２は第二の散乱Ｘ線の低周波成分の重み、Ｂ_２は第二の散乱Ｘ線の広がり、Ｃ_２は第二の散乱Ｘ線の低周波以外の成分の重みを表すパラメータである。

ここでＧ_１とＧ_２は同じであってもよいが、第一の出力データと第二の出力データをそれぞれ異なる複数の周波数帯域に分解し、重み付けして合成することが好ましい。

第一の出力データの周波数成分として式９、式１１で表される２種類の成分、第二の出力データの周波数成分として、式１３と式１５で表される２種類の成分を例示した。第一の出力データの周波数成分の２種類のうちから少なくとも１種類と、第二の出力データの周波数成分の２種類のうちから少なくとも１種類とを重み付けして合成し、散乱線成分としてもよい。

上述の処理により、第一の出力データの周波数成分Ｓ_１と第二の出力データの周波数成分Ｓ_２に基づいて、散乱線成分を得ることができる。すなわち画像処理装置１０５における散乱線モデルとは、第一の関数の出力についての特定の周波数成分と、前記第一の関数とは異なる第二の関数の出力についての特定の周波数成分と、に基づき散乱線成分を得る散乱線モデルである。

Ａ_１、Ｂ_１、Ｃ_１、Ａ_２、Ｂ_２、Ｃ_２は複数の組合せをＲＡＭ２０２に記憶し、自動で又は操作者が適宜設定できるようにしてもよい。当該放射線画像の撮影条件であるＸ線管の管電流、管電圧、Ｘ線の照射時間、撮影部位等に応じて合成の重みを設定する設定モジュールを有してもよい。例えば、設定モジュールは、操作部１０７に対する操作入力に応じて上述の撮影条件の情報を取得し、当該撮影条件に対応するパラメータセット（Ａ_１、Ｂ_１、Ｃ_１、Ａ_２、Ｂ_２、Ｃ_２）を得る。ＳＳＤ２０８に撮影条件とパラメータセットとの対応関係を示すテーブル情報を記憶させておき、当該テーブル情報を参照することにより、ＣＰＵ２０１がパラメータセットを設定する。当該テーブル情報は、例えば撮影条件ごとに予め実験的に定められる。

上述の通り、散乱Ｘ線像は低周波の成分が主であるため、低周波成分の重みを表すパラメータＡ_１、Ａ_２が、低周波以外の成分の重みを表すパラメータＣ_１、Ｃ_２よりも大きな値であることが好ましい。

なお、周波数処理モジュール２２２はフーリエ変換を用いて、低周波成分と、低周波以外の成分とを作成しても良いし、実空間で行ってもよい。また、周波数処理に用いるフィルタはガウス関数に限らない。たとえば、加重平均フィルタを用いれば対象画素の周辺を中心とした広がりを表現できるし、移動平均フィルタを用いれば対象画素から離れた画素までの広がりを表現できる。

また、Ｓ３０３、Ｓ３０４、Ｓ３０５、Ｓ３０６の順序を例示したが、Ｓ３０３、Ｓ３０５、Ｓ３０４、Ｓ３０６の順序で処理を行ってもよく、Ｓ３０３とＳ３０５の処理と、Ｓ３０４とＳ３０６の処理を並行して行うようにしてもよい。

Ｓ３０７では、最適化モジュール２２３がＣＰＵ２０１に、Ｓ３０６までに計算された第一の出力の周波数成分と、第二の出力の周波数成分を合成して散乱線成分を得て、一次Ｘ線成分の最尤推定を行わせる。最尤推定とは、たとえば式１７に表されるように、入力画像Ｍ（ｘ，ｙ）を、反復の前段階で得られた一次Ｘ線成分Ｐ^ｎ（ｘ，ｙ）と第一の出力データの周波数成分Ｓ_１ ^ｎ（ｘ，ｙ）と第二の出力データの周波数成分Ｓ_２ ^ｎ（ｘ，ｙ）の和で除算し、一次Ｘ線成分Ｐ^ｎ（ｘ，ｙ）に乗算する。ここで上添え字ｎはｎ回目の反復で得られたものであることを示す。

式１７に基づく最尤推定では、反復を繰り返すごとに

が１に収束していくことが知られている。逐次的に、式１の関係を満たすようにして一次Ｘ線成分Ｐ（ｘ，ｙ）、第一の出力データの周波数成分Ｓ_１（ｘ，ｙ）、第二の出力データの周波数成分Ｓ_２（ｘ，ｙ）を求めていく。

Ｓ３０８では、最適化モジュール２２３がＣＰＵ２０１に反復継続の判定をさせる。式１８が十分に１に近づいたか、あるいは式１の関係が十分に満たされたかを二乗残差などによって判定する。もしくは固定回数の反復が行われたかで判定してもよい。

Ｓ３０８で反復を継続する判定がなされた場合、計算されたＰ^ｎ＋１（ｘ，ｙ）を初期値としてＳ３０２に戻り処理を継続する。反復を継続しない判定がなされた場合は処理を終了し、Ｐ^ｎ＋１（ｘ，ｙ）を散乱Ｘ線成分の低減された補正画像とする。あるいは、計算されたＳ_１ ^ｎ＋１（ｘ，ｙ）とＳ_２ ^ｎ＋１（ｘ，ｙ）の和にパラメータを乗じて入力画像から減じてもよい。すなわち、散乱線像の低減効果を調節可能としても良い。

このように、画像処理装置１０５においては、第一の出力データと第二の出力データのそれぞれについて、特定の周波数成分を合成した散乱線成分Ｓ_１（ｘ，ｙ）＋Ｓ_２（ｘ，ｙ）を散乱線成分として出力する。さらにこの散乱線成分と一次Ｘ線成分との和が入力画像と一致するとの仮定の下、逐次近似的解法により放射線画像に含まれる散乱線成分を導出する。

ステップＳ３０９でＣＰＵ２０１は出力制御モジュール２１４を実行することにより、上述の処理で得られる補正画像を出力する。かかる工程では当該補正画像がＰＡＣＳ１１４に保存され、モニタ１０６に表示される。推定された散乱Ｘ線像であるＳ_１ ^ｎ＋１（ｘ，ｙ）＋Ｓ_２ ^ｎ＋１（ｘ，ｙ）を当該補正画像とは別の画像データあるいは画像ファイルとして、ＰＡＣＳ１１４に保存しても良い。

別の観点では第一の関数と第二の関数は、放射線が、放射線画像の画素値を与える放射線検出素子１０２ａに到達するまでの経路上で散乱して生じた散乱線の全強度を、当該経路上で減衰される距離が異なる複数の散乱線モデルにより当該経路毎に近似した関数である。さらに、上述の処理では、放射線検出素子からの出力に対応するそれぞれの画素値に対して一様に当該近似の処理を行って散乱線像を推定している。これにより、放射線画像に含まれる画素値毎にモデルを切り替える必要がなくなり、結果的に散乱線像の推定処理に要する時間やリソースを低減することができる。

また、式１９に示すように、散乱後に減衰される距離はパラメータによって調整できる。たとえば、散乱後に被写体１０３内で距離ｘでは減衰され、他では減衰されないと仮定するならば、

と表せる。ｘは０＜ｘ＜（Ｌ−ｚ）を満たすものとする。
散乱後に被写体１０３内で減衰されるが、距離ｙにおいては減衰されないと仮定するならば、

と表せる。ｙは０＜ｙ＜（Ｌ−ｚ）を満たすものとする。

ＸおよびＹは、パラメータＡ_１、Ｂ_１、Ｃ_１、Ａ_２、Ｂ_２、Ｃ_２に含まれ得る。したがって、本発明はこの二つのモデルのみに限られるものではなく、適宜適切なパラメータを選択することにより様々な減衰距離で近似されたモデルを含む。

パラメータＡ_１、Ｂ_１、Ｃ_１、Ａ_２、Ｂ_２、Ｃ_２は同値変形の範囲内において、いかなる項に乗じてもよい。すなわち、第一の関数及び第二の関数は、上述したような散乱線強度の近似の考え方に限定されず、他の項に各パラメータを乗じたものも本発明に含む。

また、画像処理装置１０５における第一の関数は、１−Ｐに限定されない。１−Ｐは、（１−１）正規化された放射線画像中で一次放射線成分が１すなわち最大値である場合に、散乱線成分が０すなわち最小値となる性質と、（１−２）一次放射線成分が０すなわち最小値である場合に、散乱線成分は１すなわち正の値である性質と、（１−３）一次放射線成分が最小値に近付くに従って正の値に収束する性質と、を有する関数である。また１−Ｐは、（１−４）一次放射線成分が０から１の範囲内で単調減少する性質を有する関数である。また１−Ｐは、（１−５）一次放射線成分が０である場合に、導関数が負の値をとる性質を有する関数である。すなわち、このような（１−１）乃至（１−５）の性質の一部又は全部を満たす関数であって、一次放射線成分と散乱線成分の値との対応関係を示す情報を第一の関数としてもよい。たとえば、式２１が挙げられる。

さらに、画像処理装置１０５における第二の関数は、−ＰｌｎＰに限定されない。−ＰｌｎＰは、（２−１）正規化された放射線画像中で一次放射線成分が１すなわち最大値である場合に、散乱線成分は０すなわち最小値となる性質と、（２−２）一次放射線成分が０すなわち最小値である場合に、散乱線成分は０すなわち最小値となる性質と、（２−３）一次放射線成分が０でも１でもない場合に、散乱線成分が正の値となる性質と、を有する関数である。また−ＰｌｎＰは、（２−４）一次放射線成分が０から１の範囲内で極大値をとる性質を有する関数である。また−ＰｌｎＰは、（２−５）一次放射線成分が０である場合に、導関数が正の値をとる性質を有する関数である。すなわち、このような性質の一部又は全部を満たす関数であって、一次放射線成分と散乱線成分の値との対応関係を示す情報を第二の関数としてもよい。たとえば、式２２が挙げられる。

前記第一の関数または第二の関数は、上述に示すような数式で表されていなくても良く、例えば一次Ｘ線成分Ｐの値に対する出力値を表すルックアップテーブルで表されていてもよい。

第一の関数に基づくことで、たとえば被写体に厚みがある場合など一次Ｘ線が大きく減少して散乱Ｘ線が多く発生している場合に、精度良く散乱線像を推定できる。例えば、被写体に厚みがある部分では一次Ｘ線はほとんど通過せず、発生した散乱Ｘ線もほとんど通過しない。具体的には、Ｐがゼロに近い値を取るときの散乱Ｘ線強度は、第二の関数によれば散乱Ｘ線もゼロに近い値に推定される。しかし第一の関数によれば１に近い値に推定されるため、本来発生している散乱Ｘ線を少なく見積もることがなく、精度良く推定できる。さらに、散乱Ｘ線が大きく広がることをふまえ、複数の周波数成分に分解し、それぞれの成分での広がり方を考慮したパラメータによって合成することで、被写体の存在しない領域や、被写体の厚みが薄い部分に広がった散乱Ｘ線を精度よく推定することができる。この結果、被写体が存在しない部分や、図９（ａ）に示す肺野領域９０４、スキンライン９０３等の散乱線像をより正確に推定することができる。

さらにまた別の観点では、入射したＸ線の減衰によって一次の散乱Ｘ線が生じ、第二の関数−ＰｌｎＰは、一次の散乱Ｘ線が被写体１０３内で減衰される過程において、多次の散乱Ｘ線は生じないとする仮定に基づく関数と考えてもよい。この観点によれば、第一の関数１−Ｐは、入射したＸ線の減衰によって一次の散乱Ｘ線が生じ、一次の散乱Ｘ線が被写体１０３中を通過して減衰される過程において、被写体１０３に吸収されることなく、多次の散乱Ｘ線を生じたとする仮定に基づく関数と考えることができる。実際には被写体１０３内で生じた散乱Ｘ線は多次の散乱Ｘ線を生じ得る。すなわち、第二の関数−ＰｌｎＰの出力データに対して第一の関数１−Ｐの出力データを重み付けして合成することは、被写体１０３に入射したＸ線から生じ得る多次の散乱Ｘ線と、一次の散乱Ｘ線の割合を重み付けして合成することである。よって、被写体１０３内で生じる散乱Ｘ線の減衰のされ方を適切に表現することができ、推定される散乱線成分の定量性を向上させる。すなわち、散乱線成分をより正確に推定できる。例示した２種類の関数の出力データに対して、特定の周波数成分を得てそれぞれ合成することで、減衰のされ方が異なる散乱Ｘ線の広がり方の違いを考慮することができ、さらに散乱線成分の推定精度が向上する。

第一の関数と、第二の関数に加えて、上述した様々な関数の中から第三の関数も用いることにしても良い。この場合、第一、第二及び第三それぞれの関数に入力画像若しくは式１７で計算された近似解であるＰ^ｎ（ｘ，ｙ）を代入して得られるそれぞれの出力データを用いて散乱線成分を推定する。第三の関数は第一の関数若しくは第二の関数と同じでもよいが、それぞれの出力データについて異なる周波数成分を合成するようにすることが好ましい。たとえば、第三の出力データにガウス関数Ｇ_３を用いて特定の周波数成分を得て、上述したような第一の出力データの周波数成分と第二の出力データの周波数成分と重み付けして合成してもよい。

放射線画像に適用する散乱線モデルとしては、散乱線が被写体１０３内で減衰される距離が様々であること、別の観点からすると多次の散乱線の生じ方が様々であることを表現できればよい。すなわち、関数としては単一であっても、複数の周波数成分に分解して得られる少なくとも２つの周波数成分を重み付けして合成すればよい。パラメータが異なることで、表現し得る散乱線の挙動も異なり、少なくとも低周波数側の成分と高周波数側の成分の２つの広がり方を考慮することで、散乱線の分布の変化を推定することができる。たとえば、式９と式１１で表現されるそれぞれの成分を重み付けして合成してもよい。

Ｘ線が被写体１０３内で減衰する現象は、Ｘ線が散乱したり吸収されたりすることにより、ＦＰＤ１０２に到達するＸ線強度が被写体１０３に入射したＸ線と比べて減少する現象である。被写体１０３内で生じた散乱Ｘ線であって、散乱後の透過距離に応じて減衰される散乱Ｘ線とは、一回散乱して生じた散乱Ｘ線が散乱後の透過の過程でさらに多次の散乱Ｘ線を生じたり吸収されたりするものの、多次の散乱線についてはＦＰＤ１０２に到達せず一次の散乱Ｘ線のみがＦＰＤ１０２に検出される場合を仮定した散乱Ｘ線である。このような散乱Ｘ線は被写体の厚みが大きい部分で優位となると考えられる。被写体１０３内で生じた散乱Ｘ線であって、散乱後も減衰されない散乱Ｘ線とは、被写体１０３内で放射線は吸収されず、生じた多次の散乱Ｘ線は全てＦＰＤ１０２に検出される場合を仮定した散乱Ｘ線である。被写体の厚みが小さい部分や、散乱して広がった結果素抜け領域に到達する散乱Ｘ線においては、一次の散乱Ｘ線のみではなく多次の散乱Ｘ線も含まれると考えられる。したがって、これらの２つの場合をそれぞれ仮定した関数に基づくことにより、放射線が被写体１０３内で減衰あるいは散乱する態様をより正確に表現することができる。

様々な観点から散乱線モデルについて説明したが、これまでに示した関数の変形例は、当然ながら、散乱線モデルをいかなる観点から捉えた場合にも用いることができ、本発明に含む。

以上のように、第一の関数と、第一の関数とは異なる第二の関数に基づくことにより、精度良く散乱線像を推定することができる。また、低周波数側の成分だけでなく、複数の周波数成分に分解して得られる少なくとも２つの周波数成分を重み付けして合成することにより、散乱線の分布が変化する境界を含む散乱線成分を推定できる。

次に、本発明の別の実施形態について説明する。上述した実施形態と同様に、第一の関数の出力についての特定の周波数成分と、前記第一の関数とは異なる第二の関数の出力についての特定の周波数成分と、に基づき散乱線成分を得る散乱線モデルを放射線画像に適用することにより、放射線画像に含まれる散乱線成分を推定する。ここでは、第一の出力データと第二の出力データを３つ以上の周波数帯域に分解して周波数成分を得る処理を行う。

図６は、本発明の別の実施形態に係る処理を示すフローチャートである。ステップＳ６０１、Ｓ６０２、Ｓ６０３、Ｓ６０６、Ｓ６０９、Ｓ６１０については、図３で示すＳ３０１、Ｓ３０２、Ｓ３０３、Ｓ３０５、Ｓ３０７、Ｓ３０８とそれぞれ同様であるので、詳しい説明を省略する。

ステップＳ６０３では、Ｓ３０３と同様にして第一の出力データを得る。

ステップＳ６０４ではＳ６０３で計算された第一の出力データに対し、周波数処理モジュール２２２がＣＰＵ２０１に多重周波数分解処理を行わせる。本実施形態においては汎用性があり高速なラプラシアンピラミッド分解を用いる。その他周波数帯域に分解する手法としては高速フーリエ変換、ウェーブレット分解などがあり、これらを用いてもよい。

図７（ａ）はラプラシアンピラミッド分解の概念図である。Ｓ６０３で計算された第一の出力データをＬ_１とすると

と表される。Ｌ_１に図７（ｂ）に示すような二次元ローパスフィルタＡＦでエイリアシング処理を行ってから、縦横１／２のサイズにダウンサンプリングすることで低周波画像Ｌ_２を得ることができる。

ここで↓は１／２ダウンサンプリングを示す。同様にして逐次的に低周波画像Ｌ_ｎを得ることができる。

また、高周波画像Ｈ_ｎはＬ_ｎ＋１を二倍にアップサンプリングし、図７（ｂ）に示すような二次元ローパスフィルタＡＦを用いてエイリアシング処理を行って、Ｌ_ｎから差分することによって得ることができる。

ここで↑は二倍アップサンプリングを示す。式１５、式１６の処理を逐次的に行うことで、ｎ種類の周波数帯域に分割され、Ｌ_ｎ画像およびＨ_ｎ画像を作成できる。このようにして複数の周波数帯域の画像が作成される。

図８はラプラシアンピラミッドによる周波数分解の周波数特性を示す図である。Ｓ６０４のラプラシアンピラミッド分解によって生成される各周波数帯域の画像Ｌ_ｎ、Ｈ_ｎの周波数特性の例である。たとえば縦横１０００画素の画像をＬ_１０まで分解した場合、Ｌ_１０では縦横数画素のレベルまで分解される。また、各レベルの高周波成分Ｈ_ｎは入力画像Ｌ_１を帯域制限したものとなっている。

ステップＳ６０５では、ステップＳ６０３までで得られた第一の出力データを、ステップＳ６０４で複数の周波数帯域に分解して得られた画像である各周波数分解画像Ｌ_ｎ、Ｈ_ｎを、周波数処理モジュール２２２がＣＰＵ２０１に実行されることにより再構成される。この処理は図７（ｃ）に示すように各レベルのＬ_ｎ＋１画像を縦横２倍にアップサンプリングし、図７（ｂ）に示すようなローパスフィルタＡＦでエイリアシング処理をしてＨ_ｎに加算することでＬ_ｎを再構成できる。たとえば、Ｌ_１０まで分解したとすると、Ｌ_９ ^´は

と表される。したがって、Ｌ_ｎ ^´は

と表せる。この計算を逐次的に行い、本実施形態では再構成されたＬ_１ ^´を第一の出力データの周波数成分Ｓ_１とする。ここで、α_ｎは結合定数である。α_ｎを調節することで、ｎ個の周波数成分を重み付けして合成し、Ｓ_１を表すことができる。

ステップＳ６０６では、Ｓ３０５と同様にして第二の出力データを得る。

ステップＳ６０７では、Ｓ６０６で計算された第二の出力データに対して、Ｓ６０４と同様にして、複数の周波数分解画像Ｌ_ｎ、Ｈ_ｎに分解する。

ステップＳ６０８では、Ｓ６０５と同様にして、Ｓ６０７で計算された各周波数分解画像Ｌ_ｎ、Ｈ_ｎを再構成する。すなわち、

となり、Ｌ_ｎ ^´は

と表される。この計算を逐次的に行い、本実施形態では再構成されたＬ_１ ^´を第一の散乱Ｘ線像Ｓ_２とする。ここで、β_ｎは結合定数である。β_ｎを調節することで、ｎ個の周波数成分を重み付けして合成し、Ｓ_２を表すことができる。

本実施形態でも、Ｘ線源１０１から発生した放射線が、放射線画像の画素値を与える放射線検出素子１０２ａに直線的に到達する経路上で生成した散乱線の全強度の近似手法を、第一の散乱Ｘ線モデルと、第二の散乱Ｘ線モデルの二つを用いて説明した。しかし式１９、式２０に示すように、散乱後の減衰距離の仮定はパラメータによって調整することができる。したがって、本発明はこの二つのモデルのみに限られるものではなく、適宜適切なパラメータを選択することにより様々な減衰距離で近似されたモデルを含み得る。また、ローパスフィルタとして、たとえば加重平均フィルタを用いれば対象画素の周辺を中心とした広がりを表現できるし、移動平均フィルタを用いれば対象画素から離れた画素までの広がりを表現できる。

本実施形態では、α_ｎやβ_ｎを調節することにより、複数の周波数成分で散乱線成分を表現することができる。また、第一の出力データと第二の出力データの少なくともいずれかについて、３つ以上の周波数帯域の成分に分解することとしてもよい。３つ以上の周波数帯域の成分に分解することで、体内で散乱されて大きく広がる低周波成分と、被写体１０３でのＸ線の吸収が大きく変化するスキンラインや肺野などの部分における構造を、より精度良く推定できる。

本実施形態でも、前述した様々な観点から捉えた散乱線モデルおよび関数の変形例を用いることができ、本発明に含む。

本実施形態では放射線撮影で得られた二次元画像に対して、複数の関数を用いたフィッティング処理を施すことで散乱線成分を推定する。本実施形態には、断層撮影の投影像に対して上述の散乱線像の推定処理を適用する例が含まれる。あるいは、被写体１０３に対して一方向から放射線を照射して二次元検出器を用い撮影するいわゆる一般撮影で得られる放射線画像に対して上述の散乱線像の推定処理を適用する例も本発明に含まれる。放射線の吸収が変化する部分をセグメンテーション等の解析処理によって抽出する必要がないため、よりロバストな処理により散乱線成分を推定できる。

また、撮影部位に応じて散乱線低減処理のオン／オフを切り替えられるようにすることが好ましい。

図９（ｂ）は胸部ファントムを用いて取得した放射線画像に含まれる散乱線成分を、本発明に係る画像処理を用いて取得した例を示す図である。実験的に取得した散乱線成分を例示する図である（ａ）と同様に、素抜け領域９０５と肺野領域９０８における強度分布や、スキンライン９０７が表現されている。

上述の実施形態における画像処理装置は単体の装置であったが、複数の情報処理装置を含む装置を互いに通信可能に組合せた画像処理システムで上述の処理を実行する形態も、本発明に含まれる。あるいは複数のモダリティで共通のサーバ装置あるいはサーバ群で、上述の処理を実行することとしてもよい。この場合、当該共通のサーバ装置は実施形態に係る画像処理装置に対応し、当該サーバ群は実施形態に係る画像処理システムに対応する。情報システム１２０あるいは画像処理システムを構成する複数の装置は所定の通信レートで通信可能であればよく、また同一の施設内あるいは同一の国に存在することを要しない。

また、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラムを、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータが該供給されたプログラムのコードを読みだして実行するという形態を含む。

したがって、本実施形態に係る処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明の一つである。また、コンピュータが読みだしたプログラムに含まれる指示に基づき、コンピュータで稼働しているＯＳなどが、実際の処理の一部または全部を行い、その処理によっても前述した実施形態の機能が実現される。

上述の実施形態を適宜組み合わせた形態も、本発明に含まれる。

Claims (28)

1. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、
   前記被写体内で前記放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、前記被写体内で前記放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づいて、前記画像取得手段により取得された放射線画像に含まれる散乱線成分を推定する推定手段と、
   前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力手段と、を有することを特徴とする画像処理装置。
2. 前記推定手段は、前記第一の関数の出力についての特定の周波数成分と、前記第二の関数の出力についての特定の周波数成分とを重み付けして合成することにより、前記散乱線成分を推定するように構成されていることを特徴とする請求項１に記載の画像処理装置。
3. 前記第一の関数の出力についての特定の周波数成分と、前記第二の関数の出力についての特定の周波数成分は、それぞれ周波数帯域が異なることを特徴とする請求項２に記載の画像処理装置。
4. 前記推定手段は、前記第一の関数の出力と前記第二の関数の出力について、それぞれ異なる複数の周波数成分に分解し、前記異なる複数の周波数成分を重み付けして合成することにより、前記放射線画像に含まれる散乱線成分を推定するように構成されていることを特徴とする請求項３に記載の画像処理装置。
5. 前記異なる複数の周波数成分は、前記第一の関数の出力と前記第二の関数の出力の少なくともいずれかについて、３つ以上の周波数帯域に分解された成分であることを特徴とする請求項４に記載の画像処理装置。
6. 前記第一の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す情報であって、前記一次放射線成分が、前記放射線画像中で最大値である場合に前記散乱線成分が最大値となり、前記放射線画像中で最小値である場合に前記散乱線成分が正の値となり、前記放射線画像中で最小値に近付くに従って前記散乱線成分が正の値に収束する関数であることを特徴とする請求項１乃至請求項５のいずれか一項に記載の画像処理装置。
7. 前記第一の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す情報であって、前記一次放射線成分が、前記放射線画像中で最大値である場合に前記散乱線成分が最小値となり、前記放射線画像中での前記一次放射線成分の範囲において単調減少する関数であることを特徴とする請求項１乃至請求項５のいずれか一項に記載の画像処理装置。
8. 前記第一の関数は、前記放射線が散乱して生じた散乱線は減衰されないと仮定した関数であることを特徴とする請求項１乃至請求項５のいずれか一項に記載の画像処理装置。
9. 前記放射線画像の画素値を前記被写体に対する入射線量に対応する値で正規化する正規化手段を有し、
   前記第一の関数は、Ｐを前記正規化された放射線画像における一次放射線成分としたとき、近似式１−Ｐで表されることを特徴とする請求項１乃至請求項４のいずれか一項に記載の画像処理装置。
10. 前記第二の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す情報であって、前記一次放射線成分が、前記放射線画像中で最大値である場合に前記散乱線成分が最小値となり、最大値でない場合と最小値でない場合に前記散乱線成分が正の値となる関数であることを特徴とする請求項１乃至請求項９のいずれか一項に記載の画像処理装置。
11. 前記第二の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す情報であって、前記一次放射線成分が、前記放射線画像中で最大値である場合と最小値である場合に前記散乱線成分が最小値となり、最大値でない場合と最小値でない場合に前記散乱線成分が正の値となり、前記放射線画像中での前記一次放射線成分の範囲において極大値を有する関数であることを特徴とする請求項１乃至請求項９のいずれか一項に記載の画像処理装置。
12. 前記第二の関数は、前記放射線が散乱して生じた散乱線は前記被写体内でさらに減衰されると仮定した関数であることを特徴とする請求項１乃至請求項９のいずれか一項に記載の画像処理装置。
13. 前記放射線画像の画素値を前記被写体に対する入射線量に対応する値で正規化する正規化手段を有し、
    前記第二の関数は、Ｐを前記正規化された放射線画像における一次放射線成分としたとき、近似式−ＰｌｎＰで表されることを特徴とする請求項１乃至請求項９のいずれか一項に記載の画像処理装置。
14. 前記異なる複数の周波数成分は、前記第一の関数の出力と前記第二の関数の出力について、それぞれ異なる２つの周波数成分に分解したものであって、
    前記推定手段は前記異なる２つの周波数成分を重み付けして合成することにより、前記放射線画像に含まれる散乱線成分を推定するように構成されていることを特徴とする請求項１乃至請求項１３のいずれか一項に記載の画像処理装置。
15. 前記推定手段は、前記異なる２つの周波数成分のうち低周波数側の周波数成分について重みを大きくして合成することにより、前記放射線画像に含まれる散乱線成分を推定するように構成されていることを特徴とする請求項１４に記載の画像処理装置。
16. 前記第一の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す関数であって、前記一次放射線成分が、前記放射線画像中で最小値である場合に導関数が負の値となる関数であり、
    前記第二の関数は、前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子まで直進した放射線の成分である一次放射線成分の値と、前記放射線画像の散乱線成分の値との対応関係を示す関数であって、前記一次放射線成分が、前記放射線画像中で最小値である場合に導関数が正の値となる関数であることを特徴とする請求項１乃至請求項１５のいずれか一項に記載の画像処理装置。
17. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、
    第一の関数の出力についての特定の周波数成分と、前記第一の関数とは異なる第二の関数の出力についての特定の周波数成分と、に基づき散乱線成分を得る散乱線モデルを前記放射線画像に適用することにより、前記放射線画像に含まれる散乱線成分を推定する推定手段と、
    前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力手段と、
    を有することを特徴とする画像処理装置。
18. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、
    前記放射線画像に特定の関数の出力を、複数の周波数成分に分解して得られる少なくとも２つの周波数成分を重み付けして合成することにより、前記放射線画像に含まれる散乱線成分を推定する推定手段と、
    前記放射線画像から前記散乱線成分が低減された補正画像を出力する出力手段と、
    を有することを特徴とする画像処理装置。
19. 管電流、照射時間、撮影部位のうち少なくとも一つを含む撮影条件を用いて、前記合成のパラメータを設定する設定手段を有することを特徴とする請求項１乃至請求項１８のいずれか一項に記載の画像処理装置。
20. 放射線画像の一次放射線成分を入力として、該一次放射線成分に第一の関数を適用して得られる第一の関数の出力と前記第二の関数を適用して得られる前記第二の関数の出力のそれぞれについての前記特定の周波数成分を合成することにより該放射線画像の散乱線成分を出力として得る処理手段をさらに有し、
    前記推定手段は、前記放射線画像の一次放射線成分と前記処理手段により得られる該放射線画像の散乱線成分との和が前記放射線画像に一致するとの仮定の下、前記画像取得手段により取得された放射線画像を入力として、該画像取得手段により取得された放射線画像の散乱線成分を逐次近似的解法により導出するように構成されていることを特徴とする請求項１乃至請求項１９のいずれか一項に記載の画像処理装置。
21. 前記逐次近似的解法として最尤推定法を用いることを特徴とする請求項２０に記載の画像処理装置。
22. 前記画像取得手段は、被写体に対して一の方向から放射線を照射して撮影された放射線画像を取得することを特徴とする請求項１乃至請求項２１のいずれか一項に記載の画像処理装置。
23. 前記画像取得手段は、前記放射線画像を得るための放射線検出器から、該放射線画像よりもデータ量の少ない、該放射線画像の縮小画像を先行的に受信した後、該放射線画像の受信を完了させ、
    前記推定手段は、前記縮小画像に基づいて前記放射線画像の散乱線成分を推定し、
    前記出力手段は、前記放射線画像から前記散乱線成分が低減された補正画像を出力するように構成されていることを特徴とする請求項１乃至請求項２２のいずれか一項に記載の画像処理装置。
24. 前記出力手段は、前記補正画像にＭｏｄａｌｉｔｙタグとしてＤＸ又はＲＦを付帯させたＤＩＣＯＭ画像ファイルを出力することを特徴とする請求項１乃至請求項２３のいずれか一項に記載の画像処理装置。
25. 放射線発生装置から被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、
    前記放射線画像に含まれる散乱線成分を推定する推定手段と、
    前記放射線画像から前記散乱線成分が低減された補正画像を出力する出力手段と、を有し、
    前記推定手段は、前記放射線が前記放射線発生装置から前記放射線画像の画素値を与える放射線検出素子に到達するまでの経路上で散乱して生じた散乱線の全強度を、前記経路上で減衰される距離が異なる複数の散乱線モデルにより前記経路毎に近似し、さらに、前記放射線検出素子からの出力に対応するそれぞれの画素値に対して一様に前記近似の処理を行って散乱線像を推定するように構成されていることを特徴とする画像処理装置。
26. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得手段と、
    前記被写体内で前記放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、前記被写体内で前記放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づいて、前記画像取得手段により取得された放射線画像に含まれる散乱線成分を推定する推定手段と、
    前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力手段と、を有することを特徴とする画像処理システム。
27. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得ステップと、
    前記被写体内で前記放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、前記被写体内で前記放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づいて、前記画像取得手段により取得された放射線画像に含まれる散乱線成分を推定する推定ステップと、
    前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力ステップと、を有することを特徴とする画像処理方法。
28. 被写体に放射線を照射することにより得られる放射線画像を取得する画像取得処理と、
    前記被写体内で前記放射線が複数回散乱して生じる多次散乱線を含む散乱線の強度に対応する第一の関数と、前記被写体内で前記放射線が一回散乱して生じる一次散乱線の強度に対応する第二の関数とに基づいて、前記画像取得手段により取得された放射線画像に含まれる散乱線成分を推定する推定処理と、
    前記放射線画像から前記推定手段により推定された散乱線成分が低減された補正画像を出力する出力処理と、をコンピュータに実行させるための画像処理プログラム。