JP6071144B2 - 放射線画像解析装置および方法並びにプログラム - Google Patents

Description

本発明は、被写体を撮影して得られた放射線画像を解析する画像解析装置および方法並びにプログラムに関し、特に、被写体を撮影した放射線画像を解析して放射線画像の各位置における被写体の厚さを推定する画像解析装置および方法並びにプログラムに関する。

従来、被写体を透過した放射線により被写体の放射線画像を撮影する際、被写体の厚さが大きいほど被写体内部における放射線の散乱の発生、放射線透過率の低下などの影響が大きくなり、取得される放射線画像の画質が変動することが知られている。このため、撮影条件と放射線画像の信号値、放射線画像の信号値のヒストグラム幅、被写体画像における被写体の所定方向の長さなど種々の情報によって大まかに被写体の厚さを推定し、推定した被写体の厚さに応じて、撮影された放射線画像に対する散乱線除去処理等の画像処理条件や、放射線画像の撮影に適用される撮影条件を変更する技術が提案されている。

例えば特許文献１には、予め既知の厚みを有する模擬被写体を既知の撮影条件で放射線撮影して得られた画像の画素値を測定することにより、体厚と画素値の関係を対応付けた対応付けテーブルを用意し、対応付けテーブルに基づいて、被写体画像の画素値に応じて概略的な体厚分布を推定し、被写体画像の体厚分布に応じた被写体画像の散乱線成分を推定して、被写体画像から散乱線成分を減算した処理後画像を取得する手法が開示されている。

また、非特許文献１には、人体の体厚分布に応じて放射線画像の散乱線成分を推定して除去する手法が開示されている。非特許文献１の画像処理方法によれば、被写体画像の画素値から推定した体厚分布に基づいて、入力された被写体画像に所定の関数を適用することにより被写体画像に含まれる散乱線の像を推定した推定散乱線画像を生成し、被写体画像から推定散乱線画像を減算することにより、入力された被写体画像から一次線画像を推定した推定一次線画像を生成する。さらに、生成した推定一次線画像に所定の関数を適用することによりさらなる推定散乱線画像を生成し、被写体画像からさらなる推定散乱線画像を減算してさらなる推定一次線画像を生成する処理を所定の収束条件下で収束するまで繰り返して、収束した推定散乱線画像を算出し、この推定散乱線画像を被写体画像から減算することにより最終的に散乱線成分を除去した処理後画像を取得することができる。また、非特許文献１には、被写体画像に含まれる散乱線の像を推定するための所定の関数を体厚に応じて調整する方法が開示されている。

特開平０２−２４４８８１号公報

Trotter, 他4名, "Thickness-dependent Scatter Correction Algorithm for Digital Mammography", Proc. SPIE Vol.4682, 2002年5月, p.469-478

ここで、被写体内部の肺野等の被写体の内部構造を反映した、詳細な体厚を求めるためには、実際に被写体を撮影して得られた被写体画像の画素値から被写体の厚さを算出することが好ましい。しかしながら、被写体画像には、被写体を透過して放射線検出器に直接照射された一次線の成分（一次線成分）と、被写体内で放射線が散乱した散乱線の成分（散乱線成分）が含まれている。

このため、特許文献１や非特許文献１のように、散乱線除去グリッド（グリッド）を用いないで撮影された放射線画像に画素値に基づいて体厚を推定する方法を適用した場合には、放射線画像に含まれた散乱線成分の影響によって、被写体の体厚分布を正確に推定することが難しい。また、散乱線成分の影響を避けるためにはグリッドを用いて被写体画像を撮影することが考えられるが、被写体の被曝量などの負担を低減するためにもグリッドを用いないで撮影した被写体画像から、体厚分布を正確に推定したいという要求がある。

本発明は上記事情に鑑みてなされたものであり、被写体に放射線を照射することにより撮影された放射線画像を解析して、被写体の体厚分布を正確に推定する画像解析処理を施すことができるようにすることを目的とする。

本発明による放射線画像解析装置は、被写体の放射線撮影により得られた被写体画像を解析して、被写体の体厚分布を推定する放射線画像解析装置であって、被写体画像を取得する画像取得部と、所定の体厚分布を有する被写体の仮想モデルを取得する仮想モデル取得部と、仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像とを合成した画像を、被写体の放射線撮影により得られる放射線画像を推定した画像である推定画像として生成する推定画像生成部と、推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布を修正する修正部と、修正された仮想モデルの体厚分布を被写体の体厚分布として決定する体厚分布決定部とを備えたことを特徴とする。

本発明による放射線画像解析方法は、放射線画像解析装置に実行される、被写体に放射線を照射することにより撮影された被写体画像を解析して、被写体の体厚分布を推定する放射線解析方法であって、被写体画像を取得する画像取得ステップと、所定の体厚分布を有する被写体の仮想モデルを取得する仮想モデル取得ステップと、仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像とを合成した画像を、被写体の放射線撮影により得られる放射線画像を推定した画像である推定画像として生成する推定画像生成ステップと、推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布を修正する修正ステップと、修正された仮想モデルの体厚分布を被写体の体厚分布として決定する体厚分布決定ステップとを有することを特徴とする。

なお、本発明による放射線画像解析方法をコンピュータに実行させるためのプログラムとして提供してもよい。

上記「体厚」は、照射された放射線の経路上における、空気領域を除いた被写体領域の厚さの総計を意味する。

「推定画像」は、仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像を合成した画像であると実質的に見なせるものであればよい。例えば、仮想モデルに推定一次線画像生成用の関数を適用して推定一次線画像を作成し、仮想モデルに推定散乱線画像生成用の関数を適用して推定散乱線画像とを別途生成してから合成してもよく、仮想モデルに推定画像生成用の関数を適用して推定画像を推定してもよい。

また、本発明による放射線画像解析装置において、仮想モデル取得部が、修正された体厚分布を有する仮想モデルをさらに取得し、推定画像生成部が、修正された体厚分布を有する仮想モデルから推定画像をさらに生成し、修正部が、さらに生成された推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布をさらに修正することが好ましい。

また、上記場合に、体厚分布決定部は、推定画像と被写体画像との違いが所定の閾値以下となる場合に、仮想モデルの体厚分布を被写体の体厚分布として決定することが好ましい。

また、本発明による放射線画像解析装置において、修正部は、推定画像と被写体画像の差分画像の画素値の絶対値の総和または差分画像の画素値の二乗和が小さくなるように、仮想モデルの体厚分布を修正することができる。

また、本発明による放射線画像解析装置において、修正部は、仮想モデルの一画素以上の部分領域ごとに仮想モデルの体厚分布を変動させて、推定画像と被写体画像との違いを小さくする部分領域の体厚を算出し、算出された各部分領域の体厚によって仮想モデルの体厚分布を修正することができる。

また、本発明において、「所定の体厚分布」とは、被写体の体厚分布に類似した分布である必要はなく、任意の体厚分布であってよい。例えば、所定の体厚分布を、被写体が人体の所定部位である場合に、被写体とは別の人体の同じ部位の体厚分布としてもよく、所定の体厚分布を一様分布としてもよい。

また、本発明による放射線画像解析装置において、所定の体厚分布が、被写体とは異なる比較用被写体を放射線撮影して得た比較用被写体画像と、比較用被写体を３次元撮影して得た３次元画像を取得し、取得された３次元画像の各位置において、比較用被写体画像の放射線経路に対応する直線上の比較用被写体の体厚を計測することにより作成されていてもよい。

また、本発明による放射線画像解析装置において、仮想モデルが、仮想モデルに含まれる構造物と構造物の配置と構造物の放射線に対する特性との少なくとも１つを表す特性情報をさらに有し、推定画像生成部が、特性情報に基づいて仮想モデルの各位置に対応する構造物に応じて推定画像を算出するためのパラメータを選択して推定画像を生成することが好ましい。

また、本発明による放射線画像解析装置において、推定画像生成部が、被写体画像に含まれる構造物と構造物の配置と構造物の放射線に対する特性とを表す特性情報を仮想モデルの特性情報として取得し、特性情報に基づいて、仮想モデルの各位置に対応する構造物に応じて推定画像を算出するためのパラメータを選択して推定画像を生成することが好ましい。

上記「特性情報」は、画像に含まれる構造物とその構造物の配置と、その構造物の放射線に対する特性を示すものであれば、任意の方法によって特定されたものであってよい。例えば、特性情報を被写体の肺野、骨、血管、臓器などの撮影対象である解剖学的構造物と、各解剖学的構造物の組成によって規定することができる。

また、本発明による放射線画像解析装置において、被写体の体厚分布を用いて被写体画像の散乱線を推定した散乱線情報を取得する散乱線情報取得部と、取得した散乱線情報に基づいて、被写体画像の散乱線の除去処理を行う散乱線除去部とをさらに備えることが好ましい。

また、本発明において、被写体画像は散乱線除去グリッドを用いないで撮影された画像であってもよい。

本発明によれば、被写体画像を取得し、所定の体厚分布を有する被写体の仮想モデルを取得し、仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像とを合成した画像を、被写体の放射線撮影により得られる放射線画像を推定した推定画像として生成し、推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布を修正し、修正された仮想モデルの体厚分布を被写体の体厚分布として決定する。このため、推定画像と被写体画像との違いに基づいて、推定画像が被写体画像に近似するように体厚分布を修正できるため、修正した仮想モデルの体厚分布を被写体の体厚分布とすることにより、正確に被写体画像の体厚分布を決定することができる。

また、仮想モデル取得部が、修正された体厚分布を有する仮想モデルをさらに取得し、推定画像生成部が、修正された体厚分布を有する仮想モデルから推定画像をさらに生成し、修正部が、さらに生成された推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布をさらに修正する場合には、修正された体厚分布を有する仮想モデルに基づいて、体厚分布の修正を繰り返すことにより、推定画像が被写体画像により近似するように体厚分布をより正確に修正できるため、修正した仮想モデルの体厚分布を被写体の体厚分布とすることにすることでさらに正確に被写体画像の体厚分布を決定することができる。

本発明の第１の実施形態に係る放射線画像解析装置を適用した放射線画像撮影システムの構成を示す概略ブロック図 本発明の第１の実施形態に係る放射線画像解析装置によって行われる処理を示すフローチャート 体厚分布の対応付けテーブルの一例を示す図 推定画像の生成方法の一例を説明するための図 推定画像の生成方法の別の一例を説明するための図 本発明の第３の実施形態に係る放射線画像解析装置を適用した放射線画像撮影システムの構成を示す概略ブロック図 本発明の第３の実施形態に係る放射線画像解析装置によって行われる処理を示すフローチャート 本発明の第４の実施形態に係る放射線画像解析装置によって行われる処理を示すフローチャート 本発明の第５の実施形態に係る放射線画像解析装置によって行われる処理を示すフローチャート

以下、図面を参照して本発明の実施形態について説明する。図１は本発明の第１の実施形態による放射線画像解析装置を適用した放射線画像撮影システムの構成を示す概略ブロック図である。図１に示すように、本実施形態による放射線画像撮影システムは、撮影装置１０と、システムを制御する制御装置２０と、画像解析装置３０（放射線画像解析装置）とを備える。

撮影装置１０は被写体ＫにＸ線を照射するＸ線源１２と、被写体Ｋを透過したＸ線を検出して被写体Ｋの放射線画像を取得する放射線検出器１４とを備える。なお、本実施形態においては、被写体Ｋと放射線検出器１４との間には、被写体Ｋを透過したＸ線のうち、被写体Ｋにより散乱した散乱線を除去するための散乱線除去グリッド（グリッド）は配置されない。

制御装置２０は、設定された撮影条件に従ってＸ線源１２を駆動制御する線源駆動制御部２２と、放射線検出器１４を制御し、被写体の放射線画像（被写体画像）を取得して記憶部４２に記憶する検出器制御部２４とを備える。

画像解析装置３０は、検出器制御部２４または後述の記憶部４２などから被写体Ｋを撮影して得られた被写体画像Ｉｋを取得する画像取得部３１と、初期体厚分布Ｔ_０（所定の体厚分布）を有する被写体Ｋの仮想モデルＭを取得する仮想モデル取得部３２と、仮想モデルＭに基づいて、仮想モデルの放射線撮影により得られる一次線画像を推定した推定一次線画像Ｉｐと、仮想モデルの放射線撮影により得られる散乱線画像を推定した推定散乱線画像Ｉｓとを合成した画像を、被写体Ｋの放射線撮影により得られる放射線画像を推定した推定画像Ｉｍとして生成する推定画像生成部３３と、推定画像Ｉｍと被写体画像Ｉｋに基づいて、推定画像Ｉｍと被写体画像Ｉｋの違いが小さくなるように仮想モデルＭの初期体厚分布Ｔ_０を修正する修正部３４と、修正された体厚分布Ｔ_ｎ−１（ｎは自然数）を被写体画像Ｉｋの体厚分布Ｔｋとして決定する体厚分布決定部３５と、決定した体厚分布Ｔｋ（ｘ，ｙ）に基づいて被写体画像Ｉｋに含まれるＸ線の散乱線成分を表す散乱線情報を取得する散乱線情報取得部３６と、散乱線情報取得部３６が取得した散乱線情報に基づいて、放射線検出器１４により取得された被写体画像Ｉｋの散乱線除去処理を行う散乱線除去部３７と、入力部３８と、表示部４０と、メモリおよびハードディスクなどの記憶媒体から構成され各種情報を記憶する記憶部４２とを備える。なお、入力部３８は、画像解析装置３０に対する操作者の各種入力を受け付ける。具体的には入力部３８は、キーボード、マウス、タッチパネル等からなる。表示部４０は、ＣＲＴ、液晶ディスプレイ等からなり、撮影装置１０により取得された放射線画像および後述する散乱線除去処理に必要な各種入力の補助を行う。

