一种多尺度稀疏投影数据快速CT重建方法与流程

本发明涉及CT图像重建技术领域，尤其是一种多尺度稀疏投影数据快速CT重建方法。

背景技术：

CT(Computed Tomography计算机断层摄影术)重建技术，是广泛应用于医疗辅助诊断领域，工业无损检测领域的一门重要技术，其具有无损检测，可视化等方面的有点。而稀疏角度投影数据具有扫描速度快，待检测物体受射线照射少等优点，被广泛应用。但是随着射线的照射角度减少，重建时可以使用的数据也变少，这也使得重建效果变差，重建速度变慢等，严重影响了CT重建技术的发展。为了提高重建的速度，且保证重建的质量，需要通过算法来减少重建算法的。通常重建算法主要分为两大类，一类是以Randon变换为主的解析法，另一类是以代数方法或者统计学方法为主的迭代法。由于解析法重建时需要完备的投影数据，所以在投影角度稀疏的时候，解析法并不能发挥很好的作用。

虽然迭代法在重建上表现了更优的鲁棒性，但是在重建速度上，无法同解析法相比。特别是在引入了正则化约束条件作为限制条件时，重建时间也大大增加。为了保证算法的效率，需要提供一种算法，能够保证在某一类算法的基础上都能有效的提高重建的速度。

针对上述问题，需要发明新的有效算法，在保证重建算法质量的前提下，加快重建速度。

技术实现要素：

为了解决上述技术问题，本发明的目的是：提供一种通过采集稀疏角度投影数据、基于正则化约束的凸优化算法实现的快速CT重建方法。

本发明所采用的技术方案是：一种多尺度稀疏投影数据快速CT重建方法，包括以下步骤：

A、采集稀疏角度投影数据；

B、利用原始投影数据以及原始投影数据对应的投影矩阵进行降采样计算，得到降采样的投影稀疏矩阵和降采样的投影数据；

C、使用凸优化的方法，对原始投影数据进行重建；

D、使用凸优化的迭代方法以及正则化项的约束对步骤B中得到的降采样投影数据进行重建；

E、将上述步骤C和D中的两个重建结果进行融合；

F、将融合图像作为初值，再次利用重建算法进行重建。

进一步，所述步骤A中，在投影数据的采集过程中，采集数据的间隔大于等于4°小于10°。

进一步，所述步骤B中，将原始投影矩阵与原始投影数据进行降采样计算的具体步骤为：

B1、计算原始投影数据对应的投影矩阵Aori，M为投影数据的个数，由投影角度与接收器个数相乘得到；N为重建图像按照行展开为行向量之后的长度，由原始图像的高度与宽度相乘得到，获取投影数据时，每条射线都可以获得一个投影数据；

B2、根据上一步求得的原始投影数据的投影矩阵，计算降采样的投影稀疏矩阵Adown如下：

其中，ratio为降采样倍数，M为投影数据的个数，N为重建图像按照行展开为行向量之后的长度，W为待重建图像的宽度，i、j为投影矩阵的行列坐标值；

B3、根据原始投影数据projori，计算降采样的投影数据projdown:

其中，ratio为降采样倍数，i为投影数据的序号。

进一步，所述步骤C中，使用凸优化的方法，对原始投影数据进行重建的具体步骤为：

C1、对于原始投影数据重建过程，求解目标函数：

其中，Aori为原始投影矩阵，xori为待重建图片，并被转化为列向量，图像中第i行，第j列的像素点对应于列向量中第(i-1)·W+j个像素点，projori为投影数据，W为待重建图像的宽度；

C2、初始化重建结果k＝0；

C3、根据步骤C1的目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；

C4、利用迭代计算新的重建结果zt+1；

C5、更新t＝t+1，判断t是否达到目标迭代次数，如果达到，则更新否则，跳回步骤C4继续迭代；

C6、对进行非负约束的限制：

C7、更新k＝k+1，判断k是否达到目标迭代次数，如果达到，则输出最终结果否则，跳回步骤C3，继续迭代。

进一步，所述步骤D中，使用凸优化的迭代方法以及正则化项的约束对步骤B中得到的降采样投影数据进行重建的具体算法为：

