本發(fā)明涉及一種磁共振成像(mri：magneticresonanceimaging)技術。尤其，涉及一種使用所取得的圖像計算生物體內的磁化率分布的圖像處理技術。

背景技術：

磁共振成像裝置是利用設置于靜磁場內的氫原子核(質子)與特定頻率的高頻磁場產生共振的核磁共振現象的非侵襲性醫(yī)用圖像診斷裝置。核磁共振信號根據質子密度、衰減時間等各種物性值而變化，因此可以通過mri得到的圖像描繪出生物體組織的構造或組成、細胞性狀等各種生物體信息。

近年來，作為能夠通過mri測定的物性值，關注生物體組織間的磁化率差。磁化率是表示靜磁場中的物質被磁化的程度的物性值。在生物體內，存在靜脈血中的脫氧血紅蛋白或鐵蛋白質等順磁性體、以及成為占據生物體組織的大部分的水或鈣化的基礎的鈣等反磁性體。對生物體組織間的磁化率的差即磁化率差定量地進行圖像化，從而能夠應用于腦缺血性疾病的診斷或癌癥放射治療效果預測，神經變性疾病的診斷。

使用mri對生物體組織間的磁化率差進行圖像化的方法被稱為定量磁化率匹配(qsm：quantitativesusceptibilitymapping)。qsm方法是根據測量出的mr圖像的相位信息計算出因生物體組織間的磁化率差產生的空間性的磁場變化，并使用磁場和磁化率的關系式計算出磁化率分布的方法。

然而，磁場分布是對磁化率分布進行空間性卷積積分而得的分布，因此根據磁場分布求出磁化率分布的情況下，成為卷積積分的逆問題，無法得到唯一的解。

一般，通過最小二乘法，根據磁場分布求出磁化率分布。此時，導入誤差函數，將這些最小化的值設為解。作為有代表性的方法，例如有通過磁場分布和點偶極子磁場的k空間上的運算計算出磁化率分布的tkd(truncated-basedk-spacedivision，基于階段的k空間的劃分)(例如，參照非專利文獻1)、將通過tkd方法計算出的磁化率分布和通過閾值處理提取出的微細結構而得的磁化率分布通過反復運算進行合成的迭代swim(susceptibilityweightedimagingandmapping，磁敏感加權成像以及匹配)方法(例如，參照非專利文獻2、專利文獻1)、使用附有正規(guī)化的最小二乘法的medi方法(morphologyenableddipoleinversion，形態(tài)啟用偶極反演)方法(例如，參照非專利文獻3、專利文獻2)

現有技術文獻

專利文獻

專利文獻1：美國專利申請公開第2011/0262017號說明書

專利文獻2：美國專利申請公開第2011/0044524號說明書

非專利文獻

非專利文獻1：shmuelik等，“magneticsusceptibilitymappingofbraintissueinvivousingmriphasedata”、magneticresonanceinmedicine、2009年、62卷、1510-1522頁

非專利文獻2：tangj等，“improvingsusceptibilitymappingusingathreshold-basedk-space/imagedomainiterativereconstructionapproach”，magneticresonanceinmedicine，2013年，69卷，1396-1407頁

非專利文獻3：liut等，“morphologyenableddipoleinversion(medi)fromasingle-angleacquisition:comparisonwithcosmosinhumanbrainimaging”，magneticresonanceinmedicine,2011年，66卷，777-783頁

技術實現要素：

發(fā)明要解決的課題

tkd方法將不能進行逆運算的區(qū)域取舍為0或任意值。因此，tkd方法雖然具有簡單且計算時間短的優(yōu)點，但因對魔角區(qū)域近旁進行取舍，因此存在產生偽影，或計算磁化率的準確度下降的問題。

迭代swim方法對通過tkd方法得到的磁化率分布進行可維持微細結構的掩碼處理，反復進行更新磁化率分布的處理。雖然能夠計算出維持微細結構的磁化率分布，但存在不能完全去除偽影的問題。此外，為了提取微細結構，需要進行各種閾值處理，因此存在若不對每個被檢體設定恰當的閾值，則無法恰當地提取微細結構的課題。

在medi方法中，以磁場和磁化率為制約條件，通過反復運算使被稱為正規(guī)化項的給出平滑效果的項最小化，來計算出磁化率分布。這樣的帶有正規(guī)化的最小二乘法是對與正規(guī)化參數的值對應的大小的信號值進行取舍的處理，因此需要根據混入到磁場分布中的噪聲恰當地設定正規(guī)化項。然而，通常，無法事先得到混入到磁場分布中的噪聲的信息，因此根據使正規(guī)化參數變化多次而計算出的殘差的變化曲線計算出恰當的正規(guī)化參數。因此，在medi方法中，為了設定恰當的正規(guī)化參數，需要龐大的計算時間。

本發(fā)明是鑒于上述問題提出的，其目的是提供一種使用mri計算出磁化率分布時，通過簡單處理減少偽影和噪聲，提高計算出的磁化率分布的準確率和精度的技術。

用于解決課題的手段

在本發(fā)明中，從通過mri取得的復圖像經由相位圖像獲得磁場分布。在基于磁場和磁化率的關系式的制約條件下，反復對根據磁場分布計算出的磁化率分布進行平滑化處理。具體地，進行通過反復運算使由預先決定的初始磁化率分布和磁場分布定義的誤差函數最小化的最小化處理來計算出磁化率分布。此時，在最小化處理中，每算出更新磁化率分布時，對該更新磁化率分布進行平滑化，該更新磁化率分布是根據誤差函數計算出的。

發(fā)明效果

使用mri計算出磁化率分布時，通過簡單處理減少偽影和噪聲，提高計算出的磁化率分布的精度。

附圖說明

圖1的(a)是本發(fā)明實施方式的垂直磁場方式的磁共振成像裝置的外觀圖，(b)是本發(fā)明實施方式的水平磁場方式的磁共振成像裝置的外觀圖，(c)是本發(fā)明實施方式的提高了開放感的磁共振成像裝置的外觀圖。

圖2是表示本發(fā)明實施方式的mri裝置的概要結構的框圖。