以上述べた画像取得部３１と、仮想モデル取得部３２と、推定画像生成部３３と、修正部３４と、体厚分布決定部３５と、散乱線情報取得部３６、散乱線除去部３７、入力部３８、表示部４０および記憶部４２は、例えば一般的なパーソナルコンピュータ等のコンピュータシステムから構成することができる。

画像解析装置３０は、検査対象者である被写体Ｋに放射線を照射することにより撮影された被写体画像Ｉｋを解析して、被写体Ｋの体厚分布Ｔｋを推定するものである。

画像解析装置３０の記憶部４２には、検出器制御部２４によって取得した被写体画像Ｉｋとその撮影線量、管電圧、Ｘ線源１２と放射線検出器１４の検出面との距離ＳＩＤ、線源のターゲットおよびフィルタの材質、撮影に使用される放射線検出器の種類並びにエアギャップ量（被写体から放射線検出器までの距離）などの撮影条件が記憶される。記憶部４２には、複数の撮影条件ごとに濃度値（画素値）と体厚とを対応付けた対応付けテーブルＬＵＴが予め作成されて記憶されている。また、記憶部４２には、初期体厚分布Ｔ_０（ｘ，ｙ）を有する被写体Ｋの仮想モデルＭが記憶される。また、記憶部４２には、各処理で必要とされる各種パラメータ、生成された画像（推定一次線画像、推定散乱線画像など）が適宜記憶されるものとする。なお、本明細書における体厚は、照射された放射線の経路上における空気領域を除いた被写体領域の厚さの総計を意味する。なお、ここでは検出器制御部２４より撮影条件が取得され記憶部４２に記憶されたが、線源駆動制御部２２によって取得されて記憶部４２に記憶されてもよい。

以下、図２に示すフローチャートを用いて本実施形態に係る画像解析装置３０によって実施される放射線画像解析処理の流れを説明する。

まず、画像取得部３１は、記憶部４２から、被写体Ｋである患者を放射線撮影して得られた被写体画像Ｉｋを取得する（Ｓ０１）。

次に仮想モデル取得部３２は、記憶部４２から、初期体厚分布Ｔ_０（ｘ，ｙ）を有する被写体Ｋの仮想モデルＭを取得する（Ｓ０２）。仮想モデルＭは、初期体厚分布Ｔ_０（ｘ，ｙ）に従った体厚がｘｙ平面上の各位置に対応付けられた被写体Ｋを仮想的に表すデータである。また、仮想モデルＭに含まれる構造物（ここでは肺野、骨、臓器などの解剖学的構造物）と構造物の配置と、構造物の放射線に対する特性などを示す特性情報は、比較用被写体の胸腹部の肺野、骨などの解剖学的構造物の配置および組成に基づいて設定されている。

また、仮想モデルＭの初期体厚分布Ｔ_０（ｘ，ｙ）には任意の分布とされてよいが、本実施形態においては、仮想モデル取得部３２によって初期体厚分布Ｔ_０が生成されて取得される。仮想モデル取得部３２は、被写体Ｋの撮影線量、管電圧、ＳＩＤなどの撮影条件を取得し、記憶部４２から被写体Ｋの撮影条件に応じた画素値と体厚の対応付けテーブルＬＵＴを取得する。図３に画素値と体厚の対応付けテーブルＬＵＴの一例を示す。そして、記憶部４２から比較用被写体（人体）を放射線撮影して得られた画像データを取得し、対応付けテーブルＬＵＴに基づいて、比較用被写体の画像データの各画素の画素値に対応する体厚を特定することにより比較用被写体の画像データの体厚分布を取得する。そして、仮想モデル取得部３２は、画像データの体厚分布を仮想モデルＭの初期体厚分布Ｔ_０（所定の体厚分布）として取得する。なお、初期体厚分布Ｔ_０は、本実施形態のように仮想モデルＭの取得処理の際に生成されてもよく、仮想モデルＭの取得処理に先立って予め設定されていてもよい。

次いで推定画像生成部３３は、被写体画像と同等の撮影条件で仮想モデルＭを撮影した場合に得られる推定一次線画像Ｉｐと、被写体画像と同等の撮影条件で仮想モデルＭを撮影した場合に得られる推定散乱線画像Ｉｓとを合成した推定画像Ｉｍを生成する（Ｓ０３）。図４Ａ、図４Ｂは、推定画像Ｉｍの生成方法を説明するための図である。

図４Ａに示すように、推定画像生成部３３は、仮想モデルＭを、被写体画像Ｉｋと同等の撮影条件で仮想モデルＭを撮影した場合に得られる推定一次線画像Ｉｐを、下記式（２）に従って生成し、生成した推定一次線画像Ｉｐを用いて、式（３）に従って推定散乱線画像Ｉｓを生成する。そして、推定画像生成部３３は、式（４）に示すように推定一次線画像Ｉｐと推定散乱線画像Ｉｓを合成することにより、推定画像Ｉｍを生成する（Ｓ０３）。なお、推定一次線画像Ｉｐと推定散乱線画像Ｉｓを１回目に作成する際には、推定式（２）、式（３）において初期体厚分布Ｔ_０（ｘ，ｙ）が用いられる（式（２）、（３）においてｎ＝１である）。

ここで、（ｘ，ｙ）は被写体画像Ｉｋの画素位置の座標、Ｉｐ（ｘ，ｙ）は画素位置（ｘ，ｙ）における推定一次線画像、Ｉｓ（ｘ，ｙ）は画素位置（ｘ，ｙ）における推定散乱線画像、Ｉｏ（ｘ，ｙ）は画素位置（ｘ，ｙ）における線量、Ｉｍ（ｘ，ｙ）は画素位置（ｘ，ｙ）における推定画像、μは被写体の線減弱係数、Ｋ_ｓ(ｘ，ｙ，Ｔｎ（ｘ’，ｙ’），θ_{ｘ’，ｙ’})は画素位置（ｘ，ｙ）における被写体厚に応じた点拡散関数(Point Spread Function)を表す畳みこみカーネルである。なお、線量Ｉｏ（ｘ，ｙ）は、被写体が存在しないと仮定した際に放射線が検出器で検出される線量であり、Ｘ線源１２と放射線検出器１４の検出面との距離（ＳＩＤ）、管電圧、撮影線量に応じて変化する。また、θ_{ｘ’，ｙ’}は、管電圧などの撮影条件や仮想モデルＭの特性情報によって特定されるパラメータを表している。

なお、推定画像Ｉｍは、仮想モデルＭを放射線撮影した場合に得られると推定される画像であればよく、推定一次線画像Ｉｐと推定散乱線画像Ｉｓを合成した画像と実質的に見なせるものであればよい。例えば、図３Ｂに示すように、式（２）〜（４）に替えて下記式（５）を用いて、一次線成分と散乱線成分を合わせたカーネルを畳み込み積分して推定画像Ｉｍを生成してもよい。ここで、Ｋ_ｐ＋ｓ（ｘ，ｙ，Ｔ_ｎ−１（ｘ’，ｙ’），θ_{ｘ’，ｙ’})は、一次線成分と散乱線成分を合わせた点拡散関数を表すカーネルである。また、放射線撮影により得られた画像から推定一次線画像および推定散乱線画像を合成した推定画像を生成可能であれば、任意のモデル関数を用いてよい。

図２のフローチャートに従って、続く処理を説明する。続いて、体厚分布決定部３５は、被写体画像Ｉｋと推定画像Ｉｍとの違いが終了条件を満たすか否かを判定する（Ｓ０４）。ここでは、式（６）、式（７）に示すように、下記の被写体画像Ｉｋと推定画像Ｉｍとの違いを表すエラー値Ｖ_{ｅｒｒｏｒ}を定義し、終了条件としてエラー値Ｖ_{ｅｒｒｏｒ}が閾値以下であるか否かを判定する。また、式（７）に示すように、被写体画像Ｉｋから推定画像Ｉｍを減算した差分画像Ｉｄの各画素値の２乗和をエラー関数ｆ_{ｅｒｒｏｒ}として規定する。なお、終了条件として、被写体画像Ｉｋと推定画像Ｉｍとの違いが許容可能な程度に十分小さくなったことを判定可能なあらゆる判定手法を適用可能である。

また、上記例に限定されず、エラー関数ｆ_{ｅｒｒｏｒ}を、被写体画像Ｉｋと推定画像Ｉｍとの違いを表すあらゆる方法で規定することができる。例えば、下記式（８）に示すように、被写体画像Ｉｋから推定画像Ｉｍを減算した差分画像Ｉｄの各画素値の絶対値の総和をエラー関数ｆ_{ｅｒｒｏｒ}としてもよい。

エラー値Ｖ_{ｅｒｒｏｒ}が終了条件を満たさない場合には（Ｓ０４，Ｎｏ）、修正部３４は、体厚分布Ｔ_ｎ−１（ｎ＝１の場合には、初期体厚分布Ｔ_０）を修正する修正処理を行う（Ｓ０５）。

体厚分布Ｔ_ｎ−１の修正処理を行うために、被写体画像Ｉｋと推定画像Ｉｍの違いが小さくなるように体厚分布Ｔ_ｎ−１の各位置の修正値を取得できる任意の方法を適用可能である。本実施形態では、仮想モデルＭの一画素以上の部分領域ごとに、仮想モデルＭの体厚分布Ｔ_ｎ−１を変動させて、推定画像Ｉｍと被写体画像Ｉｋとの違いを小さくする部分領域の体厚を算出する処理を実施する。そして、算出された各部分領域の体厚によって仮想モデルの体厚分布を修正する。

具体的には、本実施形態は、最急降下法を用いて体厚分布Ｔ_ｎ−１の体厚の修正値を求めるものとする。下記式（９）、（１０）を用いて、仮想モデルＭの画素のうち、Ｔ_ｎ−１（ｘ，ｙ）において１つの特定の座標の体厚のみを変動させて、エラー関数ｆ_{ｅｒｒｏｒ}の一次偏微分（勾配）に基づいて繰り返しｄＴ_ｎ−１（ｘ，ｙ）を算出することにより、エラー関数ｆ_{ｅｒｒｏｒ}の出力値を最小化することができる。そして、エラー関数ｆ_{ｅｒｒｏｒ}の出力値を最小化した際の、１つの特定の座標の体厚を、その特定の座標の体厚の修正値として決定する。また、他の画素についても同様に、それぞれ体厚の修正値を求めることにより、各画素の体厚分布を修正し、修正した体厚分布Ｔｎを取得する。

ただし、式（９）において、αは、体厚の更新速度を表すパラメータである更新係数である。式（１０）に示すＫ_ｐ＋ｓの微分値部分の算出方法の一例として、例えば、Ｔ_ｎ−１（ｘ，ｙ）に極めて小さい値ｄｔを加えたとき値の変化を式（１１）によって算出して、式（１０）のＫ_ｐ＋ｓの値とすることができる。なお、式（１）〜（１２）において、同じ要素には同じ符号を付して、説明を省略する。被写体画像Ｉｋと推定画像Ｉｍの違いを表すエラー値Ｖ_{ｅｒｒｏｒ}を最小化するあらゆる最適化手法を適用可能であり、例えば、シンプレックス法や最急降下法、共役勾配法を用いることができる。

修正された体厚分布Ｔ_ｎを取得すると、体厚分布決定部３５は、ｎの値を１つ増加して更新し（ｎ＝ｎ＋１とする）、仮想モデル取得部３２は修正された体厚分布Ｔ_ｎを取得する（Ｓ０２）。そして、取得された体厚分布Ｔ_ｎに対して、推定画像生成部３３と体厚分布決定部３５はＳ０３乃至Ｓ０４の処理をそれぞれ先述同様に実行する。そして、被写体画像Ｉｋと推定画像Ｉｍとの違いを示すエラー値Ｖ_{ｅｒｒｏｒ}が終了条件を満たすまで、上記同様に、体厚分布Ｔ_ｎの修正処理（Ｓ０５）と、修正された体厚分布Ｔ_ｎを有する仮想モデルＭの取得処理（Ｓ０２）と、体厚分布Ｔ_ｎを用いた新たな推定画像Ｉｍの生成処理（Ｓ０３）と、新たに生成された推定画像Ｉｍと被写体画像Ｉｋとの違いが終了条件を満たすかを判定する処理（Ｓ０４）とが繰り返される。

一方、体厚分布決定部３５は、エラー値Ｖ_{ｅｒｒｏｒ}が終了条件を満たしていることを判定した場合には（Ｓ０４，Ｙｅｓ）、終了条件を満たした際にエラー値Ｖ_{ｅｒｒｏｒ}に用いられた体厚分布Ｔ_ｎを被写体画像Ｉｋの体厚分布Ｔｋとして決定して本実施形態における画像解析処理を終了する（Ｓ０６）。

なお、本実施形態では、図３に示す体厚分布決定処理（Ｓ０６）が実施された後に、ユーザの被写体画像Ｉｋに対する散乱線除去処理の指示入力を受け付けた場合には、画像解析装置３０のオプションの機能として散乱線情報取得処理と散乱線除去処理を実施する。ここでは、散乱線情報取得部３６は、取得した体厚分布Ｔｋを適用して、式（２）によって被写体画像Ｉｋの推定一次線画像を取得し、式（３）に従って、被写体画像Ｉｋ（ｘ，ｙ）の推定散乱線画像Ｉｓ（ｘ，ｙ）を取得する。そして、散乱線除去部３７は、被写体画像Ｉｋの推定散乱線画像Ｉｓ（ｘ，ｙ）を被写体画像Ｉｋ（ｘ，ｙ）から減算することにより散乱線による影響を除去した処理後画像を生成して記憶部４２に記憶する。処理後画像は、ユーザ指示に応じて適宜表示部４０に表示される。なお、終了条件の判定（Ｓ０４）において、エラー値Ｖ_{ｅｒｒｏｒ}が終了条件を満たした際には、被写体画像Ｉｋと推定画像Ｉｍとの違いは許容可能な程度に小さくなっており、推定一次線画像Ｉｐと被写体画像Ｉｋの一次線画像との違いも小さくなっている。このため、被写体画像Ｉｋから散乱線を除去した処理後画像を取得したい場合には、終了条件を満たしたエラー値Ｖ_{ｅｒｒｏｒ}の算出に用いられた推定一次線画像Ｉｐを被写体画像Ｉｋの一次線画像、すなわち被写体画像Ｉｋから散乱線除去処理を行った後の処理後画像として取得してもよい。

