数字图像学笔记 —— 17. 图像退化与复原(自适应滤波之「最小二乘方滤波」)

news2024/12/25 0:27:19

文章目录

  • 维纳滤波的缺点
  • 约束最小二乘方滤波
  • 给一个实际例子吧

维纳滤波的缺点

维纳滤波(Wiener Filter),虽然是一种非常强大的退化图像还原算法,但是从实验过程我们也发现它存在着致命的缺陷,那就是要求输入退化系统的 F ( u , v ) F(u, v) F(u,v) 是已知的。分析维纳滤波的理论,我们发现它是通过给未知的退化系统 H ( u , v ) H(u, v) H(u,v) 输入原始数据 F ( u , v ) F(u, v) F(u,v) 后得到退化后的数据 G ( u , v ) G(u, v) G(u,v),然后通过梯度下降算法与约束条件 M S E = ( F ^ − F ) 2 MSE = (\hat F - F)^2 MSE=(F^F)2 找出最接近 H ( u , v ) H(u, v) H(u,v) 的近似解 H ^ ( u , v ) \hat H(u, v) H^(u,v)

所以,在求解 H ^ ( u , v ) \hat H(u, v) H^(u,v) 的过程中,它要求我们必须要预先知道输入退化系统 H ( u , v ) H(u, v) H(u,v) 前的图像 F ( u , v ) F(u, v) F(u,v) 的信息。但是我们在实际生活中,不一定都能拥有这样的条件,可以让我们顺利的找出 H ^ ( u , v ) \hat H(u, v) H^(u,v),比方说对老照片的恢复。

在这里插入图片描述

拍摄时间未详的居里夫人照,图片来自网络。

对于上述图片来说,存在很明显的退化,而且我们也没有办法回到过去拍一张高清的照片,所以要想让照片得到修复,获得较好的视觉感受,维纳滤波在这里就无法使用了。

所以后续的研究者提出了一种仅依赖均值与方差便可以估算出最好复原效果的 「约束最小二乘方滤波(Constrained Least Squares Filtering)」。接下来我们来看看它是如何展现「化腐朽为神奇」的力量。

约束最小二乘方滤波

先从退化模型的一般形式(频率空间)出发,我们有如下公式:

G ( u , v ) = H ( u , v ) F ( u , v ) + N ( u , v ) G(u, v) = H(u, v) F(u, v) + N(u, v) G(u,v)=H(u,v)F(u,v)+N(u,v)

这是因为噪声函数 D ( u , v ) D(u, v) D(u,v) 通常被认为与输入图像 F ( u , v ) F(u, v) F(u,v) 相互独立,即 N ( u , v ) N(u, v) N(u,v) 是与 F ( u , v ) F(u, v) F(u,v) 不相关的噪声函数。因此,这个模型可以被看作是信号 F ( u , v ) F(u, v) F(u,v) 通过退化系统 H ( u , v ) H(u, v) H(u,v) 得到观测结果 G ( u , v ) G(u, v) G(u,v) 后再加上噪声 N ( u , v ) N(u, v) N(u,v)

接下来,我们可以使用最小二乘法来估计未知的退化函数 H ( u , v ) H(u, v) H(u,v)。最小二乘法是一种优化方法,可以通过最小化误差平方和来得到最优解。具体来说,对于给定的观测结果 G ( u , v ) G(u, v) G(u,v) 和输入图像 F ( u , v ) F(u, v) F(u,v),我们可以计算出它们之间的误差平方和:

E = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ G ( u , v ) 2 − H ( u , v ) F ( u , v ) ∣ 2 E = \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | G(u, v)^2 - H(u,v)F(u,v)|^2 E=u=0M1v=0N1G(u,v)2H(u,v)F(u,v)2

我们的目标是找到一个 H ( u , v ) H(u,v) H(u,v),使得误差平方和 E E E 最小化。因此,我们可以通过求解下面的优化问题来得到最优的 H ( u , v ) H(u,v) H(u,v)

min ⁡ H ( u , v ) E = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ G ( u , v ) 2 − H ( u , v ) F ( u , v ) ∣ 2 \min_{H(u,v)}E = \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | G(u, v)^2 - H(u,v)F(u,v)|^2 H(u,v)minE=u=0M1v=0N1G(u,v)2H(u,v)F(u,v)2

