《数字图像处理-频域增强实验｜CSDN创作打卡》杭州电子科技大学

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

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

%% 读取图像
x = imread('cameraman.tif');
figure, imshow(x), title('原图')
 
 
%% 傅里叶变换
y = fft2(x);       % 二维傅里叶变换
y = fftshift(y);   % 将零频分量移到频谱中心
y_abs = abs(y);    % 对复数进行模运算 
z = log(y_abs+1);  % 映射到小范围区间
% 根据z中的像素值范围对显示进行转换，将z中的最小值显示为黑色，最大值显示为白色。
figure, imshow(z, []), title('空域图像');
 
%% 低通滤波器
LPF = ones(size(y));        % 初始化低通滤波器
threshold_low = 50;         % 设置截至频率
[row, col] = size(y);       % 得到图像y的高度和宽度
for i = 1: row
    for j = 1: col
        % 计算频率
        if sqrt((i-row/2)^2+(j-col/2)^2) > threshold_low
            % 将低通滤波器中高于截止频率的频率变为0
            LPF(i,j) = 0;
        end
    end
end
y_LPF = y.*LPF;             % 低通滤波
y_LPF = ifftshift(y_LPF);   % 逆零频平移
x_LPF = ifft2(y_LPF);       % 傅里叶逆变换
x_LPF = uint8(real(x_LPF)); % 转换成实数，并将double转成uint8
figure, imshow(x_LPF), title('低通滤波')
 
 
%% 计算频率，频率 = 当前点(i,j) 到 中心点(row/2, col/2)的距离
[row, col] = size(y);
distance = zeros(row, col);
for i = 1: row
    for j = 1: col
        % 计算向量的模长，即两点间距离
        distance(i,j) = norm([i-row/2 j-col/2]);
    end
end
 
 
%% 高通滤波 (HPF, High Pass Filtering)
y_HPF = y;
threshold_high = 30;
y_HPF(distance < threshold_high) = 0;
y_HPF = ifftshift(y_HPF);
x_HPF = ifft2(y_HPF);
x_HPF = uint8(real(x_HPF));
figure, imshow(x_HPF), title('高通滤波')

 

实验心得：

图像傅里叶转换后的结果是复数矩阵，因为复数没法用在图像中显示，所以在图像中我们看到的都是幅度，也就是结果实部和虚部的平方和的平方根，也就是复数的长度，但是复数还有另一个性质，也就是其相位，根据实验显示，相位存储了图像的轮廓信息，或者叫做机构信息，而幅度存储了灰度信息。本质上，理想低通滤波器和理想高通滤波器都是基于频率的，我们可以先计算出每个点对应的频率，将其存放在一个频率矩阵（频率本质上也是距离）distance中。

 

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

问题1：对输入的灰度图像进行分块，每一块图像为8*8像素的大小。对分块图像进行离散余弦变换，输出频谱图（DCT系数）；

问题2：尝试改变部分的DCT系数（例如可以自行设计量化表，或者单纯修改直流、交流系数）；

问题3：通过离散余弦逆变换，还原出图像，观察与原图像之间的区别，分析修改的系数与图像修改效果之间的关联。

I1=imread('1.jpg');
I2=rgb2gray(I1);
 mask0=[1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0];
 
 mask9=[1 1 1 1 0 0 0 0
        1 1 1 0 0 0 0 0
        1 1 0 0 0 0 0 0
        1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0];
 
mask35=[1 1 1 1 1 1 1 1
        1 1 1 1 1 1 1 0
        1 1 1 1 1 1 0 0
        1 1 1 1 1 0 0 0
        1 1 1 1 0 0 0 0
        1 1 1 0 0 0 0 0
        1 1 0 0 0 0 0 0
        1 0 0 0 0 0 0 0];
 
mask53=[1 1 1 1 1 1 1 1
        1 1 1 1 1 1 1 1
        1 1 1 1 1 1 1 1
        1 1 1 1 1 1 1 1
        1 1 1 1 1 1 1 0
        1 1 1 1 1 1 0 0
        1 1 1 1 1 0 0 0
        1 1 1 1 0 0 0 0];
 
    
 
g0=DCTgary(I2,mask0);
g9=DCTgary(I2,mask9); 
g35=DCTgary(I2,mask35); 
g53=DCTgary(I2,mask53);
 
 
 
subplot(2,3,1);imshow(I2);title('灰度图像')
subplot(2,3,2);imshow(g0);title('灰度图像DCT系数进行遮罩 0')
subplot(2,3,3);imshow(g9);title('灰度图像DCT系数进行遮罩 9')
subplot(2,3,4);imshow(g35);title('灰度图像DCT系数进行遮罩 35')
subplot(2,3,5);imshow(g53);title('灰度图像DCT系数进行遮罩 53')
 
 
 
function RGB_rec = DCTRGB( RGB,mask )
%进行彩色图像的DCT变换
%返回重构的图像
 
%亮度量化表
m=0.5*[ 16 11 10 16 24 40 51 61;
        12 12 14 19 26 58 60 55;
        14 13 16 24 40 57 69 56;
        14 17 22 29 51 87 80 62;
        18 22 37 56 68 109 103 77;
        24 35 55 64 81 104 113 92;
        49 64 78 87 103 121 120 101;
        72 92 95 98 112 100 103 99];
%RGB图分层处理 得到3个分量图
R = RGB(:,:,1);
G = RGB(:,:,2);
B = RGB(:,:,3);
%转换为双精度
IR = double(R);
IG = double(G);
IB = double(B);
%建立8*8的DCT变换矩阵
T=dctmtx(8);
%进行DCT变换
RR = blkproc(IR,[8,8],'P1*x*P2',T,T');
GG = blkproc(IG,[8,8],'P1*x*P2',T,T');
BB = blkproc(IB,[8,8],'P1*x*P2',T,T');
 
%量化
LR = blkproc(RR,[8 8], 'round(x./P1)',m);
LG = blkproc(GG,[8 8], 'round(x./P1)',m);
LB = blkproc(BB,[8 8], 'round(x./P1)',m);
 
%对DCT系数进行遮罩处理
mR=blkproc(LR,[8 8],'x.*P1.*P2',mask,m);
mG=blkproc(LG,[8 8],'x.*P1.*P2',mask,m);
mB=blkproc(LB,[8 8],'x.*P1.*P2',mask,m);
 
%反DCT变化 IDCT
YR =blkproc(mR,[8 8],'P1*x*P2',T',T);
YG =blkproc(mG,[8 8],'P1*x*P2',T',T);
YB =blkproc(mB,[8 8],'P1*x*P2',T',T);
 
%转换为uint8
YR = uint8(YR);
YG = uint8(YG);
YB = uint8(YB);
 
%重构图像
RGB_rec =cat(3,YR,YG,YB);
 
end
 
 
function gary_rec = DCTgary(gary,mask )  
%进行灰度图像的DCT变换
%返回重构的图像
 
%亮度量化表
m=0.5*[ 16 11 10 16 24 40 51 61;
        12 12 14 19 26 58 60 55;
        14 13 16 24 40 57 69 56;
        14 17 22 29 51 87 80 62;
        18 22 37 56 68 109 103 77;
        24 35 55 64 81 104 113 92;
        49 64 78 87 103 121 120 101;
        72 92 95 98 112 100 103 99];
 
%转换为双精度
I = double(gary);
%建立8*8的DCT变换矩阵
T=dctmtx(8);
%进行DCT变换
RR = blkproc(I,[8,8],'P1*x*P2',T,T');
 
%量化
LR = blkproc(RR,[8 8], 'round(x./P1)',m);
 
%对DCT系数进行遮罩处理
mR=blkproc(LR,[8 8],'x.*P1.*P2',mask,m);
 
%反DCT变化 IDCT
YR =blkproc(mR,[8 8],'P1*x*P2',T',T);
 
%转换为uint8
YR = uint8(YR);
 
%重构图像
gary_rec = YR;
 
end

 

实验心得：

离散余弦变化-DCT的应用之一是：经常被信号处理和图像处理使用，用于对信号和图像进行有损数据压缩。这是因为DCT具有很强的“能量集中”特性。对于比较平滑的图像/数据，DCT变换数据集中在图像的左上角。对于细节丰富的图像/数据，DCT变化后的数据仍然比较集中，而且从频率谱恢复原始图像，有着不错的效果。