简介：这份资源是面向计算机、电子信息工程、数学等专业学生的锥形束CT（CBCT）图像重建Matlab示例包，聚焦投影背投、FDK解析重建与迭代重建两类核心算法，可用于课程设计、期末大作业与毕业设计，也适合作为医学影像算法入门与原型验证的实践素材。压缩包共15个文件，以12个m脚本为主体，另含2个txt说明与1个mat案例数据，整体约173KB，脚本按参数设置、投影生成、滤波、加权、FDK重建及MLEM、SART、SQS等迭代重建分模块组织，便于对照理解算法流程。资源采用参数化编程，注释清晰，配套案例数据可直接运行，读者能借此掌握从投影数据到三维体积数据的完整重建链路，并对比解析法与迭代法在噪声和数据不完整情形下的表现差异。目前已有148人学习下载，适合需要快速上手CBCT重建实验的学习者。

1. 从一张投影图到三维体数据：这套 FDK 与迭代重建 Matlab 示例到底能跑出什么

手里只有一圈二维投影图，却要还原出物体内部的三维密度分布，这是锥形束 CT 重建最核心的诉求。这套 Matlab 示例把两条主流路线都摆了出来：一条是解析法里的 FDK 投影背投，另一条是迭代重建。前者快、结构清晰，适合先把几何关系和滤波背投链路跑通；后者慢，但在投影数少、噪声大、金属伪影明显的场景里往往能救回图像质量。它适合正在做 CBCT 算法验证、想拿现成脚本改参数看效果的人，也适合刚接触三维重建、需要一份能直接运行、能逐行打断点看中间变量的工程参考。整套代码围绕锥束几何、滤波函数、体素遍历和迭代更新展开，不是纯理论推导，而是能落到矩阵运算和循环里的实现。

2. FDK 解析重建：滤波、锥束加权与背投的完整链路

2.1 为什么先跑 FDK 而不是直接上迭代

FDK 是锥形束圆轨迹重建里最经典的解析方法，本质是把三维锥束问题拆成一系列二维扇束问题，再沿探测器行方向做加权和滤波，最后把滤波后的投影按锥束几何背投回三维体素。它的优势是计算量可控，一次重建基本就是“加权—滤波—背投”三步，不需要反复迭代，适合先验证投影数据、几何参数和体素网格是否对齐。如果 FDK 都跑不出结构，迭代重建大概率也救不回来，因为几何或数据本身就有问题。

常见做法是先用小尺寸体数据跑通 FDK，确认重建切片里能看到物体轮廓，再换迭代。示例里一般会提供投影生成、滤波和背投三个独立函数，方便你单独替换某一环做对比。

2.2 投影数据的组织与几何参数

CBCT 投影数据通常组织成三维矩阵：探测器行、探测器列、投影角度。假设投影数为N_proj，探测器尺寸为N_u × N_v，数据就是N_v × N_u × N_proj。几何参数包括源到旋转中心距离、源到探测器距离、探测器像素尺寸、旋转角度步长。这些参数一旦有一个对不上，重建结果就会模糊、重影甚至完全错位。

下面是一段常见的投影数据读取与几何参数初始化代码，用来把数据整理成后续滤波和背投能直接用的形式：

% 读取投影数据，假设为单精度浮点，尺寸 [N_v, N_u, N_proj] proj = load('projection_data.mat'); proj = proj.data; % 几何参数（按实际扫描设置修改） geo.SDD = 1000; % 源到探测器距离，单位 mm geo.SOD = 600; % 源到旋转中心距离，单位 mm geo.du = 0.8; % 探测器像素尺寸 u 方向，mm geo.dv = 0.8; % 探测器像素尺寸 v 方向，mm geo.angles = linspace(0, 2*pi, size(proj,3)+1); geo.angles(end) = []; % 去掉重复的 2*pi % 体素网格 vol.Nx = 256; vol.Ny = 256; vol.Nz = 256; vol.dx = 1.0; vol.dy = 1.0; vol.dz = 1.0; % 体素尺寸 mm

