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

标题

约束最小二乘方滤波
几何均值滤波

约束最小二乘方滤波

F^(u,v)=[H∗(u,v)∣H(u,v)∣2+γ∣P(u,v)∣2]G(u,v)(5.89)\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2} \bigg]G(u,v) \tag{5.89}F^(u,v)=[∣H(u,v)∣2+γ∣P(u,v)∣2H∗(u,v)]G(u,v)(5.89)

P(u,v)为函数p(x,y)=[0−10−14−10−10]P(u,v) 为函数 p(x, y) = \begin {bmatrix}0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end {bmatrix}P(u,v)为函数p(x,y)=⎣⎡0−10−14−10−10⎦⎤的傅里叶变换

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

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

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 - 3p_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 imageparam: PSF: input the PSF maskparam: eps: epsilonparam: gamma: gamma value for least square filter fuctionreturn 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).realhuv = 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()

几何均值滤波

F^(u,v)=[H∗(u,v)∣H(u,v)∣2]α[H∗(u,v)∣H(u,v)∣2+β[Sη(u,v)Sf(u,v)]]1−αG(u,v)(5.99)\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) \tag{5.99}F^(u,v)=[∣H(u,v)∣2H∗(u,v)]α[∣H(u,v)∣2+β[Sf(u,v)Sη(u,v)]H∗(u,v)]1−αG(u,v)(5.99)

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

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

β=1\beta=1β=1而α\alphaα大于1/2时，滤波器更的性能更像逆滤波器。
α\alphaα 小于1/2时，滤波器的性能更像维纳滤波器。
α\alphaα 等于1/2且β=1\beta=1β=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 imageparam: PSF  : input the PSF maskparam: eps  : epsilonparam: K    : K equal to math: \frac{S_{\eta}(u,v)}{S_f(u,v)}param: alpha: alpha value for filter fuctionparam: beta : beta value for the filter fuctionreturn 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).realhuv = 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()