基于THz-TDS图像数据集的CT算法仿真验证

Simulation Validation of CT Algorithm Based on THz-TDS Image Dataset

摘要: 本文通过仿真和实验数据验证了一种改进的计算机断层扫描(CT)重建算法。首先，利用MATLAB中的Shepp-Logan头模型对滤波反投影(Filtered Back Projection, FBP)算法进行了仿真。随后，使用Python编写了CT成像算法，并以太赫兹时域光谱(Terahertz Time-Domain Spectroscopy, THz-TDS)数据集为基础进行验证。最终，利用图像处理软件ImageJ对获得的二维断层重建图进行三维重建。

Abstract: This paper validates an improved algorithm for computed tomography (CT) reconstruction through simulation and experimental data. First, the Filtered Back Projection (FBP) algorithm was simulated using the Shepp-Logan head model in MATLAB. Next, a CT imaging algorithm was implemented in Python and validated using a dataset from Terahertz Time-Domain Spectroscopy (THz- TDS). Finally, the two-dimensional reconstructed slices were processed into a three-dimensional reconstruction using ImageJ.

1. 引言

THz CT 成像技术起源于早期的X射线计算层析成像，2002年，B．Ferguson和张希成等人首次将脉冲太赫兹波应用到透射式CT成像中，实现了对鸡骨头的三维重建，最终可以得到鸡骨的大致轮廓，而鸡骨中的精细结构则难以分辨，成像系统的分辨率较低[1]。滤波反投影算法(FBP)原理简单重建速度较快，所以在CT成像中广泛使用，但在太赫兹CT中为了减少样品的扫描时间，使得最终获取的投影数据较少，此时FBP算法重建图像效果会变差，所以迭代算法被应用到其中，法国波尔多大学的B. Recur等人在2011年进行了一项研究，比较了针对缺陷塑料样品的三种重建算法：SART (Simultaneous Algebraic Reconstruction Technique) [2]、OSEM (Ordered Subsets Expectation Maximization) [3]和传统的FBP (Filtered Back Projection) [4]算法。他们定量比较了这三种算法在充足投影角度数据和有限投影数据(少于25个)情况下的几何结构相似度。研究结果显示，在有限投影数据条件下，SART算法和OSEM算法的重建效果均优于传统的FBP算法[5]。2012年，Q. Li等人提出了一种改进的SART算法，被称为MSART算法。该算法在SART算法的基础上融合了高斯低通滤波和数学形态学中的开闭运算，旨在有效降低重建图像的边缘效应，并在较少投影数据下实现高质量图像的重建。然而，该算法的缺点在于无法很好地抑制高频噪声[6]。2014年，B. Recur等人提出了有序子集凸(OSC)算法。该算法将太赫兹波的高斯光束传播模型引入到透射式断层扫描的最大化期望(ML-TR)算法中，从而优化了重建模型。相较于FBP和SART算法，OSC算法能够有效抑制条状伪影和噪声，并且在三维重建过程中更准确地呈现扫描物体的内部细节和外观轮廓[7]。为了产生更高质量的三维图像，以及满足无损检测应用的需求，他们还提出了适用于透射式太赫兹CT的三维图像处理序列，包括分割、标记、渲染等步骤。

然而，太赫兹CT成像技术面临诸多挑战，包括较低的图像分辨率和数据处理上的复杂性。这些挑战催生了对专门针对THz数据特点优化的新型CT重建算法的需求。本文提出一种改进的CT重建算法，通过仿真研究和实际数据验证，我们评估了该算法在处理THz时域光谱(THz-TDS)数据集时的性能，还通过三维图像重建展示了其在实际应用中的潜力。

2. 算法原理及仿真

2.1. 算法原理

