CBCT重建算法解析：FDK与迭代重建的Matlab实现

简介：该资源为基于Matlab的3D锥形束CT（CBCT）投影与重建示例包，涵盖FDK解析重建及MLEM、SART、SQS等迭代重建方法，面向医学影像、电子信息工程、数学等专业学生，适用于课程设计、期末大作业与毕业设计等场景。压缩包共15个文件，主体为12个m脚本，涉及参数设置、测量投影生成、滤波反投影、Parker加权与各算法重构模块，另含2个说明性txt文档及1个可直接调用的mat数据文件，包体约173KB，结构精炼便于二次开发。项目采用参数化编程，注释清晰，读者可自由调整投影角度、噪声水平等参数，对比FDK与迭代算法在噪声抑制和稀疏数据重建方面的差异。目前已有148人学习，适合希望深入理解CBCT重建原理并快速上手Matlab算法验证的研究者与学生。

1. 为什么 CBCT 重建绕不开 FDK 和迭代法：从一张投影图说起

3D 锥形束 CT（CBCT）最常见的落地场景，是拿一圈锥束投影图还原出体数据。处理这类数据，手边至少要有两套算法：FDK 负责快速出图，迭代重建负责在低剂量、稀疏角度下救场。你拿到的 Matlab 示例不管叫什么名字，核心代码无非是前向投影、背投、滤波、迭代更新四个功能块，外加一组定义扫描几何的参数。适合谁：刚拿到 CBCT 仿真数据的研究生、做微焦点 CT 或口腔 CT 算法验证的工程师、把重建当成 Matlab 图像处理大作业想深入一点的本科生。它能帮你在半小时内跑通最小流程，然后有目标地替换模块，把重建过程从黑匣子变成可调试的算子组合。

2. 先把几何和算子说清：CBCT 前向投影与背投的 Matlab 最小实现

FDK 和迭代法都离不开两个基本算子：前向投影（把体数据变成投影图）和背投（把投影图加回体空间）。很多人一上来就找重建函数，结果几何定义错了，后面滤波再标准、迭代再耐心，出来的图也是花的。所以这一章先把几何参数钉死，再实现两个算子，最后用一个体模脚本验证算子对不对。

2.1 锥形束几何：从扇束到锥束，多了一个 z 方向锥角

扇束 CT 只在 z = 0 的平面里做计算，锥束不一样，射线在 z 方向也有锥角，探测器是一整面二维平板。Matlab 内置的iradon是平行束和扇束接口，没法直接吃 CBCT 的投影数据，这也是必须自己管理几何的原因。一个完整的锥束几何一般包含下面这些量，单位建议统一用毫米：

% 锥束几何参数结构体，单位统一为 mm geo.dsd = 1000; % 源到探测器距离 geo.dso = 600; % 源到旋转中心距离，决定几何放大率 geo.det_rows = 256; % 探测器行数（z 方向） geo.det_cols = 512; % 探测器列数（x-y 平面内） geo.du = 1.0; % 探测器列间距 mm geo.dv = 1.0; % 探测器行间距 mm geo.npx = 256; % 体数据三个方向的体素数 geo.pix = 1.0; % 体素尺寸 mm，各向同性 geo.n_angles = 360; % 一圈扫描的投影数量 geo.dtheta = 2 * pi / geo.n_angles; geo.angles = 0 : geo.dtheta : 2*pi - geo.dtheta; % 弧度列表

这套参数里最容易出问题的不是 DSD 和 DSO，而是旋转中心偏移。实际 CBCT 里探测器中心几乎不可能恰好落在旋转轴上，所以结构体里通常还要预留det_offset_x和det_offset_y，后面背投坐标映射时要用。现在先不加，等第五章讲坑的时候再补。

2.2 前向投影：把体数据“拍照”成投影图，核心是射线跟踪

前向投影在迭代重建里就是矩阵 A 乘以体数据 x，在 FDK 验证里则用来模拟“如果这个体模是真的，投影应该长什么样”。最常见的实现是 ray-driven：对每条射线，从源点出发，沿方向穿过体数据，等间隔采样并用三线性插值取体素值，最后累加得到线积分。下面这个单角度版本用interp3完成插值，适合小体模验证：

