蒲云潔，王學(xué)淵

（西南科技大學(xué) 信息工程學(xué)院，綿陽(yáng) 621000)

0 引言

隨著社會(huì)醫(yī)療技術(shù)的發(fā)展，使用計(jì)算機(jī)影像輔助醫(yī)生進(jìn)行外科手術(shù)已經(jīng)逐漸成為趨勢(shì)，而這其中最重要的一環(huán)就是圖像的配準(zhǔn)。醫(yī)學(xué)圖像以模態(tài)劃分可以分成多種：CT，MRI，PET，X光，內(nèi)窺鏡圖像等，以圖像維度劃分可以分為二維圖像和三維圖像。二維圖像較三維圖像缺乏很多有用的空間信息，因此將這兩種模態(tài)的圖像進(jìn)行信息配準(zhǔn)融合，從而定位患者手術(shù)過程中的靶區(qū)，顯得尤為重要。本文研究的內(nèi)容就是CT斷層圖像和X光圖像之間的配準(zhǔn)，這是一種多模態(tài)配準(zhǔn)模式，也是一種2D/3D配準(zhǔn)模式。目前常用的2D/3D配準(zhǔn)方法主要分為以下3種：

1）基于圖像灰度的配準(zhǔn)算法[1]，該方法一般是先用DRR算法模擬X光圖像的生成，從患者術(shù)前采集的三維體數(shù)據(jù)中得出多幅DRR圖像來和術(shù)中采集的X光照片進(jìn)行相似性判定配準(zhǔn)，此算法精度高，但耗費(fèi)的時(shí)間較長(zhǎng)，適用于剛體之間的配準(zhǔn)。

2）基于圖像特征的配準(zhǔn)算法[2]，該方法依賴在病人身上提前做好的標(biāo)記點(diǎn)或者人體內(nèi)固有的標(biāo)記點(diǎn)進(jìn)行定位，計(jì)算出相應(yīng)的空間變換，進(jìn)而計(jì)算出待配準(zhǔn)圖像中對(duì)應(yīng)標(biāo)記點(diǎn)的最短特征距離，從而達(dá)到配準(zhǔn)的目的，此方法的配準(zhǔn)速度雖然很快，但是精度與標(biāo)記點(diǎn)的選取有很大的關(guān)聯(lián)，魯棒性也比較差。

3）基于以上兩類方法結(jié)合的配準(zhǔn)算法，即既需要計(jì)算出待配準(zhǔn)圖像對(duì)應(yīng)標(biāo)記點(diǎn)的空間變換和最短特征距離，也需要利用相似性原理進(jìn)行配準(zhǔn)，配準(zhǔn)速度快精度也好，但此方法捕捉范圍較小，魯棒性較差。

本文根據(jù)以上三種配準(zhǔn)方法提出一種基于CUDA加速的多模態(tài)配準(zhǔn)方法，該方法基于圖像灰度，在DRR圖像生成上采用光線投射算法，先進(jìn)行剛體變換粗配準(zhǔn)，再采用梯度方向（Gradient orientation，GO）測(cè)度進(jìn)行相似性測(cè)度，進(jìn)而精配準(zhǔn)，在參數(shù)優(yōu)化上采用自適應(yīng)協(xié)方差矩陣進(jìn)化算法（CMA-ES），并且在DRR圖像生成和相似性測(cè)度這兩個(gè)耗時(shí)最大的過程采用支持CUDA架構(gòu)的GPU處理器進(jìn)行加速，最大化的減少配準(zhǔn)時(shí)間。

1 DRR圖像生成

DRR（Digitally Reconstructed Radiograph）數(shù)字影像重建技術(shù)是2D/3D醫(yī)學(xué)圖像配準(zhǔn)的關(guān)鍵技術(shù)，也是最耗時(shí)最影響配準(zhǔn)精度的步驟，它最終目的是為了從三維體數(shù)據(jù)中獲取一系列模擬的二維X射線圖像。獲取DRR圖像的算法主要有：拋雪球法（Splatting）、光線投射法（Ray-Casting）、剪切變形法（Shear-Warp）等。本文使用的是光線投射算法。

1.1 光線投射算法

光線投射算法是模擬X射線穿透三維體數(shù)據(jù)元（本文用的數(shù)據(jù)是CT體數(shù)據(jù)），經(jīng)過衰減和吸收后投影到DRR平面，就形成了一張DRR圖像，形成過程如圖1所示。

圖1 光線投影示意圖

光線投射算法生成DRR圖像具體過程如下：當(dāng)虛擬X光源發(fā)出射線以一定的步長(zhǎng)穿透3D體數(shù)據(jù)時(shí)，沿著其行進(jìn)方向采集CT值，再把得到的CT值轉(zhuǎn)化成相應(yīng)的衰減系數(shù)μj，累加此方向上的衰減系數(shù)，從而計(jì)算得到射線到達(dá)DRR平面上的電子強(qiáng)度I，最后形成DRR圖像。

