JP3423828B2 - X線画像作成方法およびその装置 - Google Patents

X線画像作成方法およびその装置

Info

Publication number
JP3423828B2
JP3423828B2 JP31184195A JP31184195A JP3423828B2 JP 3423828 B2 JP3423828 B2 JP 3423828B2 JP 31184195 A JP31184195 A JP 31184195A JP 31184195 A JP31184195 A JP 31184195A JP 3423828 B2 JP3423828 B2 JP 3423828B2
Authority
JP
Japan
Prior art keywords
image
diffused light
ray
component
scattered
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP31184195A
Other languages
English (en)
Other versions
JPH09149895A (ja
Inventor
理香 馬場
健 植田
広則 植木
啓二 梅谷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hitachi Healthcare Manufacturing Ltd
Original Assignee
Hitachi Medical Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hitachi Medical Corp filed Critical Hitachi Medical Corp
Priority to JP31184195A priority Critical patent/JP3423828B2/ja
Priority to US08/759,088 priority patent/US5960058A/en
Publication of JPH09149895A publication Critical patent/JPH09149895A/ja
Application granted granted Critical
Publication of JP3423828B2 publication Critical patent/JP3423828B2/ja
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/58Testing, adjusting or calibrating thereof
    • A61B6/582Calibration
    • A61B6/583Calibration using calibration phantoms
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10TECHNICAL SUBJECTS COVERED BY FORMER USPC
    • Y10STECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10S378/00X-ray or gamma ray systems or devices
    • Y10S378/901Computer tomography program or processor

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Measurement Of Radiation (AREA)

Description

【発明の詳細な説明】
【0001】
【発明の属する技術分野】本発明は、X線画像作成方法
およびその装置に係り、それらによって得られるX線透
視画像あるいはX線撮影画像のぼけを補正するX線画像
作成方法およびその装置に関する。
【0002】
【従来の技術】X線透視撮影装置を用いてそのまま得ら
れるX線透視画像あるいはX線撮影画像は通常ぼけが発
生しており、完全なX線透視画像あるいはX線撮影画像
を得るためには、このぼけをなくす補正処理を施す必要
がある。
【0003】ここで、このぼけの発生の原因について以
下説明する。すなわち、X線画像の撮影又は透視を行う
方法として、被写体を透過したX線の強度分布をX線画
像センサにより光学像に変換し、その光学像をフィルム
に記録する方法、あるいは、その光学像をテレビカメラ
で電気的に読み出し、表示及びデジタル記録する方法等
が広く用いられている。
【0004】たとえばX線画像の撮影方法としては、X
線増感紙によりX線の強度分布を光学像に変換し、これ
に密着して配置したX線フィルムに直接的に画像を記録
する直接撮影法がある。また、X線蛍光板によりX線の
強度分布を光学像に変換し、この光学像をレンズやミラ
ーから構成される光学系を用いて結像し、結像面にフィ
ルムを配置して記録する間接撮影法がある。さらに、X
線蛍光板やX線イメージインテンシファイアを用いてX
線の強度分布を光学像に変換し、この光学像をレンズ等
の光学系を用いて撮像管やCCD等の撮像素子の撮像面
上に結像し、これを電気的に読み出す方法がある。
【0005】また、X線画像のリアルタイムの透視方法
として、X線イメージインテンシファイアやX線蛍光板
等を用いてX線の強度分布を光学像に変換し、この光学
像をレンズ等の光学系を用いて撮像管やCCD等の撮像
素子の撮像面上に結像し、これを電気的に読み出す方法
が広く用いられている。
【0006】これらの方法に用いられる検出系の代表例
を図2に示す。同図において、X線管(101)から照
射されたX線のうち被検体(110)と散乱X線遮蔽グ
リッド(102)を透過した直接X線(111)による
X線の強度分布をX線イメージインテンシファイア(1
03)を用いて電子ビーム(109)による像を経てX
線イメージインテンシファイア(103)内の出力蛍光
面(104)において光学像に変換し、この光学像をレ
ンズ及びミラー(105)を用いて撮像素子(106)
の撮像面上に結像し、これを電気的に読み出すようにな
っている。
【0007】ここで、前記出力蛍光面(104)では、
直接光(107)以外に光の拡散現象によりベーリング
グレアとも呼ばれる拡散光(108)が発生し、本来の
X線強度分布による像にぼけが加わる。一般に、X線の
強度分布を光学像に変換した後に、光学像を光学的に計
測する方法では、光学像が生成される媒体において、光
の拡散現象により本来のX線像にぼけが加わるという現
象が避けられない。拡散光によるぼけは、光学像が生成
される媒体においてX線像を縮小する割合が大きい場合
ほどその影響が大きくなる。X線イメージインテンシフ
ァイア(103)における拡散光は画像のコントラスト
を低下させる原因となっている。
【0008】そして、本来のX線画像に対してぼけ成分
となる原因としては、前記拡散光(108)のほかに、
被検体(110)の内部において、X線管から出射した
X線ビーム(以下、「直接X線」と略記する)の方向と
は異なる方向に散乱された散乱X線(112)がX線画
像センサで検出されるというX線の散乱現象がある。こ
のX線の散乱現象は避けることができず、散乱X線がX
線画像センサで検出されると、X線画像センサ面におけ
るX線の強度分布にぼけが加わることになる。
【0009】散乱X線がX線画像センサに入射すること
を制限する方法としては、被検体に照射するX線ビーム
を通常用いられる2次元のX線ビーム(円錐又は角錐状
のコーンビーム)ではなく、薄い扇状のファンビームと
し、1次元のX線画像を計測する方法がある。この方法
によれば、発生する散乱線の量が少なくなり、1次元の
X線画像を計測するセンサに入射する散乱線の量が減少
する。しかし、この方法では2次元画像を撮影するため
にはX線ビームを移動走査する必要があるので、高速撮
影や透視は不可能となる。
【0010】このように、通常の方法ではX線の強度分
布を一旦光学像に変換した後に光学像を光学的に計測す
る方法では2次元のX線ビームを用いる。この場合に散
乱X線を遮蔽する方法としては散乱X線遮蔽グリッド
(102)をX線画像検出器の前面に設置する方法が広
く利用されている。しかし、グリッドに入射する散乱X
線の全てを遮蔽することは不可能であり、その結果、計
測された画像は散乱X線によるぼけを生じる。散乱X線
によりX線画像検出器面上でぼけた画像は、さらに前述
の光学像生成媒体における拡散光(108)によるぼけ
が加わることになる。
【0011】以上、説明したように、X線像を光学像に
変換し、その画像を光学的に計測する方法では、被検体
の内部で発生する散乱X線によるX線画像センサ面にお
けるぼけと、光学像生成媒体における拡散光によるぼけ
との、2つの過程により本来のX線像に対するぼけが加
わり、画質を低下させるという問題がある。
【0012】そして、これらの拡散光によるぼけと散乱
X線(以下、「散乱線」とも略記する)によるぼけを補
正する方法がいくつか知られている。以下、補正方法及
び関連従来技術について説明する。
【0013】(1)まず、X線イメージインテンシファ
イアの拡散光によるぼけを直接デコンボリューション法
により補正する技術が、たとえばSeibertらの文
献Medical Physics誌、12巻、281
−288頁(1985年)に開示されている。また、拡
散光のぼけを補正した後に引き続いて散乱X線のぼけを
直接デコンボリューション法により補正する技術が、た
とえばSeibertらの別の文献Medical P
hysics誌、15巻、567−575頁(1988
年)に開示されている。
【0014】ここで、図3に、直接デコンボリューショ
ン法により補正する方法における、画像劣化モデル(同
図(a))と、補正演算(同図(b))とを示す。この
方法の画像劣化モデルでは、散乱線による画像のぼけ過
程を、散乱線を含んだ点像分布関数202と原画像(X
線の強度分布)201とのコンボリューション203に
より表現する。さらに、拡散光によるぼけ過程を、拡散
光を含んだ点像分布関数205と散乱線によるぼけ過程
を経た画像204とのコンボリューション206により
表現する。なお、同図では、上述した点像分布関数をP
SF(ポイントスプレッドファンクションの略)と略記
している。
【0015】これに対応する補正演算は、点像分布関数
の拡散光成分を指数関数で表現し、これから導出された
拡散光に関するデコンボリューションフィルタ(20
8)を用いたデコンボリューション演算(209)によ
り、拡散光によるぼけを回復した拡散光補正画像(21
0)を生成する。さらに、点像分布関数の散乱線成分を
正規分布関数で表現し、これから導出された散乱線に関
するデコンボリューションフィルタ(211)を用いた
デコンボリューション演算(212)により、拡散光と
散乱線によるぼけを回復した補正画像(213)を生成
する。
【0016】なお、これらの補正演算におけるデコンボ
リューションフィルタとして、点像分布関数の2次元フ
ーリエ変換の逆数として周波数空間におけるデコンボリ
ューションフィルタを求める。デコンボリューション演
算は、対象とする画像の2次元フーリエ変換に対してデ
コンボリューションフィルタを乗算し、その結果を2次
元フーリエ逆変換して実行される。ここで、拡散光に関
するデコンボリューションフィルタ(208)は拡散光
の直接光に対する強度比(以下、単に「拡散光強度比」
と略記する)をパラメータとして含む。また、散乱線に
関するデコンボリューションフィルタ(211)は散乱
線の直接線に対する強度比(以下、単に「散乱線強度
比」と略記する)をパラメータとして含む。拡散光強度
比は予め、イメージインテンシファイアについて求めて
おく。散乱線強度比は、計測画像の値から、予め計測条
件毎に準備したルックアップテーブルを索表することに
より推定される。ルックアップテーブルは、計測時のX
線管電圧、被検体の厚さ、散乱X線遮蔽用グリッド、視
野径、幾何学的配置毎に予め生成される。
【0017】(2)また、ぼけ画像生成法により補正す
る技術は、Molloiらにより、Medical P
hysics誌、15巻、289−297頁(1988
年)に開示されている。また、Hondaらにより、M
edical Physics誌、20巻、59−69
頁(1993年)に開示されている。
【0018】図4に、ぼけ画像生成法により補正する方
法において、画像劣化モデル(同図(a))と、補正演
算(同図(b))とを示す。この方法の画像劣化モデル
では、散乱線による画像のぼけ過程と拡散光によるぼけ
過程を一体化した過程とみなし、その過程を拡散光成分
と散乱線成分の和の直接線の直接光に対する強度比(以
下、単に「拡散光・散乱線成分の強度比」と略記する)
(302)と拡散光と散乱線の和の点像分布関数(以
下、単に「拡散光・散乱線の点像分布関数」と略記す
る)(303)との積(304)と原画像(301)と
のコンボリューション(305)により拡散光と散乱線
を一体化した成分からなるぼけ画像(以下、単に「拡散
光・散乱線成分画像」と略記することもある)(30
6)を生成し、これと原画像(301)との加算(30
7)により計測画像(308)が生成すると表現する。
なお、拡散光・散乱線の強度比と拡散光・散乱線の点像
分布関数との積(304)を「拡散光・散乱線の強度分
布関数」と略記することもある。
【0019】これに対応する補正モデルは、拡散光・散
乱線の強度比(302)とその点像分布関数(303)
とを求め、これらの積(304)と計測画像(308)
とのデコンボリューション演算(309)により拡散光
・散乱線成分画像(310)を生成する。ぼけ画像(3
10)を計測画像(308)から減算(311)して、
補正画像(312)を得る。
【0020】デコンボリューション演算(309)とし
て、コンボリューション法が上記Molloiらの文献
に、フーリエ変換法が上記Hondaらの文献にそれぞ
れ開示されている。コンボリューション法におけるデコ
ンボリューション演算(309)は、拡散光・散乱線の
強度比(302)にその点像分布関数(303)を乗じ
て拡散光・散乱線の強度分布関数(304)を求め、こ
れと計測画像(308)との2次元コンボリューション
演算を行なう。
【0021】フーリエ変換法におけるデコンボリューシ
ョン演算(309)は、拡散光・散乱線の強度比(30
2)とその点像分布関数(303)の2次元フーリエ変
換との積の関数として周波数空間におけるデコンボリュ
ーションフィルタを求め、計測画像(308)の2次元
フーリエ変換に対してデコンボリューションフィルタを
乗算し、その結果を2次元フーリエ逆変換して実行され
る。拡散光・散乱線の強度比は、計測画像の値、計測時
のX線管電圧、被検体の厚さ、散乱X線遮蔽グリッド、
視野径、被検体とグリッドの距離(エアギャップ)等の
関数である。上記の拡散光・散乱線の強度比の決定法と
しては、計測画像の値から、予め計測条件毎に準備した
ルックアップテーブルを索表して求める方法が上記のM
olloiらの文献に開示されている。また、別の決定
法として、計測画像の値の最大値と上記計測条件から索
表と演算により算出する方法がHondaらの上記の文
献及びHondaらの別の文献Medical Phy
sics誌、18巻、219−226頁(1991年)
に開示されている。
【0022】
【発明が解決しようとする課題】しかしながら、上述し
た各X線画像作成方法において、まず、図3に示した直
接デコンボリューション法では、ぼけの回復のために計
測画像の2次元フーリエ変換に周波数空間で高周波強調
フィルタを乗算するので、ディジタル画像のサンプリン
グピッチで決定されるナイキスト周波数より高周波数の
画像成分による折り返しノイズ等の高周波ノイズが強調
されるという問題が残存されていた。
【0023】また、図4に示したぼけ画像生成法では、
計測画像の2次元フーリエ変換に周波数空間で高周波減
衰フィルタを乗算するので、周波数空間における高周波
ノイズの増加はないが、拡散光と散乱線を一体化した点
像分布関数を用いているため、これらの正確な補正はで
きないものとなっていた。
【0024】すなわち、拡散光の点像分布関数は視野サ
イズで決定され、X線管電圧、散乱線遮蔽グリッド、あ
るいは被検体の厚さには依存しないが、散乱線の点像分
布関数は拡散光の点像分布関数と異なり、視野サイズだ
けでなく、X線管電圧、散乱線遮蔽グリッド、あるいは
被検体の厚さに依存するようになる。したがって、一般
的に、被検体の厚さの変化範囲が大きな被検体に対して
は、上記のぼけ画像生成法の補正はきわめて粗い近似に
留まるを得ないという問題が残存されていた。
【0025】さらに、上述の従来技術では、撮像素子に
過大な強度の光が入射して撮像素子の出力が飽和レベル
に達し、その結果、飽和レベルに達した領域の周辺領域
において散乱線や拡散光の影響が相対的に過大になるハ
レーションへの対応が考慮されておらず、ハレーション
が発生した場合の補正は全く不十分な補正となっていた
という問題が残存されていた。
【0026】さらに、連続的に計測される透視画像に対
しても拡散光成分と散乱線成分とを個別に補正する高速
な補正処理は不可能であった。
【0027】本発明はこのような事情に基づいてなされ
たものあり、その第1の目的は、X線透視画像あるいは
X線撮影画像に対して、周波数空間で高周波強調フィル
タの乗算を行なうことによる高周波ノイズの増加なし
に、拡散光と散乱線に起因するX線画像のぼけの補正を
それぞれの点像分布関数を用いて正確に補正し、これに
より、厚さの変化範囲が大きな被検体に対しても、高精
度で補正することのできるX線画像作成方法とその装置
を提供することにある。
【0028】本発明の第2の目的は、上記の補正に必要
となる拡散光強度分布関数と散乱線強度分布関数とを、
簡便な方法により求めることのできるX線画像作成方法
とその装置を提供することにある。
【0029】本発明の第3の目的は、撮像素子の出力が
飽和してハレーションが発生した場合にも、上記の補正
を近似的に行なうことのできるX線画像作成方法とその
装置を提供することにある。
【0030】本発明の第4の目的は、上記の補正を、特
別な演算器なしのシステムでも高速で実行可能とし、特
に高速処理が要求される透視像や、例えば2000×2
000画素や4000×4000画素の超高精細画像に
対しても、また、例えば512×512画素の画像を2
40枚撮影するコーンビームCT等にも適用することの
できるX線画像作成方法とその装置を提供することにあ
る。
【0031】本発明の第5の目的は、補正演算に要する
時間が透視像の1フレームの時間より長い場合にも連続
的に変化する透視像に対する近似的な補正をすることの
できるX線画像作成方法とその装置を提供することにあ
る。
【0032】
【課題を解決するための手段】本発明の第1の目的を達
成するためには次の手段によって達成される。
【0033】すなわち、予め求めた拡散光の強度比と予
め求めた拡散光の点像分布関数との積により拡散光の強
度分布関数を求め、拡散光の強度分布関数と計測画像と
のデコンボリューション演算によって拡散光成分画像を
算出し、一方、この拡散光の強度分布関数と、別に計測
画像とその撮影条件から求めた散乱線の強度比と予め求
めた散乱線の点像分布関数との積により求めた散乱線の
強度分布関数と、計測画像とのデコンボリューション演
算により散乱線成分画像を算出し、計測画像から拡散光
成分画像及び散乱線成分画像を減算することによって補
正画像を求めることができる。
【0034】以下、この概要を図1を用いて説明する。
同図において、画像劣化モデル(同図(a))と、補正
演算(同図(b))を示す。同図に示す画像劣化モデル
では、散乱線の強度比as(402)と散乱線の点像分
布関数PSFs(403)との積(404)(以下、
「散乱線の強度分布関数」と略記する)と原画像I
p(401)とのコンボリューション(405)により
散乱線成分画像Is(406)が式1により生成され
る。以下、記号**は2次元コンボリューションを表わ
す。
【0035】
【数1】 Is=as・Ip**PSFs …(式1) 散乱線成分画像Isと原画像Ipとの加算(407)によ
り散乱線によるぼけ加算画像(408)が生成され、加
算画像(408)と拡散光の強度比av(409)と拡
散光の点像分布関数PSFv(410)との積(以下、
「拡散光の強度分布関数」と略記する)とのコンボリュ
ーション(412)により拡散光成分画像(413)が
式2により生成される。
【0036】
【数2】 Iv=av・(Ip+Is)**PSFv …(式2) IvとIpとの加算(414)により計測画像It(41
5)が式3により生成される。
【0037】
【数3】 It=Ip+Is+Iv …(式3) これに対応する補正演算として、予め求めた拡散光の強
度比(409)と予め求めた拡散光の点像分布関数(4
10)との積(411)により拡散光の強度分布関数を
求め、拡散光の強度分布関数と計測画像(415)との
デコンボリューション演算(417)によって拡散光成
分画像Iv(418)を算出する。
【0038】一方、この拡散光の強度分布関数と、別に
計測画像(415)とその撮影条件から求めた散乱線強
度比(402)と予め求めた散乱線の点像分布関数(4
03)との積(404)により散乱線の強度分布関数を
求め、この散乱線強度分布関数と計測画像(415)と
のデコンボリューション演算(419)により散乱線成
分画像Is(420)を算出し、計測画像(415)から
拡散光成分画像(418)及び散乱線成分画像(420)を
減算(421)して補正画像Ip(422)を求める。こ
の補正演算により、本発明の第1の目的が達成される。
【0039】以下、上記のデコンボリューション(41
7)及び(419)の方法を説明する。対応する補正演算
は、
【0040】
【数4】 Iv=av・(It−Iv)**PSFv …(式4)
【0041】
【数5】 Is=as・(It−Iv−Is)**PSFs …(式5)
【0042】
【数6】 Ip=It−Iv−Is …(式6) となる。2次元フーリエ変換をF2[ ]で表現する
と、式4の2次元フーリエ変換から、
【0043】
【数7】 F2[Iv]=F2[It]・{(av・F2[PSFv])/(1+av・F2[ PSFv])} …(式7) となる。式5の2次元フーリエ変換と式7から、
【0044】
【数8】 F2[Is]=F2[It]・{(as・F2[PSFs])/(1+ as・F2[PSFs])}・{1/(1+av・F2[PSFv])} …(式8) となる。従って、式7及び式8の2次元フーリエ逆変換
からIv及びIsが求まる。即ち、2次元フーリエ逆変換
をFR2{ }で表現すると、
【0045】
【数9】 Iv=FR2{F2[It]・{(av・F2[PSFv])/(1+av・F2[ PSFv])}} …(式9)
【0046】
【数10】 Is=FR2{F2[It]・{(as・F2[PSFs])/(1+ as・F2[PSFs])}・{1/(1+av・F2[PSFv])} …(式10) となる。従って、図4におけるデコンボリューション演
算(417)及び(419)はそれぞれ式9と式10と
すればよい。式6に式9と式10を代入すれば、式11
によりIpが求まり補正画像を得ることができるように
なる。
【0047】
【数11】 Ip=It−FR2{F2[It]・{(av・F2[PSFv])/ (1+av・F2[PSFv])}}−FR2{F2[It]・{(as・ F2[PSFs])/(1+as・F2[PSFs])}・{1/(1+av・ F2[PSFv])} …(式11) 周波数空間におけるフィルタは式7の右辺の中括弧内、
及び式8の右辺の中括弧内であり、いずれも高周波域で
はゼロに漸近する。従って、これらのフィルタは高周波
強調フィルタではなく高周波低下フィルタであり、周波
数空間における演算により高周波ノイズが強調されるこ
とはない。
【0048】次に、本発明の第2の目的を達成させるた
めには次のような手段によって達成される。
【0049】計測画像のマトリックスの方向、即ち縦方
向、及び横方向の各々に対し、被写体がないとき及び被
写体として均一な厚さの人体模擬被写体を配置したとき
のエッジ強度分布を計測し、これを微分して、被写体が
ないときの線像強度分布を求め、この結果に対して2成
分の関数フィッティングを行なう。この2成分のうち幅
広い分布を拡散光成分とし、その線像分布関数と拡散光
成分の強度比を求め、縦横2方向の線像分布関数の積を
拡散光成分の点像分布関数とし、縦横2方向の相対強度
の平均を拡散光成分の強度比とし、拡散光成分の点像分
布関数と拡散光成分の強度比の積を拡散光強度分布とす
る。
【0050】さらに、均一な厚さの人体模擬被写体、例
えばアクリル板、を配置したときのエッジ強度分布に対
して、拡散光の線像分布関数と拡散光成分の強度比との
積を1次元コンボリューション演算又は1次元フーリエ
変換を利用して1次元拡散光像を求め、均一な厚さの人
体模擬被写体を配置したときのエッジ強度分布から1次
元拡散光像を減算してエッジ強度分布の拡散光補正を行
なう。これを微分して、対応する被写体厚さに対する拡
散光補正線像強度分布を求め、被写体厚さに対する補正
線像強度分布に対して2成分の関数フィッティングを行
なう。2成分のうち幅広い分布を散乱線成分とし、その
線像分布関数と強度比を求め、縦横2方向の線像分布関
数の積を散乱線成分の点像分布関数とし、縦横2方向の
強度比の平均を散乱線の強度比とし、散乱線成分の点像
分布関数と散乱線成分の強度比の積を散乱線強度分布と
する。これにより目的が達成できるようになる。
【0051】また、本発明の第3の目的を達成させるた
めには次のような手段によって達成される。
【0052】図5に示すように、計測画像(X線透視画
像あるいはX線撮影画像)(ステップ501)が飽和値
を含むか否かの判別動作(ステップ502)を行なう。
飽和値を含む場合には、飽和した画素の値を別に計測さ
れた同一の位置近傍の画像の値から推定した値(ステッ
プ503)に変更(ステップ504)して飽和レベルの
補正を行ない、飽和レベルの補正後に拡散光と散乱線の
両方を補正する(ステップ505)。これにより、撮像
素子の出力が飽和してハレーションが発生した場合にも
拡散光と散乱線の両方を近似的に補正できるようにな
る。
【0053】さらに、本発明の第4の目的を達成させる
ためには次のような手段によって達成される。
【0054】すなわち、計測画像と、拡散光強度分布又
は散乱線強度分布をそれぞれ離散的サンプリングして、
それぞれのデータのマトリックスサイズを縮小した、サ
ンプリング計測画像と、サンプリング拡散光強度分布又
はサンプリング散乱線強度分布とを求め、サンプリング
計測画像とサンプリング拡散光強度分布とのコンボリュ
ーション演算、又は、サンプリング計測画像とサンプリ
ング散乱線強度分布とのコンボリューション演算を実行
し、コンボリューションの実行結果を補間法を利用して
画像拡大処理してもとのマトリックスサイズに戻し、も
とのマトリックスサイズに戻された画像を計測画像から
減算して実現できるようになる。
【0055】さらに、本発明の第5の目的を達成させる
ためには次のような手段によって達成される。
【0056】すなわち、連続的に収集された透視画像に
対し、透視画像を間歇的にサンプリングし、サンプリン
グ画像に対して、引き続いて、拡散光成分画像、及び散
乱線成分画像を算出し、最後にこれらの成分を透視画像
から減算するに際して、対応する透視画像ではなく次に
表示される最新の透視画像から減算して実現できるよう
になる。
【0057】
【発明の実施の形態】実施例1. 図6は、本発明によるX線画像作成方法およ
び装置の一実施例を示すブロック構成図である。
【0058】同図において、まず、X線透視像あるいは
X線撮影像からなる計測画像It(601)と、拡散光
点像分布関数PSFv(604)と、散乱線点像分布関
数PSFs(605)の2次元フーリエ変換(60
6)、(607)、(608)を求める。
【0059】そして、拡散光強度比av(603)と拡
散光点像分布関数(604)の2次元フーリエ変換(6
07)との積により拡散光強度分布関数の2次元フーリ
エ変換像(609)を求め、拡散光強度分布関数の2次
元フーリエ変換像(609)を拡散光強度分布関数の2
次元フーリエ変換像に1を加算した式で除算した結果
(611)を計測画像の2次元フーリエ変換像(60
6)に乗じ、その積(612)を2次元フーリエ逆変換
(613)し、拡散光成分画像(614)を算出する。
この演算は、既出の式9によって算出できる。
【0060】散乱線強度比as(602)は前記計測画
像It及び計測条件から求める。その方法の具体例は後
に詳述する。そして、散乱線強度比as(602)と散
乱線点像分布関数(605)の2次元フーリエ変換(6
08)との積により散乱線強度分布関数の2次元フーリ
エ変換像(610)を求め、散乱線強度分布関数の2次
元フーリエ変換像を散乱線強度分布関数の2次元フーリ
エ変換像に1を加算した式で除算した結果(615)と
拡散光強度分布関数の2次元フーリエ変換像に1を加算
した結果の逆数(616)を乗算し、その結果を計測画
像Itの2次元フーリエ変換像(606)に乗じる。
【0061】さらに、その結果(617)を2次元フー
リエ逆変換(618)し、散乱線成分画像(619)を
算出する。この演算は既出の式10によって算出でき
る。拡散光成分画像(614)と散乱線成分画像(61
9)とを計測画像(601)から引き算(620)す
る。この演算は既出の式6によって算出できる。この演
算により、計測画像Itに含まれる拡散光成分と散乱線
成分の両方を補正した目的の画像(621)を得ること
ができる。
【0062】実施例2.図7は、本発明によるX線画像
作成方法および装置の他の実施例を示すブロック構成図
である。この実施例では、2次元フーリエ変換の代わり
に2次元コンボリューションを用いて補正するところが
実施例1と異なっている。
【0063】まず、その原理について説明する。
【0064】コンボリューション法による真の拡散光成
分画像Ivは、既出の式4において、その右辺のIvに同
一の式4を繰返し代入することによって式12となる。
【0065】
【数12】 Iv=av・It**PSFv−av 2・It**PSFv**PSFv +av 3・It**PSFv**PSFv**PSFv−… …(式12) この結果、拡散光補正画像は式13のようになる。
【0066】
【数13】 It−Iv=It−av・It**PSFv+av 2・It**PSFv**PSFv− av 3・It**PSFv**PSFv**PSFv+… …(式13) 真の散乱線画像は、式5において右辺のIsに同一の式
5を繰返し代入することによって式14のようにして表
わされる。
【0067】
【数14】 Is=as・(It−Iv)**PSFs−as 2・(It−Iv)**PSFs **PSFs+as 3・(It−Iv)**PSFs**PSFs**PSFs−… …(式14) そして、式14に式13を代入し、さらに、その式を展
開することによって式15を得る。
【0068】
【数15】 Is=as・It**PSFs−as・av・It**PSFv**PSFs +as・av 2・It**PSFv**PSFv**PSFs−as 2・It **PSFs**PSFs+as 2・av・It**PSFv**PSFs **PSFs+as 3・It**PSFs**PSFs**PSFs+… …(式15) したがって、補正画像Ipは、式6に式13と式15を
代入し、その後、その式を整理することによって式16
を得る。
【0069】
【数16】 Ip=It−av・It**PSFv−as・It**PSFs+av 2・It **PSFv**PSFv+as・av・It**PSFv**PSFs +as 2・It**PSFs**PSFs−av 3・It**PSFv**PSFv **PSFv−as・av 2・It**PSFv**PSFv**PSFs −as 2・av・It**PSFv**PSFs**PSFs−as 3・It **PSFs**PSFs**PSFs+… (式16) ここで、式16による補正を近似式によって行う方法を
以下説明する。すなわち、拡散光成分画像を求める正確
な演算は式12であるが、avが1に比較して小さいこ
とに着目し、近似的な拡散光成分画像Iv’を式12の
第1項により求め、これにより式17を得る。
【0070】
【数17】 Iv’=av・It**PSFv …(式17) そして、拡散光成分補正画像をIq’とすると、Iq’は
式18となる。
【0071】
【数18】 Iq’=It−Iv’=It−av・It**PSFv …(式18) 散乱線画像を求める正確な演算は式14又は式15であ
るが、散乱線遮蔽グリッドを用いた場合のasは1に比
較して一般に小さいことに着目し、近似的な散乱線成分
画像Is’を式14の第1項のIvをIv’で近似した演
算により求め、式19を得る。
【0072】
【数19】 Is’=as・(It−Iv’)**PSFs …(式19) そして、近似的な散乱線成分画像Is’は、式19に式
17を代入することによって、式20となる。
【0073】
【数20】 Is’=as・It**PSFs−as・av・It**PSFv**PSFs …(式20) さらに、近似的な補正画像P’を式21の演算により求
める。
【0074】
【数21】 Ip’=Iq’−Is’ …(式21) 補正画像P’は、式21に式18と式20を代入するこ
とによって、式22となる。
【0075】
【数22】 Ip’=It−av・It**PSFv−as・It**PSFs+as・av・It **PSFv**PSFs …(式22) このようにして求めた補正画像P’の補正誤差eは、式
22のIp’と式16のIpとの差であり、その差は次式
23で表わされる。
【0076】
【数23】 e=Ip’−Ip =−av 2・It**PSFv**PSFv−as2・It**PSFs**PSFs +av 3・It**PSFv**PSFv**PSFv+as・av 2・It*PSFv **PSFv**PSFs+as 2・av・It**PSFv**PSFs** PSFs+as 3・It**PSFs**PSFs**PSFs−… …(式23) 補正誤差はav及びasに関して2次以上の項であり、一
般にav及びasが1より小さいので近似補正として成立
する。
【0077】以上により、図1における拡散光成分画像
を求めるデコンボリューション演算(417)を式17
の2次元コンボリューション演算とし、散乱線成分画像
を求めるデコンボリューション演算(419)を式19
の2次元コンボリューション演算として、拡散光成分と
散乱線成分を補正できる。
【0078】なお、式23における主要な誤差項である
2次の項、即ち第1項及び第2項の項の符号は負である
から、補正はやや過補正となる。近似度の向上が必要の
場合には、上記第1項及び第2項を演算により求めて補
正画像に加算すればよい。
【0079】図7は、以上説明したステップを具体的に
示したフロー図である。
【0080】同図において、まず、計測画像It(70
1)の取得以前に、予め、拡散光強度比av(703)
と拡散光点像分布関数PSFv(704)との積により
拡散光強度分布関数(706)を求める。
【0081】散乱線強度比as(702)は計測画像及
び計測条件から求める。その方法の具体例は後に詳述す
る。散乱線強度比as(702)と散乱線点像分布関数
PSFs(705)との積により散乱線強度分布関数
(711)を求める。計測画像It(701)に対し
て、拡散光強度分布関数(706)をコンボリューショ
ン(707)して拡散光画像成分Iv’(708)を算
出する。この演算は既出の式17によって算出できる。
【0082】そして、拡散光画像成分Iv’(708)
を計測画像It(701)から減算(709)し、拡散
光補正画像Iq’(710)を算出する。この演算は既
出の式18によって算出できる。
【0083】この拡散光補正画像Iq’(710)に対
して散乱線強度分布as(711)をコンボリューショ
ン(712)して散乱線成分画像Is’(713)を算
出する。この演算は既出の式19によって算出できる。
【0084】散乱線成分画像Is’(713)を拡散光
補正画像Iq’(710)から減算(714)する。こ
の演算は既出の式21によって算出できる。これによ
り、計測画像Itに含まれる拡散光成分Iv’と散乱線成
分Is’の両方を補正した画像Ip’(715)を得るこ
とができる。
【0085】実施例3.次に、拡散光強度比avと拡散
光の点像分布関数PSFvを実測により求める方法を、
図8、図9を用いて以下説明する。
【0086】まず、図9(A)は、拡散光の強度比及び
拡散光のエッジ分布関数を計測する計測系を示す。金属
板(902)のエッジ(903)をグリッド面(90
4)から約25cm離して視野の中央部に、X軸に直交
するよう配置する。散乱体は配置せず、エッジの計測画
像Toxを計測する(801)。図9(B)はエッジの計
測画像の概形を示す。この画像は拡散光によるボケが加
わったものとして得られる。次にエッジの計測画像の微
分画像を求める(803)。図9(C)はエッジの微分
画像の概形を示す。同図においてGauss曲線は検出
器の分解能を、exponential曲線は光拡散の
分布を示している。そして、その微分画像を2成分の関
数フィッティング(805)して、拡散光強度比a
vx(806)と線像分布関数LSFvx(807)を求め
る。以下、この方法につき詳述する。
【0087】まず、図9(A)の配置でエッジの計測画
像Toxを計測する(801)。エッジの計測画像の概形
を図9(B)の曲線906に示す。エッジの計測画像T
oxは、散乱線成分を含まないことから式3により、直接
X線成分画像Poと拡散光成分画像Voの和である式24
となる。そして、直接X線成分画像Poは式25により
表わされる。
【0088】
【数24】 Tox=Po+Vo …(式24)
【0089】
【数25】 Po=(apo・STEPx)**PSFp …(式25) 但し、STEPxはx方向の高さ1の階段関数、apo
計測時の直接X線の強度に比例する項、PSFpは直接
X線画像の点像分布関数を表している。
【0090】点像分布関数PSFvとX方向線像分布関
数LSFvx及びY方向線像分布関数LSFvyとの関係は
式26を仮定する。
【0091】
【数26】 PSFv(x,y)=LSFvx(x)・LSFvy(y) …(式26) 以下の式の導出には、線像分布関数とエッジ分布関数の
関係を表わす式27を利用する。但し、記号*は1次元
コンボリューションを表わす。
【0092】
【数27】 STEPx*LSFx=ESFx …(式27) ここで、式25に式26と式27を代入すると、直接線
成分画像Poは、次式28となる。
【0093】
【数28】 Po=apo・ESFpx …(式28) そして、拡散光成分画像Voは、計測される拡散光成分
の強度比をavxとすると、次式29となる。
【0094】
【数29】 Vo=(apo・STEPx)**(avx・PSFv) …(式29) ESFvxを拡散光成分のエッジ分布関数として、式29
に式26と式27を代入して、これにより次式30を得
る。
【0095】
【数30】 Vo=avx・apo・ESFvx …(式30) したがって、エッジの計測画像Toxは、式24に式28
と式30を代入することにより、yに依存しないxのみ
の関数となり、2成分のエッジ分布関数の1次結合によ
り表現され(802)、次式31となる。
【0096】
【数31】 To=apo・(ESFpx+avx・ESFvx) …(式31) 次にエッジの計測画像の微分画像を求める(803)。
微分画像の概形を図9(C)に示す。微分画像d
(Tox)/dxは式31により、次式32となる。
【0097】
【数32】 d(Tox)/dx=apo・(LSFpx+avx・LSFvx) …(式32) 但し、LSFvxは次式33で表わされる。
【0098】
【数33】 LSFvx≡d(ESFvx)/dx …(式33) すなわち、エッジの計測画像の微分画像d(Tox)/d
xは2成分の線像分布関数の1次結合により表現される
(804)。次に、この微分画像d(Tox)/dxにつ
いて、直接線成分の線像分布関数LSFpxを正規分布関
数、拡散光成分の線像分布関数LSFvxを指数関数と仮
定して2成分関数フィッティング(805)を行ない、
拡散光強度比avx(806)と拡散光成分の線像分布関
数LSFvx(807)を求める。微分画像は、例えば隣
同志の画素の値の差分から求められる。この場合、1ラ
インデータではX線量子による統計的なノイズが多いの
で、複数ラインのデータを加算したデータを用いる。ま
た、微分画像は左右対称形となるが、ノイズは金属板で
遮蔽されている側(図9(C)では左側)半分では小さ
く、遮蔽されていない側(図9(C)では右側)半分で
は大きい。そこで、遮蔽されている側のデータをピーク
を中心として折り返す。このようにして得られた曲線に
ついて2成分の関数フィッティングを行なう。関数フィ
ッティングでは、式32において、直接線成分の線像強
度分布を次式34と、また、拡散光成分の線像強度分布
を次式35と設定し、最小2乗法を利用したフィッティ
ングによりa1x、a2x、b1x、b2xを決定する。
【0099】
【数34】 apo・LSFpx=a1x・exp(−a2x・x2) …(式34)
【0100】
【数35】 apo・avx・LSFvx=b1x・exp(−b2x・|x|) …(式35) 拡散光強度比avxは、フィッティングにより決定された
係数を用いて、次式36で算出する。
【0101】
【数36】 avx=(2・b1x/b2x)/{a1x√(π/a2x)} …(式36) 拡散光成分の線像分布関数LSFsは、次式37で算出
する。
【0102】
【数37】 LSFvx=(b2x/2)・exp(−b2x・|x|) …(式37) Y方向についても金属板のエッジの方向をX方向につい
ての計測のときから90°回転して、視野の中央部にY
軸に直交するように配置して、X方向の場合と同様の計
測とフィッティングを行ない、avy(808)とLSF
vy(809)を、次式38、次式39により決定する。
【0103】
【数38】 avy=(2・b1y/b2y)/{a1y√(π/a2y)} …(式38)
【0104】
【数39】 LSFvy=(b2y/2)・exp(−b2y・|y|) …(式39) そして、拡散光の点像分布関数PSFvと拡散光強度比
vを、次式40、次式41で求める(810、81
1)。具体的には、次式42、次式43となる。
【0105】
【数40】 PSFv(x,y)=LSFvx(x)・LSFvy(y) …(式40)
【0106】
【数41】 av=(avx+avy)/2 …(式41)
【0107】
【数42】 PSFv(x,y)=(b2x・b2y/4)・exp(−b2x・|x| −b2y・|y|) …(式42)
【0108】
【数43】 av=(b1x/b2x)/{a1x√(π/a2x)}+(b1y/b2y)/ {a1y√(π/a2y)} …(式43) X方向とY方向の特性が近似的に等しいと仮定できる場
合には、式44となり、PSFv(x,y)は単一のパ
ラメータb2で決定される。但し、b2=b2x=b2 yであ
る。
【0109】
【数44】 PSFv(x,y)=(b2/2)2・exp(−b2・|x|−b2・|y|) …(式44) 次に、実験により一定の厚さの被検体に対する散乱線強
度比asと散乱線の点像分布関数PSFsを求める方法を
説明する。
【0110】図9(A’)は、散乱線の強度比及び散乱
線のエッジ分布関数を計測する計測系を示す。金属板9
02のエッジ903をグリッド面904から約25cm
離して視野の中央部に、X軸に直交するよう配置する。
散乱体としてたとえばアクリル板を配置し、エッジの計
測画像Tsxを計測する。図9(B’)はエッジの計測画
像の概形を示す。この画像は拡散光と散乱線によるボケ
が加わったものとして得られる。
【0111】次にこのエッジの計測画像に対して、先に
求めた拡散光の強度と線像分布関数を用いて拡散光成分
の補正を行なう。エッジの計測画像の1次元の画像プロ
ファイル(852)に対して、1次元のコンボリューシ
ョン演算(853)により1次元の拡散光成分画像プロ
ファイル(854)を求め、1次元のプロファイル画像
(852)から1次元の拡散光成分画像プロファイル
(854)を減算(855)して、拡散光補正画像の1
次元プロファイル画像(856)を求める。図9
(C’)は拡散光補正画像の1次元プロファイル画像
(856)の概形を示す。すなわち、この画像は散乱線
のみによるボケが加わったものとして得られる。
【0112】次に、この拡散光補正後の1次元プロファ
イル画像(856)を微分演算(857)し、微分画像
(858)を求める。図9(D’)は微分画像(85
8)の概形を示す。同図においてGauss曲線の幅の
狭い成分は検出器の分解能を、幅の広い成分は散乱線の
分布を示している。その微分画像(858)を2成分の
関数フィッティングして、散乱線強度比asと線像分布
関数LSFsを求め、これから点像強度分布関数を求め
る。以下、この方法につき詳述する。
【0113】まず、図9(A’)の配置でエッジの計測
画像Tsxを計測する(851)。エッジの計測画像Tsx
は散乱線成分を含み式3により、直接X線成分画像Ps
と拡散光成分画像Vsと散乱線成分画像Ssの和を示す次
式45となる。
【0114】
【数45】 Tsx=Ps+Vs+Ss …(式45) これらの各成分画像のデータは、すべてyに依存せずx
のみの関数となり、式46、式47、式48となる。
【0115】
【数46】 Ps=aps・ESFpx …(式46)
【0116】
【数47】 Ss=aps・asx・ESFsx …(式47)
【0117】
【数48】 Vs=(aps・STEPx+Ss)**(avx・PSFv) =aps・avx・ESFvx+aps・asx・avx・ESFsvx …(式48) ここで、ESFsvxは式49で定義されるエッジ分布関
数である。
【0118】
【数49】 ESFsvx=ESTsx*LSFvx …(式49) エッジの計測画像Tsは式46から式48を式45に代
入して、次式50となる。
【0119】
【数50】 Ts=aps・(ESFpx+asx・ESFsx+avx・ESFvx +asx・avx・ESFsvx) …(式50) 次に、Tsxに対する拡散光補正演算をconvolut
ion法により行なう。式37で求めたLSFvx(80
7)と式36で求めたavx(806)を利用して、以下
の演算(853)により近似的な拡散光画像Vse(85
4)を求める。
【0120】
【数51】 Vse=Tsx*(avx・LSFvx) …(式51) そして、式51に式50を代入すると、次式52を得
る。
【0121】
【数52】 Vse=aps・avx・(ESFpvx+asx・ESFsvx+avx・ESFvvx +asx・avx・ESFsvvx) …(式52) ここで、ESFpvx、ESFsvx、ESFvvx、ESF
svvxは、式53、式54、式55、式56により表わさ
れる。
【0122】
【数53】 ESFpvx=ESFpx*LSFvx …(式53)
【0123】
【数54】 ESFsvx=ESFsx*LSFvx …(式54)
【0124】
【数55】 ESFvvx=ESFvx*LSFvx …(式55)
【0125】
【数56】 ESFsvvx=ESFsvx*LSFvx …(式56) そして、拡散光補正演算を次の式57で行なう。
【0126】
【数57】Tsex=Tsx−Vse …(式57) 拡散光補正後の画像Tseは、式57に式50と式52を
代入して整理し、次式58を得る。
【0127】
【数58】 Tsex=aps・{ESFpx+asx・ESFsx+avx・(ESFvx−ESFpvx −avx・ESFvvx)−asx・avx 2・ESFsvvx} …(式58) そして、その微分画像(857)を求める。微分画像
(857)の概形を図9(D’)に示す。この微分画像
857は次式59となる。
【0128】
【数59】 d(Tsex)/dx=aps・{LSFpx+asx・LSFsx+avx・(LSFvx− LSFpvx−avx・LSFvvx)−asx・avx 2・LSFsvvx} …(式59) ここで、LSFsx、LSFpvx、LSFvvxは、式60、
式61、式62で表わされる。
【0129】
【数60】 LSFsx=d(ESFsx)/dx …(式60)
【0130】
【数61】 LSFpvx=d(ESFpvx)/dx …(式61)
【0131】
【数62】 LSFvvx=d(ESFvvx)/dx …(式62) 式59において、式63、式64、式65が成立するの
で、式66と近似できる。
【0132】
【数63】 LSFpvx≒LSFvx …(式63)
【0133】
【数64】 avx 2・LSFvvx≒0 …(式64)
【0134】
【数65】 asx・avx 2・LSFsvvx≒0 …(式65)
【0135】
【数66】 d(Tsex)/dx≒aps・(LSFpx+asx・LSFsx) …(式66) すなわち、散乱線を含み、拡散光成分を補正されたエッ
ジの微分画像(857)は2成分の線像分布関数の1次
結合により表現される(858)。
【0136】次に、この微分画像d(Tsex)/dxに
ついて、直接線成分の線像分布関数LSFpx及び散乱線
成分の線像分布関数LSFsxを正規分布関数と仮定して
2成分関数フィッティング(859)を行ない散乱線強
度比as(860)と散乱線成分の線像分布関数LSF
sx(861)を求める。
【0137】微分画像は例えば隣同志の画素の値の差分
から求められる。この場合、1ラインデータではX線量
子による統計的なノイズが多いので、複数ラインのデー
タを加算したデータを用いる。また、微分画像は左右対
称形となるが、ノイズは金属板で遮蔽されている側(図
9(D’)では左側)半分では小さく、遮蔽されていな
い側(図9(D’)では右側)半分では大きい。そこ
で、遮蔽されている側のデータをピークを中心として折
り返す。このようにして得られた曲線について2成分の
関数フィッティングを行なう。
【0138】関数フィッティングでは、式66におい
て、直接線成分の線像強度分布を式67、散乱線成分の
線像強度分布を式68と設定し、最小2乗法を利用した
フィッティングによりa3x、a4x、b3x、b4xを決定す
る。
【0139】
【数67】 aps・LSFpx=a3x・exp(−a4x・x2) …(式67)
【0140】
【数68】 aps・asx・LSFsx=b3x・exp(−b4x・x2) …(式68) 散乱線強度比asxは、フィッティングにより決定された
係数を用いて、式69で算出する。
【0141】
【数69】 asx=(b3x/a3x)・√(a4x/b4x) …(式69) そして、散乱線成分の線像分布関数は、式70で求め
る。
【0142】
【数70】 LSFsx=√(b4x/π)・exp(−b4x・x2) …(式70) Y方向についても金属板のエッジの方向をX方向につい
ての計測のときから90°回転して、視野の中央部にY
軸に直交するように配置して、X方向の場合と同様の計
測とフィッティングを行ない、式71、式72から、a
sy(862)とLSFsy(863)を決定する。
【0143】
【数71】 asy=(b3y/a3y)・√(a4y/b4y) …(式71)
【0144】
【数72】 LSFsy=√(b4y/π)・exp(−b4y・y2) …(式72) 散乱線の点像分布関数PSFsと散乱線強度比asを、式
73、式74で求める(864、865)。具体的に
は、式75、式76となる。
【0145】
【数73】 PSFs(x,y)=LSFsx(x)・LSFsy(y) …(式73)
【0146】
【数74】 as=(asx+asy)/2 …(式74)
【0147】
【数75】 PSFs(x,y)={√(b4x・b4y)/π}・exp(−b4x・x2−b4y・ y2) …(式75)
【0148】
【数76】 as={(b3x/a3x)・√(a4x/b4x)+(b3y/a3y)・√(a4y/b4y )}/2 …(式76 ) X方向とY方向の特性が近似的に等しいと仮定できる場
合には、式77となり、PSFv(x,y)は単一のパ
ラメータb4で決定される。但し、b4=b4x=b4 yであ
る。
【0149】
【数77】 PSFs(x,y)=(b4/π)・exp(−b4・x2−b4・y2) …(式77) 以上、詳述した方法により、散乱線の強度比と点像分布
関数を種々のアクリル板に対して計測を行ない、アクリ
ルの厚さと散乱線強度比の関係を求める。次に、予め種
々のアクリル厚に対して求めた散乱線強度比から、実際
の計測画像に対して用いる散乱線強度比の値を算出する
方法を説明する。
【0150】散乱線強度比asは、実験条件毎に異な
る。散乱線強度比に影響するパラメータは、前述したよ
うに、計測時の管電圧、被検体の厚さ、散乱X線遮蔽用
グリッド、視野径、被検体と上記グリッドの間の距離
(エアギャップとも呼ばれる)に依存する。パラメータ
の種類が多いので、できるだけ少数の実験データを用
い、できるだけ少数のパラメータとasとの関係をフィ
ッティングで求め、撮影条件を与えたときにasを算出
する方法を確定することが必要となる。asは位置によ
り変化しないと仮定しても、多くのパラメータの複雑な
関数となることから、以下、具体的に検討する。
【0151】Hondaらの既出の文献Medical
Physics誌(1991年)によれば、計測画像
上の直接線(1次X線)強度Ip及び散乱線強度Isは、
次の式78、式79で表わされる。
【0152】
【数78】 Ip=γ・Ep・Tp・P …(式78)
【0153】
【数79】 Is=γ・Es・Ts・f(A,F)・S …(式79) ここで、γはゲイン、Tsは散乱線遮蔽用グリッドの散
乱線透過率、Tpは散乱線遮蔽用グリッドの直接線透過
率、Aは被検体とグリッド間の距離(エアギャップ)
(cm)、FはX線照射視野のサイズ、f(A,F)は
エアギャップとX線照射視野サイズの関数であり、f
(1,∞)=1と規格化される。さらに、Es/Epは散
乱X線と1次X線の平均エネルギーの比、Pは被検体を
透過した直接X線の強度、Sは被検体から出射した散乱
X線の強度であり、次の式80、式81で表わされる。
ここで、Poは被検体への入射位置におけるX線ビーム
強度、μはX線の全線吸収係数、α1は直接X線の散乱
係数、g1は散乱線の発生確率に関し、実験的に求める
べき項である。
【0154】
【数80】 P=Po・exp(−μ・L) …(式80)
【0155】
【数81】 S=(α1・P/g1)・{exp(g1・L)−1} …(式81)
【0156】
【数82】 as=Is/Ip …(式82) そして、式82に式78、式79、式80、式81を代
入して式83を得る。
【0157】
【数83】 as=(α1/g1)・(Es/Ep)・f(A,F)・(Ts/Tp)・ {exp(g1・L)−1} …(式83) ここで、式83の一部を式84に示すようにGで置き換
え、これにより式85を得る。
【0158】
【数84】 G=(α1/g1)・(Es/Ep) …(式84)
【0159】
【数85】 as=G・f(A,F)・(Ts/Tp)・{exp(g1・L)−1} …(式85) ここで、G、Ts、Tp、及びg1は管電圧Vtの関数であ
る。また、f(A,F)は管電圧と被写体厚さに依存し
ない。そこで、まず、幾通りかの管電圧Vtに対して、
式85に現われるパラメータを決定してテーブルに記憶
しておき、テーブルにない管電圧で計測された画像に対
しては、補間演算でasを算出する。以下では、一通り
の管電圧に対する実験及び演算を説明する。まず、種々
のアクリル厚、例えば0、5、10、15、20、25
cmに対して、グリッドなし、及びグリッド有りの条件
で、種々の管電圧、例えば60、80、100、120
kVに対して、エッジの画像(2方向)を撮影する。エ
ッジの画像にグレア補正を行ない、その結果のLSFを
求める。その結果を2成分のGaussianでフィッ
ティングし、グリッド無し、及びグリッド有りの時の1
次X線強度をIpo、Ip、散乱線強度をIso、Isとして
求め、これらからasを求める。この実験及び演算の具
体的な方法は、図8及び図9を用いて既述した。次に、
グリッドの散乱線透過率Ts及び直接線透過率Tpを、次
の式86、式87で求める。
【0160】
【数86】 Ts=Is/Iso …(式86)
【0161】
【数87】 Tp=Ip/Ipo …(式87) これらの値はアクリル厚にほとんど依存しないので種々
のアクリル厚で求めた値の平均値をTs及びTpとする。
異なるグリッドごとに求める必要がある。
【0162】次に、横軸にアクリル厚、縦軸にasを取
りプロットする。このグラフを式85でフィッティング
し、2個のパラメータ、G・f(A,F)及びg1の値
を決定する。エアギャップAと照射視野サイズFについ
てはそれぞれ数通りの値につき上記の計測を行ない、デ
ータをテーブルに保存する。計測条件に応じて、これら
のA及びFの値を自動決定又は手動で選択し、内挿演算
等によりf(A,F)を決定する。さらに、Honda
らの方法に従い、式85からLを消去し、その式がLの
値にかかわらず成立すると考え、計測画像と実験条件か
らLが未知の被検体に対して適応できる式を求める。
【0163】そして、式80、式78から、次の式88
を得る。
【0164】
【数88】 L=(1/μ)・ln{(Io・Tp)/Ip} …(式88) ここで、Ioは1次X線強度であり、次の式89で表わ
される。
【0165】
【数89】 Io=γ・Ep・Po …(式89) そして、式88を式85に代入すると、次の式90を得
る(なお以下では、a↑bはべき乗、aのb乗を表わす
ものとする)。
【0166】
【数90】 Is=z1・(Ip↑(1−g))+z2・Ip …(式90) ここで、z1、z2、gは、式91、式92、式93によ
り与えられる。
【0167】
【数91】 z1=Ts・G・f(A,F)・(Io↑(g))/(Tp↑(1−g)) …(式91)
【0168】
【数92】 z2=−Ts・G・f(A,F)/Tp …(式92)
【0169】
【数93】 g=g1/μ …(式93) Ioは基準となる撮影条件における生のビームの強度の
ディジタル値として求める。具体的には、式94、式9
5により示される。
【0170】
【数94】 Io=Iref・OT・TC・ET・(100/D)2 …(式94)
【0171】
【数95】 Iref=γref・Ep・Pref …(式95) ここで、OTは光学系の絞りの伝達効率、TCはX線管
電流、ETはX線照射時間、DはX線管検出器間距離
(cm)、Irefは以下の条件における信号強度ディジ
タル値、即ち、TC・ET=1mAs、OT=1.0、
D=100cmである。式93におけるμの値は標準的
なアクリル厚、例えば15cmのときの透過率のデータ
から求める。
【0172】以上により式91のz1及び式92のz2
求まる。散乱線の補正演算を計算を簡単にできるよう、
式90において式96の近似を行なう。
【0173】
【数96】 Ip↑(1−g)=m・Ip …(式96) ここで、mは次の式97で表わされる。
【0174】
【数97】 m=Xp↑(−g) …(式97) ここで、XpはmaxIp(直接X線画像の最大値)の近
似解とする。Xpの近似値としては、計測画像の最大値
を用いる方法、計測画像から拡散光成分画像を補正した
画像Iqの最大値を用いる方法がある。Xpの近似値をよ
り精度良く求めるには、次の式98の関係を利用する。
【0175】
【数98】 Iq=Ip+Is=z1・(Ip↑(1−g))+(z2+1)・Ip …(式98) そして、maxIq(拡散光成分補正画像の最大値)を
第ゼロ近似解として、式99の近似解をNewton法
により求めると、式100となる。
【0176】
【数99】 maxIq=z1・(Xp↑(1−g))+(1+z2)・Xp …(式99)
【0177】
【数100】 Xp=(−z1・g・(maxIq↑(1−g))+maxIq)/ (z1・(1−g)・(maxIq↑(−g))+(1+z2)) …(式100) これらのいずれかの方法によりXpを求め、式97によ
りmを求める。
【0178】そして、式88を式81に代入すると、次
の式101を得る。
【0179】
【数101】 Is=(z1・m+z2)Ip …(式101) 従って、次の式102によりasを決定できる。
【0180】
【数102】 as=z1・m+z2 …(式102) 以上、Hondaらの上記文献の考え方を応用した結
果、asは予め計測されたパラメータ、画像を計測する
ときの条件、及び画像の値を用いて算出可能なことを示
した。
【0181】次に、実画像に対して適用するas、PS
s、av、PSFvを、上述したパラメータ等を用いて
算出する方法を図10及び図11を用いて説明する。
【0182】図10は、実際の画像計測前にX線透視・
撮影装置のメモリー内に記憶されているデータの内容例
を示している。メモリーは特定のグリッドを用いた時の
データを記憶したメモリー部A(1001)等と、グリ
ッドに依存しないデータを記憶したメモリー部B(10
05)がある。メモリー部A(1001)には、AとF
の値がそれぞれ標準値Ast及び標準値Fstのときに求め
たGの値、及びg1、Tp、Tsの値が幾通りかのX線管
電圧Vtに対して記憶されている(1002)。 ここ
で、Gとg1の値は、図8及び図9で説明した方法によ
り、AとFの値がそれぞれ標準値Ast及び標準値Fst
ときに種々の厚さのアクリルに対して式76によりas
を求め、横軸にアクリル厚さ、縦軸にasの計測値をプ
ロットし、式85の関数フィッティングを行なって、予
め決定されている。TpとTsは実験値からそれぞれ式8
7及び式86を用いて算出する。
【0183】また、メモリー部A(1001)には、幾
通りかのAとFの組み合わせに対して、f(A,F)の
値が記憶されている(1003)。ここで、FがFst
AがAstのときにはf(Ast,Fst)の値として、1.
0が記憶されている。f(A st,Fst)以外のf(A,
F)の値は、その条件における実験結果から式76で求
めたasの値とAst及びFstのときの実験結果から式7
6で求めたasの実測値との比の値とされる。なお、こ
の値はLに依存しない。
【0184】また、メモリー部A(1001)には、幾
通りかのFとVtの組み合わせに対してPSFsを決定す
るためのb4の値が記憶されており(1004)、これ
らの値と式77を利用してPSFsを求めることができ
る。
【0185】次に、メモリー部B(1005)には、幾
通りかのX線管電圧Vtに対して、式95により表わさ
れるIrefの実測値と、式80で定義される全線吸収係
数μの実測値が記憶されている(1006)。また、メ
モリー部B(1005)には、幾通りかの視野サイズF
に対して、拡散光成分の強度比avの実測値と拡散光成
分の点像分布関数PSFvを決定するb2の値とが記憶さ
れている(1007)。
【0186】図11は、実画像に適用するas、PS
s、av、PSFvの値を決定する具体的なフローを示
す。実画像の計測条件としては、X線管電圧Vt、視野
サイズF、エアギャップA、散乱線遮蔽用グリッド
r、光学絞りの相対効率OT、X線管電流TC、X線
照射時間ET、X線管検出器間距離D(cm)のデータを
用いる(1101)。まず、Grによりグリッドを指定
する(1102)。次に、Vtに対して、Gとg1とTp
とTsの値をテーブル1002から読み出す(110
3)。次に、FとAに対して、f(A,F)の値をテー
ブル1003から読み出す(1104)。また、一方
で、Vtに対して、Irefとμの値をテーブル1006か
ら読み出す(1105)。次に、次に示す式94の演算
を行なう(1106)。
【0187】
【数94】 Io=Iref・OT・TC・ET・(100/D)2 …(式94)
【0188】
【数93】 g=g1/μ …(式93) そして、式91、式92の演算を行なう(1109)。
【0189】
【数91】 z1=Ts・G・f(A,F)・(Io↑(g))/(Tp↑(1−g)) …(式91)
【0190】
【数92】 z2=−Ts・G・f(A,F)/Tp …(式92) 一方、計測画像It(1110)からその最大値max
tを求める(1111)。次に次式103の演算を行
なう(1112)。
【0191】
【数103】 m=(maxIt)↑(−g) …(式103) これらの結果から最後に式102により実画像に用いる
s(1114)を決定する(1113)。
【0192】
【数102】 as=z1・m+z2 …(式102) 上記演算とは別に、指定されたグリッドに対し、FとV
tに対して、PSFsを決定するパラメータb4の値をテ
ーブル1004から読みだし(1107)、式77を用
いてPSFsを算出する(1115)。
【0193】また、FとVtに対して、avの値、及びP
SFvを決定するパラメータb2の値をテーブル1007
から読み出す(1108)。avの値はそのまま用いる
(1116)。PSFvについては、式44を用いて算
出する(1117)。
【0194】実施例4.本発明の第3の目的である、撮
像素子の出力が飽和してハレーションが発生した場合に
も上記の補正を近似的に行なうX線装置の実施例を図1
2に示す。
【0195】本装置は、X線管(101)、散乱X線遮
蔽グリッド(102)、X線イメージインテンシファイ
ア(103)、集光光学系(121)、ハーフミラー
(1201)、透視及び撮影のための撮像素子(10
6)、上記撮像素子のための結像光学系(1202)、
飽和モニタ用受光素子配列(1203)、上記受光素子
のための結像光学系(1204)、撮像素子用アナログ
デジタル変換器(1205)、撮像素子用画像メモリ
(1206)、受光素子配列用アナログデジタル変換器
(1207)、受光素子配列用メモリ(1208)、飽
和レベル演算用演算回路(1209)等から構成され
る。
【0196】本実施例では撮像素子として1000×1
000画素のCCD素子、受光素子配列として5×5素
子のフォトダイオードの2次元配列を用いる。この場合
のフォトダイオードの出力は前記CCD素子が飽和する
X線条件下でも入力に比例した出力が得られるように設
定されている。受光素子用結像光学系(1204)は、
受光素子の出力が飽和しないよう絞りが調整される。撮
像素子の200×200画素の領域が受光素子の1素子
の領域に対応する。撮像素子の出力が飽和デジタル値d
sにちょうど達するときに、受光素子の出力のデジタル
変換値も同一の値dsとなるよう、各受光素子の出力電
圧を調整する。透視画像を連続的に計測する場合には、
受光素子用のアナログデジタル変換器はすべての受光素
子の出力信号を順次、繰返しデジタル変換し、すべての
受光素子について最新のX線像の値をメモリー上に保持
している。
【0197】演算回路(1209)において、撮像素子
(106)による計測画像の任意の画素値が飽和値であ
るか否かを判別する。判別以降の補正処理を図13に示
す。計測画像の画素(i,j)における値d(i,j)
(ステップ1301)が飽和値dsに等しいか否かを判
別する(ステップ1302)。飽和値に等しい場合に
は、撮像素子における画素(i,j)の位置に対応する
受光素子配列の位置座標(m,n)を算出し(ステップ
1303)、受光素子用のメモリーから画像の値a
(m,n)(1304)を読み出す(ステップ130
5)。次に、d(i,j)とa(m,n)の値を比較し
(ステップ1306)、a(m,n)がd(i,j)よ
り大の場合には、d(i,j)の真の値をa(m,n)
の値であると推定し、計測画像の値をd(i,j)から
a(m,n)に変更する(ステップ1307)。すべて
の画素について、この処理を行なって、計測画像の飽和
値の補正を行なう。この補正を行った後に、さきに説明
した拡散光と散乱線の補正(ステップ1308)を行な
う。これにより、撮像素子の出力が飽和してハレーショ
ンが発生した場合にも拡散光と散乱線の両方を近似的に
補正できる。
【0198】実施例5.本発明の第4の目的である、拡
散光と散乱線の補正を、特別な演算器なしのシステムで
も高速で実行可能とし、特に高速処理が要求される透視
像や超高精細画像やコーンビームCTなどにも適用可能
とする実施例を図14に示す。本実施例は、図7に示し
た実施例と同一の演算を、サンプリング操作を行って計
測画像とコンボリューションフィルタのマトリックスサ
イズを縮小した縮小画像と縮小コンボリューションフィ
ルタを生成し、縮小画像と縮小コンボリューションフィ
ルタを用いて拡散光と散乱線の補正演算を行なって、サ
ンプリング操作を行わない場合に比較して大幅に少ない
演算量で実行する。本サンプリングの実施例は、計測画
像を1024×1024画素、サンプリングピッチを計
測画像に対してもコンボリューションフィルタに対して
も16画素に統一し、縮小画像のサイズを64×64と
する。
【0199】図14において、計測画像(1401)の
取得以前に、予め、拡散光点像分布関数(1404)を
2次元サンプリングし(1407)、サンプリング拡散
光点像分布関数(1410)を生成する。また、散乱線
点像分布関数(1405)をサンプリングし(140
8)、サンプリング散乱線点像分布関数(1411)を
生成する。さらに、事前処理として、拡散光強度比(1
403)とサンプリング拡散光点像分布関数(141
0)との積(1412)によりサンプリング拡散光強度
分布関数(1414)を求める。
【0200】散乱線強度比(1402)は計測画像(1
401)及び計測条件から求める。その方法は前述し
た。散乱線強度比(1402)とサンプリング散乱線点
像分布関数(1411)との積(1413)によりサン
プリング散乱線強度分布関数(1415)を求める。計
測画像(1401)に対してサンプリング(1406)
を行ない、サンプリング計測画像(1409)を生成す
る。得られたサンプリング計測画像(1409)とサン
プリング拡散光強度分布関数(1414)をコンボリュ
ーション演算(1416)して、サンプリング拡散光成
分画像(1417)を求める。サンプリング計測画像
(1409)より得られたサンプリング拡散光成分画像
(1417)を減算することにより、サンプリング拡散
光補正画像(1419)を求め、この画像(1419)
とサンプリング散乱線強度分布関数(1415)とのコ
ンボリューション演算(1420)によりサンプリング
散乱線成分画像(1421)を算出する。求めたサンプ
リング散乱線成分画像(1421)とサンプリング拡散
光成分画像(1417)との加算(1422)により、
拡散光成分と散乱線成分のサンプリング合成画像(14
23)を得る。これをもとのマトリックスサイズに補間
演算により拡大し(1424)、拡散光・散乱線成分合
成画像(1425)を算出する。もとの計測画像(14
01)から、拡散光・散乱線成分合成画像(1425)
を減算して(1426)、拡散光成分と散乱線成分の補
正画像(1427)を得る。
【0201】このようにすれば、小規模の画像に対する
コンボリューション演算で補正処理が実行できるので、
特別な高速演算装置を必要としない。本実施例は、特に
高速処理が要求される透視像や、例えば2000×20
00画素や4000×4000画素の超高精細画像に対
しても、また、例えば512×512画素の画像を24
0枚撮影するコーンビームCTなどにも好適となる。
【0202】実施例6.本発明の第5の目的である、補
正演算に要する時間が透視像の1フレームの時間より長
い場合にも、透視像の近似的な補正を可能にする実施例
を図15に示す。本実施例においては、連続的に収集さ
れた透視画像に対し、透視画像を間歇的にサンプリング
し、サンプリング画像に対してのみ、拡散光と散乱線の
補正を行なう。ここで、連続する透視画像を透視画像
1、透視画像2、透視画像3、…と以下称する。図15
は、拡散光・散乱線画像を算出するのに透視画像の表示
の3フレーム分以上4フレーム分未満の時間を要する場
合の例で、透視画像1(1501)がサンプリング(1
502)された状態を示している。サンプリング画像に
対して、前述した飽和レベルの補正(1503)を行な
い、次に、前述した拡散光と散乱線の成分の演算処理
(1504)を前述したいずれかの方法を用いて行な
い、求めた拡散光・散乱線成分画像(1505)を対応
する透視画像1からではなく、次に表示される透視画像
5(1506)から減算(1507)して表示する(1
508)。
【0203】表示される画像は、計測画像そのものか、
又は計測画像と、拡散光・散乱線成分画像の差分画像と
される。本実施例では、最初の透視画像(透視画像1)
から透視画像4迄は、計測画像1から計測画像4が計測
されたままの状態で表示される。その間、透視画像1に
関して上記のサンプリングから拡散光・散乱線成分の演
算処理を実行中である。透視画像4が表示されている間
に演算処理が終了し、次の表示に際しては、透視画像5
と透視画像1の拡散光・散乱線成分画像(1505)の
差分演算画像(1507)を表示する(1508)。こ
のフレームでは、同時に透視像5がサンプリングされ、
ハレーションレベルの補正処理等が順次行われる。次の
透視画像6から透視画像8迄が表示される時刻には、こ
れらの表示される画像と差分演算を行なうための拡散光
・散乱線画像は依然として透視画像1に関するものであ
る。また、これらの時間においては、透視画像5がサン
プリングされ、拡散光と散乱線成分の画像を演算により
求めている途中である。表示画像が表示画像9に変わる
瞬間に、表示される拡散光・散乱線はそれまで処理中で
あった透視画像5が出現する。透視画像9−透視画像1
2間のフレームには差分処理を行なうための、拡散光・
散乱線の表示画像は透視画像5のまま一定である。この
間、次の拡散光・散乱線画像用として、透視像9が利用
されている。
【0204】このように、最後に補正項を透視画像から
減算するに際して、対応する透視画像ではなく次に表示
される最新の透視画像から減算して、透視画像のリアル
タイム表示を可能とするとともに、透視画像を数フレー
ム前の画像を用いて、近似的に補正する。
【0205】以上説明した実施例の効果を説明する。図
16(a)において、X線イメージインテンシファイア
103の前方の被検体を配置する個所に、階段状のアク
リル板160を配置する。すなわち、このアクリル板1
60は、階段状の形状となっていることによって、X線
管101からX線イメージインテンシファイア103へ
至るX線の前記アクリル板160を透過する厚さが場所
によって異なる(この場合、1、5、10、15、20
cm)ようになり、それに相当する出力が前記X線イメ
ージインテンシファイア103から得られ、かつ、上述
した画像補正がなされるようになっている。そして、前
記アクリル板160の階段を構成する高さのそれぞれ異
なる各平坦部には、その長手方向に沿った金属棒160
A(6mm厚、5mm幅)が被着されている。
【0206】同図(b)は、得られた補正画像の一本の
ライン上の画像の値のプロファイルを示し、また同図
(c)は前記同図(b)の一部を拡大したプロファイル
を示している。
【0207】同図(c)から明らかになるように、金属
棒160Aが配置されている個所において、レベルが0
に近づいていることが判明する。ちなみに、本発明のよ
うな補正を行わなかった場合の、図16(b)、(c)
のそれぞれに対応する図を図17(a)、(b)に示し
ている。この場合、金属棒160Aが配置されている個
所において、レベルが0から遠ざかっていることからし
て、上記補正画像との相違が顕著となる。
【0208】即ち、図16及び図17に示すように、ア
クリル厚10cmの場合を例にとると、信号(コントラ
スト差)とバイアス成分との比は補正前の画像では約
5.7、補正後の画像では約22となり、本発明により
約3.8倍改善されている。これは、散乱線や拡散光に
よる画像のコントラストの低下を効果的に補正できてい
ることを示している。
【0209】以上、本発明をX線画像作成方法として説
明したが、本発明は、このX線画像作成方法が適用され
るX線透視撮影装置をも包含されるものである。
【0210】すなわち、このX線透視撮影装置は、X線
管と、散乱線遮蔽用グリッドと、X線像を光学像に変換
する手段と、光学像を電気信号に変換する手段と、アナ
ログデジタル変換器と、デジタル画像収集装置と、デジ
タル画像処理装置等から構成され、透視又は撮影条件に
対する散乱線の強度比相対値を記憶する第1の記憶手段
と、透視又は撮影条件に対する散乱線成分の点像分布関
数を記憶する第2の記憶手段と、X線管電圧に対する基
準X線条件の画像の値を記憶する第3の記憶手段と、視
野サイズに対する拡散光成分の強度比と拡散光成分の点
像分布関数を記憶する第4の記憶手段とをもち、計測画
像の値と、上記の記憶手段に記憶されたデータを用い
て、計測画像の拡散光成分と散乱線成分を補正する手段
を有する装置である。
【0211】また、このようなX線透視撮影装置におい
て、第1の記憶手段が、X線管電圧に対する散乱線の強
度比相対値と、視野サイズとエアギャップに対する散乱
線の強度比(相対値)とを記憶する手段であり、第2の
記憶手段が、視野サイズとX線管電圧に対する散乱線成
分の点像分布関数を記憶する手段であり、第3の記憶手
段が、X線管電圧に対する基準X線条件の画像の値を記
憶する手段であり、第4の記憶手段が、視野サイズに対
する拡散光成分の強度比と拡散光成分の点像分布関数を
記憶する手段であり、第1から第3の記憶手段は相異な
るグリッド毎に個別の値を記憶する手段を有している。
【0212】さらに、X線像を光学像に変換する手段が
X線イメージインテンシファイアであり、光学像を電気
信号に変換する手段が結像光学系と撮像素子であり、結
像光学系がイメージインテンシファイアの出力像を2ヵ
所に結像する手段をもち、結像光学系により結像された
画像の一個を透視像又は撮影像を計測する撮像素子によ
り計測する手段と、結像光学系により結像された画像の
残りの一個をハレーションレベルを計測する受光素子配
列により計測する手段と、撮像素子により計測された最
新の画像信号のディジタル値及び受光素子により計測さ
れた各素子の最新の出力信号のディジタル値を記憶する
手段と、撮像素子の出力が飽和したか否かを判別する手
段と、撮像素子の出力が飽和した場合に受光素子配列の
出力信号を用いて撮像素子の飽和レベル出力を補正する
演算手段と、飽和レベルの補正後に拡散光と散乱線の両
方の補正を行なう手段を有している。
【0213】さらに、実施例で既述した別の方法を具体
化する装置は、透視画像を連続的にディジタル画像とし
て取得する手段と、透視画像を間歇的にサンプリングす
る手段と、サンプリング画像に対して、拡散光成分と散
乱線成分の合成画像を生成する手段と、次に表示される
べき計測画像から拡散光成分と散乱線成分の合成画像を
減算して表示する手段を有している。
【0214】本発明は、上述した実施例に記述した形態
だけでなく、X線撮影の広範な分野への適用が可能であ
る。その一実施例としてコーンビームCT装置に適用で
きる。すなわち、空気のみを撮影したエア画像を記憶す
る手段と、エア画像の拡散光成分を補正して感度分布画
像を生成する手段と、計測画像から拡散光成分と散乱線
成分を補正して得られた画像と感度分布画像との対数差
分演算を行なう手段と、計測系の幾何学的歪を補正する
手段と、を有するコーンビームCT装置は、拡散光と散
乱線が補正されたデータを用いて3次元画像再構成を行
うことができるため、得られるCT画像のCT値の正確
度が従来に比較して顕著に向上する。
【0215】本発明の応用の別の実施例として、X線フ
ィルムを用いるX線撮影装置とフィルムディジタイザの
組み合わせからなるX線撮影装置、あるいは、輝尽性蛍
光板を用いるX線撮影装置等がある。より具体的には、
X線像を光学像に変換する手段が、X線増感紙フィルム
系(第1の手段)と、第1の手段に記録された光学像を
光学的に読み出すフィルムディジタイザの読み出ビーム
(第2の手段)とから構成され、光学像を電気的に変換
する手段がフィルムディジタイザの光検出素子であるX
線撮影装置に適用される。
【0216】また、他のX線撮影装置として、X線像を
光学像に変換する手段として、輝尽性蛍光板を光学的に
読み出す輝尽性蛍光板読み取り装置のレーザビーム系で
あり、光学像を電気的に変換する手段として、輝尽性蛍
光板読み取り装置の光検出素子である場合等に適用でき
る。
【0217】
【発明の効果】以上説明したことから明らかなように、
本発明によるX画像作成方法およびその装置によれば、
まず、X線透視画像あるいはX線撮影画像に対して、周
波数空間で高周波強調フィルタの乗算を行なうことによ
る高周波ノイズの増加なしに、拡散光と散乱線に起因す
るX線画像のぼけをそれぞれの点像分布関数を用いて正
確に補正し、これにより、厚さの変化範囲が大きな被検
体に対しても、高精度で補正することができるようにな
る。
【0218】また、上記の補正に必要となる拡散光強度
分布関数と散乱線強度分布関数とを、簡便な方法により
求めることができるようになる。
【0219】また、撮像素子の出力が飽和してハレーシ
ョンが発生した場合にも、上記の補正を近似的に行なう
ことができるようになる。
【0220】また、上記の補正を、特別な演算器なしの
システムでも高速で実行可能とし、特に高速処理が要求
される透視像や、例えば2000×2000画素や40
00×4000画素の超高精細画像に対しても、また、
例えば512×512画素の画像を240枚撮影するコ
ーンビームCT等にも適用することができるようにな
る。
【0221】さらに、補正演算に要する時間が透視像の
1フレームの時間より長い場合にも連続的に変化する透
視像に対する近似的な補正をすることができるようにな
る。
【図面の簡単な説明】
【図1】本発明によるX線画像作成方法およびその装置
における画像劣化モデルおよび補正演算を説明する概略
フロー図である。
【図2】本発明に使用されるX線検出系の代表例を示す
構成図である。
【図3】従来法(直接デコンボリューション法)による
画像劣化モデルおよび補正演算の一例を示したフロー図
である。
【図4】従来法(ぼけ画像生成法)による画像劣化モデ
ルおよび補正演算の他の例を示したフロー図である。
【図5】本発明によるX線画像作成方法およびその装置
において適用される飽和発生時の近似的補正を示すフロ
ー図である。
【図6】本発明の第1の実施例であるフーリエ変換法を
用いる補正演算を示すフロー図である。
【図7】本発明の第2の実施例であるコンボリューショ
ン法を用いる補正演算を示すフロー図である。
【図8】本発明によるX線画像作成方法およびその装置
において適用される拡散光及び散乱線の強度比及び点像
分布関数を求めるフロー図である。
【図9】本発明によるX線画像作成方法およびその装置
において適用される拡散光及び散乱線の強度比及び点像
分布関数を求める図である。
【図10】本発明による補正演算において、実画像に適
用する拡散光、散乱線の強度比、及び点像分布関数の算
出に用いる、X線透視・撮影装置のメモリー内に記憶さ
れているデータ内容を示す図である。
【図11】本発明による補正演算において、実画像に適
用する拡散光、散乱線の強度比、及び点像分布関数を求
めるフロー図である。
【図12】本発明による実施例を示す図で、飽和発生時
に近似的補正処理を行なうX線撮影装置の構成を示す図
である。
【図13】本発明によるX線画像作成方法およびその装
置において適用される飽和発生時に近似的補正を行なう
フロー図である。
【図14】本発明によるX線画像作成方法およびその装
置において適用される特別な演算器なしで補正演算を高
速に実行するフロー図である。
【図15】本発明によるX線画像作成方法およびその装
置において適用されるものであって、補正演算に要する
時間が透視像の1フレームの時間より長い場合に、透視
像の近似的補正演算を行なうフロー図である。
【図16】本発明の効果を説明するための図である。
【図17】本発明の効果と比較される図である。
【符号の説明】
101…X線管、102…散乱X線遮蔽グリッド、10
3…X線イメージインテンシファイア、104…出力蛍
光面、105…ミラー、106…撮像素子、107…直
接光、108…拡散光成分、109…電子ビーム、11
0…被検体、111…直接X線、112…散乱X線、1
21…集光光学系、1201…ハーフミラー、1202
…撮像素子用結像光学系、1203…飽和モニタ用受光
素子、1204…受光素子用結像光学系、1205…撮
像素子用アナログディジタル変換器、1206…撮像素
子用画像メモリ、1207…受光素子配列用アナログデ
ィジタル変換器、1208…受光素子配列用メモリ、1
209…飽和レベル演算用演算回路。
───────────────────────────────────────────────────── フロントページの続き (72)発明者 梅谷 啓二 東京都国分寺市東恋ケ窪1丁目280番地 株式会社日立製作所 中央研究所内 (56)参考文献 特開 昭59−151939(JP,A) 特開 平6−14911(JP,A) 特開 平1−125679(JP,A) (58)調査した分野(Int.Cl.7,DB名) A61B 6/00 - 6/14

Claims (5)

    (57)【特許請求の範囲】
  1. 【請求項1】 被検体がない状態で計測されたエッジ強
    度分布から算出される拡散光の強度比と前記拡散光の点
    像分布関数との積により前記拡散光の強度分布関数を求
    め、前記拡散光の強度分布関数と計測された前記被検体
    の画像との演算により拡散光成分の画像を算出する工程
    と、均一な厚さの人体模擬被写体を配置した状態で計測
    されたエッジ強度分布から求めた散乱線の点像分布関数
    と、前記被検体の画像とその撮影条件から求めた散乱線
    の強度比との積により散乱線の強度分布関数を求め、前
    記被検体の画像と前記散乱線の強度分布関数との演算に
    より散乱線成分の画像を算出する工程と、前記被検体の
    画像から前記拡散光成分の画像及び前記散乱線成分の画
    像を減算する工程とを有し、前記被検体のX線像を光学
    像に変換するX線イメージインテンシファイアの出力蛍
    光面における前記拡散光成分の画像と、前記被検体によ
    る前記散乱線成分の画像とを別々に作成し、前記被検体
    の画像から前記拡散光成分及び前記散乱線成分を除く補
    正を行うことを特徴とするX線画像作成方法。
  2. 【請求項2】 予め求めた拡散光の強度比と予め求めた
    前記拡散光の点像分布関数との積により前記拡散光の強
    度分布関数を求め、前記拡散光の強度分布関数と計測さ
    れた被検体の画像とのコンボリューション演算によって
    拡散光成分の画像を算出する工程と、前記被検体の画像
    とその撮影条件から求めた散乱線の強度比と予め求めた
    前記散乱線の点像分布関数との積により求めた散乱線の
    強度分布関数と、前記被検体の画像とのコンボリューシ
    ョン演算によって散乱線成分の画像を算出する工程と、
    前記被検体の画像から前記拡散光成分の画像及び前記散
    乱線成分の画像を減算する工程とを有し、前記被検体の
    画像から前記拡散光成分及び前記散乱線成分を除く補正
    を行うことを特徴とするX線画像作成方法。
  3. 【請求項3】 被検体がない状態で計測されたエッジ強
    度分布から算出される拡散光の強度比と前記拡散光の点
    像分布関数との積により前記拡散光の強度分布関数を求
    め、前記拡散光の強度分布関数と計測された前記被検体
    の画像との演算により拡散光成分の画像を算出する処理
    と、均一な厚さの人体模擬被写体を配置した状態で計測
    されたエッジ強度分布から求めた散乱線の点像分布関数
    と、前記被検体の画像とその撮影条件から求めた散乱線
    の強度比との積により散乱線の 強度分布関数を求め、前
    記被検体の画像と前記散乱線の強度分布関数との演算に
    より散乱線成分の画像を算出する処理とを行ない、前記
    被検体の画像から前記拡散光成分の画像及び前記散乱線
    成分の画像を減算する手段を有し、前記被検体のX線像
    を光学像に変換するX線イメージインテンシファイアの
    出力蛍光面における前記拡散光成分の画像と、前記被検
    体による前記散乱線成分の画像とを別々に作成し、前記
    被検体のX線透視画像又はX線撮影画像から前記拡散光
    成分及び前記散乱線成分を除く補正を行うことを特徴と
    するX線透視撮影装置。
  4. 【請求項4】 予め求めた拡散光の強度比と予め求めた
    前記拡散光の点像分布関数との積により前記拡散光の強
    度分布関数を求め、前記拡散光の強度分布関数と計測さ
    れた被検体の画像とのコンボリューション演算によって
    拡散光成分の画像を算出する処理と、前記被検体の画像
    とその撮影条件から求めた散乱線の強度比と予め求めた
    前記散乱線の点像分布関数との積により求めた前記散乱
    線の強度分布関数と、前記被検体の画像とのコンボリュ
    ーション演算によって散乱線成分の画像を算出する処理
    とを行い、前記被検体の画像から前記拡散光成分の画像
    及び前記散乱線成分の画像を減算する手段を有し、前記
    被検体の画像から前記拡散光成分及び前記散乱線成分を
    除く補正を行うことを特徴とするX線透視撮影装置。
  5. 【請求項5】 空気のみを撮影したエア画像の拡散光成
    分を補正して感度分布画像を生成する手段と、計測され
    た被検体の画像から拡散光成分と散乱線成分を除く補正
    を行い得られた画像と前記感度分布画像との対数差分演
    算を行う手段と、予め求めた拡散光の強度比と予め求め
    た前記拡散光の点像分布関数との積により前記拡散光の
    強度分布関数を求め、前記拡散光の強度分布関数と計測
    された被検体の画像とのコンボリューション演算によっ
    て前記拡散光成分の画像を算出する処理と、前記被検体
    の画像とその撮影条件から求めた散乱線の強度比と予め
    求めた前記散乱線の点像分布関数との積により求めた前
    記散乱線の強度分布関数と、前記被検体の画像とのコン
    ボリューション演算によって前記散乱線成分の画像を算
    出する処理とを行い、前記被検体の画像から前記拡散光
    成分の画像及び前記散乱線成分の画像を減算する手段と
    を有し、前記被検体の画像から前記拡散光成分及び前記
    散乱線成分を除く補正を行ない、拡散光と散乱線が補正
    されたデー タを用いて3次元画像再構成を行なうことを
    特徴とするコーンビームCT装置。
JP31184195A 1995-11-30 1995-11-30 X線画像作成方法およびその装置 Expired - Fee Related JP3423828B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP31184195A JP3423828B2 (ja) 1995-11-30 1995-11-30 X線画像作成方法およびその装置
US08/759,088 US5960058A (en) 1995-11-30 1996-11-29 Method for generating X-ray image and apparatus therefor

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP31184195A JP3423828B2 (ja) 1995-11-30 1995-11-30 X線画像作成方法およびその装置

Publications (2)

Publication Number Publication Date
JPH09149895A JPH09149895A (ja) 1997-06-10
JP3423828B2 true JP3423828B2 (ja) 2003-07-07

Family

ID=18022053

Family Applications (1)

Application Number Title Priority Date Filing Date
JP31184195A Expired - Fee Related JP3423828B2 (ja) 1995-11-30 1995-11-30 X線画像作成方法およびその装置

Country Status (2)

Country Link
US (1) US5960058A (ja)
JP (1) JP3423828B2 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007125129A (ja) * 2005-11-02 2007-05-24 Hitachi Medical Corp X線ct装置
US9569826B2 (en) 2014-09-30 2017-02-14 Fujifilm Corporation Radiographic image processing device, radiographic image processing method, and recording medium

Families Citing this family (40)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE19803588A1 (de) * 1998-01-30 1999-08-05 Christof Dr Steiner Verfahren und Anordnung zum Auslesen von auf einer Bildplatte gespeicherter Strahlungsbildinformation
GB9904692D0 (en) 1999-03-01 1999-04-21 Isis Innovation X-ray image processing
JP3987676B2 (ja) 2000-07-10 2007-10-10 株式会社日立メディコ X線計測装置
US7283253B2 (en) * 2002-03-13 2007-10-16 Applied Precision, Llc Multi-axis integration system and method
EP1495446B1 (en) * 2002-04-05 2007-10-10 The General Hospital Corporation Systems and method for generating an image
CA2391132A1 (en) * 2002-06-21 2003-12-21 Dan Rico Method and apparatus for determining peripheral breast thickness
US6876718B2 (en) * 2003-06-27 2005-04-05 Ge Medical Systems Global Technology Company, Llc Scatter correction methods and apparatus
JP4175468B2 (ja) * 2003-07-24 2008-11-05 学校法人日本大学 画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体
US7492967B2 (en) * 2003-09-24 2009-02-17 Kabushiki Kaisha Toshiba Super-resolution processor and medical diagnostic imaging apparatus
JP4533010B2 (ja) * 2003-11-20 2010-08-25 キヤノン株式会社 放射線撮像装置、放射線撮像方法及び放射線撮像システム
US7382853B2 (en) * 2004-11-24 2008-06-03 General Electric Company Method and system of CT data correction
US20070086560A1 (en) * 2005-10-06 2007-04-19 Imaging Sciences International, Inc. Scatter correction
JP2008073208A (ja) * 2006-09-21 2008-04-03 Konica Minolta Medical & Graphic Inc 画像処理装置及び画像処理方法
JP2008073342A (ja) * 2006-09-22 2008-04-03 Konica Minolta Medical & Graphic Inc 放射線画像撮影システム及び放射線画像撮影方法
WO2008050320A2 (en) * 2006-10-23 2008-05-02 Ben Gurion University Of The Negev, Research And Development Authority Blind restoration of images degraded by isotropic blur
AR057580A1 (es) * 2006-11-14 2007-12-05 Tomografia De Hormigon Armado Metodo y disposicion para mejorar las determinaciones tomograficas por medio de radiaciones, especialmente apropiado para barras de acero en el hormigon armado
JP5220383B2 (ja) * 2007-10-29 2013-06-26 株式会社日立メディコ 放射線撮像装置
JP5416377B2 (ja) * 2008-08-28 2014-02-12 アンリツ産機システム株式会社 画像処理装置及びそれを備えたx線異物検出装置並びに画像処理方法
JP4683247B2 (ja) * 2009-07-21 2011-05-18 花王株式会社 画像形成方法
US8606032B2 (en) * 2009-12-23 2013-12-10 Sony Corporation Image processing method, device and program to process a moving image
JP2012024344A (ja) * 2010-07-23 2012-02-09 Canon Inc X線撮影装置、x線撮影方法、プログラム及びコンピュータ記録媒体
US20120163695A1 (en) * 2010-12-22 2012-06-28 General Electric Company Intra-detector scatter correction
WO2012173703A2 (en) * 2011-04-29 2012-12-20 Bp Corporation North America Inc. System and method for underwater radiography
KR102020531B1 (ko) * 2012-08-20 2019-09-10 삼성전자주식회사 선형 감마선원을 이용하여 고해상도의 pet(양전자 방출 단층 촬영) 영상을 생성하는 방법 및 장치
KR102026735B1 (ko) * 2012-10-02 2019-09-30 삼성전자주식회사 영상 촬영 장치의 검출기의 시스템 응답 및 시스템 응답을 이용하여 의료 영상을 생성하는 방법 및 장치
JP6006454B2 (ja) * 2013-03-28 2016-10-12 富士フイルム株式会社 放射線画像処理装置および方法並びにプログラム
JP6373952B2 (ja) * 2013-07-31 2018-08-15 富士フイルム株式会社 放射線画像解析装置および方法並びにプログラム
JP6128463B2 (ja) 2013-11-06 2017-05-17 富士フイルム株式会社 放射線画像処理装置および方法並びにプログラム
JP6071853B2 (ja) 2013-11-26 2017-02-01 富士フイルム株式会社 放射線画像処理装置、方法およびプログラム
JP6129125B2 (ja) * 2013-11-29 2017-05-17 富士フイルム株式会社 放射線画像解析装置および方法並びにプログラム
JP6105508B2 (ja) * 2014-03-28 2017-03-29 富士フイルム株式会社 放射線画像撮影装置、放射線画像撮影方法、及び放射線画像撮影プログラム
JP6042855B2 (ja) * 2014-03-28 2016-12-14 富士フイルム株式会社 放射線画像撮影装置、放射線画像撮影方法、及び放射線画像撮影プログラム
JP6313402B2 (ja) * 2014-03-28 2018-04-18 富士フイルム株式会社 放射線画像撮影装置、放射線画像撮影方法、及び放射線画像撮影プログラム
US9681073B1 (en) * 2014-08-25 2017-06-13 Marvell International Ltd. Method and apparatus for compensation of veiling glare in an image capturing device
JP6165695B2 (ja) * 2014-09-24 2017-07-19 富士フイルム株式会社 放射線画像解析装置および方法並びにプログラム
JP6296553B2 (ja) * 2014-09-30 2018-03-20 富士フイルム株式会社 放射線画像撮影装置および放射線画像撮影装置の作動方法
JP6525772B2 (ja) * 2015-06-30 2019-06-05 キヤノン株式会社 画像処理装置、画像処理方法、放射線撮影システムおよび画像処理プログラム
JP6653629B2 (ja) * 2016-06-21 2020-02-26 富士フイルム株式会社 放射線画像処理装置、方法およびプログラム
JP6325141B2 (ja) * 2017-03-02 2018-05-16 富士フイルム株式会社 放射線画像撮影装置、放射線画像撮影方法、及び放射線画像撮影プログラム
KR102003197B1 (ko) * 2017-10-17 2019-07-24 연세대학교 원주산학협력단 단일 방사선 영상의 산란선 교정 시스템 및 방법

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4718076A (en) * 1983-04-22 1988-01-05 Kabushiki Kaisha Toshiba X-ray imaging apparatus
JPS62191972A (ja) * 1986-02-18 1987-08-22 Toshiba Corp X線画像処理装置
US4841555A (en) * 1987-08-03 1989-06-20 University Of Chicago Method and system for removing scatter and veiling glate and other artifacts in digital radiography
DE69308350T2 (de) * 1992-04-08 1997-08-21 Philips Electronics Nv Vorrichtung zur Röntgenuntersuchung mit Korrektur der Streustrahlungseffekte in einem Röntgenbild

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007125129A (ja) * 2005-11-02 2007-05-24 Hitachi Medical Corp X線ct装置
US9569826B2 (en) 2014-09-30 2017-02-14 Fujifilm Corporation Radiographic image processing device, radiographic image processing method, and recording medium
US9679368B2 (en) 2014-09-30 2017-06-13 Fujifilm Corporation Radiographic image processing device, radiographic image processing method, and recording medium

Also Published As

Publication number Publication date
JPH09149895A (ja) 1997-06-10
US5960058A (en) 1999-09-28

Similar Documents

Publication Publication Date Title
JP3423828B2 (ja) X線画像作成方法およびその装置
US5878108A (en) Method for generating X-ray image and apparatus therefor
JP2606828B2 (ja) X線像を較正する装置
JP3112455B2 (ja) X線映像の散乱放射効果補正装置
US8009892B2 (en) X-ray image processing system
JP4020966B2 (ja) 画像内のノイズ減少方法
KR101850871B1 (ko) 방사선 영상의 처리방법 및 방사선 촬영시스템
EP0652537B1 (en) Method of and apparatus for computed tomography
JPH06259541A (ja) 画像歪み補正方法およびそのシステム
JP3540914B2 (ja) X線撮影装置
JP3583554B2 (ja) コーンビームx線断層撮影装置
EP0116941B2 (en) An x-ray diagnostic apparatus
JP2849964B2 (ja) 画像処理方法および装置
US5402463A (en) Apparatus and method for radiation imaging
JPH0448453B2 (ja)
JP4746761B2 (ja) 放射線画像処理装置、放射線画像処理方法、記憶媒体、及びプログラム
JP3349004B2 (ja) X線画像計測装置
Monnin et al. A comparison of the performance of modern screen-film and digital mammography systems
Hoeschen et al. High-spatial-resolution measurement of x-ray intensity pattern in a radiograph of the thorax
JPH0534710B2 (ja)
US20230401677A1 (en) Image processing apparatus, radiation imaging system, image processing method, and non-transitory computer-readable storage medium
Friedman et al. A method to measure the temporal MTF to determine the DQE of fluoroscopy systems
JPH05153496A (ja) X線撮影装置,x線透視装置およびこれらに用いるx線テレビカメラ装置
JP2001292324A (ja) 画像処理装置および画像処理方法
Kawahara et al. Spatial frequency components of normal radiographic anatomical features on intraoral computed radiography

Legal Events

Date Code Title Description
LAPS Cancellation because of no payment of annual fees