企业官网建设流程全解析

简介：一套面向医学成像、图像处理及CT重建研究者的Matlab实现资源，专注于锥形束CT（CBCT）投影与反投影链路中的FDK算法。FDK是圆轨迹锥束重建的经典近似方法，代码将算法拆解为加权、滤波、反投影三个独立阶段，每个阶段均有对应脚本，便于从原理层面逐模块阅读和调试，也适合用作研究生课程设计或算法对比的起点。除了FDK之外，包内还提供MLEM、SART等迭代重建及投影数据生成脚本，可辅助不同重建方法的横向对比验证。压缩包共十七个文件，以十二个m脚本为主体，涵盖投影生成、Parker加权、滤波、反投影、参数设定与重建入口等功能；另含xml配置文件、Matlab数据文件、png示意截图及txt说明文档，整体仅200KB，结构精炼实用；截至目前，已有1097人学习下载。读者可借助这份代码快速跑通FDK重建流程，理解锥束几何下加权与滤波因子设置，并在此基础上扩展自己的实验算法。

1. 项目速览：一套能跑通的 3D CBCT 锥束投影与 FDK 重建闭环

这个项目解决的是 CBCT（Cone Beam CT）三维重建里最经典的一环：给你一组平板探测器拍到的锥束投影，如何用 FDK 算法在 Matlab 里把三维体数据重建出来。它不依赖任何商用 CT 系统，自己带体模、正投影模拟器、FDK 三步重建和可视化对比，从光源几何到断层图像全覆盖。适合两类人：一类是做 CT 重建算法验证、课程设计、毕业设计的同学，另一类是刚接手 CBCT 项目、需要快速理解“投影到重建”全链路并做二次开发的工程师。全文会按“算法原理 → 文件结构 → 核心代码 → 参数调优 → 踩坑记录 → 扩展改造”的顺序讲，你可以直接用这套代码跑通你的第一套 FDK 重建。

2. FDK 算法拆解：从锥束几何到三步重建的数学动机

2.1 为什么 FDK 是“加权—滤波—反投影”三步

二维 CT 的 FBP（滤波反投影）大家都很熟：把每个角度的投影先做斜坡滤波，再沿射线方向反投影累加。但锥束 CT 的射线是发散的，每条 X 射线从点源出发覆盖一个锥体，不再是平行束，直接套 FBP 会有锥角伪影。

FDK 的核心思想是近似：当锥角不大时，把探测器上的每一行看作一个“倾斜的扇形束”，先对锥角方向做一次几何补偿（加权），再沿探测器行方向做斜坡滤波，最后逐体素反投影并做距离加权。具体到代码，就是fdkFilter里的加权与滤波，加上fdkBackProject里的反投影循环。加权补偿了射线倾斜带来的衰减差异，滤波压制了反投影带来的低频堆积，反投影的距离加权对应了锥束射线的发散程度。

这套三步分离的结构有个工程上的好处：你可以单独替换任何一步。比如把斜坡滤波从频域实现换成空间卷积核，或者把反投影的插值从双线性换成三线性，都不需要动另外两步。后面二次开发章节我会展开讲怎么改。

2.2 几何约定：坐标系、SOD/SDD、探测器坐标映射

FDK 所有公式都建立在一套明确的几何约定上，这是最容易翻车的地方。本项目采用如下约定：

• 旋转中心在三维坐标系原点，物体绕 z 轴旋转。
• X 射线源位于 y 轴正方向，坐标为 (0, SOD, 0)，SOD 是源到旋转中心的距离。
• 平板探测器位于源的对侧，源到探测器平面的垂直距离为 SDD，因此探测器平面在 y = SOD - SDD 处。
• 探测器行方向对应 v 轴（z 方向），列方向对应 u 轴（x 方向），探测器中心正对旋转中心。

对某个体素坐标 (x, y, z)，在投影角度 θ 下先做旋转坐标变换：

x' = x·cosθ + y·sinθ
y' = -x·sinθ + y·cosθ

然后计算源到体素的纵向距离 d = SOD - y'，锥束放大系数 k = SDD / d，最终探测器坐标：

