基于MATLAB 图像处理大作业

# MATLAB 图像处理 ## 作业要求 参见`docs/图像处理大作业.pdf` ## 图像处理思路 ### 3.1 图像基础知识 #### （1）在图像中画圆 首先确定圆的圆心，为测试图像的中心；然后确定圆的半径，为图像长、宽中的较小值。接下来根据圆的定义画圆：所有与圆心的距离等于半径的点构成一个圆。在这里，为了让画出的圆清晰可见，将所有与圆心距离在 0.98~1.02 倍半径之间的点全部标红，呈现在图像上。取出标红点的关键代码如下： ~~~matlab size_pic = size(hall_gray); idx = zeros([size_pic, 2]); for i = 1: 1: height for j = 1: 1: width idx(i, j, 1) = i; idx(i, j, 2) = j; end end circle_idx = ((idx(:, :, 1) - (height + 1) / 2) .^ 2 + (idx(:, :, 2) - (width + 1) / 2) .^ 2 <= (1.02 * radius) ^ 2 & ... (idx(:, :, 1) - (height + 1) / 2) .^ 2 + (idx(:, :, 2) - (width + 1) / 2) .^ 2 >= (0.98 * radius) ^ 2); ~~~ 最后效果如下：  #### （2）将测试图像涂成国际象棋状的“黑白格” 将测试图像分成 8*8 的区域，然后对每块区域根据国际象棋“黑白格”的规则染色。在黑色的地方，将原图对应区域置为0；白色的地方则保持不变。结果如下  本题代码位于`src/hw1_3_2.m`中。 ### 3.2 图像压缩编码 #### （1）在变换域预处理 图像的预处理是将每个像素灰度值减128，若要在变换域进行，可以利用DCT变换的线性性来操作。 首先对原图像作二维DCT变换，然后对相同大小的、所有元素均为128的矩阵作二维DCT变换。由于全128的矩阵DCT变换结果只有直流分量，即左上角元素，所以只需让原图的DCT变换结果的直流分量减去全128矩阵DCT变换结果的直流分量即可。具体代码如下： ~~~matlab [part_h, part_w] = size(hall_part); DC_128 = dct2(zeros(part_h, part_w) + 128); maphall_2 = dct2(hall_part); maphall_2(1, 1) = maphall_2(1, 1) - DC_128(1, 1); ~~~ 这样操作后，与先预处理再DCT变换的结果相比，计算MSE，为 7.1000e-28，非常小，说明两种算法计算结果相同。 本题代码位于`src/hw2_4_1.m`中。 #### （2）实现二维DCT变换 二维DCT变换，首先是找到DCT算子**D**。我直接根据指导书的公式（2.5），在`matlab`中实现算子**D**。关键代码如下： ~~~matlab function D_mat = DCT2_D(input_dim) row = linspace(0, input_dim - 1, input_dim)'; col = linspace(1, input_dim * 2 - 1, input_dim); cos_mat = cos(row * col * pi / (2 * input_dim)); cos_mat(1, :) = cos_mat(1, :) / sqrt(2); D_mat = sqrt(2 / input_dim) * cos_mat; end ~~~ 得到DCT算子**D**后，就可以根据输入图像的维度，计算二维DCT变换。假设灰度图像维度为 [height, width]，则其二维DCT变换结果可以用以下代码来求得 ~~~matlab [height, width] = size(input_mat); dct_coef = DCT2_D(height) * double(input_mat) * DCT2_D(width)'; ~~~ 我将我实现的DCT变换函数，与`matlab`提供的`dct2`进行比较，计算MSE，为9.9944e-26。MSE非常小，可以认为两种方法计算结果相同。 本题代码位于`src/hw2_4_2.m`中，计算DCT算子的函数封装在`src/DCT2_D.m`中，计算二维DCT变换的函数封装在`src/DCT2.m`中。 #### （3）部分DCT系数置零 将原图按 8*8 分块后，对每块作二维DCT变换，再将DCT系数矩阵右侧四列、左侧四列系数矩阵分别置零，再作逆DCT变换看结果。 理论上分析，DCT系数矩阵左上角元素代表低频分量，右下方系数表示横竖两个方向中迅速变换的高频分量。由于常见的图片通常颜色、亮度都是缓慢变换的，因此DCT系数矩阵越往左上方取值越大，越向右下方取值越小。由于低频分量值大，对图像影响大，高频分量值小，对图像影响小，所以，理论分析DCT系数矩阵右侧四列元素置零对图像影响较小，左侧四列元素置零则会对图像有很大影响。实际结果也正是如此，结果如下：  本题代码位于`src/hw2_4_3.m`中。 #### （4）DCT系数转置、逆时针旋转90°、旋转180° 若DCT系数矩阵**C**转置，考虑DCT变换的公式 $$ C=dct(P)=DPD^T\\ P=idct(C)=D^TCD\\ D^{-1}=D^T $$ 所以**C**转置后再逆变换，结果为 $$ D^TC^TD=(D^TCD)^T=P^T $$ 结果为原图的转置。 由于右上方的系数反映图像中横向变化的纹理强度，逆时针旋转90°后，原本右下方高频分量的元素来到了右上方。由于高频分量元素本身就较小，所以DCT系数矩阵逆时针旋转90°后再逆变换，得到的图像横向几乎不会变化，呈现很强的横向纹理。 若将DCT系数矩阵旋转180°，左上方低频分量将会和右下方高频分量换位，导致恢复的图像低频分量小，高频分量大，从而图像变化剧烈，呈现格子状的纹理。 实际运行时，结果与理论分析相近，如下  本题代码位于`src/hw2_4_4.m`中。 #### （5）系统特性 如果认为差分编码是一个系统，则这个系统方程为 $$ y(n)=x(n-1)-x(n) $$ 使用`freqz`函数查看该系统频率响应如下  由此可以看出，该系统可看作为高通滤波器。DC系数先进行差分编码再进行熵编码，说明DC系数含有较多的高频分量。 本题代码位于`src/hw2_4_5.m`中。 #### （6）DC预测误差与Category 记DC预测误差为`delta`，Category为`C`，则它们之间的关系为 $$ C=ceil(log2(abs(delta)+1)) $$ 由此可以推出，Category表示了DC预测误差的二进制表示位数。 #### （7）设计实现Zig-Zag扫描 由于本次作业中，只需对 8*8 的矩阵进行 Zig-Zag 扫描，因此可以直接将各索引手动编排好，然后直接用 matlab 的数组索引功能完成扫描。扫描代码如下 ~~~matlab function output = zig_zag(input) temp = reshape(input, [64, 1]); output = (temp([ 1, 9, 2, 3, 10, 17, 25, 18, ... 11, 4, 5, 12, 19, 26, 33, 41, ... 34, 27, 20, 13, 6, 7, 14, 21, ... 28, 35, 42, 49, 57, 50, 43, 36, ... 29, 22, 15, 8, 16, 23, 30, 37, ... 44, 51, 58, 59, 52, 45, 38, 31, ... 24, 32, 39, 46, 53, 60, 61, 54, ... 47, 40, 48, 55, 62, 63, 56, 64 ... ])); end ~~~ 首先，我先将输入的 `8*8` 的矩阵变为 `64*1` 的大小。然后，我根据 Zig-Zag 的规则，事先将理论的 Zig-Zag 扫描结果的各位置索引记好。最后，即可直接通过数组索引的方式完成扫描。 为了验证我这一算法是正确的，我设计了如下测试样例进行检验 ~~~matlab test_input = [ [1, 2, 6, 7, 15, 16, 28, 29]; ... [3, 5, 8, 14, 17, 27, 30, 43]; ... [4, 9, 13, 18, 26, 31, 42, 44]; ... [10, 12, 19, 25, 32, 41, 45, 54]; ... [11, 20, 24, 33, 40, 46, 53, 55]; ... [21, 23, 34, 39, 47, 52, 56, 61]; ... [22, 35, 38, 48, 51, 57, 60, 62]; ... [36, 37, 49, 50, 58, 59, 63, 64] ... ]; ~~~ 这一矩阵经过 Zig-Zag 扫描后，理论上应该输出 1~64 的序列，使用我设计的算法扫描得到结果正是如此，说明算法实现无误。 本题代码位于`src/hw2_4_7.m`中。 #### （8）对测试图像分块、DCT与量化 首先是预处理，灰度图各像素减去128；然后是将图像每个 8*8 的块进行DCT变换并点除量化矩阵；最后是将结果依次排列。 前两步的代码如下 ~~~matlab hall_pre = double(hall_gray) - 128; hall_quan = blockproc(hall_pre, [8, 8], @(mat)(zig_zag(round(dct2(mat.data) ./ QTAB)))); ~~~ 重新排列时，我按照从左至右、从上至下的方式来遍历，代码如下 ~~~matlab [height, width] = size(hall_quan); res = zeros([64, height * width / 64]); for i = 0: 1: height / 64 - 1