单目相机标定实验(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

靶标图像如下