JP7608111B2

JP7608111B2 - 磁気共鳴イメージング装置、画像処理装置、及び、画像処理方法

Info

Publication number: JP7608111B2
Application number: JP2020177601A
Authority: JP
Inventors: 眞次黒川
Original assignee: Fujifilm Corp
Current assignee: Fujifilm Corp
Priority date: 2020-10-22
Filing date: 2020-10-22
Publication date: 2025-01-06
Anticipated expiration: 2040-10-22
Also published as: CN114376554B; CN114376554A; US20220128639A1; JP2022068752A; US12241956B2

Description

本発明は、磁気共鳴イメージング装置（以下、ＭＲＩ装置という）に関し、特に不十分なデータを用いて逐次近似により画像を再構成する技術に関する。

データが不足しているときに逐次近似で画像を再構成するという技術があり、ＭＲＩにおいては圧縮センシング（ＣＳ）などで利用されている（非特許文献１、特許文献１）。撮像時に取得するデータを少なくすることで、撮像時間は短縮できるが、画像の再構成にWavelet変換やCurvelet変換とその逆変換を伴う多くの繰り返し処理が必要になり、再構成時間が長くなるという課題がある。

繰り返し処理の高速化にはいくつか手法がある。例えば、繰り返し処理においてWavelet変換やCurvelet変換の複数の変換を用いても高速な処理となるAlternating Split Bregman法という優れた手法がある（非特許文献２）。また、繰り返しの初期値において密度補正を行って収束を速めるということも行われている。密度補正は、計測データの密度分布に応じて、データを重み付けする処理である。

米国特許公開２０１５／０１２６８５０号

Michael Lustig, David Donoho, and John M. Pauly. Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging. Magnetic Resonance in Medicine 58:1182-1195 (2007). Gerlind Plonka, Jianwei Ma. Curvelet-Wavelet Regularized Split Bregman Iteration for Compressed Sensing. International Journal of Wavelets, Multiresolution and Information Processing|Vol. 09, No. 01, pp. 79-110 (2011)

非特許文献２に記載されたAlternating Split Bregman法や、事前の密度補正により、繰り返し処理の効率はかなり改善されるものの、依然として多くの繰り返しが必要で、再構成時間が長いという課題は残り、さらなる再構成時間の短縮が望まれる。

上記課題を解決するため、本発明は、繰り返しの途中結果であるデータに対して密度補正を行う。密度補正は、補正の対象であるデータの密度を前提とするが、繰り返し演算の途中結果のデータは取得されなかったデータが作成されることで密度が変化する。本発明では、変化する密度に応じて途中結果の密度補正を行う。

すなわち本発明のＭＲＩ装置は、核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備える。密度補正部は、逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行う。

本発明の画像処理装置は、上記画像作成部の機能を備えるものである。

また本発明の画像処理方法は、計測データの逐次近似演算により画像を作成する画像処理方法であって、計測データを逐次近似演算の処理空間データに変換する処理と、逐次近似演算を行った前記処理空間データを計測データに変換し、計測データを更新する処理と、更新後の計測データを用いて、前記変換する処理及び前記更新する処理を繰り返し、その際、更新された計測データに対し、密度補正を行う処理と、を含み、前記密度補正は、繰返し毎に密度補正量を異ならせることを特徴とする。

本発明によれば、計測データの密度の初期値を用いた密度補正だけでなく、毎回の繰り返し処理においても対象となるデータの密度に応じた密度補正を行うことで、逐次近似演算をより少ない繰り返し回数で収束可能にすることができる。

また本発明は、最後に密度推定して密度補正を行うことで、より正解に近づいた再構成画像にすることができる。その分、繰り返し回数を少なくし、高速化しても良い。

本発明が適用されるＭＲＩ装置の一例の全体構成を示す図。第一実施形態のＭＲＩ装置の信号処理部の機能ブロック図。実施形態のＭＲＩ装置の処理の概要を示すフロー図。第一実施形態の繰り返し再構成の手順の一例を示すフロー。第一実施形態の密度推定方法の一例を示す図。第一実施形態の効果を説明する図。第一実施形態の変形例５の繰り返し再構成の手順の一例を示すフロー。変形例５のＵＩ画面の一例を示す図。ＡｌｔｅｒｎａｔｉｎｇＳｐｌｉｔＢｒｅｇｍａｎ法による繰り返し再構成の手順を示すフロー。第二実施形態の繰り返し再構成の手順の一例を示すフロー。第二実施形態の繰り返し再構成の手順の他の例を示すフロー。第二実施形態の変形例の繰り返し再構成の手順の一例を示すフロー。第二実施形態の変形例の繰り返し再構成の手順の他の例を示すフロー。画像処理装置の構成図。第一実施形態のＵＩ画面の一例を示す図。

最初に、本発明が適用されるＭＲＩ装置の全体概要を、図面を参照して説明する。図１は、本発明に係るＭＲＩ装置の一実施形態の全体構成を示すブロック図である。このＭＲＩ装置は、ＮＭＲ現象を利用して被検体の断層画像を得るもので、静磁場発生部２と、傾斜磁場発生部３と、送信部５と、受信部６と、信号処理部７と、シーケンサ４と、中央処理装置（ＣＰＵ）８とを備えている。なお以下の説明では、静磁場発生部２、傾斜磁場発生部３、シーケンサ４、送信部５、及び受信部６を総括して撮像部ともいう。

