位图的处理算法

1、参赛队号#3335第七届“认证杯”数学中国数学建模网络挑战赛承 诺 书我们仔细阅读了第七届“认证杯”数学中国数学建模网络挑战赛的竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们接受相应处理结果。我们允许数学中国网站()公布论文，以供网友之间学习

2、交流，数学中国网站以非商业目的的论文交流不需要提前取得我们的同意。我们的参赛队号为：参赛队员 (签名) ：队员1： 队员2：队员3： 参赛队教练员 (签名)： 参赛队伍组别：第七届“认证杯”数学中国数学建模网络挑战赛编 号 专 用 页参赛队伍的参赛队号：（请各个参赛队提前填写好）：竞赛统一编号（由竞赛组委会送至评委团前编号）：竞赛评阅编号（由竞赛评委团评阅前进行编号）：2014年第七届“认证杯”数学中国数学建模网络挑战赛第一阶段论文题 目 位图的处理算法 关 键 词 灰色像素矩阵 像素坐标 最小二乘 光顺逼近 摘 要：本文通过分析题中相关要求及条件，建立数学模型解决了简单位图文件的矢量化问题。

3、 针对问题，本文首先将题中所给简单位图文件进行分辨率处理，使其分辨率为，再将图片导入Matlab软件，利用Imread函数得到图片的灰色像素矩阵，利用矩阵得到图片中心点，并以中心点为坐标原点建立直角坐标系，然后在根据矩阵确定位图文件边界线条像素坐标，最后利用B样条曲线的最小二乘光顺逼近模型得到边界线条方程，并打印了边界图像。 利用本文模型，另选分辨率相同的具有边界方程的简单位图文件进行验证，得到的边界线条方程与所知边界线条方程形式相同。从而，验证了上述模型的合理性。参赛密码 （由组委会填写）参赛队号: 3335所选题目: B 题 英文摘要【Summary】: Through analyzing

4、 the demands and conditions which are related to the problem, this thesis establishes the mathematic model to solve the problem of the simple bitmap files vectorization.Aiming at the question, first, this paper handles the simple bitmap files resolution in the problem, making it 136 Multiply 100.Sec

5、ond, get the picture into Matlab , use Imread function to get the gray pixel matrix of the picture, use the matrix to get the central point of the picture, and set up Cartesian coordinates using central point as coordinate origin.Third, ensure the bitmap file boundary line&#

6、160;pixel coordinates according to the matrix.Finally, get the boundary line equation by using the model of Fairing and Approximation Algorithm of B spline curves least square, and print the edge image.Using this model, choose another simple bitmap file which having boundary equation and t

7、he same resolution to prove, and the boundary line equation which we get is same as the one we know. Thereby, prove the rationality of this model. 【Keyword】:Gray pixel matrix Pixel coordinates Least square Fairing and Approximation Algorithm一、问题重述图形（或图像）在计算机里主要有两种存储和表示方法。矢量图是使用点、直线或多边

8、形等基于数学方程的几何对象来描述图形，位图则使用像素来描述图像。一般来说，照片等相对杂乱的图像使用位图格式较为合适，矢量图则多用于工程制图、标志、字体等场合。矢量图可以任意放缩，图形不会有任何改变。而位图一旦放大后会产生较为明显的模糊，线条也会出现锯齿边缘等现象。第一阶段问题：矢量图从本质上只是使用曲线方程对图形进行的精确描述，在以像素为基本显示单元的显示器或打印机上是无法直接表现的。将矢量图转换成以像素点阵来表示的信息，再加以显示或打印，这个过程称之为栅格化（Rasterization），见图1。栅格化的逆过程相对比较困难。假设有一个形状较为简单的图标，保存成一定分辨率的位图文件。我们希望将

9、其矢量化，请你建立合理的数学模型，尽量准确地提取出图案的边界线条，并将其用方程表示出来。二、模型假设1、假设位图文件色调均匀2、假设位图文件黑色部分没有白点3、假设位图文件边缘光滑4、假设位图文件除黑色部分外其他地方均为白色5、假设计算机是在无噪干扰下对位图文件进行处理三、符号说明表示位图文件灰度像素矩阵表示灰色像素矩阵的列数表示灰色像素矩阵的行数表示像素点数值表示某灰度像素点为黑色还是白色表示上、下端点距离坐标原点的距离表示左、右端点距离坐标原点的距离表示灰度像素矩阵某点表示点距离位图文件像素矩阵最近的水平边界距离表示点距离位图文件像素矩阵最近的垂直边界距离表示灰度像素矩阵某列表示三次均匀B

