加噪去噪初体验

本文通过python实现图像的加噪去噪:
具体代码如下(参考:https://blog.csdn.net/XITMan/article/details/105290126,含详细注释):

#import os              #import语句的作用是用来导入模块,可以出现在程序任何位置
import cv2 as cv  # 导入openCV库
import skimage  # 导入skimage模块.scikit-image是一个图像处理算法的集合。它是基于scipy的一款图像处理包,它将图片作为numpy数组进行处理,方便进行后续运算。
# 必须首先安装numpy,scipy,matplotlib
import numpy as np  # 导入numpy模块。numpy是python扩展程序库,支持数组和矩阵运算,针对数组运算提供大量数学函数库。


def boxBlur(img):
    # 使用5x5的滤波核进行平滑
    blur = cv.boxFilter(img, -1, (5, 5))
    return blur


def gaussianBlur(img):
    #     使用高斯核进行平滑
    blur = cv.GaussianBlur(img, (5, 5), 1.5)
    return blur


def main():
    # 2. 定义图片类img
    path = r"C:\Users\Administrator\Desktop\lds.jpg"
    img = cv.imread(path)
    start_t = cv.getTickCount()
    # 5. 加噪声,绘图
    ##############################################3
    # add gaussian noise

    gauss_noiseImg = skimage.util.random_noise(img, mode='gaussian')  # 添加10%的高斯噪声
    gauss_noiseImg = gauss_noiseImg
    salt_noiseImg = skimage.util.random_noise(img, mode='salt')  # 添加椒盐噪声

    lb_gauss1 = cv.medianBlur(gauss_noiseImg.astype('float32'), 1)  # 中值滤波

    lb_salt1 = cv.medianBlur(salt_noiseImg.astype('float32'), 1)  # 中值滤波
    lb_gauss3 = cv.medianBlur(gauss_noiseImg.astype('float32'), 3)
    lb_salt3 = cv.medianBlur(salt_noiseImg.astype('float32'), 3)
    lb_gauss5 = cv.medianBlur(gauss_noiseImg.astype('float32'), 5)
    lb_salt5 = cv.medianBlur(salt_noiseImg.astype('float32'), 5)
    print(gauss_noiseImg.dtype, "gaussian noisy image dtype")  # 输出一个注释
    print(gauss_noiseImg.shape, "gaussian noisy image shape")  # 输出一个注释

    print(salt_noiseImg.dtype, "salt noisy image dtype")  # 输出一个注释
    print(salt_noiseImg.shape, "salt noisy image shape")  # 输出一个注释

    cv.namedWindow("Original Image", cv.WINDOW_NORMAL)  # 输出原图片的标题
    cv.imshow('Original Image', img)  # 输出原图片

    # Gaussian noisy image
    cv.namedWindow("Added Gaussian Noise Image", cv.WINDOW_NORMAL)  # 输出高斯噪声图片的标题
    cv.imshow('Added Gaussian Noise Image', gauss_noiseImg)  # 输出高斯噪声图片

    # Salt noisy image
    cv.namedWindow("Added Salt Noise Image", cv.WINDOW_NORMAL)  # 输出椒盐噪声图片的标题
    cv.imshow('Added Salt Noise Image', salt_noiseImg)  # 输出椒盐噪声图片

    # 滤波后的图像
    cv.namedWindow("lbguass Image1", cv.WINDOW_NORMAL)  # 输出滤波后高斯噪声图片标题
    cv.imshow('lbguass Image1', lb_gauss1)  # 输出滤波后高斯噪声图片
    cv.namedWindow("lbsalt Image1", cv.WINDOW_NORMAL)  # 输出滤波后椒盐噪声图片标题
    cv.imshow('lbsalt Image1', lb_salt1)  # 输出滤波后椒盐噪声图片
    cv.namedWindow("lbguass Image3", cv.WINDOW_NORMAL)
    cv.imshow('lbguass Image3', lb_gauss3)
    cv.namedWindow("lbsalt Image3", cv.WINDOW_NORMAL)
    cv.imshow('lbsalt Image3', lb_salt3)
    cv.namedWindow("lbguass Image5", cv.WINDOW_NORMAL)
    cv.imshow('lbguass Image5', lb_gauss5)
    cv.namedWindow("lbsalt Image5", cv.WINDOW_NORMAL)
    cv.imshow('lbsalt Image5', lb_salt5)

    #####################################################

    stop_t = ((cv.getTickCount() - start_t) / cv.getTickFrequency()) * 1000  # 运行时间

    print(stop_t, "ms")  # 输出时间并加上单位

    cv.waitKey(0)
    cv.destroyAllWindows()


if __name__ == "__main__":
    main()

输出结果:
在这里插入图片描述
对比得出:经过中值滤波5的滤波后噪声比中值滤波3滤波后的图像噪声更少,但是图像也更模糊一些。
问题:在进行中值滤波7去噪时代码出错, 即如果要进行medianBlur,ksize=3,5就没有问题,只要ksize >=7 就开始报错。
可能原因:输入1-,3-,4-通道的图像; 当ksize 是 3 或 5时,图像的深度应该时 CV_8U, CV_16U, or CV_32F, 如果孔径更大, 它只可能是CV_8U.
后面需要用到OPENCV图形转换,后续将继续学习实现…
(接上)
经过不懈努力,实现了中值滤波7去噪,以下是改进的代码。

#import os              #import语句的作用是用来导入模块,可以出现在程序任何位置
import cv2 as cv  # 导入openCV库
import skimage  # 导入skimage模块.scikit-image是一个图像处理算法的集合。它是基于scipy的一款图像处理包,它将图片作为numpy数组进行处理,方便进行后续运算。
# 必须首先安装numpy,scipy,matplotlib
import numpy as np  # 导入numpy模块。numpy是python扩展程序库,支持数组和矩阵运算,针对数组运算提供大量数学函数库。


def boxBlur(img):
    
    blur = cv.boxFilter(img, -1, (7, 7))
    return blur


def gaussianBlur(img):
    #     使用高斯核进行平滑
    blur = cv.GaussianBlur(img, (7, 7), 1.5)
    return blur


def main():
    # 2. 定义图片类img
    path = r"C:\Users\Administrator\Desktop\lds.jpg"
    img = cv.imread(path)
    start_t = cv.getTickCount()
    # 5. 加噪声,绘图
    ##############################################3
    # add gaussian noise

    gauss_noiseImg = skimage.util.random_noise(img, mode='gaussian')  # 添加10%的高斯噪声
    gauss_noiseImg = gauss_noiseImg
    salt_noiseImg = skimage.util.random_noise(img, mode='salt')  # 添加椒盐噪声
    lb_gauss7 = cv.medianBlur(np.uint8(gauss_noiseImg * 255), 7)
    lb_salt7 = cv.medianBlur(np.uint8(salt_noiseImg * 255), 7)
    cv.namedWindow("Original Image", cv.WINDOW_NORMAL)  # 输出原图片的标题
    cv.imshow('Original Image', img)  # 输出原图片
    cv.namedWindow("Added Gaussian Noise Image", cv.WINDOW_NORMAL)  # 输出高斯噪声图片的标题
    cv.imshow('Added Gaussian Noise Image', gauss_noiseImg)  # 输出高斯噪声图片
    cv.namedWindow("Added Salt Noise Image", cv.WINDOW_NORMAL)  # 输出椒盐噪声图片的标题
    cv.imshow('Added Salt Noise Image', salt_noiseImg)  # 输出椒盐噪声图片
    cv.namedWindow("lbguass Image7", cv.WINDOW_NORMAL)
    cv.imshow('lbguass Image7', lb_gauss7)
    cv.namedWindow("lbsalt Image7", cv.WINDOW_NORMAL)
    cv.imshow('lbsalt Image7', lb_salt7)
    cv.waitKey(0)
    cv.destroyAllWindows()
if __name__ == "__main__":
        main()

下面是原图和加噪以及中值滤波7去噪的图像
在这里插入图片描述
不难看出此时去噪的效果更好,但图片较原图更加模糊。

下面进行了噪声检测(方法来源:https://www.cnblogs.com/ljy1227476113/p/12171379.html),分别通过FFT和MSE方法进行噪声检测(此处用的是高斯噪声图片,具体见前文)。代码及运行结果如下:

FFT

from PIL import Image
import numpy as np
import math

T = 50  # 阈值设定,大于T则判定偏离xy轴过多


# 复数类
class complex:
    def __init__(self):
        self.real = 0.0
        self.image = 0.0


# 复数乘法
def mul_ee(complex0, complex1):
    complex_ret = complex()
    complex_ret.real = complex0.real * complex1.real - complex0.image * complex1.image
    complex_ret.image = complex0.real * complex1.image + complex0.image * complex1.real
    return complex_ret


# 复数加法
def add_ee(complex0, complex1):
    complex_ret = complex()
    complex_ret.real = complex0.real + complex1.real
    complex_ret.image = complex0.image + complex1.image
    return complex_ret


# 复数减法
def sub_ee(complex0, complex1):
    complex_ret = complex()
    complex_ret.real = complex0.real - complex1.real
    complex_ret.image = complex0.image - complex1.image
    return complex_ret


# 对输入数据进行倒序排列
def forward_input_data(input_data, num):
    j = num // 2
    for i in range(1, num - 2):
        if (i < j):
            complex_tmp = input_data[i]
            input_data[i] = input_data[j]
            input_data[j] = complex_tmp
            # print "forward x[%d] <==> x[%d]" % (i, j)
        k = num // 2
        while (j >= k):
            j = j - k
            k = k // 2
        j = j + k


# 实现1D FFT
def fft_1d(in_data, num):
    PI = 3.1415926
    forward_input_data(in_data, num)  # 倒序输入数据

    # 计算蝶形级数,也就是迭代次数
    M = 1  # num = 2^m
    tmp = num // 2;
    while (tmp != 1):
        M = M + 1
        tmp = tmp // 2
    # print "FFT level:%d" % M

    complex_ret = complex()
    for L in range(1, M + 1):
        B = int(math.pow(2, L - 1))  # B为指数函数返回值,为float,需要转换integer
        for J in range(0, B):
            P = math.pow(2, M - L) * J
            for K in range(J, num, int(math.pow(2, L))):
                # print "L:%d B:%d, J:%d, K:%d, P:%f" % (L, B, J, K, P)
                complex_ret.real = math.cos((2 * PI / num) * P)
                complex_ret.image = -math.sin((2 * PI / num) * P)
                complex_mul = mul_ee(complex_ret, in_data[K + B])
                complex_add = add_ee(in_data[K], complex_mul)
                complex_sub = sub_ee(in_data[K], complex_mul)
                in_data[K] = complex_add
                in_data[K + B] = complex_sub
                # print "A[%d] real: %f, image: %f" % (K, in_data[K].real, in_data[K].image)
            # print "A[%d] real: %f, image: %f" % (K + B, in_data[K + B].real, in_data[K + B].image)


def test_fft_1d(in_data):
    # in_data = [2,3,4,5,7,9,10,11,100,12,14,11,56,12,67,12] #待测试的x点元素
    k = 1
    while (1):
        if len(in_data) > pow(2, k) and len(in_data) <= pow(2, k + 1):  # 不足的补0
            # fftlen=pow(2,k+1)
            # in_data.extend([0 for i in range(pow(2,k+1)-len(in_data))])
            fftlen = pow(2, k)
            break
        k += 1
    # 变量data为长度为x、元素为complex类实例的list,用于存储输入数据
    data = [(complex()) for i in range(len(in_data))]
    # 将8个测试点转换为complex类的形式,存储在变量data中
    for i in range(len(in_data)):
        data[i].real = in_data[i]
        data[i].image = 0.0

    ##输出FFT需要处理的数据
    # print("The input data:")
    # for i in range(len(in_data)):
    #    print("x[%d] real: %f, image: %f" % (i, data[i].real, data[i].image))

    fft_1d(data, fftlen)

    ##输出经过FFT处理后的结果
    # print("The output data:")
    # for i in range(len(in_data)):
    # print("X[%d] real: %f, image: %f" % (i, data[i].real, data[i].image))

    Tnum = 0
    for i in range(len(in_data)):  # 虚实值都大于T的才叫偏离
        if abs(data[i].real) > T and abs(data[i].image) > T:
            Tnum += 1
    print(Tnum)
    print(str(round(Tnum / len(in_data), 4) * 100) + "%")


# test the 1d fft
# in_data=[2,3,4,5,7,9,10,11]
demo = Image.open("D:\pythonproject\Added Gaussian Noise Image.jpg")
im = np.array(demo.convert('L'))  # 灰度化矩阵
in_data = []
for item in im:
    in_data.extend(item)
test_fft_1d(in_data)

FFT运行结果:

在这里插入图片描述

MSE方法:

import cv2
from PIL import Image
from PIL import ImageChops
import numpy as np
import time
import pytesseract
import warnings
import skimage
from scipy import stats
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
demo=Image.open("D:\pythonproject\Added Gaussian Noise Image.jpg")
im=np.array(demo.convert('L'))#灰度化矩阵
print(im.shape)
print(im.dtype)
#print(im)
height=im.shape[0]#尺寸
width=im.shape[1]
varlist=[]
for i in range(height):
    for j in range(width):
        for k in range(16):
            if im[i][j]>=k*16 and im[i][j]<(k+1)*16:#16级量化
                im[i][j]=8*(k*2+1)
                break
for i in range(0,height-height%3,3):
    for j in range(0,width-width%3,3):
        x=(im[i][j]+im[i+1][j]+im[i+2][j]+im[i][j+1]+im[i+1][j+1]+im[i+2][j+1]+im[i][j+2]+im[i+1][j+2]+im[i+2][j+2])/9
        x2=(pow(im[i][j],2)+pow(im[i+1][j],2)+pow(im[i+2][j],2)+pow(im[i][j+1],2)+pow(im[i+1][j+1],2)+pow(im[i+2][j+1],2)+pow(im[i][j+2],2)+pow(im[i+1][j+2],2)+pow(im[i+2][j+2],2))/9
        var=x2-pow(x,2)
        varlist.append(round(var,3))#子窗口的方差值3x3
print(im)
#print(varlist)
T=round(sum(varlist)/len(varlist),3)#保留3位小数
print(T)

MSE法运行结果:

在这里插入图片描述
可以看出FFT法比均方误差法要准确,但耗费时间更长。
现在遇到的问题是不知如何得到检测噪声的分布,接下来将攻克获取噪声分布以及绘制噪声分布图的难题,待续。
2021.11.10更新:
经过不懈努力,终于将噪声分布的图画了出来,代码如下:
MSE噪声检测图像:

import pickle

import cv2
from PIL import Image
from PIL import ImageChops
import numpy as np
import time
import pytesseract
import warnings
import skimage
from scipy import stats
import matplotlib.pyplot as plt
# warnings.filterwarnings("ignore")
# demo=Image.open("D:\pythonproject\Added Gaussian Noise Image.jpg")
# im=np.array(demo.convert('L'))#灰度化矩阵
# print(im.shape)
# print(im.dtype)
# #print(im)
# height=im.shape[0]#尺寸
# width=im.shape[1]
# varlist=[]
# for i in range(height):
#     for j in range(width):
#         for k in range(16):
#             if im[i][j]>=k*16 and im[i][j]<(k+1)*16:#16级量化
#                 im[i][j]=8*(k*2+1)
#                 break
# for i in range(0,height-height%3,3):
#     for j in range(0,width-width%3,3):
#         x=(im[i][j]+im[i+1][j]+im[i+2][j]+im[i][j+1]+im[i+1][j+1]+im[i+2][j+1]+im[i][j+2]+im[i+1][j+2]+im[i+2][j+2])/9
#         x2=(pow(im[i][j],2)+pow(im[i+1][j],2)+pow(im[i+2][j],2)+pow(im[i][j+1],2)+pow(im[i+1][j+1],2)+pow(im[i+2][j+1],2)+pow(im[i][j+2],2)+pow(im[i+1][j+2],2)+pow(im[i+2][j+2],2))/9
#         var=x2-pow(x,2)
#         varlist.append(round(var,3))#子窗口的方差值3x3
# with open('MSE.data','wb') as f:
#     pickle.dump(varlist,f)
with open('MSE.data','rb') as f:
    varlist = pickle.load(f)

#print(im)
#print(varlist)
T=round(sum(varlist)/len(varlist),3)#保留3位小数
print(T)
varlist = [i // 256 for i in varlist]
print(max(varlist))
print(min(varlist))
a=max(varlist)
b=min(varlist)
#赋值mu,sigma
mu = 0
sigma = 1
x = np.arange(b, a+1)
y = np.zeros_like(x)
for e in varlist:
        y[int(e-b)] = y[int(e-b)] + 1
print(y)
#绘图
plt.bar(x,y)
#x轴文本
plt.xlabel('随机变量:x')
#y轴文本
plt.ylabel('概率:y')
#标题
plt.title('正态分布:$\mu$=%.1f,$\sigma^2$=%.1f' % (mu,sigma))
#网格
plt.grid()
#显示图形
plt.show()

运行结果如下,这是加了高斯噪声的图片的噪声检测图,可以看出是一个近似高斯分布的样子。。。
在这里插入图片描述
然后我突发奇想检测了一下原图的噪声分布,也顺便画了一个图。
运行结果如下:
在这里插入图片描述在这里插入图片描述
由于相关知识欠缺,不知道如何分析这两个图,等后续更新!

  • 3
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值