Qbery

单目相机标定实验（MATLAB）

MATLAB标定流程

  1. 通过MATLAB生成标定用的棋盘格，打印出来或显示在固定大小的平板上，生成的代码如下：

    1
    2
    3
    K = (checkerboard(100,4,3) > 0.5);
    figure
    imshow(K);

    其中checkboard函数的参数100是格距，4是行数，3是列数

  2. Camera Calibrator APP，导入拍摄好的标定照片（这里按照MATLAB的默认设置拍摄了20张），输入测量得到的标定用棋盘格格距（其实格距在这里甚至后面的精度评估中都无影响，毕竟单目的标定丢失了投影方向的信息，不能用来测量）后进行标定，然后即可导出计算得到的标定参数。

MATLAB标定导出变量说明

ImageSize：图像尺寸；

RadialDistortion：径向畸变，即前述畸变参数里的$k_1、k_2；$

TangentialDistortion：切向畸变，即前述畸变参数里的$p_1、p_2$；

WorldPoints：世界坐标系下的理想角点坐标，即前述成像模型的$(x_w,y_w,0)$；

WorldUnits：世界坐标系下角点坐标单位；

EstimateSkew：评估偏斜标志，如果设为1，则图像的轴偏斜，设为0，则为垂直；

NumRadialDistortionCoefficients：径向畸变的阶数，一般为2或者3；

EstimateTangentialDistortion：评估切向畸变标志，如果设为1，则评估；

TranslationVectors：平移矢量（对应每个标定图案），此处为前述成像模型中的$t$的转置；

ReprojectionErrors：重投影误差（对应每张图案的每个点）；

RotationVectors：旋转矢量，本身是旋转轴，模值是旋转弧度，类似四元数，由于模型用的是正交旋转矩阵和平移向量，这里用不到此参数；

NumPatterns：标定图案数量，与平移矢量和旋转矩阵数目相当；

IntrinsicMatrix：内参矩阵，此处为前述成像模型中的A的转置；

FocalLength：内参矩阵中的$f_x,f_y$，归一化后的焦距；

PrincipalPoint：内参矩阵中的$u_0,v_0$，图像中心坐标；

Skew：内参矩阵中的$s$，偏斜系数；

MeanReprojectionError：重投影均方根误差；

ReprojectedPoints：世界坐标系下的理想角点重投影到像素坐标系下的坐标（即图像角点坐标）；

RotationMatrices：正交旋转矩阵（对应每个标定图案）

