CN111968192A - 一种ct图像的构建方法、ct设备以及存储介质 - Google Patents

Abstract

本申请提供一种CT图像构建方法、CT设备及计算机可读存储介质。该方法包括：获取CT设备进行扫描时采集的投影数据；利用投影数据构建图像重建模型，其中，图像重建模型是基于像素邻域块的全变分模型，像素邻域块为以一个像素为中心的像素区域；将预设图像数据输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像。本申请提供的CT图像构建方法通过构建基于像素邻域块的全变分模型，并利用迭代算法求解该重建模型，使得到的重建CT图像分辨率更高，更加清晰。

Description

一种CT图像的构建方法、CT设备以及存储介质

技术领域

本申请涉及医学影像的计算机断层成像技术领域，特别是涉及一种CT图像的构建方法、设备以及存储介质。

背景技术

X射线CT扫描已经广泛应用于临床医学影像诊断，但是CT扫描过程中过高的X射线辐射剂量会存在致癌风险。为了降低对使用者的损害,如何最大限度地降低X射线使用剂量已经成为医学CT成像领域研究的关键技术之一。

为了降低X射线辐射剂量，一般地，使用的最简便的途径就是降低CT扫描过程中的管电流和扫描时间。然而，由于降低了管电流和扫描时间，就会使得投影数据中含有大量的噪声，基于传统的滤波反投影方法重建的图像质量存在严重的退化现象，难以满足临床诊断需要。

为了在保证图像质量的前提下大幅降低X射线辐射剂量，诸多基于降低管电流和扫描时间的低剂量CT图像重建方法相继提出，例如基于统计模型的迭代重建方法和基于投影数据滤波的解析重建方法。其中，基于统计模型的迭代重建方法，通过对采集的投影数据的噪声以及成像系统进行图像重建模型构建，可以实现低剂量CT图像优质重建；基于投影数据滤波的解析重建方法，是通过对采集的投影数据的噪声以及成像系统进行数据滤波建模，再通过解析重建方法实现快速且优质的低剂量CT图像重建。

发明内容

本申请主要提供一种CT图像的构建方法，以提高CT图像重建质量。

为解决上述技术问题，本申请采用的一个技术方案是：提供一种CT图像的构建方法。该方法包括：获取CT设备进行扫描时采集的投影数据；利用投影数据构建图像重建模型，其中，图像重建模型是基于像素邻域块的全变分模型，像素邻域块为以一个像素为中心的像素区域；将预设图像数据输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像。

其中，将预设图像数据输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像，包括：利用解析重建算法对投影数据进行重建，获得初始CT图像，并将初始CT图像作为预设图像数据；将初始CT图像输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像。

其中，利用解析重建算法对投影数据进行重建，获得初始CT图像，包括：利用滤波反投影算法对投影数据进行重建，获得初始CT图像。

其中，基于像素邻域块的全变分模型为：

其中，f为投影数据，μ为待重建的CT图像，H＝{H_i,j}为系统矩阵，i,j分别为待重建CT图像的像素个数和CT设备探测器的探测元个数，D(f,Hμ)为数据保真项，|μ|_TV表示全变分正则化项，β表示平衡保真度和正则化项的超参数。

其中，全变分正则化项为：

其中s和t表示CT图像中衰减系数位置的指标，δ为一常数，用以保持与图像强度的可微性，μ_s,t,l表示以位于重建CT图像的第s行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s-1,t,l表示以位于重建CT图像的第s-1行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s,t-1,l表示以位于重建CT图像的第s行以及第t-1列的像素，为中心的像素邻域块的第l个像素的像素值，N_s,t表示像素邻域块包含的像素总数。

其中，全变分模型为惩罚加权最小二乘全变分算法，其中，惩罚加权最小二乘全变分算法引入基于像素邻域块的全变分正则项为惩罚项。

其中，基于像素邻域块的全变分模型为：

其中，Σ表示对角矩阵，且

表示探测器信道i处投影数据的方差。

其中，探测器信道i处投影数据的方差根据以下公式获取：

其中，I₀表示X射线入射光子强度，

表示系统电子噪声的方差，

表示探测器信道i处投影数据的均值。

其中，系统矩阵是从CT设备系统得到的基于面积加权的系统矩阵，或者基于体素加权的系统矩阵。

其中，迭代运算为梯度下降算法、共轭梯度下降算法或超松弛迭代算法中的任意一种。

为解决上述技术问题，本申请采用的又一个技术方案是，提供一种CT设备，该CT设备包括：内部总线，以及通过内部总线连接的存储器和处理器；存储器用于存储计算机程序；处理器用于执行计算机程序以实现上述CT图像构建方法的步骤。

