相机成像的过程就是世界坐标系向像素坐标系转换的过程，即：世界坐标系(3d)->相机坐标系(3d)->像平面坐标系(2d)->像素坐标系(2d)，经过这样一级一级的转换之后，物体在空间中的坐标即转换为在图像中的像素坐标。

  • 世界坐标系(Xw、Yw、Zw)
  • 相机坐标系(Xc、Yc、Zc)
  • 像平面坐标系(x、y)
  • 像素坐标系(u、v)

对于世界坐标系到相机坐标系的转换是刚体变换，是旋转动作和平移动作的结果，如下：

\[ \begin{bmatrix} X_c \\ Y_c \\ Z_c \\ 1 \end{bmatrix} = \begin{bmatrix} R & t \\ 0 & 1 \end{bmatrix} \begin{bmatrix} X_w \\ Y_w \\ Z_w \\ 1 \end{bmatrix} \]

投影模型(Projection Model)

小孔成像模型（Pinhole Model）

相机可以抽象为最简单的形式：一个小孔和一个成像平面，小孔位于成像平面和真实的三维场景之间，任何来自真实世界的光只有通过小孔才能到达成像平面。因此，在成像平面和通过小孔看到的真实三维场景存在着一种对应关系，也就是图像中的二维像点和真实三维世界的三维点存在某种变换关系。找到了这种变换关系，就可以利用图像中的二维点信息来恢复场景的三维信息。

相机的成像过程实际是将三维空间的点P=(X,Y,Z)变换到成像平面的过程，这个过程也可以被称为射影变换.

下图是小孔成像的模型，为了简化模型，将成像平面放在了小孔的前面，并且成的像也是正立的。

小孔成像实际就是将相机坐标系中的三维点变换到成像平面中的图像坐标系中的二维点。

假设，三维空间中点P, 其在相机坐标系下坐标为\(P_c = [X, Y, Z]^T\); 其像点p，在图像坐标系中坐标为\(p = [x, y]^T\).

由相似三角形原理，得到

\[ \frac{Z}{f} = \frac{X}{x} = \frac{Y}{y} \]
\[ x = f \frac{x}{Z}, \, y = f \frac{Y}{Z}, \, z = f \]

内参

相机的内参数由下面的两部分组成：

  1. 射影变换本身的参数，相机的焦点到成像平面的距离，也就是焦距f。

  2. 从成像平面坐标系到像素坐标系的变换。上面推导中使用的像点坐标p=(x,y)是成像平面坐标系下，以成像平面的中心为原点。而实际像素点的表示方法是以像素来描述，坐标原点通常是图像的左上角，X轴沿着水平方向向左，Y轴竖直向下。像素是一个矩形块，这里假设其在水平和竖直方向的长度分别为：α和β。所以像素坐标和成像平面坐标之间，相差了一个缩放和原点的平移。 假设像素坐标的水平方向的轴为μ，竖直方向的轴为ν，那么将一个成像平面的坐标(x,y)在水平方向上缩放α倍，在竖直方向上缩放β倍，同时平移(cx,cy)，就可以得到像素坐标系的坐标(μ,ν)，其公式如下：

    \[ u = \alpha \cdot x + c_x , \, v = \beta \cdot y + c_y \]
    \[ f_x = \alpha \cdot f, \, f_y = \beta \cdot f \]

带入上式：

\[ \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \frac{1}{Z} \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} \]

通过上面的的推导，就得到了相机的内参数矩阵（Camera Intrinsics）K，

\[ K = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \]

K有4个未知数和相机的构造相关，fx,fy和相机的焦距，像素的大小有关；cx,cy是平移的距离，和相机成像平面的大小有关。 求解相机内参数的过程被称为标定，在SLAM中可以假定相机的内参是已知的，而在三维重建中内参数则是未知的，需要手动的标定（比如使用标定板）。

\[ \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \begin{bmatrix} \frac{1}{dx} & 0 & c_x \\ 0 & \frac{1}{dy} & c_y \\ 0 & & 1 \end{bmatrix} \begin{bmatrix} f & 0 & 0 & 0 \\ 0 & f & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} R & t \\ 0 & 1 \end{bmatrix} \begin{bmatrix} X_w \\ Y_w \\ Z_w \\ 1 \end{bmatrix} \]

