JPWO2019230740A1 - オーバーサンプリングによる断層画像データの取得方法、取得装置、および制御プログラム - Google Patents

Abstract

計算負荷を増大させずに再構築される断層画像の再現性を高めるために、本開示の実施形態では、Ｎ個の検出素子で検出する撮像の際、（Ｎ＋ｎ）方向にわたるオーバーサンプリングにより検出する（Ｓ１２４、Ｓ１２６）。Ｎ×（Ｎ＋ｎ）個の要素をもつベクトルを求め（Ｓ１２８）、間引き順の数列｛ｋ｝３０に当たる計ｎ×Ｎ個の要素を除去するベクトルデシメーションＳ１３０を実施する。離散ラドン逆変換ステップＳ１３２で対応する離散ラドン逆変換行列ＷＳＱ−１４０を作用させ、画像データ生成Ｓ１３２でデベクトル化して断層画像データを取得する。オーバーサンプリングを利用すれば離散ラドン逆変換行列ＷＳＱ−１が得られるため、断層画像が逐次近似によることなく行列の演算により得られる。

Description

本開示はオーバーサンプリングによる断層画像データの取得方法、取得装置、および制御プログラムに関する。さらに詳細には本開示は、オーバーサンプリングに基づく離散ラドン逆変換の効率の良い代数的手法を組み込んだ断層画像データの取得方法、取得装置、および制御プログラムに関する。

Ｘ線、ガンマ線、光波などの透過性の電磁波や粒子線、または地震波やその他の任意の透過性を示す波または粒子を利用する断層撮像法（トモグラフィー）が実用化されている。その代表例が、医療、工業用途で用いられるＸ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)である。また、近年では可視光を用いた光学ＣＴ装置が開発されており、がんの放射線治療を行う前に、放射線線量の三次元分布検証のため、放射線照射後のゲル線量計を対象物とした撮像法が研究されている。これらの撮像手法で対象物の断層画像が得られる数学的原理はラドン変換である。この場合、対象物内部の各部の減衰率を撮像情報として外部に取出す処理がラドン変換によって記述され、対象物の断層情報を外界に撮像情報から再構成する処理がラドン逆変換によって記述される。すなわち、対象物で減衰した透過後の波または粒子の流れを照射方向または検出方向を変えて各位置で検出する操作がラドン変換に対応する。逆にその検出した方向別位置別の強度情報から対象物内部の各部の減衰率を推定して画像を再構築する操作がラドン逆変換に対応する。連続座標や連続方位でのラドン変換／ラドン逆変換は、現実の画像では有限画素数を扱うためサンプリングされた座標や方位で実行される。サンプリングした座標や方位の下での画像再構成のためのラドン逆変換（離散ラドン逆変換）は多元の連立一次方程式を解くことつまり行列の演算に還元される。これらの数学的側面については標準的教科書に解説されている（例えば非特許文献１、２）。

