企业官网建设流程全解析

一次CT重建翻车，我盯着正弦图看了半个小时，最后发现锅在“平行束反投影”这个最基础的环节上。很多做图像重建的朋友可能和我一样，上来就调库，iradon一行出图，根本不知道反投影在底层做了什么，也不知道为什么直接反投影会有那种墓碑一样的星形伪影。直到自己手写了一遍反投影循环，才真正把“滤波反投影”这套逻辑嵌进脑子里。这篇东西就围绕“平行束反投影重建”这个标签，从最基本的物理直觉讲起，把Radon变换、反投影公式、斜坡滤波器、代码实现和常见伪影一次说透。

1. 平行束反投影：到底在反什么

1.1 从影子说起：什么是平行束投影

先别急着上公式。平行束投影的物理场景特别朴素——想象一束平行的X射线，从某个固定方向穿过一个物体，到达对侧的探测器。射线在穿过物质时会被衰减，衰减多少取决于它走过的路径上物质密度的积分，这个积分值就是一条射线的投影值。CT扫描的原理，就是让这个平行束绕着物体转一圈，每转一个小角度记录一组一维投影，最后用所有角度下的投影数据反推出物体内部的二维密度分布。

“平行束”这个名字强调的是射线的几何关系：同一角度下的所有X射线是彼此平行的。最早的单束平移旋转扫描方式，用的就是这种几何模型。后来的扇形束CT、锥形束CT在数学上都要更复杂，但平行束仍然是理解重建算法的最佳切入点，也是很多同步辐射光源、工业检测CT系统的实际近似模型。

把投影数据画出来，横轴是探测器位置，纵轴是旋转角度，每个像素的亮度代表那一个角度、那一条射线的投影强度，这张图就叫正弦图。名字来源于单个点状物体在图中会形成一条正弦曲线。正弦图是整个重建过程的“输入”，反投影就是要把这张图里的信息“摊回”图像空间。

1.2 反投影：把影子放回去

反投影的思路可以用一句话概括：把每一个角度下测得的一维投影值，沿着采集时的射线路径，均匀地“涂抹”回图像空间，然后把所有角度的涂抹结果累加起来。

举个例子。假设物体内部只有一个亮点，在某个角度下，探测器上对应位置会记录到一个峰值。反投影会把这个峰值沿着入射方向在整个图像上涂出一条线。换一个角度再来一次，又会涂出一条线。两个角度交叉，两条亮线会在亮点位置相交。积累足够多的角度之后，各路“涂抹”在亮点位置反复叠加，就会形成一个明显的峰值。这个峰值就是重建出来的信号。

这个过程不需要什么复杂的数学直觉，就是一种很直观的“把测量值送回原路”的做法。但是从投影到重建，关键不在于单条线，而在于所有角度的线如何叠加、如何消除彼此之间的干扰。

1.3 为什么单纯反投影会模糊：星形伪影的起源

直接反投影出来的图像，肉眼一看就有一个问题——模糊，亮点周围拖着长长的放射状条纹。这是因为单个点物体的反投影，并不是一个无限窄的点响应，而是呈1/r形式衰减的拖尾，其中r是离点源的距离。离点源越远，拖尾越小，但永远不会突然降为零。

每个方向的拖尾彼此交叠，叠加之后就在点源周围形成星形伪影。角度越少，星形伪影越明显；就算角度足够多，背景上的模糊也不会消失。用信号处理的眼光看，反投影相当于一把很钝的“锅铲”——它把能量摊出去了，却没法把该集中的能量收拢到点上。这就是为什么后来要有滤波反投影，因为需要一把锋利的刀，先把投影数据“修整”一顿，再涂回去。

提示：理解星形伪影是理解滤波的钥匙。滤波反投影中的“滤波”二字，本质上就是为了抵消反投影过程里的1/r模糊。很多初学者只记FBP公式，却不知道它到底在滤什么，后面看滤波器参数就会一头雾水。

2. 数学基底：Radon变换与投影定理

2.1 Radon变换到底做了什么

上点数学，但尽量说人话。设f(x, y)是物体的二维密度分布，平行束在角度θ下的投影函数，数学上定义为：

p(θ, s) = ∫∫ f(x, y) δ(x cosθ + y sinθ - s) dx dy

其中δ是狄拉克函数，它只在x cosθ + y sinθ = s这条直线上取非零值。所以这个式子翻译过来就是：沿着一条直线对f做线积分，把所有平行于这条直线的线积分按探测器位置s排成一排，就是当前角度下的投影。

