基于特征點(diǎn)互信息的醫(yī)學(xué)圖像配準(zhǔn)*

李 敏 劉 俊

（武漢科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 武漢 430065）

1 引言

隨著科技的不斷發(fā)展，醫(yī)學(xué)圖像在臨床醫(yī)學(xué)上有著越來(lái)越多的應(yīng)用，一張圖像無(wú)法為醫(yī)生提供完整全面的信息，然而將多張醫(yī)學(xué)圖像進(jìn)行融合后，將不同的成像信息體現(xiàn)在一幅圖像上，這樣不僅可以為醫(yī)生節(jié)約時(shí)間，更可以為醫(yī)生提供更加全面的信息。在圖像融合之前，圖像配準(zhǔn)就是其中最為關(guān)鍵的一步?？梢哉f(shuō)醫(yī)學(xué)圖像配準(zhǔn)［2］技術(shù)就是圖像融合的基礎(chǔ)和關(guān)鍵。

在20 世紀(jì)80 年代，就出現(xiàn)了醫(yī)學(xué)圖像配準(zhǔn)課題的研究。此課題一經(jīng)提出就備受關(guān)注，逐漸變成了圖像處理技術(shù)的一個(gè)熱門領(lǐng)域。由于其關(guān)注度，學(xué)者們提出了許許多多的圖像配準(zhǔn)算法，通?？梢詫D像配準(zhǔn)方法分為兩類：基于圖像特征信息的配準(zhǔn)算法［3］和基于灰度信息的配準(zhǔn)算法?；趫D像特征信息的配準(zhǔn)算法是對(duì)圖像中具有穩(wěn)定性和重復(fù)性的圖像特征進(jìn)進(jìn)行分析和匹配?；诨叶刃畔⒌呐錅?zhǔn)算法［4］是充分利用圖像的灰度信息來(lái)對(duì)圖像進(jìn)行配準(zhǔn)。

每一種圖像配準(zhǔn)方法［13］都受到一些因素的影響，其中基于互信息的配準(zhǔn)方法是應(yīng)用最為廣泛的，本文根據(jù)互信息在醫(yī)學(xué)圖像配準(zhǔn)中的優(yōu)缺點(diǎn)，采用ReNet的方法提取圖像邊緣，構(gòu)建邊緣特征點(diǎn)互信息能量函數(shù)，通過(guò)改進(jìn)Powell算法對(duì)配準(zhǔn)進(jìn)行優(yōu)化。這樣不僅可以提高配準(zhǔn)的準(zhǔn)確率，減少配準(zhǔn)時(shí)間，還可以減少噪聲對(duì)配準(zhǔn)的影響。

2 特征點(diǎn)互信息配準(zhǔn)方法

2.1 圖像邊緣提取

2.1.1 圖像預(yù)處理

首先輸入圖像，使用K-SVD 方法得到冗余字典D，并通過(guò)在原有每個(gè)元素中減去元素的平均值來(lái)獲取FoE 的自適應(yīng)濾波器，之后通過(guò)使用奇異值分解（SVD）學(xué)習(xí)與當(dāng)前圖像最相關(guān)的FoE 自適應(yīng)濾波器。其中K-SVD的表達(dá)式如下：

其中，Y表示樣本集，D表示冗余字典，ε表示重構(gòu)誤差所允許的最大值。

奇異值分解（Singular Value Decomposition，SVD）是線性代數(shù)中一種重要的矩陣分析，通過(guò)SVD對(duì)數(shù)據(jù)進(jìn)行處理，可以用小的多的數(shù)據(jù)集來(lái)表示原有的數(shù)據(jù)集。就是將實(shí)際原有數(shù)據(jù)中的噪聲和冗余信息去除，以達(dá)到優(yōu)化數(shù)據(jù)、提高結(jié)果的目的。SVD的分解公式表達(dá)如下：

由于U對(duì)前一個(gè)原子的投影值較大，對(duì)后一個(gè)原子的投影值較小，所以采用U后面一個(gè)原子作為FoE 的自適應(yīng)濾波器。將基于正則化的框架應(yīng)用到恢復(fù)圖像之中來(lái)恢復(fù)圖像，表達(dá)式如下：

其中，Ri表示樣本矩陣；JK-SVD=[]J1,…,Jk表示每一列Ji對(duì)應(yīng)的自適應(yīng)濾波器。

2.1.2 空間遞歸層

空間遞歸層用來(lái)接收?qǐng)D像大小為H×W 的經(jīng)過(guò)增強(qiáng)的圖像或上一層的特征圖像，其中H和W 分別表示圖像的高度和寬度。然后將圖像分割為大小為Wp×hp的非重疊的圖像集。準(zhǔn)確地說(shuō)，使用ReNet層接收2D映射作為輸入，沿相反方向掃描網(wǎng)絡(luò)，掃描方式如下所示：