小孔成像模型中，光线通过小孔(假设小孔极小，以至于一次只能通过一束光线)时，在像平面上会呈现出物体清晰的影像。 在实际的情况下，小孔总是有一定物理尺寸的，因此像平面的一点能够接收到锥形区域内的所有光线，所以理想化的小孔成像模型描述真实的相机模型并不准确，尤其是真实的相机都是配备镜头的，使得真实的模型结构更加的复杂。尽管如此，由于小孔透视投影模型的简洁性，我们经常从这一模型开始分析。 从小孔透视投影模型可以很容易分析出以下两条显而易见的公理：

  • 成像的大小取决于物体本身的大小和物体距离小孔的距离。这也符合我们的直观感受，越近的物体看起来越大。
  • 物理空间中同一平面的两条平行线在相平面中汇集于一点。这也符合我们的视觉感受，沿着马路放眼望去，远处的路面比近处的路面要窄。

绝大多数真实的相机都是安装有镜头的，镜头主要有两个作用：收集更多的光线和聚焦，因为如在理想的成像模型小孔透视投影模型中的分析，如果小孔非常小，以至于只有一条光线可以通过，则像平面上的每个点都对应唯一的一条穿过小孔的光线，实际场景中小孔总是有一定大小的，因此像平面上的一个像点都对应一个锥形区域内的光线，小孔的物理尺寸越大，成像就越明亮，同时也就越模糊，小孔的物理尺寸越小，透过小孔的光线就越少，成像也就越暗越清晰，如果小孔的尺寸足够小，还会造成光的干涉等物理现象。为了避免以上问题，使得成像即清晰又明亮，就需要镜头捕获更多的光线。

如果不考虑光的光的干涉、衍射等其它光的物理特性，镜头的行为符合几何光学定律，即

  • 光在均匀的介质中沿直线传播;
  • 光在介质表面会发生反射行为，并且入射光线、反射光线和介质的法线的夹角相等;
  • 光从一种介质进入另一种介质，会发生折射现象，折射行为满足折射定律。

畸变模型（Distortion Model)

相机畸变的分类：

  1. 枕形畸变：又称鞍形形变，视野中边缘区域的放大率远大于光轴中心区域的放大率，常用在远摄镜头中（下图左）；
  2. 桶形畸变，与枕形畸变相反，视野中光轴中心区域的放大率远大于边缘区域的放大率，常出现在广角镜头和鱼眼镜头中（下图中）；
  3. 线性畸变：光轴与相机所拍摄的诸如建筑物类的物体的垂平面不正交，则原本应该平行的远端一侧和近端一侧，以不相同的角度汇聚产生畸变。这种畸变本质上是一种透视变换，即在某一特定角度，任何镜头都会产生相似的畸变（下图右）。

Radial Tangential

  1. 径向畸变:是由于透镜形状的制造工艺导致，且越向透镜边缘移动径向畸变越严重，实际情况中我们常用r=0处的泰勒级数展开的前几项来近似描述径向畸变，矫正径向畸变前后的坐标关系为：

    \[ x_{rad} = x(1 + k_1r^2 + k_2r^4 +k_3r^6) \]
    \[ y_{rad} = y(1 + k_1r^2 + k_2r^4 +k_3r^6) \]

  2. 切向畸变:是由于透镜和CMOS或者CCD的安装位置误差导致，切向畸变需要两个额外的畸变参数来描述，矫正前后的坐标关系为：

    \[ x_{tan} = x + 2p_1xy + p_2(r^2 + 2x^2) \]
    \[ y_{tan} = y + p_1(r^2 + 2y^2) + 2p_2xy \]

综合以上两种畸变，得到相机的畸变模型（纠正后的图像像素坐标系的坐标）：

\[ r^2 = x^2 + y^2 \]
\[ x_{distorted} = x(1 + k_1r^2 + k_2r^4 +k_3r^6) + 2p_1xy + p_2(r^2 + 2x^2) \]
\[ y_{distorted} = y(1 + k_1r^2 + k_2r^4 +k_3r^6) + p_1(r^2 + 2y^2) + 2p_2xy \]

其中x,y 是是去畸变后的图像坐标，它是归一化的坐标，以图像中心为原点。r为半径。\(x_{distorted} y_{distorted}\)是具有畸变得图像坐标。

Example

/**
******************************************************************************
* @file
* @author  maky <chengwei920412@outlook.com>
* @version
* @date    2018-12-23 16:59:58
* @brief
******************************************************************************
* @attention
*
*
******************************************************************************
*/
#include <iostream>
#include <opencv2/opencv.hpp>
#include <Eigen/Eigen>


