JP6494218B2 - 劣決定系線型方程式の求解装置、被検体情報取得装置 - Google Patents

Description

本発明は、連立線型方程式の求解方法及び求解装置に関し、特に、変数の数に対して方程式の数が少ない劣決定系線型方程式の求解技術に関する。

Ｘ線ＣＴ（Computed Tomography）やＭＲＩ（Magnetic Resonance Imaging）等の医療
機器では、人体内部の３次元構造を、２次元平面への投影観測画像から再構成する。推定したい３次元構造を未知ベクトルｘ、Ｘ線投影のような観測過程を行列Ａ、観測信号をｂとおくと、この観測過程は次のように定式化できる：

Ａは係数行列、ｂは定数ベクトルあるいは右辺ベクトルとも呼ばれる。

式（１）の方程式を求解する場合、３次元構造の推定精度を向上させるために観測の回数を多くすることが一般的である。以後、推定すべき信号ｘの次元をＮ、観測信号ｂの次元をＭとし、Ｍ＜Ｎを仮定する。式（１）は、変数の次元Ｎに対して方程式の数Ｍが少ないので、劣決定系線型方程式と呼ばれる。
劣決定系線型方程式は解が無数に存在するので、一般に正則化と呼ばれる技術を適用することによって解を一意に決定する。

近年、推定すべき信号のスパース性を仮定する（つまり、ベクトルｘの多くの要素がゼロであると仮定する）ことによって、劣決定系線型方程式を精度よく求解する技術が提案された。この技術はスパース正則化と呼ばれ、例えば、画像のノイズ除去、医用画像の再構成等に有効な技術として注目されている。
スパース正則化のアルゴリズムは大きく２種類に分類することが出来る。「Greedy algorithm」と「Ｌ１正則化」である。Greedy algorithmで最も性能が良い方法が、非特許文献１に開示されたsubspace pursuit（以下ＳＰと略記する。）であり、Ｌ１正則化で最も性能が良い方法が非特許文献２に開示されたＦＩＳＴＡである。

ＦＩＳＴＡに比べて圧倒的に高速に精度の良い解を見つけることが出来るＳＰの処理の流れを、図３を用いて説明する。
ステップＳ３０１では、係数行列Ａ、定数ベクトルbの値を設定する。また、非ゼロ要
素の数Ｋを設定する。Ｋの値はユーザが決めるものであり、Ｋ≦Ｍ／２を満たす任意の値が設定される。
以下、変数ｘの要素番号の集合をΩと表記し、全体集合をΩ_ａｌｌ：＝｛１，２，…，Ｎ｝と表記する。

ステップＳ３０２では、変数ｘ^（０）を初期化する。また、非ゼロ要素候補の要素番号の集合Ω^（０）を空集合で初期化する。反復回数ｔ＝１に設定する。なお、カッコ書きの上付き添え字は反復回数（イテレーション回数）を表す。

ステップＳ３０３では、変数ｘ^{（ｔ−１）}を用いて、残差ｒ^{（ｔ−１）}：＝ｂ−Ａｘ^{（ｔ−１）}を計算する。
ステップＳ３０４では、疑似誤差ε^{（ｔ−１）}：＝Ａ^Ｔｒ^{（ｔ−１）}を計算し、疑似誤差の絶対値が大きい順にＫ個の要素を抽出する。抽出されたＫ個の要素の要素番号の集合
をΔΩ^（ｔ）とする。
ステップＳ３０５では、集合Ω_ｐ ^（ｔ）＝Ω^{（ｔ−１）}∪ΔΩ^（ｔ）を求め、射影作用素Ｐ^（ｔ）：＝Ω_ａｌｌ→Ω_ｐ ^（ｔ）を生成する。Ω_ｐ ^（ｔ）の要素数はＫ個以上Ｍ個以下である。例えば、Ω_ａｌｌ：＝｛１，２，…，５｝でΩ_ｐ ^（ｔ）：＝｛１，３，５｝のとき、射影作用素Ｐ^（ｔ）は以下の行列になる。

