图像处理基础与图像变换

一、目的（本次实验所涉及并要求掌握的知识点）

       通过本实验加深对数字图像的理解，熟悉MATLAB中的有关函数；应用DCT对图像进行变换；熟悉图像常见的统计指标，实现图像几何变换的基本方法。

二、实验内容与设计思想（设计思路、主要代码结构、主要代码段的文字分析、输出结果截图）

（1）选择两幅图像，读入图像并显示，同时使用Matlab计算图像的大小、灰度平均值、协方差矩阵、灰度标准差和相关系数。

1. 设计思路

这段代码的主要目的是读取两张图片（一张狗的图片和一张猫的图片），将它们转换为灰度图，并确保它们具有相同的尺寸。然后，计算并显示每张图片的基本属性（如大小、平均灰度值、标准差）以及两张图片之间的协方差矩阵和相关系数。

2. 主要代码结构

• 读取图片：使用 imread 函数读取图片文件。
• 转换为灰度图：如果图片是彩色的（RGB格式），使用 rgb2gray 函数将其转换为灰度图。
• 调整图片尺寸：如果两张图片的尺寸不同，使用 imresize 函数将它们调整为相同的尺寸。
• 显示图片：使用 imshow 函数在图形窗口中显示图片。
• 计算图片属性：计算图片的大小、平均灰度值和标准差。
• 计算协方差矩阵和相关系数：使用 cov 和 corrcoef 函数计算两张图片之间的协方差矩阵和相关系数。
• 输出结果：使用 disp 函数输出图片的属性和统计信息。

3. 主要代码段的文字分析

• 显示图像

figure;

subplot(1,2,1);

imshow(imgDog);

title('Dog Image');

subplot(1,2,2);

imshow(imgCat);

title('Cat Image');

这段代码展示了：

1. 创建一个新的图形窗口；
2. 使用 subplot 函数将图形窗口分为两个子图；
3. 使用 imshow 函数显示图像，并使用 title 函数为每个子图添加标题。

• 计算统计信息

sizeDog = size(imgDog);

sizeCat = size(imgCat);

meanDog = mean2(imgDog);

meanCat = mean2(imgCat);

stdDog = std2(imgDog);

stdCat = std2(imgCat);

covDog = cov(double(imgDog(:)));

covCat = cov(double(imgCat(:)));

corrDog = corr2(double(imgDog), double(imgDog));

corrCat = corr2(double(imgCat), double(imgCat));

这段代码展示了：

1. 使用 size 函数获取图像的大小；
2. 使用 mean2 函数计算灰度平均值。
3. 使用 std2 函数计算灰度标准差。
4. 将图像展平为一维向量，然后使用 cov 函数计算协方差矩阵。
5. 使用 corr2 函数计算自相关系数。

4. 输出结果截图

（2）图像加法去除噪声

选择一幅图像，并对其添加高斯噪声

I = imread('lena.tif'); % 读取图像

J = imnoise(I,'gaussian',0,0.01); % 添加高斯噪声，均值为0，方差为0.01

重复上述步骤多次，再将所得到的图像相加取平均，检验加法去噪的效果，写出加法去噪的原理。

1. 设计思路

• 读取一张图像并将其转换为灰度图像（如果它是彩色的）。
• 向图像添加高斯噪声，生成多个噪声图像。
• 计算这些噪声图像的平均值，得到去噪图像。
• 显示原始图像、一个噪声图像和去噪图像。
• 计算并显示每个图像的平均灰度值和标准偏差

2. 主要代码结构

• 读取图像：使用 imread 函数读取图像文件。
• 图像转换：若图像是彩色的，使用 rgb2gray 函数将其转换为灰度图像。
• 添加噪声：使用 imnoise函数向图像添加高斯噪声，生成多个噪声图像。
• 计算平均图像：计算多个噪声图像的平均值，得到去噪图像。
• 显示图像：使用 imshow函数显示原始图像、一个噪声图像和去噪图像。
• 计算统计信息：使用 mean2 和 std2 函数计算每个图像的平均灰度值和标准偏差。
• 显示结果：使用 disp 函数显示计算的统计信息。

3. 主要代码段的文字分析

• 生成噪声图像

