视图重建

如果空间点在两个或多个视图中出现，如何恢复景物的具体外观呢？

两视图重建

假定某物理点出现在两个视角的图片中，基于图像匹配的算法，我们确定了这对点集： 。

三维重建的任务是确定摄像机矩阵 以及三维坐标点X。如果我们已经得到了P和P’，计算方法下一篇会做详细的阐述。

根据摄像机模型：

通过叉积构建齐次线性方程组 , 展开:

叉乘得到

所以：

考虑两个对应点(视图),得到:

这里是从每个视图中取两个方程,得到4个方程,这是一个冗余方程组,可采用最小二乘法求解即可。

多点反投影（重建）

相机的RT参数：

struct camera_params_t
{
    double R[9];     /* Rotation */
    double T[3];     /* Translation */
    double t[3];     /*-RT*/
};

反投影函数：

v3_t triangulate_Pre( vector<v2_t> &p, vector<camera_params_t>& params, double &error_out)

其中 ，有：

所以:

程序实现:

v3_t triangulate_Pre( vector<v2_t> &p, vector<camera_params_t>& params, double &error_out) 
{//T
    int num_points = p.size();
    int num_eqs = 2 * num_points;
    int num_vars = 3;

    double *A = (double *) malloc(sizeof(double) * num_eqs * num_vars);
    double *b = (double *) malloc(sizeof(double) * num_eqs);
    double *x = (double *) malloc(sizeof(double) * num_vars);

    int i;
    double error;
    v3_t r;

    for (i = 0; i < num_points; i++) {
        int row = 6 * i;
        int brow = 2 * i;

        A[row + 0] = params[i].R[0] - Vx(p[i]) * params[i].R[6];  
        A[row + 1] = params[i].R[1] - Vx(p[i]) * params[i].R[7];  
        A[row + 2] = params[i].R[2] - Vx(p[i]) * params[i].R[8];

        A[row + 3] = params[i].R[3] - Vy(p[i]) * params[i].R[6];  
        A[row + 4] = params[i].R[4] - Vy(p[i]) * params[i].R[7];  
        A[row + 5] = params[i].R[5] - Vy(p[i]) * params[i].R[8];

        b[brow + 0] = params[i].T[2] * Vx(p[i]) - params[i].T[0];
        b[brow + 1] = params[i].T[2] * Vy(p[i]) - params[i].T[1];
    }
    
    /* Find the least squares result */
    dgelsy_driver(A, b, x, num_eqs, num_vars, 1);
    error = 0.0;
    for (i = 0; i < num_points; i++) {
        double dx, dy;
        double pp[3];

        /* Compute projection error */
        matrix_product331(params[i].R , x, pp);
        pp[0] += params[i].T[0];
        pp[1] += params[i].T[1];
        pp[2] += params[i].T[2];

        dx = pp[0] / pp[2] - Vx(p[i]);
        dy = pp[1] / pp[2] - Vy(p[i]);

        error += dx * dx + dy * dy;
    }
    error = sqrt(error / num_points);
    printf("1:triangulate error is %f 
", error);
    /* Run a non-linear optimization to refine the result */
    global_num_points = num_points;
    global_p = p;
    global_params = params;
    lmdif_driver(triangulate_n_residual, num_eqs, num_vars, x, 1.0e-5);

    error = 0.0;
    for (i = 0; i < num_points; i++) {
        double dx, dy;
        double pp[3];

        /* Compute projection error */
        matrix_product331(params[i].R , x, pp);
        pp[0] += params[i].T[0];
        pp[1] += params[i].T[1];
        pp[2] += params[i].T[2];

        dx = pp[0] / pp[2] - Vx(p[i]);
        dy = pp[1] / pp[2] - Vy(p[i]);

        error += dx * dx + dy * dy;
    }
    error = sqrt(error / num_points);
    printf("2:triangulate error is %f 
", error);

    if (error_out != NULL) {
        error_out = error;
    }
    r = v3_new(x[0], x[1], x[2]);
    free(A);
    free(b);
    free(x);
    return r;
}

求优函数参考:https://github.com/TheFrenchLeaf/Bundler/blob/master/lib/matrix/matrix.h

dgelsy_driver函数: Driver for the lapack function dgelss, which finds x to minimize norm(b - A * x)

lmdif_driver函数:Driver for the minpack function lmdif, which uses Levenberg-Marquardt for non-linear least squares minimization (可参见:https://www.math.utah.edu/software/minpack/minpack/lmdif.html)

具体函数的不同见下一篇：最小二乘法

void triangulate_n_residual(const int *m, const int *n, 
    double *x, double *fvec, double *iflag) {
    int i;
    for (i = 0; i < global_num_points; i++) {
        /* Project the point into the view */
        v2_t p = project(i, x);
        fvec[2 * i + 0] = Vx(global_p[i]) - Vx(p);
        fvec[2 * i + 1] = Vy(global_p[i]) - Vy(p);
    }
}
//
v2_t project(int index, double *P) {
    /* Rigid transform */
    double tmp[3], tmp2[3];
    matrix_product331(global_params[index].R, P, tmp);
    matrix_sum(3, 1, 3, 1, tmp, global_params[index].T, tmp2);

    /* Perspective division */
    v2_t result;
    Vx(result) = tmp2[0] / tmp2[2];
    Vy(result) = tmp2[1] / tmp2[2];

    return result;
}