CT重建项目考核

gaomingyang
36
2025-07-31

seu影像组考核任务的选定如下:CT重建呼吸心率监测算法脑疾病(老年痴呆症-基于医学影像)诊断算法复现和调研脑机接口算法复现和调研脑疾病(抑郁症-基于视频语音-CV算法相关)诊断算法复现和调研,五选一

CT重建

计算机断层扫描成像(CT)是临床医学中广泛使用的一种医学图像,它可以清晰地可视化人体内部精细结构细节。在临床操作中,为防止患者暴露在高辐射X射线束下引起组织受损,通常最小化X射线以获得CT图像,但会导致成像质量严重下降。为解决上述矛盾,如何重建出符合临床需求的CT图像是国内外研究者广泛关注的、具有挑战性的难点问题。

202507231310922.jpg

1 CT成像概述

CT图像重建是将从成像设备收集的原始数据形成可解释图像的过程。在已知一组测量值的前提下,目标是确定影响接收器收集信号的原始图像结构,这一过程被称为逆问题。设y是一组原始采集的传感器测量值,在采集过程中收到一些固有噪声N的影响,目标是恢复空间域的图像x,将这一过程用下列公式描述:

y=A(F(x),N)

其中,F(·)代表成像物理结构进行建模的正向操作运算,一般包括拉东变换(Radon)或傅里叶变换,它可以是线性操作,也可以是非线性操作,具体取决于成像模式。A表示噪声和信号之间的相互作用。图像重建还是一个不确定的问题,因为测量值(M)往往比未知数(N)少得多。

CT图像重建方法是低剂量CT图像到正常剂量CT图像的映射。CT图像的稀疏视图、有限角采样、降低辐射剂量通常会减小测量信号y的大小,同时增加其稀疏性和噪声水平,从而增加重建问题的不适定性和复杂性。这就提出了对具有高特征提取能力的复杂重建算法的需求,重建算法要最大限度地利用收集的信号以及先验知识,捕捉特定于模态的成像特征。

2 传统方法

传统CT重建方法大致可分为正弦图域滤波、迭代重建和图像域恢复三类,这三类方法虽然在重建精度和伪影减少方面有了显著的改善,但它们仍然存在一些缺点,在实际临床应用中往往受到限制。

2.1 正弦图域滤波方法

基于正弦图域滤波的重建方法主要目的是从低剂量X射线束获得的CT原始数据中滤除噪声,通过直接对反投影前形成的原始投影数据进行重建,可以准确地计算噪声统计量并进行有效重建。

  • 优点:

    • 计算成本低(毫秒量级)

    • 可以在无噪声、全采样或所有角度投影的假设下生成良好的图像质量

  • 缺点:

    • 只考虑成像系统的几何形状和采样特性,而忽略系统物理特性和测量噪声的细节

    • 当处理有噪声或不完整的测量数据时,例如降低测量采样率,重建结果会随着信号变弱而严重退化,并且无法恢复信号中的缺失信息,从而导致诊断性能受损

    • 重建过程预测数据过分依赖CT设备供应商的完备数据以及CT扫描仪无法获得的原始数据,这些数据大多数不能公开访问

  • 典型方法:结构自适应滤波、双边滤波、惩罚似然法、加权最小二乘算法

2.2 迭代重建(IR)方法

迭代重建方法基于成像系统的物理、传感器和噪声统计的更复杂的模型,结合传感器域(原始测量数据)中数据的统计特性、图像域中的先验信息,有时还将成像系统的参数组合到其损失函数中。

  • 优点:

    • 对噪声和不完整数据表示问题具有更好的鲁棒性

    • 与正弦图域滤波方法相比提高了准确率并减少了伪影

  • 缺点:

    • 迭代重建技术往往是特定于供应商的,扫描仪几何形状和校正步骤的细节对用户和其他供应商不开放

    • 每次迭代都需要投影和反投影操作的负载,迭代重建技术相关的计算开销很大

    • 重建质量高度依赖于正则化函数和相关的超参数设置,需要手动调整才能达到较优的效果

  • 典型方法:基于全变分(total variation,TV)的先验、TV的变体、非局部均值、小波方法、字典学习

2.3 图像后处理方法

图像后处理重建算法(基于图像域的重建)直接应用于CT图像,而不是原始投影数据。

  • 优点:

    • 相比于迭代重建方法,重建速度快

    • 不需要供应商提供原始数据,可以很轻松地集成到CT设备的工作流程中

  • 缺点:

    • 存在模糊CT图像的结构信息等问题

    • 算法实现过程中由于噪声的非均匀性而无法计算其统计量,削弱了CT重建的准确性

  • 典型方法:非局部均值过滤方法、基于字典学习的方法、块匹配算法、基于统计的算法

3 深度学习(DL)方法

在CT图像重建领域,深度模型可以捕获高级特征,体现了它在整个数据驱动学习过程中学习CT图像上的不确定噪声分布能力,此外,数据驱动学习方法可以有效地适应任何噪声类型。因此,深度学习方法可以显著提高CT图像重建的整体性能。基于深度学习的CT图像重建方法可划分为四个子类,即投影域CT图像重建、图像域CT图像重建、双域网络CT图像重建和直接映射CT图像重建,如下图所示。

202507231400509.jpg

3.1 投影域CT图像重建

投影域CT图像重建问题可以表述为投影域图像利用深度学习方法从不完整的数据表示(低剂量、数据稀疏、有限角)到完整数据表示(正常剂量)的回归问题。这一阶段的主要目标是利用深度学习模型估计在信号采集阶段没有采集到的缺失部分,以便将更完整的信号信息输入到滤波反投影层进行重建过程。

  • 优点:

    • 减少投影域中信号损失

    • 显著提高重建图像质量

  • 缺点:

    • 对正弦图的内在一致性很敏感,任何对正弦图的不当操作都可能在整个重建图像上引入额外的伪影

    • 深度学习方法提取的特征仅限于投影域,对于完备的投影数据在图像域的重建过程仍然存在缺陷

  • 典型方法:CNN、自编码神经网络、U-Net、GAN等

3.2 图像域CT图像重建

图像域CT图像重建的任务是学习低质量重建图像和高质量重建图像之间的映射。

  • 优点:

    • 适用于处理质量相对较好的图像重建

    • 适用于去除投影域图像中噪声伪影的图像域

  • 缺点:

    • 不能保证采样的正弦图数据得到保留

    • 初始重建中丢失的信息很难通过后处理过程进行恢复

  • 典型方法:CGAN、PatchGAN、ResNet、DenseNet、VGG等

3.3 双域网络CT图像重建

基于双域网络的CT图像重建方法同时利用了投影域和图像域数据,通过原始正弦图和相应的FBP重建图像信息互补,进一步提升了重建效果。

  • 优点:

    • 具有较高的重建精度

    • 可以解决图像重建领域以外的超分辨率、灌注CT反卷积等领域的优化和逆问题

  • 缺点:

    • 训练需要更大的GPU内存

    • 重建过程更加复杂,参数相对较多

  • 典型方法:专家评估、流形学习、自学习原对偶、小波变换等

3.4 直接映射CT图像重建

直接映射CT图像重建通过学习正弦图和CT图像空间之间的映射,同时近似逆问题的基本物理模型,使用深度神经网络直接从投影域数据解码为CT图像。

  • 优点:

    • 受益于深度学习模型的多级抽象和自动特征提取能力,可以得到高定量精度重建图像

  • 缺点:

    • 对大数据极度依赖

    • 需要大量计算计算资源

    • 专用神经网络仅限于一种特定的重建几何,并不广泛适用于各种扫描仪体系结构和扫描协议

  • 典型方法:流形近似自变换、编解码网络、级联反投影等

