基于圖形處理單元加速的磁共振重建算法研究進(jìn)展

劉 燚 張 寧

在當(dāng)今磁共振成像(magnetic resonance imaging，MRI)的應(yīng)用中，其掃描速度已成為制約其發(fā)展的一大瓶頸，緩慢的掃描速度導(dǎo)致運(yùn)動(dòng)偽影的產(chǎn)生，也使動(dòng)態(tài)成像難以實(shí)施。20世紀(jì)90年代，研究者為提高掃描速度，通過(guò)高速切換梯度磁場(chǎng)來(lái)實(shí)現(xiàn)快速成像序列，但是過(guò)度密集的射頻脈沖序列導(dǎo)致人體內(nèi)射頻能量聚積，對(duì)人類的健康造成威脅，因此依賴梯度場(chǎng)切換速度來(lái)提高成像速度基本上達(dá)到極限[1]。

有研究已將提高成像速度的重點(diǎn)放在減少采樣數(shù)據(jù)，改進(jìn)重建算法的領(lǐng)域，如螺旋數(shù)據(jù)采樣網(wǎng)格化數(shù)據(jù)[2]重建，并行成像[3]等方法已經(jīng)為業(yè)內(nèi)人士所熟知。但是，這些早已成熟的算法重建速度卻難以接受，主要因?yàn)檫@些算法中包含大量的浮點(diǎn)運(yùn)算，當(dāng)前常見(jiàn)的中央處理器(central processing unit，CPU)的浮點(diǎn)運(yùn)算能力均不盡如人意，密集運(yùn)算性能的提高范圍有限，即使利用多線程并行處理技術(shù)，在微型計(jì)算機(jī)上不能帶來(lái)很明顯的速度提升。雖然超級(jí)計(jì)算機(jī)和集群計(jì)算在浮點(diǎn)運(yùn)算方面有著出色的表現(xiàn)，但是對(duì)于研究者和應(yīng)用者這種選擇過(guò)于昂貴。

1 CPU與圖形處理單元的結(jié)構(gòu)差異

由于CPU和圖形處理單元(graphics processing unit，GPU)在微結(jié)構(gòu)上的巨大差異，CPU中大量的晶體管用作高速緩存(cache)、邏輯控制單元(Control)，只有少量的用作計(jì)算的算術(shù)邏輯單元(arithmetic logic unit，ALU)。而GPU則把更多的晶體管用作了計(jì)算單元，只有少量晶體管用作了高速緩存(Cache)和邏輯控制單元(Control)，這使得GPU比CPU更適合完成密集計(jì)算任務(wù)。因此在浮點(diǎn)運(yùn)算能力方面，圖形處理器(graphics processing units，GPUs)正逐步表現(xiàn)出優(yōu)越的性能(如圖1所示)。

圖1 GPU和CPU微結(jié)構(gòu)示意圖

早在2003年，中央處理器(central processing units，CPUs)和GPUs關(guān)于浮點(diǎn)運(yùn)算的能力已經(jīng)開(kāi)始進(jìn)行了比較，2006年，美國(guó)英偉達(dá)[4]的G80 GPU在浮點(diǎn)運(yùn)算能力方面已經(jīng)最大達(dá)到英特爾Core 2 Duo的60倍。隨著GPU技術(shù)的不斷發(fā)展，在每秒浮點(diǎn)運(yùn)算的次數(shù)以及帶寬的表現(xiàn)，GPU都表現(xiàn)出很大的優(yōu)越性，如此巨大的計(jì)算效率的提升，使研究者開(kāi)始將醫(yī)學(xué)成像的算法轉(zhuǎn)向GPUs實(shí)現(xiàn)(如圖2所示)。

圖2 CPU和GPU浮點(diǎn)數(shù)運(yùn)算能力對(duì)比圖

2 基于計(jì)算統(tǒng)一設(shè)備架構(gòu)的GPU通用計(jì)算

近年來(lái)，研究者已經(jīng)嘗試將MRI的重建算法移植到DirectX和OpenGL等圖形化的語(yǔ)言，而且在整體性能上的提升也已經(jīng)表現(xiàn)出巨大的潛力[5]。為了進(jìn)一步鼓勵(lì)使用者將GPUs作為重要的計(jì)算引擎，英偉達(dá)公司在2007年發(fā)布了計(jì)算統(tǒng)一設(shè)備架構(gòu)(compute unified device architecture，CUDA)，并向GPU的開(kāi)發(fā)者提供了標(biāo)準(zhǔn)化的體系模型和框架。