ステップＳ３０６において、次式で表される部分問題を構成し、求解する。

ただし、ｚ^（ｔ）：＝Ｐ^（ｔ）ｘ^（ｔ）は部分問題での変数であり、次元はＫ以上Ｍ以下である。
次元がＭのとき、式（３）は正則であり解が一意に決まるので、通常の連立方程式の求解法、例えばJacobi法、Gauss-Seidel法、SOR法のような直接法や共役勾配法のような反
復解法を用いる。一方、次元がＭ未満のときは優決定系なので、最小二乗法の問題と捉え、例えば次式を求解する：

ステップＳ３０７では、ステップＳ３０６で得られた解ｚ^（ｔ）を用い、次式でｘ^（ｔ）の値を計算する。

ステップＳ３０８では、ｘ^（ｔ）から絶対値が大きい順にＫ個の要素を抽出し、抽出されたＫ個の要素の要素番号の集合をΩ^（ｔ）とする。このΩ^（ｔ）は次回の反復計算のステップＳ３０５で利用される。
所定の収束条件を満足するまで、ステップＳ３０３〜Ｓ３０８を繰り返す。

Wei Dai and Olgica Milenkovic: "Subspace Pursuit for Compressive Sensing Signal Reconstruction," IEEE Transactions on Information Theory, VOL.55, NO.5, pp.2230-2249, MAY 2009. Amir Beck and Marc Teboulle : "A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems," SIAM Journal on Imaging Sciences, Vol.2, No.1, pp.183-202, 2009.

従来の２つの求解法ＳＰ及びＦＩＳＴＡの性能は、図５に示す完全再構成マップで見ることが出来る。推定する信号の数Ｎ、即ち変数の数Ｎは２５６である。図５は、正規分布に従って生成した行列Ａ及び正解１０ケースに対して各求解法を用い、９ケースが相対誤差５％以下であったときを完全再構成領域として作図した。図中横軸は疎性率（＝非ゼロ要素の数／推定信号の数）、縦軸は観測率（＝観測信号の数／推定信号の数）である。グラフ中の各線はそれぞれの解法に対するものであり、線の左上の領域が完全再構成領域である。線５０１は理論限界、線５０２はＳＰの限界、線５０３はＦＩＳＴＡの限界である。同図より、いずれの方法も理論限界に対して改善の余地があることがわかる。例えば、Ｘ線ＣＴやＭＲＩ等の医療機器の場合、被験者の身体的・心理的負担を減らすために、観測回数をより少なくし測定時間をより短くすることが望まれる。よって、必要な観測信号の数を削減することは、実用上、きわめて重要な課題である。

本発明は、上記課題に鑑みなされたものであって、より小さい観測率で（より少ない観測信号数で）より高い推定精度を実現できる求解技術を提供することを目的とする。

本発明の第一態様は、求めるべき変数の要素数がＮであり線型方程式の数がＭである（Ｍ及びＮは１≦Ｍ＜Ｎを満たす整数）、劣決定系の連立線型方程式を反復計算により解くための求解装置であって、反復計算を行う計算手段を備え、前記計算手段は、各反復において、前記変数のＮ個の要素のなかから非ゼロ要素候補をＭ個以上抽出する抽出ステップと、抽出した非ゼロ要素候補のみを変数とする部分問題を設定し、前記部分問題を不良設定問題として求解する部分問題求解ステップと、前記部分問題求解ステップで得られた非
ゼロ要素候補の値に基づき、前記求めるべき変数のＮ個の要素の値を更新する更新ステップと、を実行することを特徴とする劣決定系線型方程式の求解装置である。

本発明の第二態様は、観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する被検体情報取得装置であって、前記観測系の特性を表すデータ、及び、前記観測系によって得られた観測データを取得するデータ取得手段と、本発明に係る劣決定系線型方程式の求解装置と、を備え、前記求解装置は、前記観測系の特性を表すデータを行列Ａ、前記観測系によって得られた観測データを定数ベクトルｂ、被検体の内部の情報を表す画像データを未知変数ｘとする、連立線型方程式Ａｘ＝ｂの解を計算することを特徴とする被検体情報取得装置である。