本実施形態によれば、仮想モデルＭを放射線撮影した場合に得られると推定される画像である、推定一次線画像Ｉｐと推定散乱線画像Ｉｓを合成した推定画像Ｉｍを生成して、推定画像Ｉｍと被写体画像Ｉｋとの違いが小さくなるように、仮想モデルＭの体厚分布を修正しているため、推定画像Ｉｍと被写体画像Ｉｋとの違いに基づいて、推定画像Ｉｍが被写体画像Ｉｋに近似するように正確に体厚分布Ｔ_ｎを修正でき、修正した仮想モデルＭの体厚分布Ｔ_ｎを被写体Ｋの体厚分布Ｔｋとすることにすることにより、正確に被写体画像Ｉｋの体厚分布Ｔｋを決定することができる。また、従来方法ではグリッドを用いないで撮影された画像からは散乱線成分による影響が大きく正確な体厚分布を算出することが難しかったことに対し、本実施形態の手法によれば、推定画像Ｉｍが被写体画像Ｉｋに近似するように正確に体厚分布Ｔ_ｎを修正して、被写体Ｋの体厚分布Ｔｋとして決定しているため、被写体画像Ｉｋがグリッドを用いないで撮影された画像であっても、従来方法よりも正確に体厚分布Ｔｋを得ることができる。

また、本実施形態に示すように、仮想モデル取得部３２が、修正された体厚分布Ｔ_ｎを有する仮想モデルＭをさらに取得し、推定画像生成部３３が、修正された体厚分布Ｔ_ｎを有する仮想モデルＭから推定画像Ｉｍをさらに生成し、修正部３４が、さらに生成された推定画像Ｉｍと被写体画像Ｉｋの違いが小さくなるように仮想モデルＭの体厚分布Ｔ_ｎをさらに修正する場合には、修正された体厚分布Ｔ_ｎを有する仮想モデルに基づいて、体厚分布Ｔの修正を繰り返すことにより、推定画像Ｉｍが被写体画像Ｉｋにより近似するように体厚分布Ｔをより正確に修正できるため、修正した仮想モデルＭの体厚分布Ｔ_ｎ＋１を被写体Ｋの体厚分布Ｔｋとすることで、さらに正確に被写体画像Ｉｋの体厚分布Ｔｋを決定することができる。

また、本実施形態に示すように、体厚分布決定部３５が、推定画像Ｉｍと被写体画像Ｉｋとの違いが許容可能な程度に十分小さくなった場合に、仮想モデルＭの体厚分布Ｔ_ｎを被写体Ｋの体厚分布Ｔｋとして決定している場合には、推定画像Ｉｍと被写体画像Ｉｋとが近似する体厚分布となるように、繰り返し体厚分布を修正して、非常に正確に被写体画像の体厚分布を決定することができる。また、体厚分布決定部３５は、推定画像Ｉｍと被写体画像Ｉｋとの違いが閾値以下となったか否かを判定することにより、推定画像Ｉｍと被写体画像Ｉｋとの違いが許容可能な程度に十分小さくなったか否かを好適に判別して、推定画像Ｉｍと被写体画像Ｉｋとが近似する体厚分布となるように、繰り返し体厚分布を修正して、非常に正確に被写体画像の体厚分布を決定することができる。

また、本実施形態において、修正部３４が、推定画像と被写体画像の違い画像の画素値の絶対値の総和または違い画像の画素値の二乗和が小さくなるように、仮想モデルの体厚分布を修正しているため、好適に推定画像Ｉｍと被写体画像Ｉｋとの違いの大きさを判定することができる。

また、本実施形態に示すように、修正部３４が、仮想モデルＭの１画素以上の部分領域ごとに、仮想モデルＭの体厚分布Ｔ_ｎ−１のうち１つの部分領域の体厚を変動させて、推定画像Ｉｍと被写体画像Ｉｋとの違いを最小化する場合のその１つの部分の体厚を算出し、算出された各部分の体厚によって仮想モデルＭの体厚分布を修正する場合には、各画素の体厚の修正値を正確に算出でき、好適に修正された体厚分布Ｔ_ｎを取得することができる。

本実施形態によれば、決定された被写体Ｋの体厚分布Ｔｋを用いて被写体画像の散乱線を推定した散乱線情報を取得する散乱線情報取得部３６と、取得した散乱線情報に基づいて、被写体画像の散乱線の除去処理を行う散乱線除去部３７を備えているため、精度よく散乱線除去処理を行った処理後画像を取得することができる。

また、本実施形態の第２の実施形態として、推定画像生成部３３が、被写体画像Ｉｋに含まれる構造物と、構造物の配置と、構造物の放射線に対する特性を表す特性情報を仮想モデルＭの特性情報として取得し、特性情報に基づいて、仮想モデルＭの各位置に対応する構造物に応じて推定画像Ｉｍを算出するためのパラメータを選択して推定画像Ｉｍを生成してもよい。例えば式（１２）で示すように、特性情報に基づいて、仮想モデルＭから式（２）を用いて推定一次線画像Ｉｐを作成する際の式（２）の線減弱係数を、各位置における構造物（構造物の組成）に応じて切り換えて用いることが考えられる。放射線撮影された画像において一次線成分や散乱線成分は、被写体の骨、臓器の種類、臓器などに含まれる空洞の有無などの被写体に含まれる構造物と、構造物の空間的な位置によって、放射線撮影された画像の各位置において複雑に変動する。このため、被写体画像Ｉｋの特性情報を仮想モデルＭの特性情報として取得し、仮想モデルＭの各位置に（仮想的に）含まれる構造物に応じて、推定一次線画像や推定散乱線画像などに用いるパラメータを適切に選択することにより、構造物に起因する一次線成分や散乱線成分の誤差を低減して、より正確に推定一次線画像Ｉｐ、推定散乱線画像Ｉｓを生成することができる。

なお、式（３）で示すＫ_ｓのパラメータθ_{ｘ’，ｙ’}についても、構造物ごとに異なるθ_{ｘ’，ｙ’}の値を設定し、各位置の構造物に応じて各位置に適用されるθ_{ｘ’，ｙ’}を異ならせてもよい。また、被写体画像Ｉｋと同一の被写体Ｋを撮影したＣＴ画像やＭＲＩ画像などの３次元画像を取得し、取得したＣＴ画像やＭＲＩ画像から被写体画像Ｉｋの特性情報を測定して取得してもよい。同じ被写体Ｋの３次元画像を用いて特性情報を取得した場合には、臓器や骨の空間的な位置などの情報についても正確に取得できる。

初期体厚分布（所定の体厚分布）として種々の体厚分布を用いてよい。例えば、初期体厚分布を一様分布とすることもできるが、計算負荷の観点から、ある程度被写体Ｋに近いと推測される体厚分布を初期体厚分布とすることが好ましい。この観点から、例えば、同一の被写体Ｋについて過去に測定した体厚分布Ｔを初期体厚分布Ｔ_０とすることが好ましい。この場合には、姿勢や経時変化による被写体Ｋの違いなどの微小な修正を行うだけで体厚分布の決定をすることができ、計算負荷を軽減することができる。

また、初期体厚分布Ｔ_０が、被写体とは異なる比較用被写体を放射線撮影して得た比較用被写体画像と、比較用被写体を３次元撮影して得たＣＴ画像やＭＲＩ画像などの３次元画像を取得し、取得された３次元画像の各位置において、比較用被写体画像の放射線経路に対応する直線上の比較用被写体の体厚を計測することにより作成されていてもよい。３次元画像を用いることより、初期体厚分布を正確なものとすることができるため、被写体の姿勢などによる微小な修正を行うだけで体厚分布の決定をすることができ、計算負荷を軽減することができる。

また、比較用被写体に関する体厚分布を用意する場合、複数の比較用被写体について、各比較用被写体を測定して得た体厚分布と比較用被写体の体格を表す体格情報とを互いに対応付けたセットを複数用意し、被写体の体格情報を取得し、被写体の体格情報と類似する体格情報を有する比較用被写体の体厚分布を特定し、特定された比較用被写体の体厚分布を、被写体画像の初期体厚分布Ｔ_０とすることが考えられる。この場合には、被写体と体格が類似する比較用被写体の体厚分布を初期体厚分布とすることにより、初期体厚分布を被写体の体厚分布と類似するものとできる可能性が高くなり、被写体の姿勢などによる微小な修正を行うだけで体厚分布の決定をすることができ、計算負荷を軽減することができる。また、この被写体の体格情報と比較用被写体の体格情報は、ユーザ入力により入力されたものであってもよく、被写体画像と比較用被写体の放射線画像のそれぞれの濃度ヒストグラムのヒストグラム幅（濃度値の最大値と最小値の差）など被写体画像と比較用被写体から抽出された体格情報であってもよい。また、被写体画像と比較用被写体から抽出可能な体格情報であれば、あらゆる情報を体格情報とすることができ、例えば、被写体の所定部位の長さなどを被写体の体格情報として用いてもよい。

また、被写体の放射線撮影の前後等に超音波センサまたはデジタルメジャーなどの距離を測定可能な測定装置によってＸ線源側の被写体表面と検出器との距離を測定し、測定装置またはユーザ入力などからこの測定値を取得して、取得した被写体の体表面と検出器との距離を、初期体厚分布を決定するための指標値として用いてもよい。この場合には、例えば、被写体の体表面と検出器との距離を上記被写体の体格情報として用いることができる。また、初期体厚分布を被写体の体表面と検出器との距離の一様分布としてもよい。

また、推定一次線画像Ｉｐと推定散乱線画像Ｉｓを生成できる種々の方法を適用してよい。例えば、式（２）および（３）に替えて、例えば、加藤秀起、「ディジタルX線画像の後処理による散乱線成分除去法」, 日本放射線技術学会雑誌第62巻第9号, 2006年9月, p.1359-1368に記載されたように、モンテカルロシミュレーション法を用いて推定一次線画像Ｉｐと推定散乱線画像Ｉｓを生成してもよい。また、モンテカルロシミュレーション法を用いる場合に、仮想モデルＭに含まれる構造物と構造物の配置と構造物の放射線に対する特性とを表す情報である特性情報を用いることが好ましい。この場合には、より高精度に推定一次線画像Ｉｐと推定散乱線画像Ｉｓを生成することができる。

また、上記各実施形態において、実際に使用されたと考えられる撮影条件を取得する撮影条件取得部をさらに備えてもよい。管電圧や撮影線量、Ｘ線源と検出器の面との距離ＳＩＤは、被写体の体型や診断目的や放射線撮影システムが設置される施設の環境などに応じて変更されるものである。このため、撮影条件取得部を備えて被写体画像Ｉｋの撮影条件を取得し、推定画像生成部３３が、取得した撮影条件に応じて、推定画像Ｉｍの生成に用いられる、撮影条件によって変動するパラメータ（例えば、上記式（３）、（５）におけるθ_{ｘ’，ｙ’}など）を選択し、選択したパラメータを用いて推定画像Ｉｍの生成処理（Ｓ０３）を行うことが好ましい。

撮影条件取得部は、撮影条件を取得可能であればあらゆる方法で撮影条件を取得してよく、例えば、撮影条件取得部をユーザに入力された撮影条件または被写体が存在しない位置において検出器で検出される画素値から算出した撮影条件を取得するものとしてもよい。この場合、被写体が存在しない領域である素抜け領域の濃度値と撮影線量とを対応づけたテーブルを記憶部４２に記憶しておき、素抜け領域の濃度値に基づいてこのテーブルを参照して撮影線量を取得すればよい。また、撮影条件取得部は、実際に被写体画像の撮影に適用された撮影条件を取得可能な種々の方法を適用できる。

例えば、撮影条件取得部は、Ｘ線源と検出器の面との距離ＳＩＤを任意の方法で取得してよく、例えば超音波センサまたはデジタルメジャーなどの距離を測定可能な測定装置によって測定されたＸ線源と検出器との距離をＳＩＤとして取得することができる。また、撮影条件取得部は、Ｘ線源１２と放射線検出器１４との間であって放射線検出器１４から既知の距離離間した位置に配置された３次元マーカを放射線撮影して得られた放射線画像を取得し、この放射線画像中の３次元マーカの位置や散乱線成分を解析してＳＩＤを算出してもよい。

また、撮影条件取得部は、撮影線量を任意の方法で取得してよく、例えば、面積線量計などの測定機器を用いて測定された線量を、放射線検出器１４に入射する撮影線量として取得してもよい。また、撮影条件取得部は、厚さが既知のアクリルモデルを被写体とともに撮影し、取得された放射線画像におけるアクリルモデルの部分の濃度に基づいて、撮影線量を取得するようにしてもよい。この場合、アクリルモデルの濃度と撮影線量とを対応づけたテーブルを記憶部４２に記憶しておき、アクリルモデルの濃度に基づいてこのテーブルを参照して撮影線量を取得すればよい。なお、撮影条件は、放射線画像撮影システムが設置される施設に応じて決まっていることが多い。このため、実際の撮影時の撮影条件が不明である場合は、施設に応じた撮影条件を使用すればよい。

