CN109903355B - 基于空谱双域张量自相似的能谱ct重建方法 - Google Patents

基于空谱双域张量自相似的能谱ct重建方法 Download PDF

Info

Publication number
CN109903355B
CN109903355B CN201910160075.9A CN201910160075A CN109903355B CN 109903355 B CN109903355 B CN 109903355B CN 201910160075 A CN201910160075 A CN 201910160075A CN 109903355 B CN109903355 B CN 109903355B
Authority
CN
China
Prior art keywords
image
energy spectrum
tensor
similarity
iteration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910160075.9A
Other languages
English (en)
Other versions
CN109903355A (zh
Inventor
张意
夏文军
周激流
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan University
Original Assignee
Sichuan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sichuan University filed Critical Sichuan University
Priority to CN201910160075.9A priority Critical patent/CN109903355B/zh
Publication of CN109903355A publication Critical patent/CN109903355A/zh
Application granted granted Critical
Publication of CN109903355B publication Critical patent/CN109903355B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种基于空谱双域张量自相似的能谱CT重建方法，首先充分利用能谱CT图像在空间域和能谱域的自相似性构建三阶相似张量和能量函数模型，然后采用低秩和稀疏分解将构建的三阶张量单元进行分解，再利用交替方向乘子法对能量函数模型中的目标函数进行优化，去除噪声和伪影，得到重建的CT图像；可以在数据不充分或者噪声过大时，将CT图像完美重建出来，因此可以适用于稀疏角和低剂量情况下的能谱CT重建工作，对于降低CT扫描给人体带来的辐射具有重要意义。

Description

基于空谱双域张量自相似的能谱CT重建方法
技术领域
本发明属于图像处理技术领域，涉及CT图像重建技术，具体涉及一种基于空谱双域张量自相似的能谱CT重建方法。
背景技术
CT(Computed Tomography)，即电子计算机断层扫描，它是利用精确准直的X线束、Y射线、超声波等，与灵敏度极高的探测器一同围绕人体的某一部位作一个接一个的断面扫描，具有扫描时间快，图像清晰等特点。CT已成为当今影像技术中不可缺少的临床及工业检测手段。但是在CT成像时，常常假设X射线在不同能量下光子强度和光子对物质的衰减系数均匀一致，这与实际情况相悖，会导致射束硬化，影响成像质量。
能谱CT可以对不同能量段的X光单独成像，从而避免射束硬化。相比传统CT，能谱CT作为新兴的CT技术还有许多其他优势，比如由于每个能量段内光子对物质的衰减系数基本均匀一致，从而可以实现材料分解，以获得彩色CT图像，提高对病灶的辨识度；可以利用光子对物质衰减系数的跳变进行K-edge成像，以获得高对比度图像提高诊断价值等。能谱CT的应用前景广阔，已成为CT技术发展新的方向。但是与传统CT一样，能谱CT产生的X射线辐射会对人体产生危害，导致基因、肿瘤等疾病。因此，如何在较低的辐射剂量下，获得临床可接受的图像，是能谱CT成像领域的重要课题。
压缩感知(Compressive Sensing,CS)理论被应用于医学成像领域[E.J.Candes,J.Romberg,and T.Tao,“Robust uncertainty principles:Exact signalreconstruction from highly incomplete frequency information,”IEEETransactions on Information Theory,vol.52,no.2,pp.489–509,2006]。CS理论证明，如果信号满足稀疏条件，那么在欠采样条件下能够完美地恢复信号。基于该理论，Xu等提出了利用全变分(Total Variation,TV)最小化的迭代重建算法，获得了高质量的能谱CT图像[Q.Xu,H.Yu,J.Bennett,P.He,R.Zainon,R.Doesburg,A.Opie,M.Walsh,H.Shen,andA.Butler,“Image reconstruction for hybrid true-color micro-CT,”IEEETransactions on Bio-medical Engineering,vol.59,no.6,p.1711,2012]。作为字典学习的推广，张量字典学习(Tensor Dictionary Learning,TDL)方法被Zhang等人应用于能谱CT重建中，该方法训练得到的字典原子为张量，更好地利用了张量结构的特性，在低剂量条件下，有效地去除了噪声[Y.Zhang,X.Mou,G.Wang,and H.Yu,“Tensor-based dictionarylearning for spectral CT reconstruction,”IEEE Transactions on MedicalImaging,vol.36,no.1,pp.142–154,2017]。但稀疏角导致的伪影与随机噪声不同，很难通过字典学习的方法去除。为了有效消除伪影，Wu等人将L0范数与TDL方法结合起来，取得了较好的重建效果[W.Wu,Y.Zhang,Q.Wang,F.Liu,P.Chen,and H.Yu,“Low-dose spectralCT reconstruction using image gradient l0–norm and tensor dictionary,”AppliedMathematical Modelling,vol.63,pp.538–557,2018]。
此外，为了利用能谱CT图像在能谱域中的相似性，一些研究者基于图像在谱域的相似性，提出了一些CT图像重建方法。Gao等人提出了基于PRISM模型(Prior Rank,Intensity and Sparsity Model)的能谱CT重建方法[H.Gao,H.Yu,S.Osher,and G.Wang,“Multi-energy CT based on a prior rank,intensity and sparsity model(PRISM).”Inverse Problems,vol.27,no.11,pp.115 012–115 033(22),2011]，将不同能量段的图像向量化后构建成矩阵，利用低秩矩阵分解进行重建。Li等人将PRISM方法推广到了张量形式，不对图像向量化，在张量的基础上直接约束图像的核范数，提出了tensor PRISM，避免向量化会破坏图像的结构信息，取得了更好的重建效果[L.Li,Z.Chen,G.Wang,J.Chu,andH.Gao,“A tensor prism algorithm for multi-energy CT reconstruction andcomparative studies,”Journal of X-Ray Science Technology,vol.22,no.2,pp.147–163,2014]。Rigie和La Riviere将图像在不同谱段梯度差分的相似性作为惩罚正则项，提出了总变分核范数(Total Nuclear Variation,TNV)最小化方法，有效利用了能谱图像具有相同轮廓的特性[D.S.Rigie,“Joint reconstruction of multi-channel,spectral CTdata via constrained total nuclear variation minimization.”Physics inMedicine&Biology,vol.60,no.5,p.1741,2015]。Yu等提出了能谱压缩感知先验模型，将全光谱图像作为先验图像，全光谱图像相对各谱段的图像具有更小的量子噪声和相似的结构，因此可以作为先验信息来指导各谱段图像的重建[Z.Yu,S.Leng,Z.Li,andC.H.Mccollough,“Spectral prior image constrained compressed sensing(spectralPICCS)for photon-counting computed tomography,”Physics in Medicine&Biology,vol.61,no.18,p.6707,2016]。Kazantsev等人引入了集合差分相似先验，在上一次迭代结果中选择一组谱段图像来指导本次迭代的重建，利用了谱段间的图像轮廓相似性，[D.Kazantsev,J.S.
Figure BDA0001984312420000021
M.Andersen,W.R.B.Lionheart,P.D.Lee,and P.J.Withers,“Joint image reconstruction method with correlative multi-channel prior forx-ray spectral computed tomography,”Inverse Problems,vol.34,no.6,p.064001,2018]。Niu和Kim等人将非局部思想与低秩约束结合起来，将谱域作为非局部窗口，结合非局部图像块的相似性，有效提高了能谱CT的图像质量。
虽然上述多种CT图像重建方法在一定程度上提高了能谱CT的重建图像质量。然而对于稀疏采样所导致的不完备投影数据，采样上述CT重建方法仍无法获取较好的图像质量。
发明内容
本发明的目的旨在提供一种基于空谱双域张量自相似的能谱CT重建方法，利用能谱CT图像在图像域(或称空域)和能谱域(或称谱域)内的相似性，对稀疏采用得到的图像进行鲁棒重建，并使用交替方向乘子法去除噪声和伪影，得到高质量CT重建图像。
本发明的发明思路为，同时利用图像域及能谱域的相似性构造如下的能量函数约束模型：
Figure BDA0001984312420000031
其中A(·)表示以A为***投影矩阵的投影函数，
Figure BDA0001984312420000032
表示重建能谱图像构成的张量，
Figure BDA0001984312420000033
表示投影数据组成的张量，
Figure BDA0001984312420000034
为张量单元，p＝1，2，…，P，P为从图像域选取的图像块数量，GRp(·)是利用从图像域选取的第p个图像块构造张量单元的构造算子，
Figure BDA0001984312420000035
和
Figure BDA0001984312420000036
分别为张量单元分解得到的低秩部分和稀疏部分，λ1和λ2为权重参数。
将上述约束模型转化成非约束模型：
Figure BDA0001984312420000037
Figure BDA0001984312420000038
ρ为松弛因子。
基于上述发明思路，本发明利用交替方向乘子法提出的基于空谱双域张量自相似的能谱CT重建方法，包括以下步骤：
(1)依据CT机先前扫描的投影数据，并使用滤波反投影算法(filteredback-projection,FBP)获取初始能谱CT图像；
(2)采用交替方向乘子法对输入图像进行处理得到重建能谱CT图像：以初始能谱CT图像作为输入图像时为第一次迭代过程，第一次迭代的初始值
Figure BDA0001984312420000039
(即初始能谱CT图像)，即n＝0，除第一次迭代外，每一次迭代时使用的初始值都为上一次迭代得到的结果，对于第n+1次迭代，重建能谱CT图像过程包括以下分步骤：
(21)构建三阶相似张量：当输入图像由M个谱段组成时，选取任一谱段作为参考依据，在该谱段图像上重叠提取P个图像块，然后使用图像匹配算子构建P个三阶相似张量，依次为
Figure BDA0001984312420000041
其中
Figure BDA0001984312420000042
为第n+1次迭代中构建的第p个图像块的三阶相似张量，对于初始能谱CT图像，n＝0；
(22)依据以下公式得到P个三阶相似张量的低秩分量：
Figure BDA0001984312420000043
式中，
Figure BDA0001984312420000044
为第n+1次迭代中使得(1)式中目标函数极小化时
Figure BDA0001984312420000045
的值，并将其作为第p个三阶相似张量
Figure BDA0001984312420000046
的低秩分量，
Figure BDA0001984312420000047
为第n次迭代中
Figure BDA0001984312420000048
的稀疏分量，
Figure BDA0001984312420000049
λ1为权重参数，ρ均为松弛因子，p＝1，2，…，P；
(23)依据以下公式得到P个三阶相似张量的稀疏分量：
Figure BDA00019843124200000410
式中，
Figure BDA00019843124200000411
为第n+1次迭代中
Figure BDA00019843124200000412
的低秩分量，由步骤(22)得到，
Figure BDA00019843124200000413
为第n+1次迭代中使得式(2)中目标函数极小化时
Figure BDA00019843124200000414
的值，λ2为权重参数，p＝1，2，…，P；
(24)依据以下公式得到第n+1次迭代的重建结果：
Figure BDA00019843124200000415
式中，
Figure BDA00019843124200000416
分别为第n+1次迭代中第p个三阶相似张量
Figure BDA00019843124200000417
的低秩分量和稀疏分量，分别由步骤(22)和步骤(23)得到，
Figure BDA00019843124200000418
为为第n+1次迭代中使得(3)式中目标函数极小化时
Figure BDA00019843124200000419
的值，并将其作为第n+1次迭代重建能谱CT图像的输出值，
Figure BDA00019843124200000420
为依据目标函数中待求解
Figure BDA00019843124200000421
构建的第p个三阶相似张量，本质上可以转化成由
Figure BDA00019843124200000422
表示的张量，A(·)为CT机***矩阵的投影函数，
Figure BDA00019843124200000423
为投影数据组成的张量，ρ为松弛因子，P为三阶相似张量的总数；
(3)评价图像质量，对重建能谱CT图像质量进行评价，若图像质量满足设定要求，则输出图像为最终的重建能谱CT图像；若图像质量没有达到设定要求，以重建能谱CT图像为输入图像，返回步骤(2)进行下一次迭代。
上述基于空谱双域张量自相似的能谱CT重建方法，步骤(1)的目的是使用滤波反投影算法对CT机先前扫描的投影数据进行处理，得到的CT图像作为初始能谱CT图像，使用滤波反投影算法构建CT图像的操作为本领域的常规手段，可以参见[A.C.Kak andM.Slaney,“Principles of computerized tomography imaging,”Medical Physics,vol.29,no.1,pp.107–107,2002.]。
上述基于空谱双域张量自相似的能谱CT重建方法，是首先通过图像域及能谱域的相似性构造三阶相似张量，然后采用低秩和稀疏分解，将构建的三阶张量单元分解为低秩分量和稀疏分量，再利用交替方向乘子法对构建的能量函数模型中的目标函数进行优化，得到重建的CT图像。
步骤(21)为利用图像匹配算子构建P个在图像域和能谱域(即空谱双域)均具有高度自相似性的三阶相似张量
Figure BDA0001984312420000051
具体过程包括以下分步骤：
(211)选取任一谱段作为参考依据，获取该谱段的二维CT图像，在该二维CT图像上重叠提取P个图像块；优选的实施方式中，在参考谱段图像上选取的P个图像块的重叠范围为40％～80％；为了使重建的CT图像清晰，选取的图像块总数P不少于5个；
(212)对于第p个图像块，以其中心创建搜索窗，在搜索窗内找到与第p个图像块最相似的z-1个图像块，将第p个图像块及与之相似的z-1个图像块全部转化为列向量，然后将这些列向量组成一个相似矩阵；可以采用块匹配或欧拉距离在参考谱段图像搜索窗内找到第p个图像块的z-1个相似的图像块；为了使重建的CT图像清晰，第p个图像块及与之相似的z-1个图像块的总数，即z不少于5个；所述搜索窗窗口大小可以为选取的图像块大小的2～10倍；
(213)获得第p个图像块及与之相似的z-1个图像块在参考谱段上的位置，分别提取在其余谱段图像上相同位置的图像块，并将这些图像块转化为列向量，构建相应谱段的相似矩阵；
(214)将步骤(212)和步骤(213)获得的M个相似矩阵堆叠成第p个三阶相似张量
Figure BDA0001984312420000052
(215)对于步骤(211)中提取的P个图像块，重复(212)-(214)，便获得P个三阶相似张量
Figure BDA0001984312420000053
步骤(22)中，获取三阶相似张量的低秩分量的实现方式为：为了使(1)式中目标函数极小化，定义
Figure BDA0001984312420000054
Figure BDA0001984312420000055
为
Figure BDA0001984312420000056
沿第k阶展开所得到的矩阵，
Figure BDA0001984312420000057
为
Figure BDA0001984312420000058
的核范数，即
Figure BDA0001984312420000059
Figure BDA00019843124200000510
为
Figure BDA00019843124200000511
的第l大的奇异值，
基于上述分析，
Figure BDA00019843124200000512
的最优解
Figure BDA00019843124200000513
通过构建的以下函数得到：
Figure BDA00019843124200000514
式中，
Figure BDA0001984312420000061
和
Figure BDA0001984312420000062
由
Figure BDA0001984312420000063
沿第k阶展开得到的矩阵经奇异值分解得到的第l大的奇异值和奇异向量，foldk是将矩阵还原成张量的算子；得到的
Figure BDA0001984312420000064
的最优解
Figure BDA0001984312420000065
即为
Figure BDA0001984312420000066
步骤(23)中，获取三阶相似张量的稀疏分量的实现方式为：为了使(2)式中目标函数极小化，将具有M个能谱段的稀疏分量
Figure BDA0001984312420000067
沿第三阶分离得到M个矩阵切片
Figure BDA0001984312420000068
Figure BDA0001984312420000069
是第m个矩阵切片，定义
Figure BDA00019843124200000610
Figure BDA00019843124200000611
为梯度算子，
Figure BDA00019843124200000612
是
Figure BDA00019843124200000613
的全变分，
则
Figure BDA00019843124200000614
i和j分别是矩阵
Figure BDA00019843124200000615
行和列的下标，
基于上述分析，每个切片
Figure BDA00019843124200000616
采用Chambolle投影算法求解，即：
Figure BDA00019843124200000617
式中μ＝λ2/ρ，div是散度算子，Ym是张量
Figure BDA00019843124200000618
沿第三阶分离得到的第m个切片，矩阵q通过迭代得到，其第t+1次迭代过程为
Figure BDA00019843124200000619
其中i和j分别是式中矩阵行和列的下标，τ为设定的迭代步长，t＝0，1，2，…，T，T为迭代次数，q0＝0；
将得到的
Figure BDA00019843124200000620
堆叠所得张量即为
Figure BDA00019843124200000621
步骤(24)中，对
Figure BDA00019843124200000622
采用共轭梯度法进行求解，得到重建能谱CT图像的输出值
Figure BDA00019843124200000623
Figure BDA00019843124200000624
是求解的对象，而
Figure BDA00019843124200000625
为依据
Figure BDA00019843124200000626
拟构建的P个三阶相似张量。若要求解
Figure BDA00019843124200000627
将式(3)中的目标函数对
Figure BDA00019843124200000628
求导，目标函数第一项的求导结果为
Figure BDA00019843124200000629
在目标函数第二项中，如果将
Figure BDA00019843124200000630
按位置求和放回到
Figure BDA00019843124200000631
上，因为图像块是重叠提取的，所以放回的时候，每个像素会被重复求和多次。当确定构建三阶相似张量时所使用的所有图像块时，可以计算得到每个像素被重复使用的次数，即可获得重叠系数矩阵B，B中的每一个元素都代表相应位置的像素被重复使用的次数，因此有
Figure BDA00019843124200000632
Figure BDA00019843124200000633
同理，可以定义一个变量
Figure BDA00019843124200000634
接下来，对目标函数第二项求导得到
Figure BDA00019843124200000635
因此目标函数对
Figure BDA0001984312420000071
的求导结果为
Figure BDA0001984312420000072
若要使得式(3)中的目标函数极小，需要满足导函数为0，即
Figure BDA0001984312420000073
对于该线性方程问题，可通过共轭梯度算法迭代求解得到
Figure BDA0001984312420000074
迭代次数约为5～10次。
上述基于空谱双域张量自相似的能谱CT重建方法，可以根据相邻两次获取的重建能谱CT图像的相对均方误差来对重建能谱CT图像的质量进行评价，相对均方误差如式(8)所示：
Figure BDA0001984312420000075
其中N表示图像中像素的个数，n表示第n次迭代，
Figure BDA0001984312420000076
表示第n次迭代结果的第l*个像素值。当相邻两次获取的重建能谱CT图像的相对均方误差值小于设定值时，判定为迭代收敛，否则判定为迭代不收敛，需要重返步骤(2)进行下一次的迭代，直至当相邻两次获取的重建能谱CT图像的相对均方误差值小于设定值，迭代收敛。
与现有技术相比，本发明提供的基于空谱双域张量自相似的能谱CT重建方法具有以下有益效果：
(1)本发明首先充分利用能谱CT图像在图像域和能谱域的自相似性构建三阶相似张量和能量函数模型，然后采用低秩和稀疏分解将构建的三阶张量单元进行分解，再利用交替方向乘子法对能量函数模型中的目标函数进行优化，去除噪声和伪影，得到重建的CT图像；可以在数据不充分或者噪声过大时，将CT图像完美重建出来，因此可以适用于稀疏角和低剂量情况下的能谱CT重建工作，对于降低CT扫描给人体带来的辐射具有重要意义。
(2)本发明采用的低秩和稀疏分解方法，低秩分量代表该张量结构内容一致性的成分，稀疏分量代表该张量图像块相异的空间特征和能谱特征，因此能够在去除噪声和伪影的同时，保留更多的图像细节。
(3)本发明采用交替方向乘子法优化目标函数，该方法对于大规模计算可靠有效，并且适用于移植到GPU并行计算，因此使得整个算法能够在可接受的时间范围内获得较高的图像质量，对于临床的实际应用有重要意义。
附图说明
图1为本发明基于空谱双域张量自相似的能谱CT重建方法中三阶相似张量构建流程示意图。
图2为本发明三阶相似张量的低秩部分沿三阶展开流程示意图。
图3为本发明三阶相似张量的稀疏部分沿第三阶分离流程示意图。
图4为本发明应用例1中临床采集的腹腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图5为本发明应用例1中采用FBP算法重建的腹腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图6为本发明应用例1中采用基于空谱双域张量自相似的能谱CT重建方法重建的腹腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图7为本发明应用例2中临床采集的胸腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图8为本发明应用例2中采用FBP算法重建的胸腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图9为本发明应用例2中采用基于空谱双域张量自相似的能谱CT重建方法重建的胸腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图10为本发明应用例3中采用FBP算法重建的加噪腹腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
图11为本发明应用例3中采用基于空谱双域张量自相似的能谱CT重建方法重建的加噪腹腔能谱CT图像，其中(a)表示60keV能谱段图像，(b)表示70keV能谱段图像，(c)表示80keV能谱段图像，(d)表示90keV能谱段图像，(e)表示100keV能谱段图像。
具体实施方式
以下将结合附图给出本发明实施例，并通过实施例对本发明的技术方案进行进一步的清楚、完整说明。显然，所述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明内容，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施例，都属于本发明所保护的范围
实施例
本实施例提供的基于空谱双域张量自相似的能谱CT重建方法，包括以下步骤：
(1)依据CT机先前扫描的投影数据，并使用滤波反投影算法(FBP)获取初始能谱CT图像。
本步骤使用FBP算法对CT机先前扫描的投影数据进行处理构建CT图像的操作为本领域的常规手段，可以参见[A.C.Kak and M.Slaney,“Principles of computerizedtomography imaging,”Medical Physics,vol.29,no.1,pp.107–107,2002.]。
(2)采用交替方向乘子法对输入图像进行处理得到重建能谱CT图像：以初始能谱CT图像作为输入图像时为第一次迭代过程，第一次迭代的初始值
Figure BDA0001984312420000091
(即初始能谱CT图像)，即n＝0，除第一次迭代外，每一次迭代时使用的初始值都为上一次迭代得到的结果。
本步骤首先通过图像域及能谱域的相似性构造三阶相似张量，然后采用低秩和稀疏分解，将构建的三阶张量单元分解为低秩分量和稀疏分量，再利用交替方向乘子法对构建的能量函数模型中的目标函数进行优化，得到重建的CT图像。对于第n+1次迭代，重建能谱CT图像过程包括以下分步骤(21)～(24)。
(21)构建三阶相似张量：当输入图像由M个谱段组成时，选取任一谱段作为参考依据，在该谱段图像上重叠提取P个图像块，然后使用图像匹配算子构建P个三阶相似张量，依次为
Figure BDA0001984312420000092
其中
Figure BDA0001984312420000093
为第n+1次迭代中构建的第p个图像块的三阶相似张量，对于初始能谱CT图像，n＝0。
本步骤利用图像匹配算子构建P个在图像域和能谱域(即空谱双域)均具有高度自相似性的三阶相似张量
Figure BDA0001984312420000094
三阶相似张量的具体构建过程，如图1所示，包括以下分步骤(211)～(215)。
(211)选取任一谱段作为参考依据，获取该谱段的二维CT图像XR，在该二维CT图像XR上重叠提取P个大小为w×w的图像块r1,r2,…,rp,…,rp；提取时的重叠系数约为40％～80％；为了使重建的CT图像清晰，本实施例选取的图像块总数P不少于5个。例如当图像块为6×6大小、重叠系数为50％时，每隔3个像素提取一个图像块。
(212)对于第p个图像块rp，以其中心创建搜索窗，在搜索窗内以步长为1个像素，采用欧拉距离
Figure BDA0001984312420000095
Do,s表示图像块Xo和图像块Xs之间的欧式距离)去寻找与第p个图像块rp最相似的z-1个图像块，将第p个图像块rp及与之相似的z-1个图像块全部转化成长度为w2的列向量，然后将这些列向量组成一个相似矩阵；在这个矩阵中，共有z个列向量，这z个列向量在结构和内容上都有很强的相似性。为了使重建的CT图像清晰，本实施例搜索框内相似图像块的总数z不少于5个，所述搜索窗窗口大小可以为选取的图像块大小的2～10倍。
(213)获得第p个图像块rp及与之相似的z-1个图像块在参考谱段上的位置，分别提取在其余谱段图像上相同位置的图像块，并将这些图像块转化为列向量，构建相应谱段的相似矩阵。
(214)将步骤(212)和步骤(213)获得的M个相似矩阵堆叠成第p个三阶相似张量
Figure BDA0001984312420000101
对于具有M个能谱段的能谱CT图像，共获得M个相似矩阵，将这些相似矩阵堆叠成一个三阶相似张量，获得
Figure BDA0001984312420000102
在这个三阶相似张量中，每一个列向量代表一个图像块，沿着第二阶是在一个能谱段内的局部相似图像块，沿着第三阶则是不同能谱段内但位置相同的图像块。
(215)对于步骤(211)中提取的P个图像块，重复(212)-(214)，便获得P个三阶相似张量
Figure BDA0001984312420000103
通过上述步骤(211)～(215)便获得了针对输入图像的P个三阶相似张量。对于每个三阶相似张量，具有高度内相似性，因此理论上为低秩张量。但是张量的组成成分是由不同空间位置和不同谱段所提取的图像块，图像块之间有一定的差异。因此将这些张量分解成低秩分量和稀疏分量，低秩分量代表该张量结构内容一致性的成分，稀疏分量代表该张量图像块相异的空间特征和能谱特征。该方法可以使得在去除噪声和伪影的同时，可以保留结构细节，获得更好的图像质量。去除噪声和伪影的过程本实施例采用交替方向乘子法进行，该方法利于大规模计算，对于GPU计算有很好的移植性。下面步骤(22)是对步骤(21)得到的P个三阶相似张量进行低秩分解，得到低秩分量。步骤(23)是对步骤(21)得到的P个三阶相似张量进行稀疏分解，得到稀疏分量。
(22)依据以下公式得到P个三阶相似张量的低秩分量：
Figure BDA0001984312420000104
式中，
Figure BDA0001984312420000105
为第n+1次迭代中使得(1)式中目标函数极小化时
Figure BDA0001984312420000106
的值，并将其作为第p个三阶相似张量
Figure BDA0001984312420000107
的低秩分量，
Figure BDA0001984312420000108
为第n次迭代中
Figure BDA0001984312420000109
的稀疏分量，
Figure BDA00019843124200001010
λ1为权重参数，ρ均为松弛因子，p＝1，2，…，P。
该步骤利使用张量核范数来约束低秩分量的秩，为了使(1)式中目标函数极小化，本实施例将低秩张量的核范数定义为
Figure BDA00019843124200001011
Figure BDA00019843124200001012
为
Figure BDA00019843124200001013
沿第k阶展开所得到的矩阵，展开方法如图2所示，沿第一阶展开得到若干尺寸为N1×N2的矩阵切片拼接在一起，沿第二阶展开得到若干尺寸为N2×N3的矩阵切片拼接在一起，沿第三阶展开得到若干尺寸为N3×N1的矩阵切片拼接在一起。
Figure BDA0001984312420000111
为
Figure BDA0001984312420000112
的核范数。矩阵的核范数为该矩阵的奇异值之和，即
Figure BDA0001984312420000113
Figure BDA0001984312420000114
为
Figure BDA0001984312420000115
的第l大的奇异值。
本实施例使用奇异值阈值收缩算法进行优化，
Figure BDA0001984312420000116
的最优解
Figure BDA0001984312420000117
通过构建的以下函数得到：
Figure BDA0001984312420000118
Figure BDA0001984312420000119
和
Figure BDA00019843124200001110
由
Figure BDA00019843124200001111
沿第k阶展开得到的矩阵经奇异值分解得到的第l大的奇异值和奇异向量，
Figure BDA00019843124200001112
表示
Figure BDA00019843124200001113
的转置。计算
Figure BDA00019843124200001114
沿k阶展开得到的矩阵，并经过奇异值分解，
Figure BDA00019843124200001115
得到该矩阵的奇异值和奇异向量，即
Figure BDA00019843124200001116
和
Figure BDA00019843124200001117
按照式(4)，将奇异值减去松弛因子
Figure BDA00019843124200001118
负的结果置0，再乘以奇异向量还原成矩阵。分别沿第一、二、三阶实行上述操作，再将得到的三个矩阵还原成张量，取平均值，即为
Figure BDA00019843124200001119
的最优解。foldk是将矩阵沿第k阶还原成张量的算子，是展开的逆操作；得到的
Figure BDA00019843124200001120
的最优解
Figure BDA00019843124200001121
即为
Figure BDA00019843124200001122
(23)依据以下公式得到P个三阶相似张量的稀疏分量：
Figure BDA00019843124200001123
式中，
Figure BDA00019843124200001124
为第n+1次迭代中
Figure BDA00019843124200001125
的低秩分量，由步骤(22)得到，
Figure BDA00019843124200001126
为第n+1次迭代中使得式(2)中目标函数极小化时
Figure BDA00019843124200001127
的值，λ＝和ρ均为松弛因子，p＝1，2，…，P。
本步骤利用总变分来约束稀疏分量的稀疏性，对于张量，如图3所示，将具有M个能谱段的稀疏分量
Figure BDA00019843124200001128
沿第三阶分离得到M个矩阵切片
Figure BDA00019843124200001129
Figure BDA00019843124200001130
是第m个矩阵切片，总变分定义为
Figure BDA00019843124200001131
Figure BDA00019843124200001132
为梯度算子。
Figure BDA00019843124200001133
是
Figure BDA00019843124200001134
的全变分(total variation,TV)，定义为
Figure BDA00019843124200001135
中所有像素对横向梯度和纵向梯度的平方和开根号再求和的结果，而在图像中，梯度定义为相邻两个像素的灰度差分。则
Figure BDA00019843124200001136
的表达式为：
当
Figure BDA00019843124200001137
i和j分别是矩阵行和列的下标。
基于上述分析，每个切片
Figure BDA00019843124200001138
采用Chambolle投影算法优化，即：
Figure BDA0001984312420000121
式中M＝λ2/ρ，div是散度算子，Ym是张量
Figure BDA0001984312420000122
沿第三阶分离得到的第m个切片，矩阵q通过迭代得到，其第t+1次迭代过程为
Figure BDA0001984312420000123
其中τ为设定的迭代步长，t＝0，1，2，…，T，T为迭代次数，约为10～20次，q0＝0。
将得到的
Figure BDA0001984312420000124
堆叠所得张量即为
Figure BDA0001984312420000125
(24)依据以下公式得到第n+1次迭代的重建结果，即获得重建CT能谱图像：
Figure BDA0001984312420000126
式中，
Figure BDA0001984312420000127
分别为第n+1次迭代中第p个三阶相似张量
Figure BDA0001984312420000128
的低秩分量和稀疏分量，分别由步骤(22)和步骤(23)得到，
Figure BDA0001984312420000129
为为第n+1次迭代中使得(3)式中目标函数极小化时
Figure BDA00019843124200001210
的值，并将其作为第n+1次迭代重建能谱CT图像的输出值，
Figure BDA00019843124200001211
为依据目标函数中待求解
Figure BDA00019843124200001212
构建的第p个三阶相似张量，本质上可以转化成由
Figure BDA00019843124200001213
表示的张量，A(·)为CT机***矩阵的投影函数，
Figure BDA00019843124200001214
为投影数据组成的张量，ρ为松弛因子，P为三阶相似张量的总数；
本步骤优化模型为L2范数最小化模型。
Figure BDA00019843124200001215
是求解的对象，而
Figure BDA00019843124200001216
为依据
Figure BDA00019843124200001217
拟构建的P个三阶相似张量。若要求解
Figure BDA00019843124200001218
将式(3)中的目标函数对
Figure BDA00019843124200001219
求导，目标函数第一项的求导结果为
Figure BDA00019843124200001220
在目标函数第二项中，如果将
Figure BDA00019843124200001221
Figure BDA00019843124200001222
按位置求和放回到
Figure BDA00019843124200001223
上，因为图像块是重叠提取的，所以放回的时候，每个像素会被重复求和多次。当确定构建三阶相似张量时所使用的所有图像块时，可以计算得到每个像素被重复使用的次数，即可获得重叠系数矩阵B，B中的每一个元素都代表相应位置的像素被重复使用的次数，因此有
Figure BDA00019843124200001224
Figure BDA00019843124200001225
同理，可以定义一个变量
Figure BDA00019843124200001226
接下来，对目标函数第二项求导得到
Figure BDA00019843124200001227
因此目标函数对
Figure BDA00019843124200001228
的求导结果为
Figure BDA00019843124200001229
若要使得式(3)中的目标函数极小，需要满足导函数为0，即
Figure BDA00019843124200001230
对于该线性方程问题，可通过共轭梯度算法迭代求解得到
Figure BDA0001984312420000131
迭代次数约为5～10次。
(3)评价图像质量，对重建能谱CT图像质量进行评价，若图像质量满足设定要求，则输出图像为最终的重建能谱CT图像；若图像质量没有达到设定要求，以重建能谱CT图像为输入图像，返回步骤(2)进行下一次迭代。
本步骤根据相邻两次获取的重建能谱CT图像的相对均方误差来对重建能谱CT图像的质量进行评价，相对均方误差如式(8)所示：
Figure BDA0001984312420000132
其中N表示图像中像素的个数，n表示第n次迭代，
Figure BDA0001984312420000133
表示第n次迭代结果的第l*个像素值。当相邻两次获取的重建能谱CT图像的相对均方误差值小于设定值时，判定为迭代收敛，停止迭代，将此次迭代结果作为最终的重建能谱CT图像；否则判定为迭代不收敛，需要重返步骤(2)进行下一次的迭代，直至迭代收敛。
通过上述步骤(1)～(3)，在数据不充分或者噪声过大情况下，仍能获得质量较好的CT图像。
应用例1
本应用例采用的原始数据为来自目前最先进的宝石能谱CT(GE Discovery CT750HD)，CT扫描过程中采集了数千个角度的数据。图4显示了一组腹腔的能谱CT临床数据，包括60、70、80、90、100keV五个能谱段的图像。由于为了说明本发明重建能谱CT图像的成像效果，本应用例采用的是已经扫描得到的临床数据。在实际应用中，可以直接使用CT机采集的投影数据。
本应用例中将先前扫描的腹腔全部角度投影数据通过射线驱动算法采样到65个角度，得到65个角度的投影数据进行试验。
对该65个角度的投影数据仿真生成65个角度的投影数据，仿真方法可以采用本领域已经披露的常规手段得到，本实施例采用的方法见参考文献R.L.Siddon,“Fastcalculation of the exact radiological path for a three-dimensional CT array,”Medical Physics,vol.12,no.2,pp.252–5,1985。由于该数据极度欠采样，采用传统的FBP算法对投影数据进行处理得到腹腔重建能谱CT图像(处理过程参考A.C.Kak andM.Slaney,“Principles of computerized tomography imaging,”Medical Physics,vol.29,no.1,pp.107–107,2002)，结果如图5所示，从图5中可以看出，得到的腹腔重建能谱CT图像布满伪影，难以对腹腔细节进行分辨。
采用实施例1中给出的基于空谱双域张量自相似的能谱CT重建方法对投影数据进行处理，处理过程如下：
(1)依据CT机先前扫描的投影数据，并使用FBP算法获取初始能谱CT图像。
本步骤使用传统的FBP算法对投影数据进行处理，得到初始能谱CT图像，处理过程参考A.C.Kak and M.Slaney,“Principles of computerized tomography imaging,”Medical Physics,vol.29,no.1,pp.107–107,2002。
(2)采用交替方向乘子法对输入图像进行处理得到重建能谱CT图像。以初始能谱CT图像作为输入图像时为第一次迭代过程，第一次迭代的初始值
Figure BDA0001984312420000144
(即初始能谱CT图像)，即n＝0，除第一次迭代外，每一次迭代时使用的初始值都为上一次迭代得到的结果，对于第n+1次迭代，重建能谱CT图像过程包括以下分步骤：
(21)构建三阶相似张量，本应用例中输入图像由M(＝5)个谱段组成时，选取任一谱段作为参考依据，在该谱段图像上重叠提取P(＝7225)个图像块，然后使用图像匹配算子构建P个的三阶相似张量，具体包括以下分步骤：
(211)选取任一谱段作为参考依据，获取该谱段的二维CT图像XR，在该二维CT图像XR上重叠提取P个大小为w×w(w＝6)的图像块r1,r2,…,rp,…,rP；提取时的重叠系数约为50％。
(212)对于第p个图像块rp，以其中心创建搜索窗(大小为15×15)，在搜索窗内以步长为1个像素，采用欧拉距离
Figure BDA0001984312420000141
Do,s表示图像块Xo和图像块Xs之间的欧式距离)去寻找与第p个图像块rp最相似的z-1(z＝12)个图像块，将第p个图像块rp及与之相似的11个图像块全部转化成长度为w2的列向量，然后将这些列向量组成一个相似矩阵；在这个矩阵中，共有12个列向量，这12个列向量在结构和内容上都有很强的相似性。
(213)获得第p个图像块rp及与之相似的11个图像块在参考谱段上的位置，分别提取在其余谱段图像上相同位置的图像块，并将这些图像块转化为列向量，构建相应谱段的相似矩阵。
(214)将步骤(212)和步骤(213)获得的5个相似矩阵堆叠成第p个三阶相似张量
Figure BDA0001984312420000142
对于具有5个能谱段的能谱CT图像，共获得5个相似矩阵，将这些相似矩阵堆叠成一个三阶相似张量，获得
Figure BDA0001984312420000143
在这个三阶相似张量中，每一个列向量代表一个图像块，沿着第二阶是在一个能谱段内的局部相似图像块，沿着第三阶则是不同能谱段内但位置相同的图像块。
(215)对于步骤(211)中提取的P个图像块，重复(212)-(214)，便获得P个三阶相似张量
Figure BDA0001984312420000151
其中
Figure BDA0001984312420000152
为第n+1次迭代中构建的第p个图像块的三阶相似张量。
(22)依据以下公式得到P个三阶相似张量的低秩分量：
Figure BDA0001984312420000153
式中，
Figure BDA0001984312420000154
为第n+1次迭代中使得(1)式中目标函数极小化时
Figure BDA0001984312420000155
的值，并将其作为第p个三阶相似张量
Figure BDA0001984312420000156
的低秩分量，
Figure BDA0001984312420000157
为第n次迭代中
Figure BDA0001984312420000158
的稀疏分量，
Figure BDA0001984312420000159
λ1为权重参数，ρ均为松弛因子，λ1和ρ可以根据获得的图像情况，进行设置，取值为λ1＝1×10-5,ρ＝5×10-4，p＝1，2，…，P。
本应用例采取的优化方法与实施例1中的相同，
Figure BDA00019843124200001510
的最优解
Figure BDA00019843124200001511
通过构建的以下函数得到：
Figure BDA00019843124200001512
Figure BDA00019843124200001513
和
Figure BDA00019843124200001514
由
Figure BDA00019843124200001515
沿第k阶展开得到的矩阵经奇异值分解得到的第l大的奇异值和奇异向量，
Figure BDA00019843124200001516
表示
Figure BDA00019843124200001517
的转置。计算
Figure BDA00019843124200001518
沿k阶展开得到的矩阵，并经过奇异值分解，
Figure BDA00019843124200001519
得到该矩阵的奇异值和奇异向量，即
Figure BDA00019843124200001520
和
Figure BDA00019843124200001521
按照式(4)，将奇异值减去松弛因子
Figure BDA00019843124200001522
中负的结果置0，再乘以奇异向量还原成矩阵。分别沿第一、二、三阶实行上述操作，再将得到的三个矩阵还原成张量，取平均值，即为
Figure BDA00019843124200001523
的最优解。foldk是将矩阵沿第k阶还原成张量的算子，是展开的逆操作；得到的
Figure BDA00019843124200001524
的最优解
Figure BDA00019843124200001525
即为
Figure BDA00019843124200001526
(23)依据以下公式得到P个三阶相似张量的稀疏分量：
Figure BDA00019843124200001527
式中，
Figure BDA00019843124200001528
为第n+1次迭代中
Figure BDA00019843124200001529
的低秩分量，由步骤(22)得到，
Figure BDA00019843124200001530
为第n+1次迭代中使得式(2)中目标函数极小化时
Figure BDA00019843124200001531
的值，λ2为权重参数，λ2可以根据获得的图像情况，进行设置，取值为λ2＝1.67×10-6，p＝1，2，…，P。
本应用例采取的优化方法与实施例1中的相同，将具有5个能谱段的稀疏分量
Figure BDA00019843124200001532
沿第三阶分离得到5个矩阵切片
Figure BDA00019843124200001533
Figure BDA00019843124200001534
是第m个矩阵切片，
每个切片
Figure BDA0001984312420000161
采用Chambolle投影算法优化，即：
Figure BDA0001984312420000162
式中M＝λ2/ρ，div是散度算子，Ym是张量
Figure BDA0001984312420000163
沿第三阶分离得到的第m个切片，矩阵q通过迭代得到，其第t+1次迭代过程为
Figure BDA0001984312420000164
其中i和j分别是式中矩阵行和列的下标，τ为设定的迭代步长，t＝0，1，2，…，T，T＝20，q0＝0。
将得到的
Figure BDA0001984312420000165
堆叠所得张量即为
Figure BDA0001984312420000166
(24)依据以下公式得到第n+1次迭代的重建结果，即获得重建CT能谱图像：
Figure BDA0001984312420000167
式中，
Figure BDA0001984312420000168
分别为第n+1次迭代中第p个三阶相似张量
Figure BDA0001984312420000169
的低秩分量和稀疏分量，分别由步骤(22)和步骤(23)得到，
Figure BDA00019843124200001610
为为第n+1次迭代中使得(3)式中目标函数极小化时
Figure BDA00019843124200001611
的值，并将其作为第n+1次迭代重建能谱CT图像的输出值，
Figure BDA00019843124200001612
为依据目标函数中待求解
Figure BDA00019843124200001613
构建的第p个三阶相似张量，本质上可以转化成由
Figure BDA00019843124200001614
表示的张量，A(·)为CT机***矩阵的投影函数，
Figure BDA00019843124200001615
为投影数据组成的张量，ρ为松弛因子，P为三阶相似张量的总数。
本应用例采取的优化方法与实施例1中的相同。
Figure BDA00019843124200001616
是求解的对象，而
Figure BDA00019843124200001617
为依据
Figure BDA00019843124200001618
拟构建的P个三阶相似张量。若要求解
Figure BDA00019843124200001619
将式(3)中的目标函数对
Figure BDA00019843124200001620
求导，目标函数第一项的求导结果为
Figure BDA00019843124200001621
在目标函数第二项中，如果将
Figure BDA00019843124200001622
按位置求和放回到
Figure BDA00019843124200001623
上，因为图像块是重叠提取的，所以放回的时候，每个像素会被重复求和多次。当确定构建三阶相似张量时所使用的所有图像块时，可以计算得到每个像素被重复使用的次数，即可获得重叠系数矩阵B，B中的每一个元素都代表相应位置的像素被重复使用的次数，因此有
Figure BDA00019843124200001624
Figure BDA00019843124200001625
同理，可以定义一个变量
Figure BDA00019843124200001626
接下来，对目标函数第二项求导得到
Figure BDA00019843124200001627
因此目标函数对
Figure BDA00019843124200001628
的求导结果为
Figure BDA00019843124200001629
若要使得式(3)中的目标函数极小，需要满足导函数为0，即
Figure BDA0001984312420000171
对于该线性方程问题，可通过共轭梯度算法迭代求解得到
Figure BDA0001984312420000172
迭代次数为5次。
(3)评价图像质量，对重建能谱CT图像质量进行评价，若图像质量满足设定要求，则输出图像为最终的重建能谱CT图像；若图像质量没有达到设定要求，以重建能谱CT图像为输入图像，返回步骤(2)进行下一次迭代。
本步骤根据相邻两次获取的重建能谱CT图像的相对均方误差来对重建能谱CT图像的质量进行评价，相对均方误差如式(8)所示：
Figure BDA0001984312420000173
其中N表示图像中像素的个数，n表示第n次迭代，
Figure BDA0001984312420000174
表示第n次迭代结果的第l*个像素值。当相邻两次获取的重建能谱CT图像的相对均方误差值小于设定值(本应用例取值为1×10-5)时，判定为迭代收敛，停止迭代，将此次迭代结果作为最终的重建能谱CT图像；否则判定为迭代不收敛，需要重返步骤(2)进行下一次的迭代，直至迭代收敛。
通过上述基于空谱双域张量自相似的能谱CT重建处理过程对采集的数据进行处理得到的重建能谱CT图像如图6所示，图中伪影基本去除，并且细节保护的很完整。
腹腔图像器官组织丰富，具有大量细节。在稀疏角情况下，如果采用FBP方法进行重建，得到的结果伪影严重，难以分辨细节。采用本发明提供的方法进行重建，得到的图像内容清晰，结构完整，细节丰富。65个角度的扫描相比初始数据采集的数千个角度，大大减少了对病人的X射线辐射。因此，采用本发明提供的基于空谱双域张量自相似的能谱CT重建方法，能够在保证图像质量的前提下，有效降低能谱CT扫描对人体造成的伤害。
应用例2
本应用例采用的原始数据为来自目前最先进的宝石能谱CT(GE Discovery CT750HD)，CT扫描过程中采集了数千个角度的数据。图7显示了一组胸腔的能谱CT临床数据，包括60、70、80、90、100keV五个能谱段的图像。
本应用例中将先前扫描的腹腔全部角度投影数据通过射线驱动算法采样到65个角度，得到65个角度的投影数据进行试验。
对该65个角度的投影数据仿真生成65个角度的投影数据，仿真方法可以采用本领域已经披露的常规手段得到，本实施例采用的方法见参考文献R.L.Siddon,“Fastcalculation of the exact radiological path for a three-dimensional CT array,”Medical Physics,vol.12,no.2,pp.252–5,1985。采用传统的FBP算法对仿真的投影数据进行处理得到胸腔重建能谱CT图像(处理过程参考A.C.Kak and M.Slaney,“Principles ofcomputerized tomography imaging,”Medical Physics,vol.29,no.1,pp.107–107,2002)，结果如图8所示，从图8中可以看出，得到的腹腔重建能谱CT图像布满伪影，血管结构不清晰，难以分辨。
采用与应用例1中提供的相同的基于空谱双域张量自相似的能谱CT重建方法对仿真的投影数据进行处理，得到的重建能谱CT图像如图9所示，图中伪影基本去除，并且内容清晰，血管结构完整。
应用例3
在实际临床上，因为X射线光子的发射服从泊松分布，所以得到的投影数据会含有噪声，该噪声服从泊松分布。为了验证本发明提供的图像重建方法在临床上的鲁棒性，本应用例在应用例1得到的腹腔仿真投影数据中加入泊松噪声进行实验。
根据朗伯-比尔定律，探测器接收到的光子数为
Figure BDA0001984312420000181
其中
Figure BDA0001984312420000182
为仿真投影数据，I0为X射线源发射的光子数。在仿真噪声数据的时候，设置I0＝1×106，通过无噪的仿真投影数据，可以计算得到理想情况下探测器接收到的光子数I，以I为参数利用泊松模型进行建模，使用建模后的I根据上式计算得到含有噪声的投影数据。这样泊松噪声就被加入到仿真生成65个角度的投影数据中，即得到试验用投影数据。
采用传统的FBP算法对得到的试验用投影数据进行处理，得到的腹腔重建能谱CT图像，如图10所示，图中出现了许多噪声，相比与图5，图像更加不清晰。
采用与应用例1中提供的相同的基于空谱双域张量自相似的能谱CT重建方法对试验用投影数据进行处理，得到的重建能谱CT图像如图11所示，与图6中的图像相差不大，说明本发明提供的基于空谱双域张量自相似的能谱CT重建方法对噪声具有很好的鲁棒性。

