企业官网建设流程全解析

简介：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]_{s=x\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, d=1.0) return np.abs(omega) def fbp_parallel(sinogram, angles, d_s=1.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和部分微焦斑系统用弧形探测器，相邻通道夹角相同，这是等角模型。两种模型的源点和旋转轴关系一致，差别只在探测器表面参数化不同。

这个表至少要背下前两行，后面所有公式推导都以它为起点。等间距模型里 $\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^2+t^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 重排结果检查与常见伪影排查

重排结果最常见的三个问题是角度混叠、中心偏移和边缘数据不足。

旋转中心偏移是最常见的坑。拿到新数据先跑一个均匀圆模体或金属球标定，把每个角度投影的物质中心序列算出来，均值偏移量直接补偿到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_size=1.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, d=d_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 更合适；本身就要先做视角重排、后续还要处理运动伪影的数据流，重排法更好接入现成流程。

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(f"NRMSD = {nrmsd:.6f}")

这个阈值判断写进自动化脚本后，任何几何参数改动都能被快速识别。

5.3 检查点三：视场截断边界

改变程序里扇形角上限 $\gamma_{\max}$，记录重建图中不再出现断层伪影的最大半径。该半径应该跟随 $D\sin\gamma_{\max}$ 缓慢变化。如果它明显小于理论值，说明角度方向采样不足，需要减小重建矩阵尺寸或增加源角度数；如果明显大于理论值，多半是预加权没有限制视场范围，截断伪影被当成真实信号重建了。

把这三个检查点放进回归脚本，每次改了采集几何、换了探测器间距或调了插值方法后跑一遍，比对着重建图肉眼判断可靠得多。

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