简介：这份资源是面向计算机、电子信息工程、数学等专业学生的锥形束CT图像重建Matlab实践包，聚焦投影背投、FDK解析重建与迭代重建两类核心算法，可用于课程设计、期末大作业与毕业设计。压缩包共15个文件，以12个m脚本为主，另含2个txt说明与1个mat案例数据，整体约173KB，脚本覆盖参数配置、投影生成、滤波、加权、背投及FDK、MLEM、SART、SQS等多种重建实现，数据可直接运行验证。代码采用参数化编程，参数便于修改，思路清晰且注释详细，兼容Matlab 2014、2019a与2024a。已有148人学习，读者可借此对比解析重建与迭代重建在噪声和数据不完整场景下的表现，快速上手算法验证与原型开发。1. 从一圈投影到三维体数据：FDK 与迭代重建到底在算什么手里只有一圈 X 光投影图，却想还原出物体内部的三维密度分布，这就是锥形束 CT（CBCT）重建要解决的核心问题。工业无损检测、口腔颌面成像、小型样品三维扫描，很多场景用的都是锥束几何：射线源和面阵探测器绕物体转一圈，采到几百上千张投影，然后靠算法把这堆二维投影“背投”回三维空间。标题里的 FDK 就是这类圆轨迹锥束重建最经典的解析算法，而迭代重建则是当投影稀疏、噪声大、金属伪影重时更稳的另一条路。Matlab 示例的价值在于：它把投影生成、滤波、背投、迭代更新这几步拆开，让你能改参数、看中间结果，而不是调一个黑匣子函数了事。这篇笔记面向已经懂点 CT 原理、想自己动手复现的工程师和学生，从几何参数一路讲到避坑。2. FDK 解析重建：从滤波到锥束背投的完整链路2.1 为什么圆轨迹锥束能用 FDK 近似FDK 的本质是把三维锥束问题拆成一层层二维扇束问题来近似处理。它的推导前提是圆轨迹、射线源到旋转中心距离固定、探测器平面与旋转轴平行。严格来说圆轨迹不满足 Tuy 数据充分性条件，靠近上下两端的大锥角层面会有明显误差，所以 FDK 在锥角小的时候才准。工程上判断能不能用 FDK，先看半锥角：一般控制在几度以内，重建质量还能接受；锥角一大，远端层面就会出现亮度不均和模糊。这也是为什么很多工业 CBCT 会把探测器做小、把源到探测器距离拉长，就是为了压锥角。理解这一点很关键：FDK 不是“精确重建”，而是“够用的近似”。你如果拿它去重建大锥角的平板探测器数据，然后抱怨边缘层糊，那不是代码写错了，是算法前提不满足。常见做法是先做一次 FDK 看整体，再决定要不要上迭代或者做锥角校正。2.2 投影生成：用 Matlab 造一份可复现的仿真数据要验证重建代码，先得有可控的投影数据。最稳的办法是自己造一个三维体模（phantom），前向投影得到正弦图，这样真值已知，重建对错一眼能看出来。下面这段代码生成一个含球和圆柱的简单体模，并按圆轨迹锥束几何做前向投影。% 生成三维体模并做锥束前向投影（简化版，便于理解流程） N = 128; % 体数据边长 vol = zeros(N, N, N); % 初始化体模 [X, Y, Z] = ndgrid(1:N, 1:N, 1:N); % 中心放一个球 vol((X-N/2).^2 + (Y-N/2).^2 + (Z-N/2).^2 < (N/6)^2) = 1; % 旁边放一个圆柱，用来观察边缘和伪影 vol((X-N/2-20).^2 + (Y-N/2).^2 < (N/10)^2 & abs(Z-N/2) < N/4) = 0.6; nAng = 360; % 投影角度数 ang = linspace(0, 2*pi, nAng+1); ang(end) = []; SDD = 400; % 源到探测器距离 SOD = 200; % 源到旋转中心距离 detN = 160; % 探测器像素数（面阵 detN x detN） proj = zeros(detN, detN, nAng); for k = 1:nAng % 把体模旋转到当前角度，等价于源探绕物体转 R = [cos(ang(k)) -sin(ang(k)) 0; sin(ang(k)) cos(ang(k)) 0; 0 0 1]; % 这里用简单射线驱动求和代替精确线积分，重在流程演示 proj(:,:,k) = forwardProjectSimple(vol, R, SDD, SOD, detN); end逻辑说明：vol是待重建的真值，ang是采样角度，SDD/SOD决定几何放大比。forwardProjectSimple是自定义函数，沿每条射线累加体素值。参数上，nAng越大角度采样越密，重建越稳但越慢；detN决定探测器分辨率，太小会丢细节，太大则内存吃紧。实际做仿真时，射线驱动求和精度不够会引入误差，建议用 Siddon 类精确线积分或者 Matlab 自带的radon做二维验证，再扩展到三维。2.3 滤波：斜坡滤波器为什么必须加窗FDK 的第一步是对每张投影做滤波，核心是斜坡滤波器（Ram-Lak），它补偿背投带来的低频堆积。但纯斜坡滤波器会把高频噪声放大，所以工程上几乎都要加窗，常见的有 Shepp-Logan、Hamming、Hann。窗越强，噪声越小，但空间分辨率也越低，这是一对硬矛盾。function pf = rampFilter(proj, detSpacing, windowType) % proj: 单张投影 (detN x detN)，对每行做一维滤波 [nRow, nCol] = size(proj); pf = zeros(nRow, nCol); % 构造频域滤波器 f = (-nCol/2:nCol/2-1) / (nCol * detSpacing); ramp = abs(f); switch windowType case 'ram-lak' w = ones(size(f)); case 'shepp-logan' w = sinc(f / (2*max(abs(f)))); case 'hamming' w = 0.54 + 0.46 * cos(pi * f / max(abs(f))); otherwise w = ones(size(f)); end H = ramp .* w; H = ifftshift(H); for i = 1:nRow row = proj(i, :); pf(i, :) = real(ifft(fft(row) .* H)); end end逻辑说明：滤波器在频域构造，ramp是斜坡项，w是窗函数。detSpacing是探测器像素间距，它决定频率轴刻度，填错会让滤波器尺度失真，重建出来要么过锐要么过糊。参数选择上，噪声大的数据用 Hamming 或 Hann，追求分辨率用 Ram-Lak 或 Shepp-Logan。注意滤波要沿探测器水平方向（对应扇束的径向），别滤错方向，否则重建直接翻车。2.4 锥束背投：体素驱动还是射线驱动背投是把滤波后的投影按几何关系累加回三维体素。两种实现思路：体素驱动，遍历每个体素找它投影到探测器哪个位置；射线驱动，遍历每条射线累加路径上的体素。体素驱动更适合并行、内存访问规律，射线驱动精度高但慢。FDK 里常用体素驱动加距离加权。% 体素驱动背投核心片段（简化） recon = zeros(N, N, N); for k = 1:nAng P = pf(:,:,k); % 当前角度滤波后投影 for iz = 1:N for iy = 1:N for ix = 1:N % 把体素坐标旋转到当前角度 p = R{k} * ([ix; iy; iz] - [N/2; N/2; N/2]); % 锥束投影到探测器 u = SDD * p(1) / (SOD - p(2)); v = SDD * p(3) / (SOD - p(2)); % 转成探测器索引并双线性插值 ui = u / detSpacing + detN/2 + 0.5; vi = v / detSpacing + detN/2 + 0.5; if ui>1 && ui<detN && vi>1 && vi<detN recon(ix,iy,iz) = recon(ix,iy,iz) + ... interp2(P, ui, vi, 'linear', 0) * (SOD/(SOD-p(2)))^2; end end end end end逻辑说明：(SOD/(SOD-p(2)))^2是锥束背投的距离加权因子，少了它重建亮度会随位置漂移。interp2做双线性插值，最近邻插值会让重建出现块状伪影。参数上，detSpacing必须和投影生成时一致，旋转矩阵R{k}的方向要和前向投影严格对应，差一个转置就会得到镜像或错位的结果。这段三重循环在 Matlab 里很慢，实际会用permute、向量化或者gpuArray加速。3. 迭代重建：当 FDK 不够用时怎么把它按在地上摩擦3.1 迭代重建的建模思路与代价迭代重建把重建当成一个优化问题：找一组体素值，使它的前向投影和实测投影最接近。最基础的是 ART/SIRT，进阶的有 SART、OS-SART，再往上还有带正则项的最大似然类方法。它的好处是能在投影少、噪声大、有金属的情况下通过先验约束压住伪影，代价是计算量比 FDK 大一到两个数量级，而且参数没调好会越迭代越糟。选型上，如果投影角度足够密（比如一圈几百张）、锥角小、噪声可控，FDK 又快又稳，没必要上迭代。反过来，稀疏视角（几十张）、低剂量高噪声、或者有金属植入物，迭代重建的优势才明显。别为了“高级”而用迭代，先问数据条件够不够。3.2 SIRT 的一次迭代：前向投影加残差回投SIRT 每次迭代做两件事：用当前估计做前向投影，算残差，再把残差按几何回投修正体数据。下面给出单步核心。% SIRT 单次迭代核心（简化，未含松弛因子和归一化） est = zeros(N, N, N); % 初始估计，通常全零或 FDK 结果 for iter = 1:nIter for k = 1:nAng % 前向投影当前估计 fp = forwardProjectSimple(est, R{k}, SDD, SOD, detN); % 残差 resid = proj(:,:,k) - fp; % 残差回投（用背投算子） est = est + lambda * backProject(resid, R{k}, SDD, SOD, N); end % 可选：非负约束，抑制负值伪影 est(est < 0) = 0; end逻辑说明：lambda是松弛因子，太大震荡不收敛，太小收敛慢，一般从 0.1 到 1 之间试。nIter不是越大越好，迭代到一定次数后噪声会被放大，出现“半收敛”现象，通常几十次就要停。非负约束很关键，因为密度不可能为负，加上它能明显压住伪影。实际工程里会用 OS-SART 把角度分成若干子集，每次只用一个子集更新，收敛快好几倍。3.3 用 FDK 结果做初值：省一半迭代次数纯迭代从零开始收敛慢，一个实用技巧是先用 FDK 重建一个粗糙结果，再拿它当迭代初值。FDK 已经把主要结构还原出来了，迭代只需要修伪影和补细节，迭代次数能砍掉一半以上。% 用 FDK 结果初始化迭代 est = fdkReconstruct(proj, R, SDD, SOD, detN, N); % 先跑 FDK est(est < 0) = 0; % 非负截断 for iter = 1:nIter % 后续 SIRT/OS-SART 迭代，逻辑同上 end逻辑说明：fdkReconstruct就是第 2 章那套滤波加背投的封装。用 FDK 初值的好处是起点已经接近真值，残差小，迭代方向更稳。注意 FDK 结果里的伪影如果太强，也可能把迭代带偏，所以初值最好先做一次轻度平滑或非负截断。参数上，nIter可以比零初值时少一半，lambda也可以适当调小。4. 避坑与排查：CBCT 重建最容易翻车的五个地方4.1 重建出来是镜像或旋转错位现象：重建结构和真值形状对，但左右反了或者转了个角度。原因：前向投影和背投用的旋转矩阵方向不一致，或者角度符号约定相反。解决：固定一套坐标约定，前向和背投共用同一个R{k}生成函数，先用一个不对称体模（比如偏心小球）验证方向，别用中心对称体模，那样看不出镜像。4.2 重建亮度中间亮边缘暗现象：重建切片中心亮、四周暗，像打了个聚光灯。原因：背投时漏了距离加权因子，或者加权指数写错。解决：检查(SOD/(SOD-p(2)))^2这一项，确认p(2)是体素在旋转后坐标系里沿源方向的坐标，符号别搞反。锥角越大，这个因子影响越明显。4.3 滤波后投影出现振铃和过冲现象：滤波后投影在物体边缘出现明显亮暗条纹。原因：斜坡滤波器没加窗，或者频率轴刻度detSpacing填错。解决：换 Shepp-Logan 或 Hamming 窗，核对detSpacing和投影生成时一致。如果投影本身有截断（物体超出探测器），振铃会更严重，需要先做投影延拓或截断校正。4.4 迭代重建越跑越糊、噪声越大现象：SIRT 迭代到几十次后，图像开始出现噪点，细节反而变差。原因：半收敛现象，迭代到后期噪声被放大。解决：设早停准则，比如残差下降到某个阈值就停，或者固定迭代次数在 20 到 50 之间。加非负约束和轻度正则（如全变分）也能缓解。别迷信“迭代越多越好”。4.5 内存爆掉或跑一晚上没结果现象：体数据一大，Matlab 直接 Out of memory，或者三重循环跑几小时。原因：体素驱动背投是 O(N^3 * nAng)，内存和算力都吃不消。解决：用single精度代替double，把背投向量化，或者用gpuArray上 GPU。再不行就分块重建，把体数据切成几块分别处理。参数上，先把N降到 64 或 128 验证流程，跑通再放大。5. 进阶技巧：用锥角校正和正则化把重建质量再抬一档FDK 在锥角大时上下层会糊，一个实用补救是做锥角加权再滤波，也就是在滤波前给投影乘一个与锥角相关的权重，压低大锥角射线的影响。这个权重没有统一公式，常见做法是按射线与中心平面的夹角做余弦加权，再配合探测器行方向的滤波。实现上就是在滤波前加一步逐像素乘权重，代码改动很小，但边缘层的亮度均匀性能改善不少。另一个方向是给迭代重建加正则项。最基础的全变分（TV）正则，能在压噪声的同时保边缘。做法是在每次迭代后对体数据做一次 TV 降噪，或者把 TV 梯度加进更新项。参数上，TV 权重太大图像会变“塑料感”，太小压不住噪声，一般从 0.001 到 0.01 之间试。下面是一个 TV 降噪的简化实现。function v = tvDenoise(u, weight, nStep) % 对三维体数据做简单的梯度下降 TV 降噪 v = u; for t = 1:nStep [gx, gy, gz] = gradient(v); mag = sqrt(gx.^2 + gy.^2 + gz.^2) + 1e-8; % 散度近似 div = divergence(gx./mag, gy./mag, gz./mag); v = v + weight * div; end end逻辑说明：weight控制降噪强度，nStep控制迭代步数。1e-8是防止除零。这个实现是简化版，实际会用更稳的原始-对偶算法。把它插在 SIRT 每轮迭代之后，能明显压住稀疏视角带来的条纹伪影。注意 TV 会抹掉细小结构，如果关心微细裂纹，权重别调大。验证重建质量，别只看肉眼。常用指标有 RMSE、SSIM，仿真时和真值比，实测时和标准件比。我一般会固定一个含已知缺陷的体模，每次改算法都跑同一组数据，看指标和切片，避免“感觉变好了”这种玄学判断。参数记录也要养成习惯，角度数、探测器间距、滤波窗、迭代次数、松弛因子，少记一个，下次复现就得重来。这套流程我踩过最深的坑就是旋转方向搞反，调了两天才发现是矩阵转置的问题，希望帮到你。本文还有配套的精品资源，点击获取
简介：这份资源是面向计算机、电子信息工程、数学等专业学生的锥形束CT图像重建Matlab实践包，聚焦投影背投、FDK解析重建与迭代重建两类核心算法，可用于课程设计、期末大作业与毕业设计。压缩包共15个文件，以12个m脚本为主，另含2个txt说明与1个mat案例数据，整体约173KB，脚本覆盖参数配置、投影生成、滤波、加权、背投及FDK、MLEM、SART、SQS等多种重建实现，数据可直接运行验证。代码采用参数化编程，参数便于修改，思路清晰且注释详细，兼容Matlab 2014、2019a与2024a。已有148人学习，读者可借此对比解析重建与迭代重建在噪声和数据不完整场景下的表现，快速上手算法验证与原型开发。1. 从一圈投影到三维体数据：FDK 与迭代重建到底在算什么手里只有一圈 X 光投影图，却想还原出物体内部的三维密度分布，这就是锥形束 CT（CBCT）重建要解决的核心问题。工业无损检测、口腔颌面成像、小型样品三维扫描，很多场景用的都是锥束几何：射线源和面阵探测器绕物体转一圈，采到几百上千张投影，然后靠算法把这堆二维投影“背投”回三维空间。标题里的 FDK 就是这类圆轨迹锥束重建最经典的解析算法，而迭代重建则是当投影稀疏、噪声大、金属伪影重时更稳的另一条路。Matlab 示例的价值在于：它把投影生成、滤波、背投、迭代更新这几步拆开，让你能改参数、看中间结果，而不是调一个黑匣子函数了事。这篇笔记面向已经懂点 CT 原理、想自己动手复现的工程师和学生，从几何参数一路讲到避坑。2. FDK 解析重建：从滤波到锥束背投的完整链路2.1 为什么圆轨迹锥束能用 FDK 近似FDK 的本质是把三维锥束问题拆成一层层二维扇束问题来近似处理。它的推导前提是圆轨迹、射线源到旋转中心距离固定、探测器平面与旋转轴平行。严格来说圆轨迹不满足 Tuy 数据充分性条件，靠近上下两端的大锥角层面会有明显误差，所以 FDK 在锥角小的时候才准。工程上判断能不能用 FDK，先看半锥角：一般控制在几度以内，重建质量还能接受；锥角一大，远端层面就会出现亮度不均和模糊。这也是为什么很多工业 CBCT 会把探测器做小、把源到探测器距离拉长，就是为了压锥角。理解这一点很关键：FDK 不是“精确重建”，而是“够用的近似”。你如果拿它去重建大锥角的平板探测器数据，然后抱怨边缘层糊，那不是代码写错了，是算法前提不满足。常见做法是先做一次 FDK 看整体，再决定要不要上迭代或者做锥角校正。2.2 投影生成：用 Matlab 造一份可复现的仿真数据要验证重建代码，先得有可控的投影数据。最稳的办法是自己造一个三维体模（phantom），前向投影得到正弦图，这样真值已知，重建对错一眼能看出来。下面这段代码生成一个含球和圆柱的简单体模，并按圆轨迹锥束几何做前向投影。% 生成三维体模并做锥束前向投影（简化版，便于理解流程） N = 128; % 体数据边长 vol = zeros(N, N, N); % 初始化体模 [X, Y, Z] = ndgrid(1:N, 1:N, 1:N); % 中心放一个球 vol((X-N/2).^2 + (Y-N/2).^2 + (Z-N/2).^2 < (N/6)^2) = 1; % 旁边放一个圆柱，用来观察边缘和伪影 vol((X-N/2-20).^2 + (Y-N/2).^2 < (N/10)^2 & abs(Z-N/2) < N/4) = 0.6; nAng = 360; % 投影角度数 ang = linspace(0, 2*pi, nAng+1); ang(end) = []; SDD = 400; % 源到探测器距离 SOD = 200; % 源到旋转中心距离 detN = 160; % 探测器像素数（面阵 detN x detN） proj = zeros(detN, detN, nAng); for k = 1:nAng % 把体模旋转到当前角度，等价于源探绕物体转 R = [cos(ang(k)) -sin(ang(k)) 0; sin(ang(k)) cos(ang(k)) 0; 0 0 1]; % 这里用简单射线驱动求和代替精确线积分，重在流程演示 proj(:,:,k) = forwardProjectSimple(vol, R, SDD, SOD, detN); end逻辑说明：vol是待重建的真值，ang是采样角度，SDD/SOD决定几何放大比。forwardProjectSimple是自定义函数，沿每条射线累加体素值。参数上，nAng越大角度采样越密，重建越稳但越慢；detN决定探测器分辨率，太小会丢细节，太大则内存吃紧。实际做仿真时，射线驱动求和精度不够会引入误差，建议用 Siddon 类精确线积分或者 Matlab 自带的radon做二维验证，再扩展到三维。2.3 滤波：斜坡滤波器为什么必须加窗FDK 的第一步是对每张投影做滤波，核心是斜坡滤波器（Ram-Lak），它补偿背投带来的低频堆积。但纯斜坡滤波器会把高频噪声放大，所以工程上几乎都要加窗，常见的有 Shepp-Logan、Hamming、Hann。窗越强，噪声越小，但空间分辨率也越低，这是一对硬矛盾。function pf = rampFilter(proj, detSpacing, windowType) % proj: 单张投影 (detN x detN)，对每行做一维滤波 [nRow, nCol] = size(proj); pf = zeros(nRow, nCol); % 构造频域滤波器 f = (-nCol/2:nCol/2-1) / (nCol * detSpacing); ramp = abs(f); switch windowType case 'ram-lak' w = ones(size(f)); case 'shepp-logan' w = sinc(f / (2*max(abs(f)))); case 'hamming' w = 0.54 + 0.46 * cos(pi * f / max(abs(f))); otherwise w = ones(size(f)); end H = ramp .* w; H = ifftshift(H); for i = 1:nRow row = proj(i, :); pf(i, :) = real(ifft(fft(row) .* H)); end end逻辑说明：滤波器在频域构造，ramp是斜坡项，w是窗函数。detSpacing是探测器像素间距，它决定频率轴刻度，填错会让滤波器尺度失真，重建出来要么过锐要么过糊。参数选择上，噪声大的数据用 Hamming 或 Hann，追求分辨率用 Ram-Lak 或 Shepp-Logan。注意滤波要沿探测器水平方向（对应扇束的径向），别滤错方向，否则重建直接翻车。2.4 锥束背投：体素驱动还是射线驱动背投是把滤波后的投影按几何关系累加回三维体素。两种实现思路：体素驱动，遍历每个体素找它投影到探测器哪个位置；射线驱动，遍历每条射线累加路径上的体素。体素驱动更适合并行、内存访问规律，射线驱动精度高但慢。FDK 里常用体素驱动加距离加权。% 体素驱动背投核心片段（简化） recon = zeros(N, N, N); for k = 1:nAng P = pf(:,:,k); % 当前角度滤波后投影 for iz = 1:N for iy = 1:N for ix = 1:N % 把体素坐标旋转到当前角度 p = R{k} * ([ix; iy; iz] - [N/2; N/2; N/2]); % 锥束投影到探测器 u = SDD * p(1) / (SOD - p(2)); v = SDD * p(3) / (SOD - p(2)); % 转成探测器索引并双线性插值 ui = u / detSpacing + detN/2 + 0.5; vi = v / detSpacing + detN/2 + 0.5; if ui>1 && ui<detN && vi>1 && vi<detN recon(ix,iy,iz) = recon(ix,iy,iz) + ... interp2(P, ui, vi, 'linear', 0) * (SOD/(SOD-p(2)))^2; end end end end end逻辑说明：(SOD/(SOD-p(2)))^2是锥束背投的距离加权因子，少了它重建亮度会随位置漂移。interp2做双线性插值，最近邻插值会让重建出现块状伪影。参数上，detSpacing必须和投影生成时一致，旋转矩阵R{k}的方向要和前向投影严格对应，差一个转置就会得到镜像或错位的结果。这段三重循环在 Matlab 里很慢，实际会用permute、向量化或者gpuArray加速。3. 迭代重建：当 FDK 不够用时怎么把它按在地上摩擦3.1 迭代重建的建模思路与代价迭代重建把重建当成一个优化问题：找一组体素值，使它的前向投影和实测投影最接近。最基础的是 ART/SIRT，进阶的有 SART、OS-SART，再往上还有带正则项的最大似然类方法。它的好处是能在投影少、噪声大、有金属的情况下通过先验约束压住伪影，代价是计算量比 FDK 大一到两个数量级，而且参数没调好会越迭代越糟。选型上，如果投影角度足够密（比如一圈几百张）、锥角小、噪声可控，FDK 又快又稳，没必要上迭代。反过来，稀疏视角（几十张）、低剂量高噪声、或者有金属植入物，迭代重建的优势才明显。别为了“高级”而用迭代，先问数据条件够不够。3.2 SIRT 的一次迭代：前向投影加残差回投SIRT 每次迭代做两件事：用当前估计做前向投影，算残差，再把残差按几何回投修正体数据。下面给出单步核心。% SIRT 单次迭代核心（简化，未含松弛因子和归一化） est = zeros(N, N, N); % 初始估计，通常全零或 FDK 结果 for iter = 1:nIter for k = 1:nAng % 前向投影当前估计 fp = forwardProjectSimple(est, R{k}, SDD, SOD, detN); % 残差 resid = proj(:,:,k) - fp; % 残差回投（用背投算子） est = est + lambda * backProject(resid, R{k}, SDD, SOD, N); end % 可选：非负约束，抑制负值伪影 est(est < 0) = 0; end逻辑说明：lambda是松弛因子，太大震荡不收敛，太小收敛慢，一般从 0.1 到 1 之间试。nIter不是越大越好，迭代到一定次数后噪声会被放大，出现“半收敛”现象，通常几十次就要停。非负约束很关键，因为密度不可能为负，加上它能明显压住伪影。实际工程里会用 OS-SART 把角度分成若干子集，每次只用一个子集更新，收敛快好几倍。3.3 用 FDK 结果做初值：省一半迭代次数纯迭代从零开始收敛慢，一个实用技巧是先用 FDK 重建一个粗糙结果，再拿它当迭代初值。FDK 已经把主要结构还原出来了，迭代只需要修伪影和补细节，迭代次数能砍掉一半以上。% 用 FDK 结果初始化迭代 est = fdkReconstruct(proj, R, SDD, SOD, detN, N); % 先跑 FDK est(est < 0) = 0; % 非负截断 for iter = 1:nIter % 后续 SIRT/OS-SART 迭代，逻辑同上 end逻辑说明：fdkReconstruct就是第 2 章那套滤波加背投的封装。用 FDK 初值的好处是起点已经接近真值，残差小，迭代方向更稳。注意 FDK 结果里的伪影如果太强，也可能把迭代带偏，所以初值最好先做一次轻度平滑或非负截断。参数上，nIter可以比零初值时少一半，lambda也可以适当调小。4. 避坑与排查：CBCT 重建最容易翻车的五个地方4.1 重建出来是镜像或旋转错位现象：重建结构和真值形状对，但左右反了或者转了个角度。原因：前向投影和背投用的旋转矩阵方向不一致，或者角度符号约定相反。解决：固定一套坐标约定，前向和背投共用同一个R{k}生成函数，先用一个不对称体模（比如偏心小球）验证方向，别用中心对称体模，那样看不出镜像。4.2 重建亮度中间亮边缘暗现象：重建切片中心亮、四周暗，像打了个聚光灯。原因：背投时漏了距离加权因子，或者加权指数写错。解决：检查(SOD/(SOD-p(2)))^2这一项，确认p(2)是体素在旋转后坐标系里沿源方向的坐标，符号别搞反。锥角越大，这个因子影响越明显。4.3 滤波后投影出现振铃和过冲现象：滤波后投影在物体边缘出现明显亮暗条纹。原因：斜坡滤波器没加窗，或者频率轴刻度detSpacing填错。解决：换 Shepp-Logan 或 Hamming 窗，核对detSpacing和投影生成时一致。如果投影本身有截断（物体超出探测器），振铃会更严重，需要先做投影延拓或截断校正。4.4 迭代重建越跑越糊、噪声越大现象：SIRT 迭代到几十次后，图像开始出现噪点，细节反而变差。原因：半收敛现象，迭代到后期噪声被放大。解决：设早停准则，比如残差下降到某个阈值就停，或者固定迭代次数在 20 到 50 之间。加非负约束和轻度正则（如全变分）也能缓解。别迷信“迭代越多越好”。4.5 内存爆掉或跑一晚上没结果现象：体数据一大，Matlab 直接 Out of memory，或者三重循环跑几小时。原因：体素驱动背投是 O(N^3 * nAng)，内存和算力都吃不消。解决：用single精度代替double，把背投向量化，或者用gpuArray上 GPU。再不行就分块重建，把体数据切成几块分别处理。参数上，先把N降到 64 或 128 验证流程，跑通再放大。5. 进阶技巧：用锥角校正和正则化把重建质量再抬一档FDK 在锥角大时上下层会糊，一个实用补救是做锥角加权再滤波，也就是在滤波前给投影乘一个与锥角相关的权重，压低大锥角射线的影响。这个权重没有统一公式，常见做法是按射线与中心平面的夹角做余弦加权，再配合探测器行方向的滤波。实现上就是在滤波前加一步逐像素乘权重，代码改动很小，但边缘层的亮度均匀性能改善不少。另一个方向是给迭代重建加正则项。最基础的全变分（TV）正则，能在压噪声的同时保边缘。做法是在每次迭代后对体数据做一次 TV 降噪，或者把 TV 梯度加进更新项。参数上，TV 权重太大图像会变“塑料感”，太小压不住噪声，一般从 0.001 到 0.01 之间试。下面是一个 TV 降噪的简化实现。function v = tvDenoise(u, weight, nStep) % 对三维体数据做简单的梯度下降 TV 降噪 v = u; for t = 1:nStep [gx, gy, gz] = gradient(v); mag = sqrt(gx.^2 + gy.^2 + gz.^2) + 1e-8; % 散度近似 div = divergence(gx./mag, gy./mag, gz./mag); v = v + weight * div; end end逻辑说明：weight控制降噪强度，nStep控制迭代步数。1e-8是防止除零。这个实现是简化版，实际会用更稳的原始-对偶算法。把它插在 SIRT 每轮迭代之后，能明显压住稀疏视角带来的条纹伪影。注意 TV 会抹掉细小结构，如果关心微细裂纹，权重别调大。验证重建质量，别只看肉眼。常用指标有 RMSE、SSIM，仿真时和真值比，实测时和标准件比。我一般会固定一个含已知缺陷的体模，每次改算法都跑同一组数据，看指标和切片，避免“感觉变好了”这种玄学判断。参数记录也要养成习惯，角度数、探测器间距、滤波窗、迭代次数、松弛因子，少记一个，下次复现就得重来。这套流程我踩过最深的坑就是旋转方向搞反，调了两天才发现是矩阵转置的问题，希望帮到你。本文还有配套的精品资源，点击获取