标定结果精度评价

  1. 图像重投影误差$E_{rms}$

    把世界坐标系下的已知角点理想三维坐标按照标定的外参和内参投影到像素坐标系上，计算角点重投影坐标$(\hat{x_u},\hat{y_u})$与图像上存在的角点坐标的均方根，MATLAB标定里的导出变量MeanReprojectionError即为此误差：

    $$E_{rms}=\frac 1n \sum_{n=1}^n{\sqrt{(x_{u,i}-\hat{x_{u,i}})^2 + (y_{u,i}-\hat{y_{u,i}})^2}}$$

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    function [Erms] = Cal_Erms(ReprojectedPoints, WorldPoints, RotationMatrices, TranslationVectors, IntrinsicMatrix)
    %%%%计算Erms
    [Pointcount,~,Piccount] = size(ReprojectedPoints);
    IntrinsicMatrix = IntrinsicMatrix';
    ErmsList = zeros(1,Piccount * Pointcount);
    for i = 1:Piccount
        for j = 1:Pointcount
        R = RotationMatrices(:,:,i)';
        T = TranslationVectors(i,:)';
        RT = [R(:,1) R(:,2) T];
        tmpOw = WorldPoints(j,:)';
        Ow = [tmpOw;1];
        Or0 = IntrinsicMatrix * RT * Ow;
        Or0 = Or0/Or0(3,1);
        Or = ReprojectedPoints(j,:,i)';
        dx = Or0(1,1) - Or(1,1);
        dy = Or0(2,1) - Or(2,1);   
        ErmsList(1,i * Pointcount + j - Pointcount) = sqrt(dx*dx+dy*dy);
        end
    end
    Erms = mean(ErmsList);
    end
  1. 归一化标定误差$E_{ncs}$

    由于$E_{rms}$受到图像分辨率、视场范围、空间点到相机的距离的不同程度影响，可以将整个图像投影到相机坐标系下，计算投影角点与相机坐标系下理论角点坐标形成的矩形面积的归一化均值来评价精度：

    $$E_{nce}=\frac 1n \sum_{n=1}^n{\sqrt{\frac{(x_{c,i}-\hat{x_{c,i}})^2 + (y_{c,i}-\hat{y_{c,i}})^2}{z_{c,i}^2(f_x^{-2} + f_y^{-2})/12}}}$$

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    function [Encs] = Cal_Encs(ReprojectedPoints, WorldPoints,RotationMatrices,TranslationVectors,IntrinsicMatrix)
    %%%%计算Encs
    fx = IntrinsicMatrix(1,1);
    fy = IntrinsicMatrix(2,2);
    Denominator = (fx^-2+fy^-2)/12;
    [Pointcount,~,Piccount] = size(ReprojectedPoints);
    IntrinsicMatrix = IntrinsicMatrix';
    EncsList = zeros(1, Piccount);
    for i = 1:Piccount
        for j = 1:Pointcount
        R = RotationMatrices(:,:,i)';
        T = TranslationVectors(i,:)';
        RT = [R(:,1) R(:,2) T];
        tmpOw = WorldPoints(j,:)';
        Ow = [tmpOw;1];
        Oc0 = RT * Ow;
        Oc0 = Oc0/Oc0(3,1);
        Or = ReprojectedPoints(j,:,i)';
        Orr= [Or;1];
        Oc = IntrinsicMatrix\Orr;
        dx = Oc0(1,1) - Oc(1,1);
        dy = Oc0(2,1) - Oc(2,1);   
        EncsList(1,i * Pointcount + j - Pointcount) = sqrt((dx*dx+dy*dy)/Denominator);
        end
    end
    Encs = mean(EncsList);
    end

    以下为未参与标定图像投影到世界坐标系下进行比较的精度评价方式

  2. 理想点与模型计算点距离误差$E_pt$

    计算空间角点与图像投影到靶标平面的坐标的均方根：

    $$E_{pt}=\frac 1n \sum_{n=1}^n{\sqrt{||M_{c,i}-\hat{M_{c,i}}||}}$$

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    function [RotationMatrices,TranslationVectors] = Cal_R_T(NewReprojectedPoints, NewWorldPoints, cameraParams)
    %%%%根据标定结果的内参计算未参与标定图像的外参
    [~,~,Piccount] = size(NewReprojectedPoints);
    RotationMatrices = zeros(3,3);
    TranslationVectors = zeros(1,3);
    for i =1:Piccount
        [R, T] = extrinsics(NewReprojectedPoints(:,:,i), NewWorldPoints, cameraParams);
        RotationMatrices(:,:,i) = R;
        TranslationVectors(:,:,i) = T;
    end
    end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function [Ept] = Cal_Ept(NewReprojectedPoints,NewWorldPoints,cameraParams)
%%%%计算未参与标定图像Ept，角点可以由detectCheckerboardPoints得到
[Pointcount,~,Piccount] = size(NewReprojectedPoints);
EptList = zeros(1, Piccount * Pointcount);
[RotationMatrices, TranslationVectors] = Cal_R_T(NewReprojectedPoints, NewWorldPoints,cameraParams);
for i = 1:Piccount
    Orwlist = pointsToWorld(cameraParams,RotationMatrices(:,:,i),TranslationVectors(:,:,i),NewReprojectedPoints(:,:,i));
    for j = 1:Pointcount
        Orw = Orwlist(j,:);
        dx = Orw(1,1) - NewWorldPoints(j,1);
        dy = Orw(1,2) - NewWorldPoints(j,2);
        EptList(1,i * Pointcount + j - Pointcount) = sqrt(dx*dx+dy*dy);
    end
