基于GPU的快速三維醫(yī)學(xué)圖像剛性配準(zhǔn)技術(shù)

摘 要:自動(dòng)三維配準(zhǔn)將多個(gè)圖像數(shù)據(jù)映射到同一坐標(biāo)系中，在醫(yī)學(xué)影像分析中有廣泛的應(yīng)用。但現(xiàn)有主流三維剛性配準(zhǔn)算法(如FLIRT)速度較慢，2563大小數(shù)據(jù)的剛性配準(zhǔn)需要300 s左右，不能滿足快速臨床應(yīng)用的需求。為此提出了一種基于CUDA(compute unified device architecture)架構(gòu)的快速三維配準(zhǔn)技術(shù)，利用GPU(gra-phic processing unit)并行計(jì)算實(shí)現(xiàn)配準(zhǔn)中的坐標(biāo)變換、線性插值和相似性測(cè)度計(jì)算。臨床三維醫(yī)學(xué)圖像上的實(shí)驗(yàn)表明，該技術(shù)在保持配準(zhǔn)精度的前提下將速度提高了一個(gè)數(shù)量級(jí)，非常適合海量快速的臨床應(yīng)用。

關(guān)鍵詞:三維圖像剛性配準(zhǔn);統(tǒng)一計(jì)算設(shè)備架構(gòu);三維醫(yī)學(xué)圖像;GPU

中圖法分類號(hào):TP391 文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1001-3695(2010)03-1198-03

doi:10.3969/j.issn.1001-3695.2010.03.111

Fast 3D rigid medical image registration based on GPU

QIN An， XU Jian， FENG Qian-jin， MENG Xiao-lin， CHEN Wu-fan

(Key Laboratory for Medical Image Processing of Guangdong Province， Southern Medical University， Guangzhou 510515， China)

Abstract: Automatic 3D image rigid registration technique， which maps multiple images into a common coordinate space， is an important clinical image analysis tool. The state of art rigid registration techniques like FLIRT take 300s to complete a rigid registration for a 2563 data set.This makes it hard to apply this technique in time critical clinical environment.This paper proposed a parallel accelerating registration algorithm based on the new CUDA.By using the GPU to compute the coordinate transformation， interpolation and similarity metric computation， great improvement of speed has been achieved without compromising the registration accuracy. The proposed technique is valuable for time critical clinical application.

Key words:3D medical image registration;CUDA; 3D medical image; GPU

0 引言

自動(dòng)快速的剛性三維醫(yī)學(xué)圖像配準(zhǔn)可以將多組三維醫(yī)學(xué)影像通過剛性變換(平移、旋轉(zhuǎn)等)對(duì)齊到同一坐標(biāo)空間中，是三維醫(yī)學(xué)影像診斷的一種常用工具，也是三維ASM(active shape model)訓(xùn)練集構(gòu)造的一個(gè)重要步驟。醫(yī)學(xué)圖像配準(zhǔn)技術(shù)在手術(shù)導(dǎo)航、放療、評(píng)估病情的進(jìn)展或者放療的療效方面有廣泛的應(yīng)用。

隨著螺旋CT、三維超聲等三維成像設(shè)備的普及，傳統(tǒng)的配準(zhǔn)方法在處理三維數(shù)據(jù)時(shí)速度很慢，如對(duì)2563數(shù)據(jù)，國(guó)際通用的FLIRT工具包中的剛性配準(zhǔn)模塊需要300 s左右[1]，嚴(yán)重妨礙了這類算法在要求實(shí)時(shí)或快速實(shí)施的臨床環(huán)境中的應(yīng)用，如手術(shù)中圖像配準(zhǔn)、放療計(jì)劃時(shí)序配準(zhǔn)等。

