一種改進(jìn)的病理顯微圖像亞像素快速配準(zhǔn)方法

，

(華南理工大學(xué) 機(jī)械與汽車工程學(xué)院，廣州 510641)

0 引言

圖像配準(zhǔn)是指求解兩幅或者多幅具有相同場景或者內(nèi)容的圖像之間幾何變換關(guān)系，圖像配準(zhǔn)是圖像處理、機(jī)器視覺以及醫(yī)療成像中最重要的步驟之一[1]，醫(yī)學(xué)圖像拼接融合、計算機(jī)視覺領(lǐng)域的目標(biāo)定位以及衛(wèi)星遙感圖像[2]等應(yīng)用均對圖像配準(zhǔn)的精度有較高的要求。

掃描獲取一張數(shù)字病理切片[3]通常需要進(jìn)行數(shù)以千計的病理顯微圖像配準(zhǔn)拼接，拼接過程中的圖像對的配準(zhǔn)是耗費時間的主要部分，由于病理切片對于醫(yī)療診斷有著極其重要的意義[4]，數(shù)字病理切片需要盡可能準(zhǔn)確還原病理切片的所有信息，因而病理顯微圖像的快速高精度配準(zhǔn)是獲取高質(zhì)量數(shù)字切片圖像的最為關(guān)鍵一步[5]。

常見的醫(yī)學(xué)圖像的配準(zhǔn)算法主要有以下幾種：1)Appleton采用的互信息[6]等圖像灰度區(qū)域信息配準(zhǔn)方法，精度較高但配準(zhǔn)速度慢；2)基于特征點的配準(zhǔn)方法，彭勛所采用的改進(jìn)SIFT算法[7]能夠有效配準(zhǔn)具有旋轉(zhuǎn)、畸變等情況的圖像，但特征點的篩選和匹配需要耗費大量計算，且在病理顯微圖像配準(zhǔn)過程中出現(xiàn)的空白稀疏圖像會導(dǎo)致特征點過少而配準(zhǔn)失??；3)基于相位相關(guān)法[8]，通過將圖像位移轉(zhuǎn)換為圖像頻域的相位進(jìn)行求解，相位相關(guān)法對光照變化不敏感，但是配準(zhǔn)效率不高。

針對病理顯微圖像配準(zhǔn)的應(yīng)用特點，本文提出了一種新的快速配準(zhǔn)方法。該方法利用對數(shù)差值函數(shù)對待匹配圖像進(jìn)行信息評估，并構(gòu)建模板獲取像素級粗定位，根據(jù)粗定位獲取待配準(zhǔn)子圖，由相位相關(guān)法獲取亞像素級細(xì)定位。

1 亞像素配準(zhǔn)技術(shù)

1.1 歸一化相位相關(guān)法

相位相關(guān)法主要是基于傅里葉變換中的平移定理，假設(shè)f1(x,y)和f2(x,y)是兩幅存在平移變換(xs,ys)關(guān)系的圖像，滿足：

f2(x,y)=f1(x-xs,y-ys)

(1)

它們對應(yīng)的傅里葉變換分別為F1(u，v)和F2(u，v)，則有：

F2(u,v)=e-j2π(uxs+vys)F1(u,v)

(2)

則f1(x,y)和f2(x,y)的歸一化互功率譜為：

(3)

根據(jù)變換平移理論,互相關(guān)功率譜的相位等于兩個圖像的相位差,通過求互相關(guān)功率譜的傅里葉逆變換可以得到相位相關(guān)函數(shù)：

