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

高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A
我们的参赛报名号为（如果赛区设置报名号的话）： 20002028
所属学校（请填写完整的全名）： 中南大学
参赛队员 (打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：

日期： 2008 年 9 月 22 日

赛区评阅编号（由赛区组委会评阅前进行编号）：

数码相机定位

摘 要

对于问题一，我们建立了两个模型，求当坐标系原点取在该相机的光学中心时，靶标上圆的圆心在该相机像平面的像坐标。一个是椭圆最小二乘法拟合模型（模型一），这个模型中我们利用了网格纸选点算法和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 靶标的像
要求：
（1） 建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；
（2） 对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为；
（3） 设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；
（4） 建立用此靶标给出两部固定相机相对位置的数学模型和方法。

二 问题分析

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基本假设
1) 考虑最小二乘法的要求，我们假设从图像坐标系中选取的点的坐标值的随机误差服从正态分布。
2) 图3中的椭圆都是标准椭圆，不存在畸变等因素，属于理想图形。
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）
由此建立如下数学规划模型：

我们用matlab 6.0编程并用消元法求解上述线性方程组，最终求得方程系数的值。

模型二 直线逼近模型
每个椭圆用两条水平直线和两条竖直直线逼近椭圆，使得四条直线与椭圆相切，得出每个椭圆水平、竖直边界点。Imread函数对处理的函数有特定的要求，在用imread函数调用图片之前，先用photoshop软件将图3图片格式转换为tiff格式并无任何压缩限制，然后用imread函数取得图片中椭圆区域。再在图中划分出每个椭圆的边界，如图6。

图6 椭圆边界划分

对每个区域用imread函数取得的图片数据中，分区域0HFJ,HIGF,ILKG,DKQN,JEMP在每个区域用水平线和竖直线去与椭圆相切，可以求出相切交点的坐标（以椭圆D为例）如图7。

图7 椭圆D相切交点位置

其中
由椭圆对称性可知，直线与直线交于一点且过椭圆中心，则椭圆中心可以用坐标表示，得到如下数学模型（在图像坐标系下）：
（10）

3.1.5模型的求解

模型一的求解
模型中首先要解决的问题是得到六个点，用这六个点来算出等式（7）中参数六个点的值。对于选取这六个点，给出一下两个算法；
1) 网格纸选点算法
用一个网格纸放在图形上面，建立一个坐标系，确定每一个网格对应像素，再取网格顶点正好落在每个椭圆边界处，当取不到足够的点时，可以移动网格纸位置，使得足够多的点落在网格顶点，如图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)
A
1
90.2044
39.6019
2
75.1703
44.5522
3
75.1703
54.4526
4
90.2044
59.4029
5
95.2157
44.5522
6
95.2157
54.4526
B
1
105.2384
44.5522
2
120.2725
44.5522
3
105.2384
59.4029
4
110.2498
61.8780
5
122.7782
49.5024
6
117.7668
59.4029
C
1
160.3633
54.4526
2
160.3633
59.4029
3
177.9031
49.5024
4
175.3974
47.0273
5
165.3747
47.0273
6
177.9031
59.4029
D
1
150.3406
123.7560
2
145.3293
133.6565
3
160.3633
123.7560
4
162.8690
133.6565
5
160.3633
136.1316
6
152.8463
123.7560
E
1
65.1476
133.6565
2
70.1590
123.7560
3
75.1703
141.0818
4
80.1817
138.6067
5
85.1930
128.7062
6
82.6873
123.7560

在matlab6.0中拟合出每个椭圆方程的参数值，见表2

表2 每个椭圆代数形式参数的拟合值
椭圆拟合
系数
坐标中
的点

A
-1.3953
-0.3494
-99.0371
158.7387
-0.7314
B
-1.3677
-1.8741
-159.2396
353.1355
-0.5491
C
-2.1063
-3.0887
-221.2915
686.4699
-0.4062
D
0.2270
1.7202
-337.1185
405.6205
-0.0189
E
-0.0150
0.2950
-148.0241
81.2751
-0.0887

由表2拟合的每个椭圆方程的参数值代入（2）（3），得到各椭圆圆心在图像坐标系下的坐标如下：

表3 图像坐标系下拟合各个椭圆中心坐标

椭圆
(mm)
A
B
C
D
E
X
86.9180
113.2818
167.5217
154.0275
75.0327
Y
56.6095
52.1480
53.0060
132.0610
137.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.0248
32.2151
18.7209
-60.2739
y
50.8704
49.3319
46.4739
-30.5811
-32.3639

