企业官网建设流程全解析

简介：这份资源面向医学影像重建方向的学习者与研究人员，聚焦三维锥束CT（CBCT）的投影与反投影重建，核心是FDK算法的Matlab实现。FDK流程被拆解为投影数据加权、加权后滤波、滤波后反投影三步，最终输出各方向的重建结果，适合已具备一定图像处理基础、希望动手复现经典算法的读者。压缩包共17个文件，约200KB，以12个m脚本为主体，覆盖投影、反投影、Parker加权、滤波及FDK、SART、MLEM、SQS等多套重建示例，另含mat数据、xml配置、txt说明与png截图，便于直接运行与对照。已有1098人学习下载。读者可借此理解锥束几何下的加权与滤波细节，对比解析与迭代重建的差异，并基于现有脚本快速搭建自己的实验流程。

1. 锥束 CT 重建的起点：为什么 FDK 至今仍是工程首选

如果你手头有一组平板探测器采到的锥束投影数据，想把它重建成能看的三维体数据，大概率绕不开 FDK 这个名字。Feldkamp、Davis、Kress 在 1984 年提出的这个近似算法，本质上是把二维滤波反投影推广到三维锥束几何：先对投影做余弦加权和斜坡滤波，再沿锥束射线方向做加权反投影。它不像迭代重建那样需要几十次正反投影、吃满显存，一次滤波加一次反投影就能出结果，速度优势在工程现场几乎是决定性的。

这套 MATLAB 实现的 CBCT 投影/反投影 FDK 代码，面向的正是这类场景：圆形轨迹、平板探测器、锥束射线，输入是投影图序列，输出是三维体素。适合医学影像方向的学生做课程设计，也适合做无损检测、材料内部结构分析的工程师快速搭一套能跑通的重建链路。它不追求论文级的最优图像质量，但胜在结构清晰、参数可调、能让你把 FDK 每一步都拆开看明白。下面从几何参数怎么定、滤波怎么选、反投影怎么加速，一路讲到实际跑数据时最容易翻车的地方。

2. 几何参数与投影数据准备：先把坐标系对齐

FDK 重建翻车，十次里有七次不是算法写错，而是几何参数对不上。锥束 CT 的几何比平行束复杂得多：射线源、旋转中心、探测器平面三者之间的相对位置，任何一个量的符号或量纲搞错，重建出来就是一团糊。所以这一章先把参数定义和数据读取讲透，再谈滤波和反投影。

2.1 锥束几何的五个核心参数

圆形轨迹锥束 CT 的几何可以用一组参数完整描述。常见做法是定义一个源到旋转中心距离（Source-to-Object Distance，SOD）、源到探测器距离（Source-to-Detector Distance，SDD）、探测器像素尺寸、投影图像尺寸，以及每张投影对应的旋转角度。这五个量决定了反投影时每条射线打在哪个体素上。

放大比 MAG = SDD / SOD 是个关键中间量。它决定了探测器上测到的位移对应到物体空间是多少。如果 SOD 和 SDD 的单位不统一，比如一个用毫米一个用厘米，MAG 就会差一个数量级，反投影出来的体数据尺寸完全不对。我一般会在读参数后立刻打印 MAG，确认它在 1.5 到 3 之间，超出这个范围基本可以判定参数填错了。

角度向量也要注意。有些设备导出的是 0 到 360 度均匀分布，有些是 0 到 360 度但包含首尾重复点，还有些是弧度制。代码里如果按角度索引去查投影，角度向量的长度必须和投影张数严格一致，否则反投影时会出现某几张投影被重复使用或跳过，重建结果里会出现明显的条状伪影。

2.2 投影数据读取与预处理

MATLAB 里读投影序列，常见做法是用dicomread或imread循环读入，堆成一个三维数组。这里有个内存陷阱：512×512×360 的 uint16 数组大约 180 MB，转成 double 直接翻四倍到 720 MB，如果再做滤波和反投影的中间变量，很容易把内存吃满。所以读入后先保持 uint16，只在滤波前转 double。

% 读取投影序列，假设文件名为 proj_0001.tif 到 proj_0360.tif nProj = 360; detSize = [512, 512]; proj = zeros(detSize(1), detSize(2), nProj, 'uint16'); for i = 1:nProj fname = sprintf('proj_%04d.tif', i); proj(:,:,i) = imread(fname); end % 转 double 并做暗场/亮场归一化（如果有的话） proj = double(proj); % 假设已采集暗场 dark 和亮场 flat % proj = (proj - dark) ./ (flat - dark); % proj = max(proj, 0); % 避免负值

