做CT成像仿真的朋友应该都有过这种经历手里的投影数据明明是扇束扫描出来的可最顺手的重建工具偏偏只认平行束一打开iradon的帮助文档满屏都是“parallel-beam”的字眼。我最初做实验时也是卡在这一步后来老老实实把“等角重排”算法自己实现了一遍才真正搞明白扇束几何与平行束几何之间那点微妙的关系。这篇博客就把整个推导过程、MATLAB代码和踩坑记录完整整理出来给正在做CT重建、医学物理课程作业或者想脱离商业软件自己搭重建流水线的同学一个能直接落地的参考。很多人会问扇束数据直接用扇束FBP公式不就行了吗为什么非要多一步重排其实这个问题就像“有高铁为什么还要普通列车”——扇束FBP当然能重建但平行束重建在理论推导、工具支持、算法验证上都更成熟而且把数据转换到平行束几何后你可以直接复用大量现成代码还能方便地对比各种重建算法之间的差异。等角重排这个算法本身也不算复杂但它背后涉及的坐标映射、插值策略、角度重合判断几乎是所有CT几何处理问题的缩影搞懂了它后面接触锥束重建、螺旋CT时很多思路都能迁移。1. 项目整体思路为什么“扇束转平行束”是CT重建的必备技能1.1 平行束与扇束两种扫描几何的本质差异先花两分钟把两种几何彻底说清楚。平行束几何可以想象成一台相机在很远很远的地方平移拍摄每次曝光都得到一组完全平行的X射线射线束的方向可以用一个单一的角度θ来描述每条射线到旋转中心的距离s是一个一维坐标。在这个模型下投影数据P(θ, s)天然就是一个规整的二维网格行是径向位置列是投影角度滤波反投影公式推导起来非常干净。扇束几何则完全不同。X射线源是一个点射线从点源发散出去在物体另一侧用探测器阵列接收。用旋转角β表示源在圆周上的位置每条射线与源到旋转中心连线的夹角记为γ。这时候投影数据变成P(β, γ)相当于每个β角度下都有一个关于γ的一维数组。从几何上讲平行束用的是“线”扇束用的是“束”一个是平行线族一个是点源放射线族。打个比方平行束像一群士兵排成整齐的方阵从正面推进扇束则像一盏探照灯从侧面扫射。士兵方阵的“入射方向”永远是同一个角度而探照灯的每一条光线的方向其实都不相同。CT重建算法最朴素的物理图像——把所有方向的投影叠回去——在平行束下最容易理解因为每个方向的“叠回去”就是一次标准的平移叠加换成扇束每条射线经过的路径和角度都不同反投影时还得给每条射线单独计算它在图像坐标系里的方向麻烦就多起来了。1.2 直接重建的困难与重排存在的意义理论上扇束有自己的FBP公式像MATLAB的ifanbeam就是专门干这个事的。但扇束FBP公式里多了一个和距离相关的加权因子出错时不好排查而且很多开源的CT重建工具包、科研代码都只支持平行束输入。更实际的一点是平行束重建算法在傅里叶切片定理的框架下特别优雅一个角度下的投影其傅里叶变换正好等于图像二维傅里叶变换经过原点的一根切片。这个性质在扇束下也有对应但解释起来绕很多。所以重排算法的核心价值就是把扇束投影数据通过插值重新组织成平行束投影数据让后面的重建可以完全复用成熟、简洁的平行束算法。在工业CT、医学CT的教学仿真和算法验证中这种“先转几何再走标准流程”的做法极其常见。它不是唯一的方案但一定是最容易理解、最容易调试的方案。1.3 等角重排与其他重排方式的区别重排算法按探测器几何可以分成两大类等角重排和等距重排。等角重排假设探测器是圆弧形的相邻探测通道之间对源的张角γ是恒定间隔等距重排则假设探测器是直线排列的相邻通道间距在探测器平面上均匀分布。这两者的区别直接影响映射公式的形式。等角几何下径向坐标s与探测器角度γ的关系是一个干净的D·sin(γ)推导和编程都简单等距几何下探测器坐标x_d与γ之间还隔了一个arctan关系公式多一层嵌套插值时的网格也不那么规则。我这次选等角重排一方面因为MATLAB的fanbeam函数在‘arc’几何下直接输出等角数据另一方面也因为等角公式在讲解坐标映射时最有代表性——扇束的“角度”属性是完全显式的重排的核心逻辑不会被额外的坐标变换干扰。2. 核心原理等角扇束到平行束的坐标映射推导2.1 变量约定β、γ与D到底代表什么在动手写代码之前必须把变量定义钉死否则后面所有公式都会乱成一团。我用的是这套约定DX射线源点到旋转中心的距离单位是像素β源在圆周上的旋转角也就是fanbeam返回的fanRotAngles单位通常为度γ探测器某个通道相对于中心射线的夹角中心射线指从源指向旋转中心的那条射线γ的单位也是度P(β, γ)扇束投影值表示源在β角度时与中心射线夹角为γ的那条射线的探测器读数。还有一个容易搞混的点中心射线本身的方向。假设源在β位置那么从源指向旋转中心的方向角其实是β加上180度。但如果只看“射线与物体几何的夹角”我们通常习惯把源到中心的连线视为参考轴γ是相对于这条轴的偏移量。这个参考轴选择很重要因为后面推导平行束投影角θ的时候其基准定义直接决定你重建出来的图像会不会旋转。2.2 一条射线的身份转换从(β, γ)到(θ, s)现在来推导核心公式。假设源位于β角度我们考察一条从源发出、与中心射线夹角为γ的射线。这条射线在扇束数据里由(β, γ)唯一标识。目标是找到它在平行束几何下的两个参数投影角θ和径向距离s。径向距离s的推导最经典。以旋转中心为原点建立直角坐标系源S坐标为(D·cosβ, D·sinβ)。射线从S射出方向与中心射线方向夹角为γ。这条射线是经过S点、方向向量为u的直线u可以写成(-cos(βγ), sin(βγ))这里γ的正负需要跟探测器排列的左右方向匹配但我们先按这个符号体系展开。直线到原点的有向距离s等于从原点到直线的垂足距离。用直线方程法射线方向向量u的垂直法向量n(sin(βγ), cos(βγ))过点S的直线方程可以写成n·((x,y)-S)0。把原点(0,0)代入左边得到代数距离s n·(-S) -D·sin(βγ)·cosβ - D·cos(βγ)·sinβ -D·sin(βγβ)... 这一步看似复杂但其实如果取源在β0的位置特殊化会更直观。我建议读者在推导时直接令β0验证源在(D,0)中心射线沿x轴向左γ是偏离角。此时射线方向为(-cosγ, sinγ)法向量为(sinγ, cosγ)直线方程为sinγ·(x-D)cosγ·y0即sinγ·xcosγ·yD·sinγ。原点(0,0)到该直线的有符号距离就是D·sinγ。这与β无关因为整体旋转不会改变点到直线距离的绝对值只会把角度基准旋转。因此最终公式s D·sin(γ)这个公式极其重要它说明扇束径向信息s只和探测器通道角度γ有关和源旋转角β无关。至于平行束投影角θ理论上θ应该描述射线方向在全局坐标系中的朝向。在β0、γ0的特殊情况中心射线沿x轴负方向对应平行束投影角度是0或180度取决于你以射线方向还是探测器法向为基准。不同教材、不同软件包的θ定义并不完全一致这也是我实际调试中踩坑最多的地方。普遍可行的做法是先用点源或非对称体模校准一次把重排结果直接扔给iradon看重建图像是旋转了、镜像了还是完全正常然后根据偏差调整θ公式里的常数项。我在代码里用的是一组自洽且经过点源验证的映射θ β γ - 90°角度制对应的平行束数据用iradon重建后能得到正确朝向。2.3 重排为什么需要反查与插值有人可能会想既然有了正向公式(β, γ) → (θ, s)那我遍历原始扇束数据的所有点把它们一个个填到平行束网格上不就行了理论与工程的分歧就在这里。原始扇束数据在(β, γ)网格上是均匀的但映射到(θ, s)空间后点的分布会变得不均匀——有的地方密集有的地方稀疏直接填进去会在平行束投影矩阵上留下大量空洞或重叠点。更稳妥的思路是“反查”先确定目标平行束投影矩阵的每个像素(θ_target, s_target)然后反求这应该对应原始扇束数据中的哪一点(β_source, γ_source)再从原始数据中插值读出数值。反查的好处是目标网格完全规整每条平行束投影线都能保证有一个输出值不会出现空洞。而反查过程中因为γ_sourceasin(s_target/D)往往不是探测器通道上的整数刻度β_source也未必正好落在旋转角采样点上所以必须用插值——实际中双线性插值在精度和计算量之间最平衡我用的就是interp2默认的线性插值。2.4 为什么用双线性插值而不用最近邻最近邻插值确实快但效果在CT里惨不忍睹。你可以想象把一张低分辨率照片放大到原尺寸最近邻会让每个像素的边界变成生硬的锯齿在重建图像上这些锯齿会演变成一圈圈的放射状伪影特别影响图像质量。双线性插值相当于把周围四个采样点的值按距离加权结果更平滑代价只是多几行代码和微不足道的计算时间。如果你的实验对精度要求极高还可以上三次插值但实测下来在D取得足够大、角度采样足够密的条件下双线性插值得到的重建图像和原始平行束直接重建的差距已经很小了。3. MATLAB实操从扇束投影到可用平行束投影的完整实现3.1 环境准备与测试数据生成MATLAB版本方面我用的是R2023a但下面代码用到的fanbeam、iradon、interp2都是Image Processing Toolbox里的老牌函数R2019b及之后版本基本都能直接跑。记得先确认一下自己的工具箱齐全在命令行输入ver看看有没有Image Processing Toolbox没有的话先补装。测试数据直接用Shepp-Logan体模这是CT重建领域最经典的合成图像内部有多个椭圆结构能很好地暴露伪影和形变问题。生成扇束投影代码clear; close all; clc; % 生成Shepp-Logan体模大小256x256 img phantom(256); % 设置扫描几何 D 500; % 源到旋转中心的距离单位像素 fanSpacing 0.25; % 等角探测器的通道角间距单位度 fanRotationInc 1; % 源旋转角步长单位度 fanAngles 0:fanRotationInc:359; % 完整360度扫描 % 生成扇束投影 [F, fanSensorPos, fanRotAngles] fanbeam(img, D, ... FanSensorGeometry, arc, ... FanSensorSpacing, fanSpacing, ... FanRotationIncrement, fanRotationInc, ... FanCoverage, cycle);细说一下fanbeam输出的三个量。F是扇束投影矩阵行对应探测器通道列对应旋转角F(i,j)表示第j个旋转角下第i个探测器通道的读数。fanSensorPos是每个探测器通道对应的角度位置在arc几何下它的单位直接就是度范围大概是[-20,20]这样的区间具体取决于D和体模尺寸。fanRotAngles就是每个旋转角的数值从0到359度。3.2 重排核心代码与关键步骤接下来是重排的核心部分。目标是把F矩阵重排成一个平行束投影矩阵projParallel行对应径向位置s列对应平行束投影角θ。% 将角度统一转为弧度避免反复换算 gammaVec fanSensorPos * pi/180; betaVec fanRotAngles * pi/180; % 平行束径向坐标范围 smax D * sin(max(abs(gammaVec))); M size(F, 1); % 探测器通道数 s linspace(-smax, smax, M); % 径向坐标采样 % 目标平行束投影角0~179度就足够常见FBP重建 theta_deg 0:1:179; theta theta_deg * pi/180; % 预分配输出矩阵 projParallel zeros(numel(s), numel(theta)); % 用ndgrid构造目标网格 [Theta, S] ndgrid(theta, s); % 核心反查映射gamma asin(s/D)beta theta - gamma pi/2 gammaTarget asin(S / D); betaTarget Theta - gammaTarget pi/2; % 角度周期性把beta归一到[0, 2*pi)区间 betaTarget mod(betaTarget, 2*pi); % 转换为探测器通道下标和旋转角下标 gammaIdx (gammaTarget * 180/pi - fanSensorPos(1)) / fanSpacing 1; betaIdx (betaTarget - betaVec(1)) / (betaVec(2) - betaVec(1)) 1; % 对每个目标投影角做二维插值 for k 1:numel(theta) projParallel(:, k) interp2(F, gammaIdx(:, k), betaIdx(:, k), linear, 0); end这段代码有几个点必须解释清楚。第一为什么beta的公式是θ - gamma π/2。我在前面的原理部分说过θ的基准在不同软件里可能有差异我自己这个组合是经过点源校准后的结果。如果你直接套用后发现重建图像旋转了90度不要急着改整个算法先把θ减去一个常数试试。第二mod函数处理角度周期性。扇束扫描覆盖了360度而目标平行束只需要0到179度所以反查出来的β可能落在任何区域mod一下就能保证落在采集范围内。如果只有180度加扇角的短扫描数据简单mod会引入错误这一点我在后面的常见问题部分细说。第三interp2里的extrapval参数填了0。实际重排时目标网格边缘可能会有极少数点落在探测器的物理覆盖范围之外填0是最简单的处理。但这个“0”不是真实测量值所以重建图像外圈可能出现轻微模糊后面可以通过裁剪径向范围来改善。3.3 重建验证用iradon检验重排是否正确重排做没做对空口无凭必须重建出来看效果。用iradon一把梭% 用平行束FBP重建 recon iradon(projParallel, theta_deg, Linear, Ram-Lak); % 显示对比 figure; subplot(1,3,1); imshow(img, []); title(原始体模); subplot(1,3,2); imshow(recon, []); title(重排后平行束FBP重建); subplot(1,3,3); imshow(img - recon, []); title(误差图); % 计算均方根误差 rmse sqrt(mean((img(:) - recon(:)).^2)); fprintf(RMSE %.4f\n, rmse);iradon的调用不复杂但有一个细节容易忽视projParallel的第一行是径向坐标最大的那条投影线这个顺序要和iradon内部的径向坐标定义对齐。通常我们用linspace(-smax, smax, M)生成的s是从负到正iradon默认也认为第一行对应图像下侧的径向位置所以在显示重建图后如果发现上下翻转把s倒序即可不影响数值准确性。为了进一步验证还可以生成一个点源体模做同样的流程把img设成中心区域一个单点做重排再重建观察重建点是否位于体模中心附近。这一步能在几分钟内暴露角度基准错误和镜像问题是我强烈建议做的自检步骤。3.4 参数选择的几个“不要”清单参数调参看起来小事实际上对结果影响很大。列一个我平时写代码时会在注释里标注的清单不要为了省内存把s的数量设得比探测器通道数少太多。如果s数量比M小一个数量级重排后的径向分辨率会不足重建图像会出现明显模糊。我习惯把s数量直接等于探测器通道数够用又不会太慢。不要把smax设得比D·sin(γ_max)大。超过这个范围的s在反查时算出的γ会超出探测器覆盖角度interp2只能靠extrapval强行补0重建图像外圈会多出一圈假组织。不要为了追求“重建更清晰”无限减小theta步长。平行束投影角度步长与扇形旋转角步长不一致时重排插值并不会有额外信息增益反而增加计算量。一般让theta步长和扇束旋转角步长相同即可。不要忘记检查D是否足够大。fanbeam要求D大于图像对角线的一半否则射线不能完整覆盖体模重排后的边缘全是NaN或0重建结果会出现中心区域正常、边缘大面积伪影。4. 常见问题与排查经验4.1 重建图像旋转或镜像角度基准混乱是头号元凶我做这个实验遇到的第一大坑就是重建出来的图像旋转了90度还带着左右镜像。原因就是平行束投影角θ和扇束旋转角β、探测器夹角γ之间的基准没有对齐。扇束数据里β0、γ0的那条射线是沿某个特定方向穿过物体而iradon在θ0时默认的投影方向可能是另一个方向。这两者的初始朝向一旦差了一个固定角度重建图像就会整体旋转如果γ的符号取反了图像就会镜像。排查方法很简单做一个非对称的测试体模比如在中心右上方放一个小亮点重排并重建看亮点出现在什么位置。如果出现在左上方说明图像镜像了把γ转换时的符号反转如果出现在左下方说明旋转了180度或90度给θ统一加个常数再试。这个过程不需要严谨推导两三次实验就能把基准找回来。4.2 插值结果出现NaN或边缘异常径向范围没控制好interp2在查询点超出原始数据范围时会返回NaN而NaN一旦进入iradon的输入可能就是一团糟的输出。我在早期版本里遇到过重建图像中央出现大块空洞后来才发现是smax算大了目标网格里有相当一部分点对应的γ角度超出了探测器覆盖范围。解决办法smax不要直接取经验值应该从fanSensorPos里读实际的最大通道角度再用smax D*sin(max(abs(fanSensorPos)*pi/180))计算。这样能保证径向范围被探测器物理覆盖严格约束。另外如果重排后projParallel里仍然有NaN可以在插值前用gammaIdx 1 | gammaIdx M这样的逻辑掩膜把越界点预先标出来这些位置的投影值可以不填或者用邻近有效值外推千万别留给重建函数去猜。4.3 条状伪影和模糊采样密度与插值精度的博弈如果你重建出来的图像整体还行但边缘和细节部分能看到细密的条状纹理通常是重排过程中径向采样或角度采样不够密。尤其是s的采样数少于探测器通道数时重排相当于对原始数据做了有损压缩信息丢失直接表现为重建图像的伪影。我的建议是先排查s数量再排查theta步长。如果是s数量不足把s数量设成与M相等甚至略多伪影会立刻减轻。如果你用的是三次插值还是觉得不够平滑那多半不是插值阶数的问题而是原始扫描角度步长太大两条相邻投影之间存在真实的信息缺口——这时候只能加大扇束旋转角采样密度重排算法是无能为力的。4.4 参数速查与调试顺序建议下面是我整理的一份参数速查表方便以后做类似实验时直接对照参数含义推荐值/取值范围常见错误与注意点D源到旋转中心距离500取决于体模尺寸和放大倍数必须大于图像半径否则覆盖不全fanSpacing等角探测器通道间距0.2~0.5度太小会大幅增加数据量太大则分辨率下降fanRotationInc扇束旋转角步长1度与theta步长尽量保持一致s数量重排后径向采样数等于探测器通道数少于探测器通道数会引入额外模糊theta_deg平行束投影角范围0:1:179常见FBP只需180度范围smax径向最大坐标D*sin(max(fanSensorPos))不要拍脑袋填越界会引入NaN和边缘伪影插值方法interp2插值选项linear最近邻会有马赛克伪影三次插值在精度不足时提升有限调试顺序建议从简到繁先用128×128的低分辨率体模、粗角度步长把全流程跑通确认重排后重建图像结构正确再加点源校准确认角度基准最后才换256×256乃至更高分辨率的体模做定量误差分析。这样能快速定位问题不会一上来就被各种伪影淹没。最后分享一个实际经验等角重排这个算法写起来就是几十行代码但它的价值远不止于“让扇束数据能喂给iradon”。我在实现过程中最大的收获是弄清了CT几何中“坐标基准”这件事有多重要——角度方向定义差一点重建结果就差得十万八千里。后来做锥束重建、做投影数据配准遇到几何问题时的排查思路基本都是从这里迁移过去的。再分享一个调试小技巧重排前先把扇束投影F用imagesc画出来再把重排后的projParallel用imagesc画出来对比两张图的大致形状。正常情况两者都该呈现类似“正弦图”的条纹结构只是坐标轴含义不同如果看起来完全对不上多半是坐标轴顺序搞反了这时先调轴再谈插值。我在实际项目里靠这个小检查省下了无数时间看投影图往往比直接看重建图更能快速定位几何错误。