import cv2
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.stats
%matplotlib inline
讀入我們需要的圖像
apple = cv2.imread("apple.jpg")
apple = cv2.resize(cv2.cvtColor(apple,cv2.COLOR_BGR2RGB),(200,200))
plt.imshow(apple)
plt.axis("off")
plt.show()
噪聲
高斯噪聲
簡介
高斯噪聲是指它的概率密度函數服從高斯分布（即正態分布）的一類噪聲
與椒鹽噪聲相似（Salt And Pepper Noise），高斯噪聲（gauss noise）也是數字圖像的一個常見噪聲。
椒鹽噪聲是出現在隨機位置、噪點深度基本固定的噪聲，高斯噪聲與其相反，是幾乎每個點上都出現噪聲、噪點深度隨機的噪聲。
正如上面的簡介我們只要實現一個隨機矩陣，矩陣中值總體來說符合高斯分布，與原圖像想加，就可以實現高斯噪聲了，python中的random提供了產生高斯隨機數的方法，但是numpy提供了直接生成隨機高斯矩陣的方法。
我們這里使用numpy即可
gauss = np.random.normal(mean,sigma,(row,col,ch))
因此我們可以得出產生高斯噪聲的方式
def GaussieNoisy(image,sigma):
row,col,ch= image.shape
mean = 0
gauss = np.random.normal(mean,sigma,(row,col,ch))
gauss = gauss.reshape(row,col,ch)
noisy = image + gauss
return noisy.astype(np.uint8)
plt.imshow(GaussieNoisy(apple,25))
plt.show()
上圖為施加sigma為25的高斯噪聲的效果
椒鹽噪聲
相比高斯噪聲，椒鹽噪聲的概念非常簡單，即在圖像中隨機選點，使其為0或255
def spNoisy(image,s_vs_p = 0.5,amount = 0.004):
row,col,ch = image.shape
out = np.copy(image)
num_salt = np.ceil(amount * image.size * s_vs_p)
coords = [np.random.randint(0, i - 1, int(num_salt)) for i in image.shape]
out[coords] = 1
num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
out[coords] = 0
return out
plt.imshow(spNoisy(apple))
plt.show()
濾波
算術均值濾波
算術均值濾波器即求某一范圍內圖像的均值，代替范圍中心點的值，在前面已經實現過。
def ArithmeticMeanOperator(roi):
return np.mean(roi)
def ArithmeticMeanAlogrithm(image):
new_image = np.zeros(image.shape)
image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
for i in range(1,image.shape[0]-1):
for j in range(1,image.shape[1]-1):
new_image[i-1,j-1] = ArithmeticMeanOperator(image[i-1:i+2,j-1:j+2])
new_image = (new_image-np.min(image))*(255/np.max(image))
return new_image.astype(np.uint8)
def rgbArithmeticMean(image):
r,g,b = cv2.split(image)
r = ArithmeticMeanAlogrithm(r)
g = ArithmeticMeanAlogrithm(g)
b = ArithmeticMeanAlogrithm(b)
return cv2.merge([r,g,b])
plt.imshow(rgbArithmeticMean(apple))
plt.show()
幾何均值濾波
幾何均值公式如下
\[f(x,y) = [\prod_{(s,t)\in S_{x,y}}{g(s,t)}]^{\frac 1{mn}} \]
def GeometricMeanOperator(roi):
roi = roi.astype(np.float64)
p = np.prod(roi)
return p**(1/(roi.shape[0]*roi.shape[1]))
def GeometricMeanAlogrithm(image):
new_image = np.zeros(image.shape)
image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
for i in range(1,image.shape[0]-1):
for j in range(1,image.shape[1]-1):
new_image[i-1,j-1] = GeometricMeanOperator(image[i-1:i+2,j-1:j+2])
new_image = (new_image-np.min(image))*(255/np.max(image))
return new_image.astype(np.uint8)
def rgbGemotriccMean(image):
r,g,b = cv2.split(image)
r = GeometricMeanAlogrithm(r)
g = GeometricMeanAlogrithm(g)
b = GeometricMeanAlogrithm(b)
return cv2.merge([r,g,b])
plt.imshow(rgbGemotriccMean(apple))
plt.show()
諧波均值
諧波均值公式定義如下
\[ H = \frac{n} {\frac{1}{x_1}+\frac{1}{x_2}+\frac{1}{x_3}\ldots \frac{1}{x_n}} \]
這里需要注意的是，諧波均值處理的數必須大於0,當x存在為0的數是，趨近於無窮，則H=0
因此我們此處當存在x大於0的數時，就返回0
def HMeanOperator(roi):
roi = roi.astype(np.float64)
if 0 in roi:
roi = 0
else:
roi = scipy.stats.hmean(roi.reshape(-1))
return roi
def HMeanAlogrithm(image):
new_image = np.zeros(image.shape)
image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
for i in range(1,image.shape[0]-1):
for j in range(1,image.shape[1]-1):
new_image[i-1,j-1] =HMeanOperator(image[i-1:i+2,j-1:j+2])
new_image = (new_image-np.min(image))*(255/np.max(image))
return new_image.astype(np.uint8)
def rgbHMean(image):
r,g,b = cv2.split(image)
r = HMeanAlogrithm(r)
g = HMeanAlogrithm(g)
b = HMeanAlogrithm(b)
return cv2.merge([r,g,b])
plt.imshow(rgbHMean(apple))
plt.show()
逆諧波均值
公式如下
\[f(x,y) = \frac{\sum_{(s,t)\in S_{xy}}{g(s,t)^{Q+1}}} {\sum_{(s,t)\in S_{xy}}{g(s,t)^{Q}}} \]
因此使用python實現如下
def IHMeanOperator(roi,q):
roi = roi.astype(np.float64)
return np.mean((roi)**(q+1))/np.mean((roi)**(q))
def IHMeanAlogrithm(image,q):
new_image = np.zeros(image.shape)
image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
for i in range(1,image.shape[0]-1):
for j in range(1,image.shape[1]-1):
new_image[i-1,j-1] = IHMeanOperator(image[i-1:i+2,j-1:j+2],q)
new_image = (new_image-np.min(image))*(255/np.max(image))
return new_image.astype(np.uint8)
def rgbIHMean(image,q):
r,g,b = cv2.split(image)
r = IHMeanAlogrithm(r,q)
g = IHMeanAlogrithm(g,q)
b = IHMeanAlogrithm(b,q)
return cv2.merge([r,g,b])
plt.imshow(rgbIHMean(apple,2))
plt.show()
圖像的復原
下面我們將試着對加了高斯噪聲和椒鹽噪聲的圖像進行復原
spApple = spNoisy(apple,0.5,0.1)
gaussApple = GaussieNoisy(apple,25)
plt.subplot(121)
plt.title("Salt And peper Image")
plt.imshow(spApple)
plt.axis("off")
plt.subplot(122)
plt.imshow(gaussApple)
plt.axis("off")
plt.title("Gauss noise Image")
plt.show()
arith_sp_apple = rgbArithmeticMean(spApple)
gemo_sp_apple = rgbGemotriccMean(spApple)
plt.subplot(121)
plt.title("Arithmatic to spImage")
plt.imshow(arith_sp_apple)
plt.axis("off")
plt.subplot(122)
plt.imshow(gemo_sp_apple)
plt.axis("off")
plt.title("Geomotric to spImage")
plt.show()
arith_gs_apple = rgbArithmeticMean(gaussApple)
gemo_gs_apple = rgbGemotriccMean(gaussApple)
plt.subplot(121)
plt.title("Arithmatic to gsImage")
plt.imshow(arith_gs_apple)
plt.axis("off")
plt.subplot(122)
plt.imshow(gemo_gs_apple)
plt.axis("off")
plt.title("Geomotric to gsImage")
plt.show()
算術均值能略微去除椒鹽噪聲產生的點，幾何均值效果卻有些奇怪。
對於高斯噪聲，二者的效果都非常弱
arith_sp_apple = rgbHMean(spApple)
gemo_sp_apple = rgbIHMean(spApple,3)
plt.subplot(121)
plt.title("H Mean to spImage")
plt.imshow(arith_sp_apple)
plt.axis("off")
plt.subplot(122)
plt.imshow(gemo_sp_apple)
plt.axis("off")
plt.title("IH mean to spImage")
plt.show()
arith_gs_apple = rgbHMean(gaussApple)
gemo_gs_apple = rgbIHMean(gaussApple,3)
plt.subplot(121)
plt.title("HMean to gsImage")
plt.imshow(arith_gs_apple)
plt.axis("off")
plt.subplot(122)
plt.imshow(gemo_gs_apple)
plt.axis("off")
plt.title("IHMean to gsImage")
plt.show()
如圖，IHMEAN的效果要比Hmean好很多,即使是高斯造神也能達到良好的去噪效果