p(x,y)=F-1(ej2π(uxs+vys)=δ(x-xs,y-ys)

(4)

式中，δ(x-xs,y-ys)為典型的Dirac函數(shù)，也稱為沖激函數(shù)，該函數(shù)在中心點(xs,ys)處不為零，在其他位置均為零。其坐標(biāo)位置(xs,ys)即為圖像平移量,由此得到了圖像配準(zhǔn)關(guān)系[9]。

1.2 基于相位獲取亞像素配準(zhǔn)技術(shù)

經(jīng)典的相位相關(guān)法只能獲取像素級精度的位移參數(shù)，為了進(jìn)一步獲取亞像素級精度的位移參數(shù)，常見的方法有以下幾種：1)基于擬合的方法，利用正弦函數(shù)來逼近互功率譜的傅里葉逆變換所得的Dirac函數(shù)，再對擬合的譜函數(shù)進(jìn)行上采樣，從而獲取亞像素的配準(zhǔn)估計值[10]；2)基于相位差解析的方法，兩幅存在平移變換的圖像之間的相位差是一個2D鋸齒函數(shù)，根據(jù)其角分量在兩個軸上重復(fù)的周期數(shù)與平移參數(shù)所存在的解析關(guān)系即可獲得亞像素平移參數(shù)；3)Smith提出對圖像傅里葉變換后進(jìn)行零填充(zero-padding)上采樣[11]，再進(jìn)行相位相關(guān)法求取亞像素配準(zhǔn)值，但其計算資源的耗費巨大往往難以滿足實時處理的需求。在此基礎(chǔ)上，Soummer[12]等提出了基于矩陣乘法的傅里葉變換，該算法在像素級粗匹配值的一定鄰域范圍內(nèi)進(jìn)行k倍上采樣，計算出精度為1/k的亞像素配準(zhǔn)值。

2 改進(jìn)的病理顯微圖像亞像素配準(zhǔn)算法

基于矩陣乘法的上采樣離散傅里葉變換雖然一定程度上避免了傳統(tǒng)的零填充上采樣方法在求取圖像亞像素位移時的不足，但是在利用相位相關(guān)法求取像素級位移的過程中，其計算量會隨著圖像尺寸的增大而急劇變大，這對于尺寸較大的病理顯微圖像而言，求取像素級位移的計算量及算法所耗時間將非常大。并且，由于相位相關(guān)法只提取了兩幅圖像互功率譜中的相位信息，雖然減少了對圖像本身內(nèi)容的依賴，但同時又缺乏了對于病理顯微圖像配準(zhǔn)拼接過程中對于信息量的評估，對于病理切片掃描過程中常出現(xiàn)空白的顯微圖像不能進(jìn)行有效地規(guī)避，浪費了計算資源以及時效還可能導(dǎo)致錯誤的配準(zhǔn)結(jié)果。為了充分利用局部上采樣相位相關(guān)法在亞像素配準(zhǔn)的有效性，又能使其更適用于病理顯微圖像的配準(zhǔn)應(yīng)用，本文提出采用一種由對數(shù)差值函數(shù)來評價顯微圖像信息量，由對數(shù)差值函數(shù)構(gòu)建模板匹配，從而獲取待配準(zhǔn)子圖完成相位相關(guān)亞像素配準(zhǔn)，提高配準(zhǔn)的速度。

2.1 對數(shù)差值函數(shù)

用i表示圖像的像素灰度值，則其灰度值對數(shù)差值函數(shù)[13]為：

fls(I(x,y),I(x',y'))=Kln(I(x,y)/I(x',y'))=

Kln(I(x',y'))-Kln(I(x,y))=u[I(x',y')]-u[I(x,y)]

(5)

若i表示像素灰度值，則式(5)中：

u[i]=Kln(i)

(6)

對于常用的8 bit灰度圖像，i∈[0,255]，ln(i)的值一般比較小，所以K為大于1的常數(shù)放大系數(shù)。當(dāng)i=0時，ln(i)無效，為了處理方便，用一個極小常量ε(0<ε<<1)代替，即u[0]=Klnε。

對數(shù)差值函數(shù)有具備以下特性：

1)|fls(I(x,y),I(x',y'))|=|fls(I(x',y'),I(x,y))|,即表示對數(shù)差值函數(shù)值只與像素之間的灰度變化相關(guān)，與兩像素點的先后次序無關(guān)；

2)當(dāng)I(x',y')>I(x,y)時，fls(I(x,y),I(x',y'))>0，即表示兩像素間灰度值遞增；當(dāng)I(x',y')

