一、实验目的与要求
1．了解频域变换过程，掌握频域变换特点
2．熟练掌握频域滤波中常用的平滑和锐化滤波器，能够对不同要求的图像进行滤波处理，体会并正确评价滤波效果，了解不同滤波方式的使用场合，能够从理论上作出合理的解释。
二、实验相关知识
图像增强是指按特定的需要突出一幅图像中的某些有用信息，同时消弱或去除某些不需要的信息的处理方法，其主要目的是使处理后的图像对某些特定的应用比原来的图像更加有效。图像平滑与锐化处理是图像增强的主要研究内容。
和本实验有关的几个常用Matlab函数：
(1)imnoise：用于对图像生成模拟噪声，如：
j=imnoise(i,'gaussian',0,0.02)%在图像i上叠加均值为0、方差为0.02的高斯噪声，得到含噪图像j
j=imnoise(i,'salt&pepper',0.04) %在图像i上叠加密度为0.04的椒盐噪声，得到含噪图像j
(2)fspecial：用于产生预定义滤波器，如：
h=fspecial('average',3); %产生3×3模板的均值滤波器
h=fspecial('sobel');%产生sobel水平边缘增强的滤波器
可选项还有：'gaussian'高斯低通滤波器、'laplacian'拉普拉斯滤波器、'log'高斯拉普拉斯滤波器等
(3) imfilter、filter2、conv2：均是基于卷积的图像滤波函数，都可用于图像滤波，用法类似，如：
i=imread('p1.tif');
h=fspecial('prewitt');%产生prewitt算子的水平方向模板
j1=imfilter(i,h);%或者j2=filter2(h,i);或者j3=conv2(i,h);
(4) fft2:二维快速傅里叶变换函数
(5) fftshift:中心变换函数
(6) abs:取绝对值或复数取幅值
三、实验内容
1、图像频域平滑（去噪）：使用自生成图像（包含白色区域，黑色区域，并且部分区域添加椒盐噪声），然后进行傅里叶变换，并且分别使用理想低通滤波器、巴特沃斯低通滤波器、指数低通滤波器和梯形低通滤波器（至少使用两种低通滤波器），显示滤波前后的频域能量分布图，空间图像。分析不同滤波器对噪声、边缘的处理效果及其优缺点。
2、图像频域平滑（锐化）：选择一幅图像，例如rice.png，分别使用理想高通滤波器、巴特沃斯高通滤波器、指数高通滤波器和梯形高通滤波器（至少使用两种高通滤波器），显示滤波前后的频域能量分布图，空间图像。分析不同滤波器处理效果及其优缺点。
四、实验代码与结果
1、平滑
Exp4-1.py
import mathimport timeimport cv2import numpy as npfrom matplotlib import pyplot as pltfrom numba import cudafrom utils import Utilsimport argparseparser = argparse.ArgumentParser("")parser.add_argument("--method", default=1, type=int,choices=[1, 2, 3, 4])parser.add_argument("--d0", default=40, type=int)parser.add_argument("--d1", default=80, type=int)parser.add_argument("--n", default=2, type=int)args = parser.parse_args()def generate_img(height, width):img = np.zeros((height, width))for i in range(height):for j in range(width):img[i, j] = 0for i in range(round(height / 3), 2 * round(height / 3)):for j in range(round(width / 3), 2 * round(width / 3)):img[i, j] = 255return img@cuda.jit()def ideal_LowPassFiltering(d0, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)if D <= d0:h[i, j] = 1else:h[i, j] = 0@cuda.jit()def butterWorth_LowPassFiltering(d0, n, height, width, h): # 巴特沃斯低通滤波器x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)h[i, j] = 1 / (1 + 0.414 * (D / d0) ** (2 * n))@cuda.jit()def exponential_LowPassFiltering(d0, n, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)h[i, j] = math.exp(-0.347 * (D / d0) ** n)@cuda.jit()def trapezoidal_LowPassFiltering(d0, d1, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)if D <= d0:h[i, j] = 1else:if D <= d1:h[i, j] = (D - d1) / (d0 - d1)else:h[i, j] = 0if __name__ == "__main__":plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 显示中文d0 = args.d0d1 = args.d1n = args.nimg = generate_img(4096, 4096) # 读取图像u = Utils(img)img_sp = u.sp_noise(0.1)[1]plt.subplot(231)plt.imshow(img_sp, cmap="gray")plt.axis(False)plt.title("原始图像")f = np.fft.fft2(img_sp)f_shift = np.fft.fftshift(f)plt.subplot(232)plt.imshow(np.log(1 + np.abs(f_shift)), cmap="gray")plt.axis(False)plt.title("原始图像傅里叶频谱")threadsperblock = (32, 32)blockspergrid_x = int(math.ceil(img.shape[0] / threadsperblock[0]))blockspergrid_y = int(math.ceil(img.shape[1] / threadsperblock[1]))blockspergrid = (blockspergrid_x, blockspergrid_y)h = cuda.to_device(np.zeros(img.shape))D = cuda.to_device(np.zeros(img.shape))if args.method == 1:cuda.synchronize()ideal_LowPassFiltering[blockspergrid, threadsperblock](d0, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("理想低通滤波")if args.method == 2:cuda.synchronize()butterWorth_LowPassFiltering[blockspergrid, threadsperblock](d0, n, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("巴特沃斯低通滤波", fontsize=20)if args.method == 3:cuda.synchronize()exponential_LowPassFiltering[blockspergrid, threadsperblock](d0, n, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("指数低通滤波", fontsize=20)if args.method == 4:cuda.synchronize()trapezoidal_LowPassFiltering[blockspergrid, threadsperblock](d0, d1, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("梯形低通滤波", fontsize=20)del hnew_f = np.fft.ifftshift(out)new_image = abs(np.fft.ifft2(new_f))plt.subplot(234)plt.imshow(new_image, cmap="gray")plt.title("滤波后图像")plt.axis(False)plt.subplot(233)plt.imshow(np.log(1 + np.abs(out)), cmap="gray")plt.title("滤波后傅里叶频谱")plt.axis(False)plt.subplot(235)plt.imshow(np.log(1 + abs(np.fft.fftshift(np.fft.fft2(new_image)))), "gray")plt.title("滤波后所得图像傅里叶频谱")plt.axis(False)plt.show()
utils.py(这个主要是两个加噪声的函数，高斯噪声和椒盐噪声）
import randomimport matplotlib.pyplot as pltimport numpy as npclass Utils:def __init__(self, image):self.image = image# plt.imshow(self.image,"gray")# plt.show()def sp_noise(self, prob):"""添加椒盐噪声image:原始图片prob:噪声比例"""output = np.zeros(self.image.shape, np.uint8)noise_out = np.zeros(self.image.shape, np.uint8)thresh = 1 - probfor i in range(self.image.shape[0]):for j in range(self.image.shape[1]):rdn = random.random() # 随机生成0-1之间的数字if rdn < prob: # 如果生成的随机数小于噪声比例则将该像素点添加黑点，即椒噪声output[i][j] = 0noise_out[i][j] = 0elif rdn > thresh: # 如果生成的随机数大于（1-噪声比例）则将该像素点添加白点，即盐噪声output[i][j] = 255noise_out[i][j] = 255else:output[i][j] = self.image[i][j] # 其他情况像素点不变noise_out[i][j] = 100result = [noise_out, output] # 返回椒盐噪声和加噪图像return resultdef gauss_noise(self, mean=0, var=0.001):"""添加高斯噪声image:原始图像mean : 均值var : 方差,越大，噪声越大"""image = np.array(self.image / 255, dtype=float) # 将原始图像的像素值进行归一化，除以255使得像素值在0-1之间noise = np.random.normal(mean, var ** 0.5, image.shape) # 创建一个均值为mean，方差为var呈高斯分布的图像矩阵out = image + noise # 将噪声和原始图像进行相加得到加噪后的图像if out.min() < 0:low_clip = -1.else:low_clip = 0.out = np.clip(out, low_clip, 1.0) # clip函数将元素的大小限制在了low_clip和1之间了，小于的用low_clip代替，大于1的用1代替out = np.uint8(out * 255) # 解除归一化，乘以255将加噪后的图像的像素值恢复noise = noise * 255return [noise, out]
实验结果
2、锐化
import mathimport timeimport cv2import numpy as npfrom matplotlib import pyplot as pltfrom numba import cudafrom utils import Utilsimport argparseparser = argparse.ArgumentParser("")parser.add_argument("--file", type=str)parser.add_argument("--method", default=1, type=int, choices=[1, 2, 3, 4])parser.add_argument("--d0", default=40, type=int)parser.add_argument("--d1", default=80, type=int)parser.add_argument("--n", default=2, type=int)args = parser.parse_args()def generate_img(height, width):img = np.zeros((height, width))for i in range(height):for j in range(width):img[i, j] = 0for i in range(round(height / 3), 2 * round(height / 3)):for j in range(round(width / 3), 2 * round(width / 3)):img[i, j] = 255return img@cuda.jit()def ideal_HighPassFiltering(d0, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)if D <= d0:h[i, j] = 0else:h[i, j] = 1@cuda.jit()def butterWorth_HighPassFiltering(d0, n, height, width, h): # 巴特沃斯低通滤波器x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)h[i, j] = 1 / (1 + (d0 / D) ** (2 * n))@cuda.jit()def exponential_HighPassFiltering(d0, n, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)h[i, j] = math.exp(-1 * (d0 / D) ** n)@cuda.jit()def trapezoidal_HighPassFiltering(d0, d1, height, width, h):x0 = round(height / 2)y0 = round(width / 2)i = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.xj = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.yif i < height and j < width:D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)if D < d1:h[i, j] = 0else:if D <= d0:h[i, j] = (D - d1) / (d0 - d1)else:h[i, j] = 1if __name__ == "__main__":plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 显示中文d0 = args.d0d1 = args.d1n = args.nimg = cv2.imread(args.file, 0)plt.subplot(231)plt.imshow(img, cmap="gray")plt.axis(False)plt.title("原始图像")f = np.fft.fft2(img)f_shift = np.fft.fftshift(f)plt.subplot(232)plt.imshow(np.log(1 + np.abs(f_shift)), cmap="gray")plt.axis(False)plt.title("原始图像傅里叶频谱")threadsperblock = (32, 32)blockspergrid_x = int(math.ceil(img.shape[0] / threadsperblock[0]))blockspergrid_y = int(math.ceil(img.shape[1] / threadsperblock[1]))blockspergrid = (blockspergrid_x, blockspergrid_y)h = cuda.to_device(np.zeros(img.shape))D = cuda.to_device(np.zeros(img.shape))if args.method == 1:cuda.synchronize()ideal_HighPassFiltering[blockspergrid, threadsperblock](d0, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("理想高通滤波", fontsize=20)if args.method == 2:start = time.time()cuda.synchronize()butterWorth_HighPassFiltering[blockspergrid, threadsperblock](d0, n, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()end = time.time()print(end - start)out = np.multiply(f_shift, h)plt.suptitle("巴特沃斯高通滤波", fontsize=20)if args.method == 3:cuda.synchronize()exponential_HighPassFiltering[blockspergrid, threadsperblock](d0, n, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("指数高通滤波", fontsize=20)if args.method == 4:cuda.synchronize()trapezoidal_HighPassFiltering[blockspergrid, threadsperblock](d0, d1, f_shift.shape[0], f_shift.shape[1], h)h = h.copy_to_host()out = np.multiply(f_shift, h)plt.suptitle("梯形高通滤波", fontsize=20)del hnew_f = np.fft.ifftshift(out)new_image = abs(np.fft.ifft2(new_f))plt.subplot(234)plt.imshow(new_image, cmap="gray")plt.title("滤波后图像")plt.axis(False)plt.subplot(233)plt.imshow(np.log(1 + np.abs(out)), cmap="gray")plt.title("滤波后傅里叶频谱")plt.axis(False)plt.subplot(235)plt.imshow(np.log(1 + abs(np.fft.fftshift(np.fft.fft2(new_image)))), "gray")plt.title("滤波后所得图像傅里叶频谱")plt.axis(False)plt.show()
四种低通滤波处理噪声的优缺点如下：理想低通滤波器由于高频成分包含有大量的边缘信息，因此采用该滤波器在去噪声的同时将 会导致边缘信息损失而使图像边模糊。巴特沃斯低通滤波器它的特性是连续性衰减，而不象理想滤波器那样陡峭变化，即明显的不连续 性。因此采用该滤波器滤波在抑制噪声的同时，图像边缘的模糊程度大大减小， 没有振铃效应产生。指数低通滤波器采 用 该 滤 波 器 滤 波 在 抑 制 噪 声 的 同 时 ， 图 像 边 缘 的 模 糊 程 度 较 用 Butterworth 滤波产生的大些，无明显的振铃效应。梯形低通滤波器采 用 该 滤 波 器 滤 波 在 抑 制 噪 声 的 同 时 ， 图 像 边 缘 的 模 糊 程 度 较 用 Butterworth 滤波产生的大些，无明显的振铃效应。2、图像频域锐化：分析不同滤波器对噪声、边缘的处理效果及其优缺点。理想高通滤波器：理想高通有明显振铃现象，即图像的边缘有抖动现象。 Butterworth 高通滤波器：效果较好，但计算复杂，其优点是有少量低频通过， H(u，v) 是渐变的，振铃现象不明显。 指数高通滤波器：效果比 Butterworth 差些，振铃现象不明显。 梯形高通滤波器：会产生微振铃效果，但计算简单，较常用。