目前的配準(zhǔn)方法可以分為兩類，即基于特征的配準(zhǔn)和基于灰度的配準(zhǔn)。它們主要區(qū)別在于是否包含圖像的分割或特征標(biāo)記的提取，但兩類方法的運(yùn)算量都較大，尤其是在三維醫(yī)學(xué)影像開始普及的今天。因此，如何在不損失或少損失配準(zhǔn)精度的情況下，加快大數(shù)據(jù)量的醫(yī)學(xué)配準(zhǔn)速度，一直是醫(yī)學(xué)圖像分析的研究熱點(diǎn)。許多學(xué)者提出新算法或新技術(shù)用以減少配準(zhǔn)時(shí)間，曾德賢等人[2]采用簡(jiǎn)化仿射變換模型的方式來(lái)加速圖像的配準(zhǔn)，何立新等人[3]則采用了多級(jí)分辨率的方式來(lái)加速圖像的配準(zhǔn)。

近年來(lái)，具備通用計(jì)算能力的GPU的出現(xiàn)，給普通PC上的并行計(jì)算帶來(lái)革命性的變化。其中較有代表性的是NVIDIA公司的8800及其后續(xù)系列顯卡，具備強(qiáng)大的通用計(jì)算和可編程的能力。GPU已經(jīng)在許多傳統(tǒng)的運(yùn)算量大且具備并行特點(diǎn)的圖像分析算法中得到成功的應(yīng)用，但在三維圖像配準(zhǔn)上的研究還剛剛興起。2008年，Samant等人[4]報(bào)道了GPU加速?gòu)椥耘錅?zhǔn)在自適應(yīng)放射治療中的應(yīng)用;同年Li Bo等人[5]報(bào)道了GPU加速?gòu)椥耘錅?zhǔn)在心功能評(píng)價(jià)中的應(yīng)用;國(guó)內(nèi)學(xué)者最近也開展了GPU加速的二維圖像彈性配準(zhǔn)技術(shù)的研究[6]。

本文提出了一種基于CUDA的GPU并行醫(yī)學(xué)圖像配準(zhǔn)加速技術(shù)，通過利用通用圖形顯卡的GPU并行計(jì)算能力，實(shí)現(xiàn)了三維醫(yī)學(xué)圖像配準(zhǔn)的跨數(shù)量級(jí)速度的提升，該技術(shù)在要求實(shí)時(shí)或者快速臨床環(huán)境中具有較好的臨床應(yīng)用前景。

1 基于GPU的加速醫(yī)學(xué)圖像配準(zhǔn)技術(shù)

1.1 基于CUDA的剛性配準(zhǔn)

CUDA是NVIDIA的GPU模型，它以C語(yǔ)言為基礎(chǔ)，可以按照C語(yǔ)言語(yǔ)法編寫程序[7]。CUDA的執(zhí)行模式是單指令多數(shù)據(jù)(single instruction multiple data，SIMD)模式，也就是說(shuō)CUDA進(jìn)行的是高度并行化的工作。在CUDA 的架構(gòu)下，一個(gè)程序分為兩個(gè)部分，即host 端和device端，它們之間通過PCI總線進(jìn)行數(shù)據(jù)交流。

醫(yī)學(xué)圖像配準(zhǔn)[6]就是指對(duì)一幅圖像尋求一種(一系列)空間變換，使該圖像與另一幅圖像中的對(duì)應(yīng)點(diǎn)達(dá)到空間位置上的一致，這種一致是指人體上的同一解剖點(diǎn)在兩幅圖像中有相同的空間位置[7]。例如:

S(T)=S(U(X)，V(TX))(1)

其中:S是相似性測(cè)度，T是空間變換矩陣，U(X)、V(X)是要配準(zhǔn)的兩幅圖像。配準(zhǔn)的過程可以歸結(jié)為尋求最佳空間變換T，從而使S達(dá)到最佳。

T^=argmaxT S(T)(2)