圖3是本發(fā)明實施方式的計算機的功能框圖。

圖4是本發(fā)明實施方式的qsm處理的流程圖。

圖5是用于說明rssg(rf-spoiled-steady-stateacquisitionwithrewoundgradient-echo，射頻擾相穩(wěn)態(tài)采集反轉梯度回波)序列的脈沖序列的說明圖。

圖6是本發(fā)明實施方式的磁場分布計算處理的流程圖。

圖7是本發(fā)明實施方式的磁化率分布計算處理的流程圖。

圖8是用于說明本發(fā)明實施方式的磁化率分布計算處理的流程的說明圖。

圖9是本發(fā)明實施方式的變形例的磁化率分布計算處理的流程圖。

圖10是用于說明本發(fā)明實施方式的變形例的磁化率分布計算處理的流程的說明圖。

具體實施方式

對應用本發(fā)明的實施方式的一例進行說明。以下，在用于說明實施方式的所有附圖中，只要沒有特別限定，則對具有相同功能的部件賦予相同的符號，省略其重復的說明。此外，并不通過以下的記述來限定本發(fā)明。

[計算磁化率分布的原理]

在本實施方式中，使用mri裝置在預先決定的回波時間拍攝復圖像i，根據復圖像i計算捕捉到因生物體組織間的磁化率差而產生的磁場變化的磁場分布δ，并根據磁場分布δ計算出磁化率分布χ。通過所謂的qsm方法來計算出磁化率分布。說明實現這些的mri裝置前，說明在qsm方法中根據復圖像i計算磁化率分布χ的方法。

在qsm方法中，根據通過梯度回波(gradientecho，gre)方法拍攝到的相位圖像，計算出因生物體組織間的磁化率差而產生的局部的磁場變化。

設位置向量為r時，用以下的式(1)表示因組織間的磁化率差而產生的相對的磁場變化(磁場分布)δ(r)。

[公式1]

其中，p(r)表示相位圖像，γ表示質子的核磁旋轉比，b0表示靜磁場強度，te表示回波時間。

此外，通過與靜磁場相關的麥克斯韋方程式，使用生物體內的磁化率分布χ(r)，用以下的式(2)表示磁場分布δ(r)。

[公式2]

其中，α表示由向量(r’－r)和靜磁場方向形成的角度，d(r)表示點偶極子磁場。

如式(2)所示，用磁化率分布χ(r)和點偶極子磁場d(r)的卷積積分來表示磁場分布δ(r)。因此，通過對式(2)的兩邊進行傅里葉變換，將式(2)變換為以下的式(3)。

其中，k＝(kx、ky、kz)表示k空間上的位置向量，δ(k)、x(k)、d(k)表示磁場分布δ(r)、磁化率分布χ(r)、點偶極子磁場d(r)的傅里葉分量。

如式(3)所示，能夠將磁場分布的傅里葉分量δ(k)除以點偶極子磁場的傅里葉分量d(k)來計算出磁化率分布的傅里葉分量x(k)。然而，在式(3)中，在d(k)＝0的近旁區(qū)域，其倒數發(fā)散，因此無法直接計算出x(k)。

將該成為d(k)＝0的區(qū)域稱作魔角(magicangle)，成為相對于磁場方向具有大致54.7°的2倍的頂角的反向雙圓錐區(qū)域。因魔角的存在，根據磁場分布計算出磁化率分布的qsm方法回歸到病態(tài)逆問題(ill-conditionedinverseproblem)。

作為解決該問題的有代表性的方法，如上所述，有tkd方法、迭代swim方法(以下，簡稱為swim方法)、medi方法。

在tkd方法中，取點偶極子磁場d(r)的傅里葉分量d(k)的倒數，將其魔角區(qū)域近旁舍位至0或任意值。之后，乘以所測量的磁場分布δ(r)的傅里葉分量δ(k)，在實際空間進行傅里葉逆變換來計算出磁化率分布χ(r)。

在迭代swim方法中，對通過tkd方法計算出的磁化率分布，通過閾值處理提取微細結構。接著，計算出對微細結構以外的磁化率進行了掩碼處理后的微細結構磁化率分布。在k空間上混合原本的磁化率分布和微細結構磁化率分布，并計算修正磁化率分布。對于修正磁化率分布，通過閾值處理再次提取微細結構。反復進行這樣的一系列處理，來計算出磁化率分布。

在medi方法中，如以下的式(4)所示，以磁場和磁化率的關系式為制約條件，通過反復運算使稱為正則化的項的給出平滑效果的項最小化，來計算出磁化率分布χ(r)。

[公式4]

其中，χ表示計算對象的磁化率分布向量，δ表示根據相位圖像計算出的磁場分布向量，w表示在對角分量具有基于絕對值圖像計算出的加權向量的矩陣，m表示在對角分量具有基于絕對值圖像的空間梯度將邊緣區(qū)域設為0，將其他區(qū)域設為1的二值邊緣掩碼向量的矩陣。此外，fd是表示卷積積分的算子，(納布拉，nabla)是表示三維空間梯度算子，∥a∥n是表示任意向量a的ln范數(norm)。此外，式(4)的第2項λ為取得平滑化效果和磁化率精度的平衡的正規(guī)化參數。

如式(4)所示，關于medi方法，作為正規(guī)化項，使用對計算過程中的磁化率分布的全變差(tv：totalvariation)乘以基于絕對值圖像的邊緣掩碼m而得的向量的l1范數。一般，任意向量a的帶有l(wèi)1范數正規(guī)化的最小二乘解成為稀疏(sparse：疏)的解。因此，在medi方法中，絕對值圖像的變化較少的區(qū)域被平滑化，計算出保存有變化較大的邊緣部分的磁化率分布。

可以將包含medi方法的帶有正規(guī)化的最小二乘法解釋為對與正規(guī)化參數的值對應的大小的信號值進行舍位的處理。因此，當正規(guī)化參數過大時噪聲和微細結構都被去除，并且圖像整體的磁化率減少，因此計算磁化率下降。此外，當正規(guī)化參數過小時，雖然維持微細結構或圖像整體的磁化率，但噪聲增加。

