WO2016147314A1

WO2016147314A1 - コンピュータ断層撮影方法及び装置

Info

Publication number: WO2016147314A1
Application number: PCT/JP2015/057884
Authority: WO
Inventors: 明男米山; 馬場　理香
Original assignee: 株式会社日立製作所
Priority date: 2015-03-17
Filing date: 2015-03-17
Publication date: 2016-09-22

Abstract

Ｘ線ＣＴにおいてアーチファクトのない断面像を再生する　被写体とX線撮像系を相対的に回転して異なる投影角度で得られた画像データから被写体の断面像を逆投影によって再構成するコンピュータ断層撮影について、逆投影において投影方向に一様に値を加算していくのではなく、他の投影角度（基本的には直交する方向）で得られた像データを利用して加算領域を設定し，その領域だけ値の加算を行うようにする。

Description

コンピュータ断層撮影方法及び装置

　本発明はＸ線撮像技術に係わり、特に、物体の内部を非破壊で高い定量性で検査するのに適したＸ線撮像技術に関する。

　Ｘ線ＣＴ（Computed Tomography）は、Ｘ線源に対して被写体を相対的に回転し，異なる角度から取得した透過像を用いて，計算により被写体の断面像を得る方法である。Ｘ線の高い透過能を利用して非破壊で被写体内部を三次元的に観察できることから，医療診断や製品の検査には不可欠な技術となっている。

　測定した透過像データから断面像を求める計算は「再構成」と呼ばれ，一般に以下に示す「逆投影」演算により求める。

　図１に再構成と逆投影の関係を示す。

　図１（ａ）はＣＴの測定時の操作概念であり、単純化するため、被写体が中心（中心の座標を（０，０）とする）にしか存在しないものとする。測定系においては，各投影角度θにおける投影像データはその方向への強度の積算となるので、例えば投影角度θが、白矢印で示す９０度（横方向）や、黒矢印で示す０度（縦方向）では、図１（ａ）のθ=90度の投影データ１０１、θ=0度の投影データ１０２のような、中心の値（１）のデータとなる。

　一方、再構成で用いる「逆投影」では、各投影方向で得られた像データ１０１，１０２を仮想的な被写体のデータ配列（２次元あるいは３次元）に、順次加算していくことになる。図１の例では、投影角度θ=０で取得されたデータ１０２を０度方向から、９０度で取得されたデータ１０１を９０度方向から、各画素に投影方向に沿って、図１（ｂ）の矢印の方向に加算していく。このことを各角度で取得した像データについて行えば，非破壊で被写体の断面像を取得することができる。以上が，ＣＴによる断面像の再構成の原理である。

　図１（ｂ）からわかるように、投影角度に沿って各画素に全データを加算していくために、本来は被写体のないデータ「０」の領域も値（図１（ｂ）の例では「１」）を持つことになってしまい、逆投影だけでは正確にオリジナルのデータを再構成することができない。

　この問題を解決するために、従来のCTでは逆投影の前に、投影像データにフィルターを作用（コンボリューションして）し，逆投影により全てのデータを加算すると、被写体のないところの値は０になるように各像データの値を変換している。この方法をフィルタード・バックプロジェクション（FBP）と呼ぶ。FBPについては、例えば特許文献１に概要が記載されている。

特開２０１３－１９０３３３号公報

　従来のフィルタード・バックプロジェクション（FBP）では、被写体の形状や密度によっては，被写体のないところの値が０にならないこともあり，ライン状の疑似像（アーチファクト）が生じるという問題があった。また，作用させるフィルターの形状（一般には関数を利用するので，関数式やパラメータ）によって，再生されるＣＴ値が異なるため，異なるフィルターで取得した断面像におけるＣＴ値の定量的な比較は難しいという問題があった。

　上記課題を解決するために本発明の一側面では，逆投影において投影方向に一様に値を加算していくのではなく、他の投影角度（例えば直交する方向）で得られた像データを利用して加算処理を制御するようにする。これにより、被写体が存在しない本来は０である領域の値を０のままとすることができ，アーチファクトのない断面像を再生することができる。すなわち、他の投影角度で得られた像データを用いると、被写体が存在する領域をある程度特定することができるので、この情報を用いて逆投影の処理を行う。

　本発明の他の一側面は、Ｘ線源とＸ線画像検出器と演算部を備え、Ｘ線源から被写体に対して複数の異なる方向からＸ線を照射し、Ｘ線画像検出器で前記複数の異なる方向から照射されたＸ線によって得られる複数の異なる投影角度で得られた投影像を取得し、投影像の画像データを演算部によって処理し、被写体の断面像を逆投影によって再構成する、コンピュータ断層撮影装置である。演算部は、コンピュータの処理装置によってソフトウェアを実行することで実現することができる。また、専用のハードウェアで構成してもよい。演算部は、第１の投影角度で得られた画像データからマスクデータを作成する第１の処理部と、第２の投影角度で得られた画像データをデータ配列に加算して断面像を再構成する際に、画像データをマスクデータで重み付けして前記断面像の一部を形成する処理を行う第２の処理部を有する。そして、第１の投影角度および第２の投影角度を変化させ、断面像の一部を形成する処理を繰り返すことにより、断面像を再構成する。