逻辑说明：SDD和SOD决定锥束张角，du、dv决定探测器采样密度，angles决定背投时每个投影对应的旋转位置。体素网格的Nx/Ny/Nz和dx/dy/dz要和重建目标匹配，太大跑得慢，太小看不清细节。参数改错最典型的表现是重建图像整体放大或缩小，或者边缘出现明显错层。

2.3 滤波函数的选择与实现

FDK 里的滤波通常是斜坡滤波加窗，常见的有 Ram-Lak、Shepp-Logan、Cosine、Hamming 等。Ram-Lak 分辨率最高但噪声也最明显，Shepp-Logan 和 Cosine 会压高频、降噪但牺牲一点细节。示例里一般会提供一个filter_projection函数，对每个角度的投影沿 u 方向做一维滤波。

function proj_f = filter_projection(proj, du, filter_type) % 对每个角度的投影沿探测器 u 方向做斜坡滤波 [N_v, N_u, N_proj] = size(proj); proj_f = zeros(size(proj), 'like', proj); % 构造频域滤波窗 n = 2^nextpow2(2*N_u); f = linspace(-1/(2*du), 1/(2*du), n); ramp = abs(f); switch lower(filter_type) case 'ram-lak' window = ones(size(ramp)); case 'shepp-logan' window = sinc(f / (2*max(f))); case 'cosine' window = cos(pi * f / (2*max(f))); otherwise error('未知滤波类型'); end H = ramp .* window; H = ifftshift(H); for k = 1:N_proj for j = 1:N_v row = squeeze(proj(j,:,k)); row_f = fft(row, n); row_f = row_f .* H(:); row_b = real(ifft(row_f, n)); proj_f(j,:,k) = row_b(1:N_u); end end end

逻辑说明：先把投影沿 u 方向补零到n点，转到频域乘上斜坡窗，再反变换回空域。filter_type控制窗函数，du影响频率轴刻度。参数改错时，滤波太强会让图像发糊，滤波太弱会留下明显条纹。注意补零长度要足够，否则会出现环绕伪影。

2.4 锥束加权与背投

滤波之后要做锥束加权，通常按SOD / sqrt(SDD^2 + u^2 + v^2)之类的因子对投影做逐像素修正，然后再背投。背投就是遍历每个体素，计算它在每个角度下投影到探测器的位置，做双线性插值累加。

% 锥束加权 [u_grid, v_grid] = meshgrid( ... ((1:N_u) - (N_u+1)/2) * du, ... ((1:N_v) - (N_v+1)/2) * dv); weight = geo.SOD ./ sqrt(geo.SDD^2 + u_grid.^2 + v_grid.^2); weight = reshape(weight, [N_v, N_u, 1]); proj_w = proj_f .* weight; % 背投（简化版，逐体素遍历） vol_rec = zeros(vol.Nx, vol.Ny, vol.Nz, 'single'); for iz = 1:vol.Nz z = (iz - (vol.Nz+1)/2) * vol.dz; for iy = 1:vol.Ny y = (iy - (vol.Ny+1)/2) * vol.dy; for ix = 1:vol.Nx x = (ix - (vol.Nx+1)/2) * vol.dx; acc = 0; for k = 1:N_proj theta = geo.angles(k); % 旋转到探测器坐标系 u = x*cos(theta) + y*sin(theta); v = z; % 投影到探测器像素 u_idx = u * geo.SDD / (geo.SOD + x*sin(theta) - y*cos(theta)) / geo.du + (N_u+1)/2; v_idx = v * geo.SDD / (geo.SOD + x*sin(theta) - y*cos(theta)) / geo.dv + (N_v+1)/2; if u_idx >= 1 && u_idx <= N_u && v_idx >= 1 && v_idx <= N_v acc = acc + interp2(proj_w(:,:,k), u_idx, v_idx, 'linear', 0); end end vol_rec(ix,iy,iz) = acc; end end end