また、放射線検出器１４の種類ごとに被写体画像Ｉｋの画素値が変動する可能性があるため、撮影条件に放射線検出器１４の種類を特定する情報を含むようにし、撮影線量と、管電圧と、ＳＩＤと、放射線検出器１４の種類との組合せごとに、予め測定して得られた濃度値（画素値）と放射線検出器１４に到達した線量とを対応付けた線量特定テーブルを予め作成して記憶部４２に記憶しておくことが好ましい。そして、撮影条件取得部は、被写体画像Ｉｋの撮影に用いられた放射線検出器１４の種類を取得し、取得した放射線検出器１４の種類に対応する線量特定テーブルを参照して、各位置の濃度値に対応する線量を特定して撮影線量として用いてもよい。

また、修正部３４が、取得した撮影条件に応じて、推定画像Ｉｍの生成に用いられる、撮影条件によって変動するパラメータ（例えば、上記式（１０）、（１１）におけるθ_{ｘ’，ｙ’}など）を選択し、選択したパラメータを用いて推定画像Ｉｍの体厚分布の修正処理（Ｓ０５）を行うことが好ましい。この場合には、被写体画像Ｉｋの撮影条件に応じて、撮影条件によって変動するパラメータを適切に設定して推定画像Ｉｍを生成することができるため、推定画像Ｉｍをより正確に推定して生成することができる。このため、結果として、被写体Ｋの体厚分布をより正確に決定することができる。

また、上記各実施形態において、画像取得部３１は、グリッドを用いないで撮影された被写体画像Ｉｋを取得するものとしたが、これに限定されず、画像取得部３１は、グリッドを用いて被写体Ｋを放射線撮影した放射線画像に対して、グリッドに起因する縞模様を除去する処理を施した画像を被写体画像Ｉｋとして取得してもよい。グリッドに起因する縞模様を除去する処理は、グリッドに起因する縞模様を除去可能な種々の手法によって行われたものであってよく、例えば特開２０１２−２０３５０４号公報に記載された手法などが参照できる。

また、画像取得部３１が、グリッドを用いて被写体Ｋを放射線撮影した放射線画像に対して、グリッドに起因する縞模様を除去する処理を施した画像を被写体画像Ｉｋとして取得する場合には、グリッドの種類を特定するグリッド情報ごとに散乱線透過率Ｔｓおよび一次線透過率Ｔｐを対応付けたテーブルを予め作成して記憶部４２に記憶しておくことが好ましい。そして、推定画像生成部３３は、グリッド情報ごとに散乱線透過率Ｔｓおよび一次線透過率Ｔｐを対応付けたテーブルに基づいて、被写体の撮影に用いられたグリッド情報に対応する散乱線透過率Ｔｓ^ｋおよび一次線透過率Ｔｐ^ｋを特定し、条件式（４）に換えて下記条件式（４−１）を用いて推定画像を生成することが好ましい。この場合には、仮想モデルの推定一次線画像を、被写体画像に使用されたものと同じ種類のグリッドによる一次線の吸収を反映したものとすることができ、仮想モデルの推定散乱線画像を、被写体画像に使用されたものと同じ種類のグリッドによる散乱線の吸収を反映したものとすることができる。このように、被写体画像に使用されたグリッドに基づいて、推定画像を同グリッドによる散乱線と一次線の吸収を反映したものにすることにより、被写体画像に使用されたグリッドに放射線が吸収されることに起因する体厚分布の誤差を低減することができる。なお、グリッド情報を、グリッド比、グリッド密度、収束型か平行型か、収束型の場合の集束距離およびインタースペース素材（アルミニウム、ファイバー、ベークライト等）に含まれる１つ以上の要素の任意の組合せとすることができる。

続いて、図５、図６を用いて、本発明の第３の実施形態について説明する。図５は、本発明の第３の実施形態に係る放射線画像解析装置を適用した放射線画像撮影システムの構成を示す概略ブロック図であり、図６は、本発明の第３の実施形態に係る放射線画像解析装置によって行われる処理を示すフローチャートである。

図５に示すように、第３の実施形態に係る画像解析装置３０は、グリッドを使用することなく撮影を行うことにより取得された被写体画像に対して、実際にグリッドを使用して撮影を行った場合と同様の散乱線を除去する効果を付与するように、被写体画像に対して画像処理を施す散乱線解析部４３を備える。散乱線解析部４３は、後述の特性取得部３９と、散乱線情報取得部３６と、散乱線除去部３７により構成される。第３の実施形態は、画像解析装置３０が、仮想的なグリッドの特性である仮想グリッド特性を取得する特性取得部３９を更に備え、散乱線情報取得部３６が散乱線情報として被写体画像Ｉｋに含まれるＸ線の散乱成分を表す散乱成分情報をさらに取得し、散乱線除去部が特性取得部３９に取得された仮想グリッド特性および散乱線情報取得部３６が取得した散乱成分情報に基づいて、放射線検出器１４により取得された被写体画像の散乱線除去処理を行う点が上記各実施形態と異なる。なお、第３の実施形態においては、図１に示す画像解析装置３０と共通する各構成要素を備え、上記の相違点を除いて各構成要素の機能および処理も概ね共通するため、第１の実施形態と異なる部分を中心に説明し、共通部分については適宜説明を省略する。

特性取得部３９は、操作者による入力部３８からの入力により仮想グリッド特性を取得する。第３の実施形態においては、仮想グリッド特性は、仮想グリッドについての散乱線透過率Ｔｓ、および被写体Ｋを透過して直接放射線検出器１４に照射される一次線の透過率（一次線透過率）Ｔｐとする。なお、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは０〜１の間の値をとる。

特性取得部３９は、散乱線透過率Ｔｓおよび一次線透過率Ｔｐの値の入力を直接受け付けることにより仮想グリッド特性を取得してもよいが、第３の実施形態においては、グリッドの種類を表すグリッド情報、被写体についての情報（被写体情報）、および被写体画像Ｉｋの取得時の撮影条件の少なくとも１つの指定を受け付けることにより、仮想グリッド特性、すなわち散乱線透過率Ｔｓおよび一次線透過率Ｔｐを取得する。なお、以下、特性取得部３９に用いられる撮影条件を、特性取得用撮影条件と記載する。

ここで、グリッド情報とは、グリッド比、グリッド密度、収束型か平行型か、収束型の場合の集束距離、インタースペース素材（アルミニウム、ファイバー、ベークライト等）等の、グリッドの種類を特定する情報の少なくとも１つを含む。散乱線透過率Ｔｓおよび一次線透過率Ｔｐはグリッドの種類に応じて異なるものとなる。このため、グリッド情報に関して、各種グリッド情報の少なくとも１つと仮想グリッド特性とを対応づけたテーブルが記憶部４２に記憶されている。

被写体情報は、胸部、腹部および頭部等の被写体の種類を含む。ここで、被写体画像Ｉｋの撮影時には、一般的に撮影部位に応じて使用するグリッドの種類が決められており、グリッドの種類に応じて散乱線透過率Ｔｓおよび一次線透過率Ｔｐが異なるものとなる。このため、被写体情報に関して、各種被写体情報と仮想グリッド特性とを対応づけたテーブルが記憶部４２に記憶されている。

特性取得用撮影条件は、撮影時の撮影距離（ＳＩＤ）、撮影線量、管電圧、線源のターゲットおよびフィルタの材質、並びに撮影に使用される放射線検出器の種類等のうちの少なくとも１つを含む。ここで、被写体画像Ｉｋの撮影時には、一般的に撮影条件に応じて使用するグリッドの種類が決められており、グリッドの種類に応じて散乱線透過率Ｔｓおよび一次線透過率Ｔｐが異なるものとなる。このため、特性取得用撮影条件に関して、各種特性取得用撮影条件と仮想グリッド特性とを対応づけたテーブルが記憶部４２に記憶されている。なお、特性取得用撮影条件は、仮想グリッド特性の取得に必要なパラメータを含むものであれば、体厚分布の決定に用いる撮影条件（体厚分布決定用撮影条件）や後述の散乱線情報の取得に用いる散乱線情報取得用撮影条件と同じものであってもよく、異なるものであってもよい。

特性取得部３９は、記憶部４２に記憶されたテーブルを参照して、入力部３８から入力されたグリッド情報、被写体情報および特性取得用撮影条件の少なくとも１つに基づいて、仮想グリッド特性を取得する。なお、グリッド情報、被写体情報および特性取得用撮影条件は、入力部３８の直接の入力を受け付ければよいが、各種グリッド情報、各種被写体情報および各種特性取得用撮影条件のリストを表示部４０に表示し、リストからのグリッド情報、被写体情報および特性取得用撮影条件の少なくとも１つの選択を受け付けることにより、グリッド情報、被写体情報および特性取得用撮影条件の入力を行うようにしてもよい。

また、第３の実施形態においては，散乱線除去処理は、後述するように被写体画像Ｉｋを周波数分解することにより行われる。第３の実施形態においては、仮想グリッド特性は、周波数分解による得られる被写体画像Ｉｋの複数の周波数帯域のそれぞれについて取得される。このため、上記テーブルにおける仮想グリッド特性は、複数の周波数帯域のそれぞれに対応づけられたものとなっている。

また、グリッド情報、被写体情報および特性取得用撮影条件のすべてと仮想グリッド特性とを対応づけたテーブルを記憶部４２に記憶しておき、グリッド情報、被写体情報および特性取得用撮影条件のすべてに基づいて仮想グリッド特性を取得するようにしてもよい。この場合、テーブルは、各種グリッド情報、各種被写体情報および各種特性取得用撮影条件と、仮想グリッド特性とを対応づけた少なくとも４次元のテーブルとなる。

なお、グリッドを使用することによって増加する照射線量の増加率である露出倍数、グリッドを使用した場合と使用しない場合とのコントラストの比率であるコントラスト改善係数、および一次Ｘ線透過率の散乱Ｘ線透過率に対する比率である選択度は、グリッドの特性を表す特性値であり、これらの特性値から散乱線透過率Ｔｓおよび一次線透過率Ｔｐを算出することができる。このため、特性取得部３９において、露出倍数、コントラスト改善係数および選択度の少なくとも１つの指定を受け付けることにより、仮想グリッド特性、すなわち散乱線透過率Ｔｓおよび一次線透過率Ｔｐを算出して取得するようにしてもよい。

また、第３の実施形態において画像解析装置３０は、仮想グリッド特性のみならず、散乱成分情報にも基づいて散乱線除去処理を行う。このため、散乱線情報取得部３６は散乱線情報としてさらに散乱成分情報を取得する。第３の実施形態においては、散乱成分情報は、例えば被写体Ｋが胸部であれば、縦隔が存在する被写体画像Ｉｋの中央部分ほど散乱線が多く、肺野が存在する周辺部ほど散乱線が少ないという、被写体画像Ｉｋにおける散乱線含有率分布とする。

散乱線情報取得部３６は、撮影により取得された被写体画像Ｉｋを解析することにより、散乱成分情報すなわち散乱線含有率分布を取得する。被写体画像Ｉｋの解析は、被写体画像Ｉｋの撮影時における照射野情報、被写体情報および撮影条件に基づいて行う。

照射野情報とは、照射野絞りを用いて撮影を行った場合における、被写体画像Ｉｋに含まれる照射野の位置および大きさに関する照射野分布を表す情報である。被写体情報とは、上述した胸部、腹部および頭部等の被写体の種類に加えて、被写体の被写体画像Ｉｋ上での位置、被写体の組成の分布、被写体の大きさおよび被写体の厚さ等に関する情報である。散乱線情報取得部３６に用いられる撮影条件とは、撮影線量（管電流×照射時間）、管電圧、撮影時の撮影距離（ＳＩＤ）、エアギャップ量（被写体から放射線検出器までの距離）、および放射線検出器の特性等に関する情報である。なお、以下、散乱線情報取得部３６に用いられる撮影条件を、散乱線情報取得用撮影条件と記載する。散乱線情報取得用撮影条件は、散乱線情報の取得に必要なパラメータを含むものであれば、体厚分布の決定に用いた撮影条件や特性取得用の撮影条件と同じものであってもよく、異なるものであってもよい。

これらの照射野情報、被写体情報および散乱線情報取得用撮影条件は、被写体画像Ｉｋに含まれる散乱線の分布を決める要因となっている。例えば、散乱線の大小は照射野の大きさにより左右され、被写体の厚さが大きいほど散乱線は多くなり、被写体と放射線検出器との間に空気が存在すると散乱線が減少する。したがって、これらの情報を用いることにより、より正確に散乱線含有率分布を取得することができる。

散乱線情報取得部３６は、被写体Ｋの体厚分布Ｔｋ（ｘ，ｙ）から、式（１３）、（１４）にしたがって一次線像および散乱線像を算出し、算出した一次線像および散乱線像から式（１５）に基づいて、散乱線含有率分布Ｓ（ｘ，ｙ）を算出する。なお、散乱線含有率分布Ｓ（ｘ，ｙ）は０〜１の間の値をとる。
Ip(x,y) = Io(x,y)×exp(-Tk(x,y)×μ) （１３）
Is(x,y) = Io(x,y)＊Sσ(Tk(x,y)) （１４）
S(x,y) = Is(x,y)/(Is(x,y)+Ip(x,y)) （１５）