10、样条曲线表示数据点列表示数据点列中某个数据点表示三次B样条基函数,均表示控制点向量表示数据点向量表示偏移权表示广顺权表示基函数的二阶导数四、模型建立与求解4.1问题分析整体来看，本问题要求利用数学方法提取简单位图文件边界线条，并将其用方程表示出来。可以考虑线条是由点构成的，如果可以得到位图文件边界线条的点，再利用点即可得到边界方程。具体操作，考虑所给简单位图文件整个界面仅有黑、白两种颜色色，而没有其他任何特征，而题目要求解决的又是整个图片黑、白交界处的曲线方程问题。因此，可以考虑从黑白交界处图片本身所具有的特征入手，而任何一张图片固有的特征就是图片像素。所以，首先可将题中所给的图片导入Matl

11、ab软件中，以其像素为绩点，得到图片的灰色像素矩阵，这一像素矩阵即可得到表示该图片的所以特征。为了利用简单位图文件的灰色像素矩阵得到边界线条方程，那么就需要建立相应的坐标系，得到边界线条上相应点的坐标，再利用这些坐标来得到相应的曲线方程。因此，为了简化运算，可以考虑使用简单位图文件的中心点作为坐标原点，因为不管位图文件是哪样的，其边界线条肯定会围成一个圈，这样中心作为坐标点，四个象限都具有坐标，表示曲线方程更容易。建好坐标系后，可以利用灰色像素矩阵得到黑、白交界处的点坐标，然后根据所得的点坐标，利用Matlab软件曲线拟合，即可得到黑、白交界处的拟合曲线，也就是所求的边界线条方程。并且利用该模

12、型得到的边界方程与已知的位图文件边界方程进行形式对比，若形式相似说明模型、方法合理，形式不相似那么考虑改进模型、方法或者另选思路。具体操作流程如下：模型检验求解位图文件导入Matlab灰色像素矩阵边界线条坐标模型找到位图中心点光顺逼近模型得到边界方程得到边界线条坐标建立直角坐标系模型检验位图文件已知边界方程不相似相似形式对比模型合理模型不合理比所得边界方程合理改进模型或另选方法图1 解答流程图4.2数据处理将题中所给简单位图文件导入matlab中，编写程序（具体代码见附录1），可得位图文件的灰度像素矩阵,矩阵局部像素如下:图2 矩阵局部像素图根据上述灰色像素矩阵，利用matlab编写程序得到所

13、求位图文件的长、宽距离数据，其中长为136像素，宽为100像素，并且可得位图文件页边界到边界线条的横、纵距离，其中上、下、左、右四个方位横、纵距离如下表。 方位上下左右横向距离68681111纵向距离8255050表1方位距离表4.3模型建立与求解4.3.1边界线条像素坐标模型建立将碎纸片导入matlab编程计算所得的灰色像素矩阵为：由于处理后的简单位图文件的像素为100*136，因此，上述矩阵也为100*136的，矩阵每一列数据即为位图文件相应列的像素值，其中每个像素点表示此处为黑色或白色，用表示某灰度像素点为黑色还是白色，即：由于上述矩阵为100*136，因此矩阵中心点为50*78，也就是

14、说位图文件所建立得直角坐标系原点为该矩阵中心点。假设位图文件上、下、左、右端点坐标分别为：，那么上、下端点距离坐标原点的距离为，左、右端点距离坐标原点的距离为。令为坐标系象限中第列像素矩阵中与0相近的第一个非0 像素点，令该点为。令点距离位图文件像素矩阵最近的水平边界的距离为,距离位图文件像素矩阵最近的垂直边界距离为。那么点的坐标可表示如下：其中，。若灰度像素矩阵中没有而只有的情况即：那么我们选择满足列中所有数据中间点作为所求坐标点。若灰度像素矩阵中出现与交替出现的情况即：那么像上述与0 相邻的均为所求坐标。4.3.2边界线条像素坐标模型求解根据上述边界线条像素坐标模型，利用matlab软件编