在太赫兹CT成像算法中，目前最常用的算法就是滤波反投影算法(FBP)，FBP算法首先对数据作滤波处理，之后进行反投影即可重建图像。滤波反投影算法以Radon变换[8]和傅里叶中心切片定理[9]为基本原理。Randon变换描述了光线照射样品时在该方向上形成投影数据的过程，在二维平面有一个被成像物体，当一束射线照射在物体上，透射穿透物体，被另一侧的探测器接收。物体对射线的影响体现在物体的空间分布上，定义为$f\left(x,y\right)$ ，坐标为二维空间坐标[10]。所以射线沿传播路径被衰减后的投影值可以定义为

$\text{P}\left(\theta ,\text{t}\right){\displaystyle {\int }_{L\left(\theta ,t\right)}f\left(x,y\right)\text{d}l}=\Re \left(f\left(x,y\right)\right)$ (1)

直线L为射线的路径，θ是射线的角度，t是射线距离旋转轴的水平垂直偏移量，该式子就是被成像物体的二维空间投影函数。也被称为Radon变换。对上面式子进行傅里叶变化可以得到：

$P\left(\theta ,\omega \right)={\displaystyle {\int }_{-\infty }^{+\infty }P\left(\theta ,t\right){\text{e}}^{-j\omega t}\text{d}t}$

$={\displaystyle {\int }_{-\infty }^{+\infty }{d}_t}{\displaystyle \int {\displaystyle {\int }_{-\infty }^{+\infty }f\left(x,y\right)\delta \left(x\sin \theta +y\cos \theta -t\right){\text{e}}^{-j\omega t}{\text{d}}_x{\text{d}}_y}}$

$={\displaystyle \int {\displaystyle {\int }_{-\infty }^{+\infty }f\left(x,y\right){\text{e}}^{-j\omega \left(x\sin \theta -y\cos \theta \right)}\text{d}x\text{d}y}}$ (2)

$={\displaystyle \int {\displaystyle {\int }_{-\infty }^{+\infty }{\text{e}}^{-j\left(\omega \sin {\theta }_x-\omega \cos {\theta }_y\right)}f\left(x,y\right)\text{d}x\text{d}y}}$

由上面式子可以看出，Radon变换的傅里叶变换是$f\left(x,y\right)$ 二维空间傅里叶变换。说明一维投影值的傅里叶变换等同于物体空间函数的二维傅里叶变换在探测器移动方向上过旋转轴中心的一个切片。如图1所示，傅里叶中心切片定理指出：二维函数$f\left(x,y\right)$ 在某个方向上投影$\text{P}\left(\theta ,\text{t}\right)$ 的一维傅里叶变换$P\left(\theta ,\omega \right)$ 等于该函数的二维傅里叶变换$\text{F}\left({\omega }_x,{\omega }_y\right)$ 沿垂直方向过中心点的片段。在一个方向上的反投影就是二维傅里叶变换函数$\text{F}\left({\omega }_x,{\omega }_y\right)$ 中的一条中心片段，所以当探测器与辐射源围绕被测物体旋转，得到180˚方向范围的投影数据，就能计算出完整的频域图像$\text{F}\left({\omega }_x,{\omega }_y\right)$ ，再进行傅里叶反变换就能恢复出原始切片图像。此方法称为直接傅里叶变换法，但是这种重建方法存在一定的缺陷，投影数据的采集来自各个角度不同方向。沿投影方向上采样点的密度会随着远离旋转中心而变低，在图像上会呈现出图像的低频分量较多。过多的低频分量则会导致图像发生失真变得模糊不清。滤波反投影算法就是在获得投影数据的基础之上，添加权重因子，对投影数据进行适当的滤波处理，保留尽可能多的高频分量，滤掉部分低频分量，从而使重建图像更加清晰。

Figure 1. Fourier central slice theorem

图1. 傅里叶中心切片定理

2.2. 算法仿真

在CT重建算法的研究中为了客观的评价各种算法的有效性，通常选用公认的模型作为研究对象。以经典的Shepp-Logan头模型为仿真物体，如图2所示。其是由Sheep和Logan提出的，它由10个位置、大小、方向密度各异的椭圆叠加而成，最终头模型图像的灰度值是由几个相互重叠的椭圆密度代数累加得到[11]。

