JP2015181649A - 放射線画像処理装置および方法並びにプログラム - Google Patents

Abstract

【課題】精度よくかつ効率的に、放射線画像の散乱線に起因するコントラストの低下を改善できるようにする。【解決手段】周波数分解部２２が、第１の放射線画像を周波数分解して、複数の帯域画像を取得する。散乱線除去部２３が、放射線量と線形の関係にある第２の放射線画像を用いて、第２の放射線画像に含まれる被写体による散乱線成分を算出し、散乱線成分を第２の放射線画像から除去して散乱線除去済み放射線画像を取得する。変換関数決定部２４が、第２の放射線画像に含まれる散乱線成分に応じて、少なくとも１つの帯域画像を変換する変換関数を決定する。再構成部２６が、帯域画像を変換関数により変換して変換済み帯域画像を生成する。散乱線除去済み放射線画像および変換済み帯域画像を再構成して処理済み放射線画像を取得する。【選択図】図２

Description

本発明は、放射線画像に対して散乱線除去処理を含む画像処理を施す放射線画像処理装置および方法並びに放射線画像処理方法をコンピュータに実行させるためのプログラムに関する。

従来、被写体を透過した放射線により被写体の放射線画像を撮影する際、とくに被写体の厚さが大きいと、被写体内において放射線が散乱し、この散乱放射線（以下散乱線とする）により、取得される放射線画像のコントラストが低下するという問題がある。このため、放射線画像の撮影時には、放射線を検出して放射線画像を取得するための放射線検出器に散乱線が照射されないように、被写体と放射線検出器との間に散乱線除去グリッド（以下単にグリッドとする）を配置して撮影を行うことがある。グリッドを用いて撮影を行うと被写体により散乱された放射線が放射線検出器に照射されにくくなるため、放射線画像のコントラストを向上させることができる。

一方、グリッドを用いた撮影を行うと、被写体像とともにグリッドに対応した細かな縞模様（モアレ）が放射線画像に含まれてしまうため、画像が見難いものとなってしまう。

このため、グリッドを使用することなく撮影を行うことにより取得した放射線画像に対して、散乱線の除去による画質改善の効果を、撮影条件に基づいて画像処理により付与する手法が提案されている（特許文献１，２および非特許文献１参照）。特許文献１および非特許文献１の手法は、放射線画像を複数の周波数成分に周波数分解し、散乱線の成分と見なせる低周波成分に対してコントラストを抑制する処理を行い、散乱線成分の除去によるコントラスト低下を改善するために、高周波成分に対してコントラストを強調する処理を行い、処理後の周波数成分を合成することにより、散乱線の成分が除去された放射線画像を取得するものである。

なお、上記特許文献１に記載された手法は、低周波成分の画素値に応じたゲインを低周波成分に乗算することにより、散乱線除去処理を行っている。ここで、ゲインは１未満の値であり、低周波帯域ほど、また明るい画素値ほど小さな値となっている。一方、高周波成分に対しては、コントラストの低下を補うために、１よりも大きい値のゲインを乗算することにより、放射線画像のコントラストを強調している。また、放射線画像に含まれる高コントラストの物体（骨、金属等）の過強調を抑えるために、高コントラストとなる高周波成分のゲインを非線形に抑制することも行っている。

特許文献１，２および非特許文献１の手法によれば、撮影時にグリッドが不要となるため、撮影時の患者の負担を軽減することができ、かつ濃度ムラおよびモアレによる画質の低下を防止しつつ、散乱線を除去してコントラストの改善された放射線画像を得ることができる。

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

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

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

ところで、放射線画像は、被写体内部の組成（物質およびその密度）によって異なる放射線の減弱を画像の濃淡により表したものであり、放射線画像の画素値は、被写体を透過した放射線量に相当する。ここで、放射線画像の各画素の画素値Ｉｏｕｔは、被写体に入射した放射線量をＩｉｎ、放射線の組成に応じた線減弱係数μおよび厚さｔを用いて、Ｉｏｕｔ＝Ｉｉｎ・ｅ^−μｔにより表される。一方、放射線検出器においては、各画素位置において取得した信号を対数変換している。このため、放射線画像の各画素の画素値は対数放射線量と線形の関係、すなわちｌｏｇ（Ｉｏｕｔ）＝ｌｏｇ（Ｉｉｎ）＋（−μｔ）により表すことができる。そして、対数放射線量と線形の関係にあるものとなるように対数変換した放射線画像をデジタルの信号に変換し、さらに階調処理、周波数処理およびダイナミックレンジ圧縮処理等の画質を向上させるための画像処理を行って、高画質の放射線画像を取得するようにしている。

このように取得された放射線画像の各画素の画素値は、対数変換された放射線量に対して線形の関係にあるため、撮影した被写体の厚さおよび組成の相違をそのまま画像に反映させることができる。これにより、画素値に対する線形の演算により画素値を変更することができることから、階調処理、周波数処理およびダイナミックレンジ圧縮処理等の画像処理を、放射線画像によらず、同じような強調度により再現することが可能となる。例えば、上記特許文献１および非特許文献１の手法においては、低周波成分に対して１未満のゲインを乗算することにより散乱線成分を抑制し、高周波成分に対しては１より大きいゲインを乗算することによりコントラストを強調しているが、これは放射線画像が対数放射線量と線形の関係にあるために行うことができるものである。したがって、対数変換された放射線量に対して線形の関係にある放射線画像に対して画像処理を行うことにより、撮影条件によらず、診断に適した高画質の放射線画像を取得することができる。

米国特許第８０６４６７６号明細書 特開平２０１３−２４０６９６号公報 特開平０２−２４４８８１号公報

C Fivez et al, Multi-resolution contrast amplification in digital radiography with compensation for scattered radiation, 1996 IEEE, pp339-342. Trotter, 他4名, "Thickness-dependent Scatter Correction Algorithm for Digital Mammography", Proc. SPIE Vol.4682, 2002年5月, p.469-478

一方、散乱線は、放射線が被写体内において物質を透過する際に、その一部が方向を変えることにより発生するものであり、上記特許文献３および非特許文献２等の手法を用いて、放射線画像のある画素位置における散乱線成分を推定するためには、周辺の画素の画素値からの散乱線成分の総和を算出する必要がある。ここで、散乱線成分は放射線の線量に相当するため、散乱線成分を算出するためには、画素値が放射線量に対して線形の関係にある放射線画像を用いる必要がある。

