本文的撰写参考了TUM – Computer Aided Medical Procedure I

封面来自知识共享许可：https://commons.wikimedia.org/wiki/File:Modern%C3%AD_v%C3%BDpo%C4%8Detn%C3%AD_tomografie_s_p%C5%99%C3%ADmo_digit%C3%A1ln%C3%AD_detekc%C3%AD_rentgenov%C3%A9ho_z%C3%A1%C5%99en%C3%AD.jpg

计算机断层扫描是一种通过X光旋转照射人体而获得人体内部的三维重建影像的技术，我们为了方便，只讨论1D到2D的重建。本章节的知识要点如下：

1. X光

为一种能量为50-120keV的光子。这种射线能够穿透生物组织，但是生物组织会造成X光的能量损失，具体表现在吸收（光电效应）以及散射（康普顿效应）上。

对于这一过程我们需要了解到以下的几个要点：

1. 辐射源本身会发出某种能量的射线 I_0 ；
2. 射线穿透人体的时候，人体本身对射线形成了一种阻碍，我们可以定量描述这种阻碍为某种 mu ；
3. 且其穿透的距离为某种 L ;
4. 最后打击到荧幕上，得到能量为 I 。

则可以得到以下的示意图（人体外界假设真空不损耗能量）：

均匀介质

根据这个简单的模型，我们得到的数量关系为

int_{I_0}^{I}frac{dI}{I}=-mu int_{0}^{L}dx

若推广到进入人体组织内的任意距离，则得到当前距离对应的能量为：

I(x)=I_0e^{-mu x}

这个就是很著名的Beer-Lambert定律。

但目前存在的问题是：我们所希望照射的人体组织并非一个均质实体：骨骼、血液、肌肉等等都并非具有一个同样的，更何况在这些组织中本身也不一定相等（如疏松骨骼和致密骨骼）。这导致我们可能需要一个（在这里我们考虑矢量距离）来描述不同处的衰减系数：

I=I_0text{exp}left(-int_{vec{x} in l}mu(vec{x})dsright)

此式等于构造了一个某一方向的线积分，示意图如下：

非均匀介质，CT图片来自：https://www.flickr.com/photos/pulmonary_pathology/7603361920 （知识共享许可）

对于人体而言，我们一般采用空气和水的衰减系数来将所测量的人体组织衰减系数正则化，得到Hounsfield Units（HU）：

HU=frac{mu-mu_{water}}{mu_{water}-mu_{air}}·1000

以下是HU量表：

图片来自：https://radiopaedia.org/cases/hounsfield-scale-diagram （引用：Fortin, F. Hounsfield scale (diagram). Case study, Radiopaedia.org. (accessed on 25 Jun 2022) https://doi.org/10.53347/rID-77397）

以上便是X光的基本知识。

总结：X光可以通过其透射能力以及人体各种组织间衰减系数的不同，直观地呈现人体的二维内部结构。然而我们期望获得三维结构，以期望对于病理和生理特征进行更好的可视化操作。有两种方法可以进行这种操作，一种是解析法，一种是代数法，下面将逐步介绍这两种方法并且分析其优缺点。

2. 基本的数学工具（1）

重建三维模型需要多张2D的X光照片，这是肯定的。但是如何将这些2D重组为3D模型则是一个问题。在这里简要地介绍几种数学工具，如需了解更多可以自行查阅相关资料，在这里主要给出一个大概的思路（这也将是笔者的一贯作风）。并且本章节暂时不详细介绍这些工具的用途，读者仅需要一个大致印象。

1. Radon Transformation 拉东变换

以二维空间举例：考虑以下不规则二维切片以及类似于X光的透射光束：

拉东变换

其中， f(x,y) 表示该光束路径在切片的未知密度（非均匀介质，线上的每个点密度，或者说衰减系数，都不一样），P是直线到原点的距离。我们最终的目的是希望求解这个 f(x,y) ，但我们可以先看看光束投影到接收器的过程。

因为可以转成参数方程：

xcosalpha+ysinalpha=P

按照之前的Beer-Lambert定律，我们知道这个 f(x,y) 从某种程度上就是衰减系数，因此线上存在一个投影到接收器上公式为：

Radon(P,alpha)=int_{-infin}^{+infin}f(x,y)dl

但是我们需要严谨一点，因为在人体组织外我们认为信号没有损耗，因此定义一个叫做狄拉克函数，其特点是，在定义域上值为1，否则为0，因此构造狄拉克函数以及其定义域为（大概了解含义即可，直接用结论）：

delta(P-x cos alpha-y sin alpha)= begin{cases}0, & P-x cos alpha-y sin alpha neq 0 / 1, & P-x cos alpha-y sin alpha=0end{cases}

