简介这是一份面向医学CT成像与图像重建学习者的Matlab实现资源围绕锥束CTCBCT的FDK算法展开也包含MLEM、SART、SQS等迭代重建示例。代码按FDK重建流程组织先对投影数据做加权含Parker加权再滤波最后反投影得到三维体数据同时提供投影/反投影函数、参数设置与测量数据生成脚本适合想理解解析重建与迭代重建差异的读者。压缩包共17个文件以12个m源程序为主配有两份xml环境配置、一个txt说明、一张运行截图和一份img128.mat体模数据整包仅200KB结构清晰便于快速上手。目前已有1097人浏览学习资源体积小但算法覆盖较完整可直接在Matlab中运行FDK重建及对比SART/MLEM等迭代方法是入门CBCT重建的良好参考资料。1. CBCT 投影与反投影为什么 FDK 是锥形束 CT 绕不开的基石锥形束 CT 的几何模型和传统扇形束 CT 有一个本质差异X 射线源发出的是一个锥形束探测面板上每一个像素对应的不是一条“理想直线”而是一条随空间角度变化的斜射线。如果直接用扇形束时代的滤波反投影思路去处理锥形束数据重建结果会在远离中心平面的位置出现明显的伪影越靠近锥角边缘越严重。FDK 算法之所以被广泛采用是因为它用一种近乎“聪明近似”的方式把二维扇束的滤波反投影推广到了三维锥束在工程上取得了速度和质量的平衡。搞懂 FDK 的数学假设与 Matlab 实现不只是会调一个函数而是真正理解“投影—滤波—反投影”这条链路里每一步为什么长成那样。这份资源的核心是一套完整的 3D CBCT 投影、反投影与 FDK 重建的 Matlab 实现它把数据采集几何、权重计算和重建循环拆开呈现。下面我们从几何模型开始一点点拆到代码级细节。2. 锥形束几何从扇形束到锥形束的坐标变换里藏着什么2.1 扇形束时代的一维重建逻辑在二维扇形束 CT 中X 射线源在一个平面内旋转探测器是一排线阵。某个角度下测得的一行投影数据本质上是一组沿不同射线方向的线积分。重建时先把每一行投影做一维滤波再沿着射线方向“抹回”图像域这个过程叫反投影。这里的关键在于滤波核是一维的因为每条射线都能唯一确定一个二维平面内的路径。当几何升级到锥形束探测器从线阵变成面板同一角度下采集的投影变成一幅二维图像。此时一条射线不再局限于某个中心平面而是斜着穿过三维体积。如果仍然对每一行单独做一维滤波再沿着三维射线方向反投影就会出现问题斜射线的长度和它与旋转平面的夹角都在变化直接反投影会在体积边缘留下沿 z 轴方向扩散的伪影。FDK 处理这个问题的思路是给每个体素加一个与锥角相关的权重同时让滤波仍然只发生在探测器行方向。2.2 投影坐标系的参数约定理解 FDK 的第一步是把投影几何的参数约定写清楚。常见约定是源到旋转中心的距离为 DSO源到探测器平面的距离为 DSD旋转角度记为 beta。探测器面板通常被建模为垂直于中心射线的平面面板上每个像素的位置用 (u, v) 表示其中 u 是水平方向v 是垂直方向。体素空间则用物理坐标 (x, y, z) 表示z 轴与旋转轴重合。从体素到探测器坐标的映射是 FDK 实现中最容易出错的部分。对某个旋转角 beta体素 (x, y, z) 在探测器面板上的投影位置可以用以下关系描述% 将体素旋转到探测器坐标系 xr x * cos(beta) y * sin(beta); yr -x * sin(beta) y * cos(beta); % 计算探测器坐标考虑锥形束放大率 scale DSD / (DSO - yr); u xr * scale; v z * scale;这段代码里最关键的是 scale 因子。它不是一个常数而是随 yr 变化体素离源越近投影放大率越大。很多初学实现直接用一个固定放大率代替导致远离旋转中心的结构重建后位置漂移。正确做法是在循环内对每个体素重新计算 scale这也是 FDK 反投影循环计算量大的根本原因之一。2.3 为什么锥角权重要落在滤波之后而不是之前FDK 的流程是对每个角度采集的二维投影做行方向一维滤波然后加权反投影。这里有一个容易混淆的点权重到底是加在投影上还是反投影上。标准 FDK 的权重因子是 cos 锥角的倒数它本质上是补偿斜射线路径长度变化和体素有效体积变化带来的影响。在 Matlab 实现中这个权重通常和反投影循环里的 scale 因子合并计算而不是单独对投影图像做一次乘法。我一般会这么理解滤波操作是对投影数据的空间频率成分做修正它必须发生在几何加权之前否则权重会破坏滤波核的平移不变性假设。如果你先对投影做锥角加权再做滤波重建结果会出现高频噪声放大不均匀的现象中心平面干净、边缘模糊且噪点密集。3. FDK 核心流程滤波、权重与反投影的协作方式3.1 滤波核的选择与边界效果FDK 的行方向滤波使用和扇形束类似的滤波核常用的有 Ram-Lak 核、Shepp-Logan 核和带窗函数的变体。Ram-Lak 核是理想的斜坡滤波器但会放大高频噪声Shepp-Logan 核在频域里乘了一个 sinc 窗能抑制部分噪声。Matlab 实现里滤波核通常预计算成一维数组然后用卷积或频域乘法完成滤波。一个值得注意的细节是投影数据的补齐方式。滤波器对信号边缘敏感如果直接把投影行做快速傅里叶变换再乘滤波核会产生振铃效应。常见做法是把每行投影补零到下一组 2 的幂次长度或者在滤波前对行数据做镜像延拓。镜像延拓比补零效果好因为它避免了信号在边界处突变。3.2 反投影循环的结构与性能瓶颈FDK 反投影是一个典型的多层循环遍历投影角度、遍历体素、计算投影坐标、采样投影值并累加。在 Matlab 中直接写三层 for 循环会很慢常见的优化手段是把体素网格向量化并在角度循环内部用矩阵运算处理所有体素。下面是一个简化的核心循环结构% 预分配体素空间 volume zeros(Nx, Ny, Nz); for beta_idx 1:numAngles beta angles(beta_idx); % 读取当前角度的滤波后投影 P filteredProjections(:,:,beta_idx); % 对每个体素计算探测器坐标并采样 for ix 1:Nx for iy 1:Ny x xGrid(ix, iy); y yGrid(ix, iy); xr x * cos(beta) y * sin(beta); yr -x * sin(beta) y * cos(beta); scale DSD / (DSO - yr); u xr * scale; v zGrid * scale; % zGrid 是向量 % 从投影图像采样累加到体素 volume(ix, iy, :) volume(ix, iy, :) ... interp2(uGrid, vGrid, P, u, v, linear, 0) * scale * scale; end end end这个循环里的 scale * scale 不是笔误。FDK 的权重包含两个部分一个是基于几何放大率的面积修正另一个是锥角余弦权重的简化表达。在源到旋转中心距离远大于重建视野的情况下这两个因子可以合并为一个 scale 的平方项。如果你的实现用了不同的几何参数约定这里的权重形式可能变成 DSO / (DSO - yr) 或加上锥角的余弦项具体要跟投影生成过程的约定保持一致。3.3 采样插值方式对重建质量的影响反投影时体素投影到探测器坐标往往不是整数坐标必须从投影图像中插值采样。最近邻插值最快但会产生块状伪影双线性插值是默认选择更高阶的样条插值能进一步平滑但计算代价显著上升。Matlab 的 interp2 默认不支持 GPU 加速如果数据量很大一个角度一百万体素的插值开销就很可观。实际项目中我通常会先检查插值函数对边界外坐标的处理。CBCT 重建体积的角落很容易投影到探测器视野之外此时插值结果如果是 NaN累加后整个体素都会变成 NaN。常见的解决方法是把插值范围外的值设为零或者提前判断坐标是否越界。资源里的代码在我印象里用的是 interp2 加显式边界归零这是比较稳妥的做法。4. Matlab 实现里的工程细节从数据生成到重建验证4.1 投影生成模块的几何一致性FDK 重建需要和投影生成使用完全相同的几何参数否则重建结果会出现系统性错位。资源里的实现通常包含一个 forward projection 模块用于生成模拟投影数据。这个模块的核心是与反投影互逆的坐标映射给定体素空间的一个点计算它在探测器面板上的投影位置然后根据体素强度累积到探测像素上。这里有一个常见的错误投影生成时把体素当成连续点而反投影时又把探测器坐标当成连续值两者对“体素中心”的定义不一致。正确的做法是明确体素中心坐标的计算公式让正投影和反投影访问同一套坐标网格。否则即使滤波和权重都正确重建出的图像也会有半个体素大小的偏移。4.2 FDK 重建结果的验证指标拿到重建体积后第一件事不是看图而是做定量验证。常用的做法是放置若干个已知密度和高斯形状的模拟目标重建后对比峰值位置和半高全宽。对于 CBCTz 方向的半高全宽尤其值得关注因为它直接反映锥角伪影的严重程度。验证时还要注意滤波核的归一化。FDK 重建出来的数值不是原始的衰减系数而是经过滤波反投影后的相对值。不同滤波核、不同插值方式都会影响绝对数值。如果你期望重建值和真实衰减系数成线性关系需要用均匀模体做一次定标。这个定标过程在模拟环境里很容易在真实系统里则复杂得多因为涉及射束硬化等物理效应。4.3 资源代码的模块划分与扩展点拿到这份资源后建议先读它的目录结构搞清投影、滤波、反投影分别在哪几个文件里。一个设计良好的 Matlab 实现通常把几何参数集中放在一个配置脚本里而不是散落在各个函数中。你在跑通示例后可以尝试修改探测器像素尺寸、旋转角度数量、体素分辨率观察重建质量变化这是理解 FDK 参数敏感性的最快路径。修改时有一个细节值得注意角度采样间隔和探测器像素大小必须满足采样定理。如果投影角度数太少重建会出现星芒状伪影如果探测器像素太大高频信息丢失边缘变模糊。这些现象不是 FDK 算法的缺陷而是采样不足的必然结果。在模拟环境里你可以任意提高采样率但真实系统里探测器像素尺寸和曝光时间有物理限制这正体现了重建算法和硬件设计的耦合。5. FDK 实现避坑指南三条真实踩过的坑与排查路径5.1 重建体积上下颠倒或左右镜像现象模拟一个偏中心的高密度球重建后球的位置整体偏移到镜像位置。原因投影坐标的 y 轴方向和体素坐标的 y 轴方向定义不一致。Matlab 的矩阵索引是从上到下、从右到左而物理坐标通常是左手法则或右手法则。投影生成用了右手系反投影却按矩阵索引顺序解释两者一组合就出现镜像。解决在代码开头统一坐标轴方向并用手动构造的对称模体做一次全流程验证。比如放置四个角落都有标记的模体重建后检查标记的对应关系。这个验证只需要几分钟能省下后面几小时的排查时间。5.2 重建结果在远离旋转中心处出现拉伸伪影现象中心平面图像正常但 z 轴方向越靠外的层面结构变得越模糊且强度衰减。原因锥角权重处理有误。FDK 的锥角补偿本质上是把斜射线的投影值映射回体素时乘一个与锥角相关的系数。如果权重缺失或计算位置错误远离中心平面的体素累积的投影值偏小且不均匀视觉上就是拉伸和模糊。解决把锥角权重单独提取出来打印它在不同 z 层的数值分布确认它随 |v| 增大而增大。如果权重曲线异常回到几何推导重新检查 scale 因子的定义。5.3 滤波后投影图像边缘出现黑白条纹现象滤波后的投影图像在左右边缘有竖直的亮暗条纹重建后边缘出现环形伪影。原因滤波前没有做信号延拓直接补零导致频域边缘振铃。解决改用镜像延拓或者使用 Matlab 的 filtfilt 类零相位滤波思路。需要注意 filtfilt 只适用于一维信号且不能直接处理投影数据中的 NaN 区域所以更常见的做法是手写延拓函数。6. 设计取舍近似带来的效率和效率背后的边界FDK 不是一个精确重建算法这一点在资源文档里可能不会重点强调但理解它的近似性质比记住公式更重要。精确的锥形束重建需要满足数据完备性条件即每个平面都必须有足够多的射线穿过这在标准圆轨道 CBCT 中是无法满足的。FDK 选择忽略这些缺失的数据用加权反投影直接构造一个近似解代价是远离中心平面的灰度误差和伪影收益是算法结构简单、计算量可接受且对短扫描也基本鲁棒。从工程角度看这个取舍是明智的。医学影像和工业检测场景中视觉诊断和缺陷定位往往比定量绝对密度更关键。FDK 提供的图像质量在大多数场景下够用而更复杂的迭代重建算法虽然能改善锥角伪影但计算时间从分钟级上升到小时级对实时引导类应用不现实。资源里还给出了一个值得深思的扩展方向如果扫描轨道变成螺旋或双圆数据完备性条件满足后FDK 的权重公式是否需要修改以及如何扩展。我的经验是先跑通圆轨道 FDK再扩展轨道类型比直接上手迭代重建要可靠得多。从那以后我每次做 CBCT 模拟都会先构建一个小体积的数学模体用这份资源里的实现跑通投影和重建全流程确认几何没有翻转和缩放问题再切换到大体积数据。这种习惯让我避开了很多因为一个符号错误而重算整夜的窘境。希望帮到你。本文还有配套的精品资源点击获取