数字图像处理作业：灰度图像的频域滤波器 灰度图像的离散余弦变换

实验二：数字图像频域增强实验

作业一：灰度图像的频域滤波器

问题1 对灰度图像进行DFT变换，在频域上分别使用理想的高通和低通滤波器进行滤波，显示滤波后的频域图像，以及IDFT变换后的空域图像，观察振铃现象。要求实验采用两种不同的截止频率。

DFT变换-傅里叶变换(高低通滤波)

import numpy as np
import cv2
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 使图形中的中文正常编码显示
plt.rcParams['axes.unicode_minus'] = False # 使坐标轴刻度表签正常显示正负号

# 理想低通滤波器
def LowPassFilter(img,D):

    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT) #进行傅里叶变化 cv2.DFT_COMPLEX_OUTPUT表示进行傅里叶变化的方法
    fshift = np.fft.fftshift(dft) # 进行频率变换，将图像中的低频部分移动到图像的中心

    #  cv2.magnitude(x, y) 将sqrt(x^2 + y^2) 计算矩阵维度的平方根 使用cv2.magnitude将实部和虚部转换为实部
    # 乘以20是为了使得结果更大
    img_f = 20 * np.log(cv2.magnitude(fshift[:, :, 0], fshift[:, :, 1])) # 输出傅里叶频域图


    # 设置低通滤波器
    # 定义掩模：生成的掩模中间为1周围为0
    # D:截止频率
    rows, cols = img.shape
    crow, ccol = int(rows / 2), int(cols / 2)# 计算频谱中心
    mask = np.zeros((rows, cols, 2), np.uint8) # 生成rows行cols的矩阵，数据格式为uint8
    mask[crow - D:crow + D, ccol - D:ccol + D] = 1 # 将靠近频谱中心的部分低通信息设置为1

    # 掩膜图像和频谱图像乘积
    f = fshift * mask

    img_l = 20 * np.log(cv2.magnitude(f[:, :, 0], f[:, :, 1])) # 输出滤波后的频域图像

    # 傅里叶逆变换
    ishift = np.fft.ifftshift(f) #使用np.fft.ifftshift 将低频移动到原来的位置
    iimg = cv2.idft(ishift)     #使用cv2.idft进行傅里叶的反变化
    res = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])   #使用cv2.magnitude转化为空间域内
    return img_f, img_l, res


# 理想高通滤波器
def HighPassFilter(img,D):
    # # 第一步读入图片
    # img = cv2.imread('lena.jpg', 0)
    # # 第二步：进行数据类型转换
    img_float = np.float32(img)
    # 第三步：使用cv2.dft进行傅里叶变化
    dft = cv2.dft(img_float, flags=cv2.DFT_COMPLEX_OUTPUT)
    # 第四步：使用np.fft.fftshift将低频转移到图像中心
    dft_center = np.fft.fftshift(dft)
    # 第五步：定义掩模：生成的掩模中间为0周围为1
    crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2)  # 求得图像的中心点位置
    mask = np.ones((img.shape[0], img.shape[1], 2), np.uint8)
    mask[crow - D:crow + D, ccol - D:ccol + D] = 0

    # 第六步：将掩模与傅里叶变化后图像相乘，保留中间部分
    mask_img = dft_center * mask
    img_h = 20 * np.log(cv2.magnitude(mask_img[:, :, 0], mask_img[:, :, 1])) # 输出滤波后的频域图像

    # 第七步：使用np.fft.ifftshift将低频移动到原来的位置
    img_idf = np.fft.ifftshift(mask_img)

    # 第八步：使用cv2.idft进行傅里叶的反变化
    img_idf = cv2.idft(img_idf)

    # 第九步：使用cv2.magnitude转化为空间域内
    img_idf = cv2.magnitude(img_idf[:, :, 0], img_idf[:, :, 1])

    return  img_h, img_idf

imgnum=2
img = cv2.imread('%d.jpg'%imgnum)
grayImg  = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
print('保存灰度图')
cv2.imwrite('grayImg%d.png'%imgnum, grayImg) #保存灰度图
Dl1=50   #高通滤波器截断频率
Dl2=100
Dh1=10
Dh2=50   #低通滤波器截断频率
# 把显示界面分割成2*3的网格，231将参数合起来写，2是行数，3是列数，最后一个1是标号
plt.subplot(121)
plt.axis('off')
plt.title('原始灰度图')
plt.imshow(grayImg, cmap='gray')

plt.subplot(122)
plt.axis('off')
plt.title('DFT图像')
img_f, img_l, res = LowPassFilter(grayImg,Dl1)
plt.imshow(img_f, cmap='gray')
plt.savefig('%d_dft.png'%imgnum)