このため、放射線量と線形の関係にある放射線画像を用いて散乱線成分を算出し、かつ散乱線除去処理およびコントラスト強調処理を行うことが考えられる。このような場合、放射線画像の低周波成分から散乱線成分を減算し、高周波成分にはゲインを加算することにより、散乱線を除去しつつコントラストを強調できる。

しかしながら、放射線量と線形の関係にある放射線画像においては、線形の演算により画素値を変更することが困難であるため、ゲインを一意に決定することが難しい。とくに、過強調を抑えるような非線形のゲインは、放射線量と線形の関係にある放射線画像においては決定することが困難である。このため、放射線量と線形の関係のある放射線画像を用いたのでは、精度よくコントラストを強調することができない。逆に、対数放射線量と線形の関係にある放射線画像を用いて散乱線成分を算出することが考えられるが、上述したように散乱線は放射線の一部であることから、対数放射線量と線形の関係にある放射線画像を用いて散乱線成分を算出することは困難である。

本発明は上記事情に鑑みなされたものであり、精度よくかつ効率的に、放射線画像の散乱線に起因するコントラストの低下を改善できるようにすることを目的とする。

本発明による放射線画像処理装置は、被写体に放射線を照射することにより撮影された、第１の放射線画像を周波数分解して、第１の放射線画像の複数の周波数成分を表す複数の帯域画像を取得する周波数分解手段と、
放射線の放射線量と線形の関係にある被写体の第２の放射線画像を用いて、第２の放射線画像に含まれる被写体による散乱線成分を算出し、散乱線成分を第２の放射線画像から除去して散乱線除去済み放射線画像を取得する散乱線除去手段と、
第２の放射線画像に含まれる散乱線成分に応じて、少なくとも１つの帯域画像を変換する変換関数を決定する変換関数決定手段と、
少なくとも１つの帯域画像を変換関数により変換して変換済み帯域画像を生成する変換手段と、
散乱線除去済み放射線画像および変換済み帯域画像を再構成して処理済み放射線画像を取得する再構成手段とを備えたことを特徴とするものである。

なお、本発明による放射線画像処理装置においては、変換関数決定手段を、散乱線成分と放射線に含まれる一次線成分および散乱線成分の和との比率に基づいて、変換関数を決定する手段としてもよい。

また、本発明による放射線画像処理装置においては、散乱線除去手段を、被写体の体厚分布を用いて、散乱線成分および一次線成分を算出する手段とし、
変換関数決定手段を、算出された散乱線成分および一次線成分から比率を算出する手段としてもよい。

また、本発明による放射線画像処理装置においては、第２の放射線画像に基づいて、体厚分布を推定する体厚推定手段をさらに備えるものとしてもよい。

また、本発明による放射線画像処理装置においては、第１の放射線画像を、画素値が対数放射線量と線形の関係にあるものとしたもよい。

また、本発明による放射線画像処理装置においては、第２の放射線画像を、第１の放射線画像を逆対数変換した画像としてもよい。

また、本発明による放射線画像処理装置においては、再構成手段を、散乱線除去済み放射線画像を対数変換し、対数変換された散乱線除去済み放射線画像および変換済み帯域画像を再構成して処理済み放射線画像を取得する手段としてもよい。

また、本発明による放射線画像処理装置においては、被写体の撮影時に、散乱線を除去するために使用が想定される、仮想的なグリッドの特性である仮想グリッド特性を取得する特性取得手段をさらに備えるものとし、
散乱線除去手段を、仮想グリッド特性に基づいて、散乱線成分を算出する手段とし、
変換関数決定手段を、算出された散乱線成分に応じて変換関数を決定する手段としてもよい。

また、本発明による放射線画像処理装置においては、変換関数決定手段を、帯域画像の画素値が所定のしきい値を越えると、画素値を抑制するように変換関数を決定する手段としてもよい。

本発明による放射線画像処理方法は、被写体に放射線を照射することにより撮影された、第１の放射線画像を周波数分解して、第１の放射線画像の複数の周波数成分を表す複数の帯域画像を取得し、
放射線の放射線量と線形の関係にある被写体の第２の放射線画像を用いて、第２の放射線画像に含まれる被写体による散乱線成分を算出し、
散乱線成分を第２の放射線画像から除去して散乱線除去済み放射線画像を取得し、
第２の放射線画像に含まれる散乱線成分に応じて、少なくとも１つの帯域画像を変換する変換関数を決定し、
少なくとも１つの帯域画像を変換関数により変換して変換済み帯域画像を生成し、
散乱線除去済み放射線画像および変換済み帯域画像を再構成して処理済み放射線画像を取得することを特徴とするものである。

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

本発明によれば、第１の放射線画像が周波数分解されて、第１の放射線画像の複数の周波数成分を表す複数の帯域画像が取得される。また、放射線の放射線量と線形の関係にある第２の放射線画像を用いて、第２の放射線画像に含まれる被写体による散乱線成分が算出され、散乱線成分が第２の放射線画像から除去されて散乱線除去済み放射線画像が取得される。そして、第２の放射線画像に含まれる散乱線成分に応じて、少なくとも１つの帯域画像を変換する変換関数が決定され、少なくとも１つの帯域画像が変換関数により変換されて変換済み帯域画像が生成され、散乱線除去済み放射線画像および変換済み帯域画像が再構成されて処理済み放射線画像が取得される。

ここで、第２の放射線画像は放射線量と線形の関係にあるため、精度よく散乱線成分を算出することができ、その結果、第２の放射線画像から精度よく散乱線成分を除去することができる。一方、変換関数は、放射線量と線形の関係にある第２の放射線画像の散乱線成分から算出しているため、変換関数により帯域画像を精度よく強調できる。したがって、再構成された放射線画像においては、散乱線が除去され、かつ放射線画像に含まれる高周波成分が強調されるため、散乱線に起因するコントラストの低下が改善されたものとなる。

とくに、第１の放射線画像を、対数放射線量と線形の関係にあるものとすることにより、帯域画像の画素値は対数放射線量と線形の関係にあるものとなるため、散乱線成分から変換関数を容易に決定することができ、かつ変換関数により帯域画像を変換することにより、効率よく帯域画像の周波数成分を強調できる。したがって、放射線画像の散乱線に起因するコントラストの低下を、精度よく効率的に改善することができる。