然而,我们需要加上一些约束条件,以防止退化函数 H ( u , v ) H(u,v) H(u,v) 变得过于复杂,从而导致过度拟合(overfitting)。一个常见的约束条件是 H ( u , v ) H(u,v) H(u,v) 的均值为 1 1 1,方差为 σ H 2 \sigma^2_H σH2,即:

1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 H ( u , v ) = 1 \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} H(u, v) = 1 MN1u=0M1v=0N1H(u,v)=1

1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ H ( u , v ) − 1 M N ∑ i = 0 M − 1 ∑ j = 0 N − 1 H ( i , j ) ∣ 2 = σ H 2 \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | H(u,v) - \frac{1}{MN} \sum_{i=0}^{M-1}\sum_{j=0}^{N-1}H(i, j) |^2 = \sigma^2_{H} MN1u=0M1v=0N1H(u,v)MN1i=0M1j=0N1H(i,j)2=σH2

这个约束条件的意义是,我们期望退化函数 H ( u , v ) H(u,v) H(u,v) 的均值为 1 1 1,这是因为退化函数的总能量应该保持不变。同时,方差 σ H 2 \sigma^2_H σH2 用于限制 H ( u , v ) H(u,v) H(u,v) 的复杂度,从而避免过度拟合。

根据上述约束条件,我们可以得到一个带有约束的最小二乘优化问题:

min ⁡ H ( u , v ) E = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ G ( u , v ) − H ( u , v ) F ( u , v ) ∣ 2 \min_{H(u,v)} E = \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | G(u,v) - H(u, v)F(u, v) |^2 H(u,v)minE=u=0M1v=0N1G(u,v)H(u,v)F(u,v)2

s . t . 1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 H ( u , v ) = 1 s.t. \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} H(u, v) = 1 s.t.MN1u=0M1v=0N1H(u,v)=1

1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ H ( u , v ) − 1 M N ∑ i = 0 M − 1 ∑ j = 0 N − 1 H ( i , j ) ∣ 2 = σ H 2 \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i=0}^{M-1} \sum_{j=0}^{N-1} H(i, j) |^2 = \sigma_{H}^{2} MN1u=0M1v=0N1H(u,v)MN1i=0M1j=0N1H(i,j)2=σH2

上面这个优化问题可以通过拉格朗日乘子法来求解。我们首先构造一个拉格朗日函数:

L ( H , λ 1 , λ 2 ) = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ G ( u , v ) − H ( u , v ) F ( u , v ) ∣ 2 + λ 1 ( 1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 H ( u , v ) − 1 ) + λ 2 ( 1 M N ∑ u 0 M − 1 ∑ v = 0 N − 1 ∣ H ( u , v ) − 1 M N ∑ i = 0 M − 1 ∑ j = 0 N − 1 H ( i , j ) ∣ 2 − σ H 2 ) L(H, \lambda_1, \lambda_2) = \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | G(u, v) - H(u, v)F(u, v)|^2 + \lambda_1 (\frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} H(u, v) - 1) + \lambda_2 (\frac{1}{MN} \sum_{u_0}^{M-1} \sum_{v=0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i=0}^{M-1} \sum_{j=0}^{N-1} H(i, j)|^2 - \sigma_H^2) L(H,λ1,λ2)=u=0M1v=0N1G(u,v)H(u,v)F(u,v)2+λ1(MN1u=0M1v=0N1H(u,v)1)+λ2(MN1u0M1v=0N1H(u,v)MN1i=0M1j=0N1H(i,j)2σH2)

其中 λ 1 \lambda_1 λ1 λ 2 \lambda_2 λ2 是拉格朗日乘子。接下来,我们对 L ( H , λ 1 , λ 2 ) L(H, \lambda_1, \lambda_2) L(H,λ1,λ2) 分别对 H ( u , v ) H(u,v) H(u,v) λ 1 \lambda_1 λ1 λ 2 \lambda_2 λ2 求偏导数,并令它们等于 0 0 0,可以得到下面的一组方程:

∂ L ∂ λ 1 = 1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 H ( u , v ) − 1 = 0 \frac{\partial L}{\partial \lambda_1} = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} H(u, v) - 1 =0 λ1L=MN1u=0M1v=0N1H(u,v)1=0

