用已知模糊核比较图像反卷积恢复

一张照片变模糊后,反卷积希望根据成像系统如何扩散光线,估计原来的清晰图像。本篇合并 scikit-image 0.26.x 的两个 Image Deconvolution 实例:Richardson–Lucy 迭代恢复,以及使用 Gibbs 采样自动估计正则化超参数的 Wiener–Hunt 恢复。两例都先对同一张宇航员照片施加已知模糊,再加入噪声,因此适合观察模型与参数如何影响恢复过程。

这里的点扩散函数(point spread function,PSF)由代码预先设定。它描述光学系统对一个理想点的响应。真实照片的 PSF 通常需要标定或另外估计;本教程不是从未知照片中自动猜出模糊核的盲反卷积方案。

来源与版本:根据 scikit-image team 维护的两篇官方实例编译,原页未署个人作者姓名;完整链接见文末。本文核对的是 0.26.x 文档和 v0.26.0 API 源码,核验日期为 2026-10-05。本文没有执行这些示例,以下结果图是官方页面已有的示例输出。

先固定退化模型,再讨论恢复

两篇示例均调用 data.astronaut(),再用 rgb2gray() 转为灰度。该输入为 Eileen Collins 的宇航员照片。scikit-image 数据加载器将其来源指向 NASA Great Images,并注明为公共领域、没有已知版权限制。本文保留 NASA 与 scikit-image 的归属,不把示例照片当成本次拍摄或生成的素材。

两例都使用 np.ones((5, 5)) / 25 创建 5×5 的均匀模糊核:25 个权重均为 1/25,总和为 1。scipy.signal.convolve2d(..., 'same') 返回与第一幅输入同样大小的卷积结果。示例没有显式指定边界参数,沿用该函数的默认填零边界;在边缘处观察到的恢复效果也会受这一选择影响。

后面的噪声生成与恢复方法分别成对使用。不要拿两张原页结果图直接宣布哪种算法更好:Richardson–Lucy 分支加入泊松噪声并做最大值归一化,Wiener 分支加入高斯噪声,输入尺度和显示范围也不同。

Richardson–Lucy:用迭代次数控制恢复

Richardson–Lucy 根据 PSF 逐步更新对原图的估计。原文强调迭代次数需要人工调整;scikit-image 的 API 说明也将 num_iter 视为起正则化作用的参数。增加次数不是无条件改善画质,实际使用时应同时检查细节、噪声和边缘伪影。

下面保留官方示例的完整可运行代码,未改动算法、参数或随机数初始化。首先把模糊图乘以 max_photon_count = 1000,将每个像素作为泊松分布的期望计数,再除以 1000 转回强度尺度。随后除以噪声图自身的最大值,并以 num_iter=30 调用恢复函数。

import numpy as np
import matplotlib.pyplot as plt
import skimage as ski

from scipy.signal import convolve2d as conv2


rng = np.random.default_rng()

# Convert astronaut image to grayscale
astro = ski.color.rgb2gray(ski.data.astronaut())

# Define PSF
psf = np.ones((5, 5)) / 25

# Convolve image with the PSF to simulate a blurred image
astro_blurred = conv2(astro, psf, 'same')

# Add Poisson noise to the blurred image (https://en.wikipedia.org/wiki/Shot_noise)
max_photon_count = 1000
astro_noisy = rng.poisson(astro_blurred * max_photon_count) / max_photon_count

# Normalize noisy image
astro_noisy /= np.max(astro_noisy)

# Restore image by means of deconvolution
deconvolved_RL = ski.restoration.richardson_lucy(astro_noisy, psf, num_iter=30)

fig, ax = plt.subplots(ncols=3, figsize=(8, 5))
plt.gray()

for a in (ax[0], ax[1], ax[2]):
    a.axis('off')

ax[0].imshow(astro)
ax[0].set_title('Original Data')

ax[1].imshow(astro_noisy)
ax[1].set_title('Noisy data')

ax[2].imshow(deconvolved_RL)
ax[2].set_title('Restoration using\nRichardson-Lucy')


fig.subplots_adjust(wspace=0.02, hspace=0.2, top=0.9, bottom=0.05, left=0, right=1)
plt.show()

代码中的 max_photon_count 用来设定合成泊松噪声的尺度,并不代表对真实相机进行过光子计数标定。对这张官方照片,最大值归一化用于将输入放在合适的数值范围;将示例改为任意输入时,应先检查非负性、有限性,以及最大值是否大于零,避免非法泊松参数或除零。

