只存在噪聲的復原 - 空間濾波

僅被加性噪聲退化
$$g(x,y)=f(x,y)+\eta (x,y)\text{\hspace{1em}(5.21)}$$
和
$$G(u,v)=F(u,v)+N(u,v)\text{\hspace{1em}(5.22)}$$

噪聲項通常是未知的，因此不能直接從退化圖像中減去噪聲項來得到復原的圖像。

對于周期噪聲，只是個例外而不是規律，我們可以用譜$G(u,v)$來估計$N(u,v)$，從$G(u,v)$中減去$N(u,v)$能夠得到原圖像的一個估計。

均值濾波器

算術平均濾波器

算術平均濾波器是最簡單的均值濾波器。
$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}g(r,c)\text{\hspace{1em}(5.23)}$$

$S_{xy}$表示中心為$(x,y)$、大小為$m\times n$的矩形子圖像窗口（鄰域）的一組坐標。算術平均濾波器在由$S_{xy}$定義的區域中，計算被污染圖像$g(x,y)$的平均值。復原的圖像$\hat{f}$在$(x,y)$處的值，是使用$S_{xy}$定義的鄰域中的像素算出的算術平均值。

均值濾波平滑圖像中的局部變化，它會降低圖像中的噪聲，但會模糊圖像

def arithmentic_mean(image, kernel):"""define arithmentic mean filter, math: $$\hat{f}(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} g(r,c)$$param image: input imageparam kerne: input kernel, actually use kernel shapereturn: image after arithmentic mean filter, """img_h = image.shape[0]img_w = image.shape[1]m, n = kernel.shape[:2]padding_h = int((m -1)/2)padding_w = int((n -1)/2)# 這樣的填充方式，可以奇數核或者偶數核都能正確填充image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")image_mean = image.copy()for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):temp = np.sum(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1])image_mean[i - padding_h][j - padding_w] = 1/(m * n) * tempreturn image_mean

# 算術平均濾波器
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接讀為灰度圖像mean_kernal = np.ones([3, 3])
mean_kernal = mean_kernal / mean_kernal.sizeimg_arithmentic = arithmentic_mean(img_ori, kernel=mean_kernal)img_cv2_mean = cv2.filter2D(img_ori, ddepth= -1, kernel=mean_kernal)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_arithmentic, 'gray'), plt.title('Self Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_cv2_mean, 'gray'), plt.title('CV2 mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

幾何均值濾波器

$$\hat{f}(x,y)=[\prod _{(r,c)\in S_{xy}}g(r,c){]}^{\frac{1}{mn}}\text{\hspace{1em}(5.24)}$$

幾何均值濾波器實現的平滑可與算術平均濾波器的相比，但損失的圖像細節更少

def geometric_mean(image, kernel):"""define geometric mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$param image:  input imageparam kerne:  input kernel, actually use kernel shapereturn: image after geometric mean filter, """img_h = image.shape[0]img_w = image.shape[1]m, n = kernel.shape[:2]order = 1 / (kernel.size)padding_h = int((m -1)/2)padding_w = int((n -1)/2)# 這樣的填充方式，可以奇數核或者偶數核都能正確填充image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")image_mean = image.copy()# 這里要指定數據類型，但指定是uint64或者float64，但結果不正確，反而乘以1.0，也是float64，但卻讓結果正確for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):prod = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0)image_mean[i - padding_h][j - padding_w] = np.power(prod, order)return image_mean

def geometric_mean(image, kernel):""":param image: input image:param kernel: input kernel:return: image after convolution"""img_h = image.shape[0]img_w = image.shape[1]kernel_h = kernel.shape[0]kernel_w = kernel.shape[1]# paddingpadding_h = int((kernel_h -1)/2)padding_w = int((kernel_w -1)/2)image_pad = np.pad(image.copy(), (padding_h, padding_w), mode="constant", constant_values=1)image_convol = image.copy()for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):temp = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * kernel)image_convol[i - padding_h][j - padding_w] = np.power(temp, 1/kernel.size)return image_convol

# 幾何均值濾波器
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接讀為灰度圖像mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.sizeimg_geometric = geometric_mean(img_ori, kernel=mean_kernal)img_arithmentic = arithmentic_mean(img_ori, kernel=arithmetic_kernel)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmetic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
# plt.show()

諧波平均濾波器

$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(r,c)\in S_{xy}}\frac{1}{g(r,c)}}\text{\hspace{1em}(5.25)}$$

可以處理鹽粒噪聲，又能處理類似于高斯噪聲的其他噪聲，但不能處理胡椒噪聲

