JP6300782B2 - 画像処理装置、磁気共鳴イメージング装置および画像処理方法 - Google Patents

画像処理装置、磁気共鳴イメージング装置および画像処理方法 Download PDF

Info

Publication number
JP6300782B2
JP6300782B2 JP2015502869A JP2015502869A JP6300782B2 JP 6300782 B2 JP6300782 B2 JP 6300782B2 JP 2015502869 A JP2015502869 A JP 2015502869A JP 2015502869 A JP2015502869 A JP 2015502869A JP 6300782 B2 JP6300782 B2 JP 6300782B2
Authority
JP
Japan
Prior art keywords
image
wavelength band
denoising
original
unit
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
JP2015502869A
Other languages
English (en)
Other versions
JPWO2014132830A1 (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 JPWO2014132830A1 publication Critical patent/JPWO2014132830A1/ja
Application granted granted Critical
Publication of JP6300782B2 publication Critical patent/JP6300782B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G06T5/70—
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7246—Details of waveform analysis using correlation, e.g. template matching or determination of similarity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformation in the plane of the image
    • G06T3/40—Scaling the whole image or part thereof
    • G06T3/4038—Scaling the whole image or part thereof for image mosaicing, i.e. plane images composed of plane sub-images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration by the use of more than one image, e.g. averaging, subtraction
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves  involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • 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/10088—Magnetic resonance imaging [MRI]
    • 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/20024—Filtering details
    • 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/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging
    • 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

Description

本発明は、被検体中の水素や燐等からの核磁気共鳴信号(以下、NMR信号と呼ぶ)を測定し、核の密度分布や緩和時間分布等を画像化する磁気共鳴イメージング(以下、MRIと呼ぶ)装置に関し、特に画像のフィルタリング技術に関する。
MRI装置は、被検体、特に人体の組織を構成する原子核スピンが発生するNMR信号を計測し、その頭部、腹部、四肢等の形態等を2次元的に或いは3次元的に画像化する装置である。このMRI装置において計測されるNMR信号には熱雑音と呼ばれる電子回路中の電子の不規則な動きによるノイズが含まれる。そのため、得られた画像にはノイズが重畳し、画像のSN比を低下させる。
この問題を解決する1つの方法として、フィルタリングが挙げられる。フィルタリングは計測後に行う後処理であるため、計測時間の延長を伴わないというメリットがある。しかしながら、線形フィルタではノイズだけでなく、必要とされる信号の情報も失われるため、画像のボケ、すなわち分解能の低下が問題となる。
これを避けるため、画像に含まれる必要な情報の損失を防ぐフィルタリング手法が提案されており、例えば次のような手法がある。一つは、注目画素とその周辺のテクスチャの方向性を検出するか、エッジなどの構造を検出し、その結果に基づいてフィルタの重み関数の形状を変化させることで、画像上の構造を保持する手法である(例えば、非特許文献1参照)。このとき、方向性や構造の検出は、画像の画素値を基準に行われる。
H.Takeda et.al., Kernel Regression for Image Processing and Reconstruction, IEEE Trans.Img.Process, 2008 vol.16, no.2 p349 David A. Feinberg et.al., Halving MR Imaging Time by Conjugation: Demonstration at 3.5kG, Radiology, Vol.161, no.2, p527-531, (1986)
しかしながら、非特許文献1に開示のテクスチャの方向性や構造を検出する方法では、画素値の変動をテクスチャや構造と判断するため、ノイズによる画素値の変動もまた保存される。また、画像から構造を検出するため、ノイズも構造と判断されることがある。このため、フィルタリング後に不自然な構造を生じることがある。
本発明は、上記事情に鑑みてなされたもので、ノイズによって引き起こされる構造の判断誤りを低減し、画像に重畳したノイズを、有意な情報を保持しながら、除去後の結果画像に不自然さを生じさせずに除去する技術を提供することを目的とする。
本発明は、複数の原画像から作成した参照画像と各々の原画像とを比較して類似度を算出し、ノイズ判定の指標とする。当該指標を用いて各原画像を平滑化後合成し、最終的なノイズ除去後の画像を得る。このとき、原画像を、ノイズの最大波長に基づいて、複数の波長帯域毎に分割し、上記参照画像の作成、平滑化を行ってもよい。
ノイズによって引き起こされる構造の判断誤りを低減し、画像に重畳したノイズを、有意な情報を保持しながら、除去後の結果画像に不自然さを生じさせずに除去できる。
第一の実施形態のMRI装置の全体構成図 第一の実施形態のMRI装置の制御処理系の機能ブロック図 第一の実施形態の画像処理部の機能ブロックと各機能で生成されるデータを説明するための説明図 第一の実施形態のフィルタリング処理のフローチャート 第一の実施形態の帯域分割処理のフローチャート 第一の実施形態の参照画像作成処理のフローチャート 第一の実施形態の類似度算出処理のフローチャート 第一の実施形態のデノイズ処理のフローチャート 第一の実施形態の類似度に基づいたフィルタの重み関数の変化を示すグラフ 第一の実施形態の画像合成処理のフローチャート (a)は、第一実施形態のフィルタリング処理の結果と単純加算結果とを比較した画像のプロファイルグラフであり、(b)は、(a)の一部の拡大図 第一の実施形態のGUIの一例を説明するための説明図 第二の実施形態の画像処理部の機能ブロックと各機能で生成されるデータを説明するための説明図 第二の実施形態のフィルタリング処理のフローチャート 第三の実施形態の画像処理部の機能ブロックと各機能で生成されるデータを説明するための説明図 第三の実施形態のフィルタリング処理のフローチャート 第二の実施形態と第三の実施形態とを組み合わせたフィルタリング処理のフローチャート 第四の実施形態における原画像群の作成処理を説明するための説明図
＜＜第一の実施形態＞＞
以下、添付図面に従って本発明に係る好ましい実施形態について詳説する。なお、発明の実施形態を説明するための全図において、特に明示しない限り、同一機能を有するものは同一符号を付け、その繰り返しの説明は省略する。
本実施形態では、画像取得装置としてMRI装置を用いる場合を例にあげて説明する。
まず、本実施形態のMRI装置の一例の全体概要を図1に基づいて説明する。図1は、第一の実施形態のMRI装置の全体構成を示すブロック図である。
本実施形態のMRI装置100は、NMR現象を利用して被検体101の断層画像を得るもので、図1に示すように、静磁場発生系120と、傾斜磁場発生系130と、高周波磁場発生系(以下、送信系)150と、高周波磁場検出系(以下、受信系)160と、制御処理系170と、シーケンサ140と、を備える。
静磁場発生系120は、垂直磁場方式であれば、被検体101の周りの空間にその体軸と直交する方向に、水平磁場方式であれば、体軸方向に、均一な静磁場を発生させるもので、被検体101の周りに配置される永久磁石方式、常電導方式あるいは超電導方式の静磁場発生源を備える。
傾斜磁場発生系130は、MRI装置100の座標系(装置座標系)であるX、Y、Zの3軸方向に巻かれた傾斜磁場コイル131と、それぞれの傾斜磁場コイルを駆動する傾斜磁場電源132とを備え、後述のシ−ケンサ140からの命令に従ってそれぞれの傾斜磁場コイル131の傾斜磁場電源132を駆動することにより、X、Y、Zの3軸方向に傾斜磁場Gx、Gy、Gzを印加する。
送信系150は、被検体101の生体組織を構成する原子の原子核スピンに核磁気共鳴を起こさせるために、被検体101に高周波磁場パルス(以下、「RFパルス」と呼ぶ。)を照射するもので、高周波発振器(シンセサイザ)152と変調器153と高周波増幅器154と送信側の高周波コイル(送信コイル)151とを備える。高周波発振器152はRFパルスを生成し、シーケンサ140からの指令によるタイミングで出力する。変調器153は、出力されたRFパルスを振幅変調し、高周波増幅器154は、この振幅変調されたRFパルスを増幅し、被検体101に近接して配置された送信コイル151に供給する。送信コイル151は供給されたRFパルスを被検体101に照射する。
受信系160は、被検体101の生体組織を構成する原子核スピンの核磁気共鳴により放出される核磁気共鳴信号(エコー信号、NMR信号)を検出するもので、受信側の高周波コイル(受信コイル)161と信号増幅器162と直交位相検波器163と、A/D変換器164とを備える。受信コイル161は、被検体101に近接して配置され、送信コイル151から照射された電磁波によって誘起された被検体101の応答のNMR信号を検出する。検出されたNMR信号は、信号増幅器162で増幅された後、シーケンサ140からの指令によるタイミングで直交位相検波器163により直交する二系統の信号に分割され、それぞれがA/D変換器164でディジタル量に変換されて、制御処理系170に送られる。
シーケンサ140は、制御処理系170からの指示に従って、RFパルスと傾斜磁場パルスとを印加する。具体的には、制御処理系170からの指示に従って、被検体101の断層画像のデータ収集に必要な種々の命令を送信系150、傾斜磁場発生系130、および受信系160に送信する。
制御処理系170は、MRI装置100全体の制御、各種データ処理等の演算、処理結果の表示及び保存等を行うもので、CPU171と記憶装置172と表示装置173と入力装置174とを備える。記憶装置172は、ハードディスクなどの内部記憶装置と、外付けハードディスク、光ディスク、磁気ディスクなどの外部記憶装置とにより構成される。表示装置173は、CRT、液晶などのディスプレイ装置である。入力装置174は、MRI装置100の各種制御情報や制御処理系170で行う処理の制御情報の入力のインタフェースであり、例えば、トラックボールまたはマウスとキーボードとを備える。入力装置174は、表示装置173に近接して配置される。操作者は、表示装置173を見ながら入力装置174を通してインタラクティブにMRI装置100の各種処理に必要な指示、データを入力する。
CPU171は、操作者が入力した指示に従って、記憶装置172に予め保持されるプログラムを実行することにより、MRI装置100の動作の制御、各種データ処理等の制御処理系170の各処理を実現する。例えば、受信系160からのデータが制御処理系170に入力されると、CPU171は、信号処理、画像再構成処理等を実行し、その結果である被検体101の断層像を表示装置173に表示するとともに、記憶装置172に記憶する。
送信コイル151と傾斜磁場コイル131とは、被検体101が挿入される静磁場発生系120の静磁場空間内に、垂直磁場方式であれば被検体101に対向して、水平磁場方式であれば被検体101を取り囲むようにして設置される。また、受信コイル161は、被検体101に対向して、或いは取り囲むように設置される。
現在、MRI装置の撮像対象核種で、臨床で普及しているものは、被検体101の主たる構成物質である水素原子核(プロトン)である。MRI装置100では、プロトン密度の空間分布や、励起状態の緩和時間の空間分布に関する情報を画像化することで、人体頭部、腹部、四肢等の形態または機能を、二次元もしくは三次元的に撮像する。
本実施形態の制御処理系170は、図2に示すように、パルスシーケンスに従って、静磁場発生系120、傾斜磁場発生系130、高周波磁場発生系150および高周波磁場検出系160を備える計測系の動作を制御して計測した核磁気共鳴信号をk空間に配置し、ローデータ(Raw data)を得る撮像部210と、前記ローデータから画像を再構成する画像再構成部220と、画像に対し、画像処理を施す画像処理部(画像処理装置)230と、の機能を実現する。
本実施形態の画像処理部230は、同じ撮像範囲を撮像した複数の画像を用い、ノイズを低減するフィルタリング処理を行う。画像処理部230がフィルタリング処理に用いる複数の画像を、「複数の原画像」または、「原画像群」と呼び、原画像群の中の1枚の画像を指す場合は、単に「原画像」と呼ぶ。また、フィルタリング処理後の画像を、「合成画像」と呼ぶ。
具体的には、本実施形態の画像処理部230は、複数の原画像を合成して参照画像を作成し、前記複数の原画像それぞれについて、作成した前記参照画像と比較して類似度を算出し、算出した前記類似度に基づいて平滑化し、平滑化後の前記複数の原画像を合成して合成画像を得る。
なお、原画像群は、画像再構成部220が、撮像部210が取得したローデータ群から再構成することにより得る。このローデータ群は、撮像部210が、例えば、複数の受信コイル161を用いた計測により得る。また、同じ撮像条件で複数回の計測を行って得てもよい。同じ撮像範囲を撮像した複数の画像を原画像に用いるのは、異なるランダムノイズが重畳した画像を用いることで、ノイズの抽出精度を高めるためである。
本実施形態の画像処理部230は、上記フィルタリング処理を実現するため、図3に示すように、複数の原画像410(原画像群)を、それぞれ、予め定めた波長帯域毎に分割し、複数の波長帯域原画像420(波長帯域原画像群)を生成する帯域分割部310と、複数の波長帯域原画像420を、波長帯域毎に合成して波長帯域毎の参照画像430を作成する参照画像作成部320と、複数の波長帯域原画像420を、それぞれ同一の波長帯域の参照画像430と比較して、波長帯域原画像420毎の類似度を算出し、類似度マップ440を生成する類似度算出部330と、複数の波長帯域原画像420をそれぞれ当該波長帯域原画像420の類似度を用いて平滑化することにより、各波長帯域原画像420からそれぞれデノイズ画像450(複数のデノイズ画像；デノイズ画像群)を生成するデノイズ部340と、複数のデノイズ画像450を合成して合成画像460を生成する合成部350と、を備える。
各部が実行する処理の間に生成されるデータ、画像群は、それぞれ、記憶装置172の例えば、RAMに格納される。
本実施形態の上記フィルタリング処理の流れを説明する。図4は、本実施形態のフィルタリング処理の流れを説明するための処理フローである。ここでは、原画像410はN枚、また、帯域分割部310は、各原画像410をM個の波長帯域に分割するものとする。
N、Mは、1以上の整数とする。
まず、帯域分割部310がN枚の原画像410それぞれに対し、帯域分割処理を行う(ステップS1001)。これにより、N×M枚の波長帯域原画像420を得る。
次に、参照画像作成部320が、参照画像作成処理を行い、波長帯域毎に参照画像430を得る(ステップS1002)。参照画像430は、帯域毎に1枚作成する。このため、ここでは、M枚の参照画像430を得る。
次に、類似度算出部330が、類似度算出処理を行い、類似度を算出する(ステップS1003)。本実施形態では、N×M枚の波長帯域原画像420それぞれについて、画素毎に、x方向およびy方向の2方向について類似度を算出し、類似度マップ440とする。従って、2×N×M個の類似度マップ440が得られる。
そして、デノイズ部340がデノイズ処理を行い、デノイズ画像450を生成する(ステップS1004)。本実施形態では、N×M枚の波長帯域原画像420それぞれについて、画素毎に、x方向およびy方向の2方向についてデノイズ処理を行う。従って、N×M枚のデノイズ画像を得る。
最後に合成部350が画像合成処理を行い、合成画像460を得る(ステップS1005)。ここでは、デノイズ処理を終えたN×M枚のデノイズ画像を合成し、1枚の合成画像460を生成する。
以下、各処理について詳細に説明する。
＜帯域分割処理＞
帯域分割部310は、上述のように、原画像群410を、予め定めた波長帯域毎に分割する帯域分割処理を行う。本実施形態では、フィルタにより、予め定めたカットオフ波長より大きい帯域を有する画像群とカットオフ波長以下の帯域を有する画像群との2つに分割する場合を例にあげて説明する。本実施形態の帯域分割部は、これを、同一のカットオフ波長を設定したローパスフィルタおよびハイパスフィルタを用いて実現する。なお、分割する波長帯域は、2つに限らない。
図5は、本実施形態の帯域分割処理の処理フローである。なお、各原画像410には、それぞれ、予め原画像番号が付与されているものとする。ここで、原画像番号nの原画像をIOrg(n、x、y)と表す。
帯域分割部310は、ローパスフィルタおよびハイパスフィルタのカットオフ波長を設定する(ステップS1101)。
そして、以下のステップS1103からステップS1104の処理を、原画像番号nを1から1ずつインクリメントし、複数の原画像(原画像群：IOrg(n、x、y))それぞれについて実行する(ステップS1102、S1107)。
帯域分割部310は、n番目の原画像(IOrg(n、x、y))を読み込む(ステップS1103)。そして、原画像群IOrg(n、x、y)に対し、ローパスフィルタによるローパスフィルタ処理を行う(ステップS1104)。また、原画像群に対し、ハイパスフィルタによるハイパスフィルタ処理を行う(ステップS1105)。両フィルタ処理は、いずれを先に行ってもよい。
そして、帯域分割部310は、得られたローパスフィルタ結果およびハイパスフィルタ結果を、それぞれ、波長帯域原画像420として、記憶装置172に記憶する(ステップS1106)。このとき、波長帯域分割後の各帯域に帯域番号を付与する。帯域番号mの波長帯域原画像をISep(m、n、x、y)と表す。本実施形態では、ローパスフィルタ結果には、帯域番号1を付与し、第一の波長帯域原画像ISep(1、n、x、y)と呼ぶ。また、ハイパスフィルタ結果には、帯域番号2を付与し、第二の波長帯域原画像ISep(2、n、x、y)と呼ぶ。
ここで、ステップS1101のカットオフ波長の設定の詳細について説明する。設定するカットオフ波長は、ノイズの波長に応じて決定する。すなわち、ノイズの最大波長λNoiseをJ(pixels)とすると、カットオフ波長は、1／J(1／pixels)とする。例えば、ノイズの最大波長λNoiseを3pixelsとすると、カットオフ波長は1／3(1／pixels)とする。一般に、ノイズの最大波長λNoiseは、3(pixels)が適当であるが、補間された画像を原画像に用いる場合、補間倍率に比例して点数(カットオフ波長〔pixels〕)を増やす。
なお、ローパスフィルタには、ガウシアンフィルタを用いてもよい。急峻な遮断特性を持つフィルタよりも、ガウシアンのようななだらかな遮断特性を有するフィルタを用いる方が、画像にリンギングが現れにくく、好ましい。この場合、カットオフ波長λcutoffは、以下のように設定する。
ガウシアンフィルタの振幅伝達特性H(λ)は、以下の式(1)で表される。
Figure 0006300782
 なお、αは、以下の式(2)で表される。
Figure 0006300782
 ここで、λは、波長(pixels)である。式(1)および式(2)から、ノイズの最大波長λNoise以下の信号成分の減衰率(％)をδとすると、カットオフ波長λcutoffは、以下の式(3)により計算できる。
Figure 0006300782
 この場合、例えば、ノイズの最大波長λNoiseを、上述のように3(pixels)、減衰率δを99(％)とすると、カットオフ波長λcutoffは、約7．7(pixels)となる。
また、ローパスフィルタに、ガウシアンフィルタを用いる場合、上記ステップS1104のローパスフィルタ処理は、以下の式(4)に従う。
Figure 0006300782
 ここで、h(τ)は、ガウシアンフィルタの重み関数である。
また、ステップS1105のハイパスフィルタ処理は、以下の式(5)に従う。
Figure 0006300782
 なお、ノイズの最大波長λNoiseは、ユーザが後述のGUIを介して設定する。また、予め記憶装置172に設定しておいてもよい。
＜参照画像作成処理＞
次に、参照画像作成部320による参照画像作成処理の詳細を説明する。参照画像作成処理では、波長帯域毎に、暫定的に正しいものとみなす参照画像を作成する。本実施形態では、各原画像410から作成した波長帯域毎の波長帯域原画像420を用い、参照画像430を作成する。図6は、本実施形態の参照画像作成処理の処理フローである。
参照画像作成部320は、以下のステップS1202からステップS1206の処理を、分割した波長帯域数M(Mは、2以上の整数)だけ、繰り返す(ステップS1201、S1207)。本実施形態では、カットオフ波長を用い、ローパスフィルタおよびハイパスフィルタで2つの波長帯域に分割しているため、2回繰り返す。m＝1の場合は、ローパスフィルタ結果である第一の波長帯域原画像ISep(1、n、x、y)を、m＝2の場合は、ハイパスフィルタ結果である第二の波長帯域原画像ISep(2、n、x、y)を処理する。
参照画像作成部320は、N枚の原画像410から生成された帯域番号mの波長帯域原画像420を全て読み込む(ステップS1202、S1203、S1204)。
そして、参照画像作成部320は、読み込んだ波長帯域原画像420を合成して、波長帯域毎に、1の参照画像430を得る(ステップS1205)。この時の合成処理は、例えば、以下の式(6)に従い、全ての原画像410から得た同じ波長帯域の波長帯域原画像420の平均値を算出する。
Figure 0006300782
 ここで、IRef(m、x、y)は、帯域番号mの参照画像430、Nは、原画像410の総数、ここでは、N枚の原画像410から得た、帯域番号mの波長帯域原画像数である。
参照画像作成部320は、得られた参照画像IRef(m、x、y)430を、記憶装置172に格納する(ステップS1206)。
以下、本実施形態では、第一の波長帯域原画像ISep(1、n、x、y)群から得た参照画像IRef(1、x、y)を第一の参照画像、第二の波長帯域原画像ISep(2、n、x、y)群から得た参照画像IRef(2、x、y)を第二の参照画像と呼ぶ。
なお、本実施形態では、平均値を算出することにより、波長帯域毎の参照画像を得ているが、算出手法はこれに限られない。例えば、各波長帯域原画像の2乗和の平方根をとり、参照画像を得てもよい。
＜類似度算出処理＞
次に、類似度算出部330による類似度算出処理について説明する。類似度算出処理では、参照画像430に対する類似度を、波長帯域原画像420毎に算出する。類似度の算出は、各画素について、x方向、y方向について行う。図7は、本実施形態の類似度算出処理の処理フローである。
類似度算出部330は、ステップS1302からステップS1308の処理を、分割した波長帯域数毎に、繰り返す(ステップS1301、S1309)。本実施形態では、参照画像作成処理同様、2回繰り返す。
まず、類似度算出部330は、m番目の参照画像430を読み込む(ステップS1302)。本実施形態では、m＝1の場合は、第一の参照画像IRef(1、x、y)を、m＝2の場合は、第二の参照画像IRef(2、x、y)を読み込む。
類似度算出部330は、ステップS1304からステップS1307の処理を、N枚の原画像410それぞれから得た、帯域番号mの波長帯域原画像420に対して実行する(ステップS1303、S1308)。
類似度算出部330は、まず、帯域番号mの波長帯域原画像420であって、n番目の原画像410から生成した波長帯域原画像ISep(m、n、x、y)を読み込む(ステップS1304)。
次に、類似度算出部330は、読み込んだ波長帯域原画像ISep(m、n、x、y)と参照画像IRef(m、x、y)との局所的な類似度を算出する。ここでは、まず、波長帯域原画像ISep(m、n、x、y)の各画素についてx方向の類似度を算出し、類似度マップ(x方向類似度マップ)を作成する(ステップS1305)。各画素の類似度には、例えば、相関係数を用いる。各画素のx方向の類似度(x方向の類似度マップ；Similarity X(m、n、x、y))は、例えば、以下の式(7)に従って、算出する。
Figure 0006300782
ここで、Lは、相関係数算出範囲である。相関係数算出範囲Lは、例えば事前に定義したノイズの最大波長と同じ値に設定する。上式の相関係数の算出は、類似度の表現を意図したものである。なお、本実施形態では、類似度は相関係数として算出されているため、類似度の最大値は1であり、最小値は−1である。
次に、類似度算出部330は、波長帯域原画像ISep(m、n、x、y)の各画素についてy方向の類似度を算出し、y方向の類似度マップ；Similarity Y(m、n、x、y)を作成する(ステップS1306)。y方向類似度マップSimilarity Y(m、n、x、y)は、例えば、以下の式(8)に従って算出される。
Figure 0006300782
そして、類似度算出部330は、各画素のx方向の類似度、y方向の類似度として、得られたx方向の類似度マップSimilarity X(m、n、x、y)およびy方向の類似度マップSimilarity Y(m、n、x、y)を、記憶装置172に格納する(ステップS1307)。なお、x方向、y方向を特に区別する必要がない場合は、類似度マップSimilarity(m、n、x、y)と呼ぶ。
なお、算出した類似度にノイズの影響が多く見られる場合は、類似度マップにメディアンフィルタなどのフィルタ処理を行ってもよい。フィルタ処理は、類似度は空間的に連続的に変化するという仮定の下に行う。
＜デノイズ処理＞
次に、デノイズ部340によるデノイズ処理について説明する。デノイズ処理では、波長帯域毎に分割した原画像410、すなわち、波長帯域原画像420それぞれの各画素のノイズを除去する。ノイズの除去は、画素ごとに、当該画素の類似度に応じて重み関数を作成し、それを用いて行う。重み関数は、類似度が高い画素ほど急峻な形状のものとする。これは、類似度が高い画素ほど、平滑化がされにくく、類似度が低い画素ほど平滑化されやすくするためである。
図8は、本実施形態のデノイズ処理の処理フローである。本実施形態では、各波長帯域原画像420の画素毎に、x方向、y方向それぞれに重み関数を作成し、フィルタリング処理を行う。ここでは、各画素(x、y)に、画素番号pが付与されているものとする。
全画素数は、P(Pは、1以上の整数)とする。
デノイズ部340は、ステップS1402からS1414の処理を、分割した波長帯域毎に繰り返す(ステップS1401、S1415)。本実施形態では、参照画像作成処理同様、2回繰り返す。
また、デノイズ部340は、ステップS1403からステップS1413の処理を、同一波長帯域の波長帯域原画像群420毎に繰り返す(ステップS1402、S1414)。ここでは、N回繰り返す。
また、デノイズ部340は、ステップS1404からステップS1406の処理を、全画素について順に行う(ステップS1403、S1407)。ここでは、P回繰り返す。
デノイズ部340は、画素p＝(x、y)のx方向の類似度マップSimilarity X(m，n，x、y)を読み込む(ステップS1404)。そして、読み込んだ類似度マップSimilarity X(m，n，x、y)を用い、画素pの重み関数を作成する(ステップS1405)。重み関数は、例えば、類似度マップSimilarity(m，n，x、y)およびカットオフ波長λcutoffを用い、以下の式(9)に従って作成する。
Figure 0006300782
図9は、式(9)に従って重み関数hDenoise(τ)を作成した一例である。実線501は、類似度Similarityが0以下の場合の重み関数、破線502は、類似度Similarityが0．5の場合の重み関数、点線503は、類似度Similarityが0．8の場合の重み関数である。本図に示すように、類似度が高いほど急峻な形状になる。このような重み関数を用いて、平滑化処理を行うことにより、類似度の高い画素は平滑化がされにくく、反対に類似度の低い画素は強く平滑化されることになる。
デノイズ部340は、算出した重み関数を用いて、x方向にフィルタリング処理(平滑化)を行う(ステップS1406)。
次に、デノイズ部340は、y方向について、同様の処理を行う。すなわち、画素p毎に、以下のステップS1409からS1411の処理を繰り返す(ステップS1408、S1412)。まず、画素p＝(x、y)のy方向の類似度マップSimilarity Y(m，n，x、y)を読み込む(ステップS1409)。そして、読み込んだ類似度マップSimilarity Y(m，n，x、y)を用い、重み関数を作成する(ステップS1410)。重み関数は、例えば、上記式(9)に従って作成する。そして、作成した重み関数を用いて、y方向にフィルタリング処理(平滑化)を行う(ステップS1411)。
そして、デノイズ部340は、得られた画像を、デノイズ画像deN(m、n)450として記憶装置172に格納する(ステップS1413)。
＜画像合成処理＞
次に、合成部350による、画像合成処理について説明する。図10は、本実施形態の画像合成処理の処理フローである。ここでは、同一の原画像410から生成されたデノイズ画像deN(m，n)450を合成し、N個の全波長帯域のデノイズ画像(合成デノイズ画像)deNal(n)を得る。そして、全ての合成デノイズ画像deNal(n)を合成し、最終的に1枚の合成画像460を得る。
このため、合成部350は、ステップS1502からS1505の処理を、デノイズ画像deN(m、n)450の元となった原画像410毎に実行する(ステップS1501、S1506)。
合成部350は、同一原画像410から生成された波長帯域原画像420のデノイズ画像deN(m，n)450を全て読み込み(ステップS1502、S1503、S1504)、合成し、全波長帯域のデノイズ画像deNal(n)を生成する(ステップS1505)。本実施形態では、deN(1、n)とdeN(2、n)とを読み込み、合成し、合成デノイズ画像deNal(n)を生成する。
全ての合成デノイズ画像の生成を終えると、合成部350は、得られた全ての合成デノイズ画像を合成し、合成画像460を得る(ステップS1507)。
なお、合成部350が実行する合成は、単純加算でもよいし、例えば、画像群が、複数の受信コイル161により取得されたものであれば、2乗和後に、平方根を取るなどの計算によってもよい。
ここで、複数の受信コイル161を用いて取得した原画像群410を、単純に加算して得た画像と、当該原画像群410に、本実施形態のフィルタリング処理を施して得た画像(合成画像460)とを比較する。図11(a)には、単純に加算して得た画像のプロファイル601と本実施形態によるフィルタリング処理により得た画像(合成画像460)のプロファイル602とを示す。本図において、薄いラインがプロファイル601、濃いラインがプロファイル602である。また、図11(b)に、図11(a)の一点鎖線部分603を拡大した図を示す。これらの図に示すように、原画像を単純加算した画像のプロファイル601に比べて、本実施形態によるフィルタリング処理により得た画像(合成画像460)のプロファイル602は、ノイズが低減されながらも構造を保持している。
最後に、本実施形態におけるGUI(Graphical User Interface)を説明する。本実施形態のフィルタ処理を実行するために操作者が指定すべき変数は、帯域分割部310が使用するカットオフ波長を決定するノイズの最大波長λNoiseである。
本実施形態のGUI700の一例を図12に示す。本図に示すように、本実施形態のGUI700は、ノイズの最大波長λNoiseを受け付けるボックス701と、確定の意思を受け付けるボタン702と、を備え、ノイズの最大波長λNoiseの設定を受け付ける設定受付部として機能する。帯域分割部310は、GUI700を介して、ノイズの最大波長λNoiseを受け付ける。
GUI700は、表示装置173に表示され、入力装置174を通して操作する。操作者は、ボックス701にノイズの最大波長を入力し、ボタン702を押して設定を完了する。帯域分割部310は、ボタン702の押下を受け付けると、その時点でボックス701に入力されている値を、ノイズの最大波長λNoiseとして受け付ける。
設定するノイズの最大波長λNoiseは、前述したように通常3［pixels］でよいが、補間された画像であれば、補間倍率に比例して点数を増やす。さらにノイズの最大波長λNoiseはフィルタ処理の平滑化の程度を調整する変数でもある。そのため、平滑化の程度を強くしたい場合には大きな数値を、平滑化の程度を弱くしたい場合には小さな数値を入力する。入力したノイズの最大波長以下の信号成分は除去されるため、操作者は処理後の画像に求める分解能に応じて値を設定すればよい。
なお、GUI700は、操作者によるノイズの最大波長の設定を意図したものである。
このため、構成は、上記に限られない。必ずしもテキストボックスを用いなくても、スライダバーなど他の入力方法が実現可能な構成であってもよい。
以上説明したように、本実施形態の画像処理部230は、複数の原画像410を、それぞれ、予め定めた波長帯域毎に分割し、複数の波長帯域原画像420を生成する帯域分割部310と、前記複数の波長帯域原画像420を、前記波長帯域毎に合成して前記波長帯域毎の参照画像430を作成する参照画像作成部320と、前記複数の波長帯域原画像420を、それぞれ同一の波長帯域の前記参照画像430と比較して、前記波長帯域原画像420毎の類似度マップ440を算出する類似度算出部330と、前記複数の波長帯域原画像420をそれぞれ当該波長帯域原画像420の前記類似度マップ440を用いて平滑化することにより、各波長帯域原画像420からそれぞれデノイズ画像450を生成するデノイズ部340と、前記デノイズ画像450を合成して合成画像460を生成する合成部350と、を備える。
例えば、非特許文献1等に記載されているような、注目画素とその周辺のテクスチャの方向性やエッジなどの構造を検出してフィルタの形状を変化させる方法では、単一の画像を用いて処理するため、ノイズによる画素値の変動が、テクスチャや構造として認識され保存される。しかしながら、本実施形態の手法では、参照画像430を、異なるランダムノイズが重畳した複数の画像(波長帯域原画像群420)を合成して作成する。そして、得られた参照画像430と各々の画像(波長帯域原画像群420)とを比較して算出される類似度(類似度マップ440)を構造の指標とする。従って、ランダムなノイズを構造として認識する誤りを低減することができる。
さらに、本実施形態によれば、原画像群410をノイズの最大波長から決定される波長帯域で分割する。このため、算出される類似度にノイズの波長帯域とは異なる波長帯域のデータの影響が入りにくい。これにより、ノイズの検出能を上げることができる。
また、注目画素と周囲点との画素値の違いを重みに反映させると、注目画素の画素値に近いデータだけを選んで平滑化が行われるため、フィルタリング結果の画像が平坦になりやすい。しかしながら、本実施形態の手法では、参照画像430との類似度に基づいてデノイズもしくはフィルタリングを行う。このため、先見的な仮定(同じ画素値レベルのデータを用いて平滑化を行えば構造を保存できるといった仮定)を必要としないため不自然な処理結果を生じさせにくい。しかも、算出された類似度に基づいてデノイズ処理を行うため、テクスチャの方向性やエッジを仮定する必要がない。従って、あらゆる形状を持つ画像において、良好なデノイズ処理結果を得ることができる。
加えて、本実施形態のフィルタ処理は、複雑な変数設定を必要としない。操作者は、単一の変数(ノイズの最大波長)を設定すればよいだけであるため、処理結果を想像しやすく、容易に使用することが可能である。
以上のように、本実施形態によれば、ノイズによって引き起こされる構造の判断誤りを低減し、画像に重畳したノイズを、有意な情報を保持しながら、除去後の結果画像に不自然さを生じさせずに除去できる。
なお、本実施形態では、帯域分割処理を行わなくてもよい。この場合、帯域分割部310は備えなくてもよい。
このような構成であっても、本実施形態によれば、原画像群410から作成した参照画像430と、原画像群410とを比較して類似度を算出し、算出した類似度をノイズ判定の指標とする。このため、画像に重畳したノイズを有意な信号として認識する判断誤りが少なくなる。また、算出された類似度を用いてデノイズ処理を行うため、テクスチャの方向性やエッジを仮定する必要がなくなる。このため、あらゆる形状を持つ画像において、良好なデノイズ処理結果を得ることができる。
なお、本実施形態では、MRI装置100のCPU171が、上記画像処理部230を実現する場合を例にあげて説明したが、これに限られない。画像処理部230の機能は、MRI装置100とデータの送受信可能な、他の情報処理装置上に構築されてもよい。
また、本実施形態では、画像処理部230が処理対象とする原画像410を取得する装置として、MRI装置100を用いる場合を例にあげて説明したが、画像を取得する装置は、これに限られない。他の医用画像取得装置であってもよい。例えば、超音波診断装置やCT(Computed Tomography)装置で同じ対象を繰り返し測定した画像群を原画像群として用いて同じ処理を実施することも可能である。さらに、原画像410は、医用画像取得装置で取得された医用画像でなくてもよい。
＜＜第二の実施形態＞＞
次に、本発明を適用する第二の実施形態を説明する。本実施形態では、第一の実施形態のフィルタリング処理を繰り返すことにより、柔軟かつ十分なノイズ除去効果を得る。
本実施形態のMRI装置100の構成は、基本的に第一の実施形態と同様である。但し、上述の繰り返し処理を行うため、本実施形態の画像処理部230の機能が異なる。以下、本実施形態について、第一の実施形態と異なる構成に主眼をおいて説明する。
本実施形態の画像処理部230は、図13に示すように、第一の実施形態の構成に加え、合成画像460が生成される毎に収束しているか否かを判別し、否と判別する場合、原画像410を、直前に生成された合成デノイズ画像に置き換える収束判定部360を備える。
収束判定部360は、1回前の繰り返しで得た合成画像460からの変化量で、合成画像460が収束したか否かの判定を行う。ここでは、変化量が所定より小さくなった場合、収束したと判定する。例えば、RepeatNum回目に得た合成画像IComp(RepeatNum、x、y)が、(RepeatNum−1)回目に得た合成画像IComp(RepeatNum−1、x、y)との間で、以下の式(10)を満たす場合、収束したものと判別する。
Figure 0006300782
ここで、Xはx方向の画像点数、Yはy方向の画像点数、Absは絶対値関数、RepeatNumは繰り返し回数である。式(10)は、画像全体の画素値の総和に対して繰り返し処理による変化量が1％以下であることを表す。なお、収束判定基準は、これに限られない。
収束と判定する変化量の最大値をユーザが設定可能なように構成してもよい。
なお、収束判定部360は、1回目に得た合成画像IComp(1、x、y)については、上記式(10)の計算ができないため、常に否と判定する。
以下、本実施形態の画像処理部230によるフィルタリング処理の流れを説明する。図14は、本実施形態のフィルタリング処理の流れを説明するための処理フローである。
第一の実施形態同様、帯域分割部310が帯域分割処理を行い(ステップS2101)、参照画像作成部320が参照画像作成処理を行い(ステップS2102)、類似度算出部330が類似度算出処理を行い(ステップS2103)、デノイズ部340がデノイズ処理を行い(ステップS2104)、合成部350が合成処理を行う(ステップS2105)。ただし、本実施形態では、合成部350は、合成デノイズ画像470も、記憶装置172に格納する。
合成画像460を得ると、収束判定部360は、例えば、上記式(10)を計算し、得られた合成画像460が収束しているか否かを判定する(ステップS2106)。収束していると判定されると、処理を終了する。一方、否と判定された場合、収束判定部360は、原画像410群を、直前に算出された合成デノイズ画像470群に置き換え(ステップS2107)、ステップS2101へ戻る。なお、画像の置き換えは、合成デノイズ画像生成後であればよく、収束判定の前であってもよい。
以上説明したように、本実施形態によれば、第一の実施形態同様、複数の原画像410を合成して参照画像430を作成し、複数の原画像410それぞれについて、作成した参照画像430と比較して類似度を算出し、算出した類似度をノイズ判定の指標とし、デノイズ処理を行う。このとき、参照画像の作成、類似度の算出、デノイズ処理を、予め定めた波長帯域毎に行ってもよい。従って、第一の実施形態同様、画像に重畳したノイズを有意な信号として認識する判断誤りが少なくなる。また、算出された類似度を用いてデノイズ処理を行うため、テクスチャの方向性やエッジを仮定する必要がなくなる。このため、あらゆる形状を持つ画像において、良好なデノイズ処理結果を得ることができる。
さらに、デノイズ後の画像を原画像とし、収束するまでデノイズ処理を繰り返すため、デノイズ処理のデノイズ効果(平滑化の程度)を調整する必要がない。従って、本実施形態によれば、処理パラメータの調整を行うことなく、さまざまな画像に適用することができる。
なお、本実施形態も、第一の実施形態同様、画像処理部230が処理対象とする原画像410を取得する装置は、MRI装置に限られない。他の医用画像取得装置であってもよい。また、医用画像取得装置でなくても、画像を取得可能な装置であればよい。
＜＜第三の実施形態＞＞
次に、本発明を適用する第三の実施形態を説明する。本実施形態では、帯域分割処理に用いるカットオフ波長を変化させながらフィルタリング処理を繰り返す。
本実施形態のMRI装置100の構成は、基本的に第一の実施形態と同様である。但し、上述の繰り返し処理を行うため、本実施形態の画像処理部230の機能が異なる。以下、本実施形態について、第一の実施形態と異なる構成に主眼をおいて説明する。
本実施形態の画像処理部230は、第一の実施形態の構成に加え、図15に示すように、波長制御部370を備える。また、本実施形態では、帯域分割部310は、1のカットオフ波長を用い、原画像410を2つの波長帯域原画像420に分割する。そして、波長帯域の小さい方の波長帯域原画像420についてのみデノイズ処理を行う。
本実施形態の波長制御部370は、予め定めた波長帯域内で、前記波長帯域内の最小値に設定された前記カットオフ波長を予め定めた増分ずつ増加させて更新することを繰り返すとともに、前記原画像410を置き換える。
NMR信号に含まれる熱雑音は、全波長帯域にほぼ均一な振幅スペクトルを持つため、ノイズの波長帯域は0から無限大である。しかしながら、MRI装置において特に除去したいノイズは、微細な構造を見えなくする原因となる短波長のノイズである。従って、カットオフ波長を変化させる波長帯域は、3〜12［pixels］程度と設定することが実用的であり、望ましい。なお、画像が高倍率に補間されている場合は、より広い波長帯域を設定することが望ましい。
この変化させる波長帯域の設定は、操作者が行うよう構成してもよいし、予め定め、記憶装置172に保持するよう構成してもよい。また、変化させる際の増分についても同様である。
以下、本実施形態の画像処理部230によるフィルタリング処理の流れを説明する。図16は、本実施形態のフィルタリング処理の流れを説明するための処理フローである。
まず、波長制御部370は、カットオフ波長を変化させる範囲として、ノイズの波長帯域を設定する(ステップS3101)。設定範囲は上記のとおりである。同時に増分も設定し、ループカウンタk、およびkの最大値Kを算出する。例えば、ノイズの波長帯域の最小値をλminとし、最大値をλmax、増分をΔλとする。カットオフ波長λcutoffは、ループカウンタkを用い、λcutoff＝λmin＋Δλ(k−1)で算出される。また、カウンタkの最大値Kは、K＝Int((λmax−λmin)/Δλ)＋1で算出される。なお、Int(x)は、xの整数部分を返す関数である。なお、カウンタがkの場合のカットオフ波長をλcutoff(k)と表す。
例えば、上記のように、ノイズの波長帯域として3〜12[pixels]が、増分Δλとして1pixelが設定された場合、Kは、10となる。
波長制御部370は、ループカウンタkを1からKまで1ずつ変化させて、ステップS3103からS3108の処理を繰り返す(ステップS3102、S3109)。
帯域分割部310は、カットオフ波長λcutoff(k)を用い、各原画像410に対し、帯域分割処理を行う(ステップS3103)。ここでは、カットオフ波長λcutoff(k)を用いて各原画像410を2つの波長帯域に分割する。そして、カットオフ波長λcutoff(k)以下の波長帯域の波長帯域原画像420を、第一の波長帯域原画像ISep(1、n、x、y)、他方を第二の波長帯域原画像ISep(2、n、x、y)とする。
次に、参照画像作成部320が参照画像作成処理を行う(ステップS3104)。本実施形態では、デノイズ処理を、カットオフ波長λcutoff(k)以下の波長帯域原画像420である第一の波長帯域原画像ISep(1、n、x、y)に対してのみ行うため、参照画像430についても、第一の波長帯域原画像ISep(1、n、x、y)からのみ作成する。すなわち、本実施形態の参照画像作成部320は、第一の波長帯域原画像ISep(1、n、x、y)を全て読み込み、合成し、第一の参照画像IRef(1、x、y)を作成する。
次に、類似度算出部330は、類似度算出処理を行う(ステップS3105)。本実施形態では、第一の波長帯域原画像ISep(1、n、x、y)の各画素について、第一の参照画像IRef(1、x、y)を用い、x方向およびy方向の類似度を算出し、第一の波長帯域原画像ISep(1、n、x、y)の類似度マップを作成する。
デノイズ部340は、第一の波長帯域原画像ISep(1、n、x、y)に対し、デノイズ処理を行う(ステップS3106)。デノイズ処理においても、第一の波長帯域原画像ISep(1、n、x、y)の各画素についてのみ行い、第一の波長帯域原画像ISep(1、n、x、y)のデノイズ画像(以下、第一のデノイズ画像と呼ぶ。)を生成する。
次に、合成部350は、画像合成処理を行う(ステップS3107)。本実施形態では、同一の原画像410から作成された、第一のデノイズ画像と第二の波長帯域原画像ISep(2、n、x、y)とを合成し、1の原画像の合成デノイズ画像470を得る。そして、全ての合成デノイズ画像470を合成し、合成画像460を得る。本実施形態では、合成部350は、合成デノイズ画像470も記憶装置172に格納する。
次に、波長制御部370は、原画像410群を、得られた合成デノイズ画像470群に置き換える(ステップS3108)。
なお、上記ループ処理において、ステップS3107では、合成デノイズ画像470のみを算出し、全ループ処理が終了後、合成デノイズ画像470を合成し、合成画像460を得るよう構成してもよい。
以上説明したように、本実施形態によれば、第一の実施形態同様、複数の原画像410から作成した参照画像430と、原画像410とを比較して類似度を算出し、算出した類似度をノイズ判定の指標とし、デノイズ処理を行う。このため、第一の実施形態と同様の効果を得ることができる。
さらに、本実施形態によれば、予め定めたノイズの波長帯域内でカットオフ波長を増加させながら上記処理を繰り返すため、厳密なノイズの最大波長の定義を必要としない。従って、厳密な処理パラメータの設定無しに、精度よくノイズを除去できる。
なお、本実施形態も、第一の実施形態同様、画像処理部230が処理対象とする原画像410を取得する装置は、MRI装置に限られない。他の医用画像取得装置であってもよい。また、医用画像取得装置でなくても、画像を取得可能な装置であればよい。
また、本実施形態のループ処理は、第二の実施形態の収束判定処理と組み合わせてもよい。この場合、合成画像460が生成される毎に収束しているか否かを判別し、否と判別する場合、原画像410を合成デノイズ画像470に置き変える収束判定部360をさらに備える。この場合の処理の流れを図17に示す。
本図に示すように、カットオフ波長毎に、本実施形態のステップS3103からステップS3108の処理を行い、第二の実施形態の収束判定処理(ステップS2106)を行う。収束判定処理において、否と判定された場合、また、ステップS3103からの処理を繰り返す。そして、収束判定処理において、収束したと判定された後、カットオフ波長を更新する(ステップS3102、3109)。
＜＜第四の実施形態＞＞
次に、本発明を適用する第四の実施形態を説明する。第一の実施形態、第二の実施形態では、予め複数の原画像を取得し、上記画像処理を行う。しかしながら、本実施形態では、MRI装置を用いて画像を1枚、生成可能なデータを取得し、そのデータから複数の原画像を作成する。
本実施形態のMRI装置100は、基本的に第一の実施形態のMRI装置100と同様である。ただし、本実施形態では、取得した1枚の画像用のデータから、複数の原画像を作成する。このため、本実施形態の画像再構成部220は、1枚の画像から複数の原画像を生成する。以下、本実施形態について、第一の実施形態と異なる構成に主眼をおいて説明する。
本実施形態の画像再構成部220は、図2に示すように、撮像部210が取得した1のローデータから欠損領域がそれぞれ異なる複数の欠損ローデータを生成する欠損データ生成部221と、各欠損ローデータから推定データを生成する推定データ生成部222と、各推定データからそれぞれ原画像を再構成する原画像生成部223と、を備える。
以下、本実施形態の画像再構成部220の各部による原画像410群の作成処理を、図18を用いて説明する。ローデータ801は、撮像部210が取得した、画像1枚分のNMR信号をk空間に配置したものである。
欠損データ生成部221は、取得したローデータ801から、欠損ローデータ802を複数生成する。欠損ローデータ802は、ローデータ801が配置されたk空間の一部の領域のデータをゼロにしたものである。このとき、欠損させる領域をそれぞれ変えて、複数の異なる欠損ローデータ802を生成する。
欠損させる領域は、全ローデータの半分以下の領域とする。また、中心付近を欠損させないことが望ましい。
推定データ生成部222は、公知のk空間推定技術を用いて、欠損ローデータ802群それぞれのゼロにした領域にデータを充填し、推定ローデータ803群を生成する。なお、k空間推定技術は、例えば、非特許文献2に記載の、非欠損領域のデータを用いて欠損領域のデータを充填する技術を用いる。
原画像生成部223は、それぞれの推定ローデータ803群をフーリエ変換し、画像を再構成することにより、原画像群804を生成する。
図18に示すように、本実施形態の原画像804群は、1枚の画像分のローデータ801から作成されたものであるが、各々異なる欠損領域が推定されているため、k空間上における異なる領域のノイズが強調され、異なるノイズを持つ。各原画像804は、異なるノイズを持つため、この原画像804群から作成する参照画像との比較で、ノイズと有意な信号を判別できる。
作成した原画像804を用いて行う以降の処理は、第一の実施形態、第二の実施形態および第三の実施形態のいずれかと同様である。
以上説明したように、本実施形態によれば、画像処理部230によるフィルタリング処理において、上記各実施形態同様の効果が得られる。さらに、本実施形態では、単一の受信コイル161で、または、1回の測定で、原画像群804を得ることができる。
なお、上記各実施形態では、例えば、波長帯域の分割に、式(4)、式(5)で説明したガウシアンフィルタを用いているが、これに限られない。例えば、スプラインフィルタや任意の振幅伝達特性を持ったフィルタを用いることができる。
また、上記各実施形態では、類似度の算出を、式(7)、式(8)を用いて行っているが、これに限られない。例えば、参照画像と波長帯域分割画像との偏差の2乗和などを用いてもよい。
また、上記各実施形態では、デノイズ処理において、式(9)に従って、フィルタリングを行う際に用いる重み関数hdenoise(τ)を算出しているが、これに限られない。例えば、メディアンフィルタや画像群の間の比較による外れ値処理等を行ってもよい。
また、上記各実施形態で処理対象とする原画像は、2次元画像に限定されるものではなく、3次元画像であってもよい。また、次元ごとに平滑化の程度を変化させてもよい。平滑化の程度は、ノイズの最大波長の設定を変更することにより変えることができる。さらに、平滑化の程度の指定には、ノイズの最大波長の定義を用いなくてもよい。例えば、類似度算出処理で算出した類似度Similarity X、Similarity Yに係数を乗算するなどすることで平滑化の程度を調整してもよい。
100 MRI装置、101 被検体、120 静磁場発生系、130 傾斜磁場発生系、131 傾斜磁場コイル、132 傾斜磁場電源、140 シーケンサ、150 高周波磁場発生系(送信系)、151 送信コイル、152 高周波発振器、153 変調器、154 高周波増幅器、160 高周波磁場検出系(受信系)、161 受信コイル、162 信号増幅器、163 直交位相検波器、164 A/D変換器、170 制御処理系、171 CPU、172 記憶装置、173 表示装置、174 入力装置、210 撮像部、220 画像再構成部、221 欠損データ生成部、222 推定データ生成部、223 原画像生成部、230 画像処理部、310 帯域分割部、320 参照画像作成部、330 類似度算出部、340 デノイズ部、350 合成部、360 収束判定部、370 波長制御部、410 原画像、420 波長帯域原画像、430 参照画像、440 類似度マップ、450 デノイズ画像、460 合成画像、470 合成デノイズ画像、501 実線、502 破線、503 点線、601 プロファイル、602 プロファイル、603 一点鎖線部分、700 GUI、701 ボックス、702 ボタン、801 ローデータ、802 欠損ローデータ、803 推定ローデータ、804 原画像

Claims (15)

  1. 複数の原画像を合成して参照画像を作成し、前記複数の原画像それぞれについて、作成した前記参照画像と比較して類似度を算出し、算出した前記類似度に基づいて前記複数の原画像を平滑化し、平滑化後の前記複数の原画像を合成して合成画像を得る画像処理部を備えること
    を特徴とする画像処理装置。
  2. 前記画像処理部は、
    前記複数の原画像それぞれを、予め定めた波長帯域毎に分割し、複数の波長帯域原画像を生成する帯域分割部と、
    前記複数の波長帯域原画像を、前記波長帯域毎に合成して前記波長帯域毎の前記参照画像を作成する参照画像作成部と、
    前記複数の波長帯域原画像それぞれを、同一の波長帯域の前記参照画像と比較して、前記波長帯域原画像毎の前記類似度を算出する類似度算出部と、
    前記複数の波長帯域原画像それぞれを、当該波長帯域原画像の前記類似度を用いて平滑化することにより、各前記波長帯域原画像からそれぞれデノイズ画像を生成するデノイズ部と、
    前記デノイズ画像を合成して前記合成画像を生成する合成部と、を備えること
    を特徴とする請求項1記載の画像処理装置。
  3. 前記帯域分割部は、予め定めたノイズの最大波長からカットオフ波長を算出し、当該カットオフ波長により、前記複数の原画像それぞれから、前記波長帯域原画像を生成すること を特徴とする請求項2記載の画像処理装置。
  4. 前記画像処理部は、前記合成画像が生成される毎に収束しているか否かを判別し、否と判別する場合、前記合成画像を前記原画像に置き変える収束判定部をさらに備え、
    前記帯域分割部は、前記合成画像が前記原画像に置き換えられる毎に、置き換え後の原画像を分割して新たな前記波長帯域原画像を生成し、
    前記参照画像作成部は、前記波長帯域原画像が生成される毎に前記参照画像を作成し、 前記類似度算出部は、前記参照画像が作成される毎に前記類似度を算出し、
    前記デノイズ部は、前記類似度が算出される毎に前記デノイズ画像を生成し、
    前記合成部は、前記デノイズ画像が生成される毎に、同一の前記原画像から生成されたデノイズ画像を合成して合成デノイズ画像を生成し、生成した全ての合成デノイズ画像を合成して前記合成画像を生成し、
    前記収束判定部は、否と判別する場合、前記原画像を直前に生成された前記合成デノイズ画像に置き換えること
    を特徴とする請求項2記載の画像処理装置。
  5. 前記画像処理部は、予め定めた波長帯域内で、前記波長帯域内の最小値に設定された前記カットオフ波長を予め定めた増分ずつ増加させて更新することを繰り返すとともに、前記原画像を置き換える波長制御部をさらに備え、
    前記帯域分割部は、前記カットオフ波長が設定または更新される毎に、当該カットオフ波長を用いて前記置き換え後の原画像を2つの波長帯域に分割し、前記カットオフ波長以下の波長帯域の第一の波長帯域原画像と当該カットオフ波長より大きい波長帯域の第二の波長帯域原画像とを生成し、
    前記参照画像作成部は、前記第一の波長帯域原画像が生成される毎に当該第一の波長帯域原画像から前記参照画像を作成し、
    前記類似度算出部は、前記参照画像が作成される毎に前記第一の波長帯域原画像の前記類似度を算出し、
    前記デノイズ部は、前記類似度が算出される毎に前記第一の波長帯域原画像の前記デノイズ画像を生成し、
    前記合成部は、前記デノイズ画像が生成される毎に、同一の前記原画像から生成された前記デノイズ画像と前記第二の波長帯域画像とを合成して合成デノイズ画像を生成し、生成した全ての合成デノイズ画像を合成して前記合成画像を生成し、
    前記波長制御部は、前記カットオフ波長を更新する際、前記原画像を、直前に生成された前記合成デノイズ画像に置き換えること
    を特徴とする請求項3記載の画像処理装置。
  6. 前記画像処理部は、前記合成画像が生成される毎に収束しているか否かを判別し、否と判別する場合、前記原画像を前記合成デノイズ画像に置き変える収束判定部をさらに備えること
    を特徴とする請求項5記載の画像処理装置。
  7. 前記画像処理部は、前記ノイズの最大波長の設定を受け付ける設定受付部をさらに備えること
    を特徴とする請求項3記載の画像処理装置。
  8. 前記類似度は、画素毎の、前記参照画像との局所的な類似度であって、当該類似度から算出される前記平滑化に用いる重み関数を、当該類似度が高いほど急峻な形状とするものであること
    を特徴とする請求項2記載の画像処理装置。
  9. パルスシーケンスに従って、静磁場発生系、傾斜磁場発生系、高周波磁場発生系および高周波磁場検出系を備える計測系の動作を制御して計測した核磁気共鳴信号をk空間に配置し、ローデータを得る撮像部と、
    1の前記ローデータから複数の原画像を再構成する画像再構成部と、
    請求項1記載の画像処理装置と、を備えること
    を特徴とする磁気共鳴イメージング装置。
  10. 前記画像再構成部は、
    前記1のローデータから欠損領域がそれぞれ異なる複数の欠損ローデータを生成する欠損データ生成部と、
    前記各欠損ローデータから推定データを生成する推定データ生成部と、
    各推定データから前記原画像を再構成する原画像生成部と、を備えること
    を特徴とする請求項9記載の磁気共鳴イメージング装置。
  11. パルスシーケンスに従って、静磁場発生系、傾斜磁場発生系、高周波磁場発生系および高周波磁場検出系を備える計測系の動作を制御して計測した核磁気共鳴信号をk空間に配置し、ローデータを得る撮像部と、
    1の前記ローデータから1の原画像を再構成する画像再構成部と、
    請求項1記載の画像処理装置と、を備え、
    前記高周波磁場検出系は、複数の受信コイルを備え、
    前記撮像部は、前記計測した核磁気共鳴信号を前記受信コイル毎に前記k空間に配置することにより、複数の前記ローデータを得、
    前記画像再構成部は、各前記ローデータからそれぞれ原画像を再構成することにより、前記複数の原画像を得ること
    を特徴とする磁気共鳴イメージング装置。
  12. パルスシーケンスに従って、静磁場発生系、傾斜磁場発生系、高周波磁場発生系および高周波磁場検出系を備える計測系の動作を制御して計測を行い、得られた核磁気共鳴信号をk空間に配置し、ローデータを得る撮像部と、
    1の前記ローデータから1の原画像を再構成する画像再構成部と、
    請求項1記載の画像処理装置と、を備え、
    前記撮像部は、同じ撮像条件で前記計測を繰り返すことにより複数の前記ローデータを得、
    前記画像再構成部は、各前記ローデータからそれぞれ原画像を再構成することにより、前記複数の原画像を得ること
    を特徴とする磁気共鳴イメージング装置。
  13. 複数の原画像を、それぞれ、波長帯域毎に分割し、複数の波長帯域原画像を生成する帯域分割処理と、
    前記複数の波長帯域原画像を、前記波長帯域毎に合成して前記波長帯域毎の参照画像を作成する参照画像作成処理と、
    前記複数の波長帯域原画像を、それぞれ同一の波長帯域の前記参照画像と比較して、前記波長帯域原画像毎の類似度を算出する類似度算出処理と、
    前記複数の波長帯域原画像をそれぞれ当該波長帯域原画像の前記類似度に基づいて平滑化することにより、各波長帯域原画像からそれぞれデノイズ画像を生成するデノイズ処理と、
    前記デノイズ画像を合成して合成画像を生成する合成処理と、を実行すること
    を特徴とする画像処理方法。
  14. 前記合成画像が生成される毎に収束しているか否かを判別し、否と判別する場合、前記原画像を置き換え、収束していると判定するまで前記帯域分割処理、前記参照画像作成処理、前記類似度算出処理、前記デノイズ処理、および前記合成処理を実行させる収束判定処理をさらに実行し、
    前記合成処理では、前記デノイズ画像が生成される毎に、同一の前記原画像から生成されたデノイズ画像を合成して合成デノイズ画像を生成し、生成した全ての合成デノイズ画像を合成して前記合成画像を生成し、
    前記収束判定処理では、前記原画像を前記合成デノイズ画像に置き換えること
    を特徴とする請求項13記載の画像処理方法。
  15. 予め定めた波長帯域内の最小値をカットオフ波長に設定する処理と、
    複数の原画像それぞれを、前記カットオフ波長を用いて2つの波長帯域毎に分割し、前記カットオフ波長以下の波長帯域の第一の波長帯域原画像と前記カットオフ波長より大きい波長帯域の第二の波長帯域原画像とを生成する帯域分割処理と、
    前記第一の波長帯域原画像を合成して参照画像を作成する参照画像作成処理と、
    前記第一の波長帯域原画像それぞれを、前記参照画像と比較して、類似度を算出する類似度算出処理と、
    前記第一の波長帯域原画像それぞれを、当該第一の波長帯域原画像の前記類似度に基づいて平滑化することにより、それぞれデノイズ画像を生成するデノイズ処理と、
    同一の原画像から生成された前記デノイズ画像と前記第二の波長帯域原画像と合成して前記原画像毎の合成デノイズ画像を生成し、生成した合成デノイズ画像を全て合成して合成画像を生成する合成処理と、
    前記原画像を前記合成デノイズ画像に置き換え、前記カットオフ波長を予め定めた増分だけ増加し、受け付けた範囲を超えるまで、前記帯域分割処理、参照画像作成処理、類似度算出処理、デノイズ処理および合成処理を繰り返す処理と、
    を実行することを特徴とする画像処理方法。
JP2015502869A 2013-02-28 2014-02-18 画像処理装置、磁気共鳴イメージング装置および画像処理方法 Active JP6300782B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2013038135 2013-02-28
JP2013038135 2013-02-28
PCT/JP2014/053679 WO2014132830A1 (ja) 2013-02-28 2014-02-18 画像処理装置、磁気共鳴イメージング装置および画像処理方法

Publications (2)

Publication Number Publication Date
JPWO2014132830A1 JPWO2014132830A1 (ja) 2017-02-02
JP6300782B2 true JP6300782B2 (ja) 2018-03-28

Family

ID=51428102

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2015502869A Active JP6300782B2 (ja) 2013-02-28 2014-02-18 画像処理装置、磁気共鳴イメージング装置および画像処理方法

Country Status (3)

Country Link
US (1) US9830687B2 (ja)
JP (1) JP6300782B2 (ja)
WO (1) WO2014132830A1 (ja)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9360543B2 (en) * 2012-04-19 2016-06-07 Regents Of The University Of Minnesota System and method for parallel magnetic resonance image reconstruction using digital beamforming
JPWO2017217395A1 (ja) * 2016-06-14 2019-04-11 日本電気株式会社 画像処理装置、画像処理方法及びプログラム
CN107918929B (zh) * 2016-10-08 2019-06-21 杭州海康威视数字技术股份有限公司 一种图像融合方法、装置及系统
WO2019075621A1 (zh) 2017-10-16 2019-04-25 北京深迈瑞医疗电子技术研究院有限公司 超声成像设备、系统及其超声造影成像的图像增强方法
CN107886475A (zh) * 2017-12-11 2018-04-06 奕响（大连）科技有限公司 一种单通道的图片相似判定方法
US10643313B2 (en) * 2018-01-19 2020-05-05 Bae Systems Information And Electronic Systems Integration Inc. Methods for image denoising and deblurring
CN108805840B (zh) * 2018-06-11 2021-03-26 Oppo（重庆）智能科技有限公司 图像去噪的方法、装置、终端及计算机可读存储介质
JP7186604B2 (ja) * 2018-12-25 2022-12-09 キヤノンメディカルシステムズ株式会社 医用画像処理装置、医用画像処理方法、およびプログラム
JP7302988B2 (ja) * 2019-03-07 2023-07-04 富士フイルムヘルスケア株式会社 医用撮像装置、医用画像処理装置、及び、医用画像処理プログラム
DE102019215460A1 (de) * 2019-10-09 2021-04-15 Siemens Healthcare Gmbh Verfahren und Vorrichtung zur Rauschreduktion von Bildaufnahmen
CN111049989B (zh) * 2019-12-26 2021-04-06 维沃移动通信有限公司 一种图像显示方法及电子设备
JP2023160661A (ja) * 2022-04-22 2023-11-02 キヤノン株式会社 画像処理装置、画像処理方法、及びコンピュータプログラム
CN114943639B (zh) * 2022-05-24 2023-03-28 北京瑞莱智慧科技有限公司 图像获取方法、相关装置及存储介质
CN116645661B (zh) * 2023-07-27 2023-11-14 深圳市青虹激光科技有限公司 一种防重码检测方法及系统

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3847512B2 (ja) * 2000-02-07 2006-11-22 株式会社日立メディコ 磁気共鳴イメージング装置
CN1305011C (zh) * 2001-04-19 2007-03-14 株式会社东芝 图像处理方法和图像处理设备
JP4907798B2 (ja) 2001-08-24 2012-04-04 株式会社東芝 超音波診断装置
JP4363833B2 (ja) * 2001-10-16 2009-11-11 株式会社東芝 局所血流動態に関するインデックスを演算する方法及び装置
JP4826332B2 (ja) 2006-05-11 2011-11-30 株式会社日立製作所 磁気共鳴測定装置
JP2010161521A (ja) * 2009-01-07 2010-07-22 Nec Corp 画像処理装置、撮像装置、画像ぶれの補正方法、及びプログラム
JP2010181951A (ja) 2009-02-03 2010-08-19 Mitsubishi Electric Corp 画像処理装置および画像処理プログラム
JP5917054B2 (ja) * 2011-09-12 2016-05-11 キヤノン株式会社 撮像装置、画像データ処理方法、およびプログラム

Also Published As

Publication number Publication date
US20160012569A1 (en) 2016-01-14
US9830687B2 (en) 2017-11-28
WO2014132830A1 (ja) 2014-09-04
JPWO2014132830A1 (ja) 2017-02-02

Similar Documents

Publication Publication Date Title
JP6300782B2 (ja) 画像処理装置、磁気共鳴イメージング装置および画像処理方法
CN111513716B (zh) 使用扩展灵敏度模型和深度神经网络进行磁共振图像重建的方法和系统
JP5902317B2 (ja) 磁気共鳴イメージング装置および定量的磁化率マッピング法
JP6289664B2 (ja) 磁気共鳴イメージング装置、定量的磁化率マッピング方法、計算機、磁化率分布計算方法、及び、磁化率分布計算プログラム
US7949170B2 (en) Image processing method, image processing device, computer aided detection, and method for filtering along the time axis
JP6227710B2 (ja) 磁気共鳴イメージング装置およびフリップ角決定方法
JP5843876B2 (ja) 磁気共鳴イメージング装置および磁化率強調画像生成方法
JP5815722B2 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
JP6744764B2 (ja) 画像診断装置、及び画像取得方法
WO2017221654A1 (ja) 磁気共鳴イメージング装置、画像処理装置、及び拡散強調画像計算方法
JP2018198682A (ja) 磁気共鳴イメージング装置及び磁気共鳴画像処理方法
US10132902B2 (en) Intrinsic navigation from velocity-encoding gradients in phase-contrast MRI
JPWO2016021603A1 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
WO2017038345A1 (ja) 磁気共鳴イメージング装置および撮像シーケンス生成方法
WO2016017385A1 (ja) 磁気共鳴イメージング装置および画像再構成方法
JP6757200B2 (ja) 画像診断装置、及び、磁気共鳴イメージング装置
JP4675936B2 (ja) 核磁気共鳴撮影装置
JP6783619B2 (ja) 磁気共鳴イメージング装置及び画像解析方法
WO2015170394A1 (ja) 撮像装置、画像処理装置及び画像処理方法
JP2015167567A (ja) 磁気共鳴イメージング装置
WO2016125572A1 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
JP2011031058A (ja) 核磁気共鳴撮影装置
JP2018068954A (ja) 磁気共鳴イメージング装置及び画像処理方法
JP2018079147A (ja) 磁気共鳴イメージング装置及び画像処理方法

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20170126

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20170126

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20171030

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7424

Effective date: 20171107

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: 20180130

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20180227

R150 Certificate of patent or registration of utility model

Ref document number: 6300782

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

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250