function proj = cbct_project_single(vol, geo, theta) % 单角度锥束前向投影（ray-driven） % vol : 3D 体数据 [nx, ny, nz]，单位是线性衰减系数 cm^-1 % geo : 几何结构体，见 2.1 % theta: 当前投影角度，弧度 % proj : 探测器图像 [det_rows, det_cols] [nx, ny, nz] = size(vol); % 源点：绕 z 轴旋转，距离旋转中心 dso sx = geo.dso * cos(theta); sy = geo.dso * sin(theta); sz = 0; % 探测器中心：位于射线源另一侧 dx0 = -(geo.dso + geo.dsd) * cos(theta); dy0 = -(geo.dso + geo.dsd) * sin(theta); dz0 = 0; % 探测器局部坐标基向量 ex = [sin(theta), -cos(theta), 0]; % 沿探测器列方向 ey = [0, 0, 1]; % 沿探测器行方向 proj = zeros(geo.det_rows, geo.det_cols); n_samples = 200; % 每条射线采样点数，速度与精度折中 for ud = 1:geo.det_cols for vd = 1:geo.det_rows % 探测器像素的世界坐标 u = (ud - geo.det_cols/2) * geo.du; v = (vd - geo.det_rows/2) * geo.dv; det_pos = [dx0, dy0, dz0] + u * ex + v * ey; % 射线方向（源 -> 探测器像素） dir = det_pos - [sx, sy, sz]; dir = dir / norm(dir); % 从源到探测器等间隔采样，端点各留一个保护点 L = norm(det_pos - [sx, sy, sz]); t_vals = linspace(0, L, n_samples + 2); t_vals = t_vals(2:end-1); pts = [sx, sy, sz] + t_vals' * dir; % 世界坐标转体素坐标：ix = x / pix + nx/2 xv = pts(:,1) / geo.pix + nx/2; yv = pts(:,2) / geo.pix + ny/2; zv = pts(:,3) / geo.pix + nz/2; % 剔除越界采样点 inside = xv > 1 & xv < nx & yv > 1 & yv < ny & zv > 1 & zv < nz; if ~any(inside) continue; end % 三线性插值并累加，再乘采样步长得到线积分 vals = interp3(vol, yv(inside), xv(inside), zv(inside), 'linear', 0); step = L / (n_samples + 1); proj(vd, ud) = sum(vals) * step; end end end

逻辑说明：interp3的查询坐标顺序是(Yq, Xq, Zq)，对应体数据的第一、二、三维索引，这里最容易写反。写反后投影图会转置九十度，而且这种错误在后面的滤波、背投阶段极难发现。inside做越界剔除很重要，因为锥束边缘射线一定会飞出体数据范围，interp3越界默认返回 NaN，不剔除的话整幅投影图全是 NaN。

参数说明：n_samples = 200对应体素尺寸 1mm、射线长度大约 1600mm 的情况，步长约 8mm，听起来很粗，但对 256^3 体模做验证足够。正式仿真可以加到 400 或 600，代价是耗时线性上涨。真实 CBCT 数据一般不用这种纯 Matlab 循环做前向投影，而是走 GPU 光线追踪或 Siddon 算法，但验证逻辑完全一样。

2.3 背投：它只是投影的转置，并不是逆

背投在解析重建里是 FDK 的累加器，在迭代重建里就是 A 的转置 AT。严格说，AT 不是把投影图“反推”回体数据，而是把投影图上的每个值沿射线方向涂回体空间，所以投影一次再背投一次，得到的是模糊、中心偏亮的图，这是正常现象，不是 bug。像素驱动背投的做法是：遍历每个体素，找到它沿射线打到探测器上的位置，再插值取投影值：