在彼此之上，通過(guò)疊加多層組來(lái)構(gòu)建一個(gè)簡(jiǎn)單的深層的ReNet，其中第一層組用來(lái)接收原始圖像像素信息作為輸入，而最后一個(gè)層組通過(guò)Softmax層來(lái)產(chǎn)生密集的預(yù)測(cè)，但是其產(chǎn)生的預(yù)測(cè)結(jié)果精度較低。Softmax的表達(dá)式如下：

其中，Vi表示分類器前級(jí)輸出單位的輸出；i 表示類別索引；C表示總的類別個(gè)數(shù)；Si表示當(dāng)前元素的指數(shù)與所有元素指數(shù)的比值。

2.1.3 空間遞歸層

首先我們采用VGG-16 網(wǎng)絡(luò)［10］作為全卷積網(wǎng)絡(luò)的基線，我們僅使用了VGG-16 網(wǎng)絡(luò)的前幾層。將上一節(jié)討論為ReNet 層［6］插入到所選的VGG-16網(wǎng)絡(luò)的全卷積網(wǎng)絡(luò)的頂層，然后將該層接著一個(gè)卷積層，并通過(guò)轉(zhuǎn)置卷積或分?jǐn)?shù)步長(zhǎng)卷積進(jìn)行上采樣。

由于從多個(gè)層次組合特征會(huì)導(dǎo)致更具區(qū)分性特征，因此我們將VGG-16 的不同選擇層的特征映射連接起來(lái)，在VGG-16 中系統(tǒng)地從每個(gè)層選擇特征。由于最終的組合可能由一個(gè)或兩個(gè)特征映射所支配，并且特征映射的大小也會(huì)有很大的波動(dòng)，因此我們?cè)谶B接之前進(jìn)一步對(duì)特征映射進(jìn)行規(guī)范化?？偟倪吘壧卣鹘Y(jié)構(gòu)提取的結(jié)構(gòu)體系如圖1 所示。

圖1 結(jié)構(gòu)體系圖

2.2 特征點(diǎn)互信息熵

傳統(tǒng)基于互信息的醫(yī)學(xué)圖像配準(zhǔn)方法中，互信息熵的計(jì)算量相對(duì)較大。為了在盡量不減少圖像間特征信息的前提下，減少互信息熵［16］得到計(jì)算量，提高配準(zhǔn)效率，可以將相同特征點(diǎn)集的互信息熵代替整張圖像的互信息作為配準(zhǔn)的目標(biāo)函數(shù)。假設(shè)從參考圖像F 和浮動(dòng)圖像M 中分別提取出的特征點(diǎn)集合表示為

其中，Pij表示圖像特征點(diǎn)集的聯(lián)合概率，也就是說(shuō)同時(shí)從X中選取Xi和從Y中選取Yi的概率。

2.3 配準(zhǔn)優(yōu)化策略

為了解決互信息在配準(zhǔn)過(guò)程中計(jì)算量較大，減少配準(zhǔn)所用時(shí)間的問(wèn)題，本節(jié)選用一種改進(jìn)的Powell 算法作為配準(zhǔn)測(cè)度函數(shù)的優(yōu)化策略，其優(yōu)化算法的執(zhí)行過(guò)程如下。

首先計(jì)算配準(zhǔn)圖像初始位置V0 的互信息值F0。由初始位置向左和向右平移一個(gè)像素點(diǎn)，分別計(jì)算其互信息值，選取互信息增加的一側(cè)作為尋優(yōu)方向不斷搜索，直至互信息值開始減小。若互信息值開始減少，則跳轉(zhuǎn)回前一個(gè)點(diǎn)，在垂直方向按相同方法開始所需極值點(diǎn)。水平和垂直方向搜尋出最優(yōu)點(diǎn)后，對(duì)旋轉(zhuǎn)方向求取極值點(diǎn)。計(jì)算［-180，180］范圍內(nèi)每個(gè)單元處的互信息值，選取其中的最大值，作為旋轉(zhuǎn)方向的極值點(diǎn)，將最終得到的點(diǎn)作為下一個(gè)Powell算法的初始值。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)

3.1.1 實(shí)驗(yàn)環(huán)境

本組實(shí)驗(yàn)使用的硬件配準(zhǔn)為I7 8700K處理器，GPU 型號(hào)為NVIDA GeForce GTX 1080TI，顯存11G，操作系統(tǒng)使用的是Windows 7，軟件使用的是Matlab 2016a。

