JP5978429B2 - 医用画像処理装置、医用画像処理方法 - Google Patents

医用画像処理装置、医用画像処理方法 Download PDF

Info

Publication number
JP5978429B2
JP5978429B2 JP2013511984A JP2013511984A JP5978429B2 JP 5978429 B2 JP5978429 B2 JP 5978429B2 JP 2013511984 A JP2013511984 A JP 2013511984A JP 2013511984 A JP2013511984 A JP 2013511984A JP 5978429 B2 JP5978429 B2 JP 5978429B2
Authority
JP
Japan
Prior art keywords
penalty term
weight
true subset
projection data
detector
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.)
Active
Application number
JP2013511984A
Other languages
English (en)
Other versions
JPWO2012147471A1 (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 JPWO2012147471A1 publication Critical patent/JPWO2012147471A1/ja
Application granted granted Critical
Publication of JP5978429B2 publication Critical patent/JP5978429B2/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/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • 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/70—Denoising; Smoothing
    • 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • 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)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Optics & Photonics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Pathology (AREA)
  • Biophysics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Pulmonology (AREA)
  • Quality & Reliability (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Description

本発明は、検出器の出力に応じた重み及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理装置等に関するものである。
X線CT装置は、チャネル方向及び列方向に検出素子が配列されるX線検出器を有する。そして、X線CT装置による撮影では、互いに対向するX線管とX線検出器を被検体の周囲に回転させることにより、投影データが周回方向の離散的なX線管の位置(対向するX線検出器の位置とも言える。)において収集される。以下では、各々のX線管位置における投影データの取得単位を「ビュー」と呼ぶ。投影データを収集した後、演算装置が、各ビュー及び各列について、投影データのチャネル方向に再構成フィルタを重畳してフィルタ補正投影データを得て、フィルタ補正投影データに対して、ビュー方向に重みを加重しながら逆投影することによって、被検体内部のX線減弱係数の分布図として非破壊的に断層像を画像化する。以下では、投影データから断層像を画像化する方法を「画像再構成法」と呼び、画像再構成法によって得られた画像をCT画像と呼ぶ。
画像再構成法は、解析法と逐次近似法に大別される。解析法は、投影切断面定理に基づき解析的に問題を解く方法である。逐次近似法は、投影データの取得に至る観測系を数学的にモデル化し、数学モデルに基づいて反復法により最良の画像を推定する方法である。
両者を比較すると、解析法の利点としては、投影データから直接的に再構成画像が得られるため、計算量が圧倒的に少ないことが挙げられる。一方、逐次近似法の利点としては、投影データの取得に至る物理的過程及び投影データに含まれる統計的な揺らぎをそれぞれ数学モデルや統計モデルとして考慮できるため、解析法において発生するアーチファクト(コーンビームアーチファクトなど)や画像上の量子ノイズを低減できることが挙げられる。
従来、マルチスライスCTにおける画像再構成法としては、計算量が少ないことを採用理由として、解析法であるFeldkamp法、またはFeldkamp法を改良した手法が主に用いられている。しかし、近年のコンピュータの高性能化に伴い、逐次近似法も実用化が検討され始めている。特に、逐次近似法によって画質の向上を図る技術が提案されている。
ここでは、推定される変数の違いによって、逐次近似法を次の2種類に大別する。
(1)投影データの投影値を変数として評価関数を構築する手法
(2)画像データの画素値を変数として評価関数を構築する手法
(1)の手法は、画像再構成処理の前処理である投影データの補正処理などに適用される。以下では、(1)の逐次近似法が適用される投影データの補正処理を「逐次近似投影データ補正処理」と呼ぶ。
(2)の手法は、画像再構成処理に適用される。以下では、(2)の逐次近似法が適用される画像再構成処理を「逐次近似再構成処理」と呼ぶ。
逐次近似投影データ補正処理では、投影データの評価指標を事前に設定しておき、評価指標を数値化した評価値が最大値もしくは最小値をとるように投影データを逐次更新する。評価指標は、更新の対象である暫定投影データと実投影データとの矛盾や確率的な尤もらしさなどが用いられる。一例として、非特許文献1には、次式によって表される罰則付き加重二乗誤差関数を評価関数として用いる逐次近似投影データ補正処理が提案されている。
Figure 0005978429
 (1)式において、1，・・・i，・・・，Iは、検出器のチャネル、検出器の列、及びビューの組み合わせを一意に識別する通し番号である。pi及びyiは、それぞれ通し番号iによって特定される暫定投影データ及び実投影データである。p＝｛p1，・・・pi，・・・，pI｝は、暫定投影データのベクトルである。
κiは、通し番号iによって特定される検出素子と近接する検出素子の集合である。
νikは、通し番号iによって特定される検出素子と、通し番号kによって特定される検出素子との相関性を示す定数である。非特許文献1では、νikを経験的に決定している。
Ψ(pi−pk)は、暫定投影データpiと暫定投影データpkとのコントラストを変数とするポテンシャル関数である。非特許文献1では、ポテンシャル関数として、2次関数を使用している。
diは、実投影データyiと順投影データの差分値に加重される重み係数である。CT画像の画像処理では、diは通し番号iによって特定される検出素子の検出器出力を反映した値である為、以降では、diを「検出器出力重み」と呼ぶ。尚、検出器出力は、検出フォトン数に依存する値である。
βは、任意定数である。
逐次近似再構成処理では、画像データの評価指標を事前に設定しておき、評価指標を数値化した評価値が最大値もしくは最小値をとるように画像データを逐次更新する。評価指標は、更新の過程において暫定画像データが投影データに変換されたデータである順投影データと、実投影データとの矛盾や確率的な尤もらしさなどが用いられる。一例として、非特許文献2には、次式によって表される罰則付き加重二乗誤差関数を評価関数として用いる逐次近似再構成処理が提案されている。
Figure 0005978429
 (2)式において、1，・・・j，・・・，Ｊは、画像データの画素(2次元座標)を一意に識別する通し番号である。Xjは、通し番号jによって特定される画素における画素値である。X＝｛X1，・・・Xj，・・・，XＪ｝は、画像データを表すベクトルである。
aijは、画像データと投影データを対応付ける行列の要素である。この行列は、数学モデルを介してX線CT装置の撮影系の特性を表すことから、「システムマトリクス」と呼ばれている。
Xjは、通し番号jによって特定される画素と近接する画素の集合である。wjkは、通し番号jによって特定される画素と、通し番号kによって特定される画素との相関性を示す定数である。一般的には、wjkは、通し番号jによって特定される画素の中心位置と、通し番号kによって特定される画素の中心位置との距離の逆数である。
Ψ(Xj−Xk)は、通し番号jによって特定される画素における画素値Xjと、通し番号kによって特定される画素における画素値Xkとのコントラストを変数とするポテンシャル関数である。非特許文献2では、最も簡単なポテンシャル関数として、2次関数が提案されている。また、非特許文献3には、ポテンシャル関数の形状を任意のパラメータによって変化させる手法が提案されている。
その他の定数は、(1)式と同様である。
(1)式から種々の数値解析法により導出される更新式において、更新の過程における暫定投影データpiは、(1)式の第1項により実投影データyi及び検出器出力重みdiによる制約を受けながら、第2項の罰則項により近接する暫定投影データpkとのコントラストに応じて平滑化される。従って、(1)式の任意定数βは、暫定投影データpiの平滑化の度合いを調整する為のパラメータであると言える。
同様に、(2)式の任意定数βは、画像データXjの平滑化の度合いを調整する為のパラメータと言える。
非特許文献2及び非特許文献3では、任意定数βは、逐次近似再構成処理全体を通して一定値である。
一方、非特許文献4では、βを変更する手法が提案されている。βが一定のままでは、CT画像は、画素の位置によって空間分解能及びノイズ低減特性が不均一になる。そこで、非特許文献4では、各画素においてβを変更することによって、近接する画素のコントラストに加え、空間分解能及びノイズ低減特性の不均一性を緩和している。
T. Li et. al., "Nonlinear Sinogram Smoothing for Low-Dose X-Ray CT,"IEEE. Trans. Nucl. Sci., Vol.51, No.5,pp.2505-2513, 2004 K. Sauer et. al., "A Local Update Strategy for Iterative Reconstruction form Projections,"IEEE. Trans. signal. Proc., Vol.41, No.2, pp.534-548, 1993 J. B. Thibault et. al., "A three-dimensional statistical approach to improved image quality for multislice helical ct,"Med. Phys., Vol.34, No.11, pp.4526-4544, 2007 H. R. Shi et. al., "Quadratic Regularization Design for 2-D CT,"IEEE. Trans. Med. Imag., Vol.28, No.5, pp645-656, 2009] T. Goto et. al., "Weighted-Feldkamp algorithm with selective narrowest cone angle data for cone beam CT,"Proc. Of. The Eighth International Meeting on Fully Three-dimensional Image Reconstruction in Radiology and Nuclear Medicine, pp.189-192, 2005 J. Wang et. al., "Penalized weighted least-squares approach to sonogram noise reduction and image reconstruction for low-dose X-ray computed tomography,"IEEE. Trans. Med. imag., Vol.25, No.10, pp.1272-1283, 2006
ところで、逐次近似投影データ補正処理や逐次近似再構成処理を複数の部位(胸部及び腹部等)に適用する際、CT画像の観察者にとっては、量子ノイズ低減効果及びアーチファクト低減効果(以下、両者を区別しない場合には「ノイズ低減効果」と呼ぶ。)が、部位に因らず同等であることが望ましい。
しかしながら、複数部位に逐次近似投影データ補正処理や逐次近似再構成処理を適用する場合、検出器出力重みdiの大小は、部位間におけるCT画像のノイズ低減効果の違いの原因となる。例えば、胸部から腹部にかけて、式(1)や式(2)に基づく逐次近似法を適用した場合を考える。肺(中空になっている臓器)が存在する為にX線の減衰が少なく、検出器出力重みdiの値が大きい要素iが数多く存在する胸部と比較すると、検出器出力重みdiの値が大きい要素iが少ない腹部では、(1)式及び(2)式における検出器出力重みdiの制約が弱い。すなわち、胸部と腹部に対して同じ画像処理を施すと、ノイズ低減効果が異なることになる。ひいては、観察者にとって好ましくない、部位ごとにノイズ低減効果のばらつきがあるCT画像が生成されてしまう。
また、非特許文献4のように、各画素においてβを変更する場合、部位におけるノイズ低減効果のばらつきを抑制することはできる。しかしながら、検出器出力重みdiの大小が均等化されてしまい、CT画像のノイズ低減効果が不十分、すなわち量子ノイズ及びアーチファクトの低減が不十分になってしまう。更に、非特許文献4による手法では、各画素においてβを決定する処理に膨大な時間を要するという問題もある。
本発明は、前述した問題点に鑑みてなされたものであり、その目的とすることは、複数の部位に亘る投影データに対して、計算量を大幅に増加させることなく、投影データに含まれる全ての部位において一様な量子ノイズ及びアーチファクトの低減効果が達成された医用画像を作成することができる医用画像処理装置等を提供することである。
前述した目的を達成するために第1の発明は、検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理装置であって、前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定する真部分集合決定部と、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出する罰則項重み算出部と、前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行する逐次近似法実行部と、を備える医用画像処理装置である。
第2の発明は、検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び/又は再構成画像の再構成処理を行う医用画像処理方法であって、前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定するステップと、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出するステップと、前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行するステップと、を含む医用画像処理方法である。
本発明により、複数の部位に亘る投影データに対して、計算量を大幅に増加させることなく、投影データに含まれる全ての部位において一様な量子ノイズ及びアーチファクトの低減効果が達成された医用画像を作成することができる。
X線CT装置1の全体外観図 X線CT装置1の構成図 人体を模擬したファントムの肩部の断面像、並びに断面像に対する鉛直方向のビュー及び水平方向のビューを示す模式図 人体を模擬したファントムの冠状面像、並びに冠状面像の胸部に対応するビューの範囲及び腹部に対応するビューの範囲を示す模式図 医用画像処理全体の流れを示すフローチャート 第1の罰則項重み算出処理を示すフローチャート ヒストグラムの一例を示す図 第2の罰則項重み算出処理を示すフローチャート
以下図面に基づいて、本発明の実施形態を詳細に説明する。最初に、図1、図2を参照しながら、X線CT装置1の構成を説明する。
X線CT装置1は、X線管11や検出器12が搭載されるスキャナ2、被検体10を載置する寝台4、検出器12から得られるデータの処理を行う演算装置5、マウス、トラックボール、キーボード、タッチパネルなどの入力装置6、及び再構成画像(CT画像)などを表示する表示装置7などを含む。
操作者は、入力装置6を介して、撮影条件や再構成条件などを入力する。撮影条件は、例えば、寝台送り速度、管電流、管電圧、撮影範囲(スライス位置の範囲)、周回当たりの撮影ビュー数などである。また、再構成条件は、例えば、関心領域、再構成画像サイズ(再構成画像の大きさ)、再構成フィルタ関数などである。
図2に示すように、X線CT装置1は、大きく分けて、スキャナ2、操作ユニット3、寝台4から構成される。
スキャナ2は、X線管11(X線発生装置)、検出器12、コリメータ13、駆動装置14、中央制御装置15、X線制御装置16、高電圧発生装置17、スキャナ制御装置18、寝台制御装置19、寝台移動計測装置20、コリメータ制御装置21、プリアンプ22、A/Dコンバータ23などから構成されている。
中央制御装置15は、操作ユニット3における入力装置6から撮影条件や再構成条件を入力し、撮影に必要な制御信号を、コリメータ制御装置21、X線制御装置16、スキャナ制御装置18、寝台制御装置19に送信する。
コリメータ制御装置21は、制御信号に基づいてコリメータ13の位置を制御する。
撮影スタート信号を受けて撮影が開始されると、X線制御装置16は、制御信号に基づいて高電圧発生装置17を制御する。高電圧発生装置17は、X線管11(X線発生装置)に管電圧、管電流を印加する。X線管11では、印加された管電圧に応じたエネルギーの電子が陰極から放出され、放出された電子がターゲット(陽極)に衝突することによって電子エネルギーに応じたエネルギーのX線が被検体10に照射される。
また、スキャナ制御装置18は、制御信号に基づいて駆動装置14を制御する。駆動装置14は、X線管11、検出器12、プリアンプ22等が搭載されているガントリ部を被検体10の周りに周回させる。
寝台制御装置19は、制御信号に基づいて寝台4を制御する。
X線管11から照射されるX線は、コリメータ13によって照射領域が制限され、被検体10内の各組織においてX線減弱係数に応じて吸収(減衰)され、被検体10を通過し、X線管11に対向する位置に配置された検出器12によって検出される。検出器12は、2次元方向(チャネル方向およびこれに直交する列方向)に配置された複数の検出素子によって構成される。各検出素子によって受光したX線は、実投影データに変換される。
すなわち、検出器12によって検出されるX線は、電流に変換され、プリアンプ22によって増幅され、A/Dコンバータ23によってデジタルデータに変換され、LOG変換され、キャリブレーションが行われて実投影データとして演算装置5に入力される。
このとき、互いに対向するX線管11と検出器12が、被検体10の周囲を回転するので、実投影データは、回転方向の離散的なX線管位置(対向する検出器位置とも言える。)において収集される。各々のX線管位置における実投影データの取得単位が、「ビュー」である。
演算装置5は、再構成演算装置31、画像処理装置32等から構成される。また、入出力装置9は、入力装置6(入力部)、表示装置7(表示部)、記憶装置8(記憶部)等から構成される。
再構成演算装置31は、実投影データを用いて画像再構成処理を行い、再構成画像を生成する。再構成演算装置31は、各ビューの実投影データに再構成フィルタを重畳してフィルタ補正投影データを生成し、フィルタ補正投影データに対して、ビュー方向重みを加重して逆投影処理を行うことによって、被検体10内部のX線減弱係数の分布図として非破壊的に断層像を画像化する。
再構成演算装置31は、生成される再構成画像を記憶装置8に保存する。また、再構成演算装置31は、表示装置7にCT画像として表示する。あるいは、画像処理装置32が、記憶装置8に保存される再構成画像に対して画像処理を行い、表示装置7にCT画像として表示する。
X線CT装置1は、2次元方向に検出素子が配列された検出器12を用いるマルチスライスCT、検出素子が1列すなわち1次元方向(チャネル方向のみ)に配列された検出器12を用いるシングルスライスCTに大別される。マルチスライスCTでは、検出器12に合わせてX線源であるX線管11から円錐状、もしくは角錐状に広がるX線ビームが照射される。シングルスライスCTでは、X線管11から扇状に広がるX線ビームが照射される。通常、X線CT装置1による撮影では、ガントリ部が、寝台4に載置された被検体10の周りを周回しながら、X線の照射が行われる(但し、スキャノグラム撮影は除く)。
撮影中に寝台4が固定され、X線管11が被検体10の周りを円軌道状に周回する撮影態様は、アキシャルスキャンなどと呼ばれる。特に、寝台4を固定して撮影し、寝台4を次の撮影位置に移動させることを繰り返す撮影態様は、ステップ・アンド・シュートスキャンなどと呼ばれる。アキシャルスキャンは、1回だけ寝台4を撮影位置に移動させるステップ・アンド・シュートスキャンと考えられることから、以降では、両者を纏めてステップ・アンド・シュートスキャンと表記する。
また、寝台4が連続的に移動し、X線管11が被検体10の周りをらせん軌道状に周回する撮影態様は、らせんスキャンなどと呼ばれる。
寝台制御装置19は、ステップ・アンド・シュートスキャンの場合、撮影している間、寝台4を静止した状態とする。また、寝台制御装置19は、らせんスキャンの場合、入力装置6を介して入力される撮影条件としての寝台送りの速さに応じて、撮影している間、寝台4を体軸方向に平行移動させる。
X線CT装置1は、例えば、マルチスライスCTである。また、X線CT装置1のスキャン方式は、例えば、ローテート−ローテート方式(第3世代)である。
次に、図3、図4を参照しながら、式(1)、(2)における検出器出力重みdiについて説明する。
図3では、人体を模擬したファントムの肩部について、X線CT装置1が再構成処理によって画像化した断面像41(CT画像)と、この断面像41に対する鉛直方向のビュー(a)及び水平方向のビュー(b)を模式的に示している。説明を分かり易くする為、ここでは、検出器12の列の次元を考えず、ビュー及びチャネル42の2つの次元に限定して考える。水平方向のビュー(b)では、鉛直方向のビュー(a)と比較して人体におけるX線の透過長が長く、検出器出力が小さな値を示すチャネル42が多い。つまり、水平方向のビュー(b)では、検出器出力に大きなノイズが含まれる。尚、このことにより、解析法の画像再構成法によって再構成されたCT画像では、主に水平方向にストリーク状のアーチファクトが生ずる。
式(1)及び(2)の検出器出力重みdiは、検出器出力に応じた値である為、鉛直方向のビュー(a)に対応する検出器出力重みdiよりも、水平方向のビュー(b)に対応する検出器出力重みdiの方が小さな値となる。つまり、水平方向のビュー(b)では、鉛直方向のビュー(a)と比較して、式(1)及び(2)の第1項による制約が相対的に弱くなり、第2項(罰則項)による平滑化効果を強く受けることになる。従って、式(1)及び(2)などの評価関数を用いる逐次近似法では、検出器出力に含まれるノイズを平均化し、CT画像において生ずるストリーク状のアーチファクトを効果的に低減することができる。また、同様の原理によって、逐次近似法では、CT画像において生ずる量子ノイズを効果的に低減することができる。
CT画像の量子ノイズ低減効果及びストリークアーチファクト低減効果(＝ノイズ低減効果)は、部位に因らず同程度であることが望ましい。しかし、複数部位に亘り両者を適用する場合、検出器出力重みdiの大小は、部位間におけるCT画像のノイズ低減効果の違いの原因となる。
図4では、人体を模擬したファントムの胸部及び腹部について、X線CT装置1が再構成処理によって画像化した冠状面像43(CT画像)と、この冠状面像43の胸部に対応するビューの範囲(a)及び腹部に対応するビューの範囲(b)を模式的に示している。
説明を分かり易くする為、ここでは、検出器12のチャネルの次元を考えず、ビュー及び列44の2つの次元に限定して考える。45は、スキャン方向、つまり寝台4を移動させる方向と正反対の方向を示している。
図4のように、胸部から腹部をX線CT装置1によって画像化する場合、胸部における平均的な検出器出力の値は、腹部における平均的な検出器出力の値よりも大きい。これは、胸部は肺(中空部)を含む為、腹部等と比較して、X線の減弱量が相対的に小さいからである。
つまり、式(1)及び(2)などの評価関数を用いて、胸部から腹部にかけて、部位によらず同様の逐次近似法を適用する場合、大きい値の検出器出力重みdiが多い胸部と比較して、小さい値の検出器出力重みdiが多い腹部では、式(1)及び(2)の制約が相対的に弱くなる。従って、式(1)及び(2)などの評価関数を用いる逐次近似法では、CT画像において生ずるストリーク状のアーチファクト及び量子ノイズの低減効果(＝ノイズ低減効果)が、部位によってばらついてしまう。
本発明の医用画像処理装置では、後述する医用画像処理によって、この部位によるノイズ低減効果のばらつきを一様にする。
本発明の医用画像処理装置が実行する各ステップは、実投影データが収集された後に実行される。従って、本発明の医用画像処理装置は、X線CT装置1に含まれる演算装置5でも良いし、X線CT装置1に含まれない汎用のコンピュータでも良い。更に言えば、X線CT装置1と医用画像処理装置は、ネットワークを介して接続されていなくても良い。
混乱を避ける為に、以下では、演算装置5を本発明の医用画像処理装置として説明する。
同様に、本発明の医用画像処理装置が備える入力部、表示部、及び記憶部は、X線CT装置1に含まれる入力装置6、表示装置7、及び記憶装置8でも良いし、X線CT装置1に含まれない汎用のコンピュータが備える装置又は外部装置でも良い。混乱を避ける為に、以下では、入力装置6、表示装置7、及び記憶装置8を本発明の医用画像処理装置が備える入力部、表示部、及び記憶部として説明する。
また、実投影データは、ファンビーム系で収集されるが、以降ではパラレルビーム系に変換したデータを前提として説明する。
＜第1の実施の形態＞
図5〜図8を参照しながら、第1の実施の形態について説明する。第1の実施の形態では、本発明を逐次近似投影データ補正処理に適用する場合について説明する。第1の実施の形態では、演算装置5が、検出器出力重み(検出器12の出力に応じた重み)及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理を行う。
図5に示すように、操作者は、入力装置6を介して、撮影条件及び再構成条件、並びに、逐次近似処理の所望計算時間及び所望ノイズ低減率を入力する(ステップ1)。
撮影条件については、演算装置5は、投影データとともに記憶装置8に記憶されている値を自動的に設定すれば良い。例えば、演算装置5が、記憶装置8に記憶されている撮影条件の値を表示装置7に表示し、操作者は確認作業のみ行えば良い。
また、再構成条件、並びに、逐次近似処理の所望計算時間及び所望ノイズ低減率については、演算装置5は、表示装置7にデフォルト値や選択肢などを表示し、操作者の入力作業を支援するようにしても良い。
例えば、所望計算時間の入力支援として、演算装置5は、表示装置7に「高速」、「中速」、「低速」などの選択肢を表示するようにしても良い。この場合、演算装置5は、記憶装置8に選択肢ごとの計算時間を記憶しておく。そして、演算装置5は、操作者に選択された選択肢に応じて、所望計算時間を設定する。
また、例えば、所望ノイズ低減率の入力支援として、演算装置5は、表示装置7に「高画質」、「中画質」、「低画質」などの選択肢を表示するようにしても良い。この場合、演算装置5は、記憶装置8に選択肢ごとのノイズ低減率を記憶しておく。そして、演算装置5は、操作者に選択された選択肢に応じて、所望ノイズ低減率を設定する。
次に、演算装置5は、ステップ1において入力された撮影条件及び再構成条件に基づいて、1又は複数の真部分集合｛S1、・・・、Sm、・・・、SM｝を決定する(ステップ2)。
全体集合Ωは、検出器12のチャネル、検出器12の列、及びビューの組み合わせを一意に識別する通し番号｛1，・・・i，・・・，I｝を集合要素とする。より詳細には、検出器12のチャネル数がC、検出器12の列数がR、撮影全体のビュー数がVの場合、Ω＝｛1、2、3、・・・、・・・、C×R×V−2、C×R×V−1、C×R×V｝となる。
ここで、SmがΩの真部分集合であるとは、Sm⊆Ω(SmはΩの部分集合)、かつ、Sm≠Ωが成り立つことである。更に、以下では、任意の真部分集合の組(Si、Sj)について、Si∩Sj＝φが成り立つものとして説明する。
まず、ステップ・アンド・シュートスキャンの場合について説明する。ステップ・アンド・シュートスキャンの場合、寝台4の位置が固定されている間は、同一の部位が撮影されていると仮定する。そして、演算装置5は、寝台4の位置が固定されている間に収集された投影データが、同一の真部分集合に属するように、各真部分集合の集合要素を決定する。
ステップ・アンド・シュートスキャンの場合、全体集合Ωの分割数、すなわち真部分集合Smの数Mは、寝台4の位置を移動させる回数に等しい。
ここで、寝台4の各位置において撮影されるビュー数V1、・・・、Vm、・・・、VMとする。この場合、撮影全体のビュー数V＝V1＋・・・＋Vm＋・・・＋VMである。また、各真部分集合Smの集合要素の数は、寝台4の各位置において撮影されるビュー数Vm、検出器12のチャネル数C、検出器12の列数Rを乗算した数Vm×C×Rとなる。
従って、各真部分集合｛S1、・・・、Sm、・・・、SM｝は、S1＝｛1、2、・・・、V1×C×R｝、S2＝｛V1×C×R＋1、・・・、V1×C×R＋V2×C×R｝、・・・、Sm＝｛(V1＋V2＋・・・＋Vm−1)×C×R＋1、・・・、(V1＋V2＋・・・＋Vm−1＋Vm)×C×R｝、・・・、 SM＝｛(V1＋V2＋・・・＋VM−1)×C×R＋1、・・・、(V1＋V2＋・・・＋VM−1＋VM)×C×R｝となる。
例えば、寝台4の各位置において、常に1周回(360°)分のビューが撮影され、1周回(＝360°)当たりのビュー数をV360とする。この場合、V1＝・・・＝Vm＝・・・＝VM＝V360、各真部分集合Smの集合要素の数＝V360×C×Rとなり、全ての真部分集合Smの集合要素の数が同一となる。
次に、らせんスキャンの場合について説明する。らせんスキャンの場合、演算装置5は、逐次近似法における逆投影処理の計算単位として定められるビューの数、すなわち逆投影ビュー数に基づいて、真部分集合Smに含まれる集合要素の数を決定する。これは、演算装置5は、1枚の断面像を作成するために必要な投影データを集合要素として持つように、各真部分集合Smを決定する必要があるからである。
例えば、非特許文献5に記載の方法によって逐次近似処理を行う場合、逆投影ビュー数は、逆投影処理の位相幅tw(以下、「逆投影位相幅tw」と呼ぶ。)に、撮影条件として設定される周回当たりの撮影ビュー数Vdを乗じた値となる。
逆投影位相幅twは、1周回(＝360°)分の位相を示す値が1、半周回(＝180°)分の位相を示す値が1/2、2周回(＝720°)分の位相を示す値が2、などになる。
ここで、逆投影位相幅twは任意のパラメータであり、非特許文献5では、寝台送り速度及び再構成画像サイズに応じて経験的に決定している。
そこで、演算装置5は、寝台送り速度及び再構成画像サイズごとに、逆投影位相幅twを予め記憶装置8に記憶しておくことが望ましい。そして、演算装置5は、記憶装置8から、撮影条件として設定される寝台送り速度、及び再構成条件として設定される再構成画像サイズに対応する単一の逆投影位相幅twを取得する。次に、演算装置5は、撮影条件として設定される周回当たりの撮影ビュー数Vdと、記憶装置8から取得される逆投影位相幅twとを乗算した値を逆投影ビュー数とする。
もし操作者が計算速度を優先する場合には、逆投影ビュー数を1/2周回(逆投影の最小単位)と置き換えた上で、真部分集合を決定しても良い。尚、逆投影処理は、この置き換えによらず、再構成条件の値を適用する。
らせんスキャンの場合、全体集合Ωの分割数、すなわち真部分集合Smの数Mは、撮影全体のビュー数Vを、逆投影ビュー数tw×Vdによって除したときの商となる。撮影全体のビュー数Vが、逆投影ビュー数tw×Vdによって割りきれない場合、余りのビューは、最後の真部分集合SMに含める。
各真部分集合Smの集合要素の数は、全て同一であり、逆投影位相幅tw、周回当たりの撮影ビュー数Vd、チャネル数C、列数Rを乗算した数tw×Vd×C×Rとなる。
従って、各真部分集合｛S1、・・・、Sm、・・・、SM｝は、S1＝｛1、2、・・・、tw×Vd×C×R｝、S2＝｛tw×Vd×C×R＋1、・・・、2×tw×Vd×C×R｝、・・・、Sm＝｛(m−1)×tw×Vd×C×R＋1、・・・、m×tw×Vd×C×R｝、・・・、SM＝｛(M−1)×tw×Vd×C×R＋1、・・・、M×tw×Vd×C×R｝となる。
図5の説明に戻る。次に、演算装置5は、ステップ2において決定された真部分集合Sm毎に、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdiに基づいて、罰則項に係る重み(罰則項重み)｛β1、・・・、βm、・・・、βM｝を算出する(ステップ3)。罰則項重み算出処理の詳細は、図6〜図8を参照しながら後述する。
次に、演算装置5は、i∈Smについて、真部分集合Sm毎の罰則項重みβmを用いて逐次近似法を実行する(ステップ4)。すなわち、演算装置5は、部位毎に罰則項重みβmを変更しながら、逐次近似投影データ補正処理を実行する。次式は、評価関数の一例である。
Figure 0005978429
 ここで、本発明の技術的思想は、評価関数における罰則項重みを部位毎に変更することであるため、検出器出力重み及び罰則項を含む評価関数であれば、最適化に使用する数値解析法に因らずに適用することができる。
