CT重建中的等角扇束到平行束重排：MATLAB实现与验证

开场我直接说结论：在CT重建里，从扇束到平行束的重排，是个最容易被跳过、但对后续重建质量影响极大的环节。很多人学CT重建时，一上来就背FDK公式、背滤波反投影，等到自己动手做课题时才发现，手里采集到的原始数据几乎都是扇束几何。硬着头皮套平行束FBP，出来的图像要么有伪影，要么分辨率损失，最后只能回头补课。这篇文章想用MATLAB把等角重排算法完整走一遍，从几何映射推导到代码实现，再到验证和踩坑，争取让看完的人能直接把这套流程搬到自己的仿真或真实数据集上。

1. 为什么CT重建要绕一圈从扇束变平行束

1.1 中心切片定理才是FBP的底气

平行束重建之所有那么多人用，根本原因是中心切片定理存在。这个定理做的事情非常优雅：物体某一角度下的平行投影，其傅里叶变换结果等于物体二维傅里叶变换中过原点且方向一致的那条直线切片。也就是说，只要采集足够多角度的平行投影，就能在频域里把物体整张频谱图铺满，再逆变换回来就是重建图像。

1D投影 → 1D傅里叶变换 → 填充2D频域 → 2D逆傅里叶变换，这条链路每一步都有严格数学支撑，滤波反投影也可以看作在这个框架下推导出来的实际操作形式。

然而真实CT系统几乎不会用平行束扫描。X射线源是一个点源，发出的射线天然呈扇形展开，配合对面弧形或平板探测器，构成扇束几何。不做任何处理、直接把扇形射线的投影数据当成平行束来重建，每条射线到旋转中心的距离和角度全部对不上，重建结果必然是模糊加伪影双双拉满。所以想用成熟的平行束FBP管线，第一步就是做几何重排。

1.2 重排的本质是数据重新装箱

重排算法说白了就是：把扇束几何里每一条射线，按它真实的空间位置关系，换算成"这条射线如果放在平行束坐标系里，应该属于哪个投影角度、哪个径向位置"，再对应到新的数据矩阵里。如图1，概念上完全可以想象成在同一个扫描空间里，把射线一根一根挑出来重新排队。几何关系一旦确定，剩下就是数据映射和插值。

等角扇形重排（Equiangular Rebinning）特指扇束中射线角度等间隔的情况。做CT的人接触到两类扇束：一类是等角扇束（第γ方向采样间隔均匀），一类是等间距扇束（探测器线性排列，射线在探测器上等间距分布）。前者多见于早期CT，后者在现代平板探测器CT/工业CT中更常见。不同几何对应不同的重排推导过程，本文专注等角扇束，对等间距扇束的差异会在第5节里单独说明。

2. 等角扇束重排的几何映射：不敢画图也能懂的角度换算

2.1 先把坐标系和关键参数定下来

要把重排讲清楚，参数必须先统一。本文采用如下约定（这个约定在MATLAB仿真里对应得很好）：

在源角β、扇形角γ时，这条射线在空间中的位置如图2所示（这段描述你照着参数脑补即可：源的极角是β，射线相对扇形中心线偏了γ角，中心线与源位置-旋转中心连线重合）。这条射线和旋转中心的垂直距离，几何上恰好是 R·sinγ，它在平行束坐标系里对应的投影角度是β + γ。所以映射关系非常干净：

[ t = R \sin\gamma ] [ \theta = \beta + \gamma ]

注意这个关系是完全几何等价的，不涉及任何近似。反向推导也顺手：

[ \gamma = \arcsin\left(\frac{t}{R}\right) ] [ \beta = \theta - \gamma ]

这正是等角重排的正面战场：在已知扇形数据和各采样网格的前提下，把目标平行束网格（θ, t）对应的(β, γ)算出来，再从扇束sinogram里把值取出来填回去。

2.2 为什么等角重排的t方向天生不均匀