function vol = cbct_backproject_single(proj, geo, theta) % 单角度像素驱动背投 % proj : 探测器图像 [det_rows, det_cols] % vol : 当前角度累加到体空间的贡献 [nx, ny, nz] nx = geo.npx; ny = geo.npx; nz = geo.npx; vol = zeros(nx, ny, nz); sx = geo.dso * cos(theta); sy = geo.dso * sin(theta); sz = 0; dx0 = -(geo.dso + geo.dsd) * cos(theta); dy0 = -(geo.dso + geo.dsd) * sin(theta); dz0 = 0; ex = [sin(theta), -cos(theta), 0]; ey = [0, 0, 1]; % 探测器平面的法向量：源指向探测器中心的向量 n0 = [dx0 - sx, dy0 - sy, dz0 - sz]; for ix = 1:nx for iy = 1:ny for iz = 1:nz px = (ix - nx/2) * geo.pix; py = (iy - ny/2) * geo.pix; pz = (iz - nz/2) * geo.pix; % 射线方向：源 -> 当前体素 dir = [px - sx, py - sy, pz - sz]; % 射线与探测器平面的交点参数 t denom = dot(dir, n0); if abs(denom) < 1e-12 continue; end t = dot(n0, n0) / denom; hit = [sx, sy, sz] + t * dir; % 世界坐标转探测器像素索引 u = (hit - [dx0, dy0, dz0]) * ex.'; v = (hit - [dx0, dy0, dz0]) * ey.'; ud = u / geo.du + geo.det_cols/2; vd = v / geo.dv + geo.det_rows/2; if ud > 1 && ud < geo.det_cols && vd > 1 && vd < geo.det_rows vol(ix, iy, iz) = interp2(proj, ud, vd, 'linear', 0); end end end end end

逻辑说明：这条射线从源点出发经过体素后打在哪里，用的是“射线与过探测器中心、法向量为 n0 的平面求交”，比拿体素直接算角度再查表要直观。denom接近零说明体素与探测器平面几乎平行，应跳过，否则 t 会爆炸。

参数说明：这套三重循环是纯教学实现，128^3 体数据跑一个角度大约要几十秒到几分钟，视电脑而定。验证阶段建议先用 64^3 或 96^3，把几何和插值逻辑跑通，再上全分辨率。背投里没有乘任何距离权重，这个权重留给 FDK 统一加，避免两个函数各乘一次导致亮度错乱。

2.4 验证算子对不对：投影再背投的最小脚本

算子最容易出错的不是代码，而是几何方向。最常见的问题就是投影图里椭圆方向对，但重建出来旋转了 90 度，或者上下翻转。我的习惯是先用一个椭球体模做“投影 -> 背投”闭环，只看形状不看数值：

% 椭球体模，衰减系数 0.05 cm^-1 vol = zeros(128, 128, 128); [X, Y, Z] = meshgrid(linspace(-1,1,128), linspace(-1,1,128), linspace(-1,1,128)); vol((X/0.7).^2 + (Y/0.5).^2 + (Z/0.3).^2 <= 1) = 0.05; geo.dsd = 1000; geo.dso = 600; geo.det_rows = 160; geo.det_cols = 200; geo.du = 1.2; geo.dv = 1.2; geo.npx = 128; geo.pix = 1.5; proj1 = cbct_project_single(vol, geo, 0); bp = cbct_backproject_single(proj1, geo, 0); figure; subplot(1,2,1); imagesc(proj1); axis image; title('前向投影'); subplot(1,2,2); imagesc(squeeze(bp(64,:,:))); axis image; title('背投切片');

逻辑说明：这个脚本不追求还原，只验证三点：投影图里椭球位置是否居中；背投切片里椭球是否出现在同一位置；有没有明显的 NaN 或全黑列。如果投影图椭圆被拉伸或错位，先检查ex、ey基向量方向；如果背投位置偏移，先检查denom和t的符号。几何方向最容易翻车的就是这里，后面前向投影算得越精确，错起来越离谱，所以务必要做。

3. FDK 算法：扇束 FBP 到锥束的三条近似与可运行 Matlab 骨架

FDK 是 Feldman、Davis、Kress 在 1984 年提出的锥束近似重建算法，本质上就是把扇束滤波反投影推广到锥束。理解它的三个近似，才知道参数怎么调。

3.1 FDK 为什么是近似：锥角余弦加权、z 方向近似、圆形轨迹

第一个近似：对偏离中心平面的射线，用锥角的余弦做加权补偿。这个补偿在中心平面是精确的，越靠近探测器上下边缘误差越大，所以大锥角 CBCT 在 z 方向会出现模糊和伪影。第二个近似：滤波只在探测器列方向（u 方向）做一维斜坡滤波，把 z 方向当成与中心平面一样处理，这等于假设锥角很小。第三个近似：轨迹被限定为绕 z 轴的圆形，螺旋扫描、倾斜扫描都不能直接用这套公式。理解这三点之后，你就能解释为什么口腔 CBCT 用大锥角时，下颌边缘经常有一层淡淡的环形伪影，那是 FDK 的近似误差，不是代码写错。

