约束最小二乘方滤波

$$\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 - 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()

几何均值滤波

$$\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 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()