这段代码的逻辑很直接：按文件名模式循环读入，堆成三维数组。参数nProj和detSize必须和实际数据一致，否则数组维度对不上，后面反投影索引会越界。暗场和亮场归一化是可选步骤，但如果你的数据没有做过，重建出来的体数据会带一圈明显的环形伪影，因为探测器各像素的增益不一致。归一化之后记得把负值截断到零，否则取对数时会出 NaN。

提示：如果投影数据是 16 位无符号整数，转 double 之前先确认最大值没有超过 65535，否则可能是读取时字节序搞反了，重建结果会完全不可用。

2.3 余弦加权与斜坡滤波的实现细节

FDK 的第一步是对每张投影做余弦加权。加权因子是探测器像素到源的距离与 SDD 的比值，本质上是在补偿锥束射线斜穿探测器时的路径长度差异。这一步不做，重建出来的图像中心会偏亮、边缘偏暗。

% 余弦加权 [u, v] = meshgrid(1:detSize(2), 1:detSize(1)); u = u - (detSize(2)+1)/2; % 中心化 v = v - (detSize(1)+1)/2; % 探测器像素物理坐标 u = u * detPixelSize; v = v * detPixelSize; % 加权因子 cosWeight = SDD ./ sqrt(SDD^2 + u.^2 + v.^2); projWeighted = proj .* repmat(cosWeight, [1, 1, nProj]);

中心化那两行很容易写错。(detSize+1)/2和detSize/2差 0.5 个像素，对于 512 的探测器来说，这 0.5 像素的偏移会在重建结果里表现为轻微的旋转中心偏移，图像边缘出现模糊。我一般会先用一个已知对称的模体（比如均匀圆柱）跑一遍，看重建切片是否左右对称，不对称就回来调这个中心。

滤波这一步，FDK 用的是斜坡滤波（Ram-Lak），也可以加窗变成 Shepp-Logan 或 Hann。斜坡滤波在频域是 |f|，直接离散化会有直流分量问题，常见做法是构造一个实数对称的滤波器向量，然后沿探测器行方向做一维卷积。

% 构造 Ram-Lak 滤波器 N = 2^nextpow2(2*detSize(2)); f = linspace(-0.5, 0.5, N); ramLak = abs(f); ramLak = ramLak .* 2; % 归一化因子 % 加 Hann 窗 hannWin = 0.5 + 0.5*cos(2*pi*f); ramLak = ramLak .* hannWin; % 转回空域 filterKernel = real(ifft(ifftshift(ramLak))); filterKernel = fftshift(filterKernel); % 沿列方向卷积 for i = 1:nProj projFiltered(:,:,i) = conv2(projWeighted(:,:,i), filterKernel, 'same'); end

滤波器长度取2^nextpow2(2*detSize(2))是为了避免循环卷积的混叠。加 Hann 窗会牺牲一点空间分辨率，但能明显压制高频噪声，实际数据噪声大时我一般会加上。卷积用conv2的'same'模式保持尺寸不变，方便后续反投影索引。如果数据量很大，逐张卷积会比较慢，可以考虑用fft做频域乘法，但要注意补零到足够长度。

3. 反投影与体数据生成：从投影到三维体素

滤波做完，接下来就是 FDK 最耗时的部分：反投影。二维滤波反投影是每条射线对应一个像素，三维锥束反投影是每条射线穿过整个体数据，计算量随体素数量线性增长。这一章讲反投影的坐标映射、体素遍历顺序，以及怎么用向量化把速度提上来。

3.1 反投影的坐标映射逻辑

反投影的核心是：对于体数据中的每个体素 (x, y, z)，找到它在每张投影上的对应位置 (u, v)，然后把滤波后的投影值累加到这个体素上。锥束几何下，源在旋转轨迹上，体素到源的连线与探测器平面的交点就是 (u, v)。

假设源在角度 θ 时的位置是 (SOD·sinθ, -SOD·cosθ, 0)，体素坐标是 (x, y, z)，那么从源到体素的射线方向向量是 (x - SOD·sinθ, y + SOD·cosθ, z)。这条射线与探测器平面的交点可以通过相似三角形算出来。探测器平面在旋转坐标系下始终垂直于源到旋转中心的连线，所以要先做一个旋转坐标变换。

