图像超分辨率重建--图像处理课程设计
目录
1 课程设计目的............................................................................................................ 1
2图像处理系统设计内容及要求................................................................................. 2
2.1设计内容.......................................................................................................... 2
2.2设计要求.......................................................................................................... 2
3 设计方案.................................................................................................................... 3
4 功能模块的具体实现................................................................................................ 5
4.1 空域插值放大的方法..................................................................................... 5
4.1.1 最邻近插值算法.................................................................................. 5
4.1.2 双线性插值算法.................................................................................. 6
4.1.3 双三次插值算法.................................................................................. 7
4.2 频域重建的方法............................................................................................. 8
4.2.1 DCT变换的介绍.................................................................................. 8
4.2.2 DCT放大图像放大算法原理.............................................................. 8
4.3 频域分块重建的方法................................................................................... 10
4.4 同态滤波器滤波处理................................................................................... 11
4.4.1 同态滤波器原理................................................................................ 11
4.4.2 同态滤波函数的确定........................................................................ 12
5 总结与体会.............................................................................................................. 14
参考文献...................................................................................................................... 15
附录.............................................................................................................................. 16
1课程设计目的
MATLAB7.0软件。MATLAB 是由美国mathworks 公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言的编辑模式，代表了当今国际科学计算软件的先进水平。通过用MATLAB 对图像进行处理，以实现以下目的。
1. 培养严谨的科学态度，正确的设计思想，科学的设计方法和良好的工作作风。
2. 培养独立思考的能力，独立检索资料、阅读文献、综合分析、计算机应用、数据及文字处理等能力。
3. 培养综合运用基础理论、基本知识的能力。通过课程设计得到工程设计的初步锻炼。
2图像处理系统设计内容及要求
2.1设计内容
图像超分辨率重建是利用低质量或低分辨率图像来产生高质量或高分辨率图像的技术，重建包括空域方法和频域方法。本设计要求用插值技术提高图像的分辨率。
（1）利用插值技术将原始图像在空域放大2倍。
（2）在DCT 域放大原始图像2倍，设计滤波器在DCT 域增强图像的高频信息。
（3）对图像分块进行DCT 变换，在DCT 域对子图像进行放大和滤波增强高频信息。
（4）比较上述三种图像重建结果，设计软件界面。
（5）设计方案、编写代码实现上述功能。
2.2设计要求
（1）利用数字图像处理技术，以MATLAB 为平台，建立一个实现设计主题的简易处理系统。
（2）能显示输入图像、中间图像和重建的图像。
（3）程序代码要有注释说明，调用MATLAB 函数要清楚并理解函数的功能、使用范围，在设计说明书中要写清楚函数的功能和参数意义。
（4）完成设计说明书一份。
（5）刻苦钻研，勤于思考，勇于实践，独立完成课程设计任务。
（6）遵守纪律，在指定地点进行课程设计。
（7）掌握有关课程的基本理论和基本知识。概念清楚，方案合理，数据可靠，计算正确，运行良好，图纸（图表）符合标准，设计说明书（论文）撰写规范，答辩中回答问题正确。
3 设计方案
根据课程设计题目的要求，设计界面如图3-1所示：
图3-1 设计界面
本设计分为空域放大和DCT 域放大两部分，空域放大主要工作有：在空域对现有的传统插值算法分别进行了研究与仿真实验，包括最近邻域插值，双线性插值，双三次插值等，这些插值方法均是通过低通滤波，滤除和过滤图像数据中的高频信息。所以这些插值基函数对边缘和纹理信息都比较丰富的图像的插值效果不是特别理想。DCT 域放大主要工作有：通过DCT 变换实现了由空间域到频域的转换，通过对频域处理可以方便的实现空间域较难实现的处理。而空间域与频域又存在一定的联系，为数字图像的处理提供了另一种方法。该算法在对整块图像进行处理时，尽管采用了增强系数对图像亮度效果进行补充，但对整幅图像高频部分预测采用填零方式，在图像像素位数增大即图像信息量增大时这种预测精度不如对图像分块处理后高，且基于JPEG 格式图像多采用分成8×8子块分块压缩编码，对上述算法进行改进。改进后的算法，将原始图像数据切割成接近8×8大小子块，对每一子块分别实施DCT 放大算法。改进后的算法如下：对原始图像进行分块，然后对每一子块运用DCT 图像放大算法，最后合并处理所有的子块。
系统整体设计如图3-2所示。
图3-2 系统整体设计图
4 功能模块的具体实现
4.1空域插值放大的方法
4.1.1 最邻近插值算法
这是最简单的算法，每一个输出像素都赋给输入图象中与其最邻近的采样点的值。插值核函数是：
h(x)=1 0
h(x)=0 0.5
在所有的插值方法中，这种插值方法速度最快，早期的应用比较普遍，然而当图像中包含灰度有变化的细微结构时，最近邻插值法会在图像中产生人工的痕迹。图像的边缘阶梯失真现象比较明显。其实现效果如图4-1所示。
图4-1 最邻近插值算法实现效果图
在程序中可以直接调用函数也可自行编写。如自行编写，部分源程序如下： width = K * nrows;
height = K * ncols;
J = uint8(zeros(width,height));
widthScale = nrows/width;
heightScale = ncols/height;
for x = 5:width - 5
for y = 5:height - 5
xx = x * widthScale;
yy = y * heightScale;
if (xx/double(uint16(xx)) == 1.0) & (yy/double(uint16(yy)) == 1.0)
J(x,y) = I(int16(xx),int16(yy));
else% a or b is not integer
a = double(round(xx)); % (a,b) is the base-dot
b = double(round(yy));
J(x,y) = I(a,b); % calculate J(x,y)
end
end
end
imwrite(J, 'lena2.jpg', 'jpg');
figure;
imshow(J);
4.1.2 双线性插值算法
设f ( x , y ) 为2个变量的函数, 其在单位正方形顶点的值已知。假设希望通过插值得到正方形内任意点的f ( x , y )值。可以令由双线性方程
来定义的一个双曲抛物面与4 个已知点拟合。
利用公式实现插值：
双线性插值是对待插值象素周围的4个邻近像素的灰度按照距离进行加权平均，实质上是属于一阶插值。双线性插值的平滑作用有可能会使图像的细节产生退化，当放大倍数比较大的时候，这种现象更加明显。同时，双线性插值的斜率不连续也会产生不希望的结果。其实现效果如图4-2所示。
图4-2 双线性插值算法实现效果图
在程序中可以直接调用函数也可自行编写。如自行编写，部分源程序如下： width = K * nrows;
height = K * ncols;
J = uint8(zeros(width,height));
widthScale = nrows/width;
heightScale = ncols/height;
for x = 5:width - 5
for y = 5:height - 5
xx = x * widthScale;
yy = y * heightScale;
if (xx/double(uint16(xx)) == 1.0) & (yy/double(uint16(yy)) == 1.0) % if a and b is integer,then J(x,y)
J(x,y) = I(int16(xx),int16(yy));
else
a = double(uint16(xx)); % (a,b) is the base-dot
b = double(uint16(yy));
x11 = double(I(a,b)); % x11
x12 = double(I(a,b+1)); % x12
x21 = double(I(a+1,b)); % x21
x22 = double(I(a+1,b+1)); % x22
J(x,y) = uint8( (b+1-yy) * ((xx-a)*x21 + (a+1-xx)*x11) + (yy-b) * ((xx-a)*x22 +(a+1-xx) * x12) ); % calculate J(x,y)
end
end
end
imwrite(J, 'lena2.jpg', 'jpg');
figure;
imshow(J);
4.1.3 双三次插值算法
双三次插值是用待插枝点周围16个点作为参考像素值的一种三阶插值方法，典型的双三次插值核函数是：
这里参考值a 不同文献中取值不同，a=一1，a=0．5，a=0．75等等。其实现效果如图4-3所示。
图4-3 双三次插值算法实现效果图
4.2 DCT域插值放大的方法
4.2.1 DCT变换的介绍
离散余弦变换是从二种特殊形式的傅里叶变换转化过来的，是一种性能很好的正交变换方式。离散余弦变换本质上仍然是离散傅立叶变换，二者在频域本质上是相同的。离散余弦变换因其是一种实数变换，其变换矩阵的基向量很好地描述了人类视觉的相关性，接近于最佳变换。因而DCT 在图像处理中有很广泛的应用，并成为一些静态图像和视频压缩国际标准的基本处理模块，因而采用DCT 变换可以很方便地应用于压缩域图像和视频中。数字图像可以通过傅里叶变换、离散余弦变换等由空间域转换到频域中表示，通过对频域的处理可以方便实现空间域较难实现的处理。而空间域和频域又存在一定的联系，为数字图像的处理提供了另一种方法。
尺度变换的物理含义就放大意义来讲，如果信号在时域进行扩展，即当O
4.2.2 频域重建的方法
根据图像的表示方法以及以上原理，我们可以仅在频域进行处理就可完成对图像的放大操作。对原图像进行DCT 正变换得到图像频域数据F(u，v) ，将F(u，v) 作为目标放大图像的低频部分，并与增强系数(根据放大k 倍数变化) 相乘，对处理后的频域数据进行DCT 逆变换即可得到放大后的图像。
（1）对原始图像数据进行处理，对其作DCT 正变换。
（2）对数据在频域上进行处理，得到放大后的频域。若需将图像放大k ×k 倍，则放大后的图像大小为M ×N ，M=[km]，N=[kn]([]表示取整运算) 。
（3）DCT 反变换，将频域处理的数据恢复为图像空间域原始格式。其流程图如图4-4所示。
图4-4 DCT放大图像流程图
用DCT 处理图像后的效果图如图4-5所示。
图4-5 频域重建的效果图
部分程序代码如下：
J=dct2(x);%x是读入的图像 [m,n]=size(x) k=2;
T1=[eye(m);zeros(k.*m-m,m)]; T2=[eye(n);zeros(k.*n-n,n)]; J=T1*J*T2'; J1=k*J; X=idct2(J);
figure,imshow(log(abs(J-J1)),[]),colormap(jet(64)); %figure,imshow(log(abs(J1)),[]),colormap(jet(64)); figure;
imshow(double(X)/255); [M,N]=size(X)；
4.3 频域分块重建的方法
该算法在对整块图像进行处理时，尽管采用了增强系数对图像亮度效果进行补充，但对整幅图像高频部分预测采用填零方式，在图像像素位数增大即图像信息量增大时这种预测精度不如对图像分块处理后高，且基于JPEG 格式图像多采用分成8×8子块分块压缩编码，对上述算法进行改进。改进后的算法，将原始图像数据切割成接近8×8大小子块，对每一子块分别实施DCT 放大算法。将图像分解为若干子图像，分别进行DCT 变换，大大减小了DCT 的计算量，提高了算法处理的速度。其实现效果如图4-6所示。
图4-6 用频域分块重建的效果图
部分源程序代码如下：
x=double(x)/255; figure(1); imshow(x); T=dctmtx(8);
B=blkproc(x,[8 8],'P1*x*P2',T,T'); [m,n]=size(x) k=2;
T1=[eye(8);zeros(k.*8-8,8)];
B1=blkproc(B,[8 8],'P1*x*P2',T1,T1'); B2=blkproc(B1,[8 8],'P1*x',k); T2=dctmtx(8.*k);
I=blkproc(B2,[8.*k 8.*k],'P1*x*P2',T2',T2); figure; imshow(I); size(I)
4.4同态滤波器滤波处理
4.4.1 同态滤波器原理
当光照不均匀时, 图像上照度暗的部分细节就难以分辨清楚。但是, 利用光照反射模型的同态滤波方法, 可以消除光照不足带来的影响, 同时又不损失图像的细节。同态滤波增强技术是把频率过滤和灰度变换结合起来的一种图像处理方法, 利用压缩灰度范围和增强对比度来改善图像。
图像从物理过程中产生时, 它的值正比于物理源的辐射能量, 因此, 图像f ( x , y ) 的能量一定是非零且有限的。函数f ( x , y )可由2 个分量来表示: 这2 个分量相应地称为入射分量和反射分量, 分别用i ( x , y )和r ( x , y )表示。那么最后形成的数字图像
f ( x , y )可表示为
f ( x , y ) = i ( x , y ) r ( x , y ) （1）
采用(1)式的成像模型不能直接对照射和反射的频率部分分别进行操作, 因为2 个函数乘积的傅里叶变换是不可分的。即F ( f ( x , y ) ) ≠F ( i ( x , y ) )F( r ( x , y ) ) 。反射分量r ( x , y )反映图像的内容, 它随图像细节的不同在空间作快速变化, 是频域的高频分量。入射分量i( x , y )在空间上常具有缓慢变化的特点, 是频域中的低频分量。采用同态分析方法, 就是先对上式两边取对数, 把2 个相乘的分量变为2个相加的分量, 它们分别代表了图像的高频分量和低频分量。假
设
z ( x , y ) = lnf ( x , y ) = lni ( x , y ) + lnr ( x , y ) （2） 两边取傅里叶变换, 那么有
F( z ( x , y ) ) = F( lni ( x , y ) ) + F( lnr ( x , y ) )（3） 即 Z( u, v) = Fi( u, v ) + Fr ( u, v ) （4）
（4）式表明, 照明分量的频谱可以与反射分量的频谱分离开。根据它们反映的空间变化特征, 照明分量的频谱基本上位于低频部分, 反射分量的频谱位于高频部分。这时, 增强技术都可用于求解这2 个分量。
典型的同态滤波法是: 原图先经对数变换和快速傅里叶变换, 变为频率域中的2 个分离的变量, 然后根据不同需要, 选用不同的传递函数实现不同的增强。经过频域处理的图像再经快速傅里叶逆变换及指数变换, 就可得到增强的图像。
用滤波函数H ( u, v) 对Z( u, v) 进行处理, 那么从（4）式可得： S( u, v ) = H ( u, v) Z( u, v ) =H ( u, v) Fi( u,v) + H (u,v) F r (u, v)（5） 这样, 选择适当的滤波函数H ( u, v ) 就可以对图像中照明分量的反射分量进行不同程度的处理, 使照度不均匀的图像获得良好的改善。滤波后做反变换, 得:
s( x ,y ) = F1 (S( u, v)) = F 1( H( u,v)×Fi( u,v) + H ( u, v) F r( u,v)) （6） 逆傅里叶变换到空域, 再将上式两边取指数, 最后取对数得到滤波后图像的同态滤波函数。
利用照明反射模型处理一幅图像属于同态滤波技术, 同态滤波器H ( u, v)滤波示意图如下图所示。同态滤波最关键的是用（3）式将照射分量和反射分量分开, 用同态滤波器对其进行处理。同态滤波器的过程如图4-7所示：
图4-7同态滤波器工作过程
4.4.2 同态滤波函数的确定
以Rh 代表高频增益, Rl 代表低频增益, D( u,v )表示频率( u,v)距滤波器中心( u0,v0) 的距离。当Rh> 1, Rl
频域内经常使用的高通滤波器为高斯型高通滤波器,
对高斯型高通滤波器稍加修改, 可得以下高斯型同态滤波函数: H ( u, v) = ( Rh- R l) [ 1- exp( c*(-D (u ,v ) / D0^2)))) + Rl（1） 传统的同态滤波函数如图4-8所示。
图4-8传统的同态滤波函数
该同态滤波函数如图4-9所示。
图4-9高斯型同态滤波函数图像
5 总结与体会
为期两周的课设很快接近尾声了，基于MATLAB 的设计已按计划如期全部完成，通过这次课程设计，我对课堂上所学到的理论知识的理解加深了许多，自己动脑、动手设计的能力也得到了较大提高。在这次课程设计的过程中，我对MATLAB 语言有了更深的认识。现在仔细想想，这次课程设计使得我对MATLAB 语言的理解与应用能力得到了较大的提升，也让我认识到只要深入学习，提升的空间永远是存在的。在设计的过程中我遇到了一些问题，如：编写源程序中出现了语法错误等。通过查阅书本和以前设计的程序我发现了产生错误的原因并解决了问题完成了设计。动手实践是理论知识得以灵活运用的必要前提，也是今后走上工作岗位之后能够很好的完成设计工作的技术保证。只有遇到实际问题并根据自己对课堂上获得的专业知识的理解来解决才能真正的提高自己的能力。
在这次课程设计过程中，感触很深，由于对MATLAB 图像处理的函数不熟悉，导致自己走了很多弯路，imshow 函数和image 函数都可以实现显示图像的功能，但是imshow 针对double 型的变量取值范围为0-1，而image 针对unit8型的变量，取值范围为0~255，由于对这个区别不明晰，导致生成全白的图像。由于全局变量的使用不合理导致调用的时候出现错误。还有利用小波函数进行超分辨率重建时由于忽略了对变换后的低频以及高频矩阵的大小而直接进行重建导致矩阵大小不匹配而报错，并且无法实现规定大小的缩放比例。对图像的DCT 域进行块处理时，由于对DCT 域矩阵放大后而直接采用8*8块处理导致图像出现大量黑色网格，未注意到放大后对处理块大小的影响导致了不必要的错误。
在理工科的专业应用背景下，用matlab 进行相关计算与仿真编程的优势非常突出。特定的问题处理算法，我们通常都以.M 文件的文本形式给定最终的解决方案，利用matlab 强大的运算能力，通过设计GUI 界面，使操作变得简单易行，可视化效果非常好。通过本次课程设计，使自己对MATLAB GUI设计流程有了比较深刻的体会，同时也了解了一般软件设计的过程。在设计过程中碰到了很多的问题，通过这些问题，使自己分析问题，解决问题的能力得到了较大的提高。
参考文献
[1] 章毓晋. 图像处理和分析教程．北京. 人们邮电出版社，2009 [2] 龚声蓉. 数字图像处理与分析. 北京. 清华大学出版社，2006 [3] 张强 王正林. 精通MATLAB 图像处理. 北京. 电子工业出版社，2009 [4] 姚敏. 数字图像处理. 北京 . 机械工业出版社，2006
[5] 李显宏.MATLAB7.x 界面设计与编程技巧. 北京. 电子工业出版社，2006 [6] Kenneth R.Castleman著，朱志刚等译，数字图像处理，电子工业出版社，2006
附录
空域重建程序
function[ ]=kong() %空域处理 S=imread('lena.bmp'); [row,line]=size(S);
S_new=S(1:2:row,1:2:line,:); %下采样 subplot(2,2,1);
imshow(S_new);title('低分辨图');
S1=imresize(S_new,2,'nearest'); %最邻近插值 subplot(2,2,2);imshow(S1);title('最邻近插值'); S2=imresize(S_new,2,'bilinear'); %双线性插值 subplot(2,2,3);imshow(S2);title('双线性插值');
S3=imresize(S_new,2,'bicubic'); %三次线性插值 subplot(2,2,4);imshow(S3);title('三次插值'); 频域重建程序 function[ ]=pinyu1() A=imread('lena.bmp'); [row,line]=size(A);
B=A(1:2:row,1:2:line,:); %下采样
J=dct2(B); %将图像变换到频域 subplot(221),imshow(J,colormap(jet(64)));title('128*128DCT图像'); [m,n]=size(B); %利用矩阵的乘法将图像放大 k=2;
T1=[eye(m);zeros(k.*m-m,m)]; %创建256*128的的单位阵和零阵 T2=[eye(n);zeros(k.*n-n,n)]; %创建256*128的的单位阵和零阵 J=T1*J*T2'; %将图像在DCT 域放大两倍，因为原来的图像128*128，现在通过矩阵变换将图像放大到256*256
J1=k*J;
subplot(222),imshow(J1,colormap(jet(64)));title('256*256DCT图像'); I=double(J1);%通过同态滤波器将图像加强 [M,N]=size(I);%获得图片的大小 rL=0.85;
rH=1.02; % 可根据需要效果调整参数 c=1.1;
d0=10;
I1=log(I+1); %取对数 FI=fft2(I1); %傅里叶变换
FI=fftshift(FI); %将FFT 变换后的图片中心移到图像中心 n1=floor(M/2); n2=floor(N/2); for i=1:M for j=1:N
D(i,j)=((i-n1).^2+(j-n2).^2);
H(i,j)=(rH-rL).*(1-exp(c*(-D(i,j)./(d0^2))))+rL; %高斯同态滤波 end end
FI=H.*FI; %进行滤波
FI=fftshift(FI); %将滤波后的图像中心移回来 I2=ifft2(FI); %傅里叶逆变换
I3=real(exp(I2)); %用幂变换将图象恢复过来 subplot(223),
imshow(I3,[]); title('滤波后DCT 图像');
A4=idct2(I3); %将图像变换到空域 subplot(224),
imshow(A4,[]); title('同态滤波增强后图像'); 频域分块重建程序
function pinyu2()%完成频域分块并重建 A=imread('lena.bmp'); [row,line]=size(A); B=A(1:2:row,1:2:line,:); x=double(B)/255; subplot(221),
imshow(x);title('低分辨率图像');
T=dctmtx(8); %计算dct 变换矩阵 B=blkproc(x,[8 8],'P1*x*P2',T,T'); %将图像进行dct 变换 subplot(223);imshow(B); title('8*8DCT图像');
T1=[eye(8);zeros(8,8)]; %频域放大矩阵，由单位阵和零阵构成
B1=blkproc(B,[8 8],'P1*x*P2',T1,T1');%将图像进行频域放大两倍 subplot(224);imshow(B1);title('16*16DCT图像'); k=2; %放大两倍
B2=blkproc(B1,[8 8],'P1*x',k); %增强图像输出的灰度值放大二倍 T2=dctmtx(16); %获得反变换矩阵，此时应该是16*16矩阵
I=blkproc(B2,[16 16],'P1*x*P2',T2',T2);%将图像变换到空域 h4=fspecial('gaussian',[3,3],0.5); %高斯低通滤波3*3 I=filter2(h4,I); subplot(222), imshow(I);
title('滤波后的图像'); 设计界面的程序 function[ ]=zuizhong()
set(gcf,'name','图像重建','numbertitle','off',... 'unit','normalized','position',[0.05,0.1,0.9,0.75],... 'menubar','none'); K=0;
while (K
K=menu('图像重建处理菜单',' 空域重建',' 频域重建',' 频域分块',' 五种对比','Close');
switch K % 选择滤波算法 case 1 kong();
case 2 pinyu1();
case 3 pinyu2();
case 4 sanzhongduibi(); case 5 clear;close;clc; end end