3.2 一维斜坡滤波：先在探测器行上做，别在列上做

FDK 的处理顺序是：投影图逐像素乘余弦加权，然后沿探测器列方向（也就是角度方向对应的那一个维度）做一维滤波，最后加权背投。很多人会在这里犯迷糊：滤波到底沿哪一维？记住一句话，滤波永远沿着探测器上“对应扇角”的方向，也就是投影图的行方向，而不是 z 方向。下面给一个频域滤波器：

function proj_f = fdk_filter(proj, geo, filter_name) % FDK 一维滤波：逐行处理，沿探测器列方向 % proj : 单角度投影 [det_rows, det_cols] % filter_name: 'ramp' 或 'hann' [det_rows, det_cols] = size(proj); n = det_cols; % 归一化频率轴，范围 [-1, 1) freq = (0:n-1) - floor(n/2); freq = freq / (n/2); switch filter_name case 'ramp' filt = abs(freq); case 'hann' filt = abs(freq) .* (0.5 + 0.5 * cos(pi * freq)); otherwise filt = abs(freq); end filt = ifftshift(filt); % 对齐 fft 的输入顺序 proj_f = zeros(size(proj)); for v = 1:det_rows row = proj(v, :); F = fft(row) .* filt; proj_f(v, :) = real(ifft(F)); end end

逻辑说明：abs(freq)就是 ram-lak 滤波器，低频线性升到高频，这是 FDK 保留分辨率的核心，但也放大了高频噪声。hann在斜坡上乘了一个余弦窗，牺牲一点锐度换噪声抑制。ifftshift这步经常被漏掉，漏掉后滤波核位置错半个周期，重建图会出现交替明暗条纹，属于典型翻车现场。

参数说明：这个滤波器没有处理尺度因子，因为 FDK 的整体增益通常最后统一标定。你不需要在每一行里纠结常数对不对，等重建完拿一个已知衰减系数的体模标定一个全局比例因子即可，中间乱乘只会让调试更痛苦。

3.3 加权背投：1/L^2 权重与坐标映射

FDK 和普通背投的差别在于每个体素累加投影值时，要乘上射线长度的平方反比权重。原因是锥束射线越靠近探测器边缘，射线路径越长，单位面积内的能量被摊薄了。把 2.3 的背投函数复制一份，在累加前乘一个权重即可：

function vol = fdk_backproject_weighted(proj, geo, theta) % FDK 专用加权背投，在 cbct_backproject_single 基础上加 1/L^2 权重 % 这里省略与 2.3 相同的几何初始化代码，只列出核心差异 % ... 与前一个背投函数相同的源点、探测器中心、基向量定义 ... vol = zeros(geo.npx, geo.npx, geo.npx); for ix = 1:geo.npx for iy = 1:geo.npx for iz = 1:geo.npx % ... 同上计算 hit、ud、vd ... if ud > 1 && ud < geo.det_cols && vd > 1 && vd < geo.det_rows val = interp2(proj, ud, vd, 'linear', 0); % 重点：源到探测器像素的距离 L L = norm(hit - [sx, sy, sz]); w = (geo.dso / L)^2; vol(ix, iy, iz) = val * w; end end end end end

逻辑说明：L是源点到探测器像素的实际射线长度，geo.dso是源到旋转中心的距离。权重(dso/L)^2的来源是面积放大率，它让靠近中心的射线贡献更大、边缘射线贡献被压低。如果漏掉这个权重，重建出来的图像边缘会明显过亮，像蒙了一层反光。

参数说明：interp2双线性插值在探测器行列数足够大时误差很小；如果探测器很稀疏，建议改用interp2(..., 'spline')，但耗时明显增加。对大锥角系统，这个权重公式可以进一步修正，不过对初学阶段先跑通标准版本更重要。

3.4 完整 FDK 主循环：拼装三段流程

把余弦加权、滤波、加权背投按顺序组装起来就是完整 FDK：

