【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 time

import cv2
import numpy as np
from matplotlib import pyplot as plt
from numba import cuda

from utils import Utils
import argparse

parser = 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] = 0
    for i in range(round(height / 3), 2 * round(height / 3)):
        for j in range(round(width / 3), 2 * round(width / 3)):
            img[i, j] = 255
    return 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if i < height and j < width:
        D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)
        if D <= d0:
            h[i, j] = 1
        else:
            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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if i < height and j < width:
        D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)
        if D <= d0:
            h[i, j] = 1
        else:
            if D <= d1:
                h[i, j] = (D - d1) / (d0 - d1)
            else:
                h[i, j] = 0

if __name__ == "__main__":

    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 显示中文
    d0 = args.d0
    d1 = args.d1
    n = args.n
    img = 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 h

    new_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 np


class 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 - prob
        for 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] = 0
                    noise_out[i][j] = 0
                elif rdn > thresh:  # 如果生成的随机数大于（1-噪声比例）则将该像素点添加白点，即盐噪声
                    output[i][j] = 255
                    noise_out[i][j] = 255
                else:
                    output[i][j] = self.image[i][j]  # 其他情况像素点不变
                    noise_out[i][j] = 100
        result = [noise_out, output]  # 返回椒盐噪声和加噪图像
        return result

    def 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 * 255
        return [noise, out]

实验结果

 

 

 2、锐化

import math
import time

import cv2
import numpy as np
from matplotlib import pyplot as plt
from numba import cuda

from utils import Utils
import argparse

parser = 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] = 0
    for i in range(round(height / 3), 2 * round(height / 3)):
        for j in range(round(width / 3), 2 * round(width / 3)):
            img[i, j] = 255
    return 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if i < height and j < width:
        D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)
        if D <= d0:
            h[i, j] = 0
        else:
            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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if 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.x
    j = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if i < height and j < width:
        D = math.sqrt((i - x0) ** 2 + (j - y0) ** 2)
        if D < d1:
            h[i, j] = 0
        else:
            if D <= d0:
                h[i, j] = (D - d1) / (d0 - d1)
            else:
                h[i, j] = 1


if __name__ == "__main__":

    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 显示中文
    d0 = args.d0
    d1 = args.d1
    n = args.n
    img = 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 h

    new_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 差些，振铃现象不明显。

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