陳曉艷 李健楠 王化祥（天津科技大學(xué)電子信息與自動(dòng)化學(xué)院，天津 300）

2（天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072）

引言

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

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

1 原理和方法

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

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

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

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

幅值反映該點(diǎn)處亮度值變化的大小，梯度的方向反映該點(diǎn)在哪個(gè)方向上變化最大。

由于數(shù)字圖像是離散的，所以其圖像亮度值變化的大小要用差分來表示［9］。二維圖像在x方向上的一階差分為f（x+1，y）-f（x，y），在y方向上的一階差分為f（x，y+1）-f（x，y），圖像差分時(shí)的［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ù)是通過圖像像素亮度值的一階或二階導(dǎo)數(shù)進(jìn)行判斷，當(dāng)幅值A(chǔ)大于某一閾值時(shí)則為邊緣。

由于canny算子具有較強(qiáng)的魯棒性，不容易受噪聲干擾，而且能有效地提取閉合區(qū)間，適用于確定阻抗敏感區(qū)域，所以采用它進(jìn)行邊緣提取。canny算子是先用高斯函數(shù)進(jìn)行濾波，然后計(jì)算沿著x，y兩個(gè)方向的偏導(dǎo)數(shù)，求出梯度的幅值和梯度的方向。通過梯度的方向就可以找到這個(gè)像素梯度方向的鄰接像素，然后將梯度的幅值和高低閾值進(jìn)行比較，大于高閾值的一定是邊緣，小于低閾值的一定不是邊緣，介于之間的要看這個(gè)像素的鄰接像素中有沒有梯度幅值超過高閾值的像素，如果有則為邊緣，否則就不是邊緣。筆者在用canny算子進(jìn)行邊緣提取時(shí)，高低閾值的選取是由程序自動(dòng)確定的。確定的原則就是：先通過設(shè)定非邊緣像素的比例來確定高閾值TH，非邊緣像素的比例一般選擇0.7；然后再確定低閾值TL，一般TL=0.4TH。先用disk算子進(jìn)行濾波，再用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ū)域左右肺臟及心臟的邊界，并將其導(dǎo)入Comsol平臺(tái)，如圖3所示。

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

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

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

搜索的方向?yàn)?/p>

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

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

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

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

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

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

從而有

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

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

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

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

1）對(duì)每一幅源圖像分別進(jìn)行小波分解；

2）對(duì)各分解層分別進(jìn)行融合處理，不同的分解層采用不同的融合算子；

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

由于小波變換可以把被融合的圖像分解到多個(gè)頻帶中，且具有方向性，故基于小波變換的圖像融合可以在不同頻帶、不同方向上采用不同的融合規(guī)則，以充分挖掘被融合圖像的互補(bǔ)信息，突出感興趣的特征和細(xì)節(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)小波算法進(jìn)行融合，融合效果如圖6所示。

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

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

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

3 討論和結(jié)論

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

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

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

［1］任超世，李章勇，王妍，等.電阻抗斷層成像應(yīng)用基礎(chǔ)與臨床應(yīng)用的一些研究進(jìn)展［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，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)勢(shì)及臨床應(yīng)用［J］.中國醫(yī)學(xué)物理學(xué)雜志，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］董秀珍，劉銳崗，秦明新，等.電阻抗斷層成像軟件實(shí)驗(yàn)系統(tǒng)設(shè)計(jì)［J］.第四軍醫(yī)大學(xué)學(xué)報(bào)，2000，21（11）：1381-1383.

［7］陳曉艷，王化祥，石小累，等.人體肺功能生物電阻抗成像技術(shù)［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，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］.江西：計(jì)算機(jī)與現(xiàn)代化，2003，92（4）：64-65.

［10］陳曉艷，王化祥.不同體位呼吸過程的電阻抗斷層成像研究［J］.生物醫(yī)學(xué)工程學(xué)雜志，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］李建平.小波分析與信號(hào)處理——理論、應(yīng)用及軟件實(shí)現(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］.自動(dòng)化學(xué)報(bào)，2002，28（3）：441-444.