为解决上述技术问题，本申请采用的另一个技术方案是，提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述CT图像构建方法。

本申请的有益效果是：区别于现有技术的情况，本申请提供的CT图像构建方法通过构建基于像素邻域块的全变分模型，其中，像素邻域块为以一个像素为中心的像素区域。也即，该CT图像重建模型引入基于像素邻域块的全变分正则项，且该全变分正则项不同于传统的全变分正则项，本实施例提供的全变分正则项使用与像素相关联的像素邻域块计算二维图像的该像素的梯度，能够消除阶梯效应，保留图像细节，提高重建CT图像的分辨率。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，其中：

图1是本申请提供的CT图像构建方法一实施例的流程示意图；

图2是图1中步骤S30一实施例的流程示意；

图3(a)至图3(e)描述了本申请方法与其他方法的重建效果对比；

图4为本申请提供的CT设备一实施例的结构示意图；

图5为本申请提供的存储介质一实施例的结构示意图。

具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本申请的一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

本申请实施例中的术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”、“第三”的特征可以明示或者隐含地包括至少一个该特征。本申请的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤或单元。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。

CT是用X线束对人体的某一部分一定厚度的层面进行扫描，由探测器接收透过该层面的X线，所测得的信号经过模数转换，转变为数字信息后由计算机进行处理，从而得到该层面的各个单位容积的X线吸收值即CT值，并排列成数字矩阵。这些数据信息科存储在磁光盘或磁带机中，经过数模转换后再形成模拟信号，经过计算机的一定变换后输出至显示设备上显示出图像，其密度分辨率高，可直接显示X线平面无法显示的器官和病变，它在发现病变、确定病变的位置、大小、数目方面非常敏感而可靠。CT成像本质上是人体组织的衰减系数成像，成像物理原理为通过CT扫描机构扫描获取求解衰减系数的方程组。解方程组获得人体某一体层面各个体素的衰减系数值，再将衰减系数值转换为CT值，最后将CT值变换为能视觉识别的灰度图像。CT图像重建是运用一定的物理技术测定X线在人体内的吸收系数为基础，采用一定的数学方法经计算机求解出吸收系数值在人体某剖面上的二维分布矩阵，再用电子技术把吸收系数二维分布矩阵转变诶图像面上的灰度分布，实现重建体层图像目的。每个小单元体按照扫描过程中的顺序进行排列和编号，形成一个有序得到数组，这些有序的数组在图像平面上形成图像矩阵。目前，CT图像重建算法主要包括解析法和迭代法。解析法以卷积反投影算法最为常用，解析法的优点是重建速度快，缺点是抗噪声性能较差，但是对数据的完备性要求较高。迭代法的基本思想是由测量的投影数据建立一组未知向量的代数方程式，通过方程组纠结未知图像向量。统计迭代重建算法以优化理论为基础，从投影测量过程的随机性观点出发，把图像重建看成是一个参数估计问题，通过设计合理的目标函数寻求使目标函数达到最优值的参数向量。统计迭代重建中，将提前获取的与图像相关的先验信息引入低剂量CT图像重建的目标函数中，作为正则化项，使得求解过程更加稳定，更好地恢复组织结构信息以及抑制噪声。对不完全数据也可重建图像，经典的统计迭代算法有最大似然估计算法、最小二乘算法、基于Bayesian的最大后验算法等等。

线性代数中，将不适定问题转化为适定问题称为正则化。因此可以通过引入待建图像的先验知识来作为正则化项，使高度病态的不完备投影数据获得稳定而准确的重建。其中最常见的是全变分(Total Variation，TV)正则化迭代方法。2006年已经证明一个信号如果是可稀疏表示的，则可以利用全变分最小化作为正则化项，从少量测量数据中精确重建该信号。全变分最小化的稀疏角度CT图像重建方法每次迭代都由凸集投影和梯度下降两项组成。在此基础上，很多文献提出了一些改进的全变分正则化迭代CT图像重建方法。但是当投影数据含噪声时，利用上述方法得到的重建图像中出现了较为严重的噪声干扰和条形伪影，且部分细节被噪声掩盖，不能得到满意的重建图像。本申请发明人经过长期研究发现，传统的非二次惩罚项不稳定的主要原因是图像的粗糙度是基于相邻像素之间的强度差来计算的。当图像有噪声时，像素强度的差异在区分真实边缘和噪声波动时是不可靠的。而现有的非局部正则化方法要么需要预先知道的参考图像来构造权重函数，要么涉及到非凸优化。在CT图像重建中，好的参考图像在重建前是不能直接得到的。由于μ的存在，非凸优化问题常导致图像估计不稳定。