∂ L ∂ H ( u , v ) = − 2 F ( u , v ) ( G ( u , v ) − H ( u , v ) F ( u , v ) ) + λ 1 + 2 λ 2 ( H ( u , v ) − 1 M N ∑ i = 0 M − 1 ∑ j = 0 N − 1 H ( i , j ) ) = 0 \frac{\partial L}{\partial H(u, v)} = -2 F(u, v)( G(u, v) - H(u, v)F(u, v) ) + \lambda_1 + 2 \lambda_2 ( H(u, v) - \frac{1}{MN} \sum_{i=0}^{M-1} \sum_{j=0}^{N-1} H(i, j)) = 0 H(u,v)L=2F(u,v)(G(u,v)H(u,v)F(u,v))+λ1+2λ2(H(u,v)MN1i=0M1j=0N1H(i,j))=0

∂ L ∂ λ 2 = 1 M N ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ H ( u , v ) − 1 M N ∑ i = 0 M − 1 ∑ j = 0 N − 1 H ( i , j ) ∣ 2 − σ H 2 = 0 \frac{\partial L}{ \partial \lambda_2} = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i=0}^{M-1} \sum_{j=0}^{N-1} H(i ,j) |^2 - \sigma_H^2 = 0 λ2L=MN1u=0M1v=0N1H(u,v)MN1i=0M1j=0N1H(i,j)2σH2=0

解这个方程组可以得到最优的 H ( u , v ) H(u,v) H(u,v)。具体地,我们可以先把 λ 1 \lambda_1 λ1 λ 2 \lambda_2 λ2 消去,然后将 H ( u , v ) H(u,v) H(u,v) 写成矩阵形式 H = [ h 1 , h 2 , . . . , h L ] H = [h_1, h_2, ..., h_L] H=[h1,h2,...,hL],其中 L = M N L=MN L=MN。接着,将偏导数为 0 0 0 的方程写成矩阵形式 A h = b A\textbf{h} = \textbf{b} Ah=b,其中 A A A b \textbf{b} b 可以根据上述方程计算得到, h \textbf{h} h 是一个列向量,它包含了所有的 h i h_i hi。最后,我们可以通过求解下面的矩阵方程来得到最优的 h \textbf{h} h

( A T A + λ I ) h = A T g (A^{T} A + \lambda I) \mathbf h = A^T \mathbf g (ATA+λI)h=ATg

其中 g \textbf{g} g G ( u , v ) G(u,v) G(u,v) 的列向量形式, I I I 是单位矩阵, λ \lambda λ 是一个正则化参数,用于平衡拟合误差和模型复杂度。

最终,我们可以得到最优的 H ( u , v ) H(u,v) H(u,v)

H ( u , v ) = h u N + v H(u, v) = h_{uN + v} H(u,v)=huN+v

其中 h i h_i hi 是列向量 h \textbf{h} h 中的第 i i i 个元素。根据最优的 H ( u , v ) H(u,v) H(u,v),我们可以得到修复后的图像 F r e s t o r e d ( u , v ) F_{restored}(u,v) Frestored(u,v)

F r e s t o r e d ( u , v ) = F − 1 { G ( u , v ) / H ( u , v ) } F_{restored} (u, v) = \mathcal{F}^{-1} \{ G(u,v) / H(u,v) \} Frestored(u,v)=F1{G(u,v)/H(u,v)}

其中 F − 1 \mathcal{F}^{-1} F1 是傅里叶逆变换。由于在实际应用中, H ( u , v ) H(u,v) H(u,v) 可能会变得过于复杂,从而导致过度拟合。为了避免这种情况,我们可以使用正则化技术,例如 Tikhonov 正则化,来控制模型的复杂度。具体来说,我们可以将优化问题改写为:

min ⁡ H ( u , v ) E = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ G ( u , v ) − H ( u , v ) F ( u , v ) ∣ 2 + λ ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ H ( u , v ) ∣ 2 \min_{H(u, v)} E = \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} | G(u, v) - H(u, v)F(u, v) |^2 + \lambda \sum_{u = 0}^{M-1} \sum_{v=0}^{N-1} | H(u, v)|^2 H(u,v)minE=u=0M1v=0N1G(u,v)H(u,v)F(u,v)2+λu=0M1v=0N1H(u,v)2

其中 λ \lambda λ 是一个正则化参数,用于控制 H ( u , v ) H(u,v) H(u,v) 的平滑度。这样,我们可以得到下面的矩阵方程:

( A T A + λ I ) h = A T g (A^T A + \lambda I) \textbf h = A^T \textbf g (ATA+λI)h=ATg