function recon = fdk_recon(projs, geo, filter_name) % FDK 重建入口 % projs : [det_rows, det_cols, n_angles] % recon : [nx, ny, nz] nx = geo.npx; ny = geo.npx; nz = geo.npx; recon = zeros(nx, ny, nz); for ia = 1:geo.n_angles theta = geo.angles(ia); % 第一步：锥角余弦加权 proj = projs(:, :, ia); [det_rows, det_cols] = size(proj); u = ((1:det_cols) - det_cols/2) * geo.du; v = ((1:det_rows) - det_rows/2) * geo.dv; cos_w = geo.dso ./ sqrt(geo.dso^2 + u.^2 + v.^2); proj = proj .* cos_w; % 第二步：沿行一维滤波 proj = fdk_filter(proj, geo, filter_name); % 第三步：加权背投累加 recon = recon + fdk_backproject_weighted(proj, geo, theta); end % 整体尺度归一化：让体模的已知衰减系数对得上 recon = recon / max(recon(:)) * 0.05; % 示例性标定，按实际数据改 end

逻辑说明：cos_w的公式里sqrt(dso^2 + u^2 + v^2)表示源到该探测器像素的距离，dso除以它是典型 FDK 余弦加权写法。注意u、v都是向量，Matlab 的广播机制会把它们扩张成整个探测器平面的权重。

参数说明：最后那句归一化占了很大便宜。FDK 的绝对数值受探测器尺寸、滤波尺度、角度步长影响，不同示例包给出的常数可能差几倍，与其对着论文逐项校对，不如用一个衰减系数已知的仿真体模把比例因子标定出来。等这套流程稳定了，再回头细抠尺度公式。

4. 迭代重建：从 ART 到 SART 的 Matlab 实现与收敛加速

FDK 对噪声和稀疏角度很敏感。CT 剂量越低、投影数越少，FDK 的重建质量下降得越快，这时迭代重建就顶上来了。

4.1 迭代重建的本质：解一个巨大的线性方程组

把体数据拉成一列向量 x，把全部投影值拉成一列向量 b，前向投影就等价于一个巨型线性方程组 Ax = b。A 的每一行是一条射线穿过体素的加权和，直接求逆不可能，因为 A 的维度轻松超过千万乘百万。迭代法不走求逆，而是从初始估计出发，反复做“前向投影 -> 求残差 -> 背投加回”。这个循环里，前向投影对应 A，背投对应 AT，两个算子你已经有了。

迭代重建最值钱的地方是可以加先验：稀疏角度、有限视角、低剂量噪声，都可以通过正则项或约束塞进循环。代价是慢，通常比 FDK 慢一到两个数量级。

4.2 ART：逐射线更新，简单但慢

ART 是逐射线更新的：每取一条射线，算当前估计在这条射线上的投影值，与实际值比出残差，沿这条射线把残差加回沿线体素。理解它的核心公式远比跑起来重要，因为纯 Matlab 里按射线循环 512×512×360 次，等到天荒地老也跑不完：

function x = art_update_ray(x, proj_val, geo, theta, u, v, lambda) % 单射线 ART 更新 % p_est = <a_row, x>，即这条射线穿过体素的加权和 p_est = ray_trace_forward(x, geo, theta, u, v); err = proj_val - p_est; % 背投回加：x <- x + lambda * err * a_row' / <a_row, a_row> x = ray_trace_back_add(x, err, geo, theta, u, v, lambda); end

逻辑说明：ray_trace_forward和ray_trace_back_add本质上就是第 2 章里前向投影和背投的单条射线版本，只是在更新时不需要遍历整个探测器，只处理当前这条射线穿过的体素。严格 ART 的松弛因子 lambda 通常取 0.05 到 0.3，取大了会震荡甚至发散。

参数说明：ART 的优点是原理透明，但逐射线更新导致收敛很慢，而且对射线遍历顺序很敏感，不同扫描顺序会带来条带伪影。所以实际项目里很少直接上朴素 ART，更多用它的变体 SART 或 OS-SART。这篇示例保留 ART 是为了让你理解迭代的底层机制，跑实验建议用下一节的 SART。

4.3 SART：按角度子集更新，收敛更稳