本発明の第１の実施形態による放射線画像処理装置を適用した放射線画像撮影システムの構成を示す概略ブロック図 第１の実施形態における放射線画像撮影システムのコンピュータ内部の概略構成を示すブロック図 第１の実施形態において行われる処理を示すフローチャート 周波数分解を説明するための図 第１の実施形態における散乱線除去処理のフローチャート 放射線量とその対数放射線量との関係を示すグラフ 画素値を抑制するよう帯域画像を変換する変換関数を示す図 再構成を説明するための図 第２の実施形態における放射線画像撮影システムのコンピュータ内部の概略構成を示すブロック図 体厚推定部の構成を示す概略ブロック図 体厚推定処理のフローチャート 体厚分布の対応付けテーブルの例を示す図 推定画像の生成方法の例を説明するための図 推定画像の生成方法の他の例を説明するための図 第３の実施形態における散乱線除去処理のフローチャート

以下、図面を参照して本発明の実施形態について説明する。図１は本発明の第１の実施形態による放射線画像処理装置を適用した放射線画像撮影システムの構成を示す概略ブロック図である。図１に示すように、第１の実施形態による放射線画像撮影システムは、被写体の放射線画像に対して、散乱線除去処理を含む各種画像処理を行うためのものであり、図１に示すように、撮影装置１と、本実施形態による放射線画像処理装置を内包するコンピュータ２とを備える。

撮影装置１は被写体ＭにＸ線を照射するＸ線源３と、被写体Ｍを透過したＸ線を検出して被写体Ｍの放射線画像を取得する放射線検出器５とを備える。

放射線検出器５は、放射線画像の記録と読み出しを繰り返して行うことができるものであり、放射線の照射を直接受けて電荷を発生する、いわゆる直接型の放射線検出器を用いてもよいし、放射線を一旦可視光に変換し、その可視光を電荷信号に変換する、いわゆる間接型の放射線検出器を用いるようにしてもよい。また、放射線画像信号の読出方式としては、ＴＦＴ（thin film transistor）スイッチをオン・オフさせることによって放射線画像信号が読み出される、いわゆるＴＦＴ読出方式のものや、読取り光を照射することによって放射線画像信号が読み出される、いわゆる光読出方式のものを用いることが望ましいが、これに限らずその他のものを用いるようにしてもよい。なお、放射線検出器５においては、内部で電荷信号に対して対数変換およびＡ／Ｄ変換がなされるため、放射線検出器５から出力される放射線画像は、各画素の画素値が対数放射線量と線形の関係にあるものとなっている。

コンピュータ２は、中央処理装置（ＣＰＵ）、半導体メモリ、通信インターフェースおよびハードディスクやＳＳＤ等のストレージデバイス等を備えており、これらのハードウェアによって、図２に示すように本発明の放射線画像処理装置を構成する、制御部２１、周波数分解部２２、散乱線除去部２３、変換関数決定部２４、変換部２５、再構成部２６および記憶部２７が構成されている。

制御部２１は、Ｘ線源３および放射線検出器５に対して撮影の制御を行ったり、放射線検出器５から放射線画像の読み出しを行ったり、コンピュータ２において行われる処理全体の制御を行ったりするものである。

周波数分解部２２は、放射線検出器５から取得した放射線画像（以下第１の放射線画像とする）を周波数分解して、それぞれが異なる複数の周波数帯域の周波数成分を表す複数の帯域画像を作成する。

散乱線除去部２３は、第１の放射線画像を逆対数変換し、これにより取得された放射線画像（以下、第２の放射線画像とする）を用いて、第２の放射線画像に含まれる被写体Ｍによる散乱線成分を算出し、散乱線成分を用いた散乱線除去処理を行って散乱線除去済み放射線画像を取得する。

変換関数決定部２４は、第２の放射線画像に含まれる散乱線成分に応じて、帯域画像を変換するための変換関数を決定する。

変換部２５は、帯域画像を変換関数により変換して変換済み帯域画像を生成する。

再構成部２６は、散乱線除去済み放射線画像および変換済み帯域画像を再構成して処理済み放射線画像を取得する。

記憶部２７は、コンピュータ２が行う散乱線除去処理に必要な各種情報を記憶する。

表示部６は、ＣＲＴあるいは液晶ディスプレイ等からなり、撮影により取得された放射線画像および後述する画像処理に必要な各種入力の補助を行う。入力部８は、キーボード、マウスあるいはタッチパネル等からなる。

なお、制御部２１、周波数分解部２２、散乱線除去部２３、変換関数決定部２４、変換部２５、および再構成部２６が行う処理は、記憶部２７に記憶されたコンピュータプログラムにより中央処理装置が行う。なお、各部のそれぞれの処理を行う複数の処理装置をコンピュータ２に設けるようにしてもよい。

次いで、第１の実施形態において行われる処理について説明する。図３は第１の実施形態において行われる処理を示すフローチャートである。被写体Ｍの撮影が行われてコンピュータ２の制御部２１が、画素値が対数放射線量と線形の関係にある第１の放射線画像Ｇ０を取得すると（ステップＳＴ１）、周波数分解部２２が、第１の放射線画像Ｇ０を周波数分解し、第１の放射線画像Ｇ０についての複数の周波数成分を表す複数の帯域画像を生成する（ステップＳＴ２）。図４は周波数分解を説明するための図である。まず、周波数分解部２２は、放射線画像Ｇ０に対してσ＝１のガウシアンフィルタによりフィルタリング処理を行って、放射線画像Ｇ０を１／２に縮小してガウシアン成分Ｇ１を生成する。ガウシアン成分Ｇ１は放射線画像Ｇ０を１／２に縮小したものとなる。次いで、周波数分解部２２は、例えば３次Ｂスプライン補間あるいは線形補間等の周知の補間演算を行って、ガウシアン成分Ｇ１を放射線画像Ｇ０と同一サイズに拡大し、拡大したガウシアン成分Ｇ１を放射線画像Ｇ０から減算して、最高周波数帯域の周波数成分を表す帯域画像であるラプラシアン成分Ｌ０を生成する。なお、本実施形態では、最高周波数帯域を便宜上第０の周波数帯域と称する。