% 体素网格 voxSize = 1; % mm voxDim = [256, 256, 256]; [xg, yg, zg] = meshgrid(... (1:voxDim(2)) - (voxDim(2)+1)/2, ... (1:voxDim(1)) - (voxDim(1)+1)/2, ... (1:voxDim(3)) - (voxDim(3)+1)/2); xg = xg * voxSize; yg = yg * voxSize; zg = zg * voxSize; % 反投影累加 vol = zeros(voxDim); for i = 1:nProj theta = angles(i) * pi / 180; % 旋转体素坐标到探测器坐标系 xr = xg * cos(theta) + yg * sin(theta); yr = -xg * sin(theta) + yg * cos(theta); % 放大比投影 u = xr * SDD ./ (SOD - yr); v = zg * SDD ./ (SOD - yr); % 转像素索引 uIdx = round(u / detPixelSize + (detSize(2)+1)/2); vIdx = round(v / detPixelSize + (detSize(1)+1)/2); % 边界检查 valid = uIdx >= 1 & uIdx <= detSize(2) & ... vIdx >= 1 & vIdx <= detSize(1); % 累加 projSlice = projFiltered(:,:,i); vol(valid) = vol(valid) + projSlice(sub2ind(size(projSlice), vIdx(valid), uIdx(valid))); end % 归一化 vol = vol * (pi / nProj) * (SOD / SDD)^2;

这段代码里SOD - yr是体素到源在旋转坐标系下的距离，必须为正，否则射线打到了源后面，物理上无意义。valid那个边界检查不能省，否则索引越界直接报错。最后的归一化因子(pi/nProj) * (SOD/SDD)^2是 FDK 的标准缩放，不同文献里写法略有差异，但量级要对，否则重建出来的 CT 值没有可比性。

注意：round做最近邻插值，重建结果会有块状感。如果追求质量，可以改成双线性插值，但速度会慢三到四倍。我一般先用最近邻调通流程，确认几何没问题后再换插值方式。

3.2 向量化加速与内存权衡

上面那段三重循环在 256³ 体素、360 张投影下，MATLAB 里跑一遍大概要十几分钟。瓶颈在于每张投影都要对全部体素做一次坐标计算和索引。向量化能把它压到几分钟，代价是内存占用上升。

常见做法是把体素坐标预先算好，循环里只做旋转和索引，避免重复生成网格。更进一步，可以把sub2ind换成线性索引，减少函数调用开销。如果内存够，还可以把projFiltered转成single精度，索引时自动转换，省一半内存。

% 预计算体素坐标的线性索引 voxIdx = reshape(1:numel(vol), voxDim); % 循环内用线性索引累加 vol(voxIdx(valid)) = vol(voxIdx(valid)) + projSlice(linIdx(valid));

另一个加速点是并行化。MATLAB 的parfor可以把投影循环拆到多个核上，但要注意vol是累加变量，不能直接在parfor里写，需要每张投影算一个局部体数据再求和，内存会翻几倍。我一般只在体数据小于 128³ 时用parfor，再大就老老实实串行，或者转 C mex。

3.3 重建结果的验证方法

重建跑完，怎么判断结果对不对？最直接的办法是找一个已知结构的模体，比如均匀圆柱或者带孔的圆盘，看重建切片里孔的位置和大小是否和实物一致。如果没有模体，可以拿一张投影做前向投影，再和原始投影对比，差异应该在噪声水平以内。

还有一个土办法：把重建体数据沿某个方向做最大密度投影（MIP），如果能看到清晰的边缘和结构，说明几何基本对齐了。如果 MIP 里出现重影或者拖尾，多半是角度向量和投影张数不匹配，或者旋转中心偏了。我一般会先跑一个 64³ 的小体数据，确认 MIP 正常后再放大到全尺寸，省得等半天才发现参数错了。

4. 避坑与排查：FDK 重建最常见的五类翻车

FDK 的代码结构不复杂，但参数多、坐标系绕，实际跑数据时踩坑是常态。这一章把我在调试过程中遇到最多的五类问题整理出来，每条按现象、原因、解决三步走，方便你对照排查。

4.1 重建结果全黑或全白

现象：反投影累加完，体数据要么全是零，要么全是极大值，看不到任何结构。

原因：最常见的是投影数据没做暗场亮场归一化，或者归一化时除了零。另一个可能是滤波器的直流分量没处理好，导致累加结果被一个巨大的常数淹没。

解决：先检查投影数据的直方图，确认动态范围正常。然后单独跑一张投影的滤波，看滤波后数据的均值是否在合理范围。如果均值偏离原始数据一个数量级，回去检查滤波器的构造，特别是ramLak的归一化因子。

4.2 图像中心偏移、边缘模糊

现象：重建切片里物体不在中心，或者边缘有明显的模糊带。

原因：探测器中心(detSize+1)/2和实际旋转中心不重合。平板探测器的机械安装误差、源和探测器的对准误差都会导致这个偏移。

解决：用一个已知对称的模体，重建后看切片是否左右对称。如果偏了，在uIdx计算时加一个偏移量centerOffset，手动微调，直到对称。这个偏移量通常在几个像素以内，超过十个像素说明机械安装有问题，得先调硬件。