为了克服这一问题，本申请发明人提出了一种基于像素邻域块的全变分重建方法，该方法在测量图像粗糙度时使用像素邻域块而不是单个像素。由于它比较了不同像素邻域块之间的相似性，因此基于像素邻域块的全变分重建比基于像素的全变分重建方法更具有鲁棒性。

参阅图1，图1是本申请提供的CT图像构建方法一实施例的流程示意图。该方法包括以下步骤：

S10：获取CT设备进行扫描时采集的投影数据。

利用X线CT成像设备在地毫安秒扫描协议下采集病人低剂量射线的CT医学图像的线积分投影数据，并同时获取相应的校正参数以及系统矩阵，其中，低剂量射线的射线剂量为标准剂量射线的1/7至1/20，可以由用户灵活设定，其中线积分投影数据是指对数变换后的投影数据，获取的校正参数是指X线入射光子强度I₀、系统电子噪声的方差

等，通常这些数值可以直接从测试设备读取或得到。利用siddon算法计算待测体中第j个组织块相对于第i个投影数据f_i的贡献，记为h_ij，从而得到系统矩阵H。本实施例中的投影数据可以是CT设备对病人头部扫描得到投影数据、对病人全身扫描得到的投影数据或者是对病人脊椎骨等结构细节扫描得到的投影数据，在此不做具体限制。

可选地，对扫描的投影数据进行预处理和取负对数得到处理后的投影数据。包括对投影数据进行噪声处理、扫描数据标准化处理等。

S20：利用投影数据构建图像重建模型，其中，图像重建模型是基于像素邻域块的全变分模型，像素邻域块为以一个像素为中心的像素区域。可选地，像素邻域块为以一个像素为中心，以设定像素单位为像素半径的正方形像素区域。

具体地，本申请构建的基于像素邻域块的全变分模型如下：

其中，f为投影数据，μ为待重建的CT图像，H＝{H_i,j}为系统矩阵，i,j分别为待重建CT图像的像素个数和CT探测器探测元个数，由具体的CT成像系统决定，可以采用不同的方法计算，比如，基于面积加权的系统矩阵计算，或者基于体素加权的系统矩阵计算。

D(f,Hμ)为数据保真项，|μ|_TV表示全变分正则化项，β表示平衡保真度和正则化项的超参数。经本申请发明人反复试验，发现本申请提供的图像构建方法在平衡保真度和正则化项的超参数取β＝0.02时能够取得很好的重建效果。

更具体地，全变分正则化项为：

其中s和t表示重建CT图像中衰减系数位置的指标，δ为一常数，用以保持与图像强度的可微性，μ_s,t,l表示以位于重建CT图像的第s行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s-1,t,l表示以位于重建CT图像的第s-1行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s,t-1,l表示以位于重建CT图像的第s行以及第t-1列的像素，为中心的像素邻域块的第l个像素的像素值，N_s,t表示像素邻域块包含的像素总数。

全变分算子能使重建后的图像有较高的精度，具有各项异性的平滑作用，通过梯度的计算不但能够保护图像的边缘，而且可以较好地保持图像边界原有的对比度锐度。并且受噪声污染的图形的全变分要明显比没有经过噪声污染的图像的全变分大，通过限制全变分的大小就可以限制图像的噪声。但是，其在平滑区域的处理效果相对较差，容易出现阶梯效应，并且全变分在去除噪声的过程中一些细节如问题信息容易别过滤掉，影响重建效果。

本申请提供的基于像素邻域块的全变分正则化项代替传统全变分正则化项：

也即，在对二维图像中的像素点进行梯度运算时，使用与像素相关联的像素邻域块计算二维图像的该像素的梯度。以此代替传统全变分正则项基于像素本身计算二维图像的该像素的梯度所带来的不稳定性。

可选地，本申请提供的基于像素邻域块的全变分模型的保真项为加权最小二乘算法，也就是说，在惩罚加权最小二乘全变分算法加入本申请提供的基于像素邻域块的全变分正则化项作为惩罚项，以构建惩罚加权二乘全变分重建模型。

也即，惩罚加权二乘全变分重建模型为:

其中，Σ表示对角矩阵，且

表示探测器信道i处投影数据的方差。可选地，探测器信道i处投影数据的方差

可以通过公式：

估计得到，

其中，I₀表示X射线入射光子强度，

表示系统电子噪声的方差，

