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

专栏地址:『youcans 的 OpenCV 例程 300篇 - 总目录』

【第 7 章:图像复原与重建】
106. 退化图像的逆滤波
107. 退化图像的维纳滤波
108. 约束最小二乘方滤波
109. 几何均值滤波


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


6. 退化图像复原

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

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

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


6.3 约束最小二乘方滤波(Constrained Least Squares Filtering)

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

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

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

约束最小二乘方滤波的核心是解决退化函数对噪声的敏感性问题。降低噪声敏感的一种方法是,以平滑度量的最佳复原为基础,如图像的二阶导数(拉普拉斯算子),其数学描述是求一个带约束条件的准则函数的最小值:
m i n   C = ∑ x = 0 M − 1 ∑ y = 0 N − 1 [ ▽ 2 f ( x , y ) ] 2 s . t . : ∣ ∣ g − H f ^ ∣ ∣ 2 = ∣ ∣ η ∣ ∣ 2 min \ 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 min C=x=0M1y=0N1[2f(x,y)]2s.t.:∣∣gHf^2=∣∣η2
其中, ∣ ∣ a ∣ ∣ 2 ||a||^2 ∣∣a2是欧几里德范数, f ^ \hat{f} f^ 是未退化图像的估计, ▽ 2 \triangledown ^2 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) F^(u,v)=[H(u,v)2+γP(u,v)2H(u,v)] G(u,v)
式中, γ \gamma γ 是一个必须满足约束条件的参数, P ( u , v ) P(u,v) P(u,v) 是函数 p ( x , y ) p(x,y) p(x,y) 的傅里叶变换:
p ( x , y ) = [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] p(x,y) = \begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1\\ 0 & -1 &0 \end{bmatrix} p(x,y)= 010141010
注意函数 P ( u , v ) P(u,v) P(u,v) H ( u , v ) H(u,v) H(u,v) 的大小必须相等。如果 H ( u , v ) H(u,v) H(u,v) 的大小为 M ∗ N M*N MN,则 p(x,y) 必须嵌入 M ∗ N M*N MN 零阵列的中心,并保持偶对称。


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

    # 9.21: 约束最小二乘方滤波
    def getMotionDsf(shape, angle, dist):
        xCenter = (shape[0] - 1) / 2
        yCenter = (shape[1] - 1) / 2
        sinVal = np.sin(angle * np.pi / 180)
        cosVal = np.cos(angle * np.pi / 180)
        PSF = np.zeros(shape)  # 点扩散函数
        for i in range(dist):  # 将对应角度上motion_dis个点置成1
            xOffset = round(sinVal * i)
            yOffset = round(cosVal * i)
            PSF[int(xCenter - xOffset), int(yCenter + yOffset)] = 1
        return PSF / PSF.sum()  # 归一化

    def makeBlurred(image, PSF, eps):  # 对图片进行运动模糊
        fftImg = np.fft.fft2(image)  # 进行二维数组的傅里叶变换
        fftPSF = np.fft.fft2(PSF) + eps
        fftBlur = np.fft.ifft2(fftImg * fftPSF)
        fftBlur = np.abs(np.fft.fftshift(fftBlur))
        return fftBlur

    def wienerFilter(input, PSF, eps, K=0.01):  # 维纳滤波,K=0.01
        fftImg = np.fft.fft2(input)
        fftPSF = np.fft.fft2(PSF) + eps
        fftWiener = np.conj(fftPSF) / (np.abs(fftPSF)**2 + K)
        imgWienerFilter = np.fft.ifft2(fftImg * fftWiener)
        imgWienerFilter = np.abs(np.fft.fftshift(imgWienerFilter))
        return imgWienerFilter

    def getPuv(image):
        h, w = image.shape[:2]
        hPad, wPad = h - 3, w - 3
        pxy = 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 fftPuv

    def 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) ** 2
        Huv = 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()

在这里插入图片描述

程序说明:

对于含有不同水平噪声的运动模糊图像,最小二乘方滤波的图像复原效果都较好,特别是对于高噪声和中等噪声的处理效果更好。

对于最小二乘方滤波,可以交互地或循环地调整参数 γ \gamma γ 以取得更好的复原效果,例如参数 γ \gamma γ 可以从较小的取值递增以获得最优估计。


(本节完)


版权声明:
youcans@xupt 原创作品,转载必须标注原文链接:(https://blog.csdn.net/youcans/article/details/123062903)
Copyright 2022 youcans, XUPT
Crated:2022-2-22


  • 5
    点赞
  • 34
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 4
    评论
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

youcans_

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值