JPWO2017006764A1 - 画像演算装置、画像演算方法、および、断層画像撮影装置 - Google Patents

画像演算装置、画像演算方法、および、断層画像撮影装置 Download PDF

Info

Publication number
JPWO2017006764A1
JPWO2017006764A1 JP2017527167A JP2017527167A JPWO2017006764A1 JP WO2017006764 A1 JPWO2017006764 A1 JP WO2017006764A1 JP 2017527167 A JP2017527167 A JP 2017527167A JP 2017527167 A JP2017527167 A JP 2017527167A JP WO2017006764 A1 JPWO2017006764 A1 JP WO2017006764A1
Authority
JP
Japan
Prior art keywords
image
update
successive approximation
vector
predicted
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2017527167A
Other languages
English (en)
Other versions
JP6454011B2 (ja
Inventor
高橋　悠
悠 高橋
後藤　大雅
大雅 後藤
廣川　浩一
浩一 廣川
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hitachi Ltd
Original Assignee
Hitachi Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hitachi Ltd filed Critical Hitachi Ltd
Publication of JPWO2017006764A1 publication Critical patent/JPWO2017006764A1/ja
Application granted granted Critical
Publication of JP6454011B2 publication Critical patent/JP6454011B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/58—Testing, adjusting or calibrating thereof
    • A61B6/582—Calibration
    • A61B6/583—Calibration using calibration phantoms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/4808—Multimodal MR, e.g. MR combined with positron emission tomography [PET], MR combined with ultrasound or MR combined with computed tomography [CT]
    • G01R33/4812—MR combined with X-ray or computed tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/73—Deblurring; Sharpening
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/06—Diaphragms
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/40—Arrangements for generating radiation specially adapted for radiation diagnosis
    • A61B6/4035—Arrangements for generating radiation specially adapted for radiation diagnosis the source being combined with a filter or grating
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5608—Data processing and visualization specially adapted for MR, e.g. for feature analysis and pattern recognition on the basis of measured MR data, segmentation of measured MR data, edge contour detection on the basis of measured MR data, for enhancing measured MR data in terms of signal-to-noise ratio by means of noise filtering or apodization, for enhancing measured MR data in terms of resolution by means for deblurring, windowing, zero filling, or generation of gray-scaled images, colour-coded images or images displaying vectors instead of pixels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Pathology (AREA)
  • General Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Optics & Photonics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pulmonology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Quality & Reliability (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

逐次近似法を用いながら、少ない更新回数で再構成画像を生成するために、被検体の所定の断層画像を受け取って、断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返し、生成した2回の更新処理の更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成し、予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより被検体の断層画像を生成する。

Description

本発明は、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)。
Figure 2017006764
式(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)の値が小さくなる)。
Figure 2017006764
Figure 2017006764
Figure 2017006764
ただし、式(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)で表せる。
Figure 2017006764
すなわち、式(5)により算出されるSD低減率は、初期画像x(0)のSD値を基準とした第一逐次近似法の生成画像のSD値の低減割合である。
ステップ305で更新回数が所定数Nに達したならば、ステップ307に進み、ステップ301で読み込んだ初期画像x(0)を基準として、逐次近似法が収束(または漸近)した時点での画像のSD低減率を予測する。式(4)におけるβkは逐次近似法の平滑化の強さを決定する任意のパラメータであるため、前述の再構成パラメータ(初期画像x(0)および逐次近似法の再構成画像サイズ、逆投影位相幅、再構成フィルタ関数、体軸方向の画像厚さ)や、被検体およびX線発生装置に印加する電圧および電流に応じて変化する実測投影データの値および分布といった、低減率設定のための条件が決定されれば、パラメータβkとSD低減率は一対一に対応付けられる。よって、低減率設定のための条件毎に両者の関係をテーブル化しておくことで、逐次近似法の平滑化強さに対するSD低減率の予測値を算出できる。ステップ307では、SD低減率の予測値を下式(6)により画素毎に算出する。
Figure 2017006764
式(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とすると、両者の間に次式が成り立つと仮定できる。
Figure 2017006764
Figure 2017006764
ただし、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である。
Figure 2017006764
Figure 2017006764
Figure 2017006764
ステップ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反復で更新後の再構成画像との和である。
Figure 2017006764
なお、式(12)において、ag T，af TのTは、ベクトルの転置を表す。
上記式(12)は、以下の仮定に基づいている。n3 > n2なる任意のn3について、n1回目の更新からn3回目の更新までの評価関数全体の更新ベクトルは次式(13)で表せる。
Figure 2017006764
ここで、式(10)、(11)のxg(n1,n2)とxf(n1,n2)を説明変数、式(13)のxF(n1,n3)を目的変数として、重回帰分析を適用することを考えると、結合係数agおよびafを用いて、次の式(14)近似が成り立つと仮定できる。
Figure 2017006764
式(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)に分離することも可能である。
Figure 2017006764
Figure 2017006764
OSアルゴリズムでは、式(2)式の最小化を行う代わりに式(16)の各辺の最小化を順次行う。式(16)の各辺の計算は、1/Mの実測投影データを用いて完了するため、全てのデータを更新に用いた場合には、式(1)と比較して、M倍の更新頻度が得られる。この場合n回目の反復かつm番目のサブセットにおける更新式は、下式(17)で表すことができる。
Figure 2017006764
式(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)。
Figure 2017006764
Figure 2017006764
そして、ステップ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)を求める。
Figure 2017006764
この予測更新画像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)のように書き換える。
Figure 2017006764
この式(21)を用いて、上述の式(2)を式(22)のようにおく。
Figure 2017006764
式(22)に罰則パラメータρを導入してAugmented Lagrangian methodを適用し、2つの変数xおよびuを次式に基づいて交互に更新する更新式が式(23.1)、(23.2)、(23.3)のように得られる。
Figure 2017006764
一方、OSアルゴリズムをデータフィデリティ項のみに適用すると、次式(24)、(25)が得られる。
Figure 2017006764
Figure 2017006764
式(23.1)に、式(24)および(25)のOSアルゴリズムを適用することで、次式(26.1)、(26.2)、(26.3)の更新式が得られる。
Figure 2017006764
 式(23.1)〜(23.3)および式(26.1)〜(26.3)において任意の定数L>||A||2 2かつ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)のように書き換える。