　より具体的な例では、第１の投影角度と第２の投影角度を直交する角度に設定する。また、ある具体例では、第１の処理部は、第１の投影角度で得られた画像データによって、被写体が存在しないと判定される座標が０の値となるマスクデータを作成し、第２の処理部は、第２の投影角度で得られた画像データとマスクデータを積算する処理を行う。また、第２の処理部では、第２の投影角度で得られた画像データもしくはマスクデータの少なくとも一つを座標変換する処理を加えてもよい。

　本発明の他の側面は、被写体とＸ線撮像系を相対的に回転して、異なる投影角度で得られた画像データから、被写体の断面像を逆投影加算によって再構成するコンピュータ断層撮影方法において、第１の投影方向から取得した画像データを逆投影加算する際に、これと異なる第２の投影方向から取得した画像データを用いて、再構成する断面像の加算する領域を制限することを特徴とする。

　本発明の別の側面は、被写体とＸ線撮像系を相対的に回転して異なる投影角度で得られた画像データから繰り返し計算によって被写体の断面像を再構成するコンピュータ断層撮影において、被写体の断面像を再構成するための逆等投影計算が、異なる投影方向から取得したデータ上で被写体が存在する領域だけ加算し，存在しない領域は加算しない制限を設けた逆投影計算であることを特徴とする。

　上記のより具体化された一例として計算機等の演算装置で実行する処理手順を以下に示す。なお，ここでＸ線撮像装置によって取得した各投影角度θにおける投影データI(θ)と，被写体のない背景データI_BKと，Ｘ線を照射しないで取得したダークデータI_Dとする。また，被写体は360°回転して投影データを取得したとする。

　（１）各投影データI(θ)について，
　　　I_1(θ)=(I(θ)－I_D)/(I_BK－I_D)
を計算し，規格化された投影データI_1 (θ)を求める。（θは０から３６０°）
　（２）Ｉ_1（θ）について，
　　　I_2(θ)=－ln(I_1(θ))
を計算し，線吸収係数データI_2(θ)を求める。（θは０から３６０°）
　（３）投影角度θに対して直交する角度（θ＋９０°）の線吸収係数データI_2(θ＋９０)について，各画素の位置xにおける値を元にマスクデータM(x”)を
　　　M(x”) = 1 if I_2(θ＋90，x”) ＜＞ 0
　　　　　　　　0 if I_2(θ＋90，x”) = 0
により計算する。すなわち，被写体のある位置のマスクデータを１，無いところを０とする。

　（４）仮想的な再構成像データS(x’, y’)を，I_2(θ,x)とM(x)を用いて計算する。具体的には各画素(x’, y’)について
　　　x = x’ sin θ ＋ y’cos θ
　　　x” = x’ sin (θ＋90°) ＋ y’cos （θ＋90°）
　を計算し，この値に対応したI_2（θ,x）とM(x”)を用いて
　　　S(x, y) = S(x, y) ＋ M(x”)×f(I_2(θ))
　となる加算演算を行う。なお，ここでf()はShepp－Logan関数など従来のフィルター関数のコンボリューションである。

　（５）上記（３）と（４）の過程について，θが０°から３６０°まで行う。なお，１８０°以降のM(x”)の計算は，θ－90°のデータI_2(θ－９０°)を用いれば良い。

　以上の計算によって，被写体のある領域だけ値を持ち，他の領域は０となる。すなわち，アーチファクト（疑似像）のない，断面像を再生することができる。

　本発明によれば、疑似像のない，正確なＣＴ値の断面像を求めることができる。

CTと逆投影の原理図Ｘ線撮像装置の構成図Ｘ線ＣＴの測定フローチャート図本発明の実施例のフローチャート図本発明の実施例のフローチャート図本発明の実施例１の適用例を示す原理図本発明の実施例２の適用例を示す原理図本発明の実施例３における逐次計算のフローチャート図本発明の実施例４における断面像と投影像の関係を示す原理図。本発明の実施例における照射Ｘ線と被写体の回転角度の関係を示す原理図。ラミノグラフィーにおける円盤状の被写体の観察例を示す原理図

　以下、図面を用いて本発明の実施例について説明する。以下に示す図において、同じ機能を有する部分には同じ符号を付し、重複する説明を省略することがある。

　本発明は以下に示す実施の形態の記載内容に限定して解釈されるものではない。本発明の思想ないし趣旨から逸脱しない範囲で、その具体的構成を変更し得ることは当業者であれば容易に理解される。

　本明細書等における「第１」、「第２」、「第３」などの表記は、構成要素を識別するために付するものであり、必ずしも、数または順序を限定するものではない。また、構成要素の識別のための番号は文脈毎に用いられ、一つの文脈で用いた番号が、他の文脈で必ずしも同一の構成を示すとは限らない。また、ある番号で識別された構成要素が、他の番号で識別された構成要素の機能を兼ねることを妨げるものではない。