D1、使用降采样的投影数据以及降采样的重建矩阵进行重建，目标函数为：

其中，Adown为将采样投影矩阵，xdown为降采样重建图片，并被转化为列向量，图像中第i行，第j列的像素点对应于列向量中第(i-1)·W+j个像素点，projdown为投影数据，W为待重建图像的宽度；

D2、初始化重建结果k＝0；

D3、根据步骤D1的目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；

D4、利用迭代计算新的重建结果zt+1；

D5、更新t＝t+1，判断t是否达到目标迭代次数，如果达到，则更新否则，跳回D4，继续迭代；

D6、对进行约束：

D7、对步骤D6的计算结果进行正则化项的约束，使得其满足：

D8、将从行向量转化为一个高为H，宽度为W的二维矩阵

其中i,j为图像中的行列坐标，i表示行坐标，j表示列坐标；

D9、引入水平方向和竖直方向的bregman距离并初始化t＝0，以及二维图像在水平方向和竖直方向的梯度

其中i,j分别为图像的行坐标和列坐标；

D10、利用高斯赛德尔公式将待重建图像进行更新

其中λ，μ为正则化参数；

D11、更新水平方向和竖直方向的bregman距离以及水平方向和竖直方向的梯度

其中，和分别关于水平方向和竖直方向的梯度，shrike操作如下：

shrink(a,b)＝max(abs(a)-b,0)×sign(a)

其中

D12、更新k＝k+1,判断k是否达到目标迭代次数，如果达到，则输出最终结果否则，跳回步骤D3，继续迭代。

进一步，所述步骤E中，将两个重建结果进行融合的具体方法为：

E1、对步骤D中得到的原始投影数据的重建结果以及降采样投影数据得到的重建结果进行融合处理，得到结果xnew；

E2、计算降采样投影数据的重建结果的梯度gradientx:

E3、对梯度图像gradientx进行大津法计算，计算得到梯度阈值threshhold，分别得到图像的平滑区以及边缘区；并利用降采样投影数据的重建结果的梯度gradientx的平滑区域去覆盖原始投影数据的对应区域：

其中，ratio为降采样的尺度，为向下取整，i、j为原始投影数据的行、列坐标。

进一步，所述步骤F中，将融合图像作为初值，再次利用重建算法进行重建的具体步骤为：

F1、初始化重建结果k＝0；

F2、根据目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；

F3、利用迭代计算新的重建结果zt+1；

F4、更新t＝t+1，判断t是否达到目标迭代次数，如果达到，则更新否则，跳回步骤F3，继续迭代；

F5、对进行约束：

F6、对步骤F5的计算结果进行正则化项的约束，使得其满足

F7、将从行向量转化为一个高为H，宽度为W的二维矩阵

其中i,j为图像中的行列坐标，i表示行坐标，j表示列坐标；

F8、引入水平方向和竖直方向的bregman距离并初始化t＝0，以及二维图像在水平方向和竖直方向的梯度

其中i,j分别为图像的行坐标和列坐标；

F9、利用高斯赛德尔公式将待重建图像进行更新

其中λ，μ为正则化参数；

F10、更新水平方向和竖直方向的bregman距离以及水平方向和竖直方向的梯度

其中，和分别关于水平方向和竖直方向的梯度，shrike操作如下：

shrink(a,b)＝max(abs(a)-b,0)×sign(a)

其中

F11、更新k＝k+1,判断k是否达到目标迭代次数，如果达到，则输出最终结果否则，跳回步骤B,继续迭代。

本发明的有益效果是：本发明方法通过采集稀疏角度投影数据，并在已有的基于正则化约束的凸优化算法进行重建，并且在保证原算法重建质量的基础上，通过对两个重建结果进行融合加快重建算法的速度，从而实现多尺度稀疏投影数据的快速CT重建。

附图说明

图1为本发明方法的具体步骤流程图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

参照图1，一种多尺度稀疏投影数据快速CT重建方法，包括以下步骤：

A、采集稀疏角度投影数据；

B、利用原始投影数据以及原始投影数据对应的投影矩阵进行降采样计算，得到降采样的投影稀疏矩阵和降采样的投影数据；

