WO2020151424A1 - 一种基于各向异性全变分的有限角度ct重建算法 - Google Patents

一种基于各向异性全变分的有限角度ct重建算法 Download PDF

Info

Publication number
WO2020151424A1
WO2020151424A1 PCT/CN2019/126882 CN2019126882W WO2020151424A1 WO 2020151424 A1 WO2020151424 A1 WO 2020151424A1 CN 2019126882 W CN2019126882 W CN 2019126882W WO 2020151424 A1 WO2020151424 A1 WO 2020151424A1
Authority
WO
WIPO (PCT)
Prior art keywords
reconstruction
image
angle
equation
projection data
Prior art date
Application number
PCT/CN2019/126882
Other languages
English (en)
French (fr)
Inventor
刘华锋
王婷
Original Assignee
浙江大学
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
Priority to US16/979,168 priority Critical patent/US11727609B2/en
Application filed by 浙江大学 filed Critical 浙江大学
Publication of WO2020151424A1 publication Critical patent/WO2020151424A1/zh

Links

Images

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

Definitions

  • the invention belongs to the field of CT imaging technology, and specifically relates to a limited-angle CT reconstruction method based on anisotropic total variation.
  • X-ray CT has very important applications in many fields, such as medical imaging for diagnosis and treatment, safety inspection and product quality inspection and control.
  • medical imaging for diagnosis and treatment, safety inspection and product quality inspection and control.
  • image reconstruction such as low-dose CT reconstruction and sparse view CT reconstruction
  • the problem of limited angles has also received extensive attention.
  • the X-ray exposure to the patient can be reduced as much as possible.
  • the scanning angle is limited. For example, when using an X-ray microscope to image biological samples, it is limited by the sample holder and cannot perform full-angle imaging, or some safety inspections are due to equipment limitations. It is not possible to perform full-angle imaging of objects.
  • Limited-angle CT reconstruction is a serious ill-conditioned problem, because the angle range of the projection data is less than the theoretical requirement of accurate reconstruction; in order to better solve this problem, some prior information, such as image non-negativity, contour or boundary information, Image sparsity, etc., are often used as constraints for the problem.
  • TV Total Variation minimization takes advantage of the sparse nature of the image, and is often used in the field of CT image reconstruction; the TV of an image is the L1 norm of its gradient image, which is minimized Transformation can be used as a constraint condition for the fidelity of the data obtained by CT projection. This method can obtain better reconstruction results, with fewer artifacts, and better boundary restoration.
  • ATV Anisotropic Total Variation
  • This method minimizes the data fidelity term to perform iterative reconstruction and uses the TV penalty term to solve this problem.
  • the process is constrained, and the solution idea is as follows: First, by using the Lagrangian dual method, the initial minimization problem is converted into a saddle point (primitive-dual) problem, and then the alternate projection approach in the first-order primordial-dual method is used to approach Method (Alternating Projection Proximal, APP) to solve it; because the overall structure of this method is the same as the simultaneous iterative reconstruction algorithm such as SIRT and SART algorithm, the convergence speed is very slow, so filtered back-projection (Filtered Back-Projection, FBP) The slope filter in the reconstruction algorithm is used to preprocess the iterative equation.
  • the filter is used to achieve the purpose of using FBP-type preprocessing to accelerate the first-order primal-dual method.
  • the present invention provides a limited-angle CT reconstruction method based on anisotropic total variation.
  • this method improves the quality of CT image reconstruction. It can also accurately solve the target equation with a fast convergence rate.
  • a limited-angle CT reconstruction method based on anisotropic total variation including the following steps:
  • the projection data set It is composed of projection data vectors collected in various angle directions, the dimension of the projection data vector is n and each element value in the vector is the projection data measured by the corresponding detector, and n is the number of detectors;
  • Is the CT image data set and the dimension is k
  • k is the total number of pixels of the CT image
  • Each element in the value is the X-ray absorption coefficient of the corresponding pixel in the CT image to be reconstructed
  • A is the system matrix
  • ⁇ is the given weight coefficient
  • ATV represents the anisotropic total variation
  • the preprocessing process includes introducing variables And non-singular matrix P, using Lagrangian duality to convert the reconstruction equation into a saddle point problem;
  • the preprocessing process of the CT image reconstruction equation in the low-dose mode in the step (3) is as follows:
  • ⁇ represents the inner product operation.
  • the preprocessing process of the CT image reconstruction equation in the sparse view mode in the step (3) is as follows:
  • T represents transpose
  • ⁇ is the given weight coefficient
  • F and F -1 are the one-dimensional Fourier transform operator and the inverse Fourier transform operator, respectively, ⁇ is the frequency domain variable of the Fourier transformed projection data vector collected in the corresponding angle direction, m is the number of angular directions, and ⁇ is a given parameter.
  • an alternate projection approach algorithm is used to optimize the objective function to reconstruct the CT image.
  • the parameters involved in the algorithm include the total number of iterations N iter and the number of TV denoising iterations N ATV , Iteration termination threshold ⁇ and ⁇ , TV term weight coefficient ⁇ , parameter ⁇ , parameter ⁇ , parameter ⁇ , ATV weight factor ⁇ h and ⁇ v ; among them, the value range of the total number of iterations N iter is 1 ⁇ 50,000, TV goes The number of noise iterations N ATV ranges from 1 to 10000, the iteration termination thresholds ⁇ and ⁇ range from 10 -10 to 1, the TV term weight coefficient ⁇ ranges from 10 -6 to 1, and the parameter ⁇ is in the range of 0.01 to 0.25, the range parameter ⁇ is 10 -3 to 10 3, the parameter ⁇ is in the range of 10 -3 to 10 3, the ATV weighting factors ⁇ v and ⁇ h is in the range of 10 -6 ⁇ 1.
  • the CT image reconstruction method of the present invention combines the fast iterative reconstruction algorithm with the anisotropic total variation method by using the low-dose and sparse view CT image reconstruction model, effectively solving the ubiquitous existing limited-angle CT reconstruction methods Part of the area is blurred, the convergence speed is slow, and the problem cannot be solved accurately.
  • the slope filter in the filtered back projection reconstruction algorithm is introduced to preprocess the iterative equation, and the alternate projection approach method is used to solve the problem. It is solved and iterated repeatedly until the termination condition is satisfied; the comparison with the experiment of the existing reconstruction method shows that the present invention can achieve better reconstruction effect.
  • Figure 1 is a schematic flow chart of the CT image reconstruction algorithm of the present invention
  • Figure 2 is the true value image of the Shepp&Logan phantom model, and the display gray scale range is [1.00,1.05].
  • Figure 3 is the gray scale value map corresponding to the truth image of the Shepp & Logan phantom model.
  • Figure 4 is the sine diagram corresponding to the truth image of the Shepp & Logan phantom model (that is, a diagram composed of projection data from different angles, the angle range is 180°, and the angle number is 180).
  • Figure 5 shows the variation of the root mean square error (RMSE) of the reconstructed image obtained by the two methods (the method of the present invention and the anisotropic total variation method) compared with the true value image during the reconstruction process with the number of iterations
  • RMSE root mean square error
  • Fig. 6(a) is a CT image reconstructed by the Shepp & Logan phantom model using the method of the present invention, the projection data angle range is 120°, and the number of iterations is 3000.
  • Figure 6(b) shows the CT image reconstructed by the Shepp & Logan phantom model using the anisotropic total variation method.
  • the projection data angle range is 120°, and the number of iterations is 3000.
  • FIG. 7 is a schematic diagram of the comparison between the horizontal midline profile and the true value profile of the reconstructed image obtained by two methods (the method of the present invention and the anisotropic total variation method), and the number of iterations is 3000.
  • the fast iterative reconstruction algorithm of limited angle CT based on anisotropic total variation of the present invention includes the following steps:
  • Projection data set It is composed of projection data vectors collected under corresponding angle directions.
  • the dimension of the projection data vector is n and each element value in the vector is the projection data measured by the corresponding detector, and n is the number of detectors.
  • Sparse view CT among them: Is the CT image data set and the dimension is k, k is the total number of pixels of the CT image, Each element value is the X-ray absorption coefficient of the corresponding pixel in the CT image to be reconstructed.
  • A is the system matrix, which is suitable for parallel beam projection, fan beam projection and cone beam projection; ⁇ is the weighting coefficient ; It is an anisotropic total variation item.
  • ⁇ represents the inner product operation
  • F and F -1 are the one-dimensional Fourier transform operator and the inverse Fourier transform operator, respectively
  • ⁇ is the projection data vector collected in the corresponding angle direction after Fourier transform After the frequency domain variable.
  • Equation (5) can be transformed into the following form:
  • each parameter total number of iterations N iter , TV denoising iteration number N ATV , iteration termination threshold ⁇ and ⁇ , TV term weight coefficient ⁇ , parameter ⁇ , parameter ⁇ , parameter ⁇ , ATV weight factor ⁇ h and ⁇ v ; where the total number of iterations N iter ranges from 1 to 50000, the number of TV denoising iterations N ATV ranges from 1 to 10000, and the iteration termination thresholds ⁇ and ⁇ both range from 10 -10 ⁇ 1, the value range of the TV item weight coefficient ⁇ is 10 -6 ⁇ 1, the value range of parameter ⁇ is 0.01 ⁇ 0.25, the value range of parameter ⁇ is 10 -3 ⁇ 10 3 , the value range of parameter ⁇ is 10 -3 ⁇ 10 3 , the value range of ATV weighting factors ⁇ h and ⁇ v are both 10 -6 ⁇ 1.
  • Figure 5 shows the root mean square error of the reconstructed image obtained by the two methods in the reconstruction process compared to the true value map With the change of the number of iterations, it can be seen that both methods gradually converge as the number of iterations increases.
  • the root mean square error of the reconstruction result of the method of the present invention is smaller, and the convergence speed is very fast.
  • Figure 6(a) is the CT image reconstructed using the anisotropic total variation method
  • Figure 6(b) is the reconstructed image using the method of the present invention.
  • Fig. 7 shows the comparison of the horizontal midline and the true value of the reconstructed image shown in Fig. 6(a) and Fig. 6(b).
  • the section line of the method of the present invention almost completely coincides with the true value, that is, it can obtain a reconstruction result closer to the true value.

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)
  • Pulmonology (AREA)
  • Mathematical Physics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种基于各向异性全变分的有限角度CT重建方法，其通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。

