一种电阻抗图像与CT图像融合方法研究

陈晓艳 李健楠 王化祥（天津科技大学电子信息与自动化学院，天津 300）

2（天津大学电气与自动化工程学院，天津 300072）

引言

医学图像主要包括形态图像、结构图像和功能图像。医学显微图像属于形态图像；X线图像、超声、CT、MRI等属于结构图像；PET、SPECT、FMRT及EIT等为功能图像［1-2］。由于成像原理不同，不同图像提供的人体生理、病理信息也不同。将不同原理的图像进行融合，获得信息更为丰富的复合图像，是当前临床医学图像的研究热点之一［3-4］。

电阻抗成像技术（electricalimpedance tomography，EIT）是近30年发展起来的新兴成像技术，其无创、无辐射、成本低廉，能长时间连续工作，更重要的是EIT能反映人体组织、器官的功能变化，在医学图像监护方面具有显著的优势［5-6］。笔者在电阻抗成像技术基础上，以人体呼吸过程的电阻抗图像及胸部CT图像融合为例，对EIT图像和CT图像的融合方法进行研究。在实现功能成像的同时，对阻抗变化敏感区域进行有效定位［7-8］。

1 原理和方法

1.1 提取结构信息

以人体胸部CT图像为例展开研究。让被测者躺在16螺旋CT成像设备测试床上，在吸气末时刻获取CT图像，如图1所示。其中，胸腔体表外周的圆亮点是具有透射性的EIT电极，下面是测试床截面。由图可知，CT图像不仅能提供电阻抗成像所需要的胸腔外部边界结构信息，而且也能提供丰富的内部结构信息。具体提取方法是：先采用圆形区域均值滤波disk算子进行滤波，然后再采用canny边缘提取算子提取图像轮廓。

图1 人体胸部CT图像Fig.1 The chest CT image

设图像上像素点（x，y）处的梯度f、幅值A及梯度方向φ分别为

幅值反映该点处亮度值变化的大小，梯度的方向反映该点在哪个方向上变化最大。

由于数字图像是离散的，所以其图像亮度值变化的大小要用差分来表示［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）

边缘判定的依据是通过图像像素亮度值的一阶或二阶导数进行判断，当幅值A大于某一阈值时则为边缘。

由于canny算子具有较强的鲁棒性，不容易受噪声干扰，而且能有效地提取闭合区间，适用于确定阻抗敏感区域，所以采用它进行边缘提取。canny算子是先用高斯函数进行滤波，然后计算沿着x，y两个方向的偏导数，求出梯度的幅值和梯度的方向。通过梯度的方向就可以找到这个像素梯度方向的邻接像素，然后将梯度的幅值和高低阈值进行比较，大于高阈值的一定是边缘，小于低阈值的一定不是边缘，介于之间的要看这个像素的邻接像素中有没有梯度幅值超过高阈值的像素，如果有则为边缘，否则就不是边缘。笔者在用canny算子进行边缘提取时，高低阈值的选取是由程序自动确定的。确定的原则就是：先通过设定非边缘像素的比例来确定高阈值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中提取的胸腔及肺部的边界结构，勾勒出感兴趣区域左右肺脏及心脏的边界，并将其导入Comsol平台，如图3所示。

先在Comsol软件平台上进行正问题求解，得到灵敏度系数矩阵，再将其导入Matlab，经共轭梯度算法重建EIT图像。共轭梯度算法是介于最速下降法与牛顿法之间的一个方法，仅需利用一阶导数信息，即克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hessian矩阵并求逆的缺点。共轭梯度算法经过第i次迭代得到的解xi连续逼近真解，迭代误差ri=b-Axi，搜索方向设定为Pi，迭代次数与迭代误差的更新为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）。当迭代误差小于设定阈值时，迭代结束。对于EIT成像，系数矩阵（Jacobian矩阵）往往不是对称正定矩阵，因此需要对原方程进行转换，有

由于A并非列满秩矩阵，因此ATA也不是对称正定矩阵，可以通过增加正则化矩阵的办法来解决这个问题。当正则化矩阵为单位矩阵时，式（7）可以转化为

此时，系数矩阵为对称正定矩阵，因此EIT可用上述共轭梯度算法进行求解［10］。重建后的图像如图4所示。

1.3 基于小波变换的图像融合

