本發(fā)明屬于圖像處理技術(shù)領(lǐng)域，具體涉及一種基于field map對(duì)fMRI中的幾何畸變偽影進(jìn)行校正的方法，尤其涉及通過多通道圖像重建、幅值圖分割、相位圖去卷繞、field map的生成和去噪、體素位移圖生成及幾何畸變校正實(shí)現(xiàn)對(duì)功能核磁共振圖像進(jìn)行基于field map的幾何畸變校正的方法。

背景技術(shù)：

fMRI是研究腦活動(dòng)、腦功能的主要的無創(chuàng)方法之一，具有毫米級(jí)的空間分辨率。BOLD-fMRI方法的提出和發(fā)展對(duì)腦認(rèn)知功能的研究具有突破性的進(jìn)展，它已成為神經(jīng)科學(xué)探索人類大腦神經(jīng)機(jī)制的重要工具。fMRI—般指基于血氧水平依賴(blood oxygen level-dependent,BOLD)的磁共振成像，它通過測(cè)量由神經(jīng)活動(dòng)引起的腦血流和腦血氧等成分變化而造成的磁共振信號(hào)變化來反應(yīng)腦活動(dòng)。

腦功能和神經(jīng)功能研究中通常要求把低分辨率的EPI圖像配準(zhǔn)到高分辨率的MRI結(jié)構(gòu)圖像上，將EPI圖像反映的腦部功能活動(dòng)配準(zhǔn)到MRI結(jié)構(gòu)圖像反映的大腦結(jié)構(gòu)上。但是，在實(shí)際采集EPI圖像的過程中，由于磁場(chǎng)不均勻的存在，會(huì)不可避免地引入偽影(artifact)，主要表現(xiàn)為EPI圖像的幾何畸變(geometric distortion)和信號(hào)丟失(signal dropout)。幾何畸變表現(xiàn)為圖像上的像素點(diǎn)(對(duì)應(yīng)實(shí)際數(shù)據(jù)的體素點(diǎn))偏離本來的位置、被壓縮或被拉伸。EPI圖像的幾何畸變使功能磁共振圖像難以準(zhǔn)確配準(zhǔn)到結(jié)構(gòu)像上，可能造成腦激活區(qū)域的定位錯(cuò)誤。所以，在開展后續(xù)研究之前，必須對(duì)EPI圖像進(jìn)行校正，而現(xiàn)有技術(shù)中并無相關(guān)的校正方法。

技術(shù)實(shí)現(xiàn)要素：

為了克服上述現(xiàn)有技術(shù)的缺陷，本發(fā)明的目的在于基于field map對(duì)功能磁共振成像中的幾何畸變偽影進(jìn)行校正的方法，能夠簡(jiǎn)單、穩(wěn)定、高效地對(duì)功能核磁共振圖像的幾何畸變偽影進(jìn)行校正。

為了達(dá)到上述目的，本發(fā)明的技術(shù)方案為：

基于field map對(duì)fMRI中的幾何畸變偽影進(jìn)行校正的方法，具體步驟如下：

(1)、多通道相控陣線圈圖像重建，多通道相控陣線圈是由多個(gè)基本子線圈構(gòu)成，每個(gè)子線圈都接收全部腦區(qū)或其它掃描區(qū)域產(chǎn)生的磁共振信號(hào)，最后將各個(gè)線圈接收的信號(hào)進(jìn)行重建，得到完整的腦區(qū)圖像；

(2)、幅值圖分割做掩膜，對(duì)field map做掩膜(mask)運(yùn)算，去除腦組織之外的信息，對(duì)重建后的幅值圖像進(jìn)行分割，去除頭骨部分，將分割后的幅值圖像作為field map的掩膜；

(3)、相位去卷繞，采用PRELUDE算法來進(jìn)行相位去卷繞，該算法通過一種區(qū)域生長(zhǎng)的方法優(yōu)化cost函數(shù)來判定相鄰相位區(qū)域的相位差值，最后利用該差值對(duì)各個(gè)相鄰區(qū)域進(jìn)行相位恢復(fù)去卷繞；

(4)、Field map的生成，從多通道中采集到的相位圖的取值范圍為0～4096，首先對(duì)各個(gè)通道的相位圖進(jìn)行歸一化，將取值范圍歸一化到0～2π，經(jīng)過多通道的相位重建后，相位差圖的取值范圍為-π到+π，再通過計(jì)算得到的field map；

(5)、Field map去噪，對(duì)于field map的邊緣噪聲，采取腐蝕運(yùn)算將其去除；對(duì)于field map的內(nèi)部噪聲，采用卷積濾波方式進(jìn)行去噪，具體做法是先采用去尖峰和中值濾波去除尖峰噪聲，然后再用高斯平滑濾波對(duì)含噪聲的field map做進(jìn)一步的平滑；