SART 把一个角度内的全部射线一起处理：先做一次完整前向投影，得到整幅估计投影图，减出残差图，整个背投回体空间再更新。这样噪声和条带被平均掉了，收敛稳定得多。最小可运行的 SART 可以直接复用第 2 章的两个函数：

function x = sart_recon(projs, geo, n_iter, lambda) % SART 最小实现：每个角度一个子集 nx = geo.npx; ny = geo.npx; nz = geo.npx; x = zeros(nx, ny, nz); % 权重图：用全 1 投影背投一次，得到每个体素被射线覆盖的总长度 % 这套几何固定，权重图可以只算一次 weight_map = zeros(nx, ny, nz); for ia = 1:geo.n_angles ones_proj = ones(geo.det_rows, geo.det_cols); weight_map = weight_map + cbct_backproject_single(ones_proj, geo, geo.angles(ia)); end weight_map(weight_map < 1e-6) = 1e-6; % 避免除零 for iter = 1:n_iter for ia = 1:geo.n_angles theta = geo.angles(ia); p_est = cbct_project_single(x, geo, theta); err = projs(:, :, ia) - p_est; delta = cbct_backproject_single(err, geo, theta); x = x + lambda * delta ./ weight_map; end end end

逻辑说明：weight_map是 SART 里关键的一步。它统计每个体素在所有投影中“被多少射线覆盖”，更新时除以它，相当于把残差按覆盖度平均分配，避免中心体素被过度更新、边缘体素欠更新。很多简化的 SART 示例直接除以一个常数或跳过这步，在 360 度满采样下问题不大，一旦用稀疏角度，就会看到中心亮、边缘暗的锅盖状伪影。

参数说明：n_iter对这个最小实现建议取 10 到 30 轮，超过 50 轮收益很小，还容易陷入噪声。lambda取 0.5 到 1.5 都算合理范围，建议从 1.0 开始，观测投影误差曲线再微调。注意这里的projs是原始投影图，如果源数据是强度比 I/I0，必须先在入口取负对数，否则重建出来的数值完全不可用。

4.4 超参数：迭代次数、松弛因子、正则项

迭代重建不是魔法，成败往往在三个超参数上。第一个是迭代次数：数据干净、投影角多时 10 轮就够；稀疏角度要加到 30 轮以上。第二个是松弛因子：SART 的 lambda 从 1.0 起步，ART 从 0.1 起步，出现震荡就减半。第三个是正则项：最常见的是在每次迭代后对体数据做一次总变差去噪，等价于让解更平滑，这在低剂量场景效果明显。

我一般会在 SART 里预留一个回调位置，每轮算出投影误差后决定是否继续。这样 30 轮和 10 轮之间有了依据，而不是凭感觉拍脑袋。后面第六章会给具体的误差曲线做法。

5. 重建避坑手册：几何失配、滤波过冲、内存与收敛判据

这一章是踩坑记录，每一条都是真实项目里反复出现的问题，按“现象 -> 原因 -> 解决”来写。

5.1 几何参数对不上：重影与杯状伪影

现象：FDK 重建出来的椭球边界有重影，边缘发虚，旋转中心附近出现杯状凹陷；换一组角度后伪影位置跟着转。

原因：最常见的是 DSO 和 DSD 的比值不对，或者旋转中心偏移量没有纳入几何。另一个隐蔽原因是探测器的ex、ey基向量方向写反，导致坐标映射里差一个负号。前向投影和背投用的若是一套几何，投影再背投还能对上，但重建时所有角度叠加就彻底露馅。

解决：先用单颗钢珠或仿真球体做标定。每一角度下，钢珠的投影位置应当是一个闭合椭圆或圆，如果位置跳变，说明角度步长或源点公式有问题；如果位置整体偏移，就在几何里加det_offset_x，并同步修改前向投影和背投的坐标映射。这个步骤不能省，几何不对，后面所有模块都白调。

5.2 滤波函数选错：边缘过冲与噪声放大

现象：FDK 结果里骨骼边缘有一圈黑边或白环，软组织变成颗粒状，像砂纸磨过。

原因：ram-lak 滤波器在每个高频点都线性放大，投影数据里只要有少量噪声或插值误差，就会被放大成条纹。如果探测器有坏像素，还会出现以坏像素为中心的星状伪影。

