【理论+程序】摄像机标定之图像重投误差、归一化标定误差、空间点与模型计算点的距离误差、空间点到投影射线的距离误差、模型计算点之间的距离与标准距离的误差

最新推荐文章于 2023-01-05 10:05:46 发布

什么是快乐星球！？

最新推荐文章于 2023-01-05 10:05:46 发布

阅读量6.4k

点赞数 10

分类专栏：计算机视觉文章标签：计算机视觉 matlab

本文链接：https://blog.csdn.net/snail9610/article/details/105727333

版权

计算机视觉专栏收录该内容

7 篇文章 2 订阅

订阅专栏

摄像机标定结果精度评价

文章目录

摄像机标定结果精度评价
1 图像重投误差 $E_{rms}$
2 归一化标定误差 $E_{nce}$
3 空间点与模型计算点的距离误差 $E_{pt}$
4 空间点到投影射线的距离误差 $E_{ray}$
5 模型计算点之间的距离与标准距离的误差 $E_{dis}$

对于摄像机的标定已经有了很多博客，但对其精度评价介绍的不够全面。对于已经获取的摄像机参数，需要采用合适的方法对其标定精度进行评价，这里给出五种评价方法和其Matlab代码。

目前常用的两种标定精度评价方法为：

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

将已知空间三维坐标的特征点按照已标定的摄像机参数重投回图像平面，得到重投点的像素坐标 $\hat x_u,\hat y_u)$ ，计算其与真实提取的像素坐标 $x_u,y_u)$ 之间的距离：

在这里插入图片描述
这个误差是Matlab Camera Calibrator APP自带的误差，所以这里就没有编写其代码。

2 归一化标定误差 $E_{nce}$

图像分辨率、视场范围以及空间点到摄像机的距离将对重投误差 $E_{rms}$ 有不同程度的影响，为了消除上述不利因素，将每个像素点重投回三维空间中的平面上，形成一定大小的矩形区域，其面积大小表示了像素的分辨率在该距离上的不确定度，则可以利用矩形面积的方差值来对误差作归一化处理。在摄像机坐标系下，设空间点的三维坐标为 $x_c,y_c,z_c)$ ，由图像像素坐标重投回空间平面 $Z=z_k$ 形成的点的三维坐标为 $(\hat x_u,\hat y_u,\hat z_u)$ ，有效焦距为 $f_x$ 和 $f_y$ ，则 $E_{nce}$ 可表示为：

在这里插入图片描述

function Ence = showEnce(cameraParams)
% 计算归一化标定误差Encs
% 输入参数cameraParams是一个结构体,用到的结构体中的参数有ReprojectedPoints、WorldPoints、
% RotationMatrices、TranslationVectors、IntrinsicMatrix
% 输出参数为归一化标定误差，是所有图像归一化标定误差的平均值，单位是像素
% 

ReprojectedPoints = cameraParams.ReprojectedPoints;
WorldPoints = cameraParams.WorldPoints;
RotationMatrices = cameraParams.RotationMatrices;
TranslationVectors = cameraParams.TranslationVectors;
IntrinsicMatrix = cameraParams.IntrinsicMatrix;

[~,~,NumPicture] = size(RotationMatrices);% NumPicture：输入图像的数量
[NumPoint,~,~] = size(ReprojectedPoints);% NumPoint：每幅图像重投影点的数量
IntrinsicMatrix = IntrinsicMatrix';% 转化为标准的内参矩阵
fx = IntrinsicMatrix(1,1);
fy = IntrinsicMatrix(2,2); % 获取x和y轴方向的焦距
zci = 1; % 归一化平面，投影平面是Z=1
Denominator = zci*(fx^-2+fy^-2)/12;% Encs计算公式的分母