3)I(x',y')與I(x,y)差異越小，|fls(I(x,y),I(x',y'))|越?。籌(x',y')與I(x,y)差異越大，|fls(I(x,y),I(x',y'))|越大；即|fls(I(x,y),I(x',y'))|可以反映了兩個像素之間灰度變化幅度的大?。?)u[i]的值可以預(yù)先計算并存在表格里，使用fls(I(x,y),I(x',y'))評價像素間的灰度變化幅度時只需要進(jìn)行查表和減法運算，計算量很小。

2.2 病理顯微圖像亞像素配準(zhǔn)算法的流程

第一步，構(gòu)建H型對數(shù)差值模板。

以圖像的左上角為原點O(0,0)，以水平向右為X軸正方向，垂直向下為Y軸正方向，H型對數(shù)模板由如圖所示U、H、V三條線段組成，線段U和V的長度為l，線段H的長度為s。以線段U的第一個點P(x,y)確定模板的位置，線段H的第一個點坐標(biāo)為P(x,y+m),則線段U上的像素依次為P(x,y+k),線段U上的像素依次為P(x+s,y+k),線段H上像素依次為P(x+n,y+m),k=0,1,…,l-1,n=0,1,…,s-1，下文中將線段U,V,H的各像素點的灰度值統(tǒng)一分別記為Ik,Jk,Dn,Ik,Jk,Dn∈[0,255],n=0,1,…,s-1,k=0,1,…l-1。

圖1 H型對數(shù)差值模板

使用對數(shù)差值函數(shù)分別計算線段U和線段H上相鄰兩個像素間的灰度值變化,定義特征向量：

αlsp=(a0,a1,…,al-1)T=

(fls(I0,I1),fls(I1,I2),…fls(Il-2,Il-1))

(7)

γlsp= (y0,y1,…,yl-1)T=

(fls(D0,D1),fls(D1,D2),…fls(Dn-2,Dn-1))

(8)

計算獲取線段U和V上對應(yīng)位置像素間的灰度值變化，定義特征向量：

βlsp= (b0,b1,…,bl-1)T=

(fls(I0,I1),fls(I1,I2), …fls(Il-2,Il-1))

(9)

αlsp、βlsp和γlsp，共同組成H型對數(shù)差值模板Tlsp即：

Tlsp=(αlsp,βlsp,γlsp)

(10)

在病理顯微圖像配準(zhǔn)中，圖像A,B為有一定重疊區(qū)域的待配準(zhǔn)圖像對，在圖像A重疊區(qū)域的搜索區(qū)域中，線段U為使|αlsp|獲得最大值的線段，線段V為使|βlsp|取得最大值的線段，U和V之間的距離即為s。此時便構(gòu)建了H型對數(shù)差值模板Tlsp(A)。線段U和V確定之后，線段H為平行線段U、V之間|γlsp|取得最大值的線段。

構(gòu)建模板Tlsp(A)時要先確定l的參數(shù)，H型對數(shù)差值模板長度l應(yīng)當(dāng)小于待搜索區(qū)域，增加模板長度l可以有效提高配準(zhǔn)精度，同時也會增加計算量，降低配準(zhǔn)速度。因此，長度l應(yīng)權(quán)衡速度與精度取適當(dāng)?shù)闹?。參?shù)s,l為構(gòu)建最佳匹配模板Tlsp(B(x,y))的必備參數(shù)信息。

第二步，圖像信息量評估。

在切片數(shù)字化掃描過程中，由于高倍率的物鏡成像和病理切片的特點，采集的顯微圖像時常會包含非常稀疏甚至空白的區(qū)域(如圖2例子所示)，在此區(qū)域內(nèi)進(jìn)行配準(zhǔn)會導(dǎo)致誤匹配，所以在構(gòu)建模板Tlsp(A)的過程中，通過遍歷搜索|αlsp|、|βlsp|以及|γlsp|各自的最大值，同時也是對待配準(zhǔn)區(qū)域的灰度變化信息量進(jìn)行了評估。

其中|αlsp|反映了在待配準(zhǔn)區(qū)域圖像中長度為l垂直線段相鄰像素灰度變化最大的幅度，|βlsp|反映了在待配準(zhǔn)區(qū)域圖像中長度為l平行線段對應(yīng)位置像素灰度變化最大的幅度，|γlsp|反映了模板平行線段U、V范圍內(nèi)水平線段相鄰灰度值變化的最大幅度。