空間變換矩陣包含多個(gè)參數(shù)，因此配準(zhǔn)是一個(gè)多參數(shù)優(yōu)化問題。由于GPU內(nèi)部數(shù)據(jù)傳輸?shù)膸掃h(yuǎn)快于GPU與CPU之間的帶寬，在算法設(shè)計(jì)上，大量圖像數(shù)據(jù)的交換主要通過GPU內(nèi)部總線完成;而在CPU和GPU之間傳輸?shù)闹皇亲儞Q參數(shù)和相似性測(cè)度。算法設(shè)計(jì)架構(gòu)如圖1所示，該架構(gòu)對(duì)圖像配準(zhǔn)的并行實(shí)行具有普遍意義。

1.2 空間變換、插值和相似性測(cè)度的并行計(jì)算

圖像空間變換是圖像配準(zhǔn)的基礎(chǔ)，剛性變換保持了圖像中任意兩點(diǎn)間的距離不變。配準(zhǔn)的時(shí)候，本文方法的三維數(shù)據(jù)場(chǎng)的空間變換是以其質(zhì)心為原點(diǎn)的。圖像空間變換、插值和相似性測(cè)度所用數(shù)據(jù)量且各像素完全獨(dú)立，因此可以并行進(jìn)行。

GPU中開辟的線程個(gè)數(shù)等于圖像的大小，線程和X、Y平面中的像素是一一對(duì)應(yīng)的，把線程索引作為平面像素的坐標(biāo)，用變換矩陣T進(jìn)行變換處理，就可以完成對(duì)每個(gè)像素的坐標(biāo)轉(zhuǎn)換了，由于對(duì)三維圖像的像素處理都是并行進(jìn)行的，能夠大大加快配準(zhǔn)的速度。

浮動(dòng)圖像空間坐標(biāo)點(diǎn)通過空間變換后，需要通過插值方法獲得變換點(diǎn)的灰度值，本文采用三線性插值方法。數(shù)據(jù)插值也是在GPU中進(jìn)行的，在完成像素坐標(biāo)的空間變換后，利用GPU內(nèi)建的插值機(jī)制進(jìn)行像素的插值運(yùn)算。

本文配準(zhǔn)實(shí)驗(yàn)中采用的相似性測(cè)度為SSD(sum-of-squares difference)。對(duì)于相同成像模式下的圖像來(lái)說(shuō)，對(duì)應(yīng)位置上像素的灰度值應(yīng)該是相近的，因此，當(dāng)兩幅同模態(tài)的圖像達(dá)到空間位置的一致時(shí)，它們對(duì)應(yīng)像素的灰度差值將取得最小值。設(shè)圖像F(i， j)和圖像R(i， j)分別表示浮動(dòng)圖像和參考圖像，其優(yōu)化目標(biāo)函數(shù)為

ET=∑N-1i=0∑N-1j=0(R(i，j)-FT′(i，j))2(3)

其中，FT′(i， j)是空間變換后得到的圖像，這樣使目標(biāo)函數(shù)E取得最小值的參數(shù)就是最佳的配準(zhǔn)參數(shù)。

SSD的計(jì)算也是并行的，在圖像插值之后，緊接著進(jìn)行參考圖像與浮動(dòng)圖像間對(duì)應(yīng)像素的灰度平方差的運(yùn)算，把得到的SSD傳到CPU上，供參數(shù)優(yōu)化時(shí)使用。

1.3 空間變換參數(shù)的優(yōu)化

剛性配準(zhǔn)需要優(yōu)化空間的三個(gè)平移和旋轉(zhuǎn)參數(shù)(xaix，yaix，zais，θx，θy，θz)，本文采用改進(jìn)Powell方法進(jìn)行多參數(shù)優(yōu)化[8，9]，尋找最佳空間變換參數(shù)。具體流程如下:

設(shè)置精度ε和初始參數(shù)X0，以及N個(gè)線性無(wú)關(guān)的方向向量Uk(k=1，2，…，N)。

a)設(shè)置P0=X0;

b)for k=1，2，…，N，尋找使ET(Pk-1+αk-1Uk)最小的αk-1，并設(shè)置Pk=Pk-1+αk-1Uk;