本発明の第三態様は、コンピュータを、本発明に係る劣決定系線型方程式の求解装置として機能させるためのプログラムである。

本発明によれば、より小さい観測率で（より少ない観測信号数で）より高い推定精度を実現できる求解技術を提供することができる。

本発明の実施形態に係る求解方法の処理を示すフローチャート。 本発明の実施形態に係る求解装置の構成を示すブロック図。 ＳＰの計算処理例を示すフローチャート。 ＦＩＳＴＡの計算処理例を示すフローチャート。 従来技術と本発明の実施形態の求解方法の性能を示す完全再構成マップ。 本発明の求解方法を適用した被検体情報取得装置の構成例。 従来技術（ＳＰ法）と本発明の実施形態の方法の誤差ノルムの収束の様子。

本発明は、劣決定系の連立線型方程式Ａｘ＝ｂをコンピュータによる反復計算により解くための求解技術に関し、前述したＳＰ（subspace pursuit）の改良アルゴリズムを提案するものである。

本方法と、従来のＳＰとの違いは、次のとおりである。ＳＰでは非ゼロ要素の値を求めるためにＭ次元以下の優決定系の部分問題を設定したのに対し、本方法ではＭ次元以上の劣決定系の部分問題（不良設定問題）を設定する点である。ＳＰにおいては、Ｋの値が推定精度に影響を及ぼすことが知られており、解くべき問題に応じて適切なＫの値を設定しなければならないことが実用上の課題の一つであった。この点、本方法では非ゼロ要素数を設定する必要がないので、ユーザの負担を減らせると共に、安定した推定精度が期待できる。

本方法において、部分問題を求解するための方法には、各種の正則化アルゴリズムを利用できる。好ましくは、部分問題の解法にスパース正則化を用いるとよい。前述のようにスパース正則化のアルゴリズムには大きく分けてGreedy algorithm（ＧＡ）とＬ１正則化とがあるが、ＧＡはＳＰと同様の課題を持つ場合があるので、Ｌ１正則化が好ましく利用できる。

また、本方法と従来のＬ１正則化との違いは次のとおりである。従来のＬ１正則化は元のＮ次元の問題を求解するが、本発明はＭ次元以上で２Ｍ次元以下の部分問題を求解するので、計算時間の点で劇的に有利である。

＜実施形態＞
以下、本発明に係る求解方法および求解装置の好ましい実施形態について、図面を参照して詳しく説明する。

（求解装置の構成）
図２は、本発明の実施形態に係る求解装置として機能するコンピュータの基本構成を示すブロック図である。このコンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２、入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩＣ２０６等から構成される。

ＣＰＵ（中央演算装置）２０１は、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶ディスプレイなどにより構成されており、ＣＰＵ２０１による処理結果（得られた解、再構成した画像）を文字や画像などでもって表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は，ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハード
ディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムがＲＡＭ２０４にロードされ、ＣＰＵ２０１によって実行されることで、本コンピュータが連立線型方程式を数値的に解くための求解装置として機能する。ＮＩＣ（ネットワークインターフェースコントローラ）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものである。本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。

なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。すなわち、求解装置は、クライアント−サーバ、クラウドコンピューティング、グリッドコンピューティングなどのシステム構成を採ることもできる。あるいは、医療機器や画像分析機器などに搭載される組み込み用コンピュータに求解装置の機能を実装することもできる。あるいは、求解装置の機能の全部又は一部をＡＳＩＣやＦＰＧＡなどの回路で構成してもよい。