C、使用凸优化的方法，对原始投影数据进行重建；

D、使用凸优化的迭代方法以及正则化项的约束对步骤B中得到的降采样投影数据进行重建；

E、将上述步骤C和D中的两个重建结果进行融合；

F、将融合图像作为初值，再次利用重建算法进行重建。

进一步作为优选的实施方式，所述步骤A中，在投影数据的采集过程中，采集数据的间隔大于等于4°小于10°，即在一组采集数据中，数据采集的组数小于等于90组。

进一步作为优选的实施方式，所述步骤B中，将原始投影矩阵与原始投影数据进行降采样计算的具体步骤为：

B1、计算原始投影数据对应的投影矩阵Aori，M为投影数据的个数即射线条数，由投影角度(deg)与接收器个数(rec)相乘得到；N为重建图像按照行展开为行向量之后的长度，由原始图像的高度(H)与宽度(W)相乘得到。获取投影数据时，每条射线都可以获得一个投影数据；射线可以被认为是一条宽度为0的射线，在射线穿过像素时，每条射线都会被它经过的像素方格给分成一小段，我们取投影系数为第i条射线与第j个像素相交的长度。

Aori＝{ai,j}

其中，ai,j为射线源在某个投影角度下与对应角度下某一个接收单元的连线(第i条射s)同重建图像中像素(第j个像素)的交线长度。所以

B2、根据上一步求得的原始投影数据的投影矩阵，计算降采样的投影稀疏矩阵Adown如下：

其中，ratio为降采样倍数，M为投影数据的个数，N为重建图像按照行展开为行向量之后的长度，W为待重建图像的宽度，i、j为投影矩阵的行列坐标值；

B3、根据原始投影数据projori，计算降采样的投影数据projdown:

其中，ratio为降采样倍数，i为投影数据的序号。

进一步作为优选的实施方式，所述步骤C中，使用凸优化的方法，对原始投影数据进行重建的具体步骤为：

C1、对于原始投影数据重建过程，本步骤只采取保真项约束的求取，即求解目标函数：

其中，Aori为原始投影矩阵，xori为待重建图片，并被转化为列向量，图像中第i行，第j列的像素点对应于列向量中第(i-1)·W+j个像素点，projori为投影数据，W为待重建图像的宽度；

C2、初始化重建结果k＝0；

C3、根据步骤C1的目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；首先初始化t＝0。并且计算初始梯度G0，以及初始化梯度下降步长l：

其中Aori为原始投影矩阵，为原始投影矩阵的转置，λ为梯度下降系数，本文取值20。

C4、利用迭代计算新的重建结果zt+1；

zt+1＝zt-l·Gt

其中，zt为迭代t次的重建结果，G^t为迭代第t次时的下降步长。

并且更新梯度Gt+1,以及下降步长l：

C5、更新t＝t+1，判断t是否达到目标迭代次数(Iter)，如果达到，则更新否则，跳回步骤C4继续迭代；

C6、为了保证全部落在解空间内，需要对进行非负约束的限制：

C7、更新k＝k+1，判断k是否达到目标迭代次数(Iter)，如果达到，则输出最终结果否则，跳回步骤C3，继续迭代。

进一步作为优选的实施方式，所述步骤D中，使用凸优化的迭

代方法以及正则化项的约束对步骤B中得到的降采样投影数据

进行重建的具体算法为：

D1、使用降采样的投影数据以及降采样的重建矩阵进行重建，可以在保真项的基础上加入正则化项进行约束，提高重建质量，目标函数为：

其中，Adown为将采样投影矩阵，xdown为降采样重建图片，并被转化为列向量，图像中第i行，第j列的像素点对应于列向量中第(i-1)·W+j个像素点，projdown为投影数据，W为待重建图像的宽度；

D2、初始化重建结果k＝0；

D3、根据步骤D1的目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；首先初始化t＝0。并且计算初始梯度G0，以及初始化梯度下降步长l：

其中Adown为原始投影矩阵，为原始投影矩阵的转置，λ为梯度下降系数，本文取值20。

D4、利用迭代计算新的重建结果zt+1；

Zt+1＝zt-l·Gt