Claims (9)

1.一种基于空谱双域张量自相似的能谱CT重建方法，其特征在于包括以下步骤：
(1)依据CT机先前扫描的投影数据，并使用滤波反投影算法获取初始能谱CT图像；
(2)采用交替方向乘子法对输入图像进行处理得到重建能谱CT图像：以初始能谱CT图像作为输入图像时为第一次迭代过程，第一次迭代的初始值
Figure FDA0001984312410000011
即n＝0，除第一次迭代外，每一次迭代时使用的初始值都为上一次迭代得到的结果，对于第n+1次迭代，重建能谱CT图像过程包括以下分步骤：
(21)构建三阶相似张量：当输入图像由M个谱段组成时，选取任一谱段作为参考依据，在该谱段图像上重叠提取P个图像块，然后使用图像匹配算子构建P个三阶相似张量，依次为
Figure FDA0001984312410000012
其中
Figure FDA0001984312410000013
为第n+1次迭代中构建的第p个图像块的三阶相似张量，对于初始能谱CT图像，n＝0；
(22)依据以下公式得到P个三阶相似张量的低秩分量：
Figure FDA0001984312410000014
式中，
Figure FDA0001984312410000015
为第n+1次迭代中使得(1)式中目标函数极小化时
Figure FDA0001984312410000016
的值，并将其作为第p个三阶相似张量
Figure FDA0001984312410000017
的低秩分量，
Figure FDA0001984312410000018
为第n次迭代中
Figure FDA0001984312410000019
的稀疏分量，
Figure FDA00019843124100000110
λ1为权重参数，ρ均为松弛因子，p＝1，2，…，P；
(23)依据以下公式得到P个三阶相似张量的稀疏分量：
Figure FDA00019843124100000111
式中，
Figure FDA00019843124100000112
为第n+1次迭代中
Figure FDA00019843124100000113
的低秩分量，由步骤(22)得到，
Figure FDA00019843124100000114
为第n+1次迭代中使得式(2)中目标函数极小化时
Figure FDA00019843124100000115
的值，λ2为权重参数，p＝1，2，…，P；
(24)依据以下公式得到第n+1次迭代的重建结果：
Figure FDA00019843124100000116
式中，
Figure FDA00019843124100000117
分别为第n+1次迭代中第p个三阶相似张量
Figure FDA00019843124100000118
的低秩分量和稀疏分量，分别由步骤(22)和步骤(23)得到，
Figure FDA00019843124100000119
为第n+1次迭代中使得(3)式中目标函数极小化时
Figure FDA00019843124100000120
的值，并将其作为第n+1次迭代重建能谱CT图像的输出值，
Figure FDA00019843124100000121
为依据目标函数中
Figure FDA00019843124100000122
构建的第p个三阶相似张量，A(·)为CT机***矩阵的投影函数，
Figure FDA00019843124100000123
为投影数据组成的张量，ρ为松弛因子，P为三阶相似张量的总数；
(3)评价图像质量，对重建能谱CT图像质量进行评价，若图像质量满足设定要求，则输出图像为最终的重建能谱CT图像；若图像质量没有达到设定要求，以重建能谱CT图像为输入图像，返回步骤(2)进行下一次迭代。
2.根据权利要求1所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于步骤(21)使用图像匹配算子构建P个三阶相似张量
Figure FDA00019843124100000214
的过程包括以下分步骤：
(211)选取任一谱段作为参考依据，获取该谱段的二维CT图像，在该二维CT图像上重叠提取P个图像块；
(212)对于第p个图像块，以其中心创建搜索窗，在搜索窗内找到与第p个图像块最相似的z-1个图像块，将第p个图像块及与之相似的z-1个图像块全部转化为列向量，然后将这些列向量组成一个相似矩阵；
(213)获得第p个图像块及与之相似的z-1个图像块在参考谱段上的位置，分别提取在其余谱段图像上相同位置的图像块，并将这些图像块转化为列向量，构建相应谱段的相似矩阵；
(214)将步骤(212)和步骤(213)获得的M个相似矩阵堆叠成第p个三阶相似张量
Figure FDA0001984312410000023
(215)对于步骤(211)中提取的P个图像块，重复(212)-(214)，便获得P个三阶相似张量
Figure FDA0001984312410000024
3.根据权利要求2所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于在参考谱段图像上选取的P个图像块的重叠范围为40％～80％。
4.根据权利要求2所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于采用块匹配或欧拉距离在参考谱段图像搜索窗内找到第p个图像块的z-1个相似的图像块。
5.根据权利要求2所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于所述搜索窗窗口大小为选取的图像块大小的2～10倍。
6.根据权利要求1所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于步骤(22)中，
Figure FDA0001984312410000025
的最优解
Figure FDA0001984312410000026
通过构建的以下函数得到：
Figure FDA0001984312410000027
式中，
Figure FDA0001984312410000028
和
Figure FDA0001984312410000029
由
Figure FDA00019843124100000210
沿第k阶展开得到的矩阵经奇异值分解得到的第l大的奇异值和奇异向量，foldk是将矩阵还原成张量的算子；得到的
Figure FDA00019843124100000211
的最优解
Figure FDA00019843124100000212
即为
Figure FDA00019843124100000213
7.根据权利要求1所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于步骤(23)中，将具有M个能谱段的稀疏分量
Figure FDA0001984312410000031
沿第三阶分离得到M个矩阵切片
Figure FDA0001984312410000032
Figure FDA0001984312410000033
是第m个矩阵切片，每个切片
Figure FDA0001984312410000034
采用Chambolle投影算法求解，即：
Figure FDA0001984312410000035
式中μ＝λ2/ρ，div是散度算子，Ym是张量
Figure FDA0001984312410000036
沿第三阶分离得到的第m个切片，矩阵q通过迭代得到，其第t+1次迭代过程为
Figure FDA0001984312410000037
其中i和j分别是式中矩阵行和列的下标，τ为设定的迭代步长，t＝0，1，2，…，T，T为迭代次数，q0＝0；
将得到的
Figure FDA0001984312410000038
堆叠所得张量即为
Figure FDA0001984312410000039
8.根据权利要求1所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于步骤(24)中，对
Figure FDA00019843124100000310
采用共轭梯度法进行求解，得到重建能谱CT图像的输出值
Figure FDA00019843124100000311
9.根据权利要求1所述基于空谱双域张量自相似的能谱CT重建方法，其特征在于根据相邻两次获取的重建能谱CT图像的相对均方误差来对重建能谱CT图像的质量进行评价。
CN201910160075.9A 2019-03-04 2019-03-04 基于空谱双域张量自相似的能谱ct重建方法 Active CN109903355B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910160075.9A CN109903355B (zh) 2019-03-04 2019-03-04 基于空谱双域张量自相似的能谱ct重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910160075.9A CN109903355B (zh) 2019-03-04 2019-03-04 基于空谱双域张量自相似的能谱ct重建方法

