简介：这份资源面向医学影像重建方向的学习者与研究人员，聚焦三维锥束CT（CBCT）的投影与反投影重建，核心是FDK算法的Matlab实现。FDK算法主要分为三步：先对投影数据加权，再对加权结果滤波，最后完成反投影并输出重建数据，同时包内还涉及迭代重建等对比示例，适合具备一定Matlab基础、希望理解CT成像原理与算法流程的读者。资源包共17个文件，约200KB，以12个m脚本为主体，另含xml配置、txt说明、png截图与mat数据文件，脚本覆盖投影、滤波、反投影及多种重建方法，数据文件可直接用于实验验证。目前已有1098人学习下载。通过阅读与运行这些脚本，读者能够掌握FDK重建的完整链路，理解加权、滤波、反投影各环节的作用，并借助迭代重建示例对比不同算法的效果，为后续研究或工程实现提供可复用的代码基础。
1. FDK 重建：为什么锥束 CT 的"反投影"不是简单地把像素投回去
锥束 CT 的 FDK 算法，名字听起来像三个人的缩写拼在一起，实际上它解决的是一个很具体的问题：当 X 射线源和探测器绕物体旋转一圈，采集到的是一堆二维投影图，怎么把它们还原成三维体数据。直觉上你会觉得，反投影嘛，就是把每个角度的投影值沿原路径"倒着涂"回三维空间，叠加起来就完事了。但真这么做，图像会糊得没法看——中心区域过亮、边缘发暗、密度不均匀，重建出来的东西根本没法用于诊断或测量。
FDK 的核心贡献，是在反投影之前加了一步"滤波"，并且对锥束几何做了近似处理。它把三维重建拆解成：先对每个投影图做一维滤波，再按锥束几何做加权反投影。这个"先滤波再反投影"的顺序，以及滤波器的设计，才是 FDK 能 work 的关键。搞懂这个机制，你才能理解为什么代码里要先对投影做 FFT、乘上那个奇怪的斜坡函数、再逆变换，而不是直接累加。这份 Matlab 实现的资源，把 FDK 从公式到代码的每一步都摊开了，适合想从底层搞明白锥束重建到底在算什么的人。
2. 从投影到体素：FDK 的三步拆解与几何加权
2.1 滤波反投影的基本框架
FDK 算法可以拆成三个连续步骤：投影加权、一维滤波、锥束反投影。第一步的加权因子是 ( \frac{D}{\sqrt{D^2 + u^2 + v^2}} )，其中 ( D ) 是射线源到旋转中心的距离，( u ) 和 ( v ) 是探测器上的坐标。这个加权的作用是补偿锥束射线路径长度不同带来的强度差异——靠近探测器边缘的射线穿过物体的路径更长，衰减更大，不加权的话重建出来边缘会偏暗。
第二步是滤波。对每个投影图，沿水平方向（探测器 u 方向）做一维傅里叶变换，乘以斜坡滤波器 ( |\omega| ) 加窗函数，再逆变换。斜坡滤波器的作用是锐化投影，补偿反投影过程中固有的模糊。如果不滤波直接反投影，相当于对每个体素做了低通平滑，细节全丢。
第三步是反投影。对于三维空间中的每个体素点 ( (x, y, z) )，根据当前旋转角度 ( \theta )，计算它在探测器上的投影位置 ( (u, v) )，然后取滤波后的投影值累加。这里的关键是：锥束几何下，不同高度的体素对应的射线倾斜角度不同，所以反投影时必须按实际射线路径计算 ( u ) 和 ( v )，不能简单当成平行束处理。
% 以某个旋转角度 theta 为例，对单个投影做滤波 proj = proj_raw .* weight; % 第一步：锥束加权 proj_fft = fft(proj, [], 2); % 沿探测器水平方向做 FFT filter = abs(linspace(-1, 1, size(proj,2))); % 斜坡滤波器 filter = filter .* hann(size(proj,2))'; % 加 Hann 窗抑制高频噪声 proj_filtered = real(ifft(proj_fft .* filter, [], 2)); % 逆变换回空域这段代码里，weight就是前面说的几何加权因子，filter是斜坡滤波器加窗。注意fft和ifft都是沿第二维（探测器水平方向）做的，因为滤波只在投影图的行方向进行，列方向（垂直方向）不滤波。加 Hann 窗是为了压制高频噪声，代价是稍微损失一点空间分辨率——这是 FDK 实现里常见的折中。
2.2 锥束反投影的坐标映射
反投影这一步，最容易翻车的地方是坐标计算。对于旋转角度 ( \theta )，体素 ( (x, y, z) ) 在探测器上的投影位置由以下公式给出：
[ u = \frac{D \cdot (x \cos\theta + y \sin\theta)}{D + x \sin\theta - y \cos\theta} ] [ v = \frac{D \cdot z}{D + x \sin\theta - y \cos\theta} ]
分母里的 ( D + x \sin\theta - y \cos\theta ) 是射线源到体素沿中心射线方向的距离。这个项不能忽略，否则锥束角度大的时候重建会严重变形。很多初学者写的反投影代码只用了平行束的公式，结果重建出来的物体在 z 方向被拉伸或压缩，这就是没做锥束校正的后果。
% 对每个体素 (x, y, z) 和每个角度 theta 做反投影 for iz = 1:nz for iy = 1:ny for ix = 1:nx % 当前体素坐标 x = (ix - nx/2) * dx; y = (iy - ny/2) * dy; z = (iz - nz/2) * dz; % 累加所有角度的贡献 for itheta = 1:nangles theta = angles(itheta); denom = D + x*sin(theta) - y*cos(theta); u = D * (x*cos(theta) + y*sin(theta)) / denom; v = D * z / denom; % 双线性插值取滤波后的投影值 val = bilinear_interp(proj_filtered(:,:,itheta), u, v); volume(ix, iy, iz) = volume(ix, iy, iz) + val; end end end end这段代码是 FDK 反投影的核心。denom就是前面说的距离项，u和v是探测器坐标，bilinear_interp做双线性插值因为计算出的坐标通常不是整数。最后累加时要乘以角度增量 ( \Delta\theta ) 和适当的缩放因子，具体取决于投影数据的归一化方式。实际实现中，这个四重循环会非常慢，通常会用向量化或 GPU 加速，但原理就是这些。
2.3 滤波器的选择与窗函数的影响
斜坡滤波器 ( |\omega| ) 是理论上的最优选择，但它对高频噪声的放大也很厉害。实际数据总有噪声，所以必须加窗。常见的窗函数有 Hann、Hamming、Shepp-Logan、Cosine 等。不同窗函数对重建质量的影响主要体现在两个方面：空间分辨率和噪声水平。
Hann 窗的表达式是 ( 0.5 + 0.5 \cos(\pi \omega / \omega_{\max}) )，它在高频端衰减到零，噪声抑制好，但会损失一些细节。Shepp-Logan 窗是 ( \frac{\sin(\pi \omega / \omega_{\max})}{\pi \omega / \omega_{\max}} )，在高频端衰减慢一些，分辨率保留更好，但噪声会大一点。选择哪个窗，取决于你的数据质量和应用需求——低剂量 CT 通常用 Hann，高分辨率工业 CT 可能用 Shepp-Logan 或更激进的窗。
% 不同窗函数的对比 omega = linspace(-1, 1, N); ramp = abs(omega); hann_win = 0.5 + 0.5 * cos(pi * omega); shepp_win = sin(pi * omega) ./ (pi * omega); shepp_win(isnan(shepp_win)) = 1; % 处理 omega=0 处的除零 filter_hann = ramp .* hann_win; filter_shepp = ramp .* shepp_win;这段代码展示了两种窗函数的构造方式。注意shepp_win在 ( \omega = 0 ) 处需要特殊处理，因为分母为零。实际使用时，窗函数通常还要做归一化，保证直流分量不被改变。滤波器的设计是 FDK 实现里最需要调参的部分，没有之一。
3. 避坑指南：FDK 实现中最容易翻车的五个地方
3.1 投影数据归一化不一致导致重建值偏差
现象：重建出来的体素值整体偏大或偏小，和预期密度对不上。
原因：FDK 的滤波和反投影步骤对投影数据的尺度很敏感。如果投影数据没有正确归一化（比如没除以射线源到探测器的距离平方，或者没乘以角度增量），最终体素值会差一个常数因子。
解决：在滤波之前，确保投影数据已经按照 ( \frac{D}{\sqrt{D^2 + u^2 + v^2}} ) 加权，并且在反投影累加时乘以 ( \Delta\theta )。如果用的是扇形束转锥束的近似，还要注意扇形角到锥角的转换系数。
3.2 探测器坐标原点偏移导致重建偏移
现象：重建出来的物体不在图像中心，或者出现明显的重影。
原因：探测器坐标的原点没有对准射线源到旋转中心的中心射线。实际系统中，探测器可能有横向偏移，或者旋转轴和探测器中心不重合。
解决：在计算 ( u ) 和 ( v ) 之前，先减去探测器的中心偏移量。如果不知道偏移量，可以通过重建一个已知物体（比如小球）来标定，看重建结果是否居中。
3.3 角度采样不足导致条状伪影
现象：重建图像上有明显的放射状条纹，尤其在物体边缘。
原因：FDK 要求角度采样满足 Nyquist 条件，即每个角度间隔对应的探测器像素位移不超过一个像素。如果角度采样太稀疏，高频信息会混叠，产生条状伪影。
解决：增加投影角度数，或者使用更小的角度间隔。如果数据已经采集了，可以尝试在滤波阶段加更强的窗函数，牺牲分辨率来换伪影抑制。但根本办法还是采集时保证足够的角度采样。
3.4 锥束角度过大导致 FDK 近似失效
现象：重建物体在 z 方向出现明显的模糊或变形，远离中心平面的切片质量差。
原因：FDK 本质上是把锥束近似成一系列倾斜的扇形束，当锥角很大时（比如超过 ±10°），这个近似误差会显著增大，导致重建质量下降。
解决：如果锥角确实很大，FDK 可能不是最佳选择，需要考虑迭代重建或精确的锥束重建算法。如果只能用 FDK，可以尝试对 z 方向做额外的滤波或插值来补偿，但效果有限。
3.5 内存溢出与计算速度的平衡
现象：重建一个 512×512×512 的体数据时，Matlab 报内存不足，或者跑了一晚上还没出结果。
原因：FDK 反投影是四重循环，计算量是 ( O(N^3 \times N_{\text{angles}}) )，而且需要同时存储投影数据和体数据。如果直接用三重循环加角度循环，内存和速度都会崩。
解决：把反投影向量化，一次处理一个切片或一个角度。或者用parfor并行化角度循环。如果内存还是不够，可以把体数据分块重建，最后拼接。Matlab 的gpuArray也能加速，但需要支持 CUDA 的显卡。
4. 设计取舍：FDK 为什么能在工业界活这么久
FDK 算法从 1984 年提出到现在，已经四十年了，中间出了无数种改进和替代方案，但它在工业 CT 和医学 CT 里依然是最常用的重建算法之一。原因不是它最精确，而是它在精度、速度和实现复杂度之间找到了一个极好的平衡点。
迭代重建算法（比如 SART、OSEM）理论上能处理更复杂的几何、更少的投影、更强的噪声，但计算量是 FDK 的几十倍甚至上百倍，而且需要调参、容易陷入局部最优。FDK 是一步到位的解析解，没有迭代，没有收敛问题，只要投影数据质量还行，重建结果就稳定可靠。对于产线上的工业 CT，节拍是以秒计的，FDK 的速度优势是决定性的。
另一个关键取舍是滤波器的设计。理论上最优的斜坡滤波器会放大噪声，所以必须加窗。窗函数的选择本质上是在分辨率和噪声之间做权衡。没有哪个窗是绝对最优的，只能根据具体应用调。这份 Matlab 资源把滤波器的构造过程完整写出来了，你可以自己改窗函数、看重建效果，这比直接调库函数理解深得多。
还有一个容易被忽略的点：FDK 对几何参数非常敏感。射线源到旋转中心的距离、探测器到旋转中心的距离、探测器像素尺寸、旋转轴偏移，这些参数差一点，重建质量就差很多。实际系统里，这些参数需要仔细标定。代码里如果把这些参数写死了，换一台设备就得改。好的实现应该把这些参数做成可配置的，并且有标定流程。
% 几何参数配置示例 geo.D = 1000; % 射线源到旋转中心距离 (mm) geo.DSD = 1500; % 射线源到探测器距离 (mm) geo.pixel_size = 0.2; % 探测器像素尺寸 (mm) geo.nx = 512; % 重建体数据 x 方向尺寸 geo.ny = 512; % 重建体数据 y 方向尺寸 geo.nz = 512; % 重建体数据 z 方向尺寸 geo.angles = linspace(0, 2*pi, 360); % 旋转角度这段配置代码看起来简单，但每个参数都影响重建结果。D和DSD的比值决定了锥束角度，pixel_size影响空间分辨率，angles的采样密度影响伪影水平。实际使用时，这些参数应该从扫描协议的元数据里读取，而不是手动输入。
从那以后我每次做重建之前，都会先拿一个已知几何的模拟数据跑一遍，确认参数配置和坐标映射没问题，再上真实数据。这个习惯帮我省了很多调试时间。希望帮到你。
本文还有配套的精品资源，点击获取
