企业官网建设流程全解析

简介：围绕锥束CT的FDK重建与迭代重建实验，这份Matlab资料面向医学成像课程学习者、算法入门者及需要复现课堂或科研实验的读者，帮助理清投影、加权、滤波、反投影到图像重建的完整链路，适合课程设计、实验复现和算法对照。压缩包共17个文件，约200KB，以12个.m函数与示例脚本为主，辅以2个.xml配置、1个.txt说明、1个.png截图及1个.mat图像数据，体量轻便，便于查阅、运行和二次修改。已有1100人学习，说明其在相关实验练习中具有一定参考价值。内容覆盖投影生成、Parker加权、滤波、反投影等环节，并包含FDK、MLEM、SART、SQS等重建示例、参数设置与测量数据生成脚本；读者可据此理解FDK三步流程，观察不同重建方法的实现差异，对锥束CT几何关系、滤波函数选择与反投影累加过程形成直观认识，进而完成结果验证与算法比较。

1. 从一张糊成水墨画的切片说起：这套 FDK 重建代码解决什么问题

第一次跑锥束 CT 重建，很多人拿到的是一张糊得厉害、几乎看不出结构的切片：模体边缘一圈圈晕开，致密结构外面裹着一层亮环，调窗宽窗位怎么调都不对。这套基于 Matlab 的 3D 锥束（CBCT）投影、反投影 FDK 代码，要处理的正是这件事：把一圈平板探测器上采到的二维投影，按 Feldkamp-Davis-Kress 近似算法重建出三维体数据。它适合两类人：手上有锥束几何的仿真或实测投影、想在 Matlab 里把整条重建链路完整跑一遍的学生；以及要把 FDK 当基线，后面接迭代重建、去噪或深度学习后处理的工程开发者。代码本身不长，但几何约定、滤波核、距离权重这三处只要错一个，出来的图就是典型的“能看不能用”。

2. 工程结构与几何约定：射线源、探测器和旋转轴到底怎么摆

读这类重建代码时，最耗时间的往往不是算法本身，而是搞清楚“作者心里的坐标系”。同一份投影数据，源在 +x 还是 -x、u 轴朝左还是朝右、旋转是顺时针还是逆时针，四种组合里只有一种是自洽的。翻车的表现形式还特别隐蔽：重建出来的模体不是完全错位，而是整体轻微拖影，看起来像“分辨率不够”，实际是几何反了。

2.1 文件分工：投影、滤波、反投影各在哪一层

一份典型的组织方式是这样，按“入口 → 预处理 → 滤波 → 反投影”四层切开，每个函数只做一件事，方便单独替换：

cbct_fdk/ ├── main_fdk.m % 入口：读投影、装配几何参数、调重建 ├── fdk_recon.m % FDK 主流程：加权 → 滤波 → 反投影 ├── cosine_weight.m % 锥角余弦加权 ├── ramp_filter.m % 斜坡滤波，可选 Ram-Lak / Shepp-Logan / Hann ├── ramp_kernel.m % 生成频域滤波核 ├── backproject.m % 加权反投影，按 z 分层循环 ├── forward_project.m % 体素驱动前投影，用来造仿真数据 ├── phantom_ellipsoid.m % 椭球模体生成器 ├── cfg/default_cfg.m % 几何与重建参数的默认值 └── data/proj_shepp.mat % 预生成投影 [nV nU nView]

这种分法的好处是替换成本低。想试不同的窗函数，只动ramp_kernel.m；想换 GPU 版反投影，只重写backproject.m，主流程一行不改。我不太建议把所有逻辑塞进一个fdk.m里——网上下载的单文件版本大多如此，两百行里混着几何、滤波和插值，改一个参数要在三个地方对齐，二次开发时非常难受。

forward_project.m这个模块容易被忽略，但它是验证环节的关键。手上没有实测数据时，全靠它把模体“打”成投影，再喂给 FDK 反推回去，形成闭环。闭环能对上，说明几何、滤波、权重三者自洽；对不上，至少能确定问题在代码里，而不是在数据里。

2.2 几何参数：SOD、SDD 与 u/v 轴的方向约定

FDK 是圆轨迹锥束的近似算法，它要求射线源在一个平面圆上运动，探测器绕旋转轴同步转动且始终垂直于中心射线。参数就两个距离加两个像元尺寸，但含义必须先说死，后面所有公式都建立在上面：