Figure 2. Simulation object Shepp-Logan head model

图2. 仿真物体Shepp-Logan 头模型

用MATLAB中的Shepp-Logan头模型图对CT成像算法进行仿真，校验了直接反投影法、滤波反投影算法以及两种滤波器的成像效果。仿真图的大小设为256*256个像素点，然后对其进行Radon变换，得到的正弦图如图3所示。

Figure 3. Sinusoidal image of the head model

图3. 头模型正弦图

最后采用直接反投影法和滤波反投影算法进行图像重建，如图4所示，直接反投影法得到的重建图比较模糊，而分别用了RL滤波器和SL滤波器的重建图就比较清晰，很好的还原了Shepp-Logan头模型图。因此，通过MATLAB对CT算法进行仿真，证实了太赫兹CT成像系统实验的可行性。

Figure 4. Reconstructed image of the head model

图4. 头模型重建图

3. 算法验证及三维重建

传统的CT重建算法是采用直接反投影的方法来实现样品的二维图像重建，重建流程框图如图5所示，但通过这种直接反投影得到的重建图会出现模糊和失真，使得重建效果较差。于是可以采用Ram-Lak滤波器、Shepp-Logan滤波器等滤波器对数据进行滤波处理，保留尽可能多的高频分量，滤掉部分低频分量，从而使重建图像更加清晰。

Figure 5. Algorithm flowchart

图5. 算法流程框图

在Python环境中对算法进行优化以适应THz-TDS数据的特性，由于太赫兹波穿过物体时不是像X射线那样沿直线传播，而是会出现散射、衍射和折射效应，使得重建图严重失真，本文将Ram-Lak滤波器和Shepp-Logan滤波器应用到处理THz-TDS数据的CT算法中，Ram-Lak滤波器，又称为坡度滤波器，是在计算机断层扫描(CT)图像重建中广泛使用的一种滤波器。它在频域中执行高通滤波，用来补偿图像重建过程中高频信息的损失。这种滤波器主要强调图像中的边缘和细节，从而提高重建图像的清晰度和锐度。在太赫兹CT成像中，投影数据通常经过Radon变换，得到在频域上的投影数据。为了进行反投影重建，需要对这些投影数据进行滤波处理，以增强图像质量。Shepp-Logan滤波器是一种在计算机断层扫描(CT)图像重建中使用的高通滤波器。它是由Larry Shepp和Benjamin F. Logan提出的，用于改善图像的边缘清晰度。与Ram-Lak滤波器相比，Shepp-Logan滤波器提供了更平滑的滤波效果，因此在处理时降低了噪声的影响，但同时可能会牺牲一些细节的锐度。这种滤波器在重建过程中能够更好地平衡图像细节与噪声之间的关系，适合于需要高图像质量和较低噪声水平的应用场景。

采用Weng-Tai Su等人在研究中使用的THz-TDS数据集[12]对Python实现的算法进行验证。该数据集包含七个用聚苯乙烯材料3D打印的模型在不同THz频率下的透射数据，主要有时域信号峰值和不同频域下的振幅和相位，适合用于验证CT重建算法的效果。这七个模型分别为装在不透明纸筒里的盒子、头骨、伊布、DNA、机器人、北极熊、鹿[12]。如图6所示。

Figure 6. Seven 3D printed models

图6. 七个3D打印的模型[13]

在本研究中，我们采用太赫兹时域光谱(THz-TDS)数据集中的时域信号峰值作为CT成像的投影数据。这一方法是从不同角度收集的信号峰值来构造图像的二维截面，每个模型共计获取了60个投影面，如图7所示，为某一角度下的一个投影面。

Figure 7. Frontal projection of the model

图7. 模型的正面投影图

将不同角度下获取的不同位置投影数据提取出来生成正弦图，如图8所示，以便于进一步分析和图像重建。正弦图是通过将每个角度的投影数据沿着垂直于投影方向的轴进行排列而生成的。此图形显示了数据中的傅立叶空间组成，其中每行代表一个角度的投影，而列则对应于投影数据中的位置。正弦图是计算机断层扫描(CT)中一个重要的中间产物，因为它将一维投影转换为可以应用二维逆变换的形式。