u = x' · SDD / (SOD - y')
v = z · SDD / (SOD - y')

这两组公式是整个项目的地基。正投影用它把体素映射到探测器，反投影用它把探测器上的滤波值映射回体素。符号方向一定要和代码保持一致，否则会出现镜像重建，后面避坑章节细说。

2.3 项目文件结构与数据流

整个项目是标准的函数式组织，每个环节一个文件，主脚本只负责串联。目录结构如下：

数据流是单向的：vol → forwardProject → proj → fdkFilter → pFilt → fdkBackProject → rec。也就是说，你可以随时把forwardProject替换成读真实 CBCT 投影数据的代码，后面接的三步 FDK 完全不用改。正投影模拟器在这里既是验证工具，也是算法调试的参照物，因为体模是已知的，重建结果可以直接跟原始体模做对比。

3. 核心代码逐段讲解：从体模到重建图像的完整链路

3.1 三维体模生成：用四个椭球搭一个能看出方位的测试目标

FDK 验证需要一个“你知道正确答案”的体模。Matlab 自带的phantom只有二维版本，这里直接用椭球参数写了一个简化三维版：

function vol = makeSheppLogan3D(N, L) % 生成三维椭球体模，N 为体素网格数，L 为 FOV 边长（mm） half = L / 2; xv = linspace(-half, half, N); [xx, yy, zz] = ndgrid(xv, xv, xv); vol = zeros(N, N, N); % 每行: [cx, cy, cz, a, b, c, phi(deg), density] els = [... 0, 0, 0, 0.69, 0.92, 0.81, 0, 1.0; % 外壳 0, -0.0184, 0, 0.66, 0.87, 0.78, 0, -0.9; % 内部低密度区 -0.22, 0, -0.25, 0.41, 0.16, 0.21, -18, -0.3; % 左下暗区 0.22, 0, -0.25, 0.31, 0.11, 0.22, 18, -0.3; % 右下暗区 0, 0.35, -0.25, 0.21, 0.25, 0.50, 0, 0.3]; % 顶部亮点 for k = 1:size(els, 1) phi = els(k, 7) * pi / 180; xr = (xx - els(k,1)) * cos(phi) + (yy - els(k,2)) * sin(phi); yr = -(xx - els(k,1)) * sin(phi) + (yy - els(k,2)) * cos(phi); zr = zz - els(k,3); mask = (xr/els(k,4)).^2 + (yr/els(k,5)).^2 + (zr/els(k,6)).^2 <= 1; vol(mask) = vol(mask) + els(k, 8); end % 在体模边缘放一个针尖状高亮点，用于检查重建的点扩散表现 pk = round([N*0.45, N*0.45, N*0.55]); vol(pk(1), pk(2), pk(3)) = 2; end

参数说明：每个椭球用中心坐标 (cx, cy, cz)、三个半轴 (a, b, c)、绕 z 轴的旋转角 phi 和密度值描述；椭球之间用代数叠加，密度可为负值表示低密度区域。ndgrid生成的是世界坐标网格，所有椭球判断都用相同网格。

设计逻辑：四个椭球的朝向各异，加上一个偏移的针尖亮点，这样重建后只要看任意一个剖面图，就能立刻发现左右是否镜像、轴向是否旋转、点扩散是否正常。如果只是用一个同心圆体模，很多几何错误会被对称性掩盖，这是我反复强调的一点。

3.2 锥束正投影：体素驱动的射线映射

正投影模拟器的作用是：给定三维体模和几何参数，算出每个角度下探测器上应该接收到多少信号。实现上采用体素驱动，把每个体素中心投影到探测器平面然后累加：

function proj = forwardProject(vol, geo) N = size(vol, 1); half = geo.L / 2; xv = linspace(-half, half, N); [xx, yy, zz] = ndgrid(xv, xv, xv); R = geo.detRows; C = geo.detCols; p = geo.detPixel; proj = zeros(R, C, geo.nAng, 'single'); for ia = 1:geo.nAng th = geo.angles(ia); xr = xx .* cos(th) + yy .* sin(th); yr = -xx .* sin(th) + yy .* cos(th); d = geo.SOD - yr; % 源到体素的纵向距离 k = geo.SDD ./ d; % 锥束放大倍数 uu = xr .* k; vv = zz .* k; ci = round(uu / p + (C + 1) / 2); % 探测器列索引 ri = round(vv / p + (R + 1) / 2); % 探测器行索引 valid = d > 0.1 & ci >= 1 & ci <= C & ri >= 1 & ri <= R; idx = ri(valid) + (ci(valid) - 1) * R; vals = double(vol(valid)); proj(:, :, ia) = accumarray(idx, vals, [R * C, 1]); end end

逻辑说明：每个角度先对体素坐标做旋转（xr, yr），再用d = SOD - yr算出射线纵向距离。放大系数 k 与geo.SDD/d成正比，反映锥束发散程度。最后accumarray把落在同一个探测器像素内的所有体素值累加，得到该角度投影。

参数说明：ci/ri是把毫米坐标换算到像素索引，(C+1)/2是让探测器中心对齐到索引中心；valid剔除了源后方的体素（d <= 0.1）和落在探测器视野外的体素。当前版本用最近邻散点，体素尺寸与探测器像素接近时精度够用；想要更平滑的投影可改成双线性分摊，就把round改为对四个邻域像素各分配权重。

3.3 FDK 前置处理：余弦加权与斜坡滤波

function pF = fdkFilter(proj, geo) R = geo.detRows; C = geo.detCols; p = geo.detPixel; uAxis = ((1:C) - (C + 1) / 2) * p; vAxis = ((1:R) - (R + 1) / 2) * p; [vGrid, uGrid] = ndgrid(vAxis, uAxis); W = geo.SDD ./ sqrt(geo.SDD.^2 + uGrid.^2 + vGrid.^2); pW = proj .* W; % 斜坡滤波：沿探测器列方向（u轴） f = (-C/2 : C/2 - 1) / C; % 归一化频率 [-0.5, 0.5) ramp = abs(f); hamm = 0.54 + 0.46 * cos(2 * pi * f); % Hamming 窗抑制振铃 H = ifftshift(ramp .* hamm); % 转到 fft 频率顺序 pF = zeros(size(pW), 'single'); for ia = 1:geo.nAng tmp = fft(pW(:, :, ia), [], 2); pF(:, :, ia) = real(ifft(tmp .* H, [], 2)); end end

逻辑说明：加权矩阵 W 对每个探测器像素计算射线与中心射线的夹角补偿，公式SDD / sqrt(SDD^2 + u^2 + v^2)本质是夹角的余弦倒数。越靠近探测器边缘，锥角越大，补偿越强。滤波部分采用频域乘法：对每帧投影沿列方向做 FFT，乘斜坡滤波器，再 IFFT 回来。hamm窗显著减少强边缘的振铃伪影，代价是轻微牺牲空间分辨率。

参数说明：f是归一化频率，ramp = abs(f)是理想斜坡滤波器。注意ifftshift必须用对，否则滤波器在频域的位置偏半个周期，重建结果会出现整体空间偏移。没有加窗的纯斜坡滤波在高频处放大噪声，实际项目中我基本都保留 Hamming 窗。

3.4 FDK 反投影主循环：逐角度体素累加

反投影是整个项目里计算量最大的一段，也是 FDK 公式落地的地方：

function rec = fdkBackProject(pF, geo) N = geo.N; half = geo.L / 2; xv = linspace(-half, half, N); [xx, yy, zz] = ndgrid(xv, xv, xv); R = geo.detRows; C = geo.detCols; p = geo.detPixel; uAxis = ((1:C) - (C + 1) / 2) * p; vAxis = ((1:R) - (R + 1) / 2) * p; rec = zeros(N, N, N, 'single'); dAng = abs(geo.angles(2) - geo.angles(1)); for ia = 1:geo.nAng th = geo.angles(ia); xr = xx .* cos(th) + yy .* sin(th); yr = -xx .* sin(th) + yy .* cos(th); d = geo.SOD - yr; valid = d > 0.05 * geo.SDD; % 防止源平面附近的坐标爆炸 k = geo.SDD ./ max(d, geo.SDD * 0.05); uu = xr .* k; vv = zz .* k; ci = uu / p + (C + 1) / 2; % 列坐标（允许小数） ri = vv / p + (R + 1) / 2; % 行坐标（允许小数） ci(~valid) = NaN; ri(~valid) = NaN; val = interp2(pF(:, :, ia), ci, ri, 'linear', 0); val(isnan(val)) = 0; w = (geo.SOD ./ max(d, geo.SDD * 0.05)).^2; rec = rec + single(val .* w); end rec = rec * dAng / (2 * pi); end

逻辑说明：每个角度下，先把体素坐标旋转到探测器坐标系，用与正投影相同的公式算出投影坐标，然后在滤波后的投影上用interp2双线性插值取值。累加前乘权重 w = (SOD / d)^2，这是锥束反投影的距离平方加权：离源越近的体素，射线发散角越大，单位体素在探测器上占的面积越小，需要补偿的权重越大。最后乘以dAng / (2*pi)做角度归一化。

参数说明：max(d, geo.SDD*0.05)是个 clamp，防止正好落在源平面附近的体素产生除零和 NaN。interp2的 fill value 设为 0，投影视野外的体素直接取零。这里用的是双线性插值，比最近邻平滑得多，重建质量更好，代价是插值计算耗时增加。

4. 跑通与参数调优：几何参数、扫描协议与性能取舍

4.1 主脚本与最小配置

runDemo.m只做五件事，按顺序执行就是完整链路：

%% CBCT FDK 演示：体模 → 正投影 → 滤波 → 反投影 → 显示 clear; clc; geo = configGeometry(); % 见下方参数表 vol = makeSheppLogan3D(geo.N, geo.L); proj = forwardProject(vol, geo); pFilt = fdkFilter(proj, geo); rec = fdkBackProject(pFilt, geo); viewResult(vol, rec, proj, geo);

configGeometry.m里是我推荐的一组初始参数：体素网格 128，体模边长 200 mm，SOD 400 mm，SDD 700 mm，探测器像素 1.5 mm，探测器 192 × 192，扫描角度 120 个。这套参数下，正投影加反投影在主流笔记本上大概两到三分钟跑完。首次运行时先不要改参数，确认能出图再动几何。

4.2 几何参数与图像质量的对应关系

调试建议：先固定 N=64、nAng=60 跑通全流程，确认重建图像结构正确后再放大。放大时优先增加 nAng，因为 FDK 对角度欠采样的伪影最敏感；提高 N 带来的分辨率收益会被双线性插值抹平一部分，性价比低于加角度。

4.3 性能优化：单精度、逐角度处理、GPU 加速

这套代码的性能瓶颈在反投影：体素数 N³ 乘以角度数就是插值次数，N=128、nAng=120 时约 2.5 亿次插值。我做了三个层面的优化：

• 全程使用single精度。投影、滤波结果、重建体数据都是单精度，内存减半，速度也更快。体模生成时用 double，但在进入正投影前转成 single。
• 反投影内层向量化。逐角度循环，但每个角度内的三维坐标运算全部用ndgrid数组向量化，避免三重嵌套循环。
• 可选 GPU 加速。装了 Parallel Computing Toolbox 的话，把xx/yy/zz丢到gpuArray上，反投影速度能有几十倍提升，最后gather回 CPU 即可。

提示：正投影的accumarray在 GPU 上支持不理想，建议正投影留在 CPU 跑，只把反投影的坐标运算和插值放到 GPU。

5. 避坑记录：FDK 重建翻车现场与排查方法

5.1 重建图像左右镜像，体模结构完全反过来

现象：输出的重建图里，体模左下的暗区出现在右上，整个结构沿某个轴翻转了。

原因：旋转坐标变换的符号与探测器坐标方向不一致。y' 的计算里用的是sin还是-sin，决定了旋转方向；探测器 v 轴向上还是向下，决定了行索引方向。两处符号只要错一处，就会出现镜像。

解决：把正投影和反投影的坐标公式抄在一张纸上，确认两者用的是同一套约定。我排错时会在正投影输出里看中间角度的投影图，对着体模结构看左右是否颠倒；如果正投影是对的，那问题一定在反投影的符号上。

5.2 强边缘处出现放射状亮纹

现象：高密度椭球边界往外辐射出明显的亮线或暗线，像星芒一样。

原因：斜坡滤波放大了高频噪声，强边缘本身的 Gibbs 振铃又被反投影累加。如果滤波核没有加窗，这种现象非常明显。

解决：在fdkFilter里把ramp乘上hamm窗，振铃立刻大幅缓解。代价是高频信息略有损失，边缘轻微变钝。要验证是不是窗函数的问题，把窗去掉跑一遍对比就知道。

5.3 重建图出现 NaN 或极亮的散点

现象：图像上星星点点地分布着 NaN（显示为黑点）或异常亮的像素。

原因：反投影过程中出现了除零。体素恰好落在源平面附近（d 趋近 0），或者投影坐标落在探测器边缘外但插值越界逻辑没兜住。

解决：在fdkBackProject里对 d 做 clamp，并把无效区域的插值坐标设为 NaN，插值结果里再把 NaN 置零。这两步缺一不可，只做 clamp 不处理越界坐标，边缘依然会冒点。

5.4 重建结果灰蒙蒙一片，对比度明显不对

现象：结构都在，但整体灰度被压缩，暗区和亮区的差异远小于体模。

原因：加权与反投影权重里，SOD 和 SDD 混用了。加权用的是源到探测器的距离 SDD，反投影距离加权的基准必须是源到旋转中心的距离 SOD。两处基准不一致时，锥角补偿会出现系统性误差，图像整体发灰。

解决：回到 2.2 节的公式，把加权和反投影权重的距离基准分开定义。还有个小细节：滤波核在归一化频率下没有乘以 1/像素尺寸，绝对灰度会差一个常数因子，但相对对比度不受影响。定量重建时补上这个因子即可。

5.5 内存直接溢出，N=256 时程序卡死

现象：把 N 从 128 调到 256 后，正投影或反投影阶段直接报内存错误，或者磁盘 swap 把机器拖死。

原因：double 精度下，256³ 体素是 1.3 亿个数，约 1 GB；反投影时坐标数组又额外占了几个 GB。整条链路叠加起来，16 GB 内存根本撑不住。

解决：全部改 single 精度，坐标网格只在函数开头生成一次，不要在角度循环里重复创建；反投影保持逐角度循环，单帧投影只有 192 × 192，内存占用很小。如果还想冲更高分辨率，就把反投影改成按体素分块（block-wise）处理，一次只算一个子块。

6. 二次开发：从 FDK 演示到实用重建系统的三个改造方向

第一优先级是替换数据源。把forwardProject(vol, geo)换成dicomread读入真实 CBCT 投影序列，前面体模和正投影全部删掉，后面fdkFilter和fdkBackProject原样复用。真实系统里要重点核对探测器的暗场校正和增益校正，以及投影角度表是否按顺时针采样，这两点和几何约定不一致时，重建结果会比仿真体模更难排查。

第二优先级是改轨迹。把圆周轨迹改成螺旋轨迹只需在fdkBackProject的角度循环里加一个床位移量：每个角度的体素 z 坐标减去pitch * angle / (2*pi)，再按常规 FDK 反投影。螺旋 FDK 是商用 CBCT 的常见做法，改完记得用带倾斜结构的体模验证，否则看不到螺旋伪影。

第三优先级是提速。体素驱动反投影的插值坐标只依赖几何，不依赖体数据本身，可以离线预计算一张“坐标查找表”缓存到磁盘；后续每次重建直接用查表代替三角函数和坐标变换，速度提升非常明显。更激进的做法是把反投影内核写成 MEX，或者用gpuArray包住向量化运算。

每次改完几何或轨迹，我都强制自己先跑一遍“非对称体模 + 中间角度投影目检”再进完整链路。这个习惯帮我挡掉了至少三次能把结果彻底带偏的符号错误。希望帮到你。

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