CT 图像重建¶

CT 重建（Computed Tomography Reconstruction）的过程，简单来说就是把一堆从不同角度拍摄的“影子”（投影数据），还原成一个真实物体的“三维模型”（断层图像）。

目前最常用的标准算法是 滤波后向投影 (Filtered Back Projection, FBP)。

数据采集：获取正弦图 (Sinogram)¶

CT 机架旋转，X 射线穿过人体。探测器在每一个角度记录下穿透后的射线强度。人体组织越硬（如骨骼），射线衰减越多，影子越深。所有角度的投影数据堆叠在一起，形成一张像波浪一样的图，称为 正弦图 (Sinogram)。

每个横行（平行于x轴的直线）：切面投影深浅¶

含义：代表 X 射线穿过人体后，落在探测器排上的空间位置。
物理对应：横轴上的每一个点对应探测器阵列中的一个特定像素（或对应平行束扫描中的位移距离 $s$）。

横轴：扫描角度 ($\theta$)¶

含义：代表 CT 机架转动的角度偏移。
物理对应：随着时间推移，机架从 $0^\circ$ 旋转到 $180^\circ$（或 $360^\circ$），每一行数据代表在特定角度下获取的一维投影（Projection）。

傅里叶切片定理¶

想象我们要重建一个二维物体 $f(x, y)$：

空间域： 探测器在某个角度 $\theta$ 收集到的投影数据，本质上是物体在该方向上的线积分。
频率域： 定理证明，这个角度 $\theta$ 下的投影数据的一维傅里叶变换，正好等于原物体二维傅里叶变换在角度 $\theta$ 穿过原点的一条切片（直线）。

当你旋转 C 形臂（或 CT 机架）采集 360 度的投影时，你实际上是在物体的二维频域平面上，不停地旋转这条“切片直线”：

原点附近（低频）： 所有的切片直线都会经过原点。这意味着频率域的中心被无数次重复采样，采样点极其密集。
远离原点（高频）： 随着半径 $r$（即频率 $f$）增大，切片直线之间的间距越来越大，采样点变得稀疏。

这种旋转采集方式，天然地导致了频域数据在极坐标下的分布是不均匀的——低频过剩，高频不足。

预处理：对数转换¶

探测器接收到的是光子强度，但重建需要的是“衰减系数”。根据比尔-朗伯定律（Beer-Lambert Law），对原始数据取对数，将指数级的衰减转换成线性的厚度累加，方便后续数学计算。

当一束单色光照射在物质上时，光子会被吸收或散射。随着穿透深度的增加，光强度呈指数级下降。其数学表达式为：

\[I = I_0 \cdot e^{-\mu x}\]

$I_0$：入射光强度（探测器发射出的原始强度）。
$I$：透射光强度（穿过物体后被探测器接收到的强度）。
$\mu$：线性衰减系数（Linear Attenuation Coefficient），代表物质对光的阻挡能力。在CT中，不同组织（骨骼、肌肉、脂肪）的 $\mu$ 值不同。
$x$：光线穿过物质的路径长度（厚度）。

在原始测量数据中，我们得到的是 $I$。但对于CT重建算法（如FBP或迭代重建）来说，我们需要的是物质特性的线性累加。通过数学变换：

移项：$\frac{I}{I_0} = e^{-\mu x}$
取自然对数：$\ln(\frac{I}{I_0}) = -\mu x$
整理： $$\mu x = \ln(\frac{I_0}{I}) = \ln(I_0) - \ln(I)$$

滤波 (Filtering)¶

使用 Ramp 滤波器 与每一行投影数据进行卷积。这相当于给每个投影边缘加了负值。它能抵消掉物体中心由于过度堆积而产生的模糊感。

在频率域中，Ramp 滤波器是一个理想的高通滤波器。它的频率响应 $H(f)$ 与频率 $f$ 成正比：

\[H(f) = |f|\]

在极坐标系下的积分采样会导致低频信号过度叠加。乘以 $|f|$ 恰好补偿了这种采样密度的不均匀，增强了高频边缘，削弱了低频背景。