Figure 8. Sinusoidal image of the model

图8. 模型的正弦图

随后，我们应用了I-radon变换来处理正弦图，从而实现对模型的二维重建。I-radon变换是Radon变换的逆过程，用于从投影数据重建原始图像。在本研究中，该变换通过集成不同角度的投影来恢复模型的二维断层重建图，从而观察并分析被扫描物体的内部结构。如图9所示，二维断层重建图对不透明纸筒里的七个模型进行了还原，但二维重建图出现了失真和模糊，模型越复杂重建结果越差。

Figure 9. Two-dimensional reconstruction image of the model

图9. 模型的二维重建图

完成二维重建后，我们进一步使用ImageJ软件[14]来执行三维重建。这一过程涉及将连续的二维断层重建图叠加并融合，形成完整的三维模型。ImageJ软件提供了一系列工具和插件，能够将二维图像序列转换成三维表示，使得研究者能够从不同角度和深度观察模型，更全面地理解其结构和功能。如图10所示，重建出的三维图能够反映出模型的大致轮廓及位置，结构越简单还原度越高，比如盒子、DNA和头骨模型可以重建出大致轮廓，其他复杂模型重建比较困难，说明CT算法还有待优化改进但能实现该有的功能。

Figure 10. Three-dimensional reconstruction image of the model

图10. 模型的三维重建图

通过以上方法，本研究不仅展示了利用THz-TDS数据集进行CT图像重建的有效性，也强调了高级成像软件在现代医学成像技术中的应用价值。此外，本工作的成功实施展示了跨学科技术在解决复杂科学问题中的重要性，为未来利用太赫兹成像数据进行更精确的生物医学研究奠定了基础。

最后通过自己测的实验数据测试算法的可行性，通过将一个六角扳手插入泡沫中，测试隐藏在泡沫内部的金属扳手，如图11所示

Figure 11. Physical image of a metal wrench in foam

图11. 泡沫中金属扳手实物图

样品中泡沫边长40 mm，高度为60 mm，金属扳手直径5 mm。设置X轴的扫描精度为1 mm，对样品进行X方向扫描50个数据点，每7.2度进行一次X轴方向的扫描，总共扫描一周，旋转50次，获得50个投影面。通过投影数据得到的正弦图如图12所示。

Figure 12. Sinusoidal image of a metal wrench in foam

图12. 泡沫中金属扳手正弦图

从正弦图中可以清晰地看到中心有一个暗色的阴影，由此可判断为泡沫中的金属扳手，周围较蓝色部分为泡沫，用FBP算法得到的重建图如图13所示。

Figure 13. Reconstructed image of a metal wrench in foam

图13. 泡沫中金属扳手重建图

从重建图中可以清晰地看到金属扳手的断层图，验证了该算法在应用到实际应用中的可行性。由于太赫兹波穿过样品时，不像X射线一样沿直线传播，而是发生了衍射、散射、折射等效应导致最终重建图出现一定的模糊。

4. 结论

本文使用THz-TDS数据集进行CT算法的验证。首先，利用MATLAB中的Shepp-Logan头模型图进行仿真，然后在Python环境中优化了算法，将Ram-Lak滤波器和Shepp-Logan滤波器应用到THz-TDS数据的CT算法中。通过收集60个投影面的数据，生成了正弦图，并利用I-radon变换完成了模型的二维重建。虽然在处理复杂模型时出现了一些失真和模糊，但仍然能够辨认出基本结构。进一步使用ImageJ软件进行三维重建，显示了三维成像技术的潜力与挑战。结果验证了THz-TDS数据集在医学成像领域的应用价值，并强调了使用高级成像软件的重要性。未来的工作将集中在优化算法处理更复杂模型的能力上，并探索新的算法或技术策略，以提高重建的精度和减少失真。

NOTES

^*通讯作者。