式(3)から逐次近似法における更新式を導出することは、公知の手法によって実現することができる。例えば、非特許文献6のように、Gauss-Seidel法を用いることによって更新式を導出することができる。
演算装置5は、式(3)から逐次近似法における更新式を導出し、算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、投影データの補正処理を行う。
以降では、図6〜図8を参照しながら、2種類の罰則項重み算出処理の詳細について説明する。第1の罰則項重み算出処理については、図6、図7を参照しながら説明する。第2の罰則項重み算出処理については、図8を参照しながら説明する。
操作者は、所望計算時間や罰則項重みの推定精度を考慮して、両者のうちいずれか一方を任意に選択することができる。
最初に、第1の罰則項重み算出処理について説明する。
第1の罰則項重み算出処理では、演算装置5は、ノイズ低減率ごとに、基準ファントムに対してノイズ低減率と同程度のノイズ低減効果が達成される罰則項重みの値(以下、「罰則項重み基準値」と呼ぶ。)を予め記憶装置8に記憶する。
例えば、演算装置5は、ノイズ低減率10％、20％、30％、・・・に対して、罰則項重み基準値β＊ 10、β＊ 20、β＊ 30、・・・の値が登録されているノイズ低減率テーブルを記憶装置8に記憶しておく。罰則項重み基準値β＊ 10、β＊ 20、β＊ 30、・・・の値は、予め円形の水ファントムや人体を模擬したファントムを基準ファントムとし、基準ファントムを実際に撮影し、画像再構成処理を実行して得られた値とする。
そして、演算装置5は、記憶装置8から、入力されるノイズ低減率に対応する罰則項重み基準値を取得する。次に、演算装置5は、真部分集合Sm毎に、真部分集合Sm毎内に含まれる集合要素に対応する検出器出力重みdiに基づいて代表値を算出する。次に、演算装置5は、真部分集合Sm毎の代表値と、記憶装置8から取得される罰則項重み基準値とを乗算した値を、真部分集合Sm毎の罰則項重みβmとする。
特に、演算装置5は、真部分集合Sm毎に、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdiの平均値及び中央値、並びに、検出器出力重みdiを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級の3つの値の中でいずれか1つを、真部分集合Sm毎の代表値とすることが望ましい。
図6に示すフローチャートでは、代表値として、検出器出力重みdiを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級を算出する例を示している。
また、図6に示すフローチャートでは、真部分集合Smに対する罰則項重みβmの算出処理が示されているが、他の真部分集合についても同様である。
図6に示すように、演算装置5は、記憶装置8に記憶されているノイズ低減率テーブルから、操作者が入力した所望ノイズ低減率に対応するテーブル値を選択する(ステップ11)。
次に、演算装置5は、操作者が入力した所望計算時間に応じて、真部分集合Smに含まれる集合要素から、真部分集合Smの下位集合Tmを作成する(ステップ12)。
演算装置5は、予め定められた間引きルールに基づいて、Smに含まれる集合要素を間引くことによって、下位集合Tmを作成する。
間引きルールとしては、例えば、ビューをとびとびに間引くことが考えられる。ビューをとびとびに間引く場合には、最終的に作成される断面像の画質に与える影響が少ないからである。例えば、周回方向に1°ずつビューを取得している場合、周回当たりの撮影ビュー数Vdは、360個となる。360個をとびとびに間引くことによって、180個のビューだけが残る為、下位集合Tmの集合要素の数は、真部分集合Smの集合要素の数の1/2となる。
また、間引きルールとしては、例えば、各ビューの左端及び右端のチャネルを間引くことが考えられる。各ビューの左端及び右端のチャネルは、被検体を通過しないX線を受光する場合が多い為、データを間引いたとしても、最終的に作成される断面像の画質に与える影響が少ないからである。例えば、チャネル数Cに対して、左端から10％及び右端から10％のデータを間引くことによって、下位集合Tmの集合要素の数は、真部分集合Smの集合要素の数の4/5となる。
これらの間引きルールは、組み合わせて適用することも可能である。
いずれの間引きルール(又は複数の間引きルールの組み合わせ)であっても、操作者が入力した所望計算時間に応じて間引き量を加減することが可能である。このように、真部分集合Smの下位集合Tmを作成し、下位集合Tmに対して後述する処理を実行することによって、操作者の意図に合わせて計算時間を調整することができる。
次に、演算装置5は、ステップ12において作成された下位集合Tmの集合要素に対応する検出器出力重みdiを階級として、ヒストグラムを構築する(ステップ13)。ここで、検出器出力重みdiの算出手法は、公知のあらゆる手法を使用することができる。演算装置5は、例えば、非特許文献6のように、投影データを逆対数変換した後、エアデータを乗算することによって算出するようにしても良い。
図7には、ヒストグラムの一例が示されている。図7に示すヒストグラムは、横軸(階級)が検出器出力重みdi、縦軸(頻度)が各階級の値に等しい検出器出力重みdiに対応する集合要素の数である。
ヒストグラムは、被写体の大きさや部位によって異なる分布となり、それらの特徴をよく表す指標と言える。
肩部51のヒストグラムは、階級が0に近い値のところに最大のピークが存在する。また、肩部51のヒストグラムは、全体的に、検出器出力重みdiの値が低い。これは、肩部51が、X線の減弱量が大きい骨部を多く含むからである。
また、胸部52及び腹部53のヒストグラムは、階級が1に近い値のところに最大のピークが存在する。胸部52及び腹部53のヒストグラムを比較すると、階級が1〜3の頻度は、胸部52のヒストグラム＞腹部53のヒストグラムとなり、階級が5以上の頻度は、胸部52のヒストグラム＜腹部53のヒストグラムとなっている。これは、胸部52が、X線の減弱量が小さい肺を含むからである。
尚、被写体の大きさが異なれば、最大のピークが存在する階級が異なるので、図7に示すヒストグラムは、被写体の大きさの特徴もよく表していると言える。
図6の説明に戻る。次に、演算装置5は、図7に示す各ヒストグラムの面積を算出し、ヒストグラムの面積を所定の割合に分割する階級を代表値として特定する(ステップ14)。この他、代表値は、前述した通り、i∈Smに含まれる検出器出力重みdiの中央値や平均値としても良い。操作者は、これらを任意に選択することができる。
次に、演算装置5は、ステップ11において選択されたテーブル値と、ステップ14において特定された階級(代表値)を乗算することによって、真部分集合Smの罰則項重みβmを算出する(ステップ15)。
第1の罰則項重み算出処理では複雑な計算を行わないことから、計算量を大幅に増やすことなく、真部分集合Smの罰則項重みβmを算出することができる。
次に、第2の罰則項重み算出処理について説明する。第2の罰則項重み算出処理は、以下の2つの知見に基づく。
(1)検出器出力重みを定数とし、逐次近似法を任意の罰則項重みによって実行した場合、被写体及び部位によらず、ほぼ同一のノイズ低減率を得られることが経験的に明らかである。
(2)逐次近似法による投影値(画素値)の修正量とノイズ低減率との間には、高い相関がある。
尚、これらの知見は、逐次近似投影データ補正処理及び逐次近似再構成処理の両方について言える。
前述の知見を生かして、第2の罰則項重み算出処理では、演算装置5は、逐次近似法を実行した場合の投影データの修正量を算出する修正量算出関数に対して、真部分集合Sm内に含まれる集合要素に対応する検出器出力重みdi及び投影データを代入する。ここで、修正量算出関数には、真部分集合Sm毎の罰則項重みβmが変数として含まれる。そして、演算装置5は、修正量算出関数の値(＝投影データの修正量の総和)と所定の修正量基準値(＝投影データの修正量の総和の基準値)の誤差が最小となるように、真部分集合毎Smの罰則項重みβmを決定する。
所定の修正量基準値の決定処理は、次の通りである。演算装置5は、第1の罰則項重み算出処理と同様、ノイズ低減率テーブルを記憶装置8に記憶しておき、記憶装置8から、入力されるノイズ低減率に対応する罰則項重み基準値を取得する。そして、演算装置5は、検出器出力重みを定数とし、記憶装置8から取得される罰則項重み基準値を用いて逐次近似法を実行した場合の投影データの修正量を算出し、所定の修正量基準値とする。
以下、図8を参照しながら、第2の罰則項重み算出処理の詳細を説明する。図8に示すフローチャートでは、真部分集合Smに対する罰則項重みβmの算出処理が示されているが、他の真部分集合についても同様である。
図8に示すように、演算装置5は、記憶装置8に記憶されているノイズ低減率テーブルから、操作者が入力した所望ノイズ低減率に対応するテーブル値を選択する(ステップ21)。尚、処理内容の違いから、第1の罰則項重み算出処理に用いるノイズ低減率テーブルと、第2の罰則項重み算出処理に用いるノイズ低減率テーブルとは、記憶されているテーブル値が異なる。
次に、演算装置5は、操作者が入力した所望計算時間に応じて、真部分集合Smに含まれる集合要素から、真部分集合Smの下位集合Tmを作成する(ステップ22)。下位集合Tmの作成処理は、ステップ12と同様である。
次に、演算装置5は、ステップ21において選択されたテーブル値、及び、ステップ22において作成された下位集合Tmの集合要素に対応する投影データを用いて、検出器出力重みdiを任意の定数とした上で、逐次近似投影データ補正処理の適用後の投影値を推定する。ここで、計算量及び計算に要するメモリの観点から、適用後の投影値を解析的に算出することが難しいため、近似値を代用する。注目する投影値が、近接するチャネル、列、ビューのみから算出されると仮定すると、近似値は容易に算出できる。また、公知の方法を用いて、反復後の投影値の算出に必要となる逆行列を近似逆行列に置き換えた上で、近似値を算出しても良い。そして、演算装置5は、近似値と、逐次近似投影データ補正処理の適用前の投影値との誤差を算出し、誤差の総和を修正量基準値とする(ステップ23)。ここで、誤差の算出処理には、絶対誤差や二乗誤差等公知の指標を用いることができる。
次に、演算装置5は、罰則項重みβmを変数とした修正量算出関数に、ステップ23において作成された下位集合Tmの集合要素に対応する検出器出力重みdi及び投影データを代入する。そして、演算装置5は、修正量算出関数の値が、ステップ23において算出された修正量基準値と等しくなるように、罰則項重みを決定する(ステップ24)。ここで、修正量算出関数の値と修正量基準値との誤差の最小化には、公知の数値解析法(例えば二分法)を適用することができる。
第2の罰則項重み算出処理は前述の知見に基づくものであることから、精度良く真部分集合Smの罰則項重みβmを算出することができる。
以上のように、第1の実施の形態では、逐次近似投影データ補正処理に用いられる罰則項重みβmを、真部分集合Sm内では一定値とすることによって、逐次近似投影データ補正処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
また、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
＜第2の実施の形態＞
第2の実施の形態について説明する。第2の実施の形態では、本発明を逐次近似再構成処理に適用する場合について説明する。つまり、第2の実施の形態では、演算装置5が、検出器出力重み及び罰則項を評価関数に含む逐次近似法によって、再構成画像の再構成処理を行う。
第1の実施の形態と第2の実施の形態の違いは、第1の実施の形態における式(3)が、次式になることのみである。
Figure 0005978429
 式(4)についても、式(3)と同様、式(4)から逐次近似法における更新式を導出することができる。演算装置5は、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、再構成画像の再構成処理を行う。
