武义住房和城乡建设局网站,系统优化软件哪个最好的,wordpress newsletter,网站开发学习网站文章目录维纳滤波的缺点约束最小二乘方滤波给一个实际例子吧维纳滤波的缺点
维纳滤波#xff08;Wiener Filter#xff09;#xff0c;虽然是一种非常强大的退化图像还原算法#xff0c;但是从实验过程我们也发现它存在着致命的缺陷#xff0c;那就是要求输入退化系统的 …
文章目录维纳滤波的缺点约束最小二乘方滤波给一个实际例子吧维纳滤波的缺点
维纳滤波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)然后通过梯度下降算法与约束条件 MSE(F^−F)2MSE (\hat F - F)^2MSE(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∑u0M−1∑v0N−1∣G(u,v)2−H(u,v)F(u,v)∣2E \sum_{u0}^{M-1} \sum_{v0}^{N-1} | G(u, v)^2 - H(u,v)F(u,v)|^2 Eu0∑M−1v0∑N−1∣G(u,v)2−H(u,v)F(u,v)∣2
我们的目标是找到一个 H(u,v)H(u,v)H(u,v)使得误差平方和 EEE 最小化。因此我们可以通过求解下面的优化问题来得到最优的 H(u,v)H(u,v)H(u,v)
minH(u,v)E∑u0M−1∑v0N−1∣G(u,v)2−H(u,v)F(u,v)∣2\min_{H(u,v)}E \sum_{u0}^{M-1} \sum_{v0}^{N-1} | G(u, v)^2 - H(u,v)F(u,v)|^2 H(u,v)minEu0∑M−1v0∑N−1∣G(u,v)2−H(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) 的均值为 111方差为 σH2\sigma^2_HσH2即
1MN∑u0M−1∑v0N−1H(u,v)1\frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} H(u, v) 1 MN1u0∑M−1v0∑N−1H(u,v)1
1MN∑u0M−1∑v0N−1∣H(u,v)−1MN∑i0M−1∑j0N−1H(i,j)∣2σH2\frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} | H(u,v) - \frac{1}{MN} \sum_{i0}^{M-1}\sum_{j0}^{N-1}H(i, j) |^2 \sigma^2_{H} MN1u0∑M−1v0∑N−1∣H(u,v)−MN1i0∑M−1j0∑N−1H(i,j)∣2σH2
这个约束条件的意义是我们期望退化函数 H(u,v)H(u,v)H(u,v) 的均值为 111这是因为退化函数的总能量应该保持不变。同时方差 σH2\sigma^2_HσH2 用于限制 H(u,v)H(u,v)H(u,v) 的复杂度从而避免过度拟合。
根据上述约束条件我们可以得到一个带有约束的最小二乘优化问题
minH(u,v)E∑u0M−1∑v0N−1∣G(u,v)−H(u,v)F(u,v)∣2\min_{H(u,v)} E \sum_{u0}^{M-1} \sum_{v0}^{N-1} | G(u,v) - H(u, v)F(u, v) |^2 H(u,v)minEu0∑M−1v0∑N−1∣G(u,v)−H(u,v)F(u,v)∣2
s.t.1MN∑u0M−1∑v0N−1H(u,v)1s.t. \frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} H(u, v) 1 s.t.MN1u0∑M−1v0∑N−1H(u,v)1
1MN∑u0M−1∑v0N−1∣H(u,v)−1MN∑i0M−1∑j0N−1H(i,j)∣2σH2\frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i0}^{M-1} \sum_{j0}^{N-1} H(i, j) |^2 \sigma_{H}^{2} MN1u0∑M−1v0∑N−1∣H(u,v)−MN1i0∑M−1j0∑N−1H(i,j)∣2σH2
上面这个优化问题可以通过拉格朗日乘子法来求解。我们首先构造一个拉格朗日函数
L(H,λ1,λ2)∑u0M−1∑v0N−1∣G(u,v)−H(u,v)F(u,v)∣2λ1(1MN∑u0M−1∑v0N−1H(u,v)−1)λ2(1MN∑u0M−1∑v0N−1∣H(u,v)−1MN∑i0M−1∑j0N−1H(i,j)∣2−σH2)L(H, \lambda_1, \lambda_2) \sum_{u0}^{M-1} \sum_{v0}^{N-1} | G(u, v) - H(u, v)F(u, v)|^2 \lambda_1 (\frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} H(u, v) - 1) \lambda_2 (\frac{1}{MN} \sum_{u_0}^{M-1} \sum_{v0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i0}^{M-1} \sum_{j0}^{N-1} H(i, j)|^2 - \sigma_H^2) L(H,λ1,λ2)u0∑M−1v0∑N−1∣G(u,v)−H(u,v)F(u,v)∣2λ1(MN1u0∑M−1v0∑N−1H(u,v)−1)λ2(MN1u0∑M−1v0∑N−1∣H(u,v)−MN1i0∑M−1j0∑N−1H(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 求偏导数并令它们等于 000可以得到下面的一组方程
∂L∂λ11MN∑u0M−1∑v0N−1H(u,v)−10\frac{\partial L}{\partial \lambda_1} \frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} H(u, v) - 1 0 ∂λ1∂LMN1u0∑M−1v0∑N−1H(u,v)−10
∂L∂H(u,v)−2F(u,v)(G(u,v)−H(u,v)F(u,v))λ12λ2(H(u,v)−1MN∑i0M−1∑j0N−1H(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_{i0}^{M-1} \sum_{j0}^{N-1} H(i, j)) 0 ∂H(u,v)∂L−2F(u,v)(G(u,v)−H(u,v)F(u,v))λ12λ2(H(u,v)−MN1i0∑M−1j0∑N−1H(i,j))0
∂L∂λ21MN∑u0M−1∑v0N−1∣H(u,v)−1MN∑i0M−1∑j0N−1H(i,j)∣2−σH20\frac{\partial L}{ \partial \lambda_2} \frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} | H(u, v) - \frac{1}{MN} \sum_{i0}^{M-1} \sum_{j0}^{N-1} H(i ,j) |^2 - \sigma_H^2 0 ∂λ2∂LMN1u0∑M−1v0∑N−1∣H(u,v)−MN1i0∑M−1j0∑N−1H(i,j)∣2−σH20
解这个方程组可以得到最优的 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[h1,h2,...,hL]H [h_1, h_2, ..., h_L]H[h1,h2,...,hL]其中 LMNLMNLMN。接着将偏导数为 000 的方程写成矩阵形式 AhbA\textbf{h} \textbf{b}Ahb其中 AAA 和 b\textbf{b}b 可以根据上述方程计算得到h\textbf{h}h 是一个列向量它包含了所有的 hih_ihi。最后我们可以通过求解下面的矩阵方程来得到最优的 h\textbf{h}h
(ATAλI)hATg(A^{T} A \lambda I) \mathbf h A^T \mathbf g (ATAλI)hATg
其中 g\textbf{g}g 是 G(u,v)G(u,v)G(u,v) 的列向量形式III 是单位矩阵λ\lambdaλ 是一个正则化参数用于平衡拟合误差和模型复杂度。
最终我们可以得到最优的 H(u,v)H(u,v)H(u,v)
H(u,v)huNvH(u, v) h_{uN v} H(u,v)huNv
其中 hih_ihi 是列向量 h\textbf{h}h 中的第 iii 个元素。根据最优的 H(u,v)H(u,v)H(u,v)我们可以得到修复后的图像 Frestored(u,v)F_{restored}(u,v)Frestored(u,v)
Frestored(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)F−1{G(u,v)/H(u,v)}
其中 F−1\mathcal{F}^{-1}F−1 是傅里叶逆变换。由于在实际应用中H(u,v)H(u,v)H(u,v) 可能会变得过于复杂从而导致过度拟合。为了避免这种情况我们可以使用正则化技术例如 Tikhonov 正则化来控制模型的复杂度。具体来说我们可以将优化问题改写为
minH(u,v)E∑u0M−1∑v0N−1∣G(u,v)−H(u,v)F(u,v)∣2λ∑u0M−1∑v0N−1∣H(u,v)∣2\min_{H(u, v)} E \sum_{u0}^{M-1} \sum_{v0}^{N-1} | G(u, v) - H(u, v)F(u, v) |^2 \lambda \sum_{u 0}^{M-1} \sum_{v0}^{N-1} | H(u, v)|^2 H(u,v)minEu0∑M−1v0∑N−1∣G(u,v)−H(u,v)F(u,v)∣2λu0∑M−1v0∑N−1∣H(u,v)∣2
其中 λ\lambdaλ 是一个正则化参数用于控制 H(u,v)H(u,v)H(u,v) 的平滑度。这样我们可以得到下面的矩阵方程
(ATAλI)hATg(A^T A \lambda I) \textbf h A^T \textbf g (ATAλI)hATg
解这个方程可以得到带有正则化项的最优 H(u,v)H(u,v)H(u,v)从而得到修复后的图像 Frestored(u,v)F_{restored}(u,v)Frestored(u,v)。
总之约束最小二乘方滤波是一种非常强大的图像修复技术它可以通过仅依赖于均值和方差的约束条件来估计未知的退化函数并从退化的图像中恢复出尽可能接近原始图像的图像。相比于维纳滤波它更具有实用性和适用性。
给一个实际例子吧
上述内容实在太多但是OpenCV已经帮我们减少了很多工作量所以就可以简化为几个关键函数所以我们就可以这样去调用它们
def ConstrainedLeastSquaresFiltering(image_path):# Load the image as grayscaleimg cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)# Add Gaussian noise to the imagenoisy_img img np.random.normal(0, 20, sizeimg.shape)# Define the degradation matrixdegradation_matrix np.array([[0.9, 0.1, 0], [0.1, 0.8, 0.1], [0, 0.1, 0.9]])# Define the constraint matrixconstraint_matrix np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])# Define the regularization parameteralpha 0.1# Compute the inverse of the degradation matrixdegradation_matrix_inv np.linalg.inv(degradation_matrix)# Compute the constrained least squares estimatecls_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 imagecv2.imshow(Original Image, img)cv2.imshow(Noisy Image, noisy_img)cv2.imshow(Restored Image, cls_estimate)cv2.waitKey(0)cv2.destroyAllWindows()那么自然输出的效果就变成下面这个样子了 怎么样是不是很厉害的技术呢
