python数字图像处理（五）图像的退化和复原

时间：2017-12-16 21:14:47 阅读：608 评论：0 收藏：0 [点我收藏+]

标签：icc return 生成 mat pos lin 必须 block pap

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

技术分享图片