逻辑说明：加权是为了补偿锥束路径长度差异，背投是把每个体素在所有角度下的贡献累加。u_idx、v_idx是体素投影到探测器的连续坐标，interp2做双线性插值。循环顺序可以调整，但要注意内存访问模式，Matlab 里按列优先遍历会快一些。这段代码跑小体数据可以，大体数据建议向量化或分块。

3. 迭代重建：从 ART 到 SIRT 的落地写法

3.1 迭代重建适合什么场景

当投影角度少、噪声大、或者有金属植入物时，FDK 往往会出现明显条纹和伪影。迭代重建把重建问题当成一个线性方程组Ax = b来解，A是系统矩阵，x是待重建体素，b是投影数据。常见的有 ART、SIRT、SART 等。ART 每次只用一个投影更新，收敛快但不稳定；SIRT 用所有投影的残差做平均更新，稳定但慢；SART 折中，按角度分块更新。

示例里一般会提供至少一种迭代方法的实现，方便和 FDK 做对比。迭代重建的代价是计算量大，一次重建可能要几十到几百次迭代，每次迭代都要做一次前投影和一次背投影。

3.2 SIRT 的更新公式与实现

SIRT 的更新可以写成：x^(n+1) = x^n + C A^T R (b - A x^n)，其中C和R是归一化矩阵，通常取每行或每列和倒数。实现时前投影和背投影可以用和 FDK 类似的射线驱动或体素驱动方式。

% SIRT 迭代重建（简化版） x = zeros(vol.Nx, vol.Ny, vol.Nz, 'single'); num_iter = 50; relax = 0.5; % 松弛因子 % 预计算归一化（这里用简化常数，实际应按射线长度累加） C = 1 / N_proj; R = 1 / (vol.Nx * vol.Ny * vol.Nz); for iter = 1:num_iter % 前投影：Ax Ax = forward_project(x, geo, vol); % 残差 residual = proj - Ax; % 背投影：A^T residual back = back_project(residual, geo, vol); % 更新 x = x + relax * C * R * back; % 可选：非负约束 x(x < 0) = 0; fprintf('迭代 %d，残差范数 %.4f\n', iter, norm(residual(:))); end

逻辑说明：forward_project和back_project需要和 FDK 里的几何一致，否则迭代不收敛。relax控制步长，太大震荡，太小收敛慢。非负约束在 CT 里很常用，因为密度不能为负。残差范数可以用来判断是否收敛，但不要只看它，还要看图像视觉质量。

3.3 迭代次数与正则化的取舍

迭代次数不是越多越好。次数太少，图像欠重建；次数太多，噪声和伪影会被放大。常见做法是跑 20 到 100 次，中途保存几个中间结果对比。如果噪声明显，可以加 TV 正则化或早停。示例里如果没带正则化，可以自己加一个简单的平滑项，但要注意别把细节也抹掉。

提示：迭代重建跑之前先用 FDK 确认几何参数正确，否则迭代只会把错误几何放大成更奇怪的伪影。

4. 避坑与排查：几何、滤波和内存这三处最容易翻车

4.1 重建图像整体错位或重影

现象：重建切片里物体位置偏移，或者出现双影。原因：旋转中心、源到探测器距离、角度步长有一个对不上。解决：先用一个已知位置的球体或点源做仿真，检查投影上的位置和重建位置是否一致。角度步长如果是2*pi/N_proj，注意不要多算一个重复角度。

4.2 滤波后图像发糊或条纹明显

现象：FDK 重建结果模糊，或者有明显条纹。原因：滤波窗选择不当，或者补零长度不够导致环绕。解决：换 Ram-Lak 对比，确认是不是窗函数问题；把补零长度加到2*N_u以上；检查du是否和探测器实际像素尺寸一致。

