相机与图像

习题一:图像去畸变

现实生活中的图像总存在畸变。原则上来说,针孔透视相机应该将三维世界中的直线投影成直线,但是当我们使用广角和鱼眼镜头时,由于畸变的原因,直线在图像里看起来是扭曲的。本次作业,你将尝试如何对一张图像去畸变,得到畸变前的图像。

上图为本次练习的测试图片,来自EuRoC数据集。可以明显的看到实际的柱子、箱子的直线边缘在图像中被扭曲成了曲线。这就是由相机畸变造成的。根据我们在课上的介绍,畸变前后的坐标变换为:

$$ \begin {cases} x_{distorted} = x(1+k_1r^2+k_2r^4)+2p_1xy+p_2(r^2+2x^2) \\ y_{distorted} = y(1+k_1r^2+k_2r^4)+p_1(r^2+2y^2)+2p_2xy \end {cases} $$

其中 $x,y$为去畸变后的坐标, $x_ {distorted},y_ {distroted}$为去畸变前的坐标。现给定参数:

$$ k_1 =−0.28340811,k_2=0.07395907,p_1 = 0.00019359,p_2 = 1.76187114e^{−05} $$ 以及相机内参

$$ f_x = 458.654, f_y = 457.296, c_x = 367.215, c_y = 248.375 $$ 请根据undistort_image.cpp文件中内容,完成对该图像的去畸变操作。

答案:

  1. 视觉SALM十四讲中2017年第一版中关于畸变模型的公式是错误的,这一点以这次的为主。畸变公式表示从去畸变后的坐标到去畸变前的坐标变换。
  2. 对其中$r$说明,为极坐标系下的距离,即$r^2=x^2+y^2$
    最终程序的运行结果为:

完整程序链接:https://github.com/XLMaverick/Visual-Localization-Percessing/tree/master/vslam_fourteen_lectures/undistort
代码如下所示:

    #include <opencv2/opencv.hpp>
    #include <string>
    using namespace std;
    string image_file = "../../test.png";   // 请确保路径正确
    int main(int argc, char **argv)
    {
        // 本程序需要你自己实现去畸变部分的代码。尽管我们可以调用OpenCV的去畸变,但自己实现一遍有助于理解。
        // 畸变参数
        double k1 = -0.28340811, k2 = 0.07395907, p1 = 0.00019359, p2 = 1.76187114e-05;
        // 内参
        double fx = 458.654, fy = 457.296, cx = 367.215, cy = 248.375;

        cv::Mat image = cv::imread(image_file,0);   // 图像是灰度图,CV_8UC1
        int rows = image.rows, cols = image.cols;
        cv::Mat image_undistort = cv::Mat(rows, cols, CV_8UC1);   // 去畸变以后的图

        // 计算去畸变后图像的内容
        for (int v = 0; v < rows; v++)
        {
            for (int u = 0; u < cols; u++) 
            {
                
                double u_distorted = 0, v_distorted = 0;
                // TODO 按照公式,计算点(u,v)对应到畸变图像中的坐标(u_distorted, v_distorted) (~6 lines)
                // start your code here
                double x = (u-cx)/fx;
                double y = (v-cy)/fy;

                double r2 = x*x+y*y;
                double r4 = r2*r2;

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

                u_distorted = fx*x_distorted+cx;
                v_distorted = fy*y_distorted+cy;

                // end your code here

                // 赋值 (最近邻插值)
                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;
                }
            }
        }
        // 画图去畸变后图像
        cv::imwrite("../../result.png", image_undistort);
        cv::imshow("image undistorted", image_undistort);
        cv::waitKey();
        return 0;
    }

习题二:双目视差的使用

双目相机的一大好处是可以通过左右目的视差来恢复深度。课程中我们介绍了由视差计算深度的过程。本题,你需要根据视差计算深度,进而生成点云数据。本题的数据来自 Kitti 数据集。
Kitti 中的相机部分使用了一个双目模型。双目采集到左图和右图,然后我们可以通过左右视图恢复出深度。经典双目恢复深度的算法有 BM(Block Matching),SGBM(Semi-Global Matching)等,但本题不探讨立体视觉内容(那是一个大问题)。我们假设双目计算的视差已经给定,请你根据双目模型,画出图像对应的点云,并显示到 Pangolin 中。
本题给定左右图、视差图。双目的参数如下: $$ f_x = 718.856, f_y = 718.856, c_x = 607.1928, c_y = 185.2157 $$ 且双目左右间距(即基线)为: $$ d = 0.573 m $$ 请根据以上的参数,计算相机数据对应的点云,并显示到Pangolin中。