3.1.2 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)分別使用了三組數(shù)據(jù)集，第一組為美國(guó)一實(shí)驗(yàn)室的血管超聲圖像，第二組為brain web 的腦部醫(yī)學(xué)圖像，第三組為臨床使用的肺部CT圖像。

3.2 邊緣提取實(shí)驗(yàn)

圖像邊緣時(shí)指其周圍像素灰度有階躍變化或屋頂狀變化的像素集合。常見的邊緣提取方法有Canny 算子、Sobel算子、Prewitt 算子、整體嵌套邊緣檢測(cè)（Holistically-Nested Edge Detection，HED）、CASENet等。本文使用肺部CT 圖像對(duì)邊緣提取方法進(jìn)行對(duì)比實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如圖2所示。

圖2 邊緣提取實(shí)驗(yàn)

Sobel 算子對(duì)灰度漸變?cè)肼暤膱D像檢測(cè)效果較好，但是對(duì)復(fù)雜噪聲的圖像處理效果不理想；Canny算子［11］相對(duì)于其他傳統(tǒng)邊緣檢測(cè)算子效果最佳，但是其對(duì)噪聲敏感，容易將噪聲誤判為邊界；Prewitt算子對(duì)噪聲較多、灰度漸變的圖像效果較好，但是其存在檢測(cè)出來(lái)的邊緣較粗，人為選取閾值極易造成邊緣點(diǎn)誤判的缺點(diǎn)；HED 和CASENet 都是基于深度學(xué)習(xí)的圖像邊緣輪廓提取方法，在往常的圖像邊緣提取中都取得較好的效果，但是由于醫(yī)學(xué)圖像中存在這數(shù)據(jù)據(jù)數(shù)量較少的緣故，導(dǎo)致最終的邊緣檢測(cè)效果不佳；而本節(jié)中利用特定的卷積神經(jīng)網(wǎng)絡(luò)，對(duì)圖像進(jìn)行邊緣提取，從圖2 中可以看出，基于ReNet的邊緣提取方法效果最佳。

為了更加客觀地評(píng)價(jià)［15］邊緣檢測(cè)提取的效果，我們這里用到定位誤差來(lái)對(duì)其進(jìn)行檢測(cè)。對(duì)整張圖象進(jìn)行掃描，當(dāng)遇到邊緣點(diǎn)，判斷檢測(cè)圖像中是否有對(duì)應(yīng)邊緣點(diǎn)。最終計(jì)算出整張圖像的定位誤差，定位誤差表達(dá)式如下所示：

其中，Pmis表示邊緣圖像漏檢率，Perr表示誤檢率，Na表示圖像的總像素點(diǎn)數(shù)，Nmis表示圖像中沒(méi)有被標(biāo)記的漏檢點(diǎn)個(gè)數(shù)，Nerr表示誤檢的邊緣點(diǎn)個(gè)數(shù)。

表1 客觀地表示出邊緣檢測(cè)提取的效果評(píng)測(cè)。從漏檢率上看Canny 算子最低；從誤檢率上看Prewitt算子與本文中所使用的方法最小，可以說(shuō)明本文所使用方法能較為準(zhǔn)確地判斷邊緣點(diǎn)；從定位誤差上來(lái)看本文中所使用的基于ReNet 的邊緣提取方法的定位誤差最小，可以說(shuō)明本文所使用的方法效果較好。

表1 邊緣檢測(cè)誤差判斷指標(biāo)

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

3.3.1 噪聲干擾實(shí)驗(yàn)

PC-SIFT（phase congruency based on scale invitation feature transform）是一種將相位一致性代替原先SIFT中［7］的梯度用于形象結(jié)構(gòu)描述的方法，它在特征點(diǎn)提取上有著較好的表現(xiàn)。首先比較PC-SIFT算法［9］與本文所用算法在特征提取方面的魯棒性，在特征提取之前，先對(duì)圖像進(jìn)行模擬噪聲，在對(duì)比在不同噪聲［8］程度下兩種方法的提取結(jié)果，之后再對(duì)配準(zhǔn)結(jié)果進(jìn)行分析。

圖3 為將圖像中分別加入0%、1%、3%、5%和7%的噪聲后特征點(diǎn)提取的個(gè)數(shù)，是對(duì)20 組加入噪聲的圖像進(jìn)行特征提取實(shí)驗(yàn)后，統(tǒng)計(jì)其特征點(diǎn)的平均值。從圖像中可以看出PC-SIF T 方法［14］雖然在一些情況下特征點(diǎn)提取的數(shù)量較多，但是其特征點(diǎn)提取受噪聲影響波動(dòng)較大；而本文中所使用的方法，在加入不同噪聲的情況下，特征點(diǎn)提取的數(shù)量大致在1200 左右?？偟貋?lái)說(shuō)本文所使用的方法在受到噪聲影響的情況下特征點(diǎn)數(shù)量波動(dòng)較小，所以本文所使用方法在特征提取方面魯棒性更強(qiáng)。

