JP4194124B2 - 放射線画像取得装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、物質の放射線透過分布を画像化するＸ線画像撮像装置等に用いて好適な放射線画像取得装置に関するものである。
【０００２】
【従来の技術】
放射線源から発する放射線を被写体である物体に照射し、透過した放射線分布を電気信号に変換し、ディジタル変換することによりディジタル画像を取得するような、被写体内部情報を画像化するシステムにおいては、被写体内部での放射線の散乱現象のために、得られる画像のコントラストが大きく低下するなどのエラーが発生する。このためエラーをなくすための方法について従来よりさまざまな検討、研究がなされている。
【０００３】
その方法として、検出面において鉛グリッドを使用する方法は、かなりの確率で除去できるが、グリッドを使用したためのグリッド像の除去が困難な点や検出面での検出効率の著しい低下が見られ、十分な解決方法であるとはいえない。また検出面と物体との間に空気間隙を設けて散乱線強度を低下させる方法は、完全に散乱線を除去できない場合もあり、また距離をおくための検出線量低下、放射線源が完全な点線源でないためのボケの拡大や、より大きな検出面が必要になるなどの問題がある。また散乱を点分布関数として捉え、それを演算により逆フィルタをかけて除去する方法は散乱現象そのものが画像分布に線形に依存するような安定したものではないため、その方法での解決は困難である。
【０００４】
次に放射線源から発する放射線を被写体である物体に照射し、透過した放射線分布の直接線成分を被写体内部情報として得る操作を変調の立場から考え直してみる。即ち、上記操作を放射線源から発する放射線分布をキャリアとして被写体内部情報によって変調し、伝達する行為とする。具体的にはキャリアとしての放射線分布をｃ（ｘ，ｙ）、（ｘ、ｙは分布を構成する２次元平面上での位置を表わす）、ｆ（ｘ，ｙ）を伝達すべき信号情報としての被写体内部の放射線透過率の投影情報とすれば、被写体を通過した信号は変調された信号として両者の積として得られる。変調に対する復調操作は得られた変調信号のキャリアから信号情報を取り出すことを意味するので、この場合は放射線分布ｃ（ｘ，ｙ）で除去することによって復調を行う。もしくは、対数変換により乗算関係を加算関係に変換し、キャリア成分の減算によって行う。
【０００５】
一般的には放射線分布はほぼ完全に一定値の一様分布であるので、キャリアとしての意識はされず復調としての除算を行わない場合が多い、即ちキャリアを含んだ形で信号情報とする。但し、検出系の感度分布が一定でない場合、キャリアｃ（ｘ，ｙ）に感度分布を含め、一様でないとして除算する場合もある。
【０００６】
しかしながら実際のシステムでは、図９に示すように、放射線の散乱やノイズの影響を受けることになる。図９は上記システムを原理に示すものであり、４１１は被写体、４１２は放射線分布を電気信号に変換する変換装置、４１０は放射線源から発した放射線、４１３は被写体４１１を通過した放射線、４１４は変換装置４１２で変換された電気信号である。また放射線源から発する放射線分布をｃ（ｘ，ｙ）とし、ｆ（ｘ，ｙ）を被写体内部の放射線透過率の投影情報、ｑ（ｘ，ｙ）を放射線をエネルギー粒子として考えた時の量子ノイズ、ｓ（ｘ，ｙ）を被写体内部での散乱線によって新たに発生した放射線分布、ｎ（ｘ，ｙ）を電気回路などで発する熱雑音を主体とする電気ノイズであるとする。
【０００７】
図９において、得られる電気信号としての放射線分布ａ（ｘ，ｙ）は、
ａ（ｘ，ｙ）＝ｃ（ｘ，ｙ）・ｆ（ｘ，ｙ）＋ｑ（ｘ，ｙ）＋ｓ（ｘ，ｙ）＋ｎ（ｘ，ｙ） ………（１）
と表わされる。
ここで、量子ノイズｑ（ｘ，ｙ）、電気回路などのノイズｎ（ｘ，ｙ）は量的に比較的少なく、予測も不可能であることから、ここでは無視することにしても、
ａ（ｘ，ｙ）＝ｃ（ｘ，ｙ）・ｆ（ｘ，ｙ）＋ｓ（ｘ，ｙ）………（２）
と散乱成分は残る。
上記式（２）のような信号から所望するｆ（ｘ，ｙ）を得ようとするならば、散乱成分ｓ（ｘ，ｙ）を差し引いいた後に、キャリアｃ（ｘ，ｙ）を取り除くしかない。
【０００８】
従来はキャリアｃ（ｘ，ｙ）を一様分布と考えることが支配的であったため、何らかの方法でｓ（ｘ，ｙ）を推定し、差し引くだけで解決しようとしたものが多かった。その方法として例えば、特願平３−２２１８２６号出願明細書や特開平７−１２９４８号公報にあるように、主線を外れたところに別に設けた放射線検出器によって散乱線を推測する方法や、Medical Physics 1988;15:pp289-297 にあるように、画素値のグレーレベルと散乱成分との仮説的な関係を仮定して散乱線を推測する方法などがあるが、いずれも散乱線成分を完全に推測することは不可能であり、精度に問題があった。
【０００９】
また、従来ではキャリアとしての放射線分布を一様とせずにさ散乱成分を推測する方法も検討されている。例えばMedical Physics 1988;15-732-748 にあるように、被写体と放射線源との間に小型鉛円盤を置いた場合と無い場合との２回の画像の取得により、鉛円盤下の情報は散乱によるものとして散乱成分を推定する方法があるが、２回の放射線爆射が必要なため、医療用に人体に用いる場合は被爆量が増大することや、２回の爆射間で患者の位置変動が許されないことなどの問題点がある。
【００１０】
また特開平７−２０９７７８号公報にあるように、放射線をある程度透過する既知の物質と既知の厚さを持つ複数の物体を被写体と放射線源との間に置き、被写体とともに放射線爆射して得られた既知の物質による影響が既知の物質と厚みに影響によるものと一致するような補正数値演算を行い、その補正値を散乱線成分と推定し、複数部分を補間して全体の散乱線分布を求める方法があるが、用いる物質の厚みを正確に知らなければならないことや、事前に同様な物質を用いた構成操作が必要であることなどの問題がある。
【００１１】
また、SPIE Vol.1651 Medical Imaging VI:Instrumentation(1992)pp222-233 にあるように、直接線である１次放射線分布のみを被写体の放射線源側に置かれたフィルタにより変調し、同一のフィルタを検出器側において復調することにより、直接線のみを得ることを可能にしているが、変調器と復調器を完全に整合することは実際上不可能であり、精度の高い方法であるとは言えない。
【００１２】
これらのキャリアとしての放射線分布を一様とせずに、散乱線分布を推定しようとする従来の方法は優れた観点から得られた方法ではあるが、放射線分布を一様としなかったための弊害を問題にしていない。それは、前記（２）式で示した変調の考え方に戻るが、変調信号を伝達する場合のキャリア自身の劣化である。具体的にはキャリア分布ｃ（ｘ，ｙ）のボケの発生である。即ち、（２）式はキャリアｃ（ｘ，ｙ）のボケた値ｃ′（ｘ，ｙ）により変調されて次式のようになる。
ａ（ｘ，ｙ）＝ｃ′（ｘ，ｙ）・ｆ（ｘ，ｙ）＋ｓ（ｘ，ｙ）………（３）
【００１３】
このｃ′（ｘ，ｙ）も散乱線と同様に被写体によるものであり、一意的な推測は困難である。もしボケが無ければ、ｃ（ｘ，ｙ）は被写体を置かない状態で放射線爆射して得られ、簡単に補正できる。また、キャリアｃ（ｘ，ｙ）が一様分布であるかもしくは近ければボケの影響を受けないためキャリアの劣化はなく、そのことは考えなくてよい。
【００１４】
また前記特開平７−２０９７７８号公報の方法においては、既知の物質と既知の厚さを持つ複数の物体の透過率の違いを用いて散乱線成分を求め、さらに補間によって滑らかな変動を持つ成分として全体の散乱線成分を求めた後、差し引くのみであり、キャリアとしての放射線分布の変動分のボケが考慮されていない。従って物質の境界部分の影響が残留した形でしか画像が得られない。
【００１５】
【発明が解決しようとする課題】
以上述べて来た従来の各問題点をまとめると、キャリアである放射線分布が一様な状態の場合は、散乱線成分を差し引けば、信号成分である被写体の透過率分布は求まる反面、その散乱線成分自体を精度良く求める方法がない。また、放射線分布を一様とせず、散乱線成分を際立たせるような分布形状にすれば、散乱線成分は精度良く求まる可能性はあるが、キャリアである放射線分布も散乱現象と同様な理由によりボケるため、今度はキャリアの除去が困難になるということになる。また、特に全ての従来では、被写体の無い状態で撮影した一様でない特定の分布を持った放射線分布をキャリアとして、除算を用いてキャリアを除去するという方式のものはない。
【００１６】
本発明は上記の問題点を解決するために成されたもので、一様でない放射線分布を用いた場合における散乱線の影響をなくし、良好な画質を得ることを目的としている。
【００１７】
【課題を解決するための手段】
本発明の放射線画像取得装置は、異なる厚みもしくは物質からなる少なくとも２つの部分で構成されたフィルタを透過した放射線が被写体を透過した場合に、前記被写体の散乱線分布の画像を、座標を（ｘ、ｙ）とする画像ｓ（ｘ，ｙ）として推定する放射線画像取得装置において、前記フィルタを透過した放射線の放射線分布を画像にする撮像手段と、前記フィルタと前記撮像手段との間に前記被写体がない状態で前記撮像手段において撮像された画像に基づいて、画像ｃ（ｘ，ｙ）を取得し、前記フィルタと前記撮像手段との間に前記被写体がある状態で前記撮像手段において撮像された画像に基づいて、画像ａ（ｘ，ｙ）を取得し、前記画像ｃ（ｘ，ｙ）及び前記画像ａ（ｘ，ｙ）において、前記フィルタの前記２つの部分の境界部分に対応する部分分布を、それぞれ、｛ｃ_i′（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝、及び、｛ａ_i′（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝（ただし、Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1iは画像上で前記２つの部分における境界部分を含む座標であり、ｋは部分分布の個数）として抽出し、前記ａ_i′（ｘ_i，ｙ_i）に対応する前記散乱線分布の部分をｓ_i′（ｘ_i，ｙ_i）とし、ｆ_i′（ｘ_i，ｙ_i）＝｛ａ_i′（ｘ_i，ｙ_i）−ｓ_i′（ｘ_i，ｙ_i）｝／ｃ_i′（ｘ_i，ｙ_i）で表せる式において、｛（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝の範囲で、前記ｆ_i′（ｘ_i，ｙ_i）の合計値が最小となる前記ｓ_i′（ｘｉ，ｙｉ）の分布関数を、それぞれの部分（０≦ｉ≦ｋ−１）で推定し、当該推定により得られた分布関数に応じて散乱値ｓ′ｍ（０≦ｍ≦ｋ）の値を求め、求めた前記散乱値ｓ′ｍの値に基づいて、補間により前記画像ｓ（ｘ，ｙ）を推定する処理を行う処理手段とを有する。
【００２６】
【発明の実施の形態】
以下、本発明の実施の形態を図面と共に説明する。
まず、本発明を原理的に説明する。
前述したような課題に対して、本発明ではキャリアとしての放射線分布を、予想される被写体情報と明確に分離できる分布形状にするとともに、被写体の無い状態でその分布形状をキャリアとして事前に取得し、被写体で変調された情報から信号情報のみを取り出すために、上記取得したキャリアによって除算を行う場合、第１段階として変調された信号から適当な成分を減算した後、上記取得したキャリアで除算した結果にキャリアの特徴成分が最小になる減殺成分を探索する。例えば画像であれば複数の位置でその探索を行い、それらを補間して散乱成分を推定する。
【００２７】
第２段階として、第１段階で取り去れなかったキャリアの特徴成分は、散乱による影響ではなくキャリア自身の伝送路中のボケによる劣化の影響で残留したキャリアの特徴成分として、その特徴成分がなくなるまで、画像を補修する。その補修の方法は種々考えられるが、上記取得したキャリアを適当な関数で数値的にボカした後、第１段階で得られた減算された信号成分で除算を行い、その結果中のキャリアの特徴成分が最小になる関数を探索することで実現可能である。またこの方法で行うならば、第１段階と第２段階の順序は逆転してもよく、また、第１段階第２段階を繰り返し行って最適な散乱成分およびボケ関数を見出すことも可能である。また、第１段階のみでキャリアの影響が完全に近い形で除去できるならば第１段階のみで終了しても構わない。
【００２８】
以上の操作により、散乱成分による影響を取り出すためにキャリアを一様ではなくし、散乱成分を取り出すとともに、伝送路中の伝送関数の影響を受けやすい一様ではないキャリアの劣化の影響を取り除き、信号成分である被写体の透過率の投影情報を取り出す。例えば医療用のＸ線撮像装置の場合、被写体は人体であり、キャリア分布は人体内部にほとんどありえない、例えば幾何学的形状であって急峻に変動するような分布にしておく。
【００２９】
その概略の過程は、
１）複数の急峻な空間分布を持つ放射線源の画像を被写体無しで撮影取得し、｛ｃ（ｘ，ｙ）｜０≦Ｘ≦ｎ−１，０≦Ｙ≦ｍ−１｝（Ｘ、Ｙは画素数）とする。
２）１）と同じ放射線源によって被写体の画像｛ａ（ｘ，ｙ）｜０≦ｘ≦Ｘ−１，０≦ｙ≦Ｙ−１｝を取得する。
３）ｃ（ｘ，ｙ）、ａ（ｘ，ｙ）のなかで上記急峻な空間分布の部分に対応する部分分布｛ｃ_i ′（ｘ_i ，ｙ_i ）｜Ｘ_0i≦ｘ_i ≦Ｘ_1i，Ｙ_0i≦ｙ_i ≦Ｙ_1i，０≦ｉ≦ｋ−１｝，｛ａ_i ′（ｘ_i ，ｙ_i ）｜Ｘ_0i≦ｘ_i ≦Ｘ_1i，Ｙ_0i≦ｙ_i ≦Ｙ_1i，０≦ｉ≦ｋ−１｝（ｋは部分分布の個数）を切り出す。
【００３０】
４）ａ _i ′（ｘ_i，ｙ_i）に対応する適当な分布関数ｓ _i ′（ｘ_i，ｙ_i）を定め、
ｆ _i ′（ｘ_i，ｙ_i）
＝｛ａ _i ′（ｘ_i，ｙ_i）−ｓ _i ′（ｘ_i，ｙ_i）｝／ｃ _i ′（ｘ_i，ｙ_i）
（Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１） ………（４）
ｆ _i ′（ｘ_i，ｙ_i）の中で上記急峻な変動の影響が最小となるｓ′（ｘ，ｙ）を求める。
５）複数得られたｓ′（ｘ，ｙ）に対して、補間操作によって領域全体の散乱線分布｛ｓ（ｘ，ｙ）｜０≦Ｘ≦ｎ−１，０≦Ｙ≦ｍ−１｝を推定する。
６）散乱線分布を差し引いた被写体の透過率分布ｃ（ｘ，ｙ）を下式で求める。
ｆ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ） ………（５）
１０）（４）式で除去しきれない、急峻な分布変動に起因する画像の誤差を補修する。
【００３１】
線源の放射線分布の急峻な変動が除去できない場合の理由としては以下が考えられる。
ａ）散乱線の分布がオフセットとして検出され、除算に誤差が出る。
ｂ）急峻な線源の変動が被写体の散乱によってぼかされ、検出分布の急峻さがなくなる。
ｃ）被写体なしで線源の分布を取得したときと、実際の被写体の透過率分布を取得したときとで線源の分布の位置がずれる。
【００３２】
前述の過程でオフセットとしての散乱線分布が除去できるため、ａ）の影響は取り除かれる。またｂ）の影響は、線源の分布に急峻な変動成分をおいたための副作用であり、周りの画像値などから推測してその影響を補修する必要がある。ｃ）の影響を出さないためには線源と検出面の位置関係を確実にし、被写体なしで分布を取得したときからその位置関係を不変なものにする必要がある。この放射線分布取得過程を例えばコーンビームＣＴに用いる場合には、線源と検出面との位置関係は不変であるので、この条件は満たしやすい。散乱線分布がなければ、急峻な分布変動があってもその影響は単純な除算で除去できるはずであるが、オフセットとして散乱線分布があるため、その変動分を除去仕切れない。また散乱線分布ｓ（ｘ，ｙ）は比較的滑らかな分布を示すため、補間処理によってほぼ再現が可能である。
【００３３】
上述の過程４）で急峻な変動の影響の抽出方法としては、例えば急峻な変動の前後の値の差などを抽出する方法とか、得られた部分画像の分散値を見る方法などが考えられる。
例えば、急峻な変動の前後の差を見る場合は簡単な例であるが、部分分布｛ａ_i ′（ｘ_i ，ｙ_i ）｜Ｘ_0i≦ｘ_i ≦Ｘ_1i，Ｙ_0i≦ｙ_i ≦Ｙ_1i，０≦ｉ≦ｋ−１｝において、対応する散乱線分布｛ｓ_i ′（ｘ_i ，ｙ_i ）｜Ｘ_0i≦ｘ_i ≦Ｘ_1i，Ｙ_0i≦ｙ_i ≦Ｙ_1i，０≦ｉ≦ｋ−１｝はほぼ一定値ｓ_i として求める。
【００３４】
急峻な変動の前後の画素値をそれぞれａ′_0i、ａ′_1i、ｃ′_0i、ｃ′_1iとし、
ｓ′_i ＝（ｃ′_1i・ａ′_0i−ｃ′_0i・ａ′_1i）／（ｃ′_0i−ｃ′_1i）
………（６）
としてｓ′_i （０≦ｉ≦ｋ）が求まり、それぞれを補間して全体の｛ｓ（ｘ，ｙ）｜０≦ｘ≦Ｘ，０≦ｙ≦Ｙ｝を求める。
【００３５】
現実的には、変動の前後で被写体の分布が一定であるかどうかは不明であるため、具体的な演算では、除算結果の部分分布中の分散値ε_m を最小化するｓ′_m を求める。即ち、
【００３６】
【数５】

