电阻抗成像(Electrical Impedance Tomography，EIT)是一种描述生物组织电阻抗特性的医学成像技术。电阻抗成像二维重构算法是假设流入被测物体的电流只流经电流注入点的平面，通过静态^[1-2]或动态^[3-4]的重构方法得到电导率分布。因而降低算法复杂性，减少计算时间。电阻抗成像三维重构算法^[5-6]将注入被测物体的电流形成的电流场视为三维场实现重构。三维重构算法可以得到物体完整的电导率分布，但计算时间相对二维算法较长。临床EIT监护成像需要图像实时更新，三维重构算法会影响图像更新速度，从而影响监护质量。为满足EIT监护需求，笔者对二维重构图像进行三维重建。由于二维图像层间距大，重建三维的数据少，图像层间插值是必需的。Zimmermann^[7]在提出的电阻抗成像测量系统中通过提取等值面得到了圆柱容器中聚氯乙烯圆柱体和马铃薯模型，但具体细节并未阐述。

图像间插值主要分为基于灰度^[8]、基于形状^[9]和基于小波^[10]的插值方法。近年来的发展主要有：1)根据分割区域的不同密度对对应点匹配灰度插值的改进^[11]；2)形状插值前对初始轮廓的位置校准和插值后对轮廓的光滑处理^[12-13]；3)利用数学形态学对起始边界膨胀和腐蚀使其逐渐向目标边界逼近的形态插值^[14]；4)形状与灰度结合插值；5)基于小波的插值等。

匹配插值算法在插值效果上优于线性插值，而小波插值算法复杂，插值时间长，不利于实时成像。因此笔者采用匹配点灰度插值算法对二维重构图像插值，并根据电阻抗成像重构数据的特点对梯度和梯度方向角计算公式进行修改。插值得到图像间的新图像，并对插值后的体数据重建成三维图像。

1 匹配点插值算法原理

匹配点插值算法主要分为下列步骤：1)对二维电阻抗重构数据进行预处理，以便于插值使用；2)确定最佳匹配点对，笔者结合电阻抗重构图像数据的特点，对最佳匹配点算法中的梯度和梯度方向角计算公式作了改进；3)利用找到的最佳匹配点对待插值点线性插值。

1.1 插值前数据预处理

电阻抗成像重构二维图像时使用了有限元剖分，得到的重构数据是剖分后的三角形单元顶点坐标及三角形的电导率值。为便于插值，将共用一个顶点的所有三角形电导率的平均值赋给该顶点，遍历所有三角形单元得到所有顶点的电导率值。再将三角形单元网格化，判断网格点所在的三角形单元，并用该三角形顶点的电导率值插值得到网格点的电导率值。图像插值时，网格化后的网格点对应于传统匹配插值算法^[8]中的像素，网格点的电导率值对应于传统匹配插值算法^[8]中的灰度值。

1.2 确定最佳匹配点对

图 1中S_k和S_k+1是相邻两幅图像切片，即第k层图像和第k+1层图像。在S_k和S_k+1中分别选取宽度为W×W的窗口，其中a和a′，b和b′，c和c′是匹配点对。

一般最佳匹配点对的要求如文献[11]所述，由此构造匹配度计算函数如下：

$ \begin{array}{l} \;{R_{{\rm{m, n}}}} = {(v({x_{\rm{m}}}, {y_{\rm{n}}}, {z_k})- v(x{'_{\rm{m}}}, y{'_{\rm{n}}}, {z_{k + 1}}))^2} + \\ {({D_{{\rm{m, n, }}k}}({x_{\rm{m}}}, {y_{\rm{n}}}, {z_k})- {D_{{\rm{m, n}}, k + 1}}(x{'_{\rm{m}}}, y{'_{\rm{n}}}, {z_{k + 1}}))^2} + \\ {({\theta _{{\rm{m, n}}, k}}({x_{\rm{m}}}, {y_{\rm{n}}}, {z_{\rm{k}}})- {\theta _{{\rm{m, n}}, k + 1}}(x{'_{\rm{m}}}, y{'_{\rm{n}}}, {z_{k + 1}}))^2} + \\ {\left[\sqrt{{{{{({x_{\rm{m}}}-x{'_{\rm{m}}})}^2} + {{({y_{\rm{n}}}-y{'_{\rm{n}}})}^2}}}} \right]^2}. \end{array} $

(1)