以上のように、第2の実施の形態では、逐次近似再構成処理に用いられる罰則項重みβmを、第1の実施の形態と同様、真部分集合Sm内では一定値とすることによって、逐次近似再構成処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
また、第1の実施の形態と同様、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
＜第3の実施の形態＞
第3の実施の形態について説明する。第3の実施の形態では、本発明を逐次近似投影データ補正処理及び逐次近似再構成処理の両方に適用する場合について説明する。つまり、第3の実施の形態では、演算装置5が、検出器出力重み及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び再構成画像の再構成処理を行う。
第3の実施の形態は、第1の実施の形態と第2の実施の形態の組み合わせとなる。つまり、演算装置5は、第1の実施の形態と同様、式(3)から逐次近似法における更新式を導出し、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、投影データの補正処理を行う。また、演算装置5は、第2の実施の形態と同様、式(4)から逐次近似法における更新式を導出し、第1の実施の形態と同様に算出された検出器出力重みdi及び真部分集合Sm毎の罰則項重みβmを、導出された更新式に代入し、逐次近似法によって、再構成画像の再構成処理を行う。
以上のように、第3の実施の形態では、逐次近似投影データ補正処理及び逐次近似再構成処理に用いられる罰則項重みβmを、真部分集合Sm内では一定値とすることによって、逐次近似投影データ補正処理及び逐次近似再構成処理において、従来と同様に、CT画像のノイズ低減効果を達成することができる。
また、真部分集合Smは、被写体の情報が反映された検出器出力重みdi及び投影値に基づいて真部分集合Sm毎に罰則項重みβmを算出することによって、部位間でのノイズ低減効果のばらつきを抑制することができる。
以上、添付図面を参照しながら、本発明に係る医用画像処理装置等の好適な実施形態について説明したが、本発明はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
1 X線CT装置、 4 寝台、 5 演算装置、 6 入力装置、 7 表示装置、 8 記憶装置、 11 X線管、 12 検出器、 41 断面像、 42 チャネル、 43冠状面像、 44 列、 45スキャン方向