我采用的约定是：旋转轴为 z 轴，源在第 β 个视角的位置为S = (R·cosβ, R·sinβ, 0)，探测器在源的对面、距源 D，中心射线过原点。探测器坐标 (u, v) 以中心射线与探测器平面的交点为原点，u 沿旋转切线方向，v 沿 z 轴。对体素(x, y, z)，把它投到探测器上的位置是：

U = R - x·cosβ - y·sinβ 体素在中心射线方向上的等效距离 u = D · (−x·sinβ + y·cosβ) / U v = D · z / U

U这个量在代码里出现频率最高，它同时决定了投影坐标和反投影的距离权重(R/U)²。血泪经验是：U的定义一旦和源的坐标不匹配，重建图不会整体崩掉，而是中心清晰、边缘发暗，很多人第一反应是“滤波器不对”，实际上是权重反了。所以改几何时，先把S和U写在纸上对一遍，再动代码。

实测数据还有两个坑。一是探测器偏置：有些平板为了让视野更大，把探测器横向偏移了半个探元，此时 u 要加一个常量u0，并且余弦加权的形式也要微调；不做这一步，重建会出现一对镜像的鬼影。二是探测器面外倾角，FDK 假设探测器严格垂直于中心射线，机械装配的零点几度倾斜在重建里会表现为方向性模糊，这个只能靠几何标定解决，代码层面补不了。

2.3 内存布局：为什么投影数组是 [nV nU nView]

投影数组按[行 列 视角]排列，也就是[nV nU nView]，这是有讲究的。滤波沿 u 方向（第 2 维）做，一维 FFT 连续作用于第 2 维的所有行；反投影按视角逐个取二维切片Pf(:,:,iv)。如果排成[nView nV nU]，每次取切片都要跨大步长，缓存命中率下来，Matlab 里能差出两三倍时间。

精度上建议全程single。512×512×360 的投影，single约 377 MB，double直接 754 MB；反投影阶段还要再开体数据、坐标网格、插值结果，很容易把 16 GB 的机器吃满。医学 CT 的衰减值动态范围本来就不需要双精度，single的 7 位有效数字完全够用，除非你要做定量分析并且对 1e-6 级别的误差敏感。

还有一个容易忽略的点：反投影的外层循环是 z 层，内层是视角。这样每个 z 层只需要一次性建好X, Y两张二维网格（meshgrid，尺寸 ny×nx），所有视角复用；如果把 z 放到内层，每个视角都要重建一张三维坐标网格，内存和计算量都翻几倍。这个顺序调整看起来只是两层 for 互换，实际跑 256³ 数据时差别是几分钟和几十分钟的差距。

3. 核心链路拆解：余弦加权、斜坡滤波、加权反投影的三段代码

FDK 的三段结构和二维滤波反投影几乎一样，多出来的只是锥角相关的两个修正：滤波前的余弦加权，和反投影时的1/U²距离权重。这两项是为锥束几何付出的代价——圆轨迹锥束不满足 Tuy 数据完备性条件，FDK 本质上是把倾斜射线近似成“在中心平面上”处理，锥角越大误差越大。理解这一点，后面调参和判断结果合不合理就有依据了。

3.1 余弦加权：把斜射线的路径长度补回来

平板探测器上，越靠边的像素，射线与中心射线的夹角越大，穿过的物质量也越多。余弦加权就是把这段多出来的路径折算回去：

function P = cosine_weight(P, cfg) % 输入输出 P: [nV nU nView] 的投影数据 u = ((1:cfg.nU) - (cfg.nU+1)/2) * cfg.du; % 以中心射线为 0 的坐标 v = ((1:cfg.nV) - (cfg.nV+1)/2) * cfg.dv; [U, V] = meshgrid(u, v); % [nV nU]，注意行列顺序 w = cfg.SDD ./ sqrt(cfg.SDD.^2 + U.^2 + V.^2); % 锥角的余弦 P = P .* reshape(w, cfg.nV, cfg.nU, 1); % 每个视角用同一张权重图 end

meshgrid(u, v)出来的U是 nV×nU，第 1 维随 v 变化、第 2 维随 u 变化，正好和[nV nU nView]的前两维对上。这一点必须确认：如果误用成ndgrid，权重图会转置，重建出来会出现沿对角线方向的渐变亮度，属于很难一眼看出的错误。

SDD放在分子上是个归一化选择，保证中心像素（u=v=0）权重为 1，这样加权不会改变投影的整体量级。有些实现写成SOD / sqrt(...)，只要滤波核的标定系数跟着调，结果一致；但两者混用就会让重建值的绝对大小差一个常数，做定量分析时会发现 HU 值系统性偏离。