def harmonic_mean(image, kernel):"""define harmonic mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$param image:  input imageparam kerne:  input kernel, actually use kernel shapereturn: image after harmonic mean filter, """epsilon = 1e-8img_h = image.shape[0]img_w = image.shape[1]m, n = kernel.shape[:2]order = kernel.sizepadding_h = int((m -1)/2)padding_w = int((n -1)/2)# 這樣的填充方式，可以奇數核或者偶數核都能正確填充image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")image_mean = image.copy()# 這里要指定數據類型，但指定是uint64或者float64，但結果不正確，反而乘以1.0，也是float64，但卻讓結果正確# 要加上epsilon，防止除0for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):temp = np.sum(1 / (image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0 + epsilon))image_mean[i - padding_h][j - padding_w] = order / tempreturn image_mean

# 諧波濾波處理高斯噪聲
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接讀為灰度圖像mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.sizeimg_geometric = geometric_mean(img_ori, kernel=mean_kernal)img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

# 諧波濾波處理鹽粒噪聲
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0) #直接讀為灰度圖像mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.sizeimg_geometric = geometric_mean(img_ori, kernel=mean_kernal)img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

# 諧波濾波處理胡椒噪聲
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接讀為灰度圖像mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.sizeimg_geometric = geometric_mean(img_ori, kernel=mean_kernal)img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

反（逆）諧波平均濾波器

$$\hat{f}(x,y)=\frac{{\sum }_{(r,c)\in S_{xy}}g(r,c{)}^{Q+1}}{{\sum }_{(r,c)\in S_{xy}}g(r,c{)}^Q}\text{\hspace{1em}(5.26)}$$

Q稱為濾波器的階數。這種濾波器適用于降低或消除椒鹽噪聲。Q值為正時，該濾波器消除胡椒噪聲；Q值為負時，該濾波器消除鹽粒噪聲。然而，該濾波器不能同時消除這兩種噪聲。注意當$Q=0$時，簡化為算術平均濾波器；當$Q=?1$時，簡化為諧波平均濾波器。

def inverse_harmonic_mean(image, kernel, Q=0):"""define inverse harmonic mean filter, math: $$\hat{f}(x, y) = \frac{\sum_{(r,c)\in S_{xy}} g(r,c)^{Q+1}}{\sum_{(r,c)\in S_{xy}} g(r,c)^Q}$$param image : input imageparam kernel: input kernel, actually use kernel shapeparam Q     : input order of the filter, default is 0, which equal to arithmentic mean filter, while is -1 is harmonic mean filterreturn: image after inverse harmonic mean filter, """epsilon = 1e-8img_h = image.shape[0]img_w = image.shape[1]m, n = kernel.shape[:2]padding_h = int((m -1)/2)padding_w = int((n -1)/2)# 這樣的填充方式，可以奇數核或者偶數核都能正確填充image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \(padding_w, n - 1 - padding_w)), mode="edge")image_mean = image.copy()# 這里要指定數據類型，但指定是uint64或者float64，但結果不正確，反而乘以1.0，也是float64，但卻讓結果正確# 要加上epsilon，防止除0for i in range(padding_h, img_h + padding_h):for j in range(padding_w, img_w + padding_w):temp = image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * 1.0 + epsilon# image_mean[i - padding_h][j - padding_w] = np.sum(temp**(Q+1)) / np.sum(temp**Q + epsilon)image_mean[i - padding_h][j - padding_w] = np.sum(np.power(temp, (Q+1))) / np.sum(np.power(temp, Q) + epsilon)return image_mean

# 反諧波濾波處理胡椒噪聲
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接讀為灰度圖像mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.sizeimg_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernel)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

# 反諧波濾波處理椒鹽噪聲
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0510(a)(ckt-board-saltpep-prob.pt05).tif', 0) #直接讀為灰度圖像mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.sizeimg_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)
img_arithmentic_mean = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=0)
# img_arithmentic_mean = arithmentic_mean(img_ori, kernel=mean_kernel)plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

下面是各種濾波器的對比與總結

# 算術平均濾波器和幾何均值濾波器對比
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接讀為灰度圖像
img_noise = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接讀為灰度圖像mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.sizeimg_arithmentic = arithmentic_mean(img_noise, kernel=mean_kernel)
img_geometric   = geometric_mean(img_noise, kernel=mean_kernel)plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_noise, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

反諧波濾波器可以處理胡椒與鹽粒噪聲。在處理胡椒噪聲的時候，Q值一般為正值；在處理鹽粒噪聲時，Q值一般為負值。錯誤的選擇Q值，不僅不能得到好的濾波效果，反而帶來災難性的結果。

# 反諧波濾波器分別使用不同的Q值對胡椒與鹽粒噪聲的濾波器
img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.sizeinverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=1.5)
inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=-1.5)plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

# 錯誤的使用Q值，反諧波濾波器得到的是嚴重的后果
img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.sizeinverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=-1.5)
inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=1.5)plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()