CN113706653B - 基于AwTV的CT迭代重建替代函数优化方法 - Google Patents

基于AwTV的CT迭代重建替代函数优化方法 Download PDF

Info

Publication number
CN113706653B
CN113706653B CN202111111733.9A CN202111111733A CN113706653B CN 113706653 B CN113706653 B CN 113706653B CN 202111111733 A CN202111111733 A CN 202111111733A CN 113706653 B CN113706653 B CN 113706653B
Authority
CN
China
Prior art keywords
matrix
substitution function
awtv
term
iterative
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
CN202111111733.9A
Other languages
English (en)
Other versions
CN113706653A (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.)
Minfound Medical Systems Co Ltd
Original Assignee
Minfound Medical Systems Co Ltd
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 Minfound Medical Systems Co Ltd filed Critical Minfound Medical Systems Co Ltd
Priority to CN202111111733.9A priority Critical patent/CN113706653B/zh
Publication of CN113706653A publication Critical patent/CN113706653A/zh
Application granted granted Critical
Publication of CN113706653B publication Critical patent/CN113706653B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical 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
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • Theoretical Computer Science (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Surgery (AREA)
  • Biomedical Technology (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Mathematical Physics (AREA)
  • Pulmonology (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明涉及CT图像重建技术领域，尤其涉及一种基于AwTV的CT迭代重建替代函数优化方法；该方法基于AwTV正则化模型的代价函数的凸特性，构建替代函数，实现替代函数各个像素点相互分离，从而所有像素点的优化可以同步进行，借助GPU的并行计算能力，每次迭代耗时约1分钟；同时，由于数据拟合项和正则项一起优化，保证迭代过程单调收敛到最优值点，相对于凸集投影算法，提高了收敛速度，减少了迭代次数，使总的迭代耗时大大降低。

Description

基于AwTV的CT迭代重建替代函数优化方法
技术领域
本发明涉及CT图像重建技术领域，尤其涉及一种基于AwTV的CT迭代重建替代函数优化方法。
背景技术
计算机X射线断层扫描仪(computed tomography，CT)是利用X射线旋转照射被测物体，然后通过计算机处理获得物体断层图像的设备。其中，由处理扫描数据得到物体断层图像的过程称为CT重建。现有的CT重建算法可以分为两类：滤波反投影(filtering backprojection，FBP)和迭代重建。随着CT在临床中的使用越来越广泛和频繁，被扫描患者承受的X射线辐射也越来越多，为了降低受X射线辐射所引起的病变风险，低剂量CT受到越来越多的关注。当扫描剂量变低时，FBP重建的CT图像通常会存在非常严重的噪声和伪影，影响临床诊疗，而迭代重建能够避免这一问题。
迭代重建基本过程是，首先根据CT系统的特性，构建一个代价函数，然后计算出代价函数的最优值点(通常是最小值点)，该点即为迭代重建得到CT图像。代价函数的一个常见形式，公式如下：
；
其中，A表示系统矩阵，由CT系统的光路模型构建；y表示扫描得到的投影数据；D表示权重矩阵，通常是一个对角矩阵，对角线元素代表每个投影数据的权重系数；x表示CT图像，每个元素对应CT图像每个像素点的值；U(x)表示正则化模型，称为正则项；与之对应的代价函数的另一部分称为数据拟合项。基于上述代价函数，迭代重建得到的CT图像可以用以下公式表示：
。
由于直接计算代价函数L(x)的最小值点是非常困难的，因此，迭代重建的做法是通过迭代的过程，逐步逼近代价函数的最小值点，这一过程也称为优化。但优化也并非易事，尤其是当正则化模型较为复杂的时候。
现有的正则化模型有全变分、高斯马尔科夫随机场、广义高斯马尔科夫随机场以及q-GGMRF等。其中，TV在稀疏采样CT重建以及低剂量(低kv或者低毫安秒)CT重建中具有良好的降低伪影和噪声的能力，得到了广泛应用。但是，TV在降低伪影和噪声的同时，容易导致图像边缘过平滑，丢失图像细节，进而影响医生临床诊断。为了解决这一问题，论文“Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction.”提出了自适应加权全变分正则化，AwTV，其形式如下公式：
；
其中，N是总像素个数，n是每行像素个数，wj,j+n和wj,j+1称为自适应加权系数，wj,j+n和wj,j+1的公式如下：
；
；
其中，δ表示缩放因子，能够控制自适应权重wj,j+n和wj,j+1对TV的影响，进而控制图像平滑程度，当δ→∞时，AwTV等效于TV。
由于AwTV正则化模型的非线性，直接最小化上述代价函数是非常困难的。为此，论文“Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction.”采用凸集投影算法对代价函数进行迭代优化，得到代价函数最优值，即重建出高质量CT图像。
上述所采用的凸集投影算法分为两步：凸集投影和全变分最小化。其中，凸集投影相当于对代价函数中的数据拟合项进行迭代优化；全变分最小化相当于对代价函数中的正则项进行迭代优化。凸集投影采用同步代数迭代重建(simultaneous algebraicreconstruction technique，SART)算法迭代实现；全变分最小化通过梯度下降法(gradient descent，GD)实现。由于代价函数的数据拟合项和正则项分开优化，迭代过程出现震荡，导致迭代速度较慢，迭代次数较多，总的迭代耗时较长。如图1所示，自第k次迭代开始，首先进行凸集投影得到数据拟合项迭代重建结果然后对上述数据拟合项迭代重建结果/>进行全变分最小化迭代重建，得到第k+1次的迭代重建结果，不断重复上述两个步骤，直到收敛到最优值点X*。
由于数据拟合项迭代过程中没有考虑正则项的约束，相邻两次数据拟合项的迭代重建结果可能分别位于最优值点的两侧，比如，和/>导致迭代过程出现震荡，增加迭代重建收敛所需要的次数，需要300次以上的迭代次数才能使图像收敛到满足临床需求的图像质量。
发明内容
为解决上述问题，本发明的目的在于提供一种基于AwTV的CT迭代重建替代函数优化方法，以降低基于AwTV的CT迭代重建过程的计算耗时。
为了实现上述目的，本发明的技术方案如下：
一种基于AwTV的CT迭代重建替代函数优化方法，所述CT迭代重建替代函数L(x)包括正则项U(x)和数据拟合项，公式为：
其中A表示系统矩阵，由CT系统的光路模型构建；y表示扫描得到的投影数据；D表示对角矩阵，对角线元素代表每个投影数据的权重系数；x表示CT图像，每个元素对应CT图像每个像素点的值；β表示正则项的系数，为常数值；自适应加权全变分正则化AwTV的公式为：
其中，N是总像素个数，n是每行像素个数，wj,j+n和wj,j+1称为自适应加权系数，wj,j+n和wj,j+1的公式如下：
；
；
其中，δ表示缩放因子；
利用AwTV正则化模型的凸特性，构建如下像素可分离二次型替代函数：
其中表示第k次迭代所得图像xk在像素点j处的AwTV，此时的权重系数wj,j+n和wj,j+1采用当前迭代所得图像计算/>和/>作为下次迭代的权重估计值；
利用数据拟合项的凸特性，构建如下像素可分离二次型替代函数：
；
其中，Ai表示矩阵A的第i行元素；Aij表示矩阵A的第i行、第j列元素；di表示矩阵D的第i个对角线元素；yi表示扫描所得投影数据y的第i个元素；Ai+＝ΣjAij，表示对矩阵A的第j列元素进行求和，其中矩阵中i的上限值为矩阵的总行数，j的上限值为矩阵的总列数；
由数据拟合项的替代函数和正则项的替代函数组成新的CT迭代重建替代函数，并对任意像素点xj求导，同时考虑CT图像的非负约束，得到每个像素点xj的更新公式如下：
其中，M表示采集到的投影个数。
本发明的优点在于：
1.每个像素点的更新公式只与该像素点本身以及上一次迭代的结果有关，与其它像素点的当前迭代完全无关，因此，所有的像素点可以同时进行迭代更新；
2.由于数据拟合项和正则项一起优化，每步的迭代过程都是连续单调地收敛到最优值点，不会出现震荡情况，这能够有效提高收敛速度，降低迭代次数；
3.考虑到所有像素点可以同时进行迭代，借助GPU的并行计算能力，每次迭代耗时降低到约1分钟，与现有的凸集投影方法每次迭代耗时相似，但迭代次数可由原来的300次以上降低到约5到10次，总的迭代耗时大大降低。
附图说明
图1为背景技术中采用凸集投影算法时的迭代示意图；
图2为实施例中采用优化后的算法时的迭代示意图；
图3(a)采用FBP方法获得的重建CT图像；
图3(b)采用凸集投影方法获得的重建CT图像；
图3(c)采用实施例方法获得的重建CT图像；
具体实施方式
以下结合实施例对本发明作进一步详细描述。
本实施例提出一种基于AwTV的CT迭代重建替代函数优化方法，所述CT迭代重建替代函数L(x)包括正则项U(x)和数据拟合项，公式为：
；
其中A表示系统矩阵，由CT系统的光路模型构建；y表示扫描得到的投影数据；D表示权重矩阵，通常是一个对角矩阵，对角线元素代表每个投影数据的权重系数；x表示CT图像，每个元素对应CT图像每个像素点的值；自适应加权全变分正则化AwTV的公式为：
；
其中，N是总像素个数，n是每行像素个数，wj,j+n和wj,j+1称为自适应加权系数，wj,j+n和wj,j+1的公式如下：
；
；
其中,δ表示缩放因子；
其特征在于，利用AwTV正则化模型的凸特性，构建如下像素可分离二次型替代函数：
其中,
表示第k次迭代所得图像xk在像素点j处的AwTV，此时的权重系数wj,j+n和wj,j+1采用当前迭代所得图像计算和/>作为下次迭代的权重估计值；这样做，一方面能够简化正则化函数的复杂度，另一方面，权重系数的滞后一步方法不影响迭代结果的收敛性。
利用数据拟合项的凸特性，构建如下像素可分离二次型替代函数：
其中，Ai表示矩阵A的第i行元素；Aij表示矩阵A的第i行、第j列元素；di表示矩阵D的第i个对角线元素；yi表示扫描所得投影数据y的第i个元素；Ai+＝ΣjAij；
由数据拟合项的替代函数和正则项的替代函数组成新的CT迭代重建替代函数，并对任意像素点xj求导，同时考虑CT图像的非负约束，得到每个像素点xj的更新公式如下：
其中，M表示采集到的投影个数。
采用公式(7)同时计算出所有像素点的更新值，并不断重复上述过程即可实现并行的迭代优化。
由公式(7)可见，每个像素点的更新公式只与该像素点本身以及上一次迭代的结果有关，与其它像素点的当前迭代完全无关,因此，所有的像素点可以同时进行迭代更新。并且，由于数据拟合项和正则项一起优化，每步的迭代过程都是连续单调地收敛到最优值点，不会出现震荡情况，如图2所示。这能够有效提高收敛速度，降低迭代次数。
考虑到所有像素点可以同时进行迭代，借助GPU的并行计算能力，每次迭代耗时降低到约1分钟，与现有的凸集投影方法每次迭代耗时相似。同时，由于迭代次数由原来的300次以上降低到约5到10次，总的迭代耗时大大降低。
根据上述CT迭代重建替代函数优化方法，CT图像重建包括以下步骤：
(1).用FBP的重建图像初始化迭代重建图像x0；
(2).对于每个像素点，根据公式(3a)、公式(3b)和公式(5)，同时计算相应的Uj-n(xk)、Uj-1(xk)和Uj(xk)；
(3).根据公式(7)计算像素点的更新值，得到新的估计图像；
(4).重复步骤2-3，直到满足收敛条件。
应当说明，虽然上述公式针对的是二维CT图像的重建，但本专利所提出方案可以非常容易地推广到三维CT重建的情况。
上述实施例仅用于解释说明本发明的构思，而非对本发明权利保护的限定，凡利用此构思对本发明进行非实质性的改动，均应落入本发明的保护范围。

Claims (1)

1.一种基于AwTV的CT迭代重建替代函数优化方法，所述CT迭代重建替代函数L(x)包括正则项U(x)和数据拟合项，公式为：
其中A表示系统矩阵，由CT系统的光路模型构建；y表示扫描得到的投影数据；D表示对角矩阵，对角线元素代表每个投影数据的权重系数；x表示CT图像，每个元素对应CT图像每个像素点的值；β表示正则项的系数，为常数值；
自适应加权全变分正则化AwTV的公式为：
其中，N是总像素个数，n是每行像素个数，wj,j+n和wj,j+1称为自适应加权系数，wj,j+n和wj,j+1的公式如下：
；
；
其中，δ表示缩放因子；
其特征在于，利用AwTV正则化模型的凸特性，构建如下像素可分离二次型替代函数：
其中表示第k次迭代所得图像xk在像素点j处的AwTV，此时的权重系数wj,j+n和wj,j+1采用当前迭代所得图像计算/>和/>作为下次迭代的权重估计值；
利用数据拟合项的凸特性，构建如下像素可分离二次型替代函数：
；
其中，Ai表示矩阵A的第i行元素；Aij表示矩阵A的第i行、第j列元素；di表示矩阵D的第i个对角线元素；yi表示扫描所得投影数据y的第i个元素；Ai+＝ΣjAij，表示对矩阵A的第j列元素进行求和，其中矩阵中i的上限值为矩阵的总行数，j的上限值为矩阵的总列数；
由数据拟合项的替代函数和正则项的替代函数组成新的CT迭代重建替代函数，并对任意像素点xj求导，同时考虑CT图像的非负约束，得到每个像素点xj的更新公式如下：
其中，M表示采集到的投影个数。
CN202111111733.9A 2021-09-23 2021-09-23 基于AwTV的CT迭代重建替代函数优化方法 Active CN113706653B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111111733.9A CN113706653B (zh) 2021-09-23 2021-09-23 基于AwTV的CT迭代重建替代函数优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111111733.9A CN113706653B (zh) 2021-09-23 2021-09-23 基于AwTV的CT迭代重建替代函数优化方法

Publications (2)

Publication Number Publication Date
CN113706653A CN113706653A (zh) 2021-11-26
CN113706653B true CN113706653B (zh) 2023-12-08

Family

ID=78661524

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111111733.9A Active CN113706653B (zh) 2021-09-23 2021-09-23 基于AwTV的CT迭代重建替代函数优化方法

Country Status (1)

Country Link
CN (1) CN113706653B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108765514A (zh) * 2018-06-06 2018-11-06 上海交通大学 一种ct图像重建的加速方法及装置
WO2020151424A1 (zh) * 2019-01-24 2020-07-30 浙江大学 一种基于各向异性全变分的有限角度ct重建算法
CN112634388A (zh) * 2020-11-30 2021-04-09 明峰医疗系统股份有限公司 一种ct迭代重建代价函数的优化方法及ct图像重建方法、系统及ct

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108765514A (zh) * 2018-06-06 2018-11-06 上海交通大学 一种ct图像重建的加速方法及装置
WO2020151424A1 (zh) * 2019-01-24 2020-07-30 浙江大学 一种基于各向异性全变分的有限角度ct重建算法
CN112634388A (zh) * 2020-11-30 2021-04-09 明峰医疗系统股份有限公司 一种ct迭代重建代价函数的优化方法及ct图像重建方法、系统及ct

Also Published As

Publication number Publication date
CN113706653A (zh) 2021-11-26

Similar Documents

Publication Publication Date Title
Wu et al. DRONE: Dual-domain residual-based optimization network for sparse-view CT reconstruction
WO2021077997A9 (zh) 图像去噪的多生成器生成对抗网络学习方法
Yuan et al. SIPID: A deep learning framework for sinogram interpolation and image denoising in low-dose CT reconstruction
CN107481297B (zh) 一种基于卷积神经网络的ct图像重建方法
CN112102213B (zh) 低剂量ct图像处理方法、扫描系统及计算机存储介质
US20130051516A1 (en) Noise suppression for low x-ray dose cone-beam image reconstruction
CN110223255B (zh) 一种基于残差编解码网络的低剂量ct图像去噪递归方法
Zhou et al. DuDoUFNet: dual-domain under-to-fully-complete progressive restoration network for simultaneous metal artifact reduction and low-dose CT reconstruction
US20180089864A1 (en) Cbct image processing method
Gajera et al. CT-scan denoising using a charbonnier loss generative adversarial network
CN112435164A (zh) 基于多尺度生成对抗网络的低剂量ct肺部图像的同时超分辨率和去噪方法
WO2022000192A1 (zh) 一种ct图像的构建方法、ct设备以及存储介质
CN112634388B (zh) 一种ct迭代重建代价函数的优化方法及ct图像重建方法、系统及ct
WO2022091869A1 (ja) 医用画像処理装置、医用画像処理方法及びプログラム
Zhang et al. DREAM-Net: Deep residual error iterative minimization network for sparse-view CT reconstruction
Chan et al. An attention-based deep convolutional neural network for ultra-sparse-view CT reconstruction
Trung et al. Dilated residual convolutional neural networks for low-dose CT image denoising
Choi et al. Self-supervised inter-and intra-slice correlation learning for low-dose CT image restoration without ground truth
CN113706653B (zh) 基于AwTV的CT迭代重建替代函数优化方法
Chun et al. Joint image reconstruction and nonrigid motion estimation with a simple penalty that encourages local invertibility
CN112656438A (zh) 一种基于曲面全变差的低剂量ct投影域去噪及重建方法
Xia et al. Diffusion Prior Regularized Iterative Reconstruction for Low-dose CT
CN109658464B (zh) 基于加权核范数极小的稀疏角ct图像重建方法
Ikuta et al. A deep recurrent neural network with FISTA optimization for CT metal artifact reduction
CN116894783A (zh) 基于时变约束的对抗生成网络模型的金属伪影去除方法

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