3.2.2利用模型一中计算机函数取点算法求解
用函数在每个椭圆上取出六个点坐标值如表5
表5 各点取值

x（像素）
y（像素）
椭圆A
1
283
202
2
364
182
3
310
229
4
342
153
5
326
230
6
284
206
椭圆B
1
384
192
2
384
209
3
461
206
4
454
220
5
430
235
6
462
202
椭圆C
1
606
230
2
643
176
3
652
176
4
606
232
5
623
248
6
675
202
椭圆D
1
571
475
2
606
475
3
612
515
4
582
471
5
562
481
6
585
534
椭圆E
1
604
220
2
606
204
3
607
233
4
605
229
5
604
214
6
669
232

在matlab6.0中拟合出每个椭圆方程的参数值，见表6

表6 每个椭圆代数形式参数的拟合值
椭圆拟合
系数
坐标中
的点

A
-2.0407
0.8922
-243.6945
340.9528
-1.0293
B
-1.8508
-0.2235
-470.2341
920.9710
-2.2089
C
-0.0033
0.0010
-0.5666
1.6739
-0.0003
D
-1.0074
-0.2193
-684.4368
815.3615
-0.1040
E
-0.0030
0.0004
-0.5778
1.6644
-0.0011618

以一个单位像素作为刻度，图3中的图像可以划分为个单元格，用matlab库函数中的imread函数来取得每个单元格的值，再用问题一中的算法可以取得每个椭圆心坐标值，最后转换到以光学中心为原点的坐标值，见表7

表7 图像坐标系下各椭圆圆心坐标值

A(mm)
B(mm)
C(mm)
D(mm)
E(mm)
x
-50.0000
-23.5449
33.8625
18.7831
-60.1851
y
51.4550
49.6032
45.1058
-31.4815
-31.2169

3.2.3利用模型二求解
将图片转换为tiff图片格式，用取色器在转换后图片上找出椭圆边界点的RGB值，可以得到表8值

表8 tiff格式图片上RGB值

R
G
B
椭圆
边界
内和
椭圆
上可
取值
1
51
51
51
2
41
41
41
3
34
34
34
4
28
28
28
5
22
22
22
6
12
12
12
7
4
4
4
椭圆
边界
和边
界外
可取
值
1
234
234
234
2
227
227
227
3
221
221
221
4
215
215
215
5
204
204
204
在表8中，椭圆内部的RGB值都为，椭圆外部的RGB值都为，其他的值都分布在椭圆边界附近。
以一个单位像素作为刻度，图3中的图像可以划分为个单元格，用matlab库函数中的imread函数来取得每个单元格的值，再用问题一中的算法可以取得每个椭圆心坐标值，最后转换到以光学中心为原点的坐标值（代码见附录5），见表9
表9 光心坐标系下各椭圆圆心坐标值

A
B
C
D
E
x
-50.0000
-23.5449
33.8625
18.7831
-60.1851
y
51.4550
49.6032
45.1058
-31.4815
-31.2169

3.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_index
best_elipse_par=椭圆参数
max_inliers_index=inliers_index
max_inliers=num_inliers
结束
num_inliers>max_inliers
是
否
否
是

图9 算法流程图

3.3.3算法对问题一模型的检验
现取题二中的椭圆A作为分析对象(取点代码见附录7)，样本总体上有233个样本点，则有227个样本点来计算点到拟合椭圆的距离，遍历100次。为得到精确的结果，阈值应越小越好，又因为问题二中的图片上椭圆比较多标准，阈值在这里取为2个像素单位，在matlab中运行代码，对数值进行遍历，最终得到p中保存的最优椭圆参数，见表10

表10 p中保存的最优椭圆参数

B
C
D
E
F
数值
-2.5400
0.9736
-204.7289
418.6218
-6.2642

将得到的数值

A
B
C
D
E
x
-50.0000
-23.5449
33.8625
18.7831
-60.1851
y
51.4550
49.6032
45.1058
-31.4815
-31.2169

通过与问题二的数值比较发现，问题一的假设是合理的，且用问题一种建立的模型，对比较标准的图形可以得到很精确的

3.3.4算法评价
在直接最小二乘法的情形中，并没有考虑到畸变给结果带来的误差。当给定图片存在肉眼看不见的畸变或杂质点，如果用直接最小二乘法可能会产生比较小的，或者是不易察觉的误差，分析但这些误差可能会在生产工作中慢慢的积累，最终可能产生破坏性的影响。例如在眼睛治疗中，需要对眼膜有非常高的精度要求，随机化方法可以很好的消除因设备、杂质点等带来的误差
随机化算法是通过一定量的遍历来消除畸变、杂质点等影响，得到想要的到的结果，有时这种算法能达到亚像素精度。