Figure 2017006764
式(27)において、正則化項f(z)＝0の場合には、式(27)は、式(28)で表される。
Figure 2017006764
よって、n反復目かつmサブセット目のデータフィデリティ項由来の更新ベクトルγg(n,m)は次式(29)で表される。
Figure 2017006764
一方、n反復目かつmサブセット目の正則化項由来の更新ベクトルγf(n、m)は次式で表される。
Figure 2017006764
また、式(27)の反復過程で更新ベクトルγf(n,m)を蓄積することもできる。式(27)の部分問題の更新インデックスをωとおき、暫定画像をz(ω)とおく。さらに、収束時点での反復回数をΩとおき、反復の初期値をz(0)＝x(n,m-1)−t/ps(n,m)とおくと、γf(n,m)は次式(31)で表される。
Figure 2017006764
ステップ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)で表される。
Figure 2017006764
Figure 2017006764
Figure 2017006764
Figure 2017006764
但し、式(35)では3次元画像の任意のスライスのうち、画像の縦および横のインデックスp＝｛1，…，P｝およびq＝｛1，…，Q｝と画素のインデックスjを次式で対応付けることとする。
Figure 2017006764
式(33)および式(34)は、それぞれanisotropic TVおよびisotropic TVと呼ばれ、用途に応じて選択して採用される。
ここで、式(33)における∇vxおよび∇hxを新たな変数vvおよびvhで置き換え、次式(37)のように書き換える。
Figure 2017006764
式(2)、式(32)および式(37)より次式(38)を得る。
Figure 2017006764
さらに、罰則パラメータρを導入してSplit Bregman method を適用して、次の更新式(39)、(40)を得る。
Figure 2017006764
Figure 2017006764
式(38)について、3つの変数ベクトルを交互に最小化すると、下式(41.1)〜(41.3)が得られる。
Figure 2017006764
ここで、式(41.1)のx(n) p,qは、Gauss-Seidel法より、式(42)で表される。
Figure 2017006764
式(41.2)および式(41.2)のshrink関数は次式(43)で表される。
Figure 2017006764
以上の数式のうち、式(39)、式(40)、式(41.1)、式(41.2)および式(41.3)を更新式として用いる。
本実施形態の概念は、実施形態1および2と同様であり、データフィデリティ項由来の更新ベクトルおよび正則化項由来の更新ベクトルを算出する。まず、説明の簡単のために式(41.1)の右辺をxに関して偏微分し、次式を得る。
Figure 2017006764
式(44)と対比して、式(42)におけるn反復目のデータフィデリティ項由来の更新ベクトルγg(n)は次式(45)，(46)で表される。
Figure 2017006764
Figure 2017006764
一方、n反復目の正則化項由来の更新ベクトルγf(n)は次式(47)、(48)で表される。
Figure 2017006764
Figure 2017006764
ここで、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)を用いることができる。
Figure 2017006764
Figure 2017006764
ステップ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 高速演算装置

