一种结合超声图像提高有限角度CT成像质量的方法

摘要

本发明公开了一种结合超声图像提高有限角度CT(Computed Tomography，计算机断层扫描)成像质量的方法，包括以下步骤：在目标组织上、下分别放置X射线发射器和接收器；在有限角度内使发射器绕组织移动并发射X射线，同时采集投影数据；利用同步代数重建技术进行初次重建；获取组织的三维超声图像；结合初次CT重建结果，对超声图像配准；提取超声图像的有效梯度信息；使用改进的同步代数重建技术进行再次重建，得到质量提高的有限角度CT成像结果。本发明利用超声成像在z方向边缘清晰的特性，解决了有限角度CT重建结果在z方向伪影严重、边界模糊的问题，具有创新性。

说明书

技术领域

本发明涉及生物医学及图像处理领域，特别是一种结合超声图像提高有限角 度CT(Computed Tomography，计算机断层扫描)成像质量的方法。

背景技术

传统的CT成像技术根据目标组织的不同部分对X射线的衰减系数不同这一 特性来重建物体内部的三维结构。该方法主要是利用锥束型X射线发射器环绕目 标组织360°，每隔一定角度发射一次X射线，并在每次发射时用X射线接收器在 目标组织的另一边接收，从而得到若干组投影数据。从这些投影数据重建得到目 标组织的三维图像的方法有很多，其中比较成熟的一种是同步代数重建技术 (Simultaneous Algebraic Reconstruction Technique，以下简写为SART)。该技术 可以有效重建出目标组织的内部构成，且具有分辨率高、边缘清晰、噪声小的优 点。然而在实际使用时，由于重建时间、X射线剂量等因素的制约，从360°对 目标物体进行全面扫描这一条件往往很难得到满足，取而代之的是在某一有限角 度范围内每隔一定角度进行一次X射线的发射与接收，也就是有限角度CT成像技 术。该技术使探测器平面和目标组织保持静止，在目标组织的正上方的一个有限 角度内使X射线发射器围绕目标组织弧形移动，从而采集到若干投影数据。然而 与全角度CT成像相比，由于数据量的大量缩减，有限角度CT成像结果在垂直于 接收器的方向上(以下称z方向)伪影严重，分辨率低，与实际被重建组织的差 别很大，大大影响了此技术在各领域的有效运用，亟需优化改善。

三维超声图像重建作为对X射线图像重建的一种补充和辅助，其结果在z方 向的分辨率高，边缘清晰，并且伪影比CT成像结果弱，恰恰可以弥补有限角度 CT成像在z方向的不足，有效提高有限角度CT成像结果在z方向的质量。通过 这种方法代入超声图像的边缘信息后，可以在保留X射线所成图像优点的前提 下，大大提高其图像质量，对于有限角度CT成像技术的实际运用具有重大意义。

发明内容

发明目的：本发明所要解决的技术问题是针对有限角度CT成像结果在z方 向上边缘不清晰，伪影严重的问题，提供一种结合超声图像提高有限角度CT成 像质量的万法。

为了解决上述技术问题，本发明公开了一种结合超声图像提高有限角度CT 成像质量的方法，包括以下步骤：

步骤一，在目标组织上、下分别放置X射线发射器和接收器；

步骤二，在有限角度内使发射器绕组织移动并发射X射线，同时采集投影数 据；

步骤三，利用SART算法进行初次重建；

步骤四，获取组织的三维超声图像；

步骤五，结合初次CT重建结果，对超声图像配准；

步骤六，提取超声图像的有效梯度信息；

步骤七，使用改进的SART算法进行再次重建，得到质量提高的有限角度 CT成像结果。

本发明中，优选地，所述步骤一中的X射线发射器发射锥束型X射线，目 标组织贴近X射线接收器的面的尺寸应小于接收器尺寸，使得通过物块的X射 线基本都被接收器所接收到。X射线发射器位于接收器的正上方，且与接收器达 到一定的垂直距离，使得X射线能够完全覆盖目标组织。

