MATLAB高级编程与工程应用人脸识别实验报告+源代码.docx
- 文档编号:28391537
- 上传时间:2023-07-11
- 格式:DOCX
- 页数:52
- 大小:3.91MB
MATLAB高级编程与工程应用人脸识别实验报告+源代码.docx
《MATLAB高级编程与工程应用人脸识别实验报告+源代码.docx》由会员分享,可在线阅读,更多相关《MATLAB高级编程与工程应用人脸识别实验报告+源代码.docx(52页珍藏版)》请在冰豆网上搜索。
MATLAB高级编程与工程应用人脸识别实验报告+源代码
MATLAB高级编程与工程应用
实验四图像处理
第一章基础知识
1、MATLAB提供了图像处理工具箱,在命令窗口输入helpimages可查看该工具箱内的所有函数。
请阅读并大致了解这些函数的基本功能。
大致了解。
2、利用MATLAB提供的ImagefileI/O函数分别完成以下处理:
(a)以测试图像的中心为圆心,图像的长和宽中较小值的一半为半径画一个红颜色的圆;
分析:
直接利用半径条件,满足条件的点将红色元素置为255,绿色和蓝色元素置为0,于是得到如下图像:
源代码:
load('hall_color.mat');
%首先获得三原数组
R=hall_color(:
:
1);
G=hall_color(:
:
2);
B=hall_color(:
:
3);
%将圆上的点改为红色
fori=1:
120
forj=1:
168
a=abs(i-60.5);
b=abs(j-84.5);
d=sqrt(a^2+b^2);
if(abs(d-60)<0.5)
R(i,j)=255;
G(i,j)=0;
B(i,j)=0;
end
end
end
%生成新的矩阵
hall_color1(:
:
1)=R;
hall_color1(:
:
2)=G;
hall_color1(:
:
3)=B;
imshow(hall_color1);
(b)将测试图像涂成国际象棋状的“黑白格”的样子,其中“黑”即黑色,“白”则意味着保留原图。
用一种看图软件浏览上述两个图,看是否达到了目标。
分析:
首先设置标记flag在进行循环,对不同方格实行颜色更改就行。
效果:
源代码:
clearall;
load('hall_color.mat');
R=hall_color(:
:
1);
G=hall_color(:
:
2);
B=hall_color(:
:
3);
flag=1;
fori=1:
8
flag=mod((flag+1),2);
forj=1:
8
if(flag==1)
form=15*(i-1)+1:
15*i
forn=21*(j-1)+1:
21*j
R(m,n)=0;
G(m,n)=0;
B(m,n)=0;
end
end
end
flag=mod((flag+1),2);
end
end
hall_color1(:
:
1)=R;
hall_color1(:
:
2)=G;
hall_color1(:
:
3)=B;
imshow(hall_color1);
用看图软件打开成功:
第二章图像压缩编码
1、图像的预处理是将每个像素灰度值减去128,这个步骤是否可以在变换域进行?
请在测试图像中截取一块验证你的结论。
分析:
可以在变换域进行,这个操作对应于在变换域将直流分量减去128*8*8/8=1024,于是可以得到如下图像:
原图为:
直接将灰度值减去128得到图像为:
通过改变变换域数据的方式得到图像为:
对比发现,两种变换方式得到的预处理后图像都是一样的。
另外需要说明一下,为了让预处理后的图像显示出来不全为黑色,我将原图像的灰度像素均乘以了一个2,这样可以更加方便的对比变换后的图像。
2、请编程实现二维DCT,并和MATLAB自带的库函数dct2比较是否一致。
分析:
由PDF中的补充知识我们可以知道,任意二维矩阵做DCT变换的时候都可以由它和与之大小相同的余弦序列矩阵加权表示,也就是说C=DPD¬T,那么我们先找出相应的D矩阵,于是我写出程序循环得到了相应维度的D矩阵,然后可以直接使用矩阵乘法得到相应的二维DCT变换,这样与MATLAB自带的库函数dct2虽然在代码实现上有一定的差异(dct2中调用了一维DCT函数实现二维DCT),但是结果实际上是一致的。
得到DCT系数矩阵的代码如下所示:
%生成8*8DCT系数矩阵
clc;
%参数定义及初始化
N=8;
DCT=zeros(N,N);
%计算DCT系数矩阵
DCT(1,:
)=sqrt(1/N);
fori=2:
N;
forj=1:
N
DCT(i,j)=sqrt(2/N)*cos((i-1)*(2*j-1)*pi/(2*N));
end
end
3、如果将DCT系数矩阵中右侧四列的系数全部置零,逆变换后的图像会发生什么变化?
选取一块图验证你的结论。
如果左侧的四列置零呢?
分析:
检测图像如下图所示:
如果DCT系数矩阵中右侧四列的系数全部置零,逆变后的图像高频成分将会被抹去,得到图像如下图所示:
比较两个图,可以发现在黑灰、白灰以及黑白交界处有些许变化,右下角的变化较为明显,之后我将两个图像矩阵作差,这样可以更加直观的看到除去高频成分后的差异:
很明显的我们可以看到在上面所谓的交界处差值绝对值交大,而越往相邻颜色差异不大的块走,那个差值的绝对值越小。
进一步调整。
如果DCT系数矩阵中左侧四列的系数全部置零,逆变后的图像直流分量以及低频分量都被抹去,得到的图像如下图所示:
可以看到,完全就是一片黑的。
查看逆变后的矩阵我发现,两次逆变后的矩阵相加之和就是原来的矩阵了,也就是说他们是互补的。
4、若对DCT系数分别做转置、旋转90度和旋转180度操作(rot90),逆变换后恢复的图像有何变化?
选取一块图验证你的结论。
分析:
我选择如下图像作验证:
首先,将DCT系数作转置,逆变换后恢复得到的图像如下图所示:
明显我们可以看到逆变换后恢复的图像是原图像逆时针旋转了90度后的结果,原因是对DCT系数转置,原来各横行频率信息变为各竖列频率信息,对应于原图像的横竖调换。
进一步,将DCT系数做旋转90度操作,逆变换后恢复得到的图像如下图所示:
可以看到图像发生了较大变化,分析其原因主要是DCT系数旋转90度后,频域内原来的竖列直流信息的行列最高频分量变成了图像的直流分量,这个分量较小,而原图像直流分量变成了各列最高频分量的直流分量,这个分量比较大,还原回去的新图像各竖列的波动较大,而整体灰度较暗。
下面,将DCT系数做旋转90度操作,逆变换后恢复得到的图像如下图所示:
观察发现图像的变化更大了,黑白交错得更厉害,分析其原因在于将DCT系数旋转180度后,原图像直流分量编程各横行竖列的最高频率分量,这个量是比较大的,而原图像横行竖列的最高频率分量变为了图像的直流分量,这个量是相当小的,所以说最终逆变后的图像沿横竖两个方向的灰度波动都交大,并且图像整体灰度较暗。
源代码如下:
image_2_4.m
clc;
C=dct2(P);
C1=C';
C2=rot90(C);
C3=rot90(C2);
imtool(idct2(C),[0255]);
imtool(idct2(C1),[0255]);
imtool(idct2(C2),[0255]);
imtool(idct2(C3),[0255]);
5、如果认为差分编码是一个系统,请绘出这个系统的频率响应,说明它是一个___高通__(低通、高通、带通、带阻)滤波器。
DC系数先进行差分编码再进行熵编码,说明DC系数的__低__频率分量更多。
分析:
首先写出系统的差分方程如下:
y(n)=x(n-1)–x(n)
对方程进行Z变换后得到如下传递函数:
然后画出频率相应得到:
由此可知系统呈现高通特性,说明DC系数的低频分量更多,通过差分编码后可以尽可能多的减短码长。
源代码如下:
image_2_5.m
clc,closeall;
a=[1];
b=[-11];
tf(b,a);
figure;
freqz(b,a);
6、DC预测误差的取值和Category值有何关系?
如何利用预测误差计算出其Category?
分析:
Category值就是DC预测误差的绝对值的二进制位长,可以得到一个转换公式:
7、你知道哪些实现Zig-Zag扫描的方法?
请利用MATLAB的强大功能设计一种最佳方法。
分析:
Zigzag扫描实现将每一个8*8矩阵的DCT系数按照特定的顺序存储成一列数据,但是在实现扫描的Z字型移动时,会遇到较多麻烦,在编程实现时,可以建立一个映射关系,采取顺序扫描,离散映射的方法快速获得列矢量各元素。
下表表示矩阵各系数对应编号:
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
按照二维顺序扫描此8*8图像,然后按照对应编号将每个分量映射到矢量的对应位置。
写出函数zigzag()
源代码如下所示:
zigzag.m:
functionZ=zigzag(C)
order=[126715162829;
3581417273043;
49131826314244;
1012192532414554;
1120243340465355;
2123343947525661;
2235384851576062;
3637495058596363;];
Z=zeros(64,1);
fori=1:
8
forj=1:
8
Z(order(i,j))=C(i,j);
end
end
end
8、对测试图像分块、DCT和量化,将量化后的系数写成矩阵的形式,其中每一列为一个块的DCT系数Zig-Zag扫描后形成的列矢量,第一行为各个块的DC系数。
分析:
首先将测试图像分成互补重叠的8*8小块,然后依次对每一个8*8小块进行DCT变换,得到8*8的DCT系数矩阵,系数矩阵除以量化矩阵得到根据权重量化的DCT系数,接着调用7问中的zigzag()函数将量化后的DCT系数按顺序存储成为64*1的列向量,最后将所有列向量按行并起来即可得到用来表征整个图像的DCT系数矩阵,这个整体图像DCT系数矩阵的第一行即为DC系数,其余为AC系数。
原代码如下:
image_2_8.m
clc,closeall,clearall;
loadJpegCoeff.mat;
loadhall.mat;
%实现分块
[ab]=size(hall_gray);
M=8*ones(1,a/8);
N=8*ones(1,b/8);
Block=mat2cell(hall_gray,M,N);
%DCT变换、量化、zigzag扫描
J=zeros(64,a*b/64);%a*b/64为分块数目
fori=1:
a*b/64
Block{i}=double(Block{i})-128;%预处理
C=dct2(Block{i});%DCT变换
C_qual=round(C./QTAB);%量化
J(:
i)=zigzag(C_qual);%zigzag扫描
end
9、请实现本章介绍的JPEG编码(不包括写JFIF文件),输出为DC系数的码流、AC系数的码流、图像高度和图像宽度,将这四个变量写入Jpegcodes.mat文件。
分析:
按照pdf中所给的编码方式分别写出DC编码和AC编码的函数,然后将使用8问方法处理过后的矢量调用这两个函数分别得到了DC系数的码流和AC系数的码流。
源代码如下:
clc,closeall,clearall;
loadJpegCoeff.mat;
loadhall.mat;
%实现分块
[ab]=size(hall_gray);
M=8*ones(1,a/8);
N=8*ones(1,b/8);
Block=mat2cell(hall_gray,M,N);
%DCT变换、量化、zigzag扫描
J=zeros(64,a*b/64);%a*b/64为分块的数目
fori=1:
a*b/64
Block{i}=double(Block{i})-128;%预处理
C=dct2(Block{i});%DCT变换
C_qual=round(C./QTAB);%量化
J(:
i)=zigzag(C_qual);%zigzag扫描
end
%熵编码
DC=DCencode(J(1,:
));%调用函数DCencode实现直流编码
AC=ACencode(J);%调用函数Cencode实现交流编码
%保存
saveJpegcodes.matDCACab;
DCencode.m
functionDC=DCencode(dc)
%读入DCTAB
loadJpegCoeff.mat;
%进行差分运算
[ab]=size(dc);
dc_dif=zeros(1,b);
dc_dif
(1)=dc
(1);
dc_dif(2:
end)=dc(1:
end-1)-dc(2:
end);
DC=zeros(0);
fori=1:
b%计算各category,用于定位DCTAB中对应行数
ifdc_dif(i)~=0;
category=floor(log2(abs(dc_dif(i))))+1;
else
category=0;
end
%提取category的huffman码
len=DCTAB(category+1,1);
dc_huff=DCTAB(category+1,2:
len+1);
%生成相应二进制码
dc_bin=de2bi(abs(dc_dif(i)),'left-msb');
if(dc_dif(i)<0)
dc_bin=~dc_bin;
end
%生成完整编码
if(category==0)
DC=[DC,dc_huff];
else
DC=[DC,dc_huff,dc_bin];
end
end
end
ACencode.m
functionAC=ACencode(J)
%读入ACTAB
loadJpegCoeff.mat;
%初始化
ZRL=[11111111001];
EOB=[1010];
[ab]=size(J);
AC=zeros(0);
fori=1:
b
if(any(J(2:
end,i)))%如果此子向量中存在非0系数
place=find(J(2:
end,i));%定位此序列中向量非0系数所在下标
last=place(end)+1;%校正下标值,应加1
Run=0;
forj=2:
last
if(J(j,i)==0)
Run=Run+1;
if(Run==16)
AC=[AC,ZRL];
Run=0;
end
else
Size=floor(log2(abs(J(j,i))))+1;
len=ACTAB(Run*10+Size,3);
ac_huff=ACTAB(Run*10+Size,4:
len+3);%得到对应huffman编码
ac_bin=de2bi(abs(J(j,i)),'left-msb');%转化为二进制数
if(J(j,i)<0)%对负数取反
ac_bin=~ac_bin;
end
AC=[AC,ac_huff,ac_bin];
Run=0;
end
end
end
AC=[AC,EOB];%插入块结束符
end
end
10、计算压缩比(输入文件长度/输出码流长度),注意转换为相同进制。
分析:
输出码流长度=DC码流长度+AC码流长度+高的二进制位长+宽的二进制位长,由此可以得到输出码流长度=2173+23072+8+8=25261,而输入文件长度=120*168*8=161280,最终我们可以得到压缩比=161280/25261≈6.385。
11、请实现本章介绍的JPEG解码,输入是你生成的Jpegcodes.mat文件。
分别用客观(PSNR)和主观方式评价编解码效果如何。
分析:
按照编码的逆顺序处理了自己生成的Jpegcodes.mat文件,生成如下图像:
按照课件中的方法编程计算出了PSNR系数的大小:
由此不管是从主管观察还是客观的用PSNR系数来判断,都可以看到解码出来的图像与原图像相比的效果是不太好的。
源代码如下:
image_2_11.m
clc,closeall,clearall;
%载入相关数据
loadJPEGCodes.mat;
loadJpegCoeff.mat;
%解码
J_dc=DCdecode(DC,a,b);
J_ac=ACdecode(AC,a,b,ACTAB);
J=cell2mat({J_dc;J_ac});%第一行是直流分量,之后为交流分量
%反zigzag扫描,反量化,反DCT变换
Block=cell(a/8,b/8);
fori=1:
a*b/64
Q=Izigzag(J(:
i));%调用自定义反zigzag函数,恢复量化系数矩阵
C=Q.*QTAB;%反量化
Block{i}=idct2(C);%反DCT变换
Block{i}=Block{i}+128;%反预处理
end
%拼接恢复图像
hall_gray2=cell2mat(Block);
hall_gray2=uint8(hall_gray2);
savehall_gray2.mathall_gray2;
%画出图像进行对比
loadhall.mat
figure;
subplot(1,2,1),imshow(hall_gray,[0255]);
subplot(1,2,2),imshow(hall_gray2);
izigzag.m
functionC=Izigzag(J)
order=[126715162829;
3581417273043;
49131826314244;
1012192532414554;
1120243340465355;
2123343947525661;
2235384851576062;
3637495058596364;];
C=zeros(8,8);
fori=1:
64
C(i)=J(order(i));
end
end
DCdecode.m
functionJ=DCdecode(DC,a,b)
%解码恢复预测误差序列
J_dif=ones(1,a*b/64);
i=1;
j=1;
whilei<=length(DC);
%计算category
if(DC(i)==0)
if(DC(i+1)==0)
category=0;
i=i+2;
elseif(DC(i+2)==0)
category=1;
i=i+3;
else
category=2;
i=i+3;
end
elseif(DC(i+1)==0)
if(DC(i+2)==0)
category=3;
i=i+3;
else
category=4;
i=i+3;
end
elseif(DC(i+2)==0);
category=5;
i=i+3;
elseif(DC(i+3)==0)
category=6;
i=i+4;
elseif(DC(i+4)==0)
category=7;
i=i+5;
elseif(DC(i+5)==0)
category=8;
i=i+6;
elseif(DC(i+6)==0)
category=9;
i=i+7;
elseif(DC(i+7)==0)
category=10;
i=i+8;
elseif(DC(i+8)==0)
category=11;
i=i+9;
end
%二进制转化为十进制
ifcategory==0
J_dif(j)=0;
j=j+1;
else
DC_bin=DC(i:
i+category-1);
if(DC_bin
(1)==0)
DC_bin=~DC_bin;
J_dif(j)=-bi2de(DC_bin,'left-msb');
j=j+1;
else
J_dif(j)=bi2de(DC_bin,'left-msb');
j=j+1;
end
i=i+category;
end
end
%恢复DC系数
J=zeros(1,a*b/64);
J
(1)=J_dif
(1);
fori=2:
a*b/64
J(i)=J(i-1)-J_dif(i);
end
end
ACdecode.m
functionJ=ACdecode(AC,a,b,ACTAB)
%初始化
ZRL=[11111111001];
EOB=[1010];
J=zeros(63,a*b/64);
blocknum=1;%块位置编号
i=1;%列位置编号
j=1;%熵编码位置编号
while(j
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 高级 编程 工程 用人 识别 实验 报告 源代码
![提示](https://static.bdocx.com/images/bang_tan.gif)