基于两步法的数码相机定位-毕业论文.doc

基于两步法的数码相机定位
摘要
数码相机定位在机器自动装配系统、工业视觉检测与识别、三维重建、机器人视觉导航、运动分析、海上目标跟踪、交通监管（电子警察）等诸多领域中得到了运用。
本文给出了确定靶标上圆的圆心在该相机像平面的像坐标的数学模型及确定两部固定相机相对位置的数学模型，并设计出了相应的求解算法。
首先在仅考虑单相机的情况下，在分析相机成像原理和四个坐标系之间的相互变换关系的基础上，考虑了相机径向畸变和切向畸变即非线性畸变因素，选择了一种简化模型，克服了相机内外参数未知情况下求解像坐标的困难，最终建立了基于两步法的像坐标确定模型。该模型满足牛顿迭代法的收敛条件，保证了模型解的稳定性。
其次利用该模型，针对问题1和问题2，借助于Matlab工具，计算了靶标上给定5个圆的圆心像坐标。
然后选择Canny算子对给定靶标的像的几何中心进行了精确检测，并对两种结果进行了对比，分析了误差，精度及稳定性。比较结果如表1所示：
表1 两种算法所得像坐标结果对比（单位：像素）
圆形标号
Canny算子中心定位法
两步法确定像坐标法
误差(%)
Ｕ坐标
Ｖ坐标
Ｕ坐标
Ｖ坐标
Ｕ坐标
Ｖ坐标
A
-189.0056
193.7568
-190.2873
194.7648
0.6781
0.5202
B
-90.0491
187.1209
-89.4045
185.9878
-0.7158
-0.6055
C
129.0463
170.7327
133.1423
176.1031
3.1741
3.1455
D
72.813
-120.7012
70.8932
-124.2142
-2.6366
2.9105
E
-227.2402
-117.3653
-224.9869
-115.9681
-1.0708
-1.1905
注：物坐标原点和像坐标系原点都在图片的几何中心。
Ｕ坐标的平均误差：1.6551，Ｖ坐标的平均误差：1.6754 平均误差：1.6653
最大误差点为C点，最小误差点为A点。
根据表1的数据证明了两步法确定像坐标的模型具有一定的可靠性和实用性。
然后根据给定靶标模型，并在上述模型的基础之上建立了确定两部相机相对位置的模型，同时给出基于平行线“消隐点”理论的切实可行的解法。
本文最后讨论了文中所建模型和所给算法的优缺点及改进方向。

关键词：两步法；像坐标；内外参数；边缘检测；相对位置；平行线消隐点

18

一、 问题的提出
1.1 背景说明
数码相机定位在交通监管（电子警察）等方面有广泛的应用。由于目前数字图像的处理速度越来越快，且可达0.02个像素精度[1]，因此考虑畸变系统误差的高精度标定具有重要的意义。所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。
本文是通过确定数码照相机的位置，属性参数并建立成像模型，从而确定空间坐标系中物点同它在图像平面上像点之间的对应关系，并通过对所得到的结果进行相关的处理，最终得到其在像平面理想的像坐标。
1.2 重述
有人设计靶标如下A，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以C边上距离A点30mm处的B为圆心，12mm为半径作圆，如下图1所示。

图1靶标示意图
由于图形的图像与拍摄点的位置有关，所以一下得到用一位置固定的数码相机摄得其像，如图2所示。

图2 靶标的像