對于顯微圖像重疊區(qū)域為空白或者稀疏區(qū)域時，|αlsp|、|βlsp|以及|γlsp|的值均較小，分別設(shè)置灰度變化幅度閾值THα,THβ,THγ來判定模板是否擁有足夠豐富的信息量來滿足匹配要求，當(dāng)搜索所得|αlsp|、|βlsp|、|γlsp|有：

(11)

(12)

(13)

即表示在圖像A的搜索區(qū)域范圍內(nèi)整體像素的灰度變化幅度小，可認(rèn)定為是圖像搜索區(qū)域范圍內(nèi)像素灰度值單調(diào)甚至為空白區(qū)域，無法提供足夠的灰度變化信息構(gòu)建模板Tlsp(A)完成配準(zhǔn)，即便構(gòu)建模板或者進(jìn)行相位相關(guān)配準(zhǔn)其配準(zhǔn)結(jié)果亦不具備可靠性，應(yīng)當(dāng)返回信息量不足的錯誤配準(zhǔn)信息。

第三步，進(jìn)行對數(shù)差值模板匹配獲取粗定位。

在圖像A中創(chuàng)建模板Tlsp(A)之后,并經(jīng)過信息量評估判定，便可在圖像B的待搜索區(qū)域中搜索Tlsp(A)的最佳配準(zhǔn)模板Tlsp(B)以獲取像素級別的粗配準(zhǔn)值。配準(zhǔn)示意圖如圖3所示。

圖3 匹配示意圖

在顯微圖像B的待搜索區(qū)域內(nèi)各點PB(x,y)依據(jù)模板Tlsp(A)的尺寸參數(shù)s和l構(gòu)建H型對數(shù)差值模板Tlsp(B(x,y))，計算Tlsp(A)和Tlsp(B(x,y))的差異值。

|Tlsp(B(x,y))-Tlsp(A)|=μ|αlsp(B(x,y))-αlsp(A) |+

ω|βlsp(B(x,y))-βlsp(A) |+(1-μ)|

γlsp(B(x,y))-γlsp(A) |

(14)

其中，μ(0<μ<1)、ω(0<ω<1)為各特征向量差異值的權(quán)重值。

為了進(jìn)一步減少計算量，設(shè)置閾值TOv過濾灰度變化幅度較小的匹配點,αlsp(A)中分量絕對值大于TOv的位置的集合為Cα,βlsp(A)中分量絕對值大于TOv的位置的集合為Cβ,則有：

|αlsp(B(x,y))-αlsp(A)|=∑i∈Cα|ai(B(x,y))-ai(A)|

(15)

|βlsp(B(x,y))-βlsp(A)|=∑j∈Cβ|bj(B(x,y))-bj(A)|

(16)

通過設(shè)置閾值TOv大小可以減少匹配點的數(shù)量，降低比較特征向量差異所需要進(jìn)行的減法運算次數(shù)，從而加快搜索匹配模板的速度，但是TOv過大也會造成匹配點的減少，致使誤配率升高，所以TOv的大小要在速度和精度考慮中取適中的值。

當(dāng)Tlsp(A)和Tlsp(B(x,y))的差異值取得最小值時，該位置的Tlsp(B(x,y))就是Tlsp(A)的最佳匹配。從而可以獲取顯微圖像A與B的像素級圖像偏移值Δpixel=(xpixel,ypixel)。

第四步，獲取待配準(zhǔn)子圖。

在獲取了Tlsp(A)的最佳匹配Tlsp(B(x,y))之后，為了進(jìn)一步校正以及獲取亞像素級圖像配準(zhǔn)偏移值，根據(jù)模板定位坐標(biāo)PA(x,y),PB(x,y)分別在顯微圖像A,B截取圖像大小為Na×Na、完全包含模板Tlsp(A)、Tlsp(B(x,y))的子圖a,b，l

