Qbery's blog

相机标定学习笔记(三)

Qbery

相机标定

发布于：2020年4月8日

单目相机标定实验（MATLAB）

MATLAB标定流程

通过MATLAB生成标定用的棋盘格，打印出来或显示在固定大小的平板上，生成的代码如下：
1
2
3
K = (checkerboard(100,4,3) > 0.5);
figure
imshow(K);
其中checkboard函数的参数100是格距，4是行数，3是列数
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：正交旋转矩阵（对应每个标定图案）

标定结果精度评价

图像重投影误差$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}}$$

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

归一化标定误差$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}}}$$

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

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

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

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

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

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

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

理想点到投影射线的距离误差$E_{ray}$

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

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

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

模型计算点之间的距离与理想距离的误差$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}||}}}$$

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

靶标平面测量

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

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

靶标图像如下

ARM Linux启动过程

stage1（汇编实现）

进入SVC模式，屏蔽所有的中断

设置 CPU 的速度和时钟频率

RAM 初始化

关闭 CPU 内部的指令／数据 Cache

拷贝stage2到RAM

设置堆栈指针sp

通过trampoline跳转stage2 main函数入口。

stage2（C实现）

初始化硬件设备

检测系统内存映射

加载映像

设置内核启动参数

调用内核

相机标定学习笔记(三)

单目相机标定实验（MATLAB）

MATLAB标定流程

MATLAB标定导出变量说明

标定结果精度评价

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

归一化标定误差$E_{ncs}$

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

理想点到投影射线的距离误差$E_{ray}$

模型计算点之间的距离与理想距离的误差$E_{dis}$

靶标平面测量

相机标定学习笔记(二)

相机标定的原理

非线性畸变的来源

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

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

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

建立畸变模型

求解畸变参数