次いで、周波数分解部２２は、ガウシアン成分Ｇ１に対してσ＝１のガウシアンフィルタによりフィルタリング処理を行って、ガウシアン成分Ｇ１をさらに１／２に縮小してガウシアン成分Ｇ２を生成し、ガウシアン成分Ｇ２をガウシアン成分Ｇ１と同一サイズに拡大し、拡大したガウシアン成分Ｇ２をガウシアン成分Ｇ１から減算して、第１の周波数帯域の周波数成分を表す帯域画像であるラプラシアン成分Ｌ１を生成する。さらに、所望とする周波数帯域の帯域画像であるラプラシアン成分が生成されるまで上記の処理を繰り返すことにより、複数の周波数帯域のラプラシアン成分Ｌｊ（ｊ＝０〜ｎ）を生成する。本実施形態においては、第３の周波数帯域のラプラシアン成分Ｌ３および第４の周波数帯域のガウシアン成分Ｇ４が得られるまで、上記の処理を繰り返す。

ここで、ガウシアン成分の各画素の信号値はその画素の濃度を表し、ラプラシアン成分の各画素の信号値は、その画素におけるその周波数帯域の周波数成分の大きさを表すものとなる。このため、以降の説明においては、第４の周波数帯域のガウシアン成分について、放射線画像Ｇ４と称するものとする。なお、ウェーブレット変換等の他の周波数分解の手法を用いることにより、周波数帯域が異なる複数の帯域画像を生成してもよい。また、元の放射線画像Ｇ０と同一サイズのガウシアン成分およびラプラシアン成分が取得されるように周波数分解を行ってもよい。

次いで、散乱線除去部２３が、第４の周波数帯域のガウシアン成分すなわち放射線画像Ｇ４を用いて、散乱線除去処理を行う（ステップＳＴ３）。図５は散乱線除去処理のフローチャートである。

本実施形態においては、周波数分解部２２は、画素値が対数放射線量と線形の関係にある放射線画像Ｇ０を周波数分解しているため、放射線画像Ｇ４の画素値は対数放射線量と線形の関係にある。散乱線除去部２３は、放射線画像Ｇ４を逆対数変換して、各画素の画素値が放射線量と線形の関係にある放射線画像ＧＬ４（第２の放射線画像）を生成する（ステップＳＴ１１）。ここで、放射線画像Ｇ４の画素値Ｉｏｕｔは、入射線量Ｉｉｎ、放射線の組成に応じた線減弱係数μおよび厚さｔを用いて、ｌｏｇ（Ｉｏｕｔ）＝ｌｏｇ（Ｉｉｎ）＋（−μｔ）の関係にあるが、逆対数変換を行うことにより、Ｉｏｕｔ＝Ｉｉｎ・ｅ^−μｔの関係を有するものとなる。

次いで、散乱線除去部２３は、放射線画像ＧＬ４を解析することにより、放射線画像ＧＬ４の一次線成分および散乱線成分を取得する（ステップＳＴ１２）。放射線画像ＧＬ４の解析は、例えば、放射線画像Ｇ０の撮影時における照射野情報、被写体情報および撮影条件等に基づいて行う。

照射野情報とは、照射野絞りを用いて撮影を行った場合における、放射線画像ＧＬ４に含まれる照射野の位置および大きさに関する照射野分布を表す情報である。被写体情報とは、上述した胸部、腹部および頭部等の被写体Ｍの種類に加えて、被写体Ｍの放射線画像ＧＬ４上での位置、被写体Ｍの組成の分布、被写体Ｍの大きさおよび被写体Ｍの厚さ等に関する情報である。撮影条件とは、管電圧、ｍＡｓ値、撮影距離（Ｘ線源３から被写体Ｍまでの距離と被写体Ｍから放射線検出器５までの距離との合計）に加えて、エアギャップ量（被写体Ｍから放射線検出器５までの距離）、および放射線検出器５の特性等に関する情報を含む。これらの照射野情報、被写体情報および撮影条件は、放射線画像ＧＬ４に含まれる散乱線の分布を決める要因となっている。例えば、散乱線の大小は照射野の大きさにより左右され、被写体Ｍの厚さが大きいほど散乱線は多くなり、被写体Ｍと放射線検出器５との間に空気が存在すると散乱線が減少する。したがって、これらの情報を用いることにより、より正確に一次線成分および散乱線成分を取得することができる。

散乱線除去部２３は、放射線画像ＧＬ４内の被写体厚の分布Ｔ（ｘ，ｙ）から、下記の式（１）、（２）にしたがって一次線成分および散乱線成分を算出する。