静磁場発生部２は、被検体１の周りの空間に均一な静磁場を発生させるもので、永久磁石方式、常電導方式あるいは超電導方式の静磁場発生装置を備える。発生する磁場の方向により、垂直磁場方式、水平磁場方式などがあり、本発明はいずれの方式にも適用できる。

傾斜磁場発生部３は、ＭＲＩ装置の座標系（静止座標系）であるＸ、Ｙ、Ｚの３軸方向に傾斜磁場を印加する傾斜磁場コイル９と、それぞれの傾斜磁場コイルを駆動する傾斜磁場電源１０とから成り、シーケンサ４からの命令に従ってそれぞれのコイルの傾斜磁場電源１０を駆動することにより、Ｘ、Ｙ、Ｚの３軸方向に傾斜磁場Ｇｘ、Ｇｙ、Ｇｚを印加する。これら３軸方向の傾斜磁場の組み合わせにより、任意の方向の傾斜磁場を発生することができる。例えば、撮影時には、スライス面（撮影断面）に直交する方向にスライス方向傾斜磁場パルス（Ｇｓ）を印加して被検体１に対するスライス面を設定し、そのスライス面に直交して且つ互いに直交する残りの２つの方向に位相エンコード方向傾斜磁場パルス（Ｇｐ）と周波数エンコード方向傾斜磁場パルス（Ｇｆ）を印加して、エコー信号にそれぞれの方向の位置情報をエンコードする。

シーケンサ４は、高周波磁場パルス（以下、「ＲＦパルス」という）と傾斜磁場パルスをある所定のパルスシーケンスで繰り返し印加する制御手段で、ＣＰＵ８の制御で動作し、被検体１の断層画像のデータ収集に必要な種々の命令を送信部５、傾斜磁場発生部３、および受信部６に送る。

送信部５は、被検体１の生体組織を構成する原子の原子核スピンに核磁気共鳴を起こさせるために、被検体１にＲＦパルスを照射するもので、高周波発振器１１と変調器１２と高周波増幅器１３と送信側の高周波コイル（送信コイル）１４ａとを備える。高周波発振器１１から出力されたＲＦパルスをシーケンサ４からの指令によるタイミングで変調器１２により振幅変調し、この振幅変調されたＲＦパルスを高周波増幅器１３で増幅した後に被検体１に近接して配置された送信コイル１４ａに供給することにより、ＲＦパルスが被検体１に照射される。

受信部６は、被検体１の生体組織を構成する原子核スピンの核磁気共鳴により放出されるエコー信号（ＮＭＲ信号）を検出するもので、受信側の高周波コイル（受信コイル）１４ｂと信号増幅器１５と直交位相検波器１６とＡ／Ｄ変換器１７とを備える。送信コイル１４ａから照射された電磁波によって誘起された被検体１の応答のＮＭＲ信号が被検体１に近接して配置された受信コイル１４ｂで検出され、信号増幅器１５で増幅された後、シーケンサ４からの指令によるタイミングで直交位相検波器１６により直交する二系統の信号に分割され、それぞれがＡ／Ｄ変換器１７でディジタル量に変換されて、信号処理部７に送られる。

信号処理部７は、各種データ処理と処理結果の表示及び保存等を行うもので、その機能の一部はＣＰＵ８により実現される。信号処理部７（ＣＰＵ８）が行う処理には、被検体１の画像再構成や被検体の特性を示す数値の演算、受信コイル感度等の装置の特性を用いた信号や処理結果の補正などが含まれる。

信号処理部７は、光ディスク、磁気ディスク等の外部記憶装置１８と、ディスプレイ１９と、トラックボール又はマウス、キーボード等から成る操作部２０とを備えている。受信部６からのデータがＣＰＵ８に入力されると、ＣＰＵ８が信号処理、画像再構成等の処理を実行し、その結果である被検体１の断層画像をディスプレイ１９に表示すると共に、外部記憶装置１８の磁気ディスク等に記録する。ＣＰＵ８は、上述した処理を行う計算部として機能するとともにシーケンサ４や装置全体の制御を行う制御部として機能する。

操作部２０は、ＭＲＩ装置の各種制御情報や信号処理部７で行う処理の制御情報を入力するもので、ディスプレイ１９に近接して配置され、操作者がディスプレイ１９を見ながら操作部２０を通してインタラクティブにＭＲＩ装置の各種処理を制御する。

なお、図１において、送信側の高周波コイル１４ａと傾斜磁場コイル９は、被検体１が挿入される静磁場発生部２の静磁場空間内に、垂直磁場方式であれば被検体１に対向して、水平磁場方式であれば被検体１を取り囲むようにして設置されている。また、受信側の高周波コイル１４ｂは、被検体１に対向して、或いは取り囲むように設置されている。

送信コイル１４ａと受信コイル１４ｂとは、別個の高周波コイルでもよいが、一つの高周波コイルが送信コイル１４ａと受信コイル１４ｂを兼ねる場合もある。また送受信を兼ねる全身用コイルと、それとは別の局所コイル等の高周波コイルを組み合わせて用いる場合も有る。