因此一般情况下拉东变换改写成（原式积分内部乘以一个狄拉克函数）：

operatorname{Radon}(P, alpha)=int_{-infty}^{+infty} f(x, y) delta(P-x cos alpha-y sin alpha) d l

如果我们希望获得这个切片的全貌，我们可能需要从多个角度进行多次扫描，通过一维构建二维，在这里我们可以把对应的和制作成一个二维坐标系来绘图，叫做Sinogram：

Sinogram

右图为Sinogram图，纵坐标为度数（360度等间距划分为300个间隔），横坐标自然是线到原点的距离。也就是说我们通过300次扫描（绕着切片一周）获得了Sinogram。从左边切片到有右边切片的过程我们叫做前向投影Forward Project。

2. 傅里叶切片定理 Fourier Slice Theorem

我们通过不同角度的投射得到了Sinogram。但是请注意，医疗实践中我们对于上图中的左图的信息是一无所知的，我们仅仅获得了右边的Sinogram（也就是说左边是空白）。那么我们如何通过Sinogram获得左边的图呢？

在此假设大家都明白傅里叶变换及其逆变换（如果不清楚的可以寻找相关资料，此为一信号处理方面十分重要的应用），即将时空域的信号转为频率域信号以及将频率域信号复转为时空域信号。傅里叶切片定理规定了如下几个要素：

1. F_1 和 F_2 分别是1D和2D的傅里叶算子；
2. P_1 是投影算子，将2D投影到1D（例如拉东变换）；
3. S_1 是某切片（在后面解释）。

则有如下的关系：

F_1P_1=S_1F_2

这个式子看起来比较抽象，但字面的解释为：我们知道我们希望求出 f(x,y) ，那么这个 f(x,y) 根据某种投影（如拉东变换）再进行一维傅里叶变换，等价于这个 f(x,y) 的二维傅里叶变换的切片。

在这里解释切片：我们将原函数 f(x,y) 进行二维傅里叶变换后得到频率域坐标系（可以当作将原来的XY坐标系转移到某种UV坐标系），其切片为按照参数方程在该空间上选定的值，即（该 theta 角度可以理解成拉东变换的投影角度，如前文所述）：

S_{theta}=left(F_{2} circ fright)(omega cos theta, omega sin theta)

下图是切片示意图，该图直观地展示了所谓的切片就是傅里叶空间上过原点的某一条直线，可以观察到我们将原来的时空域的 f(x,y) 转换到频率域，每一次的测量都能够对应右图上一个一维slice切片：

引用来自Wikipedia相关词条

下面的关系图更好地阐释了这些变换的关系：

拉东变换和傅里叶切片定理的关系

为了获取 f(x,y) ，基于我们已经知道的Sinogram（对应上图的 P_{theta} ），我们可以得到多角度的切片，根据这些切片重新组合较为完整的（实际上我们不可能也不需要填充整个二维傅里叶空间，只需要有限测量的几个投影即可）二维傅里叶空间，再将其进行逆变换获得到 f(x,y) 。

3. 解析法 – Filtered Backprojection (FBP)

根据上述的数学方法可以执行FBP操作来还原 f(x,y) 。首先进行Backprojection后向投影，有：

f(x, y)=int_{0}^{pi} widehat{P}_{theta}(x cos theta+y sin theta) d theta

theta 即采样角， hat{P}_{theta} 是所获得的投影的值。对于其进行180°的积分即可以恢复到原来的函数。我们在上文中提到：实际上我们不可能也不需要填充整个二维傅里叶空间，只需要有限测量的几个投影即可，那么我们需要在什么样的范围内进行什么样的测量呢？在这里引入Nyquist采样定理，即，我的采样频率应当是信号频率的两倍才可以较好地恢复原始信号：

limit=N_x frac{pi}{2}

其中， N_x 为对应方向上所要观测的像素的数量。如果不能够满足这个观测量则会造成降采样，恢复的效果不会很好：

降采样的后果

对于平行CT光束而言，对目标旋转扫描180°即可，因为：

P_{theta+180°}(t)=P(-t)

当然，目前有多种CT光束方向，如扇形CT，这种CT可能需要超过200°的扫描角，具体情况取决于设备。还有一种技术则是直接投射锥体X光进行造影（即使用2D信号接收器，传统用一维）叫做FDK，可以自行了解，解算时候会考虑到高度要素（Feldkamp, Davis, Kress: Practical cone-beam algorithm, J Opt Soc Am, A6, 612-619, 1984）。

但是单纯进行后向投影是不够的。这里考虑多种因素，如采样的噪声，频率域空间采样点均匀分布在参数方程slice上而产生笛卡尔坐标系上的不均匀，等等，需要进行滤波处理。一般会对Sinogram进行滤波，采用的滤波器有Ram-Lak，Shepp-Logan，Cosine以及Hamming等等。习惯性会把Sinogram进行二维傅里叶转换，然后乘上滤波器。请自行了解。