目前的CUDA中包括cuSPARSE、cuFFT、cuBLAS和cuDPP等常用函數(shù)庫(kù)。cuSPARSE主要針對(duì)稀疏矩陣操作的線性計(jì)算庫(kù)；cuBLAS是針對(duì)傳統(tǒng)線代計(jì)算的函數(shù)庫(kù)，只支持稠密矩陣和向量的使用；cuFFT庫(kù)是利用GPU進(jìn)行快速傅里葉變換(fast Fourier transform，F(xiàn)FT)的函數(shù)庫(kù)；cuDPP庫(kù)提供了很多基本的常用并行操作函數(shù)，如排序、搜索等基本運(yùn)算組件，用于快速搭建并行計(jì)算程序。

在CUDA中，所有的線程由線程格表示，線程格分為若干個(gè)線程塊，每個(gè)線程塊有若干個(gè)線程。在具體運(yùn)算時(shí)，每個(gè)核函數(shù)針對(duì)不同的數(shù)據(jù)同時(shí)被大量的線程執(zhí)行，這些線程被組合成具有相同大小的線程塊，具有相同大小，執(zhí)行同一個(gè)核函數(shù)的線程塊組成一維或二維的網(wǎng)格，因此CUDA非常適合需要大規(guī)模并行計(jì)算的場(chǎng)景(如圖3所示)。

3 并行采集重建算法的GPU加速

3.1 算法簡(jiǎn)介

圖3 CUDA線程結(jié)構(gòu)模型圖

全面自動(dòng)校準(zhǔn)部分并行采集(generalized autocalibrating partially parallel acquisitions，GRAPPA)[6-7]概念是Griswold在2002年所提出，相比之前的基于K空間的并行成像算法，GRAPPA采用了多區(qū)塊重建的方法，一個(gè)線圈中未采樣相位編碼由多個(gè)線圈、多個(gè)區(qū)塊中已采集相位編碼進(jìn)行線性求和得到。

GRAPPA重建算法是一種更通用的基于K空間域的圖像重建方法，該算法不再將各個(gè)線圈的數(shù)據(jù)直接擬合到全K空間，而是擬合到每個(gè)線圈的自動(dòng)校準(zhǔn)信號(hào)(auto-calibration signal，ACS)行，然后求出權(quán)系數(shù)wi(m)，S是重建通道j的第ky+mΔky行，x列的數(shù)據(jù)，其計(jì)算為公式1[8]：

式中R為加速因子，Nb是用于重建的block的大小，Nc是線圈數(shù)，W為權(quán)重系數(shù)。

在GRAPPA重建算法的整個(gè)過(guò)程中，大量的計(jì)算消耗主要在W權(quán)重系數(shù)矩陣的求取上，而W的計(jì)算步驟為：利用ACS行數(shù)據(jù)建立矩陣A，B；對(duì)矩陣A，B進(jìn)行歸一化，按照公式2、公式3、公式4進(jìn)行計(jì)算：

式中U及V均為中間計(jì)算結(jié)果；W為系數(shù)矩陣。

通過(guò)分析具體的每一個(gè)計(jì)算步驟的特點(diǎn)，對(duì)算法的結(jié)構(gòu)進(jìn)行調(diào)整，以更好的利用CUDA對(duì)其進(jìn)行加速。

3.2 GRAPPA算法的使用CUDA加速

(1)利用CUDA進(jìn)行矩陣相乘運(yùn)算。由于在GRAPPA算法中，復(fù)數(shù)矩陣的乘法運(yùn)算是非常消耗時(shí)間的運(yùn)算，而在該運(yùn)算中，可以使用GPU運(yùn)算來(lái)解決這個(gè)問(wèn)題，提高復(fù)數(shù)矩陣相乘的運(yùn)算速度。調(diào)用cuBLAS中cublasCgemm()函數(shù)，實(shí)現(xiàn)單精度復(fù)數(shù)矩陣運(yùn)算核函數(shù)。

(2)利用CUDA實(shí)現(xiàn)矩陣轉(zhuǎn)置。從算法實(shí)現(xiàn)過(guò)程可以看出，在第二步中可以使用線程同步進(jìn)行操作，其算法具體實(shí)現(xiàn)步驟如圖4所示。

圖4 矩陣轉(zhuǎn)置計(jì)算流程圖

(3)矩陣求逆。在矩陣的求逆運(yùn)算中使用的是柯列斯基分解，而該算法并不適用于并行化的計(jì)算，故該算法應(yīng)在主機(jī)CPU下進(jìn)行計(jì)算。

