Multiple View Triangulation method used by COLMAP

Question

I'm looking at how COLMAP does multi-view triangulation. I can't work out what this function is doing. I can't find any formulas which look similar. The input "proj_matricies" come from pose data, the "points" are from the mono camera measurements. Seems like it is reducing an error in the projection matrices themselves? Is it trying to find some sort of overall projection error for all the observations? Any point in the right direction would be greatly appreciated!!

Eigen::Vector3d TriangulateMultiViewPoint(
const std::vector<Eigen::Matrix3x4d>& proj_matrices,
const std::vector<Eigen::Vector2d>& points) {

    CHECK_EQ(proj_matrices.size(), points.size());

    Eigen::Matrix4d A = Eigen::Matrix4d::Zero();

    for (size_t i = 0; i < points.size(); i++) {
        const Eigen::Vector3d point = points[i].homogeneous().normalized();
        const Eigen::Matrix3x4d term =
        proj_matrices[i] - point * point.transpose() * proj_matrices[i];
        A += term.transpose() * term;
  }

  Eigen::SelfAdjointEigenSolver<Eigen::Matrix4d> eigen_solver(A);

  return eigen_solver.eigenvectors().col(0).hnormalized();
}

Answer 1

This method is minimizing the reprojection error of the triangulated point. How to see this: Let $P_i$ be the $3 \times 4$ projection matrices, and let $x_i \in \mathbb{R}^3$ be the (homogenized coordinates of the) observations. This method is computing the matrix

$$ A = \sum_{i=1}^n A_i^T A_i, $$ where $$ A_i = P_i - x_i x_i^T P_i; $$ and then it returns the eigenvector $X^* \in \mathbb{R}^4$ corresponding to the smallest eigenvalue of $A$.

It is a well-known fact in linear algebra that finding this eigenvector is equivalent to solving $$ X^* = \operatorname{argmin}_{X \in \mathbb{R}^k} X^T A X $$ where in our case $k=4$ because $A$ is $4 \times 4$. Hence \begin{align} X^* &= \operatorname{argmin}_{X \in \mathbb{R}^4} X^T A X \\ &= \operatorname{argmin}_{X \in \mathbb{R}^4} \sum_{i=1}^n (A_i X)^T (A_i X) \\ & = \operatorname{argmin}_{X \in \mathbb{R}^4} \sum_{i=1}^n \| P_i X - x_i x_i^T P_i X \|_2^2. \end{align} Hence it's a least squares minimization, where each residual is $\| P_i X - x_i x_i^T P_i X\|_2^2$. Now, $P_i X$ is the projection of the triangulated point $X$ into image coordinates via the camera $P_i$, and the map $x \mapsto x_i x_i^T x$ is the (linear) projection of $x$ onto the linear subspace spanned by $x_i$. So, the $i$-th residual is the distance in image space between $P_i X$ and $x_i$, which can be interpreted as a reprojection error.