3.2 斜坡滤波：Ram-Lak、Shepp-Logan 与窗函数

投影数据在角度方向上做反投影相当于积分，频域的 1/|ω| 会把低频放大到发散，所以要乘一个 |ω| 把 1/|ω| 抵掉，这就是斜坡滤波器。它同时放大了高频噪声，于是有了各种窗函数做折中。

function H = ramp_kernel(nPad, du, winName) % 返回已在 FFT 顺序上的频域滤波核，长度 nPad f = (-nPad/2 : nPad/2-1).' / nPad; % 归一化频率，单位 cycle/pixel H = abs(f) / du; % 斜坡 |ω|，量纲折进 du switch lower(winName) case 'ramlak' % 纯斜坡：最锐利，噪声和伪影也最猛 case 'shepp' H = H .* sinc(f); % Shepp-Logan：高频缓降，折中 case 'hann' H = H .* (0.5 + 0.5*cos(2*pi*f)); % 汉宁窗：平滑，边缘偏软 otherwise error('未知窗函数: %s', winName); end end

f用(-nPad/2 : nPad/2-1)/nPad生成，天然就是fft的输出顺序，省掉一次ifftshift，也就少了一个“shift 方向搞反”的隐患。零频落在nPad/2+1处，abs(f)在那里恰好为 0，斜坡滤波器的直流增益本来就是 0，不需要额外置零。

du折进核里是一种工程折中：连续卷积离散化后，核的幅值和采样间隔绑定。好处是换探测器时只改配置；代价是重建值的绝对量级会随du变化，做定量重建时要用已知密度的模体标定一个全局系数。我一般把这个系数单独存成cfg.calib，而不是散在代码里。

滤波本身要逐视角做，千万别一次性对[nV nPad nView]做 FFT。512×2048×360 的复数双精度数组是 6 GB 量级，多数机器会直接卡死：