解这个方程可以得到带有正则化项的最优 H ( u , v ) H(u,v) H(u,v),从而得到修复后的图像 F r e s t o r e d ( u , v ) F_{restored}(u,v) Frestored(u,v)

总之,约束最小二乘方滤波是一种非常强大的图像修复技术,它可以通过仅依赖于均值和方差的约束条件来估计未知的退化函数,并从退化的图像中恢复出尽可能接近原始图像的图像。相比于维纳滤波,它更具有实用性和适用性。

给一个实际例子吧

上述内容实在太多,但是OpenCV已经帮我们减少了很多工作量,所以就可以简化为几个关键函数,所以我们就可以这样去调用它们

def ConstrainedLeastSquaresFiltering(image_path):
    # Load the image as grayscale
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)

    # Add Gaussian noise to the image
    noisy_img = img + np.random.normal(0, 20, size=img.shape)

    # Define the degradation matrix
    degradation_matrix = np.array([[0.9, 0.1, 0], [0.1, 0.8, 0.1], [0, 0.1, 0.9]])

    # Define the constraint matrix
    constraint_matrix = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])

    # Define the regularization parameter
    alpha = 0.1

    # Compute the inverse of the degradation matrix
    degradation_matrix_inv = np.linalg.inv(degradation_matrix)

    # Compute the constrained least squares estimate
    cls_estimate = cv2.filter2D(noisy_img, -1, degradation_matrix_inv)
    cls_estimate = cls_estimate - alpha * cv2.filter2D(cls_estimate, -1, constraint_matrix)

    # Convert the image to uint8 and clip it to [0, 255]
    cls_estimate = np.clip(cls_estimate, 0, 255).astype(np.uint8)

    # Show the original image, the noisy image, and the restored image
    cv2.imshow("Original Image", img)
    cv2.imshow("Noisy Image", noisy_img)
    cv2.imshow("Restored Image", cls_estimate)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

那么自然输出的效果就变成下面这个样子了

在这里插入图片描述

怎么样,是不是很厉害的技术呢?

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/389379.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

医疗器械之模糊算法(嵌入式部分)

模糊控制 所谓模糊控制,就是对难以用已有规律描述的复杂系统,采用自然语言(如大,中,小)加以描述,借助定性的,不精确的以及模糊的条件语句来表达,模糊控制是一种基于语言的…

Java虚拟机JVM-运行时数据区域说明

及时编译器 HotSpot虚拟机中含有两个即时编译器,分别是编译耗时短但输出代码优化程度较低的客户端编译器(简称为C1)以及编译耗时长但输出代码优化质量也更高的服务端编译器(简称为C2),通常它们会在分层编译…

【Linux】手把手教你在CentOS上使用docker 安装MySQL8.0

文章目录前言一. docker的安装1.1 从阿里下载repo镜像1.2 安装docker1.3 启动docker并查看版本二. 使用docker安装MySQL8.02.1 拉取MySQL镜像2.2 创建容器2.3 操作MySQL容器2.4 远程登录测试总结前言 大家好,又见面了,我是沐风晓月,本文主要…

docker系列1:docker安装

传送门 docker官网地址: Docker: Accelerated, Containerized Application Development 安装地址:Install Docker Engine docker hub地址 docker hub:Docker 安装步骤 前置条件 由于安装docker,需要根据操作系统版本选择…

设计模式(只谈理解,没有代码)

1.什么是设计模式设计模式,是一套被反复使用、多数人知晓的、经过分类编目的、代码设计经验的总结。使用设计模式是为了可重用代码、让代码更容易被他人理解、保证代码可靠性、程序的重用性。2.为什么要学习设计模式看懂源代码:如果你不懂设计模式去看Jd…

【react全家桶】生命周期

文章目录04 【生命周期】1.简介2.初始化阶段2.1 constructor2.2 componentWillMount(即将废弃)2.3 static getDerivedStateFromProps(新钩子)2.4 render2.5 componentDidMount2.6 初始化阶段总结3.更新阶段3.1 componentWillRecei…

2023最新版本RabbitMQ的持久化和简单使用

上节讲了 RabbitMQ下载安装教程 , 本节主要介绍RabbitMQ的持久化和简单使用。 一、RabbitMQ消息持久化 当处理一个比较耗时得任务的时候,也许想知道消费者(consumers)是否运行到一半就挂掉。在当前的代码中,当RabbitM…

如何在本地跑FuzzBench的实验

