ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

用Python从零实现平行束FBP重建:RL滤波器与反投影详解

用Python从零实现平行束FBP重建:RL滤波器与反投影详解 简介平行束FBP重建Python代码RL滤波器是一份面向初学者的CT图像重建算法实现适合正在学习滤波反投影原理、希望动手验证理论的同学。压缩包共2个文件包含一个.py源码脚本和一个.txt说明文档整体仅2KB代码精简、结构清晰便于逐行阅读与修改。已有181人学习/下载。代码基于NumPy和SciPy完成了投影数据读取、RL滤波、角度变换与反向投影等关键步骤RL滤波器部分重点展示了频率域加权如何影响噪声抑制与边缘保留txt文档中可能包含参数设置、运行说明或实验结果方便对照调试。通过实际运行并调整滤波器参数可以直观观察重建图像质量变化从而加深对FBP算法和医学成像基础流程的理解为后续学习CT重建或医学影像分析打下基础。1. 平行束FBP重建的核心RL滤波器决定了你看到的是不是“真”CT图像FBP滤波反投影重建在平行束几何下有一个反直觉的结论不滤波直接反投影得到的是边缘模糊的“影子”而不是断层图像真正让图像清晰起来的是频域里那个斜坡形状的RL滤波器Ram-Lak也就是傅里叶域中的 |ω| 响应。标题里的这套Python代码之所以建议刚学FBP的同学反复对照不是因为它实现了什么新颖算法而是它把“投影模拟—滤波—反投影—参数调整”完整串在几十行代码里跑通一次就能直观看到每一个环节的作用。本文不依赖任何重建库用 numpy 和 scipy 从零搭一套可运行的平行束FBP流程并把滤波器的离散化、缩放因子、角度和探测器参数这些容易踩坑的地方逐个拆开讲清楚。2. 中心切片定理、FBP公式与RL滤波器的离散化2.1 平行束投影与正弦图先明确数据从哪来平行束FBP的输入不是一张图像而是一组“从不同方向照射得到的线积分”也就是 Radon 变换。对二维函数 f(x,y)角度 θ 下的平行束投影定义为p(s, θ) ∫ f(x, y) δ(x cos θ y sin θ − s) dx dy把不同角度的投影按列堆叠形成一张二维数组 sinogram横轴是探测器位置 s纵轴是投影角度 θ。FBP 要做的事情就是从这个 sinogram 把 f(x,y) 重建出来。在纯 Python 环境里模拟平行束投影不需要正儿八经写射线求交。常见做法是把图像旋转 θ 角后沿垂直方向求和这等价于计算一组间距为1、方向垂直于探测器阵列的平行射线积分。要注意旋转后要关闭 reshape让输出图像保持原尺寸边界补零这样探测器覆盖范围才是完整。import numpy as np from scipy import ndimage def parallel_beam_project(img, angles): 平行束前向投影旋转图像后按行求和。 img: 2D 灰度图方形angles: 单位度的旋转角数组。 proj [] for theta in angles: rot ndimage.rotate(img, theta, reshapeFalse, order1) proj.append(rot.sum(axis0)) return np.array(proj)代码逻辑不复杂对每个角度把图像旋转到该视角下再沿垂直探测器的方向求和。order1是线性插值能避免最近邻插值导致的投影边缘台阶。reshapeFalse保证旋转后画布尺寸不变超出部分补零恰好对应平行束扫描中探测器覆盖区域大于物体截面的场景。初学容易在这里把尺寸搞丢丢掉之后反投影的坐标系就对不上了。2.2 中心切片定理与 FBP 推导RL滤波器是推导出来的不是拍脑袋加的中心切片定理说图像在某方向上的平行投影做一维傅里叶变换后正好是图像二维频谱在相同方向过原点的那条线。如果采集了 0 到 π 的所有角度就能铺满整张二维频域平面理论上可以直接做二维逆傅里叶变换得到图像。但直接做有两个问题极坐标系下的频率样本不均匀——中心密、边缘疏直接插值重建会带来低频偏重、细节模糊另一个问题是插值精度差。于是把傅里叶域积分写成极坐标经过变量代换会自然出现一个 |ω| 因子这就是 RL 滤波器的由来。公式落到工程实现就三步对每个角度的投影做一维 FFT乘以斜坡 |ω|再做逆 FFT最后把所有角度下的滤波投影“涂抹回”图像平面并累加。空域离散实现时RL 滤波器有明确的收敛形式h(s) 在零点的取值由采样决定非零偶数点值为 0奇数点值为 h(k) −1/(π²k²)。用频域乘法则只需一行def ramp_filter_fft(n_det): 构造长度为 n_det 的斜坡滤波器频响。 freq np.fft.fftfreq(n_det) return np.abs(freq)np.fft.fftfreq返回的是从 -0.5 到 0.5 的归一化频率取绝对值之后得到一个 V 字形的斜坡响应。DC 分量频率0乘的是 0也就是说投影的整体基线在滤波中会被完全压掉这也解释了为什么 RL 重建出的图像没有直流偏移。2.3 滤波器不唯一为什么专门拎出 RLFBP 的滤波环节里RL 是唯一不施加任何窗函数的斜坡滤波器它的频率响应在全频段保持线性放大。好处是空间分辨率最高坏处是高频段噪声也等比例被放大带噪投影重建出的图像颗粒感极重。滤波器频域形式噪声特性学习价值RLRam-Lak|ω|完整斜坡噪声随频率线性放大能看到 FBP 最原始的数学形态Shepp-Logan|ω|·sinc(ω/2)高频压制中等伪影与分辨率的折中Hamming|ω|·[0.54 0.46·cos(πω)]高频明显压弱工程中常见的选择标题里特意注明 RL 滤波器说明代码的定位是教学——直接把最核心的斜坡响应当作“滤波器”呈现不掺入窗函数的额外光滑效果让学习者先理解 |ω| 这个因子从哪里来再去折腾改良版本。理解了 RL换成汉明窗只是一行加权乘法。3. 用Python从零搭出平行束FBP重建完整代码3.1 最小可用代码投影、滤波、反投影三块拼一起反投影的朴素实现是“把滤波后的投影退回来”把一维投影复制成二维矩阵沿探测器方向铺开再旋转回投影时对应的角度叠加到重建图上。import numpy as np from scipy import ndimage def fbp_parallel_beam(sinogram, angles, img_sizeNone): 平行束 FBP 重建 main 流程。 sinogram: shape (n_angles, n_det) 的投影数据 angles: 投影角度数组度 img_size: 重建图像边长默认等于探测器数量 if img_size is None: img_size sinogram.shape[1] # 1. RL 滤波逐角度 FFT - 乘斜坡 - IFFT freq np.fft.fftfreq(sinogram.shape[1]) ramp np.abs(freq) filtered np.zeros_like(sinogram) for i in range(sinogram.shape[0]): spec np.fft.fft(sinogram[i]) filtered[i] np.real(np.fft.ifft(spec * ramp)) # 2. 反投影累加 recon np.zeros((img_size, img_size)) dtheta np.pi / len(angles) # 角度步长单位弧度 for i, theta in enumerate(angles): # 把一维滤波投影变成二维“涂抹状”每一行都是同一根投影 column np.tile(filtered[i], (img_size, 1)) # 旋转回投影时的原方向负号补偿前向投影的正向旋转 rotated ndimage.rotate(column, -theta, reshapeFalse, order1) recon rotated return recon * dtheta逻辑分两大块。滤波部分对 n_det 个探测器造成的数组做一维 FFT乘 |ω| 后反变换得到的是与 sinogram 同形状的滤波正弦图。反投影部分从数学上看是连续积分 ∫₀^π Q_θ(s) dθ 的离散化累加后乘步长 dtheta 就是对这个积分的近似如果不乘结果绝对值会随角度采样数直线增大换一组角度数就要重新调尺度。3.2 用个假体跑通README常说的“结果已测试”到底在测什么凡标明“实验结果已测试”的代码包通常都会附一张 Shepp-Logan 或圆模的重建前后对比图。自己复现时判断跑通的标准很简单重建图像轮廓清晰、背景贴近零、中心区域灰度与原始图像在同一量级。def circle_phantom(size256): 直径约 160 像素的圆形假体灰度 1.0。 y, x np.mgrid[0:size, 0:size] return np.where((x - size/2) ** 2 (y - size/2) ** 2 80**2, 1.0, 0.0) angles np.linspace(0, 180, 180, endpointFalse) phantom circle_phantom(256) sino parallel_beam_project(phantom, angles) recon fbp_parallel_beam(sino, angles) print(f重建图像均值: {recon.mean():.3f}, 原图均值: {phantom.mean():.3f})用 180 个角度、256 个探测器、圆形假体重建均值对比如果差出数量级先检查反投影是否乘了 dtheta再检查射线方向的正负号是否配对。endpointFalse让角度覆盖 [0, 180°) 而不是 [0, 180°]因为 0° 和 180° 的平行射线集合完全重合取两者的重复投影会引入额外权重。4. RL重建与参考库对照角度数、探测器、滤波器长度怎么调4.1 用 scikit-image 的 iradon 当标尺从零写出的重建算法验证方式是对照成熟实现。scikit-image 的iradon默认 filter 就是 ramp正好等于 RL 滤波器适合当基准。对比前要统一角度数组iradon 要求 sinogram 的第一维是角度数第二维是探测器样本数角度单位默认是度和上面代码约定一致。from skimage.transform import iradon recon_ref iradon(sino, thetaangles, filter_nameramp, circleTrue) recon_mine fbp_parallel_beam(sino, angles) error np.sqrt(np.mean((recon_ref - recon_mine) ** 2)) print(f与参考实现的 RMSE: {error:.4f})由于尺度因子、插值方式和边界处理不完全一致两者会有幅度差异这是正常的。判断标准不是 RMSE 无限接近 0而是误差图中没有结构性条纹如果重建图里出现从中心向外辐射的细线多半是角度方向符号错位如果图像整体像“模糊的锅盖”说明滤波丢失或只做了半截。4.2 角度数和探测器数量的直观影响伪影怎么长出来的平行束 FBP 的欠采样伪影有非常固定的形态角度数0–180°探测器数量观测到的问题30256星芒状伪影高对比边缘重影90256有轻微扇叶状纹理边缘尚可180256与参考结果基本一致360256改善不明显重建时间翻倍180128分辨率下降圆形假体边缘发虚180512无明显改善计算量增加角度过少时信息量不足反投影的累加会让每个投影像一把“刀”多把刀片叠在一起就切出星芒探测器太少对应高频信息丢失RL 滤波也只能还原到采样可承载的分辨率。改参数的实验适合配合一个交互式命令行比如用 argparse 把角度数传进去反复跑观察伪影随角度的逐渐消除。4.3 滤波器长度和补零FFT卷积有个隐藏坑滤波长度应当等于探测器数量但实际算法里经常把投影补零到更长的长度再做 FFT防止时域卷积的循环绕回。直接对原长做 FFT 再乘斜坡会出现轻微的“滤波串扰”原因是斜坡核在空域是无限长的直接截断会产生振铃。常见做法是把每条投影先右补零到两倍探测器长度滤波后再截回原长度。def rl_filter_fft_padded(proj_line, pad_factor2): n len(proj_line) n_pad n * pad_factor freq np.fft.fftfreq(n_pad) spec np.fft.fft(proj_line, nn_pad) filtered np.real(np.fft.ifft(spec * np.abs(freq))) return filtered[:n]参数pad_factor2是空间域卷积滤波常用值把循环卷积变成线性卷积避免左边界对右边界的污染。这段代码同样适用于每条投影逐根处理。如果发现重建图左右两侧有亮度差或横向条纹优先检查是不是没做补零。5. 从RL到Hamming噪声场景下滤波器选择的验证技巧RL 滤波器在无噪声或低噪声投影上表现很好一旦投影含有噪声情况立刻翻转斜坡的高频增益连同噪声一起放大重建结果会出现密集的高频颗粒。此时实际工程会用加窗斜坡滤波器常见做法是让斜坡响应乘以一个窗函数Hamming 窗是较温和且容易实现的选择。def hamming_filter_fft(n_det): freq np.fft.fftfreq(n_det) ramp np.abs(freq) window 0.54 0.46 * np.cos(2 * np.pi * freq) return ramp * window验证时构造一个带噪声的投影对同一个 sinogram 分别用 RL 和 Hamming 重建计算噪声注入前后的 RMSE。通常 Hamming 会在分辨率上有一点损失但稳定性和视觉效果明显好得多理想的教学实验是把两幅重建图并排看再叠加点状噪声重复一次就能直观看出 RL 的放大效应。实际排查问题时还有个实用技巧滤波后的 sinogram 本身会暴露异常值。FL 滤波后的数据里如果出现单条亮线往往不是算法问题而是某个角度的投影数据有离群像素把filtered数组用 matplotlib 画成灰度图每一列对应一个角度的滤波投影横向扫描就能快速定位坏角度。这个技巧在以后接触真实螺旋CT数据时非常有用比盯着重建图找条纹高效得多。本文还有配套的精品资源点击获取
返回列表