(6)、生成體素位移圖:獲得了去噪后的field map后，依據(jù)下式計(jì)算體素位移圖(voxel shift map,VSM)：

$\frac{y_0-y}{\mathrm{\Δ}y}=(2{\mathrm{\τ}}_r+\frac{N_x}{BW})N_y\mathrm{\Δ}f(x,y)$

其中即為體素位移，Δf(x,y)即為去噪后的field map，N_y為y方向采樣點(diǎn)總個(gè)數(shù),BW讀出方向的采樣帶寬,τ_r為讀出方向轉(zhuǎn)換梯度的上升時(shí)間；

(7)、對(duì)幾何畸變圖像進(jìn)行校正:利用去噪后的field map和體素位移圖對(duì)幾何畸變的EPI圖像按體素進(jìn)行位置移動(dòng)，對(duì)移動(dòng)后產(chǎn)生的空隙進(jìn)行插值計(jì)算，得到校正的EPI圖像。

本發(fā)明的創(chuàng)新點(diǎn)在于：依據(jù)噪聲來源不同，分別采用腐蝕的幅值圖做掩膜去除邊緣噪聲和用去尖峰、中值濾波和高斯平滑方法去除內(nèi)部噪聲，兩種方法的結(jié)合有效的去除了field map的噪聲，利用field map對(duì)功能磁共振圖像的校正效果明顯，簡(jiǎn)單高效，具有很高的信噪比，能很好地改善其圖像品質(zhì)。

附圖說明

圖1是本發(fā)明的流程圖。

圖2是由八通道幅值圖重建完整的幅值圖。

圖3是由兩個(gè)回波時(shí)刻八通道相位圖重建相位差圖。

圖4利用BET對(duì)幅值圖進(jìn)行分割。左圖為分割前的幅值圖，右圖為分割后的幅值圖。

圖5對(duì)兩個(gè)TE時(shí)刻的相位圖進(jìn)行相位去卷繞的結(jié)果。左圖為帶有相位卷繞的相位圖，右圖為去卷繞后的相位圖。在相位去卷繞時(shí)利用分割的幅值圖進(jìn)行了掩膜運(yùn)算。

圖6是Field map的生成。左上圖為相位去卷繞后的相位差圖，右上圖為分割的幅值圖，利用分割的幅值圖對(duì)相位差圖做掩膜并處以兩個(gè)回波時(shí)間間隔，即獲得field map。

圖7是利用腐蝕后的幅值圖對(duì)field map做掩膜來去除field map的邊緣噪聲。左圖為原field map，可以觀察到邊緣處有明顯的噪聲；右圖為用腐蝕后的幅值圖做掩膜得到的field map。

圖8濾波去除field map的內(nèi)部噪聲結(jié)果圖。

圖9是field map生成的體素位移圖。

圖10是校正前后的EPI圖像與幅值像。左上圖為校正前的EPI圖像，右上圖為校正后的EPI圖像，下圖為沒有幾何畸變的幅值像?？梢杂^察到黃色圓圈和箭頭所示的前額葉部位的幾何畸變得到了一定的校正。

具體實(shí)施方式

下面結(jié)合附圖對(duì)本發(fā)明做詳細(xì)敘述。

本發(fā)明基于field map對(duì)功能磁共振成像中的幾何畸變偽影進(jìn)行校正方法流程如圖1所示。

(1)、首先多通道相控陣線圈圖像重建。多通道相控陣線圈是由多個(gè)基本子線圈構(gòu)成，每個(gè)子線圈都接收全部腦區(qū)(或其它掃描區(qū)域)產(chǎn)生的磁共振信號(hào)，最后將各個(gè)線圈接收的信號(hào)進(jìn)行重建，得到完整的腦區(qū)圖像。數(shù)據(jù)樣本是在SIMENSE 4T磁共振成像儀上對(duì)一名34歲的男子采集的。數(shù)據(jù)包括八通道的雙回波梯度序列(Dual Echo GR)幅值像數(shù)據(jù)和EPI數(shù)據(jù)。EPI數(shù)據(jù)是在受試者做數(shù)數(shù)任務(wù)(number task)時(shí)采集的，層面內(nèi)成像矩陣(matrix)大小為64×64，采集的層面(slice)數(shù)目為34，層厚(slice thickness)為3mm，層間距為3.45mm，翻轉(zhuǎn)角度(flip angle)為75度，回波時(shí)間TE為33ms，重復(fù)時(shí)間TR為2000ms，每像素帶寬為1445Hz/pixel，共采集170個(gè)volume。