　図面等において示す各構成の位置、大きさ、形状、範囲などは、発明の理解を容易にするため、実際の位置、大きさ、形状、範囲などを表していない場合がある。このため、本発明は、必ずしも、図面等に開示された位置、大きさ、形状、範囲などに限定されない。

　図２に本実施例で処理する測定データを取得する、Ｘ線撮像装置を示す。本装置は主にＸ線源１，被写体位置決め回転機構２，Ｘ線画像検出器３，制御部４，演算部５，及び表示部６から構成されている。被写体位置決め回転機構２は、被写体８を所定のＣＴ回転軸１０を中心に回転させることができる。Ｘ線源１から出射した照射Ｘ線７は，被写体位置決め回転機構２で位置決めされた被写体８に照射される。被写体８を透過した透過Ｘ線９は，Ｘ線画像検出器３で検出される。制御部４，演算部５，及び表示部６は、纏めて制御装置を構成する。制御装置１１は、この他、記憶装置や、オペレーターが操作するためのキーボードなどの入力装置を備えてもよい。

　図３は、図２の装置を用い、制御部４の制御によって，ＣＴの測定データを取得する処理手順を示すフロー図である。

　処理Ｓ３０１では、被写体位置決め回転機構２を用いて被写体８を照射Ｘ線の光路７から退避する。

　処理Ｓ３０２では、Ｘ線を照射せずにダークデータI_Dを取得する。

　処理Ｓ３０３ではＸ線を照射する。

　処理Ｓ３０４で背景画像データI_BKを取得する。

　処理Ｓ３０５では、被写体位置決め回転機構２を用いて，被写体８を照射Ｘ線の光路７に挿入する。この際，被写体８の回転中心と光路の中心が一致するように調整すれば，Ｘ線の光路の幅（視野）と同程度の大きさの被写体を測定することが可能になる。

　処理Ｓ３０６では、被写体８を透過した透過像データI(θ)を取得する。

　処理Ｓ３０７では、被写体８を被写体位置決め回転機構２を用いてΔθだけ回転する。

　処理Ｓ３０８では、Ｓ３０６からＳ３０７のステップを、被写体の回転角度θが０から３６０°になるまで繰り返して行う。

　処理Ｓ３０９では、被写体位置決め回転機構２を用いて被写体８を照射Ｘ線７の光路から退避する。

　処理Ｓ３１０では、Ｘ線を照射して背景画像データI_BK_2を取得し、Ｓ３１１で終了する。

　なお，ダークデータI_D及び各背景画像データI_BK及びI_BK_2は複数枚測定し，その平均値を採用することで，より揺らぎの少ない画像を得ることができる。

　以上によって取得したＣＴのデータセットを用いて，はじめに以下の前処理を行い，各投影角度θおける線吸収係数の像データを表すI_2を求める。

　図４は、制御装置１１による前処理の処理手順を示すフローである。

　処理Ｓ４０１では、測定前後の背景画像データI_BK及びI_BK_2を用いて，各投影データI(θ)における背景画像I’_BKを
　　　I’_BK = (1－a)×I_BK ＋ a×I_BK_2
により求める。ここで，aは各背景画像の重み関数で，a=θ/360 である。

　すなわち，θが０°近傍では，直近に取得したI_BKを主な背景画像データとし，θが360°近傍では，I_BK_2を主な背景画像データとする。なお，像データは１或いは２次元配列のデータであり，上記の計算は各画素における個々の計算を意味している。

　処理Ｓ４０２では、上記によって算出した背景画像データI’_BKと，ダークデータI_Dを用いて，規格化された投影データI_1(θ)を
　　　I_1(θ)=(I(θ)－I_D)/(I’_BK－I_D)
により求める。

　処理Ｓ４０３では、Ｉ_1（θ）について，
　　　I_2(θ)=－ln(I_1(θ))
を計算し，線吸収係数データI_2(θ)を求める。

　なお，被写体のサイズや密度が大きく，Ｘ線が被写体を全く透過しない場合，I_1(θ)は0となるため，上記（３）の計算はln(0)となり，エラーになってしまう。この際は，Ｘ線のエネルギーを向上したり，露光時間の延長を行った上で，再度上記のＣＴ測定を行う。

　図５に、処理装置１１による断面像の再構成の手順を示す。処理は，前処理で求めた線吸収係数データI_2(θ)を用いて，以下の手順により行う。

　処理５０１では、投影角度θに対して直交する角度（θ＋９０°）の線吸収係数データI_2(θ＋９０)について，各画素位置xにおける値に基づき、マスクデータM(x”)を以下により計算する。

　　　M(x”) = 1 if I_2(θ＋90, x”)＜＞0
　　　　　　　= 0 if I_2(θ＋90, x”)=0
　すなわち，被写体を観察した投影像データを元に，被写体の存在しない（I_2=0）領域は０，存在する(I_2＜＞0)は１とするマスクデータを作成する。ただし、「＜＞」は「等しくない」ことを示す。なお、計算の前に所定閾値以下のI_2を、0に近似して置き換える処理を行ってもよい。