因此，在medi方法中，需要設定恰當的正規(guī)化參數。設定了恰當的正規(guī)化參數的情況下，式(4)的第1項的殘差的標準偏差等于混有測量出的磁場分布的噪聲的標準偏差。然而，通常無法提前得到混入到磁場分布中的噪聲的信息，因此根據使正規(guī)化參數變化多次后計算出的殘差的變化曲線計算出恰當的正規(guī)化參數。因此，在medi方法中，為了設定恰當的正規(guī)化參數，需要大量的計算時間。

在本實施方式中，根據上述公知的各方法，提供一種使用mri計算出磁化率分布時，通過簡單的處理降低偽影(artifact)和噪聲，并且能夠計算出高準確度且高精度的磁化率分布的mri裝置以及qsm方法。

[mri裝置的外觀]

首先，對本實施方式的磁共振拍攝裝置(mri裝置)進行說明。圖1的(a)～圖1的(c)是本實施方式的mri裝置的外觀圖。圖1的(a)是使用通過電磁線圈生成靜磁場的隧道型磁鐵的水平磁場方式的mri裝置100。圖1的(b)是為了提高開放感將磁鐵上下分離的漢堡包型(開放型)垂直磁場方式的mri裝置120。此外，在圖1的(c)是使用與圖1的(a)相同的隧道型磁鐵，使磁鐵的縱深變短，且傾斜而提高了開放感的mri裝置130。

在本實施方式中，也可以使用具有這些外觀的mri裝置的任一個。

另外，這些為一例，本實施方式的mri裝置并不限定于這些方式。在本實施方式中，無論裝置的方式、類型，可以使用各種的各種mri裝置。

以下，在不需要特別區(qū)別的情況下，以mri裝置100為代表。此外，在本實施方式中，例如使用設mri裝置100的靜磁場方向為z方向，與之垂直的2個方向中與載置測定對稱的被檢體的床面平行的方向為x方向，另一方向為y方向的坐標系。此外，以下將靜磁場簡稱為磁場。

[mri裝置的功能結構]

圖2是本實施方式的mri裝置100的功能結構圖。如本圖所示，本實施方式的mri裝置100具備：靜磁場線圈102，其由在放置有被檢體101的空間生成靜磁場的、例如靜磁場線圈、靜磁場生成磁鐵等構成；勻場線圈104，其調整靜磁場分布；發(fā)送用高頻線圈105(以下，簡稱為發(fā)送線圈)，其對被檢體101的測量區(qū)域發(fā)送高頻磁場；接收用高頻線圈106(以下，簡稱為接收線圈)，其接收從被檢體101產生的核磁共振信號；梯度磁場線圈103，其為了向從被檢體101產生的核磁共振信號附加位置信息，分別向x方向、y方向、z方向施加梯度磁場；發(fā)送機107；接收機108；計算機109；梯度磁場用電源部112；勻場用電源部113；以及順序控制裝置(sequencecontroller)114。

靜磁場線圈102根據圖1的(a)、圖1的(b)、圖1的(c)所示的各mri裝置100、120、130的結構，采用各種方式。梯度磁場線圈103及勻場線圈104分別由梯度磁場用電源部112及勻場用電源部113驅動。另外，在本實施方式中，以發(fā)送線圈105和接收線圈106分別使用不同的線圈的情況為例進行說明，但也可以由兼作發(fā)送線圈105和接收線圈106的功能的1個線圈構成。發(fā)送線圈105照射的高頻磁場由發(fā)送機107生成。將接收線圈106檢測出的核磁共振信號通過接收機108發(fā)送給計算機109。

順序控制裝置114控制梯度磁場線圈103的驅動用電源即梯度磁場用電源部112、勻場線圈104的驅動用電源即勻場用電源部113、發(fā)送機107以及接收機108的動作，并控制梯度磁場、高頻磁場的施加定時以及核磁共振信號的接收定時?？刂频臅r序圖被稱作脈沖時序圖，根據測量被預先設定，存儲在后述的計算機109所具備的存儲裝置等中。

計算機109控制mri裝置100整體的動作，并且對接收的核磁共振信號進行各種運算處理。在本實施方式中，生成至少一個以上的回波時間的復圖像或相位圖像、磁場分布、磁化率分布等。計算機109是具備cpu、存儲器、存儲裝置等的信息處理裝置，顯示器110、外部存儲裝置111、輸入裝置115被連接到計算機109。

顯示器110是將通過運算處理得到的結果等顯示給操作員的接口。輸入裝置115是操作員用于輸入在本實施方式中實施的測量或運算處理所需要的條件、參數等的接口。在本實施方式的輸入裝置115中，用戶可以輸入要測量的回波數、基準的回波時間、回波間隔等測量參數。外部存儲裝置111與存儲裝置保持計算機109執(zhí)行的各種運算處理所使用的數據、通過運算處理得到的數據、所輸入的條件、參數等。

[計算機的功能框圖]

如上所述，在本實施方式中，以基于磁場和磁化率的關系式的制約條件，反復運算平滑化處理，從而能夠簡單且高精度地計算出磁化率分布。接著，對實現這些的本實施方式的計算機109的功能結構進行說明。

圖3是表示本實施方式的計算機109的功能框圖。本實施方式的計算機109具備：測量控制部310，其測量與高頻磁場脈沖的照射對應地從被檢體產生的核磁共振信號(回波信號)作為復信號；圖像重構部320，其根據由測量控制部310測量出的復信號重構像素值為復數的復圖像；磁場分布計算部330，其根據由圖像重構部320重構的復圖像生成相位圖像，并根據該相位圖像計算出磁場分布；以及磁化率分布計算部340，其根據由磁場分布計算部330計算出的磁場分布計算磁化率分布。

另外，計算機109具備cpu、存儲器和存儲裝置。并且，cpu將存儲裝置所保持的程序加載到存儲器中并執(zhí)行，由此實現計算機109實現的各種功能。將各功能處理所使用的各種數據、處理中生成的各種數據存儲到存儲裝置或外部存儲裝置111中。