int main(int argc, char * argv[])
{
    cv::Mat source = cv::imread("1545031406996771344.png");
    cv::cvtColor(source, source, cv::ColorConversionCodes::COLOR_BGR2GRAY);

    double fx = 4.329118910696126e+02;
    double s = 0;
    double cx = 3.784662176701153e+02;
    double fy = 4.273013388506971e+02;
    double cy = 2.234624281814635e+02;

    double k1 = -0.326206456728099;
    double k2 = 0.079226451508043;
    double k3 = 0;
    double p1 = 0.003256859297789;
    double p2 = -0.002936791997006;

    cv::Mat intrinsic_matrix = cv::Mat::eye(3, 3, CV_64F);
    /* matlab IntrinsicMatrix
    * --            --
    * |  fx  0   0  |
    * |  s   fy  0  |
    * |  cx  cy  1  |
    * --            --
    */
    /*
    * --            --
    * |  fx  s   cx  |
    * |  0   fy  cy  |
    * |  0   0   1   |
    * --            --
    */
    intrinsic_matrix.at<double>(0, 0)/* fx */ = fx;
    intrinsic_matrix.at<double>(0, 1)/* s  */ = s;
    intrinsic_matrix.at<double>(0, 2)/* cx */ = cx;
    intrinsic_matrix.at<double>(1, 1)/* fy */ = fy;
    intrinsic_matrix.at<double>(1, 2)/* cy */ = cy;

    cv::Mat distortion_coeffs = cv::Mat::zeros(5, 1, CV_64F);
    /* matlab RadialDistortion
    * --       --
    * | k1  k2 |
    * --       --
    */
    /* matlab TangentialDistortion
    * --       --
    * | p1  p2 |
    * --       --
    */
    /*
    * --                    --
    * |  k1, k2, p1, p2, k3  |
    * --                    --
    */
    distortion_coeffs.at<double>(0, 0) = k1;
    distortion_coeffs.at<double>(1, 0) = k2;
    distortion_coeffs.at<double>(2, 0) = p1;
    distortion_coeffs.at<double>(3, 0) = p2;
    distortion_coeffs.at<double>(4, 0) = k3;

    cv::Mat map_x, map_y;
    cv::initUndistortRectifyMap(intrinsic_matrix, distortion_coeffs, cv::Mat(), /*cv::getOptimalNewCameraMatrix(intrinsic_matrix, distortion_coeffs, cv::Size(752, 480), 1, cv::Size(752, 480), 0)*/cv::Mat() ,cv::Size(752, 480), CV_16SC2, map_x, map_y);

    {
        cv::Mat output;
        cv::remap(source, output, map_x, map_y, cv::INTER_LINEAR);
        cv::imshow("rectify map", output);
    }
    {
        cv::Mat output;
        cv::undistort(source, output, intrinsic_matrix, distortion_coeffs);
        cv::imshow("undistort", output);
    }

    {
        cv::Mat output = cv::Mat::zeros(source.rows, source.cols, CV_8UC1);
        for (int v = 0; v < source.rows; v++) {
            for (int u = 0; u < source.cols; u++) {
                auto pixel = source.at<uchar>(v, u);

                double x = (u - cx) / fx;
                double y = (v - cy) / fy;

                double r = sqrt(x*x + y*y);
                double x_distorted = x*(1 + k1*r*r + k2*r*r*r*r) + 2 * p1*x*y + p2*(r*r + 2 * x*x);
                double y_distorted = y*(1 + k1*r*r + k2*r*r*r*r) + p1*(r*r + 2 * y*y) + 2 * p2*x*y;

                double u_distorted = fx * x_distorted + cx;
                double v_distorted = fy * y_distorted + cy;
                if (u_distorted >= 0 && v_distorted >= 0 && u_distorted < source.cols && v_distorted < source.rows) {
                    output.at<uchar>(v, u) = source.at<uchar>((int)v_distorted, (int)u_distorted);
                }
            }
        }
        cv::imshow("impl", output);
    }


    while (true) {
        cv::imshow("source", source);
        cv::waitKey(1);
    }

    return 0;
}

References

J. Heikilla and O. Silven. A four-step camera calibration procedure with implicit image correction. In Computer Vision and Pattern Recognition, 1997.
https://github.com/ethz-asl/kalibr
https://github.com/hengli/camodocal
Single View Point Omnidirectional Camera Calibration from Planar Grids
https://www.uio.no/studier/emner/matnat/its/UNIK4690/v16/forelesninger/lecture_1_4-the-perspective-camera-model.pdf