幅值像的采集采用雙回波梯度序列(Dual Echo GR)，雙回波梯度序列的兩個(gè)回波時(shí)間TE分別為6ms、10ms，TR＝494ms，層面內(nèi)成像矩陣(matrix)大小為64×64，層面內(nèi)FOV為192×192mm，所以層面內(nèi)的像素大小為3×3mm，采集的層面(slice)數(shù)目為34，層厚(slice thickness)為3mm，層間距為3.45mm，翻轉(zhuǎn)角度(flip angle)為30度，采集帶寬為260Hz/pixel。

采用八通道的頭線圈采集了八通道的數(shù)據(jù)，每通道的數(shù)據(jù)包括幅值(magnitude)數(shù)據(jù)和相位(phase)數(shù)據(jù)，我們需要將八個(gè)通道中的圖像重建成完整的幅值圖和相位圖。我們只需要利用一個(gè)回波時(shí)刻的幅值圖，將各個(gè)通道的圖像按照下式進(jìn)行重建：

$A=\sqrt{\underset{i}{\Σ}{A_i}^2}$

式中A表示重建后的幅值圖，A_i表示第i個(gè)通道的幅值圖。重建過程如圖2所示，可以觀察到重建后的幅值圖各個(gè)腦區(qū)的信號(hào)強(qiáng)度基本一致。

我們需要理由兩個(gè)回波時(shí)刻的相位圖構(gòu)建field map，field map是由兩個(gè)回波時(shí)刻的相位差構(gòu)建的，所以我們直接按下式構(gòu)建相位差圖，結(jié)果如圖3所示。

$\mathrm{\Δ}\mathrm{\Φ}=\∠\underset{j=1}{\overset{n}{\Σ}}{e}^{i({\mathrm{\Φ}}_{2,j}-{\mathrm{\Φ}}_{1,j})}$

ΔΦ為由兩個(gè)回波時(shí)刻的八通道相位圖構(gòu)建的相位差圖，Φ_2,j為第2個(gè)回波時(shí)刻第j個(gè)通道的相位圖，Φ_1,j為第1個(gè)回波時(shí)刻第j個(gè)通道的相位圖，n為通道數(shù)，j為當(dāng)前通道。

(2)、幅值圖分割做掩膜。field map一般只提供大腦顱骨內(nèi)部腦組織的信息，所以我們一般會(huì)對(duì)field map做掩膜(mask)運(yùn)算，去除腦組織之外的信息。我們對(duì)重建后的幅值圖像進(jìn)行分割，去除頭骨部分，將分割后的幅值圖像作為field map的掩膜。

BET算法的實(shí)現(xiàn)流程如下：首先通過直方圖的方法找出頭部圖像中像素的灰度值上限和下限；接著對(duì)頭部圖像進(jìn)行閾值化計(jì)算，在閾值化后的圖像中計(jì)算頭部重心坐標(biāo)及半徑；在原頭部圖像的重心處定義一個(gè)半徑較小球面，球面由很小的三角形網(wǎng)格作為基本單元；每次增加一個(gè)三角形網(wǎng)格的頂點(diǎn)，使球面不斷向外擴(kuò)充，并進(jìn)行平滑運(yùn)算；當(dāng)球面達(dá)到頭部組織的外部邊界時(shí)，再進(jìn)行迭代時(shí)球面不會(huì)發(fā)生太大變化，此時(shí)的球面位置即為分割的邊界位置，分割完成。BET被集成到FMRIB Software Library工具可以直接使用，利用BET算法對(duì)重建的幅值像進(jìn)行分割結(jié)果如圖4所示，左圖為分割前的幅值圖，右圖為分割后的幅值圖。

(3)、進(jìn)行相位去卷繞。這里采用Jenkinson在2003年提出的PRELUDE(Phase Region Expanding Labeller for Unwrapping Discrete Estimates)算法來進(jìn)行相位去卷繞。該算法通過一種區(qū)域生長(zhǎng)的方法優(yōu)化cost函數(shù)來判定相鄰相位區(qū)域的相位差值，最后利用該差值對(duì)各個(gè)相鄰區(qū)域進(jìn)行相位恢復(fù)去卷繞。首先將相位圖內(nèi)部分成各個(gè)區(qū)域，每個(gè)區(qū)域內(nèi)部不存在相位卷繞；其次定義cost函數(shù)定義為相鄰區(qū)域邊界的相位差的平方和，通過使cost函數(shù)達(dá)到最小值可以獲得相鄰兩區(qū)域的相位差；獲得相位差后，將相鄰兩區(qū)域進(jìn)行合并，直至整個(gè)圖像合并成一個(gè)區(qū)域時(shí)，完成相位的去卷繞。該算法已被集成到FMRIB Software Library軟件包中。利用PRELUDE算法對(duì)兩個(gè)回波時(shí)刻相位圖進(jìn)行相位去卷繞的結(jié)果如圖5所示。

