简介CT图像重建是医学影像与生物医学工程等专业的重要知识点。这份PPT课件围绕平行束与扇形束算法转换展开从X射线投影数据采集与中心切片定理出发系统讲解滤波反投影FBP的数学原理并详细演示直角坐标系到极坐标系的雅可比变换过程。针对扇形束无对应中心切片定理的难点课件给出将扇形束射线按平行角度分组、转化为平行束问题进行重建的推导思路并完整展示等角度扇形重建算法的坐标替换、滤波与反投影步骤。资源为1个pptx文件共26页内容涵盖公式推导、几何关系化简和短扫描冗余分析适合需要深入理解CT重建数学基础或准备相关课程汇报的读者。压缩包仅697KB已有122人学习是一份精炼且逻辑清晰的算法讲解型课件。1. 平行束和扇形束CT算法转换到底在转什么CT重建的教科书推导几乎都以平行束为起点射线严格平行每个角度下的投影是一组等间隔线积分。可临床和工业CT从第一代之后就不再这样采集旋转球管加弧形或平面探测器天然产生的是扇形束锥形束更是扇形束的直接堆叠。于是平行束算法成熟而现场数据是扇形束的成了做重建的工程师都要过的一道关。平行束和扇形束算法的转换核心是把投影几何映射关系在重建算法里换成等价表达。两条常用路径重排法把扇形束投影重采样成虚拟平行束再走现成的平行束滤波反投影FBP或者直接推导扇形束自身的FBP公式把预加权和距离加权嵌进反投影。两条路径的前提相同搞清探测器坐标、旋转角、扇形角三个变量在两种几何下怎么互相表达。这篇文章写给重建模块开发、扇束截断分析或刚进入CT方向的工程师。2. 平行束FBP的数学骨架与扇形束几何参数平行束投影的定义为$$p_\theta(s) \iint f(x,y),\delta(x\cos\theta y\sin\theta - s),dx,dy$$其中 $\theta$ 是旋转角$s$ 是探测器平移坐标。重建目标是从所有角度下观测到的 $p_\theta(s)$ 恢复 $f(x,y)$。中心切片定理是这一切的纽带对固定 $\theta$ 的投影做一维傅里叶变换得到的是 $f(x,y)$ 的二维频谱中过原点且方向为 $\theta$ 的那条切片。这个定理直接决定了 FBP 实现的步骤。二维频谱在极坐标下用 $(\theta, \omega)$ 采样频域中心点被所有角度重复累计而高频在径向只有稀疏覆盖直接逆变换会出现低频过重的伪影必须用斜坡滤波器 $|\omega|$ 补偿。于是有了经典的滤波反投影公式$$f(x,y) \int_0^\pi \left[p_\theta(s) * h(s)\right]_{sx\cos\theta y\sin\theta} d\theta$$这里的 $h(s)$ 是斜坡滤波器的空域形式。工程实现里很少做空域卷积常规做法是把每行投影做 FFT、乘频域滤波器、再 IFFT 回空域之后在反投影阶段按像素坐标查表累加。2.1 平行束投影与中心切片定理的推导起点中心切片定理的价值是把“采集到的投影”和“被扫描物体本身”在频域上连接起来。重建一段已经完成的角度覆盖后理论上任何缺失角度都对应二维频谱中的楔形空洞这也是扇束截断问题在频域上的根源。理解这个结构对后续做转换非常关键无论投影来自平行束还是扇形束只要反投影时每个角度的几何关系一致频谱覆盖的问题就只跟角度采样范围相关跟探测器形态无关。这解释了为什么重排法能把扇形束数据硬生生拉回平行束框架也能解释为什么只覆盖 $180^\circ$ 扇角的短扫描重建图像会有方向性伪影。2.1.1 离散斜坡滤波器的代码实现看看频繁被引用的ram_lak_filter和平行束反投影主循环的写法。这里故意用嵌套循环不追求效率只为把几何关系展示清楚import numpy as np def ram_lak_filter(ndet): # 频域幅值 |omega|长度必须与投影行数一致 omega np.fft.fftfreq(ndet, d1.0) return np.abs(omega) def fbp_parallel(sinogram, angles, d_s1.0): # sinogram: (n_angles, ndet)每行对应一个平行束投影角度 n_angles, ndet sinogram.shape filtered np.zeros_like(sinogram) filt ram_lak_filter(ndet) # 第一步逐角度滤波 for i in range(n_angles): proj_fft np.fft.fft(sinogram[i]) filtered[i] np.fft.ifft(proj_fft * filt).real # 第二步反投影累加其中 t 是像素在旋转角上的投影坐标 recon np.zeros((ndet, ndet)) for i in range(n_angles): theta angles[i] for ix in range(ndet): x (ix - ndet / 2.0) * d_s for iy in range(ndet): y (iy - ndet / 2.0) * d_s t x * np.cos(theta) y * np.sin(theta) idx int(np.floor(t / d_s ndet / 2.0)) if 0 idx ndet: recon[ix, iy] filtered[i, idx] return recon * np.pi / n_angles这段代码有三个细节决定重建质量滤波器长度必须等于探测器通道数ndet不能随意截短滤波和角度解耦每个角度独立处理反投影里t x cosθ y sinθ是平行束下唯一的几何运算。等切到扇形束之后这行要被替换成源坐标和扇形角组合的表达式后面会专门展开。2.2 等角与等间距扇形束几何的对照进入扇形束之前先区分两种探测器模型。工业CT里旋转台加直线阵列探测器射线在平板上采样通道间距固定这是等间距模型医用CT和部分微焦斑系统用弧形探测器相邻通道夹角相同这是等角模型。两种模型的源点和旋转轴关系一致差别只在探测器表面参数化不同。项目平行束扇形束等角扇形束等间距角度变量$\theta$$0 \sim \pi$$\beta$$0 \sim 2\pi$$\beta$$0 \sim 2\pi$射线横向位置$s$$D\sin\gamma$$s D t / \sqrt{L^2 t^2}$射线方向互相平行汇聚于点源汇聚于点源滤波轴$s$ 轴$\gamma$ 轴$t$ 轴反投影附加权重无$1/L^2$$1/L^2$这个表至少要背下前两行后面所有公式推导都以它为起点。等间距模型里 $\arctan(t/L)$ 的非线性关系决定了它的重排和显式FBP都要比等角模型多一次三角变换。2.3 从平行束到扇形束的积分测度变化在平行束 FBP 公式中积分变量是 $(\theta, s)$。扇形束里要把 $(s, \theta)$ 替换成 $(\gamma, \beta)$。坐标变换 $(s, \theta) (D\sin\gamma,\ \beta\gamma)$ 对应的雅可比行列式为$$ \left|\det\begin{pmatrix} \partial s/\partial \gamma \partial s/\partial \beta \ \partial \theta/\partial \gamma \partial \theta/\partial \beta \end{pmatrix}\right| D\cos\gamma $$$D\cos\gamma$ 是扇形束重建里第一处显式出现的几何因子靠近中心射线的投影保留完整权重偏离中心越远权重越小。等间距几何里这个因子会变成 $D^3 / (L^2t^2)^{3/2}$ 一类形式。第二个显式变化是反投影到具体像素时像素到源的距离 $L$ 出现在分母上即 $1/L^2$ 权重。这两个权重加上角度积分区间从 $\pi$ 扩展到 $2\pi$ 带来的 $1/2$ 因子构成扇形束 FBP 相对平行束 FBP 的全部差别。理解了这一点后面看代码就不会被公式表面绕晕。3. 重排法转换把扇形束投影重采样成平行束重排法在工程里用得最普遍因为绝大多数重建框架都内置了平行束 FBP只需要把扇形束数据整理成平行束格式即可。代价是多一次插值。3.1 重排的核心映射关系固定源在角度 $\beta$ 处发射扇形角为 $\gamma$ 的射线这条射线在平行束几何中对应的旋转角和平移坐标是$$\theta \beta \gamma, \qquad s D\sin\gamma$$把扇形束投影 $p(\beta, \gamma)$ 写入平行束数组 $p(\theta, s)$本质上就是把散点按上面两式映射到 $(\theta, s)$ 网格。这里有一个必须处理的周期问题平行束重建只需要 $\theta$ 从 $0$ 到 $\pi$而源角度 $\beta$ 覆盖 $[0, 2\pi)$。映射后同一个 $\theta$ 可能对应多个 $\beta$ 来源取值时可以直接让插值算法决定也可以把同一 $\theta$ 区间内的多条射线值做加权平均。前者实现更简便后者数据利用更充分噪声更低代价是角度方向要多做一次合并。3.2 等间距探测器时的映射差异等间距扇形束的通道坐标沿平板均匀排列通道位置 $t$ 和扇形角之间满足$$\gamma \arctan\frac{t}{L}$$其中 $L$ 是源到探测器平面的垂直距离。对应的平行束坐标是$$s \frac{D,t}{\sqrt{L^2 t^2}}$$这条公式是非线性的通道方向不能做线性搬移必须查表插值。很多现场数据文件把探测器通道写成等间距编号却按等角公式做重排重建结果四周出现拉伸变形靠近视场边缘的细节弯成弧线。正确做法是先根据探测器物理尺寸把每个通道的 $\gamma$ 求出来之后所有步骤都按非均匀 $\gamma$ 轴处理。3.3 用MATLAB实现整段重排从目标网格倒查插值这里给出一个等角到平行的重排示例用“从目标网格反查输入网格”的写法function p_par rebin_fan_to_par(p_fan, beta_axis, gamma_axis, D, n_theta, s_range) % p_fan: 扇形束投影尺寸 [n_gamma, n_beta] % beta_axis: 源旋转角弧度[0, 2*pi) % gamma_axis: 扇形角弧度[-gamma_max, gamma_max] % D: 源到旋转中心距离 % n_theta: 重排后平行束角度个数 % s_range: 重排后平移坐标 s 的半宽 n_beta length(beta_axis); n_gamma length(gamma_axis); d_gamma gamma_axis(2) - gamma_axis(1); theta_axis linspace(0, pi, n_theta); s_axis linspace(-s_range, s_range, n_gamma); [THETA, S] meshgrid(theta_axis, s_axis); % 从目标网格反推扇形束坐标 gamma_target asin(S / D); % 由 s 反解 gamma beta_target THETA - gamma_target; % 由 theta 与 gamma 反解 beta beta_target mod(beta_target, 2*pi); % 归一到 [0, 2*pi) % 在扇形束原始网格上做双线性插值越界填 0 p_par interp2(beta_axis, gamma_axis, p_fan, ... beta_target, gamma_target, linear, 0); p_par(isnan(p_par)) 0; end逻辑说明我没有从 $(\beta, \gamma)$ 正向散点往外填而是从目标平行束网格反查扇形束坐标这是插值里的标准做法保证每个输出采样点都有确定值不会留下空洞。interp2的参数顺序是(beta_axis, gamma_axis, p_fan, ...)所以查询点矩阵beta_target和gamma_target的尺寸应当与目标网格一致。参数说明n_theta不要超过n_beta的两倍否则角度方向插值过密重建角度域出现条状伪影s_range不应超过 $D\sin\gamma_{\max}$更大只会在边缘引入空数据行。插值用linear足够cubic在投影端点会产生振铃反而污染边缘视角。提示重排法虽然实现省事但插值误差会直接进入重建结果。需要定量分析密度或做边缘精确测量的场景最好改用扇形束 FBP 直接重建。3.4 重排结果检查与常见伪影排查重排结果最常见的三个问题是角度混叠、中心偏移和边缘数据不足。表现原因处理图像边缘半月亮形伪影旋转中心偏移未标定对每个角度求投影质心取均值做偏移补偿角度方向细密条纹源角度步长过大重排后角度覆盖不足增加源角度采样或把n_theta减半视场外区域发黑发虚s_range超出 $D\sin\gamma_{\max}$缩小重建视场或增加扇形角覆盖旋转中心偏移是最常见的坑。拿到新数据先跑一个均匀圆模体或金属球标定把每个角度投影的物质中心序列算出来均值偏移量直接补偿到s_axis上再跑重排才不会出现那种看起来很对称但细节全糊的半月形伪影。4. 扇形束FBP直接重建从平行束推导的加权公式不经过重排直接在扇形束原始坐标上做滤波反投影中间省掉插值代价是公式里多出两个几何权重。4.1 从平行束到扇形束的三个改动把公式 $(s, \theta) (D\sin\gamma,\ \beta\gamma)$ 代入平行束 FBP整体会出现三处变化。第一积分变量换成 $(\beta, \gamma)$ 后产生雅可比因子 $D\cos\gamma$第二反投影到具体像素减不再请求查 $s$ 等于某个常量而是要把像素到当前位置源的距离送进 $1/L^2$ 权重第三角度积分范围从 $\pi$ 变为 $2\pi$归一化系数多一个 $1/2$。重建点 $(x,y)$ 在源角度为 $\beta$ 时对应的扇形角 $\gamma$ 计算如下$$\gamma \operatorname{atan2}\left(y\cos\beta - x\sin\beta,\ D - x\cos\beta - y\sin\beta\right)$$像素到源的距离 $L$ 用两点坐标直接算。这两个量在反投影主循环里每个像素、每个角度都要重新算计算量比重排后的平行束反投影大不少这也是当年硬件不强时重排法更受欢迎的原因之一。4.2 Python实现扇形束FBP的最小骨架def fan_fbp_direct(sinogram, beta_axis, gamma_axis, D, pixel_size1.0): # sinogram: (n_beta, n_gamma)行对应源角度列对应扇形角 n_beta, n_gamma sinogram.shape R D * np.sin(gamma_axis[-1]) # 理论重建视场半径 n_pix int(2 * R / pixel_size) recon np.zeros((n_pix, n_pix)) d_beta beta_axis[1] - beta_axis[0] d_gamma gamma_axis[1] - gamma_axis[0] # 要求等间隔 # 1) 预加权D * cos(gamma) preweight D * np.cos(gamma_axis) data sinogram * preweight[np.newaxis, :] # 2) 沿扇形角方向做斜坡滤波滤波方向换成 gamma 轴 filt np.abs(np.fft.fftfreq(n_gamma, dd_gamma)) for i in range(n_beta): data[i] np.fft.ifft(np.fft.fft(data[i]) * filt).real # 3) 反投影累加 1/L^2 权重 for i in range(n_beta): beta beta_axis[i] sx, sy D * np.cos(beta), D * np.sin(beta) # 当前源位置 for ix in range(n_pix): x (ix - n_pix / 2.0) * pixel_size for iy in range(n_pix): y (iy - n_pix / 2.0) * pixel_size num y * np.cos(beta) - x * np.sin(beta) den D - x * np.cos(beta) - y * np.sin(beta) gamma_p np.arctan2(num, den) j int(np.round((gamma_p - gamma_axis[0]) / d_gamma)) if 0 j n_gamma: L2 (x - sx) ** 2 (y - sy) ** 2 recon[ix, iy] data[i, j] / L2 # 归一化角度步长 * 扇形角步长 * 半周扩展因子 recon * d_beta * d_gamma / 2.0 return recon代码逻辑分成三步分别对应前面说的三个改动。第一步预加权是在滤波之前做的与平行束“先滤波再反投影”的顺序一致第二步滤波直接沿 $\gamma$ 轴做用的还是一维斜坡滤波器但频率轴的单位从“每像素”换成“每弧度”第三步反投影时每个像素查的是 $\gamma$ 对应的通道而不是简单的投影坐标。参数上最容易被忽略的是d_gamma。这段实现要求 $\gamma$ 轴等间隔如果探测器通道是非均匀角度采样得先做重采样或者改用分段积分。另一个容易出错的点是pixel_size与D、s_range的单位一致性工业CT数据里经常混用毫米和微米差一个数量级会让重建图像尺寸完全不对。4.3 直接扇形束FBP与重排法的参数对比对比项重排法 平行束FBP直接扇形束FBP中间插值有两维无预加权不需要$D\cos\gamma$滤波方向$s$ 轴$\gamma$ 轴反投影权重无$1/L^2$归一化系数沿用平行束$d\beta \cdot d\gamma / 2$噪声表现插值有平滑效果距离权重放大靠近源的噪声计算量插值额外开销反投影每点多次三角运算实际选型时扇束角度跨度小、探测器通道多、对分辨率有要求的场景直接扇形束 FBP 更合适本身就要先做视角重排、后续还要处理运动伪影的数据流重排法更好接入现成流程。5. 用数值模体验证平行束‑扇形束转换的三个检查点动手写转换代码不难难的是确认转换没出错。手里没有扫描仪时用 Shepp-Logan 模体做数字仿真就够。Shepp-Logan 由十几个椭圆拼成有解析投影公式能以浮点精度同时生成平行束和扇形束投影也适合做工业CT重建算法回归测试的基准。5.1 检查点一均匀圆盘的径向剖面先用单一均匀圆盘做模体分别生成平行束和扇形束投影再做两种重建沿直径取剖面。重排法如果存在角度混叠剖面会呈锯齿状扇形束 FBP 如果权重写错剖面会有明显的中心突起或凹陷。把剖面归一化后对比偏差超过 2% 就回去查 $D\cos\gamma$ 权重的符号或者查 $1/L^2$ 是不是被展开成了 $1/L$。5.2 检查点二两种算法在重叠区域的NRMSD从同一个扇形束数据集出发分别跑重排法和直接扇形束 FBP输出两张重建图计算旋转中心附近一个小区域内的归一化均方根偏差将两张图都裁剪到旋转中心半径 $R$ 内再对比。重排插值误差与 FBP 几何误差混在一起时NRMSD 一般落在 $10^{-3}$ 量级附近一旦超过 $10^{-2}$基本可以断定某个权重方向反了。import numpy as np # 假设 recon_a 是重排法结果recon_b 是直接扇形束FBP结果 cx cy recon_a.shape[0] // 2 yy, xx np.ogrid[:recon_a.shape[0], :recon_a.shape[1]] mask (xx - cx) ** 2 (yy - cy) ** 2 30 ** 2 nrmsd np.sqrt(np.sum((recon_a[mask] - recon_b[mask]) ** 2)) / \ np.sqrt(np.sum(recon_a[mask] ** 2)) print(fNRMSD {nrmsd:.6f})这个阈值判断写进自动化脚本后任何几何参数改动都能被快速识别。5.3 检查点三视场截断边界改变程序里扇形角上限 $\gamma_{\max}$记录重建图中不再出现断层伪影的最大半径。该半径应该跟随 $D\sin\gamma_{\max}$ 缓慢变化。如果它明显小于理论值说明角度方向采样不足需要减小重建矩阵尺寸或增加源角度数如果明显大于理论值多半是预加权没有限制视场范围截断伪影被当成真实信号重建了。把这三个检查点放进回归脚本每次改了采集几何、换了探测器间距或调了插值方法后跑一遍比对着重建图肉眼判断可靠得多。本文还有配套的精品资源点击获取 CT重建中平行束与扇形束算法转换：原理、重排与FBP实现
简介CT图像重建是医学影像与生物医学工程等专业的重要知识点。这份PPT课件围绕平行束与扇形束算法转换展开从X射线投影数据采集与中心切片定理出发系统讲解滤波反投影FBP的数学原理并详细演示直角坐标系到极坐标系的雅可比变换过程。针对扇形束无对应中心切片定理的难点课件给出将扇形束射线按平行角度分组、转化为平行束问题进行重建的推导思路并完整展示等角度扇形重建算法的坐标替换、滤波与反投影步骤。资源为1个pptx文件共26页内容涵盖公式推导、几何关系化简和短扫描冗余分析适合需要深入理解CT重建数学基础或准备相关课程汇报的读者。压缩包仅697KB已有122人学习是一份精炼且逻辑清晰的算法讲解型课件。1. 平行束和扇形束CT算法转换到底在转什么CT重建的教科书推导几乎都以平行束为起点射线严格平行每个角度下的投影是一组等间隔线积分。可临床和工业CT从第一代之后就不再这样采集旋转球管加弧形或平面探测器天然产生的是扇形束锥形束更是扇形束的直接堆叠。于是平行束算法成熟而现场数据是扇形束的成了做重建的工程师都要过的一道关。平行束和扇形束算法的转换核心是把投影几何映射关系在重建算法里换成等价表达。两条常用路径重排法把扇形束投影重采样成虚拟平行束再走现成的平行束滤波反投影FBP或者直接推导扇形束自身的FBP公式把预加权和距离加权嵌进反投影。两条路径的前提相同搞清探测器坐标、旋转角、扇形角三个变量在两种几何下怎么互相表达。这篇文章写给重建模块开发、扇束截断分析或刚进入CT方向的工程师。2. 平行束FBP的数学骨架与扇形束几何参数平行束投影的定义为$$p_\theta(s) \iint f(x,y),\delta(x\cos\theta y\sin\theta - s),dx,dy$$其中 $\theta$ 是旋转角$s$ 是探测器平移坐标。重建目标是从所有角度下观测到的 $p_\theta(s)$ 恢复 $f(x,y)$。中心切片定理是这一切的纽带对固定 $\theta$ 的投影做一维傅里叶变换得到的是 $f(x,y)$ 的二维频谱中过原点且方向为 $\theta$ 的那条切片。这个定理直接决定了 FBP 实现的步骤。二维频谱在极坐标下用 $(\theta, \omega)$ 采样频域中心点被所有角度重复累计而高频在径向只有稀疏覆盖直接逆变换会出现低频过重的伪影必须用斜坡滤波器 $|\omega|$ 补偿。于是有了经典的滤波反投影公式$$f(x,y) \int_0^\pi \left[p_\theta(s) * h(s)\right]_{sx\cos\theta y\sin\theta} d\theta$$这里的 $h(s)$ 是斜坡滤波器的空域形式。工程实现里很少做空域卷积常规做法是把每行投影做 FFT、乘频域滤波器、再 IFFT 回空域之后在反投影阶段按像素坐标查表累加。2.1 平行束投影与中心切片定理的推导起点中心切片定理的价值是把“采集到的投影”和“被扫描物体本身”在频域上连接起来。重建一段已经完成的角度覆盖后理论上任何缺失角度都对应二维频谱中的楔形空洞这也是扇束截断问题在频域上的根源。理解这个结构对后续做转换非常关键无论投影来自平行束还是扇形束只要反投影时每个角度的几何关系一致频谱覆盖的问题就只跟角度采样范围相关跟探测器形态无关。这解释了为什么重排法能把扇形束数据硬生生拉回平行束框架也能解释为什么只覆盖 $180^\circ$ 扇角的短扫描重建图像会有方向性伪影。2.1.1 离散斜坡滤波器的代码实现看看频繁被引用的ram_lak_filter和平行束反投影主循环的写法。这里故意用嵌套循环不追求效率只为把几何关系展示清楚import numpy as np def ram_lak_filter(ndet): # 频域幅值 |omega|长度必须与投影行数一致 omega np.fft.fftfreq(ndet, d1.0) return np.abs(omega) def fbp_parallel(sinogram, angles, d_s1.0): # sinogram: (n_angles, ndet)每行对应一个平行束投影角度 n_angles, ndet sinogram.shape filtered np.zeros_like(sinogram) filt ram_lak_filter(ndet) # 第一步逐角度滤波 for i in range(n_angles): proj_fft np.fft.fft(sinogram[i]) filtered[i] np.fft.ifft(proj_fft * filt).real # 第二步反投影累加其中 t 是像素在旋转角上的投影坐标 recon np.zeros((ndet, ndet)) for i in range(n_angles): theta angles[i] for ix in range(ndet): x (ix - ndet / 2.0) * d_s for iy in range(ndet): y (iy - ndet / 2.0) * d_s t x * np.cos(theta) y * np.sin(theta) idx int(np.floor(t / d_s ndet / 2.0)) if 0 idx ndet: recon[ix, iy] filtered[i, idx] return recon * np.pi / n_angles这段代码有三个细节决定重建质量滤波器长度必须等于探测器通道数ndet不能随意截短滤波和角度解耦每个角度独立处理反投影里t x cosθ y sinθ是平行束下唯一的几何运算。等切到扇形束之后这行要被替换成源坐标和扇形角组合的表达式后面会专门展开。2.2 等角与等间距扇形束几何的对照进入扇形束之前先区分两种探测器模型。工业CT里旋转台加直线阵列探测器射线在平板上采样通道间距固定这是等间距模型医用CT和部分微焦斑系统用弧形探测器相邻通道夹角相同这是等角模型。两种模型的源点和旋转轴关系一致差别只在探测器表面参数化不同。项目平行束扇形束等角扇形束等间距角度变量$\theta$$0 \sim \pi$$\beta$$0 \sim 2\pi$$\beta$$0 \sim 2\pi$射线横向位置$s$$D\sin\gamma$$s D t / \sqrt{L^2 t^2}$射线方向互相平行汇聚于点源汇聚于点源滤波轴$s$ 轴$\gamma$ 轴$t$ 轴反投影附加权重无$1/L^2$$1/L^2$这个表至少要背下前两行后面所有公式推导都以它为起点。等间距模型里 $\arctan(t/L)$ 的非线性关系决定了它的重排和显式FBP都要比等角模型多一次三角变换。2.3 从平行束到扇形束的积分测度变化在平行束 FBP 公式中积分变量是 $(\theta, s)$。扇形束里要把 $(s, \theta)$ 替换成 $(\gamma, \beta)$。坐标变换 $(s, \theta) (D\sin\gamma,\ \beta\gamma)$ 对应的雅可比行列式为$$ \left|\det\begin{pmatrix} \partial s/\partial \gamma \partial s/\partial \beta \ \partial \theta/\partial \gamma \partial \theta/\partial \beta \end{pmatrix}\right| D\cos\gamma $$$D\cos\gamma$ 是扇形束重建里第一处显式出现的几何因子靠近中心射线的投影保留完整权重偏离中心越远权重越小。等间距几何里这个因子会变成 $D^3 / (L^2t^2)^{3/2}$ 一类形式。第二个显式变化是反投影到具体像素时像素到源的距离 $L$ 出现在分母上即 $1/L^2$ 权重。这两个权重加上角度积分区间从 $\pi$ 扩展到 $2\pi$ 带来的 $1/2$ 因子构成扇形束 FBP 相对平行束 FBP 的全部差别。理解了这一点后面看代码就不会被公式表面绕晕。3. 重排法转换把扇形束投影重采样成平行束重排法在工程里用得最普遍因为绝大多数重建框架都内置了平行束 FBP只需要把扇形束数据整理成平行束格式即可。代价是多一次插值。3.1 重排的核心映射关系固定源在角度 $\beta$ 处发射扇形角为 $\gamma$ 的射线这条射线在平行束几何中对应的旋转角和平移坐标是$$\theta \beta \gamma, \qquad s D\sin\gamma$$把扇形束投影 $p(\beta, \gamma)$ 写入平行束数组 $p(\theta, s)$本质上就是把散点按上面两式映射到 $(\theta, s)$ 网格。这里有一个必须处理的周期问题平行束重建只需要 $\theta$ 从 $0$ 到 $\pi$而源角度 $\beta$ 覆盖 $[0, 2\pi)$。映射后同一个 $\theta$ 可能对应多个 $\beta$ 来源取值时可以直接让插值算法决定也可以把同一 $\theta$ 区间内的多条射线值做加权平均。前者实现更简便后者数据利用更充分噪声更低代价是角度方向要多做一次合并。3.2 等间距探测器时的映射差异等间距扇形束的通道坐标沿平板均匀排列通道位置 $t$ 和扇形角之间满足$$\gamma \arctan\frac{t}{L}$$其中 $L$ 是源到探测器平面的垂直距离。对应的平行束坐标是$$s \frac{D,t}{\sqrt{L^2 t^2}}$$这条公式是非线性的通道方向不能做线性搬移必须查表插值。很多现场数据文件把探测器通道写成等间距编号却按等角公式做重排重建结果四周出现拉伸变形靠近视场边缘的细节弯成弧线。正确做法是先根据探测器物理尺寸把每个通道的 $\gamma$ 求出来之后所有步骤都按非均匀 $\gamma$ 轴处理。3.3 用MATLAB实现整段重排从目标网格倒查插值这里给出一个等角到平行的重排示例用“从目标网格反查输入网格”的写法function p_par rebin_fan_to_par(p_fan, beta_axis, gamma_axis, D, n_theta, s_range) % p_fan: 扇形束投影尺寸 [n_gamma, n_beta] % beta_axis: 源旋转角弧度[0, 2*pi) % gamma_axis: 扇形角弧度[-gamma_max, gamma_max] % D: 源到旋转中心距离 % n_theta: 重排后平行束角度个数 % s_range: 重排后平移坐标 s 的半宽 n_beta length(beta_axis); n_gamma length(gamma_axis); d_gamma gamma_axis(2) - gamma_axis(1); theta_axis linspace(0, pi, n_theta); s_axis linspace(-s_range, s_range, n_gamma); [THETA, S] meshgrid(theta_axis, s_axis); % 从目标网格反推扇形束坐标 gamma_target asin(S / D); % 由 s 反解 gamma beta_target THETA - gamma_target; % 由 theta 与 gamma 反解 beta beta_target mod(beta_target, 2*pi); % 归一到 [0, 2*pi) % 在扇形束原始网格上做双线性插值越界填 0 p_par interp2(beta_axis, gamma_axis, p_fan, ... beta_target, gamma_target, linear, 0); p_par(isnan(p_par)) 0; end逻辑说明我没有从 $(\beta, \gamma)$ 正向散点往外填而是从目标平行束网格反查扇形束坐标这是插值里的标准做法保证每个输出采样点都有确定值不会留下空洞。interp2的参数顺序是(beta_axis, gamma_axis, p_fan, ...)所以查询点矩阵beta_target和gamma_target的尺寸应当与目标网格一致。参数说明n_theta不要超过n_beta的两倍否则角度方向插值过密重建角度域出现条状伪影s_range不应超过 $D\sin\gamma_{\max}$更大只会在边缘引入空数据行。插值用linear足够cubic在投影端点会产生振铃反而污染边缘视角。提示重排法虽然实现省事但插值误差会直接进入重建结果。需要定量分析密度或做边缘精确测量的场景最好改用扇形束 FBP 直接重建。3.4 重排结果检查与常见伪影排查重排结果最常见的三个问题是角度混叠、中心偏移和边缘数据不足。表现原因处理图像边缘半月亮形伪影旋转中心偏移未标定对每个角度求投影质心取均值做偏移补偿角度方向细密条纹源角度步长过大重排后角度覆盖不足增加源角度采样或把n_theta减半视场外区域发黑发虚s_range超出 $D\sin\gamma_{\max}$缩小重建视场或增加扇形角覆盖旋转中心偏移是最常见的坑。拿到新数据先跑一个均匀圆模体或金属球标定把每个角度投影的物质中心序列算出来均值偏移量直接补偿到s_axis上再跑重排才不会出现那种看起来很对称但细节全糊的半月形伪影。4. 扇形束FBP直接重建从平行束推导的加权公式不经过重排直接在扇形束原始坐标上做滤波反投影中间省掉插值代价是公式里多出两个几何权重。4.1 从平行束到扇形束的三个改动把公式 $(s, \theta) (D\sin\gamma,\ \beta\gamma)$ 代入平行束 FBP整体会出现三处变化。第一积分变量换成 $(\beta, \gamma)$ 后产生雅可比因子 $D\cos\gamma$第二反投影到具体像素减不再请求查 $s$ 等于某个常量而是要把像素到当前位置源的距离送进 $1/L^2$ 权重第三角度积分范围从 $\pi$ 变为 $2\pi$归一化系数多一个 $1/2$。重建点 $(x,y)$ 在源角度为 $\beta$ 时对应的扇形角 $\gamma$ 计算如下$$\gamma \operatorname{atan2}\left(y\cos\beta - x\sin\beta,\ D - x\cos\beta - y\sin\beta\right)$$像素到源的距离 $L$ 用两点坐标直接算。这两个量在反投影主循环里每个像素、每个角度都要重新算计算量比重排后的平行束反投影大不少这也是当年硬件不强时重排法更受欢迎的原因之一。4.2 Python实现扇形束FBP的最小骨架def fan_fbp_direct(sinogram, beta_axis, gamma_axis, D, pixel_size1.0): # sinogram: (n_beta, n_gamma)行对应源角度列对应扇形角 n_beta, n_gamma sinogram.shape R D * np.sin(gamma_axis[-1]) # 理论重建视场半径 n_pix int(2 * R / pixel_size) recon np.zeros((n_pix, n_pix)) d_beta beta_axis[1] - beta_axis[0] d_gamma gamma_axis[1] - gamma_axis[0] # 要求等间隔 # 1) 预加权D * cos(gamma) preweight D * np.cos(gamma_axis) data sinogram * preweight[np.newaxis, :] # 2) 沿扇形角方向做斜坡滤波滤波方向换成 gamma 轴 filt np.abs(np.fft.fftfreq(n_gamma, dd_gamma)) for i in range(n_beta): data[i] np.fft.ifft(np.fft.fft(data[i]) * filt).real # 3) 反投影累加 1/L^2 权重 for i in range(n_beta): beta beta_axis[i] sx, sy D * np.cos(beta), D * np.sin(beta) # 当前源位置 for ix in range(n_pix): x (ix - n_pix / 2.0) * pixel_size for iy in range(n_pix): y (iy - n_pix / 2.0) * pixel_size num y * np.cos(beta) - x * np.sin(beta) den D - x * np.cos(beta) - y * np.sin(beta) gamma_p np.arctan2(num, den) j int(np.round((gamma_p - gamma_axis[0]) / d_gamma)) if 0 j n_gamma: L2 (x - sx) ** 2 (y - sy) ** 2 recon[ix, iy] data[i, j] / L2 # 归一化角度步长 * 扇形角步长 * 半周扩展因子 recon * d_beta * d_gamma / 2.0 return recon代码逻辑分成三步分别对应前面说的三个改动。第一步预加权是在滤波之前做的与平行束“先滤波再反投影”的顺序一致第二步滤波直接沿 $\gamma$ 轴做用的还是一维斜坡滤波器但频率轴的单位从“每像素”换成“每弧度”第三步反投影时每个像素查的是 $\gamma$ 对应的通道而不是简单的投影坐标。参数上最容易被忽略的是d_gamma。这段实现要求 $\gamma$ 轴等间隔如果探测器通道是非均匀角度采样得先做重采样或者改用分段积分。另一个容易出错的点是pixel_size与D、s_range的单位一致性工业CT数据里经常混用毫米和微米差一个数量级会让重建图像尺寸完全不对。4.3 直接扇形束FBP与重排法的参数对比对比项重排法 平行束FBP直接扇形束FBP中间插值有两维无预加权不需要$D\cos\gamma$滤波方向$s$ 轴$\gamma$ 轴反投影权重无$1/L^2$归一化系数沿用平行束$d\beta \cdot d\gamma / 2$噪声表现插值有平滑效果距离权重放大靠近源的噪声计算量插值额外开销反投影每点多次三角运算实际选型时扇束角度跨度小、探测器通道多、对分辨率有要求的场景直接扇形束 FBP 更合适本身就要先做视角重排、后续还要处理运动伪影的数据流重排法更好接入现成流程。5. 用数值模体验证平行束‑扇形束转换的三个检查点动手写转换代码不难难的是确认转换没出错。手里没有扫描仪时用 Shepp-Logan 模体做数字仿真就够。Shepp-Logan 由十几个椭圆拼成有解析投影公式能以浮点精度同时生成平行束和扇形束投影也适合做工业CT重建算法回归测试的基准。5.1 检查点一均匀圆盘的径向剖面先用单一均匀圆盘做模体分别生成平行束和扇形束投影再做两种重建沿直径取剖面。重排法如果存在角度混叠剖面会呈锯齿状扇形束 FBP 如果权重写错剖面会有明显的中心突起或凹陷。把剖面归一化后对比偏差超过 2% 就回去查 $D\cos\gamma$ 权重的符号或者查 $1/L^2$ 是不是被展开成了 $1/L$。5.2 检查点二两种算法在重叠区域的NRMSD从同一个扇形束数据集出发分别跑重排法和直接扇形束 FBP输出两张重建图计算旋转中心附近一个小区域内的归一化均方根偏差将两张图都裁剪到旋转中心半径 $R$ 内再对比。重排插值误差与 FBP 几何误差混在一起时NRMSD 一般落在 $10^{-3}$ 量级附近一旦超过 $10^{-2}$基本可以断定某个权重方向反了。import numpy as np # 假设 recon_a 是重排法结果recon_b 是直接扇形束FBP结果 cx cy recon_a.shape[0] // 2 yy, xx np.ogrid[:recon_a.shape[0], :recon_a.shape[1]] mask (xx - cx) ** 2 (yy - cy) ** 2 30 ** 2 nrmsd np.sqrt(np.sum((recon_a[mask] - recon_b[mask]) ** 2)) / \ np.sqrt(np.sum(recon_a[mask] ** 2)) print(fNRMSD {nrmsd:.6f})这个阈值判断写进自动化脚本后任何几何参数改动都能被快速识别。5.3 检查点三视场截断边界改变程序里扇形角上限 $\gamma_{\max}$记录重建图中不再出现断层伪影的最大半径。该半径应该跟随 $D\sin\gamma_{\max}$ 缓慢变化。如果它明显小于理论值说明角度方向采样不足需要减小重建矩阵尺寸或增加源角度数如果明显大于理论值多半是预加权没有限制视场范围截断伪影被当成真实信号重建了。把这三个检查点放进回归脚本每次改了采集几何、换了探测器间距或调了插值方法后跑一遍比对着重建图肉眼判断可靠得多。本文还有配套的精品资源点击获取