基于靶标平面平行于图像平面的相机定位模型

基于靶标平面平行于图像平面的相机定位模型
摘要
通过对相机内部结构的分析及题目给定的条件，本文首先利用小孔成像原理建立了线性模型。

为了迅速有效的进行相机标定，我们对线性标定方法进行了改进，引入畸变因子，建立了非线性畸变模型，由坐标系之间的转换，得到了实际成像点(,)d d x y 与世界坐标系中物点()w ,,w w X Y Z 的相互关系如下：
2
12
1()(1)(cos sin )
(y )(1)(sin cos )
d w d w d z w w x d w d w d z w w y x X x Y x t k r f X Y t y X y Y t k r f X Y t βγααβγαα⎧-+++=-+⎪⎨-+++=-+⎪⎩ 该模型基于靶标平面平行于图像平面建立，从而减少了内外参数的引入，降低了求解的复杂度。

然后应用matlab7.1的图像处理功能，导入靶标的像图，利用提取黑像素法获取像坐标。

可得
的相机位置固定，因此0w Z ≠，从而造成计算结果较大的误差。

本文从模型和数值计算两方面分析了产生误差的原因，同时提出了检验精度和稳定性的方法。

模型确定后，可求得两个相机各自的像坐标，再假设两相机的相对距离确定且为B ，分下面两种方式建立相应的模型。

① 左右相机的中心点和左右像点构成的光线相交于一点。

利用几何相似三角形原理建立含有四个未知参数的数学模型：
1
12
,w x Bx x x =
-121212,w y By By x x x x =
=--12.w Bf z x x =-通过B 来确定两个相机的相对位置 ② 左右相机的中心点和左右像点构成的光线并不一定相交于一点，而是两条异面
直线。

对模型点的空间位置进行求解,基于像坐标的基础上,利用异面直线公垂线中点去逼近模型的空间点.建立的相关的模型：
121112220
0Q Q PQ Q Q P Q ∙=∙=⎧−−−→−−−→⎪⎨
−−−→−−−→⎪⎩即(()())(()())1112221111222200P a l P a l l P a l P a l l ⎧+-+⋅=⎪⎨+-+⋅=⎪
⎩ 关键字：小孔透视模型 非线性畸变模型 相机标定 双目定位 公垂线法
一、问题重述
1.背景
数码相机定位在交通监管等方面应用广泛。

即用数码相机摄制物体的相片确定物体表面某些特征点的位置。

双目定位：对于物体上一个特征点，两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。

只要知道两部相机精确的相对位置，就可用几何的方法，坐标转换确定特征点的位置。

对于双目定位，精确的确定两部相机的相对
位置的过程就叫系统标定。

标定的一种方法：在
一块平板上画若干个点，同时用这两部相机照相，
分别得到这些点在它们像平面上的像点，利用这
两组像点的几何关系求得两部相机的相对位置。

然而，无论在物平面或像平面上都无法直接得到
没有几何尺寸的“点”。

实际的做法是在物平面上
画若干个圆（称为靶标），它们的圆心就是几何的
点。

而它们的像一般会变形，如图1.1所示，所
以必须从靶标上的这些圆的像中把圆心的像精确
地找到，标定就可实现。

图 1.1 靶标上圆的像在如下靶标上，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，作为标定的特征点12mm为半径作圆。

以AC边上距A点30mm处的B为圆心，12mm为半径作圆，如图1.2所示。

图1.2 靶标示意图
用一像距（光学中心到像平面的距离）为1577个像素单位，分辨率为1024×768的数码相机。

固定位置摄得像如图1.3所示。

图1.3 靶标的像
2.提出问题
（1）以相机的光学中心为坐标原点，x-y平面平行于像平面，建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标；
（2）由图1.2和图1.3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标；
（3）设计算法检验模型的精度和稳定性
（4）借助此靶标进行系统标定，建立数学模型及方法。

二、问题分析
相机定位中，最主要的就是坐标系转换[1]，因此先建立所需坐标系。