这个变换叫Radon变换。它看起来复杂，但思想非常简单：把一个二维分布“压缩”成一堆一维函数，每个一维函数对应一个角度。CT扫描本质上就是在测量物体密度分布的Radon变换，只不过是在有限个角度、有限个探测器单元上做离散采样。

2.2 投影定理：连接时域与频域的关键桥

Radon变换和图像之间最漂亮的关系，是投影定理，也叫中心切片定理。它说的是：某角度θ下的投影p(θ, s)的一维傅里叶变换，恰好等于原图像f(x, y)的二维傅里叶变换在过原点、方向角为θ的那条直线上的切片。

这个结论的价值怎么强调都不过分。它意味着，如果采集了足够多角度的投影，理论上就测到了图像二维频谱中所有过原点的径向线上的值，也就是覆盖了整个频域平面。反过来，重建图像等价于把这些径向频谱切片“铺平”成一个完整的二维频谱，再做一次二维逆傅里叶变换。

但是注意，直接把这些频谱切片铺到频率平面上，会在低频区域重复采样，而高频区域的采样密度相对较低。反投影的模糊，本质上就和这种频谱不均匀有关。它相当于在频域中用了1/|ω|这样一个权重，导致低频被放大、高频被压低，反映在空域里就是那个1/r拖尾。

2.3 离散化：从连续公式到可计算的代码

连续公式适合理解，真正写代码时必须离散化。离散化涉及三个维度参数：

• 旋转角度个数num_angles，对应采集时转了多少个角度；
• 探测器单元个数num_detectors，对应每个角度下一维投影的采样点数；
• 重建图像的像素尺寸N × N，通常取N ≈ num_detectors，保证图像分辨率和投影采样率匹配。

离散反投影的核心循环长这样：遍历每个角度，对图像中的每个像素，计算该像素在当前角度下的射线投影位置s = x cosθ + y sinθ，然后到该角度的一维投影数组中取对应值，累加到像素上。

这个循环看着简单，写出来有不少细节。下面的代码是一个最小化的并行束反投影实现，用纯NumPy，不做任何库函数的暗箱操作：

import numpy as np def backproject(sinogram, theta_deg, img_size=None): """ 平行束直接反投影 sinogram: shape (num_angles, num_detectors) theta_deg: 各角度的度数数组，长度 num_angles img_size: 重建图像边长（像素） """ num_angles, num_detectors = sinogram.shape if img_size is None: img_size = num_detectors # 归一化坐标：像素范围 [-1, 1] x = np.linspace(-1, 1, img_size) y = np.linspace(-1, 1, img_size) X, Y = np.meshgrid(x, y) recon = np.zeros((img_size, img_size)) detector_coords = np.arange(num_detectors) for i, theta in enumerate(np.deg2rad(theta_deg)): # 每个像素对应的射线投影位置 s s = X * np.cos(theta) + Y * np.sin(theta) # 将 s 从 [-1, 1] 映射到探测器索引 [0, num_detectors-1] s_idx = (s + 1.0) * 0.5 * (num_detectors - 1) # 线性插值取投影值（用 clip 防止越界） proj = sinogram[i] s_flat = s.ravel() idx_floor = np.floor(s_idx).astype(int) idx_ceil = np.clip(idx_floor + 1, 0, num_detectors - 1) frac = s_idx - idx_floor val_floor = proj[idx_floor] val_ceil = proj[idx_ceil] value = val_floor * (1 - frac) + val_ceil * frac recon += value.reshape(img_size, img_size) # 角度归一化：把所有角度的贡献做平均 recon *= np.pi / num_angles return recon

这个实现里有两个关键细节。第一是插值：s落到两个探测器单元之间时，必须做线性插值，直接取最近邻会让重建图像出现明显的锯齿。第二是角度归一化：角度越多，累加值越大，除以角度数之后才能得到稳定的灰度尺度。我见过不少人在这一步上栽跟头，结果重建出来的图像整个偏亮，以为是窗位没调好，其实是归一化写错了。

3. 从反投影到滤波反投影：那一步关键的“修整”

3.1 反投影的频域面貌：为什么必须滤波

直接反投影的模糊可以用一个简单的观点解释：它相当于在频域对图像的真实频谱乘以了一个1/|ω|的权重。也就是说，反投影得到的图像频谱B(ω)和原始图像频谱F(ω)的关系近似为：

B(ω) ≈ F(ω) / |ω|