本发明中，优选地，所述步骤二中以X射线发射器在接收器平面上的投影点 为旋转中心，以发射器的初始位置为运动中心，过此旋转中心且平行于接收器某 一条边的直线为轴，X射线发射器在垂直于接收器的平面内绕此轴作范围有限的 弧形移动。该范围由弧形所对的圆心角表示，即旋转中心与发射器移动范围的左、 右边界点所确定的圆心角，其中左、右边界点关于运动中心对称。采集过程中， X射线发射器从运动中心开始，先向某一边界点移动，每移动一定的角度向接收 器发射X射线，接收器接收后视为采集到一组投影数据，到达该方向的边界点 后将发射器移回运动中心，再向另一边界点移动，经过同样的发射、接收过程， 发射器到达另一边界点后，采集过程结束，采集到若干组投影数据。

本发明中，优选地，所述步骤三中使用SART重建算法，对重建区域进行反 复的投影-反投影迭代，不断利用估计的投影值与步骤二中所测得的实际投影值 修正重建区域体素的值，最终使得结果满足最优化条件，获得初次有限角度CT 成像结果。

本发明中，优选地，所述步骤四中超声传感器分别在目标组织上表面与下表 面移动并向目标组织内部发射、接收超声波，从而获得上表面及下表面两组超声 成像数据。超声传感器的移动范围要覆盖整个目标组织区域，以便于有限角度 CT成像结果与超声成像结果的配准。

本发明中，优选地，所述步骤五中，由于三维超声成像的特点是沿着超声波 发射、接收方向(即z方向)进行成像，也就是所成的图像在一系列平行于xOz 或yOz的平面内(不失一般性，此处以三维超声图像重建在平行于xOz的平面内 为例)，且沿着y方向的空间分辨率很低，因此在进行配准时，首先根据有限角 度CT初次重建的结果，采用线性插值或剔除无效区域的方法，调整三维超声重 建结果中平行于xOz的平面内球的大小和位置，使得在这些平面内超声重建结果 中的肿块与相应平面内CT重建结果中的肿块大小、位置一致。然后通过线性插 值的方法扩充y方向上超声重建结果，增加三维超声结果在y方向的数据量，从 而使得平行于yOz的平面内的三维超声结果与CT成像结果的肿块位置、大小也 达到一致，最后再适当调整CT成像的成像区域，使得两种重建结果达到配准。

本发明中，优选地，所述步骤六选取超声成像结果在平行于xOz的平面内的 梯度信息，该梯度有x和z两个方向，对每个体素，其沿x方向的梯度值表示为 该体素的值减去其右方相邻体素的差值，沿z方向的梯度值表示为该体素的值减 去其上方相邻体素的差值。对于某些超声成像质量不高的图像，按照以上方法对 每张平行于xOz平面的图像求梯度后的梯度图形噪声大，肿块边缘模糊，边界线 粗，可以在不影响整体图像信息的前提下对超声成像结果使用中值滤波、全变分 优化等方法，也可以对计算得到的梯度图像进行中值滤波或者简单的形态学处 理，从而得到噪声小、边缘清晰的梯度图像。

本发明中，所述步骤七从损失函数的角度对传统的SART重建算法进行的扩 展与优化，传统的SART算法使用的损失函数为最小二乘函数其迭代公式所求解的就是使得L(x)最小的最优解 x。本发明中选取的损失函数加入步骤六中所获取的超声梯度信息，也就是利用 三维超声成像在x和z两个方向的梯度信息对有限角度CT成像在x和z两个方 向的梯度信息进行修正，其损失函数为 $L\left(x\right)=\frac{1}{2}{\left|\right|b-\mathrm{Ax}\left|\right|}_w^2+\frac{{\theta }_1}{2}\left|\right|G_{\mathrm{us}1}-G_{X1}{\left|\right|}^2+\frac{{\theta }_3}{2}\left|\right|G_{\mathrm{us}3}-G_{X3}{\left|\right|}^2,$通过求解这个最优化问 题，推导得到优化后的迭代公式。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上 述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法的流程图。

图2是一种有限角度CT成像投影采集过程的示意图。

图3是一种三维超声成像的示意图。

图4是求梯度的示意图。

图5是对目标组织体素的标号示意图。

具体实施方式：