end
Ept = mean(EptList);
end
  1. 理想点到投影射线的距离误差$E_{ray}$

    对图像进行畸变校正后，可以算出每个点投影到相机坐标系的射线$\hat{L_c}$，可以求出相机坐标系下理想角点到射线的距离的均值：

    $$E_{ray}=\frac 1n \sum_{n=1}^n{\sqrt{||M_{c,i}-\hat{L_{c,i}}||}}$$

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    function [Eray] = Cal_Eray(originimageFileNames, squareSize, cameraParams)
    [~,Piccount] = size(originimageFileNames);
    ErayPicList = zeros(1,Piccount);
    for i = 1:Piccount
        originimage = imread(originimageFileNames{i});
        %undistortImage默认OutputView为same，图像中心不变，Neworigin无须赋值
        [uimage, ~] = undistortImage(originimage,cameraParams);
        [Rpoints, boardSize] = detectCheckerboardPoints(uimage);
        worldPoints = generateCheckerboardPoints(boardSize, squareSize);
        [RotationMatrices, TranslationVectors] = extrinsics(Rpoints, worldPoints, cameraParams);
        RTMatrix = [RotationMatrices(1,:)' RotationMatrices(2,:)' TranslationVectors'];
        [Pointcount, ~] = size(worldPoints);
        ErayPointList = zeros(1,Pointcount);
        for j = 1:Pointcount
        WorldPoint = [worldPoints(j,:) 1]';
        OcWorldPoint = RTMatrix * WorldPoint;
        OcWorldPoint = OcWorldPoint/OcWorldPoint(3,1);
        Rpoint = [Rpoints(j,:) 1]';
        OcRpoint = cameraParams.IntrinsicMatrix'\Rpoint;
        d = norm(cross(OcRpoint,OcWorldPoint))/norm(OcRpoint);
        ErayPointList(1,j) = d;
        end
        ErayPicList(1,i) = mean(ErayPointList);
    end
    Eray = mean(ErayPicList);
    end
  1. 模型计算点之间的距离与理想距离的误差$E_{dis}$

    把图像角点投影到靶标平面，计算相邻点的距离与理想距离的差值的均值，这种评价方法若两计算点偏移量相近，则无法评定精度：

    $$E_{dis}=\frac 1m \sum_{n=1}^n{\sqrt{||M_{w,i}-M_{w,j}||-||\hat{M_{w,i}}-\hat{M_{w,j}||}}}$$

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    function [Edis] = Cal_Edis(originimageFileNames, squareSize, cameraParams)
    [~,Piccount] = size(originimageFileNames);
    EdisPicList = zeros(1,Piccount);
    for i = 1:Piccount
        originimage = imread(originimageFileNames{i});
        %undistortImage默认OutputView为same，图像中心不变，Neworigin无须赋值
        [uimage, ~] = undistortImage(originimage,cameraParams);
        [Rpoints, boardSize] = detectCheckerboardPoints(uimage);
        worldPoints = generateCheckerboardPoints(boardSize, squareSize);
        [RotationMatrices, TranslationVectors] = extrinsics(Rpoints, worldPoints, cameraParams);
        RWpoints = pointsToWorld(cameraParams, RotationMatrices, TranslationVectors, Rpoints);
        boardSize = boardSize - [1,1];
        rowid = boardSize(1,1);
        columnid = boardSize(1,2);
        EdisPointSum = 0;
        for j = 1:columnid
            for k = 1:rowid - 1
                P1 = RWpoints((j-1)*rowid+k,:);
                P2 = RWpoints((j-1)*rowid+k+1,:);
                d = sqrt((P1(1,1)-P2(1,1))^2 + (P1(1,2)-P2(1,2))^2);
                dd = d - squareSize;
                dd = abs(dd);
                EdisPointSum = dd + EdisPointSum;
            end
        end
        EdisPointSum = EdisPointSum/((rowid-1) * columnid);
        EdisPicList(1,i) = EdisPointSum;
    end
    Edis = mean(EdisPicList);
    end

靶标平面测量

在靶标平面上设置两个定位点（例如圆形的圆心）进行距离测量：通过拍摄含有定位点的靶标图像，根据标定结果计算外参，再把图像上定位点的像素坐标投影回世界坐标系，即可计算距离。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
function [delta, result] = Measure_dis(imagepath,squareSize,cameraParams,realdistance)
im = imread(imagepath);
[uim,~] = undistortImage(im, cameraParams);
uim_HSV = rgb2hsv(uim);
uim_Sat = uim_HSV(:,:,2);
Threshold = graythresh(uim_Sat);
Seg_im = (uim_Sat > Threshold);
[centers, ~] = imfindcircles(Seg_im,[10,60],'ObjectPolarity','bright','Sensitivity',0.75);
[Rpoints, boardSize] = detectCheckerboardPoints(uim);
worldPoints = generateCheckerboardPoints(boardSize, squareSize);
[RotationMatrices, TranslationVectors] = extrinsics(Rpoints, worldPoints, cameraParams);
RealRoundCenters = pointsToWorld(cameraParams, RotationMatrices, TranslationVectors, centers);
P1 = RealRoundCenters(1,:);
P2 = RealRoundCenters(2,:);
result = sqrt((P1(1,1)-P2(1,1))^2 + (P1(1,2)-P2(1,2))^2);
delta = abs(result - realdistance);
end

靶标图像如下

相机标定的原理

非线性畸变的来源

  1. 摄像机光学系统的透镜组不完善造成的径向畸变

    透镜系统的远光轴区域的放大率与光轴附近的放大率差异导致图像上的点向内或向外偏离光轴中心，若远光轴区域的放大率比光轴附近的大，则造成枕形畸变，若小于光轴附近的放大率，则造成桶形畸变。

    前两阶畸变量可以表示为：

    $$\begin{cases} δ_u=k_1u(u^2+v^2)+k_2u(u^2+v^2)^2 \\ δ_v=k_1v(u^2+v^2)+k_2v(u^2+v^2)^2 \end{cases}$$

  2. 不正确的镜头组合引起离心畸变

    实际相机的光学系统的光学中心和镜头各器件的光学中心不一致带来的畸变，既包含径向畸变也包含切向畸变。

    前两阶畸变量可表示为：

    $$\begin{cases} δ_u=p_1(3u^2+v^2)+2p_2uv \\ δ_v=p_2(u^2+3v^2)+2p_1uv \end{cases}$$

  3. 摄像机装配不完善造成薄透镜畸变(影响较小)

    透镜设计、生产不完善和安装误差造成的畸变，包含径向畸变和切向畸变。

    通常取一阶畸变量即可：

    $$\begin{cases} δ_u=q_1(u^2+v^2)\\ δ_v=q_2(u^2+v^2) \end{cases}$$

建立畸变模型

一般情况下建立畸变模型只需考虑主要畸变，取1个或两个主要畸变来源即可，高阶畸变量相对来说影响很小，可以忽略，因此此处建立的畸变模型仅包含一、二阶径向畸变和一、二阶离心畸变：

$$\begin{cases} x_d=x_u+ x_u(k_1r_u^2+k_2r_u^4)+p_1(r_u^2+2 x_u^2)+2p_2 x_u y_u \\ y_d=y_u+ y_u(k_1r_u^2+k_2r_u^4)+p_2(r_u^2+2 y_u^2)+2p_1 x_u y_u \end{cases}$$

其中$r_u^2=x_u^2+y_u^2$。

求解畸变参数

根据标定板数据可知世界坐标系下各个角点（靶标）的理论位置，利用无畸变模型求出的内外参，计算出各角点重投影误差，进行非线性优化（牛顿高斯法，Levenberg算法，LM等）找到畸变参数$k_1、k_2、p_1、p_2$的最优解：

$$min||m_i-\hat{m_i}(dx,dy,u_0,v_0,k_1,k_2,p_1,p_2)||$$

相机标定的目的

相机标定的最终目的是通过相机成像进行视觉测量，也就是通过图像的信息获取真实三维世界里的位置信息，而标定则是建立现实世界到图像的映射的过程，为了使测量可信，标定过程需要达成两个目的：

  1. 获取世界坐标系到像素坐标系的映射关系：获取相机的内参和外参；
  2. 获得相机的畸变参数：矫正畸变。

相机标定的原理

相机成像模型建立

相机成像是三次坐标系变换的过程：

  1. 世界坐标系到相机坐标系（刚体变换）

    $$\begin{bmatrix}x_c\\y_c\\z_c\\1\end{bmatrix}=\begin{bmatrix} R & t\\0^T & 1 \end{bmatrix} \begin{bmatrix} x_w\\y_w\\z_w\\1 \end{bmatrix}$$

    其中R是3×3 的正交旋转矩阵，t是3×1 的平移矢量，其实这个过程也可以用四元数或者欧拉角表达（这俩一般在游戏引擎里用得多），但是正交矩阵R和t在这里计算起来比较直观，可以设定世界坐标系里，标定平面与$x_w-y_w$平面重合，因此$z_w=0$，矩阵在计算时可进行以下简化：

    $$\begin{bmatrix} x_c \\ y_c \\ z_c \\ 1 \end{bmatrix} = \begin{bmatrix} R & t \\ 0^T & 1 \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ z_w \\ 1 \end{bmatrix} = \begin{bmatrix} r_1’ & r_2’ & r_3’ & t’ \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ 0 \\ 1 \end{bmatrix} = \begin{bmatrix} r_1 & r_2 & t’ \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ 1 \end{bmatrix}$$

  1. 相机坐标系到图像坐标系（中心透视投影）

    $$\begin{bmatrix}x_n\\y_n\\1\end{bmatrix}=\begin{bmatrix}f&0&0&0\\0&f&0&0\\0&0&1&0\end{bmatrix} \begin{bmatrix}x_c\\y_c\\z_c\\1\end{bmatrix}$$

    其中f为相机的焦距（这种中心投影模型仅适用于针孔相机模型，对广角不适用）。

  2. 图像坐标系到像素坐标系（仿射变换）

    $$\begin{bmatrix}x_u\\y_u\\1\end{bmatrix}=\begin{bmatrix}1/dx&0&u_0\\0&1/dy&v_0\\0&0&1\end{bmatrix}\begin{bmatrix}x_n\\y_n\\1\end{bmatrix}$$

    其中$u_0,v_0$为图像中心的像素坐标，$dx,dy$为图像每个像素的长宽。

其中1的转换矩阵$\begin{bmatrix} R & t\\0^T & 1 \end{bmatrix}$即为相机的外参，每张照片由于标定平面的位置不同，外参也不一样。

2的转换矩阵左乘上3的转换矩阵$\begin{bmatrix} f_x & 0 & u_0 & 0 \\ 0 & f_y & v_0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}$即为相机的内参，内参矩阵中(1, 2)位置的参数为扭曲参数，标准相机一般可将它设为0。