Description

一种基于各向异性全变分的有限角度CT重建算法 技术领域
本发明属于CT成像技术领域，具体涉及一种基于各向异性全变分的有限角度CT重建方法。
背景技术
X光CT在很多领域中都有很重要的应用，如用于诊断和治疗的医学成像领域、安全检查领域以及产品质量检测控制领域等。除了常规的图像重建如低剂量CT重建和稀疏视角CT重建，有限角度问题也受到了广泛的关注。在医学领域，为了病人的健康考虑，希望能够尽量减少病人受到的X光照射，同时由于病人本身往往很难长时间的保持不动配合X光照射，所以希望能够尽可能的缩短扫描时间，这一方面可以通过降低X光剂量来实现，另一方面可以通过减少X光照射时间来实现，包括以更加稀疏的视角或者限制投影视角的大小来获取投影数据进行重建。另外还有一些实际场合也会有扫描角度的限制，比如用X光显微镜对生物样品进行成像时会受到样品固定器的限制而无法进行全角度的成像，或者一些安全检查时由于设备的限制而无法对物品进行全角度的成像等。有限角度CT重建是一个严重的病态问题，因为投影数据的角度范围小于精确重建的理论要求；为了更好的解决这个问题，一些先验信息，如图像的非负性、轮廓或边界的信息、图像的稀疏性等等，常常被用来作为问题的约束条件。
对于大多数的自然图像而言，尤其是医学图像，快速变化的区域只存在于一些结构的边界处，图像的大部分区域都是局部平滑的，所以即使一幅图像其本身并不稀疏，它的梯度图像却很有可能是稀疏的。全变分(Total Variation，TV)极小化就利用了图像的这种稀疏性质，常被应用于CT图像重建领域；一幅图像的TV就是其梯度图像的L1范数，对其进行极小化可以作为由CT投影得到的数据保真项的一个约束条件，这种方法可以得到较好的重建结果，伪影较少，且边界也恢复的较好。
由于我们研究的是有限角度CT重建问题，所以除了图像的稀疏性这一先验信息之外，还有另一个可以加以利用的先验信息，即角度范围信息。1988年Quinto从理论上分析了从有限角度投影数据重建得到的图像的性质特点，给出了一个结论：与投影方向相切的方向的边界和细节信息更容易被恢复，而在另一些特定的角度则会出现一些伪影和模糊。TV极小化方法在进行图像重建时无法探测到某些特定角度的较为模糊的边界，此外由于TV是一幅图像中所有像素点的梯度幅值之和是各向同性的，它同时也会包含那些模糊边界的贡献，导致对不模糊的边界的保边能力也会下降。于是几年前各向异性全变分(Anisotropic Total 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图像的投影数据，组成投影数据集
Figure PCTCN2019126882-appb-000001
所述投影数据集
Figure PCTCN2019126882-appb-000002
由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；
(2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；
低剂量模式下的CT图像重建方程表达式为：
Figure PCTCN2019126882-appb-000003
稀疏视角模式下的CT图像重建方程表达式为：
Figure PCTCN2019126882-appb-000004
其中：
Figure PCTCN2019126882-appb-000005
为CT图像数据集且维度为k，k为CT图像的像素点总数，
Figure PCTCN2019126882-appb-000006
中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| || ATV表示各向异性全变分；
(3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量
Figure PCTCN2019126882-appb-000007
和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；
(4)根据实际情况选择对应模式下的目标函数，并对其优化求解以重建得到CT图像。
进一步地，所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：
首先，引入变量
Figure PCTCN2019126882-appb-000008
将低剂量模式下的CT图像重建方程改写为如下形式：
Figure PCTCN2019126882-appb-000009
其中：
Figure PCTCN2019126882-appb-000010
表示由CT图像数据集
Figure PCTCN2019126882-appb-000011
计算得到的正向投影数据；
然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：
Figure PCTCN2019126882-appb-000012
式(1.3)对应的拉格朗日方程
Figure PCTCN2019126882-appb-000013
如下：
Figure PCTCN2019126882-appb-000014
其中：
Figure PCTCN2019126882-appb-000015
为拉格朗日乘子向量, T表示转置；
则将式(1.3)转换成鞍点问题，具体表达形式如下：
Figure PCTCN2019126882-appb-000016
最后，通过极小化式(1.4)将变量
Figure PCTCN2019126882-appb-000017
去除，从而得到对应的目标函数如下：
Figure PCTCN2019126882-appb-000018
其中：·表示内积运算。
进一步地，所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：
首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：
Figure PCTCN2019126882-appb-000019
式(2.2)对应的拉格朗日方程
Figure PCTCN2019126882-appb-000020
如下：
Figure PCTCN2019126882-appb-000021
则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：
Figure PCTCN2019126882-appb-000022
其中：
Figure PCTCN2019126882-appb-000023
为拉格朗日乘子向量, T表示转置，β为给定的权重系数。
对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：
P＝F -1R(ω)F
Figure PCTCN2019126882-appb-000024
对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：
P＝F -1R(ω)F
Figure PCTCN2019126882-appb-000025
其中：F和F -1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
进一步地，所述步骤(4)中采用交替投影接近算法对目标函数进行优化求解以重建得到CT图像，该算法中涉及到的各项参数包括总迭代次数N iter、TV去噪迭代次数N ATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ、ATV权重因子γ h和γ v；其中，总迭代次数N iter的取值范围为1～50000，TV去噪迭代次数N ATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围为10 -10～1，TV项权重系数β的取值范围为10 -6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10 -3～10 3，参数σ的取值范围为10 -3～10 3，ATV权重因子γ h和γ v的取值范围为10 -6～1。
本发明CT图像重建方法通过利用低剂量和稀疏视角CT的图像重建模型，将快速迭代重建算法与各向异性全变分的方法相结合，有效的解决了现有有限角度CT重建方法中普遍存在的部分区域模糊、收敛速度慢及无法精确求解的问题，在对该模型的求解过程中引入滤波反投影重建算法中的斜坡滤波器来对迭代方程进行预处理，并使用交替投影接近法来对其进行求解，反复迭代直至终止条件被满足；与现有重建方法的实验比较表明，本发明能取得较好的重建效果。
附图说明
图1为本发明CT图像重建算法的流程示意图；
图2为Shepp & Logan phantom模型的真值图像，显示灰度范围为 [1.00,1.05]。
图3为Shepp & Logan phantom模型真值图像对应的灰度数值图。
图4为Shepp & Logan phantom模型真值图像对应的正弦图(即不同角度投影数据组成的图，角度范围为180°，角度数为180)。
图5为重建过程中两种方法(本发明方法和各向异性全变分方法)得到的重建图像相比于真值图像的均方根误差(root mean square error，RMSE)随迭代次数的变化示意图，投影数据的角度范围为120°，角度数为120。
图6(a)为Shepp & Logan phantom模型采用本发明方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。
图6(b)为Shepp & Logan phantom模型采用各向异性全变分方法重建后的CT图像，投影数据角度范围为120°，迭代次数为3000次。
图7为两种方法(本发明方法和各向异性全变分方法)得到的重建图像的水平中线的剖线与真值剖线的对比示意图，迭代次数为3000次。
具体实施方式
为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。
如图1所示，本发明基于各向异性全变分的有限角度CT快速迭代重建算法，包括如下步骤：
(1)采集探测器在不同方向测得的CT图像的投影数据，组成投影数据集
Figure PCTCN2019126882-appb-000026
投影数据集
Figure PCTCN2019126882-appb-000027
由对应各个角度方向下所采集得到的投影数据向量组成，投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数。
(2)分别建立低剂量和稀疏视角CT的图像重建方程式：
低剂量CT：
Figure PCTCN2019126882-appb-000028
稀疏视角CT：
Figure PCTCN2019126882-appb-000029
其中：
Figure PCTCN2019126882-appb-000030
为CT图像数据集且维度为k，k为CT图像的像素点总数，
Figure PCTCN2019126882-appb-000031
中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A是系统矩阵，它适 用于平行光束投影、扇形光束投影及锥形光束投影等各种形式；β为权重系数；
Figure PCTCN2019126882-appb-000032
为各向异性全变分项。
在低剂量CT情况下，我们假定nm＞k并且
Figure PCTCN2019126882-appb-000033
含有噪声；在稀疏视角CT情况下，我们假定nm＜k并且
Figure PCTCN2019126882-appb-000034
不含噪声，m为角度方向数量。
(3)对公式(1)和(2)分别进行预处理：
对于低剂量CT，引入变量
Figure PCTCN2019126882-appb-000035
则式(1)可改写成如下形式：
Figure PCTCN2019126882-appb-000036
其中：
Figure PCTCN2019126882-appb-000037
可以看作是从图像
Figure PCTCN2019126882-appb-000038
计算得到的正向投影。引入非奇异矩阵P 1/2，则上式可进一步改写为：
Figure PCTCN2019126882-appb-000039
上式对应的拉格朗日方程可写成如下形式：
Figure PCTCN2019126882-appb-000040
其中：
Figure PCTCN2019126882-appb-000041
是拉格朗日乘子向量，则式(3)式可转换成鞍点问题：
Figure PCTCN2019126882-appb-000042
因为
Figure PCTCN2019126882-appb-000043
可以通过极小化上式去除，于是我们可以得到：
Figure PCTCN2019126882-appb-000044
P＝F -1R(ω)F
Figure PCTCN2019126882-appb-000045
其中：·表示内积运算，F和F -1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量。
对于稀疏视角CT，同样引入非奇异矩阵P 1/2，则式(2)可改写成如下形式：
Figure PCTCN2019126882-appb-000046
同样的可写出上式对应的拉格朗日方程：
Figure PCTCN2019126882-appb-000047
则式(5)可转换成如下形式：
Figure PCTCN2019126882-appb-000048
P＝F -1R(ω)F
Figure PCTCN2019126882-appb-000049
(4)使用交替投影接近算法APP来对式(4)和(6)进行求解：
4-1初始化
Figure PCTCN2019126882-appb-000050
并设置各项参数的值：总迭代次数N iter、TV去噪迭代次数N ATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ，ATV权重因子γ h和γ v；其中，总迭代次数N iter的取值范围为1～50000，TV去噪迭代次数N ATV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围均为10 -10～1，TV项权重系数β取值范围为10 -6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10 -3～10 3，参数σ的取值范围为10 -3～10 3，ATV权重因子γ h和γ v取值范围均为10 -6～1。
4-2初始令迭代计数k＝1；
4-3对
Figure PCTCN2019126882-appb-000051
进行更新，当k为1时：
Figure PCTCN2019126882-appb-000052
当k不为1时：
低剂量CT：
Figure PCTCN2019126882-appb-000053
稀疏视角CT：
Figure PCTCN2019126882-appb-000054
4-4对
Figure PCTCN2019126882-appb-000055
进行更新：
Figure PCTCN2019126882-appb-000056
可用Chambolle的TV去噪算法来对上式进行求解，由于
Figure PCTCN2019126882-appb-000057
这一项在这一步中是常量，为方便表述引入
Figure PCTCN2019126882-appb-000058
则式(7)的求解方程为：
Figure PCTCN2019126882-appb-000059
Figure PCTCN2019126882-appb-000060
其中：
Figure PCTCN2019126882-appb-000061
为梯度算子，假定重建图像的大小为N×N像素，在本实施方式中，其定义如下：
Figure PCTCN2019126882-appb-000062
Figure PCTCN2019126882-appb-000063
Figure PCTCN2019126882-appb-000064
则各项异性全变分项
Figure PCTCN2019126882-appb-000065
的计算方式如下：
Figure PCTCN2019126882-appb-000066
对每一个
Figure PCTCN2019126882-appb-000067
div算子的计算方式如下：
Figure PCTCN2019126882-appb-000068
计算
Figure PCTCN2019126882-appb-000069
若∈ c＞∈且重复次数未达到N ATV，重复计算式(8)～式(9)。
4-5再次对
Figure PCTCN2019126882-appb-000070
进行更新：
低剂量CT：
Figure PCTCN2019126882-appb-000071
稀疏视角CT：
Figure PCTCN2019126882-appb-000072
4-6判断重建图像
Figure PCTCN2019126882-appb-000073
与重建图像
Figure PCTCN2019126882-appb-000074
的差值图像的所有像素值的平方和是否小于迭代终止阈值δ或k是否大于N iter：若是，执行步骤4-7；若否，令k＝k+1，执行步骤4-3～步骤4-5。
4-7终止迭代，得到最终重建图像
Figure PCTCN2019126882-appb-000075
以下我们通过对Shepp & Logan phantom模型的正弦图进行重建来验证实施方法的实用性和可靠性。Shepp & Logan phantom模型的真值图像如图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图像的投影数据，组成投影数据集
    Figure PCTCN2019126882-appb-100001
    所述投影数据集
    Figure PCTCN2019126882-appb-100002
    由对应各个角度方向下所采集得到的投影数据向量组成，所述投影数据向量维度为n且向量中每一元素值为对应探测器测得的投影数据，n为探测器个数；
    (2)分别建立低剂量模式和稀疏视角模式下的CT图像重建方程；
    低剂量模式下的CT图像重建方程表达式为：
    Figure PCTCN2019126882-appb-100003
    稀疏视角模式下的CT图像重建方程表达式为：
    Figure PCTCN2019126882-appb-100004
    其中：
    Figure PCTCN2019126882-appb-100005
    为CT图像数据集且维度为k，k为CT图像的像素点总数，
    Figure PCTCN2019126882-appb-100006
    中的每一元素值为待重建CT图像中对应像素点的X光线吸收系数，A为系统矩阵，β为给定的权重系数，|| || ATV表示各向异性全变分；
    (3)对两种模式下的CT图像重建方程进行预处理得到对应的目标函数，预处理过程包括引入变量
    Figure PCTCN2019126882-appb-100007
    和非奇异矩阵P、利用拉格朗日对偶将重建方程转换成鞍点问题；
    (4)根据实际情况选择对应模式下的目标函数，并对其优化求解以重建得到CT图像。
  2. 根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对低剂量模式下CT图像重建方程的预处理过程如下：
    首先，引入变量
    Figure PCTCN2019126882-appb-100008
    将低剂量模式下的CT图像重建方程改写为如下形式：
    Figure PCTCN2019126882-appb-100009
    其中：
    Figure PCTCN2019126882-appb-100010
    表示由CT图像数据集
    Figure PCTCN2019126882-appb-100011
    计算得到的正向投影数据；
    然后，引入非奇异矩阵P，将式(1.2)进一步改写为如下形式：
    Figure PCTCN2019126882-appb-100012
    式(1.3)对应的拉格朗日方程
    Figure PCTCN2019126882-appb-100013
    如下：
    Figure PCTCN2019126882-appb-100014
    其中：
    Figure PCTCN2019126882-appb-100015
    为拉格朗日乘子向量,T表示转置；
    则将式(1.3)转换成鞍点问题，具体表达形式如下：
    Figure PCTCN2019126882-appb-100016
    最后，通过极小化式(1.4)将变量
    Figure PCTCN2019126882-appb-100017
    去除，从而得到对应的目标函数如下：
    Figure PCTCN2019126882-appb-100018
    其中：·表示内积运算。
  3. 根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(3)中对稀疏视角模式下CT图像重建方程的预处理过程如下：
    首先，引入非奇异矩阵P，将稀疏视角模式下的CT图像重建方程改写为如下形式：
    Figure PCTCN2019126882-appb-100019
    式(2.2)对应的拉格朗日方程
    Figure PCTCN2019126882-appb-100020
    如下：
    Figure PCTCN2019126882-appb-100021
    则将式(2.2)转换成鞍点问题，从而得到对应的目标函数如下：
    Figure PCTCN2019126882-appb-100022
    其中：
    Figure PCTCN2019126882-appb-100023
    为拉格朗日乘子向量,T表示转置，β为给定的权重系数。
  4. 根据权利要求2所述的有限角度CT重建方法，其特征在于：对于低剂量模式下，所述非奇异矩阵P的计算表达式如下：
    P＝F -1R(ω)F
    Figure PCTCN2019126882-appb-100024
    其中：F和F -1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
  5. 根据权利要求3所述的有限角度CT重建方法，其特征在于：对于稀疏视角模式下，所述非奇异矩阵P的计算表达式如下：
    P＝F -1R(ω)F
    Figure PCTCN2019126882-appb-100025
    其中：F和F -1分别为一维傅里叶变换算子和傅里叶逆变换算子，ω为对应角度方向下所采集得到的投影数据向量经傅里叶变换后的频域变量，m为角度方向数量，τ为给定的参数。
  6. 根据权利要求1所述的有限角度CT重建方法，其特征在于：所述步骤(4)中采用交替投影接近算法对目标函数进行优化求解以重建得到CT图像，该算法中涉及到的各项参数包括总迭代次数N iter、TV去噪迭代次数N ATV、迭代终止阈值δ和∈、TV项权重系数β、参数α、参数τ、参数σ、ATV权重因子γ h和γ v；其中，总迭代次数N iter的取值范围为1～50000，TV去噪迭代次数N ARV的取值范围为1～10000，迭代终止阈值δ和∈的取值范围为10 -10～1，TV项权重系数β的取值范围为10 -6～1，参数α的取值范围为0.01～0.25，参数τ的取值范围为10 -3～10 3，参数σ的取值范围为10 -3～10 3，ATV权重因子γ h和γ v的取值范围为10 -6～1。
PCT/CN2019/126882 2019-01-24 2019-12-20 一种基于各向异性全变分的有限角度ct重建算法 WO2020151424A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/979,168 US11727609B2 (en) 2019-01-24 2019-11-20 Limited-angle CT reconstruction method based on anisotropic total variation

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910069082.8A CN109840927B (zh) 2019-01-24 2019-01-24 一种基于各向异性全变分的有限角度ct重建算法
CN201910069082.8 2019-01-24

Publications (1)

Publication Number Publication Date
WO2020151424A1 true WO2020151424A1 (zh) 2020-07-30

Family

ID=66884144

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/126882 WO2020151424A1 (zh) 2019-01-24 2019-12-20 一种基于各向异性全变分的有限角度ct重建算法

Country Status (3)

Country Link
US (1) US11727609B2 (zh)
CN (1) CN109840927B (zh)
WO (1) WO2020151424A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113706653A (zh) * 2021-09-23 2021-11-26 明峰医疗系统股份有限公司 基于AwTV的CT迭代重建替代函数优化方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109840927B (zh) * 2019-01-24 2020-11-10 浙江大学 一种基于各向异性全变分的有限角度ct重建算法
CN110717956B (zh) * 2019-09-30 2023-06-20 重庆大学 一种有限角投影超像素引导的l0范数最优化重建方法
CN111696166A (zh) * 2020-06-10 2020-09-22 浙江大学 基于fdk型预处理矩阵的圆周锥束ct快速迭代重建方法
CN112070856B (zh) * 2020-09-16 2022-08-26 重庆师范大学 基于非下采样轮廓波变换的有限角c型臂ct图像重建方法
CN113554729B (zh) * 2021-07-28 2022-06-07 中北大学 一种ct图像重建方法及系统
CN114549684B (zh) * 2022-03-05 2023-03-24 昆明理工大学 基于离散余弦变换和代数重建算法改进矿井无线电波透视成像重建方法
CN115249028B (zh) * 2022-05-13 2023-06-23 哈尔滨工业大学 一种基于稀疏正则化约束的盲解卷积无线通信信号重构方法
CN116509453B (zh) * 2023-06-29 2023-09-01 南京科进实业有限公司 一种基于超声波的胫骨骨密度检测系统及方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160063694A1 (en) * 2014-08-28 2016-03-03 Kabushiki Kaisha Toshiba Iterative reconstruction scheme for phase contrast tomography
CN106846427A (zh) * 2017-01-25 2017-06-13 浙江大学 一种基于重加权各向异性全变分的有限角度ct重建方法
CN108280859A (zh) * 2017-12-25 2018-07-13 华南理工大学 一种采样角度受限下的ct稀疏投影图像重建方法及装置
CN109840927A (zh) * 2019-01-24 2019-06-04 浙江大学 一种基于各向异性全变分的有限角度ct重建算法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102010034099B4 (de) * 2010-08-12 2017-04-06 Siemens Healthcare Gmbh Iterative Bildfilterung mit anisotropem Rauschmodell für ein CT-Bild
CN103679706B (zh) * 2013-11-26 2018-02-02 中国人民解放军第四军医大学 一种基于图像各向异性边缘检测的ct稀疏角度重建方法
CN104657950B (zh) * 2015-02-16 2017-05-03 浙江大学 一种基于Poisson TV的动态PET图像重建方法
US10140733B1 (en) * 2017-09-13 2018-11-27 Siemens Healthcare Gmbh 3-D vessel tree surface reconstruction
CN108550172B (zh) * 2018-03-07 2020-05-19 浙江大学 一种基于非局部特性和全变分联合约束的pet图像重建方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160063694A1 (en) * 2014-08-28 2016-03-03 Kabushiki Kaisha Toshiba Iterative reconstruction scheme for phase contrast tomography
CN106846427A (zh) * 2017-01-25 2017-06-13 浙江大学 一种基于重加权各向异性全变分的有限角度ct重建方法
CN108280859A (zh) * 2017-12-25 2018-07-13 华南理工大学 一种采样角度受限下的ct稀疏投影图像重建方法及装置
CN109840927A (zh) * 2019-01-24 2019-06-04 浙江大学 一种基于各向异性全变分的有限角度ct重建算法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113706653A (zh) * 2021-09-23 2021-11-26 明峰医疗系统股份有限公司 基于AwTV的CT迭代重建替代函数优化方法
CN113706653B (zh) * 2021-09-23 2023-12-08 明峰医疗系统股份有限公司 基于AwTV的CT迭代重建替代函数优化方法

Also Published As

Publication number Publication date
US11727609B2 (en) 2023-08-15
US20200402274A1 (en) 2020-12-24
CN109840927B (zh) 2020-11-10
CN109840927A (zh) 2019-06-04

Similar Documents

Publication Publication Date Title
WO2020151424A1 (zh) 一种基于各向异性全变分的有限角度ct重建算法
JP6855223B2 (ja) 医用画像処理装置、ｘ線コンピュータ断層撮像装置及び医用画像処理方法
Niu et al. Sparse-view x-ray CT reconstruction via total generalized variation regularization
Zhu et al. Improved compressed sensing‐based algorithm for sparse‐view CT image reconstruction
Ge et al. ADAPTIVE-NET: deep computed tomography reconstruction network with analytical domain transformation knowledge
Wang et al. ADMM-based deep reconstruction for limited-angle CT
Vandeghinste et al. Split-Bregman-based sparse-view CT reconstruction
Li et al. Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV)
Wang et al. An end-to-end deep network for reconstructing CT images directly from sparse sinograms
Lee et al. Interior tomography using 1D generalized total variation. Part II: Multiscale implementation
Cheslerean-Boghiu et al. WNet: A data-driven dual-domain denoising model for sparse-view computed tomography with a trainable reconstruction layer
Luo et al. Adaptive weighted total variation minimization based alternating direction method of multipliers for limited angle CT reconstruction
Cierniak An analytical iterative statistical algorithm for image reconstruction from projections
Feng et al. Dual residual convolutional neural network (DRCNN) for low-dose CT imaging
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
CN105844678A (zh) 基于全广义变分正则化的低剂量x射线ct图像重建方法
Bubba et al. Shearlet-based regularized reconstruction in region-of-interest computed tomography
Wu et al. Low dose CT reconstruction via L 1 norm dictionary learning using alternating minimization algorithm and balancing principle
Coban et al. Assessing the efficacy of tomographic reconstruction methods through physical quantification techniques
CN110853113A (zh) 基于bpf的tof-pet图像重建算法及重建系统
Gopi et al. Iterative computed tomography reconstruction from sparse-view data
Tiwari et al. A hybrid-cascaded iterative framework for positron emission tomography and single-photon emission computed tomography image reconstruction
Arya et al. Regularization based modified SART iterative method for CT image reconstruction
Adelman et al. Deformable registration and region-of-interest image reconstruction in sparse repeat CT scanning

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19912120

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19912120

Country of ref document: EP

Kind code of ref document: A1