【００３７】
である（Ｎは部分分布中の画素数）。一定の条件でｓ′_m を求めると以下の式になる。
【００３８】
【数６】

【００３９】
上式で散乱値ｓ′_m （０≦ｍ≦ｋ）を求め、各点を補間して、散乱分布｛ｓ（ｘ_i ，ｙ_i ）｜０≦ｘ_i ≦Ｘ，０≦ｙ_i ＝Ｙ｝を得る。ｓ′_m を求める方法としは、これらに限ったことではなく、除算した結果に急峻な変動の影響を最も少なくするような値を求める他のさまざまな方法が考えられる。
【００４０】
次に、１０）の｛ｓ（ｘ_i ，ｙ_i ）｜０≦ｘ_i ≦Ｘ，０≦ｙ_i ≦Ｙ｝の差分のみでは除去しきれない成分の補修は、
｛ｃ_i ′（ｘ_i ，ｙ_i ）｜Ｘ_0i≦ｘ_i ≦Ｘ_1i，Ｙ_0i≦ｙ_i ≦Ｙ_1i，０≦ｉ≦ｋ−１｝を適当な関数ｈ_m （ｘ，ｙ）との畳み込み積分（convolution ）を演算する。すなわち、
【００４１】
【数７】

【００４２】
ここで求まった、ｄ（ｘ_i ，ｙ_i ）を新たなキャリアとして再び、分数値を計算する。
【００４３】
【数８】