因此可以得到相机成像模型：

$$\begin{bmatrix} x_u \\ y_u \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & 0 & u_0 & 0 \\ 0 & f_y & v_0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} R & t \\ 0^T & 1 \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ z_w \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & 0 & u_0 & 0 \\ 0 & f_y & v_0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} r_1’ & r_2’ & t’ \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & 0 & u_0 \\ 0 & f_y & v_0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} r_1 & r_2 & t \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ 1 \end{bmatrix}$$

记内参矩阵为A，投影矩阵为H(3x3)，则有：

$$z\begin{bmatrix}x_u\\y_u\\1\end{bmatrix}=A\begin{bmatrix}r_1&r_2&t\end{bmatrix}\begin{bmatrix}x_w\\y_w\\1\end{bmatrix}=H\begin{bmatrix}x_w\\y_w\\1\end{bmatrix}=\begin{bmatrix}h_1&h_2&h_3\end{bmatrix}\begin{bmatrix}x_w\\y_w\\1\end{bmatrix}$$

其中z为尺度因子，为简化计算提取出来。

求解内外参方法

求解H矩阵

根据多视几何的结论，可以令$x=\begin{bmatrix}h_1^{-T}\\h_2^{-T}\\h_3^{-T}\end{bmatrix}$，$M=\begin{bmatrix}x_w\\y_w\\1\end{bmatrix}$，则上式可化为：