问题1：建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标，这里坐标系原点取在该相机的光心，x-y平面平行于像平面；
问题2：对由图1、图2分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标，该相机的像距（即光心到像平面的距离）是1577个像素单位(1毫米约为3.78个像素单位)，相机分辨率为1024×768；
问题3：设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；
问题4：建立用此靶标给出两部固定相机相对位置的数学模型和方法。
二、 模型假设及符号说明
2.1 模型假设：
1）假设相机坐标平面与成像坐标平面平行；
2）假设相机坐标系原点在成像平面的投影是像坐标的原点；
3）假设噪声对相机成像无影响；
2.2 主要符号说明：
1）：全局坐标系；
2）：像素坐标系；
3）：物理坐标系；
4）：数码照相机坐标系；
5）：理想坐标点在x轴上的畸变量；
6）：理想坐标点在y轴上的畸变量；
7）：相机的焦距；
注：其它所有符号在文中引用处有详细说明。
三、 问题的分析
通常，为了通过物坐标求取像坐标，需知道相机的内外参数，求取内外参数的过程本质就是相机的标定。一般相机的标定方法都采用张氏法标定，该方法能高精度地标定相机，且简便易行，可以说是相机标定方法的代表，但其算法模型仅考虑了镜头径向畸变，其对切向畸变较大的场合，如鱼眼镜头成像，则不适用。且该算法要求一个比较好的初值，若初值选择不当，则算法难以收敛或只能收敛到局部最小，从而大大降低标定精度。
为了克服计算相机内外参数的困难，提出了一种基于两步法的确定像坐标模型。首先采用一种简化的线性模型估算出像中心的近似位置，并利用该像中心点求取初值，由于该近似点畸变很小，故求取的初值能很好地逼近准确值，然后采用一种基于牛顿迭代法得到较为精确的像坐标。
相机在使用阶段其内部参数是基本保持不变的，但两台相机之间的相互位置可能会因不慎磕碰而改变，这样势必会影响测量精度。为了避免在重新确定两台相机之间的相互位置关系时要重新对相机进行标定，本文建立了仅利用空间点在两台相机像面上的坐标即可获得它们之间的相互关系的模型，给出了基于平行线消隐点理论算法。
四、 确定像坐标模型的建立和求解
4.1 数码照相机模型
如图3所示：在模板平面上，先分别以两条正交的直线为轴，并按右手定则建立全局坐标系；然后以数码照相机透视中心为坐标原点，并使其轴和轴分别平行于像面坐标系的轴和轴，且方向相同，轴为数码照相机镜头的主光轴，按右手定则建立数码照相机坐标系；再以像平面上像素坐标值的点为坐标原点，以水平像元阵列方向为方向，以垂直像元阵列方向为方向，建立像素平面坐标系；最后以主光轴与像平面的交点为坐标原点，并使坐标轴分别平行于像素坐标系的坐标轴来建立物理坐标系。模板平面上的理想成像点为，而实际透视畸变后的成像点为。

图3 透视投影模型
4.2 理想透视模型
通过数码照相机采集的数字图像，每幅数字图像在计算机内为数组，行列的图像中的每一个元素我们称之为像素，其数值即为图像点的亮度。再根据针孔成像原理，由世界坐标点到理想像素点的齐次变换如下：
(1)
其中，为内参矩阵，为数码照相机的外参矩阵，和分别为模板相对数码照相机的旋转分量和平移分量。其分别定义如下[2]：

其中分别表示在方向和方向上像点的物理坐标到图像坐标的比例系数；表示主光轴与像平面交点的图像坐标。由于全局坐标系的轴在模板平面上，因此对于模板上的点，，这样式(1)可重写为：
(2)
消去即得到以下理想透视数学模型：
(3)
4.3 实际成像模型
理想的透镜成像是针孔成像，其物与像满足相似三角形关系，而实际的透镜并不满足这个条件(如图3所示)。由于数码照相机光学系统存在加工误差和装配误差，因此物点在数码照相机像面上实际所成的像与理想成像之间主要存在径向畸变和切向畸变两种光学畸变误差。为此，本文主要考虑镜头径向畸变和切向畸变。
设理想成像点的物理坐标为，畸变后的物理坐标为，理想成像点的物理坐标与对应的世界坐标点有如下关系：
(4)
其中为相机焦距，分别是物距和相距，像素点大小（mm/像素）。
由文献[3]，可建立如下总像差模型：
(5)
为径向畸变系数，为切向畸变系数。由理想物理坐标到实际物理坐标的关系如下：