通過(guò)上述3個(gè)部分優(yōu)化，重新調(diào)整算法結(jié)構(gòu)，在計(jì)算機(jī)上對(duì)優(yōu)化前和優(yōu)化后的兩種方法分別進(jìn)行運(yùn)算，在鑫高益公司SuperScan1.5T型MRI上掃描志愿者得到采樣數(shù)據(jù)，采集矩陣為276×276，共采集8個(gè)通道19層數(shù)據(jù)，從中取出8個(gè)通道的一層數(shù)據(jù)進(jìn)行測(cè)試。測(cè)試過(guò)程是在配有8 G內(nèi)存，Intel Core i5-3230M CPU(主頻為2.60 GHz)，Navida NVS 5400 M(核心頻率為600 MHz，最大顯存為2048 MB)的計(jì)算機(jī)上進(jìn)行(如圖5所示)。

圖5 GRAPPA運(yùn)算執(zhí)行時(shí)間對(duì)比圖

圖5顯示，在GRAPPA算法中需要大量消耗計(jì)算時(shí)間的具體環(huán)節(jié)上，使用CUDA進(jìn)行計(jì)算，利用CUDA的多線程block，以及出色的浮點(diǎn)計(jì)算能力，通過(guò)在算法中加入CUDA運(yùn)算，極大的提高了GRAPPA算法的運(yùn)算速度，從而極好的解決了GRAPPA算法運(yùn)算時(shí)間過(guò)長(zhǎng)的瓶頸，而且重建效果基本一致，GPU加速為GRAPPA算法在并行重建中的應(yīng)用提供了更好的發(fā)展前景(如圖6所示)。

圖6 CPU和GPU重建后的MRI結(jié)果圖像

并行采集的重建算法中，常規(guī)的SENSE算法也可以使用GPU進(jìn)行計(jì)算，且在其他的SENSE擴(kuò)展算法中也有GPU的應(yīng)用[9]。

4 Griding算法的GPU加速

4.1 算法概述

網(wǎng)格化(Griding)算法[10]是一種當(dāng)前非常流行的基于非笛卡爾采樣軌跡的MR圖像重建技術(shù)，射線型、螺旋型采樣[11]相比笛卡爾型采樣有其他的優(yōu)勢(shì)，如采集速度快，對(duì)運(yùn)動(dòng)偽影不敏感，如果采用此種采樣就不能直接使用FFT重建圖像，因此采用網(wǎng)格化的算法可以將非笛卡爾網(wǎng)格通過(guò)插值變成笛卡爾型，然后可使用FFT進(jìn)行圖像重建[12]。對(duì)于放射型K空間采用的重建，對(duì)一些關(guān)鍵的步驟使用CUDA加速，取得比較好的效果(如圖7、圖8所示)。

圖7 笛卡爾型K空間采樣示圖

圖8 非笛卡爾型K空間采樣示圖

4.2 算法加速

在具體算法實(shí)施過(guò)程中，其中插值這一步中使用特定的插值窗函數(shù)Kaiser-Bessel窗，在一步中由于特定計(jì)算結(jié)構(gòu)，笛卡爾網(wǎng)格中的所有元素C(Kx，Ky)，可以使用GPU進(jìn)行加速計(jì)算，其卷積插值計(jì)算為公式5：

通過(guò)使用CUDA對(duì)卷積插值方法進(jìn)行優(yōu)化，以采集到的K空間軌跡上非均勻分布采集點(diǎn)為中心，CUDA中的每個(gè)線程負(fù)責(zé)一個(gè)采樣點(diǎn)或一組采樣點(diǎn)的計(jì)算，根據(jù)采樣點(diǎn)與卷積窗內(nèi)的網(wǎng)格點(diǎn)的距離，計(jì)算采樣點(diǎn)對(duì)這些網(wǎng)格點(diǎn)的貢獻(xiàn)，相當(dāng)于把采樣點(diǎn)的值“分配”到各個(gè)網(wǎng)格點(diǎn)上。測(cè)試結(jié)果也表明，網(wǎng)格化算法使用GPU計(jì)算，相比CPU計(jì)算GPU可以達(dá)到3.5倍左右的加速，且兩者的計(jì)算結(jié)果基本一致，呼吸運(yùn)動(dòng)偽影明顯減弱，表明在使用網(wǎng)格化的相關(guān)重建算法中都可以用GPU來(lái)進(jìn)行加速[13]。在對(duì)螺旋采用的K空間進(jìn)行網(wǎng)格化，且使用GPU加速，也取得很好效果(如圖8、圖9所示)。

圖9 網(wǎng)格化算法CPU重建MRI圖像

圖10 網(wǎng)格化算法GPU重建MRI圖像