#################################################################
plt.subplot(221)
plt.axis('off')
plt.title('低通滤波频域图D=%d'%Dl1)
plt.imshow(img_l, cmap='gray')

plt.subplot(222)
plt.axis('off')
plt.title('低通空域图(IDFT)D=%d'%Dl1)
plt.imshow(res, cmap='gray')

#########################
img_f, img_l, res = LowPassFilter(grayImg,Dl2)

plt.subplot(223)
plt.axis('off')
plt.title('低通滤波频域图D=%d'%Dl2)
plt.imshow(img_l, cmap='gray')

plt.subplot(224)
plt.axis('off')
plt.title('低通空域图(IDFT)D=%d'%Dl2)
plt.imshow(res, cmap='gray')

plt.savefig('%d_low.png'%imgnum)
##################################################3
plt.subplot(221)
plt.axis('off')
plt.title('高通滤波频域图D=%d'%Dh1)
img_h, res2 = HighPassFilter(grayImg,Dh1)
plt.imshow(img_h, cmap='gray')

plt.subplot(222)
plt.axis('off')
plt.title('高通空域图(IDFT)D=%d'%Dh1)
plt.imshow(res2, cmap='gray')

######################
img_h, res2 = HighPassFilter(grayImg,Dh2)

plt.subplot(223)
plt.axis('off')
plt.title('高通滤波频域图D=%d'%Dh2)
plt.imshow(img_h, cmap='gray')

plt.subplot(224)
plt.axis('off')
plt.title('高通空域图(IDFT)D=%d'%Dh2)
plt.imshow(res2, cmap='gray')

plt.savefig('%d_hig.png'%imgnum)

# plt.show()
plt.savefig('2dft.png')

作业二：灰度图像的离散余弦变换

问题1 对输入的灰度图像进行分块，每一块图像为8*8像素的大小。对分块图像进行离散余弦变换，输出频谱图（DCT系数）；
问题2 尝试改变部分的DCT系数（例如可以自行设计量化表，或者单纯修改直流、交流系数）；
问题3 通过离散余弦逆变换，还原出图像，观察与原图像之间的区别，分析修改的系数与图像修改效果之间的关联。

1.算法原理
DCT，即离散余弦变换，常用图像压缩算法，步骤如下
1）分割，首先将图像分割成8x8或16x16的小块；
2）DCT变换，对每个小块进行DCT变换；
3）舍弃高频系数（AC系数），保留低频信息（DC系数）。高频系数一般保存的是图像的边界、纹理信息，低频信息主要是保存的图像中平坦区域信息。
4）图像的低频和高频，高频区域指的是空域图像中突变程度大的区域（比如目标边界区域），通常的纹理丰富区域。

离散余弦变换（DCT)是图像频域变换的一种，它类似于离散傅里叶变换(DFT)，离散余弦变换相当于一个长度大概是它两倍的离散傅里叶变换，但是离散余弦变换只使用实数。在傅里叶级数中，如果被展开的函数是实偶函数，那么在傅里叶级数中则只包含余弦项，再将其离散化，由此便可导出离散余弦变化。

变换后DCT系数能量主要集中在左上角，其余大部分系数接近于零，因此DCT具有适用于图像压缩的特性，用于对信号和图像(包括静止图像和运动图像)进行有损数据压缩。
二维离散余弦变换及逆变换的公式分别如下：

代码

# 整张图 DCT 变换
def whole_img_dct(img_f32):
    img_dct = cv2.dct(img_f32)            # 进行离散余弦变换
    img_dct_log = np.log(abs(img_dct))    # 进行log处理
    img_idct = cv2.idct(img_dct)          # 进行离散余弦反变换
    return img_dct_log, img_idct

img_u8 = cv2.imread("3.png", 0)
img_f32 = img_u8.astype(np.float)  # 数据类型转换 转换为浮点型
img_dct_log, img_idct = whole_img_dct(img_f32)

Odd-size DCT’s are not implemented in function ‘cv::OcvDctImpl::apply’

使用函数cvDCT(…)时，它所处理的矩阵的行数和列数必须是偶数。如果行列中任何一个是奇数，则会报错，如下：

The function/feature is not implemented (Odd-size DCT’s are not
implemented)in function cvDCT.

dst = cv2.dct(src), 其中src必须为float类型，uint8不能正常运行。

error: (-215:Assertion failed) type == CV_32FC1 || type == CV_64FC1 in function ‘cv::dct’

img_dct[8*h: 8*(h+1), 8*w: 8*(w+1)] = cv2.dct(img_block)

读取的图片格式不同
直接读取cv2保存的灰度图 不行

img_u8 = cv2.imread('grayImg%d.png'%imgnum)

下面方式读灰度图 可以

img = cv2.imread('1.jpg', 0)