Claims (11)

  1. 被検体の所定の断層画像を受け取って、前記断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返す第一逐次近似処理部と、
    前記第一逐次近似処理部の生成した2回の更新処理の前記更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成する予測更新画像生成部と、
    前記予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより前記被検体の断層画像を生成する第二逐次近似処理部とを有することを特徴とする画像演算装置。
  2. 請求項1に記載の画像演算装置であって、前記予測更新画像生成部は、前記更新ベクトルを、2以上のベクトル成分に分離し、2以上のベクトル成分をそれぞれ異なる係数倍した後、加算することにより、前記予測更新ベクトルを生成することを特徴とする画像演算装置。
  3. 請求項2に記載の画像演算装置であって、前記第一逐次近似処理部は、前記逐次近似法による更新処理に用いる更新式として、データフィデリティ項と正則化項を含む式を用い、
    前記予測更新画像生成部は、前記更新ベクトルを、前記データフィデリティ項に由来するベクトル成分と、前記正則化項に由来するベクトル成分に分離することを特徴とする画像演算装置。
  4. 請求項1に記載の画像演算装置であって、前記予測更新画像生成部は、前記所定の断層画像の画素値から指標を算出し、前記指標に対応する係数を求め、求めた係数を前記更新ベクトルに掛けることを特徴とする画像演算装置。
  5. 請求項4に記載の画像演算装置であって、前記指標は、前記所定の断層画像の画素値のSD(標準偏差)低減率の予測値であることを特徴とする画像演算装置。
  6. 請求項4に記載の画像演算装置であって、前記予測更新画像生成部は、予め求めておいた前記指標と前記係数の値との関係を示す関数、または、テーブルを有し、前記算出した前記指標に対応する前記係数の値を前記関数または前記テーブルに基づいて求めることを特徴とする画像演算装置。
  7. 請求項1に記載の画像演算装置であって、前記第一逐次近似処理部は、前記更新処理を予め定めた回数行うことを特徴とする画像演算装置。
  8. 請求項1に記載の画像演算装置であって、前記第一逐次近似処理部は、前記更新画像について所定の指標を算出し、前記所定の指標が所定値以下になった場合、前記更新処理を終了することを特徴とする画像演算装置。
  9. 請求項8に記載の画像演算装置であって、前記所定の指標は、前記更新画像の画素値のSD低減率、前記SD低減率の変化率、前記更新画像の画素値の所定の評価関数の値、および、前記評価関数の値の変化率のいずれかを含むことを特徴とする画像演算装置。
  10. 被検体の所定の断層画像を初期画像として逐次近似法により更新処理して更新画像を得る処理を2回以上繰り返し、
    2回の更新処理の前記更新画像の差に相当する更新ベクトルを所定の係数倍した予測更新ベクトルを生成し、これを用いて予測更新画像を生成し、
    前記予測更新画像を新たな初期画像として、逐次近似法により更新処理して更新画像を得る処理を繰り返すことにより前記被検体の断層画像を生成することを特徴とする画像演算方法。
  11. 被検体の断層画像を撮影する撮影部と、前記撮影部の撮影した前記断層画像を用いて前記被検体の断層画像を生成する画像演算装置とを有する断層画像撮影装置であって、
    前記画像演算装置は、請求項1に記載の画像演算装置であることを特徴とする断層画像撮影装置。