本実施形態のＭＲＩ装置は、上述した信号処理部７の機能、特に画像生成の機能として、ＣＳ等の逐次近似演算を行う機能が含まれている。画像生成に関する信号処理部７の機能の概略を図２に示す。

図２に示すように、信号処理部７は、大きく分けて、撮像条件を設定してパルスシーケンスを作成する撮像設定部７１と、画像作成部７３と、を有する。画像作成部７３は、逐次近似演算部（単に演算部という場合もある）７３１と、密度推定部７３２と、密度補正部７３３と、を含む。

撮像設定部７１は、操作部２０等を介してユーザーが設定した撮像方法や撮像条件に従って、パルスシーケンスを計算し、シーケンサ４に設定する。パルスシーケンスは、デカルト座標系の軸に沿ってスキャンする撮像方法のほかラジアルスキャンなど非デカルト座標系スキャンも含まれ、このスキャンの仕方と撮像条件として設定される倍速率によって、計測空間の間引き方、計測されるデータ点の密度などが決まる。

逐次近似演算部７３１は、計測データに対し逆フーリエ変換、ウェーブレット変換などの正則化を行う空間への変換、正則化を行う空間における最小化、フーリエ変換及び正則化を行う空間への変換の逆変換等を計測データの更新を行いながら繰り返し、画像を作成する。

密度推定部７３２は繰り返し演算の途中で更新されたデータに対し、密度の推定を行う。密度の推定は、データの差分から推定する方法や予め定めた密度変化曲線に基づく方法などがある。推定方法の詳細は、後述する実施形態で説明する。密度補正部７３３は、密度推定部７３２が推定した密度を用いて、繰り返し処理途中のデータを密度補正する。

上述した信号処理部７の機能は、ＣＰＵ８に備えられたメモリ（記憶部）やＣＰＵ８に組み込まれたプログラム等のソフトウェアを実行することで実現することができる。またその一部の機能をＡＳＩＣやＦＰＧＡなどのハードウェアで実現してもよい。

次に、上述したＭＲＩ装置による撮像及び信号処理の概要を、図３のフローを参照して説明する。

まずＭＲＩ装置の静磁場空間に被検体１を配置し（Ｓ１０１）、撮像を行う（Ｓ１０２）。撮像は、撮像設定部７１が計算したパルスシーケンスに従って、撮像部が被検体にＲＦパルス及び傾斜磁場パルスを印加するとともに被検体から発生せられるＮＭＲ信号を収集することにより行われる。ＮＭＲ信号は、ｋ空間のマトリクスに配置された計測データ（ｋ空間データともいう）となる。本実施形態では、ｋ空間をすべて埋めるデータを収集するのではなく、アンダーサンプリングした計測データを得る。得られる計測データの密度は、傾斜磁場パルスによるエンコーディングステップなどの撮像パラメータや間引き率（あるいは倍速率）などの撮像条件によって決まる。

次いで画像作成部７３は、撮像により得られた計測データ（ＮＭＲ信号）を用いて、逐次近似演算部７３１の繰り返し演算によって未計測のデータを更新しながら、最終的に目的とする画像を作成する（Ｓ１０３）。目的とする画像には、プロトン密度画像、拡散強調画像等の強調画像、それら画像から派生する計算値画像などが含まれる。この画像作成ステップＳ１０３において、密度補正部７３３は、繰り返し演算の途中で変化していくデータの密度に合わせて、処理途中のデータの密度補正（以下、密度補正処理という）を行う。変化するデータの密度は、密度推定部７３２が推定した密度を用いる。密度推定の手法は、データの差分から推定する方法や予め定めた密度変化曲線に基づく方法などがある。

以下、本発明による画像作成法の具体的な実施形態を説明する。なお図１及び図２に示す装置の構成及び図３に示す動作の概要は、各実施形態に共通しているので、適宜これら図面を参照する。

＜第一実施形態＞
本実施形態は、間引かれた、即ち未取得の計測点を含む計測データに対し、繰り返し演算で未取得データを更新して再構成を行う場合に、毎回の繰り返しにおいて、密度推定部７３２は、未取得データの更新の進み具合から密度を推定する。密度補正部７３３は、推定された密度に基づき密度補正を行う。

本実施形態の画像作成処理（Ｓ１０３）における繰り返し処理の一例を図４にしめす。撮像（Ｓ１０２）で取得されたデータｓ０（ｋ）を密度補正（Ｓ１０３１）し、逆フーリエ変換（Ｓ１０３２）で画像にする。画像をウェーブレット変換（Ｓ１０３３）して、閾値処理（Ｓ１０３４）で大きな係数だけ取り出したのち、逆ウェーブレット変換（Ｓ１０３５）で画像に戻し、フーリエ変換（Ｓ１０３６）でｋ空間（取得されたデータの空間）に戻す。この更新されたデータから密度推定（Ｓ１０３７）を行う。更新されたデータのうち、撮像で取得されたデータはｓ０（ｋ）で置き換え（Ｓ１０３８）最終的な更新されたデータｓ（ｋ）とする。更新データが収束していなければ、ｓ（ｋ）に対し推定された密度で密度補正（Ｓ１０３１）をする処理から繰り返す。更新データが収束していれば繰り返しを終了し、逆フーリエ変換（Ｓ１０４０）で画像を作成して終了する。