num_images = 10;

noisy_images = cell(1, num_images);

for i = 1:num_images

    noisy_images{i} = imnoise(I, 'gaussian', 0, 0.01);

end

这段代码展示了：

定义生成噪声图像的数量 num_images。

使用 cell 函数创建一个单元数组 noisy_images 来存储噪声图像。

使用 for 循环生成 num_images 个噪声图像，并将它们存储在 noisy_images 单元数组中。

使用 imnoise 函数向图像添加高斯噪声，其中 0 是均值，0.01 是方差。

• 计算平均图像

average_image = zeros(size(I));

for i = 1:num_images

    average_image = average_image + double(noisy_images{i});

end

average_image = uint8(average_image / num_images);

这段代码展示了：

初始化一个与原始图像大小相同的零矩阵 average_image。

使用 for 循环将所有噪声图像相加，并存储在 average_image 中。

将 average_image 转换为 uint8 类型，以确保其值在有效范围内。

4. 输出结果截图

（3）FFT变换

 f=zeros(30,30);

f(5:24,13:17)=1;

figure(1)

imshow(f,'notruesize');

F=fft2(f);

F2=log(abs(F));

figure(2)

imshow(F2,[-1 5],'notruesize');

colormap(jet);

①运行程序，观察实验结果。

使用 jet 颜色映射，不同的频率成分将以不同的颜色显示，通常低频成分（接近图像中心）会显示为蓝色或绿色，高频成分（远离图像中心）会显示为红色或黄色。

将代码

F=fft2(f);

F2=log(abs(F));

修改为

F=fft2(f,256,256);

F2=fftshift(F);

②观察实验结果,写出两个结果的差别

• 零频率成分的位置：

代码1：零频率成分位于图像的左上角。

代码2：零频率成分位于图像的中心，使用 fftshift 函数将零频率成分移到图像中心。这使得频谱图的对称性更加明显，更容易观察和分析。

• 对数变换：

代码1：直接对傅里叶变换结果取对数，可能会有一些暗区，因为对数值较小的部分会被压缩到较低的亮度。

代码2：对傅里叶变换结果取绝对值后加1再取对数，避免了对数运算中的负无穷大问题，使得图像的亮度分布更加均匀，对比度更好。

③改变f的大小和f中白色矩阵的大小，说明观察结果和原来不同的原因

• 将 f 的大小改为 60x60：

图像的分辨率提高：因为 f 的大小增加了，所以傅里叶变换后的图像会有更高的分辨率，可以显示更多的细节。

频谱图的变化：由于白色矩形相对于整个图像的比例减小了，因此在频谱图中，代表该矩形的特征会变得更弱，即对应的频率分量强度降低。

• 保持 f 的大小为 60x60，但将白色矩形的大小改为更大一些：

频谱图的变化：随着白色矩形尺寸的增加，频谱图中代表矩形边缘的高频分量会减少，而低频分量会增加。这是因为较大的均匀区域对应于较低的空间频率。

对比度的变化：如果白色矩形非常大，接近整个图像的大小，那么频谱图中的对比度会降低，因为大部分能量都集中在低频部分，导致图像看起来较为平滑。

（4）DCT变换

  RGB = imread('autumn.tif');

I = rgb2gray(RGB);

fugure(1);

Imshow(I);

J = dct2(I);

Figure(2);

imshow(log(abs(J)),[]), colormap(jet(64)), colorbar

观察实验结果，说明figure2中图像色彩变换的规律，并写出其原因。

• 规律有：

1. 中心区域：在图像的中心区域，看到较亮的颜色。
2. 边缘区域：在图像的边缘区域，会看到较暗的颜色。
3. 色彩分布：jet 颜色表从蓝色逐渐过渡到绿色再到黄色和红色，最后到白色。

② 原因：

