CN107958468A - 利用空间位置不同的三个球标定中心折反射摄像机的方法 - Google Patents
利用空间位置不同的三个球标定中心折反射摄像机的方法 Download PDFInfo
- Publication number
- CN107958468A CN107958468A CN201711343962.7A CN201711343962A CN107958468A CN 107958468 A CN107958468 A CN 107958468A CN 201711343962 A CN201711343962 A CN 201711343962A CN 107958468 A CN107958468 A CN 107958468A
- Authority
- CN
- China
- Prior art keywords
- msub
- msup
- ball
- mrow
- picture
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
- G06T7/85—Stereo camera calibration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20068—Projection on vertical or horizontal image axis
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
- Image Analysis (AREA)
Abstract
本发明是利用空间位置不同的三个球标定中心折反射摄像机的方法,其特征在于仅利用球元素。首先,从拍摄的球的图像中提取图像上的点拟合出来三个球所对应的球像,再求出三个球在单位视球上投影所成的小圆的圆心所对应的像点。利用极点极线的关系求出圆心关于球像的极线。进而求出球像与极线的两个交点,即圆环点的像。由上步得到的六个圆环点的像得到绝对二次曲线,最后对绝对二次曲线进行分解并求逆得摄像机内参数。
Description
技术领域
本发明属于计算机视觉领域,涉及一种利用空间中三个球求解中心折反射摄像机内参数的方法。
背景技术
计算机视觉是指通过计算机、图像采集等设备来模拟人类的视觉对目标进行识别、检测、跟踪等。计算机视觉的中心任务就是对图像进行获取,处理,分析。而它的最终目标是使计算机具有通过二维图像认知三维世界的能力。在视觉测量中,计算机从摄像机获得二维图像信息来计算三维空间中物体的几何信息、并由此重建和识别物体,而空间物体表面点的三维几何位置与在图像中的对应点之间的相互关系是由摄像机成像的几何模型决定的,这些的几何模型就是摄像机的内参数,而确定这些参数的方法就叫做摄像机的标定。摄像机标定可分为传统标定、自标定和基于几何实体的标定。无论哪种标定方法,都旨在建立二维图像与摄像机内参数之间的约束关系,特别是线性约束关系,这是目前摄像机标定所追求的目标,也是目前计算机视觉领域研究的热点之一。
文献“Geometric properties of central catadioptric line images andtheir application in calibration”,(Barreto J.P.,Araujo H.,IEEE Transactionson Pattern Analysis and Machine Intelligence,vol.27,no.8,pp.1327-1333,2005)研究了中心折反射摄像机下直线的像的几何性质,并将这些性质应用于中心折反射摄像机的标定。文献“Calibration of central catadioptric cameras using a DLT-likeapproach”,(Puig L.,Bastanlar Y.,Sturm P.,et al.International Journal ofComputer Vision,vol.93,no.1,pp.101-114,2011)提出了一种基于三维控制点的标定方法,通过使用Veronese映射对三维点和其图像点的坐标进行了扩展,在扩展坐标的基础上基于DLT(直接线性变换)——相似方法实现了中心折反射摄像机的标定,但是这类方法需要已知三维点的位置,并且容易从图像中提取其图像点。
球作为一种常见的几何体,其最重要的优点在于无自身遮挡,从任何一个方向看空间中一个球的封闭轮廓线总是一个圆,并且它的投影轮廓线可全部提取。由于球具有丰富的视觉几何特性,因此利用球进行摄像机标定已成为近年来的一个热点。文献“Catadioptric camera calibration using geometric invariants”,(Ying X.,Hu Z.,IEEE Transactions on Pattern Analysis and Machine Intelligence,vol.26,no.10,pp. 1260-1271,2004)首次提出了利用球标定中心折反射摄像机。他们证明了球在中心折反射摄像机的单位球投影模型下的像为椭圆,并且在非退化情况下一个球的投影二次曲线提供两个不变量。他们提出了一种分步标定方法,该方法至少需要4个球的投影才能完成摄像机的标定。但是该文献提出的标定方法是非线性的,计算的复杂度较高,并且该标定方法只能标定抛物折反射摄像机的部分内参数。文献“Identical projective geometricproperties of central catadioptric line images and sphere images withapplications to calibration”,(Ying X.,Zha H.,International Journal ofComputer Vision, vol.78,no.1,pp.89-105,2008)介绍了修正绝对二次曲线的像(MIAC)在中心折反射摄像机标定中的作用。他们通过研究球在中心折反射摄像机下的像与MIAC的几何与代数关系提出了两种线性标定算法。它们得出的结论对于对偶形式也是成立的。但是这篇文献中的理论和标定方法对于抛物折反射摄像机的情况是退化的。文献“Acalibration method for paracatadioptric camera from sphere images”,(Duan H.,Wu Y.,Pattern Recognition Letters,vol.33,no.6,pp.677-684,2012)基于圆环点理论提出了一种利用对拓球像标定抛物折反射摄像机的线性方法。但是这篇文献中关于圆环点的像的选取比较复杂。文献“Camera calibration based on the common self-polartriangle of sphere images”(Huang H.,Zhang H.,Cheung Y.M.,Lecture Notes inComputer Science,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为光心的摄像机的内参数矩阵为其中fx,fy为摄像机在x,y轴方向上的焦距, (u0,v0)为像平面上的主点坐标,s为摄像机的倾斜因子。
根据斜锥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)为准线CSi上任一点;点N(xi,yi,zi) (i=1,2,3)N≠N0在锥面上的充分必要条件是:N在一条母线上,即可得斜锥TSi (i=1,2,3)的方程为,其中u为一非零参数
将上式整理可得:
其中
记即为所求斜锥TSi
(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上提取的六个已知点,
由上式求出S0i(i=1,2,3)的坐标。
三角形OCS'0iS0i相似于三角形OCOPOmi,其中OP是图像的主点连接OP和OC形成直线,过点S0i做OP和OC形成直线的垂线交于点S'0i,如图2(i=1为例)所示设点Omi的齐次坐标为设点Omi的齐次坐标为(xi,yi,1)(i=1,2,3),可得整理可得Omi (i=1,2,3)齐次坐标为
4.求解绝对二次曲线的
在求解在Omi之后可以根据极点极线关系求得lSi,从而求出lS1与CS1的两个交点(圆环点的像),lS2和CS2的两个交点(圆环点的像),lS3和CS3的两个交点(圆环点的像),这六个点同时也在绝对二次曲线的像ωC上.在这里,选取空间中位置不同的三个球得到的三个球像便可以求得绝对二次曲线的像ωC上的六个点,然后由这六个交点得ωC。
5.求解摄像机内参数
圆环点的像记为mI和mJ在绝对二次曲线的像ωC上,所以可得到下属关于ωC的约束方程:由于mI,mJ是两共轭复点,所以在上面的方程中只能提供下面两个关于绝对二次曲线的像ωC的两个实线性约束方程其中Re,Im分别表示复数的实部和虚部,由上面的方程可求出绝对二次曲线的像ωC,然后利用Cholesky分解法对ωC进行分解可唯一确定K-1,再对K-1求逆就可得到射相机的内参数矩阵K。
本发明优点:
(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为光心的摄像机的内参数矩阵为其中fx,fy为摄像机成像平面上x轴方向和y轴方向的焦距,(u0,v0)为像平面上的主点坐标,s为摄像机的倾斜因子。
如图1,根据斜锥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)为准线CSi上任一点(点N(xi,yi,zi) (i=1,2,3)N≠N0在锥面上的充分必要条件是:N在一条母线上)即可得斜锥TSi (i=1,2,3)的方程为
其中u为一非零参数,图1中i=1为例。将(2)式整理可得:
其中记:
即为所求斜锥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的投影关系满足i=(1,2,3),这里
其中i=(1,2,3)。并且易知Xi共线且在球面上的同一个大圆上,也就是说,是共面因此有:det[X1-O,X2-O,X3-O]=0,由上面的投影关系。上式可表示为
det表示矩阵行列式的值。令从而表示将图像平面的原点转化为主点p,即有:
整理有
其中的e=(0,0,1)T。令即
令
在本文中所用的中心折反射摄像机为抛物折反射摄像机所以l=1即τ=0时
可化简为
上式又可以表示为
所以分别从三个球像上取三个点时就可以得到下面方程组:
其中,δ表示为一个非零常数, 由上面的方程组可解得fy。为了得到精确的fy值,多取几组舍去接近于零的值和负值,对剩下的数用RANSAC算法选出最优的fy。
设小圆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.9 140.2 1)T,Om2的齐次坐标为(1309.9 303.4 1)T,Om3的齐次坐标为(1174.7 232.5 1)T。
3.圆心的像Omi所对应的极线
由上面所求出的Omi的齐次坐标分别代入(11)式整理可以得到
极线lS1的齐次线坐标的矩阵形式为
lS1=[-0.0053 0.0068 1.0000]T, (17)
极线lS2的齐次线坐标的矩阵形式为
lS2=[-0.0033 0.0026 1.0000]T, (18)
极线lS3的齐次线坐标的矩阵形式为
lS3=[-0.0131 0.0150 1.0000]T。 (19)
4.求解交点,即圆环点的像
分别连立(14)和(17);(15)和(18);(16)和(19)得到对应的两个交点:
lS1和CS1的交点为复点,其交点的坐标矩阵为:
m1I=[724.0418796372854+264.1304204451712i 124.2209398186203+289.3652102225750i]T, (20)
m1J=[724.0418796372854-264.1304204451712i 124.2209398186203-289.3652102225750i]T; (21)
lS2和CS2的交点为复点,其交点的坐标矩阵为:
m2I=[698.6400000000319+149.3866666666647i 226.6666666666646+266.666666666517i]T, (22)
m2J=[698.6400000000319-149.3866666666647i 226.6666666666646-266.666666666517i]T; (23)
lS3和CS3的交点为复点,其交点的坐标矩阵为:
m3I=[634.6205052005533+169.7893016344956i 208.4398216939404+209.5096582466692i]T, (24)
m3I=[634.6205052005533-169.7893016344956i 208.4398216939404-209.5096582466692i]T。 (25)
5.求解摄像机内参数
对上面的(20),(22),(22),(23),(24),(25)带入(13)可以得到关于ωC方程,使用SVD分解求解该线性方程组得到ωC的系数矩阵形式如下:
对(26)进行Cholesky分解并求逆得到:其中摄像机的5 个内参数分别为:fx=567.902,fy=507.055,s=0.1014,u0=400.000,v0=359.999。
Claims (1)
1.利用空间位置不同的三个球标定中心折反射摄像机的方法,其特征在于仅利用球元素作为靶标;所述方法的具体步骤包括:首先,从拍摄的球的图像中提取图像上的点拟合出来三个球所对应的球像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)估计小圆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为光心的摄像机的内参数矩阵为其中fx,fy为摄像机在x,y轴方向上的焦距,(u0,v0)为像平面上的主点坐标,s为摄像机的倾斜因子;
根据斜锥TSi的顶点OC和准线CSi来求出斜锥QSi的方程,对于斜锥TSi其顶点是虚拟摄像机光心OC,坐标为(0,0,l),准线是球像CSi;设点N(xi,yi,zi)(i=1,2,3)为斜锥TSi(i=1,2,3)上任一点,点为准线CSi上任一点;点N(xi,yi,zi)(i=1,2,3)N≠N0在锥面上的充分必要条件是:N在一条母线上,即得斜锥TSi(i=1,2,3)的方程为,其中u为一非零参数,
将上式整理得:
其中
记即为所求斜锥TSi(i=1,2,3)的方程;由斜锥TSi,其中i=1,2,3,的方程和单位视球方程联立得小圆Si,其中i=1,2,3,的方程:其中
(2)估计小圆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);所以分别从三个球像上取三个点时就得到下面的方程组:
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>3</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>3</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>3</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>3</mn>
</msub>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>)</mo>
</mrow>
<msub>
<mi>&delta;f</mi>
<mi>y</mi>
</msub>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>3</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<msup>
<mi>&delta;</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>6</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>5</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>4</mn>
</msub>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>4</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>4</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>5</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>5</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>6</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>6</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>6</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>4</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>5</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>5</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>4</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>6</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>6</mn>
</msub>
<msub>
<mi>x</mi>
<mn>4</mn>
</msub>
<msub>
<mi>y</mi>
<mn>4</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>5</mn>
</msub>
<msub>
<mi>x</mi>
<mn>5</mn>
</msub>
<msub>
<mi>y</mi>
<mn>5</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>4</mn>
</msub>
<msub>
<mi>x</mi>
<mn>6</mn>
</msub>
<msub>
<mi>y</mi>
<mn>6</mn>
</msub>
<mo>)</mo>
</mrow>
<msub>
<mi>&delta;f</mi>
<mi>y</mi>
</msub>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>6</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>4</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>5</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>5</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>4</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>6</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<msup>
<mi>&delta;</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>9</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>8</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>7</mn>
</msub>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>9</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>7</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>8</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>8</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>7</mn>
</msub>
<msup>
<msub>
<mi>x</mi>
<mn>9</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>9</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>7</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>8</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>8</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>7</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>9</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
<msup>
<msub>
<mi>f</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>9</mn>
</msub>
<msub>
<mi>x</mi>
<mn>7</mn>
</msub>
<msub>
<mi>y</mi>
<mn>7</mn>
</msub>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>8</mn>
</msub>
<msub>
<mi>x</mi>
<mn>8</mn>
</msub>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>7</mn>
</msub>
<msub>
<mi>x</mi>
<mn>9</mn>
</msub>
<msub>
<mi>y</mi>
<mn>9</mn>
</msub>
<mo>)</mo>
</mrow>
<msub>
<mi>&delta;f</mi>
<mi>y</mi>
</msub>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>U</mi>
<mn>9</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>7</mn>
</msub>
<mn>2</mn>
</msup>
<mo>-</mo>
<msub>
<mi>U</mi>
<mn>8</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>8</mn>
</msub>
<mn>2</mn>
</msup>
<mo>+</mo>
<msub>
<mi>U</mi>
<mn>7</mn>
</msub>
<msup>
<msub>
<mi>y</mi>
<mn>9</mn>
</msub>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<msup>
<mi>&delta;</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
其中,δ表示为一个非零常数,
由上面的方程组解得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上提取的六个已知点,
由上式求出S0i(i=1,2,3)的坐标;
三角形OCS'0iS0i相似于三角形OCOPOmi,其中OP是图像的主点连接OP和OC形成直线,过点S0i做OP和OC形成直线的垂线交于点S'0i;设点Omi的齐次坐标为设点Omi的齐次坐标为(xi,yi,1),其中i=1,2,3,得整理得Omi,其中i=1,2,3,齐次坐标为
(3)求解绝对二次曲线的
在求解在Omi之后根据极点极线关系求得lSi,从而求出lS1与CS1的两个交点,lS2和CS2的两个交点,lS3和CS3的两个交点,这六个点同时也在绝对二次曲线的像ωC上,即圆环点的像上;在这里,选取空间中位置不同的三个球得到的三个球像便求得绝对二次曲线的像ωC上的六个点,然后由这六个交点得ωC。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711343962.7A CN107958468B (zh) | 2017-12-15 | 2017-12-15 | 利用空间位置不同的三个球标定中心折反射摄像机的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711343962.7A CN107958468B (zh) | 2017-12-15 | 2017-12-15 | 利用空间位置不同的三个球标定中心折反射摄像机的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107958468A true CN107958468A (zh) | 2018-04-24 |
CN107958468B CN107958468B (zh) | 2021-06-08 |
Family
ID=61958916
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711343962.7A Expired - Fee Related CN107958468B (zh) | 2017-12-15 | 2017-12-15 | 利用空间位置不同的三个球标定中心折反射摄像机的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107958468B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109272558A (zh) * | 2018-11-28 | 2019-01-25 | 云南大学 | 分离圆的公共自极三角形及圆环点标定针孔摄像机的方法 |
CN110428472A (zh) * | 2019-07-19 | 2019-11-08 | 中北大学 | 三角形基元立体球靶标及工业相机标定方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030004694A1 (en) * | 2001-05-29 | 2003-01-02 | Daniel G. Aliaga | Camera model and calibration procedure for omnidirectional paraboloidal catadioptric cameras |
CN101414380A (zh) * | 2008-11-07 | 2009-04-22 | 浙江大学 | 全景相机的单图标定方法 |
CN102982551A (zh) * | 2012-12-14 | 2013-03-20 | 云南大学 | 空间三条不平行直线线性求解抛物折反射摄像机内参数 |
US8432435B2 (en) * | 2011-08-10 | 2013-04-30 | Seiko Epson Corporation | Ray image modeling for fast catadioptric light field rendering |
CN105303570A (zh) * | 2015-10-22 | 2016-02-03 | 云南大学 | 使用双球的相交像及圆环点的像标定拋物折反射摄像机 |
CN106558081A (zh) * | 2016-11-28 | 2017-04-05 | 云南大学 | 标定非中心轴对称系统的圆锥折反射摄像机的方法 |
CN106558082A (zh) * | 2016-11-28 | 2017-04-05 | 云南大学 | 利用Veronese映射棋盘格的投影矩阵标定中心折反射摄像机 |
-
2017
- 2017-12-15 CN CN201711343962.7A patent/CN107958468B/zh not_active Expired - Fee Related
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030004694A1 (en) * | 2001-05-29 | 2003-01-02 | Daniel G. Aliaga | Camera model and calibration procedure for omnidirectional paraboloidal catadioptric cameras |
CN101414380A (zh) * | 2008-11-07 | 2009-04-22 | 浙江大学 | 全景相机的单图标定方法 |
US8432435B2 (en) * | 2011-08-10 | 2013-04-30 | Seiko Epson Corporation | Ray image modeling for fast catadioptric light field rendering |
CN102982551A (zh) * | 2012-12-14 | 2013-03-20 | 云南大学 | 空间三条不平行直线线性求解抛物折反射摄像机内参数 |
CN105303570A (zh) * | 2015-10-22 | 2016-02-03 | 云南大学 | 使用双球的相交像及圆环点的像标定拋物折反射摄像机 |
CN106558081A (zh) * | 2016-11-28 | 2017-04-05 | 云南大学 | 标定非中心轴对称系统的圆锥折反射摄像机的方法 |
CN106558082A (zh) * | 2016-11-28 | 2017-04-05 | 云南大学 | 利用Veronese映射棋盘格的投影矩阵标定中心折反射摄像机 |
Non-Patent Citations (3)
Title |
---|
DUAN H , WU Y: "Paracatadioptric camera calibration using sphere images", 《IEEE INTERNATIONAL CONFERENCE ON IMAGE PROCESSING》 * |
ZHAO Y , WANG Y: "Intrinsic parameter determination of a paracatadioptric camera by the intersection of two sphere projections", 《JOURNAL OF THE OPTICAL SOCIETY OF AMERICA A OPTICS IMAGE SCIENCE & VISION》 * |
顾万家: "基于完全四点形利用空间球标定拋物折反射摄像机", 《中国优秀硕士学位论文全文数据库 信息科技辑》 * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109272558A (zh) * | 2018-11-28 | 2019-01-25 | 云南大学 | 分离圆的公共自极三角形及圆环点标定针孔摄像机的方法 |
CN109272558B (zh) * | 2018-11-28 | 2021-07-09 | 云南大学 | 分离圆的公共自极三角形及圆环点标定针孔摄像机的方法 |
CN110428472A (zh) * | 2019-07-19 | 2019-11-08 | 中北大学 | 三角形基元立体球靶标及工业相机标定方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107958468B (zh) | 2021-06-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107886546B (zh) | 利用球像及公共自极三角形标定抛物折反射摄像机的方法 | |
CN101216296A (zh) | 双目视觉转轴标定方法 | |
CN106327504A (zh) | 利用单个球及圆环点的像标定拋物折反射摄像机的方法 | |
CN107644445A (zh) | 利用单个球及圆切线的性质标定拋物折反射摄像机的方法 | |
CN104517291A (zh) | 基于目标同轴圆特征的位姿测量方法 | |
CN106651956A (zh) | 利用单个球及平行圆性质标定拋物折反射摄像机的方法 | |
CN107657645A (zh) | 利用直线及圆的共轭直径的性质标定拋物折反射摄像机 | |
WO2019056782A1 (zh) | 一种基于球体投影公切线的多相机标定及参数优化方法 | |
CN104200476B (zh) | 利用双平面镜装置中的圆周运动求解摄像机内参数的方法 | |
CN105279758A (zh) | 使用双球的相切像与圆环点的像标定拋物折反射摄像机 | |
CN108921889A (zh) | 一种基于扩增现实应用的室内三维定位方法 | |
CN109360248A (zh) | 利用单个球及共轭直径的性质标定拋物折反射摄像机 | |
CN109325983A (zh) | 利用无穷远点关于圆极线性质标定拋物折反射摄像机 | |
CN105321181A (zh) | 使用双球的相离像与圆环点的像标定拋物折反射摄像机 | |
CN107958468A (zh) | 利用空间位置不同的三个球标定中心折反射摄像机的方法 | |
CN105303570A (zh) | 使用双球的相交像及圆环点的像标定拋物折反射摄像机 | |
Bergamasco et al. | A robust multi-camera 3D ellipse fitting for contactless measurements | |
Chen et al. | A novel mirrored binocular vision sensor based on spherical catadioptric mirrors | |
Yang et al. | Camera calibration using projective invariants of sphere images | |
CN110148183A (zh) | 利用球及极点极线标定摄像机的方法、存储介质和系统 | |
CN109360247A (zh) | 单个球的公共自极三角形及正交消失点标定抛物摄像机 | |
CN109523598A (zh) | 单个球的公共自极三角形及圆环点标定抛物摄像机的方法 | |
CN109712195A (zh) | 利用球像的公共自极三角形进行单应性估计的方法 | |
CN109325982A (zh) | 利用单个球及平行圆切线性质标定拋物折反射摄像机 | |
CN109035342A (zh) | 利用一条直线及圆环点极线标定拋物折反射摄像机的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210608 Termination date: 20211215 |