王文賽，邢恩銘，王 韜，周 榴，王曉春，周 盛，宋學(xué)東，楊 軍

（中國(guó)醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院生物醫(yī)學(xué)工程研究所，天津 300192）

0 引言

高頻超聲三維成像可以直接提供淺表器官的立體結(jié)構(gòu)特征，顯示病灶與周圍組織的空間位置關(guān)系，從而提高臨床診斷的準(zhǔn)確性[1]。三維超聲成像分為三維數(shù)據(jù)采集、三維重建和三維圖像可視化3個(gè)階段，三維數(shù)據(jù)采集的質(zhì)量直接影響三維重建圖像的效果[2]。在三維數(shù)據(jù)采集過(guò)程中，三維超聲探頭自身的振動(dòng)、操作者和受試者不自主的運(yùn)動(dòng)均會(huì)造成圖像的抖動(dòng)，給圖像采集帶來(lái)誤差，嚴(yán)重影響三維重建的效果。本文的高頻超聲皮膚三維成像系統(tǒng)中的換能器頻率為50 MHz，其焦點(diǎn)處波束寬度約為50 μm，成像精度較高，但也易受到外部干擾的影響，一般的運(yùn)動(dòng)會(huì)帶來(lái)比普通B超更嚴(yán)重的后果。因此，在采集到高頻超聲三維數(shù)據(jù)后，首先要進(jìn)行圖像配準(zhǔn)，盡可能降低三維數(shù)據(jù)采集過(guò)程中帶來(lái)的誤差，提高三維重建圖像的精度。

圖像配準(zhǔn)是指對(duì)圖像序列中的2幅圖像，通過(guò)某種空間變換把一幅圖像映射到另一幅圖像，使得2幅圖像中同一空間位置的點(diǎn)相互對(duì)應(yīng)，即實(shí)現(xiàn)2幅圖像在空間和灰度上的映射。醫(yī)學(xué)圖像有多種配準(zhǔn)方法，大致可分為基于灰度的配準(zhǔn)方法、基于特征的配準(zhǔn)方法和基于變換域的配準(zhǔn)方法三大類?；诨叶鹊呐錅?zhǔn)方法是利用圖像灰度值的統(tǒng)計(jì)信息來(lái)度量圖像的相似程度，通常不需要對(duì)圖像進(jìn)行復(fù)雜的預(yù)處理，實(shí)現(xiàn)較為簡(jiǎn)單，應(yīng)用較廣；基于特征的配準(zhǔn)方法首先要根據(jù)待配準(zhǔn)圖像的特點(diǎn)進(jìn)行特征提取，利用提取到的特征進(jìn)行圖像配準(zhǔn)；基于變換域的配準(zhǔn)方法中最常用的是傅里葉變換域配準(zhǔn)方法，圖像的平移、旋轉(zhuǎn)等一系列變換在變換域中都有相應(yīng)的體現(xiàn)[3-5]。

基于圖像的灰度進(jìn)行配準(zhǔn)的方法能避免特征提取帶來(lái)的誤差，也無(wú)需復(fù)雜的預(yù)處理，是目前醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域廣泛應(yīng)用的配準(zhǔn)方法[6]。本文為了提高高頻超聲皮膚三維重建圖像的精度，以C++語(yǔ)言為基礎(chǔ)，結(jié)合圖像處理工具包ITK（insight segmentation and registration toolkit）中基于灰度的配準(zhǔn)方法對(duì)高頻超聲三維數(shù)據(jù)進(jìn)行配準(zhǔn)，并從配準(zhǔn)速度、配準(zhǔn)精度方面嘗試對(duì)配準(zhǔn)算法進(jìn)行評(píng)價(jià)，找出適用于高頻超聲皮膚三維成像的配準(zhǔn)算法。

1 ITK配準(zhǔn)原理

ITK源于美國(guó)的可視人計(jì)劃，是一個(gè)開源的、面向?qū)ο蟮?、跨平臺(tái)的醫(yī)學(xué)影像分析工具[7]。其封裝了大量基于C++實(shí)現(xiàn)的圖像分割和圖像配準(zhǔn)以及圖像統(tǒng)計(jì)的算法，并用CMake軟件管理編譯過(guò)程以保證編譯過(guò)程獨(dú)立于平臺(tái)。ITK開源和算法復(fù)用性強(qiáng)的特點(diǎn)，使得程序開發(fā)工作簡(jiǎn)化，因此在醫(yī)學(xué)圖像處理領(lǐng)域廣泛應(yīng)用[8]。