（求解方法）
次に、本実施形態に係る劣決定系線型方程式の求解方法について説明する。ここでは、解くべき連立線型方程式としてＡｘ＝ｂが与えられるものとする。Ａは係数行列であり、ｂが定数ベクトル（右辺ベクトル）であり、ｘが求めるべき変数ベクトルである。また、変数ベクトルの次元、つまり変数ｘの要素数はＮであり、定数ベクトルの次元、つまり方程式の数はＭであり、１≦Ｍ＜Ｎ（Ｍ，Ｎは整数）である。

図１は、求解装置によって実行される求解処理の流れを示すフローチャートである。図１に示す各処理の実行主体は、求解装置のＣＰＵ（計算手段）２０１である。

ステップＳ１０１において、ＣＰＵ２０１は、ユーザによって入力された問題、すなわち、係数行列Ａの値と定数ベクトルｂの値をそれぞれ設定する。
ステップＳ１０２において、ＣＰＵ２０１は、変数ｘのメモリ領域を確保し、初期値ｘ^（０）を設定する。また、反復回数ｔ＝１に設定する。

ステップＳ１０３において、ＣＰＵ２０１は、変数ｘ^{（ｔ−１）}を用いて、残差ｒ^{（ｔ−１）}：＝ｂ−Ａｘ^{（ｔ−１）}を計算する。
ステップＳ１０４において、ＣＰＵ２０１は、残差ｒ^{（ｔ−１）}を用いて疑似誤差ε^{（ｔ−１）}：＝Ａ^Ｔｒ^{（ｔ−１）}を計算する。そして、ＣＰＵ２０１は、疑似誤差の絶対値が大きい順にＭ個の要素を抽出し、抽出されたＭ個の要素の要素番号の集合をΔΩ^（ｔ）とする。

ステップＳ１０５において、ＣＰＵ２０１は、非ゼロ要素候補の集合Ω_ｐ ^（ｔ）＝Ω^{（ｔ−１）}∪ΔΩ^（ｔ）を求め、射影作用素Ｐ^（ｔ）：＝Ω_ａｌｌ→Ω_ｐ ^（ｔ）を生成する。Ω_ａｌｌは要素番号の全体集合である。ここで、Ω_ｐ ^（ｔ）の要素数はＭ個以上２Ｍ個以下である。例えば、Ω_ａｌｌ：＝｛１，２，…，５｝でΩ_ｐ ^（ｔ）：＝｛１，３，５｝のとき、射影作用素Ｐ^（ｔ）は以下の行列になる。

ただし、射影作用素は陽に確保する必要はない。したがって、プログラムの実装では、
非ゼロ要素候補の要素番号（つまりΩ_ｐ ^（ｔ）の要素）だけをメモリに記憶しておけばよい。以下に述べる部分問題を取り扱う場合でも、射影作用素Ｐ^（ｔ）の計算を行うよりも、非ゼロ要素候補の要素番号に対応する部分を抽出する操作を行う方が計算量的に有利である。

上記ステップＳ１０３〜Ｓ１０５の処理が、非ゼロ要素候補をＭ個以上抽出する抽出ステップに相当する。本実施形態の抽出ステップでは、連立線型方程式の残差ｒ^{（ｔ−１）}に基づいて非ゼロ要素候補の更新を行っている。詳しくは、残差ｒ^{（ｔ−１）}に基づきＮ個の要素のなかから誤差ε^{（ｔ−１）}が大きいと推定される要素を選択し（ΔΩ^（ｔ））、前回の反復にて抽出した非ゼロ要素候補（Ω^{（ｔ−１）}）に新たに選択した要素を追加する（Ω^{（ｔ−１）}∪ΔΩ^（ｔ））。かかる更新操作によって、非ゼロ要素の抽出および修正を適切に行うことが期待できる。

次に、ステップＳ１０６において、ＣＰＵ２０１は、非ゼロ要素候補のみを変数とする部分問題を設定する。非ゼロ要素候補のみの変数をｚ^（ｔ）とすると、部分問題は次式のように表せる。

