CN107958468B - 利用空间位置不同的三个球标定中心折反射摄像机的方法 - Google Patents

利用空间位置不同的三个球标定中心折反射摄像机的方法 Download PDF

Info

Publication number
CN107958468B
CN107958468B CN201711343962.7A CN201711343962A CN107958468B CN 107958468 B CN107958468 B CN 107958468B CN 201711343962 A CN201711343962 A CN 201711343962A CN 107958468 B CN107958468 B CN 107958468B
Authority
CN
China
Prior art keywords
image
points
point
camera
equation
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.)
Expired - Fee Related
Application number
CN201711343962.7A
Other languages
English (en)
Other versions
CN107958468A (zh
Inventor
赵越
张敬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yunnan University YNU
Original Assignee
Yunnan University YNU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Yunnan University YNU filed Critical Yunnan University YNU
Priority to CN201711343962.7A priority Critical patent/CN107958468B/zh
Publication of CN107958468A publication Critical patent/CN107958468A/zh
Application granted granted Critical
Publication of CN107958468B publication Critical patent/CN107958468B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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)
  • Image Analysis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明是利用空间位置不同的三个球标定中心折反射摄像机的方法，其特征在于仅利用球元素。首先，从拍摄的球的图像中提取图像上的点拟合出来三个球所对应的球像，再求出三个球在单位视球上投影所成的小圆的圆心所对应的像点。利用极点极线的关系求出圆心关于球像的极线。进而求出球像与极线的两个交点，即圆环点的像。由上步得到的六个圆环点的像得到绝对二次曲线，最后对绝对二次曲线进行