1.1 ITK圖像配準(zhǔn)框架

ITK圖像配準(zhǔn)的框架如圖1所示，主要由2幅輸入圖像、相似性度量、優(yōu)化器、校對(duì)機(jī)、變換函數(shù)組成。圖像配準(zhǔn)過(guò)程就是尋找最優(yōu)的變換，使得變換后的待配準(zhǔn)圖像和參考圖像之間的相似度達(dá)到最大[7]。

圖1 ITK圖像配準(zhǔn)的框架

1.2 ITK配準(zhǔn)流程

基于ITK進(jìn)行配準(zhǔn)時(shí)，首先，準(zhǔn)備好待輸入的2幅圖像，即參考圖像和待配準(zhǔn)圖像；其次，利用變換函數(shù)將參考圖像進(jìn)行坐標(biāo)變換得到新的區(qū)域，利用校對(duì)機(jī)得到待配準(zhǔn)圖像在上一步新區(qū)域內(nèi)的坐標(biāo)；然后，使用相似性度量計(jì)算待配準(zhǔn)圖像和參考圖像之間的相似度，相似性度量得到的結(jié)果輸入到優(yōu)化器中，優(yōu)化器通過(guò)迭代進(jìn)行最優(yōu)化計(jì)算，得到最終的變換參數(shù)；最后，輸出在最優(yōu)變換下得到的圖像，即配準(zhǔn)后圖像。

配準(zhǔn)過(guò)程就是一個(gè)尋求最優(yōu)解的過(guò)程，配準(zhǔn)過(guò)程中每迭代一次，得到一個(gè)相似性度量值，將此參數(shù)和設(shè)定的參數(shù)進(jìn)行比對(duì)，如果達(dá)到預(yù)期效果則停止迭代；若設(shè)定的參數(shù)選取不合適，迭代就會(huì)無(wú)休止進(jìn)行。一般進(jìn)行配準(zhǔn)時(shí)需要設(shè)置迭代次數(shù)，以便觀察配準(zhǔn)結(jié)果。

1.3 ITK配準(zhǔn)中的變換函數(shù)

ITK圖像配準(zhǔn)支持二維圖像配準(zhǔn)、三維圖像配準(zhǔn)以及單模態(tài)圖像配準(zhǔn)、多模態(tài)圖像配準(zhǔn)。本文在二維單模態(tài)圖像配準(zhǔn)情況下，以均方值為相似性度量函數(shù)，保持其他配準(zhǔn)組成部分不變，僅改變變換函數(shù)，選取平移變換函數(shù)、剛性變換函數(shù)、二維相似變換函數(shù)、仿射變換函數(shù)以及B樣條可變形變換函數(shù)對(duì)高頻超聲皮膚圖像進(jìn)行配準(zhǔn)實(shí)驗(yàn)。

（1）平移變換函數(shù)是最簡(jiǎn)單、有效的變換函數(shù)，表示輸入空間中點(diǎn)的簡(jiǎn)單變換，雖然僅有水平和垂直方向上的位移，但平移變換速度較快，能減少配準(zhǔn)時(shí)間。在開始進(jìn)行配準(zhǔn)實(shí)驗(yàn)時(shí)，平移變換函數(shù)是最佳變換函數(shù)。平移變換可用矩陣表示為

其中，（x0，y0）為待配準(zhǔn)圖像上的點(diǎn)，（x1，y1）為變換后坐標(biāo)，tx為水平位移量，ty為垂直位移量。

（2）剛性變換函數(shù)可以分解為圖像在一個(gè)二維平移變換后，圍繞一個(gè)固定的中心點(diǎn)進(jìn)行旋轉(zhuǎn)，這個(gè)轉(zhuǎn)動(dòng)中心一般選擇圖像的中心位置。圖像中任意兩點(diǎn)的距離在剛性變換配準(zhǔn)前后保持不變。剛性變換可用矩陣表示為

其中，θ為旋轉(zhuǎn)角度，（Cx，Cy）是轉(zhuǎn)動(dòng)中心的坐標(biāo)。

