CT成像仿真，采用滤波反投影算法，Matlab代码，主要包括生成投影信号，信号傅立叶变换，滤波...

CT成像仿真，采用滤波反投影算法，Matlab代码，主要包括生成投影信号，信号傅立叶变换，滤波，反变换等步骤，代码有详细注释

手把手实现CT滤波反投影仿真

CT成像的核心在于如何从投影数据重建出物体切片。滤波反投影（FBP）算法是经典方法之一，今天咱们用Matlab从零实现整个过程，顺便聊聊代码里的门道。

生成投影信号：Radon变换打底

CT投影的本质是获取物体在不同角度的线积分。Matlab自带的radon函数能快速生成仿真投影数据。这里我们构造一个包含椭圆和矩形的测试模型：

% 创建仿真模型（512x512像素） phantomSize = 512; model = phantom(phantomSize); % 设置0-179度共180个投影角度 theta = 0:179; % 生成投影数据（sinogram） [sinogram, sensorPos] = radon(model, theta);

这里phantom函数生成的标准模型包含不同密度的椭圆，radon输出的sinogram矩阵行代表传感器位置，列对应不同角度。生成的投影数据长这样：

![Sinogram示例图]

傅立叶变换：投影进频域

滤波反投影的关键步骤是对投影做傅立叶变换。这里有个细节：每条投影单独做一维FFT。

% 对每条投影做FFT（补零到2的整数次幂） nfft = 2^nextpow2(size(sinogram,1)); projection_fft = fft(sinogram, nfft);

补零操作（nextpow2）是为了加速FFT计算。此时projection_fft的每列存储的是对应角度投影的频谱。这一步相当于把投影数据从时域搬到频域，为后续滤波做准备。

滤波操作：Ram-Lak滤波器上场

频域滤波是为了抑制星状伪影。这里实现最常用的Ram-Lak（斜坡）滤波器：

% 生成斜坡滤波器 rampFilter = abs(linspace(-1,1,nfft)).'; % 频率轴修正（避免直流分量偏移） rampFilter = fftshift(rampFilter); % 频域相乘实现滤波 filtered_proj = projection_fft .* rampFilter;

滤波器的形状像个三角形，高频部分增益更大。fftshift的作用是让滤波器中心对应零频率点。注意这里是频域乘法，相当于时域的卷积操作。

反变换与反投影：空间域重建

滤波后的数据需要先做逆傅立叶变换，再通过反投影累加：

% 逆FFT取实部 filtered_proj_spatial = real(ifft(filtered_proj)); % 裁剪回原始长度 filtered_proj_spatial = filtered_proj_spatial(1:size(sinogram,1), :); % 反投影重建 reconImage = iradon(filtered_proj_spatial, theta, 'linear','Ram-Lak', 1, phantomSize);

这里iradon函数实现了反投影的积分操作。自己写反投影循环其实更直观：

% 手动反投影（核心代码段） recon = zeros(phantomSize); for i = 1:length(theta) proj = filtered_proj_spatial(:,i); recon += backproject(proj, theta(i), phantomSize); end

其中backproject函数需要实现将投影值沿角度反方向散布到图像网格的操作。这部分涉及坐标变换，代码略长但逻辑清晰。

结果对比与调参心得

运行完整代码后，重建图像与原始模型的对比：

![重建结果对比图]

几个踩坑经验：

1. 滤波时未做fftshift会导致相位错位，图像出现重影
2. 投影数据补零不足会引入高频振荡伪影
3. 反投影时的插值方法（线性 vs 最邻近）显著影响边缘清晰度

完整代码已在GitHub开源（地址示例：github.com/xxx/ct-fbp-demo），包含交互式参数调节模块。修改滤波器类型（试试Shepp-Logan）、调整投影角度数，立马能看到重建效果的变化——动手试试才是理解算法的最佳姿势！

后记

虽然现在深度学习在CT重建中很火，但掌握传统FBP算法依然重要。通过自己写代码，才能真正理解那些公式背后的物理意义。下次咱们可以聊聊怎么用GPU加速反投影——毕竟180个角度循环跑起来，Matlab还是挺吃力的...