c)記錄使函數(shù)下降最快的方向Umax=Uk，并記錄此時(shí)的差值max;設(shè)置UN=PN-P0尋找使ET(P0+γUN)最小的γ;

d)X1=P0+γUN，若‖X1-X0‖≤ε， 則停止計(jì)算，輸出X1為解;

e)求出ET(X1) 與ET(X0)的差值，并計(jì)算cha/max的開方s;

f)如果s<γ，設(shè)置U j=Uj(j=1，2，…，k-1)，Ui=Ui+1(i=k，…，N-1)。

g)重復(fù)步驟a)~f)，直到滿足精度ε輸出。

在成像過程中，目標(biāo)平移要比旋轉(zhuǎn)更加明顯，所以三維優(yōu)化的參數(shù)搜索順序是(xaix，yaix，θz，θx，θy，zaix)。參數(shù)的優(yōu)化過程并不能并行進(jìn)行，所以參數(shù)的優(yōu)化過程是在CPU上進(jìn)行的。

1.4 基于CUDA的算法實(shí)現(xiàn)

CUDA中每個(gè)塊中的線程數(shù)目是有一定限制的，最大不超過512，其能發(fā)揮最佳性能的線程數(shù)是64～256個(gè)，結(jié)合醫(yī)學(xué)圖像的大小，為每個(gè)block分配256個(gè)線程。

設(shè)圖像大小為(m，n，p)，使用并行的三維圖像的配準(zhǔn)，按以下步驟進(jìn)行操作:

a)將三維數(shù)據(jù)場(chǎng)讀入主機(jī)內(nèi)存并使用函數(shù)cudaMemcpy()拷貝到設(shè)備內(nèi)存中;

b)在GPU中開辟grid(m/16，n/16)和block(16，16)數(shù)量的并行線程;

c)將CPU上生成變換矩陣T傳入GPU中進(jìn)行坐標(biāo)變換;

d)在GPU上按分配好的線程進(jìn)行圖像的變換操作，計(jì)算式(3)相似性測(cè)度，并傳回CPU供優(yōu)化使用。

GPU上每個(gè)線程的執(zhí)行流程如圖2所示。首先根據(jù)其索引確定要處理的像素;然后進(jìn)行像素坐標(biāo)的空間變換，并根據(jù)插值方式確定變換后坐標(biāo)的像素灰度值;再計(jì)算對(duì)應(yīng)像素之間的灰度差值;最后輸出SSD以供CPU優(yōu)化使用，并記錄變換矩陣生成的圖像。

2 實(shí)驗(yàn)結(jié)果及討論

實(shí)驗(yàn)編程環(huán)境為Visual Studio C++.NET 2005，采用實(shí)驗(yàn)的硬件配置如表1所示。采用三維CT數(shù)據(jù)(512×512×52 16bit)作為參考圖像R;通過對(duì)其進(jìn)行人工平移、旋轉(zhuǎn)得到浮動(dòng)圖像F。本文共進(jìn)行七組配準(zhǔn)實(shí)驗(yàn)，各組實(shí)驗(yàn)的變換參數(shù)設(shè)置如表2所示，平移單位為像素，旋轉(zhuǎn)單位為度。其中第1~6組只對(duì)單個(gè)參數(shù)進(jìn)行變動(dòng)，第7組對(duì)所有參數(shù)同時(shí)變動(dòng)。

表3給出了編號(hào)1~7實(shí)驗(yàn)的定量配準(zhǔn)結(jié)果和所用的時(shí)間，對(duì)比表2可以看出其配準(zhǔn)得到的空間參數(shù)誤差很小，所需的時(shí)間在10 000 ms左右。圖3給出了第7組實(shí)驗(yàn)配準(zhǔn)前后的切面和三維可視化結(jié)果，左列是配準(zhǔn)前的，右列是配準(zhǔn)后的，視覺效果顯示兩幅圖像已經(jīng)對(duì)準(zhǔn)。本文實(shí)驗(yàn)的初始值都設(shè)置為0，實(shí)際使用中醫(yī)生通過交互設(shè)置更合適的起始點(diǎn)。

