针孔相机的成像原理

之前一直对相机的原理有些模糊，这里记一个笔记，以免遗忘。

设相机坐标系下的点$P(x,y,z)$在相机成像平面上成的像为$P^{{}^{\prime }}(X^{{}^{\prime }},Y^{{}^{\prime }},Z^{{}^{\prime }})$

根据三角形的相似原理，可以得到如下式子：
$$\begin{array}{c}\frac{Z}{f}=\frac{X}{X^{{}^{\prime }}}=\frac{Y}{Y^{{}^{\prime }}}\\ \begin{array}{rlr}{X}^{{}^{\prime }}& =& \frac{f}{Z}X\\ Y^{{}^{\prime }}& =& \frac{f}{Z}Y\end{array}\end{array}\text{\hspace{1em}(1)}$$

注意此时$X^{{}^{\prime }}$的单位是米，但是像素平面上的单位是像素。同时需要注意的一点是，成像平面的原点在成像平面的中点上，一般来说，成像平面的原点在像平面的右上角：

结合以上两点，我们对式(1)进行转化变为：
$$\begin{array}{c}\{\begin{array}{l}u=\alpha X^{{}^{\prime }}+c_x\\ v=\beta Y^{{}^{\prime }}+c_y\end{array}\end{array}$$
其中$u,v$是像素平面上的像素坐标单位是像素，$\alpha ,\beta$是统一尺度的因子，单位为像素/(毫)米。$c_x,c_y$单位为像素，是平移坐标系的偏移量。

此时可以将相机平面下的3D坐标转化为像素平面上的坐标公式写为：
$$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\frac{1}{Z}\left[\begin{array}{ccc}{f}_x,& 0,& c_x\\ 0,& f_y,& c_y\\ 0,& 0,& 1\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$$

其中$f_x$表示$\alpha f$, $f_y$表示$\beta f$。
由此我们便可以总结出将相机坐标系下的3D坐标转化到像素平面上的像素坐标需要先得到相机坐标系下3D点的归一化坐标$[X/Z,Y/Z,1{]}^T$。然后乘上相机的内参矩阵：
$$K=\left[\begin{array}{ccc}{f}_x,& 0,& c_x\\ 0,& f_y,& c_y\\ 0,& 0,& 1\end{array}\right]$$
就可以得到像素平面上的坐标值了。

相机畸变模型和校正方法

相机畸变模型和校正方法

先把公式摆在这里：
$$\begin{array}{rl}& x_{distorted}=x(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\\ & y_{distorted}=y(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\\ & r^2=x^2+y^2\end{array}$$
其中$x,y$是正常相机坐标系下的归一化坐标，$x_{distorted},y_{distorted}$是畸变后的归一化坐标。

以下是根据上文所述的相机模型以及畸变模型编写的图像去畸变矫正代码。

#include <opencv2/opencv.hpp>

using namespace std;
string image_file = "/YOUR_PATH";   // 请确保路径正确

uchar B_interpolation(cv::Mat im, float u, float v){
    if (u < 0 && v < 0 && u >= im.cols && v >= im.rows)
        return 0;
    else{
        uchar a = im.at<uchar>(floor(v), floor(u));
        uchar b = im.at<uchar>(floor(v), ceil(u));
        uchar c = im.at<uchar>(ceil(v), floor(u));
        uchar d = im.at<uchar>(ceil(v), ceil(u));
        float fx1 = (ceil(u)-u)/(ceil(u)-floor(u))*a+(u-floor(u))/(ceil(u)-floor(u))*b;
        float fx2 = (ceil(u)-u)/(ceil(u)-floor(u))*c+(u-floor(u))/(ceil(u)-floor(u))*d;
        return (ceil(v)-v)/(ceil(v)-floor(v))*fx1 + (v-floor(v))/(ceil(v)-floor(v))*fx2;
    }
}

int main(int argc, char **argv) {
 
    double k1 = -0.28340811, k2 = 0.07395907;  // 畸变参数
    double fx = 458.654, fy = 457.296, cx = 367.215, cy = 248.375;

    cv::Mat image = cv::imread(image_file, CV_8UC1);   // 图像是灰度图
    int rows = image.rows, cols = image.cols;
    cv::Mat image_undistort = cv::Mat(rows, cols, CV_8UC1);   // 去畸变以后的图
    cv::imshow("image distorted", image);
    // 计算去畸变后图像的内容
    for (int v = 0; v < rows; v++)
        for (int u = 0; u < cols; u++) {

            float u_distorted = 0, v_distorted = 0;
            float X_undistort = (u - cx)/fx;
            float Y_undistort = (v - cy)/fy;
            float r2 = X_undistort * X_undistort + Y_undistort * Y_undistort;
            float X_distort = X_undistort * (1 + k1*r2 + k2 * r2 * r2);
            float Y_distort = Y_undistort * (1 + k1*r2 + k2 * r2 * r2);
            u_distorted = fx*X_distort + cx;
            v_distorted = fy*Y_distort + cy;

//            if (u_distorted >= 0 && v_distorted >= 0 && u_distorted < cols && v_distorted < rows) {
//                image_undistort.at<uchar>(v, u) = image.at<uchar>((int) v_distorted, (int) u_distorted);
//            } else {
//                image_undistort.at<uchar>(v, u) = 0;
//            }
            image_undistort.at<uchar>(v, u) = B_interpolation(image, u_distorted, v_distorted);//双线性插值法
        }

    cv::imshow("image undistorted", image_undistort);
    cv::waitKey(0);

    return 0;
}

畸变图像：

修正后的图像：

原始图片与相机参数均取自 “计算机视觉life” 公众号。