《数码相机的定位2008年全国大学生数学建模A题一等奖(16页).doc》由会员分享，可在线阅读，更多相关《数码相机的定位2008年全国大学生数学建模A题一等奖(16页).doc（16页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-数码相机的定位2008年全国大学生数学建模A题一等奖-第 16 页数码相机的定位摘要本文运用相机标定模型确定了相机像平面的像坐标，利用本质矩阵标定双目相机，快速找出了相机的相对位置关系；利用MATAB软件和图像处理进行编程求解；通过对图像的预处理和灰度质心法对模型进行了验证，得出模型的精度。针对问题一，根据数码相机的特点，提出了一个新的标定方法，建立相机标定模型，确定了靶标上圆的圆心在该相机像平面的像坐标，为问题二的计算提供了一个好的算法。针对问题二，我们利用问题一建立的模型和方法运用MATLAB编程精确的计算了靶标上五个圆的圆心在像平面上的像坐标。ABCDEx86.2453110.8411
2、168.2104155.171273.2152y50.009950.797652.6354134.09100.2646针对问题三，我们引入了灰度质心法及像差模型对前述问题的模型的稳定性和坐标值精度进行检验后，发现两种模型的中心坐标值的误差值在03个像素区间内，说明前述模型的计算结果的精度很高，通过像差模型得出其径向畸变系数趋于无穷小，认为前述模型有很好的稳定性。针对问题四，我们提出了一种改进的的立体摄像机标定方法，通过双目匹配点，线性地求解本质矩阵，快速找出摄像机的相对位置关系。关键词：双目定位 系统定标 灰度质心法 像差模型 本质矩阵1. 问题重述数码相机定位在交通监管（电子警察）等方面有广
3、泛的应用。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像
4、平面上我们都无法直接得到没有几何尺寸的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，如图1所示，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现。 图 1 靶标上圆的像有人设计靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以AC边上距离A点30mm处的B为圆心，12mm为半径作圆，用一位置固定的数码相机摄得其像。请你们：（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；（2） 对由
5、图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024768；（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。2. 模型假设1、假设拍摄靶标的像的数码照相机的光学系统在制造和安装过程中产生的各种畸变误差很小，可忽略不计；2、假设拍摄靶标的相机的焦距是固定的，且相机的焦距等于像距；3、设像面中心相对光轴与像面的焦点在 方向和 方向的偏移为0；4、假设像面的倾斜程度为0。3. 符号说明问
6、题一：世界坐标系坐标的表示：图像坐标系的坐标表示：像坐标系的坐标表示：光心坐标系坐标表示：焦距：33阶的旋转矩阵:位移向量：图像坐标系原点在像素坐标系中的横坐标：图像坐标系原点在像素坐标系中的纵坐标：像素坐标系在方向相邻像素间的距离:像素坐标系在方向相邻像素间的距离:数码相机的内部参数矩阵：数码相机的外部参数矩阵：数码相机的投影矩阵问题二：靶标上A圆圆心在成像平面上对应的像素坐标：靶标上B圆圆心在成像平面上对应的像素坐标：靶标上C圆圆心在成像平面上对应的像素坐标：靶标上D圆圆心在成像平面上对应的像素坐标：靶标上D圆圆心在成像平面上对应的像素坐标：靶标上A圆圆心在成像平面上对应的像坐标：靶标上B
7、圆圆心在成像平面上对应的像坐标：靶标上C圆圆心在成像平面上对应的像坐标：靶标上D圆圆心在成像平面上对应的像坐标：靶标上E圆圆心在成像平面上对应的像素坐标问题三：i ： 表示坐标点底数， ： 表示第i个质点的位置坐标n ： 表示图像的像素个数 ： 表示第i个像素的坐标 ： 表示第i个像素的灰度值 ： 表示第i个质点的质量x0，y0 ： 表示全部n个离散质点的重心坐标（注：问题四的符号说明已在问题中给出）4. 问题分析计算机视觉是一门新兴的学科。随着计算机硬件、软件、图像采集、处理技术的迅速发展，计算机视觉的理论和技术已被广泛地应用于医学图像处理、机器人技术、文字识别、工业检测、军事勘察、地面勘察
8、和现场检测等。计算机视觉在多种测量中的应用是一种定量分析系统，有确定的精度要求。一般运用于现场不可到达、现场复杂不便直接测量、保留现场状态以便测量、补测和事后重构现场等场合。这主要是由于数码相机有如下几个特点：有较多像素的CCD；可以根据测量的需要选择不同的像素大小；可直接与计算机进行通讯；可选定多种焦距进行顶焦距拍摄。上述特点，给图像采集和图像处理带来了许多方便之处。在应用计算机视觉的测量系统中，无论是采用摄像机还是数码相机作为原始图像的采集设备，都必须进行系统标定。系统标定是应用计算机视觉进行现场测量和图像处理的第一步，也是关键的一步。标定矩阵的精度直接影响到最终的测量精度。我们运用数码相
9、机成像的几何模型，并将次模型分解成内、外参数矩阵。根据数码相机的特点，提出了一种新的标定方法。这种方法是先求解有较高精度的内部参数矩阵，在此基础上再来求解外部参数矩阵，由此获得高精度的内、外参数矩阵。5. 模型的建立与求解5.1 问题一 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标，这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面。5.1.1问题的分析本题要求我们通过建立合适的模型以确定靶标上的圆的圆心在相机平面的像的像坐标。数码相机是将空间三维图像转换为平面二维图像的传感设备，我们先建立空间三维物体坐标点与二维图像坐标点的对应关系，建立数码相机的几何模型。通过此模型
10、可以求出空间特征点的图像坐标值。相机的定标主要是确定数码相机的内外参数。外部参数用来把目标的世界坐标系转换为数码相机坐标系。内部参数用来确定数码相机的内部几何与光学系统。5.1.2模型的建立和求解数码相机图像拍摄实际上是一个光学成像过程。可以将此过程分为三个步骤、隶属于四个坐标系。这四个坐标系分别为：（1）世界坐标系根据自然环境所选定的坐标系，坐标值用表示。（2）图像坐标系坐标原点在CCD图像平面的中心，X轴，Y轴分别为平行于图像平面的俩条垂直边，坐标值用（x,y）表示。（3）像素坐标系坐标系在CCD图像平面的左上角，X轴，Y轴分别平行于图像坐标系的X轴和Y轴，坐标值用（u,v）来表示，且为离
11、散的整数值。（4）光心坐标系以相机的光心为坐标系点，X轴，Y轴分别为平行于图像坐标系的X轴和Y轴，相机的光轴为Z轴，坐标值用表示。3个步骤是分别将世界坐标系中的信息转换到光心光标系，再有光心坐标系转换到图像坐标系，最后有图像坐标系转换到像素坐标系。光学光学成像的理论模型为针孔模型。根据这个模型有光心坐标系向图像坐标系的过程符合中心影射或透视投影，可用齐次坐标与矩阵表示其中，是光心坐标系中空间点P的坐标，是对应图像坐标系中P点的坐标，f是相机的焦距。可用齐次坐标与矩阵到光心坐标系的转换关系为：其中，R为旋转矩阵，t为位移向量， 元素为0的列向量。由图像坐标到像素坐标系的转换关系为其中，是图像坐标
12、系圆点在像素坐标系中的坐标，（dx,dy）分别是像素坐标系在X方向和Y方向相邻像素间的距离。将式代入的可将式简化成N为式4右边地1项即相机的内部参数矩阵，M为式4右边第2项即相机外部参数，Q为投影矩阵。5.2问题二的分析5.2.1 问题分析：本题要求我们利用问题一建立的模型或方法来精确计算靶标上五个圆的圆心在靶标图上的像的坐标。由相机的像距及相机的分辨率等信息估计出相机的内部参数，然后基于圆的特征性质以及透视投影的光学原理在靶标和其像上采集六对坐标点。代入模型计算出投影矩阵后便可计算靶标上圆的圆心的像点坐标。5.2.2问题求解具体计算过程如下：首先建立相应坐标系，使靶标及其像图坐标化。如图5.
13、1所示。对于靶标建立如下的三维坐标体系。靶标上所有点的世界坐标为0.12mm /mm 100mm13mm /mm 0 图5.1 世界三维坐标系对于靶标在像平面的像图，建立如图二维坐标系为了高精度提取像图上与靶标上点对应的点的坐标值，将图像网格化处理。 图5.2 图像像素坐标图我们在A、B、C、D五个圆像上采集六对点的坐标值（其中在A 、B、C、D上采一个点，E上采集俩个点 ），六对点的坐标在附录中给出。 然后将六对点代入模型求得参数如下：f=417.1958; dx=0.2646; dy=0.2646;已知靶标上A、B、C、D、E圆的圆心坐标为：A（13，113）、B（43，113）、C（11
14、3，113）、D（113，13）、E（13，13）靶标上圆圆心在成像平面上对应的像素坐标为：靶标上圆圆心在成像平面上对应的像坐标为：5.3模型的检验5.3.1模型检验的方法说明由问题一中建立的数码相机像坐标标定的模型，我们通过数码相机的内部和外部参数计算出了问题二中各圆的像坐标值。但我们并未知道问题二中所计算得出的像坐标值的精确度到底有多高，以及数码相机的一阶畸变系数的大小。因此我们对问题一中提出的模型得到的像坐标值进行检验，通过对各种传统标定方法、自标定方法和基于主动视觉的标定方法的学习，我们采用了一种灰度质心法和像差模型相结合检验的数学检验方法。 要使用灰度质心模型对前述问题的模型进行检验
15、，首先，我们在matlab下对题目所给图骗进行了灰度处理（灰度处理程序见附件三）；其次，通过其他方式确定图像中圆或者椭圆的内中心点附近的像素点坐标值，最后即可运用灰度模型进行求解5.3.2灰度质心法在高等数学中，我们都知道平面上n个离散的质点，它们的质量重心坐标可以表示为： (1) (2)考虑到圆或者椭圆的中心对称性，其密度均匀的圆或者椭圆上的中心是重合的，其中将（1）、（2）式中的点质量用点的灰度值替换，质量重心变成了灰度重心，由此圆形或者椭圆形光斑的中心也可以通过求取光斑在像平面上上网灰度重心来得到。 原始图像结果在matlab中进行灰度处理后，把其明暗对比度加深，得到被测特征点成像的特征
16、区域成圆形或者椭圆形的光斑，设各圆形或者椭圆形光斑所占像素数量为，于是灰度质心法计算各圆形或者椭圆形光斑中心的坐标可表示为； (3) .(4)根据上面的模型，第一中方式我们采用的是：通过一种计算机软件对图像的像素坐标进行多点的坐标取值，根据这些坐标值运用matlab软件编程(程序见附件三)求出其图像中心点在像素坐标系的值，然后运用灰度质心法中 (3)、(4)关系式，把各像素坐标值带入相应的程序中，求出畸变下的像素坐标值，求的各图像的像素坐标值如下表所示： 表5.1各圆和椭圆图像的像素坐标ABCDEX325.1000418.8000641.3529584.0714282.7000Y189.263
17、2195.1500207502.8571502.8500（注：表中个坐标值的单位是像素）其对应的像坐标值有题目所购的换算关系的出： 表5.2 各图像的像坐标值ABCDEXc86.0053110.7937169.6701110.793774.7884Yc50.069651.627054.7619133.0291133.0310（注：表中个坐标值的单位是毫米）从表中，我们经过与问题二中求解的值进行对比可以得出：问题一中所建立的模型的求解的值精度较高。 但在摄像机模型中全面考虑的话，我们还需要引入镜头产生的径向畸变和且想畸变在摄像机摄像过程中对成像效果的影响，因此我们建立了一个关于畸变参数求解的模型
18、像差数学模型。其数学表达式如下： .(5) .(6)式中： -畸变坐标 -径向畸变系数 -切向畸变系数由上述数学关系式，我们通过在matlab下编程(程序见附件三)的到了各畸变系数的值，下表所示： 表5.3径向、切向畸变系数k1k2k3k4-Inf-Inf-2.5289e-006-1.4292e-006有上述各表可以知道，通过该模型对像素坐标值和像坐标值以及畸变系数的计算的到的结果值，和问题一建立模型计算得出的数值进行比较后，得出其精度很高，像素坐标的各坐标值的误差值B在像素点这样一个小区间里。所以说明摄像机的摄像精度较高。而畸变系数的值较小，有此认为摄像机在不同角度对物体进行摄像时，其图像的
19、畸变较小，图像不会产生严重失真。通过模型的检验，我们认为问题一的模型适合运用在定点摄像的场所，比如，公路电子监控、室内摄像的场所。5.4一种改进的立体数码相机标定方法对于物体上的一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是根据这一方法我们建立了双目模型并利用本质矩阵标定双目摄像机，解决了两部固定相机的相对位置。5.4.1双目模型的建立(5.12)在双目摄像机系统中，每个摄像机都分别符合图5.4.1所描述的成像模型 。而通过单个摄像机
20、模型，可以分别知道左、右摄像机的图像点与3维世界坐标的点的对应关系，即及的对应关系已分别知道。而双目标定所需要解决的问题是两个摄像机坐标系之间的对应关系 ,即Pc1=(xc1,yc1,zc1)与Pc2= ( xc2，yc2，z c2)之间的转换关系 。而由于它们分别属于3维直角坐标系，因此它们之间的转换包含有旋转以及平移的关系，形如 其中，旋转矩阵R以及平移向量T就是双目标定所要求的内容。 图5.3 双目摄像机模型5.4.2 本质矩阵的基本原理1）、符号说明1. 对左摄像机成像平面的2维坐标，对应于摄像机坐标系的点，将其改写为其次向量的形式；同理，对右摄像机，；2. 把3维的世界坐标亦记作齐次
21、向量的形式；3. 世界坐标到成像平面的 2维坐标的变换都是线性的 ，因此可以把左右摄像机的针孔成像过程写成以下的矩阵形式：，其中，A1,A2都是34的矩阵4. 为下面讨论方便起见，把Ai看作是一个方阵与一个向量的组合，即 和，其中Mi 是3的方阵，Ti 是31的向量；5. 向量，定义运算定义表示两个矩阵只相差一个比例因子，即为实数。不难看出，这样定义的符号“”满足交换性和传递性。2）、本质矩阵及其性质使用符号说明中的记号方法，把左右摄像机成像平面的坐标分别记为向量形式， 文献9指出，必然存在一个33的方阵Q，使得矩阵Q被称为本质矩阵。如何利用Q求解两台摄像机的相对位置呢？下面先考察本质矩阵的性
22、质。定理1：若左右摄像机的投影矩阵分别是=和，那么相应于这两台摄像机的本质矩阵为由这个定理知道，本质矩阵Q可以分解为一个满秩的方阵R与一个反对称矩阵S的乘积，即。以下的定理表明，这样的分解几乎是唯一的（最多相差一个尺度因子）。定理2：若一个33的方阵Q能够有两种的分解方法，其中Si是非零反对称矩阵，Ri是满秩的方阵，那么必然有，并且，如果，那么，是某个3维的向量。定理3：若左右摄像机的相对位置固定，A1,A2和A1,A2对应于相同的本质矩阵Q当且仅当存在着一个满秩的44的方阵H，使得P1HP1且P2HP2.以上三条定理的详细证明可以参考文献8和9.5.4.3利用本质矩阵标定双目摄像机上一节的定
23、理反映了一个事实：两台摄像机的相对位置关系与本质矩阵几乎是等价的。即给定了两台摄像机的相对位置关系R和T，就给定了本质矩阵Q，两台摄像机的旋转平移关系R和T也几乎确定了，最多像差一个比例常数因子，“几乎”的意思就是R和T的确定关系还需要标定一个比列因子K。 有参考文献【】给出了一个方便的通过Q求解R和T的方法，即先把Q进行特征值分解的，根据上街定理1得到， 其中, ,而t是把平移矩阵T的三个元素看成向量。 如此，我们确定了旋转矩阵和平移关系的基本形势R，T。而R，T与R，T相差一个比列因子k，即R=kR，T=kT，k的标定可以借助世界坐标中已知距离的两个点和。5.4.4 本质矩阵的求解 本小节
24、给出一种新的本质矩阵的求解方法。与已有的算法相比，新方法显得更为容易理解和实现。本质矩阵Q的定义本身就是求解的依据。假设是左右眼匹配点坐标构成的向量，记和，满足展开，得 (5.4.1)其中，q11，q12，q13,.,q33为矩阵Q的九个元素，也是要求接的未知数。一般地，可以带入九个匹配点的坐标组成方程组来求解。但根据线性代数理论，这是一个齐次方程组，当系数矩阵的行列式为零时，有无穷多解。若找到一个特解Q0。那么通解的形式是Q=。下面我们只需要找到这样的一个特解Q0令q33=1.则（5.4.1）变为 .(5.4.2)则方程组不再是齐次的，而存在唯一的解。下面是用最小二乘法来求得方程组的解。 在
25、左右眼找到N个匹配点和，i=1，k，N代入式（5.4.2），得到一个过约束方程组其中，A是的系数矩阵，q-（q11，q12q13，.,q33）是一个的解向量，过约束方程组的解 如此，就求得了Q的一个特解Q0，Q=。这里，是一个特定的常数，它的确定可以结合5.4.3的待定常数k一起进行。5.4.5具体标定步骤 1. 利用传统的方法分别对左右两台摄像机的内部参数进行标定。控制点通过精确打印的靶标图案提取，该图案沿着平行于光轴方向多次放置，精确记下各次放置的位置。 2.把已完成单目标定的两台摄像机安装到时间系统中，固定。 3.分别在左右影像中确立匹配像素度，列出方程，并按照5.4.3和5.4.4接触
26、本质矩阵Q,匹配点的位置没有严格要求，但力求均匀分布。同时，因为方程的求解用了最小二乘法，所以匹配点的数量越多，方程的解越精确。为了精确的在图像上提取匹配点，建议使用一种被广泛采用的类似于靶标的图案的摆放位置。控制点的分布根据两台摄像机的基线长度而定，尽可能的遍布摄像机的公共视角。4.根据5.4.3的原理求出两台摄像机的相对位置R和T。待定比列因子k。5.最后，标定比列因子k。具体做法是先设定两个表定点，测出他们的距离d0;其次，通过R和T求出其在摄像机坐标系下的距离d，因为R和T是k的线性函数，所以d也是k的线性函数，即。解一个一次的线性方程（f（k）=k0）即可得到比例因子k。6模型的评价
27、与改进参考文献：附件：1问题二的MATLAB算法程序：f=417.1958; %焦距dx=0.2646;%水平方向上的点距dy=0.2646;%垂直方向上的点距u0=512;v0=384;Zc=417.1958A=322,146,1;462,190,1;677,209,1;574,536,1;287,462,1;242,502,1%像素坐标值B=13,125,0,1;55,113,0,1;125,113,0,1;113,1,0,1;13,25,0,1;1,13,0,1%与像素坐标相对应的世界坐标值Q=(Zc*A)/B %系数的确定D1=Q*13,113,0,1;43,113,0,1;113,1
28、13,0,1;113,13,0,1;13,13,0,1/417.1958 %像素坐标的个像素坐标值D=Q*13,113,0,1;43,113,0,1;113,113,0,1;113,13,0,1;13,13,0,1/(417.1958*3.78) %像坐标中各坐标点的值D1=326.0413 418.9795 635.8352 586.5473 276.7535 189.0375 192.0148 198.9619 506.8640 496.9396 1.0000 1.0000 1.0000 1.0000 1.0000D =86.2543 110.8411 168.2104 155.1712
29、73.2152 50.0099 50.7976 52.6354 134.0910 131.46550.2646 0.2646 0.2646 0.2646 0.26462.灰度质心法的MATLAB程序%灰度处理附件三（1）a=imread(1.jpg);a=rgb2gray(a);imhist(a,256);imshow(a);b=imadjust(a,0/255,150/255,);figure,imshow(b);figure,imhist(b,256);%方式一附件三（2）%AX=305,319,334,348,286,303,320,336,353,283,303,319,336,351
30、,303,319,333,399,319,333;K=length(X);Pxy=1;Sx=sum(X.*Pxy);XA=Sx/(K*Pxy)Y=158,159,159,159,175,166,167,190,190,190,190,192,206,207,209,209,223,223,224;Sy=sum(Y.*Pxy);N=length(Y);Ya=Sy/(N*Pxy)XA =325.1000YA =189.2632%BX=416,429,398,415,429,449,384,401,415,432,449,385,399,415,430,447,401,416,444,422;Y=1
31、59,159,173,176,176,175,190,189,192,190,189,207,207,207,207,209,223,221,223,231;K=length(X);Pxy=1;Sx=sum(X.*Pxy);XB=Sx/(K*Pxy)Sy=sum(Y.*Pxy);N=length(Y);YB=Sy/(N*Pxy)XB=418.8000YB=195.1500%CX=638,653,624,639,655,672,608,622,639,655,670,607,622,638,653,622,639,653,642;Y=177,176,192,190,190,191,206,209
32、,209,206,206,223,224,223,221,240,240,240,248;K=length(X);Pxy=1;Sx=sum(X.*Pxy);XC=Sx/(K*Pxy)Sy=sum(Y.*Pxy);N=length(Y);YC=Sy/(N*Pxy)XC=641.3529YC=207%DX=573,605,607,560,576,591,608,560,574,591,608,557,576,591;Y=478,478,478,494,495,495,497,512,511,511,511,526,526,528;K=length(X);Pxy=1;Sx=sum(X.*Pxy);X
33、D=Sx/(K*Pxy)Sy=sum(Y.*Pxy);N=length(Y);YD=Sy/(N*Pxy)XD = 584.0714YD =502.8571%E%X=272,288,305,255,369,388,305,319,254,269,288,303,269,286,320;%Y=478,480,477,495,495,497,494,495,511,509,509,511,528,526,497;%K=length(X);%Pxy=1;%Sx=sum(X.*Pxy);%XE=Sx/(K*Pxy)%Sy=sum(Y.*Pxy);%N=length(Y);%YE=Sy/(N*Pxy)%X
34、E =282.7000%YE =502.8500X1=325.1000,418.8000,641.3529,418.8000,282.7000/3.78;Y1=189.2632,195.1500,207,502.8500,502.8571/3.78;X1 =86.0053 110.7937 169.6701 110.7937 74.7884%横的像坐标Y1 =50.0696 51.6270 54.7619 133.0291 133.0310%纵向像坐标%附件三（4）Xu=326.0413;Yu=189.0375;Xd=325.1000;Yd=189.2632;P=sqrt(Xd.*Xd)+(Y
35、d.*Yd)Wx=Xd(1)-Xu;Wy=Yd-Yu;k2=0;k4=0;%S=solve(325.1*(k1*p2+k2*p4)+k3*(3*(325.1)2+(189.2632)2)+2*k4*325.1*189.2632-Wx=0,189.2632*(k1*p2+k2*P4)+2*k3*325.1000*189.2632+k4*(325.1000)2+3*(189.2632)2)-Wy=0,418.8000*(k1*p2+k2*p4)+k3*(3*(418.8000)2+(195.1500)2)+2*k4*418.8000*195.1500-Wx=0,195.1500*(k1*p2+k2
36、*P4)+(2*k3*418.8000*195.1500)+k4*(418.8000)2+3*(195.1500)2)-Wy=0,k1,k2,k3,k4)k1=.50053831821545565380810169004052e-28*(-3261649625998908447875000.*p4*Wx-115551837114149400129533841.*p4*Wy+109168378189999454839404173.*P4*Wx+3003454391997036349100500.*P4*Wy)/p2/(-1.*p4+P4)k2=.5005383182154556538081016
37、9004052e-28*(105906728564000546391529173.*Wx-112548382722152363780433341.*Wy)/(-1.*p4+P4)k3=.26915641793591477296728378420238e-5*Wx+.20500622845529555660010816871490e-7*Wyk4=.14089603352087248947792101470719e-5*Wx-.45594569558163086224030393272402e-6*Wy%S=solve(Xd*(k1*p2+k2*p4)+k3*(3*Xd2+Yd2)+2*k4*X
38、d*Yd-Wx=0,Yd*(k1*p2+%k2*P4)+2*k3*Xd*Yd+k4*(Xd2+3*Yd2)-Wy=0,k1,k2,k3,k4)k1 = -Infk2 = -Infk3 =-2.5289e-006k4 = -1.4292e-006%方式二附件三（3）clccleardata=zeros(768,1024);data=imread(1.bmp); sumx=0; sumy=0; sum=0; for i=400:600 for j=400:1000 if data(i,j)1 sum=sum+1; sumx=sumx+j; sumy=sumy+i; end; end; end; x=sumx/sum y=sumy/sum