做 CT 重建绕不开三个词Radon 变换、正弦图、FBP。后面的迭代重建、稀疏视角、扩散模型几乎每篇论文都拿 FBP 当基线。这篇用最少的公式把它们讲清楚最后给一段可以直接运行的 Python 代码亲手看一遍投影 → 正弦图 → 重建的全过程。这是「CT 重建从零开始」系列的第一篇。一、CT 到底测了什么X 射线穿过物体会衰减。沿一条射线 L探测器读到的强度满足 Beer–Lambert 定律I I 0 exp ⁡ ( − ∫ L μ ( x , y ) d l ) I I_0 \exp\left(-\int_L \mu(x,y)\,dl\right)II0​exp(−∫L​μ(x,y)dl)其中μ \muμ是线性衰减系数也就是我们想重建的图像。两边取对数p − ln ⁡ I I 0 ∫ L μ ( x , y ) d l p -\ln\frac{I}{I_0} \int_L \mu(x,y)\,dlp−lnI0​I​∫L​μ(x,y)dl所以 CT 每一个探测器读数经过对数变换后就是图像沿一条直线的线积分。重建要回答的问题是已知大量线积分反推出图像本身。二、Radon 变换把所有线积分排好用角度θ \thetaθ和到中心的距离s ss描述一条直线平行束下的投影就是 Radon 变换p ( θ , s ) ∬ f ( x , y ) δ ( x cos ⁡ θ y sin ⁡ θ − s ) d x d y p(\theta, s) \iint f(x,y)\,\delta(x\cos\theta y\sin\theta - s)\,dx\,dyp(θ,s)∬f(x,y)δ(xcosθysinθ−s)dxdy固定一个角度θ \thetaθ扫过所有s ss得到一条一维投影角度从 0° 转到 180°平行束转半圈就够另一半是重复的就得到全部数据。三、正弦图为什么叫正弦把p ( θ , s ) p(\theta, s)p(θ,s)画成一张二维图横轴角度、纵轴探测器位置就是正弦图sinogram。原因很直观图像里一个点( x 0 , y 0 ) (x_0, y_0)(x0​,y0​)在角度θ \thetaθ下投影到的位置是s x 0 cos ⁡ θ y 0 sin ⁡ θ s x_0\cos\theta y_0\sin\thetasx0​cosθy0​sinθ这正是一条正弦曲线。整幅图像是无数个点叠加所以正弦图就是无数条正弦曲线叠在一起。离中心越远的点曲线振幅越大。四、直接反投影为什么是糊的最朴素的想法把每个角度的投影值沿原路抹回图像所有角度叠加起来叫反投影backprojection。问题在于中心附近每个角度都会被抹到越靠近中心叠加得越多。可以证明直接反投影的结果等于真实图像和1 / r 1/r1/r核做卷积所以低频被严重放大图像一片模糊见下图第三张。五、中心切片定理与斜坡滤波解决方法来自中心切片定理一个角度投影的一维傅里叶变换等于图像二维傅里叶变换中过原点、同一方向的那条切片。也就是说转半圈采集投影相当于在频域里沿放射线采样。放射状采样在中心密、在外围稀密度正比于1 / ∣ ω ∣ 1/|\omega|1/∣ω∣。要补偿这种不均匀就要给每个频率乘上∣ ω ∣ |\omega|∣ω∣这就是斜坡滤波器ramp filter。于是 FBPFiltered Back-Projection滤波反投影只有两步对每个角度的投影做一维斜坡滤波把滤波后的投影反投影回图像。斜坡滤波会放大高频也就会放大噪声。所以实际中常给它乘一个窗函数比如 Shepp–Logan、Cosine、Hamming、Hann牺牲一点分辨率换取更低的噪声。六、代码30 行跑通 FBP只依赖 numpy、matplotlib 和 scikit-imagepip install scikit-image matplotlibimportnumpyasnpimportmatplotlib.pyplotaspltfromskimage.dataimportshepp_logan_phantomfromskimage.transformimportradon,iradon,rescalefromskimage.metricsimportpeak_signal_noise_ratioaspsnr imgrescale(shepp_logan_phantom(),0.64)# 256x256 Shepp-Logan 体模thetanp.linspace(0.,180.,180,endpointFalse)# 180 个角度步长 1°sinoradon(img,thetatheta)# 正弦图探测器位置 x 角度bpiradon(sino,thetatheta,filter_nameNone)# 直接反投影不滤波fbpiradon(sino,thetatheta,filter_nameramp)# 滤波反投影sparsetheta[::6]# 只保留 30 个角度fbp30iradon(radon(img,thetasparse),thetasparse,filter_nameramp)scorelambdax:psnr(img,np.clip(x,0,1),data_range1)print(fFBP 180 views:{score(fbp):.2f}dB)print(fFBP 30 views:{score(fbp30):.2f}dB)titles[Phantom,Sinogram,Backprojection,FBP 180 views,FBP 30 views]fig,axplt.subplots(1,5,figsize(17,3.6))fora,t,xinzip(ax,titles,[img,sino,bp,fbp,fbp30]):a.imshow(x,cmapgray,aspectautoiftSinogramelseequal)a.set_title(t);a.axis(off)plt.tight_layout();plt.show()运行结果scikit-image 0.26180 个角度的 FBPPSNR 30.91 dB和原图几乎一样不滤波的直接反投影一团模糊只能看出轮廓只用 30 个角度PSNR 掉到 22.14 dB满图放射状的条纹伪影。七、换个滤波器会怎样同样 180 个角度、无噪声数据换不同窗函数滤波器PSNR (dB)Ramp30.91Shepp–Logan29.94Cosine28.19Hamming27.14Hann26.85无噪声时纯斜坡最好窗越软越模糊。再模拟一次低剂量。注意体模的灰度只是相对值直接拿它算正弦图线积分最大约 66按e − 66 e^{-66}e−66衰减几乎没有光子能到达探测器。所以先把体模缩小 10 倍让线积分最大约 6.6回到真实 CT 的量级再按入射光子数I 0 10 4 I_0 10^4I0​104加泊松噪声muimg/10# 缩到真实衰减量级线积分最大约 6.6sino_muradon(mu,thetatheta)rngnp.random.default_rng(0)# 固定随机种子结果可复现I01e4# 入射光子数countsrng.poisson(I0*np.exp(-sino_mu))# 探测器计数泊松分布sino_noisy-np.log(np.maximum(counts,1)/I0)forfin[ramp,shepp-logan]:rec10*iradon(sino_noisy,thetatheta,filter_namef)# 乘回 10print(f,f{score(rec):.2f}dB)结果Ramp27.66 dBShepp–Logan27.95 dBShepp–Logan 反超了 Ramp。换随机种子数字会有 0.1 dB 左右的浮动但排序不变。这就是窗函数存在的意义噪声越大越需要压高频。八、FBP 的天花板就是后面所有方向的起点FBP 快、稳定、可解释至今仍是临床 CT 的重要基础。但它对数据的要求很高角度不够稀疏视角条纹伪影就是上面 30 个角度那张图剂量太低噪声被斜坡滤波放大角度范围不足有限角某些方向的边缘直接丢失几何不是平行束实际 CT 多是扇束和锥束需要扇束 FBP 和 FDK 算法。这几个问题分别对应了稀疏视角 CT、低剂量 CT、有限角层析、锥束 CT 几个研究方向。这些方向的代表论文、开源代码和数据集我整理在了一个开源列表里可以按方向筛选、支持中文搜索网站https://computedtomography.org/GitHubhttps://github.com/cslht11/awesome-computed-tomography下一篇「CT 重建从零开始二」会讲迭代重建把 CT 写成线性方程组A x b Ax bAxb从 ART、SART 到加入全变分TV正则看看它们怎么在 30 个角度下把条纹压下去。如果这篇对你有帮助欢迎收藏有问题可以在评论区讨论。