第5章 Python 数字图像处理(DIP) - 图像复原与重建16 - 约束最小二乘方滤波、几何均值滤波

约束最小二乘方滤波

$$\hat{F}(u,v)=\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\gamma |P(u,v){|}^2}\right]G(u,v)\text{\hspace{1em}(5.89)}$$

$P(u,v)为函数p(x,y)=\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]$的傅里叶变换

我们认为这个函数是一个拉普拉斯核，注意式中在$\gamma =0$时会简化为逆滤波。

函数$P(u,v)和H(u,v)$的大小必须相同。意味着$p(x,y)$必须嵌入$M$ x $N$ 零阵列的中心。为保持$p(x,y)$的偶对称， $M$ 和 $N$必须是偶整数。如果由其获得$H$的一幅已知退化图像不是偶数维的，则在计算$H$之前需要酌情删除一行和/或一行，以便可以使用

2021-21-16 update, Thanks to fans id is ‘weixin_42674476’。

def get_puv(image):
    h, w = image.shape[:2]
    
    h_pad, w_pad = h - 3, w - 3
    p_xy = np.array([[0, -1, 0],
                 [-1, 4, -1],
                 [0, -1, 0]])
    p_pad = np.pad(p_xy, ((h_pad//2, h_pad - h_pad//2), (w_pad//2, w_pad - w_pad//2)), mode='constant')
    p_uv = np.fft.fft2(p_pad)
    return p_uv

def least_square_filter(image, PSF, eps, gamma=0.01):
    """
    least square filter for image denoise, math: 
    $$\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2} \bigg]G(u,v)$$
    param: image: input image
    param: PSF: input the PSF mask
    param: eps: epsilon
    param: gamma: gamma value for least square filter fuction
    return image after least square filter
    问题：为什么改变gamma值结果也没有变化，sorted 2021-12-16
    """
    fft = np.fft.fft2(image)
    PSF_fft = np.fft.fft2(PSF)
    conj = PSF_fft.conj()
    p_uv = get_puv(image)
    abs_conj = np.abs(PSF_fft) ** 2
    # abs_conj = (PSF_fft * conj).real
    huv = conj / (abs_conj + gamma * (np.abs(p_uv) ** 2))
    result = np.fft.ifft2(fft * huv)
    result = np.abs(np.fft.fftshift(result))
    return result, abs_conj

# 约束最小二乘方滤波
image = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif", 0)
PSF = get_motion_dsf(image.shape[:2], -50, 100)
blurred = make_blurred(image, PSF, 1e-5)

wiener = wiener_filter(blurred, PSF, 1e-5, K=0.03)     
least_square, abs_con = least_square_filter(blurred, PSF, 1e-5, gamma=1e-6)
least_square = np.uint8(normalize(least_square) * 255)

img_diff = image - least_square
plt.figure(figsize=(15, 5))
plt.subplot(131), plt.imshow(blurred, 'gray'), plt.title("Motion blurred"), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(wiener, 'gray'), plt.title("Wiener Filter"), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(least_square, 'gray'), plt.title("Least Square Filter"), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

几何均值滤波

$$\hat{F}(u,v)=\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2}{]}^{\alpha }\right[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\beta [\frac{S_{\eta }(u,v)}{S_f(u,v)}]}{]}^{1-\alpha }G(u,v)\text{\hspace{1em}(5.99)}$$

$\alpha 和\beta$ 是非负的实常数。

当$\alpha =1$时，几何均值滤波器简化为逆滤波器；当$\alpha =0$时，几何均值滤波器变为参数维纳滤波器，参数维纳滤波器在$\beta =1$时简化为标准维纳滤波器；当$\alpha =1/2$时几何均值滤波器变成幂次相同的两个量的乘积，这就是几何均值的定义。

$\beta =1$而$\alpha$大于1/2时，滤波器更的性能更像逆滤波器。
$\alpha$ 小于1/2时，滤波器的性能更像维纳滤波器。
$\alpha$ 等于1/2且$\beta =1$时，滤波器通常称为频谱均衡滤波器。

def geometric_mean_filter(image, PSF, eps, K=1, alpha=1, beta=1):
    """
    geometric mean filter for image denoise, math: 
    $$\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2}\bigg]^\alpha 
    \Bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \beta \big[\frac{S_{\eta}(u,v)}{S_f(u,v)} \big]}\Bigg]^{1-\alpha} G(u, v) $$
    param: image: input image
    param: PSF  : input the PSF mask
    param: eps  : epsilon
    param: K    : K equal to math: \frac{S_{\eta}(u,v)}{S_f(u,v)}
    param: alpha: alpha value for filter fuction
    param: beta : beta value for the filter fuction
    return image after least square filter
    """
    fft = np.fft.fft2(image)
    PSF_fft = np.fft.fft2(PSF)
    conj = PSF_fft.conj()
    abs_square = (PSF_fft * conj).real
    huv = np.power(conj / (abs_square), alpha) * np.power(conj / (abs_square + beta*(K)), 1 - alpha)
    result = np.fft.ifft2(fft * huv)
    result = np.abs(np.fft.fftshift(result))
    return result

#### 几何均值滤波器
image = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif", 0)
PSF = get_motion_dsf(image.shape[:2], -50, 100)
blurred = make_blurred(image, PSF, 1e-5)

wiener = wiener_filter(blurred, PSF, 1e-5, K=0.03)     
geometric_mean = geometric_mean_filter(blurred, PSF, 1e-5, K=1, alpha=1/2, beta=0)
geometric_mean = np.uint8(normalize(geometric_mean) * 255)

plt.figure(figsize=(14, 5))
plt.subplot(131), plt.imshow(blurred, 'gray'), plt.title("Motion blurred"), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(wiener, 'gray'), plt.title("Wiener Filter"), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(geometric_mean, 'gray'), plt.title("Geometric mean Filter"), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()