Figure DEST_PATH_IMAGE001
分解并求逆得摄像机内参数。

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为光心的摄像机的内参数矩阵为
Figure BDA0001508929080000031
其中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)上任一点，点
Figure BDA0001508929080000032
(i＝1,2,3)为准线CSi上任一点；点N(xi,yi,zi) (i＝1,2,3)N≠N0在锥面上的充分必要条件是：N在一条母线上，即可得斜锥TSi (i＝1,2,3)的方程为，
Figure BDA0001508929080000033
其中u为一非零参数
Figure BDA0001508929080000034
将上式整理可得：
Figure BDA0001508929080000035
其中
Figure BDA0001508929080000036
记
Figure BDA0001508929080000041
即为所求斜锥TSi
(i＝1,2,3)的方程。
由斜锥TSi(i＝1,2,3)的方程和单位视球方程联立可得小圆Si(i＝1,2,3)的方程：
Figure BDA0001508929080000042
其中
Figure BDA0001508929080000043
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)。所以分别从三个球像上取三个点时就可以得到下面的方程组：
Figure BDA0001508929080000044
其中，δ表示为一个非零常数。
Figure BDA0001508929080000045
Figure BDA0001508929080000046
由上面的方程组可解得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上提取的六个已知点，
Figure BDA0001508929080000051
由上式求出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)，可得
Figure BDA0001508929080000052
整理可得Omi (i＝1,2,3)齐次坐标为
Figure BDA0001508929080000053
4.求解绝对二次曲线的
在求解在Omi之后可以根据极点极线关系求得lSi，从而求出lS1与CS1的两个交点(圆环点的像)，lS2和CS2的两个交点(圆环点的像)，lS3和CS3的两个交点(圆环点的像)，这六个点同时也在绝对二次曲线的像ωC上.在这里，选取空间中位置不同的三个球得到的三个球像便可以求得绝对二次曲线的像ωC上的六个点，然后由这六个交点得ωC。
5.求解摄像机内参数
圆环点的像记为mI和mJ在绝对二次曲线的像ωC上，所以可得到下属关于ωC的约束方程：
Figure BDA0001508929080000054
由于mI,mJ是两共轭复点，所以在上面的方程中只能提供下面两个关于绝对二次曲线的像ωC的两个实线性约束方程
Figure BDA0001508929080000055
其中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为光心的摄像机的内参数矩阵为
Figure BDA0001508929080000061
其中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)上任一点，点
Figure BDA0001508929080000071
(i＝1,2,3)为准线CSi上任一点(点N(xi,yi,zi) (i＝1,2,3)N≠N0在锥面上的充分必要条件是：N在一条母线上)即可得斜锥TSi (i＝1,2,3)的方程为
Figure BDA0001508929080000072
其中u为一非零参数,
Figure BDA0001508929080000073
图1中i＝1为例。将(2)式整理可得：
Figure BDA0001508929080000074
其中
Figure BDA0001508929080000075
记：
Figure BDA0001508929080000076
即为所求斜锥TSi的方程。由斜锥TSi的方程和单位视球方程联立可得小圆Si的方程：
Figure BDA0001508929080000077
其中
Figure BDA0001508929080000078
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的投影关系满足
Figure BDA0001508929080000081
i＝(1,2,3)，这里
Figure BDA0001508929080000082
Figure BDA0001508929080000083
Figure BDA0001508929080000084
其中
Figure BDA0001508929080000085
i＝(1,2,3)。并且易知Xi共线且在球面上的同一个大圆上，也就是说,
Figure BDA0001508929080000086
是共面因此有：det[X1-O,X2-O,X3-O]＝0，由上面的投影关系。上式可表示为
Figure BDA0001508929080000087
det表示矩阵行列式的值。令
Figure BDA0001508929080000088
从而
Figure BDA0001508929080000089
表示将图像平面的原点转化为主点p，即有：
Figure BDA00015089290800000810
整理有
Figure BDA00015089290800000811
其中的e＝(0,0,1)T。令
Figure BDA00015089290800000812
即
Figure BDA00015089290800000813
令
Figure BDA00015089290800000814
在本文中所用的中心折反射摄像机为抛物折反射摄像机所以l＝1即τ＝0时
Figure BDA00015089290800000815
可化简为
Figure BDA00015089290800000816
上式又可以表示为
Figure BDA00015089290800000817
Figure BDA0001508929080000091
所以分别从三个球像上取三个点时就可以得到下面方程组：
Figure BDA0001508929080000092
其中，δ表示为一个非零常数，
Figure BDA0001508929080000093
Figure BDA0001508929080000094
由上面的方程组可解得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)上提取的六个已知点，
Figure BDA0001508929080000095
三角形OCS'0iS0i相似于三角形OCOPOmi，其中OP是图像的主点。连接OP和OC形成直线，过点S0i做OP和OC形成直线的垂线交于点S'0i，如图2(i＝1为例)所示。设点Omi的齐次坐标为如图2(i＝1为例)所示设点Omi的齐次坐标为(xi,yi,1)，可得
Figure BDA0001508929080000096
整理可得Omi齐次坐标为
Figure BDA0001508929080000097
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)式
Figure BDA0001508929080000101
整理可得(11)式:
Figure BDA0001508929080000102
上式即lSi的方程。由lS1的方程和CS1的方程，lS2的方程和CS2的方程，lS3的方程和CS3的方程即可求它们所对应的两个交点，共可求出6个交点。由这6个交点进行得到绝对二次曲线的像ωC。
5.求解摄像机内参数
圆环点的像在绝对二次曲线的像ωC上，所以可得到下属关于ωC的约束方程：
Figure BDA0001508929080000103
由于mI,mJ是两共轭复点，所以在上面的方程中只能提供下面两个关于绝对二次曲线的像ωC的两个实线性约束方程：
Figure BDA0001508929080000104
其中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)，结果如下：
Figure BDA0001508929080000111
Figure BDA0001508929080000112
Figure BDA0001508929080000113
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的系数矩阵形式如下：
Figure BDA0001508929080000121
对(26)进行Cholesky分解并求逆得到：
Figure BDA0001508929080000122
其中摄像机的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为光心的摄像机的内参数矩阵为
Figure FDA0002988394320000011
其中fx,fy为摄像机在x，y轴方向上的焦距，(u0,v0)为像平面上的主点坐标，s为摄像机的倾斜因子；
根据斜锥TSi的顶点OC和球像CSi来求出斜锥的方程，对于斜锥TSi其顶点是虚拟摄像机光心OC，坐标为(0,0,l)，准线是球像CSi；设点N(xi,yi,zi)(i＝1,2,3)为斜锥TSi(i＝1,2,3)上任一点，点
Figure FDA0002988394320000012
为准线CSi上任一点；点N(xi,yi,zi)(i＝1,2,3)N≠N0在锥面上的充分必要条件是：N在一条母线上，即得斜锥TSi(i＝1,2,3)的方程为，
Figure FDA0002988394320000013
其中u为一非零参数，
Figure FDA0002988394320000021
将上式整理得：
Figure FDA0002988394320000022
其中
Figure FDA0002988394320000023
记
Figure FDA0002988394320000024
即为所求斜锥TSi(i＝1,2,3)的方程；由斜锥TSi，其中i＝1,2,3，的方程和单位视球方程联立得小圆Si，其中i＝1,2,3，的方程：
Figure FDA0002988394320000025
其中
Figure FDA0002988394320000026
(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)；所以分别从三个球像上取三个点时就得到下面的方程组：
Figure FDA0002988394320000027
其中，δ表示为一个非零常数，
Figure FDA0002988394320000028
Figure FDA0002988394320000029
Figure FDA00029883943200000210
由上面的方程组解得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上提取的六个已知点，
Figure FDA0002988394320000031
由上式求出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，得
Figure FDA0002988394320000032
整理得Omi，其中i＝1,2,3，齐次坐标为
Figure FDA0002988394320000033
(3)求解绝对二次曲线的像
在求解在Omi之后根据极点极线关系求得lSi，从而求出lS1与CS1的两个交点，lS2和CS2的两个交点，lS3和CS3的两个交点，这六个点同时也在绝对二次曲线的像ωC上，即圆环点的像上；在这里，选取空间中位置不同的三个球得到的三个球像便求得绝对二次曲线的像ωC上的六个点，然后由这六个交点得ωC。
CN201711343962.7A 2017-12-15 2017-12-15 利用空间位置不同的三个球标定中心折反射摄像机的方法 Expired - Fee Related CN107958468B (zh)

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 CN107958468A (zh) 2018-04-24
CN107958468B true 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)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109272558B (zh) * 2018-11-28 2021-07-09 云南大学 分离圆的公共自极三角形及圆环点标定针孔摄像机的方法
CN110428472B (zh) * 2019-07-19 2021-11-09 中北大学 工业相机标定方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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 云南大学 使用双球的相交像及圆环点的像标定拋物折反射摄像机
CN106558082A (zh) * 2016-11-28 2017-04-05 云南大学 利用Veronese映射棋盘格的投影矩阵标定中心折反射摄像机
CN106558081A (zh) * 2016-11-28 2017-04-05 云南大学 标定非中心轴对称系统的圆锥折反射摄像机的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7362969B2 (en) * 2001-05-29 2008-04-22 Lucent Technologies Inc. Camera model and calibration procedure for omnidirectional paraboloidal catadioptric cameras

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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 云南大学 使用双球的相交像及圆环点的像标定拋物折反射摄像机
CN106558082A (zh) * 2016-11-28 2017-04-05 云南大学 利用Veronese映射棋盘格的投影矩阵标定中心折反射摄像机
CN106558081A (zh) * 2016-11-28 2017-04-05 云南大学 标定非中心轴对称系统的圆锥折反射摄像机的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Intrinsic parameter determination of a paracatadioptric camera by the intersection of two sphere projections;Zhao Y , Wang Y;《Journal of the Optical Society of America A Optics Image Science & Vision》;20151231;第32卷(第11期);第2201-2209页 *
Paracatadioptric camera calibration using sphere images;Duan H , Wu Y;《IEEE International Conference on Image Processing》;20111231;第1-4页 *
基于完全四点形利用空间球标定拋物折反射摄像机;顾万家;《中国优秀硕士学位论文全文数据库 信息科技辑》;20170215(第2期);第I138-3447页 *

