陳曉艷 李健楠 王化祥（天津科技大學電子信息與自動化學院，天津 300）

2（天津大學電氣與自動化工程學院，天津 300072）

引言

醫(yī)學圖像主要包括形態(tài)圖像、結(jié)構(gòu)圖像和功能圖像。醫(yī)學顯微圖像屬于形態(tài)圖像；X線圖像、超聲、CT、MRI等屬于結(jié)構(gòu)圖像；PET、SPECT、FMRT及EIT等為功能圖像［1-2］。由于成像原理不同，不同圖像提供的人體生理、病理信息也不同。將不同原理的圖像進行融合，獲得信息更為豐富的復(fù)合圖像，是當前臨床醫(yī)學圖像的研究熱點之一［3-4］。

電阻抗成像技術(shù)（electricalimpedance tomography，EIT）是近30年發(fā)展起來的新興成像技術(shù)，其無創(chuàng)、無輻射、成本低廉，能長時間連續(xù)工作，更重要的是EIT能反映人體組織、器官的功能變化，在醫(yī)學圖像監(jiān)護方面具有顯著的優(yōu)勢［5-6］。筆者在電阻抗成像技術(shù)基礎(chǔ)上，以人體呼吸過程的電阻抗圖像及胸部CT圖像融合為例，對EIT圖像和CT圖像的融合方法進行研究。在實現(xiàn)功能成像的同時，對阻抗變化敏感區(qū)域進行有效定位［7-8］。

1 原理和方法

1.1 提取結(jié)構(gòu)信息

以人體胸部CT圖像為例展開研究。讓被測者躺在16螺旋CT成像設(shè)備測試床上，在吸氣末時刻獲取CT圖像，如圖1所示。其中，胸腔體表外周的圓亮點是具有透射性的EIT電極，下面是測試床截面。由圖可知，CT圖像不僅能提供電阻抗成像所需要的胸腔外部邊界結(jié)構(gòu)信息，而且也能提供豐富的內(nèi)部結(jié)構(gòu)信息。具體提取方法是：先采用圓形區(qū)域均值濾波disk算子進行濾波，然后再采用canny邊緣提取算子提取圖像輪廓。

圖1 人體胸部CT圖像Fig.1 The chest CT image

設(shè)圖像上像素點（x，y）處的梯度f、幅值A(chǔ)及梯度方向φ分別為

幅值反映該點處亮度值變化的大小，梯度的方向反映該點在哪個方向上變化最大。

由于數(shù)字圖像是離散的，所以其圖像亮度值變化的大小要用差分來表示［9］。二維圖像在x方向上的一階差分為f（x+1，y）-f（x，y），在y方向上的一階差分為f（x，y+1）-f（x，y），圖像差分時的［33］鄰域模版為

f（x-1，y-1），f（x，y-1），f（x+1，y-1）

f（x-1，y），f（x，y），f（x+1，y）

f（x-1，y+1），f（x，y+1），f（x+1，y+1）

邊緣判定的依據(jù)是通過圖像像素亮度值的一階或二階導數(shù)進行判斷，當幅值A(chǔ)大于某一閾值時則為邊緣。

由于canny算子具有較強的魯棒性，不容易受噪聲干擾，而且能有效地提取閉合區(qū)間，適用于確定阻抗敏感區(qū)域，所以采用它進行邊緣提取。canny算子是先用高斯函數(shù)進行濾波，然后計算沿著x，y兩個方向的偏導數(shù)，求出梯度的幅值和梯度的方向。通過梯度的方向就可以找到這個像素梯度方向的鄰接像素，然后將梯度的幅值和高低閾值進行比較，大于高閾值的一定是邊緣，小于低閾值的一定不是邊緣，介于之間的要看這個像素的鄰接像素中有沒有梯度幅值超過高閾值的像素，如果有則為邊緣，否則就不是邊緣。筆者在用canny算子進行邊緣提取時，高低閾值的選取是由程序自動確定的。確定的原則就是：先通過設(shè)定非邊緣像素的比例來確定高閾值TH，非邊緣像素的比例一般選擇0.7；然后再確定低閾值TL，一般TL=0.4TH。先用disk算子進行濾波，再用canny算子提取的人體胸部輪廓，如圖2所示。

圖2 disk算子和canny算子處理效果Fig.2 The edge effect based on the operators of‘disk’and‘canny’

1.2 EIT圖像重建

保留圖2中提取的胸腔及肺部的邊界結(jié)構(gòu)，勾勒出感興趣區(qū)域左右肺臟及心臟的邊界，并將其導入Comsol平臺，如圖3所示。