4.3 条状伪影和环形伪影

现象：重建图像里出现放射状条纹或者同心圆环。

原因：条状伪影多半是角度向量不均匀或者某几张投影损坏。环形伪影通常是探测器像素增益不一致，暗场亮场没做好。

解决：先检查角度向量，打印相邻角度差，看是否均匀。如果有跳变，回去核对设备导出的角度文件。环形伪影的话，重新做暗场亮场归一化，或者对投影做一次沿探测器行的中值滤波，压制坏像素。

4.4 内存溢出或运行极慢

现象：MATLAB 报 Out of Memory，或者跑一晚上还没出结果。

原因：体数据太大，double 精度下 512³ 就是 1 GB，加上中间变量很容易爆。反投影的三重循环在 MATLAB 里本来就慢，没向量化的话更慢。

解决：先把体数据降到 128³ 或 64³ 调通流程，确认参数无误后再放大。精度上，投影和体数据都可以用single，内存减半，速度也快一些。反投影循环尽量向量化，或者用parfor拆到多核。

4.5 重建结果和预期方向不一致

现象：重建出来的体数据上下颠倒、左右镜像，或者旋转了 90 度。

原因：MATLAB 的数组索引是 (行, 列, 页)，而图像坐标习惯是 (x, y, z)，两者容易搞混。另外，探测器的行方向对应体数据的哪个轴，不同设备定义不一样。

解决：拿一个已知方向的模体，重建后看它在体数据里的朝向。如果不对，在反投影时交换uIdx和vIdx，或者对体数据做flip。这个没有通用解，只能对着你的数据试一次，记下来，以后就固定了。

5. 进阶技巧：用前向投影做闭环验证与参数扫描

FDK 调通之后，下一步是怎么确认重建质量真的到位了。光看切片有时候不够，因为人眼对低频的对比度不敏感，一些细微的几何误差可能被忽略。我一般会做两件事：一是用重建结果做前向投影，和原始投影对比；二是对关键参数做小范围扫描，看重建指标的敏感度。

前向投影的逻辑和反投影反过来：对于每张投影的每个像素，沿射线方向穿过体数据，累加体素值。如果几何参数完全一致，前向投影的结果应该和原始投影在噪声水平以内吻合。差异大的地方，就是几何没对齐的地方。

% 前向投影：从体数据生成投影 projFwd = zeros(detSize(1), detSize(2), nProj); for i = 1:nProj theta = angles(i) * pi / 180; for v = 1:detSize(1) for u = 1:detSize(2) % 探测器像素坐标 uu = (u - (detSize(2)+1)/2) * detPixelSize; vv = (v - (detSize(1)+1)/2) * detPixelSize; % 射线方向（旋转坐标系） dir = [uu, vv, SDD]; dir = dir / norm(dir); % 沿射线采样体数据 t = linspace(0, SOD*2, 200); xr = -SOD + t * dir(1); yr = t * dir(2); zr = t * dir(3); % 旋转回体数据坐标 x = xr * cos(theta) - yr * sin(theta); y = xr * sin(theta) + yr * cos(theta); % 三线性插值 projFwd(v, u, i) = interp3(xg, yg, zg, vol, x, y, zr, 'linear', 0); end end end

这段代码是逐像素逐投影算的，速度很慢，只适合小体数据验证。interp3的最后一个参数 0 表示超出体数据范围时返回 0，避免边界外插值引入假信号。采样点数 200 是个经验值，体数据 256³ 时够用，再大就得加。

对比前向投影和原始投影，我一般看两个指标：归一化均方误差（NMSE）和结构相似性（SSIM）。NMSE 低于 5% 说明几何基本对齐，SSIM 高于 0.9 说明结构保真度可以。如果 NMSE 偏大，先查角度向量，再查旋转中心偏移，最后查探测器像素尺寸。

参数扫描方面，最值得扫的是旋转中心偏移和角度向量。旋转中心偏移在 ±3 像素范围内扫，步长 0.5 像素，看 NMSE 的变化曲线，最低点就是最优偏移。角度向量的话，如果设备导出的角度有噪声，可以尝试用均匀角度替代，看 NMSE 是否改善。我遇到过一台设备导出的角度文件里混入了重复角度，导致重建图像里出现一对对称的条状伪影，换成均匀角度后伪影消失。

提示：参数扫描很耗时，建议在 64³ 体数据上做，找到最优参数后再放大到全尺寸跑一次确认。

从那以后我每次拿到新的投影数据，都强制走一遍前向投影闭环验证，确认 NMSE 和 SSIM 达标后才开始正式重建。这个习惯帮我省下了大量返工时间，也让我对每台设备的几何特性心里有数。希望帮到你。

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