数码相机定位数学建模论文

PAGE第2页，共46页2008高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：A 我们的参赛报名号为（如果赛区设置报名号的话）：20002028所属学校（请填写完整的全名）：中南大学参赛队员(打印并签名)：1.刘龙2.折巧梅3.黄湘龙指导教师或指导教师组负责人(打印并签名)：刘诚日期：2008年9赛区评阅编号（由赛区组委会评阅前进行编号）：2006高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：PAGE46数码相机定位摘要对于问题一，我们建立了两个模型，求当坐标系原点取在该相机的光学中心时，靶标上圆的圆心在该相机像平面的像坐标。一个是椭圆最小二乘法拟合模型（模型一），这个模型中我们利用了网格纸选点算法和matlab软件函数取点算法两种算法，这两个算法均是通过从像平面中选取六个点的坐标值，代入利用最小二乘法建立的方程组（9）式中，求得了拟合椭圆方程的参数，然后通过式（2）、式（3）的转换公式得到椭圆中心的坐标；另一个模型是直线逼近模型(模型二)，将相片转换为tiff格式后直接调用软件matlab6.0中的imread函数，取得相片中的若干个椭圆区域，然后分别在不同区域内用两条水平线和两条竖直线去与椭圆相切，求出相切交点的坐标，由椭圆对称性可知，两条过对称切点的直线会交于椭圆中心，则椭圆的中心坐标可以表示成.问题二是在问题一的基础上进行求解，将相关数据代入模型中得出了椭圆最小二乘法拟合模型和直线逼近模型的解，结果分别见表4和表9。问题三，在模型一和模型二中，由于包含误差较大样本点在内的所有样本点都参与运算，所以会对椭圆拟合的最后结果产生偏差。针对这种情况，采用随机理论的思想，先随机选取6个点拟合椭圆，然后计算与此椭圆匹配的所有样本点个数。重复此过程一定次数（一般选取100~200次），匹配样本点多的椭圆即为最优椭圆，由此构造了一种快速准确剔除误差较大样本点的改进椭圆拟合随机化算法，并通过实例验证了算法的精度和稳定性，其中精度检验时最大误差为0.73%，稳定性检验时最大误差为3.7%，说明模型能够有效地处理包含有较大比例误差点的样本空间，拟合出具有高精度的椭圆，并且该算法的速度能够满足实时性的要求。问题四属于双目定位问题范畴，可由常规的摄像机标定问题进行反推求解，常规的摄像机标定问题是将像素坐标系信息（如相片）转换到世界坐标系信息（如实物）的三维坐标，而这里是将世界坐标系信息（靶标上的圆心）向像素坐标系信息转换，从而测定两部固定的数码相机的相对位置。我们通过已知的实物和相片信息然后利用空间坐标系转换公式求解得到由旋转距阵和自由平移向量组成的两部相机的外部参数矩阵，进而得到两部数码相机的光心坐标变换模型（见式（29）），由此即得两相机的相对位置。关键词：最小二乘法、图像网格处理、取色器、标定、旋转距阵一问题重述数码相机定位在交通监管（电子警察）等方面有广泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。标定的一种做法是：在一块平板上画若干个点，同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，如图1所示，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。图1靶标上圆的像我们设计靶标如下，取1个边长为的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，为半径作圆。以AC边上距离A点处的B为圆心，为半径作圆，如图2所示。图2靶标示意图用一位置固定的数码相机摄得其像，如图3所示。图3靶标的像要求：建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标,该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为；设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；建立用此靶标给出两部固定相机相对位置的数学模型和方法。二问题分析2.1问题背景分析数码相机在计算机视觉中的应用逐渐普及和深入，在应用中数码相机的标定是相当重要的。计算机视觉是一门新兴的学科。随着计算机硬件、软件、图像采集、处理技术的迅速发展，计算机视觉的理论和技术已被广泛地应用于医学图像处理、机器人技术、文字识别、工业检侧、军事侦察、地理勘察和现场测量等。计算机视觉在多种测量中的应用是一种定量分析系统，有确定的精度要求。一般运用于现场不可到达、现场复杂不便直接测量、保留现场状态以便复测、补测和事后重构现场等场合。通常人们把摄像机和胶片相机作为获取原始图像的主要设备，随着数码相机的出现和发展，在许多场合人们已逐渐采用数码相机作为主要的图像采集设备。这主要是由于数码相机有如下几个特点:有较多像素的CCU;可以根据测量的需要选择不同的像素大小;可直接与计算机进行通讯;可选定多种焦距进行定焦距拍摄。上述特点，给图像采集和图像处理带来了许多方便之处。在应用计算机视觉的测量系统中，无论是采用摄像机还是数码相机作为原始图像的采集设备，都必须进行系统标定。系统标定是应用计算机视觉进行现场测量和图像处理的第一步，也是关键的一步。标定矩阵的精度直接影响到最终的测量精度。国内外许多学者对此做了大量的研究，提出了多种利用摄像机作为输入设备的标定方法。2.2问题要求分析对于问题一，题目要求确定靶标上圆的圆心在该相机像平面的像坐标，这里坐标原点是取在该相机的光学中心，我们若以光轴所在直线为轴，则所求靶标上圆的圆心在相机像平面上的纵坐标是相同的（在这里为该相机的像距，是一个固定值），所以只需求各和坐标即可，由实物与像之间的关系，在只考虑几乎理想状态下，我们可以试着确定相片上椭圆的中心坐标，并且认为靶标上圆的圆心在该相机像平面的像坐标应该是和椭圆的中心坐标很接近。问题二是在问题一建立的模型上赋予具体的数值进行求解，因此只需将题目中所给数据代入第一问的模型中即可有解决，值得注意的是在理想光学成像系统下焦距即是像距。问题三是要求设计一种方法来给第一问设计的模型进行检验，同时还要求这种方法具有一定的精度和稳定性，并能计算出来。问题的难点在于这里没有一个标准又不能进行真实模拟，鉴于此，我们决定采用计算机仿真的思想，即给出一个确定解同时产生一系列随机数，用该方法求得的一组值来达到分析的目的。对于问题四，由于这是属于双目定位问题。根据相机成像的原理，我们利用空间坐标系的相互转换就应该可以实现在靶标给出的条件下确定两部固定相机的相对位置。三模型的建立及求解数码相机图像拍摄实际上是一个光学成像过程，而且光学成像的理论模型是针孔模型。我们将此成像过程分为3个步骤、隶属于4个坐标系。这四个坐标系分别为:(1)世界坐标系—根据自然环境所选定的坐标系，坐标值用表示。(2)图像坐标系—坐标原点在图像平面的中心，X轴、Y轴分别为平行于图像平面的两条垂直边，坐标值用表示。(3)像素坐标系—坐标原点在图像平面的左上角，X轴、Y轴分别平行于图像坐标系的X轴和Y轴，坐标值用来表示，且为离散的整数值。(4)光心坐标系—以相机的光心为坐标原点，X轴、Y轴分别平行于图像坐标系的X轴和Y轴，相机的光轴为Z轴，坐标值用表示。由实物到相片的形成经过了3个步骤，是分别将世界坐标系中的信息转换到光心坐标系，再由光心坐标系转换到图像坐标系，最后由图像坐标系转换到像素坐标系。3.1问题一的模型3.1.1基本假设考虑最小二乘法的要求，我们假设从图像坐标系中选取的点的坐标值的随机误差服从正态分布。图3中的椭圆都是标准椭圆，不存在畸变等因素，属于理想图形。相片没有损坏，保存完好，图像清晰3.1.2符号说明——椭圆中心位置——椭圆长轴和短轴a——椭圆长轴的转角3.1.3模型的准备1）椭圆的表示方法在二维平面坐标系中，椭圆一般可以用2种形式来表示。一种是利用圆锥曲线方程的代数形式表示，即:（1）另外一种更直观的方式是用平面坐标系的几何参数表示，即椭圆中心位置，长轴和短轴，长轴的转角a。二维平面里的任意椭圆都可以用这5个参数唯一确定，参数的几何意义如图4所示。图4二维平面椭圆的表示两种表示形式的参数可用式(2)～式(6)转换。(2)(3)(4)(5)(6)2）坐标转换分析根据这种转换关系，可以先选取图3中椭圆中点的坐标值，求得等式（1）中参数，再通过转换，得到图3中各个椭圆中点在像素坐标系中的坐标值，最后由坐标平移求出各个椭圆中心在图像坐标系的的坐标值。位置如图5图5图像坐标系下各个点位置图3.1.4模型的建立建立两个模型，模型一中用到最小二乘法来拟合椭圆函数，将得到的拟合函数的参数，用式（2）、（3）转换得到题中要求坐标系情况下椭圆圆心坐标。在选取点拟合过程中，考虑对精度的要求，给出在精度要求不高情况下，用网格取点法和高精度要求中的计算机像素取点法两个算法；模型二中根据标准椭圆对称性，用matlab软件中图像处理库函数建立模型，提高运行速度。模型一椭圆最小二乘法拟合最小二乘法椭圆拟合是较常用的椭圆拟合方法。最小二乘法是在随机误差为正态分布时，由最大似然法推出的一个最优估计技术，它可使测量误差的平方和最小，因此也被视为从一组测量值中求出一组未知量的最可信赖的方法之一.最小二乘技术主要是寻找参数集合，从而最小化数据点与椭圆之间的距离度量.这里的距离度量常见的有几何距离和代数距离。几何距离表示某点到曲线最近点的距离。平面内某点到方程所代表曲线的代数距离就是。以下是以代数距离作为距离度量介绍最小二乘法。假设一般形式的椭圆方程如式(1)所示，为了避免零解，并将解的任何整数倍都视为对同一椭圆的表述，对参数做一些限制，约束条件设。显然，直接应用上述方程对边缘检测后的离散点进行最小二乘处理，就可以得到方程中的各系数，也即，求目标函数（7）的最小值来确定各系数。再由极值原理，欲值为最小，必有（8）由此建立如下数学规划模型：我们用matlab6.0编程并用消元法求解上述线性方程组，最终求得方程系数的值。模型二 直线逼近模型每个椭圆用两条水平直线和两条竖直直线逼近椭圆，使得四条直线与椭圆相切，得出每个椭圆水平、竖直边界点。Imread函数对处理的函数有特定的要求，在用imread函数调用图片之前，先用photoshop软件将图3图片格式转换为tiff格式并无任何压缩限制，然后用imread函数取得图片中椭圆区域。再在图中划分出每个椭圆的边界，如图6。图6椭圆边界划分对每个区域用imread函数取得的图片数据中，分区域0HFJ,HIGF,ILKG,DKQN,JEMP在每个区域用水平线和竖直线去与椭圆相切，可以求出相切交点的坐标（以椭圆D为例）如图7。图7椭圆D相切交点位置其中由椭圆对称性可知，直线与直线交于一点且过椭圆中心，则椭圆中心可以用坐标表示，得到如下数学模型（在图像坐标系下）：（10）3.1.5模型的求解模型一的求解模型中首先要解决的问题是得到六个点，用这六个点来算出等式（7）中参数六个点的值。对于选取这六个点，给出一下两个算法；网格纸选点算法用一个网格纸放在图形上面，建立一个坐标系，确定每一个网格对应像素，再取网格顶点正好落在每个椭圆边界处，当取不到足够的点时，可以移动网格纸位置，使得足够多的点落在网格顶点，如图8所示图8网格纸及点选取 算法具体实现步骤：Step1放置网格纸，使尽可能多的网格顶点落在椭圆边缘Step2根据相机分辨率，设置网格坐标刻度Step3选取6个点的坐标值Step4代入式（8）中，得到关于系数的线性方程组Step5用消元法求得每个椭圆对应的值Step6通过式（2），（3）求出每个椭圆圆心在下的坐标值Step7坐标变换，得到在像素坐标系下椭圆圆心坐标值2）matlab软件函数取点算法用photoshop软件将图3格式转换为tiff格式，用matlab图形工具箱中的imread函数取得图片中椭圆区域，再在图片区域找出椭圆边界上6个点的坐标值，用消元法求出最终结果 算法实现步骤：Step1用photoshop软件将图3格式转为tiff格式并无需任何压缩Step2用imread函数调用转换后图片Step3从得到的图片信息中提取出每个椭圆边界点（代码见附录1）Step4从每个椭圆边界点中取出6个点的坐标值（代码见附录2）Step5用matlab软件求出式（3）的方程组形式，并用消元法求得系数值Step6通过式（2），（3）和坐标变换，最终得到在像素坐标系下椭圆圆心坐标值[注]：消元法是以函数形式给出，代码见附录3模型二的求解在imread图像数据中，可以直接得到的值，代入式（10）便得到椭圆D的圆心在图像坐标系下的坐标。相继也可以求出靶标上其他椭圆圆心在图像坐标系下的坐标值。算法具体实现步骤：Step1用photoshop软件将图3格式转为tiff格式并无需任何压缩Step2用imread函数调用转换后图片Step3将图片按椭圆位置分块Step4在每一块中用matlab找出椭圆水平和竖直方向最大值和最小值Step5用式（10）求出每个椭圆中心在图2坐标系下的坐标值Step6坐标变换，最终得到在像素坐标系下椭圆圆心坐标值3.1.5模型的评价对于模型一，首先在图像坐标系中对图像进行图像网格处理，得到每个椭圆上的一些点的坐标值；然后根据椭圆的表达式，用最小二乘法求出每个椭圆中心在图像坐标系中对应的坐标值。给出取点的两种方法：网格取点与计算机函数取点。其中网格取点算法可用于对结果精度要求不高，或是没有高清晰摄像设备的情况下，在实际中，这种取点方法能方便快速的得到一些基本的数据，给生产工作带来便利，是一种辅助其他算法的工具，也可以对其他的算法结果正确性做简单的分析，但这种方法由于精度不高，误差大；计算机取点算法能精确的取到图片要求区域的边界坐标值，能用取到的精确坐标值进行最小二乘法拟合，利用这种取点算法得到的解精度很高，且当图形有略微畸变时，可以避开畸变区域点的选取，最终得到理想的结果，但这种方法需要高清晰度取像设备，不宜在室外作业，可以在实验室的条件下做实验数据分析。 对于模型二，解决问题的前提是利用椭圆的对称性，快速的得到椭圆上四个极值点的位置，再结合几何图形知识得出最终的结果，在处理具有对称性的图形时，能非常快速的得到精确解，并且在程序执行过程中花去的代价很低，但当图形有稍微的畸变，就容易产生很大的误差。模型一和模型二都有各自的适用范围，可以根据不同的情况需求选择不同的模型或算法来解决问题，也可以将不同的模型或算法结合起来解决待解决的问题。3.2问题二的解答由于问题二中给出了具体的数值且问题模型并没有改变，显然问题二的求解必须是要在问题一模型下进行，所以只需将具体的数值代入问题一的模型即可。利用matlab6.0的强大数值计算能力，我们求得了如下解。3.2.1利用模型一中网格取点算法求解这里只考虑用网格选点法，得到从图3中取得的点的坐标如下表1表1图像坐标系取点坐标值椭圆x(mm)y(mm)A190.204439.6019275.170344.5522375.170354.4526490.204459.4029595.215744.5522695.215754.4526B1105.238444.55222120.272544.55223105.238459.40294110.249861.87805122.778249.50246117.766859.4029C1160.363354.45262160.363359.40293177.903149.50244175.397447.02735165.374747.02736177.903159.4029D1150.3406123.75602145.3293133.65653160.3633123.75604162.8690133.65655160.3633136.13166152.8463123.7560E165.1476133.6565270.1590123.7560375.1703141.0818480.1817138.6067585.1930128.7062682.6873123.7560在matlab6.0中拟合出每个椭圆方程的参数值，见表2表2每个椭圆代数形式参数的拟合值椭圆拟合系数坐标中的点A-1.3953-0.3494-99.0371158.7387-0.7314B-1.3677-1.8741-159.2396353.1355-0.5491C-2.1063-3.0887-221.2915686.4699-0.4062D0.22701.7202-337.1185405.6205-0.0189E-0.01500.2950-148.024181.2751-0.0887由表2拟合的每个椭圆方程的参数值代入（2）（3），得到各椭圆圆心在图像坐标系下的坐标如下：表3图像坐标系下拟合各个椭圆中心坐标椭圆(mm)ABCDEX86.9180113.2818167.5217154.027575.0327Y56.609552.148053.0060132.0610137.8438表3得到的是在图像坐标系中的坐标，如图2，图片中心相对左上角的坐标为（135.3066，101.4799）(单位为毫米)进行坐标变换，用图像坐标系下的各个点X坐标值值减去135.4497（mm），得到在像素坐标系下各个点的x值，用101.5873减去在下的各个点Y坐标值，得到在像素坐标系下各个点的y值，如表4表4在像素坐标系下各个点的坐标值椭圆（mm）A(mm)B(mm)C(mm)D(mm)E(mm)x-48.3886-22.024832.215118.7209-60.2739y50.870449.331946.4739-30.5811-32.36393.2.2利用模型一中计算机函数取点算法求解 用函数在每个椭圆上取出六个点坐标值如表5表5各点取值x（像素）y（像素）椭圆A128320223641823310229434215353262306284206椭圆B138419223842093461206445422054302356462202椭圆C160623026431763652176460623256232486675202椭圆D157147526064753612515458247155624816585534椭圆E160422026062043607233460522956042146669232在matlab6.0中拟合出每个椭圆方程的参数值，见表6表6每个椭圆代数形式参数的拟合值椭圆拟合系数坐标中的点A-2.04070.8922-243.6945340.9528-1.0293B-1.8508-0.2235-470.2341920.9710-2.2089C-0.00330.0010-0.56661.6739-0.0003D-1.0074-0.2193-684.4368815.3615-0.1040E-0.00300.0004-0.57781.6644-0.0011618以一个单位像素作为刻度，图3中的图像可以划分为个单元格，用matlab库函数中的imread函数来取得每个单元格的值，再用问题一中的算法可以取得每个椭圆心坐标值，最后转换到以光学中心为原点的坐标值，见表7表7图像坐标系下各椭圆圆心坐标值A(mm)B(mm)C(mm)D(mm)E(mm)x-50.0000-23.544933.862518.7831-60.1851y51.455049.603245.1058-31.4815-31.21693.2.3利用模型二将图片转换为tiff图片格式，用取色器在转换后图片上找出椭圆边界点的RGB值，可以得到表8值表8tiff格式图片上RGB值RGB椭圆边界内和椭圆上可取值1515151241414133434344282828522222261212127444椭圆边界和边界外可取值12342342342227227227322122122142152152155204204204在表8中，椭圆内部的RGB值都为，椭圆外部的RGB值都为，其他的值都分布在椭圆边界附近。以一个单位像素作为刻度，图3中的图像可以划分为个单元格，用matlab库函数中的imread函数来取得每个单元格的值，再用问题一中的算法可以取得每个椭圆心坐标值，最后转换到以光学中心为原点的坐标值（代码见附录5），见表9表9光心坐标系下各椭圆圆心坐标值ABCDEx-50.0000-23.544933.862518.7831-60.1851y51.455049.603245.1058-31.4815-31.21693.3问题三的解答直接最小二乘法并未考虑各样本点误差的差异，它假设样本数据是零均值的，且有共同的协方差阵，这与实际情况不符，因而导致参数的有偏估计，算法结果往往并不令人满意。基于这种考虑，我们研究了一种改进的椭圆拟合的算法。这个算法主要是用来解决在取得的图像中有畸变的情况。3.3.1算法建立——基于最小二乘法和随机原理，可得到一种具有较高抗干扰能力的椭圆拟合算法。算法原理如下：1)在所有样本点(已编号)中随机选取6个样本点；2)利用最小二乘法求解椭圆参数；3)遍历所有样本点，求取各个点到已得到的椭圆之间的距离，最小二乘技术主要是寻找参数集合,从而最小化数据点与椭圆之间的距离度量。这里的距离度量常见的有几何距离和代数距离。如果小于某个自己定义点阈值，则称该样本点为匹配点，在数组a记录该样本点的编号，遍历完毕求取对于该拟合椭圆的匹配点的总个数b；4)比较匹配点总个数c与匹配点最大值d，当前者大于后者，将椭圆参数和记录匹配点编号的数组M保存下来，分别拷贝至数组p和数组q，最后将b赋值给d；5)循环执行步骤1)～步骤4)一定次数(根据运行时间、需要结果的准确度以及样本点总个数适当定义)，最后在p保存了最优椭圆参数，在数组q保存了在所有样本点中匹配点的编号，也就相应的可以得到不匹配点的编号。3.3.2算法说明算法第1)步中随机选取6个点拟合椭圆，之所以选取6个点是因为在离散数据中，往往5个点不能得到一个椭圆，本文选取6个点，用基于代数距离的最小二乘法拟合椭圆。之所以随机选取样本点进行椭圆拟合是考虑到算法的实时性和准确性。如果考虑所有的可能性，譬如总样本空间有300个点，选取6个不同的样本点一共有种可能，遍历所有可能性在时间上是不允许的，而且也不必要。要加速计算，只需选择其中的一部分子集。在样本性质并不清楚的情况下，使用随机方法抽取是一个很好的选择。算法流程图如图9所示。（代码见附录6）开始开始在所有样本点E随机选取6个点循环次数是否已到利用最小二乘法求解椭圆参数遍历所有样本点判断点到椭圆的距离是否符合要求得到num_inliers和inliers_indexbest_elipse_par=椭圆参数max_inliers_index=inliers_indexmax_inliers=num_inliers结束num_inliers>max_inliers是否否是图9算法流程图3.3.3算法对问题一模型的检验现取题二中的椭圆A作为分析对象(取点代码见附录7)，样本总体上有233个样本点，则有227个样本点来计算点到拟合椭圆的距离，遍历100次。为得到精确的结果，阈值应越小越好，又因为问题二中的图片上椭圆比较多标准，阈值在这里取为2个像素单位，在matlab中运行代码，对数值进行遍历，最终得到p中保存的最优椭圆参数，见表10表10p中保存的最优椭圆参数BCDEF数值-2.54000.9736-204.7289418.6218-6.2642将得到的数值ABCDEx-50.0000-23.544933.862518.7831-60.1851y51.455049.603245.1058-31.4815-31.2169通过与问题二的数值比较发现，问题一的假设是合理的，且用问题一种建立的模型，对比较标准的图形可以得到很精确的3.3.4算法评价 在直接最小二乘法的情形中，并没有考虑到畸变给结果带来的误差。当给定图片存在肉眼看不见的畸变或杂质点，如果用直接最小二乘法可能会产生比较小的，或者是不易察觉的误差，分析但这些误差可能会在生产工作中慢慢的积累，最终可能产生破坏性的影响。例如在眼睛治疗中，需要对眼膜有非常高的精度要求，随机化方法可以很好的消除因设备、杂质点等带来的误差 随机化算法是通过一定量的遍历来消除畸变、杂质点等影响，得到想要的到的结果，有时这种算法能达到亚像素精度。3.3.5算法精度检验1)精度概念分析当给定一个误差范围，通过拟合值与真实值之间的比较，得到误差平法和、方差等统计量来对结果进行检验2)检验方法描述给定一个图形，并知道图形所有点坐标值、参数值等；取到每个点的坐标值，并储存；对每一个图形上的点用matlab中normrnd函数生成正态（高斯）分布的随机数，将这个随机数加到图形对应点中，从而得到一个和原图形有偏差的新图形；以新图形上的点坐标作为已知值，用随机化算法找到最接近原图形的参数值用新得到的参数值与已知参数值做统计量比较，检验随机化算法拟合精度。3)检验实例假设a,b值拟合误差超过1%，认为算法不合格，否则，证明此算法在拟合方面有很高的精度给定一标准椭圆，中点在原点，a=50,b=30；取到每个点的值，以1为一个单位，总共有200个点x(i),y(i)(i=1,2,…,200)；用normrnd（0，0.1）生成400个白噪声数值nx(i,1),ny(i,2)(i=1,2,…,200)，组成一个白噪声序列，每两个用来控制一个点在x轴和y轴上移动的值;将对应的点在白噪声nx,ny的影响下产生新的位置值x1(i),y1(i);将数据x1(i),y1(i)用随机化算法，来求得a1,b1的值；对拟合值a1,b1进行误差分析其中将(x1(i),y1(i))点和原椭圆画在一个图中，如图10图10随机产生点与原椭圆用随机化算法计算得a1,b1值分别为50.1205，29.8781重新生成白噪声序列，再重复进行多次的拟合，最终得到表11中数据（代码见附录8）表11拟合值12345678910a150.120549.999850.000050.200149.973149.739150.031749.690250.100049.8973b129.878129.917330.010430.217830.078930.108330.110530.115730.030630.0338分别计算出a1,b1样本均值有a1,b1的样本方差计算每组数据的相对误差如表12表12相对误差12345678910a10.00240.00000.000000.00400.00050.00520.00060.00620.00200.0021b10.00410.00280.00030.00730.00260.00360.00370.00390.00100.0011从先对误差表中可得a1最大误差为0.4%，b1最大误差为0.73%，都没超过1%，证明此算法有很高的精度。由于产生随机序列的值可能取到很大，则会出现较大奇异点，使得在其他模型的求解过程中产生很大的误差，而随机化算法可以把这些有较大奇异点的地方，通过筛选把奇异点排除在正确解的范围内。能够有效地处理包含有较大比例误差点的样本空间，拟合出具有高精度的椭圆，并且算法的速度能够满足实时性的要求。3.3.6算法稳定性检验稳定性概念分析在一组影响变量a的因素中，当某一个因素发生一个小的变化，会使得变量发生很大的变化，则变量不稳定，如果在一定的范围内，影响因素的变化不能影响变量发生很大的变化，则这个变量是稳定的。在随机化算法中，主要解决的问题是通过循环排除畸变或奇异点等等影响，如果模型中出现较大的畸变或出现大的奇异点，而随机化算法依然能得到拟合效果很好的拟合情况。否则认为随机化算法稳定性不够。2)检验方法描述给定一个图形，并指导图形在坐标系下的所有情况和图形本身参数值；给图形制造一定范围的畸变或者奇异点；用新图形的坐标值拟合求解图形的参数值；进行参数检验，确定随机算法是否具有较高的稳定性3)检验实例假设c,d值拟合误差超过5%，认为算法不合格，否则，证明此算法在拟合方面有很高的精度给定一标准椭圆，中点在原点，c=50,d=30；取到每个点的值，以1为一个单位，总共有200个点x(i),y(i)(i=1,2,…,200)；从y=40出小区图形右边部分，用normrnd（0，0.1）生成400个白噪声数值nx(i,1),ny(i,2)(i=1,2,…,200)，组成一个白噪声序列，每两个用来控制一个点在x轴和y轴上移动的值;将对应的点在白噪声nx,ny的影响下产生新的位置值x1(i),y1(i);将数据x1(i),y1(i)用随机化算法，来求得c1,d1的值；对拟合值a1,b1进行误差分析用随机化算法计算，遍历100次，得a1,b1值分别为48.1205，29.5781重新生成白噪声序列，再重复进行多次的拟合，最终得到表13中数据表13拟合值12345678910c147.920148.299848.000049.100148.773148.837148.233747.890248.104047.8633d128.878129.735630.710429.217829.578930.207331.110530.214731.041631.0778分别计算出c1,d1样本均值有a1,b1的样本方差计算每组数据的相对误差如表14表14相对误差12345678910a10.00240.00000.000000.00400.00050.00520.00060.00620.00200.0021b10.03740.00880.02370.02610.01400.00690.03700.00720.03470.0359从先对误差表中可得a1最大误差为1.8%，b1最大误差为3.7%，都为超过5%，证明此算法有很好的稳定性。在较大比例误差点的样本空间并不会给最终的解的值带来很大的误差，依然能保持解的稳定性。模型稳定性很高。3.3.7算法评价通过以上验证了算法能够有效地处理包含有较大比例误差点的样本空间，拟合出具有高精度的椭圆，并且算法的速度能够满足实时性的要求。能够有效地处理包含有较大比例误差点的样本空间，拟合出具有高精度的椭圆，并且算法的速度能够满足实时性的要求。能很好的将引起较大误差的点排除，并寻找最优解或者近优解，将误差控制在很小的范围内。3.4问题四的模型——几何模型3.4.1基本假设=1\*GB3①设两部相机的焦距相等，各内部参数也相同②两部数码相机的坐标平面是选在光心坐标平面=3\*GB3③不考虑相机透镜的畸变，采用较好的光学镜头3.4.2符号说明——世界坐标系中一空间点的坐标——光心坐标系中一空间点的坐标——对应图像坐标系中点的坐标——图像坐标系原点在像素坐标系中的坐标——分别是像素坐标系在方向和方向相邻像素间的距离——拍摄相机的焦距——旋转矩阵——位移向量3.4.3模型的建立如图11所示,将物体所在的空间设定为世界坐标系,而相机拍摄的图像是电荷藕合器件CCD面转化的,所以,可以认为CCD平面为图像平面,定义图像平面为像素坐标系,坐标原点在CCD图像平面的左上角,u、v轴平行于CCD平面的两边,用表示,该坐标值为整数值.图11数码相机成像坐标系相机拍摄的原理是将坐标系中的物体转化成像素坐标系中的图像,转换过程可以通过建立两个中间坐标系来实现.3个步骤是分别将世界坐标系中的信息转换到光心坐标系，再由光心坐标系转换到图像坐标系，最后由图像坐标系转换到像素坐标系。（11）其中是光心坐标系中一空间点P的坐标，是对应图像坐标系中P点的坐标，f是拍摄相机的焦距。由世界坐标系到光心坐标系的转换关系为（12）其中，R为旋转矩阵，t为位移向量，元素为0的列向量。由图像坐标系到像素坐标系的转换关系为（13）其中，是图像坐标系原点在像素坐标系中的坐标，分别是像素坐标系在X方向和Y方向相邻像素间的距离。将式（11），（12）代人式（13）我们得到如下模型：（14）可将式（14）简化，得到最终模型为：（15）为式(4)右边第1项即相机的内部参数矩阵，为式(14)右边第2项即相机的外部参数，为投影矩阵。3.4.4上述是针孔模型下的成像系统的基本标定原理，相机的内、外参数都作为未知内容进行求解和标定的。在利用数码相机进行采集图像时，可以通过多种途径求解出数码相机的内部参数，然后再根据逆向投影的方法来求解旋转矩阵和位移向量。可将式(15)改写成（1６）从式(6)可见内部参数矩阵为满秩矩阵，令:，用左乘式(16)两边并整理得(17)令并可以进一步将式（17）改写成(18)其中，是世界坐标系中的第点坐标，是由对应点的图像像素坐标的计算值。将其展开可得到如下3个方程(19)将式(19)中的第3式分别代入上两式消去。可得到如下两个方程(20)若在世界坐标系中取个点，将会产生个方程，用矩阵形式写出这些方程(21)从式(21)可见，若已知世界坐标中个点的坐标，同时也知道各对应点在像素坐标系中的坐标，将式(14)两边同除以或取不影响方程的求解。式(21)可简写成(22)其中，分别为式(21)等号左侧向量和右侧向量除以后的向量，可用最小二乘求解上述线性方程组，其解为(23)则数码相机的外部参数矩阵与其有如下关系(24)因是正交矩阵，所以有(25)式(24)称为标定矩阵的隐式解，要求其显式解，多种旋转矩阵中的一种,如取（26）则有分别求出，世界坐标系相对光心坐标系的位置为3.4.5这里我们可以选取世界坐标系中靶标上的五个圆心点，并得到各对应成像点分别在这两部相机像素坐标系中的坐标，这些坐标值可由第一问的模型求解得来.利用这些坐标值再代入上述几何模型我们可得到两个相机的外部参数矩阵和代入并（12）式得：（27）（28）其中和分别表示在两部相机下光心坐标系下空间点的坐标，显然和均为可逆矩阵，所以将乘以（28）式等号两边并代入（27）式得两相机的相对位置关系，由此得到如下光心坐标变换模型：（29）3.4.7模型的评价对于该模型，分别先得出两部相机的外部参数矩阵，进而得出靶标在两部相机对应光心坐标系下的坐标转换关系，模型相对简单，且易于理解，但是并没有明确得出这两部相机的相对位置关系。因此我们可在模型的推广中进行讨论，通过坐标系的转换关系分别得出两部相机在世界坐标系下的坐标，然后由此推出其相对位置关系。4.模型的推广 问题一中，我们考虑的是实物的中在图像坐标系中的像是像图中心的情况，而在实际情况，实物的中心在图像坐标系中的像并不一定在像图中心，如果有特殊的要求，或者要求精度很高时，就要找到具体的中心对照位置，这就需要建立实物与像图位置的几何位置关系，通过求转移矩阵和旋转矩阵，逐步求出实物中心在像图中对应的中心位置。上面建立的模型都是基于针孔模型下建立的线形模型，实际中有时会遇到非针孔模型，即要建立非线性模型来对模型进行求解，可以用遗传算法来解决这类问题，遗传算法在摄像机标定方面已取得一定成效，但主要的研究工作集中在单个摄像机的参数标量采用遗传算法标定摄像机参数在收敛性精度和鲁棒性方面具有很好的性能但由于其摄像机模型中没有包含透镜畸变因素如果用于三维测量则难以保证精度，但是用于二维的条件下是适用。在问题四中，可以建立仅用一台数码相机，围绕目标进行多角度自由拍摄获得多帧图像的模型，采用一套具有身份唯一性的标记点实现多图像间点的稳定匹配；在此基础上根据多视图几何约束，精确计算各次拍摄时的相机姿态和位置；进而求解目标点的三维位置信息，适用于大型物体多视角测量的拼合、某些情况下的装配检查等，也可直接用于较规则构型物体的三维数据获取与CAD模型重建。5.参考文献[1]闫蓓王斌李媛，基于最小二乘法的椭圆拟合改进算法，北京航天航空大学学报，第34卷第三期：2008.3。[2]姜启源谢金星叶俊著，数学模型（第三版），北京：高等教育出版社，2003。[3]姜大志孙闵等著，数码相机标定方法研究，南京航天航空大学学报，第33卷第1期：P57～58，2001.2。[4]邹宏伟著，PhotoshopCS案例教程，北京：中国石化出版社，2007。[5]\o"主题标题。"MATLAB知识/1304141710/archive/2007/11/08/400066068.aspx，2008-9-22附录问题一附录：附录1从图像中选取椭圆边界坐标值m1=imread('p.tif');fori=1:768 forj=1:1024 ifm1(i,j,1)<70 p(j,i)=m1(i,j,1); else p(j,i)=0; end endendfori=2:1023 forj=2:767 ifp(i,j)~=0 ifp(i+1,j)==0||p(i-1,j)==0||p(i,j+1)==0||p(i,j-1)==0 n1(i,j)=1; else n1(i,j)=0; end else n1(i,j)=0; end endendk1=0;fori=2:1023 forj=2:767 ifn1(i,j)==1 k1=k1+1; a1(k1)=i; b1(k1)=j; end endend用0和1表示出5个椭圆形状附录2生成随机的六个数c1=randperm(233);%用于生成随机的点fori=1:6%取前6个作为随机点 x(i)=a1(c1(i)); y(i)=b1(c1(i));end附录3消元法函数代码functionx=DelGauss(a,b)%Gauss消去法[n,m]=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);fork=1:n-1fori=k+1:nifa(k,k)==0returnendm=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);fork=n:-1:1%回代forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);end问题二附录：附录4用方格取点算法带入数据求解m=[18.00008.000015.00009.000015.000011.000018.000012.000019.00009.000019.000011.000021.00009.000024.00009.000021.000012.000022.000012.500024.500010.000023.500012.000032.000011.000032.000012.000035.500010.000035.00009.500033.00009.500035.500012.000030.000025.000029.000027.000032.000025.000032.500027.000032.000027.500030.500025.000013.000027.000014.000025.000015.000028.500016.000028.000017.000026.000016.500025.0000];fori=1:30 m(i,1)=m(i,1)*1024/(54*3.784); m(i,2)=m(i,2)*768/(41*3.784);end fori=1:6 x(i)=m(i,1); y(i)=m(i,2);enda(1)=0;b(1)=0;c(1)=0;d(1)=0;e(1)=0;f(1)=0;g(1)=0;h(1)=0;j(1)=0;k(1)=0;z(1)=0;u(1)=0;v(1)=0;fori=1:6 a(i+1)=a(i)+x(i)^2*y(i)^2; b(i+1)=b(i)+x(i)^3*y(i); c(i+1)=c(i)+x(i)*y(i)^3; d(i+1)=d(i)+x(i)^2*y(i); e(i+1)=e(i)+x(i)*y(i)^2; f(i+1)=f(i)+x(i)*y(i); g(i+1)=g(i)+y(i)^4; h(i+1)=h(i)+y(i)^3; j(i+1)=j(i)+x(i)^2; k(i+1)=k(i)+x(i)^3; z(i+1)=z(i)+y(i)^2; u(i+1)=u(i)+x(i); v(i+1)=v(i)+y(i);endp(1)=a(7);p(2)=b(7);p(3)=c(7);p(4)=d(7);p(5)=e(7);p(6)=f(7);p(7)=g(7);p(8)=h(7);p(9)=j(7);p(10)=k(7);p(11)=1;p(12)=z(7);p(13)=u(7);p(14)=v(7);q(1,1)=p(1);q(1,2)=p(3);q(1,3)=p(4);q(1,4)=p(5);q(1,5)=p(6);w(1)=-p(2);q(2,1)=p(3);q(2,2)=p(7);q(2,3)=p(5);q(2,4)=p(8);q(2,5)=p(12);w(2)=-p(1);q(3,1)=p(4);q(3,2)=p(5);q(3,3)=p(9);q(3,4)=p(6);q(3,5)=p(13);w(3)=-p(10);q(4,1)=p(5);q(4,2)=p(8);q(4,3)=p(6);q(4,4)=p(12);q(4,5)=p(14);w(4)=-p(4);q(5,1)=p(6);q(5,2)=p(12);q(5,3)=p(13);q(5,4)=p(14);q(5,5)=1;w(5)=-p(9);x=DelGauss(q,w);s(1)=(x(1)*x(4)-2*x(2)*x(3))/(x(2)*4-x(1)^2);t(1)=(x(1)*x(3)-2*x(4))/(4*x(2)-x(1)^2);m=[18.00008.000015.00009.000015.000011.000018.000012.000019.00009.000019.000011.000021.00009.000024.00009.000021.000012.000022.000012.500024.500010.000023.500012.000032.000011.000032.000012.000035.500010.000035.00009.500033.00009.500035.500012.000030.000025.000029.000027.000032.000025.000032.500027.000032.000027.500030.500025.000013.000027.000014.000025.000015.000028.500016.000028.000017.000026.000016.500025.0000];fori=1:30 m(i,1)=m(i,1)*1024/(54*3.784); m(i,2)=m(i,2)*768/(41*3.784);end fori=1:6 x(i)=m(i+6,1); y(i)=m(i+6,2);enda(1)=0;b(1)=0;c(1)=0;d(1)=0;e(1)=0;f(1)=0;g(1)=0;h(1)=0;j(1)=0;k(1)=0;z(1)=0;u(1)=0;v(1)=0;fori=1:6 a(i+1)=a(i)+x(i)^2*y(i)^2; b(i+1)=b(i)+x(i)^3*y(i); c(i+1)=c(i)+x(i)*y(i)^3; d(i+1)=d(i)+x(i)^2*y(i); e(i+1)=e(i)+x(i)*y(i)^2; f(i+1)=f(i)+x(i)*y(i); g(i+1)=g(i)+y(i)^4; h(i+1)=h(i)+y(i)^3; j(i+1)=j(i)+x(i)^2; k(i+1)=k(i)+x(i)^3; z(i+1)=z(i)+y(i)^2; u(i+1)=u(i)+x(i); v(i+1)=v(i)+y(i);endp(1)=a(7);p(2)=b(7);p(3)=c(7);p(4)=d(7);p(5)=e(7);p(6)=f(7);p(7)=g(7);p(8)=h(7);p(9)=j(7);p(10)=k(7);p(11)=1;p(12)=z(7);p(13)=u(7);p(14)=v(7);q(1,1)=p(1);q(1,2)=p(3);q(1,3)=p(4);q(1,4)=p(5);q(1,5)=p(6);w(1)=-p(2);q(2,1)=p(3);q(2,2)=p(7);q(2,3)=p(5);q(2,4)=p(8);q(2,5)=p(12);w(2)=-p(1);q(3,1)=p(4);q(3,2)=p(5);q(3,3)=p(9);q(3,4)=p(6);q(3,5)=p(13);w(3)=-p(10);q(4,1)=p(5);q(4,2)=p(8);q(4,3)=p(6);q(4,4)=p(12);q(4,5)=p(14);w(4)=-p(4);q(5,1)=p(6);q(5,2)=p(12);q(5,3)=p(13);q(5,4)=p(14);q(5,5)=1;w(5)=-p(9);x=DelGauss(q,w);s(2)=(x(1)*x(4)-2*x(2)*x(3))/(x(2)*4-x(1)^2);t(2)=(x(1)*x(3)-2*x(4))/(4*x(2)-x(1)^2);m=[18.00008.000015.00009.000015.000011.000018.000012.000019.00009.000019.000011.000021.00009.000024.00009.000021.000012.000022.000012.500024.500010.000023.500012.000032.000011.000032.000012.000035.500010.000035.00009.500033.00009.500035.500012.000030.000025.000029.000027.000032.000025.000032.500027.000032.000027.500030.500025.000013.000027.000014.000025.000015.000028.500016.000028.000017.000026.000016.500025.0000];fori=1:30 m(i,1)=m(i,1)*1024/(54*3.784); m(i,2)=m(i,2)*768/(41*3.784);end fori=1:6 x(i)=m(i+24,1); y(i)=m(i+24,2);enda(1)=0;b(1)=0;c(1)=0;d(1)=0;e(1)=0;f(1)=0;g(1)=0;h(1)=0;j(1)=0;k(1)=0;z(1)=0;u(1)=0;v(1)=0;fori=1:6 a(i+1)=a(i)+x(i)^2*y(i)^2; b(i+1)=b(i)+x(i)^3*y(i); c(i+1)=c(i)+x(i)*y(i)^3; d(i+1)=d(i)+x(i)^2*y(i); e(i+1)=e(i)+x(i)*y(i)^2; f(i+1)=f(i)+x(i)*y(i); g(i+1)=g(i)+y(i)^4; h(i+1)=h(i)+y(i)^3; j(i+1)=j(i)+x(i)^2; k(i+1)=k(i)+x(i)^3; z(i+1)=z(i)+y(i)^2; u(i+1)=u(i)+x(i); v(i+1)=v(i)+y(i);endp(1)=a(7);p(2)=b(7);p(3)=c(7);p(4)=d(7);p(5)=e(7);p(6)=f(7);p(7)=g(7);p(8)=h(7);p(9)=j(7);p(10)=k(7);p(11)=1;p(12)=z(7);p(13)=u(7);p(14)=v(7);q(1,1)=p(1);q(1,2)=p(3);q(1,3)=p(4);q(1,4)=p(5);q(1,5)=p(6);w(1)=-p(2);q(2,1)=p(3);q(2,2)=p(7);q(2,3)=p(5);q(2,4)=p(8);q(2,5)=p(12);w(2)=-p(1);q(3,1)=p(4);q(3,2)=p(5);q(3,3)=p(9);q(3,4)=p(6);q(3,5)=p(13);w(3)=-p(10);q(4,1)=p(5);q(4,2)=p(8);q(4,3)=p(6);q(4,4)=p(12);q(4,5)=p(14);w(4)=-p(4);q(5,1)=p(6);q(5,2)=p(12);q(5,3)=p(13);q(5,4)=p(14);q(5,5)=1;w(5)=-p(9);x=DelGauss(q,w);s(5)=(x(1)*x(4)-2*x(2)*x(3))/(x(2)*4-x(1)^2);t(5)=(x(1)*x(3)-2*x(4))/(4*x(2)-x(1)^2);m=[18.00008.000015.00009.000015.000011.000018.000012.000019.00009.000019.000011.000021.00009.000024.00009.000021.000012.000022.000012.500024.500010.000023.500012.000032.000011.000032.000012.000035.500010.000035.00009.500033.00009.500035.500012.000030.000025.000029.000027.000032.000025.000032.500027.000032.000027.500030.500025.000013.000027.000014.000025.000015.000028.500016.000028.000017.000026.000016.500025.0000];fori=1:30 m(i,1)=m(i,1)*1024/(54*3.784); m(i,2)=m(i,2)*768/(41*3.784);endfori=1:6 x(i)=m(i+12,1); y(i)=m(i+12,2);enda(1)=0;b(1)=0;c(1)=0;d(1)=0;e(1)=0;f(1)=0;g(1)=0;h(1)=0;j(1)=0;k(1)=0;z(1)=0;u(1)=0;v(1)=0;fori=1:6 a(i+1)=a(i)+x(i)^2*y(i)^2; b(i+1)=b(i)+x(i)^3*y(i); c(i+1)=c(i)+x(i)*y(i)^3; d(i+1)=d(i)+x(i)^2*y(i); e(i+1)=e(i)+x(i)*y(i)^2; f(i+1)=f(i)+x(i)*y(i); g(i+1)=g(i)+y(i)^4; h(i+1)=h(i)+y(i)^3; j(i+1)=j(i)+x(i)^2; k(i+1)=k(i)+x(i)^3; z(i+1)=z(i)+y(i)^2; u(i+1)=u(i)+x(i); v(i+1)=v(i)+y(i);endp(1)=a(7);p(2)=b(7);p(3)=c(7);p(4)=d(7);p(5)=e(7);p(6)=f(7);p(7)=g(7);p(8)=h(7);p(9)=j(7);p(10)=k(7);p(11)=1;p(12)=z(7);p(13)=u(7);p(14)=v(7);q(1,1)=p(1);q(1,2)=p(3);q(1,3)=p(4);q(1,4)=p(5);q(1,5)=p(6);w(1)=-p(2);q(2,1)=p(3);q(2,2)=p(7);q(2,3)=p(5);q(2,4)=p(8);q(2,5)=p(12);w(2)=-p(1);q(3,1)=p(4);q(3,2)=p(5);q(3,3)=p(9);q(3,4)=p(6);q(3,5)=p(13);w(3)=-p(10);q(4,1)=p(5);q(4,2)=p(8);q(4,3)=p(6);q(4,4)=p(12);q(4,5)=p(14);w(4)=-p(4);q(5,1)=p(6);q(5,2)=p(12);q(5,3)=p(13);q(5,4)=p(14);q(5,5)=1;w(5)=-p(9);x=DelGauss(q,w);s(3)=(x(1)*x(4)-2*x(2)*x(3))/(x(2)*4-x(1)^2);t(3)=(x(1)*x(3)-2*x(4))/(4*x(2)-x(1)^2);m=[18.00008.000015.00009.000015.000011.000018.000012.000019.00009.000019.000011.000021.00009.000024.00009.000021.000012.000022.000012.500024.500010.000023.500012.000032.000011.000032.000012.000035.500010.000035.00009.500033.00009.500035.500012.000030.000025.000029.000027.000032.000025.000032.500027.000032.000027.500030.500025.000013.000027.000014.000025.000015.000028.500016.000028.000017.000026.000016.500025.0000];fori=1:30 m(i,1)=m(i,1)*1024/(54*3.784); m(i,2)=m(i,2)*768/(41*3.784);endfori=1:6 x(i)=m(i+18,1); y(i)=m(i+18,2);enda(1)=0;b(1)=0;c(1)=0;d(1)=0;e(1)=0;f(1)=0;g(1)=0;h(1)=0;j(1)=0;k(1)=0;z(1)=0;u(1)=0;v(1)=0;fori=1:6 a(i+1)=a(i)+x(i)^2*y(i)^2; b(i+1)=b(i)+x(i)^3*y(i); c(i+1)=c(i)+x(i)*y(i)^3; d(i+1)=d(i)+x(i)^2*y(i); e(i+1)=e(i)+x(i)*y(i)^2; f(i+1)=f(i)+x(i)*y(i); g(i+1)=g(i)+y(i)^4; h(i+1)=h(i)+y(i)^3; j(i+1)=j(i)+x(i)^2; k(i+1)=k(i)+x(i)^3; z(i+1)=z(i)+y(i)^2; u(i+1)=u(i)+x(i); v(i+1)=v(i)+y(i);endp(1)=a(7);p(2)=b(7);p(3)=c(7);p(4)=d(7);p(5)=e(7);p(6)=f(7);p(7)=g(7);p(8)=h(7);p(9)=j(7);p(10)=k(7);p(11)=1;p(12)=z(7);p(13)=u(7);p(14)=v(7);q(1,1)=p(1);q(1,2)=p(3);q(1,3)=p(4);q(1,4)=p(5);q(1,5)=p(6);w(1)=-p(2);q(2,1)=p(3);q(2,2)=p(7);q(2,3)=p(5);q(2,4)=p(8);q(2,5)=p(12);w(2)=-p(1);q(3,1)=p(4);q(3,2)=p(5);q(3,3)=p(9);q(3,4)=p(6);q(3,5)=p(13);w(3)=-p(10);q(4,1)=p(5);q(4,2)=p(8);q