ここで、（ｘ，ｙ）は被写体画像Ｉｋの画素位置の座標、Ｉｐ（ｘ，ｙ）は画素位置（ｘ，ｙ）における一次線像、Ｉｓ（ｘ，ｙ）は画素位置（ｘ，ｙ）における散乱線像、Ｉｏ（ｘ，ｙ）は画素位置（ｘ，ｙ）における被写体表面への入射線量、μは被写体の線減弱係数、Ｓσ（Ｔｋ（ｘ，ｙ））は画素位置（ｘ，ｙ）における被写体厚に応じた散乱の特性を表す畳みこみカーネルである。なお、先述の式（２）と式（１３）は公知の指数減弱則に基づく式であり、式（１４）は「J M Boon et al, An analytical model of the scattered radiation distribution in diagnostic radiolog, Med. Phys. 15(5), Sep/Oct 1988」（参考文献１）に記載された手法に基づく式である。なお、上記式（１３）、（１４）において、被写体表面への入射線量Ｉｏ（ｘ，ｙ）は、どのような値を定義してもＳ（ｘ，ｙ）を算出する際に除算によってキャンセルされるため、例えば値を１とする等、任意の値としてもよい。

ここで、式（１４）における＊は畳みこみ演算を表す演算子である。カーネルの性質は、被写体の厚さの他に、照射野の分布、被写体の組成の分布、撮影時の照射線量、管電圧、撮影距離、エアギャップ量、および放射線検出器の特性等によっても変化する。参考文献１に記載された手法によれば散乱線は一次線に対する位置拡張関数（point spread function、式（１４）におけるSσ（Ｔｋ（ｘ，ｙ））の畳みこみにより近似することができる。なお、Sσ（Ｔｋ（ｘ，ｙ））は、照射野情報、被写体情報および散乱線情報取得用撮影条件等に応じて実験的に求めることができる。

第３の実施形態においては、撮影時の照射野情報、被写体情報および散乱線情報取得用撮影条件に基づいてSσ（Ｔｋ（ｘ，ｙ））を算出してもよいが、各種照射野情報、各種被写体情報および各種散乱線情報取得用撮影条件とSσ（Ｔｋ（ｘ，ｙ））とを対応づけたテーブルを記憶部４２に記憶しておき、撮影時の照射野情報、被写体情報および散乱線情報取得用撮影条件に基づいて、このテーブルを参照してSσ（Ｔｋ（ｘ，ｙ））を求めるようにしてもよい。なお、Sσ（Ｔｋ（ｘ，ｙ））をＴｋ（ｘ，ｙ）にて近似するようにしてもよい。

散乱線除去部３７は、仮想グリッド特性および散乱成分情報に基づいて、被写体画像Ｉｋにおける散乱線と見なせる周波数帯域の周波数成分を低減させることにより、散乱線除去処理を行う。このため、散乱線除去部３７は、被写体画像Ｉｋを周波数分解して複数の周波数帯域毎の周波数成分を取得し、少なくとも１つの周波数成分のゲインを低減する処理を行い、処理済みの周波数成分およびこれ以外の周波数成分を合成して、散乱線除去処理済みの被写体画像Ｉｋを取得する。なお、周波数分解の手法としては、被写体画像Ｉｋを多重解像度変換する手法の他、ウェーブレット変換、フーリエ変換等、公知の任意の手法を用いることができる。

散乱線除去部３７は、仮想グリッド特性である散乱線透過率Ｔｓおよび一次線透過率Ｔｐ、並びに散乱線含有率分布Ｓ（ｘ，ｙ）から、周波数成分を変換する変換係数Ｒ（ｘ，ｙ）を下記の式（１６）により算出する。
R(x,y) = S(x,y)×Ts + (1-S(x,y))×Tp （１６）

散乱線透過率Ｔｓおよび一次線透過率Ｔｐ、並びに散乱線含有率分布Ｓ（ｘ，ｙ）は０〜１の間の値となるため、変換係数Ｒ（ｘ，ｙ）も０〜１の間の値となる。散乱線除去部３７は、変換係数Ｒ（ｘ，ｙ）を複数の周波数帯域のそれぞれについて算出する。

なお、以降の説明において、被写体画像Ｉｋの画素値をＩｋ（ｘ，ｙ）、周波数分解により得られる周波数成分画像をＩｋ（ｘ，ｙ，ｒ）、周波数合成をＩｋ（ｘ，ｙ）＝ΣｒＩｋ（ｘ，ｙ，ｒ）、周波数帯域毎の変換係数をＲ（ｘ，ｙ，ｒ）、周波数帯域毎の散乱線透過率および一次線透過率をＴｓ（ｒ）、Ｔｐ（ｒ）で表すものとする。なお、ｒは周波数帯域の階層を表し、ｒが大きいほど低周波であることを表すものとする。したがって、Ｉｋ（ｘ，ｙ，ｒ）は、ある周波数帯域の周波数成分画像となる。散乱線含有率分布Ｓ（ｘ，ｙ）は被写体画像Ｉｋについてのものをそのまま用いればよいが、散乱線透過率Ｔｓおよび一次線透過率Ｔｐと同様に周波数帯域のそれぞれについて取得するようにしてもよい。

第３の実施形態においては、周波数成分毎に変換係数Ｒ（ｘ，ｙ，ｒ）を算出し、周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）に対して対応する周波数帯域の変換係数Ｒ（ｘ，ｙ，ｒ）を乗算して、周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）の画素値を変換し、変換係数Ｒ（ｘ，ｙ，ｒ）が乗算された周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）（すなわち、Ｉｋ（ｘ，ｙ，ｒ）×Ｒ（ｘ，ｙ，ｒ））を周波数合成して処理済みの被写体画像Ｉｋ′（ｘ，ｙ）を取得する。したがって、散乱線除去部３７において行われる処理は、下記の式（１７）により表される。なお、変換係数Ｒ（ｘ，ｙ，ｒ）は０〜１の間の値となるため、周波数成分Ｉｋ（ｘ，ｙ，ｒ）に対して対応する周波数帯域の変換係数Ｒ（ｘ，ｙ，ｒ）を乗算することにより、その周波数成分の画素位置（ｘ，ｙ）における画素値すなわちゲインが低減されることとなる。
Ik′(x,y)=Σr{Ik(x,y,r)×R(x,y,r)}
=Σr{Ik(x,y,r)×(S(x,y)×Ts(r)+(1-S(x,y))×Tp(r))} （１７）

ここで、第３の実施形態においては、被写体画像Ｉｋを６つの周波数帯域に周波数分解するものとし、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは６つの周波数帯域について取得されるものとする。この場合、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは、例えば下記式（１８）に示す値となる。なお、式（１８）では右側ほど低周波数帯域の値を表すものとする。
Ts={0.7, 0.7, 0.7, 0.7, 0.3, 0.2}
Tp={0.7, 0.7, 0.7, 0.7, 0.7, 0.7} （１８）

式（１８）に示すように、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは、高周波数帯域（ｒ＝１〜４）では同一の値であるが、低周波数帯域（ｒ＝５〜６）においては、散乱線透過率Ｔｓの方が低い値となる。これは、グリッドは散乱線の周波数成分が支配的である低周波帯域ほどその除去率が高いが、一次線については除去率の周波数依存性が小さいからである。

例えば、胸部を放射線撮影して得られた被写体画像であれば、散乱線の含有率が高い縦隔部および肺野の周囲において、式（１６）、（１８）に基づいて算出した変換係数の値が小さくなり、より大きく画素値が低減されることとなる。したがって、このように算出した変換係数を用いて式（１７）に示す処理を行うことにより取得された処理済みの被写体画像Ｉｋ′（処理後画像Ｉｋ′）においては、使用が想定されるグリッドの種類に応じて散乱線成分が除去されたものとなる。

なお、散乱線除去部３７においては、下記のようにして被写体画像Ｉｋの散乱線を除去するようにしてもよい。まず、上記と同様に周波数合成をＩｋ（ｘ，ｙ）＝ΣｒＩｋ（ｘ，ｙ，ｒ）で表すとすると、散乱線除去部３７は、周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）を、下記の式（１９）により、散乱線含有率分布Ｓ（ｘ，ｙ）を用いて、散乱線像（散乱線成分）Ｉｓ（ｘ，ｙ，ｒ）と一次線像（一次線成分）Ｉｐ（ｘ，ｙ，ｒ）とに分解する。
Is(x,y,r)= S(x,y)×Ik(x,y,r)
Ip(x,y,r)=(1-S(x,y))×Ik(x,y,r) （１９）

さらに散乱線除去部３７は、下記の式（２０）により、散乱線成分Ｉｓ（ｘ，ｙ，ｒ）および一次線成分Ｉｐ（ｘ，ｙ，ｒ）のそれぞれに対して、仮想グリッド特性である散乱線透過率Ｔｓ（ｒ）および一次線透過率Ｔｐ（ｒ）を適用して画像変換し、変換された散乱線成分Ｉｓ′（ｘ，ｙ，ｒ）および一次線成分Ｉｐ′（ｘ，ｙ，ｒ）を算出する。
Is′(x,y,r)=Is(x,y,r)×Ts(r)=S(x,y)×Ik(x,y,r)×Ts(r)
Ip′(x,y,r)=Ip(x,y,r)×Tp(r)=(1-S(x,y))×Ik(x,y,r)×Tp(r) （２０）

そして下記の式（２１）により、Ｉｓ′（ｘ，ｙ，ｒ）および一次線成分Ｉｐ′（ｘ，ｙ，ｒ）を周波数合成して、処理済みの被写体画像Ｉｋ（ｘ，ｙ）′を算出する。
Ik′(x,y)=Σr{Is′(x,y,r)+Ip′(x,y,r)}
=Σr{S(x,y)×Ik(x,y,r)×Ts(r)+(1-S(x,y))×Ik(x,y,r)×Tp(r)}
=Σr{Ik(x,y,r)×(S(x,y)×Ts(r)+(1-S(x,y))×Tp(r))} （２１）

次いで、第３の実施形態において行われる処理について説明する。図６は第３の実施形態において行われる処理を示すフローチャートである。撮影装置１０において取得された被写体画像Ｉｋが画像解析装置３０に入力されると（Ｓ２１）、第１の実施形態と同様に被写体画像Ｉｋの体厚分布を決定する処理が行われる（Ｓ２２）。なおＳ２２は、図２のＳ０２〜Ｓ０６に示す処理に対応する。続いて、特性取得部３９がグリッド情報、被写体情報および撮影条件の少なくとも１つの入力部３８からの入力を受け付け（グリッド情報等入力、Ｓ２３）、仮想グリッド特性、すなわち散乱線透過率Ｔｓおよび一次線透過率Ｔｐを取得する（Ｓ２４）。

また、散乱線情報取得部３６が被写体画像Ｉｋを解析し（Ｓ２５）、散乱線成分情報、すなわち散乱線含有率分布Ｓ（ｘ，ｙ）を取得する（Ｓ２６）。一方、散乱線除去部３７は、被写体画像Ｉｋを周波数分解する（Ｓ２７）。なお、Ｓ２３，Ｓ２４の処理、Ｓ２５，Ｓ２６の処理およびＳ２７の処理は並列に行ってもよく、Ｓ２５，Ｓ２６の処理を先に行ってもよく、Ｓ２７の処理を先に行ってもよい。

そして、散乱線除去部３７は、上記式（１６）により周波数帯域毎の変換係数Ｒ（ｘ，ｙ，ｒ）を算出し（Ｓ２８）、変換係数Ｒ（ｘ，ｙ，ｒ）により周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）を変換する（Ｓ２９）。そして変換された周波数成分画像Ｉｋ′（ｘ，ｙ，ｒ）を周波数合成して、処理後画像Ｉｋ′を取得し（Ｓ３０）、処理を終了する。なお、処理後画像Ｉｋ′は表示部４０に表示されて診断に供されるか、外部の画像サーバに送信されて保存される。

このように、第３の実施形態においては、画像解析装置３０が、グリッドを使用しないで撮影された被写体画像に対して、実際にグリッドを使用して撮影を行った場合と同様の散乱線を除去する効果を付与するように、被写体画像に対して画像処理を施す散乱線解析部４３を備え、被写体画像Ｉｋの撮影時に散乱線を除去するために使用が想定されるグリッドの特性である仮想グリッド特性を取得し、さらに散乱性成分情報を取得し、仮想グリッド特性および散乱線成分情報に基づいて、被写体画像Ｉｋの散乱線除去処理を行うようにしたものである。このため、実際に使用する散乱線除去グリッドと同様の散乱線除去の効果を被写体画像Ｉｋに付与することができる。また、被写体画像Ｉｋの画質を、様々な種類の散乱線除去グリッドを使用して撮影を行うことにより取得した被写体画像の画質に近づけることができる。

また、収束型のグリッドを使用した場合、放射線の斜入により被写体画像Ｉｋに濃度ムラが発生するおそれがあるが、第３の実施形態においては、グリッドを使用しないで撮影された被写体画像に対して散乱線を除去する処理を実施することにより、放射線の斜入による濃度ムラが発生しなくなるため、より高画質の処理後画像Ｉｋ′を取得することができる。

次いで、本発明の第４の実施形態について説明する。なお、第４の実施形態においては、画像解析装置の構成は第３の実施形態と同様であり、行われる処理のみが異なるため、ここでは装置についての詳細な説明は省略する。上記第３の実施形態においは、被写体画像Ｉｋを周波数分解し、変換後の周波数成分画像を周波数合成して処理後画像Ｉｋ′を取得しているが、第４の実施形態においては、被写体画像Ｉｋから除去すべき周波数帯域の周波数成分を抽出し、抽出された周波数成分に対して散乱線除去処理を行い、処理後の周波数成分を被写体画像Ｉｋに加減算することにより、処理後画像Ｉｋ′を取得するようにした点が第３の実施形態と異なる。