（3）二維相似變換函數(shù)可以認(rèn)為是剛性變換和各向同性的縮放系數(shù)的結(jié)合，表示圍繞一個(gè)固定的中心點(diǎn)對(duì)二維圖像進(jìn)行均勻縮放、旋轉(zhuǎn)和平移。二維相似性變換可用矩陣表示為

其中，λ是縮放因子。

（4）仿射變換函數(shù)是圖像配準(zhǔn)中應(yīng)用最普遍的變換，圖像中的直線在仿射變換配準(zhǔn)之后仍為直線，能保持二維圖像的平直性。仿射變換表現(xiàn)為在各個(gè)方向尺度變換系數(shù)一致的均勻尺度變換或變換系數(shù)不一致的非均勻尺度變換和剪切變換[9-10]。仿射變換可用矩陣表示為

可以通過(guò)調(diào)節(jié)參數(shù)m00～m12實(shí)現(xiàn)圖像的平移、縮放、剪切、旋轉(zhuǎn)等操作。

（5）B樣條可變形變換函數(shù)可以解決變形配準(zhǔn)的問(wèn)題。在圖像采集過(guò)程中，由于操作者和受試者不自主的生理運(yùn)動(dòng)或者微小的移動(dòng)會(huì)造成皮膚等淺表器官或其組織內(nèi)部局部發(fā)生形變，這時(shí)剛性配準(zhǔn)不能滿足配準(zhǔn)要求[11]。為了兼顧B樣條函數(shù)配準(zhǔn)的速度和精度，本文使用三階B樣條可變形變換函數(shù)對(duì)采集到的高頻超聲皮膚三維數(shù)據(jù)進(jìn)行配準(zhǔn)實(shí)驗(yàn)。

1.4 配準(zhǔn)效果評(píng)價(jià)

本文選用相關(guān)系數(shù)（correlation coefficient，CC）、結(jié)構(gòu)相似性（structural similarity，SSIM）及余弦相似度（cosine similarity，Cos）算法嘗試對(duì)5種變換函數(shù)進(jìn)行定量評(píng)價(jià)[12-13]。

CC算法通過(guò)配準(zhǔn)后圖像M和參考圖像F之間的線性關(guān)系的強(qiáng)度和方向來(lái)評(píng)價(jià)圖像的相似度，公式如下：

其中，Mi表示配準(zhǔn)后圖像M的第i個(gè)像素值，F(xiàn)i表示參考圖像F的第i個(gè)像素值，N是圖像的像素點(diǎn)數(shù)目，μM是配準(zhǔn)后圖像M的灰度平均值，μF是參考圖像F的灰度平均值。CC取值范圍為[0，1]，其值越接近1，表明配準(zhǔn)后圖像M和參考圖像F的相似度越高，配準(zhǔn)效果越好。

SSIM算法通過(guò)比較參考圖像與配準(zhǔn)后圖像的對(duì)比度、亮度、結(jié)構(gòu)3個(gè)方面的相似程度來(lái)評(píng)價(jià)圖像的相似度，公式如下：

其中，σ2表示圖像的方差，σMF表示配準(zhǔn)后圖像M與參考圖像F的協(xié)方差，C1與C2為常數(shù)。SSIM取值范圍為[0，1]，其值接近于1，表明配準(zhǔn)后圖像M和參考圖像F的相似度越高，配準(zhǔn)效果越好。

Cos算法通過(guò)將圖片表示成向量，計(jì)算向量間的余弦距離來(lái)表示配準(zhǔn)后圖像M與參考圖像F的相似度，公式如下：

其中，Mm和Fn分別代表由配準(zhǔn)后圖像M和參考圖像F構(gòu)成的向量。Cos取值范圍為[-1，1]，其值越接近1，表明配準(zhǔn)后圖像M和參考圖像F的相似度越高，配準(zhǔn)效果越好。

本文中高頻超聲皮膚圖像配準(zhǔn)精度需要結(jié)合以上評(píng)價(jià)指標(biāo)進(jìn)行綜合評(píng)價(jià)。

2 配準(zhǔn)實(shí)驗(yàn)與結(jié)果分析