Icp(x,y) = Io(x,y)×exp(-μ×T(x,y)) …（１）
Ics(x,y) = Io(x,y)＊Sσ(T(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σ（Ｔ（ｘ，ｙ））をＴ（ｘ，ｙ）にて近似してもよい。

ここで、放射線画像ＧＬ４の各画素の画素値は、一次線成分と散乱線成分とを加算した値となる。その一方で、散乱線除去の演算は、画素値が対数放射線量と線形の関係にある放射線画像Ｇ４を用いて行うことが好ましい。このため、散乱線除去部２３は、上記式（１）および（２）に基づいて算出した一次線成分および散乱線成分から、下記式（３）に基づいて、放射線画像ＧＬ４の各画素における散乱線を除去するための変換関数Ｓ（ｘ，ｙ）を算出する。なお、変換関数Ｓ（ｘ，ｙ）は０〜１の間の値をとる。ここで、変換関数Ｓ（ｘ，ｙ）は、放射線画像ＧＬ４における散乱線の含有率分布を表すものとなる。

S(x,y) = Icp(x,y)/(Ics(x,y)+Icp(x,y)) …（３）
散乱線除去部２３は、下記の式（４）により、変換関数Ｓ（ｘ，ｙ）を，逆対数変換前の放射線画像Ｇ４の各画素の画素値Ｇ４（ｘ，ｙ）に乗算することにより、散乱線除去済みの放射線画像Ｇ４′を算出し（ステップＳＴ１３）、散乱線除去処理を終了する。

G4'(x,y)=G4(x,y)×S(x,y) …（４）
変換関数Ｓ（ｘ，ｙ）は０〜１の間の値となるため、放射線画像Ｇ４（ｘ，ｙ）に対して変換関数Ｓ（ｘ，ｙ）を乗算することにより、放射線画像Ｇ４の画素値が低減され、その結果、散乱線成分が除去されることとなる。なお、この変換関数Ｓ（ｘ，ｙ）は、変換関数決定部２４において算出してもよい。

なお、散乱線除去部２３は、下記の式（５）に示すように、放射線画像ＧＬ４の画素値から散乱線成分を減算することにより、散乱線除去済みの放射線画像ＧＬ４′を算出するようにしてもよい。この場合、散乱線除去部２３は、放射線画像ＧＬ４′を対数変換して、散乱線除去済みの放射線画像Ｇ４′を取得するものとなる。

GL4'(x,y)=GL4(x,y)-Ics(x,y) …（５）
図２に戻り、変換関数決定部２４は、ラプラシアン成分Ｌ０〜Ｌ３を変換するための変換関数を決定する（ステップＳＴ４）。以下、変換関数の決定について説明する。本実施形態においては、散乱線除去部２３において、画素値が対数線量と線形の関係にある放射線画像Ｇ４を逆対数変換して、画素値が放射線量と線形の関係にある放射線画像ＧＬ４を生成し、放射線画像ＧＬ４から一次線成分Ｉｃｐおよび散乱線成分Ｉｃｓを算出し、これを用いて散乱線除去処理を行っている。図６は放射線量とその対数放射線量との関係を示すグラフである。図６において横軸が放射線量と線形の関係にある放射線画像の画素値を、縦軸が対数放射線量と線形の関係にある放射線画像の画素値を示す。式（４）による散乱線除去処理は、式（５）に示すように、放射線量と線形の関係にある放射線画像ＧＬ４の各画素の画素値から散乱線成分Ｉｃｓを減算することにより散乱線除去処理を行っていることと同義である。ここで、ある画素の画素値をＸ１、散乱線成分Ｉｃｓを減算した後の画素値をＸ２とすると、図６に示すように対数放射線量と線形の関係にある放射線画像においては、画素値はｌｏｇＸ１からｌｏｇＸ２に小さくなる。

図６に示す関係において、画素値ｌｏｇＸ１，ｌｏｇＸ２における対数軸上での傾きは、それぞれ１／Ｘ１，１／Ｘ２となり、Ｘ１＞Ｘ２であるため、散乱線成分を減算することにより、対数軸上における画素値の変化量は大きくなる。変換関数決定部２４は、Ｘ１＝Ｉｃｐ（ｘ，ｙ）＋Ｉｓｃ（ｘ，ｙ）、Ｘ２=Ｉｃｐ（ｘ，ｙ）として、この傾きの変化を反映するように、ラプラシアン成分Ｌ０〜Ｌ３を変換する変換関数Ｒ（ｘ，ｙ）を下記の式（６）により算出する。なお、Ｘ１＞Ｘ２であるため、変換関数Ｒ（ｘ，ｙ）は１より大きい値となる。

R(x,y)=(1/X2)/(1/X1)=X1/X2=(Ics(x,y)+Icp(x,y))/Icp(x,y) …（６）
なお、変換関数決定部２４は、放射線画像Ｇ４の各画素ごとに変換関数Ｒ（ｘ，ｙ）を算出して変換関数マップを作成する。この変換関数マップは、解像度が放射線画像Ｇ４と同一であるため、変換関数マップＲ４と称する。

次いで、変換部２５は、変換関数決定部２４が決定した変換関数Ｒ（ｘ，ｙ）により、ラプラシアン成分Ｌ０〜Ｌ３を変換する（ステップＳＴ５）。ここで、変換関数決定部２４が決定した変換関数マップＲ４のサイズは、ラプラシアン成分Ｌ３のサイズの１／２となっている。このため、３次Ｂスプライン補間あるいは線形補間等の周知の補間演算により変換関数マップＲ４を２倍に拡大して、ラプラシアン成分Ｌ３と同一サイズの変換関数マップＲ３を作成する。以下同様に、ラプラシアン成分Ｌ２と同一サイズの変換関数マップＲ２、ラプラシアン成分Ｌ１と同一サイズの変換関数マップＲ１、およびラプラシアン成分Ｌ０と同一サイズの変換関数マップＲ０をそれぞれ生成する。

そして、変換部２５は、下記の式（７）により、ラプラシアン成分Ｌ０〜Ｌ３を変換して、変換済みラプラシアン成分Ｌ０′〜Ｌ３′を算出する。但し、ｊ＝１〜３である。

Lj'(x,y)=Lj(x,y)×Rj(x,y) …（７）
ところで、被写体Ｍは、肺野、骨、筋肉および脂肪等の複雑な人体構造を有するため、散乱線成分を正確に求められない場合がある。このような場合において、上述した変換関数Ｒ（ｘ，ｙ）を用いてラプラシアン成分を変換した場合、骨部、治療具に含まれる金属等、コントラストが高い構造物が過強調となりやすい。このようにコントラストが過強調されると、構造物にアーチファクトが生じる等して、放射線画像の画質が低下してしまう。このため、ラプラシアン成分の値の絶対値が所定のしきい値を超えるとその値を抑制するように変換関数Ｒ（ｘ，ｙ）を決定するようにしてもよい。例えば、図７に示すように、入力値であるラプラシアン成分の絶対値がしきい値Ｔｈ１に近づくにつれて値が１に非線形に近づくようにし、しきい値Ｔｈ１を以上となると、入力値をそのまま変換された値とするように変換関数Ｒ（ｘ，ｙ）を決定するようにしてもよい。

続いて、再構成部２６が、散乱線が除去された放射線画像Ｇ４′（すなわち最低周波数帯域のガウシアン成分）および変換済みラプラシアン成分Ｌ０′〜Ｌ３′を再構成して処理済み放射線画像Ｇ０′を生成し（ステップＳＴ６）、処理を終了する。図８は再構成を説明するための図である。まず、再構成部２６は、散乱線が除去された放射線画像が、放射線画像ＧＬ４′の場合には、放射線画像ＧＬ４′を対数変換して、放射線画像Ｇ４′を生成する。そして、放射線画像Ｇ４′を補間演算により変換済みラプラシアン成分Ｌ３′と同一サイズ（すなわち、２倍）に拡大し、拡大した放射線画像Ｇ４′と変換済みラプラシアン成分Ｌ３′とを加算して、変換済みガウシアン成分Ｇ３′を生成する。

次いで、再構成部２６は、変換済みガウシアン成分Ｇ３′を変換済みラプラシアン成分Ｌ２′と同一サイズに拡大し、拡大した放射線画像Ｇ３′と変換済みラプラシアン成分Ｌ２′とを加算して、変換済みガウシアン成分Ｇ２′を生成する。以下、最高周波数帯域の変換済みラプラシアン成分Ｇ０′、すなわち処理済み放射線画像Ｇ０′が生成されるまで上記の処理を繰り返す。

ここで、散乱線成分を算出する放射線画像ＧＬ４は、画素値が放射線量と線形の関係にあるため、精度よく散乱線成分を算出することができ、その結果、放射線画像ＧＬ４または放射線画像Ｇ４から精度よく散乱線成分を除去することができる。一方、変換関数Ｒ（ｘ，ｙ）は、放射線量と線形の関係にある放射線画像ＧＬ４の散乱線成分から算出しているため、変換関数Ｒ（ｘ，ｙ）によりラプラシアン成分を精度よく強調できる。したがって、再構成された放射線画像においては、散乱線が除去され、かつ放射線画像Ｇ０に含まれる高周波成分が強調されるため、処理済みの放射線画像Ｇ０′においては、散乱線に起因するコントラストの低下が改善されたものとなる。

とくに、放射線画像Ｇ０は、画素値が対数放射線量と線形の関係にあるため、散乱線成分から変換関数Ｒ（ｘ，ｙ）を容易に決定することができ、かつ変換関数Ｒ（ｘ，ｙ）によりラプラシアン成分を変換することにより、効率よくラプラシアン成分を強調できる。したがって、放射線画像の散乱線に起因するコントラストの低下を、精度よく効率的に改善することができる。

なお、上記第１の実施形態においては、最低周波数帯域のガウシアン成分である放射線画像ＧＬから散乱線成分を算出して、変換関数マップを作成しているが、放射線画像Ｇ０をローパスフィルタによりフィルタリング処理して、放射線画像Ｇ０をぼかし、さらにぼかした放射線画像Ｇ０を逆対数変換して、画素値が放射線量と線形の関係にある放射線画像ＧＬ０を取得し、この放射線画像ＧＬ０から一次線成分および散乱線成分を算出するようにしてもよい。この場合、作成される変換関数マップは、元の放射線画像Ｇ０すなわちラプラシアン成分Ｌ０と同一サイズとなる。このため、変換関数決定部２４は、ラプラシアン成分Ｌ０と同一サイズの変換関数マップ（上記と同様にＲ０とする）を１／２に縮小して、ラプラシアン成分Ｌ１と同一サイズの変換関数マップＲ１を作成し、以下同様にラプラシアン成分Ｌ２と同一サイズの変換関数マップＲ２、ラプラシアン成分Ｌ３と同一サイズの変換関数マップＲ３をそれぞれ生成すればよい。また、散乱線除去部２３は、変換関数Ｃを放射線画像ＧＬ０から算出した一次線成分および散乱線成分を用いて算出して、放射線画像Ｇ４に対して散乱線除去処理を行えばよい。

また、上記第１の実施形態においては、すべての周波数帯域のラプラシアン成分を変換関数により変換しているが、これらのうちの少なくとも１つのラプラシアン成分を変換関数により変換するようにしてもよい。この場合、変換関数は変換する周波数帯域のラプラシアン成分に対してのみ決定すればよい。

次いで、本発明の第２の実施形態について説明する。図９は第２の実施形態による放射線画像撮影装置のコンピュータ内部の概略構成を示すブロック図である。図９に示すように、第２の実施形態においては、放射線画像を解析して、被写体Ｍの体厚分布を推定する体厚推定部２８を備え、放射線画像ＧＬ４の散乱線成分を算出する際に使用される体厚として、体厚推定部２８が推定した体厚を用いるようにした点が第１の実施形態と異なる。

図１０は体厚推定部の構成を示す概略ブロック図である。図１０に示すように、体厚推定部２８は、仮想モデル取得部３１、推定画像生成部３２、修正部３３および体厚分布決定部３４を備える。

仮想モデル取得部３１は、初期体厚分布Ｔ０（所定の体厚分布）を有する被写体Ｍの仮想モデルＫを取得する。

推定画像生成部３２は、仮想モデルＫに基づいて、仮想モデルの放射線撮影により得られる一次線画像を推定した推定一次線画像Ｉｐと、仮想モデルの放射線撮影により得られる散乱線画像を推定した推定散乱線画像Ｉｓとを合成した画像を、被写体Ｍの放射線撮影により得られる放射線画像を推定した推定画像Ｉｍとして生成する。

修正部３３は、推定画像Ｉｍと放射線画像ＧＬ４とに基づいて、推定画像Ｉｍと放射線画像ＧＬ４の違いが小さくなるように仮想モデルＫの初期体厚分布Ｔ０を修正する。

体厚分布決定部３４は、修正された体厚分布Ｔｎ−１（ｎは自然数）を放射線画像ＧＬ４の体厚分布Ｔｋに決定する。

また、第２の実施形態においては、記憶部２７には、初期体厚分布Ｔ０（ｘ，ｙ）を有する被写体Ｍの仮想モデルＫが記憶される。なお、体厚とは、照射された放射線の経路上における空気領域を除いた被写体領域の厚さの総計を意味する。

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

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

次いで推定画像生成部３２は、放射線画像ＧＬ４を取得した放射線画像Ｇ０と同等の撮影条件で仮想モデルＫを撮影した場合に得られる推定一次線画像Ｉｐと、放射線画像Ｇ０と同等の撮影条件で仮想モデルＫを撮影した場合に得られる推定散乱線画像Ｉｓとを合成した推定画像Ｉｍを生成する（ステップＳＴ３２）。図１３および図１４は、推定画像Ｉｍの生成方法を説明するための図である。

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

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

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

なお、Ｋｓ(ｘ，ｙ，Ｔｎ（ｘ’，ｙ’），θｘ’，ｙ’)、Ｋｐ＋ｓ（ｘ，ｙ，Ｔｎ−１（ｘ’，ｙ’），θｘ’，ｙ’)は、撮影条件等に応じて実験的に求めることができる。

本実施形態においては、撮影時の撮影条件に基づいてカーネルＫｓ(ｘ，ｙ，Ｔｎ（ｘ’，ｙ’），θｘ’，ｙ’)、Ｋｐ＋ｓ（ｘ，ｙ，Ｔｎ−１（ｘ’，ｙ’），θｘ’，ｙ’)を算出してもよいが、各種撮影条件とカーネルＫｓ(ｘ，ｙ，Ｔｎ（ｘ’，ｙ’），θｘ’，ｙ’)、Ｋｐ＋ｓ（ｘ，ｙ，Ｔｎ−１（ｘ’，ｙ’），θｘ’，ｙ’)とを対応付けたテーブルを記憶部２７に記憶しておき、撮影時の照射野情報、被写体情報および撮影条件に基づいて、このテーブルを参照してカーネルＫｓ(ｘ，ｙ，Ｔｎ（ｘ’，ｙ’），θｘ’，ｙ’)、Ｋｐ＋ｓ（ｘ，ｙ，Ｔｎ−１（ｘ’，ｙ’），θｘ’，ｙ’)を求める。

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

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

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

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

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

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