表示探测器信道i处投影数据的均值。当然探测器信道i处投影数据的方差

也可以通过局部邻域方差估计等其他方式得到，在此不做具体限制。

S30：将预设图像数据输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像。

可选地，上述迭代运算为梯度下降算法、共轭梯度下降算法或超松弛迭代算法中的任意一种。

具体地，在迭代过程中，判断第n次迭代重建图像μⁿ是否满足迭代终止条件，若满足，将第n次得到的图像数据作为最终重建图像μ^*。若不满足，则基于第n次迭代重建图像μⁿ更新CT图像，得到更新迭代重建图像μⁿ⁺¹，重复迭代以得到最终重建图像μ^*。其中，迭代终止条件为相邻两次迭代重建结果的相对均方误差小于设定的阈值k，即

k为一正实数。比如，在一具体实施例中，迭代终止条件中设定的阈值k为0.001，阈值的取值根据实际需求设定，在此不做具体限制。

本实施例提供的CT图像构建方法通过构建基于像素邻域块的全变分模型，其中，像素邻域块为以一个像素为中心的像素区域。也即，该CT重建模型引入基于像素邻域块的全变分正则项，且该全变分正则项不同于传统的全变分正则项，本实施例提供的全变分正则项使用与像素相关联的像素邻域块计算二维图像的该像素的梯度，能够消除阶梯效应，保留图像细节，提高重建图像分辨率。

参阅图2，在一实施方式中，在步骤S30可以通过以下两个步骤进行：

S31：利用解析重建算法对投影数据进行重建，获得初始CT图像，并将初始CT图像作为预设图像数据。可选地，利用滤波反投影方法对投影数据进行重建，获得初始CT图像。

S32：将初始CT图像输入至图像重建模型并进行迭代运算，以得到图像重建模型输出的重建CT图像。

滤波反投影的解析重建算法是典型的解析重建算法，其优点是重建速度快，在投影数据完备的情况下重建效果好。本申请为了弥补上述迭代重建算法对CT图像进行重建可能造成的迭代时间长等问题，首先利用滤波反投影重建算法对扫描得到的投影数据进行预重建，以加快重建速度，再利用本申请提供的迭代重建算法更新初始CT图像，以增加重建CT图像细节，提高重建CT图像的分辨率。

参阅图3(a)至图3(e)，图3(a)至图3(e)描述了本申请方法与其他方法的重建效果对比。图3(a)为参考图像，图3(a)是CT设备在标准剂量射线CT医学图像。图3(b)为低剂量数据采用滤波反投影方法得到的重建图像。图3(c)为低剂量图像采用惩罚加权最小二乘方法得到的迭代重建图像。图3(d)为低剂量图像采用惩罚加权最小二乘全变分方法得到的迭代重建图像，且该方法用到的全变分惩罚项为传统全变分算子。图3(e)为采用本申请方法即采用基于像素邻域块的全变分模型得到的迭代重建图像。其中，本方法采用参数为：β＝0.02，δ＝0.001。可以对比看出本申请方法保持了良好的边缘并且分辨率更高更接近参考图像图3(a)。

参阅图4，图4是本申请提供的CT设备一实施例的结构示意图，CT设备400包括：内部总线401，以及通过内部总线401连接的存储器402、处理器403。其中存储器402用于存储计算机程序。处理器403用于执行计算机程序以实现本申请提供的CT图像构建方法的步骤。

处理器403可能是一个中央处理器CPU，或者是专用集成电路ASIC(ApplicationSpecific Integrated Circuit)，或者是被配置成实施本申请实施例的一个或多个集成电路。存储器402用于可执行的指令。存储器402可能包含高速RAM存储器，也可能包括非易失性存储器(non-volatile memory)，例如至少一个磁盘存储器。存储器402也可以是存储器阵列。存储器402还可能被分块，并且块可按一定的规则组合成虚拟卷。存储器402存储的指令可被处理器403执行，以使处理器403能够执行上述任意方法实施例中的CT图像构建的方法。

参阅图5，图5是本申请提供的计算机可读存储介质一实施例的结构示意图。该计算机可读存储介质500上存储有计算机程序501，计算机程序501被处理器执行时实现本申请提供的CT图像构建方法的步骤。计算机存储介质500可以是计算机能够存取的任何可用介质或数据存储设备，包括但不限于磁性存储器(例如软盘、硬盘、磁带、磁光盘(MO)等)、光学存储器(例如CD、DVD、BD、HVD等)、以及半导体存储器(例如ROM、EPROM、EEPROM、非易失性存储器110(NANDFLASH)、固态硬盘(SSD))等。