15、程求解可得所求简单位图136列中满足条件的所有坐标情况，其部分坐标情况如下表：（具体代码见附录2，全部坐标见附录3）(-10,4)(-9,4)(-8,4)(-7,4)(-5,4)(-4,4)(-3,4)(-2,4)(-1,4)(0,4)(1,4)(2,4)(3,4)(4,4)(5,4)(6,4)(7,4)(8,4)(9,4)(10,4)(-10,5)(-9,5)(-8,5)(-7,5)(-6,5)(-4,5)(-3,5)(-2,5)(-1,5)(0,5)(1,5)(2,5)(3,5)(4,5)(5,5)(6,5)(7,5)(8,5)(9,5)(10,5)(-23,6)(-22,6)(-21,6

16、)(-20,6)(-19,6)(-18,6)(-17,6)(-16,6)(-15,6)(-14,6)(-13,6)(11,6)(12,6)(13,6)(14,6)(15,6)(16,6)(17,6)(18,6)(19,6)(20,6)(21,6)(-22,7)(-21,7)(-20,7)(-19,7)(-18,7)(-17,7)(-16,7)(-15,7)(-14,7)(-13,7)表2 部分坐标情况表4.3.3 B样条曲线最小二乘光顺逼近模型建立1.光顺准则与常用算法 光顺概念涉及到几何外形的美观性，受人们主观因素的影响。迄今为止对光顺性还没有一个统一的标准，在不同文献中对光顺准则有不同的提

17、法。提法虽然有差异，但不同的提法中有许多共同点。可归纳如下：(1)曲线连续；(2)没有多余拐点；(3)曲率变化均匀。光顺准则只是对曲线曲面光顺性的一个大概的、定性的描述，在实际使用时还需要对其作定量的描述。为使生成的曲线具有良好的光顺性，常采用的算法有：采用良好的参数化方法；采用优良的曲线、曲面生成和表示方法；适当调整型值点(或控制顶点)的空间位置等。目前较多采用优化的方法对曲线进行光顺，如最小二乘法。2.最小二乘光顺逼近算法构造一条曲线使之在某种意义下最为接近给定的数据点，称为对这些数据点进行逼近，所构造的曲线称为逼近曲线。常用的方法是使各型值点与逼近曲线上对应点的差方和为最小，即最小二乘逼

18、近法。为了使逼近曲线比较光顺，刘鼎元在曲线逼近时引入了偏离权和光顺权的概念。以三次B样条曲线为例来介绍曲线的最小二乘光顺逼近方法。给定一组数据点列,求一条由个控制点，决定的三次均匀B样条曲线来逼近数据点列，其中为曲线的分段数，且。对每一个数据点，应该选取曲线上哪一点与之作偏差量估计的对应点，是参数曲线逼近中首先要解决的问题。本文选取规范的累加弦长参数作为值，即令：引入段序号，则；取曲线参数，记号代表小于或等于的最大整数，。这条B样条曲线的方程可以写为：其中，其表达式分别为：的矩阵表达式为：其中,为控制点向量，。并且上述公式中一共包含个方程，个未知量，属于矛盾方程组。取最小二乘逼近目标函数为：式

19、中：为数据点向量；为偏移权，越大，点越靠近逼近曲线；为广顺权，越大，逼近曲线越光滑；为基函数的二阶导数，其表达式为：那么，对最小二乘逼近目标函数关于求导，则有最小二乘解的法方程便为：求解上式即可得到逼近曲线的控制点向量,带入的方程即可得到逼近曲线的方程如下：4.3.4 B样条曲线最小二乘光顺逼近模型求解根据上述模型，结合边界线条像素坐标模型求解所得的坐标点，利用matlab软件编程求解可得所求简单位图文件边界线条方程及边界线条图像如下：（具体代码见附录4）边界方程：边界线条图像：图3边界线条图像五、模型检验为检验上述模型的合理性，下面本文利用上述模型结合圆形位图文件。具体图像如下：图4 圆形位

20、图利用上述界线条像素坐标模型可求得圆形位图文件的坐标值如下：（具体代码见附录5）表3 圆形位图文件坐标表进而利用上诉B样条曲线最小二乘光顺逼近模型即可得到上述圆形文图文件的边界线条方程为：上述结果与过原点的圆的标准方程：形式相同，故上述模型、方法合理。六、模型评价与推广5.1模型的优点使用该模型，可以很好的将简单位图文件矢量化，得到简单位图文件的边界线条方程。利用计算机提取照片灰度像素矩阵并利用矩阵得到边界线条点，进而得到边界线条方程及图像，方法科学合理，结果准确。从位图出发逐步得到将其矢量化，以点代面的做法，可使思路更加简单明；运用matlab计算及绘出图形，用数形结合的方法来进行分析，模型