經(jīng)過第二步的圖像信息量評估以及第三步的像素級粗定位，實際上是對兩幅顯微圖像待配準(zhǔn)區(qū)域的所有子圖進(jìn)行篩選，通過比值模板匹配獲取粗定位，并截取其小范圍子圖作為亞像素配準(zhǔn)的輸入圖像對，圖像尺寸的縮小可以有效減小相位相關(guān)法計算亞像素偏移值的計算量。

第五步，獲取亞像素級細(xì)定位。

對兩幅待匹配子圖a,b進(jìn)行傅里葉變換并計算其相位相關(guān)歸一化互功率譜,在兩幅大小為M×N的子圖計算所得的互動率譜的中心位置(xc,yc)附近范圍構(gòu)建大小為2c鄰域，在此鄰域區(qū)域內(nèi)采用矩陣相乘計算互功率譜矩陣M(u,v)的K倍上采樣局部互功率譜矩陣，其大小為cK×cK，矩陣形式為：

(17)

其中，X=[0,1…,cK-1]T-cK/2+xcK;Y=[0,1…,cK-1]T-cK/2+ycK;U=[0,1,…M-1]T-M/2;V=[0,1,…N-1]T-N/2；M'(U,V)為互功率譜矩陣M(U,V)的中心變換形式。

在第三步對數(shù)模板匹配正確獲取像素級的偏移值的情況下，兩幅待配準(zhǔn)的相鄰子圖僅存在亞像素級偏移，將中心位置(xc,yc)的1.5×1.5像素大小的鄰域作為細(xì)定位的搜尋區(qū)域，可以確保亞像素級的細(xì)定位峰值點在此范圍之內(nèi)。病理顯微拼接算法的上采樣倍率K=100，即細(xì)定位將在150×150區(qū)域內(nèi)獲取亞像素定位點，因而運動偏移量的配準(zhǔn)估計精度可以達(dá)到0.01像素。通過求解式(17)的傅里葉逆變換即可獲得顯微圖像A與B的亞像素級的圖像偏移值Δsub=(xsub,ysub)。

因此，改進(jìn)的病理顯微圖像亞像素配準(zhǔn)的偏移量為Δ(x,y)為：

(18)

算法流程圖如圖4所示。

圖4 病理顯微圖像配準(zhǔn)算法流程圖

3 實驗結(jié)果與分析

3.1 實驗平臺

為了驗證改進(jìn)的病理顯微圖像配準(zhǔn)方法的可行性、精度以及效率，在如圖5所示的實驗平臺對切片標(biāo)本采集圖像進(jìn)行了實驗。實驗平臺主要由機(jī)械運動平臺、圖像采集系統(tǒng)及軟件平臺組成。PC主要配置為Intel Core i5-4440,4G內(nèi)存，wndows7操作系統(tǒng)，算法測試軟件使用Microsoft Visual C++ 2013開發(fā)。

圖5 實驗平臺

3.2 實驗一

為了驗證改進(jìn)的病理顯微圖像配準(zhǔn)算法的精度，在實驗平臺上對由廣東醫(yī)學(xué)院病理教研室提供的體大息肉ESD標(biāo)本采集圖像尺寸大小為2448×2048的病理顯微圖像樣例，經(jīng)過10倍上采樣，從中截取五組的偏移量不同的待匹配圖像組，經(jīng)由降采樣后得到大小為640×480存在亞像素平移圖像組(圖6展示其中一組)。分別使用文獻(xiàn)[12]中所提出的交互相關(guān)法和本文所提出的配準(zhǔn)方法來對這五組相鄰的病理顯微圖像進(jìn)行配準(zhǔn)。得出的配準(zhǔn)平移估計值如表1所示。

圖6 精度測試待配準(zhǔn)顯微圖像組

由表1數(shù)據(jù)可知，相較于文獻(xiàn)[12]方法，本文方法在計算圖像的亞像素位移是對圖像組子圖進(jìn)行矩陣乘法相位相關(guān)法所獲取，對五組相鄰圖像配準(zhǔn)的最大偏差為0.04，配準(zhǔn)精度達(dá)到0.01像素，驗證了該方法具備亞像素級配準(zhǔn)精度，能夠很好地滿足病理顯微圖像高精度配準(zhǔn)的需求。

