2014第七届优秀第二阶段特等奖

第七届“认证杯”数学中编号专用页数学建模网络赛第二阶段 关键 本文采用领域平均法，把分成N个领域，在每个领域内中心的像素灰度值与其Robert算子、sobel算子、Robinson算子等，通过理论分析和实验结果比较了他们各自的优缺点。最后提出一种新梯度域图像锐化算（填写（填写 §1问题重述第二阶段问题：位图在放大时，图像质量常会有所下降，如容易产生较为明显的2§2模型的假设与符号的约定 所处理的是一幅连续图像g(x,Sxyf(x,*tI光点§3问题的分析，邻的值；界发使得邻看如下四种：二值图像(binaryimages)、亮度图像(intensityimages)、索引图像(indexedimages)RGB图像(RGBimages)。可以分为灰度图像和彩像；按照图像的维数，可分为二维图像，三维图像和图I(xyzt入1f(x,y,z入,t(t)(λ)（xy图像。「f(0, ⁝ ⁝ ⁝ ⁝ 对于平面图像，则与坐标z无关，故f(M1,N1)是一般图像f(x,y,z入,t在不同情况下可分为：二维图像：If(x,y,入,t)；单 像：Ifx,y,z,t；静态图像：If(x,y,z,入)顺序相反。JPEG.jpg.jpeg，一个文件多个图像，可实现动画功能。与JEPG2000图像一样，也实现渐进传输。GIFbmpjpeg§4模型的建立与求解为若干级，数字图像（DIDITALIMAGE）在空间坐标和亮度上都离散化了的图像。数字图像处理时常常是将其视为一个矩阵来处理的。首先对fx,y到一个数字矩fxy,即数字图像可以表示为一个矩阵.矩阵的元素称为数字图像上f(0,

ff

f(0,N1)f(x,y) ff f(M f(M1,N 2「f(0, ⁝ ⁝ ⁝ ⁝ f(M1,N1)可以用灰度0~255值来表示，记为fx,y，他表示图像在水平和垂直方向上的光照强若图像为彩，则是以三基色（RGB）的明亮度来作为函数来表示，即：T

fi,jf(x,y)[fR(x,y),fG(x,y),fB(x,

f(x,y)[fR(i,j),fG(i,j),fB(i,经过离散化的图像，的灰度值任是连续的；如图：连续灰度值灰度标 整数值精确度量（量化 （量化为8bite图像在时，行列点与量化的级数决定了数字图像的的质量和图像数据所占空间的大小，设图像MxN样点fxy代表(x,y)点的灰度值。在一般情况QG2Q G2m(n0,1,2,3;m0,1,2,33BMN8

bMN副N*N个像素的图像f(x,y),，平滑处理后得到一幅图像g(x,y),g(x,y)由下式决定：g(x,y)1 f(m•n) M(m•nxy0,1,2,N1;S(xy点领域中点坐标的集合，但是其中不包括(xy均由；含在(x,y的预定领域中的f(x,y的几个像素的灰度值来决定。平均算法是将原作为新图像中该灰度值。我们用如下方法来表示该操作；3*3模板 11*19

（4-「1 11

1

1 （4-1 21 1

11f(xy)ff(i1,jf(i,jf(i1,jf(i1,jf(i,f(i1,f(i1,jf(i,jf(i1,j4PAGEPAGE7由于图像放大后会出现明显的马赛克现象，使图像模糊不清，上各点像素与邻分成N24和原图局部图像图4-1-0和图4-1- 图4-1- 图4-1- 图4-1- 图4-1-通过进行多次领域平均值法，通过比较相同地方的像素灰度矩阵的值，可看出如图2-图4-2-图4-2-

jzfxyP(xyf→fi gradf(x,y) .(4-gradf(x,y)f→f→f x y zk(4-8个问题中占着很大的。fGx，Gy分别表示沿x方向和y方向的梯度[8]，那么这个梯度矢量可以表示为：f(x,y)[G,G]T[f,f xmag(f)g(x,y)(f(x,y)arctan(f

(4-(4-

x(4-f则在点x,yxy方向的一阶差分可表示为：xf(x,y)f(x1,y)f(x,yf(x,y)f(x,y1)f(x, (4-9求导图像梯度:Gx,ydxidyj;dxi,jIi1,jIi,dyi,jIi,j1Ii, (4- 其中，I是图像像素的值(如：RGB Ii1, Ii,dxi,jdyi,j

2Ii,j1Ii, (4-中分差表示图像的平均梯度(meangradient，是指图像的边界或影线两侧附近灰度程度。平均梯度即图像的清晰度definition，反映图像对细节对比的表达能力。上面说的是简单的梯度定义，其实还有更复杂的梯度。最简单的梯度近Gxf(i,j)f(i1,Gyf(i,j)f(i,j1)(4-

mag(f)Gx

(4-f(x,y变化最快的地方，当图像中存在边缘时，一定有较大的我通过计算每一个点的梯度值，利用调用程序Thegradient以上各式的偏导数需要对每个像素的位置计算，经典的图像梯度算法是考虑图像素域的用缘的阶规对原始图像中像素的某个邻域来设置梯度算子，在实际中常用小区域模板进行卷积来计算。根据模板的大小及权重的不同，人们提出了很多梯度算子，常用的算子有：Sobel算子、RobinsonLacemask波。castlman的数字图像处理对卷积讲得很详细。 D(x)e－入入F(X)Kk

(4.1-其中Dkx为k上图表明f(t)可以展开为正交的抽样函数的无穷级数且级数的系数等于抽样值f(ntifi(t)的每个抽样值上画一个峰值为f(ntiSa合成的波形就是f(t)。另外，我们知道：Sah(t)，h(t)，这些响应进行叠加便得到f(t)，从而达到恢复信号的目的。设f(t)是一带限连续信号，最高频率为,根据定理一对f(t)进行抽样，得到f(nti)，则f(nti)经过一个频率响应为如图的理想低通滤波器后便得到f(t)。由于它是讨论由离散h(x,y)f(x,y)＊g(x,f(x,y)g(xx,yy)dxdy

(4.1-h[i,j]f[i,j]＊g[i,f[kl]g[ik,jk0l

(4.1-积模板原点对应于位置E,而权重A,B,...,I是g[k,l]的值，其中 果f(x,y)和h(x,y)表示图像，则卷积就变成了对像素点的 散系统，脉冲响应函数是一个无穷权序列，系统的输出是输入序列u(t)与权序列h(t)板．对图像中每一像素点[,]ij出响应值h(x,y)是通过平移卷积模板到像素点[,]ij计算模板与像素点[i,j]邻域 板为33的示意图。卷积是线性运算，如下式，因为g[i,j]＊{a1f1[i,j]a2f2[i,j]}a1{g[i,j]＊f1[i,j]}a2{g[i,j]＊f2[i,对任何常量a1和a2都成立．换句话说，和的卷积等于卷积的和，尺度变换后的图像-2SobelSobel算子主要用作边缘检测。在技术上，它是一离散性差分算子，用来运算图像该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以AGx及Gy分别代表经横向及纵向边缘检测的图像Sobel算子，sxsy可用卷积模板来实现，其如下:「 1 「1 Sy 0,Sx 1

1 sx(a2ca3a4)(a0ca7a6sy(a0ca1a2)(a6ca5a4

Sobel提出一种将方向差分运算与局部平均相结合的方法[2]，Sobel算子将图像中 s2s2 用范数衡量梯度的幅度|G(x,y)|性max(|GX|,|Gy|) （4.1-7。Sobel算子对图像 编程Sobel算子得到的图2.3-3Robinson行卷积，用其中三个相邻点的和减去剩下的五个点的和。令3个点环绕不断Kirsch算子是一种非线性算子{10]，Kirsch算子对图像中的每个像素（xy），考虑他的8邻域灰度变化，用其中三个相邻点的和减去剩下的五个点的和。令3个临Kirsch算子使用了8个模板来确定梯度应。所有8个方向中的最大值作为边缘幅度图像的输出。最大响应掩模的序号构成了对边Kirsch算子的梯度幅度值用如下G(x,y)max(|M1|,|M2||M3|,|M4|,|M5|,|M6|,|M7|,|M8 （4.2-Gx,y)max1,max{||5sk4tkk0,17；（4.2-skakak1ak2；tkak3ak4ak7（4.2-「 1 「1 「1 「1 0 0 0211（1）90度方 （2）0度方 （3）270度方 （4）180度方「1 「2 0 「 000002112（5）45度方 （6）315度方 （7）225度方 （8）135度方 Robinson算子各方向梯度0306090120度、18045135算法的效果。同时，这里还有大小不等的圆，验证其他角度的梯度对算法的适用程度4.1.3Gxf(i1,j1)f(i,GyG可得到Roberts

f(i,j1)f(i1,

|G[f(i,j)]|f(i,j)max{|f(i1,j1)f(i,j)|,|f(i,j1)f(i1,j)|}（4.3-从可以看出，其值是与相邻像素成正比的。在图像轮廓上像素的灰度有陡然Roberts§5模型的进一步讨论Roberts边缘检测算子是一种利用局部差分算子寻找边缘的算子,Roberts「 dd d

d dd d

d d d d

1「1 11 11

使用下面的计算梯8

为了抑制噪声的影响，使用下面的分段调制一下EM（EM为相对长度单位，相对于 xxmx3 2 EMLUT(x)

x32 m2Avg

x2xxx xx mx(xx2 xx x x3对于上面的D，采用平滑算法是有一定风险的，可能噪声边界的部分模糊，这里也对于上面的D，采用平滑算法是有一定风险的，可能噪声边界的部分模糊，这里也AvgHYm,nYm,n1Ym3AvgVYm,nYm1,n3AvgLYm,nYm1,n13

AvgRYm,nYm1,n13

=abs(2*Ym,n–Ym,n-1–Ym,n+1)=abs(2*Ym,n–Ym-1,n–Ym+1,n)=abs(2*Ym,n–Ym-1,n-1–;=abs(2*Ym,n–Ym-1,n+1–Ym+1,n-;(5-对于彩像的锐化，通常的做法是对RGB三个通道分别进行锐化，这样做的一个YRGBRGBYUV使用如下：「Y 「 「RU 0.147 0.436 相反，YUV转RGB的

（5- 1.140YG 0.581

（5- B1.0 UVY通道图5-35-5-提取图5-2中梯度强度大于40的区域，如图5-5所示，该高梯度区域的梯度做适当抑5-5-5-5-度如图5-8所示。5-5-。22传统锐化方法噪声局部放大 新锐化方法噪声局部放大 原局部放大 传统锐化方法局部放大§6模型的优缺点§6.1§6.2§7参考文献张德丰，数字图像处理 张德丰，数字图像处理[M]，1 ，2007-01-signalprocessingletter,1999.6(7)onImageProcessing.2003(3)A.Agarwala,M.Dontcheva,M.Agrawala..In ctivedigitalphotomontage[J].InProc.SIGGRAPH04,2004.8.A.ZometandY.Weiss.Seamlessimagestitchinginthegradient[R].InProc.oftheEuropeanConferenceonComputerVision(ECCV04),Prague,Czech，.一种基于梯度域的彩像转灰度图像的方法[J].影像技术.数字图像处理.科学.数字图像处理学[10]RafaelcGonzalezrichardewoodsdigitalimageprocessing(secondedition)[M].peking.publishinghouseofelectronicsindustry.2003.3[11]，.实用[M]..电子工业[12]，.基于算法的图像锐化算法研究[J].电脑知识[13]吴海波,.基于算子的彩像锐化处理[J].电脑开发与应[15]，.基于局部图像锐化的自适应模糊边缘检测算法[J].电脑开发[16]齐欣.数码影像技术操作篇之六十二数码图像锐化方法[J].[17]，程义民，，.基于梯度场的拼接缝消除方法[J].计算机辅助[18]王俊文，，.利用异常边缘进行图像锐化篡改取证[J].理工[19].基于边缘检测的图像锐化算法[J].现代电子技术[20]selimesedoglu.anysisoftheperona-malikscheme[J].commpureapplonOptimalMassPreservingMap[J].IEEETransactionsonImage§8程序1-clearall; 程序1-clearall;figure, figure,imshow(I4); hx=hy';Iyimfilter(double(i),hy, %求x,yIx=imfilter(double(i),hx, %范围在pi/2 arctan(isnan(arctan)==1&Iy~=0)=max(arctan(:)); %将Ix=0的统一设成最大pi %量化成32级，以统计直方图unI %统计arctan中有多少种不同程序1-closeIimread('1灰度图.jpg');hy=fspecial('sobel');hx=Iy=imfilter(double(I),hy,'replicate');Ix=imfilter(double(I),hx,gradmag=sqrt(Ix.^2+figure('units','normalized','position',[001subplot(1,2,1);imshow(I,[]),title('灰度增强图subplot(1,2,2imshow(gradmag/255title('梯度幅值图像');imwrite(gradmag,'2梯度幅值图像.jpg')%保存图像为文件程序1-%Open[filename,pathname]=uigetfile('*.bmp','Pickupatruecolorbmpfile');ifisequal(filename,0)||isequal(pathname,0)%ReadinImage=imread(fullfile(pathname,filename));Gray=Image(:,:,1);[m,n]=%GradY=zeros(m,n);h=[­1,­2,­1;0,0,0;image90=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert90=uint8((image90-Min)*255.0/(Max­Min));title('270h=[­1,0,1;­2,0,2;­image0=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert0=uint8((image0-Min)*255.0/(Max­Min));title('0h=[0,1,2;­1,0,1;­2,­image270=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert270=uint8((image270-Min)*255.0/(Max­Min));title('45h=[1,0,­1;2,0,­2;1,0,-image180=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert180=uint8((image180-Min)*255.0/(Max­Min));title('180h=[1,2,1;0,0,0;1,-2,-image45=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c-1);

ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert45=uint8((image45-Min)*255.0/(Max­Min));title('90h=[­2,­1,0;­1,0,1;image315=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert315=uint8((image315-Min)*255.0/(Max­Min));title('315image225=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert225=uint8((image225-Min)*255.0/(Max­Min));title('135h=[0,­1,­2;1,0,­1;image135=zeros(m,n);forr=2:1:m-2forc=2:1:n-rb=int32(r­1);re=int32(r+1);cb=int32(c­1);ce=GrayLocal=Gray(rb:re,cb:ce);GrayLocalFilter=GrayLocal.*h;ImageConvert135=uint8((image135-Min)*255.0/(Max­Min));title('225程序1-%Open[filename,pathname]=uigetfile('xx.bmp','Pickupatruecolorbmpfile');ifisequal(filename,0)||isequal(pathname,0)%ReadinImage=imread(fullfile(pathname,filename));title('tzmcmR=Image(:,:,G=Image(:,:,B=Image(:,:,title('R通道图像');title('B通道图像');%Y=yuv(:,:,U=yuv(:,:,V=yuv(:,:,h=[­1,0,­1,0,­1;image=imfilter(Y,h/8);title('亮度通道梯度');imageadd=image+Y;title('亮度通道直接锐化%new% idx=find(image>40);imageGradAdd(idx)=%displeyred=image;green=image;green(idx)=0;blue=image;RGB(:,:,2)=RGB(:,:,3)=%idx=find(15<=image&image<40);title('中梯度区域加