第４の実施形態においては、散乱線除去部３７は、以下のように処理を行う。まず、第３の実施形態と同様に、被写体画像Ｉｋを周波数分解して周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）を取得し、下記の式（２２）により、除去のための変換係数Ｒ′（ｘ，ｙ，ｒ）を周波数帯域毎に算出する。
R′(x,y) = S(x,y)×(1-Ts(r)) + (1-S(x,y))×(1-Tp(r)) （２２）
そして、下記の式（２３）により、周波数帯域毎の除去成分ΔＩｋ（ｘ，ｙ，ｒ）を算出する。
ΔIk(x,y,r)=Ik(x,y,r)×R′(x,y,r)
=Ik(x,y,r)×{S(x,y)×(1-Ts(r))+(1-S(x,y))×(1-Tp(r))} （２３）
そして、除去成分ΔＩｋ（ｘ，ｙ，ｒ）を周波数合成し、被写体画像Ｉｋ（ｘ，ｙ）から周波数合成した除去成分ΣｒΔＩｋ（ｘ，ｙ，ｒ）を減算して、処理済みの被写体画像Ｉｋ′（ｘ，ｙ）を取得する。
Ik′(x,y)=Ik(x,y)-ΣrΔIk(x,y,r)
=Ik(x,y)-Σr{Ik(x,y,r)×{S(x,y)×(1-Ts(r))+(1-S(x,y))×(1-Tp(r))}}
（２４）

次いで、第４の実施形態において行われる処理について説明する。図７は第４の実施形態において行われる処理を示すフローチャートである。なお、Ｓ３１〜Ｓ３７までの処理は第３の実施形態におけるＳ２１〜Ｓ２７の処理と同一であるため、ここでは詳細な説明は省略する。

Ｓ３７に続き、散乱線除去部３７は、上記の式（２２）により除去のための変換係数Ｒ′（ｘ，ｙ，ｒ）を周波数帯域毎に算出し（Ｓ３８）、上記の式（２３）により周波数帯域毎の除去成分ΔＩｋ（ｘ，ｙ，ｒ）を算出する（Ｓ３９）。そして、除去成分ΔＩｋ（ｘ，ｙ，ｒ）を周波数合成し（Ｓ４０）、周波数合成した除去成分ΣｒΔＩｋ（ｘ，ｙ，ｒ）を被写体画像Ｉｋから減算して、処理後画像Ｉｋ′を取得し（Ｓ４１）、処理を終了する。なお、処理後画像Ｉｋ′は表示部４０に表示されて診断に供されるか、外部の画像サーバに送信されて保存される。

なお、上記第４の実施形態において、散乱線除去部３７は下記のようにして被写体画像Ｉｋの散乱線を除去するようにしてもよい。まず、上記と同様に周波数分解をＩｋ（ｘ，ｙ）＝ΣｒＩｋ（ｘ，ｙ，ｒ）で表すとすると、周波数成分画像Ｉｋ（ｘ，ｙ，ｒ）を、上記の式（１９）により、散乱線含有率分布Ｓ（ｘ，ｙ）を用いて、散乱線成分Ｉｓ（ｘ，ｙ，ｒ）と一次線成分Ｉｐ（ｘ，ｙ，ｒ）とに分解する。さらに、下記の式（２５）により、散乱線成分Ｉｓ（ｘ，ｙ，ｒ）および一次線成分Ｉｐ（ｘ，ｙ，ｒ）のそれぞれに対して、仮想グリッド性能としての散乱線透過率Ｔｓ（ｒ）、一次線透過率Ｔｐ（ｒ）を適用して画像変換し、散乱線除去成分ΔＩｓ（ｘ，ｙ，ｒ）および一次線除去成分ΔＩｐ（ｘ，ｙ，ｒ）を算出する。
ΔIs(x,y,r)=Is(x,y)×(1-Ts(r))=S(x,y)×Ik(x,y)×(1-Ts(r))
ΔIp(x,y,r)=Ip(x,y)×(1-Tp(r))=(1-S(x,y))×Ik(x,y)×(1-Tp(r)) （２５）

そして下記の式（２６）により、散乱線除去成分ΔＩｓ（ｘ，ｙ，ｒ）および一次線除去成分ΔＩｐ（ｘ，ｙ，ｒ）を周波数合成し、周波数合成した散乱線除去成分ΣｒΔＩｓ（ｘ，ｙ，ｒ）および一次線除去成分ΣｒΔＩｐ（ｘ，ｙ，ｒ）を被写体画像Ｉｋから減算して、処理済みの被写体画像Ｉｋ′（ｘ，ｙ）を算出する。
Ik′(x,y)=Ik(x,y) -Σr(ΔIs(x,y,r)+ΔIp(x,y,r)) （２６）

なお、上記第３および第４の実施形態において、被写体画像Ｉｋが放射線検出器への入射線量に比例した画素値を持つ、放射線量に対してリニアな空間において散乱線除去処理を行い、その後に対数変換を行って、人間の視覚と比例する対数リニアな空間に変換することが好ましい。

また、上記第３および第４の実施形態においては、特性取得部３９において、仮想グリッド特性として、散乱線透過率Ｔｓおよび一次線透過率Ｔｐを取得しているが、散乱線透過率Ｔｓおよび一次線透過率Ｔｐの一方のみを取得するようにしてもよい。

また、上記第３および第４の実施形態においては、グリッドを使用せずに撮影を行うことにより取得した被写体画像Ｉｋに対して散乱線除去処理を行っているが、グリッドを使用して撮影を行うことにより取得した被写体画像を処理の対象としてもよい。この場合、被写体画像に対して、グリッドに起因する縞模様を除去する処理を施し、その後、散乱線除去処理を行う。また、かかる散乱線除去処理においては、所望のグリッドである第１のグリッドが使用されて撮影された放射線画像（第１グリッド使用画像）を取得し、所望の仮想的なグリッドに対応する仮想グリッド特性を取得し、取得された第１グリッド使用画像中の散乱線量と一次線線量を、取得した仮想グリッド特性に対応するグリッド（取得した仮想グリッド特性の散乱線透過率と一次線透過率を有するグリッド）に対応する散乱線量と一次線線量になるように変換すればよい。また、第１のグリッドと仮想グリッド特性に対応するグリッドは、どちらの散乱線除去効果が大きくてもよく、目的や事情に応じて任意に選択されてよい。なお、グリッドに起因する縞模様を除去する処理としては、例えば特開２０１２−２０３５０４号公報に記載された手法を用いることができる。

また、グリッド無しで撮影された被写体画像に１つの仮想グリッド特性（第１の仮想グリッド特性）を適用して散乱線除去処理を行った処理後画像を第３および第４の実施形態における散乱線除去処理の対象としてもよい。この場合、第１の仮想グリッド特性と該第１の仮想グリッド特性が適用された処理後画像である第１処理後画像とを取得し、第１の仮想グリッド特性とは異なる所望の仮想的なグリッドに対応する第２仮想グリッド特性を取得する。そして、第２の仮想グリッド特性に基づいて、第１処理後画像中の散乱線量と一次線線量を、第２の仮想グリッド特性に対応する散乱線量と一次線線量になるように変換すればよい。また、第１の仮想グリッド特性と第２の仮想グリッド特性は、どちらの散乱線除去効果が大きくてもよく、目的や事情に応じて任意に選択されてよい。

このような処理を行うことにより、例えば、グリッド比が３：１であるグリッドを用いて撮影した被写体画像（またはグリッド無しで撮影された被写体画像に第１の仮想グリッド特性に基づいて散乱線除去処理を行った第１処理後画像）に基づいて、使用したグリッドとは異なる、１０：１のグリッド比をもつグリッドを用いて撮影したかのような処理後画像を仮想的に取得することが可能となる。また、逆に、グリッド比が１０：１であるグリッドを用いて撮影した被写体画像（またはグリッド無しで撮影された被写体画像に第１の仮想グリッド特性に基づいて散乱線除去処理を行った第１処理後画像）に基づいて、使用したグリッドとは異なる、３：１のグリッド比をもつグリッドを用いて撮影したかのような処理後画像を仮想的に取得することが可能となる。これらの場合には、被写体の撮影を繰り返さなくても、容易にグリッド比を変換した画像を取得することができるため、意図しないグリッド比をもつグリッドで撮影された被写体画像や上記の第１処理後画像から、所望のグリッド比のグリッドを使用して散乱線除去処理を行った処理後画像を得ることができる。このため、被写体の再撮影無しに、異なる度合いで散乱線除去処理を行った処理後画像を観察したいという要望に応えることができる。

具体的な方法としては、例えば第３の実施形態において、変換前のグリッドに対応する変換前グリッド情報と、変換後のグリッドに対応する変換後グリッド情報の組合せごとに、式（１４）における散乱の特性を表すＳσを対応付けたテーブルを記憶部４２に記憶しておく。なお、かかるテーブルにおけるＳσは、変換前のグリッドによる散乱の特性を変換後のグリッドによる散乱の特性に相対的に変換できるように、あらかじめ実験等により求めたものとする。そして、散乱線情報取得部３６は、実際に使用された使用グリッド（または仮想的なグリッド）に対応する第１グリッド情報を変換前グリッド情報として取得し、所望の仮想的なグリッドに対応する第２グリッド情報を変換後グリッド情報として取得し、上記テーブルに基づいて第１グリッド情報と第２グリッド情報に対応するＳσを取得する。そして、式（１３）と式（１４）を用いて、Io（ｘ，ｙ）を例えば１とし、取得したＳσを用いて一次線像Ｉｐ（ｘ，ｙ）と、散乱線像Is（ｘ，ｙ）とをそれぞれ算出する。そして、算出された一次線像Ｉｐ（ｘ，ｙ）と、散乱線像Is（ｘ，ｙ）とに基づいて、式（１５）により散乱線含有率分布Ｓ（ｘ，ｙ）を算出すればよい。

また、散乱線除去部３７は、式（１８）に示す周波数帯域毎の散乱線透過率Ｔｓおよび一次線透過率Ｔｐについて、実際に使用された使用グリッド（または仮想的なグリッド）に対応する第１のグリッド特性（一次線透過率Ｔｐ１、散乱線透過率Ｔｓ１）と、所望の仮想的なグリッドに対応する第２の仮想グリッド特性（一次線透過率Ｔｐ２、散乱線透過率Ｔｓ２）を取得し、変換前の第１のグリッドによる散乱の特性を変換後の第２のグリッドによる散乱の特性に相対的に変換するために、Ｔｐ２／Ｔｐ１を式（１８）における一次線透過率Ｔｐとして取得し、Ｔｓ２／Ｔｓ１を式（１８）における散乱線透過率Ｔｓとして取得する。そして、散乱線除去部３７は、取得した散乱線透過率Ｔｓ（=Ｔｓ２／Ｔｓ１）と一次線透過率Ｔｐ（=Ｔｐ２／Ｔｐ１）を式（１６）に適用して変換係数Ｒを求めて、この変換係数Ｒを用いて第３の実施形態と同様に散乱線除去処理を行えばよい。なお、上記式（１６）において、変換係数Ｒ（ｘ，ｙ）は、第１のグリッド特性の散乱線透過率Ｔｓ１に対して第２のグリッド特性の散乱線透過率Ｔｓ２が大きい場合には、１より大きい値となる場合がある。

また、第１のグリッド特性と第２のグリッド特性は、任意の方法で取得されてよい。例えば、グリッド情報ごとに予め実験等により求めたグリッド特性（一次線透過率Ｔｐ、散乱線透過率Ｔｓ）を対応付けたテーブルを用意して記憶部４２に記憶しておく。そして、散乱線除去部３７が、第１および第２のグリッド情報を取得し、このテーブルに基づいて、第１および第２のグリッド情報に対応する第１のグリッド特性と第２のグリッド特性を取得するようにしてもよい。また、第１および第２のグリッド特性を入力部３８からのユーザ入力に基づいて取得するようにしてもよい。なお、グリッド情報は、入力部３８からの入力により取得してもよく、例えば特開２００３−２６００５３号公報に記載されたように、グリッドの種類に応じた突起をグリッドに形成しておき、その突起を検出することによりグリッド情報を取得するようにしてもよい。

一方、撮影部位によっては、散乱線除去グリッドを使用しないで撮影した被写体画像を用いて画像の観察をしたいという要求がある。このような部位を撮影することにより取得した被写体画像に対して、上記第３および第４の実施形態の散乱線除去処理を行うことは好ましくない。このため、撮影部位に応じて、第３の実施形態の散乱線除去処理のオン／オフを撮影部位に応じて切り替えるようにすることが好ましい。なお、撮影部位の情報は、操作者が入力することにより取得してもよく、撮影フローの制御を行う周知のコンソールＰＣ（不図示）に入力される撮影依頼から自動的に取得してもよく、撮影後にシステムが被写体画像に付帯して保存する情報を利用してもよい。また、このような情報が取得できない場合には、被写体画像に対して部位認識処理を行うことにより取得してもよい。この場合、部位に応じて処理のオン／オフを対応づけたテーブルを記憶部４２に記憶しておき、このテーブルを参照して処理のオン／オフを切り替えるようにすればよい。

また、上記第３および第４の実施形態において、処理後画像Ｉｋ′および処理前の被写体画像Ｉｋの双方を表示し、いずれの被写体画像Ｉｋを診断に使用するかを選択できるようにしてもよい。

また、病気の治癒状況あるいは進行状況の診断を行うために、過去の被写体画像を用いて経時比較観察を行う場合がある。このような場合において、散乱線除去グリッドを使用せずに撮影を行うことにより取得した被写体画像（第１の被写体画像とする）と、散乱線除去グリッドを使用して撮影を行うことにより取得した被写体画像（第２の被写体画像とする）とを比較する場合には、第１の被写体画像に対してグリッドに起因する縞模様を除去する処理を施した際の処理条件に応じて、第３または第４の実施形態の散乱線除去処理の条件を修正し、第１および第２の被写体画像の画質を一致させるようにすることが好ましい。