此外，計算機109實現的各種功能中至少一個功能也可以由獨立于mri裝置100的信息處理裝置即可與mri裝置100收發(fā)數據的信息處理裝置來實現。并且，全部或一部分功能也可以由asic(applicationspecificintegratedcircuit，專用集成電路)、fpga(field-programmablegatearray，現場可編程門陣列)等硬件來實現。

[qsm處理的流程]

對基于本實施方式的計算機109的上述各功能的qsm方法的處理(qsm處理)的流程進行說明。圖4是用于說明本實施方式的qsm處理的流程的流程圖。

首先，測量控制部310按照預先決定的脈沖序列控制順序控制裝置114，測量預先決定的回波時間的回波信號(步驟s1001)。

之后，圖像重構部320將得到的回波信號配置于k空間上并進行傅里葉變換，來對復圖像i進行重構(步驟s1002)。

磁場分布計算部330進行根據復圖像i計算磁場分布δ的磁場分布計算處理(步驟s1003)。

之后，磁化率分布計算部340根據計算出的磁場分布δ，進行計算磁化率分布χ的磁化率分布計算處理(步驟s1004)，并將算出的磁化率分布χ顯示于顯示器110(步驟s1005)。

另外，在將磁化率分布χ顯示于顯示器110時，根據需要，除在步驟s1003中計算出的磁場分布δ等外，還可以將在qsm處理的過程中計算出的圖像顯示于顯示器110。

以下，對本實施方式的qsm處理的各處理的細節(jié)進行說明。

[測量：s1001]

測量控制部310按照根據用戶經由輸入裝置115輸入的參數設定的脈沖序列，使順序控制裝置114動作，進行得到預先決定的回波時間(te)的核磁共振信號(回波信號)的測量。順序控制裝置114按照來自測量控制部310的指示，控制mri裝置100的各部來進行測量。在本實施方式中，得到至少一個回波時間的回波信號。

在此，對測量控制部310在測量中所使用的脈沖序列的例子進行說明。在本實施方式中，例如使用gre(gradientecho，梯度回波)系的脈沖序列。該gre系的脈沖序列對因生物體組織間的磁化率差而生成的局部的磁場變化敏銳。

圖5中，作為gre系的脈沖序列的一例，示出了rssg(rf-spoiled-steady-stateacquisitionwithrewoundgradient-echo，射頻擾相穩(wěn)態(tài)采集反轉梯度回波)序列510的例子。在該圖中，rf、gs、gp、gr分別表示高頻磁場、切片梯度磁場、相位編碼梯度磁場、讀出梯度磁場。

在rssg序列510中，施加切片梯度磁場脈沖501，并且照射高頻磁場(rf)脈沖502，激發(fā)被檢體101內預定切片的磁化。接著，向磁化相位施加用于附加切片方向和相位編碼方向的位置信息的切片編碼梯度磁場脈沖503和相位編碼梯度磁場脈沖504。

施加使像素內的核磁化的相位分散的失相(dephase)用讀出梯度磁場脈沖505后，邊施加用于附加讀出方向的位置信息的讀出梯度磁場脈沖506，邊測量一個核磁共振信號(回波)507。然后，最后施加使被切片編碼梯度磁場脈沖503和相位編碼梯度磁場脈沖504失相的核磁化的相位收斂的重相位用切片編碼梯度磁場脈沖508和相位編碼梯度磁場脈沖509。

測量控制部310邊使切片編碼梯度磁場脈沖503、508(切片編碼數ks)和相位編碼梯度磁場脈沖504、509(相位編碼數kp)的強度、rf脈沖502的相位變化，邊在反復時間tr反復執(zhí)行以上的過程，為了得到1幅圖像測量必要的回波。另外，此時使rf脈沖502的相位例如每次增加117度。此外，在圖5中，連字符以下的數字表示是反復的第幾次。

另外，測量的各回波被配置于以kr、kp、ks為坐標軸的三維k空間上。此時，一個回波在k空間上占與kr軸平行的一條線。關于通過該rssg序列510得到的絕對值圖像，在將te設定得短時成為t1(縱衰減時間)強調圖像，在將te設定得長時成為反映像素內的相位分散的t2*強調圖像。

另外，在此示例了與k空間的坐標軸平行地取得數據的笛卡兒(cartesian)拍攝的一例即rssg序列510，但在本實施方式中，可以使用任意序列來取得任意k空間區(qū)域的數據。

例如，也可以使用通過一次的tr取得多個回波，并將這些配置于不同的k空間的序列。此外，也可以是多次進行測量，通過各測量取得不同的k空間區(qū)域的數據，對這些數據進行整合，得到1幅圖像的序列。此外，也可以使用在k空間旋轉狀地取得數據的徑向掃描等非笛卡兒拍攝。

并且，也可以使用根據多個te取得的不同的復圖像、相位圖像生成一幅復圖像、相位圖像的序列。例如，使用可實現根據通過多個te取得的信號得到圖像的dixon方法的脈沖序列的情況下，在分離水和脂肪的過程得到磁場分布δ。因此，在該情況下不需要后述的生成相位圖像p的處理。

[圖像重構：s1002]

接著，對上述步驟s1002的圖像重構處理進行說明。圖像重構部320將在步驟s1001測量出的回波時間的回波信號配置于k空間，并進行傅里葉變換。由此，圖像重構部320計算出各像素值為復數的復圖像i。

[磁場分布計算處理：s1003]

接著，對上述步驟s1003的磁場分布計算部330進行的磁場分布計算處理進行說明。本實施方式的磁場分布計算部330根據圖像重構部320重構的復圖像i計算出磁場分布δ。計算出的磁場分布δ是取了因生物體組織間的磁化率差產生的局部的磁場變化的磁場分布。

磁場分布計算部330，首先根據復圖像i計算出相位圖像p，并對計算出的相位圖像p進行用于去除背景相位的各種相位圖像處理。然后，根據相位圖像處理后的相位圖像p取得磁場分布δ。

在相位圖像處理中，去除源于磁化率以外的相位分量，提取源于磁化率的相位分量。具體地，在計算出的相位圖像p中，進行去除超過－π至π的范圍后折返的相位的相位展開(unwrap)處理、和去除因生物體形狀產生的整體(global)的相位變化的整體相位變化去除處理。