由实际物理坐标到实际像点的关系为：
(6)
综合式式即得到如下实际像坐标计算公式：
(7)
可分别由式、式确定。其中对于数码照相机的每一个位置，其外部参数就有6个；内部参数有以及畸变系数。
通过假设像素坐标的原点对应全局坐标系的原点，像素点为正方形，相机的焦距的理想情形下，得到简化的线性模型[4]，像坐标的初始值有如下关系：
(8)
求取相机内外部参数可参见文献[3]。
根据确定像坐标的模型式可知，像坐标和物坐标关系本质上是一个线性变换加上一个较小的非线性余项，完全满足牛顿迭代法的收敛要求。通过多次迭代可以确定较为精确的像坐标。
4.4 模型应用
针对问题2，通过两步法确定像坐标的模型先由Matlab[5]读入图1、图2（如图4所示），数字化图像。

图4 Matlab中的靶标中圆及其像

然后借助像坐标求取的源程序3（见附录），求出靶标上圆的圆心在像平面上的像坐标。其计算过程如下：
第一步：按照光学透镜原理（如图5）及相似三角形性质有，
图5 透镜成像原理

得到、。
第二步：然后利用简化模型在计算像坐标初值
，
其中取1/3.78像素/毫米，求得。
第三步：反复代入式，通过多次迭代可以确定较为精确的像坐标，结果见表2。

表2 问题2求解结果（单位：像素）
圆形标号
A
B
C
D
E
x坐标
-190.287
-89.4045
133.1423
70.8932
-224.987
y坐标
194.7648
185.9878
176.1031
-124.214
-115.968
注：物坐标原点和像坐标系原点都在图片的几何中心。
4.5 基于Canny算子的中心定位算法的模型检验
圆形标志经透镜成像后畸变为椭圆，为了达到对椭圆中心的高精度定位，首先精确提取标志图像的边缘点，然后对提取的边缘点进行计算从而确定标志中心的精确位置。为了得到椭圆中心的像素坐标，首先用Canny算子对椭圆的边缘进行像素级检测，最后对提取的像素边缘点进行最小二乘椭圆拟合，从而确定标志中心的精确位置，通过高精度的中心点拟合提高相机标定精度。用此算法来检定像坐标的确定模型。
4.5.1 Canny边缘检测基本原理
Canny边缘检测器是高斯函数的一阶导数，是对信噪比与定位之乘积的最优化逼近算子。它能排除低概率的错标非边缘点和低概率不标真实边缘点，使检测出来的边缘点为尽可能地靠近真实边缘中心的单值点。
它先是用高斯滤波器平滑图像，去除图像噪声，然后用一阶偏导的有限差分来计算梯度的幅值和方向，使用的梯度算子计算和方向的偏导数和，方向角，梯度幅值；接着对梯度幅值应用非极大值抑制，用双阈值算法检测并连接边缘；最后，通过对边缘点的标记和领域关系进行连接得到最后的边缘检测图。
4.5.2最小二乘椭圆拟合
对上述算法得到的每个圆区域的亚像素数据进行最小二乘拟合。椭圆的平面方程为：

拟合求得椭圆方程5个参数，椭圆中心坐标计算公式为：

4.6 本检测方法的精度分析
在相机的成像过程中，通常存在以下3种误差：
（1）由相机的镜头所造成的可以描述为一奇次多项式对称径向失真对称的径向失真或桶形失真。
（2）图像中心点偏移失真。
（3）成像面不平整失真。
为了避免以上误差，在仅考虑径向畸变和切向畸变的假设下，利用确定像坐
图6 Canny算子边缘检测靶标圆心的像坐标
标两步法的模型求得靶标上圆的圆心的像坐标。另一方面根据Canny算子边缘检测源程序，提取出靶标上圆的圆心在像平面上的像坐标，两种结果对比如图6所示。
表4 两种算法所得像坐标结果对比（单位：像素）
圆形标号
Canny算子中心定位法
两步法确定像坐标法
误差(%)
Ｕ坐标
Ｖ坐标
Ｕ坐标
Ｖ坐标
Ｕ坐标
Ｖ坐标
A
-189.0056
193.7568
-190.2873
194.7648
0.6781
0.5202
B
-90.0491
187.1209
-89.4045
185.9878
-0.7158
-0.6055
C
129.0463
170.7327
133.1423
176.1031
3.1741
3.1455
D
72.813
-120.7012
70.8932
-124.2142
-2.6366
2.9105
E
-227.2402
-117.3653
-224.9869
-115.9681
-1.0708
-1.1905
注：物坐标原点和像坐标系原点都在图片的几何中心。
由表4可知Ｕ坐标的平均误差：1.6551%，Ｖ坐标的平均误差：1.6754%平均误差：1.6653%；最大误差点为Ｃ点，最小误差点为A点。误差最大3.17%，最小0.52%，平均误差为1.67%，达到了非量测相机的精度要求[6]。
4.7 稳定性分析
由于数码照相机光学系统存在加工误差和装配误差，因此物点在数码照相机像面上实际所成的像与理想成像之间主要存在径向畸变和切向畸变两种光学畸变误差。其中径向畸变是关于数码照相机镜头的主光轴对称的而切向畸变则关于数码照相机镜头的主光轴不对称。虽然实际上还存在，如不对心像差、薄透镜像差等问题，但若考虑过多的畸变参数不仅不能提高计算精度，反而会使计算结果不稳定。为此，本文主要考虑镜头径向畸变和切向畸变。
为了提高稳定性，优化计算求精初值。在计算内外参数过程中，虽可求出