21、思路更加清晰，说服力更强。5.2模型的缺点在使用最小二乘光顺逼近时，有时因为点数量过多而不得不进行删除、调整。5.3模型的推广本文所建模型用于将简单位图文件矢量化，得到其边界方程，在实际生活中，与简单位图文件相关的矢量化问题以及与负责图像相关的相似，匹配，识别等相关问题，均可使用该模型思想进行运算。七、参考文献1 游洪.自由曲面的B样条拟合2朱心雄.白由曲线曲面造型技术，北京：科学出版社，20003刘鼎元，赵玉琦，詹廷雄等Bezier曲线和B样条曲线光顺拟合法，计算数学，1984，6(4)：360365八、附录附录1，计算图片灰色像素矩阵 使用软件：matlab fid=fopen('

22、F:Matlabyuan.txt', 'wt'); %创建记录位图的像素矩阵文件 %img为图片像素矩阵img=imread('F:Matlabpic2.bmp');%row,col为该图片像素矩阵的行和列row,col=size(img);for i=1:row for j=1:col fprintf(fid, '%dt', img(i, j); %类似于C语言的输出格式输出到记事本中 end fprintf(fid, 'n'); end fclose(fid);附录2 求解位图文件边界坐标 使用软件：matlab%im

23、g为图片像素矩阵img=imread('F:Matlabpic2.bmp');%row,col为该图片像素矩阵的行和列row,col=size(img);fid=fopen('F:Matlabt0.xls','wt');tp=0;temp=ones(1,col)*255;temp2=ones(row,1)*255;imgT=zeros(row,col);for i=1:row for j=1:col tp=(img(i,j); %fprintf(fid, '%dt', img(i, j); %类似于C语言的输出格式输出到记事本中

24、if tp>0&&tp<255 if i>row/2 k=(row/2-i); else k=i; end l=j-col/2; fprintf(fid,'%dt%dt%dn',k,l,img(i,j); %k=k+1; imgT(i,j)=imgT(i,j)+img(i,j); end endend fclose(fid);imshow(imgT)附录3 位图文件全部坐标 使用软件：matlab和excel附录4 求解边界线条方程及边界线条图像 使用软件：matlabclc; clear all; im=imread('F:Matla

25、bpic2.bmp'); im0=rgb2gray(im); im1=medfilt2(im0,3 3); BW = edge(im1,'sobel'); %finding edges imx,imy=size(BW); L = bwlabel(BW,8);% Calculating connected components mx=max(max(L) r,c = find(L=3); rc = r c; sx sy=size(rc); n1=ones(imx,imy);%zeros for i=1:sx x1=rc(i,1); y1=rc(i,2); n1(x1,y1

26、)=0; end % Storing the extracted image in an array figure,imshow(n1); title('边缘提取'); S=ellipsefit(c,r) a=S.Phi; jdg(S.Xc,S.Yc,S.A,S.B,a,n1)%plot ellipse function varargout=ellipsefit(x,y)x=x(:); % convert data to column vectors y=y(:); if numel(x)=numel(y) | numel(x)<5 error('X and Y

27、Must be the Same Length and Contain at Least 5 Values.') end D1=x.*x x.*y y.*y; % quadratic terms D2=x y ones(size(x); % linear terms S1=D1'*D1; S2=D1'*D2; Q2,R2=qr(D2,0); if condest(R2)>1.0e10 warning('ellipsefit',. 'Data is Poorly Conditioned and May Not Represent an Ell

28、ipse.') end T=-R2(R2'S2'); % -inv(S3) * S2' M=S1+S2*T; CinvM=M(3,:)/2; -M(2,:); M(1,:)/2; V,na=eig(CinvM); c=4*V(1,:).*V(3,:) - V(2,:).2; A1=V(:,c>0); P=A1; T*A1; % correct signs if needed P=sign(P(1)*P; Phi=atan(P(2)/(P(3)-P(1)/2; c=cos(Phi); s=sin(Phi); % rotate the ellipse para

29、llel to x-axis Pr=zeros(6,1); Pr(1)=P(1)*c*c - P(2)*c*s + P(3)*s*s; Pr(2)=2*(P(1)-P(3)*c*s + (c2-s2)*P(2); Pr(3)=P(1)*s*s + P(2)*s*c + P(3)*c*c; Pr(4)=P(4)*c - P(5)*s; Pr(5)=P(4)*s + P(5)*c; Pr(6)=P(6); % extract other data XcYc=c s;-s c*-Pr(4)/(2*Pr(1);-Pr(5)/(2*Pr(3); Xc=XcYc(1); Yc=XcYc(2); F=-Pr(6) + Pr(4)2/(4*Pr(1) + Pr(5)2/(4*Pr(3); AB=sqrt(F./Pr(1:2:3); A=AB(1); B=AB(2); Phi=-Phi; if A<B % x-axis not major axis, so rotate it pi/2 Phi=Phi-sign(Phi)*pi/2; A=AB(2); B=AB(1); end S.Xc=Xc; S.Yc=Yc; S.A=A; S.B=B; S.Phi=Phi; S.P=P; if nargout=1 vararg