反復処理を行う逐次近似法（iterative reconstruction, ＩＲ）によってよりノイズやアーティファクトが軽減される手法も注目されている。その反復は、例えば１００回またはそれ以上繰返すことにより解を得ている。このように代数的手法に頼って多元の連立一次方程式の近似解を得ることによって再構成像を得る手法はＡＲＴ（Algebraic Reconstruction Technique）と呼ばれている。ＡＲＴは、例えば、Maximum Likelihood - Expectation Maximization (ML-EM）法、Additive - Simultaneous Iterative Reconstruction Techniques (ASIRT) 法、Multiplicative - Simultaneous Iterative Reconstruction Techniques (MSIRT) 法、と呼ばれる手法を含んでいる。従来のＡＲＴでは、計算量が膨大となり多量の計算資源を必要とする。また、イタレーションを重ねても、真の画像に対応する解に近づかず、局所解に捕獲されてしまうことも多い。他方、代数的手法にこだわらずに計算機での近似計算の効率を高め、現実的な計算資源で対処する手法も模索されてきた。この一例がＦＢＰ（Filtered Back Projection）である（例えば特許文献１）。ＦＢＰは代数的手法に比べ計算量が少ないことから現在Ｘ線ＣＴスキャナーにおいて最も普及している。

断層撮像法の発展経過では、初期から今日まで撮像画素数増大の要請が継続しており、ＦＢＰの可能性が探求されてきた。しかしＦＢＰは、ＩＲに比して種々のアーティファクトが顕われやすく再構成像の品質では一般に不利である。計算能力の発展が続く今日では、再構成像の品質で本来優位なＡＲＴも再び注目されている（例えば特許文献２）。

国際公開第２０１４／０２１３４９号 国際公開第２０１３／１０５５８３号

Aviash C. Kak and Molcolm Slaney, "Principles of Computerized Tomographic Imaging," IEEE Press (1988) Stanley R. Deans, "The Radon Transform and Some of Its Applications," Krieger Publications (1983) Robert L. Siddon, "Fast calculation of the exact radiological path for a three-dimensional CT array," Medical Physics, vol.12, no.2, pages 252-255 (1985) Sunnegardh, J. and Danielsson, P.-E. "A New Anti-Aliased Projection Operator for Iterative CT Reconstruction," 9 th International Meeting of Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, pages 125-127 (2009)

画像の精細度の向上のための画素数増大への要請に加え、再構成像の品質に対する要請も弱まる気配がない。撮像画素数が増大した状況で従来のＡＲＴを採用しても、解くべき連立一次方程式の元（決定すべき変数）の数が増大してしまい膨大な計算資源が必要となり、その処理に長い時間を要す。このため従来のＡＲＴでの画像データ取得処理に実用性を求めるには、現状では計算能力の向上を待つしかない。特に、逐次近似に伴う反復計算を撮像の毎に再実行するような従来のＡＲＴの処理に頼る限り、例えば断層画像を多数得て三次元のボリューム像を得ることは現実的ではない。

本開示は上記問題の少なくともいくつかを解決することを課題とする。本開示は、像の再構成を担う離散ラドン逆変換のために代数的厳密解の裏付けをもつような新たなサンプリング手法を採用することにより、高品質な断層撮像法に高い実用性をもたらすものである。

本発明者は、従来のＡＲＴの問題点がサンプリング（標本化）といった離散化手法に起因していることに気づいた。従来のＡＲＴでは、多元の連立一次方程式に重複する式が多数含まれてしまい、元（決定すべき変数）の数に対し利用できる方程式の数が実質的に足りなくなっている。このために逐次近似が必要となってプロセッサーの計算能力による解決が指向されているのである。これに対し、本発明者は、離散化手法から再検討することにより、従来にないサンプリングと適切な処理とを組み合わせる上記課題への新たな解決手法を完成させた。

すなわち、本開示のある態様においては、少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置の検出範囲に対象物を配置するステップと、照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない（Ｎ＋ｎ）方向（ｎは１以上の整数）にわたって前記検出方向を変更して繰り返す検出ステップと、各行および各列を各検出方向および各検出素子に対応させた（Ｎ＋ｎ）行Ｎ列のオーバーサンプル・サイノグラムをベクトル化したものに相当するＮ×（Ｎ＋ｎ）個の要素をもつ第１ベクトルを、前記検出動作における前記検出装置による検出信号から得るベクトル化ステップと、前記第１ベクトルから、間引き順に当たる計ｎ×Ｎ個の要素を除去して、残りのＮ×Ｎ個の要素をもつ第２ベクトルを得るベクトルデシメーションステップと、前記第２ベクトルに離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第３ベクトルを得る離散ラドン逆変換ステップと、前記第３ベクトルをデベクトル化することにより、前記対象物に対して静止した画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップとを含む断層画像データ取得方法が提供される。本開示のかかる態様においては、照射（投影）方向についてオーバーサンプリングとデシメーション（間引き）を組み合わせることにより、対称性の破れを意図的に生じさせて代数的手法による厳密な解の算出を可能性にしている。

また、本開示は、断層画像データ取得装置の態様としても実施することができる。同様に、本開示は、断層画像データ取得装置のための制御プログラムの態様でも実施することができる。

オーバーサンプリングまたはオーバーサンプル（oversample）とは、デジタル信号処理理論のサンプリングの概念に基づいており、基本となるサンプリングよりも高い密度でサンプリングをする処理をいう。特に本開示の各態様におけるオーバーサンプリングは検出方向について有効な検出素子の数Ｎを超えた数となるＮ＋ｎ（ｎは１以上の整数）の検出方向のサンプル点を互いに重複しないように選択して必要な角度範囲をサンプリングすることを含む。なお、以下特に断りのない場合の検出素子の数は有効なものを指すこととする。デシメーション（間引き処理）とは、オーバーサンプルされて得た行列の要素や対応するベクトルの要素のうち、少なくともいずれかをその後の計算処理に用いない（捨てる）処理を指している。これらを含め本出願においては、不明瞭にならない限り本開示の属する技術分野の慣用に従う用語法を採用することがある。波は電磁波、音波といった伝播による減衰を観念できる任意の波動を含み、光波（赤外、可視、紫外）、振動（地震など）を含む。粒子は、流れとなった粒子線において伝播による減衰を観念できる任意の粒子を含み、中性子線やその他の放射線も含む。電子や光子のように、量子力学的に波動性と粒子性を兼ね備えるものも、これらの波または粒子のいずれかまたは両方に含まれる。同様に、古典力学的に波動や粒子と捉えられるものも、これらの波または粒子のいずれかまたは両方に含まれる。これらの波または粒子は、それらの人為的な放出源となる装置（照射装置）により生成されるものであっても、また、例えば宇宙線などの自然由来で生成されるものや、例えば対象物の各部から放射される放射線などの対象物の性質として生じるもの（本出願書面では「照射装置によらずに生じた波または粒子」と呼ぶ）であっても本開示の実施に用いることができる。

本開示のある態様では、相対的に少ない計算資源を用いる場合でも高い精細度で高品質な再構成像を高速に実現することができる。

図１は従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。 図２は従来および本開示の実施形態の離散ラドン逆変換の代数的手法の原理を説明する説明図であり、図２Ａ〜Ｃは従来の手法を示し図２Ｄは本実施形態の工夫を示す。 図３は本開示の実施形態において画像が取得される一例の断層撮像装置において、画像が取得される対象物の切断面における平面配置を含む概略構成を説明するための模式図である。 図４は本開示の実施形態の断層撮像装置のうち制御装置において実施される処理を概観するフローチャートである。 図５は本開示の実施形態の断層撮像装置における処理のフローチャート（図４）の各段階で得られる行列やベクトル、または画像についての変遷を模式的に示す説明図である。 図６は本開示の実施形態における行列デシメーション処理を示すフローチャートであり、図６Ａは乱択アルゴリズムによるもの、図６Ｂはグラム・シュミットの直交化法を採用するものである。 図７は、本開示の実施形態におけるブロックによるデシメーション処理を示す説明図であり、図７Ａは、オーバーサンプル・システム行列Ｗ_ＯＳを画像により示し、図７Ｂは、ブロックでの間引き処理によって得た正方システム行列Ｗ_ＳＱを同様に示したものである。 図８は、本開示の実施形態におけるブロックによるデシメーション処理の成果を示す説明図であり、図８Ａは離散ラドン逆変換行列Ｗ_ＳＱ ^−１であり、図８Ｂは、離散ラドン逆変換行列Ｗ_ＳＱ ^−１を正方システム行列Ｗ_ＳＱに左から作用させて得られる行列を示す。 図９は断層撮像装置によって本実施形態を実装する場合における制御装置のいくつかの特有の機能部と特有のデータ格納部を示すブロック図である。 図１０は本開示の実施形態のための実施例において、検証のために用いた検証用データ（数値ファントム）を示す図（図１０Ａ）および対応するオーバーサンプル・サイノグラム（図１０Ｂ）である。 図１１は本開示の実施形態のための実施例において得られた再構成画像であり、本実施形態により求めたもの（図１１Ａ）、従来のＦＢＰ法においてＨａｎｎフィルターを用いたもの（図１１Ｂ）、および従来の代数再構成法として逐次近似計算（ＭＬ−ＥＭ法）を３０回行って得たもの（図１１Ｃ）である。また、図１１Ｄは比較のための基準とする単純な逆ラドン変換による再構成画像である。 図１２は本開示の実施形態における追加の画像再構成手法において、マウス標本を利用して取得したＸ線ＣＴのデータによる再構成処理を行った事例を説明する説明図であり、図１２Ａはマウス標本の写真、図１２Ｂは擬似逆行列Ｗ_ｐ ^−１を利用した再構成画像、図１２ＣはＦＢＰによる再構成画像、図１２ＤはＭＬ−ＥＭによる再構成画像、図１２Ｅは単純な逆ラドン変換による再構成画像である。 図１３は本開示の実施形態における追加の画像再構成手法において擬似逆行列Ｗ_ｐ ^−１により生成した再構成画像を含んでおり、図１３Ａ〜図１３Ｊは各図に付す値にトレランス値を設定したものであり、図１３Ｋは、特異値を大きさ順に並べたグラフである。 図１４は本開示の実施形態における追加の画像再構成手法において、コンクリート標本を利用して取得した中性子線ＣＴのデータによる再構成処理を行った事例を説明する説明図であり、図１４Ａはコンクリート標本の写真、図１４Ｂは擬似逆行列Ｗ_ｐ ^−１を利用した再構成画像、図１４ＣはＦＢＰによる再構成画像、図１４ＤはＭＬ−ＥＭによる再構成画像、図１４Ｅは単純な逆ラドン変換による再構成画像である。 図１５は本開示の実施形態における追加の画像再構成手法において擬似逆行列Ｗ_ｐ ^−１により生成した再構成画像を含んでおり、図１５Ａ〜図１５Ｊは各図に付す値にトレランス値を設定したものである。 図１６は、本開示の実施形態における追加の画像再構成手法における擬似逆行列Ｗ_ｐ ^−１の例であり、図１６Ａは、図７Ｂの正方システム行列Ｗ_ＳＱに対して求めた擬似逆行列Ｗ_ｐ ^−１、図１６Ｂは、擬似逆行列Ｗ_ｐ ^−１を正方システム行列Ｗ_ＳＱに左から作用させて得られる行列を示す。

以下図面を参照し、本開示に係る断層画像撮像の実施形態を説明する。全図を通じ当該説明に際し特に言及がない限り、共通する部分または要素には共通する参照符号が付される。また、図中、各実施形態の要素のそれぞれは、必ずしも互いの縮尺比を保って示されてはいない。

１．原理
図１は、従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。この図は、波または粒子を照射し、透過後の強度データを用いて断層画像を再構築する処理の本質を示す初等的なものである。図示される縦２×横２に並ぶ正方形領域はそれぞれが断層画像の画素を意図したもので、ａ_１１〜ａ_２２の４つの値を持ちうる。これらの値は、各画素の位置での対象物によるビームの減衰の程度を示し、決定されるべき値である。複数の画素を通るように描かれた白抜き矢印は波または粒子のビームの照射と検出素子の位置に対応している。矢印が指し示す先に記す数値は、検出装置の検出素子がその矢印の先に位置している際に得られる値（例えば減衰量）を意味しており、ここでは通った画素の値の合計が数値で示されている。つまり、この矢印が波または粒子のビームを表し、矢印が通る画素と矢印の向きで照射および検出（以下「照射」とのみ記す）を示している。矢印の先の合計値は、その矢印の先に検出装置の検出素子がある場合の、各検出素子で検出された減衰成分に対応している。矢印の指し示す値つまり検出装置から得られる値のみによってａ_１１〜ａ_２２の４つの値を特定することが、検出装置の検出素子の示す値に基づいて離散ラドン逆変換を実行することにより断層画像を再構成する処理に対応する。

図１に示された例では、ａ_１１〜ａ_２２のうち矢印それぞれがつないでいるものが線形和を与えるため、連立方程式に書き直せる。その連立方程式を解くと、
ａ_１１＝４−ａ_２１
ａ_１２＝３＋ａ_２１
ａ_２１＝任意
ａ_２２＝７−ａ_２１
が求まる解となる。つまり、元（決定すべき変数）の数が４であるのに対し、矢印の示す照射によって導き出せる連立方程式の実質的な数は３である。ａ_２１にどのような値を入れても上記連立方程式は成立してしまうことから、ａ_１１〜ａ_２２の４つとも不定（indeterminate）となり値を特定できない。

図２は、従来および本実施形態の離散ラドン逆変換の代数的手法の原理を説明する説明図である。図２Ａは図１と同じもの、図２Ｂ、Ｃは別の従来の手法を例示するもの、そして図２Ｄが本実施形態の工夫を反映したものである。図２Ａ〜Ｄの解を表１にまとめている。

表１に示すように、図２Ａ（図１）の場合と同様、図２Ｂ、Ｃでもａ_１１〜ａ_２２の４つの値を特定しえない。つまり対角方向の照射を用いて変数の組み合わせを変更して連立方程式を作ったとしても、図２Ｂ、Ｃではａ_１１〜ａ_２２の４つの値が特定できるとは限らない。他方、図２Ｄの矢印で示すような照射をした場合、ａ_１１〜ａ_２２の４つの値を一意に決定することができる。注目すべきは、図２Ｄの矢印の組み合わせが図２Ａ〜Ｃにある個別の矢印（すなわち照射方位と検出素子）を取捨選択して組み合わせたものであることである。図２Ｄでは、解くべき連立方程式の数を示す矢印の数が増えているわけでもない。それにも関わらず図２Ｄでは４つの値が特定される。このような違いは対称性の現われと本発明者は考えている。つまり、図２Ａ〜Ｃでは検出方向または検出素子について対称性が高いため解が不定となるのに対し、図２Ｄでは、対称性の破れ（symmetry breaking）が生じているため代数的に値が特定できている。この違いが示唆するのは、照射の段階で適切な工夫をすれば、現実の断層画像の再構築処理のためのより一般の形をもつ多元連立一次方程式でも代数的手法によって解を求められる余地があることである。本願の発明者は、照射の段階（離散ラドン変換）においてオーバーサンプリングを行いさらにデシメーションを組み合わせれば、対称性の破れを意図的に実現することができて代数的手法によって不定ではない厳密解が得られることに気づいた。

ここで、本実施形態における例示の幾何学的配置を説明する。図３は本実施形態の画像が取得される一例の断層撮像装置１００において、画像が取得される対象物の切断面における平面配置を含む概略構成を説明するための模式図である。対象物（図３には示さない）に対し静止した空間の一平面にＮ画素×Ｎ画素を持つ画素の群が定義される。通例この整数Ｎは検出装置２０の検出素子２２の数と一致させている。このため、本実施形態において、整数Ｎを特に解像度Ｎと呼ぶこともある。波または粒子のビーム４は、その平面に含まれており、照射装置１０から照射される。検出装置２０と照射装置１０は、相対配置が固定されたまま対象物に対して回転することができる。このような制御を含む装置の動作や計算処理のために、断層撮像装置１００はコンピューター機器でありうる制御装置１５０を備えている。この回転は、上記画素の群やそこに静止している対象物に対する相対的なものである。このため必ずしも検出装置２０と照射装置１０のみが回転されるとは限らず、検出装置２０と照射装置１０が設置床面（図示しない）に対し固定されていて対象物が、画素の群とともに回転される場合もある。

照射装置１０により波または粒子のビーム４を照射し検出装置２０によりそのビーム４を対象物の反対側で検出する動作は、従来の手法では検出方向θのスキャン範囲を画素数Ｎによって分割した刻みで検出方向についてサンプリングされる。図３のようなビーム４が平行なものでは、検出方向θのスキャン範囲は例えば０°〜１８０°とされ、従来の手法では検出装置２０からの検出信号が（１８０／Ｎ）°の刻みでサンプリングされる。なお、本実施形態を含めて、断層撮像装置１００の検出方向θの典型的なスキャン動作は、一定速度で増加または減少させる定速回転をしながら、検出装置２０の値を取得するタイミングを電気的に制御することによりサンプリングが行われる。ただし、本実施形態はこのようなものに限定されるものではない。

ビーム４は上記座標にて全ての画素範囲をカバーする広がりを持ち、各検出素子２２に向かう部分はある幅τを持っている。当該検出素子２２に対するビーム４についてのその画素のもつ寄与は、典型的にはビーム４のうち、検出素子２２の一つに向かうものがもつ幅τと各画素（縦横それぞれδ）の重なりである。図３では一つの画素についてこの重なりの様子を模様で示している。Ｎ画素×Ｎ画素の全ての画素について、検出方向θ別、検出素子２２の別に寄与が定まる。この寄与を数値として保持しているのがシステム行列Ｗである。システム行列Ｗを用いると、照射と画像再構成のうち、照射を表す離散ラドン変換は、行列の積により、
Ｘ′＝ Ｗ Ｘ （１）
と表現される。ここで、システム行列Ｗは、総画素数に応じた空間についてのサンプリング点の並びを行方向に、検出方向の刻みと検出素子とに応じた照射および検出についてのサンプリング点の並びを列方向にとるのが通常である。従来の手法を図３の画素数などの設定に適用した場合には、システム行列Ｗが行方向にも列方向にもＮ×Ｎ個の要素を持つ正方行列となりその要素の総数はＮ^４となる。これは、画素がＮ画素×Ｎ画素、検出方向θのスキャン範囲がＮサンプル、検出素子２２の数がＮ素子であるためである。なお、画質を高めるためにＮの値は時代とともに増大しており、典型例では１０２４すなわち１０^３程度とされ、従来のものでは２０００程度が最大値であるが、本実施形態の場合にもＮの値が得に限定されるものではない。Ｎが１０^３程度のものではシステム行列は行方向列方向ともに１０^６程度の規模、要素数は１０^１２程度となる。

式（１）のＸは各画素の位置での減衰（吸収、散乱も含む）を表す像である減衰像を列ベクトルとして表現したものである。実際の空間での減衰像は図３に示したＮ画素×Ｎ画素の二次元画像であるが、システム行列Ｗとの演算のために、行方向を総画素数Ｎ×Ｎになるようにベクトル化したものがＸである。式（１）において減衰像Ｘにシステム行列Ｗが作用して得られる左辺の列ベクトルＸ′は、サイノグラムと呼ばれる像に対応している。システム行列Ｗにて列方向を検出方向のサンプリング点と検出素子の並びとしたことに対応して、列ベクトルＸ′も検出方向のサンプリング点と検出素子それぞれについて得られる検出強度を示している。なお、サイノグラムは検出素子の軸と検出方向刻みの軸との２次元に列ベクトルＸ′の要素を並べ替えたものである。以上のように、式（１）は照射と検出を表現したものである。図２に示した図では、列ベクトルＸはａ_１１〜ａ_２２の各値に、また列ベクトルＸ′は各矢印の先に付した値にそれぞれ対応し、システム行列Ｗは、各矢印と各画素の対応関係を反映するものといえる。

断層画像の再構成とは、式（１）が実現していることを前提として、システム行列Ｗとサイノグラムの列ベクトルＸ′から列ベクトルＸを算出し、列ベクトルＸから減衰像を得る処理である。つまり断層画像の再構成は形式的に次式により表現される。
Ｘ＝ Ｗ^−１ Ｘ′ （２）
ここで、Ｗ^−１はＷの逆行列であり、式（２）は式（１）の両辺にＷ^−１を作用させて左右辺を入れ替えたものである。実際には、正方行列であるＷが常に逆行列Ｗ^−１をもつ（正則である）とは限らない。図２Ａ〜Ｃでａ_１１〜ａ_２２の各値が定まらないのは、それらの照射と検出の仕方ではＷ^−１が存在しないことに対応し、他方、図２Ｄでａ_１１〜ａ_２２が決定できるのはＷ^−１が存在することに対応する。つまり、図２Ｄにて例示した対称性の破れを用いるという着想の実現可能性は、逆行列Ｗ^−１が存在するような正方行列であるシステム行列Ｗを十分な確かさで決定できることの可能性と等価である。本実施形態では、逆行列Ｗ^−１をもつシステム行列Ｗを確実に決定できる手法が提案され、さらにそのための照射、検出、および断層画像再構成のための十分に実用性のある具体的構成も提供される。なお上述した図３に基づく説明は平行ビームを照射するものであるが、本実施形態はファンビームやコーンビームといった平行ではないビームを利用するような照射装置１０や検出装置２０を採用する場合にも当業者には明らかなわずかな変更によって適用することができる。また、図３の画素数などの設定以外においても、本実施形態の手法を適用することができる。なお、式（１）、（２）を、数学的に等価な他の形式により表現することもでき、たとえば、システム行列Ｗのために、上述したものとは異なり、総画素数に応じた空間についてのサンプリング点の並びを列方向に、検出方向の刻みと検出素子とに応じた照射および検出についてのサンプリング点の並びを行方向にとる。そしてそれに対応させて、式（１）のＸを各画素の位置に応じて行ベクトルとして表現することもできる。その場合式（１）、（２）に相当する式は、それぞれの各辺に転置操作を施したものとなる。この転置操作では行列やベクトルの行と列を入れ替え、ベクトルに作用させる行列が左からではなく右から作用させる順序に変わるという表現上の影響がある。ベクトルが「列または行ベクトル」のように択一的に記載され、対応する行列の要素の並びが「列方向または行方向」のように択一的に記載されているときは、これらの組み合わせの選択は、それぞれの択一的記載における同順のものを組み合わせて選ぶこととする。この表現の選択は本出願の本質に影響を与えないため、本願においては特に断りの無い場合には式（１）のようにサイノグラムをベクトル化したＸ等を列ベクトルにより表現し、システム行列等の行列は左から作用させる表記にもとづいて説明する。

２．具体的構成
本実施形態では、最終的に式（２）が成立するよう、逆行列Ｗ^−１を持ちうる正則な正方行列であるシステム行列Ｗを確実に決定できる。なお、本実施形態での区別のため，正方行列のシステム行列Ｗを、特に正方システム行列と呼びＷ_ＳＱと記す。

２−１．オーバーサンプリングとデシメーション
本実施形態で正則な正方システム行列Ｗ_ＳＱが得られるのは、オーバーサンプリングとデシメーション（間引き処理）を行うためである。図４は、これらを含む本実施形態の断層撮像装置１００のうち制御装置１５０において実施される処理を概観するフローチャートであり、図５は当該フローチャートの処理にて各段階で得られる行列やベクトル、または画像についての変遷を模式的に示す説明図である。図６は行列についてのデシメーション処理を示すフローチャートであり、図６Ａは乱択アルゴリズムによるもの、図６Ｂはグラム・シュミットの直交化法を採用するものである。

２−２．処理の概要
本実施形態において実行される処理を概観すると、まず装置の構成が決定される（図４、Ｓ１０２）。後の撮像でＮ画素×Ｎ画素（Ｎは正の整数）の二次元画像を再構成する場合、検出装置２０の検出素子２２の数もＮ個とされる。検出方向θについてのサンプリングは、従来の手法では、上述したように、例えば０〜１８０°のような必要な範囲において重複しないＮ方向となるように、例えば（１８０／Ｎ）°刻みでサンプリングが行われる。これに対し本実施形態の検出方向のサンプリングは、必要な範囲において重複しないＮ＋ｎ方向（ｎは１以上の整数）でなされ、例えば（１８０／（Ｎ＋ｎ））°刻みで行われる。このｎは、装置構成の決定Ｓ１０２の段階で検出素子の素子数（断層画像の一辺の画素数）Ｎに応じて適宜に決定される。このように、Ｎ方向に代えてＮ＋ｎ方向で検出方向のサンプリングを行うことを、オーバーサンプリングと呼ぶ。またｎを特にオーバーサンプリング量とも呼ぶ。

対象物２から断層画像を取得して適切な再構成像を得るためには、オーバーサンプリングに合わせシステム行列を決定しなくてはならない。オーバーサンプリングの場合のＮ＋ｎ個の検出方向で図３に示した従来のシステム行列と同様の処理を行ってシステム行列を生成する（Ｓ１１２）。ここで求まるシステム行列は、検出方向と検出素子それぞれに対応させて列方向にＮ×（Ｎ＋ｎ）、断層画像の画素数に対応させて行方向にＮ×Ｎだけの要素をもつ。この非正方の行列を本実施形態においてオーバーサンプル・システム行列と呼びＷ_ＯＳと記す。図５に、このオーバーサンプル・システム行列Ｗ_ＯＳを模式的に示している（符号２１２）。このオーバーサンプル・システム行列Ｗ_ＯＳ２１２はオーバーサンプル・システム行列格納部５０に格納される。次の行列デシメーション処理Ｓ１１４において、オーバーサンプル・システム行列Ｗ_ＯＳ２１２から行列の適切な位置（間引き順）にあるＮ×ｎ個の行ベクトルの要素が除去され、正方システム行列Ｗ_ＳＱ（図５、符号２１４）が生成される。この正方システムＷ_ＳＱ２１４は正方システム行列格納部６０に格納される。その際には、間引き順を与える数列｛ｋ｝（以下「間引き数列｛ｋ｝」と呼ぶ）も間引き数列格納部３０に格納する。この間引き数列｛ｋ｝が適切であれば、正方システム行列Ｗ_ＳＱは正則行列となる。正方システム行列Ｗ_ＳＱからは、離散ラドン逆変換行列算出ステップＳ１１６にてその逆行列Ｗ_ＳＱ ^−１が算出され、離散ラドン逆変換行列格納部４０に格納される。本実施形態では逆行列Ｗ_ＳＱ ^−１を離散ラドン逆変換行列とも呼ぶ。ここで逆行列Ｗ_ＳＱ ^−１を算出するための処理には、例えば余因子法や掃き出し法など任意の手法を採用することができる。

一方、対象物からの断層画像を得るには、照射装置１０と各検出素子２２をもつ検出装置２０との間に、対象物２を配置する（Ｓ１２２、図５）。オーバーサンプリングでの照射および検出のために、照射装置１０からのビーム４を照射しながら、ｎを１以上の整数として重複しないＮ＋ｎ方向で検出が行われる（Ｓ１２４）。典型的な実装形態では、このビーム４は連続して照射し続けていて、対象物２と照射装置１０および検出装置２０の向き（検出方向）はなめらかに動きつつ、検出装置２０の検出素子２２からの信号が電気的にサンプリングされる。重複しないＮ＋ｎ方向での照射および検出の典型例は、（１８０／（Ｎ＋ｎ））°刻みのタイミングでのサンプリングである。得られる信号は、必要に応じて機器の物理的特性について補正やスケールの換算などの処理と配列順序の整理などをすれば、サイノグラムを生成しうるデータに整理される（Ｓ１２６）。整理したものは、検出方向と検出器の位置とに対応することにより検出方向についてオーバーサンプルしたものであるため、本実施形態ではこれをオーバーサンプル・サイノグラム（図５、符号２２６）と呼ぶ。このオーバーサンプル・サイノグラム２２６を得る処理が、照射および検出の処理であり、離散ラドン変換である。このオーバーサンプル・サイノグラム２２６からは、オーバーサンプル・システム行列Ｗ_ＯＳ２１２の列方向の並びに合わせるような順序で要素を並べ替えた列ベクトルを得ることができる（ベクトル化Ｓ１２８）。このベクトルを本実施形態では第１ベクトル２２８（図５）と呼ぶ。第１ベクトル２２８は、オーバーサンプリング・サイノグラムの値を持つのでＮ×（Ｎ＋ｎ）要素を持つ。第１ベクトル２２８は、式（１）と同様に、仮に断層画像が得られていて、それをベクトル化したもの（Ｎ×Ｎ要素の列ベクトル）にオーバーサンプル・システム行列Ｗ_ＯＳ２１２を作用させたとしたならば得られることとなる列ベクトルといえる。第１ベクトルはＸ_ＯＳ′とも記す。なお、照射および検出Ｓ１２４からベクトル化Ｓ１２８までは、信号のアナログデジタル変換による量子化、デジタル信号の処理や、必要に応じて一時的にデータを保存して処理することもできる。また、オーバーサンプル・サイノグラムは、それを生成しうるデータが得られれば十分であり、必ずしもその表示を行う必要はなく、オーバーサンプル・サイノグラムといえる明示的なデータを経由する必要もない。オーバーサンプル・サイノグラムをベクトル化したものに相当する第１ベクトルが検出装置２０からの検出信号に基づいて得られれば、照射および検出Ｓ１２４以降ベクトル化Ｓ１２８まで等価な任意の処理を行う任意の処理も本実施形態に含まれている。

正方システム行列Ｗ_ＳＱを得たのと同一の間引き順に従って、第１ベクトルをＸ_ＯＳ′の要素を間引くベクトルデシメーション処理（Ｓ１３０）を行えば、正方システム行列Ｗ_ＳＱとの間で式（１）の関係を満たすことができるＮ×Ｎ要素をもつ列ベクトルである第２ベクトル２３０を得ることができる。第２ベクトルはＸ′とも記す。第２ベクトルＸ′は正方システム行列Ｗ_ＳＱとの間で式（１）の関係を満たすため、正方システム行列Ｗ_ＳＱの逆行列である離散ラドン逆変換行列Ｗ_ＳＱ ^−１を用いて式（２）の関係が成り立ち、離散ラドン逆変換処理Ｓ１３２によって再構成像を与える第３ベクトル２３２が得られる。この第３ベクトル２３２はＸとも記す。この第３ベクトルＸの並びをベクトル化の逆操作によってＮ画素×Ｎ画素の画像に戻すデベクタライズ処理により再構成像を画像化することができる（Ｓ１３４）。これで１スライス分の断層画像が再構築されるので、必要に応じてスライス位置を変更した処理を実行することができる。

本実施形態の処理は図４、５においてＡとＢに示す２系統に分かれている。通常再構成像を撮像するユーザーは、照射および検出ステップＳ１２４〜画像化ステップＳ１３４に至る処理の系統Ｂを実行する。つまりユーザーは、装置構成の決定Ｓ１０２、システム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａを必ずしも実行する必要はない。システム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａが必要となるのは、断層撮像装置を作製したり、設置やメンテナンスしたりする場面で、いわば機器メーカーが対応するような処理である。予め決められたオーバーサンプリング量ｎに応じて間引き数列｛ｋ｝および離散ラドン逆変換行列Ｗ_ＳＱ ^−１が得られていれば、再構成像を撮像しようとするユーザーは自らの測定のためにはそれらを利用するだけでよい。系統Ａからの適切な機器情報についてのデータが得られていれば、それを使用すれば系統Ｂの範囲の処理は十分実行できる。

特筆すべきは、系統Ｂのうち離散ラドン変換Ｓ１２６〜画像化ステップＳ１３４が行列やベクトルの順序の入れ替えや行列積の演算であり、例えば逐次近似といった計算処理量が問題となる処理がその範囲に含まれていないことである。本実施形態において、ユーザーの通常の撮像において利用される系統Ｂの処理は、系統Ａの処理により準備された間引き数列格納部３０、離散ラドン逆変換行列格納部４０を使うだけで十分に軽い計算負担で実行可能であり、この点は、本実施形態に高い実用性をもたらすものといえる。断層画像を取得するという最も頻度の高い処理の計算負担がごく少ないためである。特に、三次元ボリューム像を得るために多数のスライスの撮像を繰り返す用途には本実施形態の手法が極めて有利である。その繰り返しに系統Ａの処理が不要なことに加え、系統Ｂの離散ラドン変換Ｓ１２６〜画像化ステップＳ１３４はＧＰＧＰＵ（General-purpose computing on graphics processing unit）などのメニーコアタイプのプロセッサーでの計算処理に適するため、今後の処理速度の向上も期待できるからである。

また、幾何学的な側面を反映するオーバーサンプル・システム行列Ｗ_ＯＳ２１２や正方システム行列Ｗ_ＳＱ、システム行列Ｗ_ＳＱの逆行列Ｗ_ＳＱ ^−１がシステム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａのみに含まれており、測定のための系統Ｂ（図４）と分離されていることから、代数的手法に反復を利用したときの収束の問題と測定のためのノイズの側面とが本質的に分離されている。このため、本実施形態では、ノイズとは無関係な近似解に起因した現象が生じにくく、画像処理で一般に採用される各種のノイズ軽減処理が目的通りに効果を発揮しやすいという利点もある。

２−３．正則な正方システム行列を与える間引き順の決定
概要において上述したように、離散ラドン逆変換行列Ｗ_ＳＱ ^−１が正方システム行列Ｗ_ＳＱの逆行列として求まることは、間引き数列｛ｋ｝の適切さに依存している。その間引き数列｛ｋ｝は、行列デシメーション処理Ｓ１１４において決定されている。本実施形態の行列デシメーション処理Ｓ１１４は、結果として正則な正方システム行列を与えられれば限定はされず、その具体例として次の二つの手法が挙げられる。

２−３−１．階数での判定を利用する乱択アルゴリズム
図６Ａは、乱数またはランダムな選択を利用する乱択アルゴリズム（randomized algorithm）での行列デシメーション処理の一例を示すフローチャートである。オーバーサンプル・システム行列Ｗ_ＯＳ２１２の行の総数はＮ×（Ｎ＋ｎ）である（図５）。このうち、行列を正方行列にするために残す行数はＮ×Ｎであり、間引く行の総数はＮ×ｎである。行を特定するインデックス（行番号）の整数の範囲［１，Ｎ×（Ｎ＋ｎ）］のうちから、残すＮ×Ｎ個か間引くＮ×ｎ個かいずれかの数列を決定する。この決定のために図６Ａでは乱択手法を採用する。つまり、残すＮ×Ｎ個か間引くＮ×ｎ個かのいずれかを数値生成された乱数でランダムに決定すれば、最終的に間引き数列｛ｋ｝が決定できる（Ｓ１１４−２）。間引き数列｛ｋ｝が決まれば、対応した間引き処理によってオーバーサンプル・システム行列Ｗ_ＯＳ２１２から正方システム行列Ｗ_ＳＱを得ることができる（Ｓ１１４−４）。次に、求めた正方システム行列Ｗ_ＳＱの階数（rank）を算出する（Ｓ１１４−６）。算出した階数は、正方システム行列Ｗ_ＳＱが正則であればＮ×Ｎと同じ値であり、そうでなければＮ×Ｎよりも小さい値となる。このため、階数の値は正方システム行列Ｗ_ＳＱが逆行列を持ちうるかどうかの判定となる。その判定で階数がＮ×Ｎと一致した場合には（Ｓ１１４−８、Ｙの分枝）その時点での正方システム行列Ｗ_ＳＱを与える数列｛ｋ｝をそのまま間引き数列｛ｋ｝として間引き数列格納部３０に出力し、正方システム行列Ｗ_ＳＱ２１４を正方システム行列格納部６０に出力する（Ｓ１１４−１０）。これにより行列デシメーションの処理を終える。階数がＮ×Ｎと一致しなければ（Ｓ１１４−８、Ｎの分枝）、再びＳ１１４−２から別の乱数に基づいてランダムアルゴリズムによって間引き数列｛ｋ｝の選択から再度実行する（リトライ処理）。このような乱択アルゴリズムを採用すれば、実用的な速度で行列デシメーション処理Ｓ１１４を実行することができる。なお、後述するように、ランダムに選んだ場合に正則な正方システム行列Ｗ_ＳＱが得られる確率は解像度Ｎの値とオーバーサンプリング量ｎの値の選択の影響をうける。解像度Ｎが１０００程度の値の場合、ｎはその１〜２％程度である１０〜２０としても、十分に高い確率で正則な正方システム行列Ｗ_ＳＱを得ることができる。

２−３−２．グラム・シュミットの直交化法
図６Ｂは、グラム・シュミットの直交化を用いる行列デシメーション処理の一例を示すフローチャートである。グラム・シュミットの直交化法を採用する場合、まずオーバーサンプル・システム行列Ｗ_ＯＳの対角化処理を行う（Ｓ１１４−１２）。それに伴い間引き数列｛ｋ｝が決定される。この対角化は、例えばＱＲ分解法により行うことができる。また、対角化の処理により、対応する第１ベクトル２２８の各成分も変換（混合）する必要がある。このため、決定される間引き数列｛ｋ｝とともに、この変換を規定する行列データも決定する。対角化の後にグラム・シュミットの直交化法（Ｓ１１４−１４）を実行すれば、オーバーサンプル・システム行列Ｗ_ＯＳ２１２から線形独立なＮ×Ｎ個の行ベクトルの群を決定することができる。この行ベクトルの群を並べれば、正方システム行列Ｗ_ＳＱとなる。その際採用されなかった行のインデックスを示す間引き数列｛ｋ｝、正方システム行列Ｗ_ＳＱ、および対角化による成分を混合する変換行列を出力する（Ｓ１１４−１６）。なお、当該変換行列は、ベクトルデシメーション処理Ｓ１３０にて第１ベクトル２２８を変換するために利用される。

２−４．逆行列を持つ比率（解の存在比）の確認
次に、数値シミュレーションにより調査した正方システム行列Ｗ_ＳＱが逆行列をもつ成功率を説明する。表２は、オーバーサンプル・システム行列Ｗ_ＯＳ２１２を準備し、それぞれを図６Ａの乱択アルゴリズムに従う行列デシメーションの処理を２０回行うことにより正方システム行列Ｗ_ＳＱを求めてその階数を決定したときにＮ×Ｎとなった割合（解の存在比）をまとめたものである。解の存在比がＮおよびオーバーサンプル量ｎにどのように依存するかを調査するため、ここではＮ＝３２，３０，２８、オーバーサンプル量ｎ＝０〜１０の各組み合わせを網羅して調べた。

表２に示した解の存在比の値は、同じＮおよびｎを採用したときの現実のオーバーサンプル・システム行列Ｗ_ＯＳ２１２を対象にして図６Ａの乱択アルゴリズムを実行する場面において、正方システム行列Ｗ_ＳＱの階数がＮ×Ｎとなる成功率に対する指標となる。表２に示したように、概してオーバーサンプル量ｎを増やせば解の存在比は高くなり、実際にｎが３程度まででその値が急激に高まることがわかる。しかし個別にｎの値に対する傾向をみると、解の存在比はｎに対し単調増加はしてはいない。例えばいくつかの偶数（ｎ＝２、４、８、１０）では、解の存在比が０または極小さな値となっているものがある。この理由については現時点で定かではないが、Ｎが偶数の場合にｎも偶数であれば対称性が高まってしまうことに起因していると予測している。なお、ここでは実際に採用されるＮの場合に比べて、小規模なシステム行列を利用した確認にとどまるが、例えばＮ＝１２８とした実例でも解の存在比は十分に高く（後述）、さらにＮ＝１０００程度にした場合にも同様の傾向が実現するものと予測している。その際にどの程度のオーバーサンプリング量ｎが適切であるかまでは特定できていないが、表２のＮを変更することでオーバーサンプリング量ｎへの依存性が影響を受ける。このため、オーバーサンプリング量ｎだけではなく必要に応じて解像度Ｎを調整することも有用である。なお解像度Ｎは、例えば検出装置においていずれか一方または両方の端に位置する検出素子を意図的に利用せず、対応させて断層画像の周縁の画素を使用しないといったことにより、容易に減じることができる。

２−５．ブロックによるデシメーション
本実施形態にて正則な正方システム行列Ｗ_ＳＱが得られるデシメーション（間引き処理）は別の手法によって実施することもできる。本発明者は、オーバーサンプル・システム行列Ｗ_ＯＳを対象に、検出方向を単位としてまとまったブロックでの間引き処理によっても、正則な正方システム行列Ｗ_ＳＱが得られることを確認している。

図７は、解像度Ｎを８、オーバーサンプリング量ｎを１とした例において、ブロックによるデシメーション処理を示す説明図であり、図７Ａは、オーバーサンプル・システム行列Ｗ_ＯＳを画像により示しており、数値０を白、数値１を黒、それら中間値は、対応する明度のグレーにより表現している。この数値は、図３に示したように、ビーム幅τと画素の幅δと幾何学配置によって決まる。図７Ｂは、オーバーサンプル・システム行列Ｗ_ＯＳから、ブロックによるデシメーション処理によって得た正方システム行列Ｗ_ＳＱを同様に示したものである。

各ビームに対する各画素の位置における吸収の寄与を表すオーバーサンプル・システム行列Ｗ_ＯＳは、解像度Ｎが８、オーバーサンプリング量ｎが1の場合、７２行×６４列となる。これは、図３において、解像度Ｎは８×８個の画素を持つこと、および、検出装置２０が８個の検出素子２２を持つことに対応する。ここで、説明のため、θの原点を画素の並びの水平方向（紙面上の左右方向）または垂直方向（同、上下方向）のいずれか一方に揃っている方向に選ぶこととする。本節では一例として、紙面における１２時方向にビーム４が向かう検出方向をθの原点（＝０°）とする。この場合、図７Ａのオーバーサンプル・システム行列Ｗ_ＯＳの第１〜８行の部分行列は、検出方向θを０°とした場合の各検出素子２２へ各画素が与える寄与を数値化したものである。オーバーサンプル・システム行列Ｗ_ＯＳの第９〜１６行の部分行列は、検出方向θを２０°とした場合の各検出素子２２への各画素の寄与である。同様に、第１９〜２４行、第２５〜３２行、…、第６５〜７２行の部分行列は、順に、４０°、６０°、…、１６０°の検出方向とした場合の同様のものである。このように、図７Ａは、計９方向の各検出方向についての寄与を表すオーバーサンプル・システム行列Ｗ_ＯＳである。

オーバーサンプル・システム行列Ｗ_ＯＳの各要素について詳述する。オーバーサンプル・システム行列Ｗ_ＯＳの第１〜８行の部分行列は、検出装置２０の検出方向が原点（θ＝０°）にあるときに８個備わる検出素子２２それぞれに対する８×８個の画素の各画素の寄与である。またオーバーサンプル・システム行列Ｗ_ＯＳの第１〜６４列の要素は、８×８個の各画素の寄与を、画素の第１行の８画素について第１列〜第８列に、画素の第２行の８画素について画素の第２列〜第１６列に、…、画素の第８行の８画素について画素の第５７列〜第６４列に、と順に並べたものである。

オーバーサンプル・システム行列Ｗ_ＯＳの第９〜１６行の部分行列は、θが原点の次のサンプリングの検出方向（＝２０°）となる検出方向に向けられた場合についての各画素の寄与の数値である。この数値も、その検出方向で各画素が各検出素子２２に与える寄与である。第１７行の部分行列について以下同様である。一般の検出方向θの場合には、各画素の寄与が中間的な値を持ちうる。例えば、１画素の面積δ^２のうちどれだけの面積を幅τの各ビームが通過するか、という幾何学配置によって決まる。

本実施形態のブロックによるデシメーション処理は、図７の例に即して説明すれば、例えば、オーバーサンプル・システム行列Ｗ_ＯＳの第６５〜７２の部分行列を除去して、６４行６４列の正方行列とすることである。これは、間引き数列｛ｋ｝を６５から始まる８個の要素を持つ数列、つまり、６５、６６、…、７２とすることと同等である。より一般には、デシメーション処理でオーバーサンプル・システム行列Ｗ_ＯＳから除去するオーバーサンプル量ｎ×解像度Ｎ個の要素はまとまったブロックの連続する範囲とすることができる。さらにその連続する範囲は、１、Ｎ＋１、２×Ｎ＋１、…、Ｎ×Ｎ＋１、すなわち、ｓを１以上の整数として、（ｓ−１）×Ｎ＋１から始まるものとすることができる。デシメーション処理は、オーバーサンプル・システム行列Ｗ_ＯＳに残すＮ×Ｎ個の要素を決定することにより、残りのｎ×Ｎ個の要素を間接的に特定することによっても実行可能である。その際に残すＮ×Ｎ個の要素は、ｓを１以上の整数として、ｓ×Ｎで終る連続する範囲のものとする。この場合、オーバーサンプル量ｎが２以上の場合にオーバーサンプル・システム行列Ｗ_ＯＳから除去される行の範囲のブロックは、第１行を含む一つのブロックと第（Ｎ＋ｎ）×Ｎ行（最終行）を含むもう一つのブロックとに分かれて２ブロックとなることがある。いずれにしても、除去される行の範囲のブロックとなる部分行列は、連続した範囲であるだけでなく、検出方向を単位にしたものである。

デシメーション処理で連続する範囲を除去したり残したりする処理では、かならずしも正方システム行列Ｗ_ＳＱが正則になるとは限らない。正方システム行列Ｗ_ＳＱが正則になれば同処理は一度で十分であるが、正則とならないときには条件を変更して正則となるまで同処理は再度実行される。条件の変更は例えば整数ｓの値を変更する。

図８は、本実施形態におけるブロックによるデシメーション処理の成果を示す説明図であり、図８Ａは、図７Ｂの正方システム行列Ｗ_ＳＱから算出した逆行列である離散ラドン逆変換行列Ｗ_ＳＱ ^−１であり、図８Ｂは、離散ラドン逆変換行列Ｗ_ＳＱ ^−１を正方システム行列Ｗ_ＳＱに左から作用させて得られる行列を示す。図７Ｂの正方システム行列Ｗ_ＳＱが実際に、離散ラドン逆変換行列Ｗ_ＳＱ ^−１を逆行列としてもつことは、図８Ｂに示されている行列が単位行列であることに表れている。

本発明者は、解像度Ｎ＝８、オーバーサンプル量１の場合以外にも、Ｎ＝２およびＮ＝６４それぞれでオーバーサンプル量１の場合に、最後の検出方向に対応する連続する範囲のブロックによるデシメーション処理によって正則な正方システム行列Ｗ_ＳＱが得られることを確認した。

連続する範囲によるブロックでのデシメーション処理であっても、図４、５に基づく説明は同様に成立する。特に図５からも分かるように、間引き数列｛ｋ｝が連続する範囲を除去するものや、連続する範囲を残すものである場合、本実施形態では、除去される範囲または残す範囲に対応している検出方向が特定できる。例えば、図７では、オーバーサンプル・システム行列Ｗ_ＯＳの第６５〜７２行を除去しており、１６０°の検出方向に対応する要素が除去されている。これに応じて、検出動作で１６０°の検出方向において取得した測定値も間引き数列｛ｋ｝に応じて除去される。つまり、本実施形態においてブロックによるデシメーション処理を行う場合、間引き数列｛ｋ｝にて除去する範囲を特定するｎ×Ｎ個の要素に対応する検出方向について検出動作自体が省略可能である。結果として、検出動作では、オーバーサンプリングによる（Ｎ＋ｎ）方向の検出方向のうち、ｎ方向は省略できるため、最小限の検出方向は、Ｎ方向となり、その場合にはデシメーション処理が実質不要となる。

２−６．具体的構成のまとめ
以上に示した通り、本実施形態にて用いるオーバーサンプリングとデシメーション（間引き）を組み合わせる処理は、断層画像の撮像および取得のプロセスに高い実用性をもたらすものである。またその利点の一つとして、間引き数列｛ｋ｝を決定する処理はユーザーが断層画像を取得する処理とは別個に行うことができるため、逐次近似の処理、といったが計算量の多い処理は繰り返す必要もない。

３．本実装形態の実装
本実施形態は、断層撮像装置として実装することも、また既存の断層撮像装置のための、制御および処理ソフトウエアとして実装することもできる。

３−１．断層撮像装置での実装
本実施形態の断層撮像装置１００は、上述した処理を主に制御装置１５０において行うことにより、従来の断層撮像装置において、従来のＡＲＴを用いたものと同等またはそれ以上に良質な画像の断層画像を得ることができる。その際に必要となる計算資源も十分に実用的なものであり、また処理が高速である。

図９は、断層撮像装置によって本実施形態を実装する場合における制御装置１５０のいくつかの特有の機能部と特有のデータ格納部を示すブロック図である。システム行列生成部１６２、行列デシメーション処理部１６４、および逆行列算出部１６６は、それぞれ、システム行列生成ステップＳ１１２、行列デシメーション処理Ｓ１１４、および離散ラドン逆変換行列算出ステップＳ１１６を実行する。同様に、動作制御部１７４、離散ラドン変換処理部１７６、ベクトル化処理部１７８、ベクトルデシメーション処理部１８０、離散ラドン逆変換処理部１８２、および画像データ生成部１８４は、照射および検出ステップＳ１２４、離散ラドン変換Ｓ１２６、ベクトル化Ｓ１２８、ベクトルデシメーション処理Ｓ１３０、離散ラドン逆変換処理Ｓ１３２、および画像化ステップＳ１３４をそれぞれ実行する。その動作に当たり、間引き数列格納部３０、離散ラドン逆変換行列格納部４０、オーバーサンプル・システム行列格納部５０、および正方システム行列格納部６０も利用される。また、照射および検出Ｓ１２４からベクトル化Ｓ１２８までが等価な処理で実施できることと同様に、照射および検出ステップＳ１２４、離散ラドン変換Ｓ１２６、ベクトル化Ｓ１２８は、機能上これらを組み合わせた処理を行いうる機能部により実施することができる。さらに、制御装置１５０は、系統Ａの処理と系統Ｂの処理で必要となる機能部が明瞭に分かれていることから、必要に応じ系統Ａの処理を主として行うものであったり、逆に系統Ｂの処理を主として行うものであったり、様々な実装形態をとることができる。断層撮像装置の具体例としては、Ｘ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、光学ＣＴ装置が挙げられる。

３−２．コンピュータープログラムでの実装
また、断層撮像装置１００は、例えば従来のＡＲＴを実行するために用いられる照射装置１０検出装置２０を備えている断層撮像装置において、例えばコンピュータープログラムの実装形態として、断層撮像装置１００に備わる制御装置１５０に、図４に示したフローチャートの一部または全部の制御ステップや変換、算出等の計算処理を実行させることによって、その全体を断層撮像装置１００として実施することができる。本実施形態では、さらに、ユーザーが断層画像を撮像するタイミングでは、断層撮像装置１００に備わる制御装置１５０には系統Ｂ（図４）の処理のみを実行させることもできる。

４．計算機シミュレーションによる検証（実施例）
本実施形態による再構成手法の有効性を計算機シミュレーションにより確認した。逆ラドン変換を含む処理全般は、対比のための従来の手法のものも含めて、数式処理ソフトウエアMathematica（Wolfram Research, Inc., Champaign, Illinois）および画像処理ソフトウエアＩｍａｇｅＪ（アメリカ国立衛生研究所）を用いて実現した。本実施形態の制御装置１５０の処理も含め、シミュレーション処理のために、Ｉｎｔｅｌ社製Ｘｅｏｎ Ｅ５プロセッサー（３．５ＧＨｚ、６―Ｃｏｒｅタイプ）のＣＰＵと６４ＧＢの主記憶装置とを搭載したコンピューターを採用した。図１０は、検証のために用いた検証用データ（数値ファントム）を示す図（図１０Ａ、原画像という）および対応するオーバーサンプル・サイノグラム（図１０Ｂ）である。各図には背景と区別するための外形を縁取る実線を画像の周辺に追加し、画素の座標を表す画素番号の最初および最後を示す数値を上辺および左辺に示している。原画像（数値ファントム）は、計算機シミュレーションにおいて対象物のある断面における減弱率構造を模擬するために人為的に与える数値データであり、画素が空間位置に対応したＮ＝１２８つまり１２８×１２８の画像として表現できる。図１０Ａでは吸収の強い領域を暗く、弱い領域を明るく表示している。計算機シミュレーションでは、断層撮像装置１００による照射および検出Ｓ１２４の動作も模擬し、計算によって測定値に相当する数値データを得た。これを表すオーバーサンプル・サイノグラムでは、検出強度の弱い値を暗く、強い値を明るく表示している。図１０Ｂのオーバーサンプル・サイノグラムは、測定動作にあたる離散ラドン変換Ｓ１２６（図４）に基づき得られたものである。具体的には、Ｎ＝１２８として１２８個の検出素子２２を持つ検出装置２０によって、ｎ＝１を採用して１２９（＝Ｎ＋ｎ）方向からの投影分にてオーバーサンプリング離散ラドン変換を実行している。これを反映するため、図１０Ｂのオーバーサンプル・サイノグラムは１２８×１２９画素のサイズをもつ画像として表示している。図１０Ｂのオーバーサンプル・サイノグラムは本実施形態の動作のためのものであるが、従来のＦＢＰ法やＭＬ−ＥＭ法の処理のためにも利用している。従来の手法ではオーバーサンプリングすることは特段必要ではなくＮ方向（１２８方向）で取得したサイノグラム（図示しない）を採用することもできるが、ここでは、より直接的な比較とするために本実施形態および従来のもののために（Ｎ＋ｎ）方向からのオーバーサンプル・サイノグラムを採用している。

本実施形態の検証として、図５および図６に示した処理を実施して再構成画像を取得した。行列デシメーション処理Ｓ１１４の間引き数列｛ｋ｝を決定するために採用した処理は図６Ａに示した乱択アルゴリズムである。当該アルゴリズムに基づき間引き数列｛ｋ｝を決定する処理を実際に多数回行ったところ、階数の判定処理Ｓ１１４−８で階数がＮ×Ｎ（１２８×１２８）に一致せず（達せず）Ｎの分岐のとおり、乱数による決定Ｓ１１４−２に戻ってリトライ処理が必要となることがあった。解の存在比を反映して実際に必要となったリトライ処理の割合は、２０回のアルゴリズム動作に対し３〜４回程度にとどまり、残りの１６〜１７回は乱数による決定Ｓ１１４−２を一度実行すれば十分であった。

図１１に示す再構成画像は、本実施形態により求めたもの（図１１Ａ）、従来のＦＢＰ法においてＨａｎｎフィルターを用いたもの（図１１Ｂ）、および従来の代数再構成法として逐次近似計算（ＭＬ−ＥＭ法）を３０回行って得たもの（図１１Ｃ）である。本実施形態での再構成画像は１２８×１２８画素の範囲で得られるものの、これらの再構成画像は９１×９１画素の範囲のみ示している。これは、ＦＢＰ法でシミュレーションにより再構成像が得られる範囲が１２８画素を直径とする円形の領域となり、その円に概ね内接する正方形領域が比較可能となるためである。

図１０Ａと図１１Ａを比較すれば明らかなように、本実施形態においては再現性が極めて高い再構成画像が得られた。このような高い再現性は、リトライ処理の有無にかかわらず乱択アルゴリズムにより得られたどの間引き数列｛ｋ｝を採用した場合でも同様であった。これは、従来のＨａｎｎフィルターを用いたＦＢＰ法（図１１Ｂ）、およびＭＬ−ＥＭ法（図１１Ｃ）とは異なり、離散ラドン逆変換行列Ｗ_ＳＱ ^−１が代数的な厳密解となることに裏付けられている本実施形態の利点であり、それが画像自体の高い再現性を通じ確認されている。

次に、再構成像の再現性の比較を数値的指標に基づいて実施した。具体的には、各再構成画像の再現性について式（３）のＩＳＮＲ（Improvement in Signal-to-Noise ratio）を用いて評価した。

ここで、

は、順に、図１０Ａに示す原画像、対象とする図１１Ａ〜Ｃのいずれかの再構成画像、および単純な逆ラドン変換（フィルター無しの逆投影）の再構成画像について、それぞれをベクトル表示したものである。単純な逆ラドン変換の再構成画像は、比較のための基準であり、その画像を図１１Ｄに示す。式（３）の評価値は図１１に示した９１×９１画素の範囲を対象とした。単純な逆ラドン変換の再構成画像を基準に求めた図１１Ａ〜Ｃのいずれかの再構成画像のＩＮＳＲの値を表３に示す。

各ＩＳＮＲの値から、ＦＢＰ法や逐次近似再構成法を用いる従来例に比べ、本実施形態による再構成画像の再現性が大幅に高まっていることが数値的にも確認された。

５．画像再構成手法のバリエーション 上述した画像再構成手法では、オーバーサンプル・システム行列Ｗ_ＯＳから得た正方システム行列Ｗ_ＳＱが正則であって逆行列をもつことが保証された下で、正則な正方システム行列Ｗ_ＳＱから逆行列Ｗ_ＳＱ ^−１を実際に算出して画像再構成を実施した。本開示では、追加の画像再構成手法として、同様の正則な正方システム行列Ｗ_ＳＱの擬似逆行列を利用しても画像再構成を行いうる。

本追加の画像再構成手法のための擬似逆行列は、例えば特異値分解（singular value decomposition; SVD）の手法によって算出される。例えば、ムーア・ペンローズ(Moore-Penrose)型逆行列は正方行列に対して特異値分解により生成可能であり、本追加の画像再構成手法のための擬似逆行列として有用なものの一例である。

本開示の実施形態では、上述した画像再構成手法および本追加の画像再構成手法の双方で、式（１）に相当する式（４）が利用される。
Ｘ′＝ Ｗ_ＳＱ Ｘ （４）
この際、正方システム行列Ｗ_ＳＱは正則であることすなわち逆行列をもつことが保証されている。上述した画像再構成手法は、この正則な正方システム行列Ｗ_ＳＱの逆行列である離散ラドン逆変換行列Ｗ_ＳＱ ^−１を求めて式（５）に従って画像再構成を実施するものであった。
Ｘ＝ Ｗ_ＳＱ ^−１ Ｘ′ （５）

本追加の画像再構成手法では、式（４）までは同等であるものの、式（５）に代えて、
Ｘ＝ Ｗ_ｐ ^−１ Ｘ′ （６）
に従って画像再構成を実施する。ここで、Ｗ_ｐ ^−１は正方システム行列Ｗ_ＳＱの擬似逆行列である。つまり、正則であって逆行列の存在が保証されている正方システム行列Ｗ_ＳＱの逆行列を敢えて使わず、近似的な逆行列である擬似逆行列Ｗ_ｐ ^−１を利用しても、本追加の画像再構成手法では、画像再構成を行いうる。以下、Ｘ線ＣＴのデータによる再構成処理を行った事例（図１２、図１３）および中性子線ＣＴのデータによる再構成処理を行った事例（図１４、図１５）を説明する。

図１２は、マウス標本を利用して実際に取得したＸ線ＣＴのデータによる再構成処理を行った事例を説明する説明図であり、図１２Ａはマウス標本の写真、図１２Ｂは擬似逆行列Ｗ_ｐ ^−１を利用した再構成画像、図１２Ｃは比較のためのＦＢＰによる再構成画像、図１２Ｄは比較のためのＭＬ−ＥＭによる再構成画像、図１２Ｅは比較のための単純な逆ラドン変換による再構成画像である。図１２Ａの標本は、Ｘ線ＣＴ装置（SCANXMATE-E090S、コムスキャンテクノ株式会社、横浜）により撮影した。撮影データは、本来の性能とは異なり、比較計算のために適する少ない画素数（２５６×１３０）で取得し、再構成時に６４×６４、すなわちＮ＝６４になるように粗視化した。

本追加の画像再構成手法の図１２Ｂのための処理を、必要に応じて図４の要素を付記しながら説明すると、最初に、１８０°／６５すなわち２．７６９°ステップで６５方向の検出方向についてオーバーサンプル・システム行列Ｗ_ＯＳを算出した後（図４、Ｓ１１２）、最初の方向を０°として６５番目の１７７．２３０°の検出方向に対応する連続する範囲の要素をブロックで間引くことにより行列デシメーション処理をして、６４方向に対応する連続する範囲の要素を持つ正方システム行列Ｗ_ＳＱを決定した（Ｓ１１４）。これにより、間引き数列｛ｋ｝が決定される。この正方システム行列Ｗ_ＳＱが正則であることは、階数が４０９６（＝６４×６４）であることにより確認した。さらに、当該正方システム行列Ｗ_ＳＱのムーア・ペンローズ型逆行列を算出して擬似逆行列Ｗ_ｐ ^−１を算出して、それを離散ラドン逆変換行列Ｗ_ＳＱ ^−１とした。このため、図４において、離散ラドン逆変換行列算出ステップＳ１１６では、擬似逆行列が算出され、離散ラドン逆変換行列格納部４０（図４、図５）にはその擬似逆行列が格納される。

次に、マウス標本を対象に上記６５検出方向に対応するデータとして、３６０°／１３０すなわち２．７６９°ステップで１３０方向からデータを取得した。このとき、検出器の数は２５６である。ここで得られた２５６×１３０のサイノグラム画像に対して、１８０°分に対応する１から６５行目までを残すようにトリミングし、検出器方向の２５６ピクセルに対しては、２５％画像サイズ縮小を行っている（Ｓ１２２、Ｓ１２４）。本実施例は、このような処理を行ったデータに対しても、本手法が有効であることを示している。オーバーサンプル・サイノグラムに対応する第１ベクトルＸ_ＯＳ′を生成した（Ｓ１２６、Ｓ１２８）。次に、間引き数列｛ｋ｝に基づいて６５方向６５番目に当たる１７７．２３０°の検出方向の成分を間引くことにより、ベクトルデシメーション処理（Ｓ１３０）を実行して第２ベクトルＸ′を求めた。擬似逆行列Ｗ_ｐ ^−１を離散ラドン逆変換行列Ｗ_ＳＱ ^−１のために用いることにより、第３ベクトルＸを算出し（Ｓ１３２）、さらにデベクタライズ処理により再構成画像を得た（Ｓ１３４）。擬似逆行列Ｗ_ｐ ^−１の算出は、後述するトレランスを１／１９に設定して実施した。

図１２Ｃ〜図１２Ｅのためには、１８０°／６４すなわち２．８１２５°ステップで６４方向の検出方向から取得したデータを対象にして、そのサイノグラム画像を図１２Ｂ場合と同様のサイズとなるようリサイズしたデータを使用した。図１２Ｃ〜図１２ＥはImageJ(米国アメリカ国立衛生研究所)を利用し、図１２ＣはＦＢＰ処理、図１２ＤはＭＬ−ＥＭ処理（イタレーション数１）、図１２ＥはＢＰ処理により、それぞれの再構成画像を得た。なお、ＭＬ−ＥＭ処理でのイタレーション数と画質との関係は事前に調査し、１回が最も良好となることを確認している。

図１２Ｂ〜図１２Ｅから明らかなように、擬似逆行列Ｗ_ｐ ^−１を利用する追加の画像再構成手法ではＦＢＰ、ＭＬ−ＥＭ、ＢＰとの比較で良好な品質での再構成像を得ることができることを確認した。

擬似逆行列Ｗ_ｐ ^−１の算出は、典型的には適切なトレランス値を設定して行われる。トレランス値は、特異値分解の近似を通じて擬似逆行列を生成する処理において、特異値数列において特異値を除外するための基準値を、特異値の最大値との比率で示す。例えば、１／５をトレランス値とした場合には、特異値の最大値の１／５より小さい値を持つ特異値を０で置き換えて、特異値分解した行列から擬似逆行列を算出する。図１３は、トレランスを変更して算出した擬似逆行列Ｗ_ｐ ^−１により生成した再構成画像を含んでおり、図１３Ａ〜図１３Ｊは各図に付す値にトレランス値を設定したものである。また、図１３Ｋは、特異値を大きさ順に並べたグラフであり、縦軸に特異値の値、横軸に順序としてプロットしている。グラフには、複数のトレランス値を示す複数の直線を付記している。図１３から分かるように、トレランスを１／１とすると再構成画像は認識できる像を構成しない。トレランスを小さくするにつれて再構成画像は解像度を増してゆくが、１／１３以降は比較的安定している。

図１４は、コンクリート標本（直径２５ｍｍ×高さ６０ｍｍ、７０ｇ）を利用して実際に取得した中性子線ＣＴのデータによる再構成処理を行った事例を説明する説明図であり、図１４Ａはコンクリート標本の写真、図１４Ｂは擬似逆行列Ｗ_ｐ ^−１を利用した再構成画像、図１４Ｃは比較のためのＦＢＰによる再構成画像、図１４Ｄは比較のためのＭＬ−ＥＭによる再構成画像、図１４Ｅは比較のための単純な逆ラドン変換による再構成画像である。図１４Ａの標本は、小型中性子線源加速器によるＣＴ撮像装置（RIKEN Accelerator-driven Neutron Source、理研中性子ビーム技術開発チーム、理化学研究所、和光）により撮影した。撮影データは３６０°を３°おきに１２０ステップで取得し、再構成時に６４×６４、すなわちＮ＝６４になるように粗視化した。

本追加の画像再構成手法の図１４Ｂのための処理を、必要に応じて図４の要素を付記しながら説明すると、最初に、１８０°／６５すなわち２．７６９°ステップで６５方向の検出方向についてオーバーサンプル・システム行列Ｗ_ＯＳを算出した後（図４、Ｓ１１２）、最初の方向を０°として６５番目の１７７．２３０°の検出方向に対応する連続する範囲の要素をブロックで間引くことにより行列デシメーション処理をして、６４方向に対応する連続する範囲の要素を持つ正方システム行列Ｗ_ＳＱを決定した（Ｓ１１４）。これにより、間引き数列｛ｋ｝が決定される。この正方システム行列Ｗ_ＳＱが正則であることは、階数が４０９６（＝６４×６４）であることにより確認した。さらに、当該正方システム行列Ｗ_ＳＱのムーア・ペンローズ型逆行列を算出して擬似逆行列Ｗ_ｐ ^−１を算出して、それを離散ラドン逆変換行列Ｗ_ＳＱ ^−１とした。このため、図４において、離散ラドン逆変換行列算出ステップＳ１１６では、擬似逆行列が算出され、離散ラドン逆変換行列格納部４０（図４、図５）にはその擬似逆行列が格納される。

次に、コンクリート標本を対象に上記６５検出方向に対応するデータとして、３６０°／１２０すなわち３°ステップで１２０方向からデータを取得した。ここで得られた１２０方向分のサイノグラム画像に対して、１８０°分に対応する１から６０行目までを残すようにトリミングし、画像をリサイズして６５方向分に伸長した（Ｓ１２２、Ｓ１２４）。本実施例は、このような処理を行ったデータに対しても、本手法が有効であることを示している。オーバーサンプル・サイノグラムに対応する第１ベクトルＸ_ＯＳ′を生成した（Ｓ１２６、Ｓ１２８）。次に、間引き数列｛ｋ｝に基づいて６５方向６５番目に当たる１７７．２３０°の検出方向の成分を間引くことにより、ベクトルデシメーション処理（Ｓ１３０）を実行して第２ベクトルＸ′を求めた。擬似逆行列Ｗ_ｐ ^−１を離散ラドン逆変換行列Ｗ_ＳＱ ^−１のために用いることにより、第３ベクトルＸを算出し（Ｓ１３２）、さらにデベクタライズ処理により再構成画像を得た（Ｓ１３４）。擬似逆行列Ｗ_ｐ ^−１の算出は、後述するトレランスを１／１９に設定して実施した。

図１４Ｃ〜図１４Ｅのためには、３６０°／１２０すなわち３°ステップで１２０方向の検出方向から取得したデータを対象にして、そのサイノグラム画像を図１４Ｂと同様のサイズとなるようリサイズしたデータを使用した。図１４Ｃ〜図１４ＥはImageJを利用し、図１４ＣはＦＢＰ処理、図１４ＤはＭＬ−ＥＭ処理（イタレーション数１）、図１４ＥはＢＰ処理により、それぞれの再構成画像を得た。

図１４Ｂ〜図１４Ｅから明らかなように、擬似逆行列Ｗ_ｐ ^−１を利用する追加の画像再構成手法ではＦＢＰ、ＭＬ−ＥＭ、ＢＰとの比較で良好な品質での再構成像を得ることができることを確認した。

また、図１５は、トレランスを変更して算出した擬似逆行列Ｗ_ｐ ^−１により生成した再構成画像を含んでおり、図１５Ａ〜図１５Ｊは各図に付す値にトレランス値を設定したものである。図１５から分かるように、中性子線ＣＴの場合においても、トレランスを１／１とすると再構成画像は認識できる像を構成しない。また、トレランスを小さくするにつれて再構成画像は解像度を増してゆくが、１／１３以降は比較的安定している。

図１６は、本追加の画像再構成手法における擬似逆行列Ｗ_ｐ ^−１の例であり、図１６Ａは、図７Ｂの正方システム行列Ｗ_ＳＱに対して求めた擬似逆行列Ｗ_ｐ ^−１、図１６Ｂは、擬似逆行列Ｗ_ｐ ^−１を正方システム行列Ｗ_ＳＱに左から作用させて得られる行列を示す。図８Ａと図１６Ａとの比較から、同一の正方システム行列Ｗ_ＳＱを対象にして算出されていても、擬似逆行列Ｗ_ｐ ^−１は離散ラドン逆変換行列Ｗ_ＳＱ ^−１とは一致しない。その場合であっても、図８Ｂと図１６Ｂとの比較から、逆行列Ｗ_ｐ ^−１を正方システム行列Ｗ_ＳＱに左から作用させた結果のみを見ると、単位行列に近似できる行列が算出される。これは、逆行列Ｗ_ｐ ^−１を正方システム行列Ｗ_ＳＱに右から作用させても同様である（図示しない）。

以上のように、正則で逆行列を持つことが保証されている正方システム行列Ｗ_ＳＱを対象に擬似逆行列を利用することによって、十分な品質の再構成画像が生成されることを確認した。

なお、擬似逆行列のためのトレランス値は、一般に、小さい値を設定するほど算出される擬似逆行列が逆行列に近づき、再構成画像の解像度も高くなる。しかし、現実の測定装置を用いた測定においては、小さければ良好な再構成画像が得られるわけではない。幾何学的要因に帰着できない種々の攪乱要因が存在しているためである。例えば、検出信号に含まれるノイズや、ビームの不安定性、ビームの散乱、装置自体が持っている画像特性に起因した要因が、現状の測定装置に攪乱要因として残存している。通常これらの要因をすべて反映した装置全体の撮像性能が、例えば解像特性を示すＭＴＦ（Modulation Transfer Function）や、画像の粒状度の指標となるウィーナー・スペクトルによって評価されている。本追加の画像再構成手法において０より大きな所定のトレランス値（例えば、図１３、図１５では１／１３〜１／１９）を用いることにより、幾何学的要因以外の攪乱要因を除外するフィルター動作が実現する。このフィルター動作は、トレランス値を増減することによりその効果を容易に調整できることから、解像度を重視するか、定量性を重視するかなどの着目する画像特性に応じて必要な再構成像を作成することに役立ちうる。

６．変形例
本実施形態は、断層画像データ取得方法、断層画像データの取得装置やそのためのコンピュータープログラムのいずれで実施する場合であっても、上述したもの以外の種々の形態で実施することができる。例えば、図４、５、７において系統Ａにおいて機器情報を決定する処理のためにユーザーが関与することは必ずしも必要なく、ユーザーを系統Ｂの測定情報を取得することのみに関与させるような実施態様への変形も容易である。さらに、例えば本実施形態を実施するためのコンピュータープログラムで提供される場合には、照射および検出を行うハードウエアは既存のものや既に設置済のものを採用し、制御用コンピュータープログラムに上述した処理を新たに組み込むことで本実施形態を実装することもできる。このように、本実施形態の手法は、種々の実施態様を通じて実施することができる。

また、実施形態の変形に加え、その処理の細部を変形したものも本実施形態の一部をなす。その一例が、図３を参照して説明した検出素子２２に対するビーム４についてのその画素のもつ寄与の決定手法である。図３で模様を付した部分の面積を寄与とする手法に加え、他の実用に付されている手法は本実施形態においても採用できる。例えば、面積に代え、掃く長さに応じた重みを持たせたり、その掃く長さを効率よく評価するシドンのアルゴリズム（Siddon's algorithm）を採用して上記寄与を決定することも有利である（非特許文献３）。さらに、さらなる高画質化を行うために、アンチエイリアス処理できるSunnegardh and Danielsonの方法を利用することも有利である（非特許文献４）。

さらに、図３以降の説明は、一例の断層撮像装置１００の動作により説明したが、本実施形態は、必ずしも照射装置１０を採用するものには限定されない。ＳＰＥＣＴ（単一光子放射断層撮影）など、照射装置を用いることなく対象物の断層撮影を行う他の手法においても、本実施形態を同様に適用することができる。

以上、本開示の実施形態を具体的に説明した。上述の実施形態、変形例および実施例は、本出願において開示される発明を説明するために記載されたものであり、本出願の発明の範囲は、請求の範囲の記載に基づき定められるべきものである。実施形態の他の組合せを含む本開示の範囲内に存在する変形例もまた、請求の範囲に含まれるものである。

本開示は、波または粒子が例えば対象物に照射されたりすることによって、対象物自体を波または粒子が透過するような任意の断層撮像手法のために利用可能であり、例えば、Ｘ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、光学ＣＴ装置が挙げられる。

１００ 断層撮像装置（断層画像データの取得装置）
２ 対象物
４ ビーム
１０ 照射装置
２０ 検出装置
２２ 検出素子
３０ 間引き数列格納部
４０ 離散ラドン逆変換行列格納部
５０ オーバーサンプル・システム行列格納部
６０ 正方システム行列格納部
１５０ 制御装置
１６２ システム行列生成部
１６４ 行列デシメーション処理部
１６６ 逆行列算出部
１７４ 動作制御部
１７６ 離散ラドン変換処理部
１７８ ベクトル化処理部
１８０ ベクトルデシメーション処理部
１８２ 離散ラドン逆変換処理部
１８４ 画像データ生成部
２１２ オーバーサンプル・システム行列Ｗ_ＯＳ
２１４ 正方システム行列Ｗ_ＳＱ
２２６ オーバーサンプル・サイノグラム
２２８ 第１ベクトル
２３０ 第２ベクトル
２３２ 第３ベクトル

Claims (17)

1. 少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置の検出範囲に対象物を配置するステップと、
   照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない（Ｎ＋ｎ）方向（ｎは１以上の整数）にわたって前記検出方向を変更して繰り返す検出ステップと、
   各行および各列を各検出方向および各検出素子に対応させた（Ｎ＋ｎ）行Ｎ列のオーバーサンプル・サイノグラムをベクトル化したものに相当するＮ×（Ｎ＋ｎ）個の要素をもつ第１ベクトルを、前記検出動作における前記検出装置による検出信号から得るベクトル化ステップと、
   前記第１ベクトルから、間引き順に当たる計ｎ×Ｎ個の要素を除去して、残りのＮ×Ｎ個の要素をもつ第２ベクトルを得るベクトルデシメーションステップと、
   前記第２ベクトルに離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第３ベクトルを得る離散ラドン逆変換ステップと、
   前記第３ベクトルをデベクトル化することにより、前記対象物に対して静止した画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップと
   を含む断層画像データ取得方法。
2. 前記ベクトルデシメーションステップより前に、
   （Ｎ＋ｎ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×（Ｎ＋ｎ）個の要素の列または行ベクトルを、前記二次元断層画像の画素に対応させて列または行方向に並べたオーバーサンプル・システム行列を得るステップと、
   前記オーバーサンプル・システム行列から、前記間引き順に当たる計ｎ×Ｎ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を得る行列デシメーションステップと
   をさらに含む
   請求項１に記載の断層画像データ取得方法。
3. 前記行列デシメーションステップより後でかつ前記離散ラドン逆変換ステップより前に、
   前記正方システム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得る逆行列算出ステップ
   をさらに含む
   請求項２に記載の断層画像データ取得方法。
4. 前記行列デシメーションステップより後でかつ前記離散ラドン逆変換ステップより前に、
   前記正方システム行列の擬似逆行列を算出することにより前記離散ラドン逆変換行列を得る擬似逆行列算出ステップ
   をさらに含む
   請求項２に記載の断層画像データ取得方法。
5. 前記行列デシメーションステップが、
   １以上Ｎ×（Ｎ＋ｎ）以下の整数から前記間引き順を与える数列を選択するステップであって、該数列が、１以上Ｎ×（Ｎ＋ｎ）以下の整数のなす数列の要素から、選択された前記整数列のＮ×Ｎ個の要素を除いて得たｎ×Ｎ個の要素をもつ整数列であるステップ
   を含むものである、
   請求項２に記載の断層画像データ取得方法。
6. 前記行列デシメーションステップが、
   １以上Ｎ×（Ｎ＋ｎ）以下の整数から前記間引き順を与える数列を選択するステップであって、該数列が、１以上Ｎ×（Ｎ＋ｎ）以下の整数のなす数列の要素から、ｓ×Ｎ（ｓは１以上の整数）で終る連続する範囲のＮ×Ｎ個の要素を除いて得たｎ×Ｎ個の要素をもつ整数列であるステップ
   を含むものである、
   請求項２に記載の断層画像データ取得方法。
7. 前記行列デシメーションステップが、
   前記間引き順に当たる位置の要素を除くことにより得た前記正方システム行列を対象に階数を決定するステップと、
   前記階数がＮ×Ｎに一致するかを判定するステップと
   をさらに含んでいる、
   請求項５または請求項６に記載の断層画像データ取得方法。
8. 前記行列デシメーションステップが、
   １以上Ｎ×（Ｎ＋ｎ）以下の整数から前記間引き順を与える数列を選択するステップであって、該数列がランダムに重複なく選択されたｎ×Ｎ個の要素をもつ整数列であるステップ
   を含むものである、
   請求項２に記載の断層画像データ取得方法。
9. 前記行列デシメーションステップが、
   １以上Ｎ×（Ｎ＋ｎ）以下の整数から前記間引き順を与える数列を選択するステップであって、該数列が（ｓ−１）×Ｎ＋１（ｓは１以上の整数）から始まり連続する範囲のｎ×Ｎ個の要素をもつ整数列であるステップ
   を含むものである、
   請求項２に記載の断層画像データ取得方法。
10. 前記行列デシメーションステップが、
    前記間引き順に当たる位置の要素を除くことにより得た前記正方システム行列を対象に階数を決定するステップと、
    前記階数がＮ×Ｎに一致するかを判定するステップと
    をさらに含んでいる、
    請求項８または請求項９に記載の断層画像データ取得方法。
11. 前記間引き順を与える数列の要素に対応する検出方向の前記検出ステップでの検出動作が省略されている、
    請求項６または請求項９に記載の断層画像データ取得方法。
12. 前記行列デシメーションステップが、
    前記オーバーサンプル・システム行列の要素の行ベクトルまたは列ベクトルの集合から、グラム・シュミットの直交化法により、線形独立なＮ×Ｎ個の行ベクトルまたは列ベクトルの群を生成し、前記行ベクトルまたは列ベクトルの群を列方向または行方向に並べることにより前記正方システム行列を得るステップと、
    前記オーバーサンプル・システム行列の要素の行ベクトルのうち、前記線形独立な行ベクトルの群の生成に用いられなかったものの行番号または列番号を前記間引き順のために選択するステップと
    を含むものである、
    請求項２に記載の断層画像データ取得方法。
13. 少なくとも一列に並んだＮ個（Ｎは１以上の整数）の検出素子をもつ検出装置と、
    前記検出装置の検出対象となる波または粒子の照射装置により前記検出装置に向けて前記波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記検出装置の検出範囲に配置した対象物との間の検出方向を変更し、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより検出するよう、前記検出装置を制御する制御装置と、
    を備える断層画像データの取得装置であって
    間引き順を与える数列データを格納する間引き数列格納部と、
    離散ラドン逆変換行列を格納する離散ラドン逆変換行列格納部と
    をさらに備え、
    前記制御装置は、
    前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない（Ｎ＋ｎ）方向（ｎは１以上の整数）にわたって前記検出方向を変更し、前記検出装置に、各行および各列を各検出方向および各検出素子に対応させた（Ｎ＋ｎ）行Ｎ列のオーバーサンプル・サイノグラムのための検出信号を出力させる離散ラドン変換処理部と、
    前記検出信号から、前記オーバーサンプル・サイノグラムをベクトル化したものに相当するＮ×（Ｎ＋ｎ）個の要素をもつ第１ベクトルを得るベクトル化処理部と、
    前記第１ベクトルから、前記間引き順を与える数列データに応じ計ｎ×Ｎ個の要素を除去して、残りのＮ×Ｎ個の要素をもつ第２ベクトルを得るベクトルデシメーション処理部と、
    前記第２ベクトルに前記離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第３ベクトルを得る離散ラドン逆変換処理部と、
    前記第３ベクトルをデベクトル化することにより、前記対象物に対して静止した画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成部と
    を備えている断層画像データ取得装置。
14. オーバーサンプル・システム行列格納部と、
    正方システム行列格納部と
    をさらに備え、
    前記制御装置は、
    （Ｎ＋ｎ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×（Ｎ＋ｎ）個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて列方向または行方向に並べたオーバーサンプル・システム行列を得て前記オーバーサンプル・システム行列格納部に格納するオーバーサンプル・システム行列生成部と、
    前記オーバーサンプル・システム行列から、前記間引き順に当たる計ｎ×Ｎ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を得て前記正方システム行列格納部に格納する行列デシメーション処理部と、
    前記正方システム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得て前記離散ラドン逆変換行列格納部に格納する逆行列算出部と
    をさらに備えるものである
    請求項１３に記載の断層画像データ取得装置。
15. オーバーサンプル・システム行列格納部と、
    正方システム行列格納部と
    をさらに備え、
    前記制御装置は、
    （Ｎ＋ｎ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×（Ｎ＋ｎ）個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて列方向または行方向に並べたオーバーサンプル・システム行列を得て前記オーバーサンプル・システム行列格納部に格納するオーバーサンプル・システム行列生成部と、
    前記オーバーサンプル・システム行列から、前記間引き順に当たる計ｎ×Ｎ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を得て前記正方システム行列格納部に格納する行列デシメーション処理部と、
    前記正方システム行列の擬似逆行列を算出することにより前記離散ラドン逆変換行列を得て前記離散ラドン逆変換行列格納部に格納する擬似逆行列算出部と
    をさらに備えるものである
    請求項１３に記載の断層画像データ取得装置。
16. 少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置と、
    前記検出装置の検出対象となる波または粒子の照射装置により前記検出装置に向けて前記波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記波または粒子の検出装置の検出範囲に配置した対象物に対する前記波または粒子の相対的な検出方向を変更し、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより検出するよう、前記検出装置を制御する制御装置と
    を備える断層画像データの取得装置のための制御プログラムであって
    前記制御装置に、
    前記対象物からみた前記波または粒子の相対的な検出方向における互いに重複しない（Ｎ＋ｎ）方向（ｎは１以上の整数）にわたって前記検出方向を変更し、前記検出装置に、各行および各列を各検出方向および各検出素子に対応させた（Ｎ＋ｎ）行Ｎ列のオーバーサンプル・サイノグラムのための検出信号を前記検出装置に取得させる離散ラドン変換ステップと、
    前記検出信号から、前記オーバーサンプル・サイノグラムをベクトル化したものに相当するＮ×（Ｎ＋ｎ）個の要素をもつ第１ベクトルを得るベクトル化ステップと、
    間引き数列格納部から間引き順を与える数列データを読み出して、前記第１ベクトルから、前記数列データに応じ計ｎ×Ｎ個の要素を除去して、残りのＮ×Ｎ個の要素をもつ第２ベクトルを得るベクトルデシメーションステップと、
    離散ラドン逆変換行列格納部から離散ラドン逆変換行列を読み出して、前記第２ベクトルに前記離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第３ベクトルを得る離散ラドン逆変換ステップと、
    前記第３ベクトルをデベクトル化することにより、前記対象物に対して静止した画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップと
    を実行させる制御プログラム。
17. 前記制御装置に、
    （Ｎ＋ｎ）方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×（Ｎ＋ｎ）個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて列方向または行方向に並べたオーバーサンプル・システム行列を得てオーバーサンプル・システム行列格納部に格納するステップと、
    前記オーバーサンプル・システム行列から、前記間引き順に当たる計ｎ×Ｎ個の行ベクトルまたは列ベクトルの要素を除去することにより正方システム行列を得て正方システム行列格納部に格納するステップと、
    前記正方システム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得て離散ラドン逆変換行列格納部に格納する逆行列算出ステップと
    をさらに実行させる請求項１６に記載の制御プログラム。

Images