本发明以传统的SART重建方法为基础，通过加入超声重建图像的有效梯度 信息，解决了传统有限角度CT成像在z方向图像质量不高，分辨率低的问题， 提高了重建图像的质量。

如图1所示，本发明公开了一种结合超声图像提高有限角度CT成像质量的 方法，包括以下步骤：

步骤一，在目标组织上、下分别放置X射线发射器和接收器；

步骤二，在有限角度内使发射器绕组织移动并发射X射线，同时采集投影数 据；

步骤三，利用SART算法进行初次重建；

步骤四，获取组织的三维超声图像；

步骤五，结合初次CT重建结果，对超声图像配准；

步骤六，提取超声图像的有效梯度信息；

步骤七，使用改进的SART算法进行再次重建，得到质量提高的有限角度 CT成像结果。

本发明中，步骤一，目标组织是一个长22.225cm，宽12.7cm，高5.08cm的长 方体物块，物块中包含若干大小不等的球形肿块以及菱形肿块；X射线接收器长 23.04cm，宽19.02cm，接收器由1152×951个接收传感器单元组成，每个单元为 一个边长为0.02cm的正方形，因此每次接收到的一组投影数据为一幅分辨率为 1152×951的图片。将该物块放置于X射线发射器与X射线接收器之间，贴近X 射线接收器平面，物块的底面与接收器平面之间有0.65cm的间隙。X射线发射 器位于物块正上方，与接收器的垂直距离为90cm，使得发射器发射的X射线能 基本覆盖整个目标组织。一种X射线发射器、目标组织、X射线接收器的放置 示意图如图2所示。

本发明中，步骤二，如图2所示，X射线发射器初始位置(运动中心)在点 P，P′为点P在接收器平面内的投影。点P位于目标组织的一个侧面OADG所 确定的平面内，以这个平面作为发射器的运动平面，以P′为圆心，使X射线发 射器围绕点P′作弧形移动，该弧所对的圆心角为60°，边界点L与边界点R到 运动中心所对的圆心角角度都为30°，即这两个边界点关于运动中心对称。X 射线从运动中心P开始，先向L点移动，每隔3°向X射线接收器发射一次X射 线，接收器接收到的数据即为这一角度上的一组投影值。当X射线发射器运动 点L后，将发射器移动回到运动中心P，并开始向另一边界点R移动，同样地， 每经过3°进行一次X射线的发射与接收，直到发射器到达点R，采集过程结束， 共采集到21组投影数据。

本发明中，步骤三，SART重建算法主要是通过一种迭代的方法求这样一类 方程的解x：

Ax=b       (1)

在公式(1)中，A为M×N矩阵，在重建过程中是指投影变换的系统矩阵， x是N×1列向量，由所有要重建的体素的值构成，即代表每个体素的衰减系数， N为重建体素的数目；b是M×1列向量，代表投影值，M是采集到的投影像 素的总数目。在有限角度CT成像的情况下，M的值一般远小于N，因此这个 方程为欠定方程，有无数多组解。SART算法主要的思想是投影和反投影，即先 对x给定一个初始估计值然后对这个估计值进行正投影得到投影列向量b的 一个估计值再用第一个角度下测量到的投影值b₁与之差作为修正值修 正即反投影过程，对所有21组投影值完成一次对的修正视为完成一次迭 代，最终收敛到使得||b-Ax||w²最小的x，即最小二乘解。根据SART理论， 加权矩阵W的引入可以有效地减小伪影。

从解决最优化问题的角度理解SART算法，也就是要求使得如下损失函数最 小的x的值，如公式(2)所示：

$L\left(x\right)=\frac{1}{2}{\left|\right|b-\mathrm{Ax}\left|\right|}_w^2---\left(2\right)$

利用梯度下降法解出使得L(x)最小的x即得到SART算法的迭代公式：

$x_j^{(t+1)}=x_j^{\left(t\right)}+\lambda \underset{i=1}{\overset{M}{\Sigma }}\frac{A_{\mathrm{ij}}}{A_{+,j}}·\frac{b_j-\left({\mathrm{Ax}}^{\left(t\right)}\right)}{A_{i,+}}---\left(3\right)$