なお、第３、４の実施形態における被写体画像Ｉｋの体厚分布を決定する一連の処理（図６のＳ２２、図７のＳ３２）は、少なくとも、仮想モデル取得部３２が所定の体厚分布を有する被写体の仮想モデルを取得する処理と、推定画像生成部３３が仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像とを合成した画像を、被写体の放射線撮影により得られる被写体画像を推定した画像である推定画像として生成する処理と、修正部３４が推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布を修正する処理と、体厚分布決定部３５が修正された仮想モデルの体厚分布を被写体の体厚分布として決定する処理を少なくとも含むもの（本発明に係る体厚分布を決定する処理）であれば、被写体Ｋの体厚分布Ｔｋを決定するいかなる処理としてもよい。例えば、体厚分布を決定する一連の処理（図６のＳ２２、図７のＳ３２）を、第２の実施形態に対応する処理を含むようにしてもよい。また、第３、４の実施形態における被写体画像Ｉｋの体厚分布を決定する一連の処理（図６のＳ２２、図７のＳ３２）は、散乱線情報取得部３６が被写体Ｋの体厚分布Ｔｋ（ｘ，ｙ）から、式（１３）、（１４）にしたがって一次線像および散乱線像を算出する処理に先立って実行されるものであれば、任意のタイミングで実行されてよい。

また、被写体画像Ｉｋを取得してから、本発明による被写体画像の体厚分布を決定する一連の処理（例えば図２のＳ０２〜Ｓ０６参照）を実施し、決定した被写体画像の体厚分布を用いて散乱線除去処理などのその他所望の画像処理を行って処理後画像を生成して表示する場合に、被写体画像取得から処理後画像表示までの所要時間（被写体画像を取得してから、被写体画像の体厚分布を決定する処理と、被写体画像に対する散乱線除去処理などのその他所望の画像処理とを実施し、得られた処理後画像を表示するまでの時間）をできるだけ短縮化することが好ましい。このために、上記各実施形態の被写体画像の体厚分布を決定する一連の処理において、画像取得部３１が被写体画像を所定のサイズに縮小した縮小画像を生成して取得し、仮想モデル取得部３２が縮小画像に対応するサイズの仮想モデルを取得し、推定画像生成部３３が縮小画像に対応するサイズの推定画像を生成し、修正部３４が推定画像と縮小画像の違いが小さくなるように仮想モデルの体厚分布を修正し、体厚分布決定部３５が修正された仮想モデルの体厚分布を例えば被写体画像Ｉｋのサイズなどの所望のサイズに拡大して被写体の体厚分布として決定してもよい。この場合には、被写体画像の体厚分布を決定する処理に含まれる少なくとも１つ以上の処理（特に、推定画像生成部３３が縮小画像に対応するサイズの推定画像を生成する処理、修正部３４が推定画像と縮小画像の違いが小さくなるように仮想モデルの体厚分布を修正する処理など）の処理負荷を低減できる可能性が高く、被写体画像取得から処理後画像表示までの所要時間を低減して、ユーザの待ち時間を低減することによりユーザの観察作業の効率化を支援することができる。さらに、医療機関などにおいて、複数の被写体を放射線撮影した複数の被写体画像を取得し、これらの複数の被写体画像に対してそれぞれ上記各実施形態に示すような被写体画像の体厚分布を決定する処理と散乱線除去処理などの所望の画像処理を順次実施し、得られた処理後画像をユーザが順に表示させて観察する状況下において、上記効果が大いに有益であり、ユーザの作業時間を著しく低減できることができる。なお、縮小画像は、縮小画像から取得した特徴情報とモデル画像の特徴情報を用いて類似を好適に判別可能な解像度を維持する範囲内で解像度をできるだけ小さくすることが好ましい。具体的には、撮影対象部位毎に特徴的・代表的な臓器や組織が識別・モデル選択可能な範囲で圧縮率を変えればよい。

また、医療機関などにおいて、複数の被写体を放射線撮影した複数の被写体画像を取得し、これらの複数の被写体画像に対してそれぞれ本発明における被写体画像の体厚分布を決定する処理（例えば図２のＳ０２〜Ｓ０６参照）と、散乱線除去処理などの所望の画像処理を順次実施し、得られた処理後画像をユーザが順に表示させて観察する場合には、被写体画像取得から処理後画像表示までの所要時間を有効活用することが好ましい。一例として、上記各実施形態による被写体画像取得から処理後画像表示までの所要時間中に、処理後画像の表示に先だって、画像観察の参照になる情報を表示することが考えられる。

例えば、第１の実施形態の画像解析装置の変形例である第５の実施形態として、被写体画像を取得した後、被写体画像Ｉｋの体厚分布Ｔｋを決定する一連の処理（図２のＳ０２〜Ｓ０６）の処理期間中に、散乱線情報取得部３６が予め設定された仮体厚分布を用いて被写体画像の散乱線情報を取得し、散乱線除去部３７が仮体厚分布を用いて取得した散乱線情報に基づいて被写体画像Ｉｋの散乱線除去処理を実施して仮処理後画像を生成し、表示部４０に仮処理後画像を表示させる処理（仮処理後画像の生成表示処理）の一部または全部を実施するようにしてもよい。

図８に、第５の実施形態のフローチャートを示す。第５の実施形態の画像解析装置は、被写体画像Ｉｋの体厚分布Ｔｋを決定する一連の処理（図２のＳ０２〜Ｓ０６）の処理期間中に、散乱線情報取得部３６が予め設定された仮体厚分布に基づいて被写体画像Ｉｋに含まれるＸ線の散乱線を表す散乱線情報を取得し、散乱線除去部３７が仮体厚分布に基づいて取得された散乱線情報に基づいて被写体画像の散乱線除去処理を実施して仮処理後画像を生成し、表示部４０に仮処理後画像を表示させる点において第１の実施形態と異なる。なお、その他の部分については、第５の実施形態における画像解析装置は、図１に示す画像解析装置３０と共通する各構成要素を備え、各構成要素の機能および処理も概ね共通するため、第１の実施形態と異なる部分を中心に説明し、共通部分については適宜説明を省略する。

第５の実施形態では、図８に示すように、画像取得部３１が上記第１の実施形態のＳ０１と同様に被写体画像Ｉｋを取得し（Ｓ５１）、上記第１の実施形態（図２のＳ０２〜Ｓ０６）と同様に被写体画像Ｉｋの体厚分布Ｔｋを決定する処理（Ｓ５２）を実施し、上記第１の実施形態のオプションの例と同様に条件式（１）、（２）に基づいて、決定された体厚分布Ｔｋを用いて散乱線情報（一次線像Ｉｐと散乱線像Ｉｓ）を取得し、上記第１の実施形態のオプションの例と同様に、取得した散乱線像Ｉｓを被写体画像Ｉｋから減算して散乱線除去処理（Ｓ５３）を実施する。

また、第５の実施形態では、画像取得部３１が被写体画像を取得すると（Ｓ５１）、Ｓ５２の処理と並行して、散乱線情報取得部３６は、予め設定された仮体厚分布を取得し、上記第１の実施形態のオプションの例と同様に上記条件式（１）、（２）に基づいて、取得した仮体厚分布を用いて散乱線情報（一次線像Ｉｐと散乱線像Ｉｓ）を取得し（Ｓ５５）、仮体厚分布を用いて取得した散乱線像Ｉｓ′に基づいて上記第１の実施形態のオプションの例と同様に被写体画像Ｉｋから散乱線像Ｉｓ′を減算することにより散乱線による影響を除去した仮処理後画像を生成して、記憶部４２に記憶する（Ｓ５６）。続いて、画像解析装置３０は、Ｓ５３に示す散乱線を除去する処理により得られる処理後画像が生成されているか否かを判定し、処理後画像が生成されていない場合には（Ｓ５７，Ｎｏ）、処理後画像が生成されるまで、仮処理後画像を表示部４０に表示させる（Ｓ５８）。一方、Ｓ５３に示す処理により得られる処理後画像が生成されている場合には（Ｓ５７,Ｙｅｓ）、画像解析装置は処理後画像を表示部４０に表示させる（Ｓ５４）。

第５の実施形態によれば、被写体画像Ｉｋを取得し、被写体画像Ｉｋの体厚分布Ｔｋを決定する処理を実施した後、散乱線除去処理などの所望の画像処理を実施し、得られた処理後画像を表示する放射線画像解析装置において、被写体画像Ｉｋの体厚分布Ｔｋを決定する一連の処理の処理期間中に、予め設定された仮体厚分布を用いて被写体画像の散乱線除去処理などの所望の画像処理を実施し、仮体厚分布を用いて得られた処理後画像を仮処理後画像として生成し、仮処理後画像を表示部４０に表示させる処理の一部または全部を実施するようにしたため、ユーザが処理後画像の表示までの時間を有効活用して、仮体厚分布を用いた仮処理後画像を観察することで被写体画像の注目すべき部分や、被写体画像に適用された放射線撮影の撮影条件の適否などを大まかに把握できる可能性があり、ユーザの観察を効率化するとともにユーザの観察の参考になる情報を提供することができる。

第５の実施形態において、「予め設定された体厚分布」は、被写体Ｋに概ね対応するような任意の体厚分布を用いることができ、たとえば、過去の複数の被写体画像から体厚決定処理を実施して得られた体厚分布のうち最新の体厚分布とすることができ、標準的な被写体の体厚分布とすることもできる。

また、第５の実施形態における被写体画像Ｉｋの体厚分布を決定する一連の処理（Ｓ５２）は、少なくとも、仮想モデル取得部３２が所定の体厚分布を有する被写体の仮想モデルを取得する処理と、推定画像生成部３３が仮想モデルの放射線撮影により得られる一次線画像を仮想モデルから推定した推定一次線画像と仮想モデルの放射線撮影により得られる散乱線画像を仮想モデルから推定した推定散乱線画像とを合成した画像を、被写体の放射線撮影により得られる被写体画像を推定した画像である推定画像として生成する処理と、修正部３４が推定画像と被写体画像の違いが小さくなるように仮想モデルの体厚分布を修正する処理と、体厚分布決定部３５が修正された仮想モデルの体厚分布を被写体の体厚分布として決定する処理を少なくとも含むもの（本発明に係る体厚分布を決定する処理）であれば、被写体Ｋの体厚分布Ｔｋを決定するいかなる処理としてもよい。例えば、体厚分布を決定する一連の処理（Ｓ５２）を、第２の実施形態に対応する処理を含むようにしてもよい。

また、第５の実施形態における被写体画像Ｉｋから散乱線を除去する処理（Ｓ５３、Ｓ５６）は、ここでは第１の実施形態における被写体画像Ｉｋから散乱線を除去する処理と共通するものとしたが、被写体画像Ｉｋから散乱線を除去可能な任意の方法を適用できる。例えば、第５の実施形態における被写体画像Ｉｋから散乱線除去する処理（Ｓ５３、Ｓ５６）を第３の実施形態の図６のＳ２３〜Ｓ３０に対応する処理としてもよく、または、第４の実施形態の図７のＳ３３〜Ｓ４１に対応する処理としてもよい。また、第５の実施形態において、被写体画像Ｉｋから散乱線除去処理（Ｓ５３、Ｓ５４）を行う例を説明したが、これに限られず、散乱線除去処理に換えて被写体画像Ｉｋに他の画像処理を実施してもよく、散乱線除去処理にあわせてさらに被写体画像Ｉｋに他の画像処理を行ってもよい。また、Ｓ５５〜Ｓ５８の処理は、Ｓ５２の処理と一部または全部が重複するような任意の期間に実施でき、仮処理後画像をより早く表示してユーザの観察の参考になる情報をより早く提供するためには、被写体画像Ｉｋの取得後Ｓ５５〜Ｓ５８の処理を早いタイミングで実施することが好ましい。

また、上記の各実施形態において、図２のＳ０１に示す被写体画像Ｉｋの取得処理は、Ｓ０４に示す被写体画像と推定画像の違いの判定処理に先だって行われるものであれば任意のタイミングで行ってよい。

上記実施形態において、画像解析装置３０を、散乱線情報取得部３６および散乱線除去部３７を省略した構成とし、散乱線情報取得処理および散乱線除去処理を実施しない態様としてもよい。この場合、他の装置に決定した被写体Ｋの体厚分布Ｔｋを出力し、他の装置で体厚分布Ｔｋを用いて被写体画像Ｉｋの画像処理や撮影条件を決定する処理を実施することが考えられる。

なお、本各実施形態に限定されず、本発明によって得られた被写体の体厚分布は、被写体画像に対して、被写体の体厚に応じた画像処理条件を決定するためのあらゆる処理に用いることができる。例えば、本発明によって得られた体厚分布を静止画像または動画像である被写体画像に対する、濃度、コントラストなどの階調処理、ノイズ抑制処理、ダイナミックレンジ調整処理、周波数強調処理などに用いることが考えられる。また、本発明により得られた体厚分布を、被写体画像に対して、体厚に応じた撮影条件を決定するためのあらゆる処理に用いることができる。本発明により得られた体厚分布を用いて各画像処理条件または撮影条件を決定した場合には、被写体画像に対する正確な体厚分布が適用されることにより、決定された画像処理条件または撮影条件による画質改善効果を高めることができる。

また、例えば、管電圧を異ならせて高エネルギーと低エネルギーの放射線をそれぞれ撮影して取得された２つの放射線画像の違いによって放射線画像を取得するエネルギーサブトラクション技術において、本発明によって得られた被写体の体厚分布に応じて、高エネルギー画像から低エネルギー画像を減算する際に体厚の大きい位置における高エネルギー画像の重み付けを大きくするように重み係数を決定する処理を行ってもよい。この場合には、被写体画像に対する正確な体厚分布が適用されることにより、被写体の厚さに応じて放射線の線質が変動するビームハードニング現象の発生の影響を低減して、処理後画像の画質を好適に改善することができる。