修正された体厚分布Ｔｎを取得すると、体厚分布決定部３４は、ｎの値を１つ増加して更新し（ｎ＝ｎ＋１とする）、仮想モデル取得部３１は修正された体厚分布Ｔｎを取得する（ステップＳＴ３１）。そして、取得された体厚分布Ｔｎに対して、推定画像生成部３２および体厚分布決定部３４はステップＳＴ３１〜ステップＳＴ３３の処理をそれぞれ上記と同様に実行する。そして、放射線画像ＧＬ４と推定画像Ｉｍとの違いを示すエラー値Ｖｅｒｒｏｒが終了条件を満たすまで、上記同様に、体厚分布Ｔｎの修正処理（ステップＳＴ３４）と、修正された体厚分布Ｔｎを有する仮想モデルＫの取得処理（ステップＳＴ３１）と、体厚分布Ｔｎを用いた新たな推定画像Ｉｍの生成処理（ステップＳＴ３２）と、新たに生成された推定画像Ｉｍと放射線画像ＧＬ４との違いが終了条件を満たすかを判定する処理（ステップＳＴ３３）の処理とが繰り返される。

一方、体厚分布決定部３４は、エラー値Ｖｅｒｒｏｒが終了条件を満たしていることを判定した場合には（ステップＳＴ３３：Ｙｅｓ）、終了条件を満たした際にエラー値Ｖｅｒｒｏｒに用いられた体厚分布Ｔｎを放射線画像ＧＬ４の体厚分布Ｔｋとして決定して体厚推定処理を終了する（ステップＳＴ３５）。