ここで繰り返しの１回目で行う密度補正（Ｓ１０３１）は、データｓ０（ｋ）の密度（初期値）に基づいて行う。２回目以降は、更新されたデータｓ（ｋ）の密度、即ちＳ１０３７で推定された密度に基づいて行う。密度補正は、データの密度分布に応じて、データを重み付けする処理であり、密度が高いほど重みは小さく、密度が低いほど重みは大きい。密度の逆数を密度補正量という。

ステップＳ１０３７以外の各ステップの処理は、公知の圧縮センシングなどの逐次近似演算の処理と同様であり、ここではその詳細な説明は省略する。また図４では、正則化を行う空間への変換の一例としてウェーブレット変換を示したが、ウェーブレット変換以外の基底関数、カーブレット変換などを用いてもよいし、複数の基底関数を組み合わせてもよい。

以下、繰り返し処理における密度推定（Ｓ１０３７）について、図５を参照して具体的に説明する。図５（ａ）は、ｋ空間500における計測データの存在を示す図で、取得したデータ（取得データ）の位置を黒色で示している。初期的には、撮像時のサンプリングパターンで決まる。このサンプリングパターンの密度ρ０（ｋ）は、位置ｋの周辺のデータの密度から直接求めても良いし、各位置ｋにおいてデータを取得する確率が、図５（ｂ）の等高線502で示すようなｐ０（ｋ）として与えられているのであれば、その確率ｐを密度ρとしても良い。

更新途中のデータについては、更新の進み具合に応じて密度を推定する。このため、まず、データの更新の不十分さα（ｋ）を以下のように定義する。

［数１］
α（ｋ）＝（ｂ（ｋ）―ｃ（ｋ））÷ｂ（ｋ）（１）
式（１）中、ｂ（ｋ）は撮像で取得されたデータ（図５（ａ）の黒色の領域、データ取得領域と呼ぶ）の信号値ｓ０（ｋ）の絶対値の合計を取得されたデータの点数で割った平均であり、ｃ（ｋ）は更新されたデータの信号値ｓ（ｋ）において、データが取得されていない領域（白色の領域、データ未取得領域と呼ぶ）での絶対値の合計をデータ未取得領域の点数で割った平均である。α（ｋ）はデータが更新されていない状態では１になり、データ未取得領域の平均値が、取得データの平均値の大きさまで更新されれば０になる。

次に、式（１）により求めた更新の不十分さα（ｋ）と、データ取得時の密度分布の初期値ρ０（ｋ）とを用いて、次式（２）により、データの更新により変化した密度ρ（ｋ）を推定する。

［数２］
ρ（ｋ）＝α（ｋ）×（ρ０（ｋ）－１）＋１（２）
式（２）からわかるように、α（ｋ）＝１（データが更新されていない）のときの密度はρ０（ｋ）と推定され、α（ｋ）＝０（データが完全に更新されたと推定される）のときの密度は１と推定されることになる。

なお、この例ではｂ（ｋ）、ｃ（ｋ）はデータ全体で平均をとり、α（ｋ）はｋによらず一定の値となる場合を説明したが、ｋ空間を区切った小領域毎にｂ（ｋ）及びｃ（ｋ）を定義してもよい。すなわち、図５（ａ）に小さい四角で示すように、位置ｋの周辺の領域501に含まれる、取得データの信号値ｓ０（ｋ）の絶対値の合計を、領域501で取得されたデータの点数で割った平均とし、ｃ（ｋ）の定義を、領域501内で更新されたデータの信号値ｓ（ｋ）の絶対値の合計を、領域501内で更新されたデータの点数で割った平均とする。この場合は、α（ｋ）はｋによって変化するような値となる。

この場合、ｂ（ｋ）、ｃ（ｋ）を求める際に平均をとる領域501としては、ｋの周辺の領域でなく、図５（ｃ）に示すような、ｋを通る正方形の辺付近（黒色で表す領域）503としても良いし、図５（ｂ）に示すようなデータ取得確率ｐ０（ｋ）が与えられている場合には、図５（ｄ）に示すように、確率ｐの値が近い領域504としても良く、自由に決めることができる。

密度補正部７３３は、上述したステップＳ１０３７で推定した密度を用いて、繰り返し演算途中のデータの密度補正を行う（Ｓ１０３１）。すなわち密度に応じた密度補正量で計測データの密度補正を行い、繰り返される各処理に進む。
本実施形態によれば、データ更新の進み具合をもとに密度を推定することで、繰り返し途中でも、そのときの確からしい密度を元に、高い精度で密度補正を行うことができる。

本実施形態による密度補正を行った場合の画像と、密度補正を行わない場合の画像とを比較した結果を図６に示す。図６中、（Ａ）はデータの不足がない場合（フルサンプリングした場合）の画像、（Ｂ）はデータが不足している状態で密度補正を行わずに逐次近似演算を行い、繰り返し１０回で再構成した画像、（Ｃ）は（Ｂ）と同様にデータが不足している状態で繰り返し途中で密度補正を行いながら逐次近似演算で再構成した画像で繰り返し回数は５回である。