Also Published As

Publication number Publication date
CN107958468A (zh) 2018-04-24

Similar Documents

Publication Publication Date Title
CN111414798B (zh) 基于rgb-d图像的头部姿态检测方法及系统
CN107886546B (zh) 利用球像及公共自极三角形标定抛物折反射摄像机的方法
CN107657645B (zh) 利用直线及圆的共轭直径的性质标定拋物折反射摄像机
CN109523595A (zh) 一种建筑工程直线棱角间距视觉测量方法
CN104517291A (zh) 基于目标同轴圆特征的位姿测量方法
CN105279758A (zh) 使用双球的相切像与圆环点的像标定拋物折反射摄像机
CN110675436A (zh) 基于3d特征点的激光雷达与立体视觉配准方法
CN109325983A (zh) 利用无穷远点关于圆极线性质标定拋物折反射摄像机
CN109360248A (zh) 利用单个球及共轭直径的性质标定拋物折反射摄像机
CN111080715B (zh) 三个球和无穷远点的极线性质标定摄像机内参数的方法
CN105321181A (zh) 使用双球的相离像与圆环点的像标定拋物折反射摄像机
CN107958468B (zh) 利用空间位置不同的三个球标定中心折反射摄像机的方法
CN108921904B (zh) 利用单个球及渐近线的性质标定针孔摄像机的方法
CN109360247B (zh) 单个球的公共自极三角形及正交消失点标定抛物摄像机
CN105303570A (zh) 使用双球的相交像及圆环点的像标定拋物折反射摄像机
CN109035342B (zh) 利用一条直线及圆环点极线标定拋物折反射摄像机的方法
CN109523598A (zh) 单个球的公共自极三角形及圆环点标定抛物摄像机的方法
CN113345033B (zh) 一种中心折反射摄像机内参数标定方法及系统
CN109712195A (zh) 利用球像的公共自极三角形进行单应性估计的方法
CN111429522B (zh) 利用共面圆的公共极点极线性质标定摄像机的方法和系统
CN109272558B (zh) 分离圆的公共自极三角形及圆环点标定针孔摄像机的方法
CN109191528B (zh) 利用球像与圆环点极线的性质标定针孔摄像机的方法
CN105354839A (zh) 使用相切的双球像与正交消失点标定拋物折反射摄像机
CN110428472B (zh) 工业相机标定方法
CN105354840A (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