要还原F(ω)，最自然的方式就是把反投影结果乘以|ω|。这正是滤波反投影的思路，只不过在实际实现中，不是先反投影再乘，而是先乘再反投影，即：在频域中对每个角度的一维投影做滤波，滤波核是斜坡函数|ω|，然后再做直接反投影。两种顺序在数学上是等价的，但后者实现起来更直观、效率也更高，因为只需要对一维信号做傅里叶变换，而不是对二维图像做。

3.2 斜坡滤波器：理想频率响应与其局限

理想情况下，FBP需要一个频域响应为|ω|的一维滤波器。但直接照搬这个理想响应，在离散实现中会带来一个严重问题：高频噪声会被剧烈放大。X射线探测器、电子线路、光子统计波动都会引入高频噪声，|ω|在高频端无限增大，没有上限，所以纯斜坡滤波器在实际数据上往往会产生一幅充满噪声甚至无法观看的图像。

解决办法是加窗。常见窗函数有以下几种：

在实际工程中，Ram-Lak是“母版”，其余窗都是在它上面乘一段衰减曲线。如果做的是标准体模仿真，Ram-Lak是首选；如果是真实CT数据，Shepp-Logan或Hamming更稳。

3.3 滤波步骤的代码实现

斜坡滤波可以在频域完成，也很方便。给定一维投影数组proj，用NumPy做一次FFT，乘以|ω|，再IFFT回来即可：

def ramp_filter(proj, use_window='shepp-logan'): n = proj.shape[0] freqs = np.fft.fftfreq(n) * n # 归一化到 [-n/2, n/2) ramp = np.abs(freqs) if use_window == 'shepp-logan': # 用sinc窗调制，抑制高频噪声 # 窗长度近似取n，离散实现中常用归一化的sinc函数 window = np.sinc(freqs / n) ramp = ramp * window elif use_window == 'hamming': window = 0.54 + 0.46 * np.cos(2 * np.pi * freqs / n) ramp = ramp * window return np.real(np.fft.ifft(np.fft.fft(proj) * ramp))

注意这里的频率轴处理。np.fft.fftfreq(n)默认返回的频率范围是[-0.5, 0.5]，如果直接乘以它，斜坡的增益会偏小，重建出来的图像灰度偏暗。我一般先乘以n，把频率归一化到和样本数一致的范围，这样能量关系更直觉。窗函数的截止位置也要和n对齐，否则滤波核的宽度不匹配。

做完滤波之后，把每个角度的滤波后投影送入上面的反投影函数，得到的就是FBP重建结果。

4. 实操：用Python跑通平行束FBP全流程

4.1 准备一张测试图像：Shepp-Logan体模

为了验证算法，需要一张有明确结构、有已知正确答案的图像。最经典的是Shepp-Logan体模，它是一个由多个椭圆叠加而成的模拟头部模型，灰度差异明显，内部还有高密度骨骼和低密度软组织区域。在仿真中可以直接生成这张图像，用它计算解析投影，再重建，最后对比。

有一个少有人提的经验：如果只是验证平行束重建算法，用Shepp-Logan体模甚至会“太好”了，因为它光滑、灰度层次分明，滤波参数就算调得不太对，结果也不会太难看。真正会暴露问题的是棋盘类随机图案、带有细线条纹理的测试图，这些对高频信息保留要求极高。我在实际项目里会把三类体模都跑一遍：Shepp-Logan看整体、细线看分辨率、泊松噪声合成投影看鲁棒性。

4.2 模拟投影：把图像转为正弦图

投影可以用解析公式，也可以在离散网格上做射线驱动求和。最常见也最省事的方式是用scikit-image的radon函数：

from skimage.transform import radon image = shepp_logan_phantom_128() theta = np.linspace(0, 180, 180, endpoint=False) sinogram = radon(image, theta=theta, circle=True)

这个过程就是把第二小节讲的Radon变换实际算一遍。circle=True表示在圆形支持域内做投影，可以避免方形图像四角处出现多余的零填充干扰。对这个细节不敏感的人很多，但它在边界处会产生明显差异。

4.3 直接反投影与FBP的结果对比

有了正弦图，就可以分别尝试直接反投影和FBP。直接反投影用上面手写的backproject函数，FBP则先对每一行正弦图做滤波，再反投影。

对比结果非常直观：

• 直接反投影：背景亮蒙蒙一片，高密度区域周围拖出放射状条纹，图像对比度偏低。
• FBP：背景显著变暗，边缘变锐利，高密度区域与低密度区域之间分界清晰。