本文使用高頻超聲三維數(shù)據(jù)采集系統(tǒng)采集了10組皮膚三維超聲數(shù)據(jù)，每一組由標(biāo)準(zhǔn)三維數(shù)據(jù)和增加干擾后的三維數(shù)據(jù)組成，每一個(gè)三維數(shù)據(jù)由100幀每幀間隔100 μm的皮膚二維超聲圖像組成，三維數(shù)據(jù)的體素大小為768×576×100[14]，三維數(shù)據(jù)對(duì)應(yīng)皮膚的尺寸為 15 mm×11 mm×10 mm（長(zhǎng)×寬×高）。在標(biāo)準(zhǔn)三維數(shù)據(jù)采集過(guò)程中盡可能減少三維探頭的振動(dòng)，在同一位置采集4次進(jìn)行疊加平均，以減少隨機(jī)干擾，獲取較好的標(biāo)準(zhǔn)三維數(shù)據(jù)；標(biāo)準(zhǔn)三維數(shù)據(jù)采集完成后三維探頭回到數(shù)據(jù)采集初始位置，再通過(guò)添加抖動(dòng)三維探頭干擾，采集干擾后的三維數(shù)據(jù)。

本文采集圖像所用的超聲換能器為中心頻率50 MHz的單陣元超聲換能器，高頻超聲三維數(shù)據(jù)采集系統(tǒng)主要由該換能器組成的機(jī)械定位線性三維掃描探頭構(gòu)成。如圖2所示，在進(jìn)行數(shù)據(jù)采集時(shí)將待測(cè)部位放入水中或者在待測(cè)部位增加水囊，超聲換能器不接觸待測(cè)部位，可以較好地獲得待測(cè)部位的原始形態(tài)。在工作時(shí)將三維探頭固定在機(jī)械臂上，這樣可以最大程度地降低人為影響。三維數(shù)據(jù)采集的主要干擾來(lái)源于三維探頭的振動(dòng)，該振動(dòng)產(chǎn)生的皮膚位移會(huì)給圖像質(zhì)量帶來(lái)影響，本文配準(zhǔn)實(shí)驗(yàn)的目的就是盡可能消除這種影響，增加的抖動(dòng)干擾是為了模擬實(shí)際應(yīng)用中三維探頭可能出現(xiàn)的干擾因素。

本文配準(zhǔn)實(shí)驗(yàn)以C++語(yǔ)言為基礎(chǔ)，結(jié)合ITK圖像處理工具包對(duì)采集到的高頻超聲三維數(shù)據(jù)進(jìn)行配準(zhǔn)，然后利用可視化工具包VTK（visualization toolkit）對(duì)配準(zhǔn)前后的三維數(shù)據(jù)進(jìn)行三維重建顯示。

圖2 皮膚超聲三維數(shù)據(jù)采集圖

2.1 皮膚圖像配準(zhǔn)實(shí)驗(yàn)

本文以均方值為相似性度量函數(shù)，分別使用平移變換、剛性變換、二維相似變換、仿射變換以及B樣條可變形變換5種變換算法對(duì)采集到的皮膚二維圖像序列進(jìn)行配準(zhǔn)，并比較皮膚超聲圖像中的配準(zhǔn)效果。利用VTK對(duì)配準(zhǔn)前后的高頻超聲三維數(shù)據(jù)進(jìn)行三維重建，顯示皮膚的三維立體圖像，圖3顯示了其中一組的配準(zhǔn)效果。

由圖3可以看出，高頻超聲三維數(shù)據(jù)經(jīng)過(guò)5種變換算法配準(zhǔn)后，其重建結(jié)果有相應(yīng)變化，皮膚表面的抖動(dòng)有一定程度的消除，但僅從三維圖像的表面情況不能準(zhǔn)確判斷配準(zhǔn)變換函數(shù)的優(yōu)劣?；诖?，本文在皮膚三維圖像的相同位置選取切面，通過(guò)觀察三維圖像的內(nèi)部情況來(lái)比較變換函數(shù)的優(yōu)劣，如圖4所示。

由圖4可以看出，經(jīng)過(guò)平移變換、剛性變換以及仿射變換的配準(zhǔn)后圖像質(zhì)量比配準(zhǔn)前有明顯改善，而二維相似變換和B樣條可變形變換對(duì)圖像質(zhì)量的改善并不理想。

2.2 皮膚圖像配準(zhǔn)效果評(píng)價(jià)