（Ｄ）及び（Ｅ）は、（Ａ）と（Ｂ）との差分、（Ａ）と（Ｃ）との差分である。（Ｂ）と（Ｃ）との比較或いは（Ｄ）と（Ｅ）との比較から、本実施形態では、密度補正を行わない場合の半分の繰り返し回数で、フルサンプリングして得られる画像に近づいていることがわかる。すなわち、本実施形態により画像再構成の精度が良好になること、また従前の画像再構成より短時間で元の画像に近づけることがわかる。
本実施形態の適用は予め設定しておいてもよいが、ユーザーインターフェイス（図１：操作部２０）を介して、ユーザーが繰り返し中の密度補正を行うか否かを選択するようにしてもよい。図１５にユーザーインタフェイス画面の一例を示す。この例では、繰り返し中の密度補正についてユーザーがプルダウン選択方式で、実施するか否かを選択可能にしている。

＜第一実施形態の変形例１＞
第一実施形態では、密度推定方法として、データ未取得領域の更新データから直接データ更新の進み具合を求めたが、データ未取得領域の更新データを使う代わりに、データ取得領域における更新データから間接的にデータ更新の進み具合を求めても良い。

本変形例では、データの更新の不十分さα（ｋ）を次式（３）と定義し、第一実施形態の式（２）を用いて、密度ρ（ｋ）を推定する。

［数３］
α（ｋ）＝Δｂ（ｋ）÷Δｂ０（ｋ）（３）
ここで、Δｂ（ｋ）及びΔｂ０（ｋ）は、それぞれ、次式（４）、（５）で表され、Δｂ（ｋ）はデータ取得領域内の位置ｋにおける更新データの信号値ｓ（ｋ）と同位置ｋにおける取得データの信号値ｓ０（ｋ）との差分の絶対値を、データ取得領域で総和したものを表し、Δｂ０（ｋ）は、データ取得領域内の位置ｋにおける取得データの信号値ｓ０（ｋ）をその位置の密度ρ０（ｋ）で除して密度補正したものと信号値ｓ０（ｋ）との差分の絶対値をデータ取得領域で総和したものを表す。

［数４］
Δｂ（ｋ）＝Σ｜ｓ（ｋ）－ｓ０（ｋ）｜（４）
［数５］
Δｂ０（ｋ）＝Σ｜ｓ０（ｋ）÷ρ０（ｋ）－ｓ０（ｋ）｜（５）
データ取得領域は、ｋ空間全体におけるデータ取得領域でもよいし、位置ｋを含む一部のデータ取得領域としてもよい。

式（３）においても、更新が不十分なときはｓ（ｋ）がｓ０（ｋ）からずれるため、α（ｋ）が大きくなり１に近づき、データが完全に更新されたときはｓ（ｋ）がｓ０（ｋ）からずれないため、α（ｋ）が小さくなり０に近づき、式（２）を用いて、密度を推定することができる。

＜第一実施形態の変形例２＞
第一実施形態では、式（２）で推定した密度ρ（ｋ）は、データ取得領域とデータ未取得領域とで区別されていないが、データ取得領域とデータ未取得領域とでは密度を異ならせてもよい。例えば、データ取得領域で密度を変化させず１とし、データ未取得領域のみ、式（１）で定義される更新の不十分さα（ｋ）を用いて、

［数６］
ρ（ｋ）＝１－α（ｋ）（６）
としてもよい。

または、次式（７）で示すように密度を定義し、領域に応じて密度を変化させてもよい。
［数７］
ρ（ｋ）＝ｗ（ｋ）×α（ｋ）×（（１－ｗ（ｋ））×ρ０（ｋ）－１）＋１（７）
式（７）において、ｗ（ｋ）は、領域に応じた重みで、例えば、データ取得領域ではｗ（ｋ）＝０、データ未取得領域ではｗ（ｋ）＝１とすると、
ρ（ｋ）＝１（データ取得領域）
ρ（ｋ）＝１－α（ｋ）（データ未取得領域）
となり、式（６）のケースとなる。

＜第一実施形態の変形例３＞
第一実施形態或いはその変形例１，２で推定した密度をさらに調整しても良い。調整方法としては、次のような方法が挙げられる。α（ｋ）は０以上１以下になるように調整し、ρ（ｋ）が必ずρ０（ｋ）から１の範囲となるように調整する。或いは、ρ（ｋ）が単調に増加するように調整する、すなわち、繰り返しの際にｎ番目に推定された密度をρｎとすると、ひとつ前の繰り返しにおいて推定された密度ρｎ－１に対し、ρｎ≧ρｎ－１となるようにρ（ｋ）を調整しても良い。

このような調整を行うことにより、取得されたデータや更新データの信号値にアーチファクトやノイズ等の原因により異常な値が含まれている場合にも、異常な値につられることのない密度推定とそれに基づく密度補正を行うことができる。

＜第一実施形態の変形例４＞
第一実施形態では、更新中のデータを用いて更新の進み具合αを算出し、これをもとに更新途中のデータの密度の推定を行ったが、更新の進み具合αを演算の繰り返し回数ｎを変数とする関数で定義し、それを用いて、更新途中のデータの密度を推定してもよい。関数として例えば、次式（８）を用いることができる。