Publications (2)

Publication Number Publication Date
CN109903355A CN109903355A (zh) 2019-06-18
CN109903355B true CN109903355B (zh) 2020-08-04

Family

ID=66946210

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910160075.9A Active CN109903355B (zh) 2019-03-04 2019-03-04 基于空谱双域张量自相似的能谱ct重建方法

Country Status (1)

Country Link
CN (1) CN109903355B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110428395B (zh) * 2019-06-20 2021-10-08 浙江大学 单能谱ct图像的多材料分解方法
CN116012264B (zh) * 2023-03-27 2023-06-13 山东省工业技术研究院 一种基于稀疏约束的图像恢复方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103907132A (zh) * 2011-11-03 2014-07-02 皇家飞利浦有限公司 图像数据处理
CN104103086A (zh) * 2014-06-06 2014-10-15 华南理工大学 一种稀疏采样角度下基于变分不等式的ct图像重建方法
CN104751429A (zh) * 2015-01-27 2015-07-01 南方医科大学 一种基于字典学习的低剂量能谱ct图像处理方法
WO2017088816A1 (zh) * 2015-11-27 2017-06-01 广州聚普科技有限公司 一种基于dti的颅内神经纤维束的三维重建方法
CN108010098A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种双能谱ct基材料图像迭代重建方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103907132A (zh) * 2011-11-03 2014-07-02 皇家飞利浦有限公司 图像数据处理
CN104103086A (zh) * 2014-06-06 2014-10-15 华南理工大学 一种稀疏采样角度下基于变分不等式的ct图像重建方法
CN104751429A (zh) * 2015-01-27 2015-07-01 南方医科大学 一种基于字典学习的低剂量能谱ct图像处理方法
WO2017088816A1 (zh) * 2015-11-27 2017-06-01 广州聚普科技有限公司 一种基于dti的颅内神经纤维束的三维重建方法
CN108010098A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种双能谱ct基材料图像迭代重建方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Tensor decomposition and nonlocal means based spectral CT reconstruction;Zhang Yanbo 等;《PROCEEDINGS OF SPIE》;20161231;第99670Z页 *
基于分数阶自适应正则化的统计迭代重建方法;张意 等;《第十四届中国体视学与图像分析学术会议论文集》;20150930;第113-117页 *
基于图像总变分和张量字典的多能谱CT材料识别研究;陈佩君 等;《光学学报》;20181130;第1111002页 *

