简介这套压缩资料面向CT成像与图像重建领域的学习者和研究者尤其适合正在学习扇形束扫描条件下滤波反投影FBP重建算法的读者。包内共2个文件含1个MATLAB脚本和1篇英文PDF论文压缩包仅345KB。脚本实现了FBP重建流程投影数、源到探测器距离、探测器大小、重建矩阵尺寸等关键参数均可调整这些参数共同决定重建图像的空间分辨率和伪影水平便于对比不同采样设置下的质量差异PDF论文探讨计算层析成像的实现为理解算法提供理论支撑。目前已有247人浏览学习适合医学影像、无损检测或算法课程设计、毕业课题参考。通过这份资料读者既能直接运行算法脚本又能结合论文掌握扇形束投影几何与FBP实现细节快速建立从原始投影到断层图像重建的完整知识链。1. 扇形CT的FBP重建标题里的这批数据该按什么流程吃掉看到data_扇形FBP_ct数据_扇形CT_giftcja_FBP这种命名第一反应不是去搜解析库而是先确认里面装的到底是什么。FBP滤波反投影在 CT 里是最成熟的重建路线拿到一组在不同角度下扫描的投影经过滤波和反投影叠加还原出断层上的线性衰减系数分布。扇形CT 指的是射线源和探测器组成的扇形束几何工业CT、低剂量CT实验里最常见的投影数据就是这种。后面那串giftcja_FBP大概率是某台设备或某个批处理脚本自动生成的标记真正要紧的是三件事投影数据排布、扇束几何参数、旋转中心偏移量。这篇笔记就讲这三件事怎么在动手重建前落实以及扇形束 FBP 和平行束 FBP 差在哪。2. 扇形束几何与 FBP 公式为什么这段数据不能丢给平行束算法2.1 先分清是等角还是等间距采样拿到扇形 CT 数据第一个要确认的是探测器阵列的采样方式。探测器装在弧面上每个探测单元的张角相同这种叫等角采样探测器装在平板上单元间距相同叫等间距采样。医用 CT 和大部分工业 CT 用的是等角采样因为扇形束重建公式在这种几何下最干净。但实验室自研设备、微焦斑 CT 和一部分低成本系统会用平板探测器得到的就是等间距数据。这两种几何在 FBP 里的校正方式不一样直接混用会出现边缘扭曲。等角采样下每个探测器通道对应一个固定的扇形角 γ重建时只需要在滤波前乘一个cos(γ)的预加权等间距采样下投影坐标和扇形角之间是t D * tan(γ)的关系滤波核和反投影权重都要换成更复杂的表达。采样方式探测器排列预加权反投影权重实现难度等角弧面等角排列cos(γ)1 / L²低等间距平面等距排列需重排为等角需变尺度反投影中高我一般拿到数据先看探测器数量、投影角度数量再看有没有配置文件标注 detector pitch、source-detector distance。如果完全没配置就重建一帧看看边缘扭曲方向也能反推是哪种采样。等间距数据最省事的做法是在重建前把探测器坐标按gamma atan(t / D)重新插值到等角网格再用等角 FBP 跑损失一点边缘分辨率换来整个流程稳定。2.2 FBP 的三步几何校正预加权、滤波、反投影加权重扇形束 FBP 与平行束的核心差别在反投影前多了一次几何加权。标准流程分三步第一步预加权。同一个物体位置在不同扇形角 γ 下射线穿过物体的体积元长度不同需要把投影值乘上cos(γ)做校正。这一步不只在边缘通道明显探测器张角越大越不能省。第二步滤波。沿探测器方向做一个一维斜坡滤波对应频域里乘|ω|。CT 投影在探测器方向上低频强、高频弱直接反投影会出现星形条纹斜坡滤波相当于把高频抬起来补偿角采样不足。实际工程里很少用纯 R-L 核都会叠一个 Hamming 或 Hann 窗来压高频噪声低剂量CT图像尤其依赖这一步。第三步加权反投影。对每个投影角度 β把滤波后的投影值沿射线方向铺回图像空间。扇形束里图像上某像素对应的投影位置不是由平行坐标决定而是由“源到该像素的连线与中心射线的夹角”决定。反投影时给每个像素乘上1/L²L 是射线源到该像素的距离。离源越远的区域衰减越强不校正的话图像四周会发暗。离散化之后反投影表达式可以写成f(x, y) (Δβ / 2π) * Σ_β [ 1 / L²(x, y, β) * q(γ_p, β) ]其中q是预加权并滤波后的投影γ_p是像素 P 对应的扇形角。这个公式只对全扫描0 到 2π成立如果是短扫描还需要加 Parker 加权那一套更适合放到迭代重建或者已有成熟工具库的环节来讨论手写从零实现时先按全扫描处理。2.3 从探测器计数到 sinogram负对数与空气校准FBP 的输入必须是“线积分”也就是物体对 X 射线衰减的路径积分单位与线性衰减系数 μ 一致。探测器读到的原始量是光子数衰减后的强度两者之间隔着比尔定律p -ln(I / I0)I0 是无物体时的空气投影强度I 是实际扫描强度。如果给你的数据已经是处理好的 sinogram数值范围大概在 0 到 4 左右如果原始文件里是大几百甚至几万的值那很可能是光子计数或原始探测器响应。这里最容易翻车有人把光子计数直接塞进 FBP重建出来全是负值云团。空气校准的另一个作用是消除探测器单元灵敏度差异。每个探测器的增益不可能完全一致等角排列的弧面探测器尤其明显。处理办法是采集一帧空气投影 air_scan用实际投影逐通道除以它再做负对数p(x) -ln( data[x] / air_scan[x] )这一步不做重建图里会出现一圈一圈的同心环伪影。工业CT 里很多设备自带平场校正但第三方数据往往没做过养成拿到数据先看“有没有空气帧”的习惯。3. 搭一套可复现的扇形 FBP 重建Python 代码与参数对照3.1 读入投影文件reshape 成角度探测器并做输入自检拿到一批未知的 CT 投影文件第一步是确认存储格式。常见的 .raw 是 float32 或 uint16 连续二进制文件名里一般带n_proj和n_det。下面这段代码把二进制读成二维数组并做了最基础的数值自检import numpy as np n_proj 360 # 投影数量 n_det 512 # 探测器通道数 dtype np.float32 # 多数工业CT投影文件是float32 raw np.fromfile(proj_001.raw, dtypedtype) sino raw.reshape(n_proj, n_det).astype(np.float64) # 数值自检直接判断输入是否已经是对数后的线积分 print(min/max:, sino.min(), sino.max()) print(mean:, sino.mean()) # 如果原始值是光子计数通常范围很大且全为正 # 这里需要先做空气校准再取负对数 # air np.fromfile(air.raw, dtypedtype).reshape(n_proj, n_det) # sino -np.log(sino / np.maximum(air, 1e-6))逻辑说明reshape(n_proj, n_det)的维度顺序必须确认。大多数数据存储习惯是“先角度后探测器”但有的设备文件头会反着写最好先把第一帧sino[0]画出来如果沿横轴是一条连续曲线而不是一个周期信号说明维度顺序反了。min/max是判断数据类型的快速手段线积分的值一般接近 0 到 5如果最大值在几百以上大概率还是原始计数。遇到这类情况先做负对数转换不要直接重建。3.2 最小实现等角扇形 FBP 重建函数下面是一个完整的等角扇形束 FBP 实现包含 R-L 滤波核和可选 Hamming 窗不依赖任何 CT 专用库只用 NumPy。体数据的地面是 0.01 量级、目标几何大小是几十个像素重建尺寸取 128 足够看效果。def rl_filter(proj, d_gamma, windowhamming): 对每个角度的投影沿探测器方向做一维R-L滤波。 proj: (n_proj, n_det) d_gamma: 相邻探测器通道的扇形角间距弧度 n_proj, n_det proj.shape n_idx np.arange(-(n_det - 1), n_det) kernel np.zeros_like(n_idx, dtypenp.float64) even_mask (n_idx % 2 0) kernel[even_mask] -1.0 / (np.pi**2 * n_idx[even_mask]**2) kernel[n_idx 0] 0.25 kernel / d_gamma if window hamming: kernel * np.hamming(len(kernel)) out np.zeros_like(proj) for i in range(n_proj): out[i] np.convolve(proj[i], kernel, modesame) return out * d_gamma def fan_fbp_eqang(sino, betas, gammas, R, pix_size1.0, npix128): 等角扇形束FBP重建。 sino: (n_proj, n_det)已做负对数的线积分 betas: 每个投影对应的旋转角弧度长度n_proj gammas: 探测器扇形角弧度长度n_det R: 射线源到旋转中心的距离 n_proj, n_det sino.shape d_gamma abs(gammas[1] - gammas[0]) # 1) 预权重补偿扇形角导致的路径长度差异 pw sino * np.cos(gammas)[None, :] # 2) 斜坡滤波 fproj rl_filter(pw, d_gamma, windowhamming) # 3) 加权反投影 pix_axis (np.arange(npix) - (npix - 1) / 2.0) * pix_size X, Y np.meshgrid(pix_axis, pix_axis) img np.zeros((npix, npix), dtypenp.float64) for k, beta in enumerate(betas): sx R * np.cos(beta) sy R * np.sin(beta) dx X - sx dy Y - sy L2 dx * dx dy * dy # 像素射线与中心射线的夹角扇形角 gamma_pix np.arctan2(np.sin(beta) * dx - np.cos(beta) * dy, -np.cos(beta) * dx - np.sin(beta) * dy) # 扇形角坐标 - 探测器通道坐标 pos (gamma_pix - gammas[0]) / d_gamma pos np.clip(pos, 0, n_det - 1) i0 np.floor(pos).astype(np.int32) i1 np.minimum(i0 1, n_det - 1) w pos - i0 val (1.0 - w) * fproj[k, i0] w * fproj[k, i1] img val / L2 # 离散角度积分归一化 d_beta abs(betas[1] - betas[0]) img * d_beta / (2.0 * np.pi) return img逻辑说明rl_filter里构建的是经典 R-L 离散卷积核偶数位置按-1/(π²n²)取值零点取1/4奇数位置为零。除以d_gamma是为了让卷积核的物理尺度跟上相邻探测器角间隔。Hamming 窗直接乘在核上作用是压住高频振铃低剂量 CT 数据里这行基本不能省。fan_fbp_eqang的反投影部分用了像素驱动方式对重建图像上每个像素算出从源到该像素连线的扇形角再用线性插值取对应滤波后的投影值。L2是源到像素的欧氏距离平方直接作为反投影权重。如果重建结果是左右翻转的把gamma_pix的计算取负号即可这是等角扇形几何的方向约定问题不影响算法正确性。3.3 跑通需要的三组参数角度、扇角、旋转半径没有真实的 CT 投影数据时可以用下面这段生成一个简单 phantom 投影作用等价于一个可复现的测试夹具用来验证重建代码def make_phantom(n128): ph np.zeros((n, n)) rr, cc np.ogrid[:n, :n] cx, cy n / 2.0, n / 2.0 ph[(rr - cy)**2 (cc - cx)**2 28**2] 0.02 ph[(rr - (cy 10))**2 (cc - (cx 15))**2 10**2] 0.01 return ph这个 phantom 是两个不同密度的圆重建后水区应当是平滑圆盘、高密度区灰度更高。如果你手头有真实工业 CT 数据直接读入替换sino然后按扫描协议的几何设置参数参数常见取值若不匹配会出现什么n_proj360~1440太少会出放射状条纹betas0 到 2π 均匀分布角度不均匀会出半圆错位gammas范围±0.2~±0.5 rad范围越宽 FOV 越大边缘畸变也越大R设备源到旋转中心距离mm偏大/偏小会让图像整体放大缩小pix_size0.1~1.0 mm设置过大图像糊过小只是放大噪声实际调试时我会先只用 64×64 重建尺寸跑通流程确认几何参数合理后再开全分辨率。R是最容易搞错的参数很多配置文件只给 source-to-detector distance没有给 source-to-rotation-center。后者需要标定或从数据自带的校准扫描里反推不要直接从文件名猜。4. 扇形 FBP 重建的现场排查五个高频坑与处理办法4.1 重建图像边缘重影旋转中心偏移现象重建出的圆盘边缘有两条重叠边界像是画框内又套了一个框高密度物体周围还有对称振铃。原因扇形束数据里旋转中心与重建图中心不重合。等角扇形束比平行束对中心偏移更敏感差半个探测器通道都会在边缘看到重影。很多设备虽然出厂做了标定但机械运动、热胀冷缩都会改变实际旋转中心位置。解决给反投影的探测器坐标加一个整体偏移。修改fan_fbp_eqang里pos的计算把pos减去center_offset_det这个偏移可以取负也可以取正取决于偏移方向。离线标定最土也最可靠的办法是扫描一根细金属丝重建后看它是否变为一个圆点。如果不是调整偏移量直到伪影最小。数值上可以试-2.0到2.0之间的 0.1 步长。4.2 图像出现同心环条纹探测器响应不一致现象重建图像里出现一圈一圈的同心圆尤其在空气区域最明显像树的年轮。原因探测器中某个或某几个通道的增益和别的通道不一致。这个不一致在 sinogram 里是一条固定的竖线反投影后沿不同半径扫出同心圆。工业 CT 的探测器老化、某个模块过热都是典型来源。解决先把 sinogram 画出来看是否存在固定位置的高/低亮竖线。找到坏通道后用相邻通道的均值替换bad_idx [37, 238] # 从sinogram图中定位 for i in bad_idx: lo max(0, i - 3) hi min(n_det, i 4) sino[:, i] np.median(sino[:, lo:hi], axis1)替换前先对整列做中值滤波效果更好但要注意不要误伤正常的低剂量 CT 图像边缘结构。空气扫描的平场校正是预防手段坏通道替换是事后补救两者都做才稳妥。4.3 低剂量 CT 投影噪声爆炸R-L 滤波核太尖锐现象低剂量 CT 数据重建后图像颗粒感极强像撒了一层盐边缘甚至出现彩色噪点。原因投影噪声在频域里是全频段分布的R-L 斜坡滤波把高频噪声放大。低剂量扫描的入射光子少泊松噪声占比高纯 R-L 核完全不可用。这不是算法写错了是滤波器和噪声不匹配。解决把rl_filter里的窗函数列成一个开关优先试 Hamming 和 Hann。Hamming 窗在滤波器边缘保留约 8% 增益Hann 窗衰减更彻底。还可以手动限制滤波核的有效截止频率比如把核截断到中心零点两侧各 16 个采样点。代价是空间分辨率下降但低剂量场景下人眼主观质量往往更好通常需要结合分辨率测试来决定取舍。4.4 整张图发糊投影欠采样或截止频率太低现象重建图像轮廓还在但细节全无圆盘边缘过渡带很宽。用细丝 phantom 测试时FWHM 偏大。原因两个来源。要么投影数量太少角度方向的信息量不足要么滤波核太长或窗函数太狠高频成分被压过度。区分方法是看 sinogram如果相邻角度的投影之间变化剧烈说明角度采样不够。如果 sinogram 平滑但重建糊问题在滤波参数。解决优先保证n_proj和探测器数量匹配经验公式投影角度数应不小于 π × 单角探测器覆盖范围的一半工程上常见取n_det * π / 2左右。例如 512 通道扇形束全扫描 360° 至少需要约 800 个投影角度720 到 1440 很常见。滤波核方面把窗函数换成 Hamming 且不截断核重建一张对比如果锐度回来了问题就在核截断。4.5 图像出现反白负值云团输入数据取没取负对数现象重建结果里物体区域是负的、空气区域是正的或者出现大片黑色空洞CT 值完全不正常。原因输入给 FBP 的 sinogram 不是线积分。常见坑有两种一是原始数据是光子计数直接送进了 FBP二是数据本身已经做过负对数又被重复取了一次对数。反复取对数会把弱吸收区域推到负值。解决重建前先检查数值范围。如果整个 sinogram 最大值小于 20基本都是已处理好的线积分直接进 FBP如果最大值是几万先归一化到I/I0再取-log。如果疑似重复取了对数把图像整体乘-1看看灰度关系是否恢复正常。这类问题没有通用修复只能在数据加载阶段多打几个print养成习惯。5. 重建质量验收CT 值、分辨率、噪声三个指标一起看5.1 CT 值线性空气和水模定量重建完成后先不急着调好看放三个圆形 ROI空气区域、均匀水模区域、高密度校准块。分别统计均值和标准差。空气区域的理论值接近 0或 -1000 HU看设备刻度约定水模应落在 0 左右。如果水模均值和标称值偏了几十个 HU 以上多半是滤波核常数或反投影权重没校准先把重建里那个d_beta / (2π)系数换成体模实测刻度。这个校准不只是为了“显示好看”。工业CT 里靠灰度值判断缺陷深度如果线性区间偏了后续所有定量分析都不可信。低剂量 CT 研究里更是如此散射伪影、束硬化伪影都会在 ROI 统计里暴露出来。5.2 空间分辨率的快速验证线对卡与点扩散没有线对卡时重建一根细丝 phantom计算重建图像的径向剖面取峰值半高宽来做分辨率估计。理想情况下 FWHM 应接近探测器单元投影到中心的表现尺寸。如果 FWHM 明显偏大检查像素尺寸是不是定得太粗。像素尺寸和经验值相差 2 倍以上时图像分辨率的瓶颈在重建网格而不是滤波核。注意不要为了追求高分辨率无限减小像素尺寸像素小于探测器几何极限后只会放大噪声FWHM 并不会再下降。5.3 滤波函数选择的最终权衡噪声与分辨率先定一个低剂量 CT 图像里分辨率噪声永远是矛盾关系。我的习惯是先固定噪声目标在均匀水模区域量出标准差要求它不超过某个阈值达到目标后看线对卡分辨率还能剩多少再微调 Hamming 窗截断长度。先噪声后分辨率的顺序适合临床和工业检测场景反过来做会陷入无休止的调参。这套扇形 FBP 流程我前后改过不下十次每次翻车都出在最基础的输入数据上中心偏移没标定、负对数处理错误、把等间距数据当等角处理。希望帮到你。本文还有配套的精品资源点击获取