ただし、ｚ^（ｔ）：＝Ｐ^（ｔ）ｘ^（ｔ）である。変数ｚ^（ｔ）の要素数はＭ個以上２Ｍ個以下である。方程式はＭ個なので、要素数がＭ個の場合以外は不良設定問題になる。よって、式（７）の部分問題は、なんらかの正則化法を用いて求解する。

ステップＳ１０７において、ＣＰＵ２０１は、式（７）の部分問題をスパース正則化法を用いて解くことで、変数ｚ^（ｋ）の値を求める。本実施形態では、スパース正則化法の一つであるＦＩＳＴＡを部分問題の解法に適用する。

ＦＩＳＴＡでは、以下の式（８）で定義される評価関数Ｌ_２の最小値を探索することによって、式（７）を満足し、かつ非ゼロ要素数が少ない解ベクトルｚを得る。

評価関数の第１項は式（７）の連立方程式の残差ベクトルのＬ２ノルムであり、残差ノルムと呼ばれる。また、第２項は解ベクトルｚのＬ１ノルム（ベクトルｚの全要素の絶対値和）に定数λを乗じたものである。この評価関数Ｌ_１は、解ベクトルｚが正解に近づくほど第１項が小さくなり、解ベクトルｚがスパースになるほど（ゼロ要素が増すほど）第２項が小さくなる、という性質をもつ。

ＦＩＳＴＡの処理の流れを、図４を用いて説明する。
ステップＳ４０１では、ＣＰＵ２０１が、変数ｚ^（０）とｙ^（１）に初期値を設定し、反復回数ｋ＝１、μ^（１）＝１とする。
ステップＳ４０２では、ＣＰＵ２０１が変数ｚのステップ幅βを決定する。ステップ幅βには、次式で定義される∇Ｌ_２のLipschitz定数の逆数を設定するが、次元が大きくな
ると計算時間が膨大になるので、近似値で代替してもよい。ただし、近似値はLipschitz
定数の逆数を超えてはならない。

ステップ幅βは黄金分割探索（golden section search）によって決定することもでき
る。また、計算可能であれば次式によりステップ幅βを求めることも有効である：

ステップＳ４０３では、ＣＰＵ２０１は次式で変数ｚ^（ｋ）を計算する。

ここで、ｐｒｏｘ_γは下記式で定義される関数である。

γは反復回数に対して単調減少するように、例えば以下のように設定する。

ただし、ｋ_ＭＡＸはユーザが設定した最大反復回数である。

ステップＳ４０４では、ＣＰＵ２０１は次式でμ^{（ｋ＋１）}を計算する。

ステップＳ４０５では、ＣＰＵ２０１は次式でｙ^{（ｋ＋１）}を計算する。

上記更新式のかわりに、ＩＳＴＡに従って次式を用いてもよい：

ＣＰＵ２０１は所定の収束条件を満足するまで、ステップＳ４０３〜Ｓ４０５を繰り返す（ステップＳ４０６）。収束条件としては、例えば、評価関数Ｌ_２の値が所定の閾値よりも小さくなること、式（７）の残差ノルム（評価関数Ｌ_２の第１項に相当）の値が所定の閾値よりも小さくなること、反復回数ｋが所定数に達すること、などの条件を用いることができる。あるいは、評価関数、残差ノルム、反復回数のうちの二つ以上を組み合わせた収束条件を設定してもよい。組み合わせとは、複数の条件（命題）の論理積又は論理和である。ただし、部分問題を残差ノルムが小さくなるまで求解すると、全体問題の求解が収束しないことが分かっている。部分問題求解の反復回数は、例えば１０回程度で十分である。

ステップＳ４０６において所定の収束条件を満たすと判断されたら、式（７）の部分問題の解ｚ^（ｔ）を確定し、図１のステップＳ１０８に進む。上記ステップＳ１０６〜Ｓ１０７の処理が、部分問題を不良設定問題として求解する部分問題求解ステップに相当する。

ステップＳ１０８では、ＣＰＵ２０１が、ステップＳ１０７で得られた部分問題の解ｚ^（ｔ）（つまり、非ゼロ要素候補の値）を全体問題の解ｘに展開する。具体的には、次式の計算により変数ｘ^（ｔ）のＮ個の要素の値を更新する。