3.3.5算法精度检验
1)精度概念分析
当给定一个误差范围，通过拟合值与真实值之间的比较，得到误差平法和、方差等统计量来对结果进行检验

2) 检验方法描述
a) 给定一个图形，并知道图形所有点坐标值、参数值等；
b) 取到每个点的坐标值，并储存；
c) 对每一个图形上的点用matlab中normrnd函数生成正态（高斯）分布的随机数，将这个随机数加到图形对应点中，从而得到一个和原图形有偏差的新图形；
d) 以新图形上的点坐标作为已知值，用随机化算法找到最接近原图形的参数值
e) 用新得到的参数值与已知参数值做统计量比较，检验随机化算法拟合精度。

3) 检验实例
假设a,b值拟合误差超过1%，认为算法不合格，否则，证明此算法在拟合方面有很高的精度
a) 给定一标准椭圆，中点在原点，a=50,b=30；
b) 取到每个点的值，以1为一个单位，总共有200个点x(i),y(i)(i=1,2,…,200)；
c) 用normrnd（0，0.1）生成400个白噪声数值nx(i,1),ny(i,2) (i=1,2,…,200)，组成一个白噪声序列，每两个用来控制一个点在x轴和y轴上移动的值;将对应的点在白噪声nx,ny的影响下产生新的位置值x1(i),y1(i);
d) 将数据x1(i),y1(i)用随机化算法，来求得a1,b1的值；
e) 对拟合值a1,b1进行误差分析

其中将(x1(i),y1(i))点和原椭圆画在一个图中，如图10
图10 随机产生点与原椭圆

用随机化算法计算得a1,b1值分别为50.1205，29.8781
重新生成白噪声序列，再重复进行多次的拟合，最终得到表11中数据（代码见附录8）
表11 拟合值

1
2
3
4
5
6
7
8
9
10
a1
50.1205
49.9998
50.0000
50.2001
49.9731
49.7391
50.0317
49.6902
50.1000
49.8973
b1
29.8781
29.9173
30.0104
30.2178
30.0789
30.1083
30.1105
30.1157
30.0306
30.0338

分别计算出a1,b1样本均值有
a1,b1的样本方差
计算每组数据的相对误差如表12
表12 相对误差

1
2
3
4
5
6
7
8
9
10
a1
0.0024
0.0000
0.00000
0.0040
0.0005
0.0052
0.0006
0.0062
0.0020
0.0021
b1
0.0041
0.0028
0.0003
0.0073
0.0026
0.0036
0.0037
0.0039
0.0010
0.0011
从先对误差表中可得a1最大误差为0.4%，b1最大误差为0.73%，都没超过1%，证明此算法有很高的精度。由于产生随机序列的值可能取到很大，则会出现较大奇异点，使得在其他模型的求解过程中产生很大的误差，而随机化算法可以把这些有较大奇异点的地方，通过筛选把奇异点排除在正确解的范围内。能够有效地处理包含有较大比例误差点的样本空间，拟合出具有高精度的椭圆，并且算法的速度能够满足实时性的要求。

3.3.6算法稳定性检验
1) 稳定性概念分析
在一组影响变量a的因素中，当某一个因素发生一个小的变化，会使得变量发生很大的变化，则变量不稳定，如果在一定的范围内，影响因素的变化不能影响变量发生很大的变化，则这个变量是稳定的。
在随机化算法中，主要解决的问题是通过循环排除畸变或奇异点等等影响，如果模型中出现较大的畸变或出现大的奇异点，而随机化算法依然能得到拟合效果很好的拟合情况。否则认为随机化算法稳定性不够。

2) 检验方法描述
1) 给定一个图形，并指导图形在坐标系下的所有情况和图形本身参数值；
2) 给图形制造一定范围的畸变或者奇异点；
3) 用新图形的坐标值拟合求解图形的参数值；
4) 进行参数检验，确定随机算法是否具有较高的稳定性

3) 检验实例
假设c,d值拟合误差超过5%，认为算法不合格，否则，证明此算法在拟合方面有很高的精度
a) 给定一标准椭圆，中点在原点，c=50,d=30；
b) 取到每个点的值，以1为一个单位，总共有200个点x(i),y(i)(i=1,2,…,200)；
c) 从y=40出小区图形右边部分，用normrnd（0，0.1