本发明属于计算机视觉领域，涉及一种利用空间中三个球求解中心折反射摄像机内参数的方法。
背景技术：
计算机视觉是指通过计算机、图像采集等设备来模拟人类的视觉对目标进行识别、检测、跟踪等。计算机视觉的中心任务就是对图像进行获取，处理，分析。而它的最终目标是使计算机具有通过二维图像认知三维世界的能力。在视觉测量中，计算机从摄像机获得二维图像信息来计算三维空间中物体的几何信息、并由此重建和识别物体，而空间物体表面点的三维几何位置与在图像中的对应点之间的相互关系是由摄像机成像的几何模型决定的，这些的几何模型就是摄像机的内参数，而确定这些参数的方法就叫做摄像机的标定。摄像机标定可分为传统标定、自标定和基于几何实体的标定。无论哪种标定方法，都旨在建立二维图像与摄像机内参数之间的约束关系，特别是线性约束关系，这是目前摄像机标定所追求的目标，也是目前计算机视觉领域研究的热点之一。
文献“geometricpropertiesofcentralcatadioptriclineimagesandtheirapplicationincalibration”,(barretoj.p.,araujoh.,ieeetransactionsonpatternanalysisandmachineintelligence,vol.27,no.8,pp.1327-1333,2005)研究了中心折反射摄像机下直线的像的几何性质，并将这些性质应用于中心折反射摄像机的标定。文献“calibrationofcentralcatadioptriccamerasusingadlt-likeapproach”,(puigl.,bastanlary.,sturmp.,etal.internationaljournalofcomputervision,vol.93,no.1,pp.101-114,2011)提出了一种基于三维控制点的标定方法，通过使用veronese映射对三维点和其图像点的坐标进行了扩展，在扩展坐标的基础上基于dlt(直接线性变换)——相似方法实现了中心折反射摄像机的标定，但是这类方法需要已知三维点的位置，并且容易从图像中提取其图像点。
球作为一种常见的几何体，其最重要的优点在于无自身遮挡，从任何一个方向看空间中一个球的封闭轮廓线总是一个圆，并且它的投影轮廓线可全部提取。由于球具有丰富的视觉几何特性，因此利用球进行摄像机标定已成为近年来的一个热点。文献“catadioptriccameracalibrationusinggeometricinvariants”，(yingx.,huz.,ieeetransactionsonpatternanalysisandmachineintelligence,vol.26,no.10,pp.1260-1271,2004)首次提出了利用球标定中心折反射摄像机。他们证明了球在中心折反射摄像机的单位球投影模型下的像为椭圆，并且在非退化情况下一个球的投影二次曲线提供两个不变量。他们提出了一种分步标定方法，该方法至少需要4个球的投影才能完成摄像机的标定。但是该文献提出的标定方法是非线性的，计算的复杂度较高，并且该标定方法只能标定抛物折反射摄像机的部分内参数。文献“identicalprojectivegeometricpropertiesofcentralcatadioptriclineimagesandsphereimageswithapplicationstocalibration”，(yingx.,zhah.,internationaljournalofcomputervision,vol.78,no.1,pp.89-105,2008)介绍了修正绝对二次曲线的像(miac)在中心折反射摄像机标定中的作用。他们通过研究球在中心折反射摄像机下的像与miac的几何与代数关系提出了两种线性标定算法。它们得出的结论对于对偶形式也是成立的。但是这篇文献中的理论和标定方法对于抛物折反射摄像机的情况是退化的。文献“acalibrationmethodforparacatadioptriccamerafromsphereimages”，(duanh.,wuy.,patternrecognitionletters,vol.33,no.6,pp.677-684,2012)基于圆环点理论提出了一种利用对拓球像标定抛物折反射摄像机的线性方法。但是这篇文献中关于圆环点的像的选取比较复杂。文献“cameracalibrationbasedonthecommonself-polartriangleofsphereimages”(huangh.,zhangh.,cheungy.m.,lecturenotesincomputerscience,9004:19-29,2015)利用空间球标定针孔摄像机。
技术实现要素：
本发明提供了一种制作简单，适用性广，稳定性好的利用靶标求解中心折反射摄像机内参数的方法，该靶标由空间中位置不同的三个球构成。在求解中心折反射摄像机内参数的过程中，需使用中心折反射摄像机拍摄靶标的一幅图像线性求解出中心折反射摄像机的5个内参数。本发明采用如下技术方案：
用中心折反射摄像机拍摄至少一幅含有三个球的图像。本发明是利用空间位置不同的三个球作为靶标用于求解中心折反射摄像机内参数的方法，其特征在于仅利用球元素。首先，从拍摄的球的图像中提取图像上的点拟合出来三个球所对应的球像cs1,cs2,cs3，再求出三个球在单位视球上投影所成的小圆s1,s2,s3的圆心所对应的像点om1,om2,om3。利用极点极线的关系求出om1关于cs1的极线ls1，om2关于cs2的极线ls2，om3关于cs3的极线ls3。进而求出cs1与ls1的两个交点(圆环点的像)，cs2与ls2的两个交点(圆环点的像)，cs3与ls3的两个交点(圆环点的像)。由上步得到的六个交点(圆环点的像)得到ωc，最后对ωc进行cholesky分解并求逆得摄像机内参数k。
1.拟合靶标投影方程
利用matlab程序中的edge函数提取靶标图像边缘点的像素坐标，并用最小二乘法拟合获得球像cs1,cs2,cs3的方程。
2.估计小圆s1,s2,s3的方程
空间中的球qsi(i＝1,2,3)在中心折反射摄像机的单位球模型下的投影分为两步，第一步，球qsi投影是以o为中心的单位视球上的小圆si(i＝1,2,3)，第二步，以单位视球表面上的一点oc为投影中心，这里oc看作一个摄像机的光心，将si(i＝1,2,3)分别投影为中心折反射图像平面上的二次曲线csi(i＝1,2,3)，其中称可见的二次曲线csi(i＝1,2,3)为球qsi(i＝1,2,3)的像。因为二次曲线系数可以由一个3×3的对称矩阵来表示，而且在中心折反射摄像机下球的像为二次曲线，所以可以用一个五维空间中的向量csi＝[a′b′c′d′e′f′]t(i＝1,2,3)来表示球qsi(i＝1,2,3)的像。令以oc为光心的摄像机的内参数矩阵为
根据斜锥tsi的顶点oc和准线csi来求出斜锥qsi的方程，对于斜锥tsi其顶点是虚拟摄像机光心oc，坐标为(0,0,l)，准线是球像csi。设点n(xi,yi,zi)(i＝1,2,3)为斜锥tsi(i＝1,2,3)上任一点，点
记
(i＝1,2,3)的方程。
由斜锥tsi(i＝1,2,3)的方程和单位视球方程联立可得小圆si(i＝1,2,3)的方程：
3.估计小圆s1,s2,s3圆心的像
要确定小圆s1,s2,s3圆心的像应先求出fy，首先从三个球像上分别提取3个不同点。从球像cs1中提取的三个点记为第一组，其齐次坐标的表达形式为(x1,y1,1),(x2,y2,1),(x3,y3,1)，cs2中提取的三个点记为第二组，分别为(x4,y4,1),(x5,y5,1),(x6,y6,1)；提取的三个点记为第三组，别为(x7,y7,1),(x8,y8,1),(x9,y9,1)。所以分别从三个球像上取三个点时就可以得到下面的方程组：
其中，δ表示为一个非零常数。
由上面的方程组可解得fy，为了得到精确的fy值。分别从三个球像上多取几组点由上面的方程组解出多个fy.舍去接近于零的值和负值，对剩下的数用ransac算法选出最优的fy。
设小圆si所在平面的方程为ax+by+cz+h＝0，a,b,c,h为平面方程系数。假设小圆si的圆心为s0i，坐标为(x0i,y0i,z0i)(i＝1,2,3)，则有，
x0i(xi-xj)+y0i(yi-yj)+z0i(zi-zj)＝δ，i＝1,2,3,j＝1,2,3,i≤j，其中，点(xi,yi,zi)，
(xj,yj,zj)，i＝1,2,3，j＝1,2,3为从小圆si上提取的六个已知点，
三角形ocs'0is0i相似于三角形ocopomi，其中op是图像的主点连接op和oc形成直线，过点s0i做op和oc形成直线的垂线交于点s'0i，如图2(i＝1为例)所示设点omi的齐次坐标为设点omi的齐次坐标为(xi,yi,1)(i＝1,2,3)，可得
4.求解绝对二次曲线的
在求解在omi之后可以根据极点极线关系求得lsi，从而求出ls1与cs1的两个交点(圆环点的像)，ls2和cs2的两个交点(圆环点的像)，ls3和cs3的两个交点(圆环点的像)，这六个点同时也在绝对二次曲线的像ωc上.在这里，选取空间中位置不同的三个球得到的三个球像便可以求得绝对二次曲线的像ωc上的六个点，然后由这六个交点得ωc。
5.求解摄像机内参数
圆环点的像记为mi和mj在绝对二次曲线的像ωc上，所以可得到下属关于ωc的约束方程：
本发明优点：
(1)该靶标制作简单，只需将三个不同位置的球固定在一个支架上。
(2)该靶标的图像边界点几乎可以全部提取，这样可以提高曲线拟合的精确度，从而提高标定精度。
附图说明
图1是空间球在中心折反射摄像机下的成像模型。
图2圆心以及圆心的像点之间所构成相似三角形。
图3是极线与球像的关系。
具体实施方式
本发明提供了一种利用靶标求解中心折反射摄像机内参数的方法，靶标是由空间中的三个位置不同的球构成，用此靶标完成中心折反射摄像机内参数的求解需要经过以下步骤：从中心折反射图像中提取靶标图像边缘点的像素坐标，并用最小二乘法拟合球像的方程cs1,cs2,cs3。三个空间球在镜面上的投影为三个小圆s1,s2,s3，求出三个小圆s1,s2,s3的圆心进而求出圆心所对应的像点omi(i＝1,2,3)，再求出圆心的像点omi(i＝1,2,3)关于球像的极线lsi(i＝1,2,3)，最后求出极线lsi(i＝1,2,3)与其对应的球像csi(i＝1,2,3)的交点(圆环点的像)，三条极线共对应六个位于ωc上的圆环点的像，由得到的六个圆环点的像得到ωc，然后对其进行cholesky分解，再求逆得到摄像机的5个内参数。利用本发明中的方法对用于实验的中心折反射摄像机进行标定，具体步骤如下：
1.拟合靶标投影方程
利用matlab程序中的edge函数提取靶标图像边缘点的像素坐标，并用最小二乘法拟合获得球像cs1,cs2,cs3的方程。
2.估计小圆s1,s2,s3的方程
空间中的球qsi(i＝1,2,3)在中心折反射摄像机的单位球模型下如图1(i＝1为例)投影分为两步，第一步，球qsi投影是以世界坐标系o-xwywzw的原点o为中心的单位视球上的小圆si(i＝1,2,3)，o即为单位视球的球心。第二步，以单位视球表面上的一点为摄像机坐标系oc-xcyczc，其中oc为投影中心，这里oc看作一个摄像机的光心，将si(i＝1,2,3)分别投影为中心折反射图像平面上的二次曲线csi(i＝1,2,3)，其中称可见的二次曲线csi(i＝1,2,3)为球qsi(i＝1,2,3)的像，成像面与光轴ooc垂直，zw轴和zc轴与光轴ooc重合，xw,xc轴和yw,yc轴与成像平面上x轴方向和y轴方向平行。因为二次曲线系数可以由一个3×3的对称矩阵来表示，而且在中心折反射摄像机下球的像为二次曲线，球qsi(i＝1,2,3)的像可以用一个五维空间中的向量来表示
csi＝[a′b′c′d′e′f′]t，(i＝1,2,3)。(1)
令以oc为光心的摄像机的内参数矩阵为
如图1，根据斜锥tsi的顶点oc坐标和准线csi的方程求出斜锥qsi的方程，对于斜锥tsi其顶点是虚拟摄像机光心oc(0,0,l)，准线是球像csi，设点n(xi,yi,zi)(i＝1,2,3)为斜锥tsi(i＝1,2,3)上任一点，点
其中u为一非零参数,
其中
即为所求斜锥tsi的方程。由斜锥tsi的方程和单位视球方程联立可得小圆si的方程：
其中
3.估计小圆圆心的像
要确定小圆s1,s2,s3圆心的像应先求出fy，首先从三个球像上分别提取3个不同点。从球像cs1中提取的三个点记为第一组，其齐次坐标的表达形式为(x1,y1,1),(x2,y2,1),(x3,y3,1)，cs2中提取的三个点记为第二组，分别为(x4,y4,1),(x5,y5,1),(x6,y6,1)，cs3中提取的三个点记为第三组，分别为(x7,y7,1),(x8,y8,1),(x9,y9,1)。
设mi表示球面上点xi的像点，并且xi和mi的投影关系满足
其中的e＝(0,0,1)t。令
上式又可以表示为
所以分别从三个球像上取三个点时就可以得到下面方程组：
其中，δ表示为一个非零常数，
设小圆si所在平面的方程为ax+by+cz+h＝0，a,b,c,h为平面的方程系数。假设小圆si的圆心为s0i，在成像面的投影为omi，图1中i＝1为例，坐标为(x0i,y0i,z0i)(i＝1,2,3)，则有，
x0i(xi-xj)+y0i(yi-yj)+z0i(zi-zj)＝δ,i＝1,2,3,j＝1,2,3,i≤j，(7)
其中，点(xi,yi,zi),(xj,yj,zj),i＝1,2,3,j＝1,2,3为从小圆si(i＝1,2,3)上提取的六个已知点，
三角形ocs'0is0i相似于三角形ocopomi，其中op是图像的主点。连接op和oc形成直线，过点s0i做op和oc形成直线的垂线交于点s'0i，如图2(i＝1为例)所示。设点omi的齐次坐标为如图2(i＝1为例)所示设点omi的齐次坐标为(xi,yi,1)，可得
整理可得omi齐次坐标为
4.求解绝对二次曲线的像
已知csi与ωc交于四个点，其中的两个交点的连线为lsi(i＝1,2,3)，在求得lsi(i＝1,2,3)之后，便可以求得lsi与csi的两个交点(即圆环点)，如图3(i＝1为例)这两个点同时也在ωc上.极点omi和极线lsi的关系式为(9)式:
lsi＝csiomi,(i＝1,2,3),(9)
由(9)式可以求出极线lsi的方程记为(10)式
整理可得(11)式:
上式即lsi的方程。由ls1的方程和cs1的方程，ls2的方程和cs2的方程，ls3的方程和cs3的方程即可求它们所对应的两个交点，共可求出6个交点。由这6个交点进行得到绝对二次曲线的像ωc。
5.求解摄像机内参数
圆环点的像在绝对二次曲线的像ωc上，所以可得到下属关于ωc的约束方程：
由于mi,mj是两共轭复点，所以在上面的方程中只能提供下面两个关于绝对二次曲线的像ωc的两个实线性约束方程：
其中re，im分别表示负数的实部和虚部，由上面的方程可求出绝对二次曲线的像ωc，然后利用cholesky分解法对ωc进行分解可唯一确定k-1,再对k-1求逆就可得到射相机的内参数矩阵k。
实施例
本发明提出了一种利用空间的三个球作为标定物线性求解中心折反射摄像机内参数的方法，本发明采用的实验模板结构示意图如图1所示。下面用一实施方案作出更为详细的描述。
中心折反射摄像机标定采用的模板为空间中位置不同的三个球如图1球即为qsi(i＝1,2,3)，利用本发明的方法对本实验的中心折反射摄像机进行标定，具体步骤如下：1.拟合靶标曲线方程
本发明采用的图像大小为1700×1600。用中心折反射摄像机拍摄靶标的1幅实验图像，读入图像，利用matlab中的edge函数提取图像靶标图像边缘点的像素坐标，并用最小二乘法拟合获得三个球像的方程。三个空间球的球像的方程的系数矩阵分别为csi(i＝1,2,3)，结果如下：
2.求小圆si(i＝1,2,3)的圆心的像omi的坐标
先由(6)式求出fy，再把(14)，(15)，(16)分别代入(2)式再结和(5)求得小圆si的方程，再由(7)式得到小圆si(i＝1,2,3)的圆心最后由(8)式得到圆心的像omi的坐标分别为：om1齐次坐标为(1263.9140.21)t，om2的齐次坐标为(1309.9303.41)t，om3的齐次坐标为(1174.7232.51)t。
3.圆心的像omi所对应的极线
由上面所求出的omi的齐次坐标分别代入(11)式整理可以得到
极线ls1的齐次线坐标的矩阵形式为
ls1＝[-0.00530.00681.0000]t，(17)
极线ls2的齐次线坐标的矩阵形式为
ls2＝[-0.00330.00261.0000]t，(18)
极线ls3的齐次线坐标的矩阵形式为
ls3＝[-0.01310.01501.0000]t。(19)
4.求解交点，即圆环点的像
分别连立(14)和(17)；(15)和(18)；(16)和(19)得到对应的两个交点：
ls1和cs1的交点为复点，其交点的坐标矩阵为：
m1i＝[724.0418796372854+264.1304204451712i124.2209398186203+289.3652102225750i]t，(20)
m1j＝[724.0418796372854-264.1304204451712i124.2209398186203-289.3652102225750i]t；(21)
ls2和cs2的交点为复点，其交点的坐标矩阵为：
m2i＝[698.6400000000319+149.3866666666647i226.6666666666646+266.666666666517i]t，(22)
m2j＝[698.6400000000319-149.3866666666647i226.6666666666646-266.666666666517i]t；(23)
ls3和cs3的交点为复点，其交点的坐标矩阵为：
m3i＝[634.6205052005533+169.7893016344956i208.4398216939404+209.5096582466692i]t，(24)
m3i＝[634.6205052005533-169.7893016344956i208.4398216939404-209.5096582466692i]t。(25)
5.求解摄像机内参数
对上面的(20)，(22)，(22)，(23)，(24)，(25)带入(13)可以得到关于ωc方程，使用svd分解求解该线性方程组得到ωc的系数矩阵形式如下：
对(26)进行cholesky分解并求逆得到：