DRR圖像生成過程中，若入射光線電子強(qiáng)度Ii，它與穿過3D體數(shù)據(jù)后到達(dá)DRR平面上的電子強(qiáng)度Io關(guān)系如式(1)所示：

式(1)中μ(x)是對(duì)應(yīng)組織的衰減系數(shù)，0和d分別是光線經(jīng)過3D數(shù)據(jù)的穿入和穿出點(diǎn)。然而由于人體的組織不是均勻的，因而μ(x)是不連續(xù)的，所以式(1)可以改寫成：

式(2)中μj是光線穿過不同組織的衰減系數(shù)，xj是穿過這種組織的有效長(zhǎng)度。而衰減系數(shù)μ與CT值之間又有如下的轉(zhuǎn)換關(guān)系：

式(3)中μw是純水的衰減系數(shù)，通常是定值，HU是組織的CT值。最后將CT值再轉(zhuǎn)化成灰度值就可以完成DRR圖像的生成。

1.2 DRR圖像和X射線圖像之間的剛體變換

對(duì)于膝關(guān)節(jié)圖像，無論是DRR圖像還是X射線圖像，這兩者之間的配準(zhǔn)可以看成是剛體之間的配準(zhǔn)，則兩者之間的變換就是剛體變換。剛體變換進(jìn)一步可以分解為旋轉(zhuǎn)變換和平移變換，從而可以使用六參數(shù)矢量來表示，其中θx，θy，θz分別表示物體圍繞x，y，z軸的旋轉(zhuǎn)變換，tx，ty，tz表示物體沿著x，y，z軸的平移變換。進(jìn)一步將旋轉(zhuǎn)變換和平移變換寫成矩陣的形式：

則剛體變換矩陣為：

2 相似性測(cè)度

圖像配準(zhǔn)的相似度測(cè)度是指通過幾何空間變換如剛體變換之后的兩幅圖像相似程度，常用的相似性測(cè)度方法有：互信息[3]（Mutual Information，MI）、絕對(duì)誤差和[4]（Sum of Absolute Differences，SAD）、歸一化互相關(guān)[5]（Normalized Cross Correlation，NCC），梯度相關(guān)[6]（Gradient Correlation，GC）、模式強(qiáng)度[7]（Pattern Intensity，PI）等。

梯度方向[8]測(cè)度已經(jīng)被應(yīng)用于特征匹配[9]，物體檢測(cè)[10]和圖像配準(zhǔn)[11]，但是傳統(tǒng)的梯度方向測(cè)度可能會(huì)在具有低梯度幅度的區(qū)域中受到噪聲的影響，從而影響配準(zhǔn)的成功率與精度。本文提出一種改進(jìn)方式：預(yù)先對(duì)待配準(zhǔn)的兩幅圖像分別設(shè)定兩個(gè)閾值，在計(jì)算梯度方向測(cè)度時(shí)，只計(jì)算這兩幅圖像中大于設(shè)定閾值的像素，這兩個(gè)閾值可以設(shè)定為對(duì)應(yīng)圖像所有梯度的中值，這樣就消除了50%的包含低梯度幅度的像素，減少像素遍歷數(shù)，從而降低了算法復(fù)雜度和減小了配準(zhǔn)中噪聲的影響。梯度方向測(cè)度函數(shù)如式(9)、式(10)所示：

式中，N是圖像總像素，NL是人為規(guī)定的像素總數(shù)下界，以補(bǔ)償像素?cái)?shù)低于NL的損失，?u是梯度算子，I1，I2是待配準(zhǔn)的兩幅圖像，t1，t2是設(shè)定的梯度閾值，等于對(duì)應(yīng)圖像梯度的中值，θi是?uI1(i)和?uI2(i)之間的夾角。最后計(jì)算的GO_Value值越接近1配準(zhǔn)效果越好。

3 參數(shù)優(yōu)化

參數(shù)優(yōu)化是醫(yī)學(xué)圖像配準(zhǔn)的重要過程，其實(shí)質(zhì)就是相似性測(cè)度目標(biāo)函數(shù)進(jìn)行不斷迭代優(yōu)化來尋找最佳的變換參數(shù)，以達(dá)到配準(zhǔn)成功的目的。常用的優(yōu)化算法[3]有模擬退火算法，Powell算法，梯度下降算法，遺傳算法等。本文采用的是自適應(yīng)協(xié)方差矩陣進(jìn)化算法（CMA-ES）。