5 FFT的GPU加速

FFT是一種重要的算法，幾乎在所有的磁共振圖像重建算法中都有應(yīng)用，同樣FFT算法的快慢也將對(duì)磁共振圖像重建算法的速度有著非常重要的影響。傳統(tǒng)計(jì)算FFT通常會(huì)選用FFTW庫(kù)，或是Intel MKL，這兩個(gè)函數(shù)庫(kù)都是經(jīng)過(guò)高度優(yōu)化的數(shù)學(xué)計(jì)算庫(kù)，對(duì)于計(jì)算FFTW，這兩個(gè)庫(kù)都可以充分發(fā)揮CPU的運(yùn)算能力。但是，隨著GPU的出現(xiàn)，基于GPU的浮點(diǎn)運(yùn)算表現(xiàn)出強(qiáng)大的能力，所以在目前的磁共振圖像重建中使用基于GPU的FFT已經(jīng)廣泛應(yīng)用，尤其是英偉達(dá)公司推出CUDA以后，基于CUDA的FFT庫(kù)[14](cuFFT)使得FFT計(jì)算速度有了大幅度的提升，相對(duì)于FFTW和Intel MKL，都有很大的速度優(yōu)勢(shì)。實(shí)驗(yàn)表明，cuFFT對(duì)MKL在單精度和雙精度的計(jì)算上都有著非常明顯的優(yōu)勢(shì)，其他測(cè)試也表明，在不考慮數(shù)據(jù)拷貝的情況下，cuFFT的運(yùn)算時(shí)間遠(yuǎn)低于FFTW的運(yùn)算時(shí)間[15-18](如圖11所示)。

圖11 cuFFT與MKL單精度和雙精度運(yùn)算速度對(duì)比示圖

在SuperScan1.5T型MRI上掃描志愿者得到采樣數(shù)據(jù)，采集矩陣為276×276，共掃描17層。重建過(guò)程是在配有8 G內(nèi)存，Intel Core i5-3230M CPU(主頻為2.60 GHz)，Navida NVS 5400 M(核心頻率為600 MHz，最大顯存為2048 MB)的計(jì)算機(jī)上進(jìn)行，分別使用MKL進(jìn)行FFT運(yùn)算和cuFFT進(jìn)行運(yùn)算，測(cè)試結(jié)果表明，使用MKL計(jì)算時(shí)間為112 ms，cuFFT計(jì)算時(shí)間為72 ms，運(yùn)算時(shí)間可以減少40%，加速效果比較明顯，且計(jì)算結(jié)果基本一致。在一些磁共振重建的算法中會(huì)反復(fù)大量用到FFT運(yùn)算，所以累計(jì)減少的運(yùn)算時(shí)間比較可觀，基于GPU的FFT計(jì)算所表現(xiàn)出來(lái)的強(qiáng)大優(yōu)勢(shì)，使得該技術(shù)在磁共振重建中發(fā)揮著重要的作用(如圖12、圖13所示)。

6 展望

自從GPU出現(xiàn)以后，其在科學(xué)計(jì)算方面所展現(xiàn)出強(qiáng)大的運(yùn)算能力，隨著GPU顯卡高速發(fā)展，該技術(shù)已在醫(yī)學(xué)MRI領(lǐng)域得到了廣泛的應(yīng)用，對(duì)算法的加速也取得了非常顯著的進(jìn)步，從而使得許多優(yōu)秀的重建算法得到更好的應(yīng)用，從整體上加快MRI速度，某種程度上也提高了圖像的質(zhì)量，為醫(yī)學(xué)MRI技術(shù)添加了新的發(fā)展動(dòng)力。

圖12 MKL重建結(jié)果圖像

圖13 cuFFT重建結(jié)果圖像

對(duì)3種常見(jiàn)的磁共振重建算法使用GPU加速，驗(yàn)證了GPU的加速性能，這些只是在該領(lǐng)域的嘗試，如基于深度學(xué)習(xí)的MRI圖像重建加速、基于MRI掃描圖像的藥理分析模型、MRI實(shí)時(shí)掃描重建[19]以及大數(shù)據(jù)量的3D圖形加速重建等等都可以使用GPU加速運(yùn)算[20]。由于顯卡顯存的限制，目前只能滿足二維和部分三維圖像的重建加速，對(duì)于更高維度、更大數(shù)據(jù)的計(jì)算將是巨大的挑戰(zhàn)。但是，隨著GPU技術(shù)的進(jìn)一步發(fā)展，相信該技術(shù)在磁共振領(lǐng)域也一定會(huì)有更多更廣泛的應(yīng)用。