【文献】
三上 拓也、浜本 隆之 ,“近赤外動画像を利用した低照度下におけるカラー動画像の取得”,映像メディア処理シンポジウム 第17回シンポジウム資料,日本,電子情報通信学会画像工学研究専門委員会,2012年10月24日,pp.55-56
【文献】
Dong-Hoon Lee et al.,"Sparse Sampling MR Image Reconstruction Using Bregman Iteration: A Feasibility Study at Low Tesla MRI System",2011 IEEE Nuclear Science Symposium Conference Record,米国,IEEE,2012年02月21日,pp.3446-3449
(58)【調査した分野】(Int.Cl.,DB名)
【発明を実施するための形態】
【0010】
以下、図面を参照して実施形態について説明する。なお、以下の説明において、同一の機能及び構成を有する構成要素については、共通する参照符号を付す。
【0011】
1. 実施形態
本実施形態に係る画像データ復元システムについて説明する。本実施形態に係る画像データ復元システムは、マルチスペクトルカメラにおける画像のノイズ除去、スペクトラルX線装置における画像の再構成、及びMRI装置における画像の再構成等を行うシステムに広く適用可能である。
【0012】
1.1 構成について
1.1.1 画像データ復元システムの構成について
まず、本実施形態に係る画像データ復元システムの構成について
図1を用いて説明する。
【0013】
図1は、本実施形態に係る画像データ復元システム1の構成の一例を示すブロック図である。
図1に示すように、画像データ復元システム1は、センサ2と、画像データ復元装置3と、を備えている。
【0014】
センサ2は、図示しない観測対象に対する観測を行う。センサ2は、画像データ復元装置3に対して、当該観測の結果として得られた信号を観測データyとして出力する機能を備えている。
【0015】
センサ2は、複数のチャネルに対応する観測データyを出力し得る。チャネルとは、例えば、センサ2が複数のセンサの集合として構成される場合における、当該複数のセンサの各々を意味し得る。より具体的には、例えば、センサ2がMRI装置に実装された複数のコイルに相当する場合、チャネルとは、当該複数のコイルの各々に相当する。また、チャネルとは、例えば、センサ2が複数のスペクトルにおいて観測を行う場合における、当該複数のスペクトルの各々を意味し得る。より具体的には、例えば、センサ2がマルチスペクトルカメラに実装されたセンサに相当する場合、チャネルとは、当該センサにおいて観測される複数のスペクトルの各々に相当する。なお、センサ2は、同一の観測対象について、複数のチャネルに対応する観測データyを、複数回にわたって出力してもよい。
【0016】
画像データ復元装置3は、観測データ取得部11と、記憶部12と、画像データ復元部13と、復元画像出力部14と、を備えている。観測データ取得部11、記憶部12、画像データ復元部13、及び復元画像出力部14の各々は、例えば、以下に個別に示される機能構成を実現するための回路によって構成される。観測データ取得部11、記憶部12、画像データ復元部13、及び復元画像出力部14の各々に対応する回路は、例えば、図示しないプロセッサ等によって制御されることによって協働して動作し、後述の画像データ復元動作を実現可能に構成される。
【0017】
観測データ取得部11は、センサ2で観測された複数のチャネルに対応する観測データyを取得する。観測データ取得部11は、例えば、観測データyを記憶部12に記憶させる。
【0018】
記憶部12は、例えば、HDD(Hard Disk Drive)、SSD(Solid State Drive)等の記憶媒体を含む。記憶部12は、例えば、観測データ取得部11及び画像データ復元部13から読出し/書込み可能な記憶媒体であり、観測データyを一時的に記憶する。記憶部12は、同一の観測対象を対象として観測された複数個の観測データyが存在する場合、当該複数個の観測データyを、1つの観測データフルセットy
fullとして集約して記憶する。観測データフルセットy
fullは、画像データ復元部13において使用され得る全てのチャネルに対応する観測データyを含む。観測データフルセットy
fullは、例えば、列ベクトルであり、複数のチャネルをまとめて1列にしたベクトルである。
【0019】
画像データ復元部13は、記憶部12から観測データフルセットy
fullを読出し、当該観測データフルセットy
fullに基づいて、観測対象に係る画像データ(例えば中間画像データ)を生成する。画像データとは、観測対象についての観測データyに何らかのマッピングを施すことにより、2次元又は3次元で可視化することが可能な、2次元以上のデータである。画像データ復元部13は、画像データを最適化することによって復元画像データxを復元し、復元画像出力部14に出力する。画像データ復元部13における復元画像データxを復元する処理の詳細については後述する。
【0020】
復元画像出力部14は、例えば、画像データ復元部13から復元画像データxを受けると、当該復元画像データxに基づいて、観測対象の復元画像を出力する。復元画像出力部14は、例えば、表示画面(例えば、LCD(Liquid Crystal Display)又はEL(Electroluminescence)ディスプレイ等)等を含み、当該表示画面を介してユーザに観測対象に係る画像を視覚的に提供する。なお、復元画像出力部14は、上述の例に限らず、画像データを可視化可能な機能構成を有する図示しない外部機器に対して、復元画像データxを出力してもよい。
【0021】
1.1.2 画像データ復元部の構成について
次に、本実施形態に係る画像データ復元部の構成について
図2を用いて説明する。
【0022】
図2は、本実施形態に係る画像データ復元部13の構成の一例を示すブロック図である。
図2に示すように、画像データ復元部13は、チャネル選択部131と、復元モデル132と、事前知識モデル133と、最適化演算部134と、判定部135と、を備えている。
【0023】
チャネル選択部131は、記憶部12から読み出された観測データフルセットy
fullに対応する複数のチャネルのうち、最適化演算部134に入力される観測データに対応するチャネルを選択する。チャネル選択部131は、全てのチャネルを選択する場合、観測データフルセットy
fullを最適化演算部134に出力する。また、チャネル選択部131は、全てのチャネルのうちの一部を選択する場合、当該一部のチャネルに対応する観測データを観測データサブセットy
subとして抽出し、最適化演算部134に出力する。観測データサブセットy
subは、例えば、観測データフルセットy
fullを構成する列ベクトルのうち、チャネルサブセットに対応するデータを含む新たな列ベクトルとして抽出される。
【0024】
以下の説明では、観測データフルセットy
fullに対応する複数のチャネルを「チャネルフルセット」とも言い、観測データサブセットy
subに対応する1又は複数のチャネルを「チャネルサブセット」とも言う。
【0025】
また、チャネル選択部131は、最適化演算部134において繰り返し実行される最適化演算を打切るか否かを判定した結果を含む判定結果を、判定部135から受信する。チャネル選択部131は、例えば、当該判定結果に応じて、チャネルフルセットを選択するか、又はチャネルサブセットを選択するかを決定する。具体的には、チャネル選択部131は、最適化演算を打切る旨の判定結果を受けると、チャネルフルセットを選択し、最適化演算部134に対して観測データフルセットy
fullを出力する。また、チャネル選択部131は、最適化演算を続行する(打切らない)旨の判定結果を受けると、チャネルサブセットを選択し、最適化演算部134に対して観測データサブセットy
subを出力する。
【0026】
なお、チャネル選択部131がチャネルフルセットの中からいずれのチャネルをチャネルサブセットとして選択するかは、種々の手法が適用可能である。
【0027】
例えば、チャネル選択部131は、チャネルサブセットをランダムに選択してもよい。この場合、チャネルサブセットは、最適化演算部134に観測データサブセットy
subが出力されるたびにランダムに選択されてもよいし、一旦チャネルサブセットが選択された後は固定されてもよい。このような選択方法は、チャネルサブセットの選択に要する負荷を抑えることができる。
【0028】
また、例えば、チャネル選択部131は、予め定められたチャネルの組をチャネルサブセットとして選択してもよい。この場合、当該チャネルの組は、例えば、事前学習結果に基づき、最も復元画像データへの寄与率が高いと判定されたチャネルの組が選択される。このような選択方法は、最適化演算部134の最適化演算を開始する前に事前に行うことができるため、最適化演算に要する時間を短縮することができる。また、このような選択方法は、チャネルサブセットを事前に選択しておくという性質上、事前学習時と、最適化演算の実行時と、でチャネルの性質が変化しない場合に特に有効である。
【0029】
また、例えば、チャネル選択部131は、観測データフルセットy
fullに基づいて、チャネルサブセットを選択してもよい。具体的には、例えば、チャネル選択部131は、観測データフルセットy
fullのSNR(Signal to Noise Ratio)を演算し、SNRの高いチャネルをチャネルサブセットとして選択してもよい。このような選択方法は、最適化演算部134の最適化演算に際して新たな演算を行う必要が生じるため、最適化演算に要する時間が上述の例と比較して増加する。しかしながら、チャネルサブセットを実際に観測されたデータに応じて選択できるため、事前学習時と、最適化演算の実行時と、でチャネルの性質が変化し得る場合においても有効である。
【0030】
復元モデル132及び事前知識モデル133は、最適化演算部134において実行される最適化演算に用いられる数学モデルである。復元モデル132及び事前知識モデル133は、例えば、記憶部12内にデータベースとして記憶され、最適化演算部134によって適宜読み出される。
【0031】
復元モデル132は、真の画像データ(「x
true」と呼ぶ)で表現される観測対象について、誤差のない真の観測データ(「y
true」と呼ぶ)が観測された場合に、当該真の観測データy
trueと、真の画像データx
trueとを相互に変換可能な線形変換である。復元モデル132は、例えば、線形変換Aを含む。すなわち、線形変換Aは、y
true−Ax
true=0を満たす。
【0032】
なお、センサ2によって実際に観測された観測データyには、観測誤差が含まれる。また、復元画像データxは、真の画像データx
trueに必ずしも一致しない。このため、演算y−Axは、必ずしも“0”とはならないが、“0”に近づくほど、復元画像データxが観測データyに対して適合していると考えられる。以下の説明では、復元画像データxが観測データyに対して適合しているか否かを示す指標を「データ適合度」と言い、E
d(x)と表す。データ適合度E
d(x)は、例えば、その値が小さいほど、復元画像データxが観測データyに対して適合していることを示す。
【0033】
事前知識モデル133は、真の画像データx
trueに関する既知の情報をモデル化したものである。例えば、真の画像データx
trueの画素は、完全にランダムに分布しているわけではなく、特定の特徴を有するパターンにしたがって分布している場合があることが経験的に知られている。具体的には、例えば、隣接画素同士の輝度の変化率に着目すると、真の画像データx
trueは、平坦な場合と部分のはっきりした部分とが比較的多く、ノイズ状の部分は全体的に少ない、等の特徴を有する場合がある。事前知識モデル133は、学習データを用いることによって上述のような真の画像データx
trueに関する特徴を事前に学習し、復元画像データxの関数P(x)としてモデル化したものである。関数P(x)は、実数の関数である。
【0034】
なお、関数P(x)は、例えば、線形変換Tを用いて、P(x)=Txのように更にモデル化されることができる。以下の説明では、復元画像データxが事前知識に対して適合しているか否かを示す指標を「事前知識適合度」と言い、E
p(x)と表す。事前知識適合度E
p(x)は、例えば、その値が小さいほど、復元画像データxが事前知識に対して適合していることを示す。
【0035】
最適化演算部134は、チャネル選択部131から観測データフルセットy
full又は観測データサブセットy
subを受信すると、記憶部12から復元モデル132及び事前知識モデル133を読み出す。最適化演算部134は、復元モデル132及び事前知識モデル133に基づいて、データ適合度E
d(x)及び事前知識適合度E
p(x)を最適化(最小化)するための最適化演算を繰り返し実行する。最適化演算部134は、当該繰り返し実行される最適化演算の過程において生成される中間画像データx
kを生成する。中間画像データx
kは、最適化演算がk回目(kは、自然数)に繰り返された際に生成される画像データを示す。最適化演算部134は、中間画像データx
kを判定部135に送信してもよい。
【0036】
すなわち、最適化演算部134は、k回目の最適化演算の後に判定部135から繰り返し演算要求を受けると、チャネル選択部131から観測データフルセットy
full又は観測データサブセットy
subのいずれかを受信する。そして、最適化演算部134は、(k+1)回目の最適化演算を実行して、中間画像データx
k+1を生成し、判定部135に送信する。なお、最適化演算部134において繰り返し実行される最適化演算の詳細については、後述する。
【0037】
判定部135は、中間画像データx
kを受けるたびに、最適化演算部134による最適化演算の繰り返しの打切り条件を満たすか否かを判定する。打ち切り条件は、例えば、中間画像データx
kが算出された際における評価関数E(x
k)が所定の閾値に達するか否かに基づいて設定される。この場合、判定部135は、中間画像データx
kについて評価関数E(x
k)を算出し、評価関数E(x
k)が所定の閾値を下回った場合、打切り条件を満たすと判定し、下回らなかった場合、打切り条件を満たさないと判定する。
【0038】
なお、打切り条件は、上述の例に限らず、最適化演算の繰り返し回数の上限値として設定されてもよい。この場合、判定部135は、最適化演算の繰り返し回数が所定の上限値に達した場合、打切り条件を満たすと判定し、所定の上限値に達していない場合、打切り条件を満たさないと判定する。なお、打ち切り条件が繰り返し回数の上限値として設定される場合、判定部135は、必ずしも評価関数E(x
k)を算出しなくてもよい。
【0039】
判定部135は、k回目の最適化処理の結果、打切り条件を満たさないと判定した場合、当該判定結果をチャネル選択部131に送信すると共に、(k+1)回目の最適化演算についての繰り返し演算要求を最適化演算部134に送信する。一方、判定部135は、k回目の最適化処理の結果、打切り条件を満たすと判定した場合、(k+1)回目の最適化演算が最後の最適化演算であると判定する。すなわち、判定部135は、打切り条件を満たす旨の判定結果をチャネル選択部131に送信すると共に、最後の最適化演算についての繰り返し演算要求を最適化演算部134に送信する。
【0040】
また、判定部135は、最適化演算部134から最後の最適化演算の結果として中間画像データx
kを受け取ると、当該中間画像データx
kを復元画像データxとして復元画像出力部14に出力する。
【0041】
1.1.3 最適化演算部における最適化演算について
次に、本実施形態に係る最適化演算部134において繰り返し実行される最適化演算の詳細について説明する。
【0042】
1.1.3.1 繰り返し解法について
最適化演算部134は、以下の式(1)に示す評価関数E(x)を最小化することにより、最終的に復元画像データxを算出する。すなわち、最適化演算部134における最適化演算の繰り返しは、以下の式(1)を最小化する最小化問題として定式化される。
【0043】
E(x)=E
d(x)+E
p(x) (1)
なお、画像復元問題がMAP(Maximum A Posteriori)推定における確率の最大化問題として定式化されている場合、最大化問題における評価関数は、その対数に“−1”を乗ずることによって式(1)に示す評価関数E(x)の形に変形できる。このため、いずれにしても、画像復元問題は、式(1)の最小化問題に帰着できる。
【0044】
また、評価関数E(x)の各項のうち、データ適合度E
d(x)は、以下の式(2)で定式化される。
【0046】
仮に事前知識適合度E
p(x)がデータ適合度E
d(x)と同様にL2ノルムで定式化される場合、式(1)に示される最小化問題は、評価関数E(x)の微分が“0”となるような復元画像データxを求める線形計画問題に帰着できる。このような線形計画問題は、例えば、共役勾配法等を用いることによって閉じた解を得ることができる。
【0047】
しかしながら、自然界に存在するものを観測対象とする場合、真の画像データx
trueは、スパース性を有することが事前知識として知られている場合が多い。この場合、上述の線形計画問題の解は、真の画像データx
trueを精度よく復元できない場合がある。したがって、復元精度を向上させるため、事前知識適合度E
p(x)には、スパース性を表現可能なL1ノルム又はトレースノルムの項を含むものが用いられることが多い。事前知識適合度E
p(x)は、例えば、L1ノルムを用いて、以下の式(3)のように定式化される。
【0049】
この場合、式(1)に示される最小化問題は、L2ノルムの場合と同様の手法で閉じた解を得ることは困難となる。一方、L1ノルムを含む最小化問題は、限定された条件を満たす形に変形することによって、閉じた解を得ることができることが知られている。
【0050】
具体的には、例えば、L1ノルムを含む最小化問題は、軟判定閾値法(Soft Thresholding)を適用することによって解くことができる。軟判定閾値法は、以下の式(4)で示される最小化問題に対して閉じた解を与える演算である。
【0052】
ここで、ベクトルp及びp0は列ベクトルである。設計変数αは、0以上の実数である。なお、式(4)の解を与えるベクトルpの各要素は、上付き添え字をベクトルの要素番号として、以下の式(5)で与えられる。
【0054】
そこで、一般的に、式(1)に示される最小化問題は、補助変数を導入することにより、軟判定閾値法等が適用可能な最適化問題を含む2つ以上の最小化問題に分割する手法が用いられる。
【0055】
最適化演算部134は、上述の2つ以上に分割された最小化問題の各々を1回ずつ解く処理を1回の最適化演算として実行する。そして、最適化演算部134は、当該2つ以上に分割された最小化問題を交互に解くことを繰り返すことによって、最終的に式(1)を最小化する復元画像データxを算出することができる。
【0056】
上述した繰り返し解法としては、例えば、ペナルティ法(Penalty Method)、ADMM法(Alternating Direction Method of Multipliers)、Split Bregman Method等、種々の解法が知られている。本実施形態における最適化演算部134は、上述の例に限らず、データ適合度E
d(x)の項を含む最適化問題を繰り返し解法によって解く任意の解法に対して適用可能である。
【0057】
1.1.3.2 ペナルティ法を適用した場合
以下では、最適化演算を繰り返し実行して復元画像データxを得るための繰り返し解法の一例として、最適化演算部134にペナルティ法を適用した場合について説明する。
【0058】
ペナルティ法は、内部データとして補助変数zを導入することによって、式(1)をL2ノルムにより構成される最小化問題と、L1ノルムを含む最小化問題と、の2つの最小化問題に分離する。具体的には、最適化演算部134は、(k+1)回目の最適化演算において、以下の式(6)〜(8)に示す最小化問題を解く。
【0060】
ここで、設計変数cは1以上の実数であり、設計変数δ
kは0以上の実数である。式(8)に示す通り、設計変数δ
kは、最適化処理の繰り返し数がインクリメントされるたびに小さな値から徐々に大きな値に増加していく。最適化処理の繰り返しによって設計変数δ
kが十分な大きさの値になると、式(6)〜(8)によって得られる解は、式(1)の解と等価となる。
【0061】
なお、設計変数c及びδ
kは、上述の例に限らず、種々の変形が可能である。例えば、設計変数cは、式(8)では、繰り返し数によらない定数として設定されているが、これに限らず、繰り返し数に応じて動的に設定されてもよい。また、設計変数δ
kは、設計変数cに対して、以下の式(9)を満たすように設定されてもよい。
【0063】
1.1.3.3 ADMM法を適用した場合
次に、ペナルティ法に代えてADMM法を最適化演算部134に適用した場合について説明する。
【0064】
ADMM法は、内部データとして補助変数z及びdを導入することによって、式(1)をL2ノルムにより構成される最小化問題と、L1ノルムを含む最小化問題と、の2つの最小化問題に分離する。具体的には、最適化演算部134は、(k+1)回目の最適化演算において、以下の式(10)〜(12)に示す最小化問題を解く。
【0066】
ここで、設計変数δは正の実数であり、設計変数γは実数である。
【0067】
なお、最適化演算部134において繰り返し実行される最適化演算にADMM法が適用される場合、ADMM法は、上述の式(10)〜(12)に限らず、より高速に演算可能な形に変形されてもよい。ADMM法の高速化手法については、例えば、以下に示す文献([AADMM])に示される手法が適用可能である。
【0068】
[AADMM] Mojtaba Kadkhodaie et al., “Accelerated Alternating Direction Method of Multipliers”, Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 497-506.
1.2 動作について
次に、本実施形態に係る画像データ復元装置における画像データ復元動作について説明する。
図3は、本実施形態に係る画像データ復元装置の動作の一例を示すフローチャートである。
【0069】
図3に示すように、ステップST10において、観測データ取得部11は、センサ2において観測された、複数のチャネルに対応する観測データyを取得する。記憶部12は、取得された観測データyを同一の観測対象について集約し、チャネルフルセットに対応する観測データフルセットy
fullとして記憶する。
【0070】
ステップST20において、チャネル選択部131は、チャネルフルセットのうち、チャネルサブセットとして使用するチャネルを選択する。
【0071】
以降のステップST30〜ステップST70において、最適化演算部134及び判定部135における最適化演算の繰り返し処理が実行される。以下の説明では、ステップST30〜ステップST50では、繰り返し処理が打ち切られていないk回目の最適化演算が実行され、ステップST60及びステップST60では、最後の最適化演算が実行されるものとして説明する。
【0072】
ステップST30において、チャネル選択部131は、k回目の最適化演算に際して使用する観測データyとして、観測データサブセットy
subを適用する。具体的には、チャネル選択部131は、ステップST20において選択されたチャネルサブセットに対応する観測データサブセットy
subを抽出し、最適化演算部134に送信する。
【0073】
ステップST40において、最適化演算部134は、観測データサブセットy
subを用いて、復元モデル132に基づくデータ適合度E
d(x
k)と、事前知識モデル133に基づく事前知識適合度E
p(x
k)と、を最適化するためのk回目の最適化演算を実行し、中間画像データx
kを生成する。最適化演算部134は、中間画像データx
kを判定部135に送信する。
【0074】
ステップST50において、判定部135は、k回目の最適化演算が終了した後、最適化演算の繰り返しについて、打切り条件を満たすか否かを判定する。
【0075】
打ち切り条件を満たさないと判定した場合(ステップST50;no)、判定部135は、当該判定結果をチャネル選択部131に送信すると共に、(k+1)回目の最適化演算についての繰り返し演算要求を最適化演算部134に送信し、ステップST30に戻る。ステップST30において、チャネル選択部131は、判定結果に基づき、(k+1)回目の最適化演算に観測データサブセットy
subを適用する。
【0076】
一方、打ち切り条件を満たすと判定した場合(ステップST50;yes)、判定部135は、当該判定結果をチャネル選択部131に送信すると共に、最後の最適化演算についての繰り返し演算要求を最適化演算部134に送信し、ステップST60に進む。ステップST60において、チャネル選択部131は、判定結果に基づき、(k+1)回目の最適化演算に観測データフルセットy
fullを適用する。具体的には、チャネル選択部131は、チャネルフルセットに対応する観測データフルセットy
fullを抽出し、最適化演算部134に送信する。
【0077】
ステップST70において、最適化演算部134は、観測データフルセットy
fullを用いて、復元モデル132に基づくデータ適合度と、事前知識モデル133に基づく事前知識適合度と、を最適化するための(k+1)回目の最適化演算を実行し、中間画像データx
k+1を生成する。最適化演算部134は、中間画像データx
k+1を判定部135に送信する。判定部135は、中間画像データx
k+1を復元画像データxとして復元画像出力部14に送信する。
【0078】
ステップST80において、復元画像出力部14は、復元画像データxを用いて復元画像をユーザに出力する。
【0079】
以上で、画像データ復元動作が終了する。
【0080】
1.3 本実施形態に係る効果
本実施形態によれば、画質の劣化を抑制しつつ、画像データの復元に要する時間を短縮することが出来る。本効果につき、以下説明する。
【0081】
本実施形態に係る画像データ復元装置3の最適化演算部134は、データ適合度E
d(x)及び事前知識適合度E
p(x)を最適化するための最適化演算を繰り返し実行する。そして、判定部135は、当該最適化演算の繰り返し処理が実行されるたびに、繰り返しを打切るか否かを判定し、その判定結果をチャネル選択部131に送信する。チャネル選択部131は、当該判定結果に応じて、最後の最適化演算に対してはチャネルフルセットに対応する観測データフルセットy
fullを適用し、それ以外の最適化演算に対してはチャネルサブセットに対応する観測データサブセットy
subを適用する。
【0082】
これにより、最適化演算が繰り返される間は、観測された情報のうちの一部を用いることで、最適化演算の処理量を低減することができ、ひいては処理速度を早めることができる。また、最後の最適化演算では観測された全ての情報を用いることにより、最適化演算の繰り返しの途中で一部の情報しか用いないことによる画質の劣化を抑制することができる。
【0083】
また、最適化演算部134は、内部データとして補助変数zを導入することによって、ペナルティ法又はADMM法を実行する。これにより、最適化演算部134は、復元画像データxを最適化する演算と、内部データを最適化する演算とを交互に実行することができる。このため、最適化演算部134は、当該交互に実行される演算を繰り返すことにより、最終的にデータ適合度E
d(x)及び事前知識適合度E
p(x)を最適化することができる。
【0084】
2. 変形例
なお、本実施形態に係る画像データ復元装置は、上述の例に限らず、種々の変形が適用可能である。
【0085】
2.1 第1変形例
上述の通り、データ適合度E
d(x)は、復元画像データxの観測データyに対する適合度を示し、事前知識適合度E
p(x)は、復元画像データxの事前知識に対する適合度を示す。データ適合度E
d(x)及び事前知識適合度E
p(x)を最適化するためには、最適化演算の繰り返し毎に、各々の適合度の大きさのバランスが保たれることが望ましい。
【0086】
本実施形態のように最適化演算の繰り返し回数に応じて使用するチャネル数を変化させる場合、データ適合度E
d(x)の値は、観測データサブセットy
subを用いた場合と、観測データフルセットy
fullを用いた場合とで、データセットの違いにより変化し得る。一方、事前知識適合度E
p(x)の値は、観測データフルセットy
full又は観測データサブセットy
subのどちらを用いても、データ適合度E
d(x)ほど大きく変化しない。このため、観測データサブセットy
subを用いた場合と、観測データフルセットy
fullを用いた場合とでデータ適合度E
d(x)及び事前知識適合度E
p(x)の大きさのバランスが変化してしまう可能性がある。
【0087】
したがって、各々の適合度のバランスを保つため、データ適合度E
d(x)の値は、観測データフルセットy
full又は観測データサブセットy
subのどちらを用いても、大きく変化しないようにスケーリングされることが望ましい。
【0088】
具体的には、例えば、最適化演算部134は、観測データサブセットy
subを適用する場合、式(2)に代えて、以下の式(13)を用いて、データ適合度E
d(x)をスケーリングしてもよい。
【0090】
ここで、係数β
Mは正の実数であり、データ適合度のスケーリングファクタを示す。
【0091】
より具体的には、例えば、最適化演算部134は、例えば、以下の式(14)に示すように、チャネルサブセットの数と、チャネルフルセットの数との割合に応じて、データ適合度E
d(x)をスケーリングしてもよい。
【0093】
ここで、M
subはチャネルサブセットのチャネル数であり、M
fullはチャネルフルセットのチャネル数である。
【0094】
また、最適化演算部134は、例えば、以下の式(15)に示すように、観測データサブセットy
subのノルムの大きさと、観測データフルセットy
fullのノルムの大きさとに応じて、データ適合度E
d(x)をスケーリングしてもよい。
【0096】
なお、観測データyがベクトルではなく、行列の場合は、式(10)におけるL2ノルムは、フロベニウスノルムとすればよい。
【0097】
図4は、本実施形態の第1変形例に係る画像データ復元動作の一例を示すフローチャートである。
図4では、本実施形態に係る
図3に対して、更にステップST35が追加されている。
図4の動作は、ステップST35以外の動作については、
図3と同様である。以下では、最適化演算部134による最適化演算に係る部分の動作(ステップST30、ST35、及びST40)について説明し、その他の動作(ステップST10、ST20、及びST50以降)を省略する。
【0098】
ステップST30において、チャネル選択部131は、k回目の最適化演算に際して使用する観測データyとして、観測データサブセットy
subを適用する。具体的には、チャネル選択部131は、ステップST20において選択されたチャネルサブセットに対応する観測データサブセットy
subを抽出し、最適化演算部134に送信する。
【0099】
ステップST35において、最適化演算部134は、データ適合度E
d(x)をスケーリングする。これにより、データ適合度E
d(x)及び事前知識適合度E
p(x)の値のバランスは、観測データフルセットy
fullを適用した際と同等となる。
【0100】
ステップST40において、最適化演算部134は、観測データサブセットy
subを用いて、復元モデル132に基づくデータ適合度E
d(x
k)と、事前知識モデル133に基づく事前知識適合度E
p(x
k)と、を最適化するためのk回目の最適化演算を実行し、中間画像データx
kを生成する。最適化演算部134は、中間画像データx
kを判定部135に送信する。
【0101】
以上のように動作させることにより、観測データサブセットy
subを適用した場合と、観測データフルセットy
fullを適用した場合とで、データ適合度E
d(x)及び事前知識適合度E
p(x)のバランスが大きく変化することを抑制することができる。このため、観測データサブセットy
subを適用した最適化問題は、観測データフルセットy
fullを適用した最適化問題に近似した問題とみなすことができる。したがって、最適化演算を繰り返し実行する際に、適用する観測データのチャネル数を動的に変更した場合においても、収束性を損なうことなく最適化演算を実行することができ、ひいては、繰り返し演算を安定して動作させることができる。
【0102】
2.2 第2変形例
また、本実施形態では、最適化演算の繰り返しの最後の1回に観測データフルセットy
fullを適用し、残りの最適化演算について観測データサブセットy
subを適用する場合について説明したが、これに限られない。例えば、観測データフルセットy
fullは、最適化演算の繰り返しの任意の回数において、適宜適用されてもよい。より具体的には、判定部135は、最適化演算が実行されるたびに、打切り条件を満たすか否かに加えて、チャネルサブセットを用いるか否かを更に判定してもよい。
【0103】
判定部135は、例えば、チャネルサブセットを用いるか否かを、最適化演算の繰り返し数に応じて判定してもよい。より具体的には、例えば、判定部135は、繰り返し演算の1回目には、チャネルサブセットを用いない(チャネルフルセットを用いる)旨を判定し、以降の繰り返しの際には、チャネルサブセットを用いる旨を判定してもよい。このように、繰り返し回数の初期にチャネルフルセットを用いることによって、最適化演算の初期値として復元画像データxと適合していないデータが入力された場合においても、初期に生成される中間画像データx
kの精度を高めることができる。これにより、最適化演算の収束性を安定させることができる。
【0104】
なお、判定部135は、上述の例に限らず、任意の繰り返し数における最適化演算に対して、チャネルフルセットを用いる旨を判定してもよい。具体的には、例えば、判定部135は、最初の数回の最適化演算に対してチャネルフルセットを用いるように判定を行っても良いし、最後の数回の最適化演算に対してチャネルフルセットを用いるように判定を行っても良い。
【0105】
図5は、本実施形態の第2変形例に係る画像データ復元装置における画像データ復元動作の一例を示すフローチャートである。
図5では、本実施形態に係る
図3に対して、ステップST30に代えて、ステップST32、ST34、及びST36が追加されている。
図5の動作は、ステップST32、ST34、及びST36以外の動作については、
図3と同様である。以下では、最適化演算部134による最適化演算に係る部分の動作(ステップST32、ST34、ST36、及びST40)について説明し、その他の動作(ステップST10、ST20、及びST50以降)を省略する。
【0106】
ステップST32において、判定部135は、k回目の最適化演算に際して、チャネルサブセットを用いるか否かを判定する。
【0107】
判定部135は、チャネルサブセットを用いると判定した場合(ステップST32;yes)、(k+1)回目の最適化演算に際して、チャネルサブセットを用いる旨の判定結果をチャネル選択部131に送信し、ステップST34に進む。ステップST34において、チャネル選択部131は、ステップST32における判定部135の判定結果に基づき、(k+1)回目の最適化演算に際して、観測データサブセットy
subを適用する。
【0108】
一方、判定部135は、チャネルサブセットを用いないと判定した場合(ステップST32;no)、(k+1)回目の最適化演算に際して、チャネルフルセットを用いる旨の判定結果をチャネル選択部131に送信し、ステップST36に進む。ステップST36において、チャネル選択部131は、ステップST32における判定部135の判定結果に基づき、(k+1)回目の最適化演算に際して、観測データフルセットy
fullを適用する。
【0109】
ステップST40において、最適化演算部134は、観測データサブセットy
sub又は観測データサブセットy
fullを用いて、復元モデル132に基づくデータ適合度E
d(x
k+1)と、事前知識モデル133に基づく事前知識適合度E
p(x
k+1)と、を最適化するための(k+1)回目の最適化演算を実行し、中間画像データx
k+1を生成する。最適化演算部134は、中間画像データx
k+1を判定部135に送信する。
【0110】
以上のように動作させることにより、最適化演算部134は、任意の繰り返し回数における最適化演算に対して、チャネルフルセットを用いるか、チャネルサブセットを用いるかを使い分けることができる。これにより、最適化演算部134は、復元画像データxの画質に大きな影響を与える最適化演算に対してはチャネルフルセットを用い、軽微な影響しか与えない最適化演算に対してはチャネルサブセットを用いることができる。
【0111】
より具体的には、例えば、繰り返しの序盤にチャネルフルセットを用いることにより、復元画像データxの初期値の精度が低い場合においても、収束性を高めることができる。また、例えば、繰り返しの終盤にチャネルフルセットを用いることにより、復元画像データxの最終的な画質の劣化を抑制することができる。
【0112】
2.3 第3変形例
また、本実施形態では、センサ2からの観測データyについて、実際のチャネル(センサ又はスペクトル)に基づいて最適化演算を行う場合について説明したが、これに限られない。例えば、画像データ復元部13は、実際のチャネルに対応付けられた観測データyを、仮想のチャネル(仮想チャネル)に対応付けられた観測データy’にチャネル変換し、当該観測データy’を用いて最適化演算を行ってもよい。
【0113】
図6は、本実施形態の第3変形例に係る画像データ復元部の構成の一例を示すブロック図である。
図6に示すように、画像データ復元部13は、チャネル変換部136を更に備えている。
【0114】
チャネル変換部136は、記憶部12から観測データフルセットy
fullを読み出すと、当該観測データフルセットy
fullの実際のチャネルを仮想チャネルに変換し、観測データフルセットy’
fullを生成する。具体的には、例えば、チャネル変換部136は、観測データフルセットy
fullに対して主成分分析、又は特異値分解に基づく直交変換を含む線形変換を施す。これにより、観測データフルセットy
fullは、固有値が比較的大きいチャネルと、ほぼ“0”に近い値となるチャネルとに分かれるような仮想チャネルに対応する観測データフルセットy’
fullに変換される。チャネル変換部136は、変換された観測データフルセットy’
fullをチャネル選択部131に送信する。
【0115】
チャネル選択部131は、観測データフルセットy’
fullに対応する複数の仮想チャネルのうち、最適化演算部134に入力される観測データに対応する仮想チャネルを選択する。チャネル選択部131は、複数の仮想チャネルの全て(仮想チャネルフルセット)を選択する場合、観測データフルセットy’
fullを最適化演算部134に送信する。また、チャネル選択部131は、仮想チャネルフルセットのうちの一部(仮想チャネルサブセット)を選択する場合、当該仮想チャネルサブセットに対応する観測データを観測データサブセットy’
subとして抽出し、最適化演算部134に送信する。観測データサブセットy’
subは、例えば、観測データフルセットy’
fullを構成する列ベクトルのうち、固有値の比較的大きい仮想チャネルに対応する列を含む新たな列ベクトルとして抽出される。
【0116】
このように構成することにより、チャネル数を削減した最適化演算を行うことによる画像データの復元精度の劣化を抑制することができる。
【0117】
図7は、本実施形態の第3変形例に係る画像データ復元装置における画像データ復元動作の一例を示すフローチャートである。
図7では、本実施形態に係る
図3に対して、ステップST15が追加されている。以下では、
図3と同様な動作については説明を適宜省略しつつ、全体の動作について説明する。
【0118】
図7に示すように、ステップST10において、観測データ取得部11は、センサ2において観測された、複数のチャネルに対応する観測データyを取得し、記憶部12に記憶させる。
【0119】
ステップST15において、チャネル変換部136は、観測データフルセットy
fullに対して線形変換を施し、実際のチャネルから仮想チャネルへの変換を行う。これにより、チャネル変換部136は、観測データフルセットy’を生成し、チャネル選択部131に送信する。
【0120】
ステップST20Aにおいて、チャネル選択部131は、仮想チャネルフルセットのうち、仮想チャネルサブセットとして使用する仮想チャネルを選択する。
【0121】
ステップST30Aにおいて、チャネル選択部131は、k回目の最適化演算に際して使用する観測データyとして、観測データサブセットy’
subを適用する。具体的には、チャネル選択部131は、ステップST20Aにおいて選択された仮想チャネルサブセットに対応する観測データサブセットy’
subを抽出し、最適化演算部134に送信する。
【0122】
ステップST40及びST50に係る動作は、
図3と同様であるため、その説明を省略する。
【0123】
ステップST60において、チャネル選択部131は、k回目の最適化演算に対する判定結果に基づき、(k+1)回目の最適化演算に観測データフルセットy’
fullを適用する。具体的には、チャネル選択部131は、仮想チャネルフルセットに対応する観測データフルセットy’
fullを抽出し、最適化演算部134に送信する。
【0124】
ステップST70及びST80に係る動作は、
図3と同様であるため、その説明を省略する。
【0125】
以上で、画像データ復元動作が終了する。
【0126】
このように、チャネル変換部136は、観測データyの実際のチャネルを仮想チャネルに変換する。これにより、チャネル選択部131は、有意な情報を有するチャネルの数を絞り込むことができる。チャネル選択部131は、例えば、固有値の大きい仮想チャネルを優先的に選択することにより、観測データyが本来有する情報を大きく損なうことなく、最適化演算に用いるチャネル数を減らすことができる。したがって、画質の劣化を抑制しつつ、画像データの復元に要する時間を短縮することができる。
【0127】
3. 適用例
上記実施形態及び各変形例は、種々の装置に適用することが可能である。以下に、いくつかの具体的な適用例について説明する。
【0128】
3.1 MRI装置への適用例
まず、MRI装置への適用例について説明する。
【0129】
上述の通り、本実施形態をMRI装置に適用する場合、最適化演算部134による最適化演算の繰り返しは、k空間データから再構成画像を再構成する処理に相当する。すなわち、観測データy及び復元画像データxはそれぞれ、k空間データ及び再構成画像に相当し、チャネルは、k空間データを観測するためのコイルに相当する。
【0130】
データ適合度E
d(x)に用いられる線形変換Aは、例えば、コイル感度行列Sと、フーリエ変換行列Fと、を乗じて得られる行列を含む。
【0131】
コイル感度行列Sは、復元画像データに対するコイルごとの感度を示す。コイル感度行列Sは、コイル間でSNRが一致しない場合や、コイル間の相関がある場合には、その補正項を含んでもよい。
【0132】
フーリエ変換行列Fは、センサ2がk空間データをCartesian Scanと呼ばれる手法を用いて観測する場合、離散フーリエ変換(Discrete Fourier Transform)として定義される。また、フーリエ変換行列Fは、センサ2がk空間データをNon-Cartesian Scanと呼ばれる手法を用いて観測する場合、非等方離散フーリエ変換(Nonequispaced Discrete Fourier Transform)として定義される。
【0133】
事前知識適合度E
p(x)に用いられる線形変換Tは、例えば、Total Variation項を含む。Total Variation項は、例えば、隣接画素間の輝度差又は微分に関する情報を含む。すなわち、事前知識適合度E
p(x)は、復元画像データxの空間的な滑らかさを示す制約項として機能する。
【0134】
なお、MRI装置が複数の時刻にわたって時系列に観測を行う場合、観測データy及び復元画像データxは、時刻情報を含んでもよい。具体的には、或る時刻における観測データy及び復元画像データxが列ベクトルで表現される場合、時系列の観測データy及び復元画像データxは、行列の形式で表現されてもよい。
【0135】
時系列の観測データy及び復元画像データxが行列の形式で表現される場合、データ適合度E
d(x)及び事前知識適合度E
p(x)はそれぞれ、以下の式(16)及び式(17)ように定式化することができる。
【0137】
事前知識適合度E
p(x)に対してトレースノルムを用いる場合、トレースノルムを含む式に対して最適化問題を解くことになる。これにより、事前知識適合度E
p(x)は、時間的な滑らかさを示す制約項としても機能する。
【0138】
トレースノルムを含む式に関する最適化問題は、特異値判定法(Singular Value Thresholding)を適用することにより、閉じた解を得ることができることが知られている。特異値判定法は、以下の手順(i)〜(iii)によって算出される行列を出力する演算である。
【0139】
(i)行列で表現される復元画像データxに対して、特異値分解x=UΛV
*を行う。
【0140】
(ii)対角行列Λの対角成分をベクトル化したデータに対して軟判定閾値法を適用し、係数Λ’を得る。
【0141】
(iii)係数Λ’に演算UΛ’V
*を適用することによって求められた行列を出力する。
【0142】
3.2 マルチスペクトルカメラへの適用例
次に、マルチスペクトルカメラへの適用例について説明する。
【0143】
上述の通り、本実施形態をマルチスペクトルカメラに適用する場合、最適化演算部134による最適化演算の繰り返しは、可視化に用いるチャネルに混入するノイズ除去処理に相当する。すなわち、観測データy及び復元画像データxは、ノイズを含む画像データ、及びノイズ除去された画像データに相当し、チャネルは、スペクトルに相当する。
【0144】
この場合、データ適合度E
d(x)に用いられる線形変換Aは、例えば、恒等変換行列を含む。
【0145】
事前知識適合度E
p(x)に用いられる線形変換Tは、例えば、事前学習によって得られる基底行列を含む。基底行列は、例えば、同一のマルチスペクトルカメラで事前に撮像した画像データに対して主成分分析を適用することにより算出される。事前知識適合度E
p(x)は、例えば、復元画像データxを事前学習によって得られた低ランク直交基底に射影した際の誤差として定式化される。
【0146】
3.3 スペクトラルX線装置への適用例
次に、スペクトラルX線装置への適用例について説明する。スペクトラルX線装置は、CT装置、デュアルエナジー装置、及びフォトンカウンティング装置等を含む。
【0147】
上述の通り、本実施形態をX線スペクトラム装置に適用する場合、最適化演算部134による最適化演算の繰り返しは、生データ(raw data)から物質弁別画像を再構成する処理に相当する。すなわち、観測データy及び復元画像データxはそれぞれ、生データ及び再構成画像に相当し、チャネルは、スペクトルに相当する。
【0148】
データ適合度E
d(x)に用いられる線形変換Aは、例えば、物質弁別画像と生データとの変換行列を含む。当該変換行列は、生データがサイノグラムである場合、ラドン変換を含む。
【0149】
事前知識適合度E
p(x)に用いられる線形変換Tは、例えば、Total Variation項を含む。
【0150】
4. その他
また、上述の実施形態は、以下に示す変更が適宜可能である。
【0151】
例えば、上述の実施形態では、最適化演算部134は、1つの機能ブロックとして示されているが、当該1つの機能ブロックは、複数の演算プロセッサを実装することによって実現されてもよい。具体的には、最適化演算部134は、データ適合度E
d(x)を評価するための演算部分について、当該複数の演算プロセッサによる分散処理を行うように構成されてもよい。
【0152】
すなわち、データ適合度E
d(x)を評価するための演算は、複数のチャネルごとに独立して行い得る。このため、複数の演算プロセッサの各々は、データ適合度E
d(x)の評価の際に、各々に均等に配分されたチャネルに対応する部分についての演算を担当するように構成されてもよい。これにより、データ適合度E
d(x)の評価に要する演算時間を短縮すると共に、1つの演算プロセッサにかかる負荷を低減することができる。
【0153】
また、上述した実施形態における画像データ復元動作で示された各ステップは、プログラムとして記述されたソフトウェアに基づいて実行されることが可能である。当該プログラムは、汎用の計算機システムに予め記憶され、画像データ復元装置3が当該プログラムを読み込むことにより、上述した実施形態における画像データ復元動作を実現してもよい。
【0154】
なお、上述したプログラムは、画像データ復元装置3が読み取り可能な任意の記録媒体に予め記憶されていてもよい。任意の記録媒体としては、例えば、磁気ディスク(フレキシブルディスク、ハードディスク等)、光ディスク(CD−ROM、CD−R、CD−RW、DVD−ROM、DVD±R、DVD±RW等)、半導体メモリ、又はこれに類する記録媒体等が挙げられる。また、記録媒体は、画像データ復元装置3と独立した媒体に限らず、LAN(Local Area Network)及びインターネット等を介してダウンロードしたプログラムを記憶可能な画像データ復元装置3内の記憶媒体も含まれる。画像データ復元装置3は、当該記録媒体から上述したプログラムを読み込むことによって、画像データ復元動作の各ステップを画像データ復元装置3内のプロセッサに実行させることができる。なお、画像データ復元装置3は、当該記録媒体から読み込まれたプログラムの一部を、画像データ復元装置3上で稼働しているOS(オペレーティングシステム)、データベース管理ソフト、及びMW(ミドルウェア)等に実行させてもよい。
【0155】
なお、画像データ復元装置3は、パソコン、マイコン等の1つからなる装置、複数の装置がネットワーク接続されたシステム等の何れの構成であってもよい。
【0156】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。