解决：先确认你用的到底是ramp还是hann。想看分辨率就上 ramp，噪声压不住就换 hann。另一种常用做法是给 ramp 加一个截止频率，把接近 Nyquist 频率的部分直接置零，这在 CBCT 里是常规操作。注意改正坏像素要在滤波之前完成，否则坏点被斜坡放大污染整行。

5.3 迭代发散：松弛因子与收敛判据

现象：SART 前 5 轮图像越来越清晰，到第 15 轮开始出现椒盐状噪声，数值范围越来越大，最后图像像雪花屏。

原因：松弛因子太大，或者根本没有收敛判据，迭代在追投影里的噪声。投影误差一开始下降，随后开始反弹，这个拐点就是过拟合的开始。

解决：固定一个验证角度集，每轮迭代后算一次投影误差，误差连续三轮上升就停。SART 的 lambda 从 1.0 开始，如果前几轮就震荡就把 lambda 减半。ART 同理，但起始值要降到 0.1。把迭代当作一个可控过程，而不是“越多越好”，这条血泪经验能帮你省掉大量无效调试。

5.4 内存爆炸：512^3 体数据与系统矩阵

现象：迭代重建跑到一半 Matlab 报内存不足，或者系统一直占内存直到卡死。

原因：最常见的是试图显式构造系统矩阵 A。512 列探测器、360 个角度，射线条数是亿级，再乘以 512^3 列体素，这个矩阵在物理上不可能存下。即便只用 256^3，稠密矩阵也会立刻吃光内存。

解决：永远不要显式构造 A，用算子代替矩阵。第 2 章的前向投影和背投就是算子，SART 只需要它们，不需要 A 本身。另一个缓解手段是全程用single类型，投影图和体数据的体积直接砍半。到 512^3 这个量级，纯 Matlab 循环已经不适合跑迭代，常见做法是用 ASTRA 这类 GPU 加速库做前向投影和背投，Matlab 只负责编排数据流。

5.5 灰度单位不一致：线性衰减还是强度比

现象：同一套代码，换一组投影数据后重建数值差几个数量级；FDK 和 SART 结果对不上，一个数值在 0.01 量级，另一个在 100 量级。

原因：CBCT 原始探测器数据往往是强度比 I/I0，要取-log才是线积分；而仿真体模直接给线性衰减系数，投影值就是线积分。把强度比当成线积分用，或者反过来重复取了对数，数值当然对不上。

解决：在数据入口统一做一次转换，注释写清楚当前这条路径用的是哪种约定。我习惯在结构体里加一个proj_unit字段，代码里显示判断，宁可多写一行proj = -log(proj ./ I0)，也不要让换算散落在各个函数里。这个黑白匣子问题，排查起来相当耗时间。

6. 验证与进阶：用投影误差曲线和低分辨率预跑把重建调到可信

最后分享两个每次重建前必做的习惯，它们能帮你在正式跑全分辨率之前发现大部分问题。

第一个是低分辨率预跑。先别用 512^3，把体数据降到 128^3 或 96^3，投影角度数降到 90 到 120，跑通整个 pipeline。这一步只验证代码路径和参数范围，通常几分钟内能出结果，发现几何错误再调整成本极低。等低分辨率结果形状正确、误差曲线正常，再上全分辨率，节省的是半小时起步的等待时间。

第二个是投影误差曲线。每轮迭代后对固定角度子集做一次前向投影，算归一化均方根误差 NRMSE：

err_iter(iter) = norm(p_est(:) - p_measured(:)) / norm(p_measured(:));

然后画出随迭代次数的曲线。规范的过程是这条曲线单调下降后趋平，如果在某个点开始回升，即使重建图看起来还在变干净，也要停在那里。FDK 和 SART 的对比也建议用这条曲线说话，而不是只贴两张图。

下面这张表是我常用的选择口径：

我现在拿到一套 CBCT 数据，不管是不是教学示例，都会先跑一遍 FDK 验几何，再把 FDK 的初值喂给 SART，全程盯着投影误差曲线决定停在哪。这个流程帮我避开了无数次“重建图很好看但临床数值完全不可信”的翻车。希望你也能从这套最小实现出发，把 CBCT 重建从黑匣子变成自己手里可控的模块。希望帮到你。

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