JP2017527167A 2015-07-08 2016-06-22 画像演算装置、画像演算方法、および、断層画像撮影装置 Active JP6454011B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2015137260 2015-07-08
JP2015137260 2015-07-08
PCT/JP2016/068470 WO2017006764A1 (ja) 2015-07-08 2016-06-22 画像演算装置、画像演算方法、および、断層画像撮影装置

Publications (2)

Publication Number Publication Date
JPWO2017006764A1 true JPWO2017006764A1 (ja) 2018-04-05
JP6454011B2 JP6454011B2 (ja) 2019-01-16

Family

ID=57685572

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2017527167A Active JP6454011B2 (ja) 2015-07-08 2016-06-22 画像演算装置、画像演算方法、および、断層画像撮影装置

Country Status (4)

Country Link
US (1) US10410343B2 (ja)
JP (1) JP6454011B2 (ja)
CN (1) CN107530043B (ja)
WO (1) WO2017006764A1 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10559100B2 (en) * 2017-05-22 2020-02-11 Prismatic Sensors Ab Method and devices for image reconstruction
US20200175732A1 (en) * 2017-06-02 2020-06-04 Koninklijke Philips N.V. Systems and methods to provide confidence values as a measure of quantitative assurance for iteratively reconstructed images in emission tomography
JP7114347B2 (ja) * 2018-06-04 2022-08-08 浜松ホトニクス株式会社 断層画像予測装置および断層画像予測方法
CN112215779B (zh) * 2020-10-28 2023-10-03 青岛大学 一种图像处理方法、装置、设备及计算机可读存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110044546A1 (en) * 2006-02-13 2011-02-24 Pan Xiaochuan M Image Reconstruction From Limited or Incomplete Data
WO2012050149A1 (ja) * 2010-10-14 2012-04-19 株式会社 日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
JP2014000408A (ja) * 2012-06-15 2014-01-09 Toshiba Corp 反復再構成方法およびシステムにおいて適応的に決定されるパラメータ値
WO2014030543A1 (ja) * 2012-08-22 2014-02-27 株式会社 日立メディコ Ｘ線ｃｔ装置およびｘ線ｃｔ画像生成方法
JP2015080719A (ja) * 2013-10-24 2015-04-27 株式会社東芝 Ｘ線コンピュータ断層撮影装置、医用画像処理装置および医用画像処理方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7831097B2 (en) 2006-03-16 2010-11-09 Siemens Medical Solutions Usa, Inc. System and method for image reconstruction
US8472688B2 (en) * 2008-04-17 2013-06-25 Wisconsin Alumni Research Foundation Method for image reconstruction employing sparsity-constrained iterative correction
US8620054B2 (en) * 2010-09-08 2013-12-31 SOCPRA—Sciences Sante et Humaines, s.e.c. Image reconstruction based on accelerated method using polar symmetries
JP5858928B2 (ja) * 2010-12-10 2016-02-10 株式会社日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
JP6218334B2 (ja) 2012-11-30 2017-10-25 株式会社日立製作所 Ｘ線ｃｔ装置及びその断層画像撮影方法
CN103971349B (zh) * 2013-01-30 2017-08-08 上海西门子医疗器械有限公司 计算机断层扫描图像重建方法和计算机断层扫描设备
WO2014172421A1 (en) * 2013-04-16 2014-10-23 The Research Foundation Iterative reconstruction for x-ray computed tomography using prior-image induced nonlocal regularization
US9261467B2 (en) * 2013-06-14 2016-02-16 General Electric Company System and method of iterative image reconstruction for computed tomography
US9208588B2 (en) * 2013-09-25 2015-12-08 Wisconsin Alumni Research Foundation Fast statistical imaging reconstruction via denoised ordered-subset statistically-penalized algebraic reconstruction technique
KR20160010221A (ko) * 2014-07-18 2016-01-27 삼성전자주식회사 의료 영상 촬영 장치 및 그에 따른 영상 처리 방법
CN105488823B (zh) * 2014-09-16 2019-10-18 株式会社日立制作所 Ct图像重建方法、ct图像重建装置及ct***
CN107592935A (zh) * 2015-03-20 2018-01-16 伦斯勒理工学院 X射线ct的自动***校准方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110044546A1 (en) * 2006-02-13 2011-02-24 Pan Xiaochuan M Image Reconstruction From Limited or Incomplete Data
WO2012050149A1 (ja) * 2010-10-14 2012-04-19 株式会社 日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
JP2014000408A (ja) * 2012-06-15 2014-01-09 Toshiba Corp 反復再構成方法およびシステムにおいて適応的に決定されるパラメータ値
WO2014030543A1 (ja) * 2012-08-22 2014-02-27 株式会社 日立メディコ Ｘ線ｃｔ装置およびｘ線ｃｔ画像生成方法
JP2015080719A (ja) * 2013-10-24 2015-04-27 株式会社東芝 Ｘ線コンピュータ断層撮影装置、医用画像処理装置および医用画像処理方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
今倉 禄浩，他２名: ""特性乗数によるＰＭＡＲＴ力学系の最適なパラメータの検討"", 電子情報通信学会技術研究報告, vol. 105, no. 547, JPN6018047617, 2006, JP, pages 49 - 54, ISSN: 0003931433 *

