把滤波反投影算法在Matlab里完整跑通一遍是我自己做CT成像仿真印象最深的一次。以前看教材总觉得“滤波反投影”就是个名词真到动手写代码才发现生成投影信号、对信号做傅里叶变换、频域滤波、再反投影重建每一步都有大量细节值得琢磨。这其实是一个非常适合练手的仿真项目。它不依赖昂贵的实验设备不需要真实的CT扫描仪用Matlab就能把一套完整的CT成像链路走通。无论是学医学影像、做无损检测还是研究图像重建算法这个项目都能帮你把“探测器采集什么、频域里发生什么、重建算法解决什么”这些底层问题彻底搞清楚。我会把整个项目拆开讲清楚包括算法选型、核心原理、完整代码、参数怎么调、常见的坑是什么以及后续可以往哪个方向扩展。1. 为什么我用滤波反投影算法做CT成像仿真1.1 解析重建与迭代重建的取舍CT重建的核心问题是从一组X射线穿透物体后的投影数据中反推出物体内部的衰减系数分布。解决这个问题主要有两大路线解析重建和迭代重建。滤波反投影算法Filtered Back ProjectionFBP属于解析重建是临床上绝大多数传统CT机都在用的方法。我选择FBP做这个仿真项目原因很简单一是它和傅里叶变换、Radon变换这些数学工具紧密结合理解FBP相当于同时打通了信号处理和图像重建两条线二是它的计算效率高一次完整重建在普通电脑上几秒钟就能完成迭代重建往往需要几十次甚至上百次迭代调试起来非常痛苦。当然FBP也有限制比如它对噪声敏感、投影数据不足时容易出现星形伪影。但作为教学仿真和算法验证平台FBP几乎是性价比最高的选择。你先把FBP吃透后续再理解迭代重建里的系统矩阵、正则化项会轻松很多。1.2 仿真链路从模型到图像重建整个CT成像仿真项目可以拆成四条主线构造一个数字人体断层模型常用的就是Shepp-Logan头颅模型模拟X射线从多个角度穿过模型记录投影信号这一步对应实际扫描中的数据采集对每个角度下的投影信号做一维傅里叶变换在频域使用滤波器进行补偿再变换回空间域把所有处理后的投影数据沿各自角度“拉回”到图像空间叠加得到重建图像。在这个过程中最核心的理论支撑是中心切片定理Central Slice Theorem。你不需要把公式背得多熟练但一定要理解这个定理的含义某角度下的投影信号它的一维傅里叶变换恰好等于目标断层图像二维傅里叶变换中过原点的一条径向线。有了这个关系二维图像重建问题就能被拆解成一系列一维信号处理问题来处理。仿真项目的魅力在于你能把每一个环节单独拎出来观察、打印、验证。我自己做的时候会把投影信号画出来、把滤波前后的频谱画出来、把中间每一步都可视化这样一遍走下来比死记十遍公式都管用。2. 滤波反投影怎么一步步实现核心原理解析2.1 投影信号它是怎么生成的又代表了什么在平行束CT几何中X射线从某个方向照射物体探测器记录的是这条路径上所有物质对X射线的衰减积分值。数学上这就是Radon变换( P_{\theta}(t) \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x, y) , \delta(x\cos\theta y\sin\theta - t) , dx , dy )公式看着复杂你可以把它理解成“沿着某一方向把图像压扁”。就好比你从侧面去看一个立体物体看到的是一个降了一维的轮廓不同角度看到的轮廓各不相同。CT扫描的过程就是收集多个角度下“压扁”的结果。在Matlab里radon函数可以直接完成这一步。它会返回一个二维矩阵每一列对应一个投影角度每一行对应探测器上的一个位置单元。默认探测器单元数会比图像尺寸大一些因为旋转时图像角落也会超出原始尺寸范围。这个细节在写自定义反投影时尤其重要探测器坐标轴的范围必须跟radon返回的xp变量保持一致不然反投影时插值坐标会对不上重建图会出现明显错位。2.2 中心切片定理为什么必须先滤波再反投影如果不做任何处理直接把原始投影数据原封不动地反投影回图像空间你会发现重建出来的图像糊成一片边缘完全不清晰。这是因为反投影本质上是对整个频域做了半径影响的加权叠加它会把低频分量过度放大相当于给图像加了一个和频率成反比的权重。中心切片定理给出了答案。每组角度的投影傅里叶变换只是原图二维频谱中的一条线用有限个角度去覆盖整个二维频域越远离中心的高频区域越稀疏。为了让重建图像逼近原图必须对投影信号做频域补偿补偿的权重大致与频率绝对值成正比也就是|ω|关系这个滤波器被称为斜坡滤波器ramp filter。“先滤波再反投影”这个名字就是这么来的。滤波过程就是让每个角度的投影信号在频域上乘一个绝对值频率函数反投影过程则是把所有滤波后的投影信号沿原方向回抹到图像空间。两者合在一起才能尽可能还原被模糊掉的高频细节。2.3 滤波函数选型Ram-Lak还是Hamming斜坡滤波器本身是理想化的它在频域上的形式就是一条从负频率到正频率的|ω|折线。如果直接用这个滤波器重建图像虽然清晰但一旦投影数据里含有噪声高频噪声也会被同步放大结果就是图像出现密密麻麻的颗粒感。实际工程里通常会在斜坡滤波器前面乘一个窗函数来抑制高频。这就好比手机拍照的美颜滤镜保留主要细节同时压掉一部分高频噪点。常见的组合如下滤波器名称函数形式特点适用场景Ram-Lak|ω|分辨率最高噪声也最大无噪或低噪仿真数据Shepp-Logan|ω|·sinc(ω/2)振铃小噪声抑制适中教学仿真、通用场景Cosine|ω|·cos(ω/2)平滑过渡噪声抑制较好含噪投影数据Hamming|ω|·(0.540.46cosω)噪声大幅抑制边缘稍糊真实低剂量重建我在这个项目里默认使用Ram-Lak因为仿真数据本身没有噪声能最大程度考察算法本身的重建能力。如果你往投影里加了高斯噪声或者泊松噪声建议切换成Hamming再跑一遍对比效果会非常直观。3. Matlab代码实现从投影到重建的完整流程3.1 第一步生成仿真模型与投影信号Matlab的Image Processing Toolbox自带phantom函数可以生成Shepp-Logan模型。这个模型由十几个椭圆组成模拟了头颅断层中不同组织区域的衰减差异是CT重建领域最经典的测试图像。% 生成256x256的数字断层模型 I phantom(modified shepp-logan, 256); imshow(I, []); title(原始模型); % 设置投影角度这里用0度到178度步长2度共90个角度 theta 0:2:178; % 计算投影信号P矩阵的行对应探测器单元列对应角度 [P, xp] radon(I, theta); % 画一下90度方向的投影信号感受一下 figure; plot(xp, P(:, find(theta 90))); xlabel(探测器位置); ylabel(投影值); title(90度方向的投影信号);xp是radon函数返回的探测器坐标轴。给radon加两个输出参数时它才会返回坐标轴很多人第一次用会漏掉这个然后自定义反投影时找不到正确的距离轴只能在代码里瞎猜这是没必要踩的坑。3.2 第二步信号傅里叶变换与斜坡滤波投影信号本质是一维信号要做频域滤波就得先做傅里叶变换。注意我这里用的是fft配合fftshift/ifftshift目的是把频率0放到序列中心这样构造斜坡滤波器时可以直接用与频率轴对应的一维数组。% 频域斜坡滤波器构造函数 n size(P, 1); freq linspace(-1, 1, n); % 归一化频率轴范围[-1,1] ramp abs(freq); % 理想斜坡滤波器 % 对每个角度的投影信号做傅里叶变换并滤波 Pf zeros(size(P)); for k 1:length(theta) proj P(:, k); PROJ fftshift(fft(ifftshift(proj))); % 一维傅里叶变换 PROJ_filt PROJ .* ramp(:); % 频域乘斜坡滤波器 Pf(:, k) fftshift(ifft(ifftshift(PROJ_filt))); % 逆变换回空间域 end % 看看滤波前后的差异 figure; plot(xp, P(:,1), b, xp, real(Pf(:,1)), r); legend(滤波前, 滤波后); xlabel(探测器位置); ylabel(幅度); title(第一个角度的投影信号滤波前后对比);这里有一个很容易被忽略的细节滤波之后得到的信号可能是复数值这是因为频域相乘后相位没有完全归零实际使用中应该取实部。如果你发现重建图像有细微的纹波先怀疑一下是不是忘了取real()。理论上也可以在空间域用卷积实现同样效果即把投影信号跟斜坡滤波器的空间域核函数做卷积。频域实现的好处是更直观、更容易控制滤波频率范围也方便切换不同窗函数。比如想用Hamming窗只需要在乘滤波器时再乘一个窗函数数组% 以Hamming窗为例 window 0.54 0.46 * cos(linspace(-pi, pi, n).); ramp_hamming abs(freq). .* window; PROJ_filt PROJ .* ramp_hamming;3.3 第三步反投影重建滤波做完接下来就是反投影。这一步的核心操作是对图像平面上的每一个像素点计算它在这个角度下对应的探测器位置然后从滤波后的投影数据中取出该位置的值累加到像素上。% 初始化重建图像 recon zeros(size(I)); % 以图像中心为原点建立坐标网格 N size(I, 1); half N / 2; [xg, yg] meshgrid(linspace(-half, half, N), linspace(-half, half, N)); % 角度步长弧度用于离散积分的权重 dtheta deg2rad(theta(2) - theta(1)); % 逐角度反投影累加 for k 1:length(theta) % 当前角度下每个像素对应的探测器坐标 t x*cos(theta) y*sin(theta) t_coord xg * cosd(theta(k)) yg * sind(theta(k)); % 用插值从滤波投影数据中取像素对应值 proj_interp interp1(xp, real(Pf(:, k)), t_coord(:), linear, 0); proj_interp reshape(proj_interp, size(xg)); % 累加 recon recon proj_interp; end % 乘上角度间隔完成离散积分近似 recon recon * dtheta; figure; imshow(recon, []); title(滤波反投影重建结果);这个循环就是整个项目最核心的代码段了。interp1里的最后一个参数0表示当像素点对应的探测器坐标超出范围时按0处理这样能减少边界外无效区域的干扰。如果你运行完这段代码发现重建图像和原始模型对比亮度和对比度差了一大截别急着改代码先把dtheta这个角度步长因子加上。反投影过程本质上是在做一个离散积分近似漏掉这个因子重建图像的整体幅度会错得离谱。3.4 第四步重建质量评估跑完重建不能只看图要量化评估。我用两个指标衡量重建效果一是均方误差(MSE)二是峰值信噪比(PSNR)。因为重建图像的幅度尺度和原图不完全一致评估前需要先把重建图缩放到和原图相同的数值范围。% 将重建图像缩放到[0,1]区间 recon_norm mat2gray(recon); I_norm mat2gray(I); % 计算MSE和PSNR mse_val mean((I_norm(:) - recon_norm(:)).^2); psnr_val 10 * log10(1 / mse_val); fprintf(MSE %.4f\n, mse_val); fprintf(PSNR %.2f dB\n, psnr_val);在我自己的测试里用256分辨率模型、90个角度、Ram-Lak滤波MSE大约在0.001量级PSNR在28dB附近。如果角度增加到180个PSNR能往上走3到4个dB。这里想提个醒PSNR在CT重建里只能作为参考指标因为它对整体灰度缩放非常敏感真正评价重建质量最好同时输出重建图和原图的差值图直观查看误差分布。4. 实操中遇到的问题排查与参数调试记录4.1 星形伪影严重怎么回事最典型的问题是重建图像上出现明显的放射状条纹从中心向四周扩散看起来像星星的光芒。这个现象的本质是角度采样不足。滤波反投影理论要求从0到180度连续采集投影实际中只能采集有限个角度。每个角度在频域中对应一条径向线角度太少频域覆盖就会出现缺口这些缺口在图像空间就表现为放射状伪影。解决办法很简单把角度步长从2度改成1度甚至0.5度伪影会明显减少。代价是计算量线性增长仿真环境下一般都能接受。另外探测器单元数太少也会造成类似问题。探测器分辨率不够相当于投影信号本身就被“糊”了一层反投影时再精细也没用。4.2 重建图像偏暗、灰度不对这个坑我印象太深了。第一次跑完反投影重建图像整体灰蒙蒙的对比度极低怎么看怎么不对劲最后发现是漏了角度步长因子。反投影的离散公式是( f_{recon}(x,y) \approx \sum_{k1}^{N} P_{\theta_k}^{filtered}(x\cos\theta_k y\sin\theta_k) \cdot \Delta\theta )如果你累加之后没有乘上dtheta等价于每个角度贡献的积分权重是1而不是dtheta重建结果就会偏大或偏小很多。另一个常见原因是斜坡滤波器的频率轴定义不对比如用了linspace(0, 1, n)而不是linspace(-1, 1, n)相当于滤波器没有负频率成分重建结果会差很多。4.3 投影角度和采样数怎么选这是做仿真时必问的问题。我的经验值是图像尺寸为N时探测器单元数大约取N*sqrt(2)因为旋转对角线会把模型的最长跨幅暴露出来投影角度数取和探测器单元数相近的值即可一般128到256个角度效果就不错了。可以自己做个简单实验把投影角度设为10、30、90、180、360五个档次分别重建比较PSNR。你会看到一个很明显的趋势——角度少时PSNR低、伪影重到180角度以后增长开始变缓360角度时伪影基本肉眼不可见。这个实验做完你对“角度分辨率”的理解会比看任何教材都深刻。4.4 重建速度太慢怎么办反投影循环是整个算法最耗时的地方主要是因为每个角度都用interp1对整张图像的每个像素做插值。图像尺寸256还好一旦上到512或102490个角度的循环会明显变慢。我的优化经验按性价比排序用parfor替代for多核并行处理代码改动最小提速效果明显把interp1换成内置的interp2高级用法或者用imrotate结合矩阵旋转但要注意插值精度预分配所有中间变量避免循环内动态增长数组如果追求极致性能可以把反投影核心写成MEX函数或使用GPU加速。对仿真教学来说做到第1步就够了。我自己的测试里256图像、180角度parfor开8核后耗时能从三四秒降到一秒以内。启动并行池本身有开销如果角度数少于90没必要开并行。5. 从平行束仿真到真正CT的扩展路径5.1 扇形束与锥形束的区别上面讲的是平行束几何也就是每束X射线都互相平行只适用于早期CT和教学仿真。真实临床CT用的是扇形束单排探测器或者锥形束多排探测器几何关系更复杂。从平行束换到扇形束反投影公式里要多一个距离加权因子因为扇形的每条射线路径长度不同采样密度也不同。多排探测器进一步扩展成锥形束之后就要用FDK算法这类近似重建方法。虽然公式复杂了但核心思想跟FBP完全一致频域补偿加反投影。理解了平行束FBP再往扇形束走主要工作是修几何权重。5.2 仿真代码还能怎么用这套仿真代码一旦跑通后续能玩的方向非常多。比如给投影信号加噪声研究不同滤波函数对噪声的抑制能力比如设置金属伪影场景在模型里插入高衰减区域模拟临床上常见的伪影问题比如把投影数据换成真实CT原始数据验证自己的重建流程是否还能成立。我自己最推荐的扩展是“用FBP做低剂量仿真”。你往投影里加泊松噪声然后对比Ram-Lak和Hamming窗的重建结果会直观看到噪声放大的代价和滤波带来的分辨率损失。这个实验做完你就明白为什么真实CT系统里滤波器设置那么讲究为什么低剂量扫描会催生一系列新算法。另外这个项目也很适合作为深度学习重建的基准。用FBP重建结果作为网络输入让深度学习模型去学习修正伪影和噪声是目前很多研究论文的标准做法。你想研究AI图像重建的话这套仿真代码就是最好的数据生成器。最后再说一点个人体会做这个项目最忌讳的是直接调库函数出个结果就完事。iradon确实一行代码就能重建但如果你没有亲手写过那个反投影循环你可能永远体会不到为什么滤波是必须的、为什么角度步长要乘上去、为什么插值方式会影响重建质量。我自己跑完这个项目之后再去翻CT重建的经典论文很多之前看不懂的公式突然变得非常具体。强烈建议你也动手拆一遍把每个环节都打印出来看一眼这种感觉是纯看书给不了的。
Matlab实现滤波反投影CT成像仿真：从原理到代码全解析
把滤波反投影算法在Matlab里完整跑通一遍是我自己做CT成像仿真印象最深的一次。以前看教材总觉得“滤波反投影”就是个名词真到动手写代码才发现生成投影信号、对信号做傅里叶变换、频域滤波、再反投影重建每一步都有大量细节值得琢磨。这其实是一个非常适合练手的仿真项目。它不依赖昂贵的实验设备不需要真实的CT扫描仪用Matlab就能把一套完整的CT成像链路走通。无论是学医学影像、做无损检测还是研究图像重建算法这个项目都能帮你把“探测器采集什么、频域里发生什么、重建算法解决什么”这些底层问题彻底搞清楚。我会把整个项目拆开讲清楚包括算法选型、核心原理、完整代码、参数怎么调、常见的坑是什么以及后续可以往哪个方向扩展。1. 为什么我用滤波反投影算法做CT成像仿真1.1 解析重建与迭代重建的取舍CT重建的核心问题是从一组X射线穿透物体后的投影数据中反推出物体内部的衰减系数分布。解决这个问题主要有两大路线解析重建和迭代重建。滤波反投影算法Filtered Back ProjectionFBP属于解析重建是临床上绝大多数传统CT机都在用的方法。我选择FBP做这个仿真项目原因很简单一是它和傅里叶变换、Radon变换这些数学工具紧密结合理解FBP相当于同时打通了信号处理和图像重建两条线二是它的计算效率高一次完整重建在普通电脑上几秒钟就能完成迭代重建往往需要几十次甚至上百次迭代调试起来非常痛苦。当然FBP也有限制比如它对噪声敏感、投影数据不足时容易出现星形伪影。但作为教学仿真和算法验证平台FBP几乎是性价比最高的选择。你先把FBP吃透后续再理解迭代重建里的系统矩阵、正则化项会轻松很多。1.2 仿真链路从模型到图像重建整个CT成像仿真项目可以拆成四条主线构造一个数字人体断层模型常用的就是Shepp-Logan头颅模型模拟X射线从多个角度穿过模型记录投影信号这一步对应实际扫描中的数据采集对每个角度下的投影信号做一维傅里叶变换在频域使用滤波器进行补偿再变换回空间域把所有处理后的投影数据沿各自角度“拉回”到图像空间叠加得到重建图像。在这个过程中最核心的理论支撑是中心切片定理Central Slice Theorem。你不需要把公式背得多熟练但一定要理解这个定理的含义某角度下的投影信号它的一维傅里叶变换恰好等于目标断层图像二维傅里叶变换中过原点的一条径向线。有了这个关系二维图像重建问题就能被拆解成一系列一维信号处理问题来处理。仿真项目的魅力在于你能把每一个环节单独拎出来观察、打印、验证。我自己做的时候会把投影信号画出来、把滤波前后的频谱画出来、把中间每一步都可视化这样一遍走下来比死记十遍公式都管用。2. 滤波反投影怎么一步步实现核心原理解析2.1 投影信号它是怎么生成的又代表了什么在平行束CT几何中X射线从某个方向照射物体探测器记录的是这条路径上所有物质对X射线的衰减积分值。数学上这就是Radon变换( P_{\theta}(t) \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x, y) , \delta(x\cos\theta y\sin\theta - t) , dx , dy )公式看着复杂你可以把它理解成“沿着某一方向把图像压扁”。就好比你从侧面去看一个立体物体看到的是一个降了一维的轮廓不同角度看到的轮廓各不相同。CT扫描的过程就是收集多个角度下“压扁”的结果。在Matlab里radon函数可以直接完成这一步。它会返回一个二维矩阵每一列对应一个投影角度每一行对应探测器上的一个位置单元。默认探测器单元数会比图像尺寸大一些因为旋转时图像角落也会超出原始尺寸范围。这个细节在写自定义反投影时尤其重要探测器坐标轴的范围必须跟radon返回的xp变量保持一致不然反投影时插值坐标会对不上重建图会出现明显错位。2.2 中心切片定理为什么必须先滤波再反投影如果不做任何处理直接把原始投影数据原封不动地反投影回图像空间你会发现重建出来的图像糊成一片边缘完全不清晰。这是因为反投影本质上是对整个频域做了半径影响的加权叠加它会把低频分量过度放大相当于给图像加了一个和频率成反比的权重。中心切片定理给出了答案。每组角度的投影傅里叶变换只是原图二维频谱中的一条线用有限个角度去覆盖整个二维频域越远离中心的高频区域越稀疏。为了让重建图像逼近原图必须对投影信号做频域补偿补偿的权重大致与频率绝对值成正比也就是|ω|关系这个滤波器被称为斜坡滤波器ramp filter。“先滤波再反投影”这个名字就是这么来的。滤波过程就是让每个角度的投影信号在频域上乘一个绝对值频率函数反投影过程则是把所有滤波后的投影信号沿原方向回抹到图像空间。两者合在一起才能尽可能还原被模糊掉的高频细节。2.3 滤波函数选型Ram-Lak还是Hamming斜坡滤波器本身是理想化的它在频域上的形式就是一条从负频率到正频率的|ω|折线。如果直接用这个滤波器重建图像虽然清晰但一旦投影数据里含有噪声高频噪声也会被同步放大结果就是图像出现密密麻麻的颗粒感。实际工程里通常会在斜坡滤波器前面乘一个窗函数来抑制高频。这就好比手机拍照的美颜滤镜保留主要细节同时压掉一部分高频噪点。常见的组合如下滤波器名称函数形式特点适用场景Ram-Lak|ω|分辨率最高噪声也最大无噪或低噪仿真数据Shepp-Logan|ω|·sinc(ω/2)振铃小噪声抑制适中教学仿真、通用场景Cosine|ω|·cos(ω/2)平滑过渡噪声抑制较好含噪投影数据Hamming|ω|·(0.540.46cosω)噪声大幅抑制边缘稍糊真实低剂量重建我在这个项目里默认使用Ram-Lak因为仿真数据本身没有噪声能最大程度考察算法本身的重建能力。如果你往投影里加了高斯噪声或者泊松噪声建议切换成Hamming再跑一遍对比效果会非常直观。3. Matlab代码实现从投影到重建的完整流程3.1 第一步生成仿真模型与投影信号Matlab的Image Processing Toolbox自带phantom函数可以生成Shepp-Logan模型。这个模型由十几个椭圆组成模拟了头颅断层中不同组织区域的衰减差异是CT重建领域最经典的测试图像。% 生成256x256的数字断层模型 I phantom(modified shepp-logan, 256); imshow(I, []); title(原始模型); % 设置投影角度这里用0度到178度步长2度共90个角度 theta 0:2:178; % 计算投影信号P矩阵的行对应探测器单元列对应角度 [P, xp] radon(I, theta); % 画一下90度方向的投影信号感受一下 figure; plot(xp, P(:, find(theta 90))); xlabel(探测器位置); ylabel(投影值); title(90度方向的投影信号);xp是radon函数返回的探测器坐标轴。给radon加两个输出参数时它才会返回坐标轴很多人第一次用会漏掉这个然后自定义反投影时找不到正确的距离轴只能在代码里瞎猜这是没必要踩的坑。3.2 第二步信号傅里叶变换与斜坡滤波投影信号本质是一维信号要做频域滤波就得先做傅里叶变换。注意我这里用的是fft配合fftshift/ifftshift目的是把频率0放到序列中心这样构造斜坡滤波器时可以直接用与频率轴对应的一维数组。% 频域斜坡滤波器构造函数 n size(P, 1); freq linspace(-1, 1, n); % 归一化频率轴范围[-1,1] ramp abs(freq); % 理想斜坡滤波器 % 对每个角度的投影信号做傅里叶变换并滤波 Pf zeros(size(P)); for k 1:length(theta) proj P(:, k); PROJ fftshift(fft(ifftshift(proj))); % 一维傅里叶变换 PROJ_filt PROJ .* ramp(:); % 频域乘斜坡滤波器 Pf(:, k) fftshift(ifft(ifftshift(PROJ_filt))); % 逆变换回空间域 end % 看看滤波前后的差异 figure; plot(xp, P(:,1), b, xp, real(Pf(:,1)), r); legend(滤波前, 滤波后); xlabel(探测器位置); ylabel(幅度); title(第一个角度的投影信号滤波前后对比);这里有一个很容易被忽略的细节滤波之后得到的信号可能是复数值这是因为频域相乘后相位没有完全归零实际使用中应该取实部。如果你发现重建图像有细微的纹波先怀疑一下是不是忘了取real()。理论上也可以在空间域用卷积实现同样效果即把投影信号跟斜坡滤波器的空间域核函数做卷积。频域实现的好处是更直观、更容易控制滤波频率范围也方便切换不同窗函数。比如想用Hamming窗只需要在乘滤波器时再乘一个窗函数数组% 以Hamming窗为例 window 0.54 0.46 * cos(linspace(-pi, pi, n).); ramp_hamming abs(freq). .* window; PROJ_filt PROJ .* ramp_hamming;3.3 第三步反投影重建滤波做完接下来就是反投影。这一步的核心操作是对图像平面上的每一个像素点计算它在这个角度下对应的探测器位置然后从滤波后的投影数据中取出该位置的值累加到像素上。% 初始化重建图像 recon zeros(size(I)); % 以图像中心为原点建立坐标网格 N size(I, 1); half N / 2; [xg, yg] meshgrid(linspace(-half, half, N), linspace(-half, half, N)); % 角度步长弧度用于离散积分的权重 dtheta deg2rad(theta(2) - theta(1)); % 逐角度反投影累加 for k 1:length(theta) % 当前角度下每个像素对应的探测器坐标 t x*cos(theta) y*sin(theta) t_coord xg * cosd(theta(k)) yg * sind(theta(k)); % 用插值从滤波投影数据中取像素对应值 proj_interp interp1(xp, real(Pf(:, k)), t_coord(:), linear, 0); proj_interp reshape(proj_interp, size(xg)); % 累加 recon recon proj_interp; end % 乘上角度间隔完成离散积分近似 recon recon * dtheta; figure; imshow(recon, []); title(滤波反投影重建结果);这个循环就是整个项目最核心的代码段了。interp1里的最后一个参数0表示当像素点对应的探测器坐标超出范围时按0处理这样能减少边界外无效区域的干扰。如果你运行完这段代码发现重建图像和原始模型对比亮度和对比度差了一大截别急着改代码先把dtheta这个角度步长因子加上。反投影过程本质上是在做一个离散积分近似漏掉这个因子重建图像的整体幅度会错得离谱。3.4 第四步重建质量评估跑完重建不能只看图要量化评估。我用两个指标衡量重建效果一是均方误差(MSE)二是峰值信噪比(PSNR)。因为重建图像的幅度尺度和原图不完全一致评估前需要先把重建图缩放到和原图相同的数值范围。% 将重建图像缩放到[0,1]区间 recon_norm mat2gray(recon); I_norm mat2gray(I); % 计算MSE和PSNR mse_val mean((I_norm(:) - recon_norm(:)).^2); psnr_val 10 * log10(1 / mse_val); fprintf(MSE %.4f\n, mse_val); fprintf(PSNR %.2f dB\n, psnr_val);在我自己的测试里用256分辨率模型、90个角度、Ram-Lak滤波MSE大约在0.001量级PSNR在28dB附近。如果角度增加到180个PSNR能往上走3到4个dB。这里想提个醒PSNR在CT重建里只能作为参考指标因为它对整体灰度缩放非常敏感真正评价重建质量最好同时输出重建图和原图的差值图直观查看误差分布。4. 实操中遇到的问题排查与参数调试记录4.1 星形伪影严重怎么回事最典型的问题是重建图像上出现明显的放射状条纹从中心向四周扩散看起来像星星的光芒。这个现象的本质是角度采样不足。滤波反投影理论要求从0到180度连续采集投影实际中只能采集有限个角度。每个角度在频域中对应一条径向线角度太少频域覆盖就会出现缺口这些缺口在图像空间就表现为放射状伪影。解决办法很简单把角度步长从2度改成1度甚至0.5度伪影会明显减少。代价是计算量线性增长仿真环境下一般都能接受。另外探测器单元数太少也会造成类似问题。探测器分辨率不够相当于投影信号本身就被“糊”了一层反投影时再精细也没用。4.2 重建图像偏暗、灰度不对这个坑我印象太深了。第一次跑完反投影重建图像整体灰蒙蒙的对比度极低怎么看怎么不对劲最后发现是漏了角度步长因子。反投影的离散公式是( f_{recon}(x,y) \approx \sum_{k1}^{N} P_{\theta_k}^{filtered}(x\cos\theta_k y\sin\theta_k) \cdot \Delta\theta )如果你累加之后没有乘上dtheta等价于每个角度贡献的积分权重是1而不是dtheta重建结果就会偏大或偏小很多。另一个常见原因是斜坡滤波器的频率轴定义不对比如用了linspace(0, 1, n)而不是linspace(-1, 1, n)相当于滤波器没有负频率成分重建结果会差很多。4.3 投影角度和采样数怎么选这是做仿真时必问的问题。我的经验值是图像尺寸为N时探测器单元数大约取N*sqrt(2)因为旋转对角线会把模型的最长跨幅暴露出来投影角度数取和探测器单元数相近的值即可一般128到256个角度效果就不错了。可以自己做个简单实验把投影角度设为10、30、90、180、360五个档次分别重建比较PSNR。你会看到一个很明显的趋势——角度少时PSNR低、伪影重到180角度以后增长开始变缓360角度时伪影基本肉眼不可见。这个实验做完你对“角度分辨率”的理解会比看任何教材都深刻。4.4 重建速度太慢怎么办反投影循环是整个算法最耗时的地方主要是因为每个角度都用interp1对整张图像的每个像素做插值。图像尺寸256还好一旦上到512或102490个角度的循环会明显变慢。我的优化经验按性价比排序用parfor替代for多核并行处理代码改动最小提速效果明显把interp1换成内置的interp2高级用法或者用imrotate结合矩阵旋转但要注意插值精度预分配所有中间变量避免循环内动态增长数组如果追求极致性能可以把反投影核心写成MEX函数或使用GPU加速。对仿真教学来说做到第1步就够了。我自己的测试里256图像、180角度parfor开8核后耗时能从三四秒降到一秒以内。启动并行池本身有开销如果角度数少于90没必要开并行。5. 从平行束仿真到真正CT的扩展路径5.1 扇形束与锥形束的区别上面讲的是平行束几何也就是每束X射线都互相平行只适用于早期CT和教学仿真。真实临床CT用的是扇形束单排探测器或者锥形束多排探测器几何关系更复杂。从平行束换到扇形束反投影公式里要多一个距离加权因子因为扇形的每条射线路径长度不同采样密度也不同。多排探测器进一步扩展成锥形束之后就要用FDK算法这类近似重建方法。虽然公式复杂了但核心思想跟FBP完全一致频域补偿加反投影。理解了平行束FBP再往扇形束走主要工作是修几何权重。5.2 仿真代码还能怎么用这套仿真代码一旦跑通后续能玩的方向非常多。比如给投影信号加噪声研究不同滤波函数对噪声的抑制能力比如设置金属伪影场景在模型里插入高衰减区域模拟临床上常见的伪影问题比如把投影数据换成真实CT原始数据验证自己的重建流程是否还能成立。我自己最推荐的扩展是“用FBP做低剂量仿真”。你往投影里加泊松噪声然后对比Ram-Lak和Hamming窗的重建结果会直观看到噪声放大的代价和滤波带来的分辨率损失。这个实验做完你就明白为什么真实CT系统里滤波器设置那么讲究为什么低剂量扫描会催生一系列新算法。另外这个项目也很适合作为深度学习重建的基准。用FBP重建结果作为网络输入让深度学习模型去学习修正伪影和噪声是目前很多研究论文的标准做法。你想研究AI图像重建的话这套仿真代码就是最好的数据生成器。最后再说一点个人体会做这个项目最忌讳的是直接调库函数出个结果就完事。iradon确实一行代码就能重建但如果你没有亲手写过那个反投影循环你可能永远体会不到为什么滤波是必须的、为什么角度步长要乘上去、为什么插值方式会影响重建质量。我自己跑完这个项目之后再去翻CT重建的经典论文很多之前看不懂的公式突然变得非常具体。强烈建议你也动手拆一遍把每个环节都打印出来看一眼这种感觉是纯看书给不了的。