但，不一定正交，实际上这样计算出的旋转矩阵通常不满足正交性，以下将进行初值的优化计算。
可只根据理想透视模型建立如下优化模型的评价标准：
(9)
约束等式是所有都应满足正交性的6个等式。优化目的是由初步求出的对其再求精之。包含了要求取的数码照相机所有参数，其中表示选取的第幅图中心附近的点数，表示第幅图像中第个点的实际图像坐标点，是相对应的世界坐标点，是由代入理想透视模型式(3)所求出的图像坐标点，理论上求出的图像坐标点与实际图像坐标点应越接近越好，为此可建立式的优化模型，为加快收敛可采用梯度下降法。由于其初值较准确，因此梯度法能很快收敛。
4.8 求精数码照相机内外部参数
以上利用图像中心附近点畸变量小的性质，即不考虑图像畸变，而直接选用这些点来求取内外参数的初值。为获取精确的数码照相机模型参数，应考虑所有世界坐标点和对应的图像点，并建立以下优化模型，以求精所有内外参数：
(10)
约束等式是所有都应满足正交性的6个等式。优化模型的目的是为了求精、、、、、、同模型式定义。、、、为本文考虑的畸变系数，由于其一般很小，因此其初始值可置为。本文模型的初值为第1步求取的各参数值以及初始的畸变系数。是由代入式(4)先求得再由实际成像模型式(7)求出的。由于模型参数太多，用一般的梯度法难以保证收敛，因此本文采用一种基于牛顿迭代法求解。
五、 确定两部相机相对位置的数学模型和方法
5.1 模型假设
（1）假设两部相机内部参数固定不变。
（2）假设噪声对相机成像无影响。
5.2 模型建立
设右相机坐标系到左相机坐标系的变换关系为：
(11)
空间任一点在左右象面象点分别为转换成各自相应的相机坐标分别为，；又已知两相机透视中心分别为、，这样由这点便构成一个平面，即外极平面。将这点坐标都转换到左相机坐标系下，由四点共面条件可得：
(12)
令 (13)
将个空间点的左右象点的象面坐标值代入式可得
(14)
其中