［数８］
α（ｎ）＝０．５＾ｎ（８）
この場合にも、上述した式（２）を用いて、密度ρ（ｋ）を推定することができる。なお関数は、初期値を１として漸次減少するものであれば式（８）に限定されない。

本変形例においても、上述した変形例２を適用することができ、例えば、データ取得領域の密度は１に固定し、データ未取得領域のみ密度を変化させてもよい。

本変形例によれば、更新されたデータを使わずに、予め定めた関数をもとに密度を推定できるので、密度推定を簡便に行うことができる。

＜第一実施形態の変形例５＞
第一実施形態では、繰り返し終了後にｓ（ｋ）を逆フーリエ変換（Ｓ１０４０）して最終的な画像を得る場合を示したが、ｓ（ｋ）に対し密度推定を行い、推定された密度で密度補正をしてから逆フーリエ変換して最終的な画像を得てもよい。

本変形例の処理の流れを図７に示す。Ｓ１０３８までは、図４と同様であるが、本変形例では、繰り返し演算が終了後に得られた計測データに対し、密度推定を行い、推定した密度に基づく密度補正を行う。密度推定の手法は、第一実施形態及びその変形例１～４の手法のいずれを用いてもよい。

本変形例の適用は予め設定しておいてもよいが、ユーザーインターフェイス（図１：操作部２０）を介して、ユーザーが繰り返し終了後に密度補正を追加するか否かを選択するようにしてもよい。図８にユーザーインタフェイス画面の一例を示す。この例では、（Ａ）、（Ｂ）に示すように、繰り返し中の密度補正及び繰返し終了後の密度補正それぞれについてユーザーがプルダウン選択方式で、実施するか否かを選択可能にしている。

本変形例によれば、繰り返し処理で十分収束していなくても、密度補正を行うことで画質を改善できるため、さらに少ない繰り返しで良い近似画像を得ることができる。またユーザーによる選択を可能にすることによって、画質優先か画像作成までの時間優先かなどユーザーの希望に応じた画像再構成を行うことができる。

＜第二実施形態＞
第一実施形態では、一般的な逐次近似の繰り返し演算に対して、繰り返しごとにデータｓ（ｋ）を密度補正する処理を追加する再構成を行ったが、本実施形態では、逐次近似演算部７３１が、繰り返し処理に高速なAlternating Split Bregman法を用いる。本実施形態でも繰り返し処理において密度補正を追加することは第一実施形態と同様であるが、図４のようにデータｓ（ｋ）を密度補正する処理（Ｓ１０３１）をフローに追加するのではなく、処理全体を変更する点が異なる。

まず、本実施形態で採用するAlternating Split Bregman法による逐次近似演算について説明する。

この逐次近似演算による画像再構成は、逆フーリエ変換によって画像空間データとした計測データに対し、ウェーブレット変換とカーブレット変換とを用いて、これらのＬ１ノルムの最小化を行い、再構成画像を得る処理であり、次式（９）で表される問題を解く処理である。

式（９）中、ｆは撮像データ、uが再構成画像、Φはフーリエ変換をしてデータ取得領域のみ抽出する変換、Ψcはカーブレット変換、Ψｗはウェーブレット変換、θcはカーブレット変換係数、θwはウェーブレット変換係数、Λcはカーブレット係数に設定する重み、Λwはウェーブレット係数に設定する重みである。

この問題は、Bregman法により、次式（１０）を繰り返し解くことに置き換えられる。以下の式において、上付き文字の「ｎ」は繰り返し回数を表す（以下同じ）。

式（１０）は、Split Bregman法により、さらに以下のように置き換えられる。

式（１１）は、Alternating Split Bregman法により、さらに以下のように置き換えられる。

式（１２－１）～（１２－３）で表されるｕ、θc、θwそれぞれの最適化問題は解析的に解くことができる。

このAlternating Split Bregman法による画像作成処理（図３：Ｓ１０３）は図９のようになる。まず、ｕ、θw、θc、ｂc、ｂwの初期値を設定する（Ｓ２０３１）。初期値は例えば0で良い。次に、式（１２－１）でｕを更新する（Ｓ２０３２）。次に、式（１２－２）でθcを、式（１２－３）でθwを、更新する（Ｓ２０３３）。次に、式（１２－４）でｂcを、式（１２－５）でｂwを、更新する（Ｓ２０３４）。次に、式（１２－６）でｆを更新する（Ｓ２０３５）。解ｕが収束していなければ、S２０３２に戻り処理を繰り返す。解ｕが収束してれば終了する。

本実施形態では、この処理に密度補正を導入する。このため、まず、ｕ、θc、θwを求める式（１０）を、次式（１３）ように変更する。

式（１３）において、Ｆはフーリエ変換を表し、Φはフーリエ変換を含まずデータ取得領域のみ抽出する変換と変更し、（Ω^ｎ）^－１は密度補正処理を表す。密度補正の手法は、第一実施形態で説明した手法と同様であり、n回目の繰り返しで得られたデータ（Ｆｕ^ｎ）から推定した密度或いは繰り返し回数nで決まる推定密度に基づく。Ω^ｎは任意に決めることができるので、ｕに依存しないものとする。μは項の重みであり、0より大きい値を取る。