由于人的视网膜是在不同的频率信道中进行不同的信号处理［11］，为了获得EIT图像和CT图像良好的融合效果，需要对EIT图像和CT图像在不同的频率信道上进行融合，而小波变换可实现将原始图像分解成一系列具有不同频域特性的子图像，从而为进行不同频率的图像数据融合提供了有效方法。

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

从而有

这就是说，ψ（t）与整个横坐标所围成面积的代数和为0，因此ψ（t）的图形应是在横坐标上下波动的“小波”。通常称ψ（t）为母小波，则称ψa，b（t）为小波函数，或简称为小波（wavelet）［12］。

令L2（R）为可测的、平方可积函数f（t）的矢量空间，R为实数集，对于任意的f（t）∈L2（R）的连续小波变换给出如下定义，即：

小波变换的重构（逆变换）公式为

设EIT图像和CT图像分别为A和B，F为融合后的图像，融合的基本步骤为：

1）对每一幅源图像分别进行小波分解；

2）对各分解层分别进行融合处理，不同的分解层采用不同的融合算子；

3）对融合后的数据进行小波逆变换，所得到的重构图像即为融合图像。

由于小波变换可以把被融合的图像分解到多个频带中，且具有方向性，故基于小波变换的图像融合可以在不同频带、不同方向上采用不同的融合规则，以充分挖掘被融合图像的互补信息，突出感兴趣的特征和细节信息，从而获得良好的视觉效果［13-15］。基于小波变换的图像融合流程如图5所示。

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

2 结果

将图2中胸腔边界之外无用信息滤除后的图像，与EIT重建图像（见图4）经小波算法进行融合，融合效果如图6所示。

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

通常，在电阻抗成像过程中，灵敏度矩阵的求解是在敏感区域内部各点电导率相同的假设前提下获得的。但在生物医学电阻抗成像的应用研究中，人体脑部、腹部、肺部等区域内的电导率各不相同，而其分布是十分重要的先验信息。

本研究充分利用结构和电导率的先验信息，获得更准确的灵敏度矩阵，从而达到提高电阻抗成像质量的目的。在此基础上，采用小波算法，将CT结构图像与EIT功能图像进一步融合。图6与图像融合前的图4相比，EIT图像中的电阻抗功能变化信息在结构图像中凸显出来，可以实现电阻抗变化区域的初步定位。

3 讨论和结论

电阻抗功能成像技术的优势显而易见，但是如何充分发挥其优势，服务临床、辅助诊疗是目前亟待解决的问题，笔者就是围绕如何实现病灶的定位来开展探索性的研究。借助高精度的CT结构图像，不仅可以提高电阻抗图像质量，更有望融合后的图像能提供定位性的生理和病理变化信息。

本研究探索了EIT图像与CT图像融合的有效方法，并获得了良好的融合效果，为EIT成像技术与其他技术（如超声、MRI、X线等，甚至其他功能成像技术）进行图像融合奠定了基础。但是，本研究仅限于两种静态图像的融合，对于呼吸等动态变化的过程，融合技术和方法会有所不同，有待进一步深入研究。

笔者主要研究了CT图像轮廓信息的提取方法，通过对各种滤波算子以及边缘提取算子的综合比较，表明用disk算子进行滤波，然后再用canny算子进行边缘提取，所得到的图像比较符合要求。提取的CT轮廓信息可以建立EIT重建模型，也可以与EIT图像进行融合。本研究验证了EIT图像与CT图像融合的可行性，探索了一种有效的图像融合方法，复合图像不仅可以把成像区域的结构信息与功能信息同时显示出来，实现信息互补，而且为进一步研究功能变化信息的定位问题开辟了新的途径。

［1］任超世，李章勇，王妍，等.电阻抗断层成像应用基础与临床应用的一些研究进展［J］.中国生物医学工程学报，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成像原理、优势及临床应用［J］.中国医学物理学杂志，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］董秀珍，刘锐岗，秦明新，等.电阻抗断层成像软件实验系统设计［J］.第四军医大学学报，2000，21（11）：1381-1383.

［7］陈晓艳，王化祥，石小累，等.人体肺功能生物电阻抗成像技术［J］.中国生物医学工程学报，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及其在图像处理中的应用［J］.江西：计算机与现代化，2003，92（4）：64-65.

［10］陈晓艳，王化祥.不同体位呼吸过程的电阻抗断层成像研究［J］.生物医学工程学杂志，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］李建平.小波分析与信号处理——理论、应用及软件实现［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.