圖3 噪聲干擾實(shí)驗(yàn)

3.3.2 準(zhǔn)確率實(shí)驗(yàn)

為了驗(yàn)證本文使用的醫(yī)學(xué)圖像配準(zhǔn)方法的準(zhǔn)確性，使用brain web 的腦部醫(yī)學(xué)圖像進(jìn)行實(shí)驗(yàn)，分別將本文使用的方法與基于傳統(tǒng)互信息配準(zhǔn)方法、SIFT 算法、Active Demons 算法［1］進(jìn)行對(duì)比試驗(yàn)，實(shí)驗(yàn)結(jié)果如圖4所示。

圖4 腦部配準(zhǔn)對(duì)比試驗(yàn)

從圖4 中可以直觀看出，基于傳統(tǒng)互信息配準(zhǔn)方法、SIFT 算法、Active Demons 算法［5］和本文使用的方法都可以較好地對(duì)腦部圖像進(jìn)行配準(zhǔn)。表2分別使用互信息、相關(guān)系數(shù)和算法運(yùn)行時(shí)間對(duì)腦部圖像配準(zhǔn)結(jié)果進(jìn)行定量分析。從表中可以看出，本文所使用的方法在MI、CC的值都相對(duì)較大，且運(yùn)行時(shí)間上本文使用方法略低于其他方法?？偟貋?lái)說(shuō)，本文使用的方法在提升配準(zhǔn)精度的前提下，縮短了配準(zhǔn)所需的時(shí)間。

表2 腦部圖像配準(zhǔn)對(duì)比實(shí)驗(yàn)

3.3.2 非剛性配準(zhǔn)

使用超聲血管圖像進(jìn)行配準(zhǔn)與其他圖像不同，針對(duì)超聲血管圖像進(jìn)行配準(zhǔn)僅需要對(duì)病灶區(qū)域進(jìn)行配準(zhǔn)即可，若使用整張圖像進(jìn)行配準(zhǔn)不僅會(huì)降低配準(zhǔn)的準(zhǔn)確率還會(huì)增加配準(zhǔn)所需的時(shí)間?，F(xiàn)將基于傳統(tǒng)互信息配準(zhǔn)方法、Active Demons 算法［12］、分割后使用結(jié)構(gòu)互信息配準(zhǔn)方法與本文所使用方法使用在超聲血管圖像配準(zhǔn)上，配準(zhǔn)結(jié)果如圖5 所示。

圖5 血管配準(zhǔn)對(duì)比試驗(yàn)

從圖5 可以看出，由于傳統(tǒng)互信息配準(zhǔn)方法和Active Demons 算法都是針對(duì)整張圖像進(jìn)行配準(zhǔn)，所以導(dǎo)致整張圖像發(fā)生了形變，造成最終配準(zhǔn)效果較差。而使用對(duì)病灶區(qū)域進(jìn)行分割后在結(jié)構(gòu)互信息配準(zhǔn)的方法與本文使用的方法從直觀的角度看都可以較好地對(duì)血管圖像進(jìn)行配準(zhǔn)。表3 對(duì)這兩種方法進(jìn)行定量分析，從互信息和相關(guān)系數(shù)上看，分割后結(jié)構(gòu)互信息的配準(zhǔn)準(zhǔn)確率略高于本文使用的方法，但從配準(zhǔn)所需時(shí)間上看，本文所使用方法的配準(zhǔn)效率更高。

表3 超聲血管圖像配準(zhǔn)評(píng)價(jià)指標(biāo)

4 結(jié)語(yǔ)

醫(yī)學(xué)影像在臨床醫(yī)學(xué)中有著越來(lái)越重要的作用，結(jié)合互信息在醫(yī)學(xué)圖像配準(zhǔn)中的優(yōu)缺點(diǎn)，并且結(jié)合相同醫(yī)學(xué)圖像集之間的共享配準(zhǔn)模式的特點(diǎn)，提出了基于特征點(diǎn)互信息的醫(yī)學(xué)圖像配準(zhǔn)方法。通過(guò)驗(yàn)證，該方法在一定程度上提高了配準(zhǔn)的準(zhǔn)確率，減少了配準(zhǔn)所需的時(shí)間，并且減少了圖像噪聲對(duì)配準(zhǔn)的影響。