在极坐标系 $(r, \theta)$ 下进行面积积分（还原回图像）时，面积元是 $r \, dr \, d\theta$。

\[f(x,y) = \int_0^\pi \int_{-\infty}^{\infty} P_\theta(f) \cdot |f| \cdot e^{j2\pi f(x\cos\theta + y\sin\theta)} df d\theta\]

如果不乘 $|f|$，直接把所有角度的投影叠加回去（简单反投影），中心重叠的低频分量会因为累加过多而产生巨大的能量，导致图像中心像被揉开了一样模糊。
这个 $|f|$（即 Ramp 滤波器）在数学上扮演了雅可比行列式（Jacobian）的角色，用来抵消极坐标采样带来的权重偏差：
当频率 $f$ 趋近于 0 时，权重也趋近于 0（压制过度叠加的低频）。
当频率 $f$ 增大时，权重线性增加（补偿稀疏的高频）。

2. 离散实现：空间域卷积核¶

在计算机处理中，我们通常不在连续频域操作，而是使用离散的卷积核。最著名的实现是 Ram-Lak 滤波器。

如果在空间域（距离域）表示，Ramp 滤波器的离散冲击响应 $h(n)$ 定义如下：

\[h(n) = \begin{cases} \frac{1}{4d^2}, & n=0 \\ -\frac{1}{(n\pi d)^2}, & n \text{ 为奇数} \\ 0, & n \text{ 为偶数} \end{cases}\]

其中 $d$ 是检测器单元的间距。

3. 具体代码实现步骤¶

实现 Ramp 滤波通常有两种路径：频域乘法（更常用，速度快）和时域卷积。

A. 频域实现（基于 FFT）¶

补零（Padding）： 为了防止循环卷积带来的伪影，先将投影数据（Sinogram 的每一行）长度扩展到原来的 2 倍。
生成滤波器： 创建一个从 $-f_{max}$ 到 $+f_{max}$ 的线性斜坡数组，取绝对值。
快速傅里叶变换： 将投影数据转换到频域。
相乘： 将频域信号与线性斜坡进行点乘。
逆傅里叶变换： 转回空间域，取实部，截取回原长度。

B. Python/PyTorch 伪代码示例¶

既然你在学习 PyTorch，可以参考这个逻辑：

import torch

def ramp_filter(projection_data):
    # projection_data shape: [views, detectors]
    n = projection_data.size(-1)
    n_pad = 2 ** (n - 1).bit_length() * 2 # 补零到2的幂次提高FFT效率

    # 1. 在频率域创建斜坡
    freqs = torch.fft.fftfreq(n_pad)
    h = torch.abs(freqs) # 理想斜坡 |f|

    # 2. 对投影数据做FFT
    projection_f = torch.fft.fft(projection_data, n=n_pad)

    # 3. 频域滤波
    filtered_f = projection_f * h.to(projection_data.device)

    # 4. IFFT转回空间域
    filtered_data = torch.fft.ifft(filtered_f).real
    return filtered_data[:, :n] # 截取有效部分

4. 实际应用中的改进¶

理想的 Ramp 滤波器虽然边缘增强效果好，但会急剧放大噪声（因为噪声通常存在于高频）。为了平衡画质，医学成像中常使用变形版本：

Shepp-Logan 滤波器： 使用 $sinc$ 函数对斜坡进行平滑，抑制极高频噪声。
Cosine / Hamming Window： 在斜坡后面加一个窗函数，让高频部分平滑衰减。

简单来说： 实现 Ramp 滤波就是把每一行投影数据丢进 FFT，乘以一个“V字型”的频率响应函数，再变回来。

你目前在做 BME 的转专业准备，如果涉及到医学影像处理（如 CT 原理），理解这个“高频补偿”逻辑对面试非常有帮助。需要我帮你推导一下为什么它能消除 FBP 中的 1/r 模糊吗？

后向投影 (Back-Projection)¶

将滤波后的每一个角度的影子，按照它射入时的原路径，涂回到重建矩阵中。当成百上千个角度的影子在中心重叠时，无关的干扰会互相抵消，而物体的真实形状会因为多次叠加而显现出来。