CMA-ES算法是Nikolaus Hansen等人[12]在2006年提出的一種模擬生物進(jìn)化過程的算法。該算法的特點(diǎn)是假設(shè)不論基因發(fā)生何種變化，產(chǎn)生的結(jié)果（性狀）總遵循這零均值，某一方差的高斯分布。CMA-ES算法實(shí)現(xiàn)簡(jiǎn)單，對(duì)目標(biāo)函數(shù)要求寬松，適用于解決高維度復(fù)雜的優(yōu)化問題。在實(shí)際操作中，CMA-ES算法需要先初始化參數(shù)，設(shè)置種群的相關(guān)參數(shù)：種群中心點(diǎn)（期望）、步長(zhǎng)、進(jìn)化路徑、子代數(shù)量、協(xié)方差矩陣、最大迭代次數(shù)以及一些自適應(yīng)參數(shù)。算法步驟[13]如下：

1）種群采樣操作：是指在目標(biāo)函數(shù)的解空間中隨機(jī)選擇一個(gè)解，并以此為中心生產(chǎn)正態(tài)分布種群。以式(10)進(jìn)行采樣并產(chǎn)生新解。

式中k=1，2，...λ，～表示式子兩邊服從相同分布，表示第k+1代的第n個(gè)后代，m(k)表示第k代搜索分布的均值，δ(k)表示第k代的步長(zhǎng)，N(0，C(k))表示均值為0，協(xié)方差矩陣為C(k)的正態(tài)分布。第k代的協(xié)方差矩陣為C(k)，C(0)=In*n，I是單位矩陣。協(xié)方差矩陣是CMA-ES算法進(jìn)化策略的關(guān)鍵。

2）選擇和重組操作：是指在生成的總?cè)褐羞x擇最優(yōu)子群。當(dāng)前新的最優(yōu)子群通過加權(quán)重組得到新的期望：

3）更新操作：從一代代的種群中估計(jì)協(xié)方差矩陣并且更新步長(zhǎng)。估計(jì)第代協(xié)方差矩陣的公式如下所示：

其中，c1≤1是協(xié)方差矩陣進(jìn)化路徑Pc的學(xué)習(xí)率，T表示矩陣轉(zhuǎn)置。hδ表示Heaviside函數(shù)，它的作用是控制進(jìn)化路徑Pc的過大增長(zhǎng)。通常用平滑指數(shù)來構(gòu)造進(jìn)化路徑：

步長(zhǎng)δ更新函數(shù)：

式中，e為自然底數(shù)，cδ表示步長(zhǎng)學(xué)習(xí)率，dδ為步長(zhǎng)更新的阻尼系數(shù)。E||N(0，I)||表示進(jìn)化路徑歸一化后的期望模長(zhǎng)，I為單位矩陣。

CMA-ES算法通過以上三個(gè)步驟的循環(huán)迭代，直到找到全局最優(yōu)解。

4 基于CUDA加速的多模態(tài)圖像配準(zhǔn)

傳統(tǒng)的圖像配準(zhǔn)方法通常是基于軟件層面的，耗時(shí)非常大。因此，采取基于支持CUDA架構(gòu)的GPU處理器進(jìn)行硬件加速，彌補(bǔ)傳統(tǒng)配準(zhǔn)方法在耗時(shí)上面的不足，大大提高配準(zhǔn)的速度。本文的配準(zhǔn)算法流程如圖2所示。

圖2 本文算法流程圖

5 結(jié)語(yǔ)

本文實(shí)驗(yàn)數(shù)據(jù)來源為佛山市中醫(yī)院提供的4名患者膝關(guān)節(jié)CT斷層體數(shù)據(jù)，單層圖像大小為512×512，配套有同一患者的X射線圖像。硬件配置為Intel(R) Core(TM) i7-8700K CPU，主機(jī)頻率為3.7GHz，顯卡是NVIDIA公司的GeForce GTX1060 6GB，軟件配置為Win10下的Visual Studio 2015+Python 3.6+CUDA 9.0。實(shí)驗(yàn)以裁剪后的感興趣區(qū)域X射線圖像為參考圖像，CT體數(shù)據(jù)生成的DRR圖像為浮動(dòng)圖像，在Python環(huán)境下調(diào)用VS2015和CUDA生成的動(dòng)態(tài)鏈接庫(kù)完成配準(zhǔn)，結(jié)果如圖3～圖6所示：

患者1：

圖3 患者1配準(zhǔn)圖

患者2：

圖4 患者2配準(zhǔn)圖

患者3：

圖5 患者3配準(zhǔn)圖

患者4：

圖6 患者4配準(zhǔn)圖

表1數(shù)據(jù)是每例患者數(shù)據(jù)均配準(zhǔn)50次之后取的平均值，從4組數(shù)據(jù)的配準(zhǔn)結(jié)果圖來說，也可以看出本文算法對(duì)于即使擁有噪聲（如患者3的X射線圖像）或者有其他干擾物的情況下（如患者4）仍然能夠取得很好的配準(zhǔn)結(jié)果，這證明了本文算法具有高配準(zhǔn)率和魯棒性，利用CUDA進(jìn)行加速處理后配準(zhǔn)速度也有了質(zhì)的飛躍，為臨床應(yīng)用提供了時(shí)效性。

表1 GPU與CPU配準(zhǔn)對(duì)比