4 常用数据集及损失函数

4.1 常用数据集

正常剂量和低剂量的配对CT数据集对于深度学习模型的训练和验证至关重要。在临床操作中,对患者的重复扫描是获得成对数据的唯一方式。然而,在临床实践中,这种操作是不被允许的,当患者长期暴露在X射线的辐射中,会对患者身体造成不可逆的损伤。此外,CT图像的正弦图数据是特定于供应商的,不允许从第三方提取。有监督的深度学习模型中,必须有成对的NDCT和LDCT图像,常用的解决方案是将泊松噪声和高斯噪声添加到从NDCT获得的正弦图中模拟生成不同剂量的LDCT图像。

数据集名称

解刨学部位

数据集介绍

NBIA/NCIA数据集

包括胸腔在内的很多器官

国家生物医学档案馆共有7015张256x256大小的NDCT图像

AAPM-Mayo数据集

腹部

Mayo诊所提供的真实临床AAPM低剂量CT挑战数据集,包括来自10名患者的2378张3mm全剂量和四分之一剂量CT图像,大小为512x512

猪仔数据集

全身

在4个剂量水平下获得图像,每个剂量水平由850张512x512大小的图像组成

3D-IRCADb

不同器官

10例患者的1375张NDCT图像

Luna-16

肺部

888个带有标注的临床NDCT扫描数据

COVID-19

肺部

来自9个医学中心的1141例胸部容积CT检查,312个数据被标记为PCR阳性

为了解决数据量不足的问题,研究人员常利用旋转和翻转来增加训练数据集中的样本数量,与数据增强相比,基于图像块的训练加快了网络收敛,不仅有助于增强局部区域中感知方差的检测,同时增加训练样本的数量。

4.2 损失函数

损失函数用来评估数据在特定深度学习模型中的建模效果,通常由一个或多个函数组成,对重建的最终图像质量有很大的影响。

损失函数名称

具体描述

优势

L1损失

计算两幅图像的平均绝对误差

保持图像的对称性、凸性、减少模糊,但会形成块状伪影

L2损失

计算真实图像和生成图像之间的均方误差损失

确保生成图像与真实图像之间的可微性、对称性和凸性

对抗性损失

包括生成器损失和鉴别器损失

提高训练过程的稳定性

基于Wasserstein距离的对抗性损失

基于推土机距离确定GAN损失

利用训练的收敛性,提高训练过程的稳定性

循环损失

允许噪声图像和去噪图像之间一一对应的损失

加强了循环GAN模型的逆关系

最小二乘损失

最小二乘GAN模型

提高训练过程的稳定性

结构性损失

在图像中找到最佳视觉模式的损失

可以去模糊,并保留图像结构和纹理细节

感知损失

决定特征相似性的损失,基于VGG19等预训练的网络中提取特征进行计算

保留结构和纹理信息

多尺度损失

在不同尺度上进一步计算结果和背景信息的损失

在不同比例上保存结构和上下文信息

梯度损失

决定高层特征信息的损失

利用训练过程的融合实现高层次特征的提取

细节损失

真实图像的高频子带与生成图像的高频子带之间的损失

保留边缘信息

注意力损失

注意力块产生的图像与先验/注意知识之间的损失

创建最佳的注意力映射

清晰度损失

评估待处理图像中表示的特征清晰度损失的损失

确保特征的清晰度

5 文献复现——Structure-Aware Sparse-View X-ray 3D Reconstruction

SAX-NeRF提出了一种用于稀疏视角下 X 光三维重建的NeRF方法。具体而言,主要做两个任务。一是X光的新视角合成 (Novel View Synthesis, NVS),二是CT重建,可以简单理解为体密度的重建。具体贡献如下:

  • 提出了一套全新的能够同时做X光新视角合成与CT成像的NeRF框架,名为SAX-NeRF。该框架的训练不需要用CT作为监督信号,只使用X光片即可。

  • 设计了一种新的分段式Transformer,名为Lineformer,可以捕获成像物体在三维空间中的复杂的内部结构。Lineformer是首个将 Transformer应用于X光渲染的Transformer。

  • 提出了一种新型的射线采样策略,名为MLG sampling,可以从X光片上提取出局部和全局的信息。

  • 搜集了首个大规模的X光三维重建数据集,涵盖医疗、生物、安检、工业领域。同时,设计的算法在这个数据集上取得了当前最好效果,在X光新视角合成和CT重建两大任务上比之前的最好方法要高出12.56和2.49dB。

由于时间以及后续作业原因,仅将该文献提供的示例代码环境配好并跑通,尚未在自己的数据集上进行训练。

6 参考文献

[1] Ahishakiye E, Bastiaan Van Gijzen M, Tumwiine J, et al. A survey on deep learning in medical image reconstruction[J]. Intelligent Medicine, 2021, 1(03): 118-127.

[2] Zhang S, Xia Y. CT image reconstruction algorithms: A comprehensive survey[J]. Concurrency and Computation: Practice and Experience, 2021, 33(8): e5506.

[3] Ben Yedder H, Cardoen B, Hamarneh G. Deep learning for biomedical image reconstruction: A survey[J]. Artificial intelligence review, 2021, 54(1): 215-251.

[4] McLeavy C M, Chunara M H, Gravell R J, et al. The future of CT: deep learning reconstruction[J]. Clinical radiology, 2021, 76(6): 407-415.

[5] Szczykutowicz T P, Toia G V, Dhanantwari A, et al. A review of deep learning CT reconstruction: concepts, limitations, and promise in clinical practice[J]. Current Radiology Reports, 2022, 10(9): 101-115.

[6] Li D, Ma L, Li J, et al. A comprehensive survey on deep learning techniques in CT image quality improvement[J]. Medical & Biological Engineering & Computing, 2022, 60(10): 2757-2770.

[7] 李青,李润睿,强彦,等.人工智能在医学CT图像重建中的研究进展[J].太原理工大学学报,2023,54(01):1-16.DOI:10.16355/j.cnki.issn1007-9432tyut.2023.01.001.

[8] Wang T, Xia W, Lu J, et al. A review of deep learning CT reconstruction from incomplete projection data[J]. IEEE Transactions on Radiation and Plasma Medical Sciences, 2023, 8(2): 138-152.

[9] Koetzier L R, Mastrodicasa D, Szczykutowicz T P, et al. Deep learning image reconstruction for CT: technical principles and clinical prospects[J]. Radiology, 2023, 306(3): e221257.

[10] Cai Y, Wang J, Yuille A, et al. Structure-aware sparse-view x-ray 3d reconstruction[C]//Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2024: 11174-11183.

作业1

这个试题将帮你入门断层成像(Computed Tomography,简称CT)图像重建算法的原理。

1 Radon transform(拉东变换)

CT原始数据可近似为“拉东(Radon)变换”。拉东变换是一个积分变换,它将定义在二维平面上的一个函数f(x,y) 沿着平面上的任意一条直线做线积分,相当于对函数f(x,y)做CT扫描。对于直线

x\cos\theta+y\sin\theta=\rho

沿着这一条直线的拉东变换定义如下:

R\left( \rho ,\theta \right) =\int_{l:x\cos \theta +y\sin \theta =\rho}{f\left( x,y \right) dl}

R\left( \rho ,\theta \right) =\iint_{R^2}{f\left( x,y \right) \delta \left( x\cos \theta +y\sin \theta -\rho \right) dxdy}

上述两个定义是等价的。

例题:

