简介面向计算机断层成像学习者的扇形束滤波反投影重建资源包聚焦滤波反投影算法在扇形CT数据中的参数设置与实现。共2个文件包含一个MATLAB脚本与一篇PDF论文压缩包整体约345KB。MATLAB脚本可用于实践投影、距离、探测器大小、重建矩阵大小等参数对重建结果的影响PDF文献则围绕计算层析成像的实现展开辅助理解滤波反投影在扇形束及类CT场景中的理论依据与算法细节。目前已有247人学习下载适合医学物理、无损检测及相关方向的学生和研究人员快速上手扇形CT重建实验。通过脚本调试与论文对照读者既能获得可直接运行的算法示例也能掌握从投影数据到断层图像的完整处理思路为后续研究或工程应用打下基础。1. 扇形CT重建成像明明拿到的是数据为什么第一步卡在几何标定拿到一份文件名带“giftcja”的扇形CT数据打开以后不是图像而是一堆按角度排列的探测器读数时很多人的第一反应是找现成软件直接出图。实际做过一遍就有个反直觉的结论扇形CT数据重建成像最耗时间的不是FBP算法本身而是把几何参数标对。D_so源到旋转中心、D_sd源到探测器、探测器间距、起始角这四个数里任何一个差零点几毫米或零点几度重建出来的切片上就是一层层重影和星芒算法再先进也拉不回来。这篇笔记面向的是工业CT检测、科研CT重建、低剂量CT图像处理这条链上的工程师目标是让你拿到类似giftcja这样一套扇形CT投影数据后能自己写通从数据读取、几何重排、FBP重建到验证的一条完整链路。2. 扇形FBP与平行束FBP的本质差别多出来的加权项和一套几何映射2.1 扇形CT与平行束CT的差别一个角度参数带来的整套公式变化先想清楚一件事CT重建教科书里最常见的平行束FBP输入是一个(sinogram)横轴是探测器单元纵轴是旋转角度。而扇形CT的投影形状完全不同X射线从一个点源出发覆盖一个扇形面探测器是圆弧排列或者等距直线排列。同样是转一圈采集射线并不是平行穿过物体的每条射线与旋转中心的几何关系都要用扇形角来描述。扇形束和平行束最核心的换算关系是任意一条扇形射线都可以用两个量表示——射线与中心射线的夹角γ以及这条射线到旋转中心的距离t。平行束重建时我们直接按(t, θ)来做滤波反投影而扇形CT里探测器直接给的是角度β和探测单元位置u所以多了一步把扇形投影“重排”成平行束投影或者直接在反投影公式里加距离权重。这步不是可选项。如果拿着扇形数据直接套平行束FBP公式重建出来的图像会有明显的杯状伪影和几何畸变靠近视野边缘的地方误差尤其大。原因在于平行束公式默认每条射线的路径长度权重一样而扇形射线从点源出发越远离中心射线的射线在物体内走过的路径和到达像素的距离权重都不一样。2.2 扇形数据用FBP建的三步重排、斜坡滤波、反投影扇形FBP最常见的工程做法不是硬写扇形反投影公式而是先把扇形投影重排成平行束sinogram再走标准的平行束FBP流程。这种做法稳定、好调试而且很多开源CT重建框架在二维场景下也是这么干的只是把重排藏在了底层。重排的核心是建立探测器单元与平行束坐标的映射。以等距直线探测器为例设源到旋转中心距离为D_so源到探测器距离为D_sd探测单元i到中心射线的距离为u那么这条射线与中心射线的夹角γ满足γ arctan(u / D_sd)这条射线到旋转中心的距离t为t D_so * sin(γ)如果采集投影时的源角度是β那么这条射线对应平行束sinogram的角度θ为θ β γ把每条射线按算出的(t, θ)散落到一个规则的(t, θ)网格里中间用线性插值补齐就得到了平行束格式的sinogram。这个过程叫fan-to-parallel rebinning是所有扇形FBP落地里最值得自己手写一遍的代码。重排之后做滤波滤波核仍然用斜坡核Ram-Lak就是频域里乘一个|f|。如果投影数据本身噪声偏大可以给斜坡核加窗最常见的是Hamming窗或Cosine窗。这里有个容易忽略的细节斜坡滤波是作用在探测器坐标方向的不是在角度方向滤波器方向搞反了重建结果会变成一片横向条纹状的模糊。最后是反投影。反投影这一步会把滤波后的每个角度投影值沿着对应角度的射线方向均匀铺回图像空间把所有角度的贡献叠加起来。离散实现里每个角度的投影先延拓成一张二维的“竖直条纹图”再把条纹图旋转到该投影角度逐角度累加。拼完所有角度后除以角度总数乘一个归一化系数得到的就是体密度图像。2.3 为什么工业CT、低剂量CT这些场景仍把FBP当基线现在深度学习重建、迭代重建都很多但工业CT和低剂量CT图像处理里FBP仍然是基线算法原因很实际第一FBP是线性的参数固定后同一套数据每次跑出来的结果完全一致做缺陷检测和尺寸测量时可重复性比迭代算法好第二FBP没有迭代步数、正则化强度这些玄学参数最坏情况下图像只是糊一点不太会出现迭代重建那种“硬生生收敛出假结构”的情况第三工业CT扫描的物体往往密度对比很大比如金属工件里看气孔FBP的线性特性让灰度值能直接和线衰减系数挂钩。所以在做低剂量CT图像后处理或者AI对CT超分辨率重建之前业界默认的做法是先补一版高质量的FBP重建作为参照真值。FBP出的图不好后面加什么网络都是补不完的。这不是说FBP不会被替代而是说它是链路里最稳的一环出了问题可以逐步骤排查。3. 读giftcja扇形CT数据并重建第一张图几何参数、代码和参数档3.1 拿到数据先确认四类几何参数拿到giftcja这类扇形CT投影数据不要急着写重建代码先花半小时把数据里带的几何信息翻出来。通常是三种来源数据集自带的JSON/YAML参数文件、扫描日志里的文本、或者HDF5文件里的属性。你需要确认四个参数源到旋转中心距离D_so、源到探测器距离D_sd、探测器单元间距det_pitch、投影角度序列起始角和角度步长。此外还要确认探测器是等距排列还是等角排列这决定了γ角的计算方式。如果数据集没给参数文件只能自己测量那就要用标定模体。常见做法是扫描一根已知直径的细钢针或者一个球体重建后看图像里针的位置和形状反推几何参数。这一步很磨人但必须做因为后面所有代码都建立在几何参数之上。很多公开CT数据集会附带几何标定文件类似TCIA下载的数据包里通常带扫描参数giftcja这类带编号的数据命名风格也往往能在配套说明里找到几何信息但格式不一定规范需要自己解析。3.2 最小可跑的扇形FBP重建代码重排平行FBP下面这段代码是从HDF5读投影数据、做重排、滤波、反投影的最小闭环。为了少依赖只用了NumPy和SciPy的rotate。import numpy as np import h5py from scipy.ndimage import rotate import json # ---------- 1. 读数据与几何参数 ---------- with h5py.File(giftcja_scan.h5, r) as f: sino f[projection][()] # shape (n_angle, n_det) meta json.loads(f.attrs[meta]) d_so meta[d_so] # 源到旋转中心单位mm d_sd meta[d_sd] # 源到探测器单位mm det_pitch meta[det_pitch] # 探测器单元间距mm beta_deg np.arange(sino.shape[0]) * meta[angle_step] # 每个投影角 n_angle, n_det sino.shape # ---------- 2. 扇形重排到平行束 ---------- # 输出网格t 是射线到旋转中心距离theta 是平行束角度 n_t n_det t_max d_so * np.sin(np.arctan((n_det - 1) / 2 * det_pitch / d_sd)) t_grid np.linspace(-t_max, t_max, n_t) theta_deg np.linspace(0, 180, n_angle, endpointFalse) sino_para np.zeros((n_t, n_angle)) for i in range(n_angle): beta np.deg2rad(beta_deg[i]) for j in range(n_det): u (j - (n_det - 1) / 2) * det_pitch gamma np.arctan2(u, d_sd) t d_so * np.sin(gamma) theta np.rad2deg(beta gamma) % 180 k np.argmin(np.abs(theta_deg - theta)) # t 方向用线性插值 sino_para[:, k] np.interp(t_grid, [t - 1e-6, t 1e-6], [sino[i, j], sino[i, j]])这段代码的逻辑是把每条扇形射线映射到平行束网格上。u是当前探测单元相对中心射线的位置gamma是扇形夹角t是这条射线到旋转中心的距离theta是平行束投影角。因为旋转一周采集的数据足够密角度方向用最近邻取点t方向做线性插值。注意循环里是累加因为多条投影的射线可能落在同一个输出角附近。接下来是做滤波和反投影# ---------- 3. 斜坡滤波 ---------- def ramp_filter(sino_para): n_t, n_theta sino_para.shape # 频域乘 |f|即斜坡核 freqs np.fft.rfftfreq(n_t, d1.0) ramp np.abs(freqs) * (2 * n_t) # 幅度归一化经验系数 filt np.fft.irfft(np.fft.rfft(sino_para, axis0) * ramp[:, None], nn_t, axis0) return filt sino_filt ramp_filter(sino_para) # ---------- 4. 反投影 ---------- N n_t recon np.zeros((N, N)) for k, th in enumerate(theta_deg): # 把一维投影延拓成竖直条纹图 proj_img np.tile(sino_filt[:, k], (N, 1)) # 旋转到该投影角度 recon rotate(proj_img, angle-th, reshapeFalse, order1) recon * np.pi / (2 * len(theta_deg))反投影的原理是把滤波后的投影值沿射线方向均匀铺回图像。np.tile把一维投影铺成列方向一致、行方向重复的条纹图rotate再把条纹旋转到对应的投影角度。order1是线性插值比最近邻平滑也不会像高阶样条那样过冲。最后那个归一化系数是经验值不同斜坡核幅度略有差异第一次跑完和模拟体模对比再微调即可。实际跑之前强烈建议先不处理真实数据先拿一个已知图像做正向投影生成扇形投影再跑重建确认全流程没写反。否则输出全黑或全白时几何符号错误和归一化错误根本分不清。3.3 参数怎么调探测器抽稀、角度方向、滤波核第一个要调的是探测器抽稀。如果n_det是2048甚至4096而图像只想重建到512x512不必把重排网格也开成4096。一般来说重排网格的t方向取最终图像边长即可也就是512。这样做能明显加快反投影速度也不会损失视觉分辨率。第二个要调的是角度方向符号。扇形CT围绕物体旋转如果顺时针和逆时针的定义与数据采集方向相反重排出来θ方向可能差一个负号重建图像会左右镜像。判断方法很简单重建一个椭球体模看长轴方向是否和原始放置方向一致。不一致就把theta beta gamma改成beta - gamma这一条值得写在代码注释里。第三个是滤波核选择。低噪声数据用Ram-Lak最锐利低剂量CT图像噪声大时用Hamming窗更稳。具体来说频率域乘的不是|f|而是|f| * (0.54 0.46*cos(π*f/f_max))相当于把高频部分的噪声压下去。代价是图像边缘略糊但整体看起来干净很多。后面验证章节会教你怎么定量选。4. 避坑扇形CT重建最容易翻车的四处细节4.1 旋转中心偏移现象是重影、星芒和边缘双影现象重建出的切片边缘有一圈虚影点状目标周围出现星芒状条纹左右两侧的清晰度明显不对称。很多人的第一反应是滤波核写错了其实不是。原因旋转中心没有落在探测器投影坐标的零点上。实际CT系统里旋转轴和投影坐标系原点通常有零点几个像素的偏差。重排代码里(n_det - 1) / 2假设探测器中心正对旋转中心但真实数据不满足。解决先做一个粗略标定。取0度和180度两组投影理想情况下这两组投影互为镜像。把其中一组翻转后与另一组做互相关峰值偏移量就是旋转中心偏差。更简单的做法是重建一组针孔模体改变中心偏移参数看图像最清晰时的偏移量。在重排代码里u的计算改成u (j - (n_det - 1) / 2 - offset) * det_pitch多试几个offset值一般能修好。4.2 坏道与坏探测单元一条条亮暗条纹的由来现象重建图像上出现横贯整个视野的亮条纹或暗条纹有时候是几条平行的条纹一起出现像梳子齿一样。原因探测器某些通道响应异常读数偏低或偏高。这些“坏道”在滤波阶段会被斜坡核放大成一条条的伪影。工业CT里这种情况尤其常见因为探测器老化或受辐射损伤。解决在投影域做坏道修正。先统计各探测器通道在全部角度下的均值那些均值明显偏离整体水平的通道就是坏道。修正可以用相邻通道线性插值也可以用双向平均。插值时要注意如果坏道在边缘相邻通道可能不在有效视野内这时要用内侧最近有效通道外推。更稳妥的办法是采集一次暗场和亮场先做增益校正再做坏道替换。这一步必须在滤波之前做因为滤波会把单个通道的错误扩散到整条射线。4.3 截断伪影视野不足时出现的杯状伪影现象物体超出重建视野图像边缘出现明显的杯状凹陷灰度从中心到边缘整体下降严重时还会有环状亮边。原因扫描时物体比最大成像视野大部分射线没有穿过物体就打到探测器外投影数据不完整。扇形FBP要求每个角度下物体完整覆盖射线束截断数据在重排后等于是sinogram外部补了零斜坡滤波会把截断边界上的阶跃信号变成一条条正弦状伪影。解决先确认最大视野。重排里的t_max就是最大成像视野半径如果物体半径接近或者超过这个值不应该强行重建。工程上常用的方案是裁剪重建区域只重建物体中心满足完整投影的部分边缘部分直接舍弃。如果必须看边缘就要做扩展视野重建常见做法是把截断部分做外推用正弦函数或者多项式拟合格把sinogram平滑延伸到探测器边缘再走FBP流程。这个外推参数很敏感一般要按数据噪声水平调整。4.4 低剂量CT的噪声与滤波核放大R-L核与Hamming核怎么选现象同一套投影数据Ram-Lak滤波核重建出来的图像看起来“脏”布满颗粒状噪声换Hamming核后干净了但边缘也糊了。原因低剂量CT图像投影数据是典型的泊松噪声主导噪声在频域里同样落在高频段。Ram-Lak核在高频段增益最大正好把噪声放大了。这不是重建代码的问题是滤波器选择与数据噪声不匹配的问题。解决把滤波核改成带窗函数的形式。实践里先跑一版Ram-Lak如果噪声明显再换成Hamming然后对比两版的边缘保留情况。更细的方案是调节窗函数的截止频率比如Cosine窗的可调参数比Hamming少但过渡更缓。工业CT里如果扫描的是金属工件密度对比本来就很强保留Ram-Lak更合适医疗低剂量CT图像这类场景Hamming窗是默认起点。别指望一个核走天下滤波核选择本质上是分辨率和噪声的权衡。5. 用模拟体模验证重建参数RMSE曲线和自动化标定5.1 为什么先用模拟数据再碰真实数据重建参数全凭感觉调很容易陷入“看着还行但不知道对不对”的状态。真实投影数据没有标准答案你无法判断图像里的某个细节是真实结构还是伪影。所以在跑giftcja数据之前要先生成一组模拟扇形投影用已知的图像做正向投影再用自己的重建链路把它恢复出来。如果恢复结果和已知图像接近说明重建逻辑是对的如果不对就能精准定位是几何映射、滤波核还是归一化的问题。这一点特别重要真实数据的重建经常因为旋转中心偏移、探测器响应不均匀等因素引入伪影而这些伪影和算法错误混在一起时新手几乎无法区分。模拟数据没有这些设备误差是纯净的算法自检环境。5.2 用Shepp-Logan体模生成扇形投影做回归验证Shepp-Logan体模是CT重建里最常用的标准测试模体由一组椭圆组成每个椭圆有确定的中心、半轴、旋转角和密度。生成扇形投影不需要复杂的正向投影计算直接对每条射线求它与所有椭圆的交点弦长乘密度累加即可。import numpy as np def ray_ellipse_hit(p1, p2, ell): # p1, p2: 射线上两点ell (cx, cy, a, b, phi, rho) cx, cy, a, b, phi, rho ell dx, dy p2[0] - p1[0], p2[1] - p1[1] x0, y0 p1[0] - cx, p1[1] - cy # 把射线变换到椭圆局部坐标系 c, s np.cos(phi), np.sin(phi) xt1 c * x0 s * y0 xt2 c * dx s * dy yt1 -s * x0 c * y0 yt2 -s * dx c * dy A (xt2 / a) ** 2 (yt2 / b) ** 2 B 2 * (xt1 * xt2 / a ** 2 yt1 * yt2 / b ** 2) C (xt1 / a) ** 2 (yt1 / b) ** 2 - 1 disc B * B - 4 * A * C if disc 0: return 0.0 sq np.sqrt(disc) t1 (-B - sq) / (2 * A) t2 (-B sq) / (2 * A) if t2 0 or t1 1: return 0.0 t1 max(t1, 0.0) t2 min(t2, 1.0) return rho * (t2 - t1) * np.sqrt(dx * dx dy * dy)射线用两点表示p1是源点p2是探测器单元点。这个函数返回射线在该椭圆内走过的弦长乘密度。对每个椭圆累加就得到这条射线的投影值。代码里的关键是把射线变换到椭圆坐标系再求交避免解非对称的通用圆锥曲线方程。有了射线交点函数生成投影就很直接先定义一组Shepp-Logan椭球参数然后对每个投影角度、每个探测单元计算射线源点与探测器点调用这个函数累加得到模拟sinogram。生成的模拟投影再喂给前面的重排和FBP代码重建出图像后与原始体模的像素值矩阵对比。一版代码跑通后几何参数和归一化系数就不会再错了。5.3 三个关键指标与一个自动标定旋转中心的做法重建结果对比不能只看肉眼的“像不像”要有数值指标。最常用的是RMSE和SSIMfrom skimage.metrics import structural_similarity as ssim # 重建图和真值可能整体有灰度缩放先做线性回归 a, b np.polyfit(recon.ravel(), gt.ravel(), 1) recon_cal recon * a b rmse np.sqrt(np.mean((recon_cal - gt) ** 2)) ssim_score ssim(recon_cal, gt, data_rangegt.max() - gt.min())整体灰度缩放是一个很常见又容易被忽略的坑。扇形重排后因为插值密度和归一化系数的原因重建出的绝对值往往和真实线性衰减系数差一个比例常数。直接算RMSE会把整体亮度差也当成误差掩盖局部的质量问题所以先做一阶线性回归修正再算RMSE和SSIM。SSIM超过0.9说明结构基本还原RMSE要看图像灰度范围通常归一化到0到1后小于0.05比较好。如果SSIM低于0.8多数情况下问题出在旋转中心偏移而不是滤波核。旋转中心的自动标定也可以利用RMSE把旋转中心偏移量设成一个待搜索参数在-2到2像素范围里以0.1像素步长扫描每个偏移量重建一次计算与真值的RMSERMSE最低处对应的偏移量就是估计值。这个过程看着笨但在模拟数据上几秒就能跑完真实数据也适用前提是能找到一个大致形状参考。这个技巧值得放进自己的重建工具包里以后换一台扫描设备就能用标准模体把旋转中心重新标一遍。6. 进阶从二维扇形FBP走向三维与加速二维扇形FBP跑通以后真正的工程挑战是从切片走向体积。工业CT、医疗CT的数据量都在千张投影级别的体数据纯NumPy的循环反投影速度不够用。常见的做法是把重排和FBP交给GPUASTRA Toolbox的FBP_CUDA接口直接吃扇形投影参数一条命令就能完成带几何标定的重建比自己写的循环快两个数量级。迁移时需要注意ASTRA的几何定义里旋转中心和探测器轴的零点约定与手写代码不一定一致务必用模拟数据重新验证一遍符号方向。另一个值得投入的方向是低剂量CT图像和AI后处理的衔接。FBP重建得到的图像是后续超分辨率重建、去噪网络的基础输入但网络训练前要先确认FBP的滤波核和窗函数固定下来不要在训练中途更换。我自己经历过一次换了Hamming窗之后之前调的AI模型全要重训因为高频细节分布变了。最后分享一个习惯每次改几何参数、换滤波核或迁移环境后都先跑一遍Shepp-Logan体模回归确认SSIM没有跳变再处理真实数据。这个动作花不到一分钟但能避免一整天的无效重建。CT重建这行稳定比炫技更重要。希望帮到你。本文还有配套的精品资源点击获取