iEachPointEncs = zeros(1,NumPoint);% 第i个图像每个点的归一化误差
iEachPicture = zeros(1, NumPicture);% 第i个图像的归一化误差(平均值)
for i = 1:NumPicture
    for j = 1:NumPoint
        RotationMatrix = RotationMatrices(:,:,i)';% 第i个图像的正交旋转矩阵
        TranslationVector = TranslationVectors(i,:)';% 第i个图像的平移向量
        %PositionMatrix:位姿矩阵，只需要正交旋转矩阵的前两列即可
        PositionMatrix = [RotationMatrix(:,1) RotationMatrix(:,2) TranslationVector];
        WorldPoint_temp = WorldPoints(j,:)';
        WorldPoint = [WorldPoint_temp;1];% 齐次坐标系
        CameraPoint = PositionMatrix * WorldPoint;% 摄像机坐标系
        NormalizedCameraPoint = CameraPoint/CameraPoint(3,1);% 归一化的摄像机坐标系
        PixelPoint_temp = ReprojectedPoints(j,:,i)';
        PixelPoint = [PixelPoint_temp;1];% 齐次坐标系
        NormalizedPixelPoint = inv(IntrinsicMatrix)*PixelPoint;% 图像像素坐标重投回空间平面zc=1的点
        Delta_x = NormalizedCameraPoint(1,1) - NormalizedPixelPoint(1,1);
        Delta_y = NormalizedCameraPoint(2,1) - NormalizedPixelPoint(2,1);   
        iEachPointEncs(j) = sqrt((Delta_x^2+Delta_y^2)/Denominator);
    end
    iEachPicture(i) = mean(iEachPointEncs(j));
end
Ence = mean(iEachPicture);

% 绘制条形图
figure
bar(iEachPicture);
title('Ence')
% figure
% ax = bar(iEachPicture);
% hold on
% plot(1:(NumPicture+ax.BarWidth),repmat(Encs,[1,(NumPicture+ax.BarWidth)]))
end

图像重投误差 $E_{rms}$ 和归一化标定误差 $E_{nce}$ 两种精度评价手段所采用的特征点坐标均参与了摄像机参数的计算过程，然而实际的视觉测量阶段所拍摄的图像和特征点均未参与摄像机参数的标定过程，因此，对摄像机标定结果的测量精度评价应采用未参与标定的图像和特征点坐标，并且理应在测量空间中对点的三维坐标进行比较。下面给出了三种在测量空间中评价标定结果精度的方法。

3 空间点与模型计算点的距离误差 $E_{pt}$

采用未标定的若干幅图像作为测试图像，运用参数的标定数值，计算出每幅测试图像的靶标平面位姿参数，再将每幅测试图像的图像特征点投射到靶标平面形成交点 $\hat {M_c}$ ，计算其与所有空间特征点 $\hat {M_c}$ 之间的距离的均方根误差：
在这里插入图片描述

function Ept = showEpt(NewImageFileNames,cameraParams)
% 计算空间点与模型计算点的距离误差Ept
%

% 在图像中检测棋盘格
[NewImagePoints, NewBoardSize, NewImagesUsed] = detectCheckerboardPoints(NewImageFileNames);
NewImageFileNames = NewImageFileNames(NewImagesUsed);

% 生成新的世界坐标
NewSquareSize = 10;  % 单位：mm
NewWorldPoints = generateCheckerboardPoints(NewBoardSize, NewSquareSize);

% 读取测试图像的数量
[~,NewPictureNum] = size(NewImageFileNames);

% 读取测试图像点的数量
% [NewPicturePointNum,~] = size(NewWorldPoints);

% 计算新的外参 & 映射图像指向X-Y平面上的世界坐标
NewRotationMatrix = zeros(3,3,NewPictureNum);
NewTranslationVector = zeros(1,3,NewPictureNum);
transposeMatrix = zeros(3,3,NewPictureNum);

