数学建模的格式

2009高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2009高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：双目高清定位摘要:本文是研究数码相机定位问题，先后建立了圆心定位的模型与算法和数码相机标定的模型对于问题一，考虑只有一个圆的圆心在像平面的定位情况首先，建立相机坐标系,利用边缘抽取法，得到圆的边缘像素坐标集然后基于广义霍夫变换，建立圆心定位算法，求得峰值点坐标，作为圆心在像平面的像坐标为求解问题二，先用sobel算子边缘抽取法获取五个圆的边缘像素坐标集采取聚类分割方法对其加以分离，获得五个圆各自的边缘像素坐标集沿用问题一的模型，利用Matlab软件分别求得五个圆的圆心在像平面的像坐标为 (-49.47,-51.06,417.20)，( -22.75,-49.47,417.20)，( 34.39,-45.24,417.20)，( 19.05,31.22,417.20)，(-59.52,31.22,417.20)问题三要求设计一种模型检验方法,由于圆A,B,C的圆心在同一直线上,所以可检验这三个圆心的像坐标是否保持在同一直线上,以此检验模型求解结果的精度. 基于不同的检测算子边缘抽取法,分别求解出相应的结果,并利用这种检验方法求得各个结果的精度.根据精度的偏差,可以分析该检验方法的稳定性. 问题四即是相机标定的问题，先以靶标所在平面建立世界坐标系，求得靶标五个圆的圆心以及AD与CE之间的连线交点F的世界坐标通过问题1的模型，求得五个圆心的像坐标，另外，AD与CE之间的连线交点的坐标即是点F在像平面的像坐标然后利用世界坐标系下的空间点坐标以其相应的像坐标之间的关系：，求得相机的内外参数然后通过空间任一点P分别在两部相机中的坐标关系：，得到这两部相机的相对位置关系关键词:检测算子边缘抽取法霍夫变换聚类分割法相机标定一、问题重述数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置最常用的定位方法是双目定位对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了而它们的像一般会变形，如图1所示，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现 图 1 靶标上圆的像现在，有人设计靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆以AC边上距离A点30mm处的B为圆心，12mm为半径作圆，如图2所示图 2 靶标示意图用一位置固定的数码相机摄得其像，如图3所示图3 靶标的像下面我们需解决一下4个问题：（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；（2） 对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024768；（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法二、 符号说明和基本假设2.1基本假设(1)用相机拍摄相片不会出现叠影或曝光残缺.(2)圆圈上各像素点的灰度值相差不大,相邻像素点的灰度值相差不大.(3)不考虑边缘抽取的像素点残缺较多且分散.2.2符号说明为了问题求解和阅读的方便，我们对全文所用符号进行统一：(u,v)一个点在二维像素坐标平面上的坐标；(x,y)一个点在像平面内以毫米为单位的坐标系上的坐标；O0像素坐标系的坐标原点；O1以物理单位(毫米)表示的图像坐标系得原点；(u0,v0)原点O1的坐标；(Xw,Yw,Zw)一个点在世界坐标系上的坐标；(Xc,Yc,Zc)一个点在相机坐标系的坐标；dx每个像素在x轴方向上的物理尺寸；dy每个像素在y轴方向上的物理尺寸；R世界坐标系到相机坐标系的旋转矩阵；t世界坐标系平移t个单位到相机坐标系三、模型准备为了便于建立模型和求解问题,先建立像平面的坐标系.3.1坐标系的建立1：如图4所示，在像平面上定义直角坐标系，每一像素点的坐标（u，v）分别是该像素点在数组中的列数与行数，而（u，v）是以像素为单位的图像坐标称为像素坐标；像素坐标的集合称为像素坐标集图4 像平面坐标系由于(u,v)只表示像素位于数组中的列数与行数,并没有用长度单位表示出该像素在图像中的位置，因而，为了解题的需要，我们建立一个以物理单位(毫米)表示的像平面坐标系如图4，坐标原点O1为相机光轴与像平面的交点，建立与u,v轴平行且同向的x，y轴，并以毫米为单位的xy坐标系 若（u，v）表示像素坐标，（x，y）表示以毫米为单位的像坐标系的坐标 x，y轴方向上的物理尺寸为dx，dy，则图像中任意一个像素在两个坐标系下的坐标有如下关系：用齐次坐标与矩阵形式表示上式为：逆关系可写成：四、问题1的求解4.1 问题1的分析：首先,为了便于求解,根据题目要求,建立一个相机坐标系.接着,要确定圆心的像坐标,需要建立模型,确定圆心的像坐标与圆面的其他像素点的关系. 边缘是图像最基本的特征 2 ,圆的边缘点与圆心的关系最为紧密.所以需要建立利用边缘点定位圆心坐标的模型.然后,要确定圆的边缘点位置,可用检测算子的边缘抽取法,获得边缘点的像素坐标集.最后,根据模型的算法,利用边缘点的像素坐标集,计算出圆心的像坐标.4.2 问题1的模型建立:4.2.1 建立三维坐标系首先，我们根据题目要求建立相机坐标系，如图5也就是说，以光学中心为原点，过原点垂直于像平面并以指向像平面为正方向的有向直线作为Zc轴可知，Zc轴经过像平面坐标系的原点O1 ,Xc -Yc平面平行于像平面，其中Xc 轴，Yc轴分别平行于x轴，y轴且同向图5 相机坐标系4.2.2 图像空间中圆心的GHT(广义霍夫变换)求解基本原理3 要从一幅图像中检测出半径已知的圆形来.取和图像平面一样的参数平面，以图像上每一个前景点为圆心，以已知的半径在参数平面上画圆面，并把结果进行累加.最后找出参数平面上的峰值点，这个位置就对应了图像上的圆心，其中图像平面上的每一点对应到参数平面上的一个圆面.若图像空间中圆的方程为：.对应于图像空间，参数空间中以(x，y)为圆心，半径为R的圆方程为：.在圆的半径r为常量的前提下，GHT是将图像空间中圆边界上的一个点(x，y)映射参数空间(u，v)一个圆上. GHT求取圆心坐标的算法步骤简述如下：(1)建立起一个离散的参数空间(u，v)，参数空间的量化间隔取决于所求圆心坐标参数的精度.(2)建立一个二维累加器Nf,g，并置每个元素初始值为0.(3)对参数空间的每个点(u,v)，考虑图像空间中每个边界像素点(x,y)，若成立，则相应的累加器Nf,g加1.(4)找出Nf,g的极大值，所求的圆心参数坐标为(f,g).4.2.3 基于GHT的圆心定位模型对于靶标上的圆,用相机摄得其像,在像平面上不再是一个规则的圆.所以,无法简单的用广义霍夫变换求解圆心.而且,对于圆的像,一般会变形,其边界像素点也不能精确的确定.边缘是图像最基本的特征2,基于边的分割方法,找出灰度不连续的点作为边缘.为了在图像中确定边缘点,已经有很多人设计了各种边缘检测算子,然后根据测出的边缘点作为边缘.其中Matlab软件自带的算子检测函数有sobel算子, 梯度算子, prewitt算子, log算子, canny算子.各种检测算子的边缘抽取法,在图像上所确定的边缘点会出现重叠或缺损.如果利用所有边缘点的几何关系,来确定图像的质心,误差较大,且抗干扰性差.所以,利用边缘点的几何关系确定圆心的方法不可取.如果利用图像的所有像素点来确定图像的质心,计算量非常大,不易实现.无论相机从哪个角度摄得圆的像, 靶标上圆内的每一小块(大小取决于摄像时的像距和相机的分辨率)都会在图像的像素点上获得一个灰度值.4虽然圆的像在像平面上不是一个规则的圆,但是其像素点上的灰度值并不相等.而且,对于一个平面图形的像,相邻像素点的灰度值相差不大2.所以,圆心沿着半径到边缘上各像素点的灰度值是连续(在这里,连续指灰度值相差不大)变化的.而且, 边缘上各像素点的灰度值两两彼此相差不大.这样,圆心在图像上的定位受到所有边缘点共同的限制.基于广义霍夫变换,利用边缘点定位圆心的算法如下:(1) 给像平面的每个像素点建立一个计数器,初始值为0.(2) 以检测算子抽取出来的边缘像素点(x,y)为中心,在像平面上画圆.(3) 对于像平面的每个像素点(u,v),如果满足.则像素点(u,v)的计数器Nu,v的值增加,其中为预先设置的误差值,r为半径.(4) 以每个边缘像素点为中心画圆,重复(2)(3)步.(5) 找到Nu,v的极大值,其对应的像素点(u,v)则为圆心的定位.在该模型的求解中, 取值近似于0；r为靶标上圆的半径在图像中的长度,其单位为像素单位.而且,边缘点的重叠和缺损,对模型的求解影响较小,同样对r的取值精度要求也不需要太高.五、问题2的求解5.1问题2的分析：问题2的求解可沿用问题的模型.第一，要求解圆心的像坐标,则需要先抽取出圆的边缘坐标集.我们可用Matlab软件读取图像的数据,并用sobel检测算子抽取图像的边缘坐标集.第二，要利用边缘坐标集给圆心定位,则需先把五个圆的图像的边缘像素坐标集分割出来,分别放在五个边缘像像素坐标矩阵中. 由于五个圆在像平面所成的像之间的欧式像素距离足够大，所以可用聚类的方法加以分割，获得五个圆在像平面上各自的边缘像素坐标集第三，由于边缘抽取法得到的图形不一定是个规则的椭圆,而且有重叠和残缺.所以，利用边缘点的几何关系来确定圆心在像平面的像坐标的方法不可取.因此，我们采用基于广义霍夫变换的算法,即问题1的圆心定位模型5.2 问题2的模型求解：5.2.1 坐标系的建立：我们先建立坐标系如图6,与问题一的相机坐标系和像平面坐标系一致. 图6 相机坐标系Zc轴垂直于像平面，且光学中心到像平面的距离固定.因此,点在Zc轴方向的坐标固定为417.20毫米对于点的坐标，只需要确定Xc轴，Yc轴方向的坐标，也就是要找出点在像平面xy坐标系的坐标 5.2.2 边缘像坐标集的抽取和分割：用sobel检测算子抽取边缘，得到五个圆在像平面上图形的边缘像坐标集由于五个圆的像之间的欧式像素距离足够大，并且同一个圆的边缘像素点非常接近.因此,可以用如下方法对边缘像坐标集进行分割.对于任意两个边缘像素点,如果满足，则将这两个边缘像素点,放进同一个边缘像坐标矩阵中,其中一般取值为1.如果像素点出现残缺,可能一部分像素点被遗漏,但对圆心定位的结果几乎没什么影响.根据假设3, 不考虑边缘抽取的像素点残缺较多且分散的情况.用上述方法将所有抽取出的边缘像素点类聚，获得五个圆在像平面上各自的边缘像素坐标集5.2.3 圆心的定位：根据5.1的分析,可利用基于广义霍夫变换定位圆心的算法(如4.2.3),用Matlab软件求得结果如下：图7 sobel算子边缘抽取及圆心定位结果表1 各个圆在像平面上的像素坐标(u,v)及相机坐标系的像坐标(Xc,Yc, Zc)圆心坐标圆A圆B圆C圆D圆E(u,v）(325,191)(426,197)(642,213)(584,502)(287,502)Xc(毫米）-49.47-22.7534.3919.05-59.52Yc(毫米）-51.06-49.47-45.2431.2231.22Zc(毫米）417.20417.20417.20417.20417.20由图7观察可知，圆心的定位比较准确而且，表1的结果比较符合实际六、问题的求解6.1 问题3的分析：靶标上圆A,B,C的圆心在同一直线上,用相机摄得其像后,其像也依然要保持直线.所以,可以检验这三个圆心的像坐标是否保持在同一直线上.并且,可以以此衡量模型求解结果的精度.基于不同的检测算子边缘抽取法,分别求解出相应的结果,并利用这种检验方法求得各个结果的精度.根据精度的偏差,可以分析该检验方法的稳定性.6.2 问题3的求解:6.2.1模型检验及误差衡量:一条直线,用相机拍摄后,其像依然是一条直线.在靶标上,圆A,B,C的圆心在同一条直线上.因此, 这三个圆心的像素点坐标,应该保持在同一直线上,即AB的斜率和AC的斜率应该相等.但事实上,圆心的定位有误差,而且像素点坐标系是一个离散的整数点坐标系,因此和计算结果必有偏差.设AB与AC的夹角为,满足关系式,则可利用和的计算结果求得夹角为.的值越小,则说明了模型求解的误差越小,精度越高.用Matlab软件求解得, AB与AC的夹角为,误差非常小.6.2.2稳定性的分析:分别用Matlab软件自带的算子检测函数sobel算子, 梯度(Roberts)算子, prewitt算子, log算子, canny算子抽取边缘像素点.基于不同的边缘抽取法,用基于广义霍夫变换的圆心定位和误差衡量,得到不同的结果,以此进行稳定性的分析.计算结果如下表2 问题3的检验结果(具体结果及图像见附录9.1)：圆心坐标(u,v）圆A圆B圆C圆D圆E(度)sobel算子(325,191)(426,197)(642,213)(584,502)(287,502)0.84Roberts算子(322,188)(423,197)(642,218)(582,502)(287,502)0.70prewitt算子(325,191)(426,197)(642,213)(584,502)(287,502)0.84log算子(325,191)(426,197)(642,213)(590,508)(287,502)0.84canny算子(325,191)(426,197)(642,213)(584,502)(289,497)0.84由结果可知，用不同的边缘抽取法，求得的圆心坐标相差不大，说明该模型的算法抗干扰性较好.其检验结果的误差较小,而且衡量误差的值基本保持不变,说明该检验方法的稳定性较好.七、问题4的求解：7.1 问题分析：精确地确定两部相机的相对位置，这一过程称为系统标定标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置7.2 模型假设：在标定的过程中，采用小孔相机模型，忽略相机的非线性畸变7.3 模型建立：首先，用两部相机同时观测靶标，在标定中，我们可以用单相机标定的方法分别得到，两部相机各自的内外参数，如果外参数分别用、 与、表示，则，表示相机与世界坐标系之间的相对位置，、表示相机与世界坐标系之间的相对位置对任意一点，在世界坐标系、坐标系与坐标系下的非齐次坐标分别为，则 （1）将上式中消去后得到 （2）因此，两部相机之间的几何关系可用以下和表示： （3）7.4 模型求解：第一步：建立世界坐标系以靶标的中心为原点，x,y平面与靶标所在平面重合，背向相机的那一端为z轴的正方向如图所示：图9 相机坐标系和世界坐标系图8 世界坐标系 依据题目所给数据，求出圆心点A，B，C，D，E的世界坐标值，以及AD与CE之间的连线交点F的世界坐标值，第二步：求相机内外参数 首先，根据问题1中建立的模型，求解靶标上五个圆的圆心在该相机像平面的像素坐标因为点F是AD与CE之间的连线交点，所以点F在图像坐标中所对应的点也是与连线的交点因此，利用两直线交点的计算公式可以求得点的图像坐标因为，世界坐标系表示的空间P点坐标与其投影点p的图像平面坐标之间存在着如下关系： （4）其中，是空间点P在相机坐标系下的值，是相机的焦距；（是每一个像素分别在轴方向上的物理尺寸，单位为毫米）；为34矩阵，称为投影矩阵；为相机的内部参数，为相机的外部参数所以，每个靶标上的点及其相应的投影点之间坐标关系可写成：（5）其中，为空间第个点的坐标；为第个点的图像坐标；为投影矩阵的第行列元素式（5）包含三个方程: （6）将式(6)中的第一式除第三式，第二式除第三式分别消去后，可得如下两个关于的线性方程：（）将6 个圆心及其投影点之间的关系可以写成12个关于M矩阵的线性方程，下面用矩阵形式写出这些方程：（8）由式（4）可见, M矩阵乘以任意不为零的常数,不影响与的关系，因此，在式（8）中可以指定1，从而得到关于M矩阵其他元素的12个线性方程，这些未知元素的个数为11个 ，记为11维向量，将式（8）简写成 KU （9）其中，K为式（8）左边1211矩阵；未知的11维向量；U为式（8）右边的12维向量；K、U为已知向量因为1211，所以可以用最小二乘法求出上述线性方程的解为： (10)向量与1构成了所求解的M矩阵但是现在所求解的M矩阵与式（4）所表示的矩阵M相差一个常数因子这一点可以从解方程（8）时指定1中可以看出虽然指定1不影响投影关系，但在分解M矩阵时必须考虑我们将式（4）中的M矩阵与相机内外参数的关系写成 （11）其中，为由式（8）求得的M矩阵的第行的前三个元素组成的行向量；为M矩阵的第行第4列元素；为旋转矩阵R的第行；分别为平移向量的三个向量比较上式两边可得，由于是正交单位矩阵的第三行，因此，由求得再由以下式子可以求得，至此，第一部相机的外部参数已经求解出来了同理，可求得第二部相机的外部参数第三步：求解两部相机的相对位置将代入式（2）这样，两部相机之间的几何关系可用以下和表示：八、模型的分析与推广(1)问题一和问题二的圆心定位模型和算法克服了图像变形及边缘缺损的影响，抗干扰性强(2)该模型的圆心定位算法，基于抽取的边缘点，减少了计算量，有利用计算机实现(3)圆心定位模型的算法精度较高，且计算时间较短，可推广用于快速报靶或无线传感器的图像传送(4)在问题4的求解中，用式子(10)所表示的线性方法求解这些参数时，并没有考虑这些变量间的约束条件，因此，当数据有误差的情况下，计算结果是有误差的，而且误差在个参数间的分配也没按它们之间的约束关系考虑如果对M矩阵进行适当的约束，可提高计算精度参考文献1 马颂德、张正友，计算机视觉，北京：科学出版社，19982 管伟光，体视化技术及其应用，北京：电子工业出版社，19983 张红民、何健鹰，用改进的广义Hough变换获取靶纸图像子像素圆心坐标，计算机与现代化，2003,10:43-50，20034 曹建海、路长厚，基于灰度图像和矢量的圆心定位，光电子激光，2004,15(6):714-718,2004九、附录：附录9.1 问题3根据不同边缘抽取法的圆心定位结果图10 sobel算子边缘抽取及圆心定位结果图11 Roberts算子边缘抽取及圆心定位结果图12 prewitt算子边缘抽取及圆心定位结果图13 log算子边缘抽取及圆心定位结果图14 canny算子边缘抽取及圆心定位结果附录9.2 求解问题1和问题2的Matlab程序建立M文件bianyuan.m，程序如下：。%边缘抽取clear;I=imread(image002.bmp);% 提取图像 b1=edge(I,sobel); %用SOBEL算子进行边缘检测 b2=edge(I,roberts);%用Roberts算子进行边缘检测 b3=edge(I,prewitt); %用prewitt算子进行边缘检测 b4=edge(I,log); %用log算子进行边缘检测 b5=edge(I,canny); %用canny算子进行边缘检测。建立M文件fenge.m，程序如下：。%图像分割x,y=find(b0);A=x,y;xs=sort(x);ys=sort(y);xl=length(xs);yl=length(ys);A1=0,0;A2=0,0;A3=0,0;A4=0,0;A5=0,0;A6=0,0;xx1=135;xx2=281;yy3=375;yy2=186;yy1=137;a1=zeros(xx1,yy1);a2=zeros(xx1,yy2);a3=zeros(xx1,yy3);a4=zeros(xx2,yy1);a5=zeros(xx2,yy2);a6=zeros(xx2,yy3);for i=1:xl; if A(i,1)xx1&A(i,2)yy1; a1(A(i,1),A(i,2)=1; A1=A1;A(i,:); elseif A(i,1)xx1&A(i,2)yy2; a2(A(i,1),A(i,2)=1; A2=A2;A(i,:); elseif A(i,1)xx1&A(i,2)yy3; a3(A(i,1),A(i,2)=1; A3=A3;A(i,:); elseif A(i,1)xx2&A(i,2)yy1; a4(A(i,1),A(i,2)=1; A4=A4;A(i,:); elseif A(i,1)xx2&A(i,2)yy2; a5(A(i,1),A(i,2)=1; A5=A5;A(i,:); elseif A(i,1)xx2&A(i,2)yy3; a6(A(i,1),A(i,2)=1; A6=A6;A(i,:); end end。建立M文件yuanxin.m，程序如下：。%圆心定位%第5个圆的圆心确定a6size=size(a6);aa6=a6;A6l=length(A6);r=round(12*3.78*375/1024);for n=1:A6l; for i=1:a6size(1); for j=1:a6size(2); if (i-A6(n,1).2+(j-A6(n,2).2a6max);a6row=round(mean(a6row);a6col=round(mean(a6col);aa6(a6row,a6col)=1;%第4个圆的圆心确定a4size=size(a4);aa4=a4;A4l=length(A4);for n=1:A4l; for i=1:a4size(1); for j=1:a4size(2); if (i-A4(n,1).2+(j-A4(n,2).2a4max);a4row=round(mean(a4row);a4col=round(mean(a4col);aa4(a4row,a4col)=1;%第3个圆的圆心确定a3size=size(a3);aa3=a3;A3l=length(A3);for n=1:A3l; for i=1:a3size(1); for j=1:a3size(2); if (i-A3(n,1).2+(j-A3(n,2).2a3max);a3row=round(mean(a3row);a3col=round(mean(a3col);aa3(a3row,a3col)=1;%第2个圆的圆心确定a2size=size(a2);aa2=a2;A2l=length(A2);for n=1:A2l; for i=1:a2size(1); for j=1:a2size(2); if (i-A2(n,1).2+(j-A2(n,2).2a2max);a2row=round(mean(a2row);a2col=round(mean(a2col);aa2(a2row,a2col)=1;%第1圆的圆心确定a1size=size(a1);aa1=a1;A1l=length(A1);for n=1:A1l; for i=1:a1size(1); for j=1:a1size(2); if (i-A1(n,1).2+(j-A1(n,2).2a1max);a1row=round(mean(a1row);a1col=round(mean(a1col);aa1(a1row,a1col)=1;%从375*28分辨率转化为1024*768的像素坐标%yxzb=a1col,a1row;a2col,a2row;a3col,a3row;a4col,a4row;a6col,a6row*1024/375;yxzb=a1col,a2col,a3col,a6col,a4col; a1row,a2row,a3row,a6row,a4row *1024/375;yxzb=round(yxzb);yxuvzb=yxzb;ones(1,5)uvxy=1/3.78 0 -512/3.78 0 1/3.78 -384/3.78 0 0 1577/3.78;yxxyzb=uvxy*yxuvzb。建立M文件que2.m，程序如下：。%边缘抽取和圆心定位的画图bianyuan;b=b1;fenge;yuanxin;imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(sobel算子边缘抽取及圆心定位结果);。输入：que2输出：yxuvzb = 325 426 642 584 287 191 197 213 502 502 1 1 1 1 1yxxyzb = -49.4709 -22.7513 34.3915 19.0476 -59.5238 -51.0582 -49.4709 -45.2381 31.2169 31.2169 417.1958 417.1958 417.1958 417.1958 417.1958附录9.3 求解问题3的Matlab程序建立M文件que3.m，程序如下：。%运行程序时,可以根据提示press any key按回车键,分别观看抽取出来的图像%用SOBEL算子进行边缘检测 bianyuan;b=b1;fenge;yuanxin;jianyan; imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(sobel算子边缘抽取及圆心定位结果); pause%用Roberts算子进行边缘检测 bianyuan;b=b2;fenge;yuanxin;jianyan; imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(Roberts算子边缘抽取及圆心定位结果); pause%用prewitt算子进行边缘检测 bianyuan;b=b3;fenge;yuanxin;jianyan; imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(prewitt算子边缘抽取及圆心定位结果);pause%用log算子进行边缘检测bianyuan;b=b4;fenge;yuanxin;jianyan;imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(log算子边缘抽取及圆心定位结果);pause%用canny算子进行边缘检测bianyuan;b=b5;fenge;yuanxin;jianyan;imshow(abs(aa6-1);hold onimshow(abs(aa4-1);imshow(abs(aa3-1);imshow(abs(aa2-1);imshow(abs(aa1-1);title(canny算子边缘抽取及圆心定位结果);pause。输入：que3输出：yxuvzb = 325 426 642 584 287 191 197 213 502 502 1 1 1 1 1yxxyzb = -49.4709 -22.7513 34.3915 19.0476 -59.5238 -51.0582 -49.4709 -45.2381 31.2169 31.2169 417.1958 417.1958 417.1958 417.1958 417.1958ducha = 0.8396yxuvzb = 322 423 642 582 287 188 197 218 502 502 1 1 1 1 1yxxyzb = -50.2646 -23.5450 34.3915 18.5185 -59.5238 -51.8519 -49.4709 -43.9153 31.2169 31.2169 417.1958 417.1958 417.1958 417.1958 417.1958ducha = 0.7004yxuvzb = 325 426 642 584 287 191 197 213 502 502 1 1 1 1 1yxxyzb = -49.4709 -22.7513 34.3915 19.0476 -59.5238 -51.0582 -49.4709 -45.2381 31.2169 31.2169 417.1958 417.1958 417.1958 417.1958 417.1958ducha = 0.8396yxuvzb = 325 426 642 590 287 191 197 213 508 502 1 1 1 1 1yxxyzb = -49.4709 -22.7513 34.3915 20.6349 -59.5238 -51.0582 -49.4709 -45.2381 32.8042 31.2169 417.1958 417.1958 417.1958 417.1958 417.1958ducha = 0.8396yxuvzb = 325 426 642 584 289 191 197 213 502 497 1 1 1 1 1yxxyzb = -49.4709 -22.7513 34.3915 19.0476 -58.9947 -51.0582 -49.4709 -45.2381 31.2169 29.8942 417.1958 417.1958 417.1958 417.1958 417.1958ducha = 0.8396