　処理５０２では、仮想的な再構成像データS(x’, y’)の各座標(x’, y’)について
　　　x = x’ sin θ ＋ y’cos θ
　　　x” = x’ sin (θ＋90°) ＋ y’cos （θ＋90°）
の座標変換を行う。

　処理５０３では、この座標に対応したI_2とMを用いて，
　　　S(x, y) = S(x, y) ＋ M(x”)×f(I_2(θ))
を計算する。なお，ここでf( )はShepp－Logan関数やランプ関数などのフィルター関数のコンボリューションである。

　処理５０４では、上記の手順について，投影角度θが０°から３６０°になるまで繰り返して行う。なお，１８０°以降のマスクデータM(x”)の計算は，θ－90°のデータを用いる。３６０°になると終了（Ｓ５０５）である。

　以上によって算出した再構成像データSは，オペレーターの指示に従って，表示器６で表示する。ここで、制御部４と演算部５の処理動作は、例えば計算機（コンピュータ）によるプログラムの実行を通じて実現される。プログラムは、不図示の記憶装置に格納されている。また、制御部４や演算部５で実行される処理動作の一部又は全部は、ＡＳＩＣや専用の計算モジュールとして実現されてもよい。なお、図２では、表示器６が演算部５に接続されているが、出力装置としての印刷装置が、表示器６に代えて又は表示器６とともに、演算部５に接続されていてもよい。印刷装置は、前述の再構成像を印刷する。また、表示器６に代えて又は表示器６とともに通信装置を接続し、外部装置に送信できるようにしてもよい。

　図６に図１と同様のサイズ５×５画素で、被写体が中心にしか存在しない単純な系を対象とした適用例を示す。この系で得られる各投影像データは図１（ａ）で示したようになる。

　θ＝０°のデータ１０２を逆投影で加算するとき，θ＝９０°の投影データ１０１から、中心（座標（０，０））以外には被写体がないことが分かるので、中心以外の領域には値を入れない。すなわち、マスクデータMは (0, 0, 1, 0, 0)となる。したがって，上記（２）の加算は、領域６０２に示すように中心（０，０）だけに行うことになる。

　この処理は、図１（ｂ）と図６を比較すれば明らかなように、θ＝９０°の投影データ１０１に基づいて、θ＝０°のデータ１０２を加算すべき領域を特定していることになる。データ処理上は、θ＝９０°の投影データ１０１でマスクデータを生成し、θ＝０°のデータ１０２の加算時にフィルタリング処理を行っていることと等価である。

　次にθ＝９０度のデータ１０１を加算するとき，θ＝０の投影データ１０２から中心以外に被写体がないことがわかるので，マスクデータMは同じように(0, 0, 1, 0, 0)となる。このため，領域６０１に示すように中心（０，０）だけに同様に値を加算する。以上によって，オリジナルのデータと同様に，被写体のある領域だけ値を持ち，他の領域は０となる。すなわち，アーチファクト（疑似像）のない，断面像を再生することができる。

　以上，本実施例によれば，再構成の逆投影計算において，投影するデータと直交した投影データの値が０でない領域を１，０の領域を０とするマスクデータを作成し，このマスクデータと投影データの積を再生構成データに加算することによって，アーチファクトのない断面像を再生することができる。

　実施例１では，投影データにフィルター関数をコンボリューションした後に加算を行っていた。このため，再生された断面像における各画素のＣＴ値はフィルター関数の種類に依存して変化するため一定の値を得ることができず，定量的な解析が難しかった。本実施例では，フィルター関数のコンボリューションを行う必要がなく，断面像を再生する方法について説明する。

　本実施例では実施例１のマスクデータM(x”)を計算する過程について，以下のように変更する。

　　　M’(x”) = I_2(θ＋90°, x”)/ΣI_2(θ＋90°,x_i”)
　ここで，分母はθ＋９０°の投影像データの総和で，分子は投影像データの位置x”における値である。すなわち，マスクデータM’は規格化した投影像データとなる。そして，逆投影はこのM’を用いて，
　　　S(x, y) = S(x, y) ＋ M’(x”)×I_2(θ)
　により，計算する。実施例１との違いは，フィルター関数のコンボリューションが無いことである。

　以上により求めた断面像Sは，実施例１と同様にオペレーターの指示に従って，表示器で表示する。また，その実現方法も同様に，計算機（コンピュータ）によるプログラムや，ＡＳＩＣや専用の計算モジュールで実行してもよい。

　図７に、５×５の小さい画素の測定データを対象として，本発明を適用した例を示す。なお，図１と異なり，被写体は中心部の１画素だけでなく，１×２の大きさを有している。この系で、測定時において、投影角度θ＝０°及び９０°で得られる投影データ１０２，１０１は，図１と同様に縦及び横方向の積算であるために，図７（ａ）のような値となる。