scikit-image 官方 Richardson–Lucy 示例:左为宇航员灰度原图,中为加入泊松噪声的模糊图,右为迭代30次的恢复结果。
官方 0.26.x 页面已有结果图,并非本次运行结果。由左至右为原始数据、噪声数据及 Richardson–Lucy 恢复。示例与排图归属 scikit-image team;基础照片来源 NASA,人物为 Eileen Collins。图片原页。

默认 clip=True 会把输出限制在 [-1, 1] 范围,以配合 scikit-image 的图像处理约定;这不是对任意物理强度单位都无损的操作。若处理科学成像数据,应先明确强度尺度,再决定是否保留这个默认值。本文没有据此改动原例。

Wiener 与自调正则化 Wiener–Hunt

第二篇原文先介绍普通 Wiener 滤波:它将 PSF、正则化先验以及数据拟合与先验之间的权衡结合起来。正则化先验惩罚高频变化,权衡参数需要人工调整。原文将这类基于线性模型的方法与 TV 等非线性恢复方法比较,认为其锐利边缘恢复能力较弱、速度较快;这一描述是原文的定性说明,不是本文测得的跨算法性能结论。

unsupervised_wiener() 则从数据中估计正则化相关的超参数。它使用迭代 Gibbs 采样,交替从图像、噪声功率和图像频率功率的后验条件分布中抽样。其图像结果是后验均值的估计,不能把“unsupervised”理解为“不需要 PSF”,更不能理解为参数对所有图片都会自动最优。

完整官方代码如下。它仍使用同一 5×5 模糊核,但对模糊图加入标准差为 0.1 * astro.std() 的高斯噪声。注意这里的 astro.std() 在卷积之后计算。

import numpy as np
import matplotlib.pyplot as plt

from skimage import color, data, restoration

rng = np.random.default_rng()

astro = color.rgb2gray(data.astronaut())
from scipy.signal import convolve2d as conv2

psf = np.ones((5, 5)) / 25
astro = conv2(astro, psf, 'same')
astro += 0.1 * astro.std() * rng.standard_normal(astro.shape)

deconvolved, _ = restoration.unsupervised_wiener(astro, psf)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5), sharex=True, sharey=True)

plt.gray()

ax[0].imshow(astro, vmin=deconvolved.min(), vmax=deconvolved.max())
ax[0].axis('off')
ax[0].set_title('Data')

ax[1].imshow(deconvolved)
ax[1].axis('off')
ax[1].set_title('Self tuned restoration')

fig.tight_layout()

plt.show()

返回值 deconvolved 是重建图;第二个返回值包含噪声与先验精度的采样链,原例用下划线丢弃它。进行算法诊断时可以保留这个字典,检查采样过程,而不是只看图像是否更清晰。

scikit-image 官方自调正则化 Wiener 示例:左侧为加入高斯噪声的模糊宇航员图,右侧为自调正则化恢复图。
官方 0.26.x 页面已有的 Wiener–Hunt 示例输出。左图显示时使用右图的最小值与最大值作为颜色范围;图中变化不能直接当作定量质量指标。示例与排图归属 scikit-image team;基础照片来源 NASA。图片原页。

这个页面虽然介绍了普通 Wiener 和自调正则化 Wiener,实际展示的代码只调用 unsupervised_wiener()。本文保留这一边界,没有补造一组普通 Wiener 的运行结果或给出未经测试的排名。

复现时需要分开控制的变量

原例两处 np.random.default_rng() 都未传入固定种子,每次生成的噪声可能不同。Wiener–Hunt 还在恢复函数内部进行随机采样,因此只固定噪声生成器,并不能同时固定恢复阶段的随机性。v0.26.0 的函数签名支持关键字参数 rng。若要做可复现的实验,可以显式为噪声和恢复各提供一个固定种子;这是编者建议的实验设计,不是原文已有代码,也未在本次执行验证。

若要公平比较算法,应使用相同的清晰输入、PSF、噪声 realization、边界处理和显示尺度,同时固定调参预算,并选择适合该噪声模型的评价方法。对于真实图片,通常没有可直接比较的清晰真值,更需要谨慎解释“锐化”后的新细节。视觉上更锐利不等于恢复了真实内容。

研究页给出的候选复现基线为 CPython 3.11、scikit-image 0.26.0、NumPy 2.2.6、SciPy 1.15.3、Matplotlib 3.10.3。本次只以 0.26.x 原文和 v0.26.0 源码核查代码语义,没有安装这些依赖、生成完整锁文件或验证二进制兼容性。这组版本不能作为运行通过的保证,也不表示它们是当前最新版本。