$$\begin{bmatrix}M^T & 0^T & -x_u M^T\\0^T & M^T & -y_u M^T \end{bmatrix}x=0$$

设靶标上由n个可以确定位置的点，则有n组这样的方程组，组合起来为2n×9的矩阵，因为H为齐次矩阵，9个未知数只有8个是独立的，因此至少只需要4个点即可进行求解，可由奇异值分解求出x，即可解出投影矩阵H。

求解内外参

求出H矩阵后，根据R的正交性，有$r_1•r_2=0$$和$$r_1^Tr_1=r_2^Tr_2$，即有方程组：

$$\begin{cases} h_1^TA^{-T}A^{-1}h_2=0 \\ h_1^TA^{-T}A^{-1}h_1=h_2^TA^{-T}A^{-1}h_2 \end{cases}$$

令$B_{(3×3)}=A^{-T}A^{-1}$共有6个未知量，用b表示为：

$$b=\begin{bmatrix}1/f_x^2 & 0 & 1/f_y^2 & -u_0/f_x^2 & -v_0/f_y^2 & {u_0^2/f_x^2+v_0^2/f_y^2+1}\end{bmatrix}^T$$

可推出：

$$h_i^TBh_j=v_{ij}^Tb$$

$$v_{ij}=\begin{bmatrix} h_{i1}h{j1} & h_{i1}h_{j2} + h_{i2}h_{j1} & h_{i2}h_{j2} & h_{i3}h_{j1} + h_{i1}h_{j3} & h_{i3}h_{j2} + h_{i2}h_{j3} & h_{i3}h_{j3} \end{bmatrix}^T$$

因此有：

$$\begin{bmatrix} v_{12}^T \\ (v_{11}-v_{22})^T \end{bmatrix}b=0$$

设有N幅图像，则有N组这样的方程组，组合起来为2N×6的矩阵，可通过奇异值分解求出b，再通过Cholesky分解即可得到内参矩阵A。

A求出后根据$A\begin{bmatrix}r_1&r_2&t\end{bmatrix}=H$和$r_3=r_1×r_2$即可求出外参。