对于一个边长为1的方形(内部值为1),如下图所示,沿着平行于边且与方形相交的直线的Radon变换值为1,沿着如图所示的斜着的方向变换值为\frac{\sqrt{2}}{2}

202507241440336.png

习题1

假设图像f(x,y)的定义发生改变,内部不再是均匀的,具体如下图,求解下述问题。

202507241441266.png

  1. 1、沿着图中蓝色实线的拉东变换值是多少?

  2. 2、沿着方形斜对角线的红色实线(直线\frac{\sqrt{2}}{2}x+\frac{\sqrt{2}}{2}y=0)拉东变换值是多少?

  3. 3、(此题较难)沿着图中距离原点距离为a且平行于方形斜对角线的红色虚线,其拉东变换值是?请分别为a\in \left[ \frac{\sqrt{2}}{4},\frac{\sqrt{2}}{2} \right]a\in \left[0,\frac{\sqrt{2}}{4} \right]两种情况进行讨论。

习题2

习题2.png

  1. 1、计算上图所示的f(x,y)(一个半径为1的圆,中心位于原点,圆内值为1)的拉东变换,即对于任意\rho\theta,计算R\left( \rho ,\theta \right) =\int_{l:x\cos \theta +y\sin \theta =\rho}{f\left( x,y \right) dl}。做出给定角度\thetaR(\rho,\theta)随着\rho变化的曲线;以及给定\rhoR(\rho,\theta)随着\theta变化的曲线。提示:结合f(x,y)中心对称特性思考该问题。

  2. 2、现实中我们只能在有限的数据点测量f(x,y)的拉东变换,按照下列标准对R(\rho,\theta)进行采样:\theta 的取值范围是1到360度,以1度为间隔;\rho取值范围是[-5,5],以0.01为间隔。以\rho方向为列方向,\theta为行方向,画出以此采样后所构成的R(\rho,\theta)二维图像。

2 图像重建的解析算法

在实际应用中,R(\rho,\theta)是测量获得的图像,而f(x,y)是未知的,我们需要通过算法从R(\rho,\theta)反推f(x,y),这是一个很有意思的逆问题求解过程。对于图像像素较少的情形,我们可以通过解析方法进行图像重建。

习题3

假设我们待重建的图像维度为2x2,共包含4个像素。该图像沿着如下图中5条直线的拉东变换值为已知量(已在图中标明),我们希望通过拉东变换的值求解图中的像素值(x_{11},x_{12},x_{21},x_{22})

  1. 1、将拉东变换中的线积分近似为沿着直线路径像素值的求和,写出x_{11},x_{12},x_{21},x_{22}满足的线性方程组(共五个方程)。

  2. 2、利用这五个方程,求解x_{11},x_{12},x_{21},x_{22}的值。

习题3.png

3 Fourier Slice Theorem(傅里叶中心切片定理)

在实际场景中,我们需要重建的图像通常包含超过几十万个像素,意味着需要求解一个几十万元一次方程,这个计算的开销是非常大的。一种更为高效的算法称为滤波反投影算法,其中傅里叶中心切片定理是这个算法的核心:对R(\rho,\theta)沿着\rho方向做傅里叶变换

\begin{aligned} \tilde{R}\left( k,\theta \right) &= \int R\left( \rho ,\theta \right) e^{-i2\pi \rho k} d\rho \\ &= \int \iint_{R^2} f\left( x,y \right) \delta \left( x\cos \theta +y\sin \theta -\rho \right) e^{-i2\pi \rho k} dxdy d\rho \\ &= \iint_{R^2} f\left( x,y \right) dxdy \int \delta \left( x\cos \theta +y\sin \theta -\rho \right) e^{-i2\pi \rho k} d\rho \\ &= \iint_{R^2} f\left( x,y \right) e^{-i2\pi k\left( x\cos \theta +y\sin \theta \right)} dxdy \\ &= \tilde{f}\left( k\cos \theta ,k\sin \theta \right) \end{aligned}

其中,\tilde{f}\left( k_x,k_y \right) =\iint_{R^2}{f\left( x,y \right) e^{-i2\pi \left( xk_x+yk_y \right)}dxdy}是f(x,y)的傅里叶变换。

上述推导为傅里叶中心切片定理:对R(\rho,\theta)沿着\rho方向做傅里叶变换,即是原图f(x,y)的傅里叶变换\tilde{f}\left( k_x,k_y \right)过原点沿着\theta方向的切片,如下图所示:

傅里叶中心切片定理.png

因此通过计算各个角度的\tilde{R}\left( k,\theta \right),我们可以得到\tilde{f}\left( k_x,k_y \right),进而通过二维傅里叶逆变换还原原始图像。

习题4(此题较难)

利用这个定理从你计算所得的R(\rho,\theta)中复原f(x,y)。

  1. 1、计算R(\rho,\theta)对于给定\theta的沿着\rho的傅里叶变换\tilde{R}\left( k,\theta \right)。此变换较难计算解析解,可使用问题2.2中的采样结果进行离散傅里叶变换。对于给定\theta,做出\tilde{R}\left( k,\theta \right)随着k变化的曲线(包括实部和虚部)。

  2. 2、对于给定\left( k_x,k_y \right),结合傅里叶中心切片定理以及\tilde{R}\left( k,\theta \right)的值,获得此点对应的\tilde{f}\left( k_x,k_y \right)值。此步骤涉及插值。画出\tilde{f}\left( k_x,k_y \right)对应的图像。

  3. 3、对\tilde{f}\left( k_x,k_y \right)做二维傅里叶反变换获得f(x,y),画出f(x,y)的图像,观察图像的形状和值,是否复原了原图。

提示:

  1. 1、结合f(x,y)中心对称特性可大幅度简化计算。

  2. 2、离散傅里叶变换注意计算变换后数据点对应的频率值。

  3. 3、\left( k_x,k_y \right)的取样需要取均匀的数据点,建议取值范围是-50到50(Nyquist频率),间隔为0.1。


习题解答

习题1

(1)直线方程:水平线y=\rho,对应\theta =90°\rho \in \left[ -0.5,0.5 \right]。故

R\left( \rho ,90° \right) =\int_{-0.5}^{0.5}{\left( 1-2\left| x \right| \right) dx}=2\int_0^{0.5}{\left( 1-2x \right) dx}=0.5

(2)直线方程:\frac{\sqrt{2}}{2}x+\frac{\sqrt{2}}{2}y=0\Rightarrow y=-x,对应\theta =45°\rho=0

参数化:令x=t,y=-t,则

dl=\sqrt{dx^2+dy^2}=\sqrt{2}t,t\in \left[ -0.5,0.5 \right]

R\left( 0,45° \right) =\int_{-0.5}^{0.5}{f\left( t,-t \right) \sqrt{2}tdt}=2\sqrt{2}\int_0^{0.5}{\left( 1-2t \right) dt}=\frac{\sqrt{2}}{2}

(3)直线方程:x+y=-\sqrt{2}a,对应\theta =45°\rho=a

参数化:令x=t,y=-\sqrt{2}a-t,则

dl=\sqrt{dx^2+dy^2}=\sqrt{2}t,t\in \left[ 0.5,0.5-\sqrt{2}a \right]

① 当a\in \left[0,\frac{\sqrt{2}}{4} \right]时,需分段

\begin{aligned} R\left( a,45° \right) &= \sqrt{2}\int_{-0.5}^{0.5-\sqrt{2}a}{\left( 1-2\left| t \right| \right) dt} \\ &= \sqrt{2}\left[ \int_{-0.5}^0{\left( 1+2t \right) dt}+\int_0^{0.5-\sqrt{2}a}{\left( 1-2t \right) dt} \right] \\ &= \frac{\sqrt{2}}{2}-2\sqrt{2}a^2 \end{aligned}