4. 基本的数学工具（2）

如何求解矩阵的逆？这是比较简单的事情，有的计算器就可以执行这样的操作。那如何执行大矩阵的逆，其秩可能达到几万，数十万甚至数百万？因为按照常规的求解办法，绝大多数计算机都不能承受住这样的解算压力。因此需要新的算法，在这里介绍一种ART（Algebraic Reconstruction Technique）算法。

我们此时拥有矩阵：

WF=P

其中 F 是希望求得的未知数，而 W 则是秩极大的矩阵，因此不能够使用常规的求逆方法。但我们可以分解该矩阵:

begin{aligned} &w_{11} f_{1}+w_{12} f_{2}+ldots+w_{1 N} f_{N}=vec{f} cdot vec{w}_{1}=p_{1}/ &w_{21} f_{1}+w_{22} f_{2}+ldots+w_{2 N} f_{N}=vec{f} cdot overrightarrow{w_{2}}=p_{2}/ …/ &w_{M 1} f_{1}+w_{M 2} f_{2}+ldots+w_{M N} f_{N}=vec{f} cdot overrightarrow{w_{M}}=p_{M} end{aligned}

然后我们可以重复投影 f 到各个超平面 p ，循环往复，直到残差小于某一个阈值，以此来逼近一个较为精确的解：

vec{f}^{(i)}=vec{f}^{(i-1)}+frac{p_{i}-vec{f}^{(i-1)} cdot vec{w}_{i}}{vec{w}_{i} cdot vec{w}_{i}} vec{w}_{i}

一般情况下会遇到这些问题

a. 矩阵的未知数数量多于观测数，导致不存在唯一解，ART的迭代会在某一子空间内震荡；

b. 超平面相对比较平行（即矩阵存在不满秩的倾向，近乎于重复的无用观测），会导致收敛非常缓慢。

解决的方法有很多，如，添加正则项等等，可以尝试阅读吉洪诺夫正则化的概念以获得更深的理解。目前存在各种较好的ART算法的变种，如最常使用的SART法，其将某几个超平面捆绑在一起进行计算，有兴趣可以自行了解。ART本身是一个优化问题，力图找出最小残差的解，因此也可以尝试使用其他优化方法，如高斯牛顿法及其变种等等。

5. 代数法

基于对ART的讲解，那么代数法CT的概念就相对清晰了。首先，X射线对人体组织进行一定程度的采样测量，然后该采样可以被记录在接收器上：

代数法求解CT

多次采样可以生成一个矩阵。在二维情况下，我们可以用一适当格网覆蓋该二维人体组织，那么经过多角度采样，采样射线的参数即为庞大的系数矩阵，而人体组织则为未知数矩阵，接收器所得到的就是观测值。通过矩阵的求解可以更加快捷地恢复人体组织实况，而且这种方法也更加灵活，易于理解。

ART方法存在如下特点：

其残差一般会比较快地收敛，并且随着迭代次数增加，其收敛效率也渐渐降低，因此需要一个特定的阈值来结束迭代：

残差衰减

如图所示，经过10次迭代，所获得的结果就和Ground Truth相差无几了：

Ground Truth对比

6. 解析法和代数法区别辨析

1. 解析法

优点是比较快捷，计算效率很高，最后出图效果很好，并且该方法是比较传统的方法，使用非常广泛。

缺点是只能限制于特定的几何形状以及采样射线模式，并且难以从物理方面优化其效果。

2. 代数法

优点是允许进行复杂的物理建模，且适应于各种集合模式。

缺点是具有很高的计算复杂度。

7. 伪影 Artifacts

伪影是因为CT成像过程中，受限于测量、环境等因素而造成的成像质量减损，如上文提到的降采样就是伪影的一个实例。在这里介绍几个主要伪影：

1. Beam Hardening：低能量的射线因为衰减而接收不到，导致接收器接收到的射线能量总体偏高。一种案例为人体金属植入物导致的：

引用来自American Journal of Neuroradiology July 2014, 35 (7) 1288-1292; DOI: https://doi.org/10.3174/ajnr.A3851 Add to Citavi project by DOI

2. Scattering：人体组织可能会导致X光散射，导致伪影的存在：

引用来自Proceedings of SPIE – The International Society for Optical Engineering DOI:10.1117/12.652283

3. 几何伪影：如人体的运动（呼吸），机器的不稳定，仪器使用时候的校准不精确，等等。

以上问题基本可以使用一些软件、硬件算法，先验知识，以及深度学习进行解决。