简介一份面向CT成像学习与研究者的扇形束CT重建资源包聚焦滤波反投影FBP算法在扇形束条件下的实现与理解。包内包含MATLAB脚本F1.m用于演示FBP重建流程可调整投影角度、源到物体距离、探测器大小、重建矩阵及重建尺寸等关键参数同时附有一篇关于层析成像实现的PDF文献可辅助理解从扇形束CT扩展到计算层析成像的理论背景与应用场景。资源共2个文件含1个m脚本和1个pdf文档压缩包总大小345KB轻量易用适合具备一定数值计算与图像处理基础的读者快速上手。已有247人浏览学习。通过运行脚本和对照文献读者能直观掌握扇形束投影数据获取、滤波反投影重建的核心步骤并迁移到层析成像等扩展应用中对算法细节与工程实现均具参考价值。1. 扇形FBP重建不是傻跑一个函数拿一组扇形束CT投影数据去重建横断面大多数工程师的第一反应是找现成库调iradon或者fbp接口。但真到了工业CT场景事情会变得很不舒服射线源到旋转中心的距离、探测器像元间距、旋转步进角度这些参数里只要有一个不对重建结果就是一片径向伪影叠加模糊。更麻烦的是很多工业CT导出的数据连头文件都残缺只有一串整数看起来像“ct后缀文件是啥”的答案却又不是标准DICOM。自己写一遍扇形FBP不是为了重新发明轮子而是为了在拿到任何一组数据时能直接控制从扇束几何到滤波核再到反投影的每一个环节。这篇文章讲的是一条可复现的路径把扇形束投影数据组织成正弦图推导扇形束下的FBP公式然后写一个最小Python实现用Shepp-Logan体模验证最后落在几个实际参数调试技巧上。适合至少写过平行束反投影、但没摸过扇束几何的算法工程师也适合需要把现成CT重建代码改造成自定义几何的工业检测从业者。标题里的giftcja可以理解为一组可调滤波核参数的代称后面会给出它在频域里的具体写法。2. 扇形束投影与FBP重建平行束公式如何迁移到扇形束2.1 扇束几何里的三个坐标关系平行束CT里投影数据是一个二维矩阵横轴是射线到中心的距离u纵轴是旋转角度theta。扇形束CT则不同每个角度下射线从同一个点源发出穿过物体后落在探测器上探测器上的位置用s表示。扇形束投影矩阵的行是旋转角度alpha列是探测器单元索引s也就是等角采样或等距采样。拿常见的等距扇形束来说点源到旋转中心的距离是R探测器到旋转中心的距离是D通常探测器在旋转中心另一侧但符号约定不同。某条射线的几何位置由扇形角gamma atan(s/D)决定。离轴角beta等于当前旋转角alpha加上扇形角gamma。这些几何关系直接决定了反投影时每个投影值应该放回图像平面的哪个位置。在平行束里投影数据直接放进一个(n_theta, n_s)矩阵就能做滤波和反投影但在扇束里同样的矩阵还要额外乘一个与R和s有关的加权因子。这个因子不是公式里凭空出现的它来自扇束射线的发散性本质上是换元时产生的雅可比项。2.2 FBP在扇形束下的加权、滤波、反投影三步顺序扇形FBP和平行束FBP最大的区别不只是多了加权而是处理顺序不能乱。平行束的流程是先对每个角度的投影做一维滤波再沿原路径把值平铺回图像空间。扇束下的流程变成了三步对投影数据按探测器位置做余弦加权权值是R / sqrt(R^2 s^2)对应到中心射线夹角的关系。这一步需要乘在滤波之前。对加权后的每一行投影做一维频域滤波滤波核是斜面高通的|omega|形式的截断窗也就是后面要写的giftcja滤波核。把滤波结果沿真实扇束射线路径反投影并且每条射线还要再乘1 / L^2的权重其中L是从源到重建像素的距离这是扇束反投影独有的衰减补偿。如果忽略第三步的1/L^2图像中心会偏亮边缘会偏暗。如果第一步的加权放在滤波之后虽然数值上差别不大但对于噪声放大的抑制效果会变弱所以最好保持这个顺序。下面这个参数关系表可以帮你快速对齐自己的数据参数符号典型工业CT取值影响源到旋转中心距离R200–600 mm决定扇束张开角影响边缘伪影探测器到源的距离RD300–800 mm决定几何放大比探测器单元数N_s512–2048决定采样间隔影响空间分辨率探测器像元间距ds0.1–0.5 mm与R一起决定扇形角范围旋转步进角度d_alpha0.1°–1°角度采样越多反投影越平滑3. 用Python实现扇形FBP最小可运行代码3.1 模拟一组扇形束投影数据作为输入在没有真实工业CT数据时最快验证重建算法的方式是使用解析体模。常见的做法是先构造Shepp-Logan体模然后在扇束几何下做正向投影。这里不直接调库而是对每一条射线做数值积分这样你就能完全掌控几何参数。下面的代码生成一个简单散射点体模的投影用scipy.ndimage.rotate来模拟旋转物体投影算子对每个角度计算物体图像的加和。更精确的写法是用射线穿体素的数值积分但对于验证FBP流程可以用雷达式投影加旋转的方式来近似。import numpy as np from scipy.ndimage import rotate def simulate_fan_projection(image, angles, R, D, detector_pitch, n_detector): 用旋转图像 逐射线求和来模拟扇形投影 image: 二维图像数组 angles: 旋转角度列表弧度 R: 源到旋转中心距离 D: 旋转中心到探测器距离 detector_pitch: 探测器像元间距 n_detector: 探测器单元数量 返回的形状: (len(angles), n_detector) n_angles len(angles) sino np.zeros((n_angles, n_detector)) center np.array(image.shape) / 2.0 detector_positions (np.arange(n_detector) - n_detector / 2.0) * detector_pitch for i, alpha in enumerate(angles): # 将物体旋转 -alpha使得射线方向固定为竖直方向 rotated rotate(image, anglenp.degrees(alpha), reshapeFalse, order1) # 对每一列做“源点-探测器”连线上的加权求和 for s_idx, s in enumerate(detector_positions): # 简化模拟把旋转后的图像沿竖直方向积分忽略几何发散 # 真实场景需要按扇形路径逐点求交点这里用行和代替 column_sum rotated[:, s_idx n_detector//2].sum() sino[i, s_idx] column_sum return sino这段代码的注释已经说清它是简化版真正的扇束正投影需要沿源到每个探测器单元的射线对图像做线积分而不是简单的列求和。参数里angles是旋转角度R和D在这个简化版本里没有直接用到但保留这两个参数是为了让函数签名和真实算法一致后面重建时就要用到。3.2 编写频域滤波核giftcja滤波核现在写滤波核。FBP里最经典的滤波核是Ramp滤波器也就是频域|omega|但直接截断会产生振铃。工业CT里常用某种窗函数去平滑giftcja可以理解为一组参数化的滤波核写法我把斜坡滤波和高斯窗叠在一起得到一个在截止频率附近柔和衰减的核。def giftcja_filter(n_detector, detector_pitch, cutoff1.0, beta2.0): giftcja滤波核频域斜坡 高斯窗 n_detector: 探测器单元数量 detector_pitch: 像元间距决定奈奎斯特频率 cutoff: 归一化截止频率0~11表示达到奈奎斯特频率 beta: 高斯窗宽度参数越大越平滑 返回: 时域滤波核长度为 n_detector freq np.fft.fftfreq(n_detector, ddetector_pitch) # 斜坡滤波 |f| ramp np.abs(freq) # 高斯窗在奈奎斯特频率附近衰减 nyquist 1.0 / (2.0 * detector_pitch) normalized_freq freq / nyquist window np.exp(-0.5 * (normalized_freq / (cutoff * beta))**2) filter_kernel ramp * window return filter_kernel def filter_projection(p, filter_kernel): 对单角度投影做频域滤波 p: 一维投影数组 filter_kernel: 频率域的滤波核长度与p一致 p_fft np.fft.fft(p) p_filtered np.fft.ifft(p_fft * filter_kernel).real return p_filtered逻辑说明giftcja_filter返回的滤波核是频域数组需要和np.fft.fft(p)的结果逐元素相乘再做逆FFT。cutoff控制频带宽度beta控制窗的形状。beta越大高频衰减越多重建图像越平滑但空间分辨率会下降。参数cutoff若取0.8等于把高于奈奎斯特频率80%的部分全部压低常用于噪声偏大的工业CT数据。实际使用时为了消除边界伪影滤波前要对投影做周期延拓或镜像填充。最简单的做法是把投影数组翻折成原来的三倍滤波后取中间段。这一点容易被忽略导致重建结果在图像边缘出现条纹。下面是滤波加反投影的完整主体def fan_fbp_reconstruct(sino, angles, R, D, detector_pitch, image_size, filter_kernel): sino: 形状为 (n_angles, n_detector) 的投影矩阵 angles: 每个投影角度弧度 R: 源到旋转中心距离 D: 旋转中心到探测器距离 image_size: 重建图像边长单位像素 n_angles, n_detector sino.shape img np.zeros((image_size, image_size)) # 第一步余弦加权 s_indices np.arange(n_detector) - n_detector / 2.0 fan_angle np.arctan(s_indices * detector_pitch / D) cos_weight np.cos(fan_angle) weighted_sino sino * cos_weight # 第二步逐角度滤波 filtered_sino np.zeros_like(weighted_sino) for i in range(n_angles): filtered_sino[i] filter_projection(weighted_sino[i], filter_kernel) # 第三步反投影 center image_size / 2.0 for i, alpha in enumerate(angles): for x in range(image_size): for y in range(image_size): px x - center py y - center # 旋转到探测器坐标系 x_rot px * np.cos(alpha) py * np.sin(alpha) y_rot -px * np.sin(alpha) py * np.cos(alpha) # 计算源到像素的射线在探测器上的投影位置 denominator R D - x_rot if abs(denominator) 1e-9: continue s_coord (y_rot * D) / denominator s_index s_coord / detector_pitch n_detector / 2.0 if 0 s_index n_detector: idx0 int(s_index) idx1 min(idx0 1, n_detector - 1) weight (s_index - idx0) val (filtered_sino[i, idx0] * (1 - weight) filtered_sino[i, idx1] * weight) L np.sqrt((denominator)**2 (y_rot)**2) img[x, y] val / L / L # 注意这里未乘角度步长loop结束后可统一乘d_alpha return img反投影里的denominator就是源到像素的沿光轴方向距离y_rot * D / denominator是射线在探测器平面上的落点坐标也就是投影位置。L / L的平方来自扇束反投影权重。最后在主循环外乘以d_alpha角度步长。这段代码是三重循环速度较慢实际工程中要用scipy.ndimage.map_coordinates或者先在探测器方向预积分再把反投影向量化。4. 扇形FBP的三个必调参数与ct数据文件读取策略4.1 源到旋转中心距离R的误差会怎样体现R是扇形FBP里最敏感的参数。工业CT的机械装配误差通常会导致实际R和标称值有几毫米偏差反映在重建图上就是离中心越远环形伪影和双影越明显。调试方法是用一根细金属丝放在视野边缘投影后如果重建点不是圆形而是拉长的弧就说明R偏大或偏小。调整时可以采用二分法把R从标称值加减5毫米范围各跑一次重建看边缘点扩散函数的变化。反投影中的denominator R D - x_rot直接参与除法当R变小时分母变小边缘像素的权重被放大靠近源一侧的图像会被拉伸。所以R的精确标定最好用专用的线对卡或球形标定体模完成不要依赖软件默认值。4.2 探测器像元间距与旋转角步长如何配合探测器像元间距ds和旋转角步长d_alpha共同决定了正弦图的空间采样是否均匀。若想重建图像的空间分辨率大致等于探测器分辨率角度步长应满足d_alpha ds / (2*pi*R)量级的条件。比如R300mm, ds0.2mm时角度步长要小于0.02度这在一圈360度下意味着约18000次投影很多工业CT实际是做不到的。实际遇到投影角度不足时常见的折中是先对正弦图做角向插值将角度数翻倍。插值方法用线性插值即可但更好的选择是采用带限插值因为正弦图在角向上是带限信号。插值前要先把投影数据做归一化去除射线源强度波动否则插值会引入频谱泄漏。下面的代码演示了用scipy.interpolate做角向重采样from scipy.interpolate import RegularGridInterpolator def resample_angles(sino, original_angles, target_angles): 对正弦图进行角向插值 sino: (n_angles, n_detector) original_angles: 原角度数组 target_angles: 目标角度数组 n_angles, n_detector sino.shape detector_idx np.arange(n_detector) interp RegularGridInterpolator( (original_angles, detector_idx), sino, methodlinear, bounds_errorFalse, fill_valueNone ) grid_angles, grid_det np.meshgrid(target_angles, detector_idx, indexingij) return interp((grid_angles, grid_det))这里RegularGridInterpolator期望输入网格是单调递增的而原始角度通常是从0到360度递增可以直接使用。target_angles可以用np.linspace(0, 2*np.pi, N_target, endpointFalse)生成。注意endpointFalse避免首尾角度重复否则插值边界会产生重影。4.3 读取CT数据ct后缀文件是啥当同事发来一堆没有说明的CT数据文件时先看扩展名。常见的.ct文件往往是多个文件组成的序列比如scan_0001.ct、scan_0002.ct每个文件对应一个旋转角度的投影。有些是从软件里导出的原始投影有些是重建后的切片判断方式是看单个文件大小与探测器单元数乘角度的关系。若文件字节数正好等于n_angles * n_detector * 4就是32位浮点投影。读取代码需要先确认字节序。工业CT数据常存为小端无符号短整型用以下方式读取import os, struct def read_ct_projection(filepath, n_detector, dtypeu2): 读取单个CT投影文件 filepath: 文件路径 n_detector: 探测器单元数用于判断单行长度 dtype: 字节序数据类型 返回一维数组 filesize os.path.getsize(filepath) n_values filesize // np.dtype(dtype).itemsize if n_values % n_detector ! 0: # 常见做法尝试按总数值除以n_detector得到角度数 n_angles n_values // n_detector if n_values ! n_angles * n_detector: raise ValueError(文件大小与探测器数不匹配) arr np.fromfile(filepath, dtypedtype) return arr.reshape(-1, n_detector).astype(np.float32)这里u2只适用于小端16位无符号如果数据是32位浮点改成f4。读取后要检查统计量如果投影值的最大值接近65535说明原始是16位并可能被截断如果值域在0到1之间说明已做过标定校准。遇到原始暗电流不均的先做空气归一化即每个探测器单元除以空扫投影的对应值再乘平均强度。ct后缀文件是啥这个问题的标准答案就是先按这个模板去猜而不是直接套用读取程序。5. 用Shepp-Logan体模验证扇形FBP重建的细节5.1 从体模到正弦图的最小验证命令前面已经构造了正投影函数现在用它生成一组投影再走完整重建与原始体模对比。取图像大小为256探测器数512R300D300角度从0到360度均匀取180个投影。这组参数虽然角度偏少但足以暴露FBP的几何伪影。angles np.linspace(0, 2*np.pi, 180, endpointFalse) sino simulate_fan_projection(phantom, angles, R300, D300, detector_pitch0.3, n_detector512) filter_kernel giftcja_filter(512, detector_pitch0.3, cutoff0.85, beta2.0) recon fan_fbp_reconstruct(sino, angles, R300, D300, detector_pitch0.3, image_size256, filter_kernelfilter_kernel) recon recon / recon.max() # 归一化方便显示和对比注意这里phantom是前面生成的Shepp-Logan灰度图像可以自行用skimage.data.shepp_logan_phantom()生成。重建后与原始体模做差计算均方根误差同时用skimage.metrics.structural_similarity评估结构相似度。如果重建出的图像边缘出现明显的单侧暗带说明余弦加权被重复执行或者反投影权重用错了符号。验证时还要检查图像中心值是否明显高于四周。中心过亮大概率是反投影时没有使用1/L^2权重而只用了1/L。这种误差通过肉眼不易察觉但画出沿水平中心线的剖面图可以看到中心处有明显的凸起。正确结果应该是中心值在相邻区域的平滑范围内峰值略低于原体模中的高密度椭圆中心。5.2 伪影检查与giftcja参数调整技巧重建完成后重点看图中有没有三种典型伪影第一种是条纹状伪影表现为从高密度物体边缘散射出去的亮线这是投影角度取样不足导致的可以通过增加角度数或角向插值缓解。第二种是圆环形伪影通常不是FBP本身的问题而是探测器响应不均匀或某个探测器单元损坏需要做空气校正。第三种是图像整体模糊原因是滤波核的截止频率太低调高cutoff到接近1并把beta从2.0降到1.0边缘会立刻变得更锐利。giftcja滤波核还有一个调参方向工业CT中经常遇到噪声和伪影并存的情况此时可以把高斯窗的中心频率往高频移动同时加一个很小的正则项。常见做法是在滤波核的频域表达式中增加epsilon项避免分母为零但这个参数不要超过奈奎斯特频率的千分之一否则低对比度结构会被抹平。一个可用的调试脚本是让cutoff在0.6到0.95之间扫五个值每次重建后计算重建图中已知空气区域的标准差标准差最小时的cutoff就是当前噪声水平下的最佳值。最后验证的另一个细节是反投影的旋转方向。很多真实CT的旋转角方向与标准数学坐标相反导致图像左右颠倒或镜像。把重建结果和原始体模做一次对比如果发现水平翻转就在fan_fbp_reconstruct里把alpha取反再跑一次。付出这点代价换来对整个FBP流程的完全掌控在工业CT数据格式混乱时值回票价。本文还有配套的精品资源点击获取