また、放射線透視画像の適用技術や、トモシンセシス技術など同一の被写体を複数回撮影して複数の放射線画像を取得する技術分野においては、最初に撮影された被写体画像に対して本発明の手法により被写体の体厚分布を取得し、取得した体厚分布に基づいて後続して撮影される被写体の撮影条件を決定することが好ましい。後続して撮影される被写体画像を、体厚に応じた適切な撮影条件で撮影できるため、後続して撮影される被写体画像を診断に適した画質の画像とすることができる。

また、本発明によって得られた被写体の体厚分布と撮影条件とを対応付けた線量管理情報を、被写体毎にそれぞれ記憶することが好ましい。なお、撮影条件は、線源等に設定される設定値と曝射時間を含むものとすることができ、実際に照射された線量等を検出器等で測定した実績値と曝射時間を含むものとすることもできる。被検者ごとに、該被検者の複数の線量管理情報を線量履歴情報として取得し、この線量履歴情報に基づいて、被検者の所定期間中の各位置における累積被曝量を算出することにより、累積被曝量が所定の臓器などの所定領域に対して許容できない所定の閾値以上となっていないかなど領域毎の線量管理の指標となる有用な情報を提供することができる。また、複数の異なる被検者に対する線量管理情報をそれぞれ取得して、統計解析することにより、体厚分布の傾向に応じてどのような撮影条件が用いられたかを特定することができ、新たな被写体に対する撮影条件の決定または過去の被写体に対する撮影条件の推定のために参考になる情報を提供することができる。

また、トモシンセシス撮影装置によって撮影された複数の被写体画像から、立体視画像を構成する２つの被写体画像を選択する際に、本発明によって得られた被写体の体厚分布を取得し、取得された体厚分布に応じて適切な視差量または輻輳角（立体視画像を構成する２つの被写体画像の各撮影方向がなす角度）となるように２つの被写体画像を選択することが好ましい。例えば、被写体の体厚分布から最大値、平均値、中央値など体厚分布の特徴を示す指標値を抽出し、指標値の範囲ごとに被写体の体厚が大きい程視差量（または輻輳角）が大きくなるように適切な視差量（または輻輳角）を対応付けた視差決定情報を予め作成して取得し、視差決定情報に基づいて、被写体の体厚分布から抽出された指標値に対応する視差量（または輻輳角）を、立体視画像の視差量（または輻輳角）として決定し、かかる視差量（または輻輳角）となるような２つの被写体画像を立体視画像を構成する被写体画像として選択することができる。一例として、被写体の体厚分布に応じて、体厚分布の特徴を示す指標値を算出して、該指標値により、被写体を、細め型、標準型、太め型など複数の体格的な分類に区分し、区分された分類に対応する視差量（または輻輳角）を被写体のステレオ画像の視差量として決定し、かかる視差量となるような２つの被写体画像を立体視画像を構成する被写体画像として選択することができる。このように、複数の被写体画像から、立体視画像を構成する２つの被写体画像を選択する際に、本発明によって得られた被写体の体厚分布を取得し、取得された体厚分布に応じて体厚が大きい程視差量または輻輳角が大きくなるように２つの被写体画像を適切に選択した場合には、体厚分布に応じて、観察に適した品質の立体視画像を生成することができる。なお、視差量、輻輳角については、本出願の過去出願（特開2013-198736号公報、特開2013-198508号公報、特開2013-154165号公報等）を参考にすることができる。

なお、収束型グリッドを用いて被写体画像が撮影されている場合に、放射線の斜入によって濃度むらが生じる恐れがある。このような濃度むらの発生に備えて、本発明により得られた体厚分布を取得し、人体を正面から撮影した画像などのように被写体画像の被写体が左右対称の体厚分布を有することが既知である場合には、体厚分布が実質的に左右対称であるか否かを判定し、左右対称でない場合には、表示や音声により再撮影を促すようにしてもよい。斜入による濃度むらによって観察に適した品質の処理後画像が得られないことを防ぐことができる。

上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。本発明の態様は、上述した個々の実施例（第１〜第５実施形態、その他の変形例および応用例）に限定されるものではなく、個々の実施例の各要素のいかなる組合せも本発明に含み、また、当業者が想到しうる種々の変形も含むものである。すなわち、特許請求の範囲に規定された内容およびその均等物から導き出される本発明の概念的な思想と趣旨を逸脱しない範囲で種々の追加、変更および部分的削除が可能である。

また、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。たとえば、画像解析装置の構成要素の一部または全部は、1台のワークステーションにより構成されたものであってもよく、ネットワークを介して接続された一台以上のワークステーション、サーバ、記憶装置によって構成されたものであってもよい。

また、上記実施形態においては、放射線検出器１４を用いて被写体の放射線画像を撮影する撮影装置１０において取得した放射線画像を用いて散乱線除去処理を行っているが、特開平８−２６６５２９号公報、特開平９−２４０３９号公報等に示される放射線検出体としての蓄積性蛍光体シートに被写体の放射線画像情報を蓄積記録し、蓄積性蛍光体シートから光電的に読み取ることにより取得した放射線画像を用いた場合においても、本発明を適用できることはもちろんである。

１０ 撮影装置
２０ 制御装置
３０ 画像解析装置（放射線画像解析装置）
３１ 画像取得部
３２ 仮想モデル取得部
３３ 推定画像生成部
３４ 修正部
３５ 体厚分布決定部
３６ 散乱線情報取得部
３７ 散乱線除去部
３８ 入力部
３９ 特性取得部
４０ 表示部
４２ 記憶部
Ｋ 被写体
Ｉｋ 被写体画像
Ｉｍ モデル画像
Ｍ 仮想モデル
Ｔｋ 被写体画像の体厚分布
Ｔｎ モデル画像の体厚分布
Ｔｐ 一次線透過率
Ｔｓ 散乱線透過率

Claims (16)

1. 被写体の放射線撮影により得られた被写体画像を解析して、前記被写体の体厚分布を推定する放射線画像解析装置であって、
   前記被写体画像を取得する画像取得部と、
   所定の体厚分布を有する前記被写体の仮想モデルを取得する仮想モデル取得部と、
   前記仮想モデルの放射線撮影により得られる一次線画像を前記仮想モデルから推定した推定一次線画像と前記仮想モデルの放射線撮影により得られる散乱線画像を前記仮想モデルから推定した推定散乱線画像とを合成した画像を、前記被写体の放射線撮影により得られる放射線画像を推定した画像である推定画像として生成する推定画像生成部と、
   前記推定画像と前記被写体画像の各対応する位置の画素の画素値の違いを表すエラー値が小さくなるように前記仮想モデルの前記体厚分布を修正する修正部と、
   前記修正された前記仮想モデルの前記体厚分布を前記被写体の体厚分布として決定する体厚分布決定部とを備えたことを特徴とする放射線画像解析装置。
2. 前記エラー値が、前記推定画像と前記被写体画像の各対応する位置の画素の画素値の違いの総和を表すものである請求項１記載の放射線画像解析装置。
3. 前記体厚分布決定部は、前記エラー値が閾値以下であるかを判定し、前記エラー値が前記閾値以下である場合に、前記仮想モデルの前記体厚分布を前記被写体の体厚分布として決定する請求項１又は２に記載の放射線画像解析装置。
4. 前記仮想モデル取得部が、前記修正された体厚分布を有する前記仮想モデルをさらに取得し、
   前記推定画像生成部が、前記修正された体厚分布を有する前記仮想モデルから前記推定画像をさらに生成し、
   前記修正部が、さらに生成された前記推定画像と前記被写体画像の違いが小さくなるように前記仮想モデルの前記体厚分布をさらに修正することを特徴とする請求項１から３のいずれか１項記載の放射線画像解析装置。
5. 前記修正部は、前記推定画像と前記被写体画像の差分画像の画素値の絶対値の総和または前記差分画像の画素値の二乗和が小さくなるように、前記仮想モデルの前記体厚分布を修正することを特徴とする請求項１から４のいずれか１記載の放射線画像解析装置。
6. 前記修正部は、前記仮想モデルの一画素以上の部分領域ごとに、前記仮想モデルの前記体厚分布を変動させて、前記推定画像と前記被写体画像との違いを小さくする前記部分領域の体厚を算出し、算出された各前記部分領域の体厚によって前記仮想モデルの前記体厚分布を修正することを特徴とする請求項１から５のいずれか１項記載の放射線画像解析装置。
7. 前記エラー値が前記推定画像と前記被写体画像の違いを表すエラー関数の出力値であり、
   前記修正部が、前記仮想モデルの一画素以上の部分領域ごとに、前記仮想モデルの当該部分領域の体厚を変動させて、前記体厚分布に関する前記エラー関数の一次偏微分を算出し、該エラー関数の一次偏微分に基づいて前記エラー関数の出力値を最小化し、該エラー関数の出力値が最小化された際の該部分領域の体厚を該部分領域の体厚の修正値とする請求項６記載の放射線画像解析装置。
8. 前記所定の体厚分布が、前記被写体とは異なる比較用被写体を放射線撮影して得た比較用被写体画像と、前記比較用被写体を３次元撮影して得た３次元画像を取得し、取得された該３次元画像の各位置において、前記比較用被写体画像の放射線経路に対応する直線上の前記比較用被写体の体厚を計測することにより作成されていることを特徴とする請求項１から７のいずれか１項に記載の放射線画像解析装置。
9. 前記仮想モデルが、該仮想モデルに含まれる構造物と該構造物の配置と該構造物の放射線に対する特性の少なくとも１つを表す特性情報をさらに有し、
   前記推定画像生成部が、前記特性情報に基づいて前記仮想モデルの各位置に対応する構造物に応じて前記推定画像を算出するためのパラメータを選択して前記推定画像を生成することを特徴とする請求項１から８のいずれか１項記載の放射線画像解析装置。
10. 前記推定画像生成部が、前記被写体画像に含まれる構造物と構造物の配置と該構造物の放射線に対する特性とを表す特性情報を前記仮想モデルの前記特性情報として取得し、該特性情報に基づいて、前記仮想モデルの各位置に対応する構造物に応じて前記推定画像を算出するためのパラメータを選択して前記推定画像を生成することを特徴とする請求項９に記載の放射線画像解析装置。
11. 決定された前記被写体の前記体厚分布を用いて前記被写体画像の散乱線を推定した散乱線情報を取得する散乱線情報取得部と、
    取得した前記散乱線情報に基づいて、前記被写体画像の散乱線の除去処理を行う散乱線除去部とをさらに備えたことを特徴とする請求項１から１０のいずれか１項に記載の放射線画像解析装置。
12. 前記被写体画像が散乱線除去グリッドを用いないで撮影された画像であることを特徴とする請求項１から１１のいずれか１項記載の放射線画像解析装置。
13. 前記被写体とは異なる比較用被写体を放射線撮影して得た比較用被写体画像を複数取得し、前記比較用被写体を測定して各々の前記比較用被写体の体厚分布を取得し、前記比較用被写体ごとに、該比較用被写体の体厚分布と、該比較用被写体の体格を表す第１の体格情報とを互いに対応付けし、前記被写体の体格を表す第２の体格情報を取得し、前記第２の体格情報と類似する前記第１の体格情報を有する前記比較用被写体の体厚分布を特定し、前記特定された比較用被写体の体厚分布を前記所定の体厚分布とする特定部をさらに備える、請求項１から１２のいずれか１項記載の放射線画像解析装置。
14. 前記第１の体格情報は、前記比較用被写体画像の濃度ヒストグラムであり、
    前記第２の体格情報は、前記被写体画像の濃度ヒストグラムである、請求項１３に記載の放射線画像解析装置。
15. 放射線画像解析装置に実行される、被写体に放射線を照射することにより撮影された被写体画像を解析して、前記被写体の体厚分布を推定する放射線画像解析方法であって、
    前記被写体画像を取得する画像取得ステップと、
    所定の体厚分布を有する前記被写体の仮想モデルを取得する仮想モデル取得ステップと、
    前記仮想モデルの放射線撮影により得られる一次線画像を前記仮想モデルから推定した推定一次線画像と前記仮想モデルの放射線撮影により得られる散乱線画像を前記仮想モデルから推定した推定散乱線画像とを合成した画像を、前記被写体の放射線撮影により得られる放射線画像を推定した画像である推定画像として生成する推定画像生成ステップと、
    前記推定画像と前記被写体画像の各対応する位置の画素の画素値の違いを表すエラー値が小さくなるように前記仮想モデルの前記体厚分布を修正する修正ステップと、
    前記修正された前記仮想モデルの前記体厚分布を前記被写体の体厚分布として決定する体厚分布決定ステップとを有することを特徴とする放射線画像解析方法。
16. 被写体に放射線を照射することにより撮影された被写体画像を解析して、前記被写体の体厚分布を推定させる放射線画像解析プログラムであって、
    コンピュータに、
    前記被写体画像を取得する画像取得ステップと、
    所定の体厚分布を有する前記被写体の仮想モデルを取得する仮想モデル取得ステップと、
    前記仮想モデルの放射線撮影により得られる一次線画像を前記仮想モデルから推定した推定一次線画像と前記仮想モデルの放射線撮影により得られる散乱線画像を前記仮想モデルから推定した推定散乱線画像とを合成した画像を、前記被写体の放射線撮影により得られる放射線画像を推定した画像である推定画像として生成する推定画像生成ステップと、
    前記推定画像と前記被写体画像の各対応する位置の画素の画素値の違いを表すエラー値が小さくなるように前記仮想モデルの前記体厚分布を修正する修正ステップと、
    前記修正された前記仮想モデルの前記体厚分布を前記被写体の体厚分布として決定する体厚分布決定ステップとを実行させることを特徴とする放射線画像解析プログラム。