(4)、Field map的生成。利用磁場(chǎng)的不均勻程度ΔB(x,y)可以表征field map，同樣利用ΔB(x,y)引起的拉莫爾進(jìn)動(dòng)頻率偏移Δf(x,y)也可以表征field map：

$\mathrm{\Δ}f(x,y)=\mathrm{\γ}\mathrm{\Δ}B=\frac{\mathrm{\Δ}\mathrm{\Φ}}{\mathrm{\Δ}TE}$

式中γ代表旋磁比，Δf(x,y)表示磁場(chǎng)不均勻引起的拉莫爾進(jìn)動(dòng)頻率偏移的分布，ΔB表示磁場(chǎng)不均勻的強(qiáng)度分布，ΔΦ表示兩個(gè)回波時(shí)刻的相位差，ΔTE表示兩個(gè)回波的時(shí)間間隔。進(jìn)行計(jì)算即可獲得field map Δf(x,y)，在生成field map的過程中，利用分割出的幅值圖做掩膜，去除非大腦組織的影響，其過程如圖6所示。

(5)、Field map去噪。由圖6中的field map可以觀察到，field map的邊緣和內(nèi)部存在噪聲，噪聲的存在會(huì)影響到體素位移位置的計(jì)算精度，最終影響幾何畸變校正的準(zhǔn)確度，所以必須對(duì)field map進(jìn)行去噪處理。本發(fā)明采取基于掩膜的腐蝕運(yùn)算，首先對(duì)分割過的幅值圖進(jìn)行腐蝕，幅值圖內(nèi)部基本不存在噪聲，腐蝕運(yùn)算只會(huì)腐蝕幅值圖的邊緣，再利用腐蝕后的幅值圖對(duì)field map做掩膜運(yùn)算，就將field map的邊緣噪聲去除，如圖7右圖所示。

對(duì)于field map的內(nèi)部噪聲，我們采用卷積濾波方式進(jìn)行去噪。濾波方法有三種：去尖峰方法(despike)、中值濾波方法(median filtering)和三維高斯平滑濾波方法(3D Gaussian filtering)。先采用去尖峰和中值濾波去除尖峰噪聲，然后再用高斯平滑濾波對(duì)含噪聲的field map做進(jìn)一步的平滑，結(jié)果如圖8所示。

(6)、生成體素位移圖。獲得了去噪后的field map后，就可以依據(jù)下式計(jì)算體素位移圖(voxel shift map,VSM)：

$\frac{y_0-y}{\mathrm{\Δ}y}=(2{\mathrm{\τ}}_r+\frac{N_x}{BW})N_y\mathrm{\Δ}f(x,y)$

其中即為體素位移，Δf(x,y)即為去噪后的field map，N_y為y方向采樣點(diǎn)總個(gè)數(shù),BW讀出方向的采樣帶寬,τ_r為讀出方向轉(zhuǎn)換梯度的上升時(shí)間，去噪后的field map計(jì)算得到的體素位移圖如圖9所示。

(7)、對(duì)幾何畸變圖像進(jìn)行校正。

最后我們利用去噪后的field map和體素位移圖對(duì)幾何畸變的EPI圖像按體素進(jìn)行位置移動(dòng)，對(duì)移動(dòng)后產(chǎn)生的空隙進(jìn)行插值計(jì)算，得到校正的EPI圖像。我們利用FUGUE(FMRIB's Utility for Geometrically Unwarping EPIs)實(shí)現(xiàn)幾何畸變圖像的體素位移和插值，結(jié)果如圖10所示，左上圖為校正前的EPI圖像，右上圖為校正后的EPI圖像(采用despike濾波)，下圖為沒有幾何畸變的幅值像，可以觀察到黃色圓圈和箭頭所示的前額葉部位的幾何畸變得到了一定的校正。

綜上所述，本發(fā)明的創(chuàng)新之處在于在對(duì)field map進(jìn)行去噪的過程中，針對(duì)邊緣噪聲和內(nèi)部噪聲的不同，先后采用腐蝕幅值圖對(duì)field map做掩膜濾除邊緣噪聲，采用尖峰方法(despike)、中值濾波方法(median filtering)和三維高斯平滑濾波方法(3D Gaussian filtering)濾除內(nèi)部噪聲，對(duì)邊緣噪聲和內(nèi)部噪聲濾除的結(jié)合使用達(dá)到了很好的去噪效果。利用field map對(duì)功能磁共振圖像的校正效果明顯，簡(jiǎn)單高效，具有很高的信噪比，能很好地改善其圖像品質(zhì)。