Claims (8)

  1. 検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び／又は再構成画像の再構成処理を行う医用画像処理装置であって、
    前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定する真部分集合決定部と、
    前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出する罰則項重み算出部と、
    前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行する逐次近似法実行部と、
    を備える医用画像処理装置。
  2. 前記真部分集合決定部は、らせんスキャンの場合、前記逐次近似法における逆投影処理の計算単位として定められる前記ビューの数である逆投影ビュー数に基づいて、前記真部分集合に含まれる前記集合要素の数を決定する
    請求項1に記載の医用画像処理装置。
  3. 前記真部分集合決定部は、
    寝台送り速度及び前記再構成画像の大きさごとに、前記逐次近似法における逆投影処理の位相幅である逆投影位相幅を予め記憶部に記憶し、
    前記記憶部から、前記撮影条件として設定される寝台送り速度、及び前記再構成条件として設定される前記再構成画像の大きさに対応する前記逆投影位相幅を取得し、
    前記撮影条件として設定される周回当たりの撮影ビュー数と、前記記憶部から取得される前記逆投影位相幅とを乗算した値を前記逆投影ビュー数とする
    請求項2に記載の医用画像処理装置。
  4. 前記罰則項重み算出部は、
    ノイズ低減率ごとに、基準ファントムに対して前記ノイズ低減率と同程度のノイズ低減が達成される前記罰則項重みである罰則項重み基準値を予め記憶部に記憶し、
    前記記憶部から、入力される前記ノイズ低減率に対応する前記罰則項重み基準値を取得し、
    前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて代表値を算出し、
    前記真部分集合毎の前記代表値と、前記記憶部から取得される前記罰則項重み基準値とを乗算した値を、前記真部分集合毎の前記罰則項重みとする
    請求項1に記載の医用画像処理装置。
  5. 前記罰則項重み算出部は、前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みの平均値及び中央値、並びに、前記検出器出力重みを階級とするヒストグラムにおいて全体を所定の割合によって分割する階級の3つの値の中でいずれか1つを、前記真部分集合毎の前記代表値とする
    請求項4に記載の医用画像処理装置。
  6. 前記罰則項重み算出部は、前記真部分集合毎の前記罰則項重みが変数であって、前記逐次近似法を実行した場合の前記投影データの修正量を算出する修正量算出関数に対して、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重み及び前記投影データを代入して算出される前記修正量算出関数の値と修正量基準値との誤差が最小となるように、前記真部分集合毎の前記罰則項重みを決定する
    請求項1に記載の医用画像処理装置。
  7. 前記罰則項重み算出部は、
    ノイズ低減率ごとに、基準ファントムに対して前記ノイズ低減率と同程度のノイズ低減効果が達成される前記罰則項重みである罰則項重み基準値を予め記憶部に記憶し、
    前記記憶部から、入力される前記ノイズ低減率に対応する前記罰則項重み基準値を取得し、
    前記検出器出力重みを定数とし、前記記憶部から取得される前記罰則項重み基準値を用いて前記逐次近似法を実行した場合の前記投影データの修正量を算出し、前記修正量基準値とする
    請求項6に記載の医用画像処理装置。
  8. 検出器の出力に応じた重みである検出器出力重み、及び罰則項を評価関数に含む逐次近似法によって、投影データの補正処理及び／又は再構成画像の再構成処理を行う医用画像処理方法であって、
    前記検出器のチャネル、前記検出器の列、及び前記投影データの取得単位であるビューの組み合わせを一意に識別する情報を集合要素とする全体集合から、撮影条件及び再構成条件に基づいて1又は複数の真部分集合を決定するステップと、
    前記真部分集合毎に、前記真部分集合内に含まれる前記集合要素に対応する前記検出器出力重みに基づいて、前記罰則項に係る重みである罰則項重みを算出するステップと、
    前記真部分集合毎の前記罰則項重みを用いて前記逐次近似法を実行するステップと、
    を含む医用画像処理方法。