② 当a\in \left[ \frac{\sqrt{2}}{4},\frac{\sqrt{2}}{2} \right]时,无需分段

\begin{aligned} R\left( a,45° \right) &= \sqrt{2}\int_{-0.5}^{0.5-\sqrt{2}a}{\left( 1-2\left| t \right| \right) dt} \\ &= \sqrt{2}\int_{-0.5}^{0.5-\sqrt{2}a}{\left( 1+2t \right) dt} \\ &= 2\sqrt{2}a^2-4a+\sqrt{2} \end{aligned}

习题2

(1)依题意,f(x,y)沿l的拉东变换值即为直线与圆相交线段长度,即

R\left( \rho ,\theta \right) =\left\{ \begin{array}{l} 2\sqrt{1-\rho ^2}\\ 0\\ \end{array}\begin{array}{c} ,\left| \rho \right|\le 1\\ ,\left| \rho \right|>1\\ \end{array} \right.

曲线绘制:

  • 固定\theta:由于结果与θ无关,这条曲线对于任何角度都是一样的。它是一个定义在[-1, 1]上的半圆形函数。在\rho=0(直线穿过圆心)时,取得最大值2(即圆的直径);在\rho=±1(直线与圆相切)时,值为0。

251844b0854b287147572a17f98cce9e.png

  • 固定\rho:对于一个固定的\rho,变换值R(\rho,\theta)是一个常数C,不随\theta的变化而变化。在图上表现为一条水平直线,体现了圆的旋转不变性。

3df62c787e01d0b1ce3fb545a308ee2a.png

(2)采样后所构成的R(\rho,\theta)二维图像如下:

fe8105ed7bd1c9827a16ef495bb744d4.png

习题3

(1)依题意,有

\left\{ \begin{array}{l} x_{11}+x_{12}=3\\ x_{21}+x_{22}=8\\ x_{11}+x_{21}=4\\ x_{12}+x_{22}=7\\ x_{11}+x_{22}=6\\ \end{array} \right.

(2)求解第一问方程,可得

\left\{ \begin{array}{l} x_{11}=1\\ x_{12}=2\\ x_{21}=3\\ x_{22}=5\\ \end{array} \right.

习题4

(1)由习题2可知,R\left( \rho ,\theta \right) =\left\{ \begin{array}{l} 2\sqrt{1-\rho ^2}\\ 0\\ \end{array}\begin{array}{c} ,\left| \rho \right|\le 1\\ ,\left| \rho \right|>1\\ \end{array} \right.,则

\begin{aligned} \tilde{R}\left( k,\theta \right) &= \int{R\left( \rho ,\theta \right) e^{-i2\pi \rho k}d\rho} \\ &= \int_{-1}^1{2\sqrt{1-\rho ^2}e^{-i2\pi \rho k}d\rho} \\ \end{aligned}

由于上式较难求解析解,利用习题2.2的结果进行离散傅里叶变换。

① 离散采样数据:

\rho _i=\rho _{\min}+i\varDelta \rho =-5+0.01\times i,i=0,1,\cdots ,1000\left( \text{共1001个采样点} \right)
\theta _j=j°,j=1,2,\cdots ,360\left( \text{共360个角度} \right)

故采样得到的二维数据矩阵:

R\left[ i,j \right] =R\left( \rho _i,\theta _j \right)

对每一固定角度\theta_j,对其对应的采样序列R[i,j](即固定行j的数据)做一维离散傅里叶变换:

\tilde{R}\left[ m,j \right] =\sum_{i=0}^{N_{\rho}-1}{R\left[ i,j \right] \cdot e^{-i2\pi \frac{mi}{N_{\rho}}}},m=0,1,\cdots ,N_{\rho}-1

其中,N_\rho=1001\rho方向的采样点数,m是频率索引。

② 离散频率映射:离散频率对应连续频率k_m

k_m=\frac{m}{N_{\rho}\varDelta \rho},\text{若}m\le \frac{N_{\rho}}{2}

对于m > \frac{N_{\rho}}{2},对应负频率:

k_m=\frac{m-N_{\rho}}{N_{\rho}\varDelta \rho}

因此,频率范围约为

k_m\in \left[ -\frac{1}{2\varDelta \rho},\frac{1}{2\varDelta \rho} \right] =\left[ -50,50 \right]

综上,

\tilde{R}\left( k_m,\theta _j \right) =\sum_{i=0}^{N_{\rho}-1}{R\left( \rho _i,\theta _j \right) e^{-i2\pi \frac{mi}{N_{\rho}}}},m=0,1,\cdots ,N_{\rho}-1
k_m=\left\{ \begin{array}{l} \frac{m}{N_{\rho}\varDelta \rho},0\le m\le \frac{N_{\rho}}{2}\\ \frac{m-N_{\rho}}{N_{\rho}\varDelta \rho},\frac{N_{\rho}}{2}<m < N_{\rho}\\ \end{array} \right.

绘图时可采用python中的scipy.fft函数包进行傅里叶变换:

6b0485e65e077785925013a27b0150d9.png

4b9587c0b2b9fdf40d809aa49efe40ae.png

(2)根据傅里叶中心切片定理,

\tilde{R}\left( k,\theta \right) =\tilde{f}\left( k_x,k_y \right) ,\ \text{其中,\ }k_x=k\cos \theta ,\ k_y=k\sin \theta

① 构造频率平面均匀采样点:

k_x=k_{x,p}=k_{x,\min}+p\varDelta k_x,p=0,1,\cdots ,N_x-1
k_y=k_{y,q}=k_{y,\min}+q\varDelta k_y,q=0,1,\cdots ,N_y-1

其中,k_x,k_y\in \left[ -50,50 \right] ,\varDelta k_x=\varDelta k_y=\varDelta k=\frac{1}{N_{\rho}\varDelta \rho}=\frac{1}{1001\times 0.01}\approx 0.1

② 极坐标变换

对频率平面每一点(k_x,k_y),计算其对应的极坐标:

k=\sqrt{k_{x}^{2}+k_{y}^{2}}
\theta =\arctan ^{-1}\frac{k_y}{k_x}

③ 插值操作

由于\tilde{R}\left( k,\theta \right)仅在离散的k_m\theta_j处有定义,仅需对频率和角度两个维度插值:

\tilde{f}\left( k_x,k_y \right) =\tilde{R}\left( k,\theta \right) \approx Interp_{k,\theta}\left( \tilde{R}\left( k_m,\theta _j \right) \right)

具体来说,

  • 径向频率插值:对固定角度θ_j,在离散频率k_m之间插值,得到任意k处的值。

  • 角度插值:对固定频率k_m,在离散角度\theta_j之间插值,得到任意θ处的值。

常用的插值方法有线性插值、多项式插值、双线性插值等。

令插值函数为\mathcal{I},则

\tilde{f}\left( k_x,k_y \right) =\mathcal{I}\left( \tilde{R}\left( k_m,\theta _j \right) ;k=\sqrt{k_{x}^{2}+k_{y}^{2}},\theta =\arctan ^{-1}\frac{k_y}{k_x} \right)

通过线性插值得到的\tilde{f}\left( k_x,k_y \right)对应图像:

ac70395968ba08e2b16297512a69d941.png

(3)空间域函数通过频谱的二维逆傅里叶变换得到:

f\left( x,y \right) =\iint{\tilde{f}\left( k_x,k_y \right) e^{i2\pi \left( k_xx+k_yy \right)}dk_xdk_y}

给定均匀采样的\tilde{f}\left[ p,q \right] =\tilde{f}\left( k_{x,p},k_{y,p} \right),其中

k_{x,p}=k_{x,\min}+p\varDelta k_x,\ p=0,1, \cdots ,N_x-1
k_{y,q}=k_{y,\min}+q\varDelta k_y,\ q=0,1,\cdots ,N_y-1

空间坐标对应的采样点:

x_m=x_{\min}+m\varDelta x,\ m=0,1,\cdots ,N_x-1
y_n=y_{\min}+n\varDelta y,\ n=0,1,\cdots ,N_y-1

其中,分辨率由采样频率范围和步长决定:

\varDelta x=\frac{1}{N_x\varDelta k_x},\ \varDelta y=\frac{1}{N_y\varDelta k_y}

故离散逆傅里叶变换表达式为

f\left[ m,n \right] =\sum_{p=0}^{N_x-1}{\sum_{q=0}^{N_y-1}{\tilde{f}\left[ p,q \right] e^{i2\pi \left( k_{x,p}x_m+k_{y,q}y_n \right)}\varDelta k_x\varDelta k_y}}

复原图像:

4daaa31caabefb03afc35fcdcc3239c9.png

复原图像与原始单位圆函数差别较大,分析产生差异的主要原因如下:

  • 采样的ρ∈[−5,5],但实际Radon投影非零区间是[−1,1],数据中大部分ρ超出支持区间,且值为零。采样区间过宽导致FFT时有效信号占比较小,频谱数值被稀释,进而导致复原图像幅值非常低且失真。

  • 构造的频率网格(k_x,k_y)∈[−50,50]范围较大,步长为0.1,频率采样和插值区域过大,导致频率分辨率不足,Radon数据有效频率范围较窄,插值时很多点插到零或接近零数据,进而导致频谱整体幅度较低。

  • 标准的Radon反变换(滤波反投影)需要对傅里叶变换后的数据乘以滤波函数,以补偿Radon变换中频率响应不均匀。而代码中直接对插值后的频谱做逆傅里叶变换,没有滤波,导致高频成分未被适当增强,造成复原图像模糊且幅值较低。

作业2

1、简述自然图像和医学图像的区别,以及使用深度学习算法处理自然图像和医学图像的异同。

自然图像和医学图像的区别:

  • 特征层面:自然图像大多为自然光成像,具有丰富的色彩与多样的纹理,对比度高;而医学图像以单通道灰度图为主,特征比较单一。

  • 噪声层面:在自然图像中的噪声分布在大多数情况下可以认为是均匀的,一般近似为高斯噪声;而在医学图像中由于光源单一以及探测手段等因素的影响,往往会造成噪声分布复杂,往往认为是一种泊松噪声。

  • 数据层面:自然图像数据集规模庞大,且易获取;医学图像受隐私、标注成本和设备限制,数据量小且标注稀疏。

  • 应用层面:自然图像任务侧重通用性,如分类、检测等任务;医学图像需精确识别解剖结构或病灶,对错误容忍度极低。

深度学习算法处理自然图像和医学图像的异同:

  • 相同点:

    • 均使用CNN/Transformer等基础架构

    • 均需用到大规模数据进行训练

    • 均需数据增强等数据处理操作

    • 均采用相似的损失函数或评估指标

  • 不同点:

    • 医学图像由于其低对比度等特点,需要特别设计网络结构

    • 医学图像的数据获取难度大,数据量小,部分数据处理手段也需定制化

    • 医学图像的评估指标需包含临床相关指标,且损失函数有时需要更精细化的设计

2、简要叙述原始GAN到WGAN-GP的演化历程及原因。简要叙述pix2pix GAN、Cycle GAN、Conditional GAN的原理。

原始GAN到WGAN-GP的演化历程及原因:

  • GAN

    • GAN的目标函数:

      \min_G\max_DV\left( G,D \right) =\mathbb{E}_{x\thicksim p_r\left( x \right)}\left[ \log D\left( x \right) \right] +\mathbb{E}_{z\thicksim p\left( z \right)}\left[ \log \left( 1-D\left( G\left( z \right) \right) \right) \right]
    • 生成器G:

      生成器G的目的是使上述优化函数最小,而在训练生成器G时,控制判别器D权重不变,那么只需要看函数的第二部分,即

      \min_G\mathbb{E}_{z\thicksim p\left( z \right)}\left[ \log \left( 1-D\left( G\left( z \right) \right) \right) \right]

      想让这个函数最小,那么生成器G就希望D(G(z))=1,即生成器G希望判别器D把生成的图片G(z)判别为1,也就是判别为真样本。

    • 判别器D:

      判别器D是用来判别生成器G生成的是假图片(fake)还是真图片(real)。在训练判别器D时,控制生成器G的权重不变,优化函数为

      \max_D\left( \mathbb{E}_{x\thicksim p_r\left( x \right)}\left[ \log D\left( x \right) \right] +\mathbb{E}_{z\thicksim p\left( z \right)}\left[ \log \left( 1-D\left( G\left( z \right) \right) \right) \right] \right)

      log函数的可行域是大于0的,所以必须保证判别器D的输出在0和1之间。那么优化函数最大就是D(x)=1,且D(G(z))=0,即希望判别器D判断真实样本x为真样本,生成出来的虚假样本G(z)为假样本。

  • WGAN

    原始GAN的损失函数训练十分不稳定而且难以收敛,WGAN引入了一种新的分布距离度量方法——Wasserstein距离,也称为(Earth-Mover Distance)简称EM距离,表示从一个分布变换到另一个分布的最小代价。

    • Wasserstein距离:

    W\left( P_r,P_g \right) =\inf_{\gamma \thicksim \prod{\left( P_r,P_g \right)}}\mathbb{E}_{\left( x,y \right) \thicksim \gamma}\left[ \lVert x-y \rVert \right]

    WGAN解决了GAN的主要训练问题。训练WGAN不需要维护在判别器D和生成器G的训练中保持谨慎的平衡,并且也不需要对网络架构进行仔细的设计。WGAN最引人注目的实际好处之一是能够通过训练判别器进行运算来连续地估计EM距离。

  • WGAN-GP

    • WGAN还是有问题以下两个问题:

      • 权重裁剪会导致参数基本都在限制的边界值,极大浪费了模型的参数。

      • 还是很容易梯度消失或者梯度爆炸,需要仔细的调参。

    • WGAN-GP的核心——Gradient Penalty

      Gradient Penalty:判别器相对于输入的梯度的二范数要约束在1附近,这样就能够保证Lipschitz连续。

      L=\underset{Original\ critic\ loss}{\underbrace{\_{\hat{x}\thicksim \mathbb{P}_g}\left[ D\left( \tilde{x} \right) \right] -\_{x\thicksim \mathbb{P}_r}\left[ D\left( x \right) \right] }}+\underset{Their\ gradient\ penalty}{\underbrace{\lambda \_{\hat{x}\thicksim \mathbb{P}_{\hat{x}}}\left[ \left( \lVert \nabla _{\hat{x}}D\left( \hat{x} \right) \rVert _2-1 \right) ^2 \right] }}

      WGAN-GP是对原始WGAN的改进,通过梯度惩罚(Gradient Penalty)替代权重裁剪(Weight Clipping),解决了WGAN训练不稳定、权重裁剪导致梯度消失或爆炸的问题

pix2pix GAN、Cycle GAN、Conditional GAN的原理:

  • pix2pix GAN

    基于成对映射的条件生成对抗网络(CGAN),其核心是通过U-Net结构的生成器将输入图像转换为目标图像,同时采用PatchGAN判别器对图像局部区域进行真伪,训练时结合对抗损失(迫使生成图像逼真)和L1像素损失(确保输入输出严格对齐)。

  • Cycle GAN

    核心是双向循环一致性约束,用于解决无配对数据的域迁移问题。其包含两个生成器和两个判别器,通过循环重建损失强制转换可逆性,同时配合对抗损失确保生成质量。

  • Conditional GAN

    生成器和判别器均以额外条件(如类别标签)为输入,实现可控生成。它在生成器和判别器的输入中引入额外条件信息,生成器基于噪声向量z和条件y生成目标(G(z|y)),判别器则同时评估图像真实性及与条件的匹配度(D(x|y))。

3、NeRF和3D Gaussian Splatting用于医学图像重建的优势和局限性是什么?

NeRF(神经辐射场)

NeRF使用深度学习技术,特别是一种密集的神经网络(通常是多层感知机,MLP),来建模复杂的3D场景。它通过训练一个神经网络来预测给定3D位置和观察方向下的颜色和体积密度。

  • 优势

    • NeRF基于物理模型,通过模拟光线在场景中的传播来创建逼真的图像。它的目标是从多个图像重建出一个全局一致的3D场景,并能从任意新视角进行逼真渲染

    • NeRF的隐式连续表达适合稀疏视角重建(如低剂量CT),生成高保真新视角。

  • 局限性

    • NeRF依赖于神经网络来精确捕捉和渲染复杂的场景细节,需要大量的计算资源,尤其是在训练阶段。

    • NeRF是从多张静态图像进行重建的,故难以建模动态物体(如跳动心脏)。

3D高斯溅射

3D Gaussian Splatting是一种体积渲染技术,经常用于医学影像和科学可视化。它通过将数据点表示为具有高斯权重的样本,然后将这些样本投影到视图平面上,来实现3D数据的可视化。

  • 优势

    • 3DGS更多地关注于科学数据的准确和直观表达,它强调的是数据点的直接表示和属性的清晰显示。

    • 3DGS通常计算上不如NeRF复杂,可以实时进行,适用于交互式数据探索和可视化。

  • 局限性

    • 3DGS在细节上稍逊,所以在重建时容易丢失一些细小结构(如血管)。

    • 3DGS同样需要多视角图像,并且还需要点云来进行重建,但其对动态物体具有一定适应性。

4、简要叙述Diffusion Model在医学图像重建方面的优势和局限性。并给出使用Diffusion Model进行医学图像重建时,进行数据一致性或数据保真约束的关键步骤的伪代码(任意形式均可)

优势

  • 扩散模型通过渐进式去噪过程生成图像,拥有高保真重建能力,生成质量明显优于GAN。

  • 模型训练稳定性强,避免了GAN的模式崩溃问题,且能无缝集成物理成像约束,通过投影操作强制重建结果符合实际测量值。

局限性

  • 由于模型推理需要多步采样迭代,导致计算效率低下。

  • 扩散模型的训练需要大量高质量标注数据,而医学图像标注稀缺且专业性强,较难获取。

数据一致性约束伪代码(以MRI重建为例):

def diffusion_reconstruction(y, M_mask, diffusion_model):
    # y: 欠采样k-space, M_mask: 采样掩码
    x_t = initialize_noisy_image()          # 初始化带噪图像
    for t in range(T, 0, -1):               # 逆向扩散T步
        x_t = denoise_step(x_t, t, diffusion_model)  # 去噪网络预测
        
        # 每k步执行数据一致性约束
        if t % k == 0:  
            k_pred = fft(x_t)               # 图像→k空间
            k_corrected = M_mask * y + (1 - M_mask) * k_pred  # 用实测值替换采样点
            x_t = ifft(k_corrected)         # 返回图像空间
    
    return x_t  # 输出重建图像

数据:

链接: https://pan.baidu.com/s/1_koUdpxwXhgugkH6v7WJqQ?pwd=list

提取码: list

quarter_1mm: 低剂量 CT 图像(噪声图像)

full_1mm:常规剂量CT图像(ground truth)

代码:

链接: https://pan.baidu.com/s/12N6_8HLX5Aq7gEFyVJYkvQ?pwd=list

提取码: list

制作配对数据,可针对任务调整:mk_data_list.py

5、(基础)阅读下面两篇图像后处理的论文(Low-Dose CT with a Residual Encoder-Decoder Convolutional Neural Network / Deep Convolutional Neural Network for Inverse Problems in Imaging),并复现FBPConvNet,结果放PPT汇报。参考代码: train_unet_2d.py。

Low-Dose CT with a Residual Encoder-Decoder Convolutional Neural Network:

f9b1bb553b1f73952fbff7b2ac640472.png

这篇文章提出了一种名为RED-CNN(残差编码器-解码器卷积神经网络)的深度学习方法,用于低剂量CT(LDCT)图像重建。该方法结合了自动编码器反卷积网络跳跃连接,通过对称的卷积层和反卷积层结构(无池化操作),并引入残差补偿机制,直接从图像域处理LDCT重建结果,有效抑制噪声和伪影,同时保留结构细节和低对比度病变。实验在模拟数据和真实临床数据上进行,结果表明RED-CNN在噪声抑制、结构保存和病变检测方面优于主流方法(如TV-POCS、K-SVD、BM3D),并在PSNR、RMSE和SSIM等定量指标及主观评估中表现最优,同时具有较高的计算效率。

Deep Convolutional Neural Network for Inverse Problems in Imaging:

fa663b4d92ad64c784681187ba87b338.png

这篇文章提出了一种名为FBPConvNet的新型深度学习方法,用于解决医学成像中的病态逆问题(如CT稀疏视图重建)。该方法基于关键发现:当正向模型的运算算子(H*H)具有卷积特性时,迭代重建算法可转化为卷积神经网络(CNN)结构。FBPConvNet结合物理驱动的滤波反投影(FBP)与改进的U-Net架构(含多尺度分解、残差学习和通道融合),通过直接反投影后接CNN处理,有效抑制重建伪影并保留结构细节。实验表明,在合成体模、生物医学数据和真实CT扫描中,该方法在50-143个稀疏视图下均优于TV正则化迭代重建,尤其擅长保留复杂纹理,且重建512×512图像仅需不到1秒(GPU),为低剂量成像提供了高效解决方案。

复现FBPConvNet:

  • FBPConvNet核心代码:

    class FBPConvNet(nn.Module):
        def __init__(self):
            super(FBPCONVNet, self).__init__()
            # create network model
            self.block_1_1 = None
            self.block_2_1 = None
            self.block_3_1 = None
            self.block_4_1 = None
            self.block_5 = None
            self.block_4_2 = None
            self.block_3_2 = None
            self.block_2_2 = None
            self.block_1_2 = None
            self.create_model()
    ​
        def forward(self, input):
            block_1_1_output = self.block_1_1(input)
            block_2_1_output = self.block_2_1(block_1_1_output)
            block_3_1_output = self.block_3_1(block_2_1_output)
            block_4_1_output = self.block_4_1(block_3_1_output)
            block_5_output = self.block_5(block_4_1_output)
            result = self.block_4_2(torch.cat((block_4_1_output, block_5_output), dim=1))
            result = self.block_3_2(torch.cat((block_3_1_output, result), dim=1))
            result = self.block_2_2(torch.cat((block_2_1_output, result), dim=1))
            result = self.block_1_2(torch.cat((block_1_1_output, result), dim=1))
            result = result + input
            return result
    ​
        def create_model(self):
            kernel_size = 3
            padding = kernel_size // 2
    ​
            # block_1_1
            block_1_1 = []
            block_1_1.extend(self.add_block_conv(in_channels=1, out_channels=64, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_1_1.extend(self.add_block_conv(in_channels=64, out_channels=64, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_1_1.extend(self.add_block_conv(in_channels=64, out_channels=64, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            self.block_1_1 = nn.Sequential(*block_1_1)
    ​
            # block_2_1
            block_2_1 = [nn.MaxPool2d(kernel_size=2)]
            block_2_1.extend(self.add_block_conv(in_channels=64, out_channels=128, kernel_size=kernel_size, stride=1,
                                                 padding=padding, batchOn=True, ReluOn=True))
            block_2_1.extend(self.add_block_conv(in_channels=128, out_channels=128, kernel_size=kernel_size, stride=1,
                                                 padding=padding, batchOn=True, ReluOn=True))
            self.block_2_1 = nn.Sequential(*block_2_1)
    ​
            # block_3_1
            block_3_1 = [nn.MaxPool2d(kernel_size=2)]
            block_3_1.extend(self.add_block_conv(in_channels=128, out_channels=256, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_3_1.extend(self.add_block_conv(in_channels=256, out_channels=256, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            self.block_3_1 = nn.Sequential(*block_3_1)
    ​
            # block_4_1
            block_4_1 = [nn.MaxPool2d(kernel_size=2)]
            block_4_1.extend(self.add_block_conv(in_channels=256, out_channels=512, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_4_1.extend(self.add_block_conv(in_channels=512, out_channels=512, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            self.block_4_1 = nn.Sequential(*block_4_1)
    ​
            # block_5
            block_5 = [nn.MaxPool2d(kernel_size=2)]
            block_5.extend(self.add_block_conv(in_channels=512, out_channels=1024, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_5.extend(self.add_block_conv(in_channels=1024, out_channels=1024, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_5.extend(self.add_block_conv_transpose(in_channels=1024, out_channels=512, kernel_size=kernel_size, stride=2, padding=padding, output_padding=1, batchOn=True, ReluOn=True))
            self.block_5 = nn.Sequential(*block_5)
    ​
            # block_4_2
            block_4_2 = []
            block_4_2.extend(self.add_block_conv(in_channels=1024, out_channels=512, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_4_2.extend(self.add_block_conv(in_channels=512, out_channels=512, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_4_2.extend(self.add_block_conv_transpose(in_channels=512, out_channels=256, kernel_size=kernel_size, stride=2, padding=padding, output_padding=1, batchOn=True, ReluOn=True))
            self.block_4_2 = nn.Sequential(*block_4_2)
    ​
            # block_3_2
            block_3_2 = []
            block_3_2.extend(self.add_block_conv(in_channels=512, out_channels=256, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_3_2.extend(self.add_block_conv(in_channels=256, out_channels=256, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_3_2.extend(self.add_block_conv_transpose(in_channels=256, out_channels=128, kernel_size=kernel_size, stride=2, padding=padding, output_padding=1, batchOn=True, ReluOn=True))
            self.block_3_2 = nn.Sequential(*block_3_2)
    ​
            # block_2_2
            block_2_2 = []
            block_2_2.extend(self.add_block_conv(in_channels=256, out_channels=128, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_2_2.extend(self.add_block_conv(in_channels=128, out_channels=128, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_2_2.extend(self.add_block_conv_transpose(in_channels=128, out_channels=64, kernel_size=kernel_size, stride=2, padding=padding, output_padding=1, batchOn=True, ReluOn=True))
            self.block_2_2 = nn.Sequential(*block_2_2)
    ​
            # block_1_2
            block_1_2 = []
            block_1_2.extend(self.add_block_conv(in_channels=128, out_channels=64, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_1_2.extend(self.add_block_conv(in_channels=64, out_channels=64, kernel_size=kernel_size, stride=1, padding=padding, batchOn=True, ReluOn=True))
            block_1_2.extend(self.add_block_conv(in_channels=64, out_channels=1, kernel_size=1, stride=1, padding=0, batchOn=False, ReluOn=False))
            self.block_1_2 = nn.Sequential(*block_1_2)
    ​
        @staticmethod
        def add_block_conv(in_channels, out_channels, kernel_size, stride, padding, batchOn, ReluOn):
            seq = []
            # conv layer
            conv = nn.Conv2d(in_channels=in_channels, out_channels=out_channels, kernel_size=kernel_size,
                             stride=stride, padding=padding)
            nn.init.normal_(conv.weight, 0, 0.01)
            nn.init.constant_(conv.bias, 0)
            seq.append(conv)
    ​
            # batch norm layer
            if batchOn:
                batch_norm = nn.BatchNorm2d(num_features=out_channels)
                nn.init.constant_(batch_norm.weight, 1)
                nn.init.constant_(batch_norm.bias, 0)
                seq.append(batch_norm)
    ​
            # relu layer
            if ReluOn:
                seq.append(nn.ReLU())
            return seq
    ​
        @staticmethod
        def add_block_conv_transpose(in_channels, out_channels, kernel_size, stride, padding, output_padding, batchOn, ReluOn):
            seq = []
    ​
            convt = nn.ConvTranspose2d(in_channels=in_channels, out_channels=out_channels, kernel_size=kernel_size, stride=stride, padding=padding, output_padding=output_padding)
            nn.init.normal_(convt.weight, 0, 0.01)
            nn.init.constant_(convt.bias, 0)
            seq.append(convt)
    ​
            if batchOn:
                batch_norm = nn.BatchNorm2d(num_features=out_channels)
                nn.init.constant_(batch_norm.weight, 1)
                nn.init.constant_(batch_norm.bias, 0)
                seq.append(batch_norm)
    ​
            if ReluOn:
                seq.append(nn.ReLU())
            return seq

lr=0.01结果:

  • loss曲线:

434606554bb037040d30c59f93cf564e.png

41400f9e51176d16bdfa14012740e5f8.png

lr=0.0001结果:

  • loss曲线:

a61e5acd6f0de38f224490affac6f106.png

216673d8bc0937ccde3e6fbf73acd58a.png

【分析】

我是在train_unet_2d.py基础上修改实现的FBPConvNet,参考了开源项目。一开始我设置的lr=0.01进行训练,结果发现loss直接过拟合了,后来沿用原代码设计,设置lr=0.0001,发现训练只需20个epoch即可达到收敛,且ssim性能可以达到0.9671。根据tensorboard上的重建图像可以看出,FBPConvNet的设计可以有效抑制重建伪影并保留结构细节。

6、(基础)阅读下面论文(Low-dose CT image denoising using a generative adversarial network with Wasserstein distance and perceptual loss),分别复现基于WGAN-GP的CT图像去噪任务及基于感知损失的CT图像去噪任务,结果放PPT汇报。参考代码分别为:train_unet_2d_gan.py和train_unet_2d_vgg.py。

Low-dose CT image denoising using a generative adversarial network with Wasserstein distance and perceptual loss:

0dd968f8b1a4e44dcecffb7e4e1d1b0d.png

这篇文章提出了一种结合Wasserstein生成对抗网络(WGAN)感知损失(Perceptual Loss)的低剂量CT图像去噪方法(WGAN-VGG)。针对传统MSE损失导致结构细节丢失的问题,该方法创新性地引入Wasserstein距离优化生成器与判别器的对抗训练稳定性,同时利用预训练VGG网络提取高级特征计算感知损失,在抑制噪声的同时有效保留解剖细节。实验基于临床CT数据集验证表明:相较于传统CNN-MSE、CNN-VGG等方法,WGAN-VGG在视觉质量上显著减少“蜡状伪影”,在定量评估中CT值统计特性更接近标准剂量图像,且放射科医生盲评确认其在伪影抑制(3.70分)和综合质量(3.70分)方面最优,为低剂量CT诊断提供了更可靠的图像增强方案。

基于WGAN-GP的CT图像去噪任务:

train_unet_2d_gan.py代码中已给出基于WGAN-GP的CT图像去噪任务关键代码:

  • 判别器损失:WGAN-GP的判别器损失(Wasserstein距离+梯度惩罚)

      d_loss_adv = -torch.mean(real) + torch.mean(fake) + 10 * gradient_penalty
  • 生成器损失:MSE重建损失+对抗损失

      g_loss_adv = -torch.mean(fake)
      g_loss_recon = loss_mse(predictImg, RDCTImg)
      g_loss = g_loss_recon + 0.01 * g_loss_adv
  • 梯度惩罚:

    def compute_gradient_penalty(D, real_samples, fake_samples):
        Tensor = torch.cuda.FloatTensor
        """Calculates the gradient penalty loss for WGAN GP"""
        # Random weight term for interpolation between real and fake samples
        if real_samples.ndim == 5:
            alpha = Tensor(np.random.random((real_samples.size(0), 1, 1, 1, 1)))
        else:
            alpha = Tensor(np.random.random((real_samples.size(0), 1, 1, 1)))
        # Get random interpolation between real and fake samples
        interpolates = (alpha * real_samples + ((1 - alpha) * fake_samples)).requires_grad_(True)
        d_interpolates = D(interpolates)
        fake = Variable(Tensor(real_samples.shape[0], 1).fill_(1.0), requires_grad=False)
        # Get gradient w.r.t. interpolates
        gradients = autograd.grad(
            outputs=d_interpolates,
            inputs=interpolates,
            grad_outputs=fake,
            create_graph=True,
            retain_graph=True,
            only_inputs=True,
        )[0]
        gradients = gradients.view(gradients.size(0), -1)
        gradient_penalty = ((gradients.norm(2, dim=1) - 1) ** 2).mean()
        return gradient_penalty

训练结果:

  • loss曲线:

8a679657f4fac40b963b5966de91f654.png

2fbca5afb789b10a911406e27fb2eeeb.png

39b46814f4566bb2fd5dcfb0a43799fd.png

分析:

基于WGAN-GP的CT图像去噪任务通过引入Wasserstein距离优化生成器与判别器的对抗训练稳定性,从训练曲线以及推理结果上来看,这种方法的训练在500个epoch之后趋于稳定,SSIM最高也可达0.9673,缓解了传统MSE损失导致结构细节丢失的问题。

基于感知损失的CT图像去噪任务:

train_unet_2d_vgg.py中已给出基于感知损失的CT图像去噪任务的关键代码:

  • VGG感知损失:MSE+感知损失的加权和

      loss1 = loss_mse(predictImg, RDCTImg)
      loss3 = vgg_loss_calc(gaussian_smooth(predictImg), gaussian_smooth(RDCTImg), vgg_feature_extractor)
      loss = loss3 * 0.02 + loss1

训练结果:

  • loss曲线:

83c271a9e49ea4f1c00343f7d643cdf6.png

708937c2f156bd0886139cc60bee03cf.png

分析:

基于感知损失的CT图像去噪任务利用预训练VGG网络提取高级特征计算感知损失,在抑制噪声的同时有效保留解剖细节,通过这种方法进行训练,loss曲线呈现一个完美下降并趋于平缓的趋势,从predict img可以看出,该方法确实保留了更多地细节信息。虽然单独使用感知损失最高的SSIM指标才0.9615,相较于基于WGAN-GP的CT图像去噪任务低了0.0058,但根据论文实验结果来看,WGAN-GP与感知损失相结合可以取得更优的效果。

7、(基础)另有代码train_unet_2d_mixup.py使用mixup数据增强技术,同时“model=DenseUNet2d()” 将模型从UNet2d换成了DenseUNet2d;train_unet_2d_mixup_sophia.py在数据增强的基础上用了更强的sophia优化器,跑一下感受一下简单的炼丹技巧,结果放PPT汇报。

mixup数据增强技术:

mixup是一种图像与标签同步混合的数据增强技术,它通过线性插值将两个随机样本的输入数据(如图像像素)和对应标签按比例融合,生成新的训练样本。mixup技术用于对图像进行数据增强,从而弥补训练图像数据集不足,达到对训练数据扩充的目的。需要注意的是,全部训练过程都只采用混合的新图像训练,原始图像不参与训练过程。

sophia优化器:

sophia是一种高效的大规模深度学习优化器,它通过引入轻量化的二阶信息(如损失函数的Hessian对角矩阵估计)替代传统一阶优化器(如Adam)中的动量项,并配合梯度裁剪机制,显著提升了语言模型等复杂任务的训练速度。其核心思想是利用曲率信息自适应调整参数更新步长,在平滑区域增大步长加速收敛,在陡峭区域缩小步长保持稳定,实验证明其收敛速度可达Adam的两倍,尤其在训练大型Transformer语言模型时效果突出。

DenseUNet2d的结果:

  • loss曲线:

45e3820312ae92af4abf473b51caf36f.png

23826cd240cf6dcdc60353d06db41751.png

DenseUNet2d+mixup的结果:

  • loss曲线:

2d01d63a6bfe08bb7762bf73b13f4660.png

deb78a3ae3812e60f15ccbf128681b0f.png

DenseUNet2d+mixup+Sophia的结果:

  • loss曲线:

2271ee1d9824f72203d0297d8a46365c.png

bb824c1925c07610f96767ec8cf79369.png

【分析】

综合来看,三种方法进行训练时ssim的最优值均在0.966附近,使用了mixup数据增强技术后train_mse_loss以及valid_mse_loss下降得更多,二者在原始DenseUNet2d框架下进行训练达到收敛时,分别为405、446;而在使用了mixup技术后,二者的值达到收敛时分别为296、436。使用sophia优化器可以使训练更快地达到收敛,相较于使用AdamW优化器需要大约60个epoch才能达到收敛,使用sophia优化器仅需10个epoch左右即可达到收敛,对训练速度的提升很明显,但是使用sophia优化器时,ssim_loss曲线波动较大,猜测是参数设置的问题。


由于考核时间较为紧张,下述内容尚未来得及开展。

8、(可选)阅读下面两篇双域混合处理的论文(CLEAR: comprehensive learning enabled adversarial reconstruction for subtle structure enhanced low-dose CT imaging Domain progressive 3D residual convolution network to improve low-dose CT imaging),了解双域混合处理算法的优势,然后基于FBPConvNet代码和torch-mando等深度重建工具(也可参考4中DREAM-Net的代码),复现不加GAN的CLEAR算法,以low-dose任务为例,并将FBPConvNet、RED-CNN和CLEAR结果放PPT汇报。

9、(可选)阅读下面两篇深度迭代重建的论文(LEARN: Learned experts' assessment-based reconstruction network for sparse-data CT / DREAM-Net: Deep residual error iterative minimization network for sparse-view CT reconstruction),了解深度迭代重建算法的特点,跑通DREAM-Net, 参考代码https://github.com/DarkBreakerZero/DREAM-Net,基于DREAM-Net代码复现LEARN算法,以sparse-view任务为例,并将FBPConvNet、RED-CNN、CLEAR、DERAM-Net和LEARN 结果放PPT汇报。



动物装饰