一種改進的病理顯微圖像亞像素快速配準方法

，

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

0 引言

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

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

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

針對病理顯微圖像配準的應用特點，本文提出了一種新的快速配準方法。該方法利用對數差值函數對待匹配圖像進行信息評估，并構建模板獲取像素級粗定位，根據粗定位獲取待配準子圖，由相位相關法獲取亞像素級細定位。

1 亞像素配準技術

1.1 歸一化相位相關法

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

1.2 基于相位獲取亞像素配準技術

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

2 改進的病理顯微圖像亞像素配準算法

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

2.1 對數差值函數

用i表示圖像的像素灰度值，則其灰度值對數差值函數[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的常數放大系數。當i=0時，ln(i)無效，為了處理方便，用一個極小常量ε(0<ε<<1)代替，即u[0]=Klnε。

對數差值函數有具備以下特性：

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

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

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

2.2 病理顯微圖像亞像素配準算法的流程

第一步，構建H型對數差值模板。

以圖像的左上角為原點O(0,0)，以水平向右為X軸正方向，垂直向下為Y軸正方向，H型對數模板由如圖所示U、H、V三條線段組成，線段U和V的長度為l，線段H的長度為s。以線段U的第一個點P(x,y)確定模板的位置，線段H的第一個點坐標為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的各像素點的灰度值統一分別記為Ik,Jk,Dn,Ik,Jk,Dn∈[0,255],n=0,1,…,s-1,k=0,1,…l-1。

圖1 H型對數差值模板

使用對數差值函數分別計算線段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上對應位置像素間的灰度值變化，定義特征向量：

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

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

(9)

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

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

(10)

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

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

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

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

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

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

(11)

(12)

(13)

即表示在圖像A的搜索區域范圍內整體像素的灰度變化幅度小，可認定為是圖像搜索區域范圍內像素灰度值單調甚至為空白區域，無法提供足夠的灰度變化信息構建模板Tlsp(A)完成配準，即便構建模板或者進行相位相關配準其配準結果亦不具備可靠性，應當返回信息量不足的錯誤配準信息。

第三步，進行對數差值模板匹配獲取粗定位。

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

圖3 匹配示意圖

在顯微圖像B的待搜索區域內各點PB(x,y)依據模板Tlsp(A)的尺寸參數s和l構建H型對數差值模板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)為各特征向量差異值的權重值。

為了進一步減少計算量，設置閾值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)

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

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

第四步，獲取待配準子圖。

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

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

第五步，獲取亞像素級細定位。

對兩幅待匹配子圖a,b進行傅里葉變換并計算其相位相關歸一化互功率譜,在兩幅大小為M×N的子圖計算所得的互動率譜的中心位置(xc,yc)附近范圍構建大小為2c鄰域，在此鄰域區域內采用矩陣相乘計算互功率譜矩陣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)的中心變換形式。

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

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

(18)

算法流程圖如圖4所示。

圖4 病理顯微圖像配準算法流程圖

3 實驗結果與分析

3.1 實驗平臺

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

圖5 實驗平臺

3.2 實驗一

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

圖6 精度測試待配準顯微圖像組

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

3.3 實驗二

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

圖7 兩種方法配準80對顯微圖像耗時結果

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

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

圖8 圖像對比

圖9 病理顯微圖像配準拼接結果

4 結語

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