表1 實(shí)驗(yàn)的硬件配置

處理器內(nèi)存顯卡類型GPU參數(shù)

Inter CoreTM2 DuoProcessor 2.20GHz1.0 GBGeForce980GTX+

處理器數(shù)目核的數(shù)目顯存最大線程數(shù)

168512 MB512

表2 浮動(dòng)圖像空間變換參數(shù)設(shè)置

編號(hào)123456

xaix=-10yaix=-8

θz=-5θy=-4θx=-3zaix=-30

編號(hào)7-8，-5，-4，-3，-4，-3(xaix，yaix，θz，θx，θy，zaix)

從表2和圖3的結(jié)果可以看出，本文并行配準(zhǔn)算法較好地完成了三維醫(yī)學(xué)圖像的配準(zhǔn)，誤差在許可范圍內(nèi)。在不損失精度的前提下，本文算法所用時(shí)間比主流三維配準(zhǔn)算法FLIRT縮短了一個(gè)數(shù)量級(jí)(10 s左右)，非常適合臨床實(shí)時(shí)或者快速環(huán)境中海量三維醫(yī)學(xué)圖像數(shù)據(jù)的配準(zhǔn)應(yīng)用。

3 結(jié)束語(yǔ)

本文提出了一種基于硬件加速的醫(yī)學(xué)圖配準(zhǔn)新方法，實(shí)現(xiàn)了單模三維醫(yī)學(xué)圖像的自動(dòng)快速剛性配準(zhǔn)。該技術(shù)只適合剛體(如頭部)配準(zhǔn)，而對(duì)于肝臟、肺及心臟等需要進(jìn)行彈性配準(zhǔn)的部位，筆者將進(jìn)一步進(jìn)行GPU加速的彈性配準(zhǔn)及多模醫(yī)學(xué)圖像配準(zhǔn)技術(shù)的研究。

參考文獻(xiàn):

[1]FISCHER B，MODERSITZKI J.FLIRT:a flexible image registration toolbox[C]//Proc of WBIR(Workshop Biomedical Image Registration).2003:261-270.

[2]曾德賢， 趙繼廣，曾朝陽(yáng).基于簡(jiǎn)化仿射變換模型的圖像配準(zhǔn)方法 [J]. 裝備指揮技術(shù)學(xué)院學(xué)報(bào)， 2005， 16(1):84-87.

[3]何立新，高曉暉，高琳.醫(yī)學(xué)圖像配準(zhǔn)的新方法 [J]. 四川工業(yè)學(xué)院學(xué)報(bào)， 2004， 26(增刊) :34-36.

[4]SAMANT S S，XIA Jun-yi，OWENS P J D.High performance computing for deformable image registration:towards a new paradigm in adaptive radiotherapy[J].Medical Physics， 2008， 35(8):3546-3553.

[5]LI Bo，YOUNG A A，COWAN B R.GPU accelerated non-rigid registration for the evaluation of cardiac function[C]//Proc of MICCAI.2008:880-887.

[6]張加萬(wàn)，楊甲東，孫濟(jì)洲.基于GPU的彈性圖像配準(zhǔn)方法[J].天津大學(xué)學(xué)報(bào)，2008，41(8):946-950.

[7]NVIDIA Corporation.NVIDIA CUDA computational unified device architecture programming guide version 2.0[K].2008.

[8]MATHEWS J H，F(xiàn)INK K K.Numerical methods using MATLAB[M].北京: 電子工業(yè)出版社，2004.

[9]NOCEDAL J，WRIGHT S.Numerical optimization[M].2nd ed.[S.l.]:Springer，2006.

[10]ZITOVA B，F(xiàn)LUSSER J.Image registration methods: a survey[J].Image and Vision Computing，2003， 21(11):977-1000.