CN109840927A - 一种基于各向异性全变分的有限角度ct重建算法 - Google Patents
一种基于各向异性全变分的有限角度ct重建算法 Download PDFInfo
- Publication number
- CN109840927A CN109840927A CN201910069082.8A CN201910069082A CN109840927A CN 109840927 A CN109840927 A CN 109840927A CN 201910069082 A CN201910069082 A CN 201910069082A CN 109840927 A CN109840927 A CN 109840927A
- Authority
- CN
- China
- Prior art keywords
- image
- projection
- follows
- data
- formula
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 claims abstract description 49
- 238000002247 constant time method Methods 0.000 claims abstract description 11
- 230000008569 process Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000007781 pre-processing Methods 0.000 claims description 5
- 230000017105 transposition Effects 0.000 claims description 4
- 238000010521 absorption reaction Methods 0.000 claims description 3
- VWDWKYIASSYTQR-UHFFFAOYSA-N sodium nitrate Chemical compound [Na+].[O-][N+]([O-])=O VWDWKYIASSYTQR-UHFFFAOYSA-N 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 2
- 241000208340 Araliaceae Species 0.000 claims 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims 1
- 235000003140 Panax quinquefolius Nutrition 0.000 claims 1
- 235000008434 ginseng Nutrition 0.000 claims 1
- 230000000007 visual effect Effects 0.000 abstract description 11
- 230000000694 effects Effects 0.000 abstract description 3
- 238000013459 approach Methods 0.000 abstract description 2
- 238000002474 experimental method Methods 0.000 abstract description 2
- 238000004800 variational method Methods 0.000 description 6
- 230000012447 hatching Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 102000011990 Sirtuin Human genes 0.000 description 2
- 108050002485 Sirtuin Proteins 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 239000012472 biological sample Substances 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000011478 gradient descent method Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
-
- 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/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2210/00—Indexing scheme for image generation or computer graphics
- G06T2210/41—Medical
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/436—Limited angle
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Radiology & Medical Imaging (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- High Energy & Nuclear Physics (AREA)
- Animal Behavior & Ethology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Surgery (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Algebra (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pulmonology (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种基于各向异性全变分的有限角度CT重建方法，其通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。
Description
技术领域
本发明属于CT成像技术领域，具体涉及一种基于各向异性全变分的有限角度CT重建方法。
背景技术
X光CT在很多领域中都有很重要的应用，如用于诊断和治疗的医学成像领域、安全检查领域以及产品质量检测控制领域等。除了常规的图像重建如低剂量CT重建和稀疏视角CT重建，有限角度问题也受到了广泛的关注。在医学领域，为了病人的健康考虑，希望能够尽量减少病人受到的X光照射，同时由于病人本身往往很难长时间的保持不动配合X光照射，所以希望能够尽可能的缩短扫描时间，这一方面可以通过降低X光剂量来实现，另一方面可以通过减少X光照射时间来实现，包括以更加稀疏的视角或者限制投影视角的大小来获取投影数据进行重建。另外还有一些实际场合也会有扫描角度的限制，比如用X光显微镜对生物样品进行成像时会受到样品固定器的限制而无法进行全角度的成像，或者一些安全检查时由于设备的限制而无法对物品进行全角度的成像等。有限角度CT重建是一个严重的病态问题，因为投影数据的角度范围小于精确重建的理论要求；为了更好的解决这个问题，一些先验信息，如图像的非负性、轮廓或边界的信息、图像的稀疏性等等，常常被用来作为问题的约束条件。
对于大多数的自然图像而言，尤其是医学图像，快速变化的区域只存在于一些结构的边界处，图像的大部分区域都是局部平滑的，所以即使一幅图像其本身并不稀疏，它的梯度图像却很有可能是稀疏的。全变分(Total Variation，TV)极小化就利用了图像的这种稀疏性质，常被应用于CT图像重建领域；一幅图像的TV就是其梯度图像的L1范数，对其进行极小化可以作为由CT投影得到的数据保真项的一个约束条件，这种方法可以得到较好的重建结果，伪影较少，且边界也恢复的较好。
由于我们研究的是有限角度CT重建问题，所以除了图像的稀疏性这一先验信息之外，还有另一个可以加以利用的先验信息，即角度范围信息。1988年Quinto从理论上分析了从有限角度投影数据重建得到的图像的性质特点，给出了一个结论：与投影方向相切的方向的边界和细节信息更容易被恢复，而在另一些特定的角度则会出现一些伪影和模糊。TV极小化方法在进行图像重建时无法探测到某些特定角度的较为模糊的边界，此外由于TV是一幅图像中所有像素点的梯度幅值之和是各向同性的，它同时也会包含那些模糊边界的贡献，导致对不模糊的边界的保边能力也会下降。于是几年前各向异性全变分(AnisotropicTotal Variation，ATV)的概念被提出，除了图像的稀疏性，它将投影数据的角度范围作为另一个先验信息来帮助重建，这个方法的主要思想就是尽可能的减少模糊边界的信息对边界探测的影响，从而可以得到更优的重建图像。
当前已存在很多图像重建算法，通常是基于对一些经典的迭代算法如梯度下降法(Gradient Descent，GD)、代数重建算法(Algebraic Reconstruction Technique，ART)、同步迭代重建技术(Simultaneous Iterative Reconstruction Technique，SIRT)、同步代数重建算法(Simultaneous Algebraic Reconstruction Technique，SART)等进行改进得到的，常常还会结合有序子集法(Ordered-Subsets Technique，OS)来进行加速。这些算法中的大部分主要存在两个问题：一是在从投影数据重建图像时精度不高，二是收敛速度慢导致重建所需要的时间很长。Kudo等人在2016年提出了一种快速迭代重建算法，用来求解低剂量CT和稀疏视角CT重建问题，该方法通过极小化数据保真项来进行迭代重建，并用TV惩罚项来对此过程进行约束，其求解思路如下：首先，通过使用拉格朗日对偶方法，将最初的极小化问题转换成鞍点(原始-对偶)问题，然后使用一阶原始-对偶方法中的交替投影接近法(Alternating Projection Proximal，APP)来对其进行求解；由于这种方法的整体结构和同时迭代重建算法如SIRT和SART算法相同，收敛速度很慢，于是引入了滤波反投影(Filtered Back-Projection，FBP)重建算法中的斜坡滤波器来对迭代方程进行预处理，这种处理方式并不会影响最终的迭代结果，其值与原始问题的解完全相同，最终得到的算法可以理解为通过对斜坡滤波器的使用来达到利用FBP型预处理以加速一阶原始-对偶方法的目的。
发明内容
鉴于上述，本发明提供了一种基于各向异性全变分的有限角度CT重建方法，该方法除了能够解决现有CT重建算法存在的部分边界模糊的问题，提高CT图像重建的质量之外，还能够以很快的收敛速度精确求解目标方程。
一种基于各向异性全变分的有限角度CT重建方法，包括如下步骤：
(1)利用探测器在不同角度方向采集获取CT图像的投影数据，组成投影数据集所述投影数据集由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；
(2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；
低剂量模式下的CT图像重建方程表达式为：
稀疏视角模式下的CT图像重建方程表达式为：
其中：为CT图像数据集且维度为k，k为CT图像的像素点总数，中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| ||ATV表示各向异性全变分；
(3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；
(4)根据实际情况选择对应模式下的目标函数，并对其优化求解以重建得到CT图像。
进一步地，所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：
首先，引入变量将低剂量模式下的CT图像重建方程改写为如下形式：
其中：表示由CT图像数据集计算得到的正向投影数据；
然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：
式(1.3)对应的拉格朗日方程如下：
其中：为拉格朗日乘子向量,T表示转置；
则将式(1.3)转换成鞍点问题，具体表达形式如下：
最后，通过极小化式(1.4)将变量去除，从而得到对应的目标函数如下：
其中：·表示内积运算。
进一步地，所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：
首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：
式(2.2)对应的拉格朗日方程如下：
则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：
其中：为拉格朗日乘子向量,T表示转置，β为给定的权重系数。
对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：
P＝F-1R(ω)F
对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：
P＝F-1R(ω)F
其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
进一步地，所述步骤(4)中采用交替投影接近算法对目标函数进行优化求解以重建得到CT图像，该算法中涉及到的各项参数包括总迭代次数Niter、TV去噪迭代次数NATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ、ATV权重因子γh和γv；其中，总迭代次数Niter的取值范围为1～50000，TV去噪迭代次数NATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围为10-10～1，TV项权重系数β的取值范围为10-6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10-3～103，参数σ的取值范围为10-3～103，ATV权重因子γh和γv的取值范围为10-6～1。
本发明CT图像重建方法通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。
附图说明
图1为本发明CT图像重建算法的流程示意图；
图2为Shepp&Logan phantom模型的真值图像，显示灰度范围为[1.00,1.05]。
图3为Shepp&Loganphantom模型真值图像对应的灰度数值图。
图4为Shepp&Loganphantom模型真值图像对应的正弦图(即不同角度投影数据组成的图，角度范围为180°，角度数为180)。
图5为重建过程中两种方法(本发明方法和各向异性全变分方法)得到的重建图像相比于真值图像的均方根误差(rootmean square error，RMSE)随迭代次数的变化示意图，投影数据的角度范围为120°，角度数为120。
图6(a)为Shepp&Loganphantom模型采用本发明方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。
图6(b)为Shepp&Loganphantom模型采用各向异性全变分方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。
图7为两种方法(本发明方法和各向异性全变分方法)得到的重建图像的水平中线的剖线与真值剖线的对比示意图，迭代次数为3000次。
具体实施方式
为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。
如图1所示，本发明基于各向异性全变分的有限角度CT快速迭代重建算法，包括如下步骤：
(1)采集探测器在不同方向测得的CT图像的投影数据，组成投影数据集投影数据集由对应各个角度方向下所采集得到的投影数据向量组成，投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数。
(2)分别建立低剂量和稀疏视角CT的图像重建方程式：
低剂量CT：
稀疏视角CT：其中：为CT图像数据集且维度为k，k为CT图像的像素点总数，中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A是系统矩阵，它适用于平行光束投影、扇形光束投影及锥形光束投影等各种形式；β为权重系数；为各向异性全变分项。
在低剂量CT情况下，我们假定nm＞k并且含有噪声；在稀疏视角CT情况下，我们假定nm＜k并且不含噪声，m为角度方向数量。
(3)对公式(1)和(2)分别进行预处理：
对于低剂量CT，引入变量则式(1)可改写成如下形式：
其中：可以看作是从图像计算得到的正向投影。引入非奇异矩阵P1/2，则上式可进一步改写为：
上式对应的拉格朗日方程可写成如下形式：
其中：是拉格朗日乘子向量，则式(3)式可转换成鞍点问题：
因为可以通过极小化上式去除，于是我们可以得到：
P＝F-1R(ω)F
其中：·表示内积运算，F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量。
对于稀疏视角CT，同样引入非奇异矩阵P1/2，则式(2)可改写成如下形式：
同样的可写出上式对应的拉格朗日方程：
则式(5)可转换成如下形式：
P＝F-1R(ω)F
(4)使用交替投影接近算法APP来对式(4)和(6)进行求解：
4-1初始化并设置各项参数的值：总迭代次数Niter、TV去噪迭代次数NATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ，ATV权重因子γh和γv；其中，总迭代次数Niter的取值范围为1～50000，TV去噪迭代次数NATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围均为10-10～1，TV项权重系数β取值范围为10-6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10-3～103，参数σ的取值范围为10-3～103，ATV权重因子γh和γv取值范围均为10-6～1。
4-2初始令迭代计数k＝1；
4-3对进行更新，当k为1时：
当k不为1时：
低剂量CT：
稀疏视角CT：
4-4对进行更新：
可用Chambolle的TV去噪算法来对上式进行求解，由于这一项在这一步中是常量，为方便表述引入则式(7)的求解方程为：
其中：为梯度算子，假定重建图像的大小为N×N像素，在本实施方式中，其定义如下：
则各项异性全变分项的计算方式如下：
对每一个div算子的计算方式如下：
计算若∈c＞∈且重复次数未达到NATV，重复计算式(8)～式(9)。
4-5再次对进行更新：
低剂量CT：
稀疏视角CT：
4-6判断重建图像与重建图像的差值图像的所有像素值的平方和是否小于迭代终止阈值δ或k是否大于Niter：若是，执行步骤4-7；若否，令k＝k+1，执行步骤4-3～步骤4-5。
4-7终止迭代，得到最终重建图像
以下我们通过对Shepp&Loganphantom模型的正弦图进行重建来验证实施方法的实用性和可靠性。Shepp&Loganphantom模型的真值图像如图2所示(由于其大部分结构的灰度值非常接近，为了清楚的显示其结构，这里选择灰度范围[1.00,1.05]进行显示)，其对应的灰度数值图如图3所示，正弦图(角度范围为180°，角度数为180)如图4所示。采用本发明方法(proposed algorithm)和各向异性全变分方法(ART+ATV)进行重建，图5给出了重建过程中两种方法得到的重建图像相比于真值图的均方根误差随迭代次数的变化，可以看到两种方法都随着迭代次数的增加而逐渐收敛，本发明方法的重建结果均方根误差更小，且其收敛速度非常的快。我们选取最大迭代次数(即3000次)时的重建图像加以显示，图6(a)为采用各向异性全变分方法重建后的CT图像，图6(b)为采用本发明方法重建后的CT图像，可看到本发明方法重建的图像伪影要小得多，能更好的恢复图像的细节信息，重建效果显著优于各向异性全变分方法。为了更直观的证明本发明方法的优越性，图7给出了图6(a)和图6(b)所示的重建图像的水平中线的剖线与真值剖线的对比，可以看到，本发明方法的剖线几乎和真值完全重合，即其能够得到更加接近于真值的重建结果。
上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。
Claims (6)
1.一种基于各向异性全变分的有限角度CT重建方法，包括如下步骤：
(1)利用探测器在不同角度方向采集获取CT图像的投影数据，组成投影数据集所述投影数据集由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；
(2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；
低剂量模式下的CT图像重建方程表达式为：
稀疏视角模式下的CT图像重建方程表达式为：
其中：为CT图像数据集且维度为k，k为CT图像的像素点总数，中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| ||ATV表示各向异性全变分；
(3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；
(4)根据实际情况选择对应模式下的目标函数，并对其优化求解以重建得到CT图像。
2.根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：
首先，引入变量将低剂量模式下的CT图像重建方程改写为如下形式：
其中：表示由CT图像数据集计算得到的正向投影数据；
然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：
式(1.3)对应的拉格朗日方程如下：
其中：为拉格朗日乘子向量,T表示转置；
则将式(1.3)转换成鞍点问题，具体表达形式如下：
最后，通过极小化式(1.4)将变量去除，从而得到对应的目标函数如下：
其中：·表示内积运算。
3.根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：
首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：
式(2.2)对应的拉格朗日方程如下：
则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：
其中：为拉格朗日乘子向量,T表示转置，β为给定的权重系数。
4.根据权利要求2所述的有限角度CT重建方法，其特征在于：对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：
P＝F-1R(ω)F
其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
5.根据权利要求3所述的有限角度CT重建方法，其特征在于：对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：
P＝F-1R(ω)F
其中：F和F-1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
6.根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(4)中采用交替投影接近算法对目标函数进行优化求解以重建得到CT图像，该算法中涉及到的各项参数包括总迭代次数Niter、TV去噪迭代次数NATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ、ATV权重因子γh和γv；其中，总迭代次数Niter的取值范围为1～50000，TV去噪迭代次数NATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围为10-10～1，TV项权重系数β的取值范围为10-6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10-3～103，参数σ的取值范围为10-3～103，ATV权重因子γh和γv的取值范围为10-6～1。
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910069082.8A CN109840927B (zh) | 2019-01-24 | 2019-01-24 | 一种基于各向异性全变分的有限角度ct重建算法 |
| US16/979,168 US11727609B2 (en) | 2019-01-24 | 2019-11-20 | Limited-angle CT reconstruction method based on anisotropic total variation |
| PCT/CN2019/126882 WO2020151424A1 (zh) | 2019-01-24 | 2019-12-20 | 一种基于各向异性全变分的有限角度ct重建算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910069082.8A CN109840927B (zh) | 2019-01-24 | 2019-01-24 | 一种基于各向异性全变分的有限角度ct重建算法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109840927A true CN109840927A (zh) | 2019-06-04 |
| CN109840927B CN109840927B (zh) | 2020-11-10 |
Family
ID=66884144
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910069082.8A Active CN109840927B (zh) | 2019-01-24 | 2019-01-24 | 一种基于各向异性全变分的有限角度ct重建算法 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US11727609B2 (zh) |
| CN (1) | CN109840927B (zh) |
| WO (1) | WO2020151424A1 (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110717956A (zh) * | 2019-09-30 | 2020-01-21 | 重庆大学 | 一种有限角投影超像素引导的l0范数最优化重建方法 |
| WO2020151424A1 (zh) * | 2019-01-24 | 2020-07-30 | 浙江大学 | 一种基于各向异性全变分的有限角度ct重建算法 |
| CN111696166A (zh) * | 2020-06-10 | 2020-09-22 | 浙江大学 | 基于fdk型预处理矩阵的圆周锥束ct快速迭代重建方法 |
| CN112070856A (zh) * | 2020-09-16 | 2020-12-11 | 重庆师范大学 | 基于非下采样轮廓波变换的有限角c型臂ct图像重建方法 |
| CN116509453A (zh) * | 2023-06-29 | 2023-08-01 | 南京科进实业有限公司 | 一种基于超声波的胫骨骨密度检测系统及方法 |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113554729B (zh) * | 2021-07-28 | 2022-06-07 | 中北大学 | 一种ct图像重建方法及系统 |
| CN113706653B (zh) * | 2021-09-23 | 2023-12-08 | 明峰医疗系统股份有限公司 | 基于AwTV的CT迭代重建替代函数优化方法 |
| CN114549684B (zh) * | 2022-03-05 | 2023-03-24 | 昆明理工大学 | 基于离散余弦变换和代数重建算法改进矿井无线电波透视成像重建方法 |
| CN115249028B (zh) * | 2022-05-13 | 2023-06-23 | 哈尔滨工业大学 | 一种基于稀疏正则化约束的盲解卷积无线通信信号重构方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102376084A (zh) * | 2010-08-12 | 2012-03-14 | 西门子公司 | 使用各向异性噪声模型对ct图像的迭代图像滤波 |
| CN103679706A (zh) * | 2013-11-26 | 2014-03-26 | 中国人民解放军第四军医大学 | 一种基于图像各向异性边缘检测的ct稀疏角度重建方法 |
| CN104657950A (zh) * | 2015-02-16 | 2015-05-27 | 浙江大学 | 一种基于Poisson TV的动态PET图像重建方法 |
| CN108550172A (zh) * | 2018-03-07 | 2018-09-18 | 浙江大学 | 一种基于非局部特性和全变分联合约束的pet图像重建方法 |
| US10140733B1 (en) * | 2017-09-13 | 2018-11-27 | Siemens Healthcare Gmbh | 3-D vessel tree surface reconstruction |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9626587B2 (en) * | 2014-08-28 | 2017-04-18 | Toshiba Medical Systems Corporation | Iterative reconstruction scheme for phase contrast tomography |
| CN106846427B (zh) | 2017-01-25 | 2019-06-11 | 浙江大学 | 一种基于重加权各向异性全变分的有限角度ct重建方法 |
| CN108280859B (zh) | 2017-12-25 | 2021-03-30 | 华南理工大学 | 一种采样角度受限下的ct稀疏投影图像重建方法及装置 |
| CN109840927B (zh) * | 2019-01-24 | 2020-11-10 | 浙江大学 | 一种基于各向异性全变分的有限角度ct重建算法 |
-
2019
- 2019-01-24 CN CN201910069082.8A patent/CN109840927B/zh active Active
- 2019-11-20 US US16/979,168 patent/US11727609B2/en active Active
- 2019-12-20 WO PCT/CN2019/126882 patent/WO2020151424A1/zh active Application Filing
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102376084A (zh) * | 2010-08-12 | 2012-03-14 | 西门子公司 | 使用各向异性噪声模型对ct图像的迭代图像滤波 |
| CN103679706A (zh) * | 2013-11-26 | 2014-03-26 | 中国人民解放军第四军医大学 | 一种基于图像各向异性边缘检测的ct稀疏角度重建方法 |
| CN104657950A (zh) * | 2015-02-16 | 2015-05-27 | 浙江大学 | 一种基于Poisson TV的动态PET图像重建方法 |
| US10140733B1 (en) * | 2017-09-13 | 2018-11-27 | Siemens Healthcare Gmbh | 3-D vessel tree surface reconstruction |
| CN108550172A (zh) * | 2018-03-07 | 2018-09-18 | 浙江大学 | 一种基于非局部特性和全变分联合约束的pet图像重建方法 |
Non-Patent Citations (1)
| Title |
|---|
| HONGLISHI 等: "A Novel Iterative CT Reconstruction", 《JOURNAL.PONE》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020151424A1 (zh) * | 2019-01-24 | 2020-07-30 | 浙江大学 | 一种基于各向异性全变分的有限角度ct重建算法 |
| CN110717956A (zh) * | 2019-09-30 | 2020-01-21 | 重庆大学 | 一种有限角投影超像素引导的l0范数最优化重建方法 |
| CN110717956B (zh) * | 2019-09-30 | 2023-06-20 | 重庆大学 | 一种有限角投影超像素引导的l0范数最优化重建方法 |
| CN111696166A (zh) * | 2020-06-10 | 2020-09-22 | 浙江大学 | 基于fdk型预处理矩阵的圆周锥束ct快速迭代重建方法 |
| CN112070856A (zh) * | 2020-09-16 | 2020-12-11 | 重庆师范大学 | 基于非下采样轮廓波变换的有限角c型臂ct图像重建方法 |
| CN112070856B (zh) * | 2020-09-16 | 2022-08-26 | 重庆师范大学 | 基于非下采样轮廓波变换的有限角c型臂ct图像重建方法 |
| CN116509453A (zh) * | 2023-06-29 | 2023-08-01 | 南京科进实业有限公司 | 一种基于超声波的胫骨骨密度检测系统及方法 |
| CN116509453B (zh) * | 2023-06-29 | 2023-09-01 | 南京科进实业有限公司 | 一种基于超声波的胫骨骨密度检测系统及方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| US20200402274A1 (en) | 2020-12-24 |
| CN109840927B (zh) | 2020-11-10 |
| WO2020151424A1 (zh) | 2020-07-30 |
| US11727609B2 (en) | 2023-08-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109840927A (zh) | 一种基于各向异性全变分的有限角度ct重建算法 | |
| Tang et al. | Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms | |
| Niu et al. | Sparse-view x-ray CT reconstruction via total generalized variation regularization | |
| JP6855223B2 (ja) | 医用画像処理装置、ｘ線コンピュータ断層撮像装置及び医用画像処理方法 | |
| Han et al. | Algorithm-enabled low-dose micro-CT imaging | |
| Zhang et al. | Regularization strategies in statistical image reconstruction of low‐dose x‐ray CT: A review | |
| EP2232445B1 (en) | Method for image reconstruction using sparsity-constrained correction | |
| Saba et al. | Maximizing quantitative accuracy of lung airway lumen and wall measures obtained from X-ray CT imaging | |
| CN106846427B (zh) | 一种基于重加权各向异性全变分的有限角度ct重建方法 | |
| Wang et al. | ADMM-based deep reconstruction for limited-angle CT | |
| Li et al. | Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV) | |
| Niu et al. | Iterative reconstruction for sparse-view x-ray CT using alpha-divergence constrained total generalized variation minimization | |
| Boudjelal et al. | Improved simultaneous algebraic reconstruction technique algorithm for positron-emission tomography image reconstruction via minimizing the fast total variation | |
| Mahmoudi et al. | Sparse-view statistical image reconstruction with improved total variation regularization for X-ray micro-CT imaging | |
| Xi et al. | Study of CT image reconstruction algorithm based on high order total variation | |
| Cierniak | An analytical iterative statistical algorithm for image reconstruction from projections | |
| Zhang et al. | Projection domain denoising method based on dictionary learning for low-dose CT image reconstruction | |
| Friot et al. | Iterative tomographic reconstruction with TV prior for low-dose CBCT dental imaging | |
| Bubba et al. | Shearlet-based regularized reconstruction in region-of-interest computed tomography | |
| Qi et al. | Sparse-view computed tomography image reconstruction via a combination of L 1 and SL 0 regularization | |
| Wang et al. | Guided image filtering based limited-angle CT reconstruction algorithm using wavelet frame | |
| Choi et al. | Integration of 2D iteration and a 3D CNN-based model for multi-type artifact suppression in C-arm cone-beam CT | |
| Gopi et al. | Iterative computed tomography reconstruction from sparse-view data | |
| Li et al. | Robust frame based X-ray CT reconstruction | |
| Sakhaee et al. | Gradient-based sparse approximation for computed tomography |
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 |