资讯详情

资讯详情

CT重建核心:滤波反投影(FBP)滤波器设计与调参实战

做CT、SPECT/PET重建或者工业无损检测的朋友应该都有过这种体会重建算法本身写起来并不难难的是图像质量的控制——噪声、分辨率、伪影之间的平衡最后几乎全都汇聚到一个环节上就是投影数据滤波器的设计。图像重建中常用的滤波器表面上看只是一条频率响应曲线实际上它直接决定了重建图像能不能看、能好到什么程度、能不能用于定量分析。这篇文章把我这些年在滤波反投影Filtered Backprojection, FBP里设计滤波器踩过的坑、用过的方案和排查思路整理一遍从频域原理讲到离散实现再讲到日常调参。适合刚接触重建的初学者也适合已经写了重建流程、想在质量和速度之间做取舍的工程同学。1. 重建这件事核心就是“反演”1.1 为什么先有投影后有图像很多成像模式拿到的本质都是投影数据。X射线CT采集的是X射线穿过物体后的线积分PET采集的是湮灭光子沿响应线LOR的积分工业CT则是X射线穿透工件的衰减路径积分。无论哪种方式我们手里的原始数据都不是图像本身而是一组在不同角度下测到的投影。图像重建的目的就是从这些投影中反演出衰减系数或活度分布的空间函数。这里有一个贯穿所有解析重建的基本定理叫投影切片定理Fourier Slice Theorem。它的含义很直白某角度下投影做一维傅里叶变换得到的结果正好等于原始二维图像在该角度方向上的二维傅里叶变换中的一条过中心直线。也就是说如果我们采集了足够多角度的投影就能在频域里慢慢填满整张图像的信息再做一次二维逆傅里叶变换就能得到图像。听起来很完美可真要这么做工程上有一个绕不开的问题频域里那些过原点的直线是均匀间隔采样的吗不是。低频区域靠近原点被反复覆盖采样很密高频区域远离原点只有少量角度经过采样很稀疏。如果直接把这些直线的信息拿去做逆傅里叶变换结果是一张严重模糊、带着辐射状衰减伪影的图像。这就引出了滤波器的必要性。打个比方投影数据好比用扇形的梳子去梳一个频率空间低频被多梳了几遍高频只被轻轻带过。我们要得到原图就得给高频补上适当的增益给低频适当压一压。这个增益曲线就是滤波反投影里的滤波器。1.2 滤波反投影里的“滤波”到底滤什么我们先把不带滤波的反投影结果写出来。对每个角度把投影沿着原来的方向“涂抹”回图像空间再对所有角度求和。数学上可以证明这种直接反投影得到的结果等于原始图像和一个 1/r 核做卷积其中 r 是到某个点的距离。卷积核 1/r 的频谱形状就是 1/|ω|也就是说直接反投影相当于对图像做了一次严重的低通滤波把高频信息几乎全部抹平了。滤波反投影的思路因此非常清晰在反投影之前先把每个角度的投影做一次一维傅里叶变换乘上理想的反卷积核 |ω|也就是斜坡ramp再反变换回去然后做反投影。这个 |ω| 恰好在频域里抵消掉 1/|ω| 的低通效应理论上可以恢复原始图像。到这里“滤波”这个词就落实了它滤的不是某个窄带频率而是按照频率高低重新分配权重。低频本来信息密集稍微压一点高频信息本来稀疏要放大。问题在于噪声在频域里通常是白噪声也就是各个频率都有能量。|ω| 放到高频时噪声也被同步放大而且高频放大倍数越大噪声越刺眼。所以实际使用中几乎没有谁会用纯的 |ω|都要在斜坡的基础上叠加一个低通窗函数apodization window在频率高到一定程度后把增益压下去。这个“斜坡乘窗”的组合就是图像重建里最常见的滤波器家族。1.3 从连续到离散滤波器设计为什么是工程难题公式里的 |ω| 是在连续频域定义的真正写代码时我们要处理三个问题。第一是离散化。投影数据是按探测器单元采样的频率范围从 -0.5 到 0.5 周期/像素以奈奎斯特频率为界。斜坡滤波器必须在离散频点上取值不能超出奈奎斯特频率——超出部分不存在数据强行放大只会放大混叠。第二是卷积长度。滤波本质上是一个卷积运算。直接在频域相乘等效于循环卷积而我们需要的是线性卷积。如果投影长度是 N滤波器设计长度也是 N循环卷积会把本该滑动出边界的部分绕回另一头造成所谓的“卷绕伪影”wrap-around artifact。通常做法是把投影补零到 2N 再做频域相乘。第三是归一化和零频处理。斜坡滤波器在零频处应当为零因为投影的直流分量代表物体总衰减量反投影后只贡献整体背景不需要额外放大。但离散实现时稍不注意就会在零频点留下一个错误的增益导致重建图中心出现亮点或整体背景偏移。这些看起来都是小细节实际调试起来却非常花时间。我自己的经验是滤波器设计占了重建算法里大约一半的坑而其中八成又集中在离散化、补零和归一化这几个环节上。后面我会具体讲怎么规避。2. 常用重建滤波函数逐个拆解2.1 Ram-Lak理想斜坡滤波最锋利也最“吵”Ram-Lak 是最经典的滤波器本质上就是带截至的斜坡函数H_ramp(ω) |ω|当 |ω| ≤ ω_max否则为 0。这里的 ω_max 一般取奈奎斯特频率。它在频域是一条直上直下的线没有任何平滑处理因此对高频的保留最完整图像边缘最锐利空间分辨率最高。代价也很明显由于白噪声在所有频率都有功率而 |ω| 的放大倍数随频率线性上升Ram-Lak 会把高频噪声放到最大。另外在截止点 ω_max 处频率响应是突然跳变到 0 的这种硬截断在空间域会形成旁瓣振荡表现在图像里就是高对比度边缘附近出现环状或线状振铃。Ram-Lak 的离散脉冲响应也值得记住h(0)1/4偶数点处 h0奇数点处 h-1/(π²n²)。它的能量衰减很慢属于慢拖尾型脉冲响应这也是它和“振铃”绑定的原因。什么时候用高剂量、高信噪比的CT数据或者做分辨率模体测量、想看系统极限分辨率时我才会直接用它。日常临床或普通检测场景我几乎不用原始 Ram-Lak。如果你在复查一段重建代码发现图像边缘出现一圈圈黑白条纹第一反应就该去查是不是用了 Ram-Lak 且没有叠加任何窗。2.2 Shepp-Logan给高频装上一个温顺的“刹车”Shepp-Logan 滤波器在斜坡上乘了一个 sinc 窗H_SL(ω) |ω| · [sin(π ω / (2 ω_max)) / (π ω / (2 ω_max))]这个窗函数的形状是从 1 衰减到 0 的圆弧样曲线。相比 Ram-Lak 的直上直下它在高频区平滑了很多高频分量被明显抑制因此重建图像噪声更小振铃也更轻微。代价是分辨率有所下降高频被压掉之后图像边缘的过渡会稍微变缓细微结构对比度降低。但在绝大多数CT应用里这个损失是值得的。对相同分辨率来说Shepp-Logan 的噪声方差通常比 Ram-Lak 低不少。我个人的默认配置就是如果没有特殊理由CT 重建先用 Shepp-Logan 打底。它是“性能平庸但哪儿都不出错”的代表Commercial 系统里也经常见到它的身影。你用 Shepp-Logan 重建出来的图像边缘不会像 Ram-Lak 那么透亮但也不会出现让人心烦的噪声颗粒和边缘振铃。2.3 余弦、Hamming、Hann从“陡截止”到“圆滑下降”这一组滤波器都是在斜坡上乘一个常见的谱分析窗函数差别在窗的形状。余弦窗Cosine: H_cos(ω) |ω| · cos(π ω / (2 ω_max))。它在截止点处平滑归零高频抑制介于 Shepp-Logan 和 Hann 之间振铃很小噪声水平适中。在很多CT机上叫“soft”或“smooth”类选择里都能看到它。Hann 窗也叫 Hanning: H_hann(ω) |ω| · [0.5 0.5 cos(π ω / ω_max)]。它在截止点处不仅数值归零导数也是零因此过渡最平滑对旁瓣的抑制最好噪声水平最低但分辨率也最差。适合低剂量扫描、肥胖体部扫描或噪声本来就大的数据。Hamming 窗: H_hamming(ω) |ω| · [0.54 0.46 cos(π ω / ω_max)]。它和 Hann 的区别在于截止点处留了一个小的残余值不是完全归零。这个小小的残余带来一个好处高频比 Hann 保留略多一点图像比 Hann 锐利同时又比 Ram-Lak 平稳。但它也因为这个残余在截止频率处会有轻微的频率响应跳变理论上振铃会比 Hann 稍微重一点。这三者选哪个本质上是在问你能接受多少噪声。统计上Hamming 的分辨率明显优于 Hann噪声又低于 Ram-LakHann 则是最圆滑的保守选择。对迭代重建出现之前的传统CTHann 常被用于噪声敏感部位Cosine 是折衷Hamming 则适合需要兼顾细节与噪声的一般场景。2.4 Butterworth 和高斯两个更适合“精细塑形”的选项Butterworth 滤波器在核医学SPECT和某些工业CT里很常见。它的窗函数写成W_butterworth(ω) 1 / sqrt(1 (ω / ω_c)^(2n))其中 ω_c 是截止频率n 是阶数。这个窗的频率响应非常平直通带内几乎不衰减到了截止频率附近才开始快速下降阶数越高下降越陡。它的好处是你可以把“多少频率保留”和“多快切掉”分开控制适合在噪声谱有明显形状的情况下精细设计。缺点是需要调两个参数参数之间还会耦合不熟悉的人容易调出要么太糊、要么太噪的结果。高斯窗也值得提一句W_gauss(ω) exp(-ω² / (2σ²))。高斯函数最突出的特点是它的傅里叶变换还是高斯函数也就是说频域光滑空间域也光滑不会产生振铃。它没有旁瓣频率响应单调下降非常适合低计数、统计噪声主导的数据比如PET和高噪声SPECT。但它没有严格截止频率理论上支撑无限长实际使用时要截断截断长度选得不好会重新引入瑕疵。2.5 一张表看懂这些滤波器怎么选滤波器截止处行为高频抑制程度空间分辨率噪声水平典型应用场景Ram-Lak硬截止无额外抑制最高最大高剂量CT、分辨率模体、算法验证Shepp-Logan平滑衰减中等高中低常规CT默认选项Cosine平滑归零中等偏强中高中低低剂量CT、需要均衡场景Hann零值且导数为零强中低低剂量、低信噪比数据Hamming保留小残余中强中高中平衡分辨与噪声Butterworth可控陡降可调可调可调SPECT、需要独立控制截止和陡度的场合Gaussian无旁瓣单调下降强中低很低PET/SPECT、强噪声数据这张表并不是标准答案。实际选型一定要结合你自己的探测器响应、剂量和临床需求我后面会讲怎么量化验证。3. 滤波器设计里的核心参数与取舍逻辑3.1 截止频率比选窗函数更关键很多初学者把注意力全放在“用哪个窗函数”上却忽略了截止频率。其实对最终图像影响最大的往往就是截止频率设多高。先算一个基础量奈奎斯特频率 f_N 1 / (2Δ)其中 Δ 是探测器采样间距换算到图像域的尺寸。如果探测器中心距是 0.5 mm那么奈奎斯特频率就是 1 lp/mm。比它更高的空间频率根本不存在不用考虑。所谓截止频率一般写成奈奎斯特频率的百分比。比如 cutoff0.8意思是把 0.8 lp/mm 以上的频率全部置零。数据噪声大的时候把 cutoff 降到 0.5~0.7相当于主动放弃一部分高频信息换取更平滑的图像。高剂量、高信噪比数据才值得用 0.9~1.0。我给一个非常实用的经验起点普通CT用 Shepp-Logan cutoff 0.7-0.8低剂量扫描用 Hann cutoff 0.5-0.6PET 用高斯窗 cutoff 0.4-0.5按奈奎斯特计。这几个组合在多数系统上都能跑出一个“可看”的结果再根据模体测量微调。这里有个反直觉的点把 cutoff 从 1.0 降到 0.6空间分辨率并没有想象中掉得那么快因为大部分现实系统的探测器光学模糊、焦点尺寸已经天然限制了高频灵敏度你保留 0.8 以上的频率放大最多的其实是噪声和重建伪影。所以别心疼那点高频先降下来看噪声改善多少。3.2 窗函数的旁瓣能量与振铃窗函数选型本质上是在调“旁瓣能量”。频域的窗函数形状越陡峭脉冲响应的旁瓣越高振铃越明显形状越圆滑旁瓣越低但主瓣会被展宽分辨率下降。这个关系可以用一个生活化的类比理解你把一个箱子推下台阶落差越大箱子落地后弹跳越厉害如果你让它顺着斜坡滑下去落地就安静很多。频率响应在截止点处的“落差”决定了脉冲响应在空间域“弹跳”的幅度。Ram-Lak 是直接从台阶上摔下去Hann 则是顺着平滑斜坡滑到零弹跳自然小。做实际设计时我建议画一下频率响应和脉冲响应两条曲线。频率响应看截止处倾斜是否可接受脉冲响应看旁瓣振幅是否会导致振铃。如果脉冲响应的第一旁瓣已经超过主瓣的 10%重建高对比度物体如骨与软组织的界面时大概率会看到人工纹路。3.3 定量权衡分辨率、噪声和对比度到了这一步就要把“感觉”变成“数字”。两个最常用的评价量是调制传递函数MTF和噪声功率谱NPS。MTF 描述系统对各个频率的传递能力由滤波器的频率响应、探测器孔径、焦点尺寸等共同决定。你可以用细线或钨丝模体测点扩散函数再做傅里叶变换得到 MTF。常用指标是 50% MTF 或 10% MTF 对应的频率比如 10% MTF 0.6 lp/mm就说明 0.6 lp/mm 以上的细节基本看不出来了。噪声功率谱描述噪声在频域的分布。白噪声经过 |ω| 斜坡放大后NPS 会呈现随频率平方上升的形状因此高频噪声能量非常大。评价噪声不能只看像素标准差要看 NPS 的形态如果 NPS 在低频段明显抬高说明出现了低频伪影或DC漂移如果仍然在高频段翘起说明还得多压高频。滤波器的噪声放大可以通过计算窗函数平方在通带内的积分近似估计。同样截止频率下Hann 的噪声放大通常比 Shepp-Logan 低 20%~40%这个差异在低剂量扫描里体感非常明显。选型的时候我会先测一组同一数据下不同滤波器的 NPS 和 MTF然后用“10% MTF 频率除以噪声标准差”这样一个简单比值做综合评分分数高的组合优先。3.4 别忘了系统和滤波器的联合作用滤波器从来不是孤立存在的。探测器有限孔径带来几何模糊X射线焦点尺寸、探测器像素串扰、散射都会进一步压低高频响应。也就是说即使你把滤波器 cutoff 设成 1.0系统物理极限可能早就把 0.8 lp/mm 以上的信号磨没了。因此一个严谨的滤波器设计流程应该先把系统自身的 MTF 测出来然后问一个问题我们期望的输出图像分辨率是多少如果目标是在 10% MTF 处达到 0.7 lp/mm而系统在 0.7 lp/mm 处已经衰减到 0.4那么滤波器就必须在这个频率附近补偿一部分增益。不过要小心补偿的只是信号噪声也在同样频率上被放大补偿越多噪声越重最终可能得不偿失。实战中我更倾向于做“频带匹配”把滤波器的截止频率放在系统 MTF 开始明显下降的位置附近也就是信噪比尚有正向余量的频带里。再往上数据里已经几乎没有可恢复的信号纯粹是噪声区果断切掉。这个思路比单纯调窗函数更系统也更容易解释为什么某种参数组合在某个机型上成立、在另一台机器上不成立。4. 实操把滤波器真正写进重建流程4.1 一维频域滤波的标准流程滤波反投影里的一维滤波可以放在重投影后的频域里完成也可以用空间卷积完成。频域实现通常更快但要注意循环卷积问题。标准流程是这样读入某一角度的投影序列长度 N。对投影做补零长度扩展到 2N目的是把线性卷积变成“伪线性”避免卷绕伪影。对补零后的序列做 FFT得到频谱 P(ω)。构造滤波函数 F(ω)长度同样为 2N或在频域直接采样 2N 个点与 P(ω) 逐点相乘。做 IFFT取前 N 个点作为该角度的滤波后投影。后续按常规反投影流程处理。这里有一个经验细节滤波器可以在空间域设计然后和投影做空间卷积也可以直接在频域设计。我推荐用频域设计因为窗函数表达式在频域写起来更直观而且 N 变化时不需要重新推导卷积核。4.2 一个可直接改的滤波函数实现下面这段代码是我自己项目里的精简版用 Python 风格写的。它根据输入长度和类型生成一维频域滤波响应。import numpy as np def ramp_window(n, kindshepp_logan, cutoff1.0): 生成频域滤波响应。 n: 投影长度实际用 2N 传入 cutoff: 截止频率1.0 表示截止到奈奎斯特频率0.5 cycles/pixel u np.fft.fftfreq(n) # 频率单位 cycles/pixel范围约 [-0.5, 0.5] u_abs np.abs(u) ramp u_abs.copy() mask u_abs cutoff * 0.5 1e-12 if kind ramp: filt ramp elif kind shepp_logan: # sin(pi*u/(2*c)) / (pi*u/(2*c))c 为归一化截止 c cutoff * 0.5 filt ramp * np.sinc(u / (2 * c)) elif kind cosine: c cutoff * 0.5 filt ramp * np.cos(np.pi * u / (2 * c)) elif kind hann: c cutoff * 0.5 filt ramp * (0.5 0.5 * np.cos(np.pi * np.clip(u / c, -1, 1))) elif kind hamming: c cutoff * 0.5 filt ramp * (0.54 0.46 * np.cos(np.pi * np.clip(u / c, -1, 1))) elif kind gaussian: sigma 0.3 * cutoff filt ramp * np.exp(-(u ** 2) / (2 * sigma ** 2)) else: raise ValueError(unknown filter type) # 超过截止频率的部分清成 0避免噪声放大 filt[~mask] 0.0 return filt用的时候注意传入的 n 应该是补零后的长度 2N。fftfreq 返回的单位是 cycles/pixel奈奎斯特是 0.5所以 cutoff1.0 表示保留到奈奎斯特极限。如果你希望截止频率对应 0.8 lp/mm而像素尺寸是 0.5 mm那么奈奎斯特是 1.0 lp/mmcutoff 就填 0.8。这段代码有几个我故意留下的处理一是 np.clip 防止 cos 参数越界导致 NaN二是先把超截止区域置零三是高斯窗里用 sigma 来控制窗宽实际项目里建议把 sigma 单独暴露为参数不要和 cutoff 绑定。你可以根据平台要求改写成 C/C 或 CUDA逻辑完全一样。4.3 三个最容易出错的细节第一个是补零长度。如果你直接在长度为 N 的投影上做 FFT 再乘长度为 N 的滤波函数得到的是循环卷积结果投影首尾会互相污染。这种污染在重建图像上表现为从物体边界往外辐射的细线伪影很容易被误认为数据有问题。至少要补到 2N我用过 4N效果没有本质提升但速度和内存开销更差。第二个是零频和负频率处理。fftfreq 返回的数组包括零频、正频率和负频率。斜坡的绝对值运算必须对每个频点逐个做不能让负频率段变成负数。如果滤波函数在某处出现负值意味着该频率分量被反相图像上会出现非常奇怪的“挖空”感。我曾经在移植代码时把 abs 写漏了一次重建出来的图像中央直接出现一个黑洞排查很久才发现是负频率段乘错了。第三个是归一化因子。FFT 的缩放约定在不同库里有差异。numpy 的 ifft 自带 1/N 缩放但有些库不做。你需要在流程里明确一次缩放并写进注释否则升级依赖库后图像整体对比度可能莫名其妙地变化。常用的方式是滤波后做 ifft然后乘 2因为前面补零再取前半段再在反投影累加时统一做 1/角度数 的归一化。具体比例以你的模体标定为准关键是全流程只做一次统一缩放。4.4 到扇束、锥束重建里要注意什么如果你做的是扇束CT或锥束FDK重建一维滤波的方向和几何权重会变。扇束重建中用的是“加权滤波反投影”滤波作用在探测器投影序列上位置不变但反投影前要做射线几何权重补偿。锥束FDK里每一行探测器投影用同样的斜坡窗滤波滤波后还要做体素相关的锥角权重。在实际写FDK时我有一个反复踩过的坑滤波函数必须沿着“探测器单元索引”方向做不能沿着图像行方向做。探测器单元间距可能是等角采样的对应的频域含义和图像像素间距完全不同。如果你把图像域的频率当成滤波器的频率单位重建图像的分辨率和噪声会一起坏掉而且是那种找不到逻辑的坏——因为它在不同位置表现不一样。另一个注意事项是3D重建里的轴向z方向通常不推荐做同样的斜坡滤波。斜坡滤波是为了补偿二维扇束投影的频域采样不均匀z方向是直线采样直接放大高频会把层间噪声拉高。我的习惯是 z 方向最多做一个轻度的平滑比如高斯窗或者干脆不处理。5. 常见问题与排查心得5.1 图像边缘出现一圈圈振铃这是最经典的滤波器问题。原因通常是截止频率处响应太陡。先确认是否用了 Ram-Lak再用 Hann 或 Shepp-Logan 替换对比。如果振铃还在可能是滤波器长度不够导致吉布斯现象试着把滤波器过渡带加宽或者降低 cutoff。另一个隐蔽原因是窗函数里用了硬截断比如对窗函数本身直接 ifft 截断这种截断本身就会引入旁瓣正确的做法是在频域里连续采样、不要对窗函数做空间截断。5.2 噪声颗粒大到无法接受看到满屏雪花点第一反应不是降低重建参数而是测一下你实际的有效截止频率。很多人代码里写了一个 cutoff0.6但因为频率数组单位不对实际乘到了 1.0 甚至更高。我建议在滤波函数里加一个简单的自检把 cutoff 对应的频点打印出来对照像素尺寸换算成 lp/mm和你的系统预期值对比。排除单位问题后再尝试 Hann 或高斯窗并逐步降低 cutoff。如果噪声仍然集中在高频检查是不是探测器坏像素太多坏像素会产生脉冲状高频噪声滤波器再平滑也救不回来。5.3 图像发雾、发糊边缘过渡变缓糊的原因一般是低通过头要么 cutoff 太低要么窗函数过于圆滑要么滤波器的归一化导致信号被整体压低。先看 MTF如果 10% MTF 频率低于你的目标就提高 cutoff 或换成陡度更大的窗Shepp-Logan 换 Ram-Lak。还有一个我遇过两次的原因滤波后取了 2N 的后半段而不是前半段等效于把投影平移了半个窗口图像整体被额外的空间低通模糊掉。检查取区间逻辑确认取的是卷积结果的有效部分。5.4 图像中央出现亮点或整体背景漂移这几乎一定是零频处理问题。斜坡滤波器在零频处应为 0但如果你构造滤波函数时用 np.linspace 生成频率数组后在零频点不小心赋了一个非零值或者 FFT 结果没做 fftshift 导致零频位置错位投影的直流分量就会被错误放大。排查方法是把滤波函数单独打印出来看 zero 频点是不是 0再看滤波后投影的均值正常应该是接近 0 的因为斜坡滤波会把直流偏置滤掉。如果均值不为 0问题就出在频域零频点处理上。5.5 从物体边缘向外辐射的细线伪影这类伪影优先怀疑补零不足。把补零长度从 N 改成 2N 试试如果显著改善就是循环卷积的问题。另外一个来源是角度采样不足但那种伪影是“放射状星芒”不是“从物体边缘往外拖尾”。区分两者的办法是旋转物体方向重测如果伪影方向始终跟随物体边缘多半是滤波器卷积问题如果伪影方向固定多半是系统几何或角度采样问题。5.6 一张问题速查表现象首要怀疑处理方向边缘振铃截止过陡换HHann/Shepp-Logan或降cutoff高频噪声重截止过高或窗太平降cutoff换高斯/Hann图像发糊低通过度提cutoff换Shepp-Logan/Ram-Lak中央亮点/背景漂移零频处理错误检查零频点置零边缘拖尾细线补零不足投影补零到2N固定方向星芒角度采样不足/几何错误检查扫描角度数和几何标定6. 我自己在选型上的经验习惯最后分享一点个人习惯。接到一个新重建项目我不会一上来就精调窗口参数而是先定三件事数据信噪比大概什么水平、目标空间分辨率是多少、重建时间允许多久。然后从 Shepp-Logan cutoff0.75 这个“安全区”开始跑出一版图像看两个指标高对比界面的振铃强度和均匀区域的噪声水平。再根据这两个指标在“更锐利”和“更平滑”之间做一次横向切换通常两次就能锁定最终参数。我还习惯把滤波器响应、补零长度、归一化因子写成一个独立函数并加断言零频点为 0、cutoff 外全 0、响应范围在 0 到奈奎斯特之间。这样无论是自己调试还是后来人接手都不至于在看不见的地方被某个细节坑掉。图像重建里的滤波器看似只是一个数学公式实际上它是连接物理系统、噪声模型和图像质量的桥梁。把这条曲线设计明白了重建算法里最难啃的那块骨头也就啃下来了。
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →