CN108280859B - 一种采样角度受限下的ct稀疏投影图像重建方法及装置 - Google Patents

一种采样角度受限下的ct稀疏投影图像重建方法及装置 Download PDF

Info

Publication number
CN108280859B
CN108280859B CN201711420601.8A CN201711420601A CN108280859B CN 108280859 B CN108280859 B CN 108280859B CN 201711420601 A CN201711420601 A CN 201711420601A CN 108280859 B CN108280859 B CN 108280859B
Authority
CN
China
Prior art keywords
solution
random
optimal solution
projection
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
CN201711420601.8A
Other languages
English (en)
Other versions
CN108280859A (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.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
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 South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN201711420601.8A priority Critical patent/CN108280859B/zh
Publication of CN108280859A publication Critical patent/CN108280859A/zh
Application granted granted Critical
Publication of CN108280859B publication Critical patent/CN108280859B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种采样角度受限下的CT稀疏投影图像重建方法及装置。所述方法包括根据投影数据，获得投影方程的解的伪逆矩阵，根据伪逆矩阵，生成随机解集，对当前随机解集中的各个解进行相应的保留或替换，迭代次数达预设最大值时从当前随机解集中选择出最优解作为所需获得的重建结果等步骤。所述装置包括用于存储至少一个程序的存储器和用于加载至少一个程序以执行本发明方法的处理器。本发明通过求离散化投影重建方程的伪逆作为算法的初始值，保证了初始解的质量，通过随机游走生成一组解并分别进行迭代优化，保证了优化路径的多样性，从而克服了传统重建方法由初始解和和迭代路径的缺陷带来的困难。本发明应用于图像处理技术领域。

Description

一种采样角度受限下的CT稀疏投影图像重建方法及装置
技术领域
本发明涉及图像处理技术领域，尤其是一种采样角度受限下的CT稀疏投影图像重建方法及装置。
背景技术
术语解释：
CT：Computed Tomography(计算机断层成像)
TV：Total Variation(全变分)
POCS-TV：Projection On Convex Sets-Total Variation(凸集投影-全变分最小化算法)
CT成像是通过利用X射线穿透待检测物体时发生的强度衰减信息，重建得到物体的断层图像。而提高角度受限下稀疏投影数据的重建质量是CT重建领域十分重要且亟需解决的难题。如医学中，在对人体大脑进行扫描时，过大的射线照射剂量可能会对大脑造成二次伤害，因此需要减小CT射线的照射剂量，只能得到角度受限的稀疏投影数据；在工业中，为了对大型的不规则印制电路板或在役的管道进行无损探伤，受CT成像系统几何位置限制，只能在有限角度下对待重建物体进行扫描，无法得到完整角度投影数据。在采样角度受限情况下，因缺失连续角度范围的投影数据，不满足解析法对获得精确重建结果的数据完整性要求(Tuy-Smith条件)，直接采用解析法重建存在严重伪影。针对角度受限投影数据的重建，现有技术通常在TV约束的正则化框架下利用基于凸优化的迭代法进行求解。现有技术在迭代重建法的原目标函数
Figure BDA0001522874350000013
(其中A为投影系数矩阵，x为待重建图像的列向量形式，b为投影数据)基础上，根据CT图像的分片光滑特性，添加待重建图像的TV项作为正则化约束，得到目标函数
Figure BDA0001522874350000012
然后通过基于凸优化的迭代法求解得到重建结果。
但是，由于采样角度受限下的稀疏投影数据对称性、完备性极低，重建算法的搜索空间巨大，传统的基于凸优化的迭代重建方法从单一的初始解出发，按照目标函数在当前解的负梯度方向路径进行搜索，一方面重建质量受初始解质量的严重影响，另一方面，单一的固定迭代路径使得优化过程无法在有限时间内遍历巨大的搜索空间，易陷入局部最优且无法跳离，重建图像质量差。由此可见，现有的重建方法面临的困难，主要是由初始解和迭代路径引起的。
发明内容
为了解决上述技术问题，本发明的第一目的在于提供一种采样角度受限下的CT稀疏投影图像重建方法，第二目的在于提供一种采样角度受限下的CT稀疏投影图像重建装置。
本发明所采取的第一技术方案是：
一种采样角度受限下的CT稀疏投影图像重建方法，包括：
S1.获取投影方程的解的伪逆矩阵，所述投影方程根据投影数据而建立；
S2.根据伪逆矩阵，生成与本轮迭代相对应的随机解集；
S3.若当前为第一轮迭代过程，则返回步骤S2开始执行，反之，将当前随机解集与上一轮迭代过程的随机解集进行对应解的比较，根据比较结果，对当前随机解集中的各个解进行相应的保留或替换；
S4.判断当前迭代次数是否已达预设最大值，若是，则根据适应度评价，从当前随机解集中选择出最优解，然后执行步骤S5，反之，返回从步骤S2开始执行；
S5.根据最优解重建得到CT图像。
进一步地，所述步骤S2具体为：根据伪逆矩阵，通过以下公式获得多个随机解，将获得的所有解构成与该轮迭代相对应的随机解集：
Figure BDA0001522874350000021
式中，t为迭代轮次，i表示序号，N表示自然数集，
Figure BDA0001522874350000022
为第t轮迭代所得的随机解集中第i个解，xini为投影方程的解的伪逆矩阵，ε是一维度与xini相同的服从(0,σ)分布的随机变量。
进一步地，步骤S3具体包括：
S31.对当前随机解集中的每一个解进行有效性约束；
S32.计算当前随机解集中每一个解的后验概率密度；
S33.根据当前随机解集中每一个解的后验概率密度与上一轮迭代过程的随机解集中每一个解的后验概率密度，计算梅特罗波利斯比率；
S34.判断梅特罗波利斯比率是否大于第一阈值，若是，则保留当前随机解集中的解，反之，将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解。
进一步地，步骤S31具体为：利用下式，对当前随机解集中的每一个解进行有效性约束：
Figure BDA0001522874350000031
式中，
Figure BDA0001522874350000032
为
Figure BDA0001522874350000033
的第j个元素；
步骤S32具体为：利用下式，计算当前随机解集中每一个解的后验概率密度：
Figure BDA0001522874350000034
式中，
Figure BDA0001522874350000035
为
Figure BDA0001522874350000036
的后验概率密度，C为
Figure BDA0001522874350000037
的协方差矩阵；
步骤S33具体为：利用下式，计算梅特罗波利斯比率：
Figure BDA0001522874350000038
式中，α为梅特罗波利斯比率；
步骤S34具体为：利用下式，对梅特罗波利斯比率是否大于第一阈值进行判断，从而保留当前随机解集中的解，或者将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解：
Figure BDA0001522874350000039
式中，rand()为一值域为[0,1]的随机函数，rand()的值为第一阈值。
进一步地，步骤S4中，选择最优解所用的公式为：
Figure BDA00015228743500000310
式中，
Figure BDA00015228743500000311
为所要获取的最优解，F(·)为适应度评价函数，
Figure BDA00015228743500000312
A为投影系数矩阵，b为投影数据，
Figure BDA00015228743500000313
为
Figure BDA00015228743500000314
的TV正则化约束项。
进一步地，步骤S4与步骤S5之间还包括以下步骤：
S4A.利用基于凸优化的迭代重建算法，对最优解进行更新。
进一步地，步骤S4A具体包括：
S4A1.采用ART迭代算法对最优解进行处理，从而得到一致最优解；
S4A2.对一致最优解进行图像非负性约束，从而得到非负最优解；
S4A3.根据非负最优解的TV项的梯度，以非负最优解作为迭代变量进行迭代，将得到的迭代结果作为更新后的最优解。
进一步地，所述步骤S4A1具体为：使用以下公式对最优解进行处理，从而得到一致最优解：
Figure BDA0001522874350000041
式中，M为投影数据总数，bm为第m个投影数据，Am为投影系数矩阵A的第m行，根据公式求得的fM为一致最优解；
所述步骤S4A2具体为：使用以下公式对一致最优解进行图像非负性约束，从而得到非负最优解：
Figure BDA0001522874350000046
式中，g0(j)为g0的第j个元素，根据公式求得的g0为非负最优解。
进一步地，所述步骤S4A3具体为：使用以下公式进行迭代，将得到的迭代结果作为更新后的最优解：
Figure BDA0001522874350000042
式中，Ngrad为梯度下降迭代总次数，
Figure BDA0001522874350000043
Figure BDA0001522874350000044
表示gn中坐标为(p,q)的点的像素值，
根据公式计算求得的
Figure BDA0001522874350000045
为更新后的最优解。
本发明所采取的第二技术方案是：
一种采样角度受限下的CT稀疏投影图像重建装置，包括：
存储器，用于存储至少一个程序；
处理器，用于加载所述至少一个程序以执行第一技术方案所述一种采样角度受限下的CT稀疏投影图像重建方法。
本发明的有益效果是：本发明通过求离散化投影重建方程的伪逆作为算法的初始值，保证了初始解的质量，通过随机游走生成一组解并分别进行迭代优化，保证了优化路径的多样性，从而克服了传统重建方法由初始解和迭代路径的缺陷带来的困难。
进一步地，本发明还利用高斯马尔科夫模型对解的进化能力进行预测，基于解的进化潜能对其进行随机更新，使得算法能够在提高全局搜索能力的同时保证效率；对基于适应度评价下的最优解使用基于凸优化的迭代更新，并将更新后的最优解作为最终的重建结果，进一步提高了包含随机因素下算法的效率。本发明全局搜索能力强，重建图像不存在严重伪影，显著提高了重建图像的质量，降低重建误差，能够在数据缺失严重的情况下实现图像重建，特别在处理含噪数据时表现出很强的鲁棒性和适应性。
附图说明
图1为本发明方法的流程图；
图2为Shepp-Logan模型原图；
图3为投影角度间隔为3°时应用POCS-TV算法的重建效果图；
图4为投影角度间隔为3°时应用本发明方法的重建效果图；
图5为投影角度间隔为3°时重建结果第64行的像素值曲线图；
图6为投影角度间隔为3°时重建结果第64列的像素值曲线图；
图7为投影角度间隔为6°时应用POCS-TV算法的重建效果图；
图8为投影角度间隔为6°时应用本发明方法的重建效果图；
图9为投影角度间隔为6°时重建结果第64行的像素值曲线图；
图10为投影角度间隔为6°时重建结果第64列的像素值曲线图；
图11为本发明装置的结构框图。
具体实施方式
实施例1
一种采样角度受限下的CT稀疏投影图像重建方法，如图1所示，包括：
S1.获取投影方程的解的伪逆矩阵，所述投影方程根据投影数据而建立；
S2.根据伪逆矩阵，生成与本轮迭代相对应的随机解集；
S3.若当前为第一轮迭代过程，则返回步骤S2开始执行，反之，将当前随机解集与上一轮迭代过程的随机解集进行对应解的比较，根据比较结果，对当前随机解集中的各个解进行相应的保留或替换；
S4.判断当前迭代次数是否已达预设最大值，若是，则根据适应度评价，从当前随机解集中选择出最优解，然后执行步骤S5，反之，返回从步骤S2开始执行；
S5.根据最优解重建得到CT图像。
在应用本发明方法之前，可以使用各种方法采集角度受限下的稀疏投影数据。所述角度受限下的稀疏投影数据，可以是投影角度范围为(0°,90°)，投影间隔大于等于3°、小于等于6°所得的投影数据。根据投影数据重建图像，实质是要求解离散化的投影方程Ax＝b，其中A为投影系数矩阵，x为待重建图像的列向量形式，b为投影数据。求得方程Ax＝b的解x，便能根据x重建图像。
x的伪逆矩阵可以使用以下公式xini＝(ATA)-1ATb求得。由于在本发明方法中，x的伪逆矩阵用于作为后续一系列迭代运算的初值，因此下标ini；本发明方法中只使用到了x的伪逆矩阵，因此本发明方法中“伪逆矩阵”专指x的伪逆矩阵。
根据伪逆矩阵，可以使用随机游走等算法产生一系列矩阵，由于这些矩阵都有可能是方程Ax＝b的解或者解的近似，对这些矩阵使用迭代算法，将收敛得到方程Ax＝b的解，因此将一系列矩阵的集合作为随机解集。
在以下方法中，需要使用到迭代算法或者循环的处理步骤，在一轮迭代过程当中，可能会对随机解集或者随机解集当中特定的解做处理，但在此轮迭代过程当中，处理前和处理后的随机解集都称作当前随机解集。
在执行第一轮迭代的情况下，由于不存在上一轮迭代，因此步骤S3将无法直接执行。可以在执行S1和S2从而得到一个随机解集之后，再次执行S2，在执行后一次S2时便进入第二轮迭代，此时便可以执行步骤S3。预设的迭代次数最大值可以为1，这样，可以在执行S1和S2后，跳过S3，然后执行S4中根据适应度评价，从当前随机解集中选择出最优解的步骤，然后跳到步骤S5执行。
步骤S3中，将本轮迭代的随机解集，也就是当前随机解集，与上一轮迭代过程的随机解集一起针对各对对应解分别进行比较。对应解指在两个解集中索引号、序号、编号相同或者地位相同的两个解。经过各个对应解的比较，按照进化能力的标准，在当前随机解集中的解保留下来，或者将当前随机解集中的部分或者全部解替换成上一轮迭代过程的随机解集中的对应解，所得的随机解集仍称作当前随机解集。上述操作使得当前随机解集中的各个解具有更佳的进化能力。
步骤S4中判断当前迭代次数是否已达预设最大值，如果没有，那么就返回步骤S2重新执行，使得当前随机解集不断完善。如果已达预设最大值，那么就使用适应度评价的方法，根据对当前随机解集的适应度评价结果，在当前随机解集中选择出一个最优解作为方程Ax＝b的解，也就是所需获得的重建结果。
进一步作为优选的实施方式，所述步骤S2具体为：根据伪逆矩阵，通过以下公式获得多个随机解，将获得的所有解构成与该轮迭代相对应的随机解集：
Figure BDA0001522874350000071
式中，t为迭代轮次，i为解的编号，
Figure BDA0001522874350000072
为第t轮迭代所得的随机解集中第i个解，xini为投影方程的解的伪逆矩阵，ε是一维度与xini相同的服从(0,σ)分布的随机变量。
由于ε是随机变量，每一次执行S2，ε的取值都是随机的，其与xini叠加后所得的结果也是随机的，由此每一次执行S2都可以得到一个新的随机解集。可以根据实际计算的需要，确定随机解集中解的总个数。
下面使用高斯马尔科夫模型对解的进化能力进行定量评价，然后基于解的进化能力对当前随机解集进行更新。
进一步作为优选的实施方式，步骤S3具体为：
S31.对当前随机解集中的每一个解进行有效性约束；
S32.计算当前随机解集中每一个解的后验概率密度；
S33.根据当前随机解集中每一个解的后验概率密度与上一轮迭代过程的随机解集中每一个解的后验概率密度，计算梅特罗波利斯比率；
S34.判断梅特罗波利斯比率是否大于第一阈值，若是，则保留当前随机解集中的解，反之，将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解。
步骤S31-S33可以计算出梅特罗波利斯(Metropolis)比率，梅特罗波利斯比率可以反映出当前随机解集与上一轮迭代过程的随机解集的进化能力孰大孰小。根据实际需要，确定第一阈值，即梅特罗波利斯比率的对比值，如果梅特罗波利斯比率小于或等于第一阈值，表明当前随机解集的进化能力不如上一轮迭代过程的随机解集，那么就将当前随机解集中的部分解或者全部解替换成上一轮迭代过程的随机解集的对应解；如果梅特罗波利斯比率大于第一阈值，表明当前随机解集的进化能力比上一轮迭代过程的随机解集好，可以将当前随机解集中的部分解或全部解保留下来。经过这样处理后的当前随机解集仍称作当前随机解集。反复执行步骤S34，直至对当前随机解集中的所有解进行保留或替换。
进一步作为优选的实施方式，步骤S31具体为：利用下式，对当前随机解集中的每一个解进行有效性约束：
Figure BDA0001522874350000081
式中，
Figure BDA0001522874350000082
为
Figure BDA0001522874350000083
的第j个元素。
由公式可知，有效性约束可以将当前随机解集中的各个解的各元素都约束在(0,1)的范围内。即使当前随机解集中各个解是经过随机过程得到的，通过有效性约束，也可以保证各个解中每个元素的值对应于图像像素值都是有意义的。
步骤S32具体为：利用下式，计算当前随机解集中每一个解的后验概率密度：
Figure BDA0001522874350000084
式中，
Figure BDA0001522874350000085
为
Figure BDA0001522874350000086
的后验概率密度，C为
Figure BDA0001522874350000087
的协方差矩阵。
步骤S33具体为：利用下式，计算梅特罗波利斯比率：
Figure BDA0001522874350000088
式中，α为梅特罗波利斯比率。
步骤S34具体为：利用下式，对梅特罗波利斯比率是否大于第一阈值进行判断，从而保留当前随机解集中的解，或者将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解：
Figure BDA0001522874350000089
式中，rand()为一值域为[0,1]的随机函数，rand()的值为第一阈值。
在步骤S34中，可以分别对每对对应解进行判断，而且第一阈值不是一个固定不变的值，由于其是随机函数，其对于不同的对应解可以取不同的值。由于梅特罗波利斯比率进行对比的阈值不同，也就是对每对对应解的进化能力的评价标准是不同的，于是可能出现当前随机解集中部分解被保留、部分解被替换的情况。而且，由于阈值的随机性，还可能出现当
Figure BDA00015228743500000810
小于
Figure BDA00015228743500000811
也就是
Figure BDA00015228743500000812
的进化能力不及
Figure BDA00015228743500000813
时，
Figure BDA00015228743500000814
仍有可能被保留下来的情况，这样，本发明方法可以对更广阔的解空间进行探索，提高全局搜索能力，可以缓解陷入局部最优的情况。
进一步作为优选的实施方式，步骤S4中，选择最优解所用的公式为：
Figure BDA0001522874350000091
式中，
Figure BDA0001522874350000092
为所要获取的最优解，F(·)为适应度评价函数，
Figure BDA0001522874350000093
A为投影系数矩阵，b为投影数据，
Figure BDA0001522874350000094
表示
Figure BDA0001522874350000095
的TV正则化约束项。
适应度评价函数中，
Figure BDA0001522874350000096
为数据保真项，评价重建图像与投影数据的一致性，
Figure BDA0001522874350000097
为TV正则化约束项，评价重建图像的分片光滑性，适应度评价函数值越小，表明解的质量越好。
在经过预定轮次的迭代后，便可以从当前随机解集中选择出最优解作为所需获得的重建结果。优选地，选择最优解的标准是，使得适应度评价函数取得最小值时的xi t的值。
进一步作为优选的实施方式，步骤S4与步骤S5之间还包括以下步骤：
S4A.利用基于凸优化的迭代重建算法，对最优解进行更新。
步骤S4中得到的最优解，可以直接执行步骤S5，也就是直接用于作为所需获得的重建结果。为了取得更好的效果，可以在执行步骤S4得到最优解后，接着执行步骤S4A，对最优解进行更新，然后执行步骤S5，根据经过步骤S4A处理的最优解，也就是更新后的最优解，来重建得到CT图像。
进一步作为优选的实施方式，步骤S4A具体包括：
S4A1.采用ART迭代算法对最优解进行处理，从而得到一致最优解；
S4A2.对一致最优解进行图像非负性约束，从而得到非负最优解；
S4A3.根据非负最优解的TV项的梯度，以非负最优解作为迭代变量进行迭代，将得到的迭代结果作为更新后的最优解。
进一步作为优选的实施方式，所述步骤S4A1具体为：使用以下公式对最优解进行处理，从而得到一致最优解：
Figure BDA0001522874350000098
式中，M为投影数据总数，bm为第m个投影数据，Am为投影系数矩阵A的第m行，根据公式求得的fM为一致最优解。
上述公式是ART迭代算法的一种实现方法，最终得到的结果fM为一致最优解。
进一步作为优选的实施方式，所述步骤S4A2具体为：使用以下公式对一致最优解进行图像非负性约束，从而得到非负最优解：
Figure BDA0001522874350000101
式中，g0(j)为g0的第j个元素，根据公式求得的g0为非负最优解。
进一步作为优选的实施方式，所述步骤S4A3具体为：使用以下公式进行迭代，将得到的迭代结果作为更新后的最优解：
Figure BDA0001522874350000102
式中，Ngrad为梯度下降迭代总次数，
Figure BDA0001522874350000103
Figure BDA0001522874350000104
表示gn中坐标为(p,q)的点的像素值，根据公式计算求得的
Figure BDA0001522874350000105
为更新后的最优解。经过步骤S4A对最优解使用基于凸优化的更新处理，可以保证包含随机因素下算法的效率。
实施例2
在本实施例中，使用实施例1的方法，对二维Shepp-Logan模型(128×128)的模拟投影数据进行重建，并在相同的迭代次数下与POCS-TV算法进行比较。假设模拟投影数据获取及重建过程的参数如下：
(1)X射线源与待重建物体中心的距离为256mm；
(2)X射线检测器与待重建物体中心的距离为256mm；
(3)X射线检测器单元的个数为256，每单元的宽度为0.5mm；
(4)投影角度范围为(0°,90°)，投影角度间隔为3°和6°；
(5)为说明本发明在含噪数据的重建上具有更强的鲁棒性和适应性，在投影数据上添加0均值，0.15方差的高斯噪声；
(6)设置总的迭代次数为5000次。
Shepp-Logan模型原图如图2所示。
投影角度间隔为3°时，用POCS-TV算法的重建结果如图3所示，用本发明方法的重建结果如图4所示，重建结果第64行的像素值曲线如图5所示，第64列的像素值曲线如图6所示。
投影角度间隔为6°时，用POCS-TV算法的重建结果如图7所示，用本发明方法的重建结果如图8所示，重建结果第64行的像素值曲线如图9所示，第64列的像素值曲线如图10所示。
从图3至图10可以看出，本发明方法的重建结果轮廓完整，细节更加清晰，更接近于模型原图，重建质量远胜于POCS-TV算法。
实施例3
一种采样角度受限下的CT稀疏投影图像重建装置，如图11所示，包括：
存储器，用于存储至少一个程序；
处理器，用于加载所述至少一个程序以执行实施例1和实施例2所述一种采样角度受限下的CT稀疏投影图像重建方法。
以上是对本发明的较佳实施进行了具体说明，但对本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做作出种种的等同变形或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。

Claims (6)

1.一种采样角度受限下的CT稀疏投影图像重建方法，其特征在于，包括：
S1.获取投影方程的解的伪逆矩阵，所述投影方程根据投影数据而建立；
S2.根据伪逆矩阵，生成与本轮迭代相对应的随机解集；
S3.若当前为第一轮迭代过程，则返回步骤S2开始执行，反之，将当前随机解集与上一轮迭代过程的随机解集进行对应解的比较，根据比较结果，对当前随机解集中的各个解进行相应的保留或替换；
S4.判断当前迭代次数是否已达预设最大值，若是，则根据适应度评价，从当前随机解集中选择出最优解，然后执行步骤S5，反之，返回从步骤S2开始执行；
S5.根据最优解重建得到CT图像；
所述步骤S2具体为：根据伪逆矩阵，通过以下公式获得多个随机解，将获得的所有解构成与该轮迭代相对应的随机解集：
Figure FDA0002907407820000011
式中，t为迭代轮次，i表示序号，N表示自然数集，
Figure FDA0002907407820000012
为第t轮迭代所得的随机解集中第i个解，xini为投影方程的解的伪逆矩阵，ε是一维度与xini相同的服从(0,σ)分布的随机变量；
步骤S3具体包括：
S31.对当前随机解集中的每一个解进行有效性约束；
S32.计算当前随机解集中每一个解的后验概率密度；
S33.根据当前随机解集中每一个解的后验概率密度与上一轮迭代过程的随机解集中每一个解的后验概率密度，计算梅特罗波利斯比率；
S34.判断梅特罗波利斯比率是否大于第一阈值，若是，则保留当前随机解集中的解，反之，将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解；
步骤S31具体为：利用下式，对当前随机解集中的每一个解进行有效性约束：
Figure FDA0002907407820000013
式中，
Figure FDA0002907407820000014
为
Figure FDA0002907407820000015
的第j个元素；
步骤S32具体为：利用下式，计算当前随机解集中每一个解的后验概率密度：
Figure FDA0002907407820000021
式中，
Figure FDA0002907407820000022
为
Figure FDA0002907407820000023
的后验概率密度，C为
Figure FDA0002907407820000024
的协方差矩阵；
步骤S33具体为：利用下式，计算梅特罗波利斯比率：
Figure FDA0002907407820000025
式中，α为梅特罗波利斯比率；
步骤S34具体为：利用下式，对梅特罗波利斯比率是否大于第一阈值进行判断，从而保留当前随机解集中的解，或者将当前随机解集中的解替换成上一轮迭代过程的随机解集中的对应解：
Figure FDA0002907407820000026
式中，rand()为一值域为[0,1]的随机函数，rand()的值为第一阈值；
步骤S4中，选择最优解所用的公式为：
Figure FDA0002907407820000027
式中，
Figure FDA0002907407820000028
为所要获取的最优解，F(·)为适应度评价函数，
Figure FDA0002907407820000029
A为投影系数矩阵，b为投影数据，
Figure FDA00029074078200000210
为
Figure FDA00029074078200000211
的TV正则化约束项。
2.根据权利要求1所述的一种采样角度受限下的CT稀疏投影图像重建方法，其特征在于，步骤S4与步骤S5之间还包括以下步骤：
S4A.利用基于凸优化的迭代重建算法，对最优解进行更新。
3.根据权利要求2所述的一种采样角度受限下的CT稀疏投影图像重建方法，其特征在于，步骤S4A具体包括：
S4A1.采用ART迭代算法对最优解进行处理，从而得到一致最优解；
S4A2.对一致最优解进行图像非负性约束，从而得到非负最优解；
S4A3.根据非负最优解的TV项的梯度，以非负最优解作为迭代变量进行迭代，将得到的迭代结果作为更新后的最优解。
4.根据权利要求3所述的一种采样角度受限下的CT稀疏投影图像重建方法，其特征在于：
所述步骤S4A1具体为：使用以下公式对最优解进行处理，从而得到一致最优解：
Figure FDA0002907407820000031
式中，M为投影数据总数，bm为第m个投影数据，Am为投影系数矩阵A的第m行，根据公式求得的fM为一致最优解；
所述步骤S4A2具体为：使用以下公式对一致最优解进行图像非负性约束，从而得到非负最优解：
Figure FDA0002907407820000032
式中，g0(j)为g0的第j个元素，根据公式求得的g0为非负最优解。
5.根据权利要求4所述的一种采样角度受限下的CT稀疏投影图像重建方法，其特征在于，所述步骤S4A3具体为：使用以下公式进行迭代，将得到的迭代结果作为更新后的最优解：
Figure FDA0002907407820000033
式中，Ngrad为梯度下降迭代总次数，
Figure FDA0002907407820000034
Figure FDA0002907407820000035
表示gn中坐标为(p,q)的点的像素值，
根据公式计算求得的
Figure FDA0002907407820000036
为更新后的最优解。
6.一种采样角度受限下的CT稀疏投影图像重建装置，其特征在于，包括：
存储器，用于存储至少一个程序；
处理器，用于加载所述至少一个程序以执行权利要求1-5任一项所述一种采样角度受限下的CT稀疏投影图像重建方法。
CN201711420601.8A 2017-12-25 2017-12-25 一种采样角度受限下的ct稀疏投影图像重建方法及装置 Active CN108280859B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711420601.8A CN108280859B (zh) 2017-12-25 2017-12-25 一种采样角度受限下的ct稀疏投影图像重建方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711420601.8A CN108280859B (zh) 2017-12-25 2017-12-25 一种采样角度受限下的ct稀疏投影图像重建方法及装置

Publications (2)

Publication Number Publication Date
CN108280859A CN108280859A (zh) 2018-07-13
CN108280859B true CN108280859B (zh) 2021-03-30

Family

ID=62802211

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711420601.8A Active CN108280859B (zh) 2017-12-25 2017-12-25 一种采样角度受限下的ct稀疏投影图像重建方法及装置

Country Status (1)

Country Link
CN (1) CN108280859B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109840927B (zh) * 2019-01-24 2020-11-10 浙江大学 一种基于各向异性全变分的有限角度ct重建算法
CN109903353B (zh) * 2019-01-28 2023-02-14 华南理工大学 一种迭代演化模型的ct图像稀疏重建方法
CN111956180B (zh) * 2019-05-20 2023-06-27 华北电力大学（保定） 一种重建光声内窥层析图像的方法
CN110210443B (zh) * 2019-06-11 2022-03-15 西北工业大学 一种优化投影对称性近似稀疏分类的手势识别方法
CN112288832B (zh) * 2020-12-24 2021-03-23 中国人民解放军国防科技大学 一种角度受限及稀疏采样的层析成像图像重构方法
CN115619890B (zh) * 2022-12-05 2023-04-07 山东省计算中心（国家超级计算济南中心） 基于并行随机迭代求解线性方程组的断层成像方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103824295A (zh) * 2014-03-03 2014-05-28 天津医科大学 一种肺部ct图像中粘连血管型肺结节的分割方法
CN104091355A (zh) * 2014-06-06 2014-10-08 华南理工大学 一种采样角度受限下的ct图像重建方法
CN104103086A (zh) * 2014-06-06 2014-10-15 华南理工大学 一种稀疏采样角度下基于变分不等式的ct图像重建方法
CN105118078A (zh) * 2015-09-24 2015-12-02 中国科学院苏州生物医学工程技术研究所 欠采样的ct图像重建方法
CN106163403A (zh) * 2013-12-18 2016-11-23 伊利克塔股份有限公司 Ct图像中的目标特定剂量和散射估计
CN107122725A (zh) * 2017-04-18 2017-09-01 深圳大学 一种基于联合稀疏判别分析的人脸识别方法及其系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160242721A1 (en) * 2015-02-20 2016-08-25 Kabushiki Kaisha Toshiba Apparatus and method for fast iterative reconstruction in computed tomography

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106163403A (zh) * 2013-12-18 2016-11-23 伊利克塔股份有限公司 Ct图像中的目标特定剂量和散射估计
CN103824295A (zh) * 2014-03-03 2014-05-28 天津医科大学 一种肺部ct图像中粘连血管型肺结节的分割方法
CN104091355A (zh) * 2014-06-06 2014-10-08 华南理工大学 一种采样角度受限下的ct图像重建方法
CN104103086A (zh) * 2014-06-06 2014-10-15 华南理工大学 一种稀疏采样角度下基于变分不等式的ct图像重建方法
CN105118078A (zh) * 2015-09-24 2015-12-02 中国科学院苏州生物医学工程技术研究所 欠采样的ct图像重建方法
CN107122725A (zh) * 2017-04-18 2017-09-01 深圳大学 一种基于联合稀疏判别分析的人脸识别方法及其系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于稀疏优化的计算机断层成像图像不完全角度重建综述;王林元，等;《物理学报》;20141231;第208702-1至208702-10页 *
基于缺陷投影的稀疏角度数据CT图像迭代重建仿真研究;陈海，等;《医疗卫生装备》;20170430;第29-31、37页 *

Also Published As

Publication number Publication date
CN108280859A (zh) 2018-07-13

Similar Documents

Publication Publication Date Title
CN108280859B (zh) 一种采样角度受限下的ct稀疏投影图像重建方法及装置
Anirudh et al. Lose the views: Limited angle CT reconstruction via implicit sinogram completion
US11769277B2 (en) Deep learning based scatter correction
Zhuge et al. TVR-DART: A more robust algorithm for discrete tomography from limited projection data with automated gray value estimation
CN110298447B (zh) 用于处理机器学习方法的参数的方法以及重建方法
CN109472841B (zh) 基于混合高斯/泊松最大似然函数的cbct三维重建方法
Xu et al. Efficient low‐dose CT artifact mitigation using an artifact‐matched prior scan
Ertas et al. An iterative tomosynthesis reconstruction using total variation combined with non-local means filtering
Jin et al. Interior tomography with continuous singular value decomposition
Haque et al. Adaptive projection selection for computed tomography
Abascal et al. A novel prior-and motion-based compressed sensing method for small-animal respiratory gated CT
Cierniak An analytical iterative statistical algorithm for image reconstruction from projections
Kim et al. An effective post-filtering framework for 3-D PET image denoising based on noise and sensitivity characteristics
Abascal et al. Investigation of different sparsity transforms for the PICCS algorithm in small-animal respiratory gated CT
Langet et al. Compressed‐sensing‐based content‐driven hierarchical reconstruction: Theory and application to C‐arm cone‐beam tomography
JP2017221339A (ja) Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム
CN111080740A (zh) 一种图像校正方法、装置、设备及介质
Guo et al. Iterative image reconstruction for limited-angle CT using optimized initial image
CN110827370B (zh) 一种非等厚构件的多能ct循环迭代重建方法
JP2023546208A (ja) ノイズ除去のためのニューラルネットワークモデルをトレーニング及びチューニングする方法並びにシステム
Zheng et al. Few-view computed tomography image reconstruction using mean curvature model with curvature smoothing and surface fitting
Sun et al. Quantification of local reconstruction accuracy for helical CT with motion correction
Gopi et al. Iterative computed tomography reconstruction from sparse-view data
Guhathakurta et al. Reducing Computed Tomography Reconstruction and Beam Hardening Artefacts by Data Fusion
CN109903353B (zh) 一种迭代演化模型的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