JP2013511984A 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法 Active JP5978429B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2011100571 2011-04-28
JP2011100571 2011-04-28
PCT/JP2012/059136 WO2012147471A1 (ja) 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法

Publications (2)

Publication Number Publication Date
JPWO2012147471A1 JPWO2012147471A1 (ja) 2014-07-28
JP5978429B2 true JP5978429B2 (ja) 2016-08-24

Family

ID=47071996

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013511984A Active JP5978429B2 (ja) 2011-04-28 2012-04-04 医用画像処理装置、医用画像処理方法

Country Status (4)

Country Link
US (1) US9123098B2 (ja)
JP (1) JP5978429B2 (ja)
CN (1) CN103501702B (ja)
WO (1) WO2012147471A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11145094B2 (en) 2018-11-12 2021-10-12 Hitachi, Ltd. Image reconstruction apparatus and image reconstruction method
US11348226B2 (en) * 2018-06-07 2022-05-31 Canon Medical Systems Corporation Medical image diagnostic apparatus

Families Citing this family (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012147471A1 (ja) * 2011-04-28 2012-11-01 株式会社 日立メディコ 医用画像処理装置、医用画像処理方法
US9159146B2 (en) * 2011-07-08 2015-10-13 Hitachi Medical Corporation Image reconstruction device and image reconstruction method configured to perform iteratively reconstructed image using weight coefficient
US9406121B2 (en) * 2012-04-24 2016-08-02 Hitachi Medical Corporation X-ray CT apparatus and image reconstruction method
JP6312401B2 (ja) * 2012-11-30 2018-04-18 キヤノン株式会社 画像処理装置、画像処理方法、及びプログラム
JP6218334B2 (ja) * 2012-11-30 2017-10-25 株式会社日立製作所 Ｘ線ｃｔ装置及びその断層画像撮影方法
JP6329490B2 (ja) * 2013-02-05 2018-05-23 株式会社日立製作所 Ｘ線ｃｔ装置及び画像再構成方法
JP6492005B2 (ja) * 2013-04-08 2019-03-27 株式会社日立製作所 Ｘ線ｃｔ装置、再構成演算装置、及び再構成演算方法
US9895128B2 (en) * 2013-08-08 2018-02-20 Hitachi, Ltd. X-ray CT apparatus and correction processing device
US9974495B2 (en) 2014-01-20 2018-05-22 Hitachi, Ltd. X-ray CT apparatus, image processing device, and image reconstruction method
WO2015137011A1 (ja) * 2014-03-14 2015-09-17 株式会社日立メディコ Ｘ線ｃｔ装置、及び処理装置
CN104318536B (zh) * 2014-10-21 2018-03-20 沈阳东软医疗系统有限公司 Ct图像的校正方法及装置
JPWO2016132880A1 (ja) * 2015-02-16 2017-11-30 株式会社日立製作所 演算装置、ｘ線ｃｔ装置、及び画像再構成方法
JP6327188B2 (ja) * 2015-03-30 2018-05-23 三菱電機株式会社 窓関数決定装置、パルス圧縮装置、レーダ信号解析装置、レーダ装置、窓関数決定方法およびプログラム
JP6470837B2 (ja) * 2015-06-12 2019-02-13 株式会社日立製作所 Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
JP2017122705A (ja) * 2016-01-06 2017-07-13 三菱電機株式会社 算出方法、判定方法、選別方法および選別装置
CN106994021B (zh) * 2016-01-22 2022-10-14 通用电气公司 一种计算ct影像上的噪声的方法及装置
WO2017154217A1 (ja) * 2016-03-11 2017-09-14 株式会社島津製作所 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置
US10650621B1 (en) 2016-09-13 2020-05-12 Iocurrents, Inc. Interfacing with a vehicular controller area network
US10789738B2 (en) * 2017-11-03 2020-09-29 The University Of Chicago Method and apparatus to reduce artifacts in a computed-tomography (CT) image by iterative reconstruction (IR) using a cost function with a de-emphasis operator
US10679385B1 (en) 2018-12-17 2020-06-09 General Electric Company System and method for statistical iterative reconstruction and material decomposition
CN110070588B (zh) 2019-04-24 2023-01-31 上海联影医疗科技股份有限公司 Pet图像重建方法、系统、可读存储介质和设备

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6017568A (ja) * 1983-07-11 1985-01-29 Hitachi Ltd 画像処理方法および装置
WO2005077278A1 (ja) * 2004-02-16 2005-08-25 Hitachi Medical Corporation 断層撮影像の再構成方法及び断層撮影装置
JP2009540966A (ja) * 2006-06-22 2009-11-26 ゼネラル・エレクトリック・カンパニイ 画像の分解能を高めるシステム及び方法
WO2010062956A1 (en) * 2008-11-26 2010-06-03 Wisconsin Alumni Research Foundation Method for prior image constrained image reconstruction in cardiac cone beam computed tomography
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品
US20110052023A1 (en) * 2009-08-28 2011-03-03 International Business Machines Corporation Reconstruction of Images Using Sparse Representation

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5909476A (en) * 1997-09-22 1999-06-01 University Of Iowa Research Foundation Iterative process for reconstructing cone-beam tomographic images
US8204173B2 (en) * 2003-04-25 2012-06-19 Rapiscan Systems, Inc. System and method for image reconstruction by using multi-sheet surface rebinning
US7596204B2 (en) * 2005-03-17 2009-09-29 Koninklijke Philips Electronics N.V. Method and device for the iterative reconstruction of cardiac images
US8737711B2 (en) * 2008-08-07 2014-05-27 Hitachi Medical Corporation X-ray CT image forming method and X-ray CT apparatus using the same
CN101901472B (zh) * 2010-07-07 2012-12-19 清华大学 一种基于矩阵秩最小化的非刚性鲁棒批量图像对齐方法
JP5828841B2 (ja) * 2010-10-14 2015-12-09 株式会社日立メディコ Ｘ線ｃｔ装置及び画像再構成方法
WO2012147471A1 (ja) * 2011-04-28 2012-11-01 株式会社 日立メディコ 医用画像処理装置、医用画像処理方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS6017568A (ja) * 1983-07-11 1985-01-29 Hitachi Ltd 画像処理方法および装置
WO2005077278A1 (ja) * 2004-02-16 2005-08-25 Hitachi Medical Corporation 断層撮影像の再構成方法及び断層撮影装置
JP2009540966A (ja) * 2006-06-22 2009-11-26 ゼネラル・エレクトリック・カンパニイ 画像の分解能を高めるシステム及び方法
WO2010062956A1 (en) * 2008-11-26 2010-06-03 Wisconsin Alumni Research Foundation Method for prior image constrained image reconstruction in cardiac cone beam computed tomography
JP2010136958A (ja) * 2008-12-13 2010-06-24 Univ Of Tokushima Ｃｔ装置、ｃｔ装置における画像再構成方法、及び電子回路部品
US20110052023A1 (en) * 2009-08-28 2011-03-03 International Business Machines Corporation Reconstruction of Images Using Sparse Representation

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
JPN6012034527; 後藤大雅、他: '逐次近似再構成アルゴリズムの開発' 日本放射線技術学会総会学術大会予稿集 第６７巻, 20110225, 第１８４頁 *
JPN6012034528; Ｊ．Ｗａｎｇ，　ｅｔ．ａｌ: 'Ｍｕｌｔｉｓｃａｌｅ　Ｐｅｎａｌｉｚｅｄ　Ｗｅｉｇｈｔｅｄ　Ｌｅａｓｔ-Ｓｑｕａｒｅｓ　Ｓｉｎｏｇｒａｍ　Ｒｅｓｔｏｒａｔｉｏｎ　ｆｏｒ　Ｌｏｗ-Ｄｏｓｅ　Ｘ-Ｒａｙ　Ｃｏｍｐｕｔｅｄ　Ｔｏｍｏｇｒａ' Ｃｏｎｆ．Ｐｒｏｃ．ＩＥＥＥ　Ｅｎｇ．Ｍｅｄ．Ｂｉｏｌ．Ｓｏｃ． Ｖｏｌ．５５，　Ｉｓｓｕｅ３, 200803, ｐｐ．１０２２-１０３１, ＩＥＥＥ *
JPN6012034529; Ｈ．Ｒ．Ｓｈｉ，　ｅｔ．ａｌ．: 'Ｑｕａｄｒａｔｉｃ　Ｒｅｇｕｌａｒｉｚａｔｉｏｎ　Ｄｅｓｉｇｎ　ｆｏｒ　２-Ｄ　ＣＴ' ＩＥＥＥ　Ｔｒａｎｓ．Ｍｅｄ．Ｉｍｇ． Ｖｏｌ．２８，Ｎｏ．５, 200905, ｐｐ．６４５-６５６, ＩＥＥＥ *
JPN7012002539; Ｚ．Ｔｉａｎ，　ｅｔ．ａｌ．: 'Ｌｏｗ-ｄｏｓｅ　ＣＴ　ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ　ｖｉａ　ｅｄｇｅ-ｐｒｅｓｅｒｖｉｎｇ　ｔｏｔａｌ　ｖａｒｉａｔｉｏｎ　ｒｅｇｕｌａｒｉｚａｔｉｏｎ．' Ｐｈｙｓ．Ｍｅｄ．Ｂｉｏｌ． Ｖｏｌ．５６，Ｎｏ．１８, 201109, ｐｐ．５９４９-５９６７ *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11348226B2 (en) * 2018-06-07 2022-05-31 Canon Medical Systems Corporation Medical image diagnostic apparatus
US11145094B2 (en) 2018-11-12 2021-10-12 Hitachi, Ltd. Image reconstruction apparatus and image reconstruction method

Also Published As

Publication number Publication date
WO2012147471A1 (ja) 2012-11-01
US9123098B2 (en) 2015-09-01
CN103501702B (zh) 2015-09-02
JPWO2012147471A1 (ja) 2014-07-28
US20140226887A1 (en) 2014-08-14
CN103501702A (zh) 2014-01-08

Similar Documents

Publication Publication Date Title
JP5978429B2 (ja) 医用画像処理装置、医用画像処理方法
JP6937157B2 (ja) 放射線画像診断装置及び医用画像処理装置
JP5530637B2 (ja) 画像再構成の方法及びシステム
Liu et al. Total variation-stokes strategy for sparse-view X-ray CT image reconstruction
JP5828841B2 (ja) Ｘ線ｃｔ装置及び画像再構成方法
RU2471204C2 (ru) Локальная позитронная эмиссионная томография
JP5968316B2 (ja) 画像再構成装置及び画像再構成方法
JP6492005B2 (ja) Ｘ線ｃｔ装置、再構成演算装置、及び再構成演算方法
JP6470837B2 (ja) Ｘ線ｃｔ装置および逐次修正パラメータ決定方法
WO2017096609A1 (en) System and method for image reconstruction
JP6329490B2 (ja) Ｘ線ｃｔ装置及び画像再構成方法
WO2014041889A1 (ja) Ｘ線ｃｔ装置およびｘ線ｃｔ画像の処理方法
JP2016152916A (ja) Ｘ線コンピュータ断層撮像装置及び医用画像処理装置
US10813616B2 (en) Variance reduction for monte carlo-based scatter estimation
JP6124868B2 (ja) 画像処理装置及び画像処理方法
JP5637768B2 (ja) コンピュータ断層撮影画像の生成方法およびコンピュータ断層撮影装置
JP2018110867A (ja) 医用画像生成装置及び医用画像生成方法
CN114387359A (zh) 一种三维x射线低剂量成像方法及装置
US9858688B2 (en) Methods and systems for computed tomography motion compensation
WO2016132880A1 (ja) 演算装置、ｘ線ｃｔ装置、及び画像再構成方法
US20220375038A1 (en) Systems and methods for computed tomography image denoising with a bias-reducing loss function
US9508164B2 (en) Fast iterative image reconstruction method for 3D computed tomography
US11353411B2 (en) Methods and systems for multi-material decomposition

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20150402

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20160330

TRDD Decision of grant or rejection written
A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A712

Effective date: 20160427

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20160509

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160530

R150 Certificate of patent or registration of utility model

Ref document number: 5978429

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

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