【００４４】
このε_m が最小になる（１０）式内のｈ（ｘ，ｙ）を探索する。具体的にはｈ（ｘ，ｙ）をバイリニア関数でもしくはガウス分布などで近似し、関数のスパンを除々に変化させながら探索する。もしキャリアの変動が一次元のみの方向性しかなければ、上記畳み込み積分も一次元の演算でよく演算速度も増す。
この方法以外にも、１０）において単純に画像中にキャリアにある急峻な変動による影響をフィルタリングなどの画像処理で補修する方法も考えられる。
【００４５】
さらに、他の方法としてキャリアが被写体によって影響を受け難い、すなわちボケの影響が出難い分布にしておき、前述の例であるキャリアのボケによる影響を補修を行わない次のような、方法も考えられる。
１）２次元的に直線的に傾斜変動する放射線源分布ｃ（ｘ，ｙ）を構成する。例えば
ｃ（ｘ，ｙ）＝ｋ１・（ｘ＋ｙ）＋ｋ２
ｋ１、ｋ２は定数
ｘ、ｙは２次元的な位置
であり、放射線源の画像を被写体無しで撮影取得する。
２）１）と同じ放射線源によって被写体の画像｛ａ（ｘ，ｙ）｜０≦ｘ≦Ｘ−１，０≦ｙ≦Ｙ−１｝を取得する。
３）ｃ（ｘ，ｙ）、ａ（ｘ，ｙ）を小画像の部分画像に分割する。
【００４６】
４）各部分画像ａ（ｘ_i ，ｙ_i ）に対応する適当なオフセット値ｓを定め
ｆ（ｘ_i ，ｙ_i ）＝｛ａ（ｘ_i ，ｙ_i ）−ｓ｝／ｃ（ｘ_i ，ｙ_i ）
の中で、分布の変動分が最も小さくなるｓの値を見出す。
５）複数得られたｓに対して、補間操作によって画像全体の散乱線分布｛ｓ（ｘ，ｙ）｜０≦Ｘ≦ｎ−１，０≦Ｙ≦ｍ−１｝を推定する。
６）散乱線分布を差し引いた被写体の透過率分布ｃ（ｘ，ｙ）を下式で求める
ｆ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）
【００４７】
この方式の利点は、キャリアの分布に急峻に変動する部分がなく、キャリアのボケの影響が少ないことである。また、直線的にキャリア強度を分布させると、キャリア取得位置による影響がなくなり（即ち、位置によるキャリアの違いはオフセット値のみで傾きは不変）、過程１）のキャリア取得時と２）の被写体取得時の影響位置が多少ずれても補正が行えることにある。問題点としては、小画像に分割したときに、小画像の大きさが小さいと、過程４）の分布の変動分の検出精度が鈍くなり、大きいと散乱線分布ｓ（ｘ，ｙ）のサンプリング精度が鈍くなる。また小画像の中の被写体画像自身の変動があまりにも大きいと散乱線成分の検出精度が鈍くなる可能性がある。また、
【００４８】
【数９】

【００４９】
を用いるのが簡便である。
【００５０】
次に上述した原理に基づく本発明の実施の形態について説明する。
まず第１の実施の形態を説明する。
図１は本発明を実施したＸ線画像撮影装置の実施の形態のブロック図であり、１はＸ線発生装置を示し、２は被写体である人体、３はキャリアであるＸ線量分布を一様にしないためのＸ線を１部透過するフィルタであり後述する。４はＸ線の放射分布を電気信号に変換する撮像素子であり、５は変換された電気信号をディジタル値に変換するＡ／Ｄ変換器、８は変換されたディジタル値を一旦画像の形態で記憶するフレームメモリ、９はキャリアとしてのＸ線放射分布を記録するフレームメモリ、１０は撮像素子４の暗出力値を補正するための暗出力値を画像として記憶するフレームメモリ、１１は計算によって求まった散乱線分布を記憶するフレームメモリである。
【００５１】
６はシステム全体を管理する中央処理装置であり、記憶媒体としてのメモリ７に格納されたプログラムによって動作する。このプログラムを格納する記憶媒体としては、ＲＯＭ、ＦＤ、ＣＤ−ＲＯＭ、ＨＤ、メモリーカード、光磁気ディスクなどを用いることができる。
【００５２】
図２は図１で示したＸ線を１部透過するフィルタ３を模式的に示したものであり、アルミニュウムもしくは銅もしくは樹脂などで構成され、ハッチングで示した部分５２とそうでない部分５１とがあり、これらの部分５１、５２は異なる厚みもしくは物質で構成されている。このフィルタ３を透過したＸ線放射分布はこれらの部分５１、５２に合わせた強度分布となる。この２つの部分５１、５２の境界部分での変動は急峻である。
次に動作について図３のフローチャートのステップＳ１〜Ｓ８と共に説明する。
【００５３】
Ｓ１：Ｘ線発生装置１を動作させず、撮像素子４からの画像データをディジタル値に変換した後フレームメモリ８に記憶すると共に、暗出力ｄ（ｘ，ｙ）としてフレームメモリ１０に記憶する。
Ｓ２：被写体２が無い状態でＸ線発生装置１を動作させ、撮像素子４からの画像データをディジタル値に変換した後、フレームメモリ８に記憶し、Ｓ１で取得したｄ（ｘ，ｙ）をそれぞれの画素値から差し引き、その値をキャリアｃ（ｘ，ｙ）としてフレームメモリ９に記憶する。
Ｓ３：被写体２をフィルタ３と撮像素子４との間に入れ、Ｘ線発生装置１を動作させ、撮像素子４からの画像データををディジタル値に変換した後フレームメモリ８に記憶し、上述の暗出力ｄ（ｘ，ｙ）を差し引いて、再びフレームメモリ８にａ（ｘ，ｙ）として記憶する。
【００５４】
Ｓ４：フレームメモリ９のｃ（ｘ，ｙ）の画像から例えば図２の５３に示すような、強度分布が急峻に変動する部分｛Ｘ_0i≦ｘ≦Ｘ_1i，Ｙ_0i≦ｙ≦Ｙ_1i｜０≦ｉ≦ｋ−１｝を抽出、もしくは予め計算によって求めておく。図２のように縦方向に一定の分布であるならば、縦方向には考えず、横方向のみの一次元の演算でも構わない。
また上記部分５３に対応する画像としてフレームメモリ９の｛ｃ_i （ｘ，ｙ）｜Ｘ_0i≦ｘ≦Ｘ_1i，Ｙ_0i≦ｙ≦Ｙ_1i｜０≦ｉ≦ｋ−１｝と、フレームメモリ８の｛ａ_i （ｘ，ｙ）｜Ｘ_0i≦ｘ≦Ｘ_1i，Ｙ_0i≦ｙ≦Ｙ_1i｜０≦ｉ≦ｋ−１｝を抽出する。
Ｓ５：それらの部分分布に対して、（６）式又は（８）式を用いて各部分に対応する散乱値ｓ_i （０＜ｉ＜ｋ−１）を求め、各部分の中央における散乱値として対応するフレームメモリ１１中のアドレスに記憶する。
【００５５】
Ｓ６：Ｓ５で記憶されなかった位置のデータを各散乱値の補間によって計算し、全体の散乱分布ｓ（ｘ，ｙ）を求め、フレームメモリ１１に記憶する。
Ｓ７：後述。
Ｓ８：上記Ｓ５の演算を実行しても急峻な分布の変動影響は例えば図４に示すような形で残る。図４は部分分布を１次元で模式的に示したもので、６１がキャリアｃ（ｘ，ｙ）の変動を表し、６２が画像ａ（ｘ，ｙ）を表す。そして６３が
ｂ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）………（１２）
である。
【００５６】
６４で示す点が急峻に変動する点であるが、６１に比して６２の波形の方がボケているため、除算結果のｂ（ｘ，ｙ）には、図で示すような急峻に変動する位置に影響された微分波形に振動が発生している。この振動は、はやり散乱によって発生したものであるが、キャリアであるＸ線放射分布が変動しなければ、発生しないものであり、被写体の情報をボカしたものではない。そこで、ステップＳ７においては、上記取得したキャリアを関数でボカして分散値が最小になる関数を探索する。
【００５７】
図５は上記ステップＳ７を詳細に示すフローチャートである。ある部分にｉにおいて、
Ｓ７１：ε_-1＝０
Ｓ７２：ｐ＝０とする
Ｓ７３：ボケ関数をスパン２ｐ＋１のバイリニア関数で近似して選択（本実施の形態では、キャリアの変動方向は一次元のものであるので、ボケ関数も一次元で考えられる。）
【００５８】
【数１０】

【００５９】
Ｓ７４：畳み込み積分
【００６０】
【数１１】

【００６１】
Ｓ７５：除算演算および分散値
【００６２】
【数１２】

【００６３】
Ｓ７６：もしε_p ＞ε_p-1 ならば、Ｓ７７へそうでなければＳ７８へ。
Ｓ７７：Ｐ←Ｐ＋１としてＳ７３へ。
Ｓ７８：Ｐ←Ｐ−１とする。
Ｓ７９：ｈ_p （ｘ）を再び求め、ｄ_p （ｘ，ｙ）を計算し、フレームメモリ８の値の内容ａ（ｘ，ｙ）と部分分布において計算された内容ｄ_p （ｘ，ｙ）とのフレームメモリ１１内の内容ｓ（ｘ，ｙ）とを用いてフレームメモリ９の対応する部分の内容ｃ_i （ｘ，ｙ）を書き換え、ｂ（ｘ，ｙ）とする。
ｃ_i （ｘ，ｙ）←ｄ_p （ｘ，ｙ）
これらをすべての部分について行う。この時、ｈ_p （ｘ）を別な関数、例えばガウス分布などを用いても同様な結果が得られる。
Ｓ８０：今まで求められたデータｓ（ｘ，ｙ）、ｃ（ｘ，ｙ）を用いてフレームメモリ８の内容ａ（ｘ，ｙ）を下記の演算によってｂ（ｘ，ｙ）に書きかえる。
ｂ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）………（１６）
【００６４】
上記の方法により、散乱成分を除去し、さらにキャリアの変動に影響を取り除いた画像ｂ（ｘ，ｙ）をフレームメモリ８に記憶できる。
尚、上記過程においてＳ６で求められたｓ（ｘ，ｙ）を用いて、急峻な変動による影響が無視できるような場合、それ以降の過程を省略し、
｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ） ………（１７）
をそのままフレームメモリ８に記憶させて終了することも可能である。
【００６５】
次に本発明の第２の実施の形態を図６のフローチャートと共に説明する。
本実施の形態は、図１の構成において、第１の実施の形態の場合と逆の演算順次で行うようにしたものである。即ち、図６のＳ１１〜Ｓ１４は図３のＳ１〜Ｓと同じ処理が行われるが、次のステップＳ１５で図３のＳ７に対応する処理を先に行うようにしたものである。ステップＳ１５においては、図５のフローチャートに沿って行われ、Ｓ７１〜Ｓ７８は前述と同様に行われる。
【００６６】
次にｈ_p （ｘ）を再び求め、ｄ_p （ｘ，ｙ）を計算する。そしてフレームメモリ８の内容ａ（ｘ，ｙ）と部分分布において計算された内容ｄ_p （ｘ，ｙ）でフレームメモリ９の部分分布の内容ｃ_i （ｘ，ｙ）を書き換える。
ｃ_i （ｘ，ｙ）←ｄ_p （ｘ，ｙ） ………（１８）
これらをすべての部分について行う。
【００６７】
次に、図６のＳ１６で、それらの部分分布に対して、（６）式もしくは（８）式を用いて各部分に対応する散乱値ｓ_i （０＜ｉ＜ｋ−１）を求め、各部分の中央における散乱値として対応するフレームメモリ１１中のアドレスに記憶する。Ｓ１７では、Ｓ１６で記憶されなかった位置のデータを各散乱値の補間によって計算し、全体の散乱分布ｓ（ｘ，ｙ）を求め、フレームメモリ１１に記憶する。
そしてＳ１８でフレームメモリ８の内容ａ（ｘ，ｙ）とフレームメモリ９の内容ｃ（ｘ，ｙ）とフレームメモリ１１の内容ｓ（ｘ，ｙ）とを用いてフレームメモリ８の内容ａ（ｘ，ｙ）を書き換え、ｂ（ｘ，ｙ）とする。
ｂ（ｘ，ｙ）←｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）………（１９）
【００６８】
上記の方法により散乱成分を除去し、さらにキャリアの変動に影響を取り除いた画像ｂ（ｘ，ｙ）をフレームメモリ８に記憶できる。
【００６９】
次に第３の実施の形態を説明する。
本実施の形態は、上述の第１又は第２の実施の形態において散乱分布ｓ（ｘ，ｙ）とボケ関数ｈ（ｘ，ｙ）を求める分散値が最小になるまで収束させるようにしたものである。これによってさらに精度を上げることができる。
【００７０】
次に第４の実施の形態を説明する。
前述までの各実施の形態では、キャリアのボケを仮定した関数での畳み込み積分で近似したが本実施の形態では、散乱線成分を抽出した後、フィルタリングで画像を補修するようにしている構成は図１と同じである。
【００７１】
フローチャートは省略するが、図３のＳ１〜Ｓ６が実質的に同様に行われる。次のステップでは、フレームメモリ８の内容ａ（ｘ，ｙ）を下記の演算でｂ（ｘ，ｙ）に書き換える。
ｂ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）………（２０）
次にｂ（ｘ，ｙ）において、ｃ（ｘ，ｙ）が急峻に変動する部分に対応する影響がなくなるように空間フィルタリングをかける。
【００７２】
次に第５の実施の形態を説明する。
本実施の形態は、放射線分布を直線的に傾斜変動させることにより、撮影位置の変動にも影響されにくく、さらに急峻な変動もないようにして、それらの影響を取り除く過程を省略したことを特徴とする。
【００７３】
構成は図１と同様なものを用いるが、被写体とＸ線発生装置１との間に入れるフィルタ３を均一な物質で作り、その２次元的な厚みの関数ｔ（ｘ，ｙ）を
ｔ（ｘ，ｙ）＝ｋ３−（１／μ）ｌｎ（ｋ１・（ｘ＋ｙ）＋ｋ２）………（２１）
但し、ｋ１、ｋ２、ｋ３は定数、μは物質の放射線透過係数
となるように構成することにより、放射線分布が直線的な傾斜変動になるようにする。そして、前述と同様に図７のフローチャートに示す過程で計算する。この場合の放射線分布の概略図を図８に示す。
【００７４】
図７においてＳ２１〜Ｓ２３は図３のＳ１〜Ｓ３と同様である。
Ｓ２４：フレームメモリ８、９全体を小画像部分（５×５など）に分割する。Ｓ２５：それら部分画像に対して、前記式（１１）を用いて各部分に対応する散乱値ｓを求め、各部分の中央における散乱値として対応するフレームメモリ１１中のアドレスに記憶する。
Ｓ２６：Ｓ２５で記憶されなかった位置のデータを各散乱値の補間によって計算し、全体の散乱分布ｓ（ｘ，ｙ）を求め、フレームメモリ１１に記憶する。
Ｓ２７：フレームメモリ８の内容ａ（ｘ，ｙ）を下記の演算でｂ（ｘ，ｙ）に書き換える。
ｂ（ｘ，ｙ）＝｛ａ（ｘ，ｙ）−ｓ（ｘ，ｙ）｝／ｃ（ｘ，ｙ）
以上の過程により、より簡便に散乱線分布を除去した画像がフレームメモリ８に得られる。
【００７５】
本実施の形態では、撮影位置の移動にたいして、計算誤差が増えないように直線的傾斜した放射線分布を用いたが、キャリアのボケがあまり起こらない変動ならどのような変動でもかまわず、散乱線成分を仮定してその変動を最小にするものを求める（例えば正弦波てきな変動など）。
【００７６】
尚、本発明のすべての各実施の形態において、前述した各プログラムを実装したメモリ７としての記憶媒体７（図１）を他の装置に用いても同様の動作が可能である。
【００７７】
【発明の効果】
以上説明したように、本発明によれば被写体を透過した放射線分布を画像として取得する場合において、放射線源の分布を一様とせずに特徴を持たせた分布とし、取得画像からその特徴をなくすように、適当な分布の減算および前置的に取得した放射線分布を加工することによって、散乱線による影響をなくし、画質の良好な被写体の画像を取得することができる。
【図面の簡単な説明】
【図１】本発明の実施の形態を示すブロック図である。
【図２】本発明で用いるＸ線フィルタの一例を示す構成図である。
【図３】第１の実施の形態の動作を示すフローチャートである。
【図４】急峻な変動が残留することを説明するための構成図である。
【図５】ボケ関数を推定する部分のフローチャートである。
【図６】第２の実施の形態の動作を示すフローチャートである。
【図７】第５の実施の形態の動作を示すフローチャートである。
【図８】第５の実施の形態で用いる放射線分布を示す構成図である。
【図９】従来の放射線撮影システムを摸式的に示す構成図である。
【符号の説明】
１ Ｘ線発生装置
３ 放射線分布を一様にしないためのフィルタ
４ 放射線分布を電気信号に変換する撮像素子
６ 中央処理装置
７ プログラムを記憶する記憶媒体
８ 画像を記憶するフレームメモリ
９ 放射線分布を記憶するフレームメモリ
１０ 暗出力値を記憶するフレームメモリ
１１ 散乱分布を記憶するフレームメモリ