由于旋转矩阵为正交矩阵，所以式是一个含有个未知数的非线性齐次方程组。
基于平行线消隐点理论的“两步法”首先利用多于组平行线方向矢量的左右象面消隐点坐标，线性求解一个纯旋转关系的方程组，得到矩阵；然后将矩阵代入式，再利用多于两个空间点的匹配坐标值及左右相机坐标系之间的先验知识，可以线性求解矩阵[7]。
5.3模型求解
1）消隐点原理及矩阵的求解
设在左相机坐标系下单位方向矢量为且过点的空间直线的参数方程为：
(15)
设同一直线在右相机坐标系下的参数方程为：
(16)
其中为此直线在右相机坐标系下的单位方向矢量。根据左右相机坐标系之间的转换关系可知：
(17)
空间平行直线只要不与相机象面平行(即)，经过透视变换后在象面上应相交于一点，即“消隐点”。
对式和式进行透视变换可得：
(18)
消隐点是空间直线(或)时在像平面上的透视投影，由此可得：
(19)
其中()和()分别为左右像面消隐点坐标，和分别为左右相机有效焦距。
设条空间平行直线在左像面透视投影为：
(20)
这组像面直线相交于一消隐点()，同样可求得左像面内的其它消隐点坐标。
同理可得右像面内各消隐点坐标。由式~可得：
(21)
将多于对消隐点的坐标值代入式，并结合矩阵的正交约束条件，解此未知数线性方程组可得旋转矩阵的初值，然后再利用迭代法可更精确地求取矩阵。
2） 矩阵的求解
将矩阵代入式，给定多于两个空间点的左右像面坐标值，并利用两台摄相机之间的先验知识来判断的符号，就可线性求得矩阵。
由双目视觉传感器模型[8]可以得到以下方程：
(22)
由此可得 (23)
将，连同矩阵的值一起代入式，计算两个空间点之间的距离，然后将式算出的含参数的测量值和实际距离值进行比较，便可计算出值，从而求得矩阵。
六、 模型和算法的优缺点分析及改进方向
由于数码照相机模型的参数太多，而且是一个复杂的非线性方程，因此一步求解不仅繁冗，而且稳定性差；本文所采用的两步法确定像坐标和两相机相对位值的方法，不仅考虑了径向畸变，还考虑了切向畸变因素，求得较为准确的像坐标和相机相对位置，精度高，稳定性好，解的收敛速度快。
由于考虑了相机的径向畸变和切向畸变即非线性畸变，计算程序略显复杂，计算量稍大。
由于BP神经网络在逼近能力，收敛速度等方面都优于其它的非线性优化算法，很适宜于相机的标定。将神经网络技术用于相机的标定，不需要预先确定摄像机具体内外参数及有关模型的先验知识，利用BP神经网络对复杂非线性映射关系强大的逼近能力来对视觉系统中的摄像机进行隐式标定，算法标定精度较高，可靠性好，实用性较强，是相机定位的一种趋势。

参考文献
[1] 刘红宁. 基于双目立体视觉拟人机器人定位关键技术的研究[D]. 河北工业大学,2006.
[2] 马颂德,张正友. 计算机视觉[M]. 北京:科学技术出版社,1998.
[3] 毛剑飞,邹细勇,诸　静. 改进的平面模板两步法标定摄像机[J]. 中国图象图形学报,2004,9(7):847~852.
[4] 张　磊,刘小力,罗　翔,田梦倩,史金飞. 一种考虑畸变的摄像机线性标定简化方法[J]. 工业仪表与自动化装置,2007(1):3~7.
[5] 陈杰. Matlab宝典[M]. 北京:电子工业出版社,2007.
[6] 李海滨, 单文军, 刘彬. 双目立体视觉测距系统误差模型的研究[J]. 光学技术,2006(01) :13~16.
[7] 张健新,段发阶,叶声华. ”两步法“求取双目视觉传感器中摄象机位置关系[J]. 光电工程,1998,25(5):37~41.
[8] 李雪鹃, 蔡怀宇, 黄战华. 基于双目视差的测距系统设计[J]. 电子测量技术,2007(08) :31~34.

附 录
基于Canny算子边界检测源程序1
%基于Canny算子求边界坐标
%参数name为图片名
function [x,y]=tp(name)
I=imread(name);%读入图片
I=double(I);
I=I(1:768,1:1024);
G=edge(I,'Canny'); %基于Canny算子求边界
imshow(G);
I=uint8(I);
G=uint8(G);
subplot(1,2,1)
imshow(I);%输出图像
subplot(1,2,2)
imshow(G);%输出图像
x=[];
y=[];
for i=100:600 %求边界坐标
for j=200:700
if G(i,j)==0）
y=[y,i];
x=[x,j];
end
end
end

基于简化模型求取像坐标初值源程序2
%基于简化模型求解像圆心
%参数x，y为物的圆心坐标u为物距
function [x0,y0]=jiandan(x,y,u)
x0=[];
y0=[];
v=1577;%给像距赋初值
h1=24;
f=ff(h1,h2,v) %调用求解焦距的函数，
n=length(u);
for j=1:n;
k=f/(u(j)+v)*3.78;%求比例系数
a=x(j)*k;%求像坐标并保存
b=y(j)*k;
x0=[x0,a];
y0=[y0,b];
end