其中，x_j^(t+1)表示第(t+1)次迭代后体素中第j点的值，i为所取的投影值的序号， A_+，j、A_i，+的值分别为：

$A_{+,j}=\underset{i=1}{\overset{M}{\Sigma }}\left|A_{\mathrm{ij}}\right|,A_{i,+}=\underset{j=1}{\overset{N}{\Sigma }}\left|A_{\mathrm{ij}}\right|(i=\mathrm{1,2},...,M;j=\mathrm{1,2},...,N)---\left(4\right)$

b_i为第i个投影点的值，(Ax^(t))_i为第i个估计的投影点的值，λ为松弛因子，影 响到收敛速度与重建结果。为了保证收敛，本例中，λ的值取为0.1，总迭代次 数取为3次。

如图5所示，通过使用SART算法，初次重建得到有限角度CT成像结果， 该重建结果的三维尺寸为m×n×l，本例中相应尺寸为1008×524×430，此重建 空间的每个体素在x和y方向的长度都为0.02cm，在z方向高度为0.01cm，即相 应的物理尺寸为20.16cm×10.48cm×4.3cm，其底面距离接收器平面1.57cm，该 尺寸略小于目标组织的实际尺寸，但是重建区域中能够包含物块中所有的肿块信 息。

本发明中，步骤四，三维超声成像示意图如图3所示，其中的箭头表示超声 传感器发射、接收超声波的过程，通过超声传感器分别在目标组织上、下两个表 面的平移、扫描过程，得到目标组织的两组三维超声重建结果，这样的两组结果 在各自靠近传感器的方向重建分辨率高，肿块清晰。由于这两组结果有一部分区 域会发生重叠，所以可以通过一系列的标记与融合处理，将这两组结果进行融合 从而得到完整的三维超声重建结果。

三维超声重建得到的结果在平行于xOz的平面内分辨率高，边缘清晰，但是 沿着y方向分辨率低。在本例中，超声成像结果的尺寸为192×489×498，其中 192为y方向尺寸，489为x方向尺寸，498为z方向尺寸。

本发明中，步骤五，由于有限角度CT成像与超声成像是两种不同的成像模 式，所以得到的结果必然会存在一定的差别，尤其是在二者的尺寸方面。然而由 于这两种方法又是对同一种目标组织的重建，那么这两种结果必定都与原先的目 标组织有一定的线性比例关系，相应的，这两种结果之间也必然存在着线性比例 关系，这就为配准工作的进行提供了理论依据。

配准的核心目的是调准两种重建结果中对应肿块的大小和位置。参照步骤三 与步骤四可知，本例中，有限角度CT重建结果的区域尺寸为1008×524×430， 而三维超声成像结果的尺寸为192×489×498，由此可见，在CT结果与超声结果 沿y方向物理长度差不多的情况下，超声重建结果在y方向的层数比CT图像少 很多，因此超声结果在y方向分辨率不高，需要对其进行插值，使得层数得到扩 充。配准工作的第一步是完成与xOz平面相平行的平面的配准，由于CT重建过 程中，每个体素的高度是体素在x方向和y方向长度的一半，因此首先要将超声 重建结果在z方向进行线性插值，使其高度变为原先的两倍；然后通过有限角度 CT与超声图像中相对应肿块的大小和位置，确定两种重建结果的一个线性比例， 对超声成像结果每一平行于xOz的平面按照这个比例线性拉伸。由于超声成像对 整个物块都有进行重建，因此重建区域大于X射线重建区域，所以需要对拉伸 后的结果剔除四周的无效区域(即CT成像没有重建的区域)，从而得到尺寸为 192×524×430的超声结果。第二步是对三维超声结果在y方向的进行扩充，根 据有限角度CT成像中y方向上肿块的大小和位置，确定超声结果在y方向的拉 伸比例，按照这个比例对超声结果在y方向进行拉伸；同样地，拉伸结束后适当 剔除两边的无效区域，获得完全配准好的尺寸为1008×524×430的三维超声图 像。