為了對(duì)變換函數(shù)進(jìn)行定量評(píng)價(jià)，在使用每種變換函數(shù)配準(zhǔn)后，統(tǒng)計(jì)CC、SSIM、Cos及配準(zhǔn)時(shí)間。10組三維數(shù)據(jù)配準(zhǔn)時(shí)間的平均值以及相似性度量的平均值取值見(jiàn)表1。

從表1中的3種相似性度量（CC、SSIM和Cos）的取值可以看出，與配準(zhǔn)前比較，幾種變換函數(shù)配準(zhǔn)方法都使圖像質(zhì)量有了改善，且剛性變換函數(shù)的配準(zhǔn)效果最好，平移變換以及仿射變換的配準(zhǔn)效果相似，二維相似變換和B樣條可變形變換配準(zhǔn)效果較差。由此可見(jiàn)，雖然皮膚是軟組織，在數(shù)據(jù)采集過(guò)程中可能存在變形，但由于皮膚受呼吸等不自主生理運(yùn)動(dòng)的影響較小，可近似看作剛體。綜上，在高頻超聲皮膚三維重建過(guò)程中，三維超聲數(shù)據(jù)配準(zhǔn)最理想的變換函數(shù)是剛性變換。

圖3 高頻超聲皮膚圖像配準(zhǔn)前和5種變換算法配準(zhǔn)后的三維重建圖

雖然剛性變換函數(shù)的配準(zhǔn)效果最好，但是其配準(zhǔn)皮膚超聲三維數(shù)據(jù)的時(shí)間較長(zhǎng)，不利于高頻超聲皮膚三維成像系統(tǒng)的實(shí)時(shí)性。從表1中可以看出，平移變換的效果和剛性變換很接近，但是平移變換的時(shí)間還不到剛性變換的十分之一。考慮到高頻超聲皮膚三維成像系統(tǒng)的實(shí)時(shí)性，在臨床使用中選用平移變換函數(shù)進(jìn)行圖像配準(zhǔn)是較好的選擇。

3 討論

本文基于C++語(yǔ)言與ITK醫(yī)學(xué)圖像處理工具包在配準(zhǔn)框架的其他組成部分一致的情況下，通過(guò)改變變換函數(shù)進(jìn)行圖像配準(zhǔn)實(shí)驗(yàn)。實(shí)驗(yàn)分別選取了平移變換函數(shù)、剛性變換函數(shù)、二維相似變換函數(shù)、仿射變換函數(shù)以及B樣條可變形變換函數(shù)對(duì)高頻超聲皮膚三維數(shù)據(jù)進(jìn)行配準(zhǔn)，結(jié)果發(fā)現(xiàn)人體皮膚在數(shù)據(jù)采集過(guò)程中可以近似看作剛體，5種變換函數(shù)中剛性變換的配準(zhǔn)效果最佳，但是考慮到其配準(zhǔn)時(shí)間較長(zhǎng)，為了兼顧配準(zhǔn)速度與精度，建議在高頻超聲皮膚三維成像系統(tǒng)中使用與剛性變換配準(zhǔn)效果接近且配準(zhǔn)時(shí)間較短的平移變換。

圖4 高頻超聲皮膚圖像配準(zhǔn)前和5種變換算法配準(zhǔn)后三維重建的切面圖

表1 變換函數(shù)配準(zhǔn)效果比較

本文中基于圖像的灰度值進(jìn)行配準(zhǔn)的方法能避免特征提取帶來(lái)的誤差，且無(wú)需復(fù)雜的圖像預(yù)處理。在皮膚超聲三維數(shù)據(jù)的配準(zhǔn)中剛性變換函數(shù)的配準(zhǔn)效果較好，但是所需時(shí)間較長(zhǎng)，為了平衡時(shí)間與精度要求建議選用平移變換函數(shù)。未來(lái)研究?jī)?nèi)容主要是配準(zhǔn)算法的加速研究，其中GPU（graphics processing unit）加速可以有效提升配準(zhǔn)的運(yùn)行速度，大大節(jié)省配準(zhǔn)時(shí)間，達(dá)到配準(zhǔn)方法的實(shí)時(shí)性要求，更好地應(yīng)用于高頻超聲皮膚三維成像系統(tǒng)，提高臨床診斷的準(zhǔn)確率。