4.3 迭代重建不收敛或震荡

现象：残差范数不降反升，或者图像出现棋盘格。原因：前投影和背投影的几何不一致，或者松弛因子太大。解决：先用一个简单仿真数据验证forward_project和back_project是否互为转置；把relax降到 0.1 到 0.3 再试；检查归一化矩阵是否合理。

4.4 内存爆掉或跑得太慢

现象：Matlab 报内存不足，或者循环跑几个小时。原因：体素网格太大，或者背投用了逐体素四重循环。解决：先把Nx/Ny/Nz降到 128 或 64 跑通；把背投向量化，或者用parfor按切片并行；投影数据用single而不是double。

4.5 投影数据和体素网格单位不统一

现象：重建结果尺寸明显不对。原因：du、dv、dx单位混用，比如一个用 mm 一个用 cm。解决：把所有几何参数统一到同一单位，通常用 mm；检查SOD和SDD的比例是否合理。

5. 进阶技巧：用仿真数据验证重建链路，再换真实投影

5.1 先造一个已知体模，跑通全链路

拿真实投影直接调重建，出了问题很难判断是数据问题还是代码问题。我一般会先用 Shepp-Logan 体模或几个球体造仿真投影，跑一遍 FDK 和 SIRT，看重建结果和原始体模的差异。这样几何、滤波、背投、迭代每一环都能单独验证。

% 造一个简单体模：中心球 + 两个偏心球 vol_true = zeros(vol.Nx, vol.Ny, vol.Nz, 'single'); [xx, yy, zz] = ndgrid( ... (1:vol.Nx) - (vol.Nx+1)/2, ... (1:vol.Ny) - (vol.Ny+1)/2, ... (1:vol.Nz) - (vol.Nz+1)/2); vol_true(sqrt(xx.^2 + yy.^2 + zz.^2) < 20) = 1.0; vol_true(sqrt((xx-30).^2 + yy.^2 + zz.^2) < 10) = 0.8; vol_true(sqrt((xx+30).^2 + (yy-20).^2 + zz.^2) < 8) = 0.6; % 前投影生成仿真数据 proj_sim = forward_project(vol_true, geo, vol); % 用仿真数据跑 FDK 和 SIRT，对比 vol_fdk = fdk_reconstruct(proj_sim, geo, vol); vol_sirt = sirt_reconstruct(proj_sim, geo, vol, 50, 0.3); % 计算均方误差 mse_fdk = mean((vol_fdk(:) - vol_true(:)).^2); mse_sirt = mean((vol_sirt(:) - vol_true(:)).^2); fprintf('FDK MSE: %.6f, SIRT MSE: %.6f\n', mse_fdk, mse_sirt);

逻辑说明：vol_true是已知真值，forward_project生成仿真投影，再分别用 FDK 和 SIRT 重建，最后算 MSE。MSE 不是唯一指标，但能快速判断链路是否跑通。如果 FDK 的 MSE 比 SIRT 还大很多，说明几何或滤波有问题。

5.2 换真实投影时的检查清单

仿真跑通后换真实投影，按这个顺序检查：投影数据是否做了暗场和亮场校正；角度范围是否覆盖 360 度；探测器像素尺寸和几何参数是否来自同一份标定文件；投影数据是否需要取负对数。常见做法是先重建一个中心切片，看有没有明显伪影，再决定要不要调滤波或加迭代。

5.3 参数对比表

5.4 我踩过的一个坑

有一次换了一套投影数据，FDK 重建出来全是条纹，查了半天以为是滤波问题，最后发现是角度步长用了2*pi/(N_proj-1)，而数据实际是2*pi/N_proj。从那以后我每次换数据都强制先跑一遍仿真体模，确认几何参数和角度定义一致，再动真实投影。希望帮到你。

本文还有配套的精品资源，点击获取