　この投影データを用いた再構成計算の逆投影は，上記に従いはじめに＋９０°の投影データ１０１を用いてマスクデータM’を計算する。本例では，＋９０°の投影データI_2(90°)は
　　　I_2(90°) = (0, 1, 2, 0, 0)
となっているので，総和は3となる。このため，マスクデータM’は
　　　M’ = (0, 1/3, 2/3, 0, 0)
となる。次にこのマスクデータを用いて，θ=0°の投影データ１０２を縦方向に加算することになるが， θ＝０°における投影像データI_2（０）は
　　　I_2(0) = (0, 0, 3, 0, 0)
のため，実際の加算は中央の縦方向だけに行う。加算するデータはマスクデータM’との積であることから
　　　S(0, y) = S(0, y) ＋ M’(y)×I_2(0, 0)
　　　　　　　= S(0, y) ＋ (0, 1, 2, 0, )
となる。すなわち

となる。なお，θ = 0°では座標の変換
　　　x” = x’ sin (θ＋90°) ＋ y’cos （θ＋90°）
から，M’の座標はxとyが入れ替わった形になっている。次にθ＝９０°の計算では，上記と同様にはじめにθ＝0°の投影データ１０２を用いてマスクデータM’を作成する。０°のデータI_2(0)は
　　　I_2(0) = (0, 0, 3, 0, 0)
であることから，
　　　M’ = (0, 0, 1, 0, 0)
となる。θ=90°の逆投影は横方向への投影データの加算になるので，
　　　S(x, y) = S(x, y) ＋ M’(x)×I_2(90, y)
となり，M’(x)×I_2(90, y)の演算結果である以下のデータをS(x, y)に加算することになる。

したがって，最終的なS(x, y)は

となり，投影データの加算回数2で除算すると，

が得られる。この結果は被写体の分布と一致しており，本法により被写体の像を正確に再現可能なことがわかる。

　以上，本実施例によれば，再構成の逆投影計算において，投影するデータI_2(θ)と直交した投影データI_2(θ＋９０°)の各位置における値をその総和で除算したマスクデータを作成し，このマスクデータと投影データの積を再生構成データに加算することによって，フィルターを用いることなく，アーチファクトのない断面像を再生することができる。

　実施例１および２では，１回の計算だけで断面像を再構成していた。このために，実施例１ではフィルターをコンボリューションし，実施例２では逆投影するデータと直交するデータを用いてその加算の比率を決めていた。本実施例では繰り返しの計算（逐次計算）によってこのような計算を行うことなく，断面像を再生する方法について説明する。

　逐次計算における１回の計算過程は図５に示した実施例１と同様に，
　処理Ｓ５０１で、投影角度θに対して直交する角度（θ＋９０°）の線吸収係数データI_2(θ＋９０)について，各画素xにおける値を元にマスクデータM(x”)を以下により計算する。

　　　M(x”) = 1 if I_2(θ＋90, x”) ＜＞ 0
　　　　　　　= 0 if I_2(θ＋90, x”) = 0
すなわち，横から画像データを元に，被写体の存在しない（I_2=0）領域は０，存在する(I_2＜＞0)領域は１とする。

　処理Ｓ５０２で、仮想的な再構成像データS(x’, y’)の各座標について
　　　x = x’ sin θ ＋ y’cos θ
　　　x” = x’ sin (θ＋90°) ＋ y’cos （θ＋90°）
の座標変換を行う。

　処理Ｓ５０３で、この座標に対応したI_2とMを用いて，
　　　S(x, y) = S(x, y) ＋ M(x”)×f(I_2(θ))
を計算する。

　処理Ｓ５０４で、上記の手順について，投影角度θが０°から３６０°になるまで繰り返して行う。なお，１８０°以降のマスクデータM(x”)の計算は，θ－90°のデータを用いる。

　但し，処理Ｓ５０３において加算するのはMとI_2の単純な積であり，フィルターのコンボリューションは行わないことが実施例１との大きな違いである。

　逐次計算は上記の計算によって得られた再構成像データSを初期の再構成像データS_0として以下の手順により行う。

　最初に、上述のように、再構成像データS_0を再構成像データの初期値に設定する。
（Ａ）S_0を用いて，投影像データI_k（θr）を以下により求める。なお，ここでθrは投影角度で，計算はθrが０から３６０°まで，測定データと同じ角度間隔で行う。

　　　I_k(θr, xr) = ΣS(xr_i, yr_i)
　　　xr_i = i cos θr ＋ xr sin θr
　　　yr_i = i sin θr ＋ xr cos θr
（Ｂ）得られたI_k（θr）について，図５に示した逆投影計算を行い，再構成像データS_kを求める。
（Ｃ）S_0とS_kの差dS_k
　　　dS_k = 2(S_0 - S_k)
を用いて，次ステップにおける再構成像データS_k＋1を
　　　S_k＋1 = S_k ＋ dS_k
から求め，さらに次の投影像データI_k＋1を上記と同様の計算で算出する。

　以上の一連の計算過程（Ａ）から（Ｂ）を
　　　E = Σ|(I_2－I_k)|
で与えられる誤差Eが予め設定した値以下となるまで行う。これにより，フィルターのコンボリューション等を用いることなく，定量性の高い断面像を再構成することができる。

　図８には，以上の計算過程のフローチャートを示す。なお，上記処理Ｓ５０２の計算過程において，ＣＴ値の上限や非負など既知の情報を制約条件として組み込むことで，さらに正確なＣＴ値の断面像を再構成することができる。

　処理Ｓ８０１では、測定データから初期再構成像データS_0を計算し、ｋを１に設定する。

　処理Ｓ８０２では、S_k－1の投影から投影像データI_kを計算する。

　処理Ｓ８０３では、E = Σ|(I_2－I_k)|で与えられる誤差Eを計算する。

　処理Ｓ８０４では、誤差Eが設定値か否かを判定する。設定値を満たせば終了し、そのときの再構成データ（再生画像）が結果となる（終了処理８０５）。誤差が設定値を満たさない場合は、処理Ｓ８０６に進む。

　処理Ｓ８０６では，I_kの逆投影によりS’_k－1を計算し再構成像データを求める。

　処理Ｓ８０７では，dS = 2(S_0 － S’_k－1)を計算して、初期再構成像データと再構成データの差分dSを求める。

　処理Ｓ８０８では、S_k = S_k－1 ＋ dSを計算して、次ステップにおける再構成像データを求め、k = k＋1に設定して、処理Ｓ８０２に戻る。

　以上により求めた断面像Sは，実施例１と同様にオペレーターの指示に従って，表示器６で表示する。また，その実現方法も同様に，計算機（コンピュータ）によるプログラムや，ＡＳＩＣや専用の計算モジュールで実行してもよい。さらに，繰り返しの計算過程を表示器６で毎回表示することによって，収束の過程や制約条件の妥当性を計算の途中でも判断することできるようになる。

　以上，本実施例によれば，繰り返し計算で断面像を再構成する方法の逆投影計算において，投影するデータと直交した投影データの値が０でない領域を１，０の領域を０とするマスクデータを作成し，このマスクデータと投影データの積を再生構成データに加算することによって，アーチファクトのない定量性の高い断面像を再生することができる。

　実施例１から３では，逆投影において，投影するデータI_2（θ）と直交する投影データO_2(θ＋９０°)を用いて，マスクデータを作成していた。本実施例では，９０°以外の投影データを用いてマスクデータを作成し，これを用いて断面像を再構成する方法について説明する。

　図９に投影角度θ＝０°における投影データI_2(0)と、任意の投影角度αにおける投影像データI_2(θ＋α)の関係を示した模式図を示す。ここで，矩形の領域９０１が角度αにおいて，被写体９０２が存在していることになる領域，すなわちマスクデータが０でない領域である。従って、この領域では、投影角度αのデータからの加算を許すことになる。

　このことから，α＝９０°のとき（実施例１から３に該当），マスクデータM”はI_2の横方向の位置x’に対して不変となるが，９０°以外の角度ではx’に依存してマスクデータの位置x”をシフトする必要があることがわかる。

　本実施例では，実施例１のマスクデータMを作成する過程（処理Ｓ５０１）を、投影角度θに対して角度α（θ＋α）の線吸収係数データI_2(θ＋α)について，各画素の位置xにおける値を元にマスクデータM”(x”)を以下により計算する。

　　　M”(x”) = 1 if I_2(θ＋α, x_m) ＜＞ 0
　　　　　　　　= 0 if I_2(θ＋α, x_m) = 1
　ただし、x_m = x” sin αとする。なお，αが小さくなるとx”に対応するI_2(x_m)がデータ範囲外となってしまうので，あまり小さくすることはできない。そして，上記によって作成したマスクデータM”を用いて，本実施例では座標変換処理Ｓ５０２と再構成像データ作成処理Ｓ５０３において、以下のように再構成の加算を行う。

　仮想的な再構成像データS(x’, y’)の各座標について以下の変換
　　　x = x’ sin θ ＋ y’cos θ
　　　x” = x’ sin (θ＋α) ＋ y’cos （θ＋α）
を行い，このx及びx”に対応する，I_2(θ,x)にフィルター関数をコンボリューションしたf(I_2)と，マスクデータM(x”)の積をSに以下のように加算する。

　　　S(x, y) = S(x, y) ＋ M(x”)×I_2(θ)
　αが小さくなるに従って，マスクデータと投影データが重なる領域が大きくなり，その効果が減少して被写体周辺でアーチファクトが発生しやすくなるが，被写体より離れた領域では依然として有効に機能する。すなわち，直交からずれた投影データを使用してもアーチファクトの発生を抑えることができる。また，上記のマスクデータM”を作成する過程において，実施例２のように，I_2(θ＋α)のデータを用いて，M”の値を0 or１でなく，加算の比率としても良い。この場合，フィルター関数のコンボリューションが不要となるために，より定量的なＣＴ値の断面像を求めることができる。さらに，上記処理Ｓ５０２～Ｓ５０３の計算過程においてコンボリューションを除き，実施例３のように繰り返し計算により，像を再生してもよい。この場合も，より定量的なＣＴ値の断面像を求めることができる。

　以上，本実施例によれば，再構成の逆投影計算において，投影するデータと投影角度がαだけ異なる投影データの値に基づいてマスクデータを作成し，このマスクデータと投影データの積を再生構成データに加算することによって，アーチファクトのない断面像を再生することができる。