ここで（Ω^ｎ）^－１の密度補正量を決めるための密度推定において、あるｎ以降、密度が１で密度補正をしない値になるとすれば、式（９）と同じ解に収束することになり、式（１３）で式（９）を解くことができる。どんなにnを大きくしても密度が1にならない場合でも、あるｎ以降Ω^ｎがΩ∞になるとすれば、式（９）の変形である式（１４）の解に収束することになる。つまり、第一実施形態の変形例３のようにΩ^ｎが収束するようにしさえすれば、式（１３）の解は収束することが保証される。

本実施形態では、さらに、式（１３）を以下の式（１５）に置き換えることで、式（１３）と同じ解の組み合わせ（シーケンス）｛（ｕ¹，θ¹c，θ¹w）（ｕ²，θ²c，θ²w）（ｕ³，θ³c，θ³w）・・・｝を得ることが可能となる。なお、ｕ, θc, θw, ｆ, ｂc, ｂwは実数として、ベクトルとしても良いし、複素数として、複素数の行列表現を使ってｎ行２列の行列としても良い。

なお、この繰り返しによる不動点は確かに式（１３）の解となっている。

さらに、式（１５－１）を式（１６－１）～（１６－３）のように変更し、ｕ, θc, θwで順番に最適化する。式（１６－４）～（１６－６）は式（１５－２）～（１５－４）と同じである。

式（１６－１）は具体的には、

式（１６－２）、（１６－３）は具体的には、

上述した密度補正処理（式（１６）に基づく処理）を導入した本実施形態の画像作成処理（Ｓ１０３）の流れを図１０に示す。

まず、ｕ、θw、θc、ｂc、ｂwの初期値を設定する（Ｓ３０３１）。このステップは、図９のステップＳ２０３１と同じであり、初期値は例えば０で良い。次に、密度補正を含む式（１６－１）でｕを更新する（Ｓ３０３２）。次に、密度補正を含む式（１６－２）でθcを、式（１６－３）でθwを、それぞれ更新する（Ｓ３０３３）。次に、密度補正を含む式（１６－４）でｂcを、式（１６－５）でｂwを、それぞれ更新する（Ｓ３０３４）。次に、密度補正を含む式（１６－６）でｆを更新する（Ｓ３０３５）。解ｕが収束していなければ、Ｓ３０３２に戻り処理を繰り返す。解ｕが収束してれば終了する。

各繰り返しで更新する変数（θc,θw,ｂc,ｂw）は、以下のように変えても良い。第三式及び第四式の左辺は、本文中、ｂ（～）^ｎ _ｃ、ｂ（～）^ｎ _ｗと記す。

その場合、式（１６－１）～（１６－６）は以下のようになる。

また式（１７）は以下のようになる。

また式（１８）は以下のようになる。

この場合の、画像作成処理（Ｓ１０３）の流れを図１１に示す。まず、ｕ、ｕw、ｕc、ｂ（～）^ｎ _ｃ、ｂ（～）^ｎ _wの初期値を設定する（Ｓ４０３１）。初期値は例えば０で良い。次に、密度補正を含む式（２０－１）でｕを更新する（Ｓ４０３２）。次に、密度補正を含む式（２０－２）でｕcを、密度補正を含む式（２０－３）でｕwを、更新する（Ｓ４０３３）。次に、密度補正を含む式（２０－４）でｂ（～）^ｎ _ｃを、密度補正を含む式（２０－５）でｂ（～）^ｎ _wを、更新する（Ｓ４０３４）。次に、密度補正を含む式（２０－６）でｆを更新する（Ｓ４０３５）。解ｕが収束していなければ、Ｓ４０３２に戻り処理を繰り返す。解ｕが収束してれば終了する。

＜第二実施形態の変形例１＞
第二実施形態では、式（１４）で示すように、ウェーブレット変換やカーブレット変換（すなわちL１ノルムを最小化する空間への変換）を行う前に密度補正（Ω）^－１を行ったが、計測データに一致させるという条件ｆ＝Φｕにおいて計測データに密度補正を行ってもよい。この場合、第二実施形態の式（１３）における密度補正Ωをカーブレット変換及びウェーブレット変換の項に追加する代わりに、次式のように、||ΦＦｕ－ｆ||の項に追加する。

この場合も処理の流れは図１０、図１１に示す流れと同じである。

＜第二実施形態の変形例２＞
第二実施形態においても、第一実施形態の変形例５と同様に、繰り返し終了後の計測データから最終的な画像を得る前に、繰り返し終了時に得られた計測データに対し密度補正をしてから逆フーリエ変換して最終的な画像を得てもよい。

この場合は、図１２、図１３に示すように、図１０、図１１の繰り返し終了後にｕをフーリエ変換する処理（Ｓ５０３８）、密度補正（Ｓ１０３９）、及び逆フーリエ変換（Ｓ１０４０）が追加される。密度推定の手法は、第一実施形態及びその変形例１～４の手法のいずれを用いてもよい。

本変形例によれば、繰り返し処理で十分収束していなくても、密度補正を行うことで画質を改善できるため、さらに少ない繰り返しで良い近似画像を得ることができる。