ステップＳ１０９では、ＣＰＵ２０１がｘ^（ｔ）から絶対値が大きい順にＭ個の要素を抽出し、抽出されたＭ個の要素の要素番号の集合をΩ^（ｔ）とする。このΩ^（ｔ）は次回の反復計算における非ゼロ要素候補の抽出ステップ（Ｓ１０３〜Ｓ１０５）で利用される。

ステップＳ１１０は、反復計算を終了するか否かを判断するステップである。ＣＰＵ２０１は、反復計算の状態が所定の収束条件を満足するか否かを評価し、収束条件を満たさない場合は反復回数ｔをインクリメントしてステップＳ１０３に戻り、収束条件を満たす場合は反復計算を抜けてステップＳ１１１に進む。収束条件としては、例えば、所定の評価関数の値が所定の閾値よりも小さくなること、残差ノルム（残差ベクトルｂ−ＡｘのＬ２ノルム）の値が所定の閾値よりも小さくなること、反復回数ｔが所定数に達すること、などの条件を用いることができる。評価関数としては、求めるべき変数ｘの関数として定義されるものであればいかなる関数を用いてもよい。あるいは、単独の条件（命題）ではなく、評価関数、残差ノルム、反復回数のうちの二つ以上を組み合わせた収束条件を設定してもよい。組み合わせとは、複数の条件（命題）の論理積又は論理和である。例えば、「反復回数が所定数を超え、且つ、評価関数の値が閾値よりも小さくなること」を収束条件としたり、「評価関数の値が第１の閾値よりも小さくなるか、又は、残差ノルムの値が第２の閾値よりも小さくなること」を収束条件としてもよい。あるいは、評価関数や残差ノルムの値が閾値を下回る状態が複数回続くこと、を収束条件としてもよい。

ステップＳ１１１は、反復計算終了後に、連立線型方程式の解を決定するステップである。簡単には、最後の反復計算で得られた変数ｘの値を方程式の解に選べばよい。あるいは、ＣＰＵ２０１が、ステップＳ１０３〜Ｓ１０９の反復計算の中で、評価関数の値又は
残差ノルムを最も小さくする変数ｘをメモリに記憶しておき、そのときの変数ｘの値を最適解として出力してもよい。

以上述べた処理によって、劣決定系の連立線型方程式の解を精度よく求めることができる。本実施形態の求解方法の性能を図５の完全再構成マップの線５０４に示す。従来方法のＳＰ（５０２）及びＦＩＳＴＡ（５０３）に比べて、本実施形態の方法の性能（５０４）は理論限界（５０１）に近づいている。特に疎性率が０．２以上の広い領域において、従来方法では不可能だった問題の再構成が可能になっていることが分かる。

次にノイズ除去の問題に適用した結果を図７に示す。対象画像は６４×６４画素であり、方程式数、即ち観測回数は画素数と同じであるが、ノイズが重畳している。図中、横軸が計算時間、縦軸が正解と近似解との誤差ノルムである。また、ＳＰ法の結果７０１、７０２、７０３に対して本実施形態の方法の結果は７０４、７０５、７０６であり、それぞれ計算パラメタを１０２４、１６００、２０４８に設定したものである。正解の非零要素数は８００であり、ＳＰ法も本実施形態の方法もパラメタ値１０２４（７０１と７０４）および１６００（７０２と７０５）では正解に収束した。しかし、パラメタ値が正解の非零要素数と大きく異なる場合、即ち２０４８のとき、本実施形態の方法は正解に収束する（７０６）が、ＳＰ法は収束しなかった（７０３）。このように、本実施形態の方法は、ＳＰ法の計算パラメタ、即ちユーザが設定すべき非零要素数に対して頑健であることがわかる。また、本実施形態の方法では部分問題を完全に求解しないので、ＳＰ法に比べて計算時間が短縮できる。これは、７０１と７０４、７０２と７０５の比較によって確認することが出来る。