為了實現這些，本實施方式的磁場分布計算部330，如圖3所示，具備：根據復圖像i計算出相位圖像p的相位圖像計算部331、進行相位展開處理的相位展開處理部332、進行整體相位變化去除處理的整體相位變化去除部333、和將相位處理后的相位圖像p變換為磁場分布δ的磁場分布變換部334。

相位圖像計算部331求出復圖像i的各像素的像素值的偏角，并將該偏角設為各個像素值來計算出相位圖像p。即，相位圖像p是以復圖像i的各像素的復數的相位分量為像素值的圖像。此處，使用復圖像i的像素值i(r)，通過式(5)來計算出像素r中的相位圖像p的像素值p(r)。

p(r)＝arg{i(r)}···(5)

相位展開處理部332對相位圖像p進行用于去除超過－π至π的范圍后折返的相位的相位展開處理。在相位圖像p中的部分區(qū)域，超過－π至π的范圍的相位值在－π至π的范圍內折返。因此，為了在拍攝部位(例如，頭部等)的全部區(qū)域求出準確的相位，需要對這些折返進行修正。在本實施方式中，例如使用區(qū)域擴大法等，對在－π至π的范圍折返的相位值進行修正。

整體相位變化去除部333在相位展開處理后的相位圖像p中去除因生物體形狀而產生的整體的相位變化。

相位圖像p的各像素值是取決于攝像部位的形狀等而產生的靜磁場不均勻而引起的整體的相位變化、組織間的磁化率變化引起的局部的相位變化的和。整體的相位變化和局部的相位變化分別對應于相位圖像p的k空間上的低頻分量和高頻分量。通過本處理，去除整體的相位變化，并提取組織間的磁化率變化引起的局部的相位變化。

本實施方式的整體相位變化去除部333，首先對拍攝的三維圖像(原圖像)，針對每個二維圖像分別實施低通濾波處理，并計算出低分辨率圖像，作為整體相位變化去除處理。然后，將原圖像用低分辨率圖像進行復數除法運算。由此，從原圖像去除低分辨率圖像所包含的整體的相位變化。

另外，去除整體的相位變化的方法包括各種方法。例如，有將三維圖像以低階多項式進行擬合來提取整體的相位變化，并將這些從原圖像中減去的方法等。在整體相位變化去除中，也可以使用這樣的其他的方法。

另外，相位圖像處理的相位展開處理和整體相位變化去除處理僅為一例，并不限定于此。此外，也可以省略這2個處理中的任一個處理。此外，各處理的處理順序是任意的。此外，也可以不進行這2個處理。

磁場分布變換部334使用式(1)，將相位圖像處理后的相位圖像p變換成磁場分布δ。

[磁場分布計算處理的流程]

對基于本實施方式的磁場分布計算部330的各功能的磁場分布計算處理的流程的一例進行說明。圖6是本發(fā)明實施方式的磁場分布計算處理的處理流程。

相位圖像計算部331根據復圖像i計算出相位圖像p(步驟s1101)。

然后，磁場分布計算部330對相位圖像p進行各種相位圖像處理，來得到相位處理后的相位圖像。在本實施方式中，例如相位展開處理部332對計算出的相位圖像p實施相位展開處理(步驟s1102)，整體相位變化去除部333對相位展開處理后的相位圖像p進行整體相位變化去除處理(步驟s1103)。如上所述，無論這些相位圖像處理的處理順序。此外，相位圖像處理可以進行任一個，也可以完全不進行。并且，并不限定于這些。

最終，磁場分布變換部334將位相處理后的相位圖像p變換為磁場分布(步驟s1104)，結束處理。

[磁化率分布計算處理：s1004]

接著，對上述步驟s1004的磁化率分布計算部340進行的磁化率分布計算處理進行說明。本實施方式的磁化率分布計算部340使用由磁場分布計算部330計算出的磁場分布δ，來計算出磁化率分布χ。

磁化率分布計算部340利用通過反復運算使誤差函數最小化的最小化處理來計算出磁化率分布χ，其中，誤差函數是由預先決定的初始磁化率分布χ0和磁場分布δ定義的函數。此時，在最小化處理中，對于根據誤差函數計算出的更新磁化率分布χupd在每次算出該更新磁化率分布時進行平滑化。另外，更新磁化率分布χupd是將初始磁化率分布χ0根據誤差函數更新后的磁化率分布。

為了實現這些，本實施方式的磁化率分布計算部340具備：誤差函數設定部341，其設定誤差函數；更新磁化率分布計算部342，其使用所設定的誤差函數，根據磁場分布δ和初始磁化率分布χ0計算更新磁化率分布χupd；平滑化處理部343，其對更新磁化率分布χupd進行平滑化，并生成合成磁化率分布χadd；以及判定部346，其進行反復運算的結束判定，并且將初始磁化率分布χ0更新為計算出的合成磁化率分布χadd。

此外，平滑化處理部343具備：濾波處理部344，其對更新磁化率分布χupd應用預先決定的平滑化濾波器，來生成濾波處理后的磁化率分布(處理后磁化率分布)χsmt；以及合成部345，其對更新磁化率分布χupd和處理后磁化率分布χsmt進行合成，生成合成磁化率分布χadd。

[誤差函數設定部]

如式(2)所示，磁場分布δ是在空間上對磁化率分布χ進行卷積積分而得的分布。以相位圖像p內的全部像素為對象，因此通過向量和矩陣表現式(2)時，表現為以下的式(6)。

δ＝fdχ···(6)

其中，δ為具有全部像素數n的大小的磁場分布的列向量，χ為計算出的磁化率分布的列向量，此外，fd為具有n×n大小，相當于對磁化率分布χ的卷積運算的矩陣。

本實施方式的更新磁化率分布計算部342根據式(6)，通過最小二乘法來求出磁化率分布χ。此時，導入誤差函數，并使其最小化來求出解。

本實施方式的誤差函數設定部341設定更新磁化率分布計算部342計算更新磁化率分布時所使用的誤差函數。例如，通過以下的式(7)來表示所設定的誤差函數e(χ0)。