以上、本発明の画像処理をＭＲＩ装置の信号処理部７（画像作成部）で実施する実施形態を説明したが、図１４に示すように、ＭＲＩ装置１００とは独立した画像処理装置２００で画像作成部７３の機能を実施することも可能である。画像処理装置２００は、ＭＲＩ装置１００と有線、無線或いはインターネット回線等で接続されているか、可搬媒体によって、ＭＲＩ装置１００の撮像部が取得した計測したデータを受け取り、画像再構成演算を行う。作成した画像は、画像処理装置２００に付属するディスプレイに表示したり記憶装置に格納したりしてもよいし、ＭＲＩ装置１００に送ってもよい。

２：静磁場発生部、３：傾斜磁場発生部、４：シーケンサ、５：送信部、６：受信部、７：信号処理部、８：中央処理装置（ＣＰＵ）、１４ａ：送信コイル、１４ｂ：受信コイル、１８：外部記憶装置、１９：ディスプレイ、２０：操作部、７１：撮像設定部、７３：画像作成部、７３１：逐次近似演算部、７３２：密度推定部、７３３：密度補正部、１００：ＭＲＩ装置、２００：画像処理装置。

Claims

核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記密度補正部は、前記計測データを複数の領域に分割し、領域毎に求めた密度から密度補正量を決定して、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行うことを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行い、その際、補正対象となる計測データのうち、前記撮像部が取得した計測点のデータと、逐次近似演算によって推定された計測点のデータと、を異なる重みで密度補正することを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行い、その際、予め設定された密度の変化に基づき密度補正を行うことを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行い、その際、逐次近似演算によって更新された計測データの更新の不十分さが、０以上１以下であるという条件のもとで、前記密度補正量を調整して密度補正を行うことを特徴とする磁気共鳴イメージング装置。
請求項１ないし３のいずれか1項に記載の磁気共鳴イメージング装置であって、
前記密度補正部は、密度の変化が単調増加するという条件のもとで、前記密度補正量を変化させて密度補正を行うことを特徴とする磁気共鳴イメージング装置。
請求項１ないし４のいずれか１項に記載の磁気共鳴イメージング装置であって、
前記密度補正部は、密度補正量の初期値を前記撮像部におけるデータ収集方法に応じて決定することを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行うとともに、
前記密度補正部は、前記逐次近似演算部による逐次近似演算の繰り返し終了後に、密度補正を行うことを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記逐次近似演算部は、前記計測データから実空間データの変換及びその逆変換を行うフーリエ変換部と、前記実空間データから正則化を行う空間データへの変換及びその逆変換を行う空間変換部とを備え、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行うことを特徴とする磁気共鳴イメージング装置。
請求項１ないし４のいずれか１項に記載の磁気共鳴イメージング装置であって、
前記逐次近似演算部は、前記計測データから実空間データの変換及びその逆変換を行うフーリエ変換部と、前記実空間データから正則化を行う空間データへの変換及びその逆変換を行う空間変換部とを備えることを特徴とする磁気共鳴イメージング装置。
核磁気共鳴信号からなる計測データを取得する撮像部と、前記撮像部がアンダーサンプリングして取得した計測データを用いて逐次近似演算により画像を作成する画像作成部と、を備え、
前記画像作成部は、前記計測データに対し逐次近似演算によって前記計測データを更新する逐次近似演算部と、前記計測データの密度補正を行う密度補正部と、を備え、
前記逐次近似演算部は、前記逐次近似演算として、Alternating Split Bregman法に基づく演算を行い、
前記密度補正部は、前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら密度補正を行うことを特徴とする磁気共鳴イメージング装置。
請求項１ないし４のいずれか１項に記載の磁気共鳴イメージング装置であって、
前記逐次近似演算部は、前記逐次近似演算として、Alternating Split Bregman法に基づく演算を行うことを特徴とする磁気共鳴イメージング装置。
計測データを入力し、逐次近似演算により画像を作成する演算部を備えた画像処理装置であって、
前記逐次近似演算の繰返しにおいて、繰り返し毎に密度補正量を変化させながら前記計測データの密度補正を行う密度補正を行う密度補正部を備え、前記密度補正部は、請求項１ないし４のいずれか１項に記載の密度補正部であることを特徴とする画像処理装置。
請求項１２に記載の画像処理装置であって、
前記計測データは、磁気共鳴イメージング装置で計測された核磁気共鳴信号からなる計測データであることを特徴とする画像処理装置。
計測データの逐次近似演算により画像を作成する画像処理方法であって、
前記逐次近似演算は、前記計測データから実空間データの変換及びその逆変換、並びに、前記実空間データから正則化を行う空間データへの変換及びその逆変換を含む変換処理、
逐次近似演算を行った処理空間のデータを計測データに変換し、計測データを更新する処理と、
更新後の計測データを用いて、前記変換処理及び前記更新する処理を繰り返し、
その際、更新された計測データに対し、密度補正を行う処理と、を含み
前記密度補正は、繰返し毎に密度補正量を異ならせることを特徴とする画像処理方法。
計測データの逐次近似演算により画像を作成する画像処理方法であって、
前記逐次近似演算は、Alternating Split Bregman法に基づく演算を行い、計測データの更新を繰り返し、
その際、更新された計測データに対し、密度補正を行う処理、を含み
前記密度補正は、繰返し毎に密度補正量を異ならせることを特徴とする画像処理方法。