概述 FuzzBench是谷歌做的一个评估模糊测试的benchmark,主要目的是为了评测覆盖率导向模糊测试工具能达到多少覆盖率,能够发现多少漏洞。本文主要介绍如何在本地运行fuzzbench的实验。当然也可以用谷歌的云服务来跑,具体教程在这&#xff1a…

计算机底层:BDC码

计算机底层:BDC码 BDC码的作用: 人类喜欢十进制,而机器适合二进制,因此当机器要翻译二进制给人看时,就会进行二进制和十进制的转换,而常规的转换法(k*位权)太麻烦。因此就出现了不同…

基本类型、包装类型、引用类型、String等作为实参传递后值会不会改变?

看了半天帖子,讲得乱七八糟,坑死了 [1] 先说结论 基本类型、包装类型、String类型作为参数传递之后,在方法里面修改他们的值,原值不会改变!引用类型不一定,要看是怎么修改它的。 [2] 为什么基本类型、包装类…

程序设计语言-软件设计(二十一)

数据结构与算法(二十)快速排序、堆排序(四)https://blog.csdn.net/ke1ying/article/details/129269655 这篇主要讲的是 编译与解释、文法、正规式、有限自动机、表达式、传值与传址、多种程序语言特点。 编译的过程 解释型 和 编译型 编译型过程&#…

【IDEA】IDEA使用有道翻译引擎—详细配置步骤

目录 前言 步骤一:下载翻译工具Translate 步骤二:注册登录有道云平台 步骤三:配置有道翻译 前言 2022年10月 谷歌翻译已经不在中国了,所以IDEA配置谷歌翻译会出错。 步骤一:下载翻译工具Translate 打开idea设置set…

c语言经典例题-选择结构程序设计进阶

(创作不易,感谢有你,你的支持,就是我前行的最大动力,如果看完对你有帮助,请留下您的足迹) 目录 快递费用计算: 题目: 代码思路: 代码表示: 计算一元二…

Java奠基】数组的讲解与使用

目录 数组概述 数组的定义与初始化 数值遍历 数组的常见操作 数组内存图 数组概述 数组是一种容器,可以用来存储同种数据类型的多个值,数组容器在存储数据的时候,需要结合隐式转换考虑。例如:int类型的数组容器不能存放取值…

191、【动态规划】AcWing ——AcWing 900. 整数划分:完全背包解法+加减1解法(C++版本)

题目描述 参考文章:900. 整数划分 解题思路 因为本题中规定了数字从大到小,其实也就是不论是1 2 1 4,还是2 1 1 4,都会被看作是2 1 1 4这一种情况,因此本题是在遍历中不考虑结果顺序。 背包问题中只需考虑…

mysql一条语句的写入原理

mysql写入原理 我们知道在mysql数据库最核心的大脑就是执行引擎; 其中的默认引擎Innodb在可靠执行和性能中做出来平衡; innodb支持在事务控制、读写效率,多用户并发,索引搜索方面都表现不俗; innodb如何进行数据写入…

MYSQL性能分析,Explain,Show Profile

文章目录一、MYSQL常见瓶颈二、ExplainExplain是什么三、Show Profile四、慢查询日志和全局查询日志一、MYSQL常见瓶颈 CPU: CPU饱和IO:磁盘IO速度过慢。服务器的硬件性能瓶颈。 二、Explain Explain是什么 使用explain关键字可以模拟优化器执行sql查…

YApi分析从NoSQL注入到RCE远程命令执行.md

0x00 前提 这个是前几个月的漏洞,之前爆出来发现没人分析就看了一下,也写了一片 Nosql注入的文章,最近生病在家,把这个写一半的完善一下发出来吧。 0x01 介绍 YApi是一个可本地部署的、打通前后端及QA的、可视化的接口管理平台…

【C++】命名空间

🏖️作者:malloc不出对象 ⛺专栏:C的学习之路 👦个人简介:一名双非本科院校大二在读的科班编程菜鸟,努力编程只为赶上各位大佬的步伐🙈🙈 目录前言一、命名空间产生的背景二、命名空…

【一天时间|简历模板】世界上最好的前端简历

一天时间系列文章是博主精心整理的面试热点问题和难点问题,吸收了大量的技术博客与面试文章,总结多年的面试经历,带你快速并高效地审视前端面试知识。直击技术痛点,主动出击,精密打击,这才是面试拿到高薪的…