代码与数据边界

静态阅读这两段示例,没有发现命令执行、SQL 拼接、硬编码秘密或破坏性文件操作;这不等于完成漏洞审计。代码会导入第三方库并加载数据,大尺寸外部图片、恶意文件和依赖本身的风险不在这两段短例子的保护范围内。移植到上传服务时,需要另行限制图像尺寸、资源用量和输入格式。

data.astronaut() 经由 scikit-image 的数据加载器读取图像:先检查缓存或随发行版提供的数据,缺失时可能下载。研究页曾核验过某一发行包包含该输入,但这不能推及每一种安装方式;离线运行前应确认所需样本已存在。本次没有执行任何数据下载函数或文章算法代码,只下载了官方网页与图片用于核验和交付。

来源、参考文献与许可

本文依据以下两篇完整官方实例合并编译,保留原始算法示例,并增加版本、复现与风险说明。原站和上游代码的署名及许可仍需随稿保留。

保留的完整许可声明

以下为原项目适用许可文本,原作者与文档归属及本文改动说明见正文。

Files: *
Copyright: 2009-2022 the scikit-image team
License: BSD-3-Clause

Files: doc/source/themes/scikit-image/layout.html
Copyright: 2007-2010 the Sphinx team
License: BSD-3-Clause
Files: skimage/feature/_canny.py
       skimage/filters/edges.py
       skimage/filters/_rank_order.py
       skimage/morphology/_skeletonize.py
       skimage/morphology/tests/test_watershed.py
       skimage/morphology/watershed.py
       skimage/segmentation/heap_general.pxi
       skimage/segmentation/heap_watershed.pxi
       skimage/segmentation/_watershed.py
       skimage/segmentation/_watershed_cy.pyx
Copyright: 2003-2009 Massachusetts Institute of Technology
           2009-2011 Broad Institute
           2003 Lee Kamentsky
           2003-2005 Peter J. Verveer
License: BSD-3-Clause
Files: skimage/filters/thresholding.py
       skimage/graph/_mcp.pyx
       skimage/graph/heap.pyx
       skimage/morphology/grayreconstruct.py
Copyright: 2009-2015 Board of Regents of the University of
           Wisconsin-Madison, Broad Institute of MIT and Harvard,
           and Max Planck Institute of Molecular Cell Biology
           and Genetics
           2009 Zachary Pincus
           2009 Almar Klein
License: BSD-2-Clause
File: skimage/morphology/grayreconstruct.py
Copyright: 2003-2009 Massachusetts Institute of Technology
           2009-2011 Broad Institute
           2003 Lee Kamentsky
License: BSD-3-Clause
File: skimage/morphology/_grayreconstruct.pyx
Copyright: 2003-2009 Massachusetts Institute of Technology
           2009-2011 Broad Institute
           2003 Lee Kamentsky
           2022 Gregory Lee (added a 64-bit integer variant for large images)
License: BSD-3-Clause

File: skimage/segmentation/_expand_labels.py
Copyright: 2020 Broad Institute
           2020 CellProfiler team
License: BSD-3-Clause

File: skimage/exposure/_adapthist.py
Copyright: 1994 Karel Zuiderveld
License: BSD-3-Clause
Function: skimage/morphology/_skeletonize_various_cy.pyx:_skeletonize_loop
Copyright: 2003-2009 Massachusetts Institute of Technology
           2009-2011 Broad Institute
           2003 Lee Kamentsky
License: BSD-3-Clause

Function: skimage/_shared/version_requirements.py:_check_version
Copyright: 2013 The IPython Development Team
License: BSD-3-Clause

Function: skimage/_shared/version_requirements.py:is_installed
Copyright: 2009-2011 Pierre Raybaut
License: MIT
File: skimage/feature/_fisher_vector.py
Copyright: 2014 2014 Dan Oneata
License: MIT

File: skimage/_vendored/numpy_lookfor.py
Copyright: 2005-2023, NumPy Developers
License: BSD-3-Clause

File: skimage/transform/_thin_plate_splines.py
Copyright: 2007 Zachary Pincus
License: BSD-3-Clause

License: BSD-2-Clause
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE HOLDERS OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
License: BSD-3-Clause
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
3. Neither the name of the University nor the names of its contributors
   may be used to endorse or promote products derived from this software
   without specific prior written permission.
.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE HOLDERS OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
License: MIT

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容