(57)【要約】２次元検出器を使用し３ＰＩアルゴリズムを用いたスパイラル走査される移動する対象物の画像を復元する方法及びシステムである。移動する対象物は、従来のシステムよりも最大約３倍ゆっくりした割合で走査することができる。好ましい実施形態では、本発明は、従来のスパイラル走査システムのものと同じ大きさの検出器アレイを使用したままで、スパイラル走査装置の中で寝台に載せた患者を最大約３倍ゆっくりと移動させることができる。
2006524059.tif 000002
【特許請求の範囲】
【請求項１】
少なくとも１つの検出器から得られたデータから画像を復元する方法であって、
軸方向に、現在の放射源位置近傍におけるスパイラル４回転分の投影よりも幅広のコーンビーム投影を少なくとも１つ検出する検出器を少なくとも１つ使用して、対象物をスパイラル方式で走査するステップと、
畳み込みに基づくＦＢＰ（フィルタ補正逆投影）アルゴリズムにより、走査した対象物の精密な画像を効率良く復元するステップとを含むことを特徴とする方法。
【請求項２】
請求項１に記載の方法であって、
前記走査ステップは、前記検出器を使用して前記対象物の２次元コーンビーム（ＣＢ）投影データを取得することを特徴とする方法。
【請求項３】
請求項２に記載の方法であって、
１Ｐ１アルゴリズムに必要とされる検出器と実質的に同様の検出器を使用するステップをさらに含むことを特徴とする方法。
【請求項４】
請求項２に記載の方法であって、
前記走査ステップは、１Ｐ１アルゴリズムで使用されるコーンビーム投影よりも軸方向に幅広のコーンビーム投影を検出するステップをさらに含むことを特徴とする方法。
【請求項５】
請求項４に記載の方法であって、
前記走査ステップは、１Ｐ１アルゴリズムで使用されるコーンビーム投影よりも軸方向に少なくとも３倍幅広のコーンビーム投影を検出するステップをさらに含むことを特徴とする方法。
【請求項６】
請求項１に記載の方法であって、前記対象物は人を含むことを特徴とする方法。
【請求項７】
２次元検出器を使用したスパイラル・コンピュータ断層撮影走査によって得られた精密な画像をコンピュータ処理する方法であって、
（ａ）１Ｐ１アルゴリズムに必要とされる検出器よりも幅広い検出器から、現在の放射源位置近傍におけるスパイラル４回転分の投影の範囲をカバーするコーンビーム（ＣＢ）投影データを収集するステップと、
（ｂ）コーンビーム投影と交差する平面П上の線群を特定するステップと、
（ｃ）ＣＢ投影データを前処理するステップと、
（ｄ）前記線に沿って、前記前処理したＣＢ投影データを畳み込みフィルタ処理するステップと、
（ｅ）前記画像の前段階と称すべきものを作成すべく、前記フィルタ処理したデータを逆投影するステップと、
（ｆ）前記対象物の精密な画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ），（ｄ）及び（ｅ）を繰り返すステップとを含むことを特徴とする方法。
【請求項８】
請求項７に記載の方法であって、前記走査は前記対象物のＸ線暴露を含むことを特徴とする方法。
【請求項９】
請求項７に記載の方法であって、前記ステップ（ａ）〜（ｆ）は３ＰＩアルゴリズムを含むことを特徴とする方法。
【請求項１０】
検出器を使用したコンピュータ断層撮影走査によって得られた画像をコンピュータ処理する方法であって、
（ａ）対象物の走査中に検出器からコーンビーム（ＣＢ）データを収集するステップと、
（ｂ）前記コーンビームと交差する平面ＤＰ（ｓ）上における、次の３つの線群を特定するステップと、
（ただし、ｓは、走査経路を表し、現在の放射源位置に相当するパラメータ値である）
（ｂｉ）第１の線群は、現在の放射源位置における走査接線の方向である
2006524059.tif 000031
に対して平行である。
（ｂｉｉ）第２の線群は、ｓ＜ｑ＜ｓ＋２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ１、及びｓ−２π＜ｑ＜ｓによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ−１に対して接する。
（ただし、ｑは、走査経路に沿ったパラメータであり、投影される点を表す）
（ｂｉｉｉ）第３の線群は、ｓ＋２π＜ｑ＜ｓ＋４πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ２、及びｓ−４π＜ｑ＜ｓ−２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ−２に対して接する。
（ｃ）前記３つの群の前記線に沿って、前記データを前処理及び推移不変フィルタ処理するステップと、
（ｄ）前記画像の前段階と称すべきもの作成すべく、前記フィルタ処理したデータを逆投影するステップと、
（ｅ）前記対象物の画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ）及び（ｄ）を繰り返すステップとを含むことを特徴とする方法。
【請求項１１】
請求項１０に記載の方法であって、
前記前処理は、現在の放射源位置についてのＣＢデータの微分計算を含むことを特徴とする方法。
【請求項１２】
請求項１０に記載の方法であって、
前記推移不変フィルタ処理は、前記前処理されたデータを１／ｓｉｎγフィルタにより畳み込むことを特徴とする方法。
【請求項１３】
請求項１０に記載の方法であって、
前記第１の線群からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、
2006524059.tif 000032
に平行であり、Γ２に接するＬ２ｃｒと、
2006524059.tif 000033
に平行であり、Γ−２に接するＬ２−ｃｒとの間に位置する場合に、前記フィルタ処理されたデータに係数ｃｍ＝２／３を乗じることを特徴とする方法。
【請求項１４】
請求項１０に記載の方法であって、
前記第１の線群の線からの前記フィルタ処理されたデータの逆投影は、
ｘのＤＰ（ｓ）への投影が、Ｌ２ｃｒの上側に位置する、又はＬ２−ｃｒの下側に位置する場合に、前記フィルタ処理されたデータに係数ｃｍ＝１／３を乗じることを特徴とする方法。
【請求項１５】
請求項１０に記載の方法であって、
前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
ｘのＤＰ（ｓ）への投影が、Γ１とΓ−１との間に位置し、線がΓ１∪Γ−１と接する点がｘの１ＰＩパラメータ区間の内側にある場合に、前記フィルタ処理されたデータに係数ｃｍ＝２／３を乗じることを特徴とする方法。
【請求項１６】
請求項１０に記載の方法であって、
前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
ＤＰ（ｓ）へのｘの投影が、Γ１とΓ−１との間に位置し、線がΓ１∪Γ−１と接する点がｘの１ＰＩパラメータ区間の外側にある場合に、前記フィルタ処理されたデータに係数ｃｍ＝−２／３を乗じることを特徴とする方法。
【請求項１７】
請求項１０に記載の方法であって、
前記第３の線群の線からの前記フィルタ処理されたデータの逆投影は、
前記フィルタ処理されたデータに係数ｃｍ＝１／３を乗じることを含むことを特徴とする方法。
【請求項１８】
検出器を使用したコンピュータ断層撮影走査によって得られた画像をコンピュータ処理する方法であって、
（ａ）対象物の走査中に検出器からコーンビーム・データを収集するステップと、
（ｂ）前記コーンビームと交差する平面ＤＰ（ｓ）上における、次の３つの線群を特定するステップと、
（ただし、ｓは、走査経路を表し、現在の放射源位置に相当するパラメータ値である）
（ｂｉ）第１の線群は、現在の放射源位置における走査接線の方向である
2006524059.tif 000034
に対して平行である。
（ｂｉｉ）第２の線群は、ｓ＜ｑ＜ｓ＋２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ１、及びｓ−２π＜ｑ＜ｓによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ−１に対して接する。
（ただし、ｑは、走査経路に沿ったパラメータであり、投影される点を表す）
（ｂｉｉｉ）第３の線群は、ｓ＋２π＜ｑ＜ｓ＋４πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ２、及びｓ−４π＜ｑ＜ｓ−２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ−２に対して、少なくとも３つの交点ｓ１，ｓ２，ｓ３を有する。
（ｃ）前記３つの群の前記線に沿って、前記データを前処理及び推移不変フィルタ処理するステップと、
（ｄ）前記画像の前段階と称すべきもの作成すべく、前記フィルタ処理したデータを逆投影するステップと、
（ｅ）前記対象物の画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ）及び（ｄ）を繰り返すステップとを含むことを特徴とする方法。
【請求項１９】
請求項１８に記載の方法であって、
前記交点ｓ１，ｓ２，ｓ３は、次の規則に従って決定されることを特徴とする方法。
ｓ＋２π＜ｓ３＜ｓ＋４πならば、ｓ１−ｓ＝Ψ（ｓ３−ｓ２）
ｓ−４π＜ｓ３＜ｓ−２πならば、ｓ３−ｓ２＝Ψ（ｓ１−ｓ）
（ただし、Ψ（ｔ）は、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという性質を有する関数である）
【請求項２０】
請求項１８に記載の方法であって、
前記前処理は、現在の放射源位置に関するＣＢデータの微分計算を含むことを特徴とする方法。
【請求項２１】
請求項１８に記載の方法であって、
前記推移不変フィルタ処理は、前記前処理されたデータを１／ｓｉｎγフィルタにより畳み込むことを含むことを特徴とする方法。
【請求項２２】
請求項１８に記載の方法であって、
前記第１の線群からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、Ｌ２ｃｒとＬ２−ｃｒとの間に位置する場合に、前記フィルタ処理されたデータに係数ｃｍ＝２／３を乗じることを特徴とする方法。
【請求項２３】
請求項１８に記載の方法であって、
前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
ＤＰ（ｓ）へのｘの投影が、Γ１とΓ−１との間に位置し、線がΓ１∪Γ−１と接する点がｘの１ＰＩパラメータ区間の内側にある場合に、前記フィルタ処理されたデータに係数ｃｍ＝２／３を乗じることを特徴とする方法。
【請求項２４】
請求項１８に記載の方法であって、
前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
ｘのＤＰ（ｓ）への投影が、Γ１とΓ−１との間に位置し、線がΓ１∪Γ−１と接する点がｘの１ＰＩパラメータ区間の外側にある場合に、前記フィルタ処理されたデータに係数ｃｍ＝２／３を乗じることを特徴とする方法。
【請求項２５】
請求項１８に記載の方法であって、
前記第３の線群の線からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、Ｌ２ｃｒ及びＬ−２ｃｒの上側に位置する場合に、前記フィルタ処理されたデータに係数ｃｍ＝−２／３を乗じることを特徴とする方法。
【請求項２６】
コンピュータ断層撮影システムにおいて、対象物の走査に基づいて画像復元するのに使用される線群を特定する方法であって、
（ｉ）現在の放射源位置ｓを決定するステップと、
（ただし、ｓは、走査経路を表すパラメータである）
（ｉｉ）コーンビーム投影と交差する平面ＤＰ（ｓ）を選ぶステップと、
（ｉｉｉ）パラメータ値を用いて、次のｓ１，ｓ２，ｓ３として表される走査経路上の３点を選択するステップと、
（ｉｉｉａ）｜ｓ−ｓ１｜＜２π
（ｉｉｉｂ）２π＜ｓ２−ｓ＜４π及び２π＜ｓ３−ｓ＜４π、又は、−４π＜ｓ２−ｓ＜−２π及び−４π＜ｓ３−ｓ＜−２π
（ｉｉｉｃ）２π＜ｓ３−ｓ＜４πならばｓ１−ｓ＝Ψ（ｓ３−ｓ２）、又は−４π＜ｓ３−ｓ＜−２πならばｓ３−ｓ２＝Ψ（ｓ１−ｓ）、
（ただし、Ψ（ｔ）は、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという関数である）
（ｉｖ）前記３点をＤＰ（ｓ）上に投影するステップと、
（ｖ）前記投影を通る境界線を引くステップとを含むことを特徴とする方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、米国特許仮出願第６０／４３０，８０２号（出願日：２００２年１２月４日）の優先権の利益を主張するものである。また、本発明は、米国特許出願第１０／３８９，５３４号（出願日：２００３年３月１４日）の一部継続出願である。なお、米国特許出願第１０／３８９，５３４号は、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）の一部継続出願である、米国特許出願第１０／３８９，０９０号（出願日：２００３年３月１４日）の一部継続出願である。
【０００２】
本発明は、コンピュータ断層撮影（computer tomography：CT）に関し、特に、３ＰＩアルゴリズムを使用したスパイラル走査によって取得したデータから３次元画像を復元するプロセス、方法及びシステムに関する。
【背景技術】
【０００３】
この３０年間、コンピュータ断層撮影（ＣＴ）は、スライス・バイ・スライス走査に基づく画像復元から、スパイラル走査による画像復元まで行われてきた。１９７０年代から１９８０年代までは、スライス・バイ・スライス走査が用いられた。この方法では、寝台に載せた患者をガントリーの中で少しずつ移動させるのと、ガントリーの回転とが交互に行われる。ガントリーの回転中、患者は静止しているので、患者の周りを回るＸ線源の軌道は円形になる。事前に選択された患者を横切るスライスは、そのような円形走査によって取得されたデータを使用して復元される。１９８０年代の中頃から今日まで、ＣＴでのデータ収集には、スパイラル型の走査が好ましい方法となった。スパイラル走査では、患者を載せた寝台は、寝台の周りを連続的に回転するガントリーの中を連続的に移動する。初めは、スパイラル走査は、データを１次元で受け取る１次元検出器を使用していた（検出器は一列に配置される）。その後、検出器が複数列（２列以上）で互いに隣接して配置される２次元検出器が導入された。ＣＴでは、特に２次元検出器を使用する場合に、画像の復元について重大な問題があった。以下では、２次元検出器で作成されるデータを、コーンビーム（cone-beam：CB）データ又はＣＢ投影と呼ぶ。
【０００４】
２次元検出器を使用したスパイラス走査により得られたデータからの３次元（別名：ボリュームメトリック）画像復元用としては、厳密アルゴリズムと近似アルゴリズムの２種類のアルゴリズムが知られている。厳密アルゴリズムと近似アルゴリズムは、どちらも既知の問題を有している。理想的な環境下では、厳密アルゴリズムは、精密な画像を復元することができる。したがって、厳密アルゴリズムには、非理想的な（すなわち、現実的な）環境下でも良質な画像を作成することが期待される。また一方、厳密アルゴリズムは、画像を復元するのに長時間かかることが知られている。また、厳密アルゴリズムは、その使用時にコンピュータの処理能力を多大に必要とする。これらのアルゴリズムは、大量のコーンビーム投影をメモリに保存する必要がある。さらに、ある厳密アルゴリズムでは、実施するために大きな検出アレイを必要としたり、走査する患者の大きさに制限があったりする。
【０００５】
近似アルゴリズムは、フィルタ補正逆投影（filtered back projection：FBP）の構造を有しているので、画像をとても効率良く作成することができる。また、厳密アルゴリズムよりも、使用するコンピュータの処理能力が少ない。また一方、近似アルゴリズムは、理想的な環境下でも、近似画像（精密な画像とは似ているものの、精密な画像とは異なっている）を作成することができる。特に、近似アルゴリズムは、画像における偽の特徴であるアーチファクトが作成され得る。ある環境下では、このようなアーチファクトはかなり強まる。
【０００６】
理論的には厳密であり推移不変（shift invariant）のＦＢＰ構造を有する前者のアルゴリズムは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）に開示されている。これらのアルゴリズムの短所は、１ＰＩモードで動作することである。このことは、検出器アレイが軸方向に大きい場合は、検出器すべてを使用するために、ガントリー内で患者を急いで平行移動させる必要があることを意味する。しかしながら、患者を急いで移動させることは、明白な理由により実用的ではない。とは言え、患者をゆっくりと移動させると、検出器の一部しか使用できない。この場合、Ｘ線照射の一部が無駄になる、打切りアーチファクト（discretization artifact）が強まる、ノイズ安定性が減少するなどの望ましくない結果をもたらす。
【発明の開示】
【課題を解決するための手段】
【０００７】
本発明の主な目的は、２次元検出器を使用したスパイラル方式により走査された対象物の画像を復元するための３ＰＩアルゴリズムを提供することにある。任意のボクセルで画像を復元するためには、これらのアルゴリズムは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許出願第６，５７４，２９９号）の１ＰＩのアルゴリズムよりも、スパイラルの断面が長いことが必要とされる。したがって、新規のアルゴリズムは、同じ大きさの検出器アレイを使用したままで、患者を約３倍ゆっくりと移動させることを可能にする。
【０００８】
本発明の第１の好ましい実施形態では、スパイラル走査された対象物の画像を復元するのに、５つの全体的なステップから成るプロセスを使用する。第１のステップでは、現在のＣＢ投影が測定される。次に、新規のアルゴリズムに従って、検出器上の線群が特定される。次に、隣接する投影間で導関数のコンピュータ処理が行われ、その後、選択された線群からの線に沿って、フィルタ処理によって導関数の畳み込みが行われる。次に、フィルタ処理されたデータを使用して逆投影を行うことにより、画像は更新される。最終的に、この前処理ステップは、全ての対象物が走査されるまで、各ＣＢ投影について繰り返される。この実施形態は、いくつかの（約２〜４）ＣＢ投影を同時にメモリに保存し、１つの線群を使用する。
【０００９】
第２の実施形態では、いくつかのＣＢ投影をメモリに保存するのと共に、異なる線群を利用することができる。
【００１０】
これらの実施形態を変更して、一度に、１つのＣＢ投影だけをコンピュータ・メモリに保存することが可能である。これは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許出願第６，５７４，２９９号）（この参照により本発明に含まれるものとする）で行われたのと同じようにして実現することができる。
【００１１】
したがって、新規のアルゴリズムは、同じ大きさの検出器アレイを使用したままで、患者を約３倍ゆっくりと移動させることを可能にする。
【００１２】
本発明のさらなる目的及び利点は、以下の添付図面に概略的に図示した、現時点での好ましい実施形態についての詳細な説明によって明らかになるであろう。
【発明を実施するための最良の形態】
【００１３】
本発明の開示された実施形態を詳細に説明する前に、本発明はここで説明される詳細によって限定されるものではなく、他の実施形態をとることもできることを理解されたい。また、ここで使用される専門用語は、説明目的のためのものであって、本発明を限定するためのものではない。
【００１４】
本発明は、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）（この参照により本発明に含まれるものとする）の一部継続出願である。
【００１５】
以下、本発明について詳細に説明する。
【００１６】
《理論背景》
まず、必要な表記について説明する。
2006524059.tif 000003
と仮定する。ただし、ｈ＞０はスパイラル（らせん）であり、Ｕはスパイラルの完全に内側の開集合である。
2006524059.tif 000004
Ｓ２は、Ｒ３における単位球面である。そして、
2006524059.tif 000005
Ｄｆ（ｙ，β）は、ｆのＣＢ変換である。
2006524059.tif 000006
とすると、
2006524059.tif 000007
ｊ＝１，２，…，は、平面
2006524059.tif 000008
とＣとの交点を表すものとする。また、
2006524059.tif 000009
である。
β∈Ｓ２に関しては、β⊥は大円｛α∈Ｓ２：α・β＝０｝を表す。コンピュータ処理が必要なｆについて、任意のｘ∈Ｒ３を決定する。ｆ（ｘ）をコンピュータ処理するために、後に決定される、有限の範囲のスパイラルの断面を使用する。とりあえず、今はＣ（ｘ）で示す。対応するパラメータ区間（parametric interval）はＩ（ｘ）で示される。Ｃ（ｘ）についての主な仮定は、下記の通りである。
【００１７】
〈性質Ｃ１〉（完全状態）
ｘを通る任意の平面は、少なくとも１点でＣと交差する。
【００１８】
反転公式の構成における重要な要素は加重関数ｎ（ｓ，ｘ，α），α∈β⊥（ｓ，ｘ）である。関数ｎは、次のように理解される。ｘとαは平面П（ｘ，α）を決定し、加重関数ｎは、ｘの位置によって決まるｙ（ｓ）∈П（ｘ，α）ПＣに割り当てられる。従って、この解釈から、ｎ（ｓ，ｘ，α）＝ｎ（ｓ，ｘ，−α）であると推定する。ｎについての主な仮定は、下記の通りである。
【００１９】
〈性質Ｗ１〉（通常状態）
2006524059.tif 000010
〈性質Ｗ２〉
ｎ（ｓ，ｘ，α）が任意の（ｓ，α）の近傍で局地的に一定である、有限の多くのＣ１関数αｋ（ｓ，ｘ）∈β⊥（ｓ，ｘ），ｓ∈Ｉ（ｘ）が存在する。ただし、ｓ∈Ｉ（ｘ）であり、α∈β⊥（ｓ，ｘ），
2006524059.tif 000011
である。
2006524059.tif 000012
を示す。
【００２０】
Ｃｌ、Ｗ１、及びＷ２の仮定の下では、次の一般的な反転公式は、A. Katsevichの論文「A general scheme for constructing inversion algorithms for cone beam CT," International Journal of Mathematics and Mathematical Sciences, Vol. 21, pp. 1305-1321 (2003)」から導かれる。
2006524059.tif 000013
ただし、θｍ∈（０，π）は、φ（ｓ，ｘ，θ）が不連続になる点であり、ｃｍ（ｓ，ｘ）はジャンプ値である。
2006524059.tif 000014
ｎが適切に選択される限り、反転公式は必ずしもＦＢＰ型を必要としない。
【００２１】
《３ＰＩ線及びその性質》
あるｓ０∈Ｒについて、Ｘ線源がｙ（ｓ０）に固定されていると仮定する。検出器アレイはＸ線源と共に回転するので、検出器平面（detector plane）はｓ０よって決まり、ＤＰ（ｓ０）で示される。ＤＰ（ｓ０）はスパイラルの軸と平行であり、Ｘ線源と反対側の点において円柱：ｙ１２＋ｙ２２＝Ｒ２（（１）参照）に接すると仮定する。したがって、ｙ（ｓ０）と検出器平面との間の距離は、２Ｒである（図３参照）。検出器平面内の座標は、次の通りである。ｄ１軸がスパイラル軸と垂直であり、ｄ２軸がスパイラル軸と平行であり、起点はｙ（ｓ０）の投影と一致すると仮定する。スパイラルの上側及び下側の回転を検出器表面上に立体投影したのを図３に示す。これにより、次のパラメータ曲線が与えられる。
2006524059.tif 000015
2006524059.tif 000016
ただし、ρは、患者の支持体（support）の半径によって決まる。
ρ＝２ｃｏｓ−１（ｒ／Ｒ）（（２）参照）。これらの曲線は、Γｊ，ｊ＝±１，±２,…,で示される（図４参照）。
2006524059.tif 000017
はｘの投影を示す。
（なお、
2006524059.tif 000018
は、明細書中では、「ｘ´」とも表記する）
任意のＸ線源位置を、スパイラル上の全ての点ｙ（ｓ）に結び付ける。ただし、
2006524059.tif 000019
である。
【００２２】
これにより、ＳＵ３ＰＩ（ｓ０）とＳL３ＰＩ（ｓ０）で表される２つの表面が得られる。この２つの表面と、円柱：ｘ１２＋ｘ２２＝Ｒ２とによって囲まれた領域は、Ｖ３ＰＩ（ｓ０）で示される。ｘは、固定されていると仮定する。ｓ０が十分に小さければ、ＳＵ３ＰＩ（ｓ０）はｘよりも下側になる。ｓ０が増大すると、２つの場合が有り得る。１つは、連続的照度（continuous illumination）として知られており、ｘはＳＵ３ＰＩ（ｓ０）を通ってＶ３ＰＩ（ｓ０）に入り、ＳL３ＰＩ（ｓ０）を通ってＶ３ＰＩ（ｓ０）から出る。明らかに、上記の方法は、一意的な３ＰＩ区間［ｂ０（ｘ），ｔ０（ｘ）］と、対応する一意的な線Ｌ０３ＰＩ（ｘ）を与える。２番目の場合は、遮断照度（interrupted illumination）として知られており、ｘは、Ｖ３ＰＩ（ｓ０）に何回も出入りする。より正確には、ｘは、ＳＵ３ＰＩ（ｓ０）とＳＬ３ＰＩ（ｓ０）の各表面と、正確に３回づつ交差する。したがって、上記の方法は、３つの３ＰＩ線を与える。Ｌｉ３ＰＩ（ｘ），ｉ＝１，２，３。対応するパラメータの値は、ｂｉ（ｘ），ｔｉ（ｘ），ｉ＝１，２，３と表される。スパイラルは対称であるため、ｓ０=ｂｉ（ｘ）の場合、ｘはＶ３ＰＩ（ｓ０）から出る（に入る）。そして、ｓ０=ｔｉ（ｘ），ｉ＝１，２，３の場合、ｘはＶ３ＰＩ（ｓ０）に入る（から出る）。平面Ｘ３＝０と見なす。この２つの場合の間の境界線を図５に示す。図５から、曲線は自分自身とは交差しておらず、
円：ｘ１２＋ｘ２２＜１を、Ｘ１とＸ３の領域に分割していることが分かる。Ｘ１を示す中央の領域では、各ｘに対して１つの３ＰＩ線がある。Ｘ１を示す外側の領域では、各ｘに対して３つの３ＰＩ線がある。
【００２３】
次の特性を確立できる。点ｘに対して３ＰＩが１つしかない場合は、
2006524059.tif 000020
である。
【００２４】
点ｘに対して３つの３ＰＩがある場合は、
2006524059.tif 000021
である。
【００２５】
検出器平面に関しては、方程式（１３ａ）及び（１３ｂ）は、ｓ∈Ｉ３ＰＩ（ｘ）の時かつその時に限り、ｘ´はΓ２とΓ−２との間にあることを示す。したがって、１ＰＩの場合との類推から、この検出器上の領域を３ＰＩウインドウと呼ぶことができる。１ＰＩの場合では、ｘの１ＰＩ線の終点によって囲まれるパラメータ区間は、ｘの１ＰＩパラメータ区間と呼ばれる。同様に、Ｉ３ＰＩ（ｘ）は、ｘの３ＰＩパラメータ区間と呼ばれる。
【００２６】
《補助構成》
まず、連続的照度が行われると仮定する。単位球面Ｓ２の表面上にある次の２つの曲線を検討する。一方の曲線は、Ｌ０Ｉ３ＰＩ（ｘ）に対して垂直である全ての単位ベクトルから成り、Ａで示される。他方の曲線は、ベクトル
2006524059.tif 000022
から成り、Ｔで示される。これらの曲線は、
2006524059.tif 000023
により定義される球座標（θ１，θ２）を使用して、平面上で示されるので、これらの曲線をＳ２上に直接的に示すのは、都合が良くない。αと−αは同一の平面を定義するので、θ１を任意の長さπの区間に限定することができる。識別することにより、反対側の境界を「接着（glue）」することができる。
2006524059.tif 000024
ｘ／Ｒ＝（０，０．２５，０）の典型的な状態を図６に示す。単位ベクトルとしてだけではなく、平面についても、Ｓ２上の点について考えることは都合がよい。各α∈Ｓ２は、法線ベクトルαを有するｘを通る一意的な平面と対応している。この対応は、反対方向のベクトルが特定できる場合は、１対１である。ＡとＴとで、Ｓ２を複数の領域Ｄ１，Ｄ２,…,に分割する。固定されたｉの構成によって、Ｃ３ＰＩ（ｘ）∩П（ｘ，α）における点の数は、全てのα∈Ｄｉについて同じである。Ｄｉがｋ交点（ＩＰｓ）を含む場合は、ｋＩＰ領域を呼ばれる。
【００２７】
次に、遮断照度が行われる場合を考える。同様の方法で、Ｓ２の表面上に、いくつかの曲線を定義する。初めの３つの曲線は、それぞれ３つの３ＰＩ軸に垂直な単位ベクトルを考慮しながら取得される。それらは、Ａｋ、ｋ＝１，２，３と示される。曲線の第２の集合は、式（１５）中のｓを、区間［ｂ１，ｂ２］，［ｂ３，ｔ１］，［ｔ２，ｔ３］に限定することにより得られる。これらの曲線は、それぞれＴ１２、Ｔ３１、及びＴ２３と示される。ｘ／Ｒ=（−０．５，０，０）についての典型的な状況を図７に示す。
【００２８】
ｘ１及びｘ２が互いに近接していれば、それらは同じダイアグラムを共有することは明らかである。したがって、数列を滑らかなに変化させることにより、あるダイアグラムは他のダイアグラムに変換される。そして対応する領域では、ＩＰの数と、Ｉ３ＰＩでの部分区間におけるそれらの分配は、同じままである。本質的な変化は、３つの境界がＳ２上の一点で互いに交差するという、「重大な事象」が起こるｘの近傍でのみ可能性がある。これらの点は、数的に発見することができる。その結果を図８に示す。これらの点から回転の中心までの最少距離は、
2006524059.tif 000025
である。したがって、医療用に関心のある状況下では（ｒFOV５Ｒ）、図６に示した連続的照度と、図７に示した遮断照度との、２つの場合が起こり得る。
【００２９】
《加重関数ｎの構成》
任意のｓ∈Ｉ３ＰＩ（ｘ）について、次の規則に従って、加重関数ｎ（ｓ，ｘ，α）、α∈β⊥（ｓ，ｘ）を決定する。
【００３０】
・１ＩＰ領域では、唯一のＩＰは重量ｎ＝１を得る。
【００３１】
・３ＩＰ領域では、各ＩＰは重量ｎ＝１／３を得る。
【００３２】
・５ＩＰ領域では、２つのＩＰは重量ｎ＝０を得、残り全てのＩｐはｎ＝１／３を得る（図９及び図１０参照）。
【００３３】
《第１の復元アルゴリズム》
重量が求められると、アルゴリズムを構成する全ての要素が得られる。まず第１に、ｘの画像を復元するために、ｘ：Ｃ（ｘ）＝Ｃ３ＰＩ（ｘ）とＩ（ｘ）＝Ｉ３ＰＩ（ｘ）の３ＩＰスパイラル分割を行う（Section 1参照）。（６）や（７）から、各ｓ∈Ｉ３ＰＩ（ｘ）について、フィルタ処理する方向は、
2006524059.tif 000026
の切れ目を発見することにより決定されることが分かる。これらの切れ目の研究は、次の３つの線群をもたらす。
【００３４】
第１の群は、スパイラル接線に平行な線から成る。この群は、Ｌ０と表される。この線と関連する係数を、図１１に示す。明らかに、任意のｘ´について、ｘ´を含む１つのＬ∈Ｌ０だけがある。
【００３５】
図１２は、フィルタ処理する線の第２の群を示す。この群は、Γ±１に接する線から成り、Ｌ１と表される。１ＰＩウインドウ内の各ｘ´に対して、Ｌ１から２つの線を引く。これらの線は、図１３に示した規則に従って決定する。また、この図は、関連した係数ｃｍを示す。
【００３６】
図１４（左側の図）は、フィルタ処理したΓ±１に接する線の群を示す。この群は、Ｌ２と表される。任意のｘ´∈Ｆ１∪Ｆ２∪Ｆ４∪Ｆ５では、１つ以上のｘ´を含む線Ｌ∈Ｌ２がある。一意的な線は、次の規則に基づいて決定される（図１４の右側の図を参照されたい）。ｘ´がｘ´∈Ｆ１∪Ｆ４である場合は、接点はｘ´の右側となる。また、ｘ´がｘ´∈Ｆ２∪Ｆ４である場合は、接点はｘ´の右側となる。全ての場合において、ｃｍ＝１／３である。
【００３７】
図１５は、図１１〜１４に含まれる情報を要約しており、ｘ´となり得る全ての可能性がある場合と、全ての関連したフィルタ処理方向及び定数ｃｍを示している。全ての場合において、フィルタ処理の方向は左から右である仮定する。この構成の結果として（ｃｆ．図１５）、一般的な復元式（７）、（８）中のθｍ（ｓ，ｘ）及びｃｍ（ｓ，ｘ）は、β（ｓ，ｘ）だけによって、ｘによってのみ決まる。したがって、ｃｍとα⊥の理由で、ｘをβ（ｓ，ｘ）によって置き換えることができ、式（７）を次の形式で書き直すことができる。
2006524059.tif 000027
また、ここで留意すべきは、フィルタ処理する線Ｌが定められると、その投影がＬに含まれるかつ上記した規則を満たす全てのｘは、同じフィルタ処理線Ｌを共有する（図１５参照）。したがって、（１９）は畳み込まれ、（１８）は逆投影となり、アルゴリズム（１８）・（１９）は畳み込みに基づくＦＢＰ型となる。このことは、米国特許第１０／１４３，１６０号（出願日：２００２年５月１０日）にも同様に見られる。逆投影法に関しては、他の方法及び技術を用いることもできる。さらに、例えば、Ｋａｔｓｅｖｉｃｈによる米国特許第６，５７４，２９９号を参照されたい（この参照により本発明に含まれるものとする）。
【００３８】
《改善された復元アルゴリズム》
Ψ（ｔ）をＲ上に定義される、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという性質を有する滑らかな関数と仮定する。Γ±１∪Γ±２：ｓ１，ｓ２，ｓ３と３つの交差点を有し、それらの点が
2006524059.tif 000028
を満たす任意の既知の線Ｌ∈Ｌ２´を求めることにより、新しい線群Ｌ２´を定義する。ｓがｘ線現の現在位置であることを思い出されたい。この線Ｌ∈Ｌ´２は、例えば、
ｓ３，２π＜｜ｓ３｜＜４πによりパラメータ化される。ｓ１及びｓ２の位置は、ｓ３の位置によって決まる（図１６に示す）。
【００３９】
Ψの性質を利用して、それぞれのｘ´∈Ｆ１∪Ｆ２∪Ｆ４∪Ｆ５について、ｘ´を含む一意的なＬ∈Ｌ´２を容易に決定することができる。また、
ｘ´→Ｌ２ｃｒであれば、ｓ２，ｓ３→２Δであり、ｓ１→ｓである。同様に、
ｘ´→Ｌ−２ｃｒであれば、ｓ２，ｓ３→−２Δであり、ｓ１→ｓである。
【００４０】
さらに重要なことには、ｘ´∈Ｆ１∪Ｆ２∪Ｆ４∪Ｆ５の場合は、フィルタ処理方向を２つから１つに減らすことができる。この付加的な利益は、検出器の使用法を改善する。したがって、図１５に２つの上側の図面は、図１７に示したダイアグラムと置き換えられる。ｘ´がＬ０（即ち、ｘ´∈Ｆ４∪Ｆ５）の下側に現れる場合は、検出器平面の基点に対する対称により、図１７から得られる。
【００４１】
反転公式（１８）、（１９）の形式は、変更されない。第１の違いは、ｘ´∈Ｆ１∪Ｆ２∪Ｆ４∪Ｆ５の場合は、Ｍ（ｓ，β）＝１となることである（２ではなく）。第２の違いは、フィルタ処理の方向であり、ここでは式（２１）、（２２）で表される。
【００４２】
《復元アルゴリズムの一般的な概要》
図２は、本発明の基本処理ステップ１０、２０、３０、４０、５０の概要を示す。以下、各ステップについて説明する。
【００４３】
ステップ１０：現在のＣＢ（コーンビーム）投影を、コンピュータ・メモリにロードする。現在メモリに記憶されているＣＢ投影の中央点を、ｙ（ｓ０）と仮定する。
【００４４】
ステップ２０：フィルタ処理する線群を求める。図１５を使用して、又は改善されたアルゴリズムの場合はさらに図１７を使用して、必要な線群を特定する。その後、各群から線の離散集合を選択する。
【００４５】
ステップ３０：フィルタ処理の準備。ステップ２０で選択された前記線上の点を、極角によりパラメータ化する。補間法を用いて、前記極角の離散値に対応する前記線上の点βにおけるＣＢデータ
2006524059.tif 000029
の導関数を計算する。
【００４６】
ステップ４０：フィルタ処理。ステップ２０で特定された各線について、フィルタ１／sinγによってステップ３０で計算されたデータを畳み込む。
【００４７】
ステップ５０：逆投影。各復元点ｘについて、ステップ４０で求められたフィルタ処理されたデータを、方程式（１８）に従って逆投影する。その後、処理すべき新しいＣＢが無くなるまで、又は全ての必要な点ｘにおいて画像復元が完成するまで、ステップ１０に進む。
【００４８】
次に、図２に示した５つのステップ１０〜５０のアルゴリズムについて、詳細に説明する。
【００４９】
〈ステップ１０〉
現在のＣＢ（コーンビーム）投影を、コンピュータ・メモリにロードする。現在メモリに記憶されているＣＢ投影の中央点を、ｙ（ｓ０）と仮定する。
【００５０】
〈ステップ２０〉
フィルタ処理する線群を求める。線群は、次のサブステップ２１、２２、及び２３により選択することができる。
【００５１】
ステップ２１：線群Ｌ０から、スパイラル接線と平行であり、Γ２とΓ−２との間に位置する検出器の表面上における関心領域の投影を覆う線の等距離集合を選択する（図１１を参照）。
【００５２】
ステップ２２：線群Ｌ１から、Γ１とΓ−１に接する線の離散集合を選択する（図１２参照）。Γ１の接線における左端の点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。同様に、Γ−１の接線における右端の点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。
【００５３】
ステップ２３：線群Ｌ２から、Γ２とΓ−２に接する線の離散集合を選択する（図１４の左側の図を参照されたい）。どちらの場合も、接点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。改善されたアルゴリズムの場合は、Γ２とΓ−２に接する線の代わりに、Γ２とΓ−２曲線上のｓ３に対する離散（等距離）集合の値を選択し、方程式（２１）、（２２）を解くことにより、線Ｌ∈Ｌ２´を決定する。両方の曲線上で、点ｓ３は、検出器の表面上における関心領域の投影を超えて広がる必要はない。
【００５４】
〈ステップ３０〉
フィルタ処理の準備。
【００５５】
ステップ３１：ステップ２０で得られた線の集合から、線Ｌを決定する。
【００５６】
ステップ３２：ｙ（ｓ０）とＬを通る平面内で、前記線上の点を極角γによりパラメータ化する。
【００５７】
ステップ３３：後にステップ４０で離散フィルタ処理に使用される、等距離値γｊの離散集合を選択する。
【００５８】
ステップ３４：各γｊについて、ｙ（ｓ０）から、γｋに対応するＬ上の点に向う単位ベクトルβｊを求める。
【００５９】
ステップ３５：ｓ０に近い少しの値ｑについてのＣＢ投影データＤｆ（ｙ（ｑ），Θ）を使用して、全てのΘ=βｊについての導関数
2006524059.tif 000030
を数的に求める。
【００６０】
ステップ３６：計算された導関数の値をコンピュータ・メモリに記憶させる。
【００６１】
ステップ３７：ステップ２０で特定されたすべての線Ｌについて、ステップ３１〜３６を繰り返す。このようにして、ｙ（ｓ０）に位置するＸ線源に対応する処理されたＣＢ投影データを作成する。
【００６２】
〈ステップ４０〉
フィルタ処理。
【００６３】
ステップ４１：ステップ２０で特定された前記線群Ｌｍの内の１つから、線Ｌを決定する。
【００６４】
ステップ４２：前記線に沿って、ステップ３０で計算された前記処理されたＣＢデータの値のＦＦＴを計算する。
【００６５】
ステップ４３：フィルタ１／sinγのＦＦＴを計算する。
【００６６】
ステップ４４：フィルタ１／sinγのＦＦＴ（ステップ４３の結果）と、前記処理されたＣＢデータの値のＦＦＴ（ステップ４２の結果）とを乗算する。
【００６７】
ステップ４５：ステップ４４の結果結果の逆ＦＦＴを行う。
【００６８】
ステップ４６：ステップ４５の結果をコンピュータ・メモリに記憶させる。
【００６９】
ステップ４７：前記線群の全ての線について、ステップ４１〜４６を繰り返す。このことにより、フィルタ処理されたＣＢデータΨｍ（ｓ０，βｊ）が得られる。ただし、ｍは、Ｌが選択される線群の番号を表す。ｍ＝０，１，２。
【００７０】
フィルタ処理は、当該技術分野では周知であり、例えば、タム（Tam）による米国特許第５，８８１，１２３号（この参照により本発明に含まれるものとする）に示され説明されているようにして実施される。
【００７１】
〈ステップ５０〉
逆投影。
【００７２】
ステップ５１：画像復元が求められている患者の内側の点である復元点ｘを求める。
【００７３】
ステップ５２：ｓ０がＩ３ＰＩ（ｘ）に属する場合は、前記フィルタ処理されたＣＢデータは、ｘにおける画像に影響を与える。そして、ステップ５３〜５８を実行する。ｓ０がＩ３ＰＩ（ｘ）区間の内側にない場合は、前記フィルタ処理されたＣＢデータはｘにおける画像復元に使用されない。この場合は、ステップ５１に戻り、他の復元点を選択する。
【００７４】
ステップ５３：検出器の平面ＤＰ（ｓ０）上にｘの投影ｘ´を求める。そして、ｙ（ｓ０）からｘに向う、単位ベクトルβ（ｓ０，ｘ）を求める。
【００７５】
ステップ５４：図１１、図１３、及び図１４の右側の図を使用して、前記線群から線を特定する。また、前記投影ｘ´に近い前記線上の点を特定する。ｘ´がＬ２ｃｒよりも上側、又はＬ２−ｃｒよりも下側であれば、改善されたアルゴリズムを使用する場合は、前記投影ｘ´に近いフィルタ処理された線を見つけるために、方程式（２１）、（２２）、及び図１６を使用する。このことにより、β（ｓ０，ｘ）に近いβｊについての少しの値Ψｍ（ｓ０，βｊ）が得られる。
【００７６】
ステップ５５：補間法を使用して、前記したβ（ｓ０，ｘ）に近いβｊについての値Ψｍ（ｓ０，βｊ）の値から、Ψｍ（ｓ０，β（ｓ０，ｘ））の値を推定する。
【００７７】
ステップ５６：Ψｍ（ｓ０，β（ｓ０，ｘ））に−ｃｍ（ｓ，β（ｓ０，ｘ）／［４π２｜ｘ−ｙ（ｓ０）｜］を掛けることにより、前記フィルタ処理されたＣＢデータから、点ｘで復元される画像に対する寄与率（contribution）を計算する。適切な逆投影率ｃｍは、図１１、図１３、及び図１４の右側の図を使用して選択される（概要については図１５を参照されたい）。ｘ´がＬ２ｃｒよりも上側、又はＬ２−ｃｒよりも下側であれば、改善されたアルゴリズムの場合は、適切な逆投影率ｃｍを見つけるために図１７を使用する。
【００７８】
ステップ５７：方程式（１８）での積分の近似値を求めるために、前記寄与率を、点ｘにおいて復元される前記画像に加える。なお、前記画像は、予め選択された方法（例えば、台形法）に従って復元される
ステップ５８：ステップ５１に進み、異なる復元点ｘを選択する。全ての復元点ｘが処理されれば、ステップ５９に進む。
【００７９】
ステップ５９：ステップ１０へ進み、次のＣＢ投影を、コンピュータ・メモリに読み込む。
【００８０】
画像データは、画像復元処理が完了した全ての復元点ｘで表示することができる（すなわち、それらの点では、画像復元するためのそれに続くＣＢ投影は全て必要としない）。コンピュータ・メモリからは、画像復元処理が完了していない点での画像復元に不要な全てのＣＢ投影を捨て去る。アルゴリズムは、走査が終了したとき又は必要な点での画像復元が全て完了したときに、結論を出す。
【００８１】
例えば、検出器が３ＰＩアルゴリズムに必要なデータを取得しなかった場合は、欠けているデータを推定するために、様々な方法を採ることができる。この場合は、検出器は、各コーンビーム投影の測定データ及び推測データを取得する仮想の検出器となる。欠けているデータを正確に推測することができれば、３ＰＩアルゴリズムは正確に復元することができる。ここでは、依然として、欠けているデータが近似的に求められる現実の環境下における正確な復元について述べている。
【００８２】
以上、本発明の実施の形態について説明したが、本発明は前記した実施の形態のみに限定されるものではなく、本発明の技術的思想に基づく限りにおいて、種々の変形が可能である。
【図面の簡単な説明】
【００８３】
【図１】Ｘ線源及び検出器アレイを有する回転するガントリー内を移動する寝台に載せた患者の一般的な配置を示す。コーンビーム投影データ集合は、Ｘ線検出器によって受け取られる。そして、コンピュータにより画像復元処理が行われる。復元された画像はモニタに表示される。
【図２】本発明の基本処理ステップの概要を示す。
【図３】検出器平面上におけるスパイラルの立体投影を示す。
【図４】検出器平面上の様々な投影及び重要な線を示す。
【図５】境界曲線を示す。
【図６】連続的照度の場合を示す。ｘ／Ｒ＝（０，０．２５，０）
【図７】遮断照度の場合を示す。ｘ／Ｒ＝（−０．５，０，０）
【図８】様々な重要な事柄が生じる点である。
【図９】連続的照度の場合における５ＰＩ領域の内側の重量分布を示す。
【図１０】遮断照度の場合における５ＰＩ領域の内側の重量分布を示す。
【図１１】Ｌ０に平行なフィルタ処理された線の群を示す。
【図１２】Γ±１に接するフィルタ処理された線の群を示す。
【図１３】フィルタ処理された線を、ｘ´の位置によって選択する方法を示す。
【図１４】Γ±２に接するフィルタ処理された線の群を示す。
【図１５】様々な場合における、フィルタ処理された線及び関連する定数ｃｍを示す。左上の図：ｘ´∈Ｆ１、右上の図：ｘ´∈Ｆ２、中央の図：ｘ´∈Ｆ３、及びΓ１以上、左下の図：ｘ´∈Ｇ１、右下の図：ｘ´∈Ｇ２。
【図１６】点ｓ１，ｓ２，ｓ３が取り得る位置を示す。
【図１７】様々な場合における、フィルタ処理された線及び関連する定数ｃｍ。左の図：ｘ´∈Ｆ１、右の図：ｘ´∈Ｆ２、
【図１８】図２に示したステップ２０に対応する、フィルタ処理する線群を求める３つのサブステップを表わすフローチャートである。
【図１９】図２に示したステップ３０に対応する、フィルタ処理準備の７つのサブステップを表わすフローチャートである。
【図２０】図２に示したステップ４０に対応する、フィルタ処理の７つのサブステップを表わすフローチャートである。
【図２１】図２に示したステップ５０に対応する、逆投影の５つのサブステップを表わすフローチャートである。
【図２２】図２に示したステップ５０に対応する、逆投影の後半部分の３つのサブステップを表わすフローチャートである。
【図１】
2006524059.tif 000035
【図２】
2006524059.tif 000036
【図３】
2006524059.tif 000037
【図４】
2006524059.tif 000038
【図５】
2006524059.tif 000039
【図６】
2006524059.tif 000040
【図７】
2006524059.tif 000041
【図８】
2006524059.tif 000042
【図９】
2006524059.tif 000043
【図１０】
2006524059.tif 000044
【図１１】
2006524059.tif 000045
【図１２】
2006524059.tif 000046
【図１３】
2006524059.tif 000047
【図１４】
2006524059.tif 000048
【図１５】
2006524059.tif 000049
【図１６】
2006524059.tif 000050
【図１７】
2006524059.tif 000051
【図１８】
2006524059.tif 000052
【図１９】
2006524059.tif 000053
【図２０】
2006524059.tif 000054
【図２１】
2006524059.tif 000055
【図２２】
2006524059.tif 000056
【国際調査報告】
2006524059.tif 000057
2006524059.tif 000058