for i=1:NewPictureNum
    [NewRotationMatrix(:,:,i), NewTranslationVector(:,:,i)] = extrinsics(...
    NewImagePoints(:,:,i),NewWorldPoints,cameraParams);
    % 计算靶标平面位姿参数
    % NewPositionMatrix:位姿矩阵，只需要正交旋转矩阵的前两列即可
    % 这个地方需要对旋转矩阵进行转置
    transposeMatrix(:,:,i) = NewRotationMatrix(:,:,i)';
    % 可以输出每个图像的位姿矩阵NewPositionMatrix
    NewPositionMatrix(:,:,i) = [transposeMatrix(:,1,i) transposeMatrix(:,2,i) NewTranslationVector(:,:,i)'];

    % NewRotationMatrix = NewRotationMatrix';
    % NewTranslationVector = NewTranslationVector';

    % 映射图像指向X-Y平面上的世界坐标
    ReproductedNewWorldPoints(:,:,i) = pointsToWorld(cameraParams,NewRotationMatrix(:,:,i),...
        NewTranslationVector(:,:,i),NewImagePoints(:,:,i));
    
    % 计算空间点与模型计算点的距离误差Ept
    DeltaMatrix(:,:,i) = NewWorldPoints-ReproductedNewWorldPoints(:,:,i);
    SquareMatrix(:,:,i) = DeltaMatrix(:,:,i).^2;
    ErrorVector(i) = mean(sqrt(SquareMatrix(:,1,i)+SquareMatrix(:,2,i)));
end
Ept = mean(ErrorVector);

% 绘制条形图
figure
bar(ErrorVector)
title('Ept')
end

4 空间点到投影射线的距离误差 $E_{ray}$

运用获取的摄像机参数，对图像点进行畸变校正并投射到三维空间形成射线 $\hat L_c$ 。在摄像机坐标系下，计算已知空间点 $M_c$ 到对应射线 $L_c$ 的距离：
在这里插入图片描述

function Eray = showEray(NewImageFileNames,cameraParams)
% 计算空间点到投影射线的距离误差

squareSize=10;

[~,NumPicture] = size(NewImageFileNames);
% str = strings(3,10);
% strs = strings(3,10,5);% 存储每个直线方程
NumPoint = 121;
iEachPointEray = zeros(1,NumPoint);% 第i个图像每个点的Eray
iEachPicture = zeros(1, NumPicture);% 第i个图像的Eray(平均值)
for i = 1:NumPicture
    OriginImage = imread(NewImageFileNames{i});
    [Uimage,newOrigin] = undistortImage(OriginImage,cameraParams);
    [ImagePoints, boardSize] = detectCheckerboardPoints(Uimage);
    worldPoints = generateCheckerboardPoints(boardSize, squareSize);
    [RotationMatrices, TranslationVectors] = extrinsics(ImagePoints, worldPoints, cameraParams);
    PositionMatrix = [RotationMatrices(1,:)' RotationMatrices(2,:)' TranslationVectors'];
    [PointNum, ~] = size(worldPoints);
    ErayPointList = zeros(1,PointNum);
    for j = 1:PointNum
        WorldPoint = [worldPoints(j,:) 1]';
        WorldtoCameraPoint = PositionMatrix * WorldPoint;
        PixelPoint = [ImagePoints(j,:) 1]';
        NormalizedPixelPoint = cameraParams.IntrinsicMatrix'\PixelPoint;
        % https://www.ilovematlab.cn/thread-102421-1-1.html
        % 这里利用了直线上的原点坐标和NormalizedPixelPoint坐标
        distance = norm(cross(NormalizedPixelPoint,WorldtoCameraPoint))/norm(NormalizedPixelPoint);
        iEachPointEray(j) = distance;
    end
    iEachPicture(i) = mean(iEachPointEray);
end
Eray = mean(iEachPicture);
% 绘制条形图
figure
bar(iEachPicture);
title('Eray');
end

5 模型计算点之间的距离与标准距离的误差 $E_{dis}$

运用获取的摄像机参数，计算每幅测试图像的图像特征点投射到靶标平面形成的交点 $\hat {M_c}$ 。在摄像机坐标系下，计算相邻格点之间的三维距离 $\parallel \hat{M_{c,i}}-\hat{M_{c,j}} \parallel$ ，并与已知的靶标最小格点距离 ${M_{w,i}}-{M_{w,j}} \parallel$ 进行比较：
在这里插入图片描述
模型计算点之间的距离与标准距离的误差 $E_{dis}$ 不适合用于测量精度的评价。这是由于当摄像机参数出现偏差时，计算出的相邻两点的 $\hat {M_c}$ 可能同时出现一个较为接近的偏移，导致计算出的距离也是比较准确的，从而无法准确评价摄像机标定结果。

function Edis = showEdis(NewImageFileNames,cameraParams)
% 计算空间点与模型计算点的距离误差Ept

% 在图像中检测棋盘格
[NewImagePoints, NewBoardSize, NewImagesUsed] = detectCheckerboardPoints(NewImageFileNames);
NewImageFileNames = NewImageFileNames(NewImagesUsed);

% 生成新的世界坐标
NewSquareSize = 10;  % 单位：mm
NewWorldPoints = generateCheckerboardPoints(NewBoardSize, NewSquareSize);

% 读取测试图像的数量
[~,NewPictureNum] = size(NewImageFileNames);

% 读取测试图像点的数量
[NewPicturePointNum,~] = size(NewWorldPoints);

% 计算新的外参 & 映射图像指向X-Y平面上的世界坐标
NewRotationMatrix = zeros(3,3,NewPictureNum);
NewTranslationVector = zeros(1,3,NewPictureNum);
% transposeMatrix = zeros(3,3,NewPictureNum);

for i=1:NewPictureNum

    [NewRotationMatrix(:,:,i), NewTranslationVector(:,:,i)] = extrinsics(...
    NewImagePoints(:,:,i),NewWorldPoints,cameraParams);
    
    % NewRotationMatrix = NewRotationMatrix';
    % NewTranslationVector = NewTranslationVector';

    % 映射图像指向X-Y平面上的世界坐标
    ReproductedNewWorldPoints(:,:,i) = pointsToWorld(cameraParams,NewRotationMatrix(:,:,i),...
        NewTranslationVector(:,:,i),NewImagePoints(:,:,i));
       
    % 计算Edis：y轴方向上的距离（误差）
    k = 1;
    for j=1:(NewPicturePointNum-1)
        DeltaMatrix(j,:) = ReproductedNewWorldPoints(j,:,i)-ReproductedNewWorldPoints(j+1,:,i);
        SquareMatrix(j,:) = DeltaMatrix(j,:).^2;
        Distance = sqrt(SquareMatrix(j,1)+SquareMatrix(j,2));
        if Distance<40
            yDistanceVector(k) = Distance;
            k = k+1;
        end
    end
    
    % 计算Edis：x轴方向上的距离（误差）
    
    k = 1;
    % 交换ReproductedNewWorldPoints纵横坐标(不是转置)
    % ReproductedNewWorldPoints(:,[1 2])=ReproductedNewWorldPoints(:,[2 1]);
    % xReproductedNewWorldPoints = permute(ReproductedNewWorldPoints,[2 1 3]);
    tempReproductedNewWorldPoints(:,:) = ReproductedNewWorldPoints(:,:,i);
    tempReproductedNewWorldPoints(:,[1 2]) = tempReproductedNewWorldPoints(:,[2 1]);
    xReproductedNewWorldPoints(:,:,i) = tempReproductedNewWorldPoints;
    for j=1:(NewPicturePointNum-1)
        DeltaMatrix(j,:) = xReproductedNewWorldPoints(j,:,i)-xReproductedNewWorldPoints(j+1,:,i);
        SquareMatrix(j,:) = DeltaMatrix(j,:).^2;
        Distance = sqrt(SquareMatrix(j,1)+SquareMatrix(j,2));
        if Distance<10
            xDistanceVector(k) = Distance;
            k = k+1;
        end
    end
    
    DistanceVector = [yDistanceVector,xDistanceVector];
   
    iPictureEdis(i) = mean(sqrt((10 - DistanceVector).^2));

end
Edis = mean(iPictureEdis);

% 绘制条形图
figure
bar(iPictureEdis)
title('Edis')
end