3.3 實驗二

為了驗證改進(jìn)的病理顯微圖像配準(zhǔn)算法在實際配準(zhǔn)應(yīng)用的配準(zhǔn)速度以及可行性，在實驗平臺上分別對體大息肉ESD標(biāo)本采集兩種不同圖像尺寸相鄰顯微圖像共80對，40對顯微圖像尺寸大小為800×600，水平方向上重疊區(qū)域大小在[520,544](像素)范圍內(nèi)波動；另外40對顯微圖像尺寸大小為2448×2048，水平方向上重疊區(qū)域大小在[1636,1650](像素)范圍內(nèi)波動。分別采用文獻(xiàn)[12]算法與本文方法對80對顯微圖像進(jìn)行配準(zhǔn)，并將兩種算法80次水平配準(zhǔn)耗時進(jìn)行對比，以此來比較兩種算法的配準(zhǔn)速度。兩種算法分別進(jìn)行80對病理顯微圖像的配準(zhǔn)耗時如圖7所示。圖8為其中一組相鄰顯微圖像以及它們的進(jìn)行相位相關(guān)計算的兩幅子圖。

圖7 兩種方法配準(zhǔn)80對顯微圖像耗時結(jié)果

從圖7數(shù)據(jù)可以得知，對于尺寸大小為800×600的40對顯微圖像，文獻(xiàn)[12]的方法配準(zhǔn)平均耗時230.125 ms，本文方法配準(zhǔn)平均耗時89.225 ms，速度提升1.5倍；對于尺寸大小為2448×2048的40對顯微圖像，文獻(xiàn)[12]的方法配準(zhǔn)平均耗時988.4 ms，本文方法配準(zhǔn)平均耗時261.45 ms，速度提升2.7倍；本文算法速度相較于文獻(xiàn)[12]方法有著明顯提高?？梢姡S著圖像尺寸的增大，文獻(xiàn)[12]算法配準(zhǔn)耗時急劇增大，而本文方法得益于模板匹配與子圖的應(yīng)用，配準(zhǔn)效率有著明顯的提升，這對于病理切片掃描需要經(jīng)過數(shù)千次圖像配準(zhǔn)的應(yīng)用需求有著重要的意義。選取的子圖表明經(jīng)過信息評估，H型模板避開了顯微圖像的空白區(qū)域提高了配準(zhǔn)可靠性。

并將這80對顯微圖像按照計算所得的配準(zhǔn)值進(jìn)行配準(zhǔn)拼接，檢驗改進(jìn)算法的有效性。圖9為其中一對顯微圖像的拼接效果圖，拼接無明顯錯位，配準(zhǔn)融合效果良好，黑邊為兩幅顯微圖像在垂直方向的波動。

圖8 圖像對比

圖9 病理顯微圖像配準(zhǔn)拼接結(jié)果

4 結(jié)語

針對病理顯微圖像快速高精度配準(zhǔn)的應(yīng)用，本文提出一種結(jié)合對數(shù)差值模板匹配和局部相位相關(guān)法的顯微圖像亞像素配準(zhǔn)方法。并且針對病理切片掃描過程中時常出現(xiàn)的空白區(qū)域，通過信息量評估來進(jìn)行規(guī)避。本文方法雖然在創(chuàng)建H型對數(shù)差值模板和匹配過程中增加了一定的計算量，但是由于對數(shù)差值模板的特性，僅是增加了查表和減法運算，并通過對關(guān)鍵點的篩選，進(jìn)一步減少了計算量。其增加的計算量小于由圖像尺寸增加而帶來的相位相關(guān)法所增加的運算，因此改進(jìn)的病理顯微圖像配準(zhǔn)算法速度有明顯提升。經(jīng)實驗驗證，本文算法配準(zhǔn)精度可達(dá)0.01像素，速度為相位相關(guān)法的3.7倍(圖像大小為2448×2048)，且圖像尺寸越大速度優(yōu)勢越明顯。因此，改進(jìn)的配準(zhǔn)方法更適用于病理切片掃描過程中海量的大尺寸顯微圖像配準(zhǔn)的應(yīng)用需求。