第２の実施形態においては、体厚推定部２８が推定した体厚分布Ｔｎを用いて、上記式（１）、（２）により、一次線成分Ｉｃｐおよび散乱線成分Ｉｃｓを算出し、算出した一次線成分Ｉｃｐおよび散乱線成分Ｉｃｓを用いて、上記第１の実施形態と同様に、散乱線除去処理、変換関数の決定、ラプラシアン成分の変換および再構成を行うことにより、処理済み放射線画像Ｇ０′を取得するようにしたものである。このため、推定された体厚分布に応じてより適切に散乱線成分を除去することができるため、より高画質の放射線画像を取得することができる。

なお、上記第１および第２の実施形態において、散乱線除去処理を行う際に、放射線画像の撮影時に散乱線を除去するために使用が想定される仮想的なグリッドの特性である仮想グリッド特性を取得し、仮想グリッド特性に基づいて、放射線画像の散乱線除去処理を行うようにしてもよい。以下、仮想グリッド特性を用いた散乱線除去処理を第３の実施形態として説明する。図１５は第３の実施形態において行われる処理を示すフローチャートである。散乱線除去部２３は、第１の実施形態と同様に放射線画像Ｇ４を逆対数変換して、各画素の画素値が放射線量と線形の関係にある放射線画像ＧＬ４を生成する（ステップＳＴ４１）。次いで、散乱線除去部２３は、放射線画像ＧＬ４から散乱線の影響を除去するために、被写体の撮影時に散乱線を除去するために使用が想定される仮想的なグリッドの特性である仮想グリッド特性を取得する（ステップＳＴ４２）。本実施形態においては、仮想グリッド特性は、仮想グリッドについての散乱線透過率Ｔｓ、および被写体Ｍを透過して直接放射線検出器５に照射される一次線の透過率（一次線透過率）Ｔｐとする。なお、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは０〜１の間の値をとる。

散乱線除去部２３は、入力部８から散乱線透過率Ｔｓおよび一次線透過率Ｔｐの値の入力を直接受け付けることにより仮想グリッド特性を取得してもよいが、本実施形態においては、撮影条件に基づいて仮想グリッド特性、すなわち散乱線透過率Ｔｓおよび一次線透過率Ｔｐを取得する。前者の場合は入力部８が、後者の場合は制御部２１が，それぞれ特性取得手段に相当する。

仮想グリッド特性を取得するための撮影条件は、撮影時の撮影距離（ＳＩＤ）、管電圧、ｍＡｓ値（管電流時間積）、線源のターゲットおよびフィルタの材質、並びに撮影に使用される放射線検出器の種類等のうちの少なくとも１つを含む。ここで、放射線画像の撮影時には、撮影条件に応じて使用するグリッドの種類が定められている。また、グリッドは、グリッド比、グリッド密度、収束型か平行型か、収束型の場合の集束距離、インタースペース素材（アルミニウム、ファイバー、ベークライト等）等に応じて種類が異なり、グリッドの種類に応じて、散乱線透過率Ｔｓおよび一次線透過率Ｔｐは異なるものとなる。このため、本実施形態においては、各種撮影条件と仮想グリッド特性とを対応付けたテーブルが記憶部２７に記憶されている。散乱線除去部２３は、このテーブルを参照して、撮影条件から仮想グリッド特性を取得する。

そして、散乱線除去部２３は、第１の実施形態と同様に、放射線画像ＧＬ４の一次線成分および散乱線成分を取得する（ステップＳＴ４３）。次いで、散乱線除去部２３は、仮想グリッド特性である散乱線透過率Ｔｓおよび一次線透過率Ｔｐ、並びに上記式（３）により算出した変換関数、すなわち散乱線含有率Ｓ（ｘ，ｙ）から、逆対数変換前の放射線画像Ｇ４の各画素の変換関数Ｑ（ｘ，ｙ）を下記の式（２２）により算出する。