先在Comsol軟件平臺上進行正問題求解，得到靈敏度系數(shù)矩陣，再將其導入Matlab，經(jīng)共軛梯度算法重建EIT圖像。共軛梯度算法是介于最速下降法與牛頓法之間的一個方法，僅需利用一階導數(shù)信息，即克服了最速下降法收斂慢的缺點，又避免了牛頓法需要存儲和計算Hessian矩陣并求逆的缺點。共軛梯度算法經(jīng)過第i次迭代得到的解xi連續(xù)逼近真解，迭代誤差ri=b-Axi，搜索方向設(shè)定為Pi，迭代次數(shù)與迭代誤差的更新為qi=Api，從而有

圖3 EIT正問題模型Fig.3 The forward problem model of EIT

式中，αi是標量，使得r（α）*A-1r（α）最小，其中r（α）=ri-1-αipi-1。所以有

搜索的方向為

保證Pi與所有Api正交，ri與所有的rj都正交（j＜i）。當?shù)`差小于設(shè)定閾值時，迭代結(jié)束。對于EIT成像，系數(shù)矩陣（Jacobian矩陣）往往不是對稱正定矩陣，因此需要對原方程進行轉(zhuǎn)換，有

由于A并非列滿秩矩陣，因此ATA也不是對稱正定矩陣，可以通過增加正則化矩陣的辦法來解決這個問題。當正則化矩陣為單位矩陣時，式（7）可以轉(zhuǎn)化為

此時，系數(shù)矩陣為對稱正定矩陣，因此EIT可用上述共軛梯度算法進行求解［10］。重建后的圖像如圖4所示。

1.3 基于小波變換的圖像融合

由于人的視網(wǎng)膜是在不同的頻率信道中進行不同的信號處理［11］，為了獲得EIT圖像和CT圖像良好的融合效果，需要對EIT圖像和CT圖像在不同的頻率信道上進行融合，而小波變換可實現(xiàn)將原始圖像分解成一系列具有不同頻域特性的子圖像，從而為進行不同頻率的圖像數(shù)據(jù)融合提供了有效方法。

圖4 EIT重建圖像Fig.4 EIT reconstruction image

從而有

這就是說，ψ（t）與整個橫坐標所圍成面積的代數(shù)和為0，因此ψ（t）的圖形應(yīng)是在橫坐標上下波動的“小波”。通常稱ψ（t）為母小波，則稱ψa，b（t）為小波函數(shù)，或簡稱為小波（wavelet）［12］。

令L2（R）為可測的、平方可積函數(shù)f（t）的矢量空間，R為實數(shù)集，對于任意的f（t）∈L2（R）的連續(xù)小波變換給出如下定義，即：

小波變換的重構(gòu)（逆變換）公式為

設(shè)EIT圖像和CT圖像分別為A和B，F(xiàn)為融合后的圖像，融合的基本步驟為：

1）對每一幅源圖像分別進行小波分解；

2）對各分解層分別進行融合處理，不同的分解層采用不同的融合算子；

3）對融合后的數(shù)據(jù)進行小波逆變換，所得到的重構(gòu)圖像即為融合圖像。

由于小波變換可以把被融合的圖像分解到多個頻帶中，且具有方向性，故基于小波變換的圖像融合可以在不同頻帶、不同方向上采用不同的融合規(guī)則，以充分挖掘被融合圖像的互補信息，突出感興趣的特征和細節(jié)信息，從而獲得良好的視覺效果［13-15］?；谛〔ㄗ儞Q的圖像融合流程如圖5所示。

圖5 基于小波變換的圖像融合流程圖Fig.5 The flow chart of image fusion based on the wavelet decomposition

2 結(jié)果

將圖2中胸腔邊界之外無用信息濾除后的圖像，與EIT重建圖像（見圖4）經(jīng)小波算法進行融合，融合效果如圖6所示。

圖6 EIT/CT融合圖像Fig.6 EIT/CT fusion image

通常，在電阻抗成像過程中，靈敏度矩陣的求解是在敏感區(qū)域內(nèi)部各點電導率相同的假設(shè)前提下獲得的。但在生物醫(yī)學電阻抗成像的應(yīng)用研究中，人體腦部、腹部、肺部等區(qū)域內(nèi)的電導率各不相同，而其分布是十分重要的先驗信息。

本研究充分利用結(jié)構(gòu)和電導率的先驗信息，獲得更準確的靈敏度矩陣，從而達到提高電阻抗成像質(zhì)量的目的。在此基礎(chǔ)上，采用小波算法，將CT結(jié)構(gòu)圖像與EIT功能圖像進一步融合。圖6與圖像融合前的圖4相比，EIT圖像中的電阻抗功能變化信息在結(jié)構(gòu)圖像中凸顯出來，可以實現(xiàn)電阻抗變化區(qū)域的初步定位。

3 討論和結(jié)論

電阻抗功能成像技術(shù)的優(yōu)勢顯而易見，但是如何充分發(fā)揮其優(yōu)勢，服務(wù)臨床、輔助診療是目前亟待解決的問題，筆者就是圍繞如何實現(xiàn)病灶的定位來開展探索性的研究。借助高精度的CT結(jié)構(gòu)圖像，不僅可以提高電阻抗圖像質(zhì)量，更有望融合后的圖像能提供定位性的生理和病理變化信息。

本研究探索了EIT圖像與CT圖像融合的有效方法，并獲得了良好的融合效果，為EIT成像技術(shù)與其他技術(shù)（如超聲、MRI、X線等，甚至其他功能成像技術(shù)）進行圖像融合奠定了基礎(chǔ)。但是，本研究僅限于兩種靜態(tài)圖像的融合，對于呼吸等動態(tài)變化的過程，融合技術(shù)和方法會有所不同，有待進一步深入研究。

筆者主要研究了CT圖像輪廓信息的提取方法，通過對各種濾波算子以及邊緣提取算子的綜合比較，表明用disk算子進行濾波，然后再用canny算子進行邊緣提取，所得到的圖像比較符合要求。提取的CT輪廓信息可以建立EIT重建模型，也可以與EIT圖像進行融合。本研究驗證了EIT圖像與CT圖像融合的可行性，探索了一種有效的圖像融合方法，復(fù)合圖像不僅可以把成像區(qū)域的結(jié)構(gòu)信息與功能信息同時顯示出來，實現(xiàn)信息互補，而且為進一步研究功能變化信息的定位問題開辟了新的途徑。

［1］任超世，李章勇，王妍，等.電阻抗斷層成像應(yīng)用基礎(chǔ)與臨床應(yīng)用的一些研究進展［J］.中國生物醫(yī)學工程學報，2010，29（2）：300-303.

［2］Townsend D W，Beyer T.A combined PET/CT scanner：the path to true image fusion［J］.The British Journal of Radiology，2002，75：24-30.

［3］孫濤，韓善清，汪家旺.PET/CT成像原理、優(yōu)勢及臨床應(yīng)用［J］.中國醫(yī)學物理學雜志，2010，7（1）：1581-1582.

［4］Lardinois D，Weder W，Hary，et al.Staging of non-small-cell lung cancer with integrated positron-emission tomography and comupted tomography［J］.N-Engl-J-Med，2003，348（25）：2500-2507.

［5］Bayford RH.Bioimpedance tomography（electrical impedance tomography）［J］.Annu Rev Biomed Eng，2006，8：63-91.

［6］董秀珍，劉銳崗，秦明新，等.電阻抗斷層成像軟件實驗系統(tǒng)設(shè)計［J］.第四軍醫(yī)大學學報，2000，21（11）：1381-1383.

［7］陳曉艷，王化祥，石小累，等.人體肺功能生物電阻抗成像技術(shù)［J］.中國生物醫(yī)學工程學報，2008，27（5）：663-668.

［8］Hahn G，Dittmar J，Just A，et al.Different approaches for quantifying ventilation distribution and lung tissue properties by functional EIT［C］//Proceedings of Electrical Impedance Tomography International Conference.Manchester：Institute of Physics Publishing，2009：73-84.

［9］許志影，李晉平.Matlab及其在圖像處理中的應(yīng)用［J］.江西：計算機與現(xiàn)代化，2003，92（4）：64-65.

［10］陳曉艷，王化祥.不同體位呼吸過程的電阻抗斷層成像研究［J］.生物醫(yī)學工程學雜志，2010，5（27）：1004-1007.

［11］Goodman TNT，Lee SL.Wavelets in wandering subspace［J］.Tran Amer Math Soc，1993，338（2）：639-654.

［12］李建平.小波分析與信號處理——理論、應(yīng)用及軟件實現(xiàn)［M］.重慶：重慶出版社，1997.

［13］Chibani Houacine A.Redundant versus orthogonal wavelet decomposition for multisensor image fusion［J］.Pattern Recognition，2003，36（4）：879-887.

［14］Garg Shruti，Ushah K，Mohan Ran，et al.Multilevel medical image fusion using segmented image by level set evolution with region competition［C］//Proceedings of the IEEE Engineering in Medicine and Biology 27thAnnual Conference.Shanghai：IEEE，2005：7680-7683.

［15］楊恒，裴繼紅，楊萬海.基于邊緣信息的多光譜高分辨圖像融合方法［J］.自動化學報，2002，28（3）：441-444.