(57)【要約】  計算負荷を増大させずに再構築される断層画像の再現性を高めるために、本開示の実施形態では、角度スキャン範囲を画定する二つの境界に近接する二つの方向のいずれかが、Ｎ画素×Ｎ画素の二次元断層画像の座標軸のいずれかからオフセットされており、そのオフセット角を０度を超え９０度未満または−９０度を超え０度未満とする。Ｎ個の検出素子をもつ検出装置で各検出方向にて検出し、Ｎ×Ｎ個の要素をもつ第１ベクトルを検出動作における検出装置による検出信号から得る。第１ベクトルに離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第２ベクトルを得る。第２ベクトルをデベクトル化してＮ画素×Ｎ画素の二次元断層画像のための画像データを得る。オフセットのためにシステム行列の逆行列が求まり、離散ラドン逆変換行列として用いる。
2019230741.tif 000002
【特許請求の範囲】
【請求項１】
少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置の検出範囲に対象物を配置するステップと、
照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、前記対象物からみた前記波または粒子の相対的な各検出方向で行う検出ステップと、
各行および各列を各検出方向および各検出素子に対応させたＮ行Ｎ列のサイノグラムをベクトル化したものに相当するＮ×Ｎ個の要素をもつ第１ベクトルを、前記検出動作における前記検出装置による検出信号から得るベクトル化ステップと、
前記第１ベクトルに離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第２ベクトルを得る離散ラドン逆変換ステップと、
前記第２ベクトルをデベクトル化することにより、前記対象物に対して定められる二つの座標軸を持つ二次元座標で各画素が指定され、両座標軸に共通のピッチで縦横に並んだ画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップと
を含み、
前記検出ステップの前記各検出方向は、前記二次元断層画像を再構成するための角度スキャン範囲にわたり、該角度スキャン範囲を均等にＮ分割した各角度刻みに対応する角度サンプリング方向であり、前記角度スキャン範囲を規定する二つの境界に近接する二つの方向のいずれかが、前記座標軸のいずれかからオフセットされており、そのオフセット角が０度を超え９０度未満または−９０度を超え０度未満である、断層画像データ取得方法。
【請求項２】
前記各検出方向は、前記座標軸のいずれとも平行でも垂直でもない、
請求項１に記載の断層画像データ取得方法。
【請求項３】
Ｎ方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×Ｎ個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて行方向または列方向に並べたシステム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得る逆行列算出ステップ
をさらに含む
請求項２に記載の断層画像データ取得方法。
【請求項４】
Ｎ方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×Ｎ個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて行方向または列方向に並べたシステム行列の擬似逆行列を算出することにより前記離散ラドン逆変換行列を得る擬似逆行列算出ステップ
をさらに含む
請求項２に記載の断層画像データ取得方法。
【請求項５】
前記システム行列が正則行列であることを判定する正則判定ステップ
をさらに含む、
請求項３または請求項４に記載の断層画像データ取得方法。
【請求項６】
前記正則判定ステップと、
前記システム行列が正則行列でない場合に前記オフセット角を変更するオフセット角再設定ステップと
を前記システム行列が正則行列となるまで繰り返すステップ
を含む、
請求項５に記載の断層画像データ取得方法。
【請求項７】
前記正則判定ステップが、前記システム行列を対象に階数を決定し、該階数がＮ×Ｎに一致するかどうかを判定するステップである、
請求項５に記載の断層画像データ取得方法。
【請求項８】
少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置と、
前記検出装置の検出対象となる波または粒子の照射装置により前記検出装置に向けて前記波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記検出装置の検出範囲に配置した対象物との間の検出方向を変更し、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより検出するよう、前記検出装置を制御する制御装置と
を備える断層画像データの取得装置であって
離散ラドン逆変換行列を格納する離散ラドン逆変換行列格納部と
をさらに備え、
前記制御装置は、
前記対象物からみた前記波または粒子の相対的な検出方向を変更し、前記検出装置に、各行および各列を各検出方向および各検出素子に対応させたＮ行Ｎ列のサイノグラムのための検出信号を前記検出装置に取得させる離散ラドン変換処理部と、
前記検出信号から、前記サイノグラムをベクトル化したものに相当するＮ×Ｎ個の要素をもつ第１ベクトルを得るベクトル化処理部と、
前記第１ベクトルに前記離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第２ベクトルを得る離散ラドン逆変換処理部と、
前記第２ベクトルをデベクトル化することにより、前記対象物に対して定められる二つの座標軸を持つ二次元座標で各画素が指定され、両座標軸に共通のピッチで縦横に並んだ画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成部と
を備えており、
前記制御装置は、前記二次元断層画像を再構成するための角度スキャン範囲にわたり、該角度スキャン範囲を均等にＮ分割した各角度刻みに対応する角度サンプリング方向である前記各検出方向で検出するよう前記検出装置を制御するものであり、
前記角度スキャン範囲を規定する二つの境界に近接する二つの方向のいずれかが、前記座標軸のいずれかからオフセットされており、そのオフセット角が０度を超え９０度未満または−９０度を超え０度未満である、断層画像データの取得装置。
【請求項９】
前記各検出方向は、前記座標軸のいずれとも平行でも垂直でもない、
請求項８に記載の断層画像データの取得装置。
【請求項１０】
システム行列格納部
をさらに備え、
前記制御装置は、
Ｎ方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×Ｎ個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて行方向または列方向に並べたシステム行列を得て前記システム行列格納部に格納するシステム行列生成部と、
前記システム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得て前記離散ラドン逆変換行列格納部に格納する逆行列算出部と
をさらに備えるものである
請求項８に記載の断層画像データ取得装置。
【請求項１１】
少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置と、
前記検出装置の検出対象となる波または粒子の照射装置により前記検出装置に向けて前記波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記波または粒子の検出装置の検出範囲に配置した対象物に対する前記波または粒子の相対的な検出方向を変更し、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより検出するよう、前記検出装置を制御する制御装置と
を備える断層画像データの取得装置のための制御プログラムであって
前記制御装置に、
前記対象物からみた前記波または粒子の相対的な検出方向を変更し、前記検出装置に、各行および各列を各検出方向および各検出素子に対応させたＮ行Ｎ列のサイノグラムのための検出信号を前記検出装置に取得させる離散ラドン変換ステップと、
前記検出信号から、前記サイノグラムをベクトル化したものに相当するＮ×Ｎ個の要素をもつ第１ベクトルを得るベクトル化ステップと、
離散ラドン逆変換行列格納部から離散ラドン逆変換行列を読み出して、前記第１ベクトルに前記離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第２ベクトルを得る離散ラドン逆変換ステップと、
前記第２ベクトルをデベクトル化することにより、前記対象物に対して定められる二つの座標軸を持つ二次元座標で各画素が指定され、両座標軸に共通のピッチで縦横に並んだ画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップと
を実行させ、
前記検出ステップの前記角検出方向は、前記二次元断層画像を再構成するための角度スキャン範囲にわたり、該角度スキャン範囲を均等にＮ分割した各角度刻みに対応する角度サンプリング方向であり、
前記角度スキャン範囲を規定する二つの境界に近接する二つの方向のいずれかが、前記座標軸のいずれかからオフセットされており、そのオフセット角が０度を超え９０度未満または−９０度を超え０度未満である、断層画像データの取得装置のための制御プログラム。
【請求項１２】
前記各検出方向は、前記座標軸のいずれとも平行でも垂直でもない、
請求項１１に記載の制御プログラム。
【請求項１３】
前記制御装置に、
Ｎ方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×Ｎ個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて行方向または列方向に並べたシステム行列の逆行列を算出することにより前記離散ラドン逆変換行列を得て離散ラドン逆変換行列格納部に格納する逆行列算出ステップ
をさらに実行させる、
請求項１２に記載の制御プログラム。
【請求項１４】
前記制御装置に、
Ｎ方向の前記検出方向それぞれに対応させて、前記波または粒子のうちＮ個の前記検出素子のそれぞれに向かう部分の行路に対する前記二次元断層画像の各画素の寄与を示す重み値を要素にもつＮ×Ｎ個の要素の列ベクトルまたは行ベクトルを、前記二次元断層画像の画素に対応させて行方向または列方向に並べたシステム行列の擬似逆行列を算出することにより前記離散ラドン逆変換行列を得て離散ラドン逆変換行列格納部に格納する擬似逆行列算出ステップ
をさらに実行させる、
請求項１２に記載の制御プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本開示は角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラムに関する。さらに詳細には本開示は、角度オフセットに基づく離散ラドン逆変換の効率の良い代数的手法を組み込んだ断層画像データの取得方法、取得装置、および制御プログラムに関する。
【背景技術】
【０００２】
Ｘ線、ガンマ線、光波などの透過性の電磁波や粒子線、または地震波やその他の任意の透過性を示す波または粒子を利用する断層撮像法（トモグラフィー）が実用化されている。その代表例が、医療、工業用途で用いられるＸ線ＣＴスキャナー、ＳＰＥＣＴ（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)
【発明の概要】
【発明が解決しようとする課題】
【０００７】
画像の精細度の向上のための画素数増大への要請に加え、再構成像の品質に対する要請も弱まる気配がない。撮像画素数が増大した状況で従来のＡＲＴを採用しても、解くべき連立一次方程式の元（決定すべき変数）の数が増大してしまい膨大な計算資源が必要となり、その処理に長い時間を要す。このため従来のＡＲＴでの画像データ取得処理に実用性を求めるには、現状では計算能力の向上を待つしかない。特に、逐次近似に伴う反復計算を撮像の毎に再実行するような従来のＡＲＴの処理に頼る限り、例えば断層画像を多数得て三次元のボリューム像を得ることは現実的ではない。【０００８】
本開示は上記問題の少なくともいくつかを解決することを課題とする。本開示は、像の再構成を担う離散ラドン逆変換のために代数的厳密解の裏付けをもつような新たなサンプリング手法を採用することにより、高品質な断層撮像法に高い実用性をもたらすものである。
【課題を解決するための手段】
【０００９】
本発明者は、従来のＡＲＴの問題点がサンプリング（標本化）といった離散化手法に起因していることに気づいた。従来のＡＲＴでは、多元の連立一次方程式に重複する式が多数含まれてしまい、元（決定すべき変数）の数に対し利用できる方程式の数が実質的に足りなくなっている。このために逐次近似が必要となってプロセッサーの計算能力による解決が指向されているのである。これに対し、本発明者は、離散化手法から再検討することにより、従来にないサンプリングと適切な処理とを組み合わせる上記課題への新たな解決手法を完成させた。【００１０】
すなわち、本開示のある態様においては、少なくとも一列に並んだＮ個（Ｎは正の整数）の検出素子をもつ検出装置の検出範囲に対象物を配置するステップと、照射装置により前記検出装置に向けて前記検出装置の検出対象となる波または粒子を照射しながら、または照射装置によらずに生じた前記波または粒子を対象物の各部に透過させながら、前記対象物の各部を透過した後の前記波または粒子を前記検出素子それぞれにより受けて前記検出素子別の強度値を得る検出動作を、前記対象物からみた前記波または粒子の相対的な各検出方向で行う検出ステップと、各行および各列を各検出方向および各検出素子に対応させたＮ行Ｎ列のサイノグラムをベクトル化したものに相当するＮ×Ｎ個の要素をもつ第１ベクトルを、前記検出動作における前記検出装置による検出信号から得るベクトル化ステップと、前記第１ベクトルに離散ラドン逆変換行列を作用させて、Ｎ×Ｎ個の要素をもつ第２ベクトルを得る離散ラドン逆変換ステップと、前記第２ベクトルをデベクトル化することにより、前記対象物に対して定められる二つの座標軸を持つ二次元座標で各画素が指定され、両座標軸に共通のピッチで縦横に並んだ画素配置をもつＮ画素×Ｎ画素の二次元断層画像のための画像データを得る画像データ生成ステップとを含み、前記検出ステップの前記各検出方向は、前記二次元断層画像を再構成するための角度スキャン範囲にわたり、該角度スキャン範囲を均等にＮ分割した各角度刻みに対応する角度サンプリング方向であり、前記角度スキャン範囲を規定する二つの境界に近接する二つの方向のいずれかが、前記座標軸のいずれかからオフセットされており、そのオフセット角が０度を超え９０度未満または−９０度を超え０度未満である、断層画像データ取得方法が提供される。本開示のかかる態様においては、検出方向を走査する範囲（角度スキャン範囲）内を均等にＮ分割した各検出方向で検出が行われる。その各検出方向は、二次元断層画像の画素位置を縦横で特定する二つの座標軸のいずれとも平行でも垂直でもない。このような画像取得方法を採用することにより、代数的手法による厳密な解の算出が可能となって、近似を含まない離散ラドン逆変換が可能となる。【００１１】
また、本開示は、断層画像データ取得装置の態様としても実施することができる。同様に、本開示は、断層画像データ取得装置のための制御プログラムの態様でも実施することができる。【００１２】
さらに本開示のいずれかの態様においては、角度スキャン範囲を画定する二つの境界方向のいずれかが、座標軸から０度を超え９０度未満のオフセット角だけオフセットされていると好ましい。オフセット角は、たとえば４５度といった角度に設定されるが、整数Ｎの値に応じて異なる値が採用されることもある。また、同一の整数Ｎに対して複数のオフセット角が適する場合もある。なお、以下特に断りのない場合の検出素子の数Ｎは、検出装置に含まれる検出素子のうち有効なものの数を指すこととする。これらを含め本出願においては、不明瞭にならない限り本開示の属する技術分野の慣用に従う用語法を採用することがある。波は電磁波、音波といった伝播による減衰を観念できる任意の波動を含み、光波（赤外、可視、紫外）、振動（地震など）を含む。粒子は、流れとなった粒子線において伝播による減衰を観念できる任意の粒子を含み、中性子線やその他の放射線も含む。電子や光子のように、量子力学的に波動性と粒子性を兼ね備えるものも、これらの波または粒子のいずれかまたは両方に含まれる。同様に、古典力学的に波動や粒子と捉えられるものも、これらの波または粒子のいずれかまたは両方に含まれる。これらの波または粒子は、それらの人為的な放出源となる装置（照射装置）により生成されるものであっても、また、例えば宇宙線などの自然由来で生成されるものや、例えば対象物の各部から放射される放射線などの対象物の性質として生じるもの（本出願書面では「照射装置によらずに生じた波または粒子」と呼ぶ）であっても本開示の実施に用いることができる。
【発明の効果】
【００１３】
本開示のある態様では、相対的に少ない計算資源を用いる場合でも高い精細度で高品質な再構成像を高速に実現することができる。
【図面の簡単な説明】
【００１４】
【図１】図１は本開示の実施形態において画像が取得される一例の断層撮像装置において、画像が取得される対象物と検出装置の関係を示す斜視図（図１Ａ）、および画像が取得される対象物の切断面における平面配置を含む概略構成を説明するための模式図（図１Ｂ）である。
【図２】図２は従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。
【図３】図３は本開示の実施形態の離散ラドン逆変換の代数的手法の原理を説明する説明図であり、図３Ａは解像度Ｎが２の場合の画素と検出方向の関係を示し、図３Ｂは一般の解像度Ｎの場合の画像の座標軸と検出方向の関係を示す。
【図４】図４は、本開示の実施形態において、オフセット角に対してシステム行列のランク（階数）を計算した計算例であり、図３に示した解像度Ｎが２の場合のものである。
【図５】図５は、本開示の実施形態において、オフセット角に対してシステム行列のランク（階数）を計算した計算例であり、図５Ａ〜５Ｃの順に、解像度Ｎが４、６、８の場合のものである。
【図６】図６は、本開示の実施形態において、オフセット角に対してシステム行列のランク（階数）を計算した計算例であり、図６Ａ〜６Ｃの順に、解像度Ｎが１０、１２、１４の場合のものである。
【図７】図７は、本開示の実施形態において、ランクがＮ２となるオフセット角の例示の値を、各Ｎに対してプロットするグラフであり、奇数の二倍に当たる整数の解像度Ｎである場合（図７Ａ）、および偶数の二倍に当たる整数の解像度Ｎである場合（図７Ｂ）、および図７Ｂの各点に解析関数によるフィッティングカーブを重ねて表示するグラフである。
【図８】図８は本開示の実施形態の断層撮像装置のうち制御装置において実施される処理を概観するフローチャートである。
【図９】図９は本開示の実施形態の断層撮像装置においてオフセット角を再設定する場合の処理を概観するフローチャートである。
【図１０】図１０は断層撮像装置によって本実施形態を実装する場合における制御装置のいくつかの特有の機能部と特有のデータ格納部を示すブロック図である。
【図１１】図８は本開示の実施形態のための実施例において、検証のために用いた検証用データ（数値ファントム）を示す図である。
【図１２】図１２は本開示の実施形態のための実施例において得られた再構成画像であり、本実施形態により求めたもの（図１２Ａ）、従来のＦＢＰ法においてＨａｎｎフィルターを用いたもの（図１２Ｂ）、および従来の代数再構成法として逐次近似計算（ＭＬ−ＥＭ法）を３０回行って得たもの（図１２Ｃ）である。
【発明を実施するための形態】
【００１５】
以下図面を参照し、本開示に係る断層画像撮像の実施形態を説明する。全図を通じ当該説明に際し特に言及がない限り、共通する部分または要素には共通する参照符号が付される。また、図中、各実施形態の要素のそれぞれは、必ずしも互いの縮尺比を保って示されてはいない。
【００１６】
１．原理
本実施形態における例示の幾何学的配置を説明する。図１は本実施形態の画像が取得される一例の断層撮像装置１００において、画像が取得される対象物の切断面における平面配置を含む概略構成を説明するための模式図である。対象物（図１には示さない）に対し静止した空間の一平面にＮ画素×Ｎ画素を持つ画素の群が定義される。通例この整数Ｎは検出装置２０の検出素子２２の数と一致させている。このため、本実施形態において、整数Ｎを特に解像度Ｎと呼ぶこともある。波または粒子のビーム４は、その平面に含まれており、照射装置１０から照射される。検出装置２０と照射装置１０は、相対配置が固定されたまま対象物に対して回転することができる。このような制御を含む装置の動作や計算処理のために、断層撮像装置１００はコンピューター機器でありうる制御装置１５０を備えている。この回転は、上記画素の群やそこに静止している対象物に対する相対的なものである。このため必ずしも検出装置２０と照射装置１０のみが回転されるとは限らず、検出装置２０と照射装置１０が設置床面（図示しない）に対し固定されていて対象物が、画素の群とともに回転される場合もある。
【００１７】
照射装置１０により波または粒子のビーム４を照射し検出装置２０によりそのビーム４を対象物の反対側で検出する動作は、検出方向θのスキャン範囲を画素数Ｎによって分割した刻みで検出方向についてサンプリングされる。図１のようなビーム４が平行なものでは、従来は検出方向θでカバーする角度スキャン範囲は、例えば０〜１８０°とされ、検出装置２０からの検出信号が（１８０／Ｎ）°の刻みでサンプリングされる。その際、角度スキャン範囲の境界に近接するような検出方向は、たとえば再構成される画像の画素を指定するための座標軸に平行にしたり垂直にされる。検出方向がいずれかの座標軸に平行または垂直となっている場合、角度にはオフセットがない。本実施形態を含めて、断層撮像装置１００の検出方向ｄ１〜ｄＮの典型的なスキャン動作は、一定速度で増加または減少させる定速回転をしながら、検出装置２０の値を取得するタイミングを電気的に制御することによりサンプリングが行われる。ただし、本実施形態はこのようなものに限定されるものではない。
【００１８】
図２は、従来の離散ラドン逆変換の代数的手法の原理を説明する説明図である。この図は波または粒子を照射し、透過後の強度データを用いて断層画像を再構築する処理の本質を示す初等的なものである。図示される縦２×横２に並ぶ正方形領域は、解像度Ｎ＝２に対応させてそれぞれが断層画像の画素を意図したもので、ｆ１〜ｆ4の４つの値を持ちうる。これらの値は、各画素の位置での対象物によるビームの減衰の程度を示し、決定されるべき値である。複数の画素を通り、解像度Ｎ＝２に対応させて１８０°／Ｎつまり９０°向きを変えて描かれた白抜き矢印は波または粒子のビームの照射と検出素子の位置に対応している。矢印が指し示す先に記す数値は、検出装置の検出素子がその矢印の先に位置している際に得られる値（例えば減衰量）を意味しており、ここでは通った画素の値の合計が数値で示されている。つまり、この矢印が波または粒子のビームを表し、矢印が通る画素と矢印の向きで照射および検出（以下「照射」とのみ記す）を示している。矢印の先の合計値は、その矢印の先に検出装置の検出素子がある場合の、各検出素子で検出された減衰成分に対応している。矢印の指し示す値つまり検出装置から得られる値のみによってｆ１〜ｆ4の４つの値を決定することが、検出装置の検出素子の示す値に基づいて離散ラドン逆変換を実行することにより断層画像を再構成する処理に対応する。
【００１９】
図２に示された例では、ｆ１〜ｆ4のうち矢印それぞれがつないでいるものが線形和を与えるため、次の連立方程式に書き直せる。
【数１】
2019230741.tif 000003
ここで、ｘ線ＣＴの場合について例示すれば、ｐ１〜ｐ4は図２のｘ線の各画素および各検出方向での測定値であり、ｆ１〜ｆ4はｘ線吸収係数である。この連立方程式は一意に値を定めることはできない。ここで、
【数２】
2019230741.tif 000004
のように任意の数Ａを用いてシフトさせた新たなｆ′１〜ｆ′4を代入しても式（１）は不変である。つまり、式（１）は元（決定すべき変数）の数が４であるのに対し、矢印の示す検出によって導き出せる連立方程式の実質的な数は３である。この場合、ｆ１〜ｆ4の４つの値はいずれも不定（indeterminate）となり値を特定できない。
【００２０】
この関係を別の数学的関係で整理すると、上述式（１）は、
【数３】
2019230741.tif 000005
のように、より一般に行列の形式に表現される。ｗｎｍはＮ２行Ｎ２列の正方行列である行列Ｗ（「システム行列」と呼ぶ）の要素である。一般に、行列Ｗが逆行列Ｗ−１を持ちうることは、行列Ｗが正則であるとも呼ばれ、種々の数学的手法により判定することができる。そして、それら手法の一つの典型が階数の値による判定である。すなわち、行列Ｗが正則であることは、Ｎ２行Ｎ２列の正方行列Ｗの階数（ｒａｎｋ（ｗｎｍ））の値がＮ×Ｎに達することと同値である。行列Ｗが正則であれば式（１）に対応する連立方程式を一意に解くことができる。上記式（３）の行列では、階数を計算するとｒａｎｋ（ｗｎｍ）＜４となってあり、Ｎ＝２に対応して正則であるために必要な値に達していない。
【００２１】
図３は、本実施形態の離散ラドン逆変換の代数的手法の原理を説明する説明図である。図３Ａは解像度Ｎが２の場合の画素と検出方向の関係を示し、図３Ｂは一般の解像度Ｎの場合の画像の座標軸と検出方向の関係を示す。本実施形態は、最も端的には、図３Ａに示すように図２から検出方向を全体にオフセットさせてずらすよう設定して検出を行う。図３Ａでは、検出方向ｄ１、ｄ２が白抜き矢印で示される複数の画素の検出方向を代表しており、−４５°のオフセット角を設定している。
【００２２】
図３Ａの本実施形態の例でみると、上述式（１）に相当する式が
【数４】
2019230741.tif 000006
のように書くことができ、同様に階数を求めると、ｒａｎｋ（ｗｎｍ）＝４となって、Ｎ×Ｎに達する。つまり、本実施形態ではｆ１〜ｆ4を一意に決定することができる。実際、式（４）は
【数５】
2019230741.tif 000007
と整理され、その逆行列Ｗ−１（｛ｗｎｍ−１｝とも記載する）が
【数６】
2019230741.tif 000008
と算出される。つまり、測定値であるｐ１〜ｐ４にこの逆行列を作用させるという代数的に厳密な手法のみによって、ｆ１〜ｆ4を一意に決定することが可能となる。本発明者は、このような検出方向をずらすのみで逆行列を持つようなシステム行列を決定しうる現象を発見し、イニャシャル・フェーズ・イフェクト（ＩＰＥ）と命名した。ＩＰＥは、実際の測定処理では角度のサンプリングのタイミングをずらすなどにより、二次元再構成画像の画素を指定する座標軸に対してオフセットされた方向での検出動作を行いうることを意味しており、極めて実用性が高い手法である。
【００２３】
つづけて、ＩＰＥがより一般の解像度Ｎに対して成立するかどうかを確認した。図３Ｂは一般の解像度Ｎの場合の画像の座標軸と検出方向の関係を示しており、解像度Ｎにおける理解の容易さのために方向（座標軸、検出方向）のみに着目している。一般の解像度Ｎでは、図３Ａにおけるｆ１〜ｆ4のような決定されるべき値は、ｆ１〜ｆＮ×Ｎと表される（図示しない）。これらの値ｆ１〜ｆＮ×Ｎは、たとえば二次元直交座標で指定されたＮ画素×Ｎ画素の二次元断層画像の各画素に対応している。この画像は、図示しない対象物に対し、たとえば静止するように定められるため、その画像は座標軸であるｘ軸およびｙ軸に沿って広がりをもち、画素をその座標軸の値によって特定できる。本実施形態では、検出方向も、検出範囲をＮ分割した検出方向ｄ１〜ｄＮの方向で行われる（図３Ｂ）。本実施形態では、検出方向ｄ１〜ｄＮのいずれもが、ｘ軸およびｙ軸のいずれとも、平行でも垂直でもないようにして行われる。上述したオフセット角は、より一般には、検出スキャン範囲を画定する二つの境界方向に（角度において）最近接する二つの検出方向と、いずれかの座標軸との間で決定される。再び図３Ａを例にすれば、オフセット角Δθは、検出スキャン範囲がｘ軸の負の方向から時計回りにｘ軸の正の方向に向かうまでの１８０°の範囲とすることができ、ｘ軸の負の方向に最近接する検出方向ｄ１と座標軸であるｙ軸の負の方向（紙面上の上方向）との間で−４５°とされている。留意すべきは、検出スキャン範囲の方向の範囲内に検出方向が含まれているものの、オフセット角は、検出スキャン範囲よりはむしろ、検出方向と画像の座標軸との間に対して決定されることである。これは、検出スキャン範囲の取り方がある程度の許容幅があることと関係している。たとえば、図３Ａの検出方向ｄ１、ｄ２を固定したままでも、検出スキャン範囲は、ｘ軸の負の方向から時計回りにｘ軸の正の方向に向かう１８０°の範囲（図３Ａの紙面上方を１２時方向とした時計文字盤の方向で９時から３時の範囲）とみることもできるし、これとは異なり、９時方向と１２時方向の二等分線の方向（１０．５時つまり１０時３０分方向）から、３時方向と６時方向の二等分線の方向（４．５時つまり４時３０分の方向）までの１８０°の範囲とみることもできる。一般に、検出方向ｄ１〜ｄＮは、検出スキャン範囲を均等にＮ分割（Ｎ等分）する角度刻みそれぞれに対応するサンプリング方向となっていればよい。図３Ｂに示すように、一般の解像度Ｎの場合には検出素子数Ｎに応じ検出方向ｄ１〜ｄＮが多数設定される。本実施形態では、検出方向ｄ１〜ｄＮのいずれもがｘ軸およびｙ軸のいずれとも平行でも垂直でもないように設定される。この場合、オフセット角Δθは検出方向ｄ１または検出方向ｄＮに対して、ｘ軸またはｙ軸のいずれかに対して決定される。なお、本実施形態では、画像のための座標軸と検出方向との間に相対的にオフセット角Δθが定められる限り、画像のための座標軸や検出方向それ自体は種々の都合で定めても良い。たとえば、装置のもつ座標軸に画像のための座標軸が一致されている場合には、検出方向を調整してオフセット角Δθが実現される。逆に、装置のもつ座標軸に検出方向が一致されている��合に、画像のための座標軸を回転させることによりオフセット角Δθを実現することもできる。
【００２４】
解像度Ｎに一般化してＩＰＥが有効であるかについて、数値計算によって確認した。図４は、解像度Ｎ＝２の場合において、オフセット角に対してシステム行列のランク（階数）を計算した計算例である。横軸はオフセット角を、ラジアン単位（rad）で取っており、縦軸は各オフセット角の値のときの階数である。計算は、数式処理ソフトウエアMathematica（Wolfram Research, Inc., Champaign, Illinois）を採用し、オフセット角を０〜π／２（rad）（０〜９０°）の範囲でステップπ／３２０（rad）ずつ変化させて行った。これらの条件は、図の上方に明示している。図２に対応するのが横軸０またはπ／２（約１．５７）（rad）の位置である。また、図３Ａに対応するのが、π／４（rad）の位置である。オフセット角π／４（rad）では式（４）〜（６）に示されるようにこの場合にはシステム行列Ｗが逆行列Ｗ−１をもち、階数が４（つまりＮ×Ｎ）となっている。この数値的な確認をより一般の解像度Ｎにより行ない、図５、図６の結果を得た。図５Ａ〜５Ｃの順に、解像度Ｎが４、６、８のものであり、図６Ａ〜６Ｃの順に、解像度Ｎが１０、１２、１４のものである。図５Ａは、オフセット角をπ／３２０００（rad）刻みで計算したため、計算したオフセット角の範囲を７．５π／３２〜π／２（rad）に限定した。
【００２５】
図示した結果以外も含め、本発明者が確認したところ、ＩＰＥの有効性、つまり各解像度Ｎにおけるは、次のとおりである。
・解像度Ｎ＝４においては、階数＝Ｎ×Ｎとならないようである。
・次の解像度Ｎの場合には、階数＝Ｎ×Ｎとなることが確認された。
Ｎ＝２，６，８，１０，１２，１４，１６，１８，２０，２２，２４，２６，２８，３０，３２，３４，３６，３８，４０，４２，４４，４６，４８，６４，１２６，１２８，１３０
すなわち、オフセット角を導入することにより、多くの場合において、階数＝Ｎ×Ｎの条件が達成され、図２の場合に見られたような階数の低下（ランク落ち）が解消してしまうことを確認した。こうして多くの実用的な解像度Ｎにおいて、オフセット角を導入することにより、ＩＰＥが有効であることを確認した。
【００２６】
本発明者は、数値計算によって確認した結果に基づいて、解像度Ｎに応じて変化する階数の規則性に着目した。具体的には、オフセット角がπ／４（rad）（４５°）の近傍でランク落ちが解消していることを検討した。まず、対称性は比較的明確であり、オフセット角を０〜π／２（rad）の範囲で変化させた場合、図４〜６のグラフがπ／４（rad）を中心に左右対称の関係がある。また、解像度Ｎが２×奇数であるものと２×偶数であるものとの間に異なる傾向が見られることに気づいた。つまりＮが２×奇数である場合、階数が回復するオフセット角の範囲の一つは、２×偶数の場合と比べて相対的にπ／４（rad）に近いところに現われるようである。図７は、本実施形態において、ランクがＮ２となるオフセット角の例示の値を、各Ｎに対してプロットするグラフであり、奇数の二倍に当たる整数の解像度Ｎである場合（図７Ａ）、および偶数の二倍に当たる整数の解像度Ｎである場合（図７Ｂ）、および図７Ｂの各点に解析関数によるフィッティングカーブを重ねて表示するグラフである。
【００２７】
現時点では、計算資源上の制約により、ランク落ちが解消するＩＰＥの有効性が確認できた解像度Ｎは最大でもＮ＝１３０である。
表１は、階数＝Ｎ×Ｎとなるオフセット角が、オフセット角を０からπ／４（rad）と変化させて行く間に複数の値または範囲にて観察されるとき、最もπ／４（rad）に近い値となるようなオフセット角を解像度Ｎに対して示すものである。
【表１】
2019230741.tif 000009
【００２８】
より大きなサイズの解像度ＮでのＩＰＥの有効なオフセット角の出現位置を予測するために、図７Ｃのように計算した範囲にて良好に計算結果を説明できる数式
【数７】
2019230741.tif 000010
を用いると、たとえば
Ｎ= 512について、オフセット角＝0.779506（rad）(44.662°)
Ｎ=1024について、オフセット角＝0.779526（rad）(44.663°)
との値が求まる。
現時点で本発明者は、ＩＰＥの有効性の原理的側面についての以下の知見を得ている。
・オフセット角を調整することにより、値ｆ１〜ｆＮ×Ｎについての不定性は解消し値ｆ１〜ｆＮ×Ｎを一意に決定することができる。
・他方、一意に値ｆ１〜ｆＮ×Ｎを決定できない可能性の高い解像度Ｎも、一例（Ｎ＝４）見つかっている。
以上の原理的側面に従えば、適切なオフセット角を採用する本実施形態では、代数的解法による厳密な画像再構成が可能となる。
【００２９】
２．具体的構成
次に、上記原理を採用する本実施形態の画像再構成法の具体的構成について説明する。図１に示したとおり、ビーム４は上記座標にて全ての画素範囲をカバーする広がりを持ち、各検出素子２２に向かう部分はある幅τを持っている。当該検出素子２２に対するビーム４についてのその画素のもつ寄与は、典型的にはビーム４のうち、検出素子２２の一つに向かうものがもつ幅τと各画素（縦横それぞれδ）の重なりである。図１では一つの画素についてこの重なりの様子を模様で示している。Ｎ画素×Ｎ画素の全ての画素について、検出方向ｄ１〜ｄＮ別、検出素子２２の別に寄与が定まる。この寄与を数値として保持しているのが上述したシステム行列Ｗである。式（３）をシステム行列Ｗ２より明示的に表現すると、照射と画像再構成のうち、照射を表す離散ラドン変換が次の行列の積により表現される。
Ｘ′＝  Ｗ  Ｘ            （８）
と表現される。ここで、システム行列Ｗは、総画素数に応じた空間についてのサンプリング点の並びを行方向に、検出方向の刻みと検出素子とに応じた照射および検出についてのサンプリング点の並びを列方向にとるのが通常である。図１の画素数などの設定に適用した場合には、システム行列Ｗが行方向にも列方向にもＮ×Ｎ個の要素を持つ正方行列となりその要素の総数はＮ４となる。これは、画素がＮ画素×Ｎ画素、検出方向θのスキャン範囲がＮサンプル、検出素子２２の数がＮ素子であるためである。なお、画質を高めるためにＮの値は時代とともに増大しており、典型例では１０２４すなわち１０３程度とされ、従来のものでは２０００程度が最大値である。本実施形態の場合にもＮの値が得に限定されるものではない。Ｎが１０３程度のものではシステム行列は行方向列方向ともに１０６程度の規模、要素数は１０１２程度となる。
【００３０】
式（８）のＸは各画素の位置での減衰（吸収、散乱も含む）を表す像である減衰像を列ベクトルとして表現したものであり、式（３）におけるｆｎｍである。実際の空間での減衰像は図１に示したＮ画素×Ｎ画素の二次元画像であるが、システム行列Ｗとの演算のために、行方向を総画素数Ｎ×Ｎになるようにベクトル化したものがＸである。式（１）において減衰像Ｘにシステム行列Ｗが作用して得られる左辺の列ベクトルＸ′は、サイノグラムと呼ばれる像に対応している。システム行列Ｗにて列方向を検出方向のサンプリング点と検出素子の並びとしたことに対応して、列ベクトルＸ′も検出方向のサンプリング点と検出素子それぞれについて得られる検出強度を示している。なお、サイノグラムは検出素子の軸と検出方向刻みの軸との２次元に列ベクトルＸ′の要素を並べ替えたものである。以上のように、式（８）は照射と検出を表現したものである。
【００３１】
断層画像の再構成とは、式（８）が実現していることを前提として、システム行列Ｗとサイノグラムの列ベクトルＸ′から列ベクトルＸを算出し、列ベクトルＸから減衰像を得る処理である。つまり断層画像の再構成は形式的に次式により表現される。
Ｘ＝  Ｗ−１  Ｘ′            （９）
ここで、Ｗ−１はＷの逆行列であり、式（９）は式（８）の両辺にＷ−１を作用させて左右辺を入れ替えたものである。ここの式（８）と（９）との関係は、両立した式（５）と（６）の関係である。そして本実施形態においては適切にオフセット角を設定することにより、正方行列であるＷがほとんどのＮの場合に逆行列Ｗ−１をもち正則であるようにすることができる。本実施形態では、最終的に式（９）が成立するよう、逆行列Ｗ−１を持ちうる正則な正方行列であるシステム行列Ｗを十分な確かさで決定できる手法が提案され、さらにそのための照射、検出、および断層画像再構成のための十分に実用性のある具体的構成も提供される。なお上述した図１に基づく説明は平行ビームを照射するものであるが、本実施形態はファンビームやコーンビームといった平行ではないビームを利用するような照射装置１０や検出装置２０を採用する場合にも当業者には明らかなわずかな変更によって適用することができる。また、図１の画素数などの設定以外においても、本実施形態の手法を適用することができる。なお、式（８）、（９）を、数学的に等価な他の形式により表現することもでき、たとえば、システム行列Ｗのために、上述したものとは異なり、総画素数に応じた空間についてのサンプリング点の並びを列方向に、検出方向の刻みと検出素子とに応じた照射および検出についてのサンプリング点の並びを行方向にとる。そしてそれに対応させて、式（８）のＸを各画素の位置に応じて行ベクトルとして表現することもできる。その場合式（８）、（９）に相当する式は、それぞれの各辺に転置操作を施したものとなる。この転置操作では行列やベクトルの行と列を入れ替え、ベクトルに作用させる行列が左からではなく右から作用させる順序に変わるという表現上の影響がある。ベクトルが「列または行ベクトル」のように択一的に記載され、対応する行列の要素の並びが「列方向または行方向」のように択一的に記載されているときは、これらの組み合わせの選択は、それぞれの択一的記載における同順のものを組み合わせて選ぶこととする。この表現の選択は本出願の本質に影響を与えないため、本願においては特に断りの無い場合には式（８）のようにサイノグラムをベクトル化したＸ等を列ベクトルにより表現し、システム行列等の行列は左から作用させる表記にもとづいて説明する。
【００３２】
２−１．処理フロー
本実施形態で正則なシステム行列Ｗが得られるのは、適切なオフセット角を採用するためである。図８は、これらを含む本実施形態の断層撮像装置１００のうち制御装置１５０において実施される処理を概観するフローチャートであり、図９は、本実施形態の断層撮像装置においてオフセット角を再設定する場合の処理を概観するフローチャートである。また、図１０は、断層撮像装置によって本実施形態を実装する場合における制御装置１５０のいくつかの特有の機能部と特有のデータ格納部を示すブロック図である。
【００３３】
２−２．処理の概要
本実施形態において実行される処理を概観すると、まず装置の構成が決定される（図８、Ｓ１０２）。後の撮像でＮ画素×Ｎ画素（Ｎは正の整数）の二次元画像を再構成する場合、検出装置２０の検出素子２２の数もＮ個とされる。サンプリングされる検出方向ｄ１〜ｄＮは、従来の手法では、角度スキャン範囲が例えば（１８０／Ｎ）°刻みで、例えばｘ軸またはｙ軸（図３Ａ、図３Ｂ）の正または負の向きをいっぽうの境界に、他方の境界を、そこから時計回りに１８０°回ったむき、といったように、画像の軸にあわせて、０〜１８０°といった範囲に設定される。これに対し本実施形態の検出方向のサンプリングは、例えば（１８０／Ｎ）°刻みで、オフセット角Δθ（たとえば±４５度に近い角度）だけｙ軸（図３Ｂ）の負の向きからオフセットされた、Δθ〜Δθ＋１８０°の範囲、のように設定される。検出方向ｄ１〜ｄＮそれぞれは、（１８０／Ｎ）°ずつ方向がずれており、上記オフセット角は、検出方向ｄ１〜ｄＮの最初または最後の検出方向ｄ１または検出方向ｄＮが、ｘ軸またはｙ軸とのいずれかの間で作る角度として定義される。本実施形態において上記オフセット角は、−９０°を超え０°未満、または０°を超え９０°未満となるようにされる。
【００３４】
対象物２から断層画像を取得して適切な再構成像を得るためには、システム行列を決定しなくてはならない。このために、オフセット角のために値が異なる以外は従来のシステム行列と同様の処理を行ってシステム行列Ｗを生成する（Ｓ１１２）。ここで求まるシステム行列Ｗは、検出方向と検出素子それぞれに対応させて列方向にＮ×Ｎ、断層画像の画素数に対応させて行方向にＮ×Ｎだけの要素をもつ。システム行列はシステム行列格納部６０（図１０）に格納される。この段階でオフセット角Δθが適切であれば、システム行列Ｗは正則行列となる。システム行列Ｗからは、離散ラドン逆変換行列算出ステップＳ１１６にてその逆行列Ｗ−１が算出され、離散ラドン逆変換行列格納部４０に格納される。本実施形態では逆行列Ｗ−１を離散ラドン逆変換行列とも呼ぶ。ここで逆行列Ｗ−１を算出するための処理には、例えば余因子法や掃き出し法など任意の手法を採用することができる。
【００３５】
一方、対象物からの断層画像を得るには、照射装置１０と各検出素子２２をもつ検出装置２０との間に、対象物２を配置する。オフセット角を利用した照射および検出のために、照射装置１０からのビーム４を照射しながら、Ｎ方向で検出が行われる（Ｓ１２４）。典型的な実装形態では、このビーム４は連続して照射し続けていて、対象物２と照射装置１０および検出装置２０の向き（検出方向）はなめらかに動きつつ、検出装置２０の検出素子２２からの信号が電気的にサンプリングされる。Ｎ方向での照射および検出の典型例は、画像の座標軸に対してオフセット角Δθを持つようなタイミングが調整された（１８０／Ｎ）°刻みの繰り返しのサンプリングである。得られる信号は、必要に応じて機器の物理的特性について補正やスケールの換算などの処理と配列順序の整理などをすれば、サイノグラムを生成しうるデータに整理される（Ｓ１２６）。整理したものは、オフセット角Δθが反映された検出方向と検出器の位置とに対応したものである。このサイノグラムを得る処理が、照射および検出の処理であり、離散ラドン変換である。このサイノグラムからは、システム行列Ｗの列方向の並びに合わせるような順序で要素を並べ替えた列ベクトルを得ることができる（ベクトル化Ｓ１２８）。このベクトルを本実施形態では第１ベクトルと呼ぶ。第１ベクトルは、サイノグラムの値を持つのでＮ×Ｎ要素を持ち、式（１）と同様に、断層画像が得られていて、それをベクトル化したもの（Ｎ×Ｎ要素の列ベクトル）にシステム行列Ｗを作用させたとしたならば得られる列ベクトルといえる。第１ベクトルはＸ′とも記す。なお、照射および検出Ｓ１２４からベクトル化Ｓ１２８までは、信号のアナログデジタル変換による量子化、デジタル信号の処理や、必要に応じて一時的にデータを保存して処理することもできる。また、サイノグラムは、それを生成しうるデータが得られれば十分であり、必ずしもその表示を行う必要はなく、サイノグラムといえる明示的なデータを経由する必要もない。サイノグラムをベクトル化したものに相当する第１ベクトルが検出装置２０からの検出信号に基づいて得られれば、照射および検出Ｓ１２４以降ベクトル化Ｓ１２８まで等価な任意の処理を行う任意の処理も本実施形態に含まれている。なお、上記オフセット角が反映された検出方向を実現するためには、オフセット角の分だけ通常の角度（たとえば装置が基準として持つ方向）からずれるように装置の検出動作を設定したり制御することは必ずしも要さず、数学的に等価な処理により本実施形態を実施することもできる。たとえば、装置が基準として持つ方向に対して検出方向は一致させておき、むしろ画像のための座標軸に必要な回転変換を施してオフセット角Δθを実現して再構成までを実行し、最後に再構成画像にその分の逆回転変換を施すことも本実施形態のために有用である。
【００３６】
第１ベクトルＸ′はシステム行列Ｗとの間で式（８）の関係を満たすため、システム行列Ｗの逆行列である離散ラドン逆変換行列Ｗ−１を用いて式（９）の関係が成り立ち、離散ラドン逆変換処理Ｓ１３２によって再構成像を与える第２ベクトルが得られる。この第２ベクトルはＸとも記す。この第２ベクトルＸの並びをベクトル化の逆操作によってＮ画素×Ｎ画素の画像に戻すデベクタライズ処理により再構成像を画像化することができる（Ｓ１３４）。これで１スライス分の断層画像が再構築されるので、必要に応じてスライス位置を変更した処理を実行することができる。
【００３７】
本実施形態の処理は図８においてＡとＢに示す２系統に分かれている。通常再構成像を撮像するユーザーは、照射および検出ステップＳ１２４〜画像化ステップＳ１３４に至る処理の系統Ｂを実行する。つまりユーザーは、装置構成の決定Ｓ１０２、システム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａを必ずしも実行する必要はない。システム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａが必要となるのは、断層撮像装置を作製したり、設置やメンテナンスしたりする場面で、いわば機器メーカーが対応するような処理である。予め決められたオフセット角Δθが得られていれば、再構成像を撮像しようとするユーザーは自らの測定のためにはそれらを利用するだけでよい。系統Ａからの適切な機器情報についてのデータが得られていれば、それを使用すれば系統Ｂの範囲の処理は十分実行できる。
【００３８】
特筆すべきは、系統Ｂのうち離散ラドン変換Ｓ１２６〜画像化ステップＳ１３４が行列やベクトルの順序の入れ替えや行列積の演算であり、例えば逐次近似といった計算処理量が問題となる処理がその範囲に含まれていないことである。本実施形態において、ユーザーの通常の撮像において利用される系統Ｂの処理は、系統Ａの処理により準備された離散ラドン逆変換行列格納部４０を使うだけで十分に軽い計算負担で実行可能であり、この点は、本実施形態に高い実用性をもたらすものといえる。断層画像を取得するという最も頻度の高い処理の計算負担がごく少ないためである。特に、三次元ボリューム像を得るために多数のスライスの撮像を繰り返す用途には本実施形態の手法が極めて有利である。その繰り返しに系統Ａの処理が不要なことに加え、系統Ｂの離散ラドン変換Ｓ１２６〜画像化ステップＳ１３４はＧＰＧＰＵ（General-purpose computing on graphics processing unit）などのメニーコアタイプのプロセッサーでの計算処理に適するため、今後の処理速度の向上も期待できるからである。
【００３９】
また、幾何学的な側面を反映するシステム行列Ｗやシステム行列Ｗの逆行列Ｗ−１がシステム行列生成ステップＳ１１２〜離散ラドン逆変換行列算出ステップＳ１１６に至る処理の系統Ａのみに含まれており、測定のための系統Ｂ（図４）と分離されていることから、代数的手法に反復を利用したときの収束の問題と測定のためのノイズの側面とが本質的に分離されている。このため、本実施形態では、ノイズとは無関係な近似解に起因した現象が生じにくく、画像処理で一般に採用される各種のノイズ軽減処理が目的通りに効果を発揮しやすいという利点もある。
【００４０】
２−３．正則なシステム行列を与えるオフセット角Δθの決定
概要において上述したように、離散ラドン逆変換行列Ｗ−１がシステム行列Ｗの逆行列として求まることは、オフセット角Δθの適切さに依存している。そのオフセット角Δθは、図８には図示しない処理において決定されている。本実施形態のオフセット角Δθ決定の処理は、結果として正則なシステム行列を与えるものである限り限定はされない。
【００４１】
オフセット角Δθを探索する手順は図９のように行われる。まず図７Ｃを通じ説明したように、現実の解像度Ｎの値の場合には±π／４（rad）（±４５°）付近を初期値に用いて探索すれば少ない手順で見つかることが期待できる。オフセット角Δθの初期値を設定した後、システム行列Ｗを算出する（Ｓ４２）。そのシステム行列Ｗを対象にして計算した階数が解像度Ｎの値から求めたＮ×Ｎと一致するかを判定する（Ｓ４４）。一致しない場合、オフセット角Δθの値を変更し（Ｓ４６）、一致した場合には、その時点での値をオフセット角Δθに設定し（Ｓ４８）をそのまま採用する。そしてシステム行列Ｗの生成（Ｓ４２）に戻る。このような処理により、確実に逆行列を求めうるシステム行列Ｗとそのためのオフセット角Δθを決定することができる。なお、図９に示した処理は、検出装置２０の検出素子２２の数である解像度Ｎを決めオフセット角Δθを一度決定すれば、その後は不要となる。このためユーザーは撮影の度にオフセット角Δθを決定したりする必要は無い。解像度Ｎごとに与えられるオフセット角Δθの値を何らかの手段により得てそれをそのまま利用したとしても、本実施形態の実施には何ら支障は生じない
【００４２】
２−４．具体的構成のまとめ
以上に示した通り、本実施形態にて用いるオフセット角Δθを利用する処理は、断層画像の撮像および取得のプロセスに高い実用性をもたらすものである。またその利点の一つとして、オフセット角Δθを決定したり、それに応じたシステム行列やその逆行列を生成する処理はユーザーが断層画像を取得する処理とは別個に行うことができる。このため、本開示の実施形態においては、逐次近似の処理、といったが計算量の多い処理は繰り返す必要もない。
【００４３】
３．本実装形態の実装
本実施形態は、断層撮像装置として実装することも、また既存の断層撮像装置のための、制御および処理ソフトウエアとして実装することもできる。
【００４４】
３−１．断層撮像装置での実装
本実施形態を断層撮像装置１００は、上述した処理を主に制御装置１５０において行うことにより、断層撮像装置１００として従来の断層撮像装置において、従来のＡＲＴを用いたものと同等またはそれ以上に良質な画像の断層画像を得ることができる。その際に必要となる計算資源も十分に実用的なものであり、また処理が高速である。
【００４５】
図１０に示したようにシステム行列生成部１６２および逆行列算出部１６６は、それぞれ、システム行列生成ステップＳ１１２、および離散ラドン逆変換行列算出ステップＳ１１６を実行する。同様に、動作制御部１７４、離散ラドン変換処理部１７６、ベクトル化処理部１７８、離散ラドン逆変換処理部１８２、および画像データ生成部１８４は、照射および検出ステップＳ１２４、離散ラドン変換Ｓ１２６、ベクトル化Ｓ１２８、離散ラドン逆変換処理Ｓ１３２、および画像化ステップＳ１３４をそれぞれ実行する。その動作に当たり、離散ラドン逆変換行列格納部４０、およびシステム行列格納部６０も利用される。また、照射および検出Ｓ１２４からベクトル化Ｓ１２８までが等価な処理で実施できることと同様に、照射および検出ステップＳ１２４、離散ラドン変換Ｓ１２６、ベクトル化Ｓ１２８は、機能上これらを組み合わせた処理を行いうる機能部により実施することができる。さらに、制御装置１５０は、系統Ａの処理と系統Ｂの処理で必要となる機能部が明瞭に分かれていることから、必要に応じ系統Ａの処理を主として行うものであったり、逆に系統Ｂの処理を主として行うものであったり、様々な実装形態をとることができる。断層撮像装置の具体例としては、Ｘ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、光学ＣＴ装置が挙げられる。
【００４６】
３−２．コンピュータープログラムでの実装
また、断層撮像装置１００は、例えば従来のＡＲＴを実行するために用いられる照射装置１０検出装置２０を備えている断層撮像装置において、例えばコンピュータープログラムの実装形態として、断層撮像装置１００に備わる制御装置１５０に、図８、図９に示したフローチャートの一部または全部の制御ステップや変換、算出等の計算処理を実行させることによって、その全体を断層撮像装置１００として実施することができる。本実施形態では、さらに、ユーザーが断層画像を撮像するタイミングでは、断層撮像装置１００に備わる制御装置１５０には系統Ｂ（図８）の処理のみを実行させることもできる。
【００４７】
４．計算機シミュレーションによる検証（実施例）
本実施形態による再構成手法の有効性を計算機シミュレーションにより確認した。逆ラドン変換を含む処理全般は、対比のための従来の手法のものも含めて、数式処理ソフトウエアMathematica（Wolfram Research, Inc., Champaign, Illinois）および画像処理ソフトウエアＩｍａｇｅＪ（アメリカ国立衛生研究所）を用いて実現した。本実施形態の制御装置１５０の処理も含め、シミュレーション処理のために、Ｉｎｔｅｌ社製Ｘｅｏｎ  Ｅ５プロセッサー（３．５ＧＨｚ、６―Ｃｏｒｅタイプ）のＣＰＵと６４ＧＢの主記憶装置とを搭載したコンピューターを採用した。図１１は、検証のために用いた検証用データ（数値ファントム）を示す図（図１１、原画像という）である。各図には背景と区別するための外形を縁取る実線を画像の周辺に追加し、画素の座標を表す画素番号の最初および最後を示す数値を上辺および左辺に示している。原画像（数値ファントム）は、計算機シミュレーションにおいて対象物のある断面における減弱率構造を模擬するために人為的に与える数値データであり、画素が空間位置に対応したＮ＝１２８つまり１２８×１２８の画像として表現できる。図１１では吸収の強い領域を暗く、弱い領域を明るく表示している。計算機シミュレーションでは、断層撮像装置１００による照射および検出Ｓ１２４の動作も模擬し、計算によって測定値に相当する数値データを得た。本手法ではＮ方向（１２８方向）で取得したサイノグラム（図示しない）が取得される。
【００４８】
本実施形態の検証として、解像度Ｎを１２８とした場合のオフセット角Δθは、表１に示すとおり０．７７７０（rad）(４４．５１８８°)とした。図１２に示す再構成画像は、本実施形態により求めたもの（図１２Ａ）、従来のＦＢＰ法においてＨａｎｎフィルターを用いたもの（図１２Ｂ）、および従来の代数再構成法として逐次近似計算（ＭＬ−ＥＭ法）を３０回行って得たもの（図１２Ｃ）である。ＦＢＰ法、およびＭＬ−ＥＭ法を適用するにあたり、本手法で得られるサイノグラムに対しオフセット角Δθに相当する投影数分（３２投影）をシフトし、最終的に得る再構成画像に生じる画像の回転（角度Δθ）を打ち消した。本実施形態での再構成画像は１２８×１２８画素の範囲で得られるものの、これらの再構成画像は９１×９１画素の範囲のみ示している。これは、ＦＢＰ法でシミュレーションにより再構成像が得られる範囲が１２８画素を直径とする円形の領域となり、その円に概ね内接する正方形領域が比較可能となるためである。
【００４９】
図１１と図１２Ａを比較すれば明らかなように、本実施形態においては再現性が極めて高い再構成画像が得られた。これは、従来のＨａｎｎフィルターを用いたＦＢＰ法（図１２Ｂ）、およびＭＬ−ＥＭ法（図１２Ｃ）とは異なり、離散ラドン逆変換行列Ｗ−１が代数的な厳密解となることに裏付けられている本実施形態の利点であり、それが画像自体の高い再現性を通じ確認されている。
【００５０】
５．画像再構成手法のバリエーション
上述した画像再構成手法では、適切なオフセット角を採用することによりシステム行列Ｗが正則であって逆行列をもつことが保証された下で、正則なシステム行列Ｗから逆行列Ｗ−１を実際に算出して画像再構成を実施した。本開示では、追加の画像再構成手法として、同様の適切なオフセット角を採用して正則なシステム行列Ｗの擬似逆行列を利用しても画像再構成を行いうる。
【００５１】
本追加の画像再構成手法のための擬似逆行列は、例えば特異値分解（singular value decomposition; SVD）の手法によって算出される。例えば、ムーア・ペンローズ(Moore-Penrose)型逆行列は正方行列に対して特異値分解により生成可能であり、本追加の画像再構成手法のための擬似逆行列として有用なものの一例である。
【００５２】
本開示の実施形態では、上述した画像再構成手法および本追加の画像再構成手法の双方で、上記式（８）が利用される。この際、適切なオフセット角を採用することによりシステム行列Ｗは正則であることすなわち逆行列をもつことが保証されている。上述した画像再構成手法は、この正則なシステム行列Ｗの逆行列である離散ラドン逆変換行列Ｗ−１を求めて式（９）に従って画像再構成を実施するものであった。
【００５３】
本追加の画像再構成手法では、式（８）までは同等であるものの、式（９）に代えて、
Ｘ＝  Ｗｐ−１  Ｘ′      （１０）
に従って画像再構成を実施する。ここで、Ｗｐ−１はシステム行列Ｗの擬似逆行列である。つまり、適切なオフセット角を採用することにより逆行列の存在が保証されているシステム行列Ｗの逆行列を敢えて使わず、近似的な逆行列である擬似逆行列Ｗｐ−１を利用しても、本追加の画像再構成手法では、画像再構成を行いうる。
【００５４】
本追加の画像再構成手法のための処理を、必要に応じて図８〜図１０の要素を付記しながら説明すると、最初に、Ｎ方向の検出方向について適切なオフセット角を採用する（図８、Ｓ１０４）。ついで、システム行列Ｗを算出する（Ｓ１１２）。このシステム行列Ｗが正則であることは、階数の計算により確認する（図９）。さらに、当該システム行列Ｗのムーア・ペンローズ型逆行列を算出して擬似逆行列Ｗｐ−１を算出して、それを離散ラドン逆変換行列Ｗ−１とする。このため、図８において、離散ラドン逆変換行列算出ステップＳ１１６では、擬似逆行列が算出され、離散ラドン逆変換行列格納部４０（図８、図１０）にはその擬似逆行列が格納される。
【００５５】
次に、上記Ｎ検出方向で目的とする対象物からデータを取得する（Ｓ１２２、Ｓ１２４）。ついでサイノグラムに対応する第１ベクトルＸ′を生成する（Ｓ１２６、Ｓ１２８）。次に、擬似逆行列Ｗｐ−１により、第２ベクトルＸを算出し（Ｓ１３２）、さらにデベクタライズ処理により再構成画像を得る（Ｓ１３４）。擬似逆行列Ｗｐ−１の算出は、後述するトレランスを適切に設定して実施することができる。このように擬似逆行列Ｗｐ−１を利用する追加の画像再構成手法で良好な品質での再構成像を得ることができる。
【００５６】
擬似逆行列Ｗｐ−１の算出は、典型的には適切なトレランス値を設定して行われる。トレランス値は、特異値分解の近似を通じて擬似逆行列を生成する処理において、特異値数列において特異値を除外するための基準値を、特異値の最大値との比率で示す。例えば、１／５をトレランス値とした場合には、特異値の最大値の１／５より小さい値を持つ特異値を０で置き換えて、特異値分解した行列から擬似逆行列を算出する。本発明者がＮ＝６４の例で検討したところ、トレランスを１／１とすると再構成画像は認識できる像を構成しない。トレランスを小さくするにつれて再構成画像は解像度を増してゆくが、１／１３以降は比較的安定している。
【００５７】
以上のように、適切なオフセット角を採用することにより正則で逆行列を持つことが保証されているシステム行列Ｗを対象にした場合であっても、擬似逆行列を利用することによって、十分な品質の再構成画像が生成される。
【００５８】
なお、擬似逆行列のためのトレランス値は、一般に、小さい値を設定するほど算出される擬似逆行列が逆行列に近づき、再構成画像の解像度も高くなる。しかし、現実の測定装置を用いた測定においては、小さければ良好な再構成画像が得られるわけではない。幾何学的要因に帰着できない種々の攪乱要因が存在しているためである。例えば、検出信号に含まれるノイズや、ビームの不安定性、ビームの散乱、装置自体が持っている画像特性に起因した要因が、現状の測定装置に攪乱要因として残存している。通常これらの要因をすべて反映した装置全体の撮像性能が、例えば解像特性を示すＭＴＦ（Modulation Transfer Function）や、画像の粒状度の指標となるウィーナー・スペクトルによって評価されている。本追加の画像再構成手法において０より大きな所定のトレランス値を用いることにより、幾何学的要因以外の攪乱要因を除外するフィルター動作が実現する。このフィルター動作は、トレランス値を増減することによりその効果を容易に調整できることから、解像度を重視するか、定量性を重視するかなどの着目する画像特性に応じて必要な再構成像を作成することに役立ちうる。
【００５９】
６．変形例
本実施形態は、断層画像データ取得方法、断層画像データの取得装置やそのためのコンピュータープログラムのいずれで実施する場合であっても、上述したもの以外の種々の形態で実施することができる。例えば、図８、１０において系統Ａにおいて機器情報を決定する処理のためにユーザーが関与することは必ずしも必要なく、ユーザーを系統Ｂの測定情報を取得することのみに関与させるような実施態様への変形も容易である。さらに、例えばコンピュータープログラムで提供される場合には、照射および検出を行うハードウエアは既存のものや既に設置済のものを採用し、制御用コンピュータープログラムに上述した処理を新たに組み込むことで本実施形態を実装することもできる。このように、本実施形態の手法は、種々の実施態様を通じて実施することができる。
【００６０】
また、実施形態の変形に加え、その処理の細部を変形したものも本実施形態の一部をなす。その一例が、図１を参照して説明した検出素子２２に対するビーム４についてのその画素のもつ寄与の決定手法である。図１で模様を付した部分の面積を寄与とする手法に加え、他の実用に付されている手法は本実施形態においても採用できる。例えば、面積に代え、掃く長さに応じた重みを持たせたり、その掃く長さを効率よく評価するシドンのアルゴリズム（Siddon's algorithm）を採用して上記寄与を決定することも有利である（非特許文献３）。さらに、さらなる高画質化を行うために、アンチエイリアス処理できるSunnegardh and Danielsonの方法を利用することも有利である（非特許文献４）。
【００６１】
さらに、図１以降の説明は、一例の断層撮像装置１００の動作により説明したが、本実施形態は、必ずしも照射装置１０を採用するものには限定されない。ＳＰＥＣＴ（単一光子放射断層撮影）など、照射装置を用いることなく対象物の断層撮影を行う他の手法においても、本実施形態を同様に適用することができる。
【００６２】
以上、本開示の実施形態を具体的に説明した。上述の実施形態、変形例および実施例は、本出願において開示される発明を説明するために記載されたものであり、本出願の発明の範囲は、請求の範囲の記載に基づき定められるべきものである。実施形態の他の組合せを含む本開示の範囲内に存在する変形例もまた、請求の範囲に含まれるものである。
【産業上の利用可能性】
【００６３】
本開示は、波または粒子が例えば対象物に照射されたりすることによって、対象物自体を波または粒子が透過するような任意の断層撮像手法のために利用可能であり、例えば、Ｘ線ＣＴスキャナー、ＳＰＥＣＴ（Single Photon Emission Computed Tomography）、光トモグラフィー(Optical Tomography)、光学ＣＴ装置が挙げられる。
【符号の説明】
【００６４】
１００  断層撮像装置（断層画像データの取得装置）
２  対象物
４  ビーム
１０  照射装置
２０  検出装置
２２  検出素子
４０  離散ラドン逆変換行列格納部
６０  システム行列格納部
１５０  制御装置
１６２  システム行列生成部
１６６  逆行列算出部
１７４  動作制御部
１７６  離散ラドン変換処理部
１７８  ベクトル化処理部
１８２  離散ラドン逆変換処理部
１８４  画像データ生成部
【図１】
2019230741.tif 000011
【図２】
2019230741.tif 000012
【図３】
2019230741.tif 000013
【図４】
2019230741.tif 000014
【図５】
2019230741.tif 000015
【図６】
2019230741.tif 000016
【図７】
2019230741.tif 000017
【図８】
2019230741.tif 000018
【図９】
2019230741.tif 000019
【図１０】
2019230741.tif 000020
【図１１】
2019230741.tif 000021
【図１２】
2019230741.tif 000022
【国際調査報告】
2019230741.tif 000023
2019230741.tif 000024
2019230741.tif 000025