本实施例提供的CT图像构建方法通过构建基于像素邻域块的全变分模型。也即，该CT重建模型引入基于像素邻域块的全变分正则项，且该全变分正则项不同于传统的全变分正则项，本实施例提供的全变分正则项使用与像素相关联的像素邻域块计算二维图像的该像素的梯度，能够消除阶梯效应，保留图像细节，提高重建CT图像的分辨率。

另外本申请为了弥补上述迭代重建算法对CT图像进行重建可能存在的迭代时间长等问题，首先利用滤波反投影重建算法对扫描得到的投影数据进行预重建，以加快重建速度，再利用本申请提供的迭代重建算法更新初始CT图像，以增加最终CT图像细节，使得到的CT图像更清晰。综上，本申请提供的CT图像构建方法既能加快重建速度，又能提高重建质量。

以上，仅为本申请中的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉该技术的人在本申请所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本申请的包含范围之内，因此，本申请的保护范围应该以权利要求书的保护范围为准。

Claims (12)

1.一种CT图像的构建方法，其特征在于，所述方法包括：

获取CT设备进行扫描时采集的投影数据；

利用所述投影数据构建图像重建模型，其中，所述图像重建模型是基于像素邻域块的全变分模型，所述像素邻域块为以一个像素为中心的像素区域；

将预设图像数据输入至所述图像重建模型并进行迭代运算，以得到所述图像重建模型输出的重建CT图像。

2.根据权利要求1所述的方法，其特征在于，所述将预设图像数据输入至所述图像重建模型并进行迭代运算，以得到所述图像重建模型输出的重建CT图像，包括：

利用解析重建算法对所述投影数据进行重建，获得初始CT图像，并将所述初始CT图像作为所述预设图像数据；

将所述初始CT图像输入至所述图像重建模型并进行迭代运算，以得到所述图像重建模型输出的重建CT图像。

3.根据权利要求2所述的方法，其特征在于，所述利用解析重建算法对所述投影数据进行重建，获得初始CT图像，包括：

利用滤波反投影算法对所述投影数据进行重建，获得初始CT图像。

4.根据权利要求1或2所述的方法，其特征在于，所述基于像素邻域块的全变分模型为：

其中，f为所述投影数据，μ为待重建的所述重建CT图像，H＝{H_i,j}为系统矩阵，i,j分别为待重建所述重建CT图像的像素个数和所述CT设备探测器的探测元个数，D(f,Hμ)为数据保真项，|μ|_TV表示全变分正则化项，β表示平衡保真度和正则化项的超参数。

5.根据权利要求4所述的方法，其特征在于，所述全变分正则化项为：

其中s和t表示所述重建CT图像中衰减系数位置的指标，δ为一常数，用以保持与图像强度的可微性，μ_s,t,l表示以位于所述重建CT图像的第s行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s-1,t,l表示以位于所述重建CT图像的第s-1行以及第t列的像素，为中心的像素邻域块的第l个像素的像素值，μ_s,t-1,l表示以位于所述重建CT图像的第s行以及第t-1列的像素，为中心的像素邻域块的第l个像素的像素值，N_s,t表示所述像素邻域块包含的像素总数。

6.根据权利要求5所述的方法，其特征在于，所述基于像素邻域块的全变分模型为惩罚加权最小二乘全变分算法，其中，所述惩罚加权最小二乘全变分算法引入所述基于像素邻域块的全变分正则项为惩罚项。

7.根据权利要求6所述的方法，其特征在于，所述基于像素邻域块的全变分模型为：

其中，Σ表示对角矩阵，且

表示探测器信道i处的所述投影数据的方差。

8.根据权利要求7所述的方法，其特征在于，所述探测器信道i处的所述投影数据的方差根据以下公式获取：

其中，I₀表示X射线入射光子强度，

表示系统电子噪声的方差，

表示所述探测器信道i处投影数据的均值。

9.根据权利要求4所述的方法，其特征在于，所述系统矩阵是从所述CT设备系统得到的基于面积加权的系统矩阵，或者基于体素加权的系统矩阵。

10.根据权利要求1所述的方法，其特征在于，所述迭代运算为梯度下降算法、共轭梯度下降算法或超松弛迭代算法中的任意一种。

11.一种CT设备，其特征在于，所述CT设备包括：

内部总线，以及通过内部总线连接的存储器和处理器；

所述存储器用于存储计算机程序；

所述处理器用于执行所述计算机程序以实现如权利要求1-10任意一项所述的CT图像构建方法的步骤。

12.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1-10任意一项所述的CT图像构建方法的步骤。

Images