视图重建
如果空间点在两个或多个视图中出现，如何恢复景物的具体外观呢？
两视图重建
假定某物理点出现在两个视角的图片中，基于图像匹配的算法，我们确定了这对点集： 。
三维重建的任务是确定摄像机矩阵 以及三维坐标点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 \n", 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 \n", 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; }
