【Python CUDA版】河北工业大学计算机图像处理实验四：频域平滑与锐化

一、实验目的与要求

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 math
import timeimport cv2
import numpy as np
from matplotlib import pyplot as plt
from numba import cudafrom utils import Utils
import 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 random
import matplotlib.pyplot as plt
import 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 math
import timeimport cv2
import numpy as np
from matplotlib import pyplot as plt
from numba import cudafrom utils import Utils
import 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 差些，振铃现象不明显。

梯形高通滤波器：会产生微振铃效果，但计算简单，较常用。

【Python CUDA版】河北工业大学计算机图像处理实验四：频域平滑与锐化相关推荐

【Python CUDA版】河北工业大学计算机图像处理实验二：图像直方图及灰度变换
一.实验目的与要求 1．掌握图像灰度直方图的概念及其计算方法,编写灰度直方图统计程序. 2．通过对图像直方图的分析,学习应用直方图法解决诸如图像二值化等具体问题. 3．熟悉直方图均衡化的计算过程及其应 ...
【Python CUDA版】河北工业大学计算机图像处理实验五：图像分割
一.实验目的与要求 1．了解图像分割的意义和常用方法 2．掌握常用的图像分割方法二.实验相关知识图像分割就是把图像分成互不重叠的区域并提取出感兴趣目标的技术,是由图像处理到图像分析的关键步骤. 和 ...
中南民族大学计算机图像处理实验报告,中南民族大学数字图像处理程序及图像...
数字图像处理实验课 matlab实验程序以及图像结果附带有命名图仅供参考定有错误木有绝对正确实验1 1-1 close all; clear; f=zeros(40,40);% f(10: ...
中南民族大学计算机图像处理实验报告,中南民族大学数字图像处理实验报告.pdf...
WORD格式院系:计算机科学学院专业:计算机科学与技术年级: 2011 级课程名称:数字图像处理组号: 02 姓名 ( 学号 ) :谢枫石小飞黄煜柳卫平李春豪指导教师:徐胜舟 20 ...
中南民族大学计算机图像处理实验报告,中南民族大学数字图像处理实验报告.docx...
院系:计算机科学学院专业: 计算机科学与技术年级: 2011 级课程名称: 数字图像处理组号: 02 姓名(学号): 谢枫石小飞黄煜柳卫平李春豪指导教师:徐胜舟 2014年4月29日 ...
计算机图像处理实验四：图像对比度调整（MATLAB）
[实验名称] 图像对比度调整 [实验目的] 通过本次实验掌握图像对比度自动调整技术: .熟悉使用库函数imhist.histeq: [实验内容] 利用库函数imhist计算并分别显示下图(图1)中图 ...
大学计算机基础实践教程第二版素材,920730-大学计算机基础-实验素材-大学计算机基础教程2...
<920730-大学计算机基础-实验素材-大学计算机基础教程2>由会员分享,可在线阅读,更多相关<920730-大学计算机基础-实验素材-大学计算机基础教程2(7页珍藏版)>请 ...
华北电力大学计算机图形学实验报告,华北电力大学计算机图形学实验报告分析.doc...
华北电力大学计算机图形学实验报告分析科技学院课程设计(综合实验)报告 ( 2013 -- 2014 年度第 2 学期) 实验名称 OpenGL基本图元绘制实验课程名称计算机图形学 | | ...
计算机改考408的学校,2020考研河北工业大学计算机改考408
2020考研河北工业大学计算机改考408~今年,改考408的学校真的是越来越多了!又有一所211大学,计算机专业改考408了,这回是--河北工业大学. 河北工业大学虽然名字上带"河北&qu ...

【Python CUDA版】河北工业大学计算机图像处理实验四：频域平滑与锐化

【Python CUDA版】河北工业大学计算机图像处理实验四：频域平滑与锐化相关推荐

最新文章

热门文章