為獲得大視野高分辨率腎小球圖像, 需對多個小視野高分辨率圖像進行拼接。由于灰度不均勻和幾何形變兩種畸變的影響, 導致圖像拼接的效果不理想, 對病理診斷難以起到實質性的協助, 為此本文提出了腎小球透射電鏡(TEM)圖像的灰度與形變校正方法。本文首先描述了兩種畸變對圖像的影響, 并根據TEM成像系統原理提出TEM圖像獲取模型和畸變校正模型, 其次根據該模型分別利用灰度不均勻校正和幾何形變校正的方法有效減小了兩種畸變的影響。實驗結果表明在拼接前用本文方法進行預處理后, 拼接效果明顯提升, 有效減少了由圖像畸變引起的條紋、模糊、偽影等。本文方法能有效解決腎小球TEM圖像的畸變問題, 有助于提高多個腎小球TEM圖像的拼接效果。
引用本文： 鄔卓彬, 李穆, 路艷蒙, 曹蕾. 腎小球透射電鏡圖像的灰度與形變畸變校正. 生物醫學工程學雜志, 2016, 33(4): 755-761. doi: 10.7507/1001-5515.20160123 復制
引言
透射電鏡(transmission electron microscope，TEM)能觀察到腎小球的超微結構，絕大多數腎臟疾病的確診和新疾病的發現，都需要在TEM的協助下進行。為協助病理專家對腎小球進行全面正確的病理診斷，需要獲得高放大倍數下腎小球的完整圖像，即需要對多張高分辨率小視野的TEM腎小球圖像進行拼接。由于成像系統的幾何形變和灰度不均勻等缺陷會引起拼接圖像產生條紋、模糊、偽影等問題，直接拼接圖像并不能滿足臨床診斷的要求。但是目前相關研究僅對一種畸變進行處理，如Tasdizen等[1]僅關注圖像灰度不均勻，而Kaynig等[2]僅針對幾何形變提出了相關校正模型和方法。本文在此基礎上根據TEM的成像原理提出圖像獲取模型與校正模型，并據此模型結合灰度不均勻校正和幾何形變校正的方法以改善畸變對圖像的影響，提高腎小球圖像拼接的質量。
Meek等[3]提出TEM燈絲的老化、錯誤的參考電壓、電鏡內部的污染以及電子束等都會導致TEM圖像的灰度不均勻。在表征上，這些現象和在核磁共振成像(magnetic resonance image, MRI)圖像中觀察到的灰度不均勻極為相似。目前，針對MRI圖像的灰度不均勻已經做了大量的研究[4]。Styner等[5]利用了MRI圖像的腦組織灰度信息對圖像進行處理，Samsonov等[6]提出了先使用各向異性擴散濾波平滑和閾值去除邊緣像素進行預處理，再根據光照模型擬合灰度不均勻場。由于TEM圖像沒有腦組織種類信息，且是一種復雜的紋理圖像而并非分塊常量圖像，并不適用于各向異性擴散濾波平滑。因而Tasdizen等[1]在Samsonov等[6]提出的模型和擬合方法基礎上，使用高斯濾波取代各向異性擴散濾波平滑，用邊緣像素權重取代閾值，以解決TEM圖像中灰度不均勻問題。Scherzer等[7]則提出TEM的電子透鏡中存在軸向畸變、橢圓畸變等幾何形變。現在普遍認同且成熟的幾何形變校正方法是利用已知結構的規則晶體樣品的圖像進行誤差計算，以校正TEM的幾何形變[8-10]，該方法準確度極高，但是也存在操作難度高、出錯率高、校正頻繁等缺點。Koshevoy等[11]提出利用多張有序的TEM圖像重疊部分的Harris角點信息得出圖像形變校正場以校正圖像幾何形變。該方法簡單高效，但校正精度不高。Kaynig等[2]對Koshevoy的方法進行改進，構建了一個簡單的TEM成像系統模型，并使用尺度不變特征變換(scale-invariant feature transform，SIFT)特征點代替Harris角點估計校正場，繼而提高了校正精度，同時省去了人工確認圖像順序的步驟。
實際應用中，TEM圖像往往是灰度不均勻和幾何形變兩種畸變同時存在，因此本文參考Kaynig等[2]文中的TEM成像系統模型，構建出結合灰度不均勻和幾何形變的TEM圖像獲取模型與校正模型，并基于此模型，對成像過程中的兩類圖像畸變進行畸變校正，從而減少圖像畸變所造成的影響，改善圖像拼接的效果。本文首先詳細介紹了圖像畸變的現象、系統的建模及畸變校正的方法，再通過實驗結果和參數的討論證明了本文方法的有效性，并在最后對全文進行總結。
1 模型與方法
1.1 圖像畸變現象
本文用于拼接的腎小球橫切面圖來源于南方醫科大學電子顯微鏡室。這些實驗數據具有20%~30%重疊區域，主要有兩類畸變，分別為針對圖像像素點灰度值的灰度畸變和圖像像素點位置的形變畸變。這兩種畸變對單幅圖像所構成的影響有限，但是當需要獲得高分辨率大視野圖像的時候，這些畸變對拼接效果的影響就不能被忽略。對于腎小球TEM圖像，單幅圖像的灰度不均勻會造成拼接后圖像產生深色條紋。如圖 1所示，拼接圖像中顯示了6張圖像拼接的效果，增強圖像分別為拼接圖像中虛線框部分的原圖通過直方圖均衡增強后的圖像。如圖 1所示，黑色的小箭頭標示了拼接后圖像產生的深色條紋，通過觀察可以發現，原增強圖像的灰度值從中部往左右兩端緩慢下降。這些灰度不均勻問題在多張圖像拼接以后變得更明顯，影響整體圖像視覺效果。
圖1
原圖增強后的圖像和拼接圖像
Figure1.
Intensified images and mosaic image
單幅圖像形變畸變表現在圖像重疊部分像素點位置不能對齊。如圖 2所示，顯示了兩張相鄰圖像重疊部分的拼接效果。從虛線邊框區域內可以看出不能對齊的像素點使拼接的圖像產生模糊與偽影。
圖2
兩張相鄰圖像重疊部分的拼接效果
Figure2.
Mosaic outcome of overlaping area from two neighbouring images
1.2 圖像獲取模型
本文根據TEM成像系統的原理提出了TEM圖像獲取模型以描述圖像獲取過程中畸變所構成的影響。TEM圖像是由原始樣片通過成像系統一系列變換所得，為了有效描述系統中的變換，建立了如圖 3所示的TEM圖像獲取、校正模型。
圖3
TEM圖像獲取模型和校正模型
Figure3.
Acquisition model and correction model of TEM images
首先將整塊樣片S所在的平面定義為原始樣片坐標系統。從TEM系統中獲取圖像的過程相當于將樣片S的原始圖像坐標通過變換A1投影到圖像坐標系統，成為圖像P1。由于TEM圖像視野的偏移會引起圖像的平移，而電子的螺旋運動會造成圖像的旋轉變換，放大倍數的設定對應于圖像的尺度變化，因此從原始坐標系統到圖像坐標系統的變化，A1近似于仿射變換。其次，由于電子透鏡的軸向畸變以及橢圓畸變等對圖像的影響使得圖像產生一定的形變畸變，P1經過形變畸變場D的影響后得到圖像P2。上述變換可以用式(1)表示：
| ${P_{2L}} = D\left( {{A_1}{S_L}} \right)$ |
在上式中，SL表示樣品S的所有圖像像素點坐標，SL=[…xs… …ys… …1…]T，xs是像素點的橫坐標，ys是像素點的縱坐標，A1是3×3的仿射變換矩陣，D是一個非線性形變畸變場，P2L則是通過上述轉換的像素點坐標。
由于TEM燈絲的老化、錯誤的參考電壓等的影響會使得圖像產生一定的灰度不均勻，其中灰度不均勻畸變場由I表示，再加上加性噪聲N，得到最終的圖像F。該變換主要引起像素點灰度值的改變，可用式(2)表示：
| ${F_I}\left( {x,y} \right) = {P_{2I}}\left( {x,y} \right)I\left( {x,y} \right) + N\left( {x,y} \right)$ |
在式(2)中P2I(x, y)表示在圖像坐標點(x, y)處的灰度值，I(x, y)表示灰度不均勻畸變場，它會引起圖像坐標點(x, y)灰度的變化，N(x, y)表示加性系統噪聲。
本文圖像的校正過程如圖 3所示，根據所獲得的圖像F，先利用高斯低通濾波器去除圖像的系統噪聲，再結合所估計的灰度不均勻畸變場得到灰度校正后的圖像，最后通過形變校正場結合仿射變換矩陣得到最終校正后的圖像，下文將分別予以介紹。
1.3 圖像灰度不均勻畸變
為有效解決圖像中的灰度不均勻問題，本文根據圖像梯度，采用低維度多項式擬合乘性灰度不均勻畸變場的方法處理TEM圖像中的非線性灰度不均勻問題。TEM圖像梯度由以下的3個部分構成：①結構清晰的邊緣，表現為具有較大梯度、空間不平滑但有一定規律的部分；②系統噪聲，表現為具有較大梯度、空間不平滑但沒有規律的部分；③灰度不均勻畸變場I，表現為梯度小、變化緩慢(低頻)的部分。由此，本文對圖像的梯度進行擬合，根據式(2)和圖 3的圖像校正模型，對樣品結構邊緣和系統噪聲用高斯低通濾波器濾除(“*”表示卷積)，然后利用對數的性質，分離乘性灰度不均勻畸變場，最后通過計算圖像的梯度得到式(3)：
| $\begin{gathered} \nabla \left\{ {\ln \left[{{F_1}\left( {x,y} \right) * {G_\sigma }} \right]} \right\} = \hfill \\ \nabla \left\{ {\ln \left[{{P_{2I}}\left( {x,y} \right)I\left( {x,y} \right) + N\left( {x,y} \right) * {G_\sigma }} \right]} \right\} \hfill \\ \end{gathered} $ |
由于低通濾波對低頻的灰度不均勻畸變場I影響幾乎可以忽略，而高頻的噪聲場被濾除，再利用對數分離乘性灰度不均勻畸變場可得：
| $\nabla \ln {F_{1\sigma }}\left( {x,y} \right) \approx \nabla \ln {P_{2I\sigma }}\left( {x,y} \right) + \nabla \ln I\left( {x,y} \right)$ |
本文使用一個乘性K1維雙變量多項式灰度不均勻畸變場來模擬TEM圖像灰度不均勻畸變場I[1]，為保證乘性的灰度不均勻畸變場總是一個非負的形式，(x, y)表示為：
| $\hat I\left( {x,y;\Gamma } \right) = e\sum\limits_{i = 0}^{{K_1}} {\sum\limits_{j = 0}^i {{\gamma _{i - j,{j^{{x^{i - j}}{y^i}}}}}} } $ |
其中(x, y)表示任意像素點的坐標，Г是多項式參數γi-j, j的向量表示形式，γi-j, j是多項式中每一項的權重，K1表示多項式的維度。
結合式(4)和式(5)，用如下罰函數求出灰度不均勻畸變場的參數Г：
| $E\left( \Gamma \right) = \sum\limits_{y = 1}^M {\sum\limits_{x = 1}^M {{{\left| {\nabla \ln {F_{1\sigma }}\left( {x,y} \right) - \nabla \ln \hat I\left( {x,y;\Gamma } \right)} \right|}^2}} } $ |
其中M，N分別為圖像的行和列。
對式(6)進行求解就能求出對應圖像的灰度不均勻畸變場(x, y;Γ)。最后根據式(2)和圖 3的圖像校正模型，在忽略系統噪聲的情況下，可通過下式求得灰度校正后的圖像：
| ${P_{2I}}\left( {x,y} \right) \approx \widehat {{P_{2I}}}\left( {x,y} \right) = \frac{{F\left( {x,y} \right)}}{{\hat I\left( {x,y;\Gamma } \right)}}$ |
1.4 圖像幾何形變畸變
對于圖像幾何形變畸變，本文并沒有使用已知規則晶體樣品的圖像計算誤差[8-10]，而是使用多個待校正的，具有較大重疊面積的圖像獲取形變校正場。因為相鄰圖像重疊部分的重疊點(簡稱“重疊點”)的坐標在忽略形變畸變的情況下能近似對齊，可以據此估算仿射變換，同時這也是相鄰圖像之間的仿射變換。通過對齊圖像，然后利用對齊時重疊點的歐式距離及多項式擬合的形變校正場構建能量函數，求出能量最小時的形變校正場，最后用形變校正場校正所有圖像的幾何形變畸變。
進行圖像像素點的坐標轉換，如式(8)所示：
| $\begin{gathered} \left[\begin{gathered} {x_{\widehat {{P_2}'}}} \hfill \\ {y_{\widehat {{P_2}'}}} \hfill \\ 1 \hfill \\ \end{gathered} \right] = \hat C\left( {\left[\begin{gathered} {x_{\widehat {{P_2}}}} \hfill \\ {y_{\widehat {{P_2}}}} \hfill \\ 1 \hfill \\ \end{gathered} \right]} \right) = {\left[{{\Phi _{{K_2}}}\left( {\left[\begin{gathered} {x_{\widehat {{P_2}}}} \hfill \\ {y_{\widehat {{P_2}}}} \hfill \\ 1 \hfill \\ \end{gathered} \right]} \right)\beta } \right]^T} \hfill \\ = \left[\begin{gathered} {\beta _{1,1}} + {\beta _{1,2}}{x_{\widehat {{P_2}}}} + {\beta _{1,3}}{y_{\widehat {{P_2}}}} + {\beta _{1,4}}{x_{{{\widehat {{P_2}}}^2}}} + {\beta _{1,5}}{x_{\widehat {{P_2}}}}{y_{\widehat {{P_2}}}} \hfill \\ + {\beta _{1,6}}{y_{\widehat {{P_2}}}}^2 + \cdots + {\beta _{1,\frac{{\left( {{K_2} + 1} \right)\left( {{K_2} + 2} \right)}}{2}}}{y_{\widehat {{P_2}}}}{K_2} \hfill \\ {\beta _{2,1}} + {\beta _{2,2}}{x_{\widehat {{P_2}}}} + {\beta _{2,3}}{y_{\widehat {{P_2}}}} + {\beta _{2,4}}{x_{{{\widehat {{P_2}}}^2}}} + {\beta _{2,5}}{x_{\widehat {{P_2}}}}{y_{\widehat {{P_2}}}} \hfill \\ + {\beta _{2,6}}{y_{\widehat {{P_2}}}}^2 + \cdots + {\beta _{2,\frac{{\left( {{K_2} + 1} \right)\left( {{K_2} + 2} \right)}}{2}}}{y_{\widehat {{P_2}}}}{K_2} \hfill \\ \end{gathered} \right] \hfill \\ \end{gathered} $ |
其中為校正后的坐標，為校正前的坐標，代表校正場，代表K2維的多項式展開核，展開后有1×[(K2+1)(K2+2)/2]項，β代表多項式的參數矩陣，共有[(K2+1)(K2+2)/2]×3項。
如圖 3所示，圖像獲取模型中，若無形變畸變場D的影響，多張相鄰且含重疊區域的圖像P1經過仿射變換后便能在樣品坐標系準確對齊。而對于圖像校正模型中，圖像已受到形變畸變場D影響，可以利用相鄰圖像中的重疊點對的歐式距離作為限制，同時為保證圖像重疊點校正前后坐標不至于偏離太遠，通過式(9)的能量函數估計形變校正場的系數β：
| $\begin{gathered} \mathop {\min }\limits_{\beta ,R,T} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {\left( {{{\left\| \begin{gathered} \left( {\hat C\left( {{L^{i,j}}} \right){R^i} + {T^i}} \right) - \hfill \\ \left( {\hat C\left( {{L^{i,j}}} \right){R^i} + {T^i}} \right) \hfill \\ \end{gathered} \right\|}^2} + \lambda {{\left\| {\hat C\left( {{L^{i,j}}} \right) - {L^{i,j}}} \right\|}^2}} \right)} \hfill \\ = \mathop {\min }\limits_{\beta ,R,T} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {\left( \begin{gathered} {\left\| {\left( {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta {R^i} + {T^i}} \right) - \left( {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta {R^i} + {T^j}} \right)} \right\|^2} \hfill \\ + \lambda {\left\| {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta - {L^{i,j}}} \right\|^2} \hfill \\ \end{gathered} \right)} \hfill \\ \end{gathered} $ |
其中，圖像的數量是Num，λ控制圖像重疊點校正前后坐標偏離程度，Li, j表示圖i中與圖j存在重疊點的坐標，其在圖j中的對應重疊點為Li, jA2i，Ri和Ti分別為對應于第i個圖像的仿射變換分離出來的旋轉矩陣和平移矩陣，表示圖像i從圖像坐標系統投影到原始樣品坐標系統所對應的仿射矩陣的估計值，可以認為是從原始樣品坐標系統到圖像坐標系統的反投影變換。
而Ri和Ti對應的求解過程正如開始所說：假設圖像不存在幾何形變畸變，則重疊點對在同一坐標下能近似對齊，即滿足式(10)：
| $\mathop {\min }\limits_{{A_2}} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {{{\left\| {{L_{i,j}}\widehat {A_2^i} - {L_{j,i}}\widehat {A_2^i}} \right\|}^2}} $ |
由于重疊部分的重疊點對數量極多，如果將所有的重疊點進行匹配并運算并不實際。由于SIFT特征點具有仿射不變性、光照不變性等特點[4, 13]，且對于有鏡頭畸變及噪聲影響的TEM圖像所匹配的特征點對正確率也高，所以本文使用了SIFT特征點作為部分重疊點進行匹配。
通過式(9)求出參數β得到校正場，最后根據式(11)得到樣片圖像的估計值：
| $\widehat {{S_L}} = \widehat C\left( {\widehat {{P_2}}} \right)\widehat {{A_2}}$ |
2 結果與討論
首先，本實驗對比了使用本文方法前后的腎小球橫切面拼接圖像的結果。如圖 4所示，直接拼接圖像中展示的是6張TEM腎小球圖像直接拼接的結果，圖像上的黑色箭頭標記了單幅圖像的灰度不均勻導致的深色條紋。使用本文方法進行校正處理后，黑色箭頭的位置與直接拼接圖像中保持一致，但是相比直接拼接的結果，本文方法拼接后深色條紋明顯減弱，不存在比較明顯的灰度不均勻。局部放大圖(×10倍)分別對應于完整圖片中白色方框區域，其中直接拼接的局部放大圖中的結構出現了明顯的偽影，這是由于圖像幾何形變畸變造成圖像對齊誤差過大而導致的，與本文方法的局部放大圖對比可以得出用本文方法處理后圖中的偽影得到明顯改善，圖像結構清晰可見，圖像拼接效果明顯提高。
圖4
圖像畸變校正前后的結果及局部放大圖(×10)
Figure4.
Outcomes of images bofore and after distortion correction and images of zoom-in (×10)
2.1 灰度不均勻校正
2.1.1 灰度不均勻畸變場維度K1
圖像的灰度不均勻畸變場是一個緩慢變化的低頻場，多項式的自由度能決定該場是否穩定[1]。在式(5)中的維度K1過高或過低都會使灰度不均勻畸變場對圖像的一些高對比度結構產生錯誤校正。取不同的維度K1時，灰度不均勻校正的結果，如圖 5所示。其中左邊第一幅為原圖像，往右依次為K1=2、K1=3和K1=5的處理效果。與原圖比較K1=2、K1=5雖然減弱了圖中部與四周的灰度差異，但是右上角都出現了不同程度的過度校正。而K1=3的處理效果中部灰度和右上角的灰度與整幅圖像灰度一致且比較均勻。因此，本文建議使用K1=3來擬合多項式的灰度不均勻畸變場。
圖5
不同K1值的灰度校正結果
Figure5.
Outcomes of Different value of K1 for gray scale distortion correction
2.1.2 梯度算子?f
在式(3)中，圖像的梯度需要求圖像的一階導數，而不同的圖像梯度算子也會對灰度不均勻的效果造成影響。本文中使用了如下三種一階導數算子：
算子一：Gx=[1-1]；
算子二：Gx=[1 0-1]；
算子三：
| $ {G_x} = \left[\begin{gathered} 1\;\;\;0\;\;\; - 1 \hfill \\ 2\;\;\;0\;\;\; - 2 \hfill \\ 1\;\;\;0\;\;\; - 1 \hfill \\ \end{gathered} \right], $ |
所有的Gy=GxT。如圖 6所示，比較了上述一階導數算子對灰度不均勻校正效果的影響，圖像都通過直方圖均衡增強以提高對比度。
圖6
不同算子的灰度校正結果
Figure6.
Outcomes of different operator in gray scale distortion correction
通過觀察可以發現，使用算子二以及算子三的方法計算梯度時，圖像的左上角和右上角都因為過度校正，而使圖像灰度值偏高，但使用算子一的方法則沒有出現這個現象，因此本文中建議使用算子一作為圖像的梯度算子。
2.2 幾何形變校正
2.2.1 形變校正場維度K2
幾何形變校正中同樣需要對多項式的維度K2進行選取，在式(8)中，增大維度K2能給予多項式更高的自由度，能表示更復雜的圖像形變校正場。本文用圖像重疊部分像素灰度的相關系數和梯度的相關系數來評價維度K2對幾何形變校正的效果，灰度相關系數越大說明重疊部分圖像差異越小，重疊部分更準確，而梯度的相關系數越大則說明人眼更為敏感的邊緣部分差異更小。
如圖 7所示，顯示了不同的維度K2下幾何形變校正的圖像相關系數和梯度的相關系數。
圖7
維度K2與相關系數、梯度相關系數的關系曲線
Figure7.
Relation curves of K2-Correlation and K2-Correlation of Gradient
當維度K2大于3以后，變化趨向平穩，而過高的維度雖然在一定程度上提高圖像重疊部分的相關性，但同時會增加計算量，綜合考慮計算時間和圖像重疊部分的相關性和幾何形變校正的效果，本文建議使用維度K2=5。
3 結論
TEM的高分辨率圖像能提供腎小球內部組織和結構信息，而大視野圖像卻能提供全面的信息。為了同時滿足以上2個苛刻的要求，最簡單而直接的方法就是對多個高分辨率圖像進行拼接。然而，TEM圖像畸變嚴重的影響了拼接圖像信息的準確性，為此本文根據TEM成像系統提出了圖像獲取模型，并針對畸變在成像模型中的影響，與之相對的提出了一個圖像校正模型。本文的方法能有效減少圖像的模糊、偽影，最終得到的拼接圖像，可以提供大視野、完整的腎小球結構，為腎小球病理診斷提供便利。目前該方法僅僅在腎小球圖像拼接中進行了應用，而本文方法的有效性不僅僅如此，由于本文方法進行校正后能很好地保持位置和灰度信息的連續性，保留圖像的細節信息，有效針TEM圖像成像系統中存在的畸變。因此該方法可以作為TEM圖像分割、配準及三維重建等研究領域的預處理，其應用范圍廣泛。期待將來能將該方法推廣到其他TEM圖像中進行應用，還原出真實完整的圖像。
引言
透射電鏡(transmission electron microscope，TEM)能觀察到腎小球的超微結構，絕大多數腎臟疾病的確診和新疾病的發現，都需要在TEM的協助下進行。為協助病理專家對腎小球進行全面正確的病理診斷，需要獲得高放大倍數下腎小球的完整圖像，即需要對多張高分辨率小視野的TEM腎小球圖像進行拼接。由于成像系統的幾何形變和灰度不均勻等缺陷會引起拼接圖像產生條紋、模糊、偽影等問題，直接拼接圖像并不能滿足臨床診斷的要求。但是目前相關研究僅對一種畸變進行處理，如Tasdizen等[1]僅關注圖像灰度不均勻，而Kaynig等[2]僅針對幾何形變提出了相關校正模型和方法。本文在此基礎上根據TEM的成像原理提出圖像獲取模型與校正模型，并據此模型結合灰度不均勻校正和幾何形變校正的方法以改善畸變對圖像的影響，提高腎小球圖像拼接的質量。
Meek等[3]提出TEM燈絲的老化、錯誤的參考電壓、電鏡內部的污染以及電子束等都會導致TEM圖像的灰度不均勻。在表征上，這些現象和在核磁共振成像(magnetic resonance image, MRI)圖像中觀察到的灰度不均勻極為相似。目前，針對MRI圖像的灰度不均勻已經做了大量的研究[4]。Styner等[5]利用了MRI圖像的腦組織灰度信息對圖像進行處理，Samsonov等[6]提出了先使用各向異性擴散濾波平滑和閾值去除邊緣像素進行預處理，再根據光照模型擬合灰度不均勻場。由于TEM圖像沒有腦組織種類信息，且是一種復雜的紋理圖像而并非分塊常量圖像，并不適用于各向異性擴散濾波平滑。因而Tasdizen等[1]在Samsonov等[6]提出的模型和擬合方法基礎上，使用高斯濾波取代各向異性擴散濾波平滑，用邊緣像素權重取代閾值，以解決TEM圖像中灰度不均勻問題。Scherzer等[7]則提出TEM的電子透鏡中存在軸向畸變、橢圓畸變等幾何形變。現在普遍認同且成熟的幾何形變校正方法是利用已知結構的規則晶體樣品的圖像進行誤差計算，以校正TEM的幾何形變[8-10]，該方法準確度極高，但是也存在操作難度高、出錯率高、校正頻繁等缺點。Koshevoy等[11]提出利用多張有序的TEM圖像重疊部分的Harris角點信息得出圖像形變校正場以校正圖像幾何形變。該方法簡單高效，但校正精度不高。Kaynig等[2]對Koshevoy的方法進行改進，構建了一個簡單的TEM成像系統模型，并使用尺度不變特征變換(scale-invariant feature transform，SIFT)特征點代替Harris角點估計校正場，繼而提高了校正精度，同時省去了人工確認圖像順序的步驟。
實際應用中，TEM圖像往往是灰度不均勻和幾何形變兩種畸變同時存在，因此本文參考Kaynig等[2]文中的TEM成像系統模型，構建出結合灰度不均勻和幾何形變的TEM圖像獲取模型與校正模型，并基于此模型，對成像過程中的兩類圖像畸變進行畸變校正，從而減少圖像畸變所造成的影響，改善圖像拼接的效果。本文首先詳細介紹了圖像畸變的現象、系統的建模及畸變校正的方法，再通過實驗結果和參數的討論證明了本文方法的有效性，并在最后對全文進行總結。
1 模型與方法
1.1 圖像畸變現象
本文用于拼接的腎小球橫切面圖來源于南方醫科大學電子顯微鏡室。這些實驗數據具有20%~30%重疊區域，主要有兩類畸變，分別為針對圖像像素點灰度值的灰度畸變和圖像像素點位置的形變畸變。這兩種畸變對單幅圖像所構成的影響有限，但是當需要獲得高分辨率大視野圖像的時候，這些畸變對拼接效果的影響就不能被忽略。對于腎小球TEM圖像，單幅圖像的灰度不均勻會造成拼接后圖像產生深色條紋。如圖 1所示，拼接圖像中顯示了6張圖像拼接的效果，增強圖像分別為拼接圖像中虛線框部分的原圖通過直方圖均衡增強后的圖像。如圖 1所示，黑色的小箭頭標示了拼接后圖像產生的深色條紋，通過觀察可以發現，原增強圖像的灰度值從中部往左右兩端緩慢下降。這些灰度不均勻問題在多張圖像拼接以后變得更明顯，影響整體圖像視覺效果。
圖1
原圖增強后的圖像和拼接圖像
Figure1.
Intensified images and mosaic image
單幅圖像形變畸變表現在圖像重疊部分像素點位置不能對齊。如圖 2所示，顯示了兩張相鄰圖像重疊部分的拼接效果。從虛線邊框區域內可以看出不能對齊的像素點使拼接的圖像產生模糊與偽影。
圖2
兩張相鄰圖像重疊部分的拼接效果
Figure2.
Mosaic outcome of overlaping area from two neighbouring images
1.2 圖像獲取模型
本文根據TEM成像系統的原理提出了TEM圖像獲取模型以描述圖像獲取過程中畸變所構成的影響。TEM圖像是由原始樣片通過成像系統一系列變換所得，為了有效描述系統中的變換，建立了如圖 3所示的TEM圖像獲取、校正模型。
圖3
TEM圖像獲取模型和校正模型
Figure3.
Acquisition model and correction model of TEM images
首先將整塊樣片S所在的平面定義為原始樣片坐標系統。從TEM系統中獲取圖像的過程相當于將樣片S的原始圖像坐標通過變換A1投影到圖像坐標系統，成為圖像P1。由于TEM圖像視野的偏移會引起圖像的平移，而電子的螺旋運動會造成圖像的旋轉變換，放大倍數的設定對應于圖像的尺度變化，因此從原始坐標系統到圖像坐標系統的變化，A1近似于仿射變換。其次，由于電子透鏡的軸向畸變以及橢圓畸變等對圖像的影響使得圖像產生一定的形變畸變，P1經過形變畸變場D的影響后得到圖像P2。上述變換可以用式(1)表示：
| ${P_{2L}} = D\left( {{A_1}{S_L}} \right)$ |
在上式中，SL表示樣品S的所有圖像像素點坐標，SL=[…xs… …ys… …1…]T，xs是像素點的橫坐標，ys是像素點的縱坐標，A1是3×3的仿射變換矩陣，D是一個非線性形變畸變場，P2L則是通過上述轉換的像素點坐標。
由于TEM燈絲的老化、錯誤的參考電壓等的影響會使得圖像產生一定的灰度不均勻，其中灰度不均勻畸變場由I表示，再加上加性噪聲N，得到最終的圖像F。該變換主要引起像素點灰度值的改變，可用式(2)表示：
| ${F_I}\left( {x,y} \right) = {P_{2I}}\left( {x,y} \right)I\left( {x,y} \right) + N\left( {x,y} \right)$ |
在式(2)中P2I(x, y)表示在圖像坐標點(x, y)處的灰度值，I(x, y)表示灰度不均勻畸變場，它會引起圖像坐標點(x, y)灰度的變化，N(x, y)表示加性系統噪聲。
本文圖像的校正過程如圖 3所示，根據所獲得的圖像F，先利用高斯低通濾波器去除圖像的系統噪聲，再結合所估計的灰度不均勻畸變場得到灰度校正后的圖像，最后通過形變校正場結合仿射變換矩陣得到最終校正后的圖像，下文將分別予以介紹。
1.3 圖像灰度不均勻畸變
為有效解決圖像中的灰度不均勻問題，本文根據圖像梯度，采用低維度多項式擬合乘性灰度不均勻畸變場的方法處理TEM圖像中的非線性灰度不均勻問題。TEM圖像梯度由以下的3個部分構成：①結構清晰的邊緣，表現為具有較大梯度、空間不平滑但有一定規律的部分；②系統噪聲，表現為具有較大梯度、空間不平滑但沒有規律的部分；③灰度不均勻畸變場I，表現為梯度小、變化緩慢(低頻)的部分。由此，本文對圖像的梯度進行擬合，根據式(2)和圖 3的圖像校正模型，對樣品結構邊緣和系統噪聲用高斯低通濾波器濾除(“*”表示卷積)，然后利用對數的性質，分離乘性灰度不均勻畸變場，最后通過計算圖像的梯度得到式(3)：
| $\begin{gathered} \nabla \left\{ {\ln \left[{{F_1}\left( {x,y} \right) * {G_\sigma }} \right]} \right\} = \hfill \\ \nabla \left\{ {\ln \left[{{P_{2I}}\left( {x,y} \right)I\left( {x,y} \right) + N\left( {x,y} \right) * {G_\sigma }} \right]} \right\} \hfill \\ \end{gathered} $ |
由于低通濾波對低頻的灰度不均勻畸變場I影響幾乎可以忽略，而高頻的噪聲場被濾除，再利用對數分離乘性灰度不均勻畸變場可得：
| $\nabla \ln {F_{1\sigma }}\left( {x,y} \right) \approx \nabla \ln {P_{2I\sigma }}\left( {x,y} \right) + \nabla \ln I\left( {x,y} \right)$ |
本文使用一個乘性K1維雙變量多項式灰度不均勻畸變場來模擬TEM圖像灰度不均勻畸變場I[1]，為保證乘性的灰度不均勻畸變場總是一個非負的形式，(x, y)表示為：
| $\hat I\left( {x,y;\Gamma } \right) = e\sum\limits_{i = 0}^{{K_1}} {\sum\limits_{j = 0}^i {{\gamma _{i - j,{j^{{x^{i - j}}{y^i}}}}}} } $ |
其中(x, y)表示任意像素點的坐標，Г是多項式參數γi-j, j的向量表示形式，γi-j, j是多項式中每一項的權重，K1表示多項式的維度。
結合式(4)和式(5)，用如下罰函數求出灰度不均勻畸變場的參數Г：
| $E\left( \Gamma \right) = \sum\limits_{y = 1}^M {\sum\limits_{x = 1}^M {{{\left| {\nabla \ln {F_{1\sigma }}\left( {x,y} \right) - \nabla \ln \hat I\left( {x,y;\Gamma } \right)} \right|}^2}} } $ |
其中M，N分別為圖像的行和列。
對式(6)進行求解就能求出對應圖像的灰度不均勻畸變場(x, y;Γ)。最后根據式(2)和圖 3的圖像校正模型，在忽略系統噪聲的情況下，可通過下式求得灰度校正后的圖像：
| ${P_{2I}}\left( {x,y} \right) \approx \widehat {{P_{2I}}}\left( {x,y} \right) = \frac{{F\left( {x,y} \right)}}{{\hat I\left( {x,y;\Gamma } \right)}}$ |
1.4 圖像幾何形變畸變
對于圖像幾何形變畸變，本文并沒有使用已知規則晶體樣品的圖像計算誤差[8-10]，而是使用多個待校正的，具有較大重疊面積的圖像獲取形變校正場。因為相鄰圖像重疊部分的重疊點(簡稱“重疊點”)的坐標在忽略形變畸變的情況下能近似對齊，可以據此估算仿射變換，同時這也是相鄰圖像之間的仿射變換。通過對齊圖像，然后利用對齊時重疊點的歐式距離及多項式擬合的形變校正場構建能量函數，求出能量最小時的形變校正場，最后用形變校正場校正所有圖像的幾何形變畸變。
進行圖像像素點的坐標轉換，如式(8)所示：
| $\begin{gathered} \left[\begin{gathered} {x_{\widehat {{P_2}'}}} \hfill \\ {y_{\widehat {{P_2}'}}} \hfill \\ 1 \hfill \\ \end{gathered} \right] = \hat C\left( {\left[\begin{gathered} {x_{\widehat {{P_2}}}} \hfill \\ {y_{\widehat {{P_2}}}} \hfill \\ 1 \hfill \\ \end{gathered} \right]} \right) = {\left[{{\Phi _{{K_2}}}\left( {\left[\begin{gathered} {x_{\widehat {{P_2}}}} \hfill \\ {y_{\widehat {{P_2}}}} \hfill \\ 1 \hfill \\ \end{gathered} \right]} \right)\beta } \right]^T} \hfill \\ = \left[\begin{gathered} {\beta _{1,1}} + {\beta _{1,2}}{x_{\widehat {{P_2}}}} + {\beta _{1,3}}{y_{\widehat {{P_2}}}} + {\beta _{1,4}}{x_{{{\widehat {{P_2}}}^2}}} + {\beta _{1,5}}{x_{\widehat {{P_2}}}}{y_{\widehat {{P_2}}}} \hfill \\ + {\beta _{1,6}}{y_{\widehat {{P_2}}}}^2 + \cdots + {\beta _{1,\frac{{\left( {{K_2} + 1} \right)\left( {{K_2} + 2} \right)}}{2}}}{y_{\widehat {{P_2}}}}{K_2} \hfill \\ {\beta _{2,1}} + {\beta _{2,2}}{x_{\widehat {{P_2}}}} + {\beta _{2,3}}{y_{\widehat {{P_2}}}} + {\beta _{2,4}}{x_{{{\widehat {{P_2}}}^2}}} + {\beta _{2,5}}{x_{\widehat {{P_2}}}}{y_{\widehat {{P_2}}}} \hfill \\ + {\beta _{2,6}}{y_{\widehat {{P_2}}}}^2 + \cdots + {\beta _{2,\frac{{\left( {{K_2} + 1} \right)\left( {{K_2} + 2} \right)}}{2}}}{y_{\widehat {{P_2}}}}{K_2} \hfill \\ \end{gathered} \right] \hfill \\ \end{gathered} $ |
其中為校正后的坐標，為校正前的坐標，代表校正場，代表K2維的多項式展開核，展開后有1×[(K2+1)(K2+2)/2]項，β代表多項式的參數矩陣，共有[(K2+1)(K2+2)/2]×3項。
如圖 3所示，圖像獲取模型中，若無形變畸變場D的影響，多張相鄰且含重疊區域的圖像P1經過仿射變換后便能在樣品坐標系準確對齊。而對于圖像校正模型中，圖像已受到形變畸變場D影響，可以利用相鄰圖像中的重疊點對的歐式距離作為限制，同時為保證圖像重疊點校正前后坐標不至于偏離太遠，通過式(9)的能量函數估計形變校正場的系數β：
| $\begin{gathered} \mathop {\min }\limits_{\beta ,R,T} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {\left( {{{\left\| \begin{gathered} \left( {\hat C\left( {{L^{i,j}}} \right){R^i} + {T^i}} \right) - \hfill \\ \left( {\hat C\left( {{L^{i,j}}} \right){R^i} + {T^i}} \right) \hfill \\ \end{gathered} \right\|}^2} + \lambda {{\left\| {\hat C\left( {{L^{i,j}}} \right) - {L^{i,j}}} \right\|}^2}} \right)} \hfill \\ = \mathop {\min }\limits_{\beta ,R,T} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {\left( \begin{gathered} {\left\| {\left( {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta {R^i} + {T^i}} \right) - \left( {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta {R^i} + {T^j}} \right)} \right\|^2} \hfill \\ + \lambda {\left\| {\Phi {K_2}\left( {{L^{i,j}}} \right)\beta - {L^{i,j}}} \right\|^2} \hfill \\ \end{gathered} \right)} \hfill \\ \end{gathered} $ |
其中，圖像的數量是Num，λ控制圖像重疊點校正前后坐標偏離程度，Li, j表示圖i中與圖j存在重疊點的坐標，其在圖j中的對應重疊點為Li, jA2i，Ri和Ti分別為對應于第i個圖像的仿射變換分離出來的旋轉矩陣和平移矩陣，表示圖像i從圖像坐標系統投影到原始樣品坐標系統所對應的仿射矩陣的估計值，可以認為是從原始樣品坐標系統到圖像坐標系統的反投影變換。
而Ri和Ti對應的求解過程正如開始所說：假設圖像不存在幾何形變畸變，則重疊點對在同一坐標下能近似對齊，即滿足式(10)：
| $\mathop {\min }\limits_{{A_2}} \sum\limits_{i,j = 1\;\;j \ne i}^{Num} {{{\left\| {{L_{i,j}}\widehat {A_2^i} - {L_{j,i}}\widehat {A_2^i}} \right\|}^2}} $ |
由于重疊部分的重疊點對數量極多，如果將所有的重疊點進行匹配并運算并不實際。由于SIFT特征點具有仿射不變性、光照不變性等特點[4, 13]，且對于有鏡頭畸變及噪聲影響的TEM圖像所匹配的特征點對正確率也高，所以本文使用了SIFT特征點作為部分重疊點進行匹配。
通過式(9)求出參數β得到校正場，最后根據式(11)得到樣片圖像的估計值：
| $\widehat {{S_L}} = \widehat C\left( {\widehat {{P_2}}} \right)\widehat {{A_2}}$ |
2 結果與討論
首先，本實驗對比了使用本文方法前后的腎小球橫切面拼接圖像的結果。如圖 4所示，直接拼接圖像中展示的是6張TEM腎小球圖像直接拼接的結果，圖像上的黑色箭頭標記了單幅圖像的灰度不均勻導致的深色條紋。使用本文方法進行校正處理后，黑色箭頭的位置與直接拼接圖像中保持一致，但是相比直接拼接的結果，本文方法拼接后深色條紋明顯減弱，不存在比較明顯的灰度不均勻。局部放大圖(×10倍)分別對應于完整圖片中白色方框區域，其中直接拼接的局部放大圖中的結構出現了明顯的偽影，這是由于圖像幾何形變畸變造成圖像對齊誤差過大而導致的，與本文方法的局部放大圖對比可以得出用本文方法處理后圖中的偽影得到明顯改善，圖像結構清晰可見，圖像拼接效果明顯提高。
圖4
圖像畸變校正前后的結果及局部放大圖(×10)
Figure4.
Outcomes of images bofore and after distortion correction and images of zoom-in (×10)
2.1 灰度不均勻校正
2.1.1 灰度不均勻畸變場維度K1
圖像的灰度不均勻畸變場是一個緩慢變化的低頻場，多項式的自由度能決定該場是否穩定[1]。在式(5)中的維度K1過高或過低都會使灰度不均勻畸變場對圖像的一些高對比度結構產生錯誤校正。取不同的維度K1時，灰度不均勻校正的結果，如圖 5所示。其中左邊第一幅為原圖像，往右依次為K1=2、K1=3和K1=5的處理效果。與原圖比較K1=2、K1=5雖然減弱了圖中部與四周的灰度差異，但是右上角都出現了不同程度的過度校正。而K1=3的處理效果中部灰度和右上角的灰度與整幅圖像灰度一致且比較均勻。因此，本文建議使用K1=3來擬合多項式的灰度不均勻畸變場。
圖5
不同K1值的灰度校正結果
Figure5.
Outcomes of Different value of K1 for gray scale distortion correction
2.1.2 梯度算子?f
在式(3)中，圖像的梯度需要求圖像的一階導數，而不同的圖像梯度算子也會對灰度不均勻的效果造成影響。本文中使用了如下三種一階導數算子：
算子一：Gx=[1-1]；
算子二：Gx=[1 0-1]；
算子三：
| $ {G_x} = \left[\begin{gathered} 1\;\;\;0\;\;\; - 1 \hfill \\ 2\;\;\;0\;\;\; - 2 \hfill \\ 1\;\;\;0\;\;\; - 1 \hfill \\ \end{gathered} \right], $ |
所有的Gy=GxT。如圖 6所示，比較了上述一階導數算子對灰度不均勻校正效果的影響，圖像都通過直方圖均衡增強以提高對比度。
圖6
不同算子的灰度校正結果
Figure6.
Outcomes of different operator in gray scale distortion correction
通過觀察可以發現，使用算子二以及算子三的方法計算梯度時，圖像的左上角和右上角都因為過度校正，而使圖像灰度值偏高，但使用算子一的方法則沒有出現這個現象，因此本文中建議使用算子一作為圖像的梯度算子。
2.2 幾何形變校正
2.2.1 形變校正場維度K2
幾何形變校正中同樣需要對多項式的維度K2進行選取，在式(8)中，增大維度K2能給予多項式更高的自由度，能表示更復雜的圖像形變校正場。本文用圖像重疊部分像素灰度的相關系數和梯度的相關系數來評價維度K2對幾何形變校正的效果，灰度相關系數越大說明重疊部分圖像差異越小，重疊部分更準確，而梯度的相關系數越大則說明人眼更為敏感的邊緣部分差異更小。
如圖 7所示，顯示了不同的維度K2下幾何形變校正的圖像相關系數和梯度的相關系數。
圖7
維度K2與相關系數、梯度相關系數的關系曲線
Figure7.
Relation curves of K2-Correlation and K2-Correlation of Gradient
當維度K2大于3以后，變化趨向平穩，而過高的維度雖然在一定程度上提高圖像重疊部分的相關性，但同時會增加計算量，綜合考慮計算時間和圖像重疊部分的相關性和幾何形變校正的效果，本文建議使用維度K2=5。
3 結論
TEM的高分辨率圖像能提供腎小球內部組織和結構信息，而大視野圖像卻能提供全面的信息。為了同時滿足以上2個苛刻的要求，最簡單而直接的方法就是對多個高分辨率圖像進行拼接。然而，TEM圖像畸變嚴重的影響了拼接圖像信息的準確性，為此本文根據TEM成像系統提出了圖像獲取模型，并針對畸變在成像模型中的影響，與之相對的提出了一個圖像校正模型。本文的方法能有效減少圖像的模糊、偽影，最終得到的拼接圖像，可以提供大視野、完整的腎小球結構，為腎小球病理診斷提供便利。目前該方法僅僅在腎小球圖像拼接中進行了應用，而本文方法的有效性不僅僅如此，由于本文方法進行校正后能很好地保持位置和灰度信息的連續性，保留圖像的細節信息，有效針TEM圖像成像系統中存在的畸變。因此該方法可以作為TEM圖像分割、配準及三維重建等研究領域的預處理，其應用范圍廣泛。期待將來能將該方法推廣到其他TEM圖像中進行應用，還原出真實完整的圖像。
