約束最小二乘方濾波

$$\hat{F}(u,v)=\left[\frac{H^{?}(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)=?\begin{array}{ccc}0& ?1& 0\\ ?1& 4& ?1\\ 0& ?1& 0\end{array}?$的傅里葉變換

我們認為這個函數是一個拉普拉斯核，注意式中在$\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^{?}(u,v)}{|H(u,v){|}^2}{]}^{\alpha }\right[\frac{H^{?}(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()