其中，zt为迭代t次的重建结果，G^t为迭代第t次时的下降步长。

并且更新梯度Gt+1,以及下降步长l：

D5、更新t＝t+1，判断t是否达到目标迭代次数(Iter)，如果达到，则更新否则，跳回D4，继续迭代；

D6、为了保证步骤(5)中得到的解完全落在正整数的解空间内，需要对进行约束：

D7、对步骤D6的计算结果进行正则化项的约束，使得其满足：

D8、在得到了之后，需要将从行向量转化为一个高为H，宽度为W的二维矩阵

其中i,j为图像中的行列坐标，i表示行坐标，j表示列坐标；

D9、引入水平方向和竖直方向的bregman距离并初始化t＝0，以及二维图像在水平方向和竖直方向的梯度

其中i,j分别为图像的行坐标和列坐标；

D10、利用高斯赛德尔公式将待重建图像进行更新

其中λ，μ为正则化参数，在本算法中，λ的取值为0.001，u的取值为0.0005；

D11、更新水平方向和竖直方向的bregman距离以及水平方向和竖直方向的梯度

其中，和分别关于水平方向和竖直方向的梯度，shrike操作如下：

shrink(a,b)＝max(abs(a)-b,0)×sign(a)

其中

D12、更新k＝k+1,判断k是否达到目标迭代次数(Iter)，如果达到，则输出最终结果否则，跳回步骤D3，继续迭代。

进一步作为优选的实施方式，所述步骤E中，将两个重建结果进

行融合的具体方法为：

E1、对步骤D中得到的原始投影数据的重建结果以及降采样投影数据得到的重建结果进行融合处理，得到新的算法结果xnew；

E2、计算降采样投影数据的重建结果的梯度gradientx:

E3、对梯度图像gradientx进行大津法计算，计算得到梯度阈值threshhold，分别得到图像的平滑区以及边缘区；并利用降采样投影数据的重建结果的梯度gradientx的平滑区域去覆盖原始投影数据的对应区域：

其中，ratio为降采样的尺度，为向下取整，i、j为原始投影数据的行、列坐标。

进一步作为优选的实施方式，所述步骤F中，将融合图像作为初

值，再次利用重建算法进行重建的具体步骤为：

F1、初始化重建结果k＝0；

F2、根据目标函数，利用梯度下降法来计算待重建图像其中k为迭代的次数；首先初始化t＝0。并且计算初始梯度G0，以及初始化梯度下降步长l：

其中Aori为原始投影矩阵，为原始投影矩阵的转置，λ为梯度下降系数，本文取值20。

F3、利用迭代计算新的重建结果zt+1；

zt+1＝zt-l·Gt

其中，zt为迭代t次的重建结果，G^t为迭代第t次时的下降步长。

并且更新梯度Gt+1,以及下降步长l：

F4、更新t＝t+1，判断t是否达到目标迭代次数(Iter)，如果达到，则更新否则，跳回步骤F3，继续迭代；

F5、为了保证步骤F4中得到的解完全落在正整数的解空间内，需要对进行约束：

F6、对步骤F5的计算结果进行正则化项的约束，使得其满足

F7、在得到了之后，需要将从行向量转化为一个高为H，宽度为W的二维矩阵

其中i,j为图像中的行列坐标，i表示行坐标，j表示列坐标；

F8、引入水平方向和竖直方向的bregman距离并初始化t＝0，以及二维图像在水平方向和竖直方向的梯度

其中i,j分别为图像的行坐标和列坐标；

F9、利用高斯赛德尔公式将待重建图像进行更新

其中λ，μ为正则化参数；在本算法中，λ的取值为0.001，u的取值为0.0005。

F10、更新水平方向和竖直方向的bregman距离以及水平方向和竖直方向的梯度

其中，和分别关于水平方向和竖直方向的梯度，shrike操作如下：

shrink(a,b)＝max(abs(a)-b,0)×sign(a)

其中

F11、更新k＝k+1,判断k是否达到目标迭代次数，如果达到，则输出最终结果否则，跳回步骤B,继续迭代。

以上是对本发明的较佳实施进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可以作出种种的等同变换或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。