　図１０は、本発明の実施例における照射Ｘ線１００１と被写体１００２の回転角度の関係を示す原理図である。

　実施例１から４では，図１０（ａ）に示すように被写体１００２の回転軸１００３に対して，直交する断面像を求めていた。本実施例では，図１０（ｂ）のように被写体の回転軸１００３が入射Ｘ線１００１と直交せず，βだけ傾いた系における断面像の計算について説明する。なお，一般に回転軸とＸ線光路が傾いた系をラミノグラフィーと呼び、その傾き角βをラミノグラフィーアングル１００４と呼ぶ。このように回転軸を傾けることによって，板状の縦と横の長さと厚さが大きく異なる被写体も投影角度によらずに積分した厚さは一定となるので，通常のＣＴと同様に測定することが可能になる。

　ラミノグラフィーでは，投影角度によって被写体の垂直方向（Ｚ軸方向）の位置もビーム（画像検出器）に対して変動することになるため，再構成を３次元で考慮する必要がある。

　図１１には厚みのある円盤状の被写体を真横から測定した形状（a）および斜めから測定した形状(b)の例を示す。なお，説明を簡略化するために，内部の吸収係数は非常に大きく，Ｘ線は透過しないものとして扱っている。斜め方向から観察した像では，両画像共に上下に丸みを帯びた形状をしていることがわかる。このため，マスクデータとしてＺ軸方向に対する制限はどうしても甘くなるが，水平面内に対しては従来通り有効に機能する。マスクデータM_3の具体的な作成方法は実施例１を２次元に拡張する。

　具体的には、実施例１のマスクデータの作成処理Ｓ５０１において、投影角度θに対して直交する角度（θ＋９０°）の線吸収係数データI_2(θ＋９０, x”, y”)について，各画素(x”, y”)における値を元にマスクデータM_3(x”, y”)を以下により計算する。

　　　M_3(x”, y”) = 1 if I_2(θ＋90, x”, y”) ＜＞ 0
　　　　　　　　　　= 0 if I_2(θ＋90, x”, y”) = 0
　すなわち，横から被写体を観察した画像データを元に，被写体の存在しない（I_2=0）領域は０，存在する(I_2＜＞0)は１とする。

　そして，座標変換処理Ｓ５０２において、このマスクデータM_3を用いて投影角度θにおける逆投影は仮想的な再構成像データS(x’, y’, z’)の各座標について以下の変換を行う。

　　　x = x’ sin θ ＋ y’cos θ
　　　y = x’ cos θsin β－ y’sin θsin β＋z’ cos β
　　　x” = x’ sin (θ＋α) ＋ y’cos (θ＋α)
　　　y” = x’ cos (θ＋α)sin β－ y’sin (θ＋α)sin β＋z’ cos β
　そして，これらの座標に対応する，I_2(θ,x, y)にフィルター関数をコンボリューションしたf(I_2)と，マスクデータM_3(x”, y”)の積をSに以下のように加算する。
S(x’, y’, z’) = S(x’, y’, z’) ＋ M_3(x”, y”)×I_2(θ, x, y)
　以上により，水平面内方向にはアーチファクトが無い断面像を再生することができる。なお，フィルター関数のコンボリューションは水平面内（X－Y面）に行えば像を再構成することができるが，より正確なＣＴ値を有した像を再生するためには，Z方向も含めた２次元的なフィルター関数（Shepp－Loganの２次元関数など）を２次元的にコンボリューションする必要がある。また，実施例３のように、上記の過程を含めた繰り返しの計算によって像を再生することができる。この計算において，被写体が円盤状であることなどの既知の情報を制約条件として組み込むことができればより正確に像を再生することができる。なお，本実施例ではマスクデータを直交したI_2(θ＋９０°)を利用しているが，当然のことながら実施例４と同様に９０°と異なる角度のデータを利用して作成してもよい。

　以上により求めた断面像Sは，実施例１から４と同様にオペレーターの指示に従って，表示器６で表示する。また，その実現方法も同様に，計算機（コンピュータ）によるプログラムや，ＡＳＩＣや専用の計算モジュールで実行してもよい。また，特に繰り返し計算で求める場合は，３次元的な計算を行うために，膨大な計算量となる。このため，超高速な演算が可能な計算機や専用のハードウェアを実装した演算部５を用いると良い。

　以上，本実施例によれば，入射Ｘ線と被写体の回転軸が直角でない測定系においても，再構成の逆投影計算において，投影するデータと投影角度がαだけ異なる投影データの値に基づいてマスクデータを作成し，このマスクデータと投影データの積を再生構成データに加算することによって，アーチファクトのない断面像を再生することができる。

　本明細書において単数形で表される構成要素は、特段文脈で明らかに示されない限り、複数形を含むものとする。

　以上の構成においてデータ処理は、単体のコンピュータで構成してもよいし、あるいは、入力装置、出力装置、処理装置、記憶装置の任意の部分が、ネットワークで接続された他のコンピュータで構成されてもよい。発明の思想としては等価であり、変わるところがない。

