简介三维锥束CT在医学成像中应用广泛这份源代码与Matlab示例为投影、反投影、FDK滤波反投影及MLEM迭代重建提供了完整的可运行实现。资源面向医学影像相关专业的学生、科研人员和算法初学者既可用于理解Radon变换、锥束几何校正、滤波与反投影的数学原理也能通过实际Demo观察不同重建算法的效果差异。压缩包共13个文件主体为10个m脚本分别演示投影数据生成、FBP重建和MLEM迭代过程另有2个mat文件保存预生成的Phantom体模及投影数据1个txt为使用说明整体体积仅73KB轻量且易于下载。目前已有884人浏览学习无论是课程实验还是科研预研都可基于这些代码快速搭建CBCT重建流程进一步对比FDK与MLEM在噪声抑制、伪影处理及收敛速度上的表现是一份兼具教学与参考价值的实践资料。1. 从投影到重建为什么这套CBCT示例值得先跑一遍拿到一套三维锥束CT重建源码最容易踩的坑是先跑Demo3_MLEM.m然后盯着迭代过程发呆。真正合理的顺序是先把投影和反投影模型吃透再去看FDK和MLEM。这套Matlab示例包恰好按这个链路拆分Demo1_projection.m负责从Phantom64.mat这类数字体模生成锥束投影Demo2_FBP.m实现FDK滤波反投影Demo3_MLEM.m走统计迭代重建。每个demo都带.mat体模数据和注释学生可以直接跑通工程师也能用它验证投影算子、反投影算子是否自洽。我拆这套代码时最深的感受是FDK和MLEM的差距不在数学公式而在对锥束几何和采样密度的敏感度。下面按三个demo展开最后给出替换成自己数据时的调试建议。2. 投影与反投影先把锥束几何和Matlab数据组织讲透2.1 从Radon变换到锥束投影三个几何参数决定一切二维CT的投影可以用Radon变换描述三维锥束CT则是Radon变换在锥形束几何下的扩展。X射线源绕z轴旋转每个角度下探测器采集一幅二维投影把所有角度堆叠起来就是三维投影数据。Matlab示例里的Phantom64.mat是64×64×64的体模对应一组离散的衰减系数投影数据则是n_angle个二维角度视图每个视图尺寸由探测器行列数决定。实际编程前需要对齐三组参数旋转角度序列、射线源到旋转中心的距离source distance、探测器像素间距。常见做法是把角度序列定义为0到360度均匀分布间隔越小重建质量越高但投影数据量和计算时间线性增长。资源里的Demo1_projection.m会生成一组投影保存成多维数组或.mat文件供后续重建调用。下面是我的调用模板% 从体模生成锥束投影 load(Phantom64.mat); % image_3d: 64x64x64 衰减系数体模 n_angle 360; % 总投影角度数 det_rows 64; % 探测器行数竖直方向 det_cols 64; % 探测器列数水平方向 sad 100; % source-axis distance单位像素 proj zeros(det_rows, det_cols, n_angle); for ia 1:n_angle theta (ia - 1) * 360 / n_angle; % 调用项目自带的投影函数输入体模和几何返回一个角度视图 proj(:, :, ia) projectConeBeam(image_3d, theta, sad, det_rows, det_cols); end这段代码里sad是射线源到旋转中心的距离它会直接影响锥束权重和反投影坐标映射det_rows和det_cols不一定要等于体模尺寸但采样密度低于体模分辨率时FDK重建出来的边缘会明显模糊。参数说明角度序列用0到360度而非0到180度是因为锥束扇形覆盖一圈才能获得完整采样投影函数名并不一定叫projectConeBeam实际以readme.txt中的函数名为准我这里只是说明常见的参数组织形式。2.2 Demo1_projection.m投影循环里到底发生了什么打开Demo1_projection.m你会发现它并不是简单地调用现成函数而是把“计算每条X射线穿过体模的线积分”这件事显式拆成了若干步。核心逻辑是对每个角度将体模旋转到源-探测器坐标系下沿着射线方向做插值并累加衰减系数。这一步如果写不好后面FDK和MLEM都会跟着错。2.2.1 正向投影的四种实现方式教学代码通常不是用最简洁的向量化写法而是用三重循环方便初学者对照数学公式。常见实现有四条路坐标旋转法、射线遍历法类似Bresenham、体素驱动法、基于傅里叶切片定理的近似法。Demo1_projection.m大概率会采用坐标旋转法因为思路最直白function p projectOneView(vol, theta, sad, nu, nv) % 将体模旋转到当前角度下的源-探测器坐标系 % vol: 三维体模; theta: 旋转角度; sad: 源到旋转中心距离 % nu, nv: 探测器行列数 vol_rot imrotate3(vol, theta, [0 0 1], linear, crop); % 绕z轴旋转 p zeros(nu, nv); for v 1:nv for u 1:nu % 源位置在球坐标外用sad控制 % 射线方向由(u,v)和源位置确定 % 这里简化为沿x方向线积分 p(v, u) sum(vol_rot(:, round(u), round(v))); end end这个写法非常简化真实计算要考虑探测器像素坐标到体素坐标的映射以及衰减系数沿射线长度的积分权重。imrotate3是Matlab图像处理工具箱的函数不是资源自带的这里只是为了说明旋转步骤。实际项目里会把旋转改成直接在射线追踪时计算坐标变换避免三维插值带来的额外模糊。2.2.2 反投影的边界效应与伪影来源反投影是重建逆过程的直观体现把每一个角度视图沿原射线方向“涂回”体模空间。如果直接把投影数据做反投影而不滤波得到的是模糊重叠的物体轮廓而且边界处会出现明显的星形伪影。这是因为反投影本质上是把高频分量放大却没有补偿采样密度。几乎所有CBCT教材都会强调解析重建的核心是“先滤波、再反投影”滤波器的选择直接决定图像锐度和噪声水平。对于这套Matlab示例反投影部分通常被封装在FDK和MLEM的迭代里单独运行Demo1只生成投影。但我习惯在跑重建前做一个“圆环闭合测试”把投影数据先反投影再重新投影看两次投影残差是否在合理范围。这个测试能快速暴露坐标系方向不一致、角度单位错误、探测器行序颠倒等低级问题。3. FDK重建Demo2_FBP.m 里的滤波反投影流水线3.1 FDK三步骤锥束加权、滤波、三维反投影FDK算法是Feldkamp、Davis和Kress在1984年提出的是目前商用CBCT最常用的解析重建方法。它不是精确重建而是在扇束FBP基础上做了锥角近似当锥角较小时重建误差可控。核心分三步对每个角度视图做锥束加权沿探测器水平方向做一维滤波再把滤波后的视图三维反投影到体模空间。锥束加权的意义在于校正X射线因倾斜路径带来的强度衰减差异。权重因子与源到探测器的距离、像素偏离中心射线的程度有关。滤波步骤绕不开Ram-Lak滤波器斜坡滤波和它的窗函数变体。下面是滤波器对照表也是我在实验中常用的组合滤波器频率响应特点适用场景伪影表现Ram-Lak全频段线性放大理想无噪声投影噪声放大明显分辨率高Shepp-Logan高频衰减低剂量投影边缘过冲减小噪声降低Cosine平滑衰减临床软组织成像图像更平滑分辨损失可接受Hamming/Hann中高频平滑一般教学示例噪声和分辨率折中Demo2_FBP.m里大概率用的是Ram-Lak或Shepp-Logan因为这两种实现最简单只要在频域乘上一条斜坡线。注意二维探测器只需要沿水平方向对应扇束内角度方向滤波竖直方向不能随意滤波否则会破坏锥束几何关系。3.2 把FDK写成Matlab循环权重、滤波器与累加我见过多个版本的FDK实现最稳定的结构是先预计算权重表再对每个视图做一维FFT滤波最后三重循环反投影累加。结构如下function vol fdk_recon(proj, angles, sad, sdd, vol_size) % proj: [nu, nv, n_angle] 投影数据 % angles: 每个投影对应的角度(度) % sad: 源到旋转中心距离 % sdd: 源到探测器距离通常sdd sad % vol_size: 重建体模尺寸例如[64 64 64] vol zeros(vol_size); N size(proj, 2); % 探测器水平采样数 for ia 1:length(angles) view double(proj(:, :, ia)); % Step 1: 锥束权重 for u 1:N gamma atan((u - N/2 - 1) / sdd); % 像素相对中心射线的角度 weight cos(gamma); view(:, u) view(:, u) .* weight; end % Step 2: 水平方向斜坡滤波 filt abs(fft((-N/2:N/2-1)/N)); % 生成Ram-Lak频域响应 view_f fft(view, [], 2); view_f view_f .* filt; view_f real(ifft(view_f, [], 2)); % Step 3: 反投影到体积 for z 1:vol_size(3) for y 1:vol_size(2) for x 1:vol_size(1) % 计算该体素在探测器上的投影位置 % 然后从view_f中插值取得贡献 vol(x, y, z) vol(x, y, z) ...; end end end end这段代码的过滤频响写法是为了展示原理实际使用时应该用ifftshift调整直流分量位置。权重生成里sdd必须与投影生成时的几何保持一致否则重建出的图像会有杯状伪影。反投影那部分省略了坐标变换完整实现要先把体素坐标转换到全局坐标再投影到探测器平面。Demo2_FBP.m里的循环顺序不一定相同但如果你打开代码发现用interp2做探测器插值那是正常的双线性插值能明显减少角度混叠。3.3 重建质量检查从横断面到锥角伪影FDK跑完后第一件事不是看三维渲染而是切三个正交面。用Matlab的squeeze和imshow直接看x/y/z三个方向切片figure; slice_xy squeeze(recon(:, :, 32)); % z32层横断面 subplot(1,3,1); imshow(slice_xy, []); title(Transverse); slice_xz squeeze(recon(:, 32, :)); % y32层矢状面 subplot(1,3,2); imshow(slice_xz, []); title(Sagittal); slice_yz squeeze(recon(32, :, :)); % x32层冠状面 subplot(1,3,3); imshow(slice_yz, []); title(Coronal);注意观察两个地方边缘是否有拖尾和负值远离旋转中心的地方是否灰度下陷。拖尾说明滤波器设计有偏负值可能是滤波后反投影叠加产生的截断伪影这在FDK里很常见。锥角伪影通常出现在体模上下两端图像会变淡或产生弧线原因是FDK对大锥角重建不精确。资源里的Phantom64是64体素尺寸锥角不大伪影不会太明显如果换成128体素差异会更明显同时计算时间翻数倍。4. MLEM迭代重建Demo3_MLEM.m 的统计模型与收敛控制4.1 为什么用MLEM噪声模型与迭代更新方程FDK是解析算法一次滤波反投影出结果速度快但对投影噪声敏感。MLEM的全称是Maximum Likelihood Expectation Maximization它把投影过程建模为泊松随机过程通过迭代寻找最符合测量投影的体模衰减系数分布。数学上迭代公式是图像新估计 图像旧估计 / 灵敏度 × 反投影(测量投影 / 前向投影(图像旧估计))。这个更新公式看起来简单但每一步都涉及大量计算前向投影要重新走一遍射线追踪反投影也要把误差分布回体模。这也是为什么MLEM在Demo3_MLEM.m里跑起来明显比FDK慢。好处是迭代重建可以建模噪声统计特性在低剂量投影下能保持相对干净的图像还能加入先验约束。对只想跑通示例的学生来说先明白迭代公式就可以了真正的坑在收敛判断和灵敏度矩阵。4.2 Demo3_MLEM.m 的迭代主体变量、循环与收敛判断教学代码里的MLEM实现通常会预先算好一个灵敏度矩阵sensitivity map表示每个体素在所有角度下被射线穿过的总权重避免在每次迭代里重复计算。迭代主体如下function vol mlem_recon(proj, geom, n_iter) % proj: 投影数据 [nu, nv, n_angle] % geom: 几何参数结构体包含sad, sdd, angles, det_size % n_iter: 预设迭代次数 sensitivity computeSensitivity(geom); % 全1体模前向投影 vol ones(geom.vol_size); % 初始化为全1 for iter 1:n_iter % 前向投影模拟当前估计产生的投影 proj_est forwardProject(vol, geom); % 比值测量投影 / 估计投影处理分母为0 ratio proj ./ (proj_est eps); % 反投影误差 correction backProject(ratio, geom); % 更新体模 vol vol ./ sensitivity .* correction; % 计算投影差异便于观察收敛 residual sum((proj_est(:) - proj(:)).^2); fprintf(Iter %d: residual %.6f\n, iter, residual); end这里的computeSensitivity、forwardProject、backProject在资源里分别对应不同的函数名但逻辑一致。sensitivity矩阵中接近0的位置在重建体外更新时除以0会产生NaN所以加上eps或掩膜处理。迭代次数不是越多越好通常20到50次视觉质量提升明显100次后主要是高频噪声被放大。我看Demo3_MLEM.m时习惯把fprintf改成记录residual数组画一条收敛曲线比直接看每一轮图像方便得多。4.3 FDK与MLEM对比从sinc伪影到剂量优化两种方法在同一个Phantom64数据上跑出来的结果差异比很多入门书描述的更大。FDK的优点是速度快、参数少、结果可预期缺点是解析反投影会放大高频噪声且锥角大于5度时边缘结构产生几何失真。MLEM的优点是在低剂量投影下依然能保持轮廓清晰缺点是需要调迭代次数并且每次迭代都相当于一次FDK的工作量。对比项FDKMLEM计算时间单轮循环秒级数十轮迭代分钟级噪声抑制依赖滤波器窗函数迭代自然抑制统计噪声几何误差锥角大时明显可建模校正收敛判断无需需要残差或图像差异阈值教学难度低中高在我自己的实验里对于Phantom128.mat这批数据FDK配合Shepp-Logan滤波器在视觉上接近MLEM迭代20次的结果但MLEM在低对比度的细节上更稳。如果你做剂量优化实验用FDK做初始估计再切换到MLEM能省掉前几次迭代这是一个实用技巧。5. 把示例改造成自己的CBCT数据参数替换与调试技巧5.1 几何参数对齐探测器间距、旋转中心与角度的坑换数据时最容易错的是几何参数不匹配。Demo系列用的是Phantom64.mat和Phantom128.mat体模坐标轴是等间隔体素。真实CBCT投影数据通常来自平板探测器几何会涉及探测器像素物理尺寸、源到探测器的距离、旋转中心在探测器上的投影位置。常见做法是先做几何校准用一根细针或钢珠在不同角度投影反算几何参数。如果跳过这一步FDK重建的切片会有重影或同心圆伪影。Matlab里调整旋转中心偏移量很简单在反投影索引上加一个offset即可但这个offset差一个像素都会让重建图像边缘发毛。5.2 用投影误差曲线判断迭代是否收敛修改Demo3_MLEM.m时不要只看重建图片要把每轮前后的投影差异记录下来。我一般会在代码里加一个向量err_history存每个iteration的均方误差跑完后画semilogy。正常情况曲线单调下降并趋于平缓如果曲线震荡说明投影几何里存在不一致比如角度顺序颠倒或探测器方向翻转。若曲线下降缓慢检查灵敏度矩阵是否包含0值导致体模内部产生空洞。这里给出简易判断代码err_history []; for iter 1:n_iter % ... 原有迭代代码 ... err_history(end1) mean((proj_est(:) - proj(:)).^2); end semilogy(err_history); grid on; xlabel(Iteration); ylabel(MSE of projection);如果MSE在第10次已经不再下降再迭代只是过拟合投影噪声可以提前终止。批量试验时用这个曲线作为每次运行的存档比保存全部重建体积省空间。5.3 性能优化从三重循环到GPU并行处理这套示例的原始循环适合教学不适合大规模实验。我在实际项目里的做法是先用单一组几何参数跑通再针对热点函数做三件事把不依赖角度数的权重矩阵提前计算将反投影中的interp2替换为查表式索引最后用parfor按角度并行。Matlab的Parallel Computing Toolbox可以这样处理parfor ia 1:numel(angles) view_proc processView(proj(:, :, ia), geom); vol_contrib backProjectSingle(view_proc, geom, angles(ia)); % 每个角度对体积的贡献单独保存最后累加 end注意parfor里不能直接更新同一个vol变量必须把每个角度的贡献存成cell数组再累加否则Matlab会报错或产生不可预测的结果。GPU加速MLEM时前向投影和反投影都换成gpuArray但注意矩阵大小要适合显存Phantom128已经接近常见显卡上限我之前用8GB显存跑有点紧张降低数据类型到single能缓解。本文还有配套的精品资源点击获取