Also Published As

Publication number Publication date
CN109903355A (zh) 2019-06-18

Similar Documents

Publication Publication Date Title
Yin et al. Domain progressive 3D residual convolution network to improve low-dose CT imaging
Hu et al. Hybrid-domain neural network processing for sparse-view CT reconstruction
Kang et al. Deep convolutional framelet denosing for low-dose CT via wavelet residual network
Wu et al. Low-dose spectral CT reconstruction using image gradient ℓ0–norm and tensor dictionary
Huang et al. CaGAN: A cycle-consistent generative adversarial network with attention for low-dose CT imaging
Wang et al. FBP-Net for direct reconstruction of dynamic PET images
Zhang et al. Tensor-based dictionary learning for spectral CT reconstruction
Yuan et al. SIPID: A deep learning framework for sinogram interpolation and image denoising in low-dose CT reconstruction
CN110415307B (zh) 一种基于张量补全的多能ct成像方法、装置及其存储设备
Xue et al. LCPR-Net: low-count PET image reconstruction using the domain transform and cycle-consistent generative adversarial networks
Cheslerean-Boghiu et al. WNet: A data-driven dual-domain denoising model for sparse-view computed tomography with a trainable reconstruction layer
CN109903355B (zh) 基于空谱双域张量自相似的能谱ct重建方法
Zhang et al. Spectral CT image-domain material decomposition via sparsity residual prior and dictionary learning
Liu et al. Deep residual constrained reconstruction via learned convolutional sparse coding for low-dose CT imaging
CN112656438B (zh) 一种基于曲面全变差的低剂量ct投影域去噪及重建方法
CN109658464B (zh) 基于加权核范数极小的稀疏角ct图像重建方法
Zhang et al. Deep generalized learning model for PET image reconstruction
Yang et al. Improve 3D cone-beam CT reconstruction by slice-wise deep learning
Wu et al. Deep learning-based low-dose tomography reconstruction with hybrid-dose measurements
Hu et al. CROSS: Cross-domain residual-optimization-based structure strengthening reconstruction for limited-angle CT
Kim et al. CNN-based CT denoising with an accurate image domain noise insertion technique
Wang et al. Hybrid-Domain Integrative Transformer Iterative Network for Spectral CT Imaging
Wang et al. Helical ct reconstruction from sparse-view data through exploiting the 3d anatomical structure sparsity
Chen et al. Fourth-order nonlocal tensor decomposition model for spectral computed tomography
CN113920216A (zh) 基于张量核范数和变换Lp范数的能谱CT重建方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant