(57)【特許請求の範囲】
【請求項１】
対象物を囲む円周上に設けたＸ線源および複数の検出器からの透過データによりＣＴ画像を再構成するＣＴ画像再構成システムであって、
前記透過データから、角度θのｒ軸上における不等間隔の投影データｐ（ｒm，θn）を得る手段と、
前記不等間隔の投影データｐ（ｒm，θn）から、直接的に極座標（ρ，θ）における等間隔の空間周波数スペクトルＧ（ρk，θn）を求める１次元フーリエ変換を行う１次元フーリエ変換手段と、
前記空間周波数スペクトルＧ（ρk，θn）から、直接２次元フーリエ逆変換を行い、直交座標系（ｘ，ｙ）のＣＴ画像であるｆ（ｘl，ｙm）を求める２次元フーリエ逆変換手段と
を備え、
前記１次元フーリエ変換手段は、Δmを投影データの検出点が占める微少範囲とすると、
【数１】
0003751788.tif 000023
により空間周波数スペクトルを求め、
前記２次元フーリエ逆変換手段は、空間周波数スペクトル標本値が極座標（ρ，θ）上でになう微少面積をΔＳkとすると、
【数２】
0003751788.tif 000024
によりＣＴ再構成画像を求めることを特徴とするＣＴ画像再構成システム。
【請求項２】
請求項１に記載のＣＴ画像再構成システムにおいて、
透過データから投影データｐ（ｒm，θn）を得る前記手段は、前記複数の検出器からの透過データをメモリに記憶してから、角度θに対応したｒ軸上の投影データを選択することにより、投影データｐ（ｒm，θn）を得ることを特徴とするＣＴ画像再構成システム。
【請求項３】
対象物を囲む円周上に設けたＸ線源および複数の検出器からの透過データによりＣＴ画像を再構成するために、
前記透過データから、角度θのｒ軸上における不等間隔の投影データｐ（ｒm，θn）を得る手段と、
前記不等間隔の投影データｐ（ｒm，θn）から、直接的に極座標（ρ，θ）における等間隔の空間周波数スペクトルＧ（ρk，θn）を求める１次元フーリエ変換を行う１次元フーリエ変換手段と、
前記空間周波数スペクトルＧ（ρk，θn）から、直接２次元フーリエ逆変換を行い、直交座標系（ｘ，ｙ）のＣＴ画像であるｆ（ｘl，ｙm）を求める２次元フーリエ逆変換手段として、コンピュータに機能させるプログラムであって、
前記１次元フーリエ変換手段は、Δmを投影データの検出点が占める微少範囲とすると、
【数３】
0003751788.tif 000025
により空間周波数スペクトルを求め、
前記２次元フーリエ逆変換手段は、空間周波数スペクトル標本値が極座標（ρ，θ）上でになう微少面積をΔＳkとすると、
【数４】
0003751788.tif 000026
により、ＣＴ再構成画像を求めるプログラムを格納した記録媒体。  
【請求項４】
透過データから投影データｐ（ｒm，θn）を得る前記手段は、前記複数の検出器からの透過データをメモリに記憶してから、角度θに対応したｒ軸上の投影データを選択することにより、投影データｐ（ｒm，θn）を得ることを特徴とする、請求項３に記載のＣＴ再構成画像を求めるプログラムを格納した記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、ＣＴ画像再構成処理に関し、特に不等間隔に得られたデータに基づく高精度ＣＴ画像再構成処理に関する。
【０００２】
【技術的背景】
コンピュータＣＴ等の画像再構成の分野では、投影データの採取の高速化に対して工夫が行われている。この１例を図１で説明する。図１に示す構成では、検出器を物体を囲む円周４０上に多数設置し、Ｘ線源２０から扇（ファン）状に広がるビームを対象物を覆うように広角度で照射して、この範囲のデータを同時に収集している。これにより平行ビーム投影のように、Ｘ線源で一直線上に走査する操作は不要となり、投影データを迅速に収集することができる。画像再構成を行うためには、原理的には物体の任意の点に対して全ての方向からのＸ線ビームがあればよい。このため、図１に示した構成において、図示のようにＸ線源２０を対象物体３０の周りを１周させてデータを採取することにより、再構成に必要なすべてのデータを短時間で採取することができる。
【０００３】
得られた採取データから画像を再構成するためには、平行ビームと比較すると計算が複雑となる（例えば、斉藤恒雄「画像処理アルゴリズム」近代科学社  １９９３年  ｐ１１１〜ｐ１２２参照のこと）。
本発明の発明者は、図１のような構成でデータを採取した場合、得られたデータが不等間隔データであり、この不等間隔データを適切に処理することができれば、計算が容易で高精度に画像再構成が行えることを見いだした。
【０００４】
【発明が解決しようとする課題】
本発明の目的は、ＣＴ画像処理において、不等間隔で得られたデータから、高精度に画像を再構成できる処理の提供である。
【０００５】
【課題を解決するための手段】
上記目的を達成するために、本発明は、対象物を囲む円周上に設けたＸ線源および複数の検出器からの透過データによりＣＴ画像を再構成するＣＴ画像再構成システムであって、前記透過データから、角度θのｒ軸上における不等間隔の投影データｐ（ｒm，θn）を得る手段と、前記不等間隔の投影データｐ（ｒm，θn）から、直接的に極座標（ρ，θ）における等間隔の空間周波数スペクトルＧ（ρk，θn）を求める１次元フーリエ変換を行う１次元フーリエ変換手段と、前記空間周波数スペクトルＧ（ρk，θn）から、直接２次元フーリエ逆変換を行い、直交座標系（ｘ，ｙ）のＣＴ画像であるｆ（ｘl，ｙm）を求める２次元フーリエ逆変換手段とを備え、前記１次元フーリエ変換手段は、Δmを投影データの検出点が占める微少範囲とすると、
【数５】
0003751788.tif 000002
により空間周波数スペクトルを求め、前記２次元フーリエ逆変換手段は、空間周波数スペクトル標本値が極座標（ρ，θ）上でになう微少面積をΔＳkとすると、
【数６】
0003751788.tif 000003
によりＣＴ再構成画像を求めることを特徴とする。
これにより、不等間隔で得られたデータから、内挿や補間をすることなく、高精度に画像を再構成できる。
【０００６】
また、透過データから投影データｐ（ｒm，θn）を得る前記手段は、前記複数の検出器からの透過データをメモリに記憶してから、角度θに対応したｒ軸上の投影データを選択することにより、投影データｐ（ｒm，θn）を得ることができる。
対象物を囲む円周上に設けたＸ線源および複数の検出器からの透過データによりＣＴ画像を再構成するために、上述した手段としてコンピュータに機能させるプログラムを格納した記録媒体も本発明である。
【０００７】
【発明の実施の形態】
本発明の実施形態を、図面を参照して詳細に説明する。
まず、本発明の基礎である不等間隔標本化定理について説明する。
＜不等間隔標本化定理＞
説明のために、実数値信号ｘ（ｔ）を考え、そのフーリエ変換として定義される、フーリエ・スペクトルＸ（ｆ）は［−ｆmax，ｆmax］に帯域制限されていると仮定する。即ち、
【数７】
Ｘ（ｆ）≡０  ｜ｆ｜＞ｆmax    （１）
である。この仮定の下に、実数値信号ｘ（ｔ）の不等間隔サンプル値系列ｘ（ｔn）（ｎ＝−∞，・・・，＋∞）が与えられたとする。ただし、標本化間隔
【数８】
0003751788.tif 000004
の最大値Δtmaxは次の条件を満たすと仮定する。
【数９】
0003751788.tif 000005
このとき
【数１０】
0003751788.tif 000006
として、次の定理が成り立つ。
１．ｘ（ｔ）が式（１）のように帯域制限されているとき、式（２）を満たす間隔で不等間隔標本化してもｘ（ｔ）の情報は何も失われない。
２．実際、サンプル値系列ｘ（ｔn）（ｎ＝−∞，・・・，＋∞）から、ｘ（ｔ）は、次の内挿公式によって復元される。
【数１１】
0003751788.tif 000007
式（３）の内挿公式は、等間隔標本化の場合の内挿公式の一般化になっている。全てのｎに対してｔ＝ｎＴとおくと、
【数１２】
0003751788.tif 000008
となり、式（３）は、等間隔標本化における内挿公式となる。さらに、式（２）の条件を満たす不等間隔標本化は、等間隔標本化よりも冗長な標本化である。また、式（３）でｔ＝ｔnにとると、ｘ（ｔn）が他の標本値の線形結合で与えられ、標本値の独立性が失われることも分かる。
【０００８】
＜データ採取＞
さて、本発明のＣＴ画像処理に用いるデータの採取について、図２〜図５を用いて説明する。
図２は、本発明のおけるＣＴ画像のためのデータを採取するための構成を示している。図２において、Ｘ線源１２０および円周上に設置された複数の検出器１４０により、対象物１３０のＸ線の透過データを採取している。円周上に設置された複数の検出器１４０は、角度Δθの等間隔で設置されている。この様な構成で、対象物との相対回転運動により、Ｘ線源からのＸ線に対する対象物１３０の全周方向のＸ線透過データを検出器１４０から得る。図２ではΔθ＝２２．５°で検出器を設置した場合を示しているが、実際はもっと多数の検出器を設置する。なお、図２ではＸ線源１２０と検出器１４０とを同じ円周上に設置しているが、別々の円周上に設置し、Ｘ線源のみが回転するようにしてもよい。また、相対運動を行えばよいので、対象物１３０が回転するようにしてもよい。
【０００９】
図３に、例えば、対象物の周りを角度３０°づつ回転して、各検出器からデータを得る場合を示している。図３（ａ）〜図３（ｃ）では、それぞれ実線で示したＸ線は、Ｘ軸に対する角度θ＝０°のｒ軸上の投影データを得るものである。すなわち、図３のようにＸ線源を円周上で回転した場合、図４に示すようにｘ軸と垂直になるＸ線を受ける検出器からのデータを集めることにより、ｒθ=0軸上の投影データｐ（ｒm，θn）を求めることができる。
このようにして、Ｘ線源を対象物と相対的に一回転することにより、投影データｐ（ｒm，θn）を、各θnに対応したｒ軸に対して求めることができる。
【００１０】
相対的回転運動をΔθ（ＮΔθ＝３６０°）ごとに行い、θn＝ｎΔθの所でのｒ軸上のＭ個の投影データをそれぞれ、
【数１３】
0003751788.tif 000009
とする。すなわち、図５に図示するような、角度θのｒ軸上のＭs個の投影データｐを考えることができる。この場合、すべてのθnに対して、ｒm（ｍ＝０，１，・・・Ｍs−１）は
【数１４】
ｒm＝Ｒ0ｃｏｓ（ｍΔθ−π）    ｍ＝０，１，・・・Ｍs−１
となり、不等間隔となる。投影データｐとしては、中心部が粗く、周辺部ほど細かくなるｒ軸上の点のみのデータが採取されることを意味する。
【００１１】
図３に示したように、Ｘ線源を対象物の周囲を回転させて検出器によりデータを採取し、採取されたデータから図４に示すように平行のＸ線透過データを得ることを、図６および図７を用いて詳しく説明する。
図６は説明のために、検出器を３０°ごとに設け、また、Ｘ線源を３０°づつ回転させて、検出器からデータを収集することを想定して図示したものである。図６において、対象物を取り囲む円周上に設けられたＸ線源および検出器の位置により、各検出器から採集できるデータをＤ（φs，φd）とする。このとき、φsはＸ線源の位置（角度）、φdは検出器の位置（角度）を示す（０°≦φs≦３６０°，０°≦φd≦３６０°）。ｓおよびｄは、図６に示すように、φs＝０°のときｓ＝０、φs＝３０°のときｓ＝１となり、φs＝３３０°のときｓ＝１１となる。また、φd＝０°のときｄ＝０、φd＝３０°のときｄ＝１となり、φd＝３３０°のときｄ＝１１となる。
図６は、ｘ軸からφs＝０（即ち、Ｓ＝０）の所にＸ線源をおいた場合、Ｘ線の照射で得られる各検出器からデータＤ（φs，φd）を取り込むところを示している。３０°ごと（検出器１つごと）に採取された採取データは、Ｄ（φ0，φ0），Ｄ（φ0，φ1），・・・，Ｄ（φ0，φ11）となる。これを、円周上でＸ線源を３０°づつずらして、１周して各検出器からのデータを取り、メモリに記憶して置いて、図４に示したような平行のＸ線で照射した、ｘ軸からθn傾斜しているｒ軸上の投影データｐ（ｒm，θn）として取り出す。
図６のようにして採取したデータから投影データｐ（ｒm，θn）の取り出し方を、図７を用いて説明する。図７において、θn＝３０°（ｎ＝１）におけるｒ軸上の投影データｐ（ｒm，θ1）は、
【数１５】
ｐ（ｒ0，θ1）＝Ｄ（φ7，φ7）
ｐ（ｒ1，θ1）＝Ｄ（φ6，φ8）＝Ｄ（φ8，φ6）
ｐ（ｒ2，θ1）＝Ｄ（φ5，φ9）＝Ｄ（φ9，φ5）
ｐ（ｒ3，θ1）＝Ｄ（φ4，φ10）＝Ｄ（φ10，φ4）
ｐ（ｒ4，θ1）＝Ｄ（φ3，φ11）＝Ｄ（φ11，φ3）
ｐ（ｒ5，θ1）＝Ｄ（φ2，φ0）＝Ｄ（φ0，φ2）
ｐ（ｒ6，θ1）＝Ｄ（φ1，φ1）
の様に、採取データを取り出すことにより得ることができる。
【００１２】
＜空間周波数スペクトルの推定＞
さて、従来のＣＴ画像再構成では、得られた投影データｐから一次元フーリエ変換を行って、θn方向の空間周波数スペクトル情報を算出する必要がある。しかし、この様な不等間隔のデータから空間周波数スペクトル情報を直接算出することができず、このため、従来の手法では、不等間隔の投影データｐを、ｒmに対して補間、内挿し、ｒ軸に関して等間隔な射影データを得る必要がある。
本発明の画像再構成では、上述のようにして得られた不等間隔の投影データｐから、補間、内挿のような演算を行わず、空間周波数スペクトルを直接得ることができる。これを以下で説明する。
【００１３】
空間周波数スペクトルの推定は、２次元離散フーリエ逆変換で画像再構成を行うため、各角度ごとに、上述の採取データ（投影データ）から、１次元離散フーリエ変換を行い、空間周波数スペクトルを得ることである。
まず、θ＝θnのｒ軸上の投影データｐ（ｒm，θn）を図示した図８を参照して、検出点１個の占める微少範囲Δmは、
【数１６】
0003751788.tif 000010
としたとき、
【数１７】
0003751788.tif 000011
である。これを考慮すると、投影データｐ（ｒm，θn）に対する周波数空間（μ，ν）の各θn毎の極座標ρk軸（図９参照）上の周波数スペクトルの推定値Ｇ（ρk，θn）は、Ｍをρk軸上でのスペクトル推定を行う点の数、Δｒは周波数領域におけるスペクトルの推定値を求める間隔１／（ＭΔｒ）を定めるための定数とすると、
【数１８】
0003751788.tif 000012
により、求めることができる。なお、ＭおよびΔｒは自由に選定することができる。そのため、ＭΔｒを大きくとれば、ρ軸の計算間隔を細くとることになり、高精度な再生を行うことができる。
【００１４】
上述の式（５）の計算を行うときにまず、図１０に示すように、
【数１９】
0003751788.tif 000013
の範囲にある周波数スペクトルの推定値Ｇを、
【数２０】
0003751788.tif 000014
の範囲に移動させる。これは、実数データのフーリエ・スペクトルの負の周波数に対応する部分を正規の位置に戻すためである。このため、移動後の範囲では、上述の式（５）におけるｋの範囲は、
【数２１】
0003751788.tif 000015
となる。なお、Ｍ0については後で説明する。
【００１５】
＜画像再構成＞
上述のように、極座標に関して等間隔の空間周波数スペクトルが求まるので、直接２次元フーリエ逆変換を行うことができる。これを行うために、図６と同様に、まず、図１１に示すように、空間周波数スペクトルの標本値が（ρ，θ）座標上でになう微少面積ΔＳkを考えると、
【数２２】
0003751788.tif 000016
である。これを考慮して、直交座標（ｘ，ｙ）領域における等間隔の再生点（ｌΔｘ，ｍΔｙ）における再構成画像ｆ（ｘl，ｙm）は、極座標領域での等間隔周波数スペクトルの２次元離散フーリエ逆変換の式である
【数２３】
0003751788.tif 000017
により、直接求めることができる。ここで、Ｍ0は、空間周波数スペクトルの極座標（ρ，θ）領域の上限により、以下に説明するように定めることができる定数である。
式（５）は、上述の不等間隔標本化定理によるものであるから、式（２）で示した条件を満足しないとエリアシングが生じる。このため、式（５）において、空間周波数の最大値ρmaxは、
【数２４】
0003751788.tif 000018
の範囲内とする必要がある。
このとき、式（６）において、ｋに関する和を
【数２５】
0003751788.tif 000019
の範囲内で取らなければならない。すなわち、上述のＭ0は、
【数２６】
0003751788.tif 000020
とする必要がある。
【００１６】
＜数値シミュレーション結果＞
上述の画像再構成の処理が、従来の画像再構成の処理と比較して有効であることを、数値シミュレーション結果を用いて説明する。
数値シミュレーションに用いた物体の吸収率分布をｆ（ｘ、ｙ）として、以下に示すように、ｓｉｎｃ関数にテーパリング関数ｗ（ｘ、ｙ）を掛け、原点を中心とする半径Ｒ0の円の範囲に領域を制限した。
【数２７】
0003751788.tif 000021
これを用いて、ラドン変換を行うことで投影データを求める。そして、従来の画像再構成の処理（１次元のフーリエ変換、２次元のフーリエ逆変換）を用いて画像を求めることにより、数値シミュレーションを行う。この場合、従来の画像再構成の処理を行う場合、１次元、２次元の内挿や近似が必要となる。
【００１７】
また、上述の吸収率分布に対する投影データを、同様のラドン変換を行って求め、本発明の処理を用いて画像を得る画像再構成のシミュレーションを行う。これを従来の画像再構成処理により求めた画像と比較した結果を以下に述べる。
なお、検出器設置角度あるいは投影データ採取角度間隔は、等間隔Δθ＝１°単位とし、検出器半径はＲ0＝０．３ｍに、また、投影データのｒ方向採取点数Ｍsは５１２に固定していると仮定している。
【表１】
0003751788.tif 000022
なお、上述の表は、再生精度を相対ＲＭＳ（２乗平均）誤差で示している。
また、２−ＤＦＴ法は２次元離散フーリエ変換法、ＦＢＰ法はフィルタ補正逆投影法であり、これらの従来の画像再構成法に関しては、例えば斉藤恒雄著「画像処理アルゴリズム」（近代科学社  １９９３年発行）「第６章  計算機トモグラフィ」および「第７章  ラドン変換に基づく計算機トモグラフィの解析」を参照されたい。
このシミュレーションの結果に示されるように、本発明の画像再構成においては、従来の処理による画像再構成と比較すると、より高精度の画像再構成を行うことができる。
【００１８】
本発明の画像再構成システムを構築するシステムとしては、スタンド・アローンのコンピュータ・システムばかりではなく、複数のシステムから構成される例えばクライアント・サーバ・システム等でもよい。
本発明に関するプログラムを格納した記憶媒体から、プログラムをシステムで読み出して実行することにより、本発明の構成を実現することができる。この記録媒体には、フロッピー・ディスク、ＣＤ−ＲＯＭ、磁気テープ、ＲＯＭカセット等がある。
【００１９】
【発明の効果】
本発明のＣＴ画像再構成処理を用いると、データ補間などを行なわずに画像再構成できるので、高精度の画像再構成を行うことができる。
【図面の簡単な説明】
【図１】従来のファンビームによるＣＴデータ採取法を説明する図である。
【図２】本発明のＣＴデータ採取法を説明するための図である。
【図３】Ｘ線源と検出器との関係を示す図である。
【図４】本発明における不等間隔の採取データを示す図である。
【図５】本発明における座標関係を示す図である。
【図６】本発明のＸ線源と検出器によるデータ採取を説明する図である。
【図７】採取されたデータから投影データを得ることを説明する図である。
【図８】ｒ軸上の検出点１個の占める微少範囲を示す図である。
【図９】周波数空間における座標関係を示す図である。
【図１０】データの並び替えを説明する図である。
【図１１】極座標で周波数スペクトルの占める微少面積を示す図である。
【図１】
0003751788.tif 000027
【図２】
0003751788.tif 000028
【図３】
0003751788.tif 000029
【図４】
0003751788.tif 000030
【図５】
0003751788.tif 000031
【図６】
0003751788.tif 000032
【図７】
0003751788.tif 000033
【図８】
0003751788.tif 000034
【図９】
0003751788.tif 000035
【図１０】
0003751788.tif 000036
【図１１】
0003751788.tif 000037