(57)【要約】  逐次近似法を用いながら、少ない更新回数で再構成画像を生成するために、被検体の所定の断層画像を受け取って、断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返し、生成した2回の更新処理の更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成し、予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより被検体の断層画像を生成する。
2017006764.tif 000002
【特許請求の範囲】
【請求項１】
被検体の所定の断層画像を受け取って、前記断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返す第一逐次近似処理部と、
前記第一逐次近似処理部の生成した2回の更新処理の前記更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成する予測更新画像生成部と、
前記予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより前記被検体の断層画像を生成する第二逐次近似処理部とを有することを特徴とする画像演算装置。
【請求項２】
請求項1に記載の画像演算装置であって、前記予測更新画像生成部は、前記更新ベクトルを、2以上のベクトル成分に分離し、2以上のベクトル成分をそれぞれ異なる係数倍した後、加算することにより、前記予測更新ベクトルを生成することを特徴とする画像演算装置。
【請求項３】
請求項2に記載の画像演算装置であって、前記第一逐次近似処理部は、前記逐次近似法による更新処理に用いる更新式として、データフィデリティ項と正則化項を含む式を用い、
前記予測更新画像生成部は、前記更新ベクトルを、前記データフィデリティ項に由来するベクトル成分と、前記正則化項に由来するベクトル成分に分離することを特徴とする画像演算装置。
【請求項４】
請求項1に記載の画像演算装置であって、前記予測更新画像生成部は、前記所定の断層画像の画素値から指標を算出し、前記指標に対応する係数を求め、求めた係数を前記更新ベクトルに掛けることを特徴とする画像演算装置。
【請求項５】
請求項4に記載の画像演算装置であって、前記指標は、前記所定の断層画像の画素値のSD(標準偏差)低減率の予測値であることを特徴とする画像演算装置。
【請求項６】
請求項4に記載の画像演算装置であって、前記予測更新画像生成部は、予め求めておいた前記指標と前記係数の値との関係を示す関数、または、テーブルを有し、前記算出した前記指標に対応する前記係数の値を前記関数または前記テーブルに基づいて求めることを特徴とする画像演算装置。
【請求項７】
請求項1に記載の画像演算装置であって、前記第一逐次近似処理部は、前記更新処理を予め定めた回数行うことを特徴とする画像演算装置。
【請求項８】
請求項1に記載の画像演算装置であって、前記第一逐次近似処理部は、前記更新画像について所定の指標を算出し、前記所定の指標が所定値以下になった場合、前記更新処理を終了することを特徴とする画像演算装置。
【請求項９】
請求項8に記載の画像演算装置であって、前記所定の指標は、前記更新画像の画素値のSD低減率、前記SD低減率の変化率、前記更新画像の画素値の所定の評価関数の値、および、前記評価関数の値の変化率のいずれかを含むことを特徴とする画像演算装置。
【請求項１０】
被検体の所定の断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返し、
2回の更新処理の前記更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成し、
前記予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより前記被検体の断層画像を生成することを特徴とする画像演算方法。
【請求項１１】
被検体の断層画像を撮影する撮影部と、前記撮影部の撮影した前記断層画像を用いて前記被検体の断層画像を生成する画像演算装置とを有する断層画像撮影装置であって、
前記画像演算装置は、請求項1に記載の画像演算装置であることを特徴とする断層画像撮影装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、CT(Computed Tomography)装置や磁気共鳴イメージング(以下、MRIと呼ぶ)装置等の画像演算装置における画質改善技術に関し、特に、高速な逐次近似法を提供するための画像演算装置に関する。
【背景技術】
【０００２】
人体を非侵襲的に画像化する医療用の断層画像撮像装置として、CT装置やMRI装置が広く利用されている。撮像の高精度化および診断可能な症例の拡大を目指し、CT装置やMRI装置等の断層画像撮像装置の要素技術の研究開発が行われている。CT装置やMRI装置等の画像の再構成方法には、大別して、解析的再構成法と代数的再構成法とがあることが知られている。解析的再構成法は、例えばフーリエ変換法やフィルタ補正逆投影法や重畳積分法である。一方、代数的再構成法は、例えば逐次近似再構成法である。【０００３】
近年、画質改善を目的として、逐次近似法を画像再構成アルゴリズムに適用する方法が盛んに提案されている。例えば、特許文献1の段落0005、0014〜0019には、逐次近似法(反復法)として、測定したサイノグラムと再構成画像の順投影との間の類似性を強制するデータフィデリティ項と、信号についての事前知識(例えば、鋭いエッジを有する滑らかな曲面)を強制する正則化項とを含む目的関数(評価関数)を用いる方法が開示されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開2007−244871号公報
【非特許文献】
【０００５】
【非特許文献１】S. Ahn et. al., “Globally Convergent Image Reconstruction for Emission Tomography Using Relaxed Ordered Subsets Algorithms,” IEEE. Trans. Med. Imag., vol.22, no.5, pp.613-626, 2003
【非特許文献２】A. Beck et. al., “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM J. Imaging Sci., vol.2, no.1, pp.183-202, 2009
【非特許文献３】D. Kim et. al., “Combining Ordered Subsets and Momentum for Accelerated X-ray CT Image Reconstruction,” IEEE. Trans. Med. Imag., vol.34, no.1, pp.167-178, 2015
【非特許文献４】H. Nien et. al., “Fast X-ray CT Image Reconstruction Using the Linearized Augmented Lagrangian Method with Ordered Subsets,” Technical report, [online] , 18 Feb 2014 , arXiv e-print, [平成27年7月8日検索] , インターネット<URL:http://arxiv.org/abs/1402.4381>
【非特許文献５】T. Goldstein et. al., “The split Bregman method for L1-regularized problem,” SIAM J. Imaging Sci., vol.2, no.2, pp.323-343, 2009
【非特許文献６】S. Ramani et. al., “A Splitting-Based Iterative Algorithm for Accelerated Statistical X-Ray CT Reconstruction,” IEEE. Trans. Med. Imag., vol.31, no.3, pp.677-688, 2012
【非特許文献７】H. Nien et. al., “Fast Splitting-Based Ordered-Subsets X-Ray CT Image Reconstruction,” Proc of The third international conference on image formation in X-ray computed tomography, pp.291-294, 2014
【発明の概要】
【発明が解決しようとする課題】
【０００６】
従来の逐次近似法を用いた画像再構成アルゴリズムは、多数回反復(更新)しなければ、目的関数(評価関数)が収束せず、例えば数十回(場合によっては数百回)を超える回数の反復(更新)が必要であった。そのため、一般的なCT装置やMRI装置等の断層撮像装置の演算装置で、逐次近似法を用いた場合、実用的な計算時間で再構成画像を生成することが困難であるという問題があった。【０００７】
本発明の目的は、逐次近似法を用いながら、少ない更新回数で再構成画像を生成することにある。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、本発明の画像演算装置は、被検体の所定の断層画像を受け取って、断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返す第一逐次近似処理部と、第一逐次近似処理部の生成した2回の更新処理の更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成する予測更新画像生成部と、予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより被検体の断層画像を生成する第二逐次近似処理部とを有する。
【発明の効果】
【０００９】
本発明によれば、逐次近似法を用いた画像再構成アルゴリズムの収束に要する更新回数を削減し、少ない計算時間で再構成画像を生成できる。
【図面の簡単な説明】
【００１０】
【図１】実施形態1のCT装置の外観を示す説明図
【図２】実施形態1のCT装置の構成を示すブロック図
【図３】実施形態1のCT装置の高速演算装置の動作を示すフローチャート
【図４】実施形態1のCT装置の高速演算装置の動作を示すフローチャート
【図５】(a)実施形態1の第一逐次近似処理により求めた画像(画像ベクトル)と、その更新ベクトルを示す説明図、(b)更新ベクトルをベクトル成分に分解することを示す説明図
【図６】(a)更新ベクトルを係数倍して予測更新ベクトルを求めることを示す説明図、(b)予測更新ベクトルを用いて求めた画像(画像ベクトル)から第二逐次近似処理により求めた画像(画像ベクトル)と、その更新ベクトルを示す説明図
【図７】比較例の逐次近似処理により求めた画像(画像ベクトル)と、その更新ベクトルを示す説明図
【図８】実施形態1のサブセットを用いる場合の第一逐次近似処理の動作を示すフローチャート
【図９】実施形態1のサブセットを用いる場合の第二逐次近似処理の動作を示すフローチャート
【図１０】実施形態5のCT装置の高速演算装置の動作を示すフローチャート
【図１１】実施形態5のCT装置の高速演算装置の動作を示すフローチャート
【発明を実施するための形態】
【００１１】
以下に、添付の図を用い、本発明の実施の形態の例について詳説する。また、実施形態を説明するための全図において、同一機能を有するものは同一符号を付け、その繰り返しの説明は省略する。
【００１２】
＜＜実施形態1＞＞
図1に、本実施形態に係るCT装置の外観図を、図2に、CT装置の全体構成を示す。
図1のように、CT装置は、撮影用に用いるスキャナ101、被検体をのせて移動するための寝台102、マウスやキーボードなどで構成され寝台移動速度情報や再構成位置など計測や再構成に用いるパラメータの入力を操作者から受け付ける入力装置103、スキャナ101が出力する計測データを処理する演算装置104、及び、再構成画像を表示する表示装置105を備えている。
【００１３】
図2のように、スキャナ101は、X線発生装置111、高電圧発生装置112、X線制御装置113、データ収集系129、スキャナ制御装置115および中央制御装置126を含む。高電圧発生装置112は、X線制御装置113の制御下で所定の電流及び高電圧を発生し、X線発生装置111に供給する。これにより、X線発生装置111は、X線を発生する。
【００１４】
X線発生装置111と検出器121は、中央に被検体117を挿入するための開口を備えた円盤114に搭載されている。円盤114には、円盤114を回転駆動する駆動装置116が備えられている。また、円盤114には、X線発生装置111が発生したX線が透過する位置に、Bow−tieフィルタ118とコリメータ119も搭載されている。
【００１５】
コリメータ119には、コリメータ制御装置120が接続されている。スキャナ制御装置115は、駆動装置116およびコリメータ制御装置120に接続され、円盤114の回転および停止、ならびに、コリメータ119の開口を制御する。
【００１６】
検出器121には、プリアンプ122とA/Dコンバータ123が順に接続されている。プリアンプ122は、検出器121の出力を増幅し、A/Dコンバータ123は、デジタル信号に変換し、演算装置104に受け渡す。プリアンプ122およびA/Dコンバータ123も円盤114に搭載されている。検出器121、プリアンプおよびA/Dコンバータ123は、データ収集系129を構成している。
【００１７】
寝台102は、円盤114に対して寝台102を移動させる寝台駆動部を内蔵している。寝台駆動部には、その駆動量を制御する寝台制御装置124と、寝台移動計測装置125が接続されている。
【００１８】
表示装置105と入力装置103は、入出力装置107を構成している。入出力装置107には、記憶装置108も配置されている。一方、演算装置104は、補正処理装置131と、再構成演算装置132と、画像処理装置133と、高速演算装置134とを備えている。入出力装置107と演算装置104は、操作ユニット106を構成している。
【００１９】
各部の動作を説明する。操作者が、操作ユニット106における入力装置103から撮影条件(寝台移動速度、管電流、管電圧、スライス位置など)や再構成パラメータ(関心領域、再構成画像サイズ、逆投影位相幅、再構成フィルタ関数、体軸方向の画像厚さなど)を入力すると、その指示に基づいて、中央制御装置126は、撮影に必要な制御信号を、X線制御装置113、寝台制御装置124、スキャナ制御装置115に出力する。
【００２０】
これにより、操作者が入力装置103を操作して、撮影スタート信号が出力したならば、撮影が開始される。撮影が開始されるとX線制御装置113により高電圧発生装置112に制御信号が送られ、高電圧がX線発生装置111に印加され、X線発生装置111からX線が被検体117へ照射される。同時に、スキャナ制御装置115から駆動装置116に制御信号が送られ、円盤114を回転させる。これにより、X線発生装置111、検出器121、プリアンプ122等が被検体の周りを周回する。一方、寝台制御装置124の制御により、被検体117を乗せた寝台102は、体軸方向に平行移動したり、静止したりする。
【００２１】
X線発生装置111から照射されたX線は、Bow−tieフィルタ118によってX線ビームの形状を整形された後、コリメータ119によって照射領域を制限され、被検体117に照射される。X線は、被検体117内の各組織で吸収(減衰)され、被検体117を通過し、回転方向について定められたサンプリング間隔において検出器121で検出される。この回転方向のデータ収集単位をビューと呼ぶ。なお、検出器121は、検出素子を2次元に配置した構成であり、回転方向の素子の並びをチャネル、それに直交する方向を列と呼び、収集するデータはビュー、チャネルおよび列によって識別される。
【００２２】
検出器121の各々の検出素子で検出されたX線は、電流に変換され、プリアンプ122で増幅され、A/Dコンバータ123でデジタル信号に変換されて演算装置104に出力される。
【００２３】
演算装置104の補正処理装置131は、A/Dコンバータ123からのデジタル信号について、検出素子の暗電流による出力のオフセットを補正するオフセット補正処理、エアー補正処理、リファレンス補正処理、対数変換処理、ビームハードニング効果を抑制するためのファントム補正処理等を行う。補正後のデータは、実測投影データとして入出力装置107内の記憶装置108に保存される。
【００２４】
再構成演算装置132は、実測投影データを画像再構成処理する。再構成演算装置132は、画像再構成処理を、一般的に短時間で演算が完了する解析的再構成法(例えば、フーリエ変換法やフィルタ補正逆投影法や重畳積分法等)に基づき実行する。再構成画像は、入出力装置107内の記憶装置108に保存され、表示装置105にCT画像として表示される。さらに、表示されたCT画像に対して表示断面の変更処理等を操作者が、入力装置103を介して指示した場合には、画像処理装置133が、指示された処理を実行する。
【００２５】
高速演算装置134は、再構成演算装置132が解析的再構成法に基づいて生成した再構成画像を、代数的再構成法である逐次近似法によりさらに処理し、精度の高い画像を生成する。
【００２６】
具体的には、高速演算装置134は、CPU (Central Processing Unit) とメモリ(いずれも不図示)を内蔵し、メモリに予め格納されたプログラムをCPUが読み込んで実行することにより、図3および図4のフローのように動作して逐次近似法により画像を生成する。あるいは、CPUの代わりにGPU (Graphics Processing Unit) チップを用いる構成や、FPGA (Field-Programmable Gate Array)やASIC (Application Specific Integrated Circuit) 等の集積回路(ハードウエア)が逐次近似法の演算を担う構成であっても良い。
【００２７】
まず、高速演算装置134の動作の概要を、図3、図4のフローと、図5および図6の概念図を用いて説明する。
【００２８】
高速演算装置134は、再構成演算装置132が生成した再構成画像を初期画像x(0)として、逐次近似法により画像の更新(反復)演算を所定回数N(N≧2、例えばN=3)行って、更新後の画像x(1)〜x(N)(例えばx(1)、x(2)、x(3))を得る(ステップ301〜306、図5(a)参照)。これを第一逐次近似処理と呼ぶ。
【００２９】
第一逐次近似処理では、更新(反復)演算に用いる更新式として、データフィデリティ項と、正則化項とが両方含まれるものを用いる。データフィデリティ項は、データ収集系129が収集した実測投影データと、再構成画像を順投影した順投影データとの間の類似性を評価する項である。正則化項は、生成される画像についての事前知識に基づいて再構成画像を評価する項であり、ここでは画素値の空間的な変化が少ない平坦な画像であることを強制するように作用する。
【００３０】
つぎに、高速演算装置134は、得られた2枚の更新後の画像x(N)とx(N−1)(例えばx(3)とx(2))の間の差分ベクトルに相当する更新ベクトルxφ(N-1,N)を、所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成する。具体的には、更新ベクトルxφ(N-1,N)のうち、データフィデリティ項由来のベクトル成分xg(N-1,N)と、正則化項由来のベクトル成分xf(N-1,N)を算出し、それぞれのベクトル成分に予め求めておいた所定の係数(結合係数)ag，afをかけて加算することにより、予測更新ベクトルを求め、これを用いて予測更新画像xφ(N)を生成する(ステップ307〜309、図5(b)および図6(a)参照)。これにより、複数回の更新演算を行った後の画像に近似した予測更新画像xφ(N)を、1回の演算で得ることができる。
【００３１】
つぎに、高速演算装置134は、算出した予測更新画像xφ(N)を初期画像として、所定回数だけ逐次近似法により画像の更新(反復)演算を行って更新後の画像x(P)(被検体の断層画像)を得る。(ステップ310〜315、図6(b))。これを第二逐次近似処理と呼ぶ。更新後の画像x(N)を表示等する(ステップ316)。
【００３２】
本実施形態では、更新ベクトルを、結合係数倍して予測更新ベクトルを求め、これを用いて予測更新画像を演算により求めることにより、複数回の更新演算を省略することができる。特に、データフィデリティ項由来のベクトル成分xg(N-1,N)と、正則化項由来のベクトル成分xf(N-1,N)を算出して、これらをそれぞれ結合係数倍することにより、比較例として図7に示した逐次近似処理のみを行う場合と比較して、図6(b)のように少ない更新演算回数で最適解x＊に収束(または漸近)することができる。よって、逐次近似法を用いた精度の高い画像を、高速で生成することが可能になる。
【００３３】
以下、図3、図4のフローを詳しく説明する。
【００３４】
ステップ301において、高速演算装置134は、再構成演算装置132が生成した再構成画像を記憶装置108から読み込み、初期画像x(0)とする。そして、初期画像x(0)に予め定めた更新式を適用し、更新後の画像x(1)を生成する(ステップ302、303)。更新式として、データフィデリティ項と、正則化項とが両方含まれるものを用いる。より具体的には、ここでは下記(A)および(B)の条件を満たす更新式を用いる。
【００３５】
(A)  データフィデリティ項と正則化項を評価関数に含む
(B)  画像の更新に際し、(A)の両者を同時に(同じタイミングで)用いる
一例として、本実施形態では、非特許文献1に記載されている更新式を用いる。
【００３６】
具体的な更新式は、式(1)により与えられる。ただし、x(n)は、n回更新後の画像であり、x(n−1)は、n−1回目の更新画像である。以下、x(n)を画像ベクトルとも呼ぶ。画像ベクトルの始点は、全画素の値が0となる点である(図5、6、7)。
【００３７】
【数１】
2017006764.tif 000003
【００３８】
式(1)において、∇φ(x)は、評価関数φ(x)の一次導関数である。λnは、更新の度合いを調整する緩和係数であり、更新の回数nごとに予め設定されている。Dは、評価関数φ(x)のヘッセ行列であるが、φ(x)のリプシッツ定数を用いることも可能である。φ(x)のリプシッツ定数を用いるDは、例えば、非特許文献2に記載されているものを用いることができる。
【００３９】
ここでは、評価関数φ(x)を、データフィデリティ項g(x)と、正則化項f(x)との和により式(2)のように表す。データフィデリティ項g(x)は、実測投影データから画像の状態を制約する、すなわち、データ収集系129が収集した実測投影データと、再構成画像を順投影した順投影データとの間の類似性を評価する式(3)を用いる。正則化項f(x)は、近接画素間における滑らかさを評価する式(4)を用いる。データフィデリティ項g(x)と、正則化項f(x)は、値が小さいほど最適な画像xに近く、評価が高い(すなわち評価関数φ(x)の値が小さくなる)。
【００４０】
【数２】
2017006764.tif 000004
【００４１】
【数３】
2017006764.tif 000005
【００４２】
【数４】
2017006764.tif 000006
【００４３】
ただし、式(3)において、yは、データ収集系129が収集した実測投影データのベクトルである。Aは、画像から投影データへの変換行列である。Wは実測投影データyの統計揺らぎを数値化した重み係数を対角成分に持つ行列である。
【００４４】
式(4)において、kは、画像xの対象画素から見た近接画素の方向を表すインデックスであり、Kは方向の総数である。例えば、3次元の画像で近接26画素を対象とする場合には、重複を避けるとKは13となる。βkは、k方向の正則化強さを決定する任意のパラメータであり、Ckは、画像ベクトルの近接要素を抽出するための行列である。また、j＝｛1，…，J｝は、画素のインデックスである。ξjkは、画素jとkの間の重み係数であって、一例として画素間の距離の逆数によって決定される。ψは、ポテンシャル関数であり、Huber関数やFair関数など至る所で微分可能な連続関数が用いられる。
【００４５】
なお、ここでは、データフィデリティ項g(x)と、正則化項f(x)を、上記式(3)、(4)で表したが、本実施形態は式(3)、(4)に限定されるものではなく、一次導関数∇g(x)、∇f(x)が算出可能な関数であれば、所望の関数を用いることができる。また、更新式も、上記式(1)に限らず、所望の式を用いることができる。
【００４６】
ステップ303により、画像x(0)から1回更新後の画像x(1)が生成されたならば、更新後の画像x(1)を記憶装置108に格納する。ステップ303、304の更新動作を、更新回数が所定数Nに到達するまで繰り返す(ステップ305，306)。これにより、記憶装置108には、更新後の画像x(1)〜x(N)(例えばx(1)、x(2)、x(3))が格納される。以上により第一逐次近似処理(ステップ301〜305)が完了する。
【００４７】
なお、図5、図6の例では、更新回数の所定数Nを3回としているが、この回数に限定されるものではなく、予め求めておいた回数に設定することができる。例えば、評価関数φ(x)または評価関数の変化率∇φ(x)が所定値以下に多くの場合到達する回数を、予め求めておき、その回数をNとして設定することができる。
【００４８】
また例えば、逐次近似法のノイズ低減効果の定量値の一つとして知られるSD (Standard Deviation) 低減率を算出し、SD低減率が所定値以下に多くの場合到達する回数を予め求めておき、その回数をNとして設定することもできる。画像上の画素値がほぼ一様な小領域において計測されるSD値に関し、SD低減率は、前述の解析的再構成法によって生成された再構成画像(本実施形態においては初期画像x(0)に等しい)のSD値を基準として、逐次近似法によって生成された画像のSD値が低減した割合であり、式(5)で表せる。
【００４９】
【数５】
2017006764.tif 000007
【００５０】
すなわち、式(5)により算出されるSD低減率は、初期画像x(0)のSD値を基準とした第一逐次近似法の生成画像のSD値の低減割合である。
【００５１】
ステップ305で更新回数が所定数Nに達したならば、ステップ307に進み、ステップ301で読み込んだ初期画像x(0)を基準として、逐次近似法が収束(または漸近)した時点での画像のSD低減率を予測する。式(4)におけるβkは逐次近似法の平滑化の強さを決定する任意のパラメータであるため、前述の再構成パラメータ(初期画像x(0)および逐次近似法の再構成画像サイズ、逆投影位相幅、再構成フィルタ関数、体軸方向の画像厚さ)や、被検体およびX線発生装置に印加する電圧および電流に応じて変化する実測投影データの値および分布といった、低減率設定のための条件が決定されれば、パラメータβkとSD低減率は一対一に対応付けられる。よって、低減率設定のための条件毎に両者の関係をテーブル化しておくことで、逐次近似法の平滑化強さに対するSD低減率の予測値を算出できる。ステップ307では、SD低減率の予測値を下式(6)により画素毎に算出する。
【００５２】
【数６】
2017006764.tif 000008
【００５３】
式(6)において、Fは低減率設定のための条件毎に予め準備しておいたテーブルである。
【００５４】
上記テーブルFの作成方法の一例について説明する。上述の再構成パラメータ(初期画像x(0)および逐次近似法の再構成画像サイズ、逆投影位相幅、再構成フィルタ関数、体軸方向の画像厚さ)等の低減率設定のための条件を固定し、パラメータβkを変化させながらSD低減率をプロットし、最小二乗法によって関数をフィッティングする。これを、低減率設定のための条件、例えば印加電圧を80kV、100kV、120kV、140kVというように変化させながら、それぞれの条件で関数を導出する。
【００５５】
各条件のうち、実測投影データの値や分布は撮影毎に変化するため、被検体の代わりに円筒型の水ファントムを撮影した実測投影データを用い、直径が異なるいくつかのファントムについて関数を算出しておく。SD低減率の予測値の算出には、被検体の実測投影データの面積に応じて最も面積が近いファントムの関数をテーブルFとして採用する。
【００５６】
なお、ステップ307のSD低減率の予測方法は上述の一例に限定されるものではなく、公知のあらゆる方法を適用できる。
【００５７】
画像x(0)の任意の画素jのSD低減率をbj、結合係数agのj番目の要素をag，j、結合係数afのj番目の要素をaf，jとすると、両者の間に次式が成り立つと仮定できる。
【００５８】
【数７】
2017006764.tif 000009
【００５９】
【数８】
2017006764.tif 000010
【００６０】
ただし、hg(b)およびhf(b)はそれぞれ変数bに対して一意に決定される関数であり、臨床データもしくは画質定量ファントム等についての実測投影データについて、下式(14)の近似誤差が最小となるような関数hg(b)およびhf(b)を最小二乗法等によって予め算出したものである。
【００６１】
そして、記憶装置108に予め格納しておいた、SD低減率の予測値と、結合係数ag、afとの関係を示す関数hg(b)およびhf(b)、またはテーブルに基づいて、ステップ307で算出したSD低減率の予測値に対応する結合係数ag、afを画素毎に算出、または読み出す。
【００６２】
N−1回目の更新画像からN回目の更新画像の差分(更新)ベクトルのうち、データフィデリティ項由来のベクトル成分をxg(N-1,N)、正則化項由来のベクトル成分をxf(N-1,N)で表すと、これらは、更新式である上述の式(1)および下式(9)から式(10)、(11)により表される。ただし、n1＝N−1、n2＝Nである。
【００６３】
【数９】
2017006764.tif 000011
【００６４】
【数１０】
2017006764.tif 000012
【００６５】
【数１１】
2017006764.tif 000013
【００６６】
ステップ309において、高速演算装置134は、式(10)、(11)より、更新ベクトルのうちデータフィデリティ項由来のベクトル成分xg(N-1,N)、正則化項由来のベクトル成分xf(N-1,N)を求め、ステップ308で求めた結合係数ag，afを掛けて予測更新ベクトルを求め、これと式(12)から、予測更新画像xφ(N)＝xφ(n1,n2)を求める。ただし、n1＝N−1、n2＝Nである。予測更新画像xφ(N)＝xφ(n1,n2)は、結合係数倍された2種類の更新ベクトル(予測更新ベクトル)と、n1反復で更新後の再構成画像との和である。
【００６７】
【数１２】
2017006764.tif 000014
【００６８】
なお、式(12)において、agT，afTのTは、ベクトルの転置を表す。
【００６９】
上記式(12)は、以下の仮定に基づいている。n3 > n2なる任意のn3について、n1回目の更新からn3回目の更新までの評価関数全体の更新ベクトルは次式(13)で表せる。
【００７０】
【数１３】
2017006764.tif 000015
【００７１】
ここで、式(10)、(11)のxg(n1,n2)とxf(n1,n2)を説明変数、式(13)のxF(n1,n3)を目的変数として、重回帰分析を適用することを考えると、結合係数agおよびafを用いて、次の式(14)近似が成り立つと仮定できる。
【００７２】
【数１４】
2017006764.tif 000016
【００７３】
式(14)の近似によって、ちょうど(n3 - n2)回分の反復にあたる計算を省略することができることがわかる。但し、結合係数によっては予測更新ベクトルの近似精度が低く、xφ(n1,n3)との間で誤差が大きい場合が考えられるので、近似誤差を低減するため、予測更新ベクトルにより算出した画像(予測更新画像)を初期値とし、以下のステップ310〜315でさらに逐次近似処理(第二逐次近似処理)を行う。
【００７４】
具体的には、高速演算装置134は、算出した予測更新画像xφ(N)を初期画像として記憶装置108から読み込み(ステップ310)、所定の更新式により、更新後の画像x(n＋1)を生成する(ステップ311、312)。更新式としては、上述の式(1)を用いることも可能であるし、式(1)とは異なる式を用いることももちろん可能である。生成した画像x(n＋1)は、記憶装置108に格納する(ステップ313)。これを、トータルの更新回数が所定数Pになるまで繰り返す(ステップ314，315)。更新後の画像x(P)は、記憶装置108に格納するとともに、表示装置105に表示させる(ステップ316)。
【００７５】
本実施形態では、逐次近似法を用い、最適解に近い精度の高い画像を、少ない更新回数で高速で生成することが可能である。
【００７６】
上記図3、図4のフローに、画像再構成における逐次近似法の高速化として知られるOrdered Subset (OS)アルゴリズムを上述の式(2)に適用して、次式のように評価関数を複数(M個)のサブセット(サブセット番号m)に分離することも可能である。
【００７７】
【数１５】
2017006764.tif 000017
【００７８】
【数１６】
2017006764.tif 000018
【００７９】
OSアルゴリズムでは、式(2)式の最小化を行う代わりに式(16)の各辺の最小化を順次行う。式(16)の各辺の計算は、1/Mの実測投影データを用いて完了するため、全てのデータを更新に用いた場合には、式(1)と比較して、M倍の更新頻度が得られる。この場合n回目の反復かつm番目のサブセットにおける更新式は、下式(17)で表すことができる。
【００８０】
【数１７】
2017006764.tif 000019
【００８１】
式(17)において、x(n，m)は、n回目の反復かつm番目のサブセットにおける暫定の画像ベクトルであり、x(n，M)はx(n＋1，0)と同値であって、n反復目かつMサブセット目の画像をn＋1反復目かつ0サブセット目の画像として使用することを表す。
【００８２】
式(17)の更新式を用いる場合、再構成演算装置132は、データ収集系129から受け取った実測投影データを、M個のサブセットに分割する。例えば、1000ビューの実測投影データに対しM＝5とした場合、ビュー方向について実測投影データを5ビューごとに選択し、1番目(m＝1)のサブセットに1、6、11、…、996ビュー、2番目(m＝2)のサブセットに2、7、12、…、997ビュー、…5番目(m＝5)のサブセットに5、10、15、…、1000ビューを割り当てることにより、ビュー方向に均等な5つのサブセットに分割する。
【００８３】
図3のステップ301では、高速演算装置134は、再構成演算装置132が生成した再構成画像を初期画像x(0)として読み込む。図3のステップ302、303、305、306は、図8のフローのように実行する。
【００８４】
図8のフローにおいて、まずステップS4-1およびS4-2では、反復(更新)のインデックスnおよびサブセットのインデックスmを初期化(n＝0、m＝1)する。ステップS4-3では、m(＝1)番目のサブセットの画像についてn(＝0)回目の更新を更新式(17)に従って行う。サブセットのインデックスを1増加させ(ステップS4-6)、ステップS4-3に戻って、m(＝2)番目のサブセットの画像についてn(＝0)回目の更新を更新式(17)に従って行う。これを、ステップS4-7において、サブセットのインデックスmがMに到達するまで繰り返し、mがMを超えたならば、ステップS4-8に進んで、反復のインデックスnを1増加させて、さらに繰り返す。
【００８５】
ステップS4-4で現在の反復のインデックスnが上述のn1回(n1=N−1)を超えた場合には、ステップS4-5において式(18)および式(19)に従ってγg(n,m)およびγf(n,m)を算出して、記憶装置108に格納する。さらに、サブセットのインデックスmを更新して画像更新を続ける(ステップS4-6,S4-7,S4-3)。
【００８６】
【数１８】
2017006764.tif 000020
【００８７】
【数１９】
2017006764.tif 000021
【００８８】
そして、ステップS4-8で増加させた反復のインデックスnが、上述のn2回(n2=N)を超えた場合には、第一逐次近似を終了し、図3のステップ307に進む。そして、高速演算装置134は、上述した通りステップ307，308を行った後、ステップ309において、ステップS4-5で算出したγg(n,m)およびγf(n,m)から式(18)および式(19)により、更新ベクトルのうちデータフィデリティ項由来のベクトル成分xg(n1,n2)および正則化項由来のベクトル成分xf(n1,n2)を算出し、ステップ308で求めた結合係数ag，afと式(20)から、予測更新画像xφ(N)＝xφ(n1,n2)を求める。
【００８９】
【数２０】
2017006764.tif 000022
【００９０】
この予測更新画像xφ(N)を初期画像とし(図4のステップ310)、第二逐次近似処理を行う。具体的には、図4のステップ311〜315を図9のステップS5-1〜S5-7により行う。反復開始の反復のインデックスnは、n2+1、かつ反復終了のインデックスnをn4(n4＝P)とする。まず、ステップS5-1およびS5-2では、反復のインデックスnおよびサブセットのインデックスmを初期化する。ステップS5-3では、更新式(17)に従って、画像の更新を行う。これをサブセットのインデックスmがMに到達するまで繰り返す(ステップS5-4,S5-5)。
【００９１】
ステップS5-5でサブセットのインデックスmが一巡してMに到達した場合、ステップS5-6で反復のインデックスnを更新し、さらに画像の更新を行う。ステップS5-7で、反復のインデックスnがn4(n4＝P)に達したならば、反復を終了し、ステップ316へ進んで、再構成画像を格納および表示する。
【００９２】
本実施形態では、第一および第二逐次近似処理の更新式として、式(1)、(17)を用いる例について説明したが、上述の条件(A)および(B)を満たす更新式であれば、式(1)、(17)に限定されるものではない。例えば、非特許文献3に開示されている、momentumアルゴリズムを用いた3から5ステップの更新式で構成される3種類の方法を用いることができる。
【００９３】
非特許文献3に開示されている方法は、上述の条件(A)および(B)を満たし、画像を更新するステップにおいてデータフィデリティ項と正則化項の更新ベクトルをそれぞれ蓄積することができる。よって、本実施形態と同様に、予測更新画像を生成することができる。
【００９４】
＜＜実施形態2＞＞
実施形態2について説明する。
【００９５】
実施形態2は、実施形態1と同様に、第一逐次近似処理(図3のステップ302〜305)を行ったのち、更新画像の差分ベクトル(更新ベクトル)のデータフィデリティ項のベクトル成分xg(N-1,N)と、正則化項由来のベクトル成分xf(N-1,N)を算出して、これらを結合係数倍した予測更新ベクトルを求め、これを用いて予測更新画像を演算により求める。このとき、実施形態2は、実施形態1とは異なり、第一逐次近似処理の更新式として、下記(A)および(C)の条件を満たす更新式を用いる。
【００９６】
(A)  データフィデリティ項と正則化項を評価関数に含む
(C)  画像の更新に際し、(A)のデータフィデリティ項と正則化項を交互に更新する(すなわち、データフィデリティ項を定数として正則化項を更新し、その後に正則化項を定数としてデータフィデリティ項を更新する)
具体的には、(A)および(C)の条件を満たす更新式として、非特許文献4に記載の更新式を用いる。
【００９７】
非特許文献4に記載の更新式について説明する。本実施形態の更新式は、上述の式(2)〜(4)の最適化に際して、データフィデリティ項に新たな変数を導入し、導入した変数の関係式を制約条件として定式化した上で最適化を行う。
【００９８】
まず、式(3)のAxを新たな変数uで置き換え、式(3)を次式(21)のように書き換える。
【００９９】
【数２１】
2017006764.tif 000023
【０１００】
この式(21)を用いて、上述の式(2)を式(22)のようにおく。
【０１０１】
【数２２】
2017006764.tif 000024
【０１０２】
式(22)に罰則パラメータρを導入してAugmented Lagrangian methodを適用し、2つの変数xおよびuを次式に基づいて交互に更新する更新式が式(23.1)、(23.2)、(23.3)のように得られる。
【０１０３】
【数２３】
2017006764.tif 000025
【０１０４】
一方、OSアルゴリズムをデータフィデリティ項のみに適用すると、次式(24)、(25)が得られる。
【０１０５】
【数２４】
2017006764.tif 000026
【０１０６】
【数２５】
2017006764.tif 000027
【０１０７】
式(23.1)に、式(24)および(25)のOSアルゴリズムを適用することで、次式(26.1)、(26.2)、(26.3)の更新式が得られる。
【０１０８】
【数２６】
2017006764.tif 000028
式(23.1)〜(23.3)および式(26.1)〜(26.3)において任意の定数L>||A||22かつt=1/Lである。
【０１０９】
式(26.1)、(26.2)および(26.3)は、n反復目かつmサブセット目における一連の更新手順であるが、正則化項f(x)の性質に依存して式(26.2)式におけるx(n，m)は反復的に算出する場合がある。そこで、式(26.2)は、公知のFast Iterative Shrinkage Thresholding Algorithm (FISTA、非特許文献2)を利用して、x(n，m)を反復的に算出する部分問題を解く。
【０１１０】
実施形態2においても実施形態1の式(18)、(19)と同様に、図3のステップ309においてデータフィデリティ項由来の更新ベクトルおよび正則化項由来の更新ベクトルを算出する必要がある。そのため、式(26.2)式の上記部分問題を次式(27)のように書き換える。
【０１１１】
【数２７】
2017006764.tif 000029
【０１１２】
式(27)において、正則化項f(z)＝0の場合には、式(27)は、式(28)で表される。
【０１１３】
【数２８】
2017006764.tif 000030
【０１１４】
よって、n反復目かつmサブセット目のデータフィデリティ項由来の更新ベクトルγg(n,m)は次式(29)で表される。
【０１１５】
【数２９】
2017006764.tif 000031
【０１１６】
一方、n反復目かつmサブセット目の正則化項由来の更新ベクトルγf(n、m)は次式で表される。
【０１１７】
【数３０】
2017006764.tif 000032
【０１１８】
また、式(27)の反復過程で更新ベクトルγf(n,m)を蓄積することもできる。式(27)の部分問題の更新インデックスをωとおき、暫定画像をz(ω)とおく。さらに、収束時点での反復回数をΩとおき、反復の初期値をz(0)＝x(n,m-1)−t/ps(n,m)とおくと、γf(n,m)は次式(31)で表される。
【０１１９】
【数３１】
2017006764.tif 000033
【０１２０】
ステップ309において、式(29)のγg(n,m)を式(18)のγg(n,m)に、式(30)あるいは式(31)のγf(n,m)を式(19)のγf(n,m)に代入して、更新ベクトルのベクトル成分xg(n1、n2)、xf(n1,n2)を算出することができる。
【０１２１】
以上のように、データフィデリティ項と正則化項を交互に更新する逐次近似法を用いる場合においても、更新式中で互いの項の更新成分に注目することで、それぞれの項に由来する更新ベクトルを容易に分離することができる。
【０１２２】
実施形態2では、上述した他、図3の他のステップは、第1の実施形態のサブセットを用いる処理動作を同様であるので説明を省略する。また、図4の第二逐次近似処理においても、本実施形態で説明した更新式を用いることが可能である。
【０１２３】
なお、上述の条件(A)および(C)を満たしていれば、本実施形態で説明した更新式以外の更新式を用いることももちろん可能である。
【０１２４】
＜＜実施形態3＞＞
実施形態3について説明する。
【０１２５】
実施形態3は、実施形態2と同様に、上述の条件(A)および(C)を満たし、なおかつ正則化項に新たな変数を導入した更新式を用いる。
【０１２６】
実施形態3で用いる更新式についてまず説明する。
【０１２７】
本実施形態では、条件(A)および(C)を満たし、かつ正則化項に微分不可能な連続関数を用いる。微分不可能な連続関数はTotal Variation (TV)に代表され、次式(32)〜(35)で表される。
【０１２８】
【数３２】
2017006764.tif 000034
【０１２９】
【数３３】
2017006764.tif 000035
【０１３０】
【数３４】
2017006764.tif 000036
【０１３１】
【数３５】
2017006764.tif 000037
【０１３２】
但し、式(35)では3次元画像の任意のスライスのうち、画像の縦および横のインデックスp＝｛1，…，P｝およびq＝｛1，…，Q｝と画素のインデックスjを次式で対応付けることとする。
【０１３３】
【数３６】
2017006764.tif 000038
【０１３４】
式(33)および式(34)は、それぞれanisotropic TVおよびisotropic TVと呼ばれ、用途に応じて選択して採用される。
【０１３５】
ここで、式(33)における∇vxおよび∇hxを新たな変数vvおよびvhで置き換え、次式(37)のように書き換える。
【０１３６】
【数３７】
2017006764.tif 000039
【０１３７】
式(2)、式(32)および式(37)より次式(38)を得る。
【０１３８】
【数３８】
2017006764.tif 000040
【０１３９】
さらに、罰則パラメータρを導入してSplit Bregman method を適用して、次の更新式(39)、(40)を得る。
【０１４０】
【数３９】
2017006764.tif 000041
【０１４１】
【数４０】
2017006764.tif 000042
【０１４２】
式(38)について、3つの変数ベクトルを交互に最小化すると、下式(41.1)〜(41.3)が得られる。
【０１４３】
【数４１】
2017006764.tif 000043
【０１４４】
ここで、式(41.1)のx(n)p,qは、Gauss-Seidel法より、式(42)で表される。
【０１４５】
【数４２】
2017006764.tif 000044
【０１４６】
式(41.2)および式(41.2)のshrink関数は次式(43)で表される。
【０１４７】
【数４３】
2017006764.tif 000045
【０１４８】
以上の数式のうち、式(39)、式(40)、式(41.1)、式(41.2)および式(41.3)を更新式として用いる。
【０１４９】
本実施形態の概念は、実施形態1および2と同様であり、データフィデリティ項由来の更新ベクトルおよび正則化項由来の更新ベクトルを算出する。まず、説明の簡単のために式(41.1)の右辺をxに関して偏微分し、次式を得る。
【０１５０】
【数４４】
2017006764.tif 000046
【０１５１】
式(44)と対比して、式(42)におけるn反復目のデータフィデリティ項由来の更新ベクトルγg(n)は次式(45)，(46)で表される。
【０１５２】
【数４５】
2017006764.tif 000047
【０１５３】
【数４６】
2017006764.tif 000048
【０１５４】
一方、n反復目の正則化項由来の更新ベクトルγf(n)は次式(47)、(48)で表される。
【０１５５】
【数４７】
2017006764.tif 000049
【０１５６】
【数４８】
2017006764.tif 000050
【０１５７】
ここで、M＝1として、式(44)を式(17)のγg(n，m)に、式(47)を式(19)式のγf(n，m)に代入して、更新ベクトルのベクトル成分xg(n1、n2)、xf(n1,n2)を算出することができる。
【０１５８】
以上より、データフィデリティ項と正則化項を交互に更新する逐次近似法においても、更新式中で互いの項の更新成分に注目することで、それぞれの項に由来する予測ベクトルを容易に分離することができる。
【０１５９】
なお、本実施形態は、非特許文献5に記載のアルゴリズムを使用する更新式を示したが、上述の条件(A)および(C)を満たしていれば更新式は、上述の式に限定されるものではない。
【０１６０】
例えば、非特許文献6では、データフィデリティ項だけでなく正則化項の変分離を行った上でAlternating Direction Method of Multipliers (ADMM)アルゴリズムを適用する方法およびSplit Bregman法を適用するアルゴリズムが記載されている。また、非特許文献7でも同様にデータフィデリティ項と正則化項の変数分離を行うアルゴリズムが述べられている。これらに基づく、更新式を用いることも可能である。
【０１６１】
＜＜実施形態4＞＞
実施形態4について説明する。
【０１６２】
実施形態1〜3では、図3、図4に示したように、第一逐次近似処理(ステップ302〜306)では、予め定めたN回の反復を行い、第二逐次近似処理(ステップ311〜315)では、予め定めた(P−N)回の反復処理を行う構成であった。本発明はこの構成に限られるものではない。
【０１６３】
例えば、高速演算装置134は、操作者から、第一および第二逐次近似処理の反復処理の合計回数Pの値を受け付けるように構成することができる。
【０１６４】
この場合、第一逐次近似処理と第二逐次近似処理のそれぞれの反復回数N、(P−N)は、前者が多く後者が少なければ、予測更新画像の精度はおおむね高くなるが、第二逐次近似処理の反復回数が少なすぎるため、最適解に十分に漸近できなくなる可能性がある。
【０１６５】
一方、後者が多く前者が少なければ、予測更新画像の精度が十分得られず、本発明による高速化の効果が不十分になる可能性がある。両者を鑑み、第一逐次近似処理と第二逐次近似処理のそれぞれの反復回数のバランスがとれるように、予め実験によりまたは計算によりPとNとの関係を求めておくことが望ましい。求めたPとNとの関係は、テーブルまたは関数として記憶装置108に格納しておくことが可能である。これにより、高速演算装置134は、操作者が設定したPに対応するNの値をテーブルから読み出して、または算出して用いることができる。
【０１６６】
＜＜実施形態5＞＞
実施形態5について説明する。
【０１６７】
本実施形態は、図10に示したように、図3の第一逐次近似処理(ステップ302〜306)において、反復処理の回数が所定数Nに到達したかどうかを判定するステップ305の前に、更新後の画像x(n＋1)がステップ309において予測更新画像Xφ(N)を生成するのに適した程度にまで最適解に近づいているかどうかを判断するために、指標Eを算出するステップ501，502を挿入している。この点が、実施形態1とは異なっている。図10のフローの他のステップは、図3と同じであるので説明を省略する。
【０１６８】
指標Eとしては、更新後の画像x(n＋1)の全体もしくは注目画素周辺のSD低減率、今回の更新画像x(n＋1)と前回の更新画像x(n)のSD低減率の変化率(例えば変化量)、更新後の画像x(n＋1)の評価関数φの値、今回の更新画像x(n＋1)と前回の更新画像x(n)の評価関数φの値の変化率(例えば変化量)のいずれかを用いることができる。
【０１６９】
ステップ501では、高速演算装置134は、ステップ304で生成した更新後の画像x(n＋1)について、SD低減率または評価関数φを所定の数式を用いて算出し、記憶装置108に格納する。変化量を指標Eとする場合には、前回格納したSD低減率または評価関数φの値との差分をさらに算出する。所定の数式としては、例えば、式(49)，(50)を用いることができる。
【０１７０】
【数４９】
2017006764.tif 000051
【０１７１】
【数５０】
2017006764.tif 000052
【０１７２】
ステップ502では、ステップ501で算出した指標Eを、予め定めておいた値と比較し、指標Eが予め定めた値よりも小さい場合には、更新後の画像x(n＋1)が予測更新画像Xφ(N)を生成するのに適した程度にまで最適解に近づいていると判断し、第一逐次近似処理の反復を終了して、ステップ307に進む。
【０１７３】
これにより、高速演算装置134は、反復処理の回数nがNに到達する前に、第一逐次近似処理を終了し、ステップ307〜309で予測更新画像xφ(N)を生成することができる。
【０１７４】
一方、ステップ502において、指標Eが予め定めた値よりも小さくなっていない場合には、ステップ305に進み、反復処理の回数nがNを超えていなければ、反復処理を継続する。指標Eが予め定めた値よりも小さくなる前に、反復処理の回数nがNを超えた場合には、ステップ307に進んで、反復処理を終了する。よって、反復処理の回数nの上限は、N回に維持できる。
【０１７５】
このように、本実施形態では、更新後の画像の状態から、反復処理を継続するかどうかを判断できるため、第一逐次近似処理において、必要以上に反復処理を行うことがなく、より高速に予測更新画像を生成できる。
【０１７６】
同様に、図11に示したフローのように、図4の第二逐次近似処理(ステップ311〜315)においても、反復処理の回数が所定数P
に到達したかどうかを判定するステップ314の前に、更新後の画像x(n＋1)が最適解に近づいているかどうかを判断するために、指標Eを算出するステップ511，512を挿入することが可能である。図11のフローの他のステップは、図4と同じであるので説明を省略する。
【０１７７】
ステップ511では、ステップ312で生成した更新後の画像x(n＋1)について、ステップ501と同様に、いずれかの指標Eを算出する。
【０１７８】
ステップ512では、ステップ501で算出した指標Eを、予め定めておいた値と比較し、指標Eが予め定めた値よりも小さい場合には、十分最適解に近づいていると判断し、第二逐次近似処理の反復を終了して、ステップ316に進む。なお、ステップ512の指標Eを判定する所定値は、ステップ502の所定値とは異なる値を用いることができる。
【０１７９】
これにより、高速演算装置134は、反復処理の回数nがPに到達する前に、第一逐次近似処理を終了させることができる。
【０１８０】
一方、ステップ502において、指標Eが予め定めた値よりも小さくなっていない場合には、ステップ314に進み、反復処理の回数nがPを超えていなければ、反復処理を継続する。ステップ511で求めた指標Eが予め定めた値よりも小さくなる前に、反復処理の回数nがPを超えた場合には、ステップ316に進んで反復処理を終了する。よって、第一および第二逐次近似処理の合計の反復処理の回数nの上限は、P回に維持できる。
【０１８１】
このように、第二逐次近似処理においても、更新後の画像の状態から、反復処理を継続するかどうかを判断できるため、第二逐次近似処理において、必要以上に反復処理を行うことがなく、最適解に漸近した画像をより高速に生成できる。
【０１８２】
上述してきた実施形態では、CT装置について説明したが、本実施形態の高速演算装置134の処理は、MRI装置等の他の画像撮影装置の画像生成に用いることができる。
【０１８３】
なお、実施形態1から3では評価関数がデータフィデリティ項と正則化項の2つからなる場合を説明したが、3項以上で構成される評価関数に対しては3つ以上の予測更新ベクトルを対応させるようにしても良い。
【符号の説明】
【０１８４】
101  スキャナ、102  寝台、103  入力装置、104  演算装置、105  表示装置、108  記憶装置、111  X線発生装置、121  検出器、132  再構成演算装置、134  高速演算装置
【図１】
2017006764.tif 000053
【図２】
2017006764.tif 000054
【図３】
2017006764.tif 000055
【図４】
2017006764.tif 000056
【図５】
2017006764.tif 000057
【図６】
2017006764.tif 000058
【図７】
2017006764.tif 000059
【図８】
2017006764.tif 000060
【図９】
2017006764.tif 000061
【図１０】
2017006764.tif 000062
【図１１】
2017006764.tif 000063
【手続補正書】
【提出日】平成29年12月14日(2017.12.14)
【手続補正１】
【補正対象書類名】図面
【補正対象項目名】図８
【補正方法】変更
【補正の内容】
【図８】
2017006764.tif 000064
【手続補正２】
【補正対象書類名】図面
【補正対象項目名】図９
【補正方法】変更
【補正の内容】
【図９】
2017006764.tif 000065
【国際調査報告】
2017006764.tif 000066
2017006764.tif 000067