（適用例）
本実施形態の求解方法を被検体情報取得装置における画像再構成処理に適用した例について説明する。被検体情報取得装置とは、観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する装置であり、例えば、医用画像診断装置（Ｘ線ＣＴ、ＭＲＩ、光音響トモグラフィ、超音波診断装置など）や非破壊検査装置が該当する。

図６は、本実施形態の求解装置を搭載した医用画像診断装置の構成を模式的に示す図である。図に示すように、医用画像診断装置は、観測系（測定装置）６０１、投影データ生成部６０２、データ取得部６０３、再構成部６０４、表示装置６０５などを備える。

観測系６０１は、測定装置筐体６０６、信号発生源６０８、６０９、６１０、検出器６１１、６１２、６１３を有している。信号発生源６０８、６０９、６１０は、筐体６０６の内部に置かれた被検体６０７に対し、測定用信号を照射する装置である。例えば、Ｘ線ＣＴであればＸ線源、ＭＲＩであれば磁界発生源、光音響トモグラフィであればレーザ光源、超音波診断装置であれば超音波送信機が該当する。被検体６０７を透過した信号、又は被検体６０７の内部で発生もしくは反射した信号が、検出器６１１、６１２、６１３で検出される。

投影データ生成部６０２は、複数の検出器６１１、６１２、６１３で検出された信号を所定のフォーマットに従って一つの観測データにまとめる。例えば２次元画像のラスタースキャンによるベクトル表現のようにである。これによって、それぞれの検出信号ｂ_１、ｂ_２、ｂ_３が一つのベクトルｂとして統合される。

データ取得部６０３は、投影データ生成部６０２から取得した観測データ（ベクトルｂ）と、予め記憶されている観測系６０１の特性を表すデータ（係数行例Ａ）とから、連立方程式Ａｘ＝ｂを作成する。ｘは未知変数であり、Ｎ次元のベクトルである。係数行列Ａ
は、観測系６０１を構成する信号発生源や検出器の物理的特性などを元に理論的に作成してもよいし、内部構造が既知のサンプルを観測系６０１で実測した結果から作成してもよい。図６のＡ_１、Ａ_２、Ａ_３はそれぞれ、信号発生源６０８と検出器６１１の特性、信号発生源６０９と検出器６１２の特性、信号発生源６１０と検出器６１３の特性を示している。

再構成部６０４は、図１に示した求解方法を用いて未知変数ｘを求める。そして、再構成部６０４は、得られた未知変数ｘの解に基づき、被検体６０７の内部構造を表す画像データを生成する。生成されたデータは表示装置６０５に表示される。このとき、被検体の断面（２次元画像）を表示してもよいし、被検体内部の３次元構造を表示してもよい。その際、ユーザーインターフェイスによって、画像の回転・並進などの操作を指示できるとよい。

ところで、本実施形態の求解方法は、未知変数ｘのうちの多くがゼロ要素とみなせるとの仮定をおいている。それゆえ、未知変数ｘの疎性が十分でない場合には、解の推定精度が低下するなど、期待する性能が得られない可能性がある。そこで、データ取得部６０３が、係数行列Ａの代わりに係数行列（ＡΦ^−１）を用いて、連立方程式（ＡΦ^−１）ｗ＝ｂを作成し、再構成部６０４が図１に示した求解方法を用いて未知変数ｗを求めた後、ｘ＝Φ^−１ｗにより本来の未知変数ｘを計算してもよい。Φは、未知変数ｘの疎性を高める作用をもつ特徴量変換行列である。つまり、未知変数ｘに変換行列Φを作用させた変数であるｗ＝Φｘは、変数ｘに比べて高い疎性をもつので、本実施形態の方法による求解が容易となるのである。変換行列Φとしては、被検体の内部の情報を表す画像データ（変数ｘ）のゼロ要素を増加させるフィルタ、例えば、微分フィルタ、２次微分フィルタ、ハイパスフィルタ、カットオフフィルタなどを用いることができる。このような方法により、画像再構成の精度をより高めることが期待できる。