本发明中，步骤六，超声成像在与xOz平面平行的平面内边缘清晰(参照图 5)，这也是超声成像的优势所在，因此取这些平面内的肿块梯度信息作为优化信 息加入改进后的SART算法当中。选取梯度信息时，在不影响超声成像中肿块主 要信息的前提下对超声成像结果进行一些图像处理(如中值滤波、全变分优化处 理)，对计算得到的边缘图像进行了中值滤波和形态学处理，从而使得所提取的 边缘清晰、平滑。

计算超声重建图像在x方向与z方向的梯度时，其计算方法要与CT重建结 果中计算梯度的方法一致。本例中，梯度的计算方法为：某一体素点的梯度值表 示为这一点的体素值减去其右边(x方向)或上边(z方向)点的体素值，如图 3所示，点f(x₀，z₀)在x方向的梯度值表示为：

G_x(x₀，z₀)=f(x₀，z₀)-f(x₀+1，z₀)    (5)

类似的，点f(x₀，z₀)在z方向的梯度值可以表示为：

G_z(x₀，z₀)=f(x₀，z₀)-f(x₀，z₀+1)    (6)

公式(5)、(6)中，G表示单个体素值的梯度信息，是一个标量，下标x与 z表示x方向与z方向的梯度值。

对于每张平行于xOz平面的图片的最右面一列像素和最上面一行像素，需要 将其梯度值置为零。

本发明中，步骤七，对传统的SART算法进行了优化，在其损失函数(2) 的基础上加上了包含超声成像结果在x方向和z方向信息的损失项，如公式(7)：

$L\left(x\right)=\frac{1}{2}{\left|\right|b-\mathrm{Ax}\left|\right|}_w^2+\frac{{\theta }_1}{2}\left|\right|G_{\mathrm{us}1}-G_{X1}{\left|\right|}^2+\frac{{\theta }_3}{2}\left|\right|G_{\mathrm{us}3}-G_{X3}{\left|\right|}^2---\left(7\right)$

公式(7)中，θ₁、θ₃用来控制梯度信息对重建结果的影响大小。G表示梯度， 加粗表示为向量。为避免混淆，下标中的1、3分别表示x方向和z方向，下标us 表示超声成像，下标X表示有限角度CT成像。G_us1、G_us3即步骤六中所计算得 到的超声重建结果中所有点在x方向和z方向的梯度值，G_x1、G_x3并未在事先 计算得到，而是在迭代过程中计算得到，可以表示为公式(8)、(9)：

G_x1=Bx      (8)

G_x3=Cx      (9)

其中，x是待重建的未知体素的体素值，是N维列向量，N即所要重建的 体素数。对所重建的目标组织中的体素按照以下规则进行标号：以原点处的体素 为第一个体素x₁，从原点(x=0，y=0，z=0)开始，先沿y方向增加体素的序号， 达到末尾时转到下一行(x=1)继续计数，以此类推，计完一个xOy平面后移 至下一个平面(z=1)继续计数，这样计完所有体素点，从而将目标组织中每个 体素的序号从三维坐标转化为一维坐标，如图5所示，其中m表示重建区域在y 方向的尺寸，n表示重建区域在x方向的尺寸，l表示重建区域在z方向的尺寸。

公式(8)、(9)中的矩阵B、C用来计算x方向和z方向的梯度，与步骤六 中求梯度的方法一样，此处体素的梯度值也表示为该体素与其右边(x方向)或上 边(z方向)相邻体素值之差，因此，矩阵B、C可分别表示为：

$B=\left(\begin{array}{c}{B}_1\\ B_2\\ ·\\ ·\\ ·\\ B_j\\ ·\\ ·\\ ·\\ B_N\end{array}\right)---\left(10\right)$

$C=\left(\begin{array}{c}{C}_1\\ C_2\\ ·\\ ·\\ ·\\ C_j\\ ·\\ ·\\ ·\\ C_N\end{array}\right)---\left(11\right)$

其中，B_j(j=1，2，...，N)，C_j(j=1，2，...，N)都是N维行向量，分别满足以下 两个条件：

1、B_j：当j∈平面CBEF(如图5)时，B_j=0；j平面CBEF时，B_j的第j 个元素为1，第j+(m+1)个元素为-1，其余元素都为0，即使得B_jx=x_j-x_j+(m+1)， 也就是求得体素x_j在x方向的梯度。