上式计算匹配点对P_k(x_m, y_n, z_k)和P_k+1(x′_m, y′_n, z_k+1)间的匹配程度。P_k(x_m, y_n, z_k)是根据传统匹配算法^[8]中的匹配点选取公式获得，P_k+1(x′_m, y′_n, z_k+1)是依据文献[15]中计算匹配点位置公式得到。v(·)是网格点的电导率值，D_{m, n, k}是网格点的电导率梯度值，θ_{m, n, k}是网格点的电导率梯度方向，$\sqrt {{{\left( {{x_{\rm{m}}}-x{'_{\rm{m}}}} \right)}^2} + {{({y_{\rm{n}}}-y{'_{\rm{n}}})}^2}} $是匹配点对P_k(x_m, y_n, z_k)和P_k+1(x′_m, y′_n, z_k+1)在xy平面上的投影距离。

1.2.1 对梯度计算公式的改进

传统的灰度梯度计算公式^[11]如下：

$ \begin{array}{l} D\left( {x, y, z} \right) = \\ \sqrt {{{\left[{\frac{{v\left( {x + 1, y, z} \right)-v\left( {x-1, y, z} \right)}}{2}} \right]}^2} + } \\ \overline {{{\left[{\frac{{v\left( {x, y + 1, z} \right)-v\left( {x, y-1, z} \right)}}{2}} \right]}^2} + } \\ \overline {{{\left[{\frac{{v\left( {x, y, z + 1} \right)-v\left( {x, y, z-1} \right)}}{2}} \right]}^2}} 。\end{array} $

(2)

式(2)中v(·)为点的灰度值，对应于网格点的电导率值，该式计算x方向的梯度变化和y方向的梯度变化时，分别使用待计算点相邻两侧沿x方向和y方向的点。如图 1中计算b点梯度值时使用x方向的c点和d点，y方向的e点和f点。

如果待计算点相邻两侧中的一侧在x方向或者y方向上没有相邻点则视为边界点，当待插值点在上下两层对应的窗口中存在边界点时，该边界点的梯度将无法计算，同时也无法计算匹配函数值，该插值方法就不能插值此待插值点。电阻抗成像重构数据网格化后的网格点数据有限，如果继续使用式(2)则会减少真正用于插值的数据，不利于后续的三维重建。因此，再结合网格点x，y，z方向的间距将式(2)修改如下：

$ \begin{array}{l} D = \sqrt {{{\left[{\frac{{v\left( {x + 1, y, z} \right)-v\left( {x-1, y, z} \right)}}{{{d_{\rm{x}}}}}} \right]}^2} + } \\ \overline {{{\left[{\frac{{v\left( {x, y + 1, z} \right)-v\left( {x, y, z} \right)}}{{{d_{\rm{y}}}}}} \right]}^2} + } \\ \overline {{{\left[{\frac{{v\left( {x, y, z + 1} \right)-v\left( {x, y, z} \right)}}{{{d_{\rm{z}}}}}} \right]}^2}} 。\end{array} $

(3)

如图 1所示，d_x, d_y分别是两个邻近网格点间x和y方向的距离，d_z是图像层间距离。当待插值点对应的上下两层窗口中有边界点时，只利用边界点和边界内与其相邻的点计算梯度值。当窗口中的点落在右边界或下边界时，用该边界点以及其相邻左边或上边的点计算梯度值。如计算图 1中点a′的梯度值时，使用点m和点n。窗口中的点落在左边界或上边界也做同样的处理。当待插值点为边界点时，使用其在上下层中沿z轴的直接对应点线性插值得到。

1.2.2 梯度方向角计算公式的改进

θ_{m, n, k}是网格点的电导率梯度方向角，如图 1中的θ角，点a和点a′是一对匹配点。文献[15]使用反正切函数求得，并没有考虑z轴的影响，笔者是计算匹配点对的连线与z轴的夹角，公式如下：

$ \begin{array}{l} \;\;\;{\theta _{{\rm{m, n, }}k}}{\rm{ = }}\\ ar\cos \left( {\frac{{\left| {v({x_{\rm{m}}}, {y_{\rm{n}}}, {z_k})-v(x{'_{\rm{m}}}, y{'_{\rm{n}}}, {z_{k + 1}})} \right|}}{{{d_{\rm{z}}}{D_{{\rm{m, n}}, k}}}}} \right)。\end{array} $

(4)

其中D_{m, n, k}是式(3)得到的电导率梯度值，点(x_m, y_n, z_k)和点(x′_m, y′_n, z_k+1)是一对匹配点。当电导率梯度值为0时，电导率梯度方向角设为0。

得到每对匹配点的匹配度函数值后，选取匹配度函数值最小的一对对应点作为最佳匹配点对。

1.3 插值待插值点

找到最佳匹配点对后，利用线性插值对待插值点插值得到其电导率值，从而插值出整幅图像。

$ v({x_i}, {y_j}, {z_{k + {{\rm{d}}_1}}}) = \frac{{{d_2}}}{d}v({x_{\rm{m}}}, {y_{\rm{n}}}, {z_k}) + \frac{{{d_1}}}{d}v(x{'_{\rm{m}}}, y{'_{\rm{n}}}, {z_{k + 1}})。$

(5)

其中P_k(x_m, y_n, z_k)和P_k+1(x′_m, y′_n, z_k+1)为最佳匹配点对；d₁为待插值图像到上层图像的距离；d₂为待插值图像到下层图像的距离。

2 插值结果及算法验证 2.1 插值结果

实验开发环境为Visual C++6.0，实验模型是浸泡在圆柱体水槽中电导率为2.1 S/m的圆柱体异物，如图 2所示，圆柱形水槽的上表面是电极阵列，底部是背电极。图 3是圆柱体异物在圆柱形水槽中的位置和尺寸，圆柱体上表面距水槽上表面0.5 mm，圆柱体异物的半径为0.8 mm，高度为0.68 mm。沿圆柱体直径方向在每排电极阵列正下方获得8片切片，由一步牛顿法重构得到二维图像。插值后共有15片切片数据，即每相邻两层之间插值出一层新图像。

图 5(a)、(c)是图 4所示模型中的两相邻原始切片，(b)是对(a)和(c)应用笔者所述插值算法插值得到的新图像。其中黑色部分是原模型中的圆柱体异物，图像的灰度值对应网格点的电导率值。可以看出，从形状轮廓上新图像(b)呈现出(a)到(c)的过渡，基本无变形，而且也保持了电导率值。

2.2 插值算法验证

为验证匹配点灰度插值算法的有效性，取连续的3幅原始图像记为I₁、I₂、I₃。由两幅原始图像I₁和I₃插值出新的第2幅图像，并以原图I₂为标准和线性插值方法及文献[16]中提到的形状插值得到的结果进行比较。

笔者采用均方差和算法的运行时间^[17]来评判算法，均方差由$\sigma = \frac{1}{N}\sum {{{\left( {v\left( {x, y, z} \right)-v' \left( {x, y, z} \right)} \right)}^2}} $得到，其中N是图像中的点数，v(x, y, z)是插值后图像中点的电导率值，v′(x, y, z)是原图像中点的电导率值。算法的比较结果见表 1。

由表 1可知，笔者算法明显优于线性算法，而文献[16]中形状插值算法的均方差比笔者算法好，但耗时长。所以笔者用到的基于匹配点的插值算法兼顾插值精度和计算效率，符合电阻抗成像监护要求。

3 三维重建

对原始的8幅电阻抗图像插值，每两层插值出一层新图像，得到的插值后的15幅大小为81×86的电阻抗成像仿真重构切片作为三维重建的数据来源。在MATLAB R2008a环境下，调用isosurface(x, y, z, D, isovalue)函数^[18]提取等值面绘制得到, 其中x, y, z分别是体数据的x, y, z坐标值，大小为81×86×15的矩阵，D是体数据中点的电导率值，大小为81×86×15的矩阵，isovalue是根据电导率值设置的所要提取的等值面值，三维图像如图 6，7所示。

图 6是三维重建后得到的三维模型，由图 7可知重建结果很好地反映了原物体的位置特点，这对诊断病变位置有极大意义。但在模型形状上，原模型中是圆柱体，三维重建的结果却是稍微扁圆的物体。主要是以下原因造成的，在图 2的模型中，沿上表面电极注入电流，由背电极接收，由于上表面电极和背电极的位置关系，流经圆柱体异物的电流并非垂直于圆柱底面，而是倾斜的。二维重构的图像是与圆柱底面垂直的，忽略了这一倾斜角度，因而造成二维重构图像与原模型的圆柱体切面有差异，三维重建后的物体也与原模型的圆柱体异物有差异。

4 结语

笔者对电阻抗成像的二维重构断层图像序列进行图像间插值，并对插值后的体数据进行三维重建。结合了电阻抗二维重构计算时间短和三维重构可以得到体素电导率的优点，同时便于临床EIT监护。针对电阻抗重构数据的特点对插值算法的梯度和梯度方向角计算公式作了改进。插值出的新图像体现了图像层间缺失的信息，在电导率数值信息和形状上都有体现。并通过与线性插值和文献[17]中形状插值算法的对比可知，笔者算法可以兼顾插值精度和时间的要求，适合EIT图像插值。三维重建结果可以较好地反映原模型中圆柱体异物的位置特点。为提高精度，呈现更多的图像层间信息，笔者还可以在相邻图像间多插值几幅图像。