　本実施例中のデータ処理は、ソフトウエアで構成することができる。また、同等の機能は、FPGA（Field Programmable Gate Array）、ASIC（Application Specific Integrated Circuit）などのハードウエアでも実現できる。そのような態様も本願発明の範囲に含まれる。

　本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。例えば、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることが可能である。また、各実施例の構成の一部について、他の実施例の構成の追加・削除・置換をすることが可能である。

１　Ｘ線源
２　被写体位置決め回転機構
３　Ｘ線画像検出器
４　制御部
５　演算部
６　表示器
７　照射Ｘ線
８　被写体
９　透過Ｘ線

Claims

　Ｘ線源とＸ線画像検出器と演算部を備え、
　前記Ｘ線源から被写体に対して複数の異なる方向からＸ線を照射し、前記Ｘ線画像検出器で前記複数の異なる方向から照射されたＸ線によって得られる複数の異なる投影角度で得られた投影像を取得し、前記投影像の画像データを演算部によって処理し、前記被写体の断面像を逆投影によって再構成する、コンピュータ断層撮影装置において、
　前記演算部は、
　第１の投影角度で得られた前記画像データからマスクデータを作成する第１の処理部、
　第２の投影角度で得られた前記画像データをデータ配列に加算して前記断面像を再構成する際に、前記画像データを前記マスクデータで重み付けして前記断面像の一部を形成する処理を行う第２の処理部、
　を有し、
　前記第１の投影角度および第２の投影角度を変化させ、前記断面像の一部を形成する処理を繰り返すことにより、前記断面像を再構成するコンピュータ断層撮影装置。
　前記第１の投影角度と前記第２の投影角度が直交する角度である、請求項１記載のコンピュータ断層撮影装置。
　前記第１の処理部は、
　前記第１の投影角度で得られた画像データによって、前記被写体が存在しないと判定される座標が０の値となる前記マスクデータを作成し、
　前記第２の処理部は、
　前記第２の投影角度で得られた前記画像データと前記マスクデータを積算する処理を行う、請求項１記載のコンピュータ断層撮影装置。
　前記第２の処理部は、
　前記第２の投影角度で得られた画像データもしくは前記マスクデータの少なくとも一つを座標変換して演算を行う、請求項１記載のコンピュータ断層撮影装置。
　前記第２の処理部は、
　前記第２の投影角度で得られた前記画像データをフィルター関数により演算処理する、請求項１記載のコンピュータ断層撮影装置。
　前記第１の処理部は、
　前記第１の投影角度で得られた前記画像データを規格化して前記マスクデータを作成する、請求項１記載のコンピュータ断層撮影装置。
　前記逆投影によって再構成された第１の断面像から第１の投影像データを計算し、前記第１の投影像データから逆投影によって第２の断面像を計算し、前記第１の断面像と第２の断面像の差分を計算し、前記差分により前記第１の断面像を修正する、第３の処理部を有する、請求項１記載のコンピュータ断層撮影装置。
　前記第１の投影角度と前記第２の投影角度が直交する角度でない場合、前記マスクデータの座標をシフトする処理を行う、請求項１記載のコンピュータ断層撮影装置。
　前記Ｘ線源から被写体に対して複数の異なる方向からＸ線を照射するために、
　前記被写体の回転軸とＸ線光路が傾いた系を用い、前記マスクデータとして２次元の値を用いる、請求項１記載のコンピュータ断層撮影装置。
　被写体とＸ線撮像系を相対的に回転して、異なる投影角度で得られた画像データから、被写体の断面像を逆投影加算によって再構成するコンピュータ断層撮影方法において、
　第１の投影方向から取得した画像データを前記逆投影加算する際に、これと異なる第２の投影方向から取得した画像データを用いて、再構成する断面像の加算する領域を制限することを特徴とするコンピュータ断層撮影方法。
　請求項１０記載のコンピュータ断層撮影方法において、上記制限は前記第２の投影方向から取得した画像データ上で前記被写体が存在する領域だけ加算し，存在しない領域は加算しない制限であることを特徴とするコンピュータ断層撮影方法。
　請求項１０記載のコンピュータ断層撮影方法において、上記制限は前記第２の投影方向から取得した画像データの各位置の信号強度に比例した係数を、前記第１の投影方向から取得した画像データに乗算することによって行われることを特徴とするコンピュータ断層撮影方法。
　請求項１０記載のコンピュータ断層撮影方法において、前記第１の投影方向と第２の投影方向は、直交することを特徴とする。コンピュータ断層撮影方法。
　請求項１０記載のコンピュータ断層撮影方法において、被写体を照射するＸ線と被写体の回転軸が直交しないことを特徴とするコンピュータ断層撮影方法。
　被写体とＸ線撮像系を相対的に回転して異なる投影角度で得られた画像データから繰り返し計算によって被写体の断面像を再構成するコンピュータ断層撮影において、
　前記被写体の断面像を再構成するための逆等投影計算が、異なる投影方向から取得したデータ上で被写体が存在する領域だけ加算し，存在しない領域は加算しない制限を設けた逆投影計算であることを特徴とするコンピュータ断層撮影方法。