[公式7]

其中，w表示對角分量具有根據復圖像i的相位值分量計算出的加權向量的矩陣(以下，稱為加權w),fd是表示卷積積分的算子，∥a∥2²表示任意向量a的l2范數(歐幾里得范數)的平方。

即，誤差函數對由磁場分布δ和初始磁化率分布χ0進行卷積運算而得的結果決定的函數賦予加權w后的函數。如上所述根據復圖像的相位值分量(相位圖像p)，計算出該加權w。

具體地，誤差函數設定部341根據相位圖像p計算出各像素的相位偏差，根據所計算出的相位偏差的大小來計算出加權w。例如，相位偏差越大，則將加權w計算得越小。

對相位偏差例如使用標準偏差等統(tǒng)計指標，該標準偏差為使用計算相位偏差的對象像素的邊緣區(qū)域的多個像素的像素值(相位值)而計算出的偏差。

另外，使用的相位圖像p既可以是根據復圖像i得到的圖像，也可以是對根據復圖像i得到的相位圖像p實施上述各種相位圖像處理而得到的圖像。

[更新磁化率計算部]

如上所述，更新磁化率分布計算部342使用誤差函數e(χ0)，從初始磁化率分布χ0計算出更新磁化率分布χupd。在此，更新磁化率分布計算部342計算出誤差函數e(χ0)的梯度使用計算出的勾配通過以下的式(8)計算出更新磁化率分布χupd。

另外，在本實施方式中，最初設定的初始磁化率分布χ0可以被設定為任意的分布。例如，可以將全部值為0的分布設為磁化率分布χ0。

[平滑化處理部；濾波處理部]

平滑化處理部343使用平滑化濾波器對更新磁化率分布χupd進行平滑化，得到處理后磁化率分布χsmt。其中，濾波處理部344使用邊緣保存平滑化濾波器作為平滑化濾波器，與更新磁化率分布χupd中的各像素值的局部的磁化率變動對應地適當地進行線性平滑化處理。

具體地，將更新磁化率分布χupd中的任意的像素位置設為r時，若將其周邊的n×n(n為1以上的整數)的局部區(qū)域的平均值設為μ(r)，將方差設為σ(r)²時，位置r處的邊緣保存平滑化濾波器應用后的磁化率分布即處理后磁化率分布χsmt的像素值χsmt(r)用以下的式(9)來表示

[公式9]

其中，ν²為在磁化率分布χ混入了相當于噪聲方差的參數。在本實施方式中，采用更新磁化率分布χupd的各像素中的局部方差σ(r)²的平均值作為該參數。

根據在本實施方式中應用的邊緣保存平滑化濾波器，如式(9)所示，在局部方差小的區(qū)域平滑化較強，成為接近局部平均值，因此能夠維持磁化率的定量值。此外，在局部方差大的區(qū)域平滑化變弱，保存包含微細結構的邊緣部分。因此，通過向更新磁化率分布χupd應用該邊緣保存平滑化濾波器，能夠降低磁化率分布的噪聲，邊維持定量值，邊自動提取并保存細微靜脈或磁化率變化大的組織間的邊緣。

另外，在本實施方式中，使用了上述的邊緣保存平滑化濾波器，但在平滑化處理中使用的濾波器并不限定于此。例如，也可以使用中值濾波器或雙邊濾波器等公知的邊緣保存平滑化濾波器。

[平滑化處理部；合成部]

合成部345對更新磁化率分布χupd和處理后磁化率分布χsmt進行合成來計算出合成磁化率分布χadd。為了防止噪聲的增加，且維持磁化率分布的微細結構而進行合成。例如，以向噪聲的影響大的k空間的魔角近旁插入被邊緣保存平滑化后的處理后磁化率分布χsmt的傅里葉分量的方式進行合成。

合成部345將更新磁化率分布χupd和處理后磁化率分布χsmt變換為k空間數據，分別對變換后的更新磁化率分布χupd和處理后磁化率分布χsmt的k空間數據加權并進行加法運算，來得到合成磁化率分布χadd。