等角扇束中γ等间隔分布，但t和γ的关系不是线性的，而是正弦关系。这意味着如果直接把每个角度β下的射线数据按t填进平行束网格，靠近中心的位置（γ接近0）密度高，往两侧（γ接近±γ_max）密度变稀疏。这就带来一个潜在问题：如果目标平行束要求在t方向均匀采样，采样点对应的γ反算回去并不是等间隔的，所以必须做插值，而不是一对一地搬数据。

反过来说，目标平行束网格的t范围也受限。因为γ最大只有γ_max，所以t的取值范围理论上限是R·sinγ_max。超过这个半径，扇形几何根本没有射线能覆盖，重排结果也必然是空的，这点在设置输出网格边界时必须卡死。

2.3 重排数据的角度覆盖范围要比π大一圈

平行束FBP的投影角度只需要0到π，但扇束重排之后，目标平行束的有效角度覆盖不是正好π。由于θ = β + γ，扇束在β从0到2π范围内提供的θ范围是[-γ_max, 2π+γ_max]，这比理想平行束的需求宽得多，多出来的部分包含冗余信息。

实际操作中，通常把重排输出目标θ放在0到π+2γ_max范围内，其中0到π这一段用于FBP重建。直接取0到π也凑合，但会损失两侧的小角度数据，尤其在γ_max较大时，边缘位置的数据会被硬生生截掉，带来重建边缘伪影。仿真如果是从零模拟投影数据再重排，那就更不应该截断了，要做到能利用的全利用。

3. MATLAB实现重排：从sinogram到平行束数据

3.1 先准备一套可以自证正确性的仿真数据

为了验证等角重排算法的正确性，直接用真实CT数据是最麻烦的，因为真实数据没有标准平行束投影可以对照。更稳的做法是：先造一个已知的体模（比如经典的Shepp-Logan），分别用解析的方法计算它的扇束投影和平行束投影。扇束投影进重排算法，平行束投影当作标准答案，两者对比就知道重排做得对不对。

我用一个简化的数值体模来实现。真正的Shepp-Logan需要用椭圆叠加解析投影公式，写起来略长。这里把思路演示清楚即可——先定义一组椭圆参数，然后写一个通用函数，输入坐标即可得到体模的吸收系数；投影时沿着射线方向数值积分。这样做的好处是扇束投影和平行束投影可以共用一套体模API，不用维护两份几何代码。

3.2 扇束投影模拟函数

模拟扇束投影时，遍历每一个源角β和每一个扇形角γ，确定射线的起点（源位置）和方向，然后从射线起点到探测器位置做等间隔采样，累加体模衰减系数乘以步长，得到投影值。核心代码如下，简化理解：

function sino = fanbeam_forward(model, R, D, beta_list, gamma_list, sample_n) % 等角扇束前向投影数值仿真 % model: 体模结构体 % R: 源到旋转中心距离 % D: 旋转中心到探测器距离 % beta_list: 源角度列表(rad)，长度N_beta % gamma_list: 扇形角列表(rad)，长度N_gamma % sample_n: 每条射线采样点数 N_beta = length(beta_list); N_gamma = length(gamma_list); sino = zeros(N_beta, N_gamma); for ib = 1:N_beta bx = R * cos(beta_list(ib)); by = R * sin(beta_list(ib)); for ig = 1:N_gamma % 射线方向: 从源指向探测器上某点 gamma = gamma_list(ig); dir_x = -cos(beta_list(ib) + gamma); dir_y = -sin(beta_list(ib) + gamma); % 沿射线方向采样体模衰减系数并积分 t_vals = linspace(0, D + sqrt(R^2 - (R*sin(gamma))^2), sample_n); px = bx + dir_x * t_vals; py = by + dir_y * t_vals; vals = arrayfun(@(px,py) phantom_value(model, px, py), px, py); sino(ib, ig) = trapz(t_vals, vals); end end end

沿射线方向采样的时候，起始点就是源位置，终点应该越过旋转中心和探测器平面，这样积分才能把射线路径上的体模覆盖完。上面代码里终点的计算略粗糙，真实实现里应精确算出射线与探测器平面的交点。这里说一点，这种逐射线积分的做法速度一般，只适合小尺寸验证。你也完全可以换用radon变换配合几何变化来模拟，但几何直观性会差一些。

3.3 平行束投影模拟函数