这个对比几乎做一百次一百次是同样的结论。滤波反投影之所以叫“滤波反投影”，就是因为在反投影之前对投影数据做了高频补偿，抵消了反投影自带的低频放大效应。

4.4 关键参数怎么选：角度数、探测器数、像素数

三个核心参数要一起考虑。

角度数太少，正弦图上沿角度方向采样不足，重建图像会有明显的放射状条纹，这叫有限角度伪影。经验值是：num_angles要大于等于π * N，N是图像边长。用Nyquist采样定理推一下就能理解。比如128×128的图像，大概需要180个角度以上才靠谱。

探测器数和图像像素数要保持一个合理的比例。探测器数太少，投影分辨率不足，重建图像就会糊；探测器数太大于像素数，只会徒增计算量，并不会带来额外细节。一般取N * sqrt(2)或直接取N，在工程上效果都可以接受。

像素数决定图像尺寸，同时也决定计算开销。反投影循环的复杂度是O(num_angles * N^2)，所以图像从256升到512，计算量直接翻四倍。在写原型验证时，用128或256先跑通流程，确认伪影和参数都正常，再放大尺寸。

5. 常见问题与排查技巧实录

5.1 重建图像整体偏移或错位

这一般不是算法问题，而是投影几何参数不一致。常见的原因包括：

• 角度数组写错，比如有的函数要求弧度，你传了度数；
• 探测器坐标中心没有对齐图像中心，导致反投影时所有角度的射线都偏离了真实路径；
• 图像坐标的x/y方向定义和投影生成时不一致，也会让重建结果旋转了90度。

排查看似简单但最容易忽略。我习惯先重建一个非常简单的点源投影，只有一小段锐利的峰值，如果反投影后点源位置正确且没有整体偏移，说明几何没问题，再上真实体模。

5.2 星形伪影非常严重，怎么压

两种情况要分开看。一是在角度数严重不足的情况下，星形伪影是原理性的，再怎么调滤波器都没用，唯一的办法是增加角度数，或者改用迭代重建算法。二是在角度数足够但伪影趋势依然明显的情况下，检查是不是直接反投影没有滤波。如果确认已经做了斜坡滤波但伪影依然重，可以换用Hamming/Hann窗，牺牲一点分辨率换取更干净的背景。

5.3 噪声被放大、图像颗粒感重

斜坡滤波器的高频放大效应在低剂量CT、短曝光工业CT中非常常见。实际项目中我不会看都不看就上Ram-Lak，而是会根据投影信噪比选择窗函数。信噪比较高时用Shepp-Logan，信噪比偏低时用Hamming。如果噪声还是压不住，先看投影数据本身，不要急着上更贵的重建算法。很多时候是探测器响应没有做平坦场校正，或者坏像素没有剔除，这些对断层图的污染会被重建算法放大好几倍。

表格式速查：

5.4 计算太慢，怎么优化

如果只是做原型验证，纯NumPy循环没问题。但如果要处理大尺寸图像或者批量数据，反投影的嵌套循环会成为瓶颈。几个实用优化方向：

• 用内存换时间：把每个角度的s索引预先算好，存成查找表，不要在循环里重复计算三角函数。
• 用并行：把角度维度拆开，分给多核或GPU并行处理。numba加@njit装饰器就能把纯Python反投影加速几十倍。
• 频域重建：投影定理告诉我们，所有角度的投影频谱拼起来就是图像的二维频谱，可以用雷登空间插值方式直接重建，绕开反投影循环。工程上这叫直接傅里叶重建，速度比FBP快，但插值伪影需要小心处理。

我个人最常用的是numba加速。代码改动量极小，就是把backproject函数加一个@njit，把NumPy操作改成显式循环。实测128×128×180角度的重建从秒级降到毫秒级，完全够用。

6. 一点个人体会

写过一遍手写反投影之后，再去看不同成像厂商的重建文档，很多东西一下就顺了。那些所谓的迭代重建、深度学习重建，底层也还是在回答同一个问题：如何在已知投影数据的前提下，尽可能少地被噪声误导、被模糊干扰地恢复物体内部结构。平行束反投影是这一切的起点。

如果你现在也是在调库做重建，我真心建议抽半小时把反投影循环手写一遍。你觉得最没技术含量、最底层的一步，恰恰是理解所有后续算法的地基。踩过一次坑，比看十篇综述都管用。