具體地，合成部345首先將更新磁化率分布χupd和處理后磁化率分布χsmt通過傅里葉變換變換為k空間數據。然后，對變換后的各k空間數據(傅里葉分量進行加權以便達成上述目的，對加權后的各k空間數據進行加法運算。然后，計算后，通過傅里葉逆變換得到合成磁化率分布χadd。

通過生成權重圖像g，并對k空間數據乘以權重圖像g來進行k空間的加權。此時，在與根據更新磁化率分布χupd得到的k空間數據相乘的權重圖像gupd中，例如將向魔角區(qū)域近旁的數據賦予的權重設定得比向其他區(qū)域的數據賦予的權重小。此外，在與根據處理后磁化率分布χsmt得到的k空間數據相乘的權重圖像gsmt中，相反，將向魔角區(qū)域近旁的數據賦予的權重設定得比向其他區(qū)域的數據賦予的權重大。

這是因為更新磁化率分布χupd的噪聲比較多，因此降低魔角區(qū)域近旁的數據的貢獻，來防止噪聲的增加。此外，處理后磁化率分布χsmt因噪聲降低，使魔角區(qū)域近旁的數據的貢獻增加，來維持微細結構。

作為滿足這些條件的權重圖像gupd(k)、gsmt(k)，可以使用點偶極子磁場的傅里葉分量d(k)和調整魔角區(qū)域大小的調整參數b，例如可以使用由以下的式(10)定義的權重圖像。另外，k表示k空間的數據位置。

[公式10]

在上述式(10)中，在更新磁化率分布χupd中，將傅里葉分量的值小于預定值(b)的區(qū)域設為魔角區(qū)域。

另外，在調整參數b為0的情況下，gupd(k)成為1，gsmt(k)成為0。因此，不實施邊緣保存平滑化處理，而成為使式(6)的誤差函數最小化的簡單的最小二乘法。

[判定部]

在通過反復運算使誤差函數e(χ0)最小化的最小化處理中，判定部346進行反復運算的結束判定。在判定為繼續(xù)進行處理，即反復的情況下，將在該時間點得到的最新的合成磁化率分布χadd設定為初始磁化率分布χ0，來更新初始磁化率分布χ0。另外，在判定為結束時，輸出該時間點的最新的合成磁化率分布χadd作為磁化率分布。

在本實施方式中，判定部346在計算更新磁化率分布χupd后，使用該更新磁化率分布χupd和該時間點的初始磁化率分布χ0，來判定是否結束。判定部346例如使用以下的式(11)所示的評價函數f(χupd)，在其值小于預先決定的閾值ε的情況下，判定為結束。

[公式11]

[磁化率分布計算處理的流程]

接著，對本實施方式的磁化率分布計算部340進行的磁化率分布計算處理的流程進行說明。圖7是本實施方式的磁化率分布計算處理的處理流程。此外，圖8是用于說明磁化率分布計算處理的說明圖。

在本實施方式中，通過反復運算進行使利用預先決定的初始磁化率分布χ0和磁場分布δ定義的誤差函數最小化的最小化處理來計算出磁化率分布χ。此時，在最小化處理中，對根據誤差函數計算出的更新磁化率分布進行平滑化處理，得到合成磁化率分布，通過該合成磁化率分布更新初始磁化率分布，反復進行處理。具體地，進行以下的處理。

磁化率分布計算部340設定初始磁化率分布χ0410(步驟s1201)。然后，誤差函數設定部341利用初始磁化率分布χ0410和磁場分布δ400設定誤差函數e(χ0)(步驟s1202)。

然后，更新磁化率分布計算部342計算出誤差函數的梯度(步驟s1203)，按照上述式(8)，對初始磁化率分布χ0進行更新，并計算出更新磁化率分布χupd420(步驟s1204)。

判定部346使用上述式(11)的評價函數f(χupd)來進行結束判定(步驟s1205)。

在步驟s1205中，判定為繼續(xù)的情況下，濾波處理部344對更新磁化率分布χupd420應用預先決定的平滑化濾波器，計算出處理后磁化率分布χsmt430(步驟s1206)。

合成部345對更新磁化率分布χupd420和處理后磁化率分布χsmt430進行合成來得到合成磁化率分布χadd440(步驟s1207)。對更新磁化率分布χupd420和處理后磁化率分布χsmt430分別進行傅里葉變換(ft)而得到的k空間數據421、431乘以上述式(10)的權重圖像422、432并進行加法運算，對其加權相加的結果進行傅里葉逆變換(ift)來進行合成。

判定部346將合成磁化率分布χadd440設定為初始磁化率分布χ0410(步驟s1208)，返回到步驟s1202，結束處理。

在此，在步驟s1205中，在判定為結束的情況下，判定部346將在該時間點得到的最新的合成磁化率分布χadd設為磁化率分布計算部340輸出的磁化率分布χ(步驟s1209)，結束處理。

另外，在步驟s1205中，在沒有算出合成磁化率分布χadd的情況下，即，通過第1次的運算判斷為結束的情況下，通過上述方法根據該時間點的更新磁化率分布χupd計算出合成磁化率分布χadd，并將計算出的合成磁化率分布χadd設為要輸出的磁化率分布χ。

＜磁化率分布計算處理流程的變形例＞

另外，判定部346也可以使用合成磁化率分布χadd來判定是否結束。即，計算出合成磁化率分布χadd時，也可以使用該合成磁化率分布χadd和該時間點的初始磁化率分布χ0，來判定是否結束。

此時，作為判定結束的評價函數，可以使用以下的f(χadd)。

[公式12]

圖9表示此時的磁化率分布計算部340的各部的磁化率分布計算處理流程。此外，圖10是用于說明本處理的說明圖。

磁化率分布計算部340設定初始磁化率分布χ0410(步驟s1301)。然后，誤差函數設定部341利用初始磁化率分布χ0410和磁場分布δ400設定誤差函數e(χ0)(步驟s1302)。

然后，更新磁化率分布計算部342計算出誤差函數的梯度(步驟s1303)，按照上述式(8)，計算出更新磁化率分布χupd420(步驟s1304)。

濾波處理部344對更新磁化率分布χupd420應用預先決定的平滑化濾波器，計算出處理后磁化率分布χsmt430(步驟s1305)。

合成部345對更新磁化率分布χupd420和處理后磁化率分布χsmt430進行合成來得到合成磁化率分布χadd440(步驟s1306)。對更新磁化率分布χupd420和處理后磁化率分布χsmt430分別進行傅里葉變換(ft)而得到的k空間數據421、431乘以上述式(10)的權重圖像422、432并進行加法運算，對其加權相加的結果進行傅里葉逆變換(ift)來進行合成。

接著，判定部346使用上述式(12)的評價函數f(χadd)來進行結束判定(步驟s1307)。

在步驟s1307中，在判定為處理繼續(xù)的情況下，磁化率分布計算部340將合成磁化率分布χadd440設定為初始磁化率分布χ0410(步驟s1308)，返回到步驟s1302，重復處理。

在此，在步驟s1307中，在判定為結束的情況下，將在該時間點得到的最新的合成磁化率分布χadd設為磁化率分布計算部340輸出的磁化率分布χ(步驟s1309)，結束處理。

通過合成磁化率分布χadd進行判定時，與通過更新磁化率分布χupd進行判定時相比，能夠得到濾波處理的效果大的圖像。另一方面，通過更新磁化率分布χupd進行判定時，與通過合成磁化率分布χadd進行判定時相比，能夠以與實際的磁化率分布接近的數據進行判定，因此能夠在更確切的范圍內進行判定。

如上所述，本實施方式的mri裝置具備：測量控制部310，其根據高頻磁場脈沖的照射將從被檢體101發(fā)生的核磁共振信號測量為復信號；圖像重構部320，其根據上述復信號重構像素值為復數的復圖像i；磁場分布計算部330，其根據上述復圖像i生成相位圖像p，根據該相位圖像p計算出磁場分布δ；以及磁化率分布計算部340，其根據上述磁場分布δ計算出磁化率分布χ，上述磁化率分布計算部340進行通過反復運算使利用預先決定的初始磁化率分布χ0和上述磁場分布δ定義的誤差函數e(χ0)最小化的最小化處理來計算出上述磁化率分布χ，并且，在該最小化處理中，每次算出更新磁化率分布χupd時，對根據上述誤差函數e(χ0)計算出的該更新磁化率分布χupd進行平滑化。

上述磁化率分布計算部340具備：更新磁化率分布342，其使用上述誤差函數e(χ0)，根據上述初始磁化率分布χ0計算出上述更新磁化率分布χupd；平滑化處理部343，其對上述更新磁化率分布χupd進行平滑化來生成合成磁化率分布χadd；以及判定部346，其進行上述反復運算的結束判定，并且，將上述合成磁化率分布χadd設定為上述初始磁化率分布χ0，上述平滑化處理部343具備：濾波處理部344，其對上述更新磁化率分布χupd應用預先決定的平滑化濾波器，得到處理后磁化率分布χsmt；以及合成部345，其對上述更新磁化率分布χupd和上述處理后磁化率分布χsmt進行合成，來生成合成磁化率分布χadd，上述判定部346也可以將判定為結束時的上述合成磁化率分布χadd設為磁化率分布χ。

并且，上述誤差函數e(χ0)被設定為向對根據上述磁場分布δ和上述初始磁化率分布χ0進行卷積運算而得的結果決定的函數賦予加權w，上述加權w也可以基于根據上述復圖像i獲得的相位圖像p計算出。

這樣，在本實施方式中，以基于磁場和磁化率的關系式的制約條件，對平滑化處理進行反復運算，從而能夠計算出高精度的磁化率分布。

即，在本實施方式中，在使誤差函數最小化的最小化處理的過程中，對所獲得的更新磁化率分布χupd進行平滑化，反復進行處理。在平滑化過程中，在k空間上進行加權來合成通過平滑化濾波器進行平滑化而得的處理后磁化率分布χsmt和使用誤差函數的梯度計算出的更新磁化率分布χupd，以便能夠降低噪聲且維持微細結構。

這樣，根據本實施方式，通過平滑化濾波器進行的濾波，能夠使噪聲部分的數據平均化，因此能夠抑制所計算出的磁化率分布χ中的噪聲。并且，通過合成處理，在該磁化率分布χ中能夠維持微細結構。

根據本實施方式，如現有方法那樣，不使用正規(guī)化項，因此可以通過簡單的處理獲得抑制噪聲且維持微細結構的磁化率分布χ。即，根據本實施方式，不需要進行用于設定恰當的正規(guī)化參數的龐大的計算。

此外，根據本實施方式，當解決不良條件問題時，不需要對魔角區(qū)域近旁的值進行舍位。因此，還能夠抑制因舍位而產生的偽影，還可以提高計算的磁化率分布的精度。

并且，為了將更新磁化率分布χupd設為能夠抑制噪聲且維持微細結構的磁化率分布而使用的處理后磁化率分布χsmt、以及進行反復運算時更新初始磁化率分布χ0的合成磁化率分布χadd均是根據更新磁化率分布χupd計算出的。即，是取決于相位信息的數據。

在本實施方式中，使用這樣的取決于相位信息的數據來更新在最小化處理的反復運算中計算出的磁化率分布，因此具有難以產生所測量的相位信息和所計算出的磁化率分布的不匹配的問題。

并且，根據本實施方式，根據相位圖像計算出誤差函數內的加權w。相位圖像與絕對值圖像相比相位偏差的反應精度高，因此能夠獲得高精度地反應了該相位偏差的加權w。然后，通過使用了該加權w的誤差函數計算出磁化率分布χ，因此能夠更高精度地計算出磁化率分布。

如以上說明的那樣，根據本實施方式，通過簡單的處理，邊維持磁場和磁化率的關系，邊反復進行保存微細結構等邊緣信息的平滑化處理，從而能夠降低偽影和噪聲，并且提高計算的磁化率分布的準確度和精度。

＜合成時權重圖像的變形例＞

另外，合成部345對更新磁化率分布χupd和處理后磁化率分布χsmt進行合成時所使用的權重圖像g并不限定于上述的權重圖像。例如，也可以分別將gupd(k)和gsmt(k)設為固定值。即，也可以對各k空間數據內的全部數據設定一樣的權重。

通過設定一樣的權重，容易控制最終得到的圖像的snr。

＜誤差函數內的權重的變形例＞

此外，在上述實施方式中，根據復圖像i的相位值分量計算出由誤差函數設定部341設定的誤差函數e(χ0)內的加權w，但并不限定于此。例如，可以根據復圖像i的絕對值分量計算。根據復圖像i的絕對值分量計算，能夠以較少的處理簡單地計算出權重。

附圖標記說明

100：mri裝置、101：被檢體、102：靜磁場線圈、103：梯度磁場線圈、104：勻場線圈、105：發(fā)送線圈、106：接收線圈、107：發(fā)送機、108：接收機、109：計算機、110：顯示器、111：外部存儲裝置、112：梯度磁場用電源部、113：勻場用電源部、114：順序控制裝置、115：輸入裝置、120：mri裝置、130：mri裝置、310：測量控制部、320：圖像重構部、330：磁場分布計算部、331：相位圖像計算部、332：相位展開處理部、333：整體相位變化去除部、334：磁場分布變換部、340：磁化率分布計算部、341：誤差函數設定部、342：更新磁化率分布計算部、343：平滑化處理部、344：濾波處理部、345：合成部、346：判定部、400：磁場分布、410：初始磁化率分布、420：更新磁化率分布、421：更新磁化率分布的k空間數據、422：對更新磁化率分布進行乘法運算的權重圖像、430：處理后磁化率分布、431：處理后磁化率分布的k空間數據、432：對處理后磁化率分布進行乘法運算的權重圖像、440：合成磁化率分布、501：切片梯度磁場脈沖、502：rf脈沖、503：切片編碼梯度磁場脈沖、504：相位編碼梯度磁場脈沖、505：讀出梯度磁場脈沖、506：讀出梯度磁場脈沖、507：回波、508：切片編碼梯度磁場脈沖、509：相位編碼梯度磁場脈沖、510：rssg序列。