平行束前向投影就简单了。给定投影角度θ，平行射线族的法向量方向就是θ方向，沿法向量方向的偏移距离是t。对每条射线，同样沿射线方向积分体模衰减系数：

function sino_para = parallelbeam_forward(model, theta_list, t_list) % 平行束前向投影数值仿真 % theta_list: 投影角度列表 % t_list: 径向偏移距离列表 N_theta = length(theta_list); N_t = length(t_list); sino_para = zeros(N_theta, N_t); for it = 1:N_theta theta = theta_list(it); % 射线单位方向向量, 与法向量垂直 unit_dir = [-sin(theta); cos(theta)]; unit_norm = [cos(theta); sin(theta)]; for it_ = 1:N_t t = t_list(it_); % 射线过点: t*unit_norm + s*unit_dir % 沿射线采样, s范围取足够大(-1.5~1.5), 体模半径控制在1.0内 s_vals = linspace(-2, 2, sample_n); px = t * unit_norm(1) + s_vals * unit_dir(1); py = t * unit_norm(2) + s_vals * unit_dir(2); vals = arrayfun(@(px,py) phantom_value(model, px, py), px, py); sino_para(it, it_) = trapz(s_vals, vals); end end end

这段代码虽然是双循环，但好处是每一步对应关系都很清楚，出了bug容易查。

3.4 核心：等角重排函数的实现

重排函数是核心。输入为扇束数据和对应的β、γ坐标网格，输出为目标θ、t网格上的平行束投影数据。实现里我采用双循环遍历目标平行束网格点，对每个点反算(β, γ)，再从扇束数据插值。这种实现方式速度不是最优，但逻辑直白、可读性强，非常适合理解和调试。等你确认算法无误，再优化成向量化或预先计算索引矩阵也不迟。

function para_sino = fan2para_rebinning(fan_sino, beta_list, gamma_list, R, theta_list, t_list) % 等角扇束到平行束重排 % fan_sino: N_beta x N_gamma 的扇束投影数据 % beta_list: 源角度(rad) % gamma_list: 扇形角(rad) % R: 源到旋转中心距离 % theta_list: 目标平行束投影角度(rad) % t_list: 目标平行束径向偏移 N_theta = length(theta_list); N_t = length(t_list); N_beta = length(beta_list); N_gamma = length(gamma_list); % 构造二维插值所需的网格. 注意beta/gamma都需要转成矩阵网格形式 [Beta_grid, Gamma_grid] = meshgrid(beta_list, gamma_list); Beta_grid = Beta_grid'; % N_beta x N_gamma Gamma_grid = Gamma_grid'; para_sino = zeros(N_theta, N_t); for it = 1:N_theta theta = theta_list(it); for jt = 1:N_t t = t_list(jt); % 反向几何计算 gamma_target = asin(clamp(t / R, -1, 1)); beta_target = theta - gamma_target; % 角度的周期延拓处理: beta需要落在[0, 2*pi)区间 beta_target = mod(beta_target, 2*pi); % 插值: 线性插值, 越界取NaN val = interp2(Gamma_grid, Beta_grid, fan_sino, ... gamma_target, beta_target, 'linear', NaN); para_sino(it, jt) = val; end end end

这段代码里有两个特别容易出错的点。第一个是meshgrid出来的维度方向问题，interp2第一个参数是要查询的x坐标网格，第二个是y坐标网格，顺序反了结果直接全NaN，而且非常难排查。第二个是角度周期延拓，mod操作确保超出一个扫描周期的β映射回有效区间，但这种操作只对完整360°扫描成立。如果扫描范围本身小于2π，比如某些快速扫描只有200°甚至180°+扇角，那重排出来的平行束数据会有大块空洞，这时候就不是插值能解决的了，必须从扫描协议层面保证覆盖度。

3.5 向量化版本：大项目跑得动才有意义

fan2para_rebinning对每个目标网格点做一次插值，Θ×T如果是400×401，循环次数是16万次，每次调用interp2还有不小的开销。在实际项目里，这种实现跑起来人会等疯掉。好在重排的插值逻辑完全是逐点独立的，天然适合向量化或预计算。

预计算的思路：既然β_target和γ_target只依赖于目标网格(θ, t)和固定的R，那就在重排之前把所有网格点对应的β、γ矩阵算出来，一次性传给interp2：

function para_sino = fan2para_rebinning_vec(fan_sino, beta_list, gamma_list, R, theta_list, t_list) [Theta_grid, T_grid] = meshgrid(theta_list, t_list); % Theta_grid, T_grid(N_t, N_theta), 需要转置成与输出维度一致 theta_flat = Theta_grid(:)'; t_flat = T_grid(:)'; gamma_target = asin(clamp(t_flat / R, -1, 1)); beta_target = mod(theta_flat - gamma_target, 2*pi); % 用griddata或者interp2插值, interp2要求输入网格单调, 这里直接用scattered插值 para_sino = griddata(beta_list, gamma_list, fan_sino', ... beta_target, gamma_target, 'linear'); para_sino = reshape(para_sino, N_t, N_theta)'; end

griddata在速度上不如interp2经过优化的网格插值，但它接受散点坐标查询。如果你的数据量不太大，这是最省事的向量化写法。更大规模场景建议用interp2，先把β_grid、γ_grid的meshgrid存下来，再一次传递所有查询点。这个优化可以留到算法跑通后再做，前期求稳，后期求快。

4. 验证重排结果：用Shepp-Logan体模当裁判

4.1 实验参数设计

我采用如下参数做验证：

• 模型：修改过的Shepp-Logan（尺寸归一化到半径1）
• R = 3，D = 3
• γ_max = π/6（即30°扇形角范围）
• 扇束投影：β从0到2π共720个角度，γ从-30°到30°共201条射线
• 平行束投影：θ从0到π共360个角度（重排后用），t范围[-1.5, 1.5]共301个径向位置

理论上重排出来的平行束数据θ覆盖范围为γ_min~π+γ_max，我这里取0到π，中间大部分均有效。为了对比，先用解析法算一遍真平行束投影，再用扇束投影重排出平行束投影。

4.2 对比维度一：sinogram可视化

把真平行束投影和重排平行束投影并排显示，第一眼就能看出结构是否一致。平行束sinogram的形状是"中央横向亮带、两端弯月形暗区"，重排结果如果正确，会和真平行束投影有高度相似的纹理结构。如果角度对应错了，sinogram会沿轴向错位；如果t方向对应错了，横向结构会变形。

4.3 对比维度二：像素级差值与剖面线

可视化只能算定性判断，定量还是要看数值。将两组sinogram的相对误差画出来，中心区域误差应当在1%以内，边缘由于插值精度和截断效应会略差一些。比较时要注意，真平行束投影的角度采样范围是θ∈[0, π)，而重排数据有效θ范围其实比这个大。对比时要确保两者用的是同一θ子区间，比如都切成θ∈[10°, 170°]来比，避开边缘处受边界效应影响大的部分。

更细的做法是在某个固定θ角度画一条t方向剖面曲线，把真值和重排值叠在一起，二者应当高度重合。剖面线可以直观看出插值有没有带来系统性的偏移或振荡——如果曲线出现齿状波纹，基本可以断定插值选择或γ_max边界裁剪有问题。

4.4 对比维度三：重建图像对比（终极大考）

sinogram数值接近还不够，重排的最终目的毕竟是重建。把真平行束投影和重排投影分别送入FBP重建，重建网格256×256，结果放在一起看：

• 真平行束重建：作为参考，图像质量最好
• 重排后平行束重建：图像应当与参考几乎一致，边缘可看到轻微模糊（因为插值有平滑效应）
• 未经重排直接扇束FBP重建：拿来当"反面教材"，图像有明显变形或伪影，用于反衬重排的必要性

如果两种平行束结果差异过大，但有使用原始数据重建校验，问题大概率出在γ_target和β_target换算出错，或者插值方向用反了。

5. 等角重排实战中的那些坑

5.1 β角度循环边界：mod之后的世界并不完美

扇束扫描β是从0转到2π，重排目标β_target算出来可能落到负值或超过2π。用mod(beta_target, 2π)处理后，逻辑上确实对应同一个物理位置，但插值时要把扇束数据当成循环数组处理。interp2本身不知道β维是循环的，在0附近和2π附近插值时会出问题：如果查询点在2π-0.01，插值函数会在β≈2π的邻近网格点找数据，但2π和0是相邻的，数据却存到了数组两端，插值结果自然不对。

有两个办法：一是把扇束sinogram沿β维度手动首尾拼接，复制一份开场角度的数据接到结尾，扩大插值范围；二是在interp2之前把β_target用线性映射的方式做一次平滑偏移处理。第一种最常用，实现简单，代价是内存翻倍但一般完全可以接受。

5.2 γ_max截断：扇形覆盖范围之外的t要怎么处理

等角重排有个硬边界：t_max = R·sinγ_max。目标平行束网格里|t|超过这个值的点，几何上根本不成立，反算γ时会超过arcsin定义域。如果不对t_list做限制，asin传参会出NaN，后续插值全被污染。正确做法是重排前先计算t_max，把目标t_list裁剪到[-t_max×0.98, +t_max×0.98]，留出插值边缘余量。这同样意味着如果你想要更大的平行束覆盖面积，就得在扫描设计时提高γ_max或者增大R，不能指望重排算法凭空变出数据。

5.3 插值方法：线性还是三次，真的要考虑场合

双线性插值速度最快，对数值积分得到的投影数据来说，精度其实已经足够。因为CT投影数据本身高频成分有限，重排又是做数据重采样，不是最终成像，引入的微小平滑会被后续FBP滤波器放大，但不至于危害图像质量。三次插值（spline或cubic）减少平滑效果，但会产生过冲振荡，对边缘锋利的体模可能引入振铃伪影。实测下来，线性插值是常用且稳妥的选择，除非有特殊精度要求，不建议轻易上cubic。

5.4 等角重排和等间距重排：千万别混用

开头提过还有一类等间距扇束几何（探测器线性排列）。等间距扇束中，t与探测器像素位置的关系不是R·sinγ，而是t = u·R/√(R²+u²)（u为探测器线性坐标），推导过程和等角情况完全不同。如果你拿到的原始数据是平板探测器扫描出来的，却套用等角重排的公式，误差会直接反映在重建图像的几何失真上，难以靠调参解决。判断几何类型最简单的方法是看投影数据的采样坐标：如果每角度下探测器维等间隔对应角度的sin/cos关系，很可能是等间距；如果探测器维等间隔对应等角度γ，那就是等角。动手前一定要确认清楚。

5.5 数据精度陷阱：单精度改双精度，重排结果焕然一新

这个点极其容易被忽视。工业CT或医疗CT导出的投影数据很多是16位整型，但进行角度换算时，如果沿用单精度浮点存储网格坐标，累计几个重排步骤之后精度损失会明显影响重建质量。在MATLAB里尤其要养成习惯：重排时用double全程计算，只在最后输出时才允许转回需要的类型。别问为什么你的重排老是出现横向条纹，可能就是精度不够。

5.6 角度采样密度：重排后的有效角度分辨率低于直觉

重排不是无损数据变换。扇束投影有Nβ个投影角、Nγ条射线，重排到平行束后，θ方向的有效采样数不是简单的Nβ，而是Nβ+Nγ；t方向的采样密度也受到γ采样间隔和R共同影响。如果目标平行束网格设置得比原始数据可支撑的密度更细，插值结果只是把原始信息平滑放大，不会带来更高分辨率。这点在设计成像协议时就要留意，网格密度是否可以降低。比如θ方向太密集反而拉长计算时间，还增加插值伪影。

6. 完整示例代码与使用说明

把上面所有内容整理成一个可直接运行的MATLAB示例脚本，建议直接复制保存，边看注释边跑。脚本结构为：生成体模 → 扇束投影 → 等角重排 → 与真平行束对比 → FBP重建。

%% 参数定义 clear; close all; clc; R = 3; % 源到旋转中心距离 D = 3; % 旋转中心到探测器中心距离 gamma_max = pi/6; % 最大扇形角 30度 N_beta = 720; % 源角度采样数 N_gamma = 201; % 每角度下射线数 N_theta = 360; % 目标平行束角度数 N_t = 301; % 目标平行束径向采样数 beta_list = linspace(0, 2*pi, N_beta); gamma_list = linspace(-gamma_max, gamma_max, N_gamma); theta_list = linspace(0, pi, N_theta); t_max = R * sin(gamma_max); % 理论最大可覆盖半径 t_list = linspace(-t_max*0.98, t_max*0.98, N_t); %% 体模生成 model = define_shepp_logan(); % 显示体模 [xg, yg] = meshgrid(linspace(-1.5, 1.5, 256)); img = arrayfun(@(x,y) phantom_value(model, x, y), xg, yg); figure; imshow(img, []); title('Shepp-Logan体模'); %% 扇束与平行束前向投影 fan_sino = fanbeam_forward(model, R, D, beta_list, gamma_list, 200); para_sino_ref = parallelbeam_forward(model, theta_list, t_list, 200); figure; subplot(121); imshow(fan_sino, []); title('扇束投影 sinogram'); subplot(122); imshow(para_sino_ref, []); title('真平行束投影 sinogram'); %% 等角重排 para_sino = fan2para_rebinning(fan_sino, beta_list, gamma_list, R, theta_list, t_list); figure; subplot(121); imshow(para_sino_ref, []); title('参考平行束'); subplot(122); imshow(para_sino, []); title('重排平行束'); % 相对误差 mask = para_sino_ref > max(para_sino_ref(:))*0.1; err = abs(para_sino - para_sino_ref) ./ max(abs(para_sino_ref(:)), 1e-8); fprintf('中心区域平均相对误差: %.4f%%\n', mean(err(mask(:)), 'all')*100); %% FBP重建对比 (用简单斜坡滤波) recon_ref = iradon(para_sino_ref, theta_list*180/pi, 'linear', 'Ram-Lak', 1, 256); recon_rebin = iradon(para_sino, theta_list*180/pi, 'linear', 'Ram-Lak', 1, 256); figure; subplot(121); imshow(recon_ref, []); title('参考平行束FBP重建'); subplot(122); imshow(recon_rebin, []); title('重排平行束FBP重建');

脚本跑完后如果一切正常，你会看到重排后的sinogram和参考版本视觉上几乎一样，相对误差集中在1%上下，重建图像之间的差异更是微乎其微。反之，如果发现横向错位、纵向拉伸或者大面积NaN，那就按第5节提到的几个坑逐一排查。

7. 从重排到实战：给研究者和工程师的几点建议

等角重排只是CT重建链路上的一环，但它决定了后续所有步骤的输入质量。在做科研或项目预研时，挺建议把这套流程作为标准动作固定下来。队里新同学上手时，先让他们跑通重排和验证流程，再碰实际采集数据，踩坑率会下降不少。

实际落地时，还有几个值得延伸的方向：

• 多排CT/锥束CT场景，等角重排的思想可以推广到锥束数据的重排，但要额外考虑层间z方向的权重和倾斜射线，复杂度和二维情况不在一个量级。
• 如果投影数据角度范围不足360°（比如短扫描），重排中利用冗余投影信息时就要引入加权策略，否则边界处会明显不连续。
• 在GPU加速环境里，重排算法的高度并行特性非常适合用CUDA实现，每个目标网格点的插值完全独立，映射表预计算之后直接用texture memory读取，能比CPU快一到两数量级。

我过去给一位做CT安全检查设备的同行调过类似的算法，当时发现他们的工业软件里重排过程用的是最近邻插值，图快，结果重建的边缘分辨率一直在指标线上徘徊，换成双线性之后整个图像质量立刻上了一个台阶。这种问题不到实际落地，很难从论文里体会到。

最后提一个小心得：无论算法写得多么花哨，永远先用仿真数据验证几何推导。等角重排这种纯几何算法，几何对了就基本对了，剩下的都是插值精度问题。任何时候发现重排结果不对，第一反应不要调插值参数，而是回去检查R、γ_max、β范围这些几何量有没有保持一致，十有八九问题出在几何上。