你有没有想过CT 图像并不是传统意义上“拍”出来的反而是“算”出来的探测器真正接收到的只是 X 射线穿过物体后的强度衰减数据至于物体内部哪里密度高、哪里密度低全都要靠数学模型从这些投影数据里反推。CT 技术与图像重建的核心就是把“重建”变成一个求解代数方程组 Axb 的问题而在迭代求解过程中又天然地引入了差分方程递推的思想。这篇文章我结合一个纯 Python 的模拟实验把从投影采集、滤波反投影到代数迭代重建的完整链路拆开讲清楚适合数学建模竞赛选手、图像处理入门者和想理解 CT 重建底层逻辑的开发者参考。1. 从射线到方程CT 重建的第一块积木1.1 投影数据到底是什么先回到物理层。一束强度为 I₀ 的 X 射线穿过物体后强度会衰减为 I这个过程可以用朗伯-比尔定律近似描述I I₀ · exp(-∫ μ(l) dl)其中 μ(l) 是路径上每个位置的线衰减系数积分沿射线路径进行。CT 探测器记录的是 I而我们真正关心的其实是投影值p ln(I₀ / I) ∫ μ(l) dl这个 p 是一个标量吗不是。针对每一个入射角度探测器上会有成百上千个单元每个单元都测量一条独立射线的投影值。把所有角度的投影值按顺序拼在一起就得到一张二维的“正弦图”英文叫 sinogram。之所以叫这个名字是因为任意一个内部点在不同角度下的投影位置会画出一条正弦曲线。所以CT 重建并不是“看图恢复”而是要回答这样一个问题给了我一系列角度下的一维投影值 p如何反推出二维空间里每个位置的 μ 值1.2 离散化把积分变成加权和连续积分在计算机里没法直接求解必须离散化。把待重建的二维区域划分成 n×n 的像素网格每个像素里的 μ 值就是未知数把它们按行展开成一个长度为 N n² 的列向量 x。对应地每条射线沿路径的积分可以近似为路径穿过的每个像素的 μ 值与穿行长度的加权和。于是整条射线对应一个方程∑ w_ij · x_j ≈ p_i把所有角度的所有射线方程堆在一起就得到标准的线性方程组Ax b这里 A 是投影矩阵行数是“角度数 × 探测器单元数”列数是像素数 Nb 是实测投影向量。从纯数学模型角度看图像重建问题到这里已经变成了一个普通的代数方程求解问题。A 的规模通常大得吓人比如一个 128×128 的图像配上 180 个角度、185 个探测器单元A 就有 33300 行、16384 列。如果行数小于列数方程组欠定反过来则是超定。实际中为了稳定性和抗噪声能力通常会布置尽量多的测量数据使系统接近“满定”或“超定”然后通过最小二乘或迭代方法求近似解。1.3 差分方程模型在这里扮演什么角色标题里的“代数方程与差分方程模型”不是生硬拼凑。代数方程模型显然对应 Axb差分方程模型则体现在迭代重建的过程中。无论是 ART 还是 SIRT迭代重建都不是一步到位的求解而是构造一个离散递推序列x_{k1} x_k Δx_k每一次修正量 Δx_k 都是基于当前残差 r_k b - Ax_k 计算出来的。这个公式本质就是一个一阶差分方程。更广义地说对连续 Radon 变换在角度方向和探测器方向上的采样本身也是一种离散化和数值解偏微分方程时构造差分格式的思路完全一致。把图像重建同时理解为“代数的静态求解”和“差分的动态逼近”其实很有价值后面讲迭代算法时你会发现很多反直觉的现象都能从这层理解推出来。2. 三类主流重建思路原理才是根本2.1 中心切片定理所有线性重建方法的地基如果只记一个 CT 理论选中心切片定理。它的表述非常漂亮物体在角度 θ 下的一维投影做傅里叶变换结果等于物体二维傅里叶变换在过原点、方向角 θ 的直线上的切片。翻译成人话就是每一个角度的投影频率信息都只是二维频谱里的一根线。想要完整还原物体理论上需要无穷多个角度的投影把频域平面铺满。实际只能采集有限角度所以频域某些区域信息不足这就为伪影埋下了伏笔。这个定理同时解释了为什么重建问题本质上是困难问题因为我们永远只能拿到有损采样后的频域信息。2.2 直接反投影为什么图像是糊的最直觉的重建想法是反投影把某个角度的投影值沿着入射方向均匀地“抹回”图像平面把所有角度叠加起来。这个过程数学上相当于对每个像素把所有穿过它的射线投影值求和。听起来合理结果却不忍直视重建出来的图像被严重模糊物体边缘会有放射状星形伪影均匀区域出现杯状畸变。原因是反投影本身并不是 Radon 变换的逆运算它其实是在频域里给频谱乘了一个 1/|ω| 的权重因子也就是说高频信息被系统性削弱了。用一个生活类比想这件事反投影像把几张“月光投影”叠在一起每张投影都只是模糊的大致轮廓叠完了还是模糊而且多叠几次还会产生杂乱的星芒。2.3 滤波反投影补偿被削掉的高频既然反投影在频域乘了 1/|ω|那重建前先把投影在频域乘回 |ω| 不就行了这就是滤波反投影 FBP 的核心思路。具体分三步对每个角度的一维投影做傅里叶变换在频域乘一个斜坡滤波器也就是幅值响应为 |ω| 的滤波器逆傅里叶变换回空间域得到滤波后的一维投影再做反投影。斜坡滤波器是理论上的精确选择但直接把 |ω| 用起来会把高频噪声放大到离谱。所以工程上几乎总是要加窗函数做截断最常用的是 Hamming 窗和 Hann 窗。加了窗之后高频部分被压下来图像更光滑代价是分辨率略有下降。FBP 的优势是快、稳定只要有 sinogram 就能在线性时间内重建所以传统 CT 设备里长期都是主力算法。2.4 迭代重建绕开频率域直接解方程FBP 建立在频域分析之上而迭代重建选择直接面对代数方程组 Axb不搞傅里叶变换。最经典的代数重建技术 ART 是 Kaczmarz 方法的直接应用。思路是每次只拿一条射线的方程做修正x ← x λ · (p_i - aᵢᵀx) / ‖aᵢ‖² · aᵢ其中 aᵢ 是 A 的第 i 行λ 是松弛因子。翻译成直觉语言看当前解对这条射线的预测值和实际测量差多少然后把误差按照该射线穿过像素的权重大小“分摊”回像素上去。SIRT 则是 ART 的温和版本每轮迭代把所有射线的修正量累积起来取平均后再统一更新一次。因为不逐线突变SIRT 的抗噪声能力更强收敛也更平稳但需要更多轮迭代。这类方法最大的优势在于即使投影角度很少、测量不完全迭代法仍然能给出比 FBP 更合理的结果因为频率域里缺失的信息在迭代过程中会被空间域约束补回来一些。这也是低剂量 CT 和稀疏角度 CT 里迭代方法大受欢迎的根本原因。3. 用 Python 实现一个完整的最小 CT 重建闭环纸上谈兵没意思直接上代码。我用 numpy 和 scipy 造一个极简但完整的模拟实验生成模拟体模、构造投影矩阵、计算正弦图、跑 FBP 和迭代重建、最后对比质量。整个流程大概一百来行。3.1 准备模拟数据叠加椭圆的体模真正的 CT 验证常用某个公开的合成体模里面用一堆椭圆模拟不同组织的密度。这里我手动生成一个简化版直接叠三个区域。它不涉及任何现实病人数据纯粹是数学实验。import numpy as np from scipy.sparse import csr_matrix from scipy.sparse.linalg import lsqr def make_phantom(size128): yy, xx np.mgrid[0:size, 0:size] cx, cy size / 2.0, size / 2.0 img np.zeros((size, size), dtypenp.float32) # 背景 img[(((xx - cx) / (size * 0.45)) ** 2 ((yy - cy) / (size * 0.45)) ** 2) 1.0] 0.2 # 一个大椭圆 img[(((xx - cx) / (size * 0.22)) ** 2 ((yy - cy - size * 0.14) / (size * 0.16)) ** 2) 1.0] 0.9 # 一个小圆 img[(xx - cx - size * 0.16) ** 2 (yy - cy size * 0.12) ** 2 (size * 0.04) ** 2] 0.6 return img生成出来是这样一块“灰色背景 亮色椭圆 亮色圆点”的合成物。内部灰度差异足够大重建质量一眼就能看出来。3.2 构造投影矩阵最简单的权重模型要跑迭代重建光有正弦图还不够我需要 A 矩阵。正规工程里会用射线逐像素求交的方式精确计算权重典型代表是 Siddon 算法。演示代码里我用一个更粗暴的近似把一个像素中心点投影到探测器的一维坐标上落在某个探测器 bin 的宽度范围内就认为这条射线对该像素有权重权重直接取 1。def build_weight_matrix(size, angles_deg, ndet, det_len): rows len(angles_deg) * ndet cols size * size cy cx size / 2.0 bins np.linspace(-det_len / 2.0, det_len / 2.0, ndet) bin_w (det_len) / ndet ri, ci, val [], [], [] for ai, ang in enumerate(angles_deg): theta np.deg2rad(ang) for di, bpos in enumerate(bins): r ai * ndet di for j in range(size): for i in range(size): x i 0.5 - cx y j 0.5 - cy s x * np.cos(theta) y * np.sin(theta) if abs(s - bpos) bin_w / 2.0: ri.append(r) ci.append(j * size i) val.append(1.0) return csr_matrix((val, (ri, ci)), shape(rows, cols)), bins注意两点。第一这个三重循环在 128×128 下能跑但已经偏慢实际工程中一定会向量化或用射线追踪。第二权重值统一取 1 会损失定量精度但定性演示完全够用后面你可以看到伪影特征和理论预期完全一致。3.3 正投影与带噪正弦图有了 A 矩阵正投影就是一次稀疏矩阵乘向量。size 128 angles np.arange(0, 180, 2) # 90 个角度 ndet 185 det_len size * 1.5 phantom make_phantom(size) A, bins build_weight_matrix(size, angles, ndet, det_len) x_true phantom.ravel() b_clean A x_true rng np.random.default_rng(42) noise rng.normal(0, 0.03 * b_clean.std(), sizeb_clean.shape) b b_clean noise这里我加了 3% 标准差的高斯噪声。不加噪声也能重建但加了才真实因为真实探测数据一定是有噪声的而且 FBP 对噪声的放大特性很典型。3.4 手写滤波反投影直接用某个图像库的现成函数当然可以但既然我们要聊原理最好自己写。滤波过程很简单def ramp_filter(p, det_pitch): n p.size freq np.fft.rfftfreq(n, ddet_pitch) pf np.fft.rfft(p) * np.abs(freq) return np.fft.irfft(pf, nn)反投影部分我直接用 A.T 来实现。理解这张矩阵A 的行是“射线”列是“像素”A.T 就是把一条射线的值加到它经过的那些像素上正是反投影操作。把每个角度的投影滤波后乘到对应行上再由 A.T 聚合回图像空间。def fbp_reconstruct(A, b_sino, angles, ndet, det_len): ncol A.shape[1] recon np.zeros(ncol) det_pitch det_len / ndet for r in range(len(angles) * ndet): row_vec A.getrow(r).toarray().ravel() if row_vec.sum() 0: continue di r % ndet filtered ramp_filter(b_sino[r - di : r - di ndet], det_pitch) recon filtered[di] * row_vec return recon.reshape(size, size)严格来说这里每个角度内部逐线滤波再反投影是可行的只是效率不高。跑出来的结果会看到典型问题边缘比直接反投影清楚得多但噪声被抬高了圆点边缘出现亮暗交替的振铃。3.5 ART 和 SIRT 迭代求解ART 的核心循环def art_reconstruct(A, b, size, n_epoch20, lam0.1): x np.zeros(A.shape[1]) m A.shape[0] for _ in range(n_epoch): for i in range(m): ai A.getrow(i).toarray().ravel() denom ai ai if denom 1e-12: continue resid b[i] - ai x x x lam * resid / denom * ai return x.reshape(size, size)逐线更新的好处是每一轮都在快速修正前几轮误差下降极快。缺点是噪声大的时候容易过拟合λ 调太大还会发散。我试过 0.3 以上时结果剧烈震荡0.05 到 0.1 比较稳。SIRT 则换成整轮更新def sirt_reconstruct(A, b, size, n_iter200, lam0.5): x np.zeros(A.shape[1]) row_sum np.asarray(A.sum(axis1)).ravel() col_sum np.asarray(A.sum(axis0)).ravel() row_sum[row_sum 1e-12] 1.0 col_sum[col_sum 1e-12] 1.0 for _ in range(n_iter): r b - A x update A.T (r / row_sum) x x lam * update / col_sum return x.reshape(size, size)SIRT 每轮的计算量比 ART 大但每一轮都不是“过冲式”修正而是把各路信息平均后再走一小步所以曲线更平滑。在噪声比较大的情况下SIRT 重建出来的图像往往最干净。3.6 结果对比不要只用眼睛看肉眼判断太主观我再算两个常用指标相对误差和结构相似度 SSIM。用一个小函数def rel_err(pred, true): diff pred - true return np.sqrt(np.sum(diff ** 2)) / np.sqrt(np.sum(true ** 2))在我的模拟配置下典型结果是这样方法相对误差SSIM直观感受直接反投影偏高偏低明显模糊有星形伪影FBP较低较高边缘清晰但有噪声振铃ART 20 轮中间较高边缘尚可噪声被放大SIRT 200 轮最低最高平滑伪影少实际跑的时候你会注意到直接反投影虽然难看但它其实是 FBP 的“未滤波版本”两者的计算速度都极快。迭代算法在 CPU 上跑几百轮要好几秒到十几秒换回来越强的抗噪声能力和稀疏角度适应性。这里没有绝对最优只有取舍。4. 从模拟到工程几个必须注意的坑4.1 角度采样不够一切伪影的根源模拟中我把角度步长从 2 度改到 6 度FBP 图像立刻出现大面积的放射状条纹看起来像从几个方向向四周散开的光线。这是因为频域里的采样不够密高频信息沿某些方向完全缺失。实际项目里角度数和探测器单元数是直接决定重建上限的硬约束算法再好也补不出信息都没采集进来的频率内容。所以拿到一个重建任务先问一句投影数据质量如何、角度覆盖够不够。忽略这个前提去调算法纯属缘木求鱼。4.2 A 矩阵的内存爆炸与稀疏化如果不注意数据结构直接用稠密矩阵存储 A128×128 图像、33300 条射线时矩阵尺寸是 33300×16384每个 float64 占 8 字节总内存约 4.4GB。这还只是一个小规模的演示。但理论上 A 里绝大多数元素都是 0因为一条射线最多穿过几十个像素所以正确做法是压缩稀疏行格式。上面代码里用的稀疏矩阵就是这个原因。工程上投影矩阵的构建还会继续优化用射线追踪遍历、只在像素边界处求交、只生成有效射线对应的行能再省一个数量级的时间和内存。4.3 迭代发散松弛因子不是越大越快我最初跑 ART 时为了“快速收敛”把 λ 直接设成 0.5结果前几轮误差确实降得很快但后面突然开始震荡最终图像出现黑白相间的棋盘状噪声。这是典型的发散现象。松弛因子的作用相当于差分方程递推里的步长。步长太大残差修正会跨过真解步长太小收敛慢到让人失去耐心。经验值上ART 的 λ 在 0.050.3 之间比较安全SIRT 的 λ 可以到 0.51.0因为它是平均更新天然更稳。判断是否发散不要只看误差绝对值要看每轮误差曲线的末端趋势是下降还是反弹。4.4 滤波器选择噪声放大与分辨率损失的博弈FBP 里斜坡滤波器的窗函数选择是个经典权衡。用纯斜坡滤波图像最锐利但投影里有噪声时会出现密密麻麻的高频毛刺。用 Hann 窗噪声被压住了但小细节也跟着糊了。工程上的做法通常是先估算投影数据的噪声水平噪声高就选更保守的窗噪声低可以激进一些。还有一个技巧是同时显示滤波前后的一维投影曲线如果投影曲线在物体边缘处看着已经有明显抖动那 FBP 出来基本没法看这时候应该转身考虑迭代法。4.5 常见问题速查表症状可能原因处理方案图像整体模糊、边界发虚高频补偿不足可能是滤波窗口太保守改用更宽的滤波器或减少平滑星形放射状伪影角度数不够频域缺方向采样增加角度数或用迭代法图像布满细密噪声斜坡滤波器放大了高频噪声加窗滤波或改用 SIRT迭代后期误差反弹松弛因子过大发散调小 λ或换 SIRT重建结果区域亮度偏暗反投影或迭代缺少归一化检查权重矩阵行和列归一化内存溢出A 用了稠密矩阵改用稀疏矩阵或流式逐角度重建我自己的习惯是小规模实验先跑全流程确认正投影和反投影配对正确正式处理带噪声数据时优先上 SIRT用 200 轮左右的迭代对比 FBP 的结果再下结论。经典方程模型的价值在于它让整条链路保持透明哪里效果好、哪里效果差全都可以追溯到矩阵和频率的细节上。做完这一次完整模拟我对“CT 重建是一个病态反问题”这句话算有了实打实的理解。这一期的“代数方程与差分方程模型一”先聊到这儿CT 图像重建里还有动态序列和时空递推的场景那就是差分方程模型进一步延伸的方向了。