行业知识
单目相机标定实现--张正友标定法
文章目录
- 一：相机坐标系，像素平面坐标系，世界坐标系，归一化坐标系介绍
- 1：概述
- 公式
- 二:实现
- 1：整体流程
- 4：求出每张图像的单应性矩阵并用LMA优化
- 5：求解理想无畸变情况下的摄像机的内参数和外参数
- 6：应用最小二乘求出实际的畸变系数
- 7：综合内参，外参，畸变系数，使用极大似然法(LMA)，优化估计，提升估计精度
- 8：得出相机的内参，外参和畸变系数
- 9：OpenCV模型
- 三：畸变修复（去畸变）
- 四：���结
原文链接：地址
个人笔记：
本次介绍针对于单目相机标定，实现方法：张正友标定法。
一：相机坐标系，像素平面坐标系，世界坐标系，归一化坐标系介绍
1：概述
如图，现实世界中有一个P点和一个相机（光心），描述这个P点的空间坐标首先得有一个坐标系，那么以光心为原点O建一个坐标系，叫相机坐标系。
那么就可以在相机坐标系下，设 P 坐标 ( X , Y , Z ) P坐标(X,Y,Z) P坐标(X,Y,Z)和P的投影点 P ′ ( x ′ , y ′ , z ′ ) P'(x',y',z') P′(x′,y′,z′)。值得一提的是， P ′ ( x ′ , y ′ , z ′ ) P'(x',y',z') P′(x′,y′,z′)坐落在物理成像平面和像素平面。
物理成像平面，像素平面是二维的，他们的坐标系并不一样：
物理成像平面在 O ′ ( x ′ , y ′ ) O'(x',y') O′(x′,y′)平面上；
像素平面的原点在那个黑灰色图的左上角(图片的左上角)，横轴向右称为 u u u轴，纵轴向下称为 v v v轴。
这样就得到了 P ′ P' P′的像素坐标 P ( u , v ) P(u,v) P(u,v)，称为 P u v Puv Puv。
所谓的归一化的成像平面，就是将三维空间点的坐标都除以Z，在相机坐标系下，P有X, Y, Z 三个量，如果把它们投影到归一化平面 Z = 1 上，就会得到P的归一化坐标P(X/Z, Y/Z, 1)。
公式
[ X Y Z 1 ] − > \left[\begin{array}{c} X \\ Y \\ Z \\ 1 \end{array}\right]-> XYZ1 −>物体坐标.
[ R T 0 1 ] − > \left[\begin{array}{cc} R & T \\ 0 & 1 \end{array}\right]-> [R0T1]−>外参
[ α γ u 0 0 0 β v 0 0 0 0 1 0 ] − > \left[\begin{array}{cccc} \alpha & \gamma & u_{0} & 0 \\ 0 & \beta & v_{0} & 0 \\ 0 & 0 & 1 & 0 \end{array}\right]-> α00γβ0u0v01000 −>内参
[ u v 1 ] − > \left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]-> uv1 −>像素坐标
其中外参 T T T是平移向量 ( t 1 , t 2 , t 3 ) T (t1,t2,t3)^T (t1,t2,t3)T.
R R R旋转矩阵如下图：
二:实现
1：整体流程
第1步，第2步，第3步 暂不介绍了（可以用halcon算子块或者OpenCV获取特征点坐标），主要介绍获取到特征点以后，优化获取参数部分。
4：求出每张图像的单应性矩阵并用LMA优化
如何计算单应性矩阵：
[ x b y b w b ] = [ h 1 h 2 h 3 h 4 h 5 h 6 h 7 h 8 h 9 ] [ x a y a w a ] x b w b = h 1 x a + h 2 y a + h 3 w a h 7 x a + h 8 y a + h 9 w a = h 1 x a / w a + h 2 y a / w a + h 3 h 7 x a / w a + h 8 y a / w a + h 9 y b w b = h 4 x a + h 5 y a + h 6 w a h 7 x a + h 8 y a + h 9 w a = h 4 x a / w a + h 5 y a / w a + h 6 h 7 x a / w a + h 8 y a / w a + h 9 \begin{array}{c} \left[\begin{array}{l} x_{b} \\ y_{b} \\ w_{b} \end{array}\right]=\left[\begin{array}{ccc} \mathrm{h}_{1} & \mathrm{~h}_{2} & \mathrm{~h}_{3} \\ \mathrm{~h}_{4} & \mathrm{~h}_{5} & \mathrm{~h}_{6} \\ \mathrm{~h}_{7} & \mathrm{~h}_{8} & \mathrm{~h}_{9} \end{array}\right]\left[\begin{array}{l} x_{a} \\ y_{a} \\ w_{a} \end{array}\right] \\\\ \frac{x_{b}}{w_{b}}=\frac{h_{1} x_{a}+\mathrm{h}_{2} y_{a}+\mathrm{h}_{3} w_{a}}{h_{7} x_{a}+\mathrm{h}_{8} y_{a}+\mathrm{h}_{9} w_{a}}=\frac{h_{1} x_{a} / w_{a}+\mathrm{h}_{2} y_{a} / w_{a}+\mathrm{h}_{3}}{h_{7} x_{a} / w_{a}+\mathrm{h}_{8} y_{a} / w_{a}+\mathrm{h}_{9}} \\ \frac{y_{b}}{w_{b}}=\frac{h_{4} x_{a}+\mathrm{h}_{5} y_{a}+\mathrm{h}_{6} w_{a}}{h_{7} x_{a}+\mathrm{h}_{8} y_{a}+\mathrm{h}_{9} w_{a}}=\frac{h_{4} x_{a} / w_{a}+\mathrm{h}_{5} y_{a} / w_{a}+\mathrm{h}_{6}}{h_{7} x_{a} / w_{a}+\mathrm{h}_{8} y_{a} / w_{a}+\mathrm{h}_{9}} \end{array} xbybwb = h1 h4 h7 h2 h5 h8 h3 h6 h9 xayawa wbxb=h7xa+h8ya+h9wah1xa+h2ya+h3wa=h7xa/wa+h8ya/wa+h9h1xa/wa+h2ya/wa+h3wbyb=h7xa+h8ya+h9wah4xa+h5ya+h6wa=h7xa/wa+h8ya/wa+h9h4xa/wa+h5ya/wa+h6
令 u a = x a w a , v a = y a w a , u b = x b w b , v b = y b w b , 上式化简为 \text { 令 } u_{a}=\frac{x_{a}}{w_{a}}, v_{a}=\frac{y_{a}}{w_{a}}, u_{b}=\frac{x_{b}}{w_{b}}, v_{b}=\frac{y_{b}}{w_{b}} \text {, 上式化简为 } 令 ua=waxa,va=waya,ub=wbxb,vb=wbyb, 上式化简为
u b = h 1 u a + h 2 v a + h 3 h 7 u a + h 8 v a + h 9 v b = h 4 u a + h 5 v a + h 6 h 7 u a + h 8 v a + h 9 \begin{array}{l} u_{b}=\frac{h_{1} u_{a}+\mathrm{h}_{2} v_{a}+\mathrm{h}_{3}}{h_{7} u_{a}+\mathrm{h}_{8} v_{a}+\mathrm{h}_{9}} \\ v_{b}=\frac{h_{4} u_{a}+\mathrm{h}_{5} v_{a}+\mathrm{h}_{6}}{h_{7} u_{a}+\mathrm{h}_{8} v_{a}+\mathrm{h}_{9}} \end{array} ub=h7ua+h8va+h9h1ua+h2va+h3vb=h7ua+h8va+h9h4ua+h5va+h6
h 1 u a + h 2 v a + h 3 − h 7 u a u b − h 8 v a u b − h 9 u b = 0 h 4 u a + h 5 v a + h 6 − h 7 u a v b − h 8 v a v b − h 9 v b = 0 \begin{array}{c} h_{1} u_{a}+h_{2} v_{a}+h_{3}-h_{7} u_{a} u_{b}-h_{8} v_{a} u_{b}-h_{9} u_{b}=0 \\ h_{4} u_{a}+h_{5} v_{a}+h_{6}-h_{7} u_{a} v_{b}-h_{8} v_{a} v_{b}-h_{9} v_{b}=0 \end{array} h1ua+h2va+h3−h7uaub−h8vaub−h9ub=0h4ua+h5va+h6−h7uavb−h8vavb−h9vb=0
可以直接设 ∥ h ∥ 2 2 = 1 ，此时仍然可以固定住尺度，且有: \text { 可以直接设 }\|h\|_{2}^{2}=1 \text { ，此时仍然可以固定住尺度，且有: } 可以直接设 ∥h∥22=1 ，此时仍然可以固定住尺度，且有:
此时系数矩阵秩为8，有线性空间解，解的自由度为1，满足齐次性，又由于限制单位长度，有唯一解，此时仍可以通过SVD分解求解出 h，从而得到单应矩阵。
代码实现：
//获取标准差double CameraCalibration::StdDiffer(const Eigen::VectorXd& data){ //获取平均值 double mean = data.mean(); //std::sqrt((Σ(x-_x)²) / n) return std::sqrt((data.array() - mean).pow(2).sum() / data.rows());}// 归一化Eigen::Matrix3d CameraCalibration::Normalization (const Eigen::MatrixXd& P){ Eigen::Matrix3d T; double cx = P.col ( 0 ).mean(); double cy = P.col ( 1 ).mean(); double stdx = StdDiffer(P.col(0)); double stdy = StdDiffer(P.col(1));; double sqrt_2 = std::sqrt ( 2 ); double scalex = sqrt_2 / stdx; double scaley = sqrt_2 / stdy; T Eigen::Matrix3d realNormMat = Normalization(srcNormal); Eigen::Matrix3d picNormMat = Normalization(dstNormal); int n = srcNormal.rows(); // 2. DLT Eigen::MatrixXd input ( 2*n, 9 ); for ( int i=0; i //转换齐次坐标 Eigen::MatrixXd singleSrcCoor(3,1),singleDstCoor(3,1); singleSrcCoor(0,0) = srcNormal(i,0); singleSrcCoor(1,0) = srcNormal(i,1); singleSrcCoor(2,0) = 1; singleDstCoor(0,0) = dstNormal(i,0); singleDstCoor(1,0) = dstNormal(i,1); singleDstCoor(2,0) = 1; //坐标归一化 Eigen::MatrixXd realNorm(3,1),picNorm(3,1); realNorm = realNormMat * singleSrcCoor; picNorm = picNormMat * singleDstCoor; input ( 2*i, 0 ) = realNorm ( 0, 0 ); input ( 2*i, 1 ) = realNorm ( 1, 0 ); input ( 2*i, 2 ) = 1.; input ( 2*i, 3 ) = 0.; input ( 2*i, 4 ) = 0.; input ( 2*i, 5 ) = 0.; input ( 2*i, 6 ) = -picNorm ( 0, 0 ) * realNorm ( 0, 0 ); input ( 2*i, 7 ) = -picNorm ( 0, 0 ) * realNorm ( 1, 0 ); input ( 2*i, 8 ) = -picNorm ( 0, 0 ); input ( 2*i+1, 0 ) = 0.; input ( 2*i+1, 1 ) = 0.; input ( 2*i+1, 2 ) = 0.; input ( 2*i+1, 3 ) = realNorm ( 0, 0 ); input ( 2*i+1, 4 ) = realNorm ( 1, 0 ); input ( 2*i+1, 5 ) = 1.; input ( 2*i+1, 6 ) = -picNorm ( 1, 0 ) * realNorm ( 0, 0 ); input ( 2*i+1, 7 ) = -picNorm ( 1, 0 ) * realNorm ( 1, 0 ); input ( 2*i+1, 8 ) = -picNorm ( 1, 0 ); } // 3. SVD分解 JacobiSVDpublic: Eigen::VectorXd operator()(const Eigen::VectorXd& parameter,const QList Eigen::MatrixXd inValue = otherArgs.at(0); Eigen::MatrixXd outValue = otherArgs.at(1); int dataCount = inValue.rows(); //保存残差值 Eigen::VectorXd residual(dataCount*2) ,residualNew(dataCount*2); Eigen::Matrix3d HH = parameter.reshaped //转换齐次坐标 Eigen::VectorXd singleRealCoor(3),U(3); singleRealCoor(0,0) = inValue(i,0); singleRealCoor(1,0) = inValue(i,1); singleRealCoor(2,0) = 1; U = HH * singleRealCoor; U /= U(2); residual(i*2) = U(0); residual(i*2+1) = U(1); residualNew(i*2) = outValue(i,0); residualNew(i*2+1) = outValue(i,1); } residual -= residualNew; return residual; }}; //求偏导1 double PartialDeriv_1(const Eigen::VectorXd& parameter,int paraIndex,const Eigen::MatrixXd &inValue,int objIndex) { Eigen::VectorXd para1 = parameter; Eigen::VectorXd para2 = parameter; para1(paraIndex) -= DERIV_STEP; para2(paraIndex) += DERIV_STEP; //逻辑 double obj1 = ((para1(0))*inValue(objIndex,0) + para1(1)*inValue(objIndex,1) + para1(2)) / (para1(6)*inValue(objIndex,0) + para1(7)*inValue(objIndex,1) + para1(8)); double obj2 = ((para2(0))*inValue(objIndex,0) + para2(1)*inValue(objIndex,1) + para2(2)) / (para2(6)*inValue(objIndex,0) + para2(7)*inValue(objIndex,1) + para2(8));; return (obj2 - obj1) / (2 * DERIV_STEP); } //求偏导2 double PartialDeriv_2(const Eigen::VectorXd& parameter,int paraIndex,const Eigen::MatrixXd &inValue,int objIndex) { Eigen::VectorXd para1 = parameter; Eigen::VectorXd para2 = parameter; para1(paraIndex) -= DERIV_STEP; para2(paraIndex) += DERIV_STEP; //逻辑 double obj1 = ((para1(3))*inValue(objIndex,0) + para1(4)*inValue(objIndex,1) + para1(5)) / (para1(6)*inValue(objIndex,0) + para1(7)*inValue(objIndex,1) + para1(8)); double obj2 = ((para2(3))*inValue(objIndex,0) + para2(4)*inValue(objIndex,1) + para2(5)) / (para2(6)*inValue(objIndex,0) + para2(7)*inValue(objIndex,1) + para2(8));; return (obj2 - obj1) / (2 * DERIV_STEP); }public: Eigen::MatrixXd operator()(const Eigen::VectorXd& parameter,const QList Eigen::MatrixXd inValue = otherArgs.at(0); int rowNum = inValue.rows(); Eigen::MatrixXd Jac(rowNum*2, parameter.rows()); for (int i = 0; i
上一篇：揭秘JS打造万年历的绝招！