1. DCT 的性质：DCT 将图像从空间域转换到频率域。在频率域中，低频分量表示图像的整体结构和平均亮度，而高频分量表示图像的细节和边缘信息。低频分量通常具有较高的幅度值，因此在对数尺度下表现为较亮的颜色。高频分量通常具有较低的幅度值，因此在对数尺度下表现为较暗的颜色。
2. 对数变换：log(abs(J)) 是对 DCT 结果进行对数变换，目的是压缩动态范围，使图像中的细节更加明显。对数变换可以将大的幅度值压缩到较小的范围内，从而使图像中的高频和低频分量都能清晰可见。
3. 颜色映射：colormap(jet(64)) 使用了一种常见的颜色映射表，从蓝色到红色渐变。这种颜色映射有助于区分不同频率分量的强度，使得图像中的频率分布更加直观。

（5）DCT与图像压缩

   I=imread('cameraman.tif');

I=im2double(I);

T=dctmtx(8);

B=blkproc(I,[8 8],'P1*x*P2',T,T');

mask=[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];

B2=blkproc(B,[8 8],'P1.*x',mask);

I2=blkproc(B2,[8 8],'P1*x*P2',T',T);

subplot(1,2,1),imshow(I);

subplot(1,2,2),imshow(I2);

• 运行代码，观察实验结果，写出图像压缩的基本原理。

观察到了一个原始图像和一个压缩后的图像，基本原理是：利用图像数据中的冗余信息，通过数学变换和量化来减少数据量，从而达到压缩的目的。

• 计算原图像和压缩后图像的差（相减）

1. 设计思路

1. 读取和预处理图像：读取图像并将其转换为双精度类型，以便进行数学运算。
2. 离散余弦变换 (DCT)：将图像从空间域转换到频率域，以便进行频域处理。
3. 量化：通过掩码矩阵去除高频分量，保留低频分量，实现图像压缩。
4. 逆离散余弦变换 (IDCT)：将处理后的频域图像转换回空间域，得到压缩后的图像。
5. 计算和显示差值图像：计算原始图像和压缩后图像的差值，以便观察压缩效果。

2. 主要代码结构

1. 读取和预处理图像
2. DCT 变换
3. 量化处理
4. IDCT 变换
5. 计算和显示差值图像

3. 主要代码段的文字分析

diff_image = I - I2 计算原始图像 I 和压缩后图像 I2 的差值。

• 修改mask矩阵中非0元素的个数，写出修改后的压缩的图像和原图像之间的差和非0元素个数的关系，并比较变换后系数的重要性

mask中非0系数的数量直接影响了图像压缩后的质量，在非0系数大于10时，压缩后的图像和原始图像差别肉眼观察不出，在只有3个非0系数时，图像失真严重，细节损失最多，可见在非0系数较少（小于10）时，非0系数的个数直接影响了压缩后图像的质量。

非0系数越少，压缩前后图像的差越大，并且观察图像可以知道，差值在图像的边缘处更加明显，说明图像压缩对于物体边缘细节的影响最大。

• 使用别的图像进行上述实验，验证结论正确与否

（6）使用C++与OpenCV实现（5）

1. 设计思路

1. 读取图像：从指定路径读取灰度图像，并将其转换为双精度浮点数类型。
2. 计算DCT矩阵：生成8x8的离散余弦变换矩阵。
3. 块处理：将图像分成8x8的块，对每个块应用DCT变换。
4. 应用掩码：对DCT变换后的系数应用掩码，保留低频部分，去除高频部分。
5. 逆DCT变换：对应用掩码后的块进行逆DCT变换，恢复图像。
6. 显示结果：显示原始图像和处理后的图像。

2. 主要代码结构

1. 头文件和命名空间：引入必要的OpenCV和标准库头文件，并使用相应的命名空间。
2. DCT矩阵计算函数：定义 dctmtx 函数计算8x8的DCT矩阵。
3. 块处理函数：定义 blkproc 函数对图像进行块处理，包括DCT变换和逆DCT变换。
4. 主函数：

读取图像并转换为双精度浮点数。

计算8x8的DCT矩阵。

对图像进行块处理，应用DCT变换。

定义并应用掩码。

对处理后的图像进行逆DCT变换。

归一化处理后的图像。

显示原始图像和处理后的图像。

3. 主要代码段的文字分析

cv::Mat dctmtx(int n) {

    cv::Mat T(n, n, CV_64F);

    double c = sqrt(1.0 / n);

    for (int i = 0; i < n; ++i) {

        for (int j = 0; j < n; ++j) {

            if (i == 0) {

                T.at<double>(i, j) = c;

            } else {

                T.at<double>(i, j) = sqrt(2.0 / n) * cos(M_PI * (2 * j + 1) * i / (2 * n));

            }

        }

    }

    return T;

}

定义 dctmtx 函数，计算8x8的DCT矩阵。

使用 cv::Mat 创建一个8x8的矩阵。

根据DCT公式计算每个元素的值。

4. 输出结果截图

（7）编写程序实现图像的水平垂直和中心对称（注:不能用系统的旋转函数）

     

     

1. 设计思路

• 读取图像：首先读取图像文件 cameraman.tif。
• 显示原始图像：使用 subplot 创建一个 2x2 的子图布局，并在第一个位置显示原始图像。
• 图像变换：分别调用 horizontalFlip、verticalFlip 和 centerSymmetry 函数对图像进行水平翻转、垂直翻转和中心对称变换。
• 显示变换后的图像：在子图布局的其他位置显示变换后的图像。
• 定义局部函数：在脚本的末尾定义 horizontalFlip、verticalFlip 和 centerSymmetry 函数，分别实现水平翻转、垂直翻转和中心对称变换。

2. 主要代码结构

• 主脚本部分：

1. 读取图像。
2. 显示原始图像。
3. 调用变换函数并显示变换后的图像。

② 局部函数部分：

1. horizontalFlip：实现水平翻转。
2. verticalFlip：实现垂直翻转。
3. centerSymmetry：实现中心对称。

1. 主要代码段的文字分析

• I_horizontal = horizontalFlip(I);

subplot(2, 2, 2);

imshow(I_horizontal);

title('水平翻转');

功能：

1. 调用 horizontalFlip 函数对图像进行水平翻转，并在子图布局的第二个位置显示结果。
2. I_horizontal = horizontalFlip(I)：调用 horizontalFlip 函数，传入图像 I，并将结果存储在 I_horizontal 中。
3. subplot(2, 2, 2)：选择子图布局的第二个位置。
4. imshow(I_horizontal)：显示水平翻转后的图像。
5. title('水平翻转')：为当前子图添加标题。

  ② function flippedImage = horizontalFlip(image)

    [rows, cols] = size(image);

    flippedImage = uint8(zeros(rows, cols));

    for r = 1:rows

        for c = 1:cols

            flippedImage(r, cols - c + 1) = image(r, c);

        end

    end

end

功能：实现图像的水平翻转。

[rows, cols] = size(image)：获取图像的行数和列数。

flippedImage = uint8(zeros(rows, cols))：创建一个与原图像大小相同的零矩阵，并将其数据类型设置为 uint8。

for r = 1:rows 和 for c = 1:cols：使用双重循环遍历图像的每一个像素。

flippedImage(r, cols - c + 1) = image(r, c)：将原图像的第 c 列像素值复制到新图像的第 cols - c + 1 列位置，实现水平翻转。

③ function symmetricImage = centerSymmetry(image)

    [rows, cols] = size(image);

    symmetricImage = uint8(zeros(rows, cols));

    for r = 1:rows

        for c = 1:cols

            symmetricImage(rows - r + 1, cols - c + 1) = image(r, c);

        end

    end

end

功能：实现图像的中心对称变换。

[rows, cols] = size(image)：获取图像的行数和列数。

symmetricImage = uint8(zeros(rows, cols))：创建一个与原图像大小相同的零矩阵，并将其数据类型设置为 uint8。

for r = 1:rows 和 for c = 1:cols：使用双重循环遍历图像的每一个像素。

symmetricImage(rows - r + 1, cols - c + 1) = image(r, c)：将原图像的第 r 行第 c 列像素值复制到新图像的第 rows - r + 1 行第 cols - c + 1 列位置，实现中心对称变换。

3. 输出结果截图

三、实验使用环境（本次实验所使用的平台和相关软件）

MATLAB R2022b

四、实验小结（实验中遇到的问题及解决过程、实验中产生的错误及原因分析、实验体会和收获）

1. 通过实验，我更加深入地理解了数字图像的本质——即图像通过数字化的方式表示为矩阵，矩阵中的每个元素表示图像的一个像素值。通过MATLAB的各种图像处理函数，如imread、imshow、rgb2gray、imresize等，我能够方便地加载、显示、转换、调整图像大小等操作。对MATLAB中图像矩阵的操作也让我更加熟悉了图像数据的存储方式和处理流程。
2. 在实验中，我通过实现离散余弦变换（DCT），学习了如何将图像从时域转换到频域。DCT是一种常用的图像变换技术，能够有效地对图像进行压缩，因为它能够将大部分能量集中在图像的少数几个频率分量中。通过MATLAB中的dct2函数，我能够方便地对图像进行DCT变换，查看图像的频域表示，并通过阈值压缩频域信息，实现图像压缩。通过这个实验，我更好地理解了图像在频域中的特性以及如何通过频域变换来达到图像处理的目的。
3. 题（3）中的代码一开始在MATLAB R2022b中无法运行，后来我去查找资料发现是因为imshow 函数中的 'notruesize' 参数不再被支持。在较新的 MATLAB 版本中，imshow 函数的参数有所变化。于是我的修改为：

移除无效参数 'notruesize'：原代码中，imshow 函数调用时使用了 'notruesize' 参数，这是无效的参数，会导致 MATLAB 报错。修改后的代码中，去掉了这两个地方的 'notruesize' 参数。

添加 axis image 控制图像比例：在每个 imshow 调用之后，添加了 axis image; 命令。axis image; 确保图像按实际比例显示，即每个像素都是正方形。这有助于避免图像变形。

五、参考文献（书籍、网址链接等）

1. CSND
2. 《数字信号与图像处理》

六、源代码附件（完整的程序源代码，注意排版紧凑，如果代码较长，字体、字间距选择小一点的）

%% 题（1）

imagePathDog = 'D:\大三上\数字信号与图像处理\实验\实验三\dog1.jpg';

imagePathCat = 'D:\大三上\数字信号与图像处理\实验\实验三\cat1.jpg';



imgDog = imread(imagePathDog);

imgCat = imread(imagePathCat);



figure;

subplot(1,2,1);

imshow(imgDog);

title('Dog Image');



subplot(1,2,2);

imshow(imgCat);

title('Cat Image');



if size(imgDog, 3) == 3

imgDog = rgb2gray(imgDog);

end



if size(imgCat, 3) == 3

imgCat = rgb2gray(imgCat);

end



sizeDog = size(imgDog);

sizeCat = size(imgCat);



meanDog = mean2(imgDog);

meanCat = mean2(imgCat);



stdDog = std2(imgDog);

stdCat = std2(imgCat);



covDog = cov(double(imgDog(:)));

covCat = cov(double(imgCat(:)));



corrDog = corr2(double(imgDog), double(imgDog));

corrCat = corr2(double(imgCat), double(imgCat));



disp('狗狗图像:');

disp(['大小: ', num2str(sizeDog)]);

disp(['灰度平均值: ', num2str(meanDog)]);

disp(['灰度标准差: ', num2str(stdDog)]);

disp(['协方差矩阵: ', num2str(covDog)]);

disp(['相关系数: ', num2str(corrDog)]);



disp('猫猫图像:');

disp(['大小: ', num2str(sizeCat)]);

disp(['灰度平均值: ', num2str(meanCat)]);

disp(['灰度标准差: ', num2str(stdCat)]);

disp(['协方差矩阵: ', num2str(covCat)]);

disp(['相关系数: ', num2str(corrCat)]);



%% 题（2）

image_path = 'D:\大三上\数字信号与图像处理\实验\实验三\labixiaoxin.jpg';



I = imread(image_path);

if ndims(I) == 3 && size(I, 3) == 3

I = rgb2gray(I);

end



num_images = 10;

noisy_images = cell(1, num_images);



for i = 1:num_images

noisy_images{i} = imnoise(I, 'gaussian', 0, 0.01);

end



average_image = zeros(size(I));

for i = 1:num_images

average_image = average_image + double(noisy_images{i});

end

average_image = uint8(average_image / num_images);



figure;

subplot(1,3,1); imshow(I); title('原始图像属性');

subplot(1,3,2); imshow(noisy_images{1}); title('噪声图像属性');

subplot(1,3,3); imshow(average_image); title('去噪图像属性');



disp('原始图像属性:');

disp(['平均灰度值: ', num2str(mean2(I))]);

disp(['标准偏差: ', num2str(std2(I))]);



disp('噪声图像属性:');

disp(['平均灰度值: ', num2str(mean2(noisy_images{1}))]);

disp(['标准偏差: ', num2str(std2(noisy_images{1}))]);



disp('去噪图像属性:');

disp(['平均灰度值: ', num2str(mean2(average_image))]);

disp(['标准偏差: ', num2str(std2(average_image))]);



%% 题（3）①

f = zeros(30, 30);

f(5:24, 13:17) = 1;

figure(1);

imshow(f);

axis image;

F = fft2(f);

F2=log(abs(F));

figure(2);

imshow(F2, [-1,5]);

axis image;

colormap(jet);



%% 题（3）②

f = zeros(30, 30);

f(5:24, 13:17) = 1;

figure(1);

imshow(f);

axis image;

F = fft2(f, 256, 256);

F2 = fftshift(F);

F2 = log(abs(F2) + 1);

figure(2);

imshow(F2, [-1, 5]);

axis image;

colormap(jet);



%% 题（3）③

f = zeros(60, 60);

f(10:48, 26:34) = 1;

figure(1);

imshow(f);

axis image;

F = fft2(f, 256, 256);

F2 = fftshift(F);

F2 = log(abs(F2) + 1);

figure(2);

imshow(F2, [-1, 5]);

axis image;

colormap(jet);



f = zeros(60, 60);

f(10:50, 25:35) = 1;

figure(1);

imshow(f);

axis image;

F = fft2(f, 256, 256);

F2 = fftshift(F);

F2 = log(abs(F2) + 1);

figure(2);

imshow(F2, [-1, 5]);

axis image;

colormap(jet);



%% 题（4）

RGB = imread('autumn.tif');

I = rgb2gray(RGB);

figure(1);

imshow(I);

J = dct2(I);

figure(2);

imshow(log(abs(J) + 1), []), colormap(jet(64)), colorbar



%% 题（5）①

I=imread('cameraman.tif');

I=im2double(I);

T=dctmtx(8);

B=blkproc(I,[8 8],'P1*x*P2',T,T');

mask=[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];

B2=blkproc(B,[8 8],'P1.*x',mask);

I2=blkproc(B2,[8 8],'P1*x*P2',T',T);

subplot(1,2,1),imshow(I);

subplot(1,2,2),imshow(I2);



%% 题（5）②

I = imread('cameraman.tif');

I = im2double(I);

T = dctmtx(8);

B = blkproc(I, [8 8], 'P1*x*P2', T, T');

mask = [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];

B2 = blkproc(B, [8 8], 'P1.*x', mask);

I2 = blkproc(B2, [8 8], 'P1*x*P2', T', T);

diff_image = I - I2;

figure;

subplot(1, 3, 1);

imshow(I);

title('原始图像');

subplot(1, 3, 2);

imshow(I2);

title('压缩后的图像');

subplot(1, 3, 3);

imshow(diff_image, []);

title('差值图像');

colorbar;



%% 题（5）③

I = imread('cameraman.tif');

I = im2double(I);

T = dctmtx(8);

B = blkproc(I, [8 8], 'P1*x*P2', T, T');

mask1 = [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];



mask2 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



mask3 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



B2_1 = blkproc(B, [8 8], 'P1.*x', mask1);

B2_2 = blkproc(B, [8 8], 'P1.*x', mask2);

B2_3 = blkproc(B, [8 8], 'P1.*x', mask3);



I2_1 = blkproc(B2_1, [8 8], 'P1*x*P2', T', T);

I2_2 = blkproc(B2_2, [8 8], 'P1*x*P2', T', T);

I2_3 = blkproc(B2_3, [8 8], 'P1*x*P2', T', T);



diff_image_1 = I - I2_1;

diff_image_2 = I - I2_2;

diff_image_3 = I - I2_3;



figure;

subplot(3, 3, 1);

imshow(I);

title('原始图像');



subplot(3, 3, 2);

imshow(I2_1);

title('压缩后的图像 (10 非零元素)');



subplot(3, 3, 3);

imshow(diff_image_1, []);

title('差值图像 (10 非零元素)');

colorbar;



subplot(3, 3, 4);

imshow(I);

title('原始图像');



subplot(3, 3, 5);

imshow(I2_2);

title('压缩后的图像 (3 非零元素)');



subplot(3, 3, 6);

imshow(diff_image_2, []);

title('差值图像 (3 非零元素)');

colorbar;



subplot(3, 3, 7);

imshow(I);

title('原始图像');



subplot(3, 3, 8);

imshow(I2_3);

title('压缩后的图像 (15 非零元素)');



subplot(3, 3, 9);

imshow(diff_image_3, []);

title('差值图像 (15 非零元素)');

colorbar;



%% 题（5）④

% 读取图像

I = imread('D:\大三上\数字信号与图像处理\实验\实验三\cat1.jpg');

I = im2double(I);



% 获取 8x8 的 DCT 矩阵

T = dctmtx(8);



% 定义不同的掩码矩阵

mask1 = [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];



mask2 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



mask3 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



% 对每个通道进行处理

I2_1 = zeros(size(I));

I2_2 = zeros(size(I));

I2_3 = zeros(size(I));



for channel = 1:3

% 对当前通道进行 DCT 变换

B = blkproc(I(:, :, channel), [8 8], 'P1*x*P2', T, T');

% 对 DCT 变换后的图像应用不同的掩码

B2_1 = blkproc(B, [8 8], 'P1.*x', mask1);

B2_2 = blkproc(B, [8 8], 'P1.*x', mask2);

B2_3 = blkproc(B, [8 8], 'P1.*x', mask3);

% 对应用掩码后的图像进行逆 DCT 变换

I2_1(:, :, channel) = blkproc(B2_1, [8 8], 'P1*x*P2', T', T);

I2_2(:, :, channel) = blkproc(B2_2, [8 8], 'P1*x*P2', T', T);

I2_3(:, :, channel) = blkproc(B2_3, [8 8], 'P1*x*P2', T', T);

end



% 计算差值图像

diff_image_1 = I - I2_1;

diff_image_2 = I - I2_2;

diff_image_3 = I - I2_3;



% 显示结果

figure;

subplot(3, 3, 1);

imshow(I);

title('原始图像');



subplot(3, 3, 2);

imshow(I2_1);

title('压缩后的图像 (10 非零元素)');



subplot(3, 3, 3);

imshow(diff_image_1, []);

title('差值图像 (10 非零元素)');

colorbar;



subplot(3, 3, 4);

imshow(I);

title('原始图像');



subplot(3, 3, 5);

imshow(I2_2);

title('压缩后的图像 (3 非零元素)');



subplot(3, 3, 6);

imshow(diff_image_2, []);

title('差值图像 (3 非零元素)');

colorbar;



subplot(3, 3, 7);

imshow(I);

title('原始图像');



subplot(3, 3, 8);

imshow(I2_3);

title('压缩后的图像 (15 非零元素)');



subplot(3, 3, 9);

imshow(diff_image_3, []);

title('差值图像 (15 非零元素)');

colorbar;



%% 题（5） ⑤

% 读取图像

I = imread('D:\大三上\数字信号与图像处理\实验\实验三\dog1.jpg');

I = im2double(I);



% 获取 8x8 的 DCT 矩阵

T = dctmtx(8);



% 定义不同的掩码矩阵

mask1 = [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];



mask2 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



mask3 = [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

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0];



% 对每个通道进行处理

I2_1 = zeros(size(I));

I2_2 = zeros(size(I));

I2_3 = zeros(size(I));



for channel = 1:3

% 对当前通道进行 DCT 变换

B = blkproc(I(:, :, channel), [8 8], 'P1*x*P2', T, T');

% 对 DCT 变换后的图像应用不同的掩码

B2_1 = blkproc(B, [8 8], 'P1.*x', mask1);

B2_2 = blkproc(B, [8 8], 'P1.*x', mask2);

B2_3 = blkproc(B, [8 8], 'P1.*x', mask3);

% 对应用掩码后的图像进行逆 DCT 变换

I2_1(:, :, channel) = blkproc(B2_1, [8 8], 'P1*x*P2', T', T);

I2_2(:, :, channel) = blkproc(B2_2, [8 8], 'P1*x*P2', T', T);

I2_3(:, :, channel) = blkproc(B2_3, [8 8], 'P1*x*P2', T', T);

end



% 计算差值图像

diff_image_1 = I - I2_1;

diff_image_2 = I - I2_2;

diff_image_3 = I - I2_3;



% 显示结果

figure;

subplot(3, 3, 1);

imshow(I);

title('原始图像');



subplot(3, 3, 2);

imshow(I2_1);

title('压缩后的图像 (10 非零元素)');



subplot(3, 3, 3);

imshow(diff_image_1, []);

title('差值图像 (10 非零元素)');

colorbar;



subplot(3, 3, 4);

imshow(I);

title('原始图像');



subplot(3, 3, 5);

imshow(I2_2);

title('压缩后的图像 (3 非零元素)');



subplot(3, 3, 6);

imshow(diff_image_2, []);

title('差值图像 (3 非零元素)');

colorbar;



subplot(3, 3, 7);

imshow(I);

title('原始图像');



subplot(3, 3, 8);

imshow(I2_3);

title('压缩后的图像 (15 非零元素)');



subplot(3, 3, 9);

imshow(diff_image_3, []);

title('差值图像 (15 非零元素)');

colorbar;



题（6）

#include <opencv2/opencv.hpp>

#include <iostream>

#include <cmath>

using namespace cv;

using namespace std;

cv::Mat dctmtx(int n) {

    cv::Mat T(n, n, CV_64F);

    double c = sqrt(1.0 / n);

    for (int i = 0; i < n; ++i) {

        for (int j = 0; j < n; ++j) {

            if (i == 0) {

                T.at<double>(i, j) = c;

            } else {

                T.at<double>(i, j) = sqrt(2.0 / n) * cos(M_PI * (2 * j + 1) * i / (2 * n));

            }

        }

    }

    return T;

}

void blkproc(const cv::Mat& src, cv::Mat& dst, const cv::Mat& P1, const cv::Mat& P2) {

    int block_size = 8;

    dst.create(src.size(), CV_64F);

    for (int y = 0; y < src.rows; y += block_size) {

        for (int x = 0; x < src.cols; x += block_size) {

            cv::Mat block = src(cv::Rect(x, y, block_size, block_size));

            cv::Mat transformed_block;

            cv::gemm(P1, block, 1.0, cv::noArray(), 0.0, transformed_block);

            cv::gemm(transformed_block, P2, 1.0, cv::noArray(), 0.0, transformed_block);

            transformed_block.copyTo(dst(cv::Rect(x, y, block_size, block_size)));

        }

    }

}

int main() {

    cv::Mat I = cv::imread("cameraman.tif", cv::IMREAD_GRAYSCALE);

    if (I.empty()) {

        std::cerr << "Error: Could not open or find the image!" << std::endl;

        return -1;

    }

    I.convertTo(I, CV_64F);

I /= 255.0;

cv::Mat T = dctmtx(8);

    cv::Mat B;

    blkproc(I, B, T, T.t());

    cv::Mat mask(8, 8, CV_64F, cv::Scalar(0));

    mask.at<double>(0, 0) = 1; mask.at<double>(0, 1) = 1; mask.at<double>(0, 2) = 1; mask.at<double>(0, 3) = 1;

    mask.at<double>(1, 0) = 1; mask.at<double>(1, 1) = 1; mask.at<double>(1, 2) = 1; mask.at<double>(1, 3) = 1;

    mask.at<double>(2, 0) = 1; mask.at<double>(2, 1) = 1; mask.at<double>(2, 2) = 1;

    mask.at<double>(3, 0) = 1; mask.at<double>(3, 1) = 1;

    cv::Mat B2;

    blkproc(B, B2, mask, cv::noArray());

    cv::Mat I2;

    blkproc(B2, I2, T.t(), T);

    cv::normalize(I2, I2, 0, 1, cv::NORM_MINMAX);

    cv::imshow("Original Image", I);

    cv::imshow("Processed Image", I2);

    cv::waitKey(0);

    cv::destroyAllWindows();

    return 0;

}