CV-图像降噪-基于高斯的掩模

# -*- coding:utf-8 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt
global gauss_map


s1 = cv2.imread("original-image.jpg", 0)
s2 = cv2.imread("noisy-image.jpg", 0)
subimg = s1-s2

plt.hist(subimg.ravel(), 256, [0, 256])  # 画灰度分布图
plt.show()


# take a Fast Fourier Transform of subimg
dft = cv2.dft(np.float32(subimg), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

# convert the result into log-magnitude image
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))

# plot the input image and its log-magnitude image
plt.subplot(121), plt.imshow(subimg, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])

plt.subplot(122), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

plt.imshow(magnitude_spectrum, cmap='gray')
plt.show()

# design the gaussion based filter
img = cv2.imread('magnitude_spectrum.png')
sz = img.shape

# the matrix "gauss_map" is used as a mask
gauss_map = np.zeros((sz[0], sz[1]))


def gauss(center_y, center_x, sigma):
    sum = 0
    pi = 3.1415926
    for i in range(center_x - 5, center_x + 6):
        for j in range(center_y - 5, center_y + 6):
            x2 = pow(i-center_x, 2)
            y2 = pow(j-center_y, 2)
            g = np.exp(-(x2+y2)/(2*sigma*sigma))
            g /= 2*pi*sigma*sigma
            sum += g
            gauss_map[i][j] = g
    for i in range(center_x - 5, center_x + 6):
        for j in range(center_y - 5, center_y + 6):
            gauss_map[i][j] = gauss_map[i][j]/sum*30


# aplly the gauss_map on the four peaks of the
# log-magnitude image
gauss(43, 52, 2.8)
gauss(106, 40, 2.8)
gauss(22, 88, 2.8)
gauss(85, 76, 2.8)

plt.figure()
plt.imshow(gauss_map, plt.cm.gray)
plt.show()

# get the filter as 1 - (G1 + G2 + ... + Gn)
gauss_map = 1-gauss_map


# plot the filter
plt.figure()
plt.imshow(gauss_map, plt.cm.gray)
plt.show()


# take a Fast Fourier Transform of noisy-image
f = np.fft.fft2(s2)
fshift = np.fft.fftshift(f)


# multiply the filter with the FFT of the noisy image
fshift = fshift*gauss_map  # 两个矩阵对应元素相乘

# using an inverse Fast Fourier Transform of fshift
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

# plot the final result
plt.subplot(121), plt.imshow(s2, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img_back, cmap='gray')
plt.title('Image after removing noise'), plt.xticks([]), plt.yticks([])
plt.show()

diff = s1 - img_back
cv2.imshow("diff", diff)
cv2.waitKey(0)

dft = cv2.dft(np.float32(diff), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
ms = 20*np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
plt.imshow(ms, cmap='gray')
plt.show()