Q(x,y) = S(x,y)×Ts + (1-S(x,y))×Tp …（２２）
散乱線透過率Ｔｓおよび一次線透過率Ｔｐ、並びに散乱線含有率分布Ｓ（ｘ，ｙ）は０〜１の間の値となるため、変換関数Ｑ（ｘ，ｙ）も０〜１の間の値となる。散乱線除去部２３は、下記の式（２３）により、変換関数Ｑ（ｘ，ｙ）を放射線画像Ｇ４の各画素の画素値Ｇ４（ｘ，ｙ）に乗算することにより散乱線除去処理を行い、散乱線が除去された放射線画像Ｇ４′を算出し（ステップＳＴ４４）、散乱線除去処理を終了する。

G4'(x,y)=G4(x,y)×Q(x,y) …（２３）
このように、第３の実施形態においては、仮想グリッド特性を用いて散乱線除去処理を行うようにしたため、使用が想定されるグリッドを用いて撮影を行った場合と同様の散乱線除去の効果を得ることができる。

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

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

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

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

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

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

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

また、上記第１から第３の実施形態において、処理済みの放射線画像および処理前の放射線画像の双方を表示し、いずれの放射線画像を診断に使用するかを選択できるようにしてもよい。

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

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

１ 撮影装置
２ コンピュータ
３ Ｘ線源
５ 放射線検出器
６ 表示部
８ 入力部
２１ 制御部
２２ 周波数分解部
２３ 散乱線除去部
２４ 変換関数決定部
２５ 変換部
２６ 再構成部
２７ 記憶部
２８ 体厚推定部
３１ 仮想モデル取得部
３２ 推定画像生成部
３３ 修正部
３４ 体厚分布決定部

Claims (11)

1. 被写体に放射線を照射することにより撮影された、第１の放射線画像を周波数分解して、該第１の放射線画像の複数の周波数成分を表す複数の帯域画像を取得する周波数分解手段と、
   前記放射線の放射線量と線形の関係にある前記被写体の第２の放射線画像を用いて、該第２の放射線画像に含まれる前記被写体による散乱線成分を算出し、該散乱線成分を該第２の放射線画像から除去して散乱線除去済み放射線画像を取得する散乱線除去手段と、
   前記第２の放射線画像に含まれる前記散乱線成分に応じて、少なくとも１つの前記帯域画像を変換する変換関数を決定する変換関数決定手段と、
   前記少なくとも１つの帯域画像を前記変換関数により変換して変換済み帯域画像を生成する変換手段と、
   前記散乱線除去済み放射線画像および前記変換済み帯域画像を再構成して処理済み放射線画像を取得する再構成手段とを備えたことを特徴とする放射線画像処理装置。
2. 前記変換関数決定手段は、前記散乱線成分と前記放射線に含まれる一次線成分および前記散乱線成分の和との比率に基づいて、前記変換関数を決定する手段であることを特徴とする請求項１記載の放射線画像処理装置。
3. 前記散乱線除去手段は、前記被写体の体厚分布を用いて、前記散乱線成分および前記一次線成分を算出し、
   前記変換関数決定手段は、該算出された散乱線成分および一次線成分から前記比率を算出する手段であることを特徴とする請求項２記載の放射線画像処理装置。
4. 前記第２の放射線画像に基づいて、前記体厚分布を推定する体厚推定手段をさらに備えたことを特徴とする請求項３記載の放射線画像処理装置。
5. 前記第１の放射線画像は、画素値が対数放射線量と線形の関係にあることを特徴とする請求項１から４のいずれか１項記載の放射線画像処理装置。
6. 前記第２の放射線画像は、前記第１の放射線画像を逆対数変換した画像であることを特徴とする請求項５記載の放射線画像処理装置。
7. 前記再構成手段は、前記散乱線除去済み放射線画像を対数変換し、該対数変換された散乱線除去済み放射線画像および前記変換済み帯域画像を再構成して前記処理済み放射線画像を取得する手段であることを特徴とする請求項１から６のいずれか１項記載の放射線画像処理装置。
8. 前記被写体の撮影時に、散乱線を除去するために使用が想定される、仮想的なグリッドの特性である仮想グリッド特性を取得する特性取得手段をさらに備え、
   前記散乱線除去手段は、前記仮想グリッド特性に基づいて、前記散乱線成分を算出する手段であり、
   前記変換関数決定手段は、前記算出された散乱線成分に応じて前記変換関数を決定する手段であることを特徴とする請求項１から７のいずれか１項記載の放射線画像処理装置。
9. 前記変換関数決定手段は、前記帯域画像の画素値が所定のしきい値を越えると、該画素値を抑制する非線形関数となるように前記変換関数を決定する手段であることを特徴とする請求項１から８のいずれか１項記載の放射線画像処理装置。
10. 被写体に放射線を照射することにより撮影された、第１の放射線画像を周波数分解して、該第１の放射線画像の複数の周波数成分を表す複数の帯域画像を取得し、
    前記放射線の放射線量と線形の関係にある前記被写体の第２の放射線画像を用いて、該第２の放射線画像に含まれる前記被写体による散乱線成分を算出し、
    該散乱線成分を該第２の放射線画像から除去して散乱線除去済み放射線画像を取得し、
    前記第２の放射線画像に含まれる前記散乱線成分に応じて、少なくとも１つの前記帯域画像を変換する変換関数を決定し、
    前記少なくとも１つの帯域画像を前記変換関数により変換して変換済み帯域画像を生成し、
    前記散乱線除去済み放射線画像および前記変換済み帯域画像を再構成して処理済み放射線画像を取得することを特徴とする放射線画像処理方法。
11. 被写体に放射線を照射することにより撮影された、第１の放射線画像を周波数分解して、該第１の放射線画像の複数の周波数成分を表す複数の帯域画像を取得する手順と、
    前記放射線の放射線量と線形の関係にある前記被写体の第２の放射線画像を用いて、該第２の放射線画像に含まれる前記被写体による散乱線成分を算出する手順と、
    該散乱線成分を該第２の放射線画像から除去して散乱線除去済み放射線画像を取得する手順と、
    前記第２の放射線画像に含まれる前記散乱線成分に応じて、少なくとも１つの前記帯域画像を変換する変換関数を決定する手順と、
    前記少なくとも１つの帯域画像を前記変換関数により変換して変換済み帯域画像を生成する手順と、
    前記散乱線除去済み放射線画像および前記変換済み帯域画像を再構成して処理済み放射線画像を取得する手順とをコンピュータに実行させることを特徴とする放射線画像処理プログラム。