蒲云潔，王學(xué)淵

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

0 引言

隨著社會(huì)醫(yī)療技術(shù)的發(fā)展，使用計(jì)算機(jī)影像輔助醫(yī)生進(jìn)行外科手術(shù)已經(jīng)逐漸成為趨勢，而這其中最重要的一環(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ǔ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)圖像中對應(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)圖像對應(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）測度進(jìn)行相似性測度，進(jìn)而精配準(zhǔn)，在參數(shù)優(yōu)化上采用自適應(yīng)協(xié)方差矩陣進(jìn)化算法（CMA-ES），并且在DRR圖像生成和相似性測度這兩個(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ā)出射線以一定的步長穿透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)是對應(yīng)組織的衰減系數(shù)，0和d分別是光線經(jīng)過3D數(shù)據(jù)的穿入和穿出點(diǎn)。然而由于人體的組織不是均勻的，因而μ(x)是不連續(xù)的，所以式(1)可以改寫成：

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

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

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

對于膝關(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 相似性測度

圖像配準(zhǔn)的相似度測度是指通過幾何空間變換如剛體變換之后的兩幅圖像相似程度，常用的相似性測度方法有：互信息[3]（Mutual Information，MI）、絕對誤差和[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]測度已經(jīng)被應(yīng)用于特征匹配[9]，物體檢測[10]和圖像配準(zhǔn)[11]，但是傳統(tǒng)的梯度方向測度可能會(huì)在具有低梯度幅度的區(qū)域中受到噪聲的影響，從而影響配準(zhǔn)的成功率與精度。本文提出一種改進(jìn)方式：預(yù)先對待配準(zhǔn)的兩幅圖像分別設(shè)定兩個(gè)閾值，在計(jì)算梯度方向測度時(shí)，只計(jì)算這兩幅圖像中大于設(shè)定閾值的像素，這兩個(gè)閾值可以設(shè)定為對應(yīng)圖像所有梯度的中值，這樣就消除了50%的包含低梯度幅度的像素，減少像素遍歷數(shù)，從而降低了算法復(fù)雜度和減小了配準(zhǔn)中噪聲的影響。梯度方向測度函數(shù)如式(9)、式(10)所示：

式中，N是圖像總像素，NL是人為規(guī)定的像素總數(shù)下界，以補(bǔ)償像素?cái)?shù)低于NL的損失，?u是梯度算子，I1，I2是待配準(zhǔn)的兩幅圖像，t1，t2是設(shè)定的梯度閾值，等于對應(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ì)就是相似性測度目標(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)簡單，對目標(biāo)函數(shù)要求寬松，適用于解決高維度復(fù)雜的優(yōu)化問題。在實(shí)際操作中，CMA-ES算法需要先初始化參數(shù)，設(shè)置種群的相關(guān)參數(shù)：種群中心點(diǎn)（期望）、步長、進(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代的步長，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é)方差矩陣并且更新步長。估計(jì)第代協(xié)方差矩陣的公式如下所示：

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

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

式中，e為自然底數(shù)，cδ表示步長學(xué)習(xí)率，dδ為步長更新的阻尼系數(shù)。E||N(0，I)||表示進(jìn)化路徑歸一化后的期望模長，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é)語

本文實(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)鏈接庫完成配準(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é)果圖來說，也可以看出本文算法對于即使擁有噪聲（如患者3的X射線圖像）或者有其他干擾物的情況下（如患者4）仍然能夠取得很好的配準(zhǔn)結(jié)果，這證明了本文算法具有高配準(zhǔn)率和魯棒性，利用CUDA進(jìn)行加速處理后配準(zhǔn)速度也有了質(zhì)的飛躍，為臨床應(yīng)用提供了時(shí)效性。

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