答案:

  1. 注意代码中vector<Vector4d, Eigen::aligned_allocator<Vector4d>> pointcloud,其中Eigen::aligned_allocator<Vector4d>是描述vector中的Allocator type
  2. 相机模型中使用的是坐标,而相应坐标里面的数值并没有考虑。
  3. 相机的模型中我们共用到四种坐标:世界坐标、相机坐标、归一化相机坐标和像素坐标。其中畸变模型是在归一化相机坐标上讨论的。对于具体的问题,需要考虑合适的正确的坐标对应关系。
  4. 视差图转为深度图信息,公式当中的量纲,双目图片两侧小部分没有视差信息。

最终程序运行的结果(结果中是一个三维的图,前后重叠在一起了):
左右图如下所示:

视差图如下所示:
最终生成的点云图如下所示:

完整程序链接:https://github.com/XLMaverick/Visual-Localization-Percessing/tree/master/vslam_fourteen_lectures/disparity
代码如下:

    #include <opencv2/opencv.hpp>
    #include <string>
    #include <Eigen/Core>
    #include <pangolin/pangolin.h>
    #include <unistd.h>

    using namespace std;
    using namespace Eigen;

    // 文件路径,如果不对,请调整
    string left_file = "../../left.png";
    string right_file = "../../right.png";
    string disparity_file = "../../disparity.png";

    // 在panglin中画图,已写好,无需调整
    void showPointCloud(const vector<Vector4d, Eigen::aligned_allocator<Vector4d>> &pointcloud);

    int main(int argc, char **argv) {

        // 内参
        double fx = 718.856, fy = 718.856, cx = 607.1928, cy = 185.2157;
        // 间距
        double b = 0.573;

        // 读取图像
        cv::Mat left = cv::imread(left_file, 0);
        cv::Mat right = cv::imread(right_file, 0);
        cv::Mat disparity = cv::imread(disparity_file, 0); // disparty 为CV_8U,单位为像素

        // 生成点云
        vector<Vector4d, Eigen::aligned_allocator<Vector4d>> pointcloud;

        // TODO 根据双目模型计算点云
        // 如果你的机器慢,请把后面的v++和u++改成v+=2, u+=2
        for (int v = 0; v < left.rows; v++)
            for (int u = 0; u < left.cols; u++) {

                Vector4d point(0, 0, 0, left.at<uchar>(v, u) / 255.0); // 前三维为xyz,第四维为颜色

                // start your code here (~6 lines)
                // 根据双目模型计算 point 的位置
                uchar  d= disparity.at<uchar>(v,u);
                point[2] = (fx*b)/d;
                point[1] = point[2]*(v-cy)/fy;
                point[0] = point[2]*(u-cx)/fx;
                pointcloud.push_back(point);
                // end your code here
            }
        cout<<"pointcloud size: "<<pointcloud.size()<<endl;
        // 画出点云
        showPointCloud(pointcloud);
        return 0;
    }

    void showPointCloud(const vector<Vector4d, Eigen::aligned_allocator<Vector4d>> &pointcloud) {

        if (pointcloud.empty()) {
            cerr << "Point cloud is empty!" << endl;
            return;
        }

        pangolin::CreateWindowAndBind("Point Cloud Viewer", 1024, 768);
        glEnable(GL_DEPTH_TEST);
        glEnable(GL_BLEND);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

        pangolin::OpenGlRenderState s_cam(
                pangolin::ProjectionMatrix(1024, 768, 500, 500, 512, 389, 0.1, 1000),
                pangolin::ModelViewLookAt(0, -0.1, -1.8, 0, 0, 0, 0.0, -1.0, 0.0)
        );

        pangolin::View &d_cam = pangolin::CreateDisplay()
                .SetBounds(0.0, 1.0, pangolin::Attach::Pix(175), 1.0, -1024.0f / 768.0f)
                .SetHandler(new pangolin::Handler3D(s_cam));

        while (pangolin::ShouldQuit() == false) {
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

            d_cam.Activate(s_cam);
            glClearColor(1.0f, 1.0f, 1.0f, 1.0f);

            glPointSize(2);
            glBegin(GL_POINTS);
            for (auto &p: pointcloud) {
                glColor3f(p[3], p[3], p[3]);
                glVertex3d(p[0], p[1], p[2]);
            }
            glEnd();
            pangolin::FinishFrame();
            usleep(5000);   // sleep 5 ms
        }
        return;
    }