Also Published As

Publication number Publication date
CN107530043A (zh) 2018-01-02
US10410343B2 (en) 2019-09-10
WO2017006764A1 (ja) 2017-01-12
CN107530043B (zh) 2020-07-10
JP6454011B2 (ja) 2019-01-16
US20180165807A1 (en) 2018-06-14

Similar Documents

Publication Publication Date Title
JP6208271B2 (ja) 医用画像処理装置
Schmidt et al. A spectral CT method to directly estimate basis material maps from experimental photon-counting data
JP5978429B2 (ja) 医用画像処理装置、医用画像処理方法
US9672638B2 (en) Spectral X-ray computed tomography reconstruction using a vectorial total variation
JP6454011B2 (ja) 画像演算装置、画像演算方法、および、断層画像撮影装置
EP3050028B1 (en) Joint reconstruction of electron density images.
CN102667852B (zh) 增强图像数据/剂量减小
US9406154B2 (en) Iterative reconstruction in image formation
JP2020168353A (ja) 医用装置及びプログラム
US9261467B2 (en) System and method of iterative image reconstruction for computed tomography
US10314556B2 (en) Optimal energy weighting of dark field signal in differential phase contrast X-ray imaging
JP6753798B2 (ja) 医用撮像装置、画像処理方法及びプログラム
JP2016049455A (ja) Ｘ線コンピュータ断層撮影装置及び画像再構成装置
JP6044046B2 (ja) 動き追従ｘ線ｃｔ画像処理方法および動き追従ｘ線ｃｔ画像処理装置
JP2019535405A (ja) 最適化された再帰法によって拡張された高効率のコンピュータ断層撮影及び適用
US10970885B2 (en) Iterative image reconstruction
JP6744764B2 (ja) 画像診断装置、及び画像取得方法
Zhang et al. PET image reconstruction using a cascading back-projection neural network
JP2014526299A (ja) 制限角度トモグラフィーにおけるフィルターバックプロジェクションのための画像再構成方法
JP6569733B2 (ja) 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置
CN114387359A (zh) 一种三维x射线低剂量成像方法及装置
CN113424227B (zh) 锥束计算机断层扫描(cbct)中的运动估计和补偿
Shu et al. A Global Constraint to Improve CT Reconstruction Under Non-Ideal Conditions
Vijayalakshmi et al. Comparison of algebraic reconstruction methods in computed tomography
JP2019013268A (ja) 画像生成装置、画像生成方法及びプログラム

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20171214

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20171214

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20181204

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20181213

R150 Certificate of patent or registration of utility model

Ref document number: 6454011

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313111

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250