图像畸变标定与矫正

合集 - 图像处理(3)

1.图像畸变标定与矫正2023-06-21

2.张正友标定法2023-06-25 3.基础 | 详解3D结构光如何标定2023-08-02

畸变类型

相机畸变主要分为两种：

径向畸变（枕形、桶形）：光线在远离透镜中心的地方比靠近透镜中心的地方更加弯曲。
切向畸变：透镜不完全平行于图像平面，即sensor在安装时与镜头之间的角度不准。

1.摄像机成像原理简述

成像的过程实质上是几个坐标系的转换。首先空间中的一点由世界坐标系转换到摄像机坐标系，然后再将其投影到成像平面 ( 图像物理坐标系 ) ，最后再将成像平面上的数据转换到图像平面 ( 图像像素坐标系 ) 。

详细的可以参考我之前的博客[图像]摄像机标定(1) 标定中的四个坐标系

图像像素坐标系 (uOv坐标系) 下的无畸变坐标 (U, V)，经过经向畸变和切向畸变后落在了uOv坐标系的 (Ud, Vd) 上。即就是说，真实图像 imgR 与畸变图像 imgD 之间的关系为: imgR(U, V) = imgD(Ud, Vd) 。

2.成像畸变

2.1. 畸变数学模型

径向畸变

成像仪光轴中心的畸变为 0 ，沿着镜头半径方向越靠近边缘，畸变越严重。

畸变的数学模型可以用主点（principle point）周围的泰勒级数展开式的前几项来进行描述。

通常使用前两项和，即 k1 和 k2
对于畸变很大的镜头，如鱼眼镜头，需要增加第三项 k3 来进行描述。

其中（u‘，v'）是畸变后的归一化图像坐标，（u，v）是无畸变的归一化图像坐标。

　　非常值得注意的是：畸变是发生在世界点通过透镜投影到相机坐标系这个过程中的，所以去畸变操作是在相机坐标系下进行的，由于我们经常把世界点投影到归一化相机平面，所以我们在归一化相机平面上对图像去畸变！

2.2. 公式推导

公式推导:

3.畸变校正

3.1. 理论推导

我们已知的是畸变后的图像，要得到没有畸变的图像就要通过畸变模型推导其映射关系。真实图像 imgR 与畸变图像 imgD 之间的关系为: imgR(U, V) = imgD(Ud, Vd) 。通过这个关系，找出所有的 imgR(U, V) 。(U, V) 映射到 (Ud, Vd) 中的 (Ud, Vd) 往往不是整数 (U和V是整数，因为它是我们要组成图像的像素坐标位置，以这正常图像的坐标位置去求在畸变图像中的坐标位置，取出对应的像素值，这也是正常图像的像素值)。但是畸变的像素往往不是整数，所以需要通过插值来进行求解，详细见我之前的博客 [图像]图像缩放算法-双线性内插法。

3.2. 代码实现

clear;
    clc;
    A =[5.9418398977142772e+002 0 3.1950000000000000e+002;
        0    5.941839897714e+002 2.3950000000000000e+002;
        0 0 1];%相机内参
    D = [6.7442695842244271e-002 2.4180872220967919e-001 0 0 -3.3970575589699975e-001];%畸变向量
    fx = A(1,1);
    fy = A(2,2);
    cx = A(1,3);
    cy = A(2,3);
    k1 = D(1);
    k2 = D(2);
    k3 = D(5);
    p1 = D(3);
    p2 = D(4);
  
    K = A;
    I_d = imread('img_d.png');
    I_d = rgb2gray(I_d);
    I_d = im2double(I_d);
  
    I_r = zeros(size(I_d));
  
    %图像坐标系和矩阵的表示是相反的
    %[row,col] = find(X)，坐标按照列的顺序排列，这样好和reshape()匹配出响应的图像矩阵
    [v u] = find(~isnan(I_r));
  
    % XYZc 摄像机坐标系的值，但是已经归一化了，因为没有乘比例因子
    %公式 s[u v 1]' = A*[Xc Yc Zc]' ，其中s为比例因子，不加比例因子，Zc就为1，所以此时的Xc相对于( Xc/Zc )
    XYZc= inv(A)*[u v ones(length(u),1)]';
    % 此时的x和y是没有畸变的
    r2 = XYZc(1,:).^2+XYZc(2,:).^2;
    x = XYZc(1,:);
    y = XYZc(2,:);
    % x_1和y_1进行畸变的
    x_1 = x.*(1+k1*r2 + k2*r2.^2) + 2*p1.*x.*y + p2*(r2 + 2*x.^2);
    y_1 = y.*(1+k1*r2 + k2*r2.^2) + 2*p2.*x.*y + p1*(r2 + 2*y.^2);
    % (u, v) 对应的畸变坐标 (u_d, v_d)
    u_d = reshape(fx*x_1 + cx,size(I_r));
    v_d = reshape(fy*y_1 + cy,size(I_r));
    % 线性插值出非畸变的图像
    I_r = interp2(I_d, u_d, v_d);
    %对比图像
    subplot(121);     
        imagesc(I_d);
        title('畸变原图像');
    subplot(122);
        imagesc(I_r);
        title('校正后图像');