JP3698167B2

JP3698167B2 - Ｘ線画像形成装置の幾何学的較正を行うための装置と、それを自動的に行う方法

Info

Publication number: JP3698167B2
Application number: JP02481894A
Authority: JP
Inventors: ピカールカトリーヌ; ルージュアンヌ; トルッセイヴ
Original assignee: ジーイーメディカルシステムズソシエテアノニム
Priority date: 1993-01-27
Filing date: 1994-01-27
Publication date: 2005-09-21
Anticipated expiration: 2020-09-21
Also published as: FR2700909A1; JPH0716220A; FR2700909B1; US5442674A

Description

【０００１】
【産業上の利用分野】
本発明はＸ線画像形成装置の幾何学的較正を行う装置と、それを自動的に行う方法とに関するものである。
本発明の装置は主として医療分野で人体の内部構造、特に血管構造を解析するために使用される。
本発明の目的はＸ線画像形成装置の較正（キャリブレーション）を容易にすることによって、人体の内部構造を定性的に知るだけでなく定量的に分析できるようにすることにある。
本発明は主として２次元放射線検出器、例えば放射線平面検出器を用いた装置に関するが、本発明方法はＸ線画像形成装置の回転装置にも関連する。
【０００２】
【従来の技術】
Ｘ線源と２次元検出器とで構成されるＸ線画像形成装置（回転式と非回転式がある）で画像形成に使われるものは、３次元空間内での分析対象物の２次元空間（これは検出平面に対応する）上に幾何学的に広がる円錐形の投影画像である。この円錐形投影画像を完全に記述するのに必要な幾何学パラメータの数は９つである。これらの９つのパラメータはパラメータ化方法の違いによって種々の形で表される。パラメータ化方法の１例を以下に説明する。
(1) 最初の３つのパラメータは、分析対象物が配置される３次元空間の基準系Ｒ内でのＸ線源Ｓの３つの座標に関するものである。
(2) 次の３つのパラメータは３次元空間の基準系Ｒと円錐形投影画像に関する基準系Ｒ’すなわち射影面Ｐの２次元基準系とこの射影面と直角な第３の軸とから成る基準系との間の基準系の変化に関連する３つのオイラー(Euler) 角度に関するものである。
(3) 最後の３つのパラメータは射影面Ｐの任意の点Ｉ（例えば、基準系Ｒの画像の中心）の３つの座標に関するものである。
【０００３】
本発明では計算を単純化するために、以下で説明するように、最初の３つの座標および最後の３つの座標として円筒座標を用いる。従って、基準系Ｒも円筒基準系である。
これらのパラメータの全てまたは一部を知ることはＸ線撮影法(radiology) では多くの場合有用であり、特に３次元の物体に関する定量データを物体の２次元投射像または複数の２次元投射像から取る場合に有用である。
この場合、各視野(viewpoint) のパラメータが分かっている必要があるが、この視野は物体に対する画像形成装置の向きに関連している。例えば１枚または複数の２次元射投影像から血管断面を測定するという標準的な問題を考えれば分かるように、これらのパラメータに直接アクセスすることは不可能 (例えば画像形成装置でＸ線源と検出器との間の距離を直接測定することはできない) か、例えできたとしても精度が著しく悪い。
【０００４】
「画像形成装置の幾何学的な較正」とは、画像形成で役割を果たす幾何学的なパラメータに関する間接的な知識を正確なものにするための操作を意味する。その原理自体は公知であり、３次元空間的寸法が分かっている幾何学的ファントム(phantom) を用いてその２次元投影像を得るという原理である。これで用いられる操作シーケンスは下記段階で構成される：
(1) 寸法が既に分かっている物体、すなわち、この物体に固有な基準系に対して測定された座標によって空間内での位置が分かっている一定数の特定点を有する較正ファントムを用い、
(2) 較正すべき視野（または入射角) の幾何学的条件下でこの較正ファントムの画像を取り、
(3) この画像の特性点の投影像を知るために、物体の各特性点を取った画像すなわち投影像での物体の軌跡と関連付け、
(4) 投影像を記述する方程式を数学的に反転し、
(5) 所定の視野での投影像の全てのパラメータが得られる。
【０００５】
これまでは較正ファントムの投影像での特性点の確認はオペレータが「手で」行っている。一般に用いられる較正ファントムの幾何学形状は立方体で、この立方体の８つの隅に光線を通さない金属のビーズまたは小球が取付けられている。較正精度を大きくするために追加のビーズが付いている場合もある。Ｘ線管の方向に沿った投影像は物体の「透明」像を作るが、この像はＸ線管の向きによって変わるので、ビーズの２次元軌跡を対応する３次元軌跡と間違いなく組み合わせることは極めて困難である。
従来の幾何学的較正に関しては下記論文を参照されたい。
【０００６】
(1) パーカー(D. L. Parker)、ウー(J. Wu) 、ポプ(D. L. Pope)、ブリー(R. Va n Bree) 、カプト(G. R. Caputo)およびマーシャル(H. W. Marshall)著『マルチビューデジタル血管造影法を用いた冠状動脈の３次元再構成と流量測定 (Three-Dimensional Reconstruction And Flow Measurements of Coronary Arteries Using Multiview Digital Angiography) 』〔レイバー(J. C. Rei ber)およびサーリュイ(P. W. Serrruys)編『量的な冠状動脈造影法の新しい展開(New Developments in Quantitative Coronary Angiography) 』、クリュワーアカデミックパブリシャーズ(Kluwer Academic Publishers)、1988年、 247頁〕
(2) ホークス(D. J. Hawkes)、コルチェスター(A. C. F. Colchester) およびモル (C. R. Mol)著『Ｘ線記録からの血管ツリーの幾何学的形態の正確な３次元再構成(The Accurate 3D Reconstruction of the Geometric Configurati on of the Vascular Trees from X-Ray Recordings) 』〔グッツアルディ(R . Guzzardi) 編『医学画像化の物理学及びエンジニアリング(Physics and E ngineerings of Medical Imagings)』、ニジョフ(Nijoff)、1987年〕
(3) ガロー(M. Garreau)、コートリュー(J-L Coatrieux) 、コロレック(R. Coll orec) およびシャルドノン(C. Chardenon)著『２平面血管造営法射影からの血管網の３次元再構成とラベル付けの知識に基づく方法(A Knowledge-Based Approach for 3D Reconstruction and Labellings of Vascular Networks FromBiplane Angiographic Projections) 』〔「アイイーイーイーメディカルイメージングス（IEEE Medical Imagings)」、第10巻、第２号、1991年６月、122 〜131 頁〕
【０００７】
人間が介入する必要があるということは大きな欠点であり、場合によって、特に、較正すべき視野の数が多い場合には禁止しなければならない場合もある。このことは３次元の再構成のためにＸ線による２次元投影像を取る装置の場合には常にいえることである。この種の画像形成装置は例えば1989年３月20日に出願されたフランス国特許第 2,644,590号に記載されている。いずれにせよ、どのような３次元再構成方法を採用するにせよ、各投影像を特徴付ける幾何学的パラメータを完全に知ることが最初に必要である。大抵の場合は間接的較正法が必要であり、その場合には較正を自動的に行うことが必要不可欠であろう。
【０００８】
【発明が解決しようとする課題】
本発明の目的は任意の２次元Ｘ線画像形成装置、特に３次元再構成用の２次元Ｘ線画像形成装置で使用可能な自動較正が可能な装置を提案することにある。
【０００９】
【課題を解決するための手段】
本発明は、幾何学的較正装置を備え、この幾何学的較正装置は寸法が分かっている３次元ファントムと、２次元検出器上に投射されたこの３次元ファントムの画像からＸ線画像形成装置の較正係数を算出するための測定手段とを有し、２次元検出器はＸ線管に対応して配置され、検出器およびＸ線管は必要に応じて軸線を中心として一緒に回転することができるようなＸ線画像形成装置において、３次元ファントムがＸ線吸収度がバックグラウンドとは異なるレベルすなわち高レベルまたは低レベルを有する段階的な一連のセル構造を有し、各セル構造はこのセル構造を同定するため順番付けられた特性によって自動的に同定されることを特徴とする装置を提供する。
【００１０】
【作用】
本発明装置ではビーズとその軌跡との間の関連が自動的に分かる。これはビーズを段階または程度に従って一定の順序で分布させることによって、画像形成装置の回転軸線、特にファントムの軸線に沿って測定したビーズの位置が上記順序でのビーズの順番と共に単調に大きくなる（または小さくなる）ようにしたファントムを用いることによって実現できる。
【００１１】
本発明で解決すべき別の問題は較正の回数である。すなわち、この種の回転式画像形成装置では実際に１回の回転で1000またはそれ以上の角度位置を使用することができる。全ての角度位置での較正は例え較正を自動的に実施したとしても極めて長時間かかるので、毎週または毎月定期的に実行することはできない。実際、この較正のために１日が費やされることになる。
【００１２】
本発明方法は、画像形成装置の構造の理論上の寸法が実際の装置に合っていないしても、画像形成装置が完全な円回転またはほぼ完全な円回転をすると仮定することによって上記の別の問題を解決するものである。この方法の考えでは結局のところ画像形成装置の固有パラメータすなわちその回転軸線の向きを与えるパラメータを求めることになる。これらのパラメータは極めて高い精度で知ることができので、これらの固有パラメータから各角度位置に対する較正係数を計算することができる。例えば約 100回の較正でこれらの固有のパラメータを極めて正確に知ることができる。この知識から残りの 900回の「較正」を計算することができる。換言すれば、100 回の較正をしただけで、1000回の全ての較正値が分かるので、得られる時間の利益は極めて大きなものになる。
【００１３】
本発明の別の対象は、下記の段階、
(1) 寸法が分かっているファントムをＸ線管と２次元検出器との間に配置し、
(2) ファントムに対するＸ線管／２次元検出器組立体の回転位置に対して、２次元検出器上に投射されたファントムのＸ線画像でファントムの特性点の画像の座標位置を測定し、
(3) Ｘ線管／２次元検出器組立体の回転位置に対して、ファントムの位置と組み合わされた基準系内でＸ線管の放射線の集束点および２次元検出器からの位置と関連した較正係数を算出し、
(4) Ｘ線管／２次元検出器組立体の所望の各回転位置で、(2) と(3) の２つの段階を反復する
を含む、１つの軸線を中心として回転するＸ線管とそれに対向した２次元検出器とを有するＸ線画像形成装置を幾何学的に較正する方法において、
(5) Ｘ線管／２次元検出器組立体の所定数の回転位置に関する較正係数を処理して、Ｘ線管／２次元検出器組立体の回転位置とは無関係なＸ線画像形成装置に固有なパラメータを抽出し、
(6) この固有のパラメータは、Ｘ線管／２次元検出器組立体の各回転位置での角度値に対してこれらの固有なパラメータの分析関数値を評価することによって、Ｘ線管／２次元検出器組立体の任意の回転位置に関する較正係数から算出される
ことを特徴とする方法にある。
本発明は添付図面を参照した下記実施例の説明から明らかになろう。しかし、本発明が下記実施例に限定されるものではない。
【００１４】
【実施例】
図１は本発明で使用される画像形成装置を示している。この画像形成装置はＸ線源（焦点）Ｓを有するＸ線管１を有し、このＸ線管１に対向してＸ線検出器２が配置されている。Ｘ線検出器２はＸ線撮影で使用される例えば輝度増幅器型のものである。Ｘ線検出器２の検出平面Ｐへの法線はほぼ焦点ＳとＸ線検出器２の中心Ｉとを通っている。これは従来の画像形成では必要であるが、較正の問題には直接関係がない。特に、本発明では相似平行移動型のトモグラフを較正することができる。装置全体はトモ濃度計(tomodensitometer)用の公知の機構またはボウ(bow, 弓形アーム)(図示せず）によって軸線３を中心として回転される。ここでは軸線３は鉛直になっているが、通常は水平である。トモグラフの場合には平行移動と１つの相似中心(hometic center)だけがあり、軸線は必ずしもない。静止装置の場合も同じである。画像形成装置はさらに、Ｘ線検出器２が検出した電気信号をプログラムメモリ６に記憶されたプログラム５で演算処理するプロセッサ４を有している。演算処理で得られた画像または統計結果はモニタ７にディスプレイされる。
【００１５】
本発明では、較正ファントム(calibration phantom) として図２ａ〜図２ｄおよび図３に示したファントムを使用することによって２次元／３次元点を自動的に知るという問題を解決することができる。
本発明の１実施例ではこのファントムは、一連の小さい金属ビーズを支持する支持部材の役目をするプレキシガラスシリンダ８で構成されている。各金属ビーズはこのシリンダ８の壁に等間隔に螺旋状に配置されている。各ビーズは全て同じ直径、例えば１mmの直径を有するが、中央のビーズだけは他のビースより大きい直径、例えば３mmの直径を有している。螺旋パラメータはその射影が検出器上でできる限り大きくなると同時に画像内に含まれるように画像形成装置に合せて作られている。同様に、螺旋のピッチおよび螺旋上のビーズのサンプリング用ピッチは射影内に十分な数のビーズが重ならずに現われるように選択される。ファントムの変形、特にプレキシガラスの老化に起因する変形を防止するために、シリンダ８の両側を蓋９で閉じて剛性を持たせてある。
【００１６】
このファントム８は、例えばＸ線源Ｓと、それから約135cm の距離に配置された直径が37cmの円形Ｘ線検出器とで構成された画像形成装置（検査対象物は焦点Ｓから約100cm の所に配置される）に適している。上記寸法は単なる例で、本発明がこのファントムの寸法に限定されるものではない。図示した実施例では61個のビーズがある。ビースの座標の１つ (例えば縦軸Ｚ) の方向で測定した高さは１つのビーズから次のビーズまで５mmずつ変化し、その精度は５／100 である。螺旋が１回転するまでに20個のビーズが等間隔に分布している。
【００１７】
ファントム８の主要な３つの特徴は以下の通りである。
先ず第１に、ファントム８は較正用に十分な数のビーズを有している。各ビーズは３次元空間に分布しており、それらの射影像は平面Ｐ内に分布している（これは３次元較正の必須条件である）。しかし、分布の均等性は重要ではない。各ビーズはＺ軸に沿った長さの関数として容易に番号付けすることができる。ここでは下記の極めて単純な〔式１〕で与えられる：
【００１８】
【式１】
ｚ＝ｈ・ｎ＋ｚ０
（ここで、ｚ＝Ｚ軸に沿った座標
ｚ０＝０番目のビーズ座標
ｈ＝Ｚ軸に沿ったビーズのサンプリングピッチ
ｎ＝ビーズの番号）
【００１９】
ここで重要なのは順番の関係すなわち順番付けである：
ｚ１＞ｚ２の時はｎ１＞ｎ２
（ここで、ｚ１、ｚ２＝Ｚ軸に沿ったビーズ１、２の座標
ｎ１、ｎ２＝ビーズ番号）
【００２０】
ファントム８がＸ線検出器の前に正しく配置され（すなわちＺ軸がＸ線検出器２の平面Ｐとほぼ平行に配置されており）ており且つビームの円錐度(conicity)があまり大きくなければ、この順番関係が２次元射影でも見出される。
いま、図３でＶをＺ軸の射影に対応する射影平面Ｐの軸線であるとすると、下記の関係式が成り立つ：
ｖ１＞ｖ２ならば、ｎ１＞ｎ２
（ここで、ｖ１、ｖ２＝Ｖ軸に沿ったビーズ１、２の座標ｎ１、ｎ２＝ビーズ番号）
【００２１】
以下では説明を簡単にするため、Ｖ軸は射影画像の鉛直軸とほぼ一致する（すなわち画像の列に平行である）と仮定する。そうではない場合には射影平面を回転させてこの状態にする。
また、下記では説明を簡単にするために、３次元空間の基準系を螺旋と関係付けて定義する。すなわち、螺旋のＺ軸が基準系Ｒの軸の１つを定め、他の２つの軸はＺに垂直である。この基準系ではビーズの座標は極めて簡単に表される。このことは、較正を行った後の基準系は較正時にファントムの軸が占めていた鉛直軸を有するということを意味する。従って、各較正に関する３次元再構成を直接比較することはできない。しかし、本発明の方法では以下で説明するようにこれが可能である。
図３図は、Ｘ線源ＳとＸ線検出器Ｄとの間に配置されたファントム８を示している。このファントム８は、その主軸（Ｚ軸）は平面にほぼ平行で且つＸ線検出器のＶ軸に沿って投射されるように配置されている。
【００２２】
ビーズは絶対的に番号付けできる。螺旋の射影図の頭と底とが万一切れた場合（すなわちファントムの最初のビーズも最後のビーズも見えない場合）でも、中間の特別に大きなビーズを絶対基準にすることができる。この特別な１つのビーズまたは２つの特別なビーズやその他を視界内にくるであろうファントムの位置すなわちファントムの中央に配置するのが好ましい。
【００２３】
以上、本発明の較正方法の原理の理解を助けるために本発明に適したファントムの１実施例を説明したが、本発明では下記特徴を有する全てのファントムを自動較正することができる：
(1) ファントムはセル構造で、記述を単純化するためのＸ線で見えるビーズが３次元空間内に良好に分布されていなければならない。ビーズとバックグランド (環境) とのＸ線コントラストは高くなければならない。ビーズの射影は少なくとも中央ビーズの３次元画像内に含まれていなければならない。ビーズの支持部材は中空なプレキシガラスシリンダである必要はない。
(2) 画像形成装置内のファントムの向きとは無関係に、ビーズに割り当てられた数とＺ軸に沿った３次元空間内および１つの軸、例えばＶ軸に沿った射影内のビーズの座標との間には順番付けすなわち順番関係が存在しなければならない。ファントムの軸は検出器の平面とほぼ平行でなければならない。上記ファントムのように、螺旋状に均等に分布したビースを使用する必要は必ずしもない。
(3) 射影内の最初のビーズと最後のビーズが見えない時でも、少なくとも１つのビーズを他のビーズの中で認識してそれを絶対的に番号付けできなければならない。基準ビーズはその直径以外のファクター、例えば形状、Ｘ線に対する不透明度、不存在、分布ピッチの不均一度、その他任意のファクターで区別できればよい。
【００２４】
例えば、ファントムの別の実施例では、ビーズを円筒状に配置するのではなくプレキシガラス球にに沿って「変形螺旋」状に分布させる例を挙げることができる。以下では、好ましい「螺旋」状ファントムについて説明するが、下記の説明は本発明の範囲に入る他の方法にも適用することができる。
【００２５】
先ず、画像形成装置の較正するための操作順序を説明する。
１．上記ファントム８を、その主軸（Ｚ軸）が検出器の平面とほぼ平行となるように配置する。例えばファントムを血管造影装置の検査台上に載置し、Ｘ線管および検出器を支持したボウ (弓形アーム) を検査台に対して垂直な平面内に配置する。この位置でビーズの射影像から上記順番関係を確認する。
２．較正すべき入射条件下でファントムのデジタル化射影像を取る。
３．デジタル化射影像を画像プロセッサ４へ送る。
４．デジタル化射影像の幾何学的歪みを例えば1988年７月１日出願のフランス国特許出願第 2,633,793号に記載の方法で補正する。
５．補正済み画像Ｉを以下で説明するビーズの検出、位置決めおよび自動ラベル付け用アルゴリズムで処理する。
６．３次元基準系Ｒ内のビーズの座標と平面Ｐの２次元基準系内のビーズの射影の座標とを展開して、射影の幾何学的パラメータすなわち較正係数を算出する。以下、共役勾配(conjugated gradient) 法によって誤差を最少にするアルゴリズムを説明する。
【００２６】
以下、上記の較正ファントムの投射画像（図４ａ、画像Ｉ）から機械に固有な幾何学的パラメータを得るための較正アルゴリズムの１実施例を記載する。以下の説明では各種の処理演算した結果の画像を説明する：
１．先ず最初に画像Ｉをセグメント化する。ここでの目的は、画像Ｉからビーズを支持するプレキシガラスシリンダの投射信号を除去することである。そのためには数学的形態の道具 (すなわち、画像Ｉ内のビーズの寸法より大きい寸法の円形構造要素を備えたトップハット変換(top-hat transformation)) を用いる。得られる信号Ｓ１はビーズの吸収に対応する信号以外は何も含まない。数学的形態の道具は公知で、多くの研究発表がある。
２．シリンダのエッジ部から来る人工的な残留信号を画像Ｓ１から除く必要がある。そのために再び数学的形態の変換を用いて、画像内の鉛直構造を除去する。この結果得られる画像は図４ｂにネガで示した画像Ｓ２である。この画像Ｓ２にはノイズが除かれたバックグランドをゼロとした全てのビーズ（正の信号）が含まれている。
３．次に、画像Ｓ２をセグメント化する。ここでの目的はファントムの大きな中央ビーズから来る信号のみを取ることにある。そのため小さいビーズの寸法より大きく且つ大きいビーズの寸法より小さいものを画像Ｓ２から除く。従って、大きいビーズのみを含む画像Ｓ３が得られる。
【００２７】
４．次に、画像Ｓ２の全てのビーズを検出する。ここでは、画像内の全てのビーズの射影を最も近い画素 (ピクセル) に位置決めする。画像の０行（検出器の頭部）から演算を開始される。画像全体を上から下へ横断して演算する。画像内の正の信号は、その信号の所定閾値１より高い時にはビーズであることが明らかにされる。番号０は最初に出会うビーズに与えられ、画像の底ままで遭遇するビーズに番号をインクリメントして与える。上記の螺旋の特性からこの順番が正しいことは確かである。
５．次に、ビーズの正確な座標を計算する。ここでの目的は画素以下の精度(sub -pixcel precision)で画像内の各ビーズの座標を知ることにある。画像Ｓ２内の各ビーズについて、検出されたビーズの周囲の画素の組の重心をグレーレベルで計算し、このグレーレベルが第２の閾値 (閾値Ｔ２)(閾値Ｔ２＜閾値Ｔ１）より高いか否かを調べる。この計算は例えば上記フランス国特許出願第 2,633,793号に記載されている。その原理は、ビーズは十分に小さいのでＸ線を全て吸収することはないという事実に基づいている。すなわち、吸収飽和がないと、画素が得られる。
６．次に、画像Ｓ３内で大きいビーズを検出する。この座標は最も近い画素に対して分かる。この座標と画像Ｓ２で検出された全てのビーズの座標とを比較することによって、全てのビーズリストの中で最大のビーズを知ることができる。
７．投影像内の最大ビーズが分かったので、全てのビーズを絶対的に番号付けることができ、従って、各ビーズについてその番号とその３次元座標（分かっている）との間の正確な対応関係およびその射影の２次元座標との間の正確な対応関係が分かる。
【００２８】
各ビーズについて３次元座標とその投影像の２次元座標が分かれば、円錐形投影画像の９つのパラメータに対して選択した任意の値に起因する投影誤差を計算することができる。この誤差は任意の値を用いて計算した位置と較正ファントムの画像での位置との間の差であるので、検出された全てのビーズについて投影誤差を最小にするようなパラメータを選択することによって、各視野すなわち較正すべき回転位置に関係付けられたパラメータまたは較正係数を計算することができる。
【００２９】
パラメータ計算アルゴリズムとしては、平均平方エラーに対して共役勾配最小化(conjugate gradient minimization) 法を用いることができる。これは標準的な方法で、パラメータの関数として平均平方エラーと一次導関数とを用いて計算する。このアルゴリズムはプレス(W. H. Press) 、フランネリ(B. P. Flannery)、テウコルスキィ(S. A. Teukolsky) およびベッタリング(W. T. Vetterling)の『Ｃでの計算方法(Numerical Recipes in C)』ケンブリッジユニヴァーシティプレス(Cambridge University Press)、1988年、第10・６章、第 317頁以降に記載されている。この文献の同じ頁に記載されているシンプレックス(Simplex) やニュートン(Newton)とよばれる別のアルゴリズムを使用することもできる。シンプレックスとよばれる第１のアルゴリズムはパラメータ化の任意のコヒーレント系に適用できるという利点があるが、容易に収束しない。共役勾配アルゴリズムはより速く収束するがパラメータの解析表現ができるようなパラメータ化する必要がある。これは上記パラメータ化の場合である。このアルゴリズムに基づく計算は1992年11月２〜４日にランヌ(Rennes)で開催された第14回ＩＥＥＥＥＭＢＳ年次総会で発表されたルージュ(Anne Rougee) 、ピカール(Catherine Picard)、ポンシュ(Cyril Ponchut) およびトルッセ(Yves Trousset) の『医学における３次元画像処理(3D Advanced Imaging Processing In Medicine)』に記載されている。
【００３０】
このアルゴリズムのステップは以下の通り：
１．初期化：パラメータを実際の値で初期化する。それによって解に近い領域で最小化を開始することができる。
２．標準化：この初期値の近くで、パラメータの空間の各軸上での１単位分の移動に対応する基準(criterion) 値を計算し、次に、各軸上の目盛りを変えることによってこの空間内の変数を変えて、各軸上で１単位分の移動をした時に基準に与える作用を同じにする。
３．最小化：初期値から開始して、標準化段階での目盛りの変化を考えながら共役勾配法を繰り返して、検出された全てのビーズの平均平方誤差を最小化する。各反復時に基準の導関数の関数としてパラメータの新しい予測値を計算する。現在の予測値と次の予測値との間の相対差が所定閾値より小さくなった時に計算を止める。
【００３１】
「３次元再構成」をする装置とは、軸線を中心として回転することによって入射角を変えて対象物体をＸ線撮影することができる装置である。こうして得られたＸ線撮影像を３次元情報の再構成で用いる。この３次元再構成装置は例えばＸ線管／検出器組立体が剛性リング上を回転するような構造にすることができる。この場合には回転軸は一定で、時間が経過しても全く変化しない。あるいは、Ｘ線管／検出器組立体を傾斜台上に支持させる構造にすることもできる。この場合回転軸は別の方法で空間内に位置決めされる。
【００３２】
本発明では、Ｘ線管／検出器組立体の各視野すなわち各回転位置で較正係数とよばれるパラメータを計算する代わりに回転装置の固有のパラメータを計算することができる。この固有のパラメータはＸ線管／検出器組立体の回転装置の回転位置の解析関数である。位置に関する較正パラメータを得るには、この回転装置の所定の向きの値に対する固有パラメータ値を計算すればよい。従来の装置では空間内での向きを表す測定値（大抵は電気信号）を出しているので、この測定値を用いて向きに関する較正係数を得ることができる。すなわち測定値を固有パラメータの演算表現にする。
【００３３】
３次元再構成装置の全体的幾何形状を凝縮した形で記述する上記の固有パラメータの算出には以下の手順を用いることができる。すなわち、先ず第１段階で各視野毎に円錐形投影画像のパラメータを算出し、次に第２段階でこれらパラメータに基づいて固有のパラメータを計算する。
【００３４】
第１段階
３次元再構成装置には回転中に各視野毎に較正係数を算出する自動較正方法が特に適している。
螺旋体を３次元再構成装置の回転軸とほぼ同軸なＺ軸を有する検査台上に設置する。回転中に３次元再構成装置が占める各視野で螺旋体の投影像を撮る。得られた全ての映像を上記の方法で処理する。
本発明で提案するパラメータ化法（円筒形の３次元座標）はＺ軸を中心とした円形軌跡に沿ったパラメータの算出に適している。初期化と標準化段階とを最初の視野だけで行う。次に、最初の視野で得られたパラメータ値を用いて次の視野のパラメータ算出の初期化を行い、以下、最後の視野までこの操作を反復する。これによって計算を加速することができる。
こうして各視野に対する幾何学的パラメータが得られる。これらの幾何学的パラメータは３次元再構成装置の磁気ディスクに記憶され、各３次元再構成段階で使用される。
【００３５】
第２段階
本発明では、各視野ごとの較正係数を知ることによって３次元再構成装置全体の幾何学的配置を凝縮した形で記述する固有パラメータを知ることができる。
今、Ｎを較正済みの視野の数とする (例えばＮ＝100)とする。各視野はその近傍と無関係に較正しなければならないので、画像(aquisition)の幾何学的配置を記述するためには９Ｎ個のパラメータが存在する。しかし、Ｘ線管と検出器とが結合されて固定軸線を中心として回転する系の場合にはこれら９Ｎ個のパラメータは冗長過ぎる。３次元再構成装置の画像の幾何学的配置は一定数のパラメータでパラメータ化することができ、下記パラメータが好ましい：
(1) ３次元再構成装置のアイソセンターＯの座標（３つのパラメータ）
(2) アイソセンターＯを通る回転軸線Ｚ’の向き (２つの角度で与えられる)(２つのパラメータ)
(3) 射影面上の点ＩとＸ線源Ｓの軌跡の半径ｒ_I、ｒ_s（２つのパラメータ）
(4) Ｘ線源Ｓおよび射影面上の点Ｉの軌跡の距離の差 (回転軸Ｚ’方向の)(１つのパラメータ)
(5) Ｚ’に対するＳの極角度とＺ’に対するＩの極角度との間の角度差ε（１つのパラメータ）
(6) 検出器のＺ’軸に対する方向角度である検出器の傾斜α₁、α₂（２つのパラメータ）
(7) 画像のカラムと画像中でのＺ’軸の射影像との間の角度α₃（１つのパラメータ）
【００３６】
全部で12のパラメータがある。これ以外のパラメータを選択することも可能である。その選択は上記の場合と同じ理由すなわち計算の単純化のために行うのが好ましい。この場合、これらの固有パラメータ値が表されている基準系はファントムの位置とは無関係であることは理解できよう。従って、１つの較正から別の較正に交換できる。これらの12個のパラメータを９Ｎ個の初期パラメータの関数として用いるには例えば下記の方法を使うことができる：
【００３７】
１．回転軸Ｚ’の向きを求める
Ｉの軌跡を形成するＮ個の点Ｉ_n（ｎ−１、Ｎ）を考える。この軌跡は理論上は平らであり、Ｚ’に直角であるので、最小平方で表した点Ｉ_nの組に最も近い平面を求め、て、それからこの平面の法線であるＺ’の方向を求める。
点Ｉ_nの座標を（Ｘ_n、Ｙ_n、Ｚ_n）とする。点Ｉ_nの組に「最も近い平面」の式ｚ＝ａｘ＋ｂｙ＋ｃを求める。平面ｚ＝ａｘ＋ｂｙ＋ｃと点Ｉ_nの組との間の距離をｄ（ａ、ｂ、ｃ）とする。この距離ｄ（ａ、ｂ、ｃ）を最小にする係数（ａ、ｂ、ｃ）が求める。ｄ（ａ、ｂ、ｃ）を最小にする３つの数の組（ａ、ｂ、ｃ）は下記：
ｄ（ａ、ｂ、ｃ）＝Σ（ａｘ_n＋ｂｙ_n＋ｃ−ｚ_n）²
と定義すると、一次回帰によって極めて容易に見つけることができる。
【００３８】
２．アイソセンターの座標とＩの軌跡の半径を求める
アイソセンターはＩの軌跡およびＳの軌跡の両方の中心である。従って、この中心は、Ｉの軌跡を解析すれは求められる。
（注）この段階ではＩの軌跡のパラメータ（中心と半径）を求め、次の段階でＳの軌跡の半径を求める。論理的には全く逆に行ってもよい（Ｓの軌跡の中心と半径を求め、次いでそれからＩの半径と軌跡を演繹する）が、この方法に利点はない。すなわち、点Ｓは点Ｉより回転軸から離れているので、Ｓの位置の情報要素はＩの位置の情報要素より誤差が多い。従って、本発明方法では「確実な」情報要素（Ｉの情報要素）を基礎として用い、不確実なＳについての情報要素は用いない。また、ＩおよびＳの軌跡の半径と中心を同時に求める方法も同じ理由で採用しない。
【００３９】
アイソセンターの座標すなわちＩの軌跡の中心座標を見つける方法
Ｚ’に直交し且つＺ’に沿った点Ｉ_nの平均ディメンジョンに等しいディメンジョンを有する平面πを定義する。全ての点Ｉ_nは点Ｉ_nでこの平面πに直交するように投射される。Ｉの軌跡を円と仮定すると、点Ｉ’_nは原則として円上にくる。点Ｉ’_nの組に最も近い円を最小平方で求めると、以下のものが分かる：
▲１▼ Ｉの軌跡の半径ｒ₁（すなわち求める円の半径）
▲２▼ アイソセンターの座標 (πと回転軸との交点として定義される)(求める円の中心)
【００４０】
最も近い円を定義し、発見する方法
中心がＡで半径がｒの円をＣ（Ａ、ｒ）とする。先ず、円Ｃ（Ａ、ｒ）と点Ｉ’_nの組との間の距離ｄ（Ａ、ｒ）を求める。これはいくつかの方法で得られる。
方法１ (図５(a))
Ｉ’_nの座標を（Ｘ’_n、Ｙ’_n）とする。下記式を定義する：
ｄ（Ａ、ｒ）＝ΣＩ’_nＪ’_n
（ここで、点Ｊ’_nは円Ｃ（Ａ、ｒ）ｐと直線ＡＩ’_nとの交点）
【００４１】
方法２ (図５(b))
３次元再構成装置が等間隔な角度位置すなわちδΘ毎（δΘは一定）にＸ線を放射するとすると、この知識を円のパラメータを求めるのに用いることができる。
ｄ（Ａ、ｒ）はｄ（Ａ、ｒ）＝ΣＩ’_nＪ’_nで定義される。
点Ｊ’_nは基準方向に対して極角度：Θ＝Θ_o＋（ｎ−１）＊δΘを有する円Ｃ（Ａ、ｒ）の点として定義される。
【００４２】
円と点群との間の距離の測定に何れの方法（方法１、方法２）を用いた場合でも、後の手順は同じであり、任意の最小化法（例えば共役勾配法）を用いて、円への距離を最小にする円のパラメータＡとｒ（方法２ではパラメータａ、ｒおよびΘ₀）を求める。これらのパラメータが分かればこの段階は終了する。すなわち、
▲１▼ Ｊの軌跡の半径が分かる。すなわち、それはｒである。
▲２▼ 回転軸Ｚ’の向き前記段階で求めてあり、この軸Ｚ’がＡを通ることは分かっているので、回転軸Ｚ’が完全に分かる。
【００４３】
３．Ｓの軌跡の半径を求める
定義から、Ｓの軌跡の半径は前の段階で求めたＩの軌跡の半径と同じであるので、Ｓの軌跡の中心が分かる。従って、半径を求めることができる。平面π上の点ｓの投影である点Ｓ’_nを考えると、Ｓの軌跡の半径はアイソセンターＯと点Ｓ’_nとの間の距離の平均として定義される：
Ｉの半径＝１／Ｎ＊ΣＯＳ’_n
【００４４】
４．Ｘ線源の軌跡とＩの軌跡との間の軸Ｚ’方向の距離差を求める
これは、点Ｓ_nの軸Ｚ’方向の平均距離である（図５ｃを参照）。この距離差は平面πまでの距離から計算される。
【００４５】
５．角度εを求める
この演算では平面πの投影像Ｉ’_nとＳ’_nとを用いる。Ｉ’_nはＳ’_nの極角度（０に対する）を各ｎに対して計算し、その差をε_nとする（図５ (d) 参照) このεは値ε_nの平均として算出される。
【００４６】
６．検出器の傾斜α ₁ とα ₂ とを求める
螺旋の軸線に関連付けた基準系で、各視野で検出器に対して直角な座標を計算する。次に、螺旋の軸線に関連付けた基準系とアイソセンターおよび回転軸Ｚ’に関連付けた基準系との間で基準系を変化させる式を用いて、アイソセンターおよび回転軸Ｚ’に関連付けた基準系の検出器に直角な座標を求める。上記の基準系を変化させる式は螺旋に関連付けた基準系の知識から演繹される。アイソセンターは上記の円の中心であり分かっているので、Ｚ’軸の向きも分かっている。この基準系の他の２つの軸は１回転分以外は固定される。この固定のために、新しい基準系で第１の視野に３次元再構成装置の好ましい位置に関連付けた角度Θ₀を割り当てるのが好ましい。この新しい基準系の座標は２つの角度α₁、α₂のみで決まる。この２つの角度α_1n、α_2nを各ｎについて計算する。α₁はα_1nの平均値として定義される。α _2nについても同じ。
【００４７】
角度α ₃ を求める
螺旋の軸に関する基準系で、各視野について検出器の列に平行な単位ベクトルの座礁を計算する。上記の基準変化式を用いて、アイソセンターと回転軸Ｚ’に関する基準系での座標を演繹する。この新しい基準系ではこれらの座標は角度α₃のみに依存する。各ｎについてこの角度（α_3nとよぶ）を計算する。α₃はα_3nの平均値として定義される。
【００４８】
第３段階
以上のようにして、３次元再構成装置の幾何学配列を決定する12個の固有パラメータが分かる。これらの12個の固有のパラメータの知識に基づいて、各視野ｎについて、この視野を定義する９つの局所較正パラメータを計算するのは極めて容易である。この計算は３次元再構成装置の固有基準系すなわちアイソセンターと回転軸とに関連付けた基準系で行った計算と同様に簡単である。視野と関係付けた９つの局所パラメータは、レンヌ会議の前記文献に記載のように (ｒｓ_n、ａｓ_n、ｚｓ_n、theta _n、ｐｈｉ_n、ｐｓｉ_n、ｒｉ_n、ａｉ_n、ｚｉ_n）とよぶ。これから下記〔式２〕が得られる：
【００４９】
【式２】
ｒｓ_n＝Ｘ線源の軌跡の半径（＝定数）
ａｓ_n＝Θ₀＋180 ＋（ｎ−１）＊δΘ
ｚｓ_n＝上記で計算した距離の差（＝定数）
theta _n＝Θ₀＋（ｎ−１）＊δΘ＋α１
ｐｈｉ_n＝α₂
ｐｓｉ_n＝α₃
ｒｉ_n＝Ｓの軌跡の半径（＝定数）
ａｉ_n＝Θ₀＋（ｎ−１）＊δΘ
ｚｉ_n＝０（＝定数）
【００５０】
すなわち、これらの９つの較正係数またはパラメータはΘ₀＋（ｎ−１）δΘから演繹される入射角から決定される。これらは高い精度で分かっているので、入射角を正確に測定することによって、較正係数を極めて正確に計算することができる。入射角それ自体はインデックスｎ（３次元再構成装置の外周を所定の予め決められた数で割った数で、分かっている）から計算される。
【００５１】
上記のＮ個の較正された視野に関する較正係数を、固有係数に従った解析的演算法に基づいて、再計算するのが好ましいということは理解できよう。それによって、残りの視野（他の900 個の視野) をＮ(100) 個の視野の較正係数と同じ精度で調べることができ、全体の結果が改良される。
【００５２】
別の可能の方法は、螺旋体の全ての射影像で検出されたビーズの全座標に基づいて３次元再構成装置の全パラメータを直接測定する方法であるが、この方法は時間がかかる方法である。
【図面の簡単な説明】
【図１】本発明を実施するのに使用される画像を別の基準系に作る画像形成装置の概念図。
【図２】本発明のファントムを示す図で、 (a)はファントムの側面図、(b) はその正面図、(c) は投影図でのビーズの軌跡、(d) はビーズの取付け方法 (位置決め方法) を示す詳細拡大図。
【図３】画像形成装置内でのファントムの配置を示す図。
【図４】 (a)および(b) はファントムの特定の形状の投影図を処理した後の画像を示す図。
【図５】 (a)〜(d) は本発明方法の理解を助けるための幾何学的な解説図。
【符号の説明】
１Ｘ線管２検出器
４プロセッサ６プログラムメモリ
７モニター８ファントム

Claims

(1) 画像形成装置の全て回転位置で各セル構造の順番を順次自動的に同定可能な予め寸法が分かっている一連のセル構造からなり、バックグラウンドとはＸ線吸収度が異なる一つの軸線を有する３次元ファントムを、Ｘ線管と２次元検出器との間に３次元ファントムの軸線が３次元Ｘ線画像形成用装置の軸線に対して平行となるように配置し、
(2) ファントムに対するＸ線管／２次元検出器組立体の回転位置に対して、２次元検出器上に投射されたファントムのＸ線画像でファントムの特性点の画像の座標位置を測定し、
(3) Ｘ線管／２次元検出器組立体の回転位置に対して、ファントムの位置と組み合わされた基準系内で、Ｘ線管の放射線の集束点および２次元検出器からの位置に関連した較正係数を算出し、
(4) Ｘ線管／２次元検出器組立体の所望の各回転位置で、上記(2)と(3)の２つの段階を反復する、
ことを含む、１つの軸線を中心として回転するＸ線管とそれに対向した２次元検出器とを有するＸ線画像形成装置を幾何学的に較正する方法において、
(5) Ｘ線管／２次元検出器組立体の所定数の回転位置に関する較正係数を用いてＸ線管／２次元検出器組立体の回転位置とは無関係なＸ線画像形成装置に固有なパラメータを抽出し、
(6) この固有のパラメータはＸ線管／２次元検出器組立体の各回転位置での角度値に対してこれらの固有なパラメータの分析関数値を評価することによってＸ線管／２次元検出器組立体の任意の回転位置に関する較正係数から算出する、
ことを特徴とする方法。
(a) ２次元検出器上のファントムの特性点の画像である理論上の位置の座標を任意の較正係数を用いてシミュレートし、
(b) この理論上の位置の座標と各ファントムの画像で測定された位置の座標とを比較し、
(c) その結果に応じて上記の任意の較正係数を変え、
(d) 両者の差が無視できるようになるまで (a)段階と (b)段階とを反復し、この場合、変更した最後の較正係数をＸ線管／２次元検出器組立体のこの場合の回転位置の較正係数として選択する、
ことによって、Ｘ線管／２次元検出器組立体の所定数の回転位置に関する較正係数を計算する請求項１に記載の方法。
(a) Ｘ線画像形成装置の物理的構造から演繹された較正係数を用いてシミュレーションを行い、
(b) この理論上の位置の座標とファントムの画像で測定された位置の座標とを比較し、
(c) 較正係数を変え、
(d) 両者の差が無視できるようになるまで(a)段階と(b)段階とを反復し、この場合、変更した最後の較正係数をＸ線管／２次元検出器組立体のこの場合の回転位置の較正係数として選択する、
ことによって、Ｘ線管／２次元検出器組立体の所定数の回転位置に関する較正係数を計算する請求項１に記載の方法。
前回の回転位置に対してずれているＸ線管／２次元検出器組立体の新しい所望の回転位置に関する較正係数を演繹するために、この新しい回転位置に関する任意の較正係数を使用する代わりにその前の回転位置で得られた最後に修正された較正係数を用いる請求項２または３に記載の方法。
較正係数が円筒形基準系で記載された固有パラメータの解析関数で表される請求項１に記載の方法。
固有のパラメータを、画像形成装置の各向きでの所定数のファントム画像を統計処理して得る請求項１〜５のいずれか一項に記載の方法。
固有パラメータに基づいて、ファントムの特性点の画像位置の座標を測定し且つ較正係数を演繹したＸ線管／２次元検出器組立体に関する較正係数を演繹する請求項１〜６のいずれか一項に記載の方法。
下記の（ 1 ）〜（ 6 ）とを有することを特徴とする請求項１〜７のいずれか一項に記載の方法を用いて３次元Ｘ線画像形成装置の回転系を幾何学的に較正する自動較正システム：
（ 1 ）１つの軸線を中心として回転するＸ線管とそれに対向した２次元検出器とを有するＸ線画像形成装置、
（ 2 ）画像形成装置の全て回転位置で各セル構造の順番を順次自動的に同定可能な予め寸法が分かっている一連のセル構造からなり、バックグラウンドとはＸ線吸収度が異なる、Ｘ線管と２次元検出器との間に３次元Ｘ線画像形成用装置の軸線に対して平行となるように配置される３次元ファントム、
（ 3 ）ファントムに対するＸ線管／２次元検出器組立体の回転位置に対して、２次元検出器上に投射されたファントムのＸ線画像でファントムの特性点の画像の座標位置を測定する測定手段、
（ 4 ）Ｘ線管／２次元検出器組立体の回転位置に対して、ファントムの位置と組み合わされた基準系内で、Ｘ線管の放射線の集束点および２次元検出器からの位置に関連した較正係数を算出し、Ｘ線管／２次元検出器組立体の所望の各回転位置で、上記（ 2 ）と（ 3 ）の２つの段階を反復するする第１の計算手段、
（ 5 ）Ｘ線管／２次元検出器組立体の所定数の回転位置に関する較正係数を用いてＸ線管／２次元検出器組立体の回転位置とは無関係なＸ線画像形成装置に固有なパラメータを抽出する抽出手段、および
（ 6 ）Ｘ線管／２次元検出器組立体の各回転位置での角度値に対してこれらの固有なパラメータの分析関数値を評価することによってＸ線管／２次元検出器組立体の任意の回転位置に関する較正係数から上記の固有のパラメータを算出する第２の計算手段。
３次元ファントムが１つの軸線と、支持部材と、この支持部材上に分布した複数のセル構造とから成り、各セル構造の順番はファントムの上記軸線に沿って測定した時に単調に増加し、この順番はファントムの射影画像でも同じ順番である請求項８に記載の自動較正システム。
支持部材が円筒形シリンダで、この円筒形シリンダの母線は上記２次元検出器の軸線と平行であり、一連のセル構造は支持部材上で螺旋状に並んでいる請求項８または９に記載の自動較正システム。
円筒形の支持部材の端部に変形を防ぐ固定板を有する請求項９に記載の自動較正システム。
３次元ファントムの少なくとも１つの特定セル構造に対して各セルを参照することで各セル構造の順番を順次自動的に同定する請求項８〜１０のいずれか一項に記載の自動較正システム。
上記の特定セル構造が他のセル構造とは異なる寸法または放射線吸収能を有し、セル構造群の中心に配置されている請求項１２に記載の自動較正システム。
各セル構造がＸ線の通過を完全に防ぐには不十分な放射線吸収能を有する請求項８〜１３のいずれか一項に記載の自動較正システム。
各セル構造が支持部材に形成された凹部の底に配置された金属ビーズである請求項８〜１４のいずれか一項に記載の自動較正システム。