Claims (1)

1. 異なる厚みもしくは物質からなる少なくとも２つの部分で構成されたフィルタを透過した放射線が被写体を透過した場合に、前記被写体の散乱線分布の画像を、座標を（ｘ、ｙ）とする画像ｓ（ｘ，ｙ）として推定する放射線画像取得装置において、
   前記フィルタを透過した放射線の放射線分布を画像にする撮像手段と、
   前記フィルタと前記撮像手段との間に前記被写体がない状態で前記撮像手段において撮像された画像に基づいて、画像ｃ（ｘ，ｙ）を取得し、
   前記フィルタと前記撮像手段との間に前記被写体がある状態で前記撮像手段において撮像された画像に基づいて、画像ａ（ｘ，ｙ）を取得し、
   前記画像ｃ（ｘ，ｙ）及び前記画像ａ（ｘ，ｙ）において、前記フィルタの前記２つの部分の境界部分に対応する部分分布を、それぞれ、
   ｛ｃ_i′（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝、
   及び、
   ｛ａ_i′（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝
   （ただし、Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1iは画像上で前記２つの部分における境界部分を含む座標であり、ｋは部分分布の個数）として抽出し、
   前記ａ_i′（ｘ_i，ｙ_i）に対応する前記散乱線分布の部分をｓ_i′（ｘ_i，ｙ_i）とし、
   ｆ_i′（ｘ_i，ｙ_i）＝｛ａ_i′（ｘ_i，ｙ_i）−ｓ_i′（ｘ_i，ｙ_i）｝／ｃ_i′（ｘ_i，ｙ_i）で表せる式において、｛（ｘ_i，ｙ_i）｜Ｘ_0i≦ｘ_i≦Ｘ_1i，Ｙ_0i≦ｙ_i≦Ｙ_1i，０≦ｉ≦ｋ−１｝の範囲で、
   前記ｆ_i′（ｘ_i，ｙ_i）の合計値が最小となる前記ｓ_i′（ｘｉ，ｙｉ）の分布関数を、それぞれの部分（０≦ｉ≦ｋ−１）で推定し、
   当該推定により得られた分布関数に応じて散乱値ｓ′ｍ（０≦ｍ≦ｋ）の値を求め、求めた前記散乱値ｓ′ｍの値に基づいて、補間により前記画像ｓ（ｘ，ｙ）を推定する処理を行う処理手段と
   を有することを特徴とする放射線画像取得装置。

Images