2、C_j：当j∈平面DEFG时(如图5)时，C_j=0；j平面DEFG时，C_j的 第j个元素为1，第j+(m+1)(n+1)个元素为-1，其余元素都为0，即使得 C_jx=x_j-x_j+(m+1)(n+1)，也就是求得体素x_j在z方向的梯度。

将公式(8)、(9)代入公式(7)，得到如下损失函数：

$L\left(x\right)=\frac{1}{2}{\left|\right|b-\mathrm{Ax}\left|\right|}_w^2+\frac{{\theta }_1}{2}\left|\right|G_{\mathrm{us}1}-\mathrm{Bx}{\left|\right|}^2+\frac{{\theta }_3}{2}{\left|\right|G_{\mathrm{us}3}-\mathrm{Cx}\left|\right|}^2---\left(12\right)$

利用梯度下降法求在满足Ax=b的条件下使L(x)最小的解x，对公式(12) 求梯度得到：

根据梯度下降法可知，解x可由以下公式迭代得到：

根据SART算法，ω取为λ，并加入系数V^-1，其中，V是一个对角矩阵：

加权矩阵W也是对角矩阵：

W与V中的A_i，+A_+，j如公式(4)所示。

本发明中，将ωθ₁、ωθ₃设定为常数λ₁、λ₃，为了保证收敛，本例中将这 两个参数设定为0.03，在其它实施条件中还应根据具体的梯度值大小对λ₁、λ₃进 行一些测试与调整。因此，可以得到如下公式：

x^(t+1)=x^(t)+λV^-1A^TW(b-Ax^(t))+λ₁B^T(G_us1-Bx^(t))+λ₃C^T(G_us3-Cx^(t))(17)

将公式(17)展开，便得到对于每个体素点的迭代公式：

$x_j^{(t+1)}=x_j^{\left(t\right)}+\lambda \underset{i=1}{\overset{M}{\Sigma }}\frac{A_{\mathrm{ij}}}{A_{+,j}}·\frac{b_i-\left({\mathrm{Ax}}^{\left(t\right)}\right)}{A_{i,+}}+{\lambda }_1{P}_j\left(x_j^{\left(t\right)}\right)+{\lambda }_3{Q}_j\left(x_j^{\left(t\right)}\right)---\left(18\right)$

其中，如式(19)、(20)所示：

其中，Gu_s1，j表示超声成像中第j体素点沿x方向的梯度值，G_us3，j表示超声成 像中第j体素点沿z方向的梯度值。G_X1，j表示第t次迭代后重建结果的第j体素 点沿x方向的梯度值，G_X3，j表示第t次迭代后重建结果的第j体素点沿z方向的 梯度值，可以表示为：

$G_{X1,j}=x_j^{\left(t\right)}-x_{j+(m+1)}^{\left(t\right)}---\left(21\right)$

$G_{X3,j}=x_j^{\left(t\right)}-x_{j+(m+1)(n+1)}^{\left(t\right)}---\left(22\right)$

将梯度信息代入到改进后的SART算法当中，再次进行重建就可以获得z方 向图像质量提高的有限角度CT重建结果。需要指出的是，再次重建的过程中， 由于有限角度CT重建需要进行3次迭代，每次迭代过程中对每个体素点进行21 次更新，也就是对每个体素点一共进行63次更新，这不能够使得加入的梯度信 息收敛。要使梯度信息收敛，需要进行约1000次的迭代。因此，在本发明中将 梯度信息的迭代过程与SART迭代过程分离，即在对每一组投影信息进行SART 迭代以后，以这个迭代结果作为初始值再进行15次梯度信息的迭代，这样作为 对一组投影数据的迭代过程，对所有21组投影数据完成计算作为一次完整的迭 代过程结束。这样进行3次迭代后，每一个投影点进行了63×15=945次梯度信 息的迭代，可以保证重建结果的收敛。

本发明提供了一种结合超声图像提高有限角度CT成像质量的方法，应当指 出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可 以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。另外，本 实施例中未明确的各组成部分均可用现有技术加以实现。