function Pf = ramp_filter(P, cfg, winName) nPad = 2^nextpow2(2*cfg.nU); % 至少两倍补零，压住循环卷积 H = ramp_kernel(nPad, cfg.du, winName); Pf = zeros(size(P), 'like', P); buf = zeros(cfg.nV, nPad); % 单视角缓冲区，循环复用 for iv = 1:cfg.nView buf(:) = 0; % 补零段清零，避免上一轮残留 buf(:, 1:cfg.nU) = P(:, :, iv); F = fft(buf, [], 2); tmp = real(ifft(F .* H.', [], 2)); Pf(:, :, iv) = tmp(:, 1:cfg.nU); % 裁掉补零段 end end

补零到两倍以上是必须的。斜坡核在空间域是无限长、衰减很慢的，不补零就等于做了循环卷积，投影两端的数据会“绕”到另一端，重建出来在视野边缘出现一条弧形的亮带。补到2^nextpow2(2*nU)是个经验值，再大收益有限。

H.'这个转置容易写错。H是列向量，F是 nV×nPad 的矩阵，要按行广播必须转成行向量。写成H直接乘会触发隐式扩展报错，写成'在实数下没问题，但一旦哪天改成复数核就会变成共轭，建议统一用.'。

3.3 加权反投影：U² 权重和双线性插值

反投影是整条链路里最耗时的部分，也是精度损失最大的部分。每个体素要沿所有视角把滤波后的投影值“取”回来，取的坐标是小数，必须插值：

function vol = backproject(Pf, cfg) xu = ((1:cfg.nxyz(1)) - (cfg.nxyz(1)+1)/2) * cfg.dxyz; yu = ((1:cfg.nxyz(2)) - (cfg.nxyz(2)+1)/2) * cfg.dxyz; zu = ((1:cfg.nxyz(3)) - (cfg.nxyz(3)+1)/2) * cfg.dxyz; [X, Y] = meshgrid(xu, yu); % 固定一层，所有视角复用 uAxis = ((1:cfg.nU) - (cfg.nU+1)/2) * cfg.du; vAxis = ((1:cfg.nV) - (cfg.nV+1)/2) * cfg.dv; dbeta = cfg.Range / cfg.nView; vol = zeros(cfg.nxyz, 'single'); for iz = 1:cfg.nxyz(3) z = zu(iz); acc = zeros(cfg.nxyz(1:2), 'single'); for iv = 1:cfg.nView b = (iv-1) * dbeta; cb = cos(b); sb = sin(b); U = cfg.SOD - X*cb - Y*sb; % 中心射线方向上的等效距离 uq = cfg.SDD * (-X*sb + Y*cb) ./ U; vq = cfg.SDD * z ./ U; val = interp2(uAxis, vAxis, Pf(:,:,iv), uq, vq, 'linear', 0); acc = acc + (cfg.SOD ./ U).^2 .* val; % 距离平方反比权重 end vol(:,:,iz) = acc * dbeta; % 角采样间隔 end end

interp2的最后两个参数是外插值 0，超出探测器范围的点补零。这一步不能省：体素网格是立方体，探测器是矩形，锥束投影在边界处一定会有体素落在探测器之外，不给外插值会返回 NaN，一个 NaN 就能污染整层切片。

(SOD./U).^2是 FDK 里最容易被忽略的一项。距离源近的体素，同样的投影值对应更小的实际体积，反投影时必须按距离平方反比加权，否则重建结果会呈现以旋转中心为原点的径向渐变——中心偏暗或者偏亮，取决于权重是乘还是除。这条现象和“余弦加权漏了”很像，区分方法很简单：漏余弦加权的图，误差随 u 和 v 增大，在探测器边角最明显；错U²权重的图，误差随体素到旋转轴的距离变化，是同心圆状的。

最后的* dbeta是角度积分转求和的系数。用2π/nView还是π/nView取决于扫描范围：全扫描 360° 用2π，半扫描 180° 加扇形加权用π。这一项错了，重建值的绝对大小会差一个整数倍，但图像形状完全正常，属于“看着对、测着不对”的典型。

3.4 正向投影：用椭球体走一遍闭环

验证环节我一般先造数据再做重建，不依赖外部文件。体素驱动的前投影写起来最短，把每个体素按双线性“泼”到探测器上：

function proj = forward_project(vol, cfg) % 体素驱动前投影：逐层把体素累加到探测器 bin proj = zeros(cfg.nV, cfg.nU, cfg.nView, 'single'); for iz = 1:size(vol, 3) z = ((iz) - (size(vol,3)+1)/2) * cfg.dxyz; [vi, ui] = find(vol(:,:,iz)); % 只处理非零体素，稀疏更省 if isempty(vi), continue; end val = vol(sub2ind(size(vol), vi, ui, iz*ones(size(vi)))); x = (ui - (size(vol,1)+1)/2) * cfg.dxyz; y = (vi - (size(vol,2)+1)/2) * cfg.dxyz; for iv = 1:cfg.nView b = (iv-1) * cfg.Range / cfg.nView; cb = cos(b); sb = sin(b); U = cfg.SOD - x*cb - y*sb; iu = cfg.SDD * (-x*sb + y*cb) ./ U / cfg.du + (cfg.nU+1)/2; jv = cfg.SDD * z ./ U / cfg.dv + (cfg.nV+1)/2; proj(:,:,iv) = deposit(proj(:,:,iv), iu, jv, val); % 双线性累加 end end end

find那两行只取非零体素，是因为模体通常很稀疏：一个 256³ 的椭球模体，非零体素不到总量的 5%，跳过背景能省下大部分时间。deposit是双线性累加，把每个体素的值按距离权重分到四个相邻探测器 bin 上，和interp2的线性插值互为伴随操作，这也是它能和 FDK 闭环对上的原因。

要说清楚一点：体素驱动前投影本身有模糊，体素尺寸和探测器像元尺寸不匹配时更明显。用它造的投影重建出来，边缘会额外发虚，这部分不是 FDK 的问题。要做精度对比，换成解析投影（对椭球求射线与曲面的交点算弦长）更干净，代价是代码长一些，而且只对解析形状有效。我在做算法对比时会用解析版，做流程联调时用体素版，够快也够用。

4. 跑通一次完整重建：环境、参数表与分步执行

代码读懂了不等于跑得起来。这一章按“装什么 → 配什么 → 看什么”的顺序走一遍，中间结果该检查哪几项也一并说清，省得跑完二十分钟才发现参数填错。

4.1 MATLAB 版本与依赖

基础链路只依赖 MATLAB 本体，fft、interp2、meshgrid都是核心函数，不需要任何工具箱。如果要用现成的 Shepp-Logan 模体做验证，会用到 Image Processing Toolbox 里的phantom，但那只是二维版本，三维椭球模体自己写十几行更省事。

版本上建议 R2019b 及以上。原因有两个：一是zeros(..., 'like', P)这种类型跟随的写法在旧版本支持不全；二是隐式扩展（不用bsxfun）从 R2016b 才开始，代码里虽然已经避开了，但新版跑起来更省心。

内存是硬门槛。按 512×512×360 的投影、256³ 的重建算，峰值占用大概这样：

峰值在 1 GB 以内，16 GB 的机器跑起来很宽松。真正吃内存的是把single写成double，或者把滤波改成整块 FFT，那样轻松上到 6 GB。

4.2 参数表与配置文件写法

几何参数建议集中成一个结构体，入口脚本只负责填值，主流程只读不写。这样换数据时改一个文件就够了：

% main_fdk.m —— 一次完整重建的入口 clear; clc; cfg = default_cfg(); cfg.SOD = 500; % 源到旋转轴 mm cfg.SDD = 1000; % 源到探测器 mm cfg.du = 0.4; % 探测器列像元 mm cfg.dv = 0.4; % 探测器行像元 mm cfg.nU = 512; % 探测器列数 cfg.nV = 512; % 探测器行数 cfg.nView = 360; % 一圈投影张数 cfg.Range = 2*pi; % 全扫描；半扫描改 pi 并加 Parker 权重 cfg.nxyz = [256 256 256]; % 重建矩阵 cfg.dxyz = 0.5; % 体素边长 mm，一般取 du/2 ~ du 之间 S = load('data/proj_shepp.mat'); % 变量名 proj，[nV nU nView] proj = single(S.proj); proj = cosine_weight(proj, cfg); Pf = ramp_filter(proj, cfg, 'shepp'); vol = backproject(Pf, cfg); save('vol_fdk.mat', 'vol', '-v7.3');

dxyz怎么取是个实际问题。体素边长取du/2是最常见的做法，重建矩阵随之变成探测器尺寸的两倍，看起来“更清晰”，但那是插值出来的，不是真的分辨率。探测器像元 0.4 mm、放大率约 2（SDD/SOD）时，等中心处的等效探测器采样约 0.2 mm，dxyz = 0.2~0.3 mm是合理区间；再细就是给自己添内存。

Range这一项要特别小心。全扫描用2π，但很多实测数据是“短扫描”，转了 180° 加扇角，此时直接用π做反投影会出现严重的条状伪影，必须给每个视角乘 Parker 权重做冗余数据归一化。这是半扫描的必修课，代码里如果没有对应开关，宁可先按全扫描处理。

4.3 分步执行与中间结果检查

分三步跑，每步都看一眼中间结果，比一口气跑完再调试省时间：

第一步，只做余弦加权，看投影的数值分布。中心区域的均值应该基本不变，四个角的数值应该略微下降。如果角上反而变大，说明权重取成了倒数。

第二步，跑完滤波，看单个视角的滤波后投影。斜坡滤波会让数据出现负值，这是正常的——投影的高频被增强，边界附近必然过冲到负区。如果滤波后全为正，八成是滤波核没生效（比如H全是零，或者fft的方向维写成了 3）。

% 快速自检：滤波前后的值域对比 fprintf('raw : [%.4f, %.4f]\n', min(proj(:)), max(proj(:))); fprintf('filt : [%.4f, %.4f]\n', min(Pf(:)), max(Pf(:)));

第三步，跑完反投影，看一张中心切片。用imagesc(vol(:,:,128))加colormap gray、axis image，正常结果应该是背景接近均匀的暗区、模体边界清楚、没有明显的放射状条纹。再过一遍三维的vol(:)直方图，背景应该集中在 0 附近，形成一个尖峰，模体部分是一个稍宽的峰。如果背景峰的宽度和模体峰的宽度差不多，多半是角采样不够或者滤波窗开得太宽。

时间上有个感性的预期：256³ 体素、360 个视角，纯 Matlab 循环的反投影在普通笔记本上大概几分钟到十几分钟，取决于interp2的调用次数。如果跑了一小时还没完，先检查是不是把nxyz写成了[512 512 512]又用了double。

5. 避坑与排查：重建糊、环状伪影、内存爆掉这几类问题

5.1 图像质量问题：三个最常见的现象

现象一：整体发灰，中心亮边缘暗，像套了一层径向渐变。

原因基本锁定在距离权重上。(SOD./U).^2这一项要么漏了，要么写成了除法，要么U的定义和源的坐标不匹配（比如源写成了(R·sinβ, -R·cosβ, 0)，U还是R - x·cosβ - y·sinβ，符号就翻了）。解决方法是拿一个已知均匀的圆柱模体做单视角测试：只反投影一个视角，看图像是不是应该近似均匀。如果单视角就不均匀，问题一定在几何或权重，跟滤波无关。

现象二：模体沿某个方向拖尾、重影，像被抹了一把。

这类问题的根源几乎都是角度方向和探测器 u 轴方向不自洽。具体说，如果源按逆时针转，而 u 轴按探测器正面朝向来定义，那么 u 的正方向在实验室坐标系里是顺时针的切向，两者符号相反。表现就是重建图沿旋转方向被“拉”开，越靠视野边缘越明显。解决办法是把 u 取反再重建一次，两次结果里只有一次是锐利的。这个坑我用“玄学”形容过很多次——因为它在小视野、高对称模体上几乎看不出来，换到真实数据才暴露。

现象三：切片上出现一圈圈同心亮环。

环状伪影在 CT 里指向性很强：通常是探测器通道增益不一致或者有坏点。锥束平板探测器几千个通道，增益差异千分之几就足以形成可见的环。解决分两步：先做暗场和亮场归一化（(raw - dark) ./ (flat - dark)），再对明显偏离邻域的通道做插值修补。代码层面不要试图用滤波窗去压环状伪影，压不住，只会把图像一起糊掉。

现象四：锥角大的数据，越靠 z 两端越糊。

这不是 bug，是 FDK 的固有近似误差。圆轨迹锥束的数据不完备，FDK 把倾斜射线当成中心平面射线处理，锥角越大偏差越大。表现是远离中心平面的切片分辨率下降、出现纵向的模糊和伪影。可行的缓解手段是缩小锥角（增大 SOD 或者减小探测器 v 向覆盖），或者换迭代重建。我见过有人反复调滤波窗想解决这个问题，纯属白费功夫。

5.2 内存与性能问题：两条血泪经验

现象：跑到滤波那一步 MATLAB 直接卡死或者报 out of memory。

原因是把整个[nV nPad nView]数组一次性做了fft。512×2048×360 的复数双精度是 6 GB 量级，加上ifft时的临时副本，很容易翻倍。解决就是前面代码里的写法：开一个[nV nPad]的缓冲区，逐视角做完再写回。内存从 GB 级降到几十 MB，速度反而更快，因为大数组的缓存命中率本来就差。

现象：反投影慢得离谱，任务管理器显示内存占用一路飙升。

多半是在循环里反复创建全尺寸坐标网格。ndgrid出 256³ 的三张double网格就是 400 MB，放在视角循环里每轮新建一次，GC 根本跟不上。解决是坐标网格提到最外层、只用二维meshgrid、并且全程single。另外，如果只是想先看效果，把nxyz降到[128 128 128]跑一遍，确认几何没问题再上全尺寸，这个习惯能省下大量等待。

现象：并行化开了以后结果反而不对。

parfor里对acc这类累加变量做归约，Matlab 会自动处理，但如果切片里出现了对vol(:,:,iz)的写操作又没有被正确识别成切片输出，就会出现数据竞争。稳妥的写法是把每个 z 层的结果先存到 cell 或者临时三维数组里，parfor结束后再拼。这个坑排查起来很费劲，因为错误是概率性的，跑十次对九次。

6. 二次开发：把 FDK 搬到 GPU、加迭代、接深度学习的改造路径

FDK 在工程上最有价值的定位是“基线”。它快、确定性强、没有超参数，只要几何对，出来的图就一定是那个样子。所以二次开发通常不是改 FDK 本身，而是在它前后做文章。这里给三条我实际走过的路径。

第一条是把反投影搬上 GPU。这段代码天生适合并行：每个 z 层的每个体素互相独立，只是插值查表。改成gpuArray后，把X, Y, Pf都放到显存，interp2在 GPU 上有对应实现，256³×360 的反投影能从几分钟压到几十秒：

function vol = backproject_gpu(Pf, cfg) X = gpuArray(single(X)); Y = gpuArray(single(Y)); Pf = gpuArray(Pf); % 投影数据一次性上显存，约 377 MB uAxis = gpuArray(single(uAxis)); vAxis = gpuArray(single(vAxis)); % 循环体和 CPU 版完全一致，只把单精度类型统一 vol = gather(vol); % 只在最后取回主机内存 end

要注意的是显

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