【youcans 的 OpenCV 例程 300 篇】108. 约束最小二乘方滤波

6. 退化图像复原

图像复原是对图像退化的过程进行估计，并补偿退化过程造成的失真，以便获得未经退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。

典型的图像复原方法是根据图像退化的先验知识建立退化模型，以退化模型为基础采用滤波等手段进行处理，使复原后的图像符合一定的准则，达到改善图像质量的目的。

因此，图像复原是沿着质量降低的逆过程来重现真实的原始图像，通过去模糊函数而去除图像模糊。

6.3 约束最小二乘方滤波（Constrained Least Squares Filtering）

维纳滤波建立在退化函数和信噪比已知的前提下，这在实践中并不容易满足。

约束最小二乘方滤波仅要求噪声方差和均值的知识或估计，这些参数通常可以由一幅给定的退化图像算出，因而具有更为广泛的应用。而且，维纳滤波是以最小化一个统计准则为基础，因此是平均意义上的最优，而约束最小二乘方滤波则对每幅图像都会产生最优估计。

与维纳滤波相比，最小二乘方滤波对于高噪声和中等噪声处理效果更好，对于低噪声二者处理结果基本相同。当手工选择参数以取得更好的视觉效果时，约束最小二乘方滤波的效果有可能比维纳滤波的效果更好。

约束最小二乘方滤波的核心是解决退化函数对噪声的敏感性问题。降低噪声敏感的一种方法是，以平滑度量的最佳复原为基础，如图像的二阶导数（拉普拉斯算子），其数学描述是求一个带约束条件的准则函数的最小值：
$minC=∑x=0M−1∑y=0N−1[▽2f(x,y)]2s.t.:∣∣g−Hf^∣∣2=∣∣η∣∣2min \ C = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} [\triangledown ^2 f(x,y)]^2 \\ s.t.: ||g - H \hat{f}||^2 = ||\eta||^2$
其中， $a||^2$ 是欧几里德范数， $f^\hat{f}$ 是未退化图像的估计， $▽2\triangledown ^2$ 是拉普拉斯算子。求准则函数 C 的最小值，便得到最好的平滑效果即退化图像的最佳复原。

该约束最小二乘方问题的解是：
$F^(u,v)=[H∗(u,v)∣H(u,v)∣2+γ∣P(u,v)∣2]G(u,v)\hat{F}(u,v) = [\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2}] \ G(u,v)$
式中， $γ\gamma$ 是一个必须满足约束条件的参数， $P (u, v)$ 是函数 $p (x, y)$ 的傅里叶变换：
$\begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1\\ 0 & -1 &0 \end{bmatrix}$
注意函数 $P (u, v)$ 和 $H (u, v)$ 的大小必须相等。如果 $H (u, v)$ 的大小为 $M * N$ ，则 p(x,y) 必须嵌入 $M * N$ 零阵列的中心，并保持偶对称。

例程 9.21: 约束最小二乘方滤波

    # 9.21: 约束最小二乘方滤波def getMotionDsf(shape, angle, dist):xCenter = (shape[0] - 1) / 2yCenter = (shape[1] - 1) / 2sinVal = np.sin(angle * np.pi / 180)cosVal = np.cos(angle * np.pi / 180)PSF = np.zeros(shape)  # 点扩散函数for i in range(dist):  # 将对应角度上motion_dis个点置成1xOffset = round(sinVal * i)yOffset = round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter + yOffset)] = 1return PSF / PSF.sum()  # 归一化def makeBlurred(image, PSF, eps):  # 对图片进行运动模糊fftImg = np.fft.fft2(image)  # 进行二维数组的傅里叶变换fftPSF = np.fft.fft2(PSF) + epsfftBlur = np.fft.ifft2(fftImg * fftPSF)fftBlur = np.abs(np.fft.fftshift(fftBlur))return fftBlurdef wienerFilter(input, PSF, eps, K=0.01):  # 维纳滤波，K=0.01fftImg = np.fft.fft2(input)fftPSF = np.fft.fft2(PSF) + epsfftWiener = np.conj(fftPSF) / (np.abs(fftPSF)**2 + K)imgWienerFilter = np.fft.ifft2(fftImg * fftWiener)imgWienerFilter = np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilterdef getPuv(image):h, w = image.shape[:2]hPad, wPad = h - 3, w - 3pxy = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]])pxyPad = np.pad(pxy, ((hPad//2, hPad - hPad//2), (wPad//2, wPad - wPad//2)), mode='constant')fftPuv = np.fft.fft2(pxyPad)return fftPuvdef leastSquareFilter(image, PSF, eps, gamma=0.01):  # 约束最小二乘方滤波fftImg = np.fft.fft2(image)fftPSF = np.fft.fft2(PSF)conj = fftPSF.conj()fftPuv = getPuv(image)# absConj = np.abs(fftPSF) ** 2Huv = conj / (np.abs(fftPSF)**2 + gamma * (np.abs(fftPuv)**2))ifftImg = np.fft.ifft2(fftImg * Huv)ifftShift = np.abs(np.fft.fftshift(ifftImg))imgLSFilter = np.uint8(cv2.normalize(np.abs(ifftShift), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgLSFilter# # 读取原始图像img = cv2.imread("../images/Fig0526a.tif", 0)  # flags=0 读取为灰度图像hImg, wImg = img.shape[:2]# 带有噪声的运动模糊PSF = getMotionDsf((hImg, wImg), 45, 100)  # 运动模糊函数imgBlurred = np.abs(makeBlurred(img, PSF, 1e-6))  # 生成不含噪声的运动模糊图像scale = 0.01  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 对含有噪声的模糊图像进行维纳滤波imgLSFilter = leastSquareFilter(imgBlurNoisy, PSF, scale, gamma=0.01)  # 约束最小二乘方滤波plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("blurred image (dev=0.01)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(232), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(233), plt.title("least square filter"), plt.axis('off'), plt.imshow(imgLSFilter, 'gray')scale = 0.1  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 维纳滤波imgLSFilter = leastSquareFilter(imgBlurNoisy, PSF, scale, gamma=0.1)  # 约束最小二乘方滤波plt.subplot(234), plt.title("blurred image (dev=0.1)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(235), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(236), plt.title("least square filter"), plt.axis('off'), plt.imshow(imgLSFilter, 'gray')plt.tight_layout()plt.show()