２０１…ＣＰＵ、２０２…表示装置、２０３…入力装置、２０４…ＲＡＭ、２０５…ハードディスク装置、２０６…ＮＩＣ、２０７…バス

Claims (10)

1. 求めるべき変数の要素数がＮであり線型方程式の数がＭである（Ｍ及びＮは１≦Ｍ＜Ｎを満たす整数）、劣決定系の連立線型方程式を反復計算により解くための求解装置であって、
   反復計算を行う計算手段を備え、
   前記計算手段は、各反復において、
   前記変数のＮ個の要素のなかから非ゼロ要素候補をＭ個以上抽出する抽出ステップと、
   抽出した非ゼロ要素候補のみを変数とする部分問題を設定し、前記部分問題を不良設定問題として求解する部分問題求解ステップと、
   前記部分問題求解ステップで得られた非ゼロ要素候補の値に基づき、前記求めるべき変数のＮ個の要素の値を更新する更新ステップと、を実行することを特徴とする劣決定系線型方程式の求解装置。
2. 前記計算手段は、前記部分問題求解ステップでは、スパース正則化法を用いて前記部分問題を求解することを特徴とする請求項１に記載の劣決定系線型方程式の求解装置。
3. 前記計算手段は、前記抽出ステップでは、前記連立線型方程式に前記変数の値を代入したときの左辺と右辺の差である残差に基づいて、抽出する非ゼロ要素候補を決定することを特徴とする請求項１又は２に記載の劣決定系線型方程式の求解装置。
4. 前記計算手段は、前記抽出ステップでは、前記残差に基づき前記Ｎ個の要素のなかから誤差の絶対値が大きい順に所定個の要素を選択し、前回の反復において抽出した非ゼロ要素候補に対し前記選択した要素を追加することを特徴とする請求項３に記載の劣決定系線型方程式の求解装置。
5. 前記計算手段は、各反復において、前記求めるべき変数の関数として定義される評価関数の値、前記連立線型方程式に前記変数の値を代入したときの左辺と右辺の差で定義されるベクトルのノルムである残差ノルム、反復計算の回数、又は、これらのうちの二つ以上の組み合わせにより構成される収束条件を満足するか否かで、反復計算の終了を判断する
   ステップ、を実行することを特徴とする請求項１〜４のうちいずれか１項に記載の劣決定系線型方程式の求解装置。
6. 前記計算手段は、反復計算の終了後、各反復で得られた前記変数の値のうち、前記評価関数の値又は前記残差ノルムを最も小さくするものを、前記連立線型方程式の解に決定するステップ、を実行することを特徴とする請求項５に記載の劣決定系線型方程式の求解装置。
7. 観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する被検体情報取得装置であって、
   前記観測系の特性を表すデータ、及び、前記観測系によって得られた観測データを取得するデータ取得手段と、
   請求項１〜６のいずれか１項に記載の劣決定系線型方程式の求解装置と、を備え、
   前記求解装置は、前記観測系の特性を表すデータを行列Ａ、前記観測系によって得られた観測データを定数ベクトルｂ、被検体の内部の情報を表す画像データを未知変数ｘとする、連立線型方程式Ａｘ＝ｂの解を計算することを特徴とする被検体情報取得装置。
8. 前記求解装置は、
   予め与えられた変換行列Φを用いて連立線型方程式（ＡΦ^−１）ｗ＝ｂにおける変数ｗを求解した後、
   ｘ＝Φ^−１ｗにより、変数ｘを求めることを特徴とする請求項７に記載の被検体情報取得装置。
9. 前記変換行列Φは、被検体の内部の情報を表す画像データのゼロ要素を増加させるフィルタであることを特徴とする請求項８に記載の被検体情報取得装置。
10. コンピュータを、請求項１〜６のうちいずれか１項に記載の劣決定系線型方程式の求解装置として機能させるためのプログラム。

Images