问题一为确定空间物体表面某点通过相机投影之后，所在空间的三维几何位置与像平面所成的像点之间的关系，对相机成像原理加以分析，设定内外参数，采用“两步法”标定，首先忽略畸变因子建立线性模型，求得理想成像坐标，再考虑畸变，建立非线性模型。

问题二应用matlab的图像处理功能，编程读取靶标上圆在像平面上的坐标值。

以所读取的点为特征点，返回模型对非线性模型求解相机内外参数，从而进行精度和稳
定性的检验。

图2.1 成像步骤
双目测定法标定物坐标的位置,然后再通过物坐标反求出另一部相机的位置.因此只要假定两个机的相位位置固定,又由于之前求出的像坐标,利用几何关系求解其中的参数值,从而确定了另外一个相机的相对位置。

三、基本假设
标定模板平面与图像平面平行
双目测量所需的两数码相机型号一致。

左摄站和左像点构成的光线及右摄站和右像点构成的光线相交于一点。

假定I 和Ⅱ在同一个平面内，2个CCD 的像距1o 1o ＝2o 2o =f 。

四、符号说明
m ：物距（物体与光学中心的距离） n ：像距（像平面到光学中心的距离） f ：焦距（焦点到光学中心的距离）
w P ()w ,,w w X Y Z ：世界坐标系中的物点
(),,c c c c P X Y Z ：摄像机坐标系中的特征点
)(u,v ：图像像素坐标系中的像点 )(x,y ：图像平面坐标系中的像点
s ：横纵坐标像素比例系数， R ：旋转矩阵 1k ：畸变因子
五、模型的建立与求解
坐标系的建立
图像平面坐标系：在图5.1中，以坐标轴x,y 组成图像平面直角坐标，以光轴为z 轴，光轴z 与平面x-y 的交点为坐标原点1o 建立坐标系1o xyz -。

图像像素坐标系：根据数码相机的存储原理[2]，在图5.1中，以图像平面直角坐标系x-y 的左上角为坐标原点o ，以存储单元（像素）为单位建立坐标轴为u,v 的图像像素坐标系o-uv 。

即)(u,v 是以像素为单位的图像坐标系坐标，每一像素的坐标)(u,v 分别是该像素在数组中的列数和行数。

摄相机坐标系：在图5.1中，以光学中心为坐标原点c z ，以坐标轴c c c x ,y ,z 组成摄
像机坐标系。

其中z 轴为光轴，c c x -y 平面平行于图像平面x-y 。

世界坐标系：以一个基准坐标系描述摄像机安放在现实世界的位置，并用它描述三维环境中任何物体的位置。

在图5.1中，由基准观测原点和,,w w w X Y Z 轴组成世界坐标系w P ()w ,,w w X Y Z 。

图5.1坐标系示意图
问题一
模型I ：小孔透视理想模型
由相机的成像原理，即透镜的成像原理，如图5.2所示，其满足：
111f n m
=+ 在相机构造中，图像成像于焦点所在位置的CCD ，物距m 又远远大于焦距f ，即
,n f m f ≈>>，因此可将透镜模型简单近似为小孔成像模型。

图5.2小孔成像透视原理
小孔透视模型是一种最常用的理性模型，其成像关系是线性的。

摄像机的成像只与三维物体点的空间位置、摄像机焦距以及物体或摄像机相对运
动参数有关，即如图2所示，摄像机成像模型的建立从世界坐标轴中的（Xw,Yw,Zw ）到图像像素坐标轴(u ，v)的过程。

三维空间中的物点投影到图像平面上，有如下相互关系：
(1) 利用几何学三角形相似原理，摄相机坐标系和图像平面坐标系变换存在如下关系式：
c c c
c X x f Z Y y f Z ⎧=⎪⎪
⎨
⎪=⎪⎩
（5.1） (,)x y 用齐次坐标表示为(,,1)x y ，用矩阵表示(5-1)为：
000010
1c c c x f
X y f
Y Z λ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
（5.2） 其中λ为常数因子，c Z λ=。

（2）图像平面坐标系(x,y)和图像像素坐标系(u,v)之间存在如下关系式：
00x y u f x u v f y v =+=+ （5.3）
其中
/1/x y
f s dx f dy =⎧⎪⎨
=⎪⎩ 其中s 为横纵坐标像素比例系数，,x y d d 为像素点在图像平面坐标系x,y 轴方向上的宽度，00(,)u v 为图像坐标的原点在图像平面上的坐标值。

用矩阵表示为：
00000
011x y f u x u f v y v ⎡⎤⎡⎤
⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
（5.4） （3）世界坐标系w P 和摄相机坐标系c P 的转换，引入旋转矩阵R 和三维平移向量t ，有如下关系：
c w P RP t =+ （5.5） 其中t 为三维平移向量，R 为33⨯的正交单位矩阵[3]，R 用欧拉角表示如下。

cos cos cos sin sin sin cos cos sin cos sin sin sin cos sin sin sin cos cos sin sin cos cos sin sin cos sin cos cos R αβαβγαγαβγαγαβαβγαγαβγαγββγβγ-+⎡⎤
⎢⎥=+-⎢⎥⎢⎥-⎣⎦ （5.6）
其中,,αβγ为欧拉角。

用矩阵表示为
'0111c w x c w y c w z X X t Y R t Y t t Z Z t ⎡⎤⎡⎤
⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦
⎣⎦⎣⎦
， (5.7)
（4）世界坐标系和图像坐标系之间的变换关系： 由（5-2），（5-4），（5-7）进行矩阵运算可得 '1000010001100101w w w X u R t Y v A Z λ⎡⎤⎡⎤⎡⎤⎢⎥
⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥
⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
⎣⎦
（5.8） 由于对A 的影响因子均为相机基本构造，称A 为相机内部参数矩阵，结构为
000
1x y
f s u A f v ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦
模型Ⅱ 畸变影响的非线性模型
由于相机光学系统的加工误差和装配误差，在很多情况下摄像机并不是理想的针孔透视模型，会使实际成像产生畸变。

因此需要在线性模型的基础上，加以改进建立非线性模型。

畸变主要为径向畸变和切向畸变，而径向畸变为产生误差的主要因素，在此仅考虑径向畸变。

则畸变量
22
122
1()
()
x d d d y d d d k x x y k y x y δδ⎧=+⎪⎨=+⎪⎩ （5.9） (,)d d x y 在畸变影响模型下的坐标。

理想坐标(,)x y 即为实际成像点(,)d d x y 与畸变量之和。

22122
1()
()
d x d d d d d y d d d d x x x k x x y y y y k y x y δδ⎧=+=++⎪⎨=+=++⎪⎩ （5.10） 而（5-12）即为要求解的模型，其中的相机内外部参数需要标定。

模型求解：
（1）相机内外参数的求解
因为标定控制点共面，可取0w Z =,由于控制点所在平面与传感器平面平行，则
cos 1,cos 1,sin ,sin βγββγγ≈≈≈≈，代入式（5.6）有
cos sin cos sin sin cos sin cos 1
R αααβαγα
α
αβαγβ
γ
-+=--⎡⎤
⎢⎥⎢⎥
⎢⎥⎣⎦
（5.11）
由（5.7）（5.10）（5.11）计算可得
22
1221cos sin ()sin cos ()w w x d d d
d w w z w w y
d d d d w w z X Y t x k x x y f X Y t X Y t y k y x y X Y t ααβγααβγ-+++=-++++++=-++⎧⎪⎪
⎨
⎪⎪⎩
（5.12） 由（5-12）中两式相比得
cos sin sin cos d w w x
d w w y
x X Y t y X Y t αααα-+=
-+ （5.13） 假定世界坐标系的原点远离相机坐标系的Z 轴，则,x y t t 不同时为0，设0y t ≠，有：
111111
11cos sin sin cos 0cos []sin y w y w d y d x w d y w d y d y w d w d w d w d d y d
x y
Y t Y y t y t t X x t Y x t x t X y Y y Y y X x y t x t t αααααα--------⎧-+---=⎪
⎪⎡⎤
⎨⎢⎥
---=⎪⎢⎥⎪⎢⎥⎣⎦⎩
（5.14） 令1
112113cos sin y y x y C t C C t C t t αα----⎡⎤
⎡⎤⎢⎥⎢⎥==⎢
⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
（5.15） 因为0y t ≠，所以12,C C 不可能同时为0，若10C ≠，则
2
31
1
arctan
cos
y x y C t t C t C C α
α=== （5.16）
再进行转换，确定欧拉角,,z t βγ及畸变系数1k 由（5.12）式可得
2
1
2
1()(1)(cos sin )(y )(1)(sin cos )d w d w d z w w x d w d w d z w w y x X x Y x t k r f X Y t y X y Y t k r f X Y t βγααβγαα⎧-+++=-+⎪⎨-+++=-+⎪⎩ （5.17） 即为参数模型，代入特征点坐标解非线性方程组即可求得。

问题二
应用matlab 读取图像平面上所成图像中心的像素坐标[4,5]。

采用imread 函数读取所给靶标的像，用im2bw 函数将图形转换为二值图像，如图5.3：
图5.3经过matlab 处理后的二值图像
依据像点的RGB 值不同，搜索图像获得像点()'''''A B C D E ，，，，的像素坐标值，取中心点为像点上所有点在x ，y 轴上的坐标值的中值。

为避开图像平面上多个像点造成的存储不便，依次截取像点所在的区域。

以A 为例，截取图像A 的区域以（240，130）为起始点，（x ，y ）长分别为（135，110）。

源代码见附录一：
经过坐标平移，可得A 点在图像像素坐标系的中心'A 坐标为（319.6580，190.6814），同理，可得其他点在图像像素坐标系上的中心，如表5.1
问题三
据坐标系的建立关系，由几何关系有
()00,(0,768)u v =，相机提供的分辨率，可得1024768s =，由于
1577417.1958mm f n ≈==像素
经过坐标变换，可得在像平面坐标系中的坐标分别为
'(239.7435dx,-577.3186dy)A ,'(318.0056dx, -570.2056dy)B , '(480.6746dx,-553.5478dy)C ,'(437.7997dx,-264.0176dy)D , '(214.2493dx,-265.2269dy)E
所有表达式整合，为matlab 求解方便，令cos sin m n αα==，
，满足 22sin cos 1αα+=，具体为：
()()22122
12
2()(1)()
()(1)()1d w d w d z d d w w x d w d w d z d d w w y x X x Y x t k x y f X m Y n t y X y Y y t k x y f X n Y m t m n βγβγ⎧-++++=-+⎪⎪-++++=-+⎨⎪+=⎪⎩
在此靶标及其所成像中，世界坐标系以左下角E 点为坐标原点，据之前的坐标系的确定及假设建立靶标所在的平面坐标系俯视图如图5.4：
图5.4 世界坐标系俯视图
由题所给条件有(0,100),(30,100),(100,100),(100,0),(0,0)w w w w w A B C D E ，编写matlab 程序求解。

源代码见附录二。

调运程序得到 x =
-0.00000000094455 -0.01116086703969 -0.00468153802128 -0.02285754934409 -63.09276196630432 -52.57710142718334 18.65170644943549 0.00000003082829 0.06702869359774 0.99775104822763
-0.00000011629373 显然，所求解明显错误，，分析产生错误的原因。

（1）模型求解过程中，有假设控制点所在平面与传感器平面平行，而题中所提供的靶标与传感器平面不一定平行；
（2）数值计算过程中，多次取近似值，造成误差。

按照标定方法。

借助题中所给定的靶标，及matlab 所确定的靶标上圆的圆心在像平面和世界平面中的坐标，求解得到非线性模型中的相机内外参数值，求得参数后，添加标注点，利用所求模型求出目标坐标，与实际坐标比较分析，同样可进行数据模拟或实验，以检验精度和稳定性。

问题四
模型Ⅰ 双目视觉测量原理模型：
双目视觉直径测量是模拟人类眼睛对物体直径进行测量的方法，由两个不同的相
机对同一物体的成像差异来恢复被测物体的三维尺寸[6]。

有双目视觉原理可得图5.5
图5.5双目视觉原理图
图中I 和Ⅱ分别为2个CCD 成像平面，1z 和2z 分别为I 和Ⅱ的光轴，1o 和2o 分别为I 和Ⅱ的光心，1o 和2o 分别为2个CCD 的中心。

如图所示，假定在同一个平面内的I 和Ⅱ分别为2个CCD 成像平面且，2个CCD 的像距1o 1o ＝2o 2o =f ，1o 到2o 的距离1o 2o 为B ，由双目视觉原理可得如下坐标系：
（1）以I 中CCD 的中心1o 为原点，I 的光轴为1z 轴，以xy 的直角坐标系得
11O xyz 的现实世界坐标系；
（2）以I 中CCD 的中心1o 为原点，建立1o 1x 1y 成像平面I 坐标系； （3）以II 中CCD 的中心2o 为原点，建立2o 2x 2y 成像平面Ⅱ坐标系。

P 为现实世界中的一点，在I 中的成像点为1P ，在1o 1x 1y 坐标下的坐标为)(1,1x y ；在Ⅱ中的成像点为2P ，在2o 2x 2y 坐标系下的坐标为)(2,2x y 。

已知)(1,1x y 和)(2,2x y 可以求解出P 点在现实世界中的坐标为()w ,,w w X Y Z
利用相似三角形原理，P 点的坐标可以由以下的表达式求得：
1
12
12121212
,
,.w w w x y Bx x x By By x x x x Bf
z x x ==⎧⎪-⎪⎪=⎨--⎪⎪=⎪-⎩
借助已经确定的相机内外参数，像坐标的模型()11x ,y 和()22x ,y 可转化为含
()w ,,w w X Y Z 的参数表达式，那么利用上述的四个方程就可以解出w x ，w y ，w z 及B 的
值，即确定了两数码相机的相对位置。

模型Ⅱ
双目视觉测量中，为了求取线激光投影到物体上的投影点所在的空间坐标，需要将物体上的一个投影点在两个CCD 摄像机同一时刻所拍摄的两副图像上进行匹配即是将一副激光条纹图像上的目标像素点与另一副图像的目标像素点对准。

同名像点进行匹配时，左摄站和左像点构成的光线及右摄站和右像点构成的光线，由于摄像机成像模型并非理想的小孔成像模型以及CCD 成像是噪音的影响和标定参数的微小误差。

实际上这两条光线在空间可能无交点，因此由两条光线构成的方程组可能无解，因此交会时无法确定模型点的空间坐标。

为了确定该模型点的具体空间位置，我们采用公垂线法。

所谓公垂线是同时垂直于两异面直线又同时与该两异面直线有交点的直线。

这样就解决了对应像点匹配有误差时模型点的空间坐标确定间题。

由于对应像点的匹配有误差，采用上述方案确定的模型点坐标位置与其真实空间坐标位置是存在误差的，我们将这一方案称为公垂线法[7]。

由于误差的影响，这两条线可能不会相交与一个点上，而且由于相机位置的影响，这两条线不可能平行，所以这两条直线是异面直线。

(1)当公垂线的长度为零，那么两条光线相交于同一点。

(2)当公垂线的长度较小时，那么模型点的位置取公垂线的中点。

(3)当公垂线的长度较长时，那么投射点在图像匹配是误差太大，所求点可以忽略。

模型的求解
图2 空间点坐标求取模型
如图2所示，假设两条异面直线的1L 和2L 的单位向量为1l 、2l ，1Q 、2Q 为在两条异面直线上滑动的点，当12Q Q 值最小时，12Q Q 即为两条异面的公垂线，()122Q Q +为
公垂线中点的坐标。

由于异面直线的公垂线与异面直线垂直容易得到：
121112220
0Q Q PQ Q Q P Q ∙=∙=⎧−−−→−−−→⎪⎨
−−−→−−−→⎪⎩即(()())(()())1112221111222200P a l P a l l P a l P a l l ⎧+-+⋅=⎪⎨+-+⋅=⎪⎩
求解方程得：
)()())((
))((
)))(())(
())((
21121212112122111221221212p p l p p l l l a 1l l l l p p l l l p p l a 1l l l l ⎧-⋅--⋅⋅⎪=⎪-⋅⋅⎪
⎨
-⋅⋅⋅--⋅⎪
=
⎪-⋅⋅⎪⎩
将求出的位置系数1a 、2a 代入到公垂线所在的直线方程中，得到12Q Q 的大小，如果12Q Q 小于假设某个特定的阀值，就可以近似的认为空间的坐标为()12Q Q +。

当空间坐标确定以后在由上述给定的四个表达式中的任何一个表达式都可以求出B 的
值。

六、模型评价与改进
问题一，采用的两步法，其中的线性标定方法，不考虑即便因子的影响，求解线性方程组，取得理想坐标，为非线性方程组的迭代提供初值，考虑即便因子的非线性模型，能较为准确的求得相机的内外参数，在一定程度上弥补了线性模型的不足；但是在此，对于畸变因子，仅考虑径向畸变，且为一次。

改进：由于径向畸变和切向畸变均影响成像，因此同时考虑多个径向畸变和切向畸变因子，改进模型如下：
径向畸变量为
24
1224
12(...)
(...)x d y d x k r k r y k r k r δδ⎧=++⎪⎨=++⎪⎩ 切向畸变量
()()
'2212'
22
212222x
y d d d d d d p x y p r x p x y p r y δδ⎧=++⎪⎨=++⎪⎩ 最终模型为：
()()
'24221212'2422
1221(...)22(...)22x
y d x d d d d d d y d d d d x x x x k r k r p x y p r x y y y k r k r p x y p r y δδδδ⎧=++=++++++⎪⎨=++=+++++⎪⎩
异面直线公垂线中点逼近待测点的方法，能够提高测量数据的精度。

减少噪声点，因此是恢复模型点空间位置的一个比较好的解决办法，但是它并没有考虑双目视觉中匹配误差的统计分布特征，因此公垂线法不是统计意义上误差最小的空间位置恢复方案，只有考虑匹配误差的统计分布特征，才能提出统计上最优的模型点空间坐标恢复方案。

七、模型推广
利用双目定位原理，如每幅图像中相对静止的三个标志为参考，实时计算目标图像的雅可比短阵，从而预测出目标下一步运动方向，实现了对动方式未知的目标的自适应跟踪。

该系统仅要求两幅图像中都有静止的参考标志，无需摄像机参数。

而传统的视觉跟踪伺服系统需事先知道摄像机的运动、光学等参数和目标的运动方式。

双目立体视觉和机器人整体姿态信息集成，开发了仿真机器人动态行长导航系统。

该系统实现分两个步骤：首先，利用平面分割算法分离所拍摄图像对中的地面与障碍物，再结合机器人身体姿态的信息，将图像从摄像机的二维平面坐标系转换到描述躯体姿态的世界坐标系，建立机器人周围区域的地图；基次根据实时建立的地图进行障碍物检测，从而确定机器人的行走方向。

参考文献
[1] 马颂德,张正友，计算机视觉[M]，北京:科学出版社，1998。

[2]（美）赫恩（Hearn,D.）等著，蔡士杰，宋继强，蔡敏译，计算机图形学（第三版）[M]，北京:电子工业出版社，2005。

[3] M A Penna, Determining camera parameters from the perspective projection of a quadrilateral[J],Pattern Recognition,1991,24(6):533-541
[4] 王沫然，MATLAB与科学计算（第2版）[M]，北京:电子工业出版社，2003。

[5] 苏金明，王永利，MATLAB7.0实用指南（下册）[M]，北京:电子工业出版社，2004。

[6]李德良,刘艳玉,张飞龙,尹志勇，双目视觉测径仪的单CCD实现[J]，军械工程学院学报，2007，19（5）：55-56
[7] 罗世民，李茂西，双目视觉测量中三维坐标的求取方法研究[J]，计算机工程与设计，2006，27（19）：3622-3624。

附录一：
I=imread('c.bmp');%读取原图像
BW=im2bw(I);%将图像转化为二值图像
I2=imcrop(BW,[240 130 135 110]); %截取区域
dim=size(I2);
m = zeros(1,3);
k=0;%记录像点个数
x=0;
y=0;
i=1;
j=1;
%循环搜索成像的像素点，根据其RGB进行判断
for i=1:1:dim(1)
for j=1:1:dim(2)
m=impixel(I2,i,j);%读取(i,j)点的RGB值
if m(1)==0
x=i+x;
y=j+y;
k=k+1;
end
end
end
centerx=x/k %中心坐标的x值
centery=y/k %中心坐标的y值
运行：
centerx =79.6580，centery =60.6814
附录二
myfun.m
function F=myfun(x)
%dx-x(1),a-x(2),b-x(3),k1-x(4),tx-x(5),ty-x(6),tz=x(7),dy=x(8),sin-x(9),co s
%-x(10),f-x(11)
F=[(-239.7435*x(1)*100*x(3)+239.7435*x(1)*x(7))*(1+x(4)*((239.7435*x(1))^2 +(-577.3186*x(8))^2))-x(11)*(-100*x(9)+x(5));
(577.3186*x(1)*100*x(3)-577.3186*x(1)*x(7))*(1+x(4)*((239.7435*x(1))^2+(-5 77.3186*x(8))^2))-x(11)*(-100*x(10)+x(6));
(-318.0056*x(1)*30*x(2)+318.0056*x(1)*100*x(3)+318.0056*x(1)*x(7))*(1+x(4) *((318.0056*x(1))^2+(-570.2056*x(8))^2))-x(11)*(30*x(10)-100*x(9)+x(5));
(570.2056*x(1)*30*x(2)-570.2056*x(1)*100*x(3)-570.2056*x(1)*x(7))*(1+x(4)* ((318.0056*x(1))^2+(-570.2056*x(8))^2))-x(11)*(30*x(9)-100*x(10)+x(6));
(-480.6746*x(1)*100*x(2)+480.6746*x(1)*100*x(3)+480.6746*x(1)*x(7))*(1+x(4 )*((480.6746*x(1))^2+(-553.5478*x(8))^2))-x(11)*(100*x(10)-100*x(9)+x(5));
(553.5478*x(1)*100*x(2)-553.5478*x(1)*100*x(3)-553.5478*x(1)*x(7))*(1+x(4) *((480.6746*x(1))^2+(-553.5748*x(8))^2))-x(11)*(100*x(9)-100*x(10)+x(6));
(-437.7997*x(1)*100*x(2)+437.7997*x(1)*x(7))*(1+x(4)*((437.7997*x(1))^2+(-264.0176*x(8))^2))-x(11)*(100*x(10)+x(5));
(264.0176*x(1)*100*x(2)-264.0176*x(1)*x(7))*(1+x(4)*((437.7997*x(1))^2+(-2 64.0176*x(8))^2))-x(11)*(100*x(9)+x(6));
(214.2493*x(1)*x(7))*(1+x(4)*((214.2493*x(1))^2+(-265.2269*x(8))^2))-x(11) *x(5);
(-265.2269*x(1)*x(7))*(1+x(4)*((214.2493*x(1))^2+(-265.2269*x(8))^2))-x(11 )*x(6);
x(9)^2+x(10)^2-1]
命令程序
qiujie.m
x0=[0.001;-0.12;-0.05;0.00002;-60;-50;200;0.0005;0.05;0.99;-250]
options=optimset('Display','iter');
[x,fval]=fsolve(@myfun,x0,options)。