两步法求解像坐标源程序3
%两步法求解像坐标
%读入数据并整理成可处理的形式
%M为靶标数据
%m为靶标的像数据
%求解焦距
function f=ff(h1,h2,v)
n=length(h2); for k=1:n
u(k)=h1*v/h2(k);%求每一个像高对应的物距
end
u=u*3.78;
for k=1:n
F(k)=1/(1/v+1/u(k));%求每一个物距对应的焦距
end
f=mean(F);%求所以焦距的算术平平均
clear all;close all;
M=imread('tu.bmp');%读入图片
m1=imread('tu1.bmp');
m2=imread('tu2.bmp');
m3=imread('tu3.bmp');
m4=imread('tu4.bmp');
m5=imread('tu5.bmp');
MM=[M(:,1:2);M(:,3:4);M(:,5:6);M(:,7:8)];%提取数据
mm1=[m1(:,1:2);m1(:,3:4);m1(:,5:6);m1(:,7:8)];
mm2=[m2(:,1:2);m2(:,3:4);m2(:,5:6);m2(:,7:8)];
mm3=[m3(:,1:2);m3(:,3:4);m3(:,5:6);m3(:,7:8)];
mm4=[m4(:,1:2);m4(:,3:4);m4(:,5:6);m4(:,7:8)];
mm5=[m5(:,1:2);m5(:,3:4);m5(:,5:6);m5(:,7:8)];
M=[M(:,1:2) ; M(:,3:4) ; M(:,5:6) ; M(:,7:8)];
m1=[m1(:,1:2) ; m1(:,3:4) ; m1(:,5:6) ; m1(:,7:8)];
m2=[m2(:,1:2) ; m2(:,3:4) ; m2(:,5:6) ; m2(:,7:8)];
m3=[m3(:,1:2) ; m3(:,3:4) ; m3(:,5:6) ; m3(:,7:8)];
m4=[m4(:,1:2) ; m4(:,3:4) ; m4(:,5:6) ; m4(:,7:8)];
m5=[m5(:,1:2) ; m5(:,3:4) ; m5(:,5:6) ; m5(:,7:8)];
M=M';
m(:,:,1)=m1';
m(:,:,2)=m2';
m(:,:,3)=m3';
m(:,:,4)=m4';
m(:,:,5)=m5';
%M为靶标数据
%m为靶标的像数据
function liangbufa(M,m)
%M=[X,Y]’转化为M=[X,Y,1]'
%m=[u,v]'转化为m=[u,v,1]'
[rows,npts]=size(M);
matrixone=ones(1,npts);
M=[M;matrixone];
num=size(m,3);
for i=1:num
m(3,:,i)=matrixone;
end
H=jiandan(x,y,u)%调用简单模型求解初值
V=[];
for flag=1:num %迭代
v12(:,:,flag)=[H(1,1,flag)*H(2,1,flag),H(1,1,flag)*H(2,2,flag)+H(1,2,flag)*H(2,1,flag),H(1,2,flag)*H(2,2,flag),H(1,3,flag)*H(2,1,flag)+H(1,1,flag)*H(2,3,flag),H(1,3,flag)*H(2,2,flag)+H(1,2,flag)*H(2,3,flag),H(1,3,flag)*H(2,3,flag)];
v11(:,:,flag)=[H(1,1,flag)*H(1,1,flag),H(1,1,flag)*H(1,2,flag)+H(1,2,flag)*H(1,1,flag),H(1,2,flag)*H(1,2,flag),H(1,3,flag)*H(1,1,flag)+H(1,1,flag)*H(1,3,flag),H(1,3,flag)*H(1,2,flag)+H(1,2,flag)*H(1,3,flag), H(1,3,flag)*H(1,3,flag)];
v22(:,:,flag)=[H(2,1,flag)*H(2,1,flag),H(2,1,flag)*H(2,2,flag)+H(2,2,flag)*H(2,1,flag),H(2,2,flag)*H(2,2,flag),H(2,3,flag)*H(2,1,flag)+H(2,1,flag)*H(2,3,flag),H(2,3,flag)*H(2,2,flag)+H(2,2,flag)*H(2,3,flag), H(2,3,flag)*H(2,3,flag)];
V=[V;v12(:,:,flag);v11(:,:,flag)-v22(:,:,flag)];
end
k=V'*V;%将Vb=0转化为kb=V'*V*b=0,这样k就是一个对称正定矩阵，
[u,v,d]=svd(k);
[e,d2]=eig(k);
b=d(:,6);
v0=(b(2)*b(4)-b(1)*b(5))/(b(1)*b(3)-b(2)^2);
s=b(6)-(b(4)^2+v0*(b(2)*b(4)-b(1)*b(5)))/b(1);
alpha_u=sqrt(s/b(1));
alpha_v=sqrt(s*b(1)/(b(1)*b(3)-b(2)^2));
skewness=-b(2)*alpha_u*alpha_u*alpha_v/s;
u0=skewness*v0/alpha_u-b(4)*alpha_u*alpha_u/s;
A=[alpha_u skewness u0
0 alpha_v v0
0 0 1];
D=[];%求解所有外部参数
d=[];
Rm=[];
for flag=1:num
s=(1/norm(inv(A)*H(1,:,flag)')+1/norm(inv(A)*H(2,:,flag)'))/2;
rl1=s*inv(A)*H(1,:,flag)';
rl2=s*inv(A)*H(2,:,flag)';
rl3=cross(rl1,rl2);
RL=[rl1,rl2,rl3];
[U,S,V] = svd(RL);
RL=U*V';
TL=s*inv(A)*H(3,:,flag)'; TL是位移矩阵t(外参)
RT=[rl1,rl2,TL];%H=A[r1 r2 t]
XY=RT*M;% M是model plane 点的坐标
UV=A*XY;% sm=A[R t]M,UV是等式的右边
UV=[UV(1,:)./UV(3,:); UV(2,:)./UV(3,:); UV(3,:)./UV(3,:)];
XY=[XY(1,:)./XY(3,:); XY(2,:)./XY(3,:); XY(3,:)./XY(3,:)];
end
k=inv(D'*D)*D'*d;
para=[Rm,k(1),k(2),alpha_u,skewness,u0,alpha_v,v0];
options = optimset('LargeScale','off','LevenbergMarquardt','on');% 用最小二乘法求解
[x,resnorm,residual,exitflag,output] = lsqnonlin( @simon_HHH, para, [],[],options, m, M);
k1=x(num*6+1)% 显示结果
k2=x(num*6+2)
A=[x(num*6+3) x(num*6+4) x(num*6+5); 0 x(num*6+6) x(num*6+7); 0,0,1]

基于Canny算子边界检测确定圆心源程序4
%基于Canny算子求出的边界坐标来拟合圆心
%参数x，y为边界坐标
function [x0,y0]=yuanxin(x,y)
x1=(min(x)+max(x))/2;%估计圆心的坐标
y1=(min(y)+max(y))/2;
X=[];
Y=[];
for a=1:length(x)-2 %任意取三个边界点确定一个圆心的坐标
for b=a+1:length(x)-1
for c=b+1:length(x)
x0=((x(b)^2-x(a)^2+y(b)^2-y(a)^2)*(y(b)-y(a))/2-(x(a)^2-x(c)^2+y(a)^2-y(c)^2)*(y(a)-y(b))/2)/((x(a)-x(b))*(y(a)-y(c))-(x(a)-x(c))*(y(a)-y(b)));
y0=((x(a)^2-x(c)^2+y(a)^2-y(c)^2)*(x(a)-x(b))/2-(x(a)^2-x(b)^2+y(a)^2-y(b)^2)*(x(a)-x(c))/2)/...
((x(a)-x(b))*(y(a)-y(c))-(x(a)-x(c))*(y(a)-y(b)));
if (x1-x0)^2+(y1-y0)^2<100
X=[X,x0];
Y=[Y,y0];
end
end
end
end
x0=mean(X);%求满足条件的坐标均值
y0=mean(Y);