decomposeProjectionMatrix gives unexpected result

Question

I have the following projection matrix P:

-375   0    2000  262500
-375  2000    0   262500
 -1    0      0    700

This projection matrix projects 3D points in mm on a detector in px (with 1px equals to 0.5mm) and is built from the intrinsic matrix K and the extrinsic matrix [R|t] (where R is a rotation matrix and t a translation vector) according the relation P = K [R|t].

     2000    0    375             0   0   1              0
K =    0   2000   375       R =   0   1   0        t =   0
       0     0     1             -1   0   0             700

For some reasons I need to decompose P back into these matrices. When I use decomposeProjectionMatrix I get this as a rotation matrix:

 0   0   0
 0   0   0
-1   0   0

Which doesn't look like a rotation matrix to me.

Moreover when I build back the projection matrix from the Open CV decomposition I get this matrix:

-375   0      0   262500
-375   0      0   262500
 -1    0      0    700

Looks similar but it is not the same.

I'm wondering if I'm doing something wrong or if I'm unlucky and that was one of the rare cases where this function fails.

Note that I did the decomposition by myself and I get coherent results but I would rather use Open CV functions as much as possible.

Answer 1

The problem seems to be in the RQ decomposition used by decomposeProjectionMatrix. Even though the first square of the matrix P is non singular, the RQDecomp3x3 function gives incorrect results:

    0   0  375              0  0  0
R = 0   0  375         Q =  0  0  0
    0   0   1              -1  0  0

So a work around is to use a homemade function (here written in Python) based on the section 2.2 of Peter Sturm's lectures:

def decomposeP(P):
    import numpy as np
    M = P[0:3,0:3]
    Q = np.eye(3)[::-1]
    P_b = Q @ M @ M.T @ Q
    K_h = Q @ np.linalg.cholesky(P_b) @ Q
    K = K_h / K_h[2,2]
    A = np.linalg.inv(K) @ M
    l = (1/np.linalg.det(A)) ** (1/3)
    R = l * A
    t = l * np.linalg.inv(K) @ P[0:3,3]
    return K, R, t

I use the anti-identity matrix Q to build the non conventional Cholesky decomposition U U* where U is upper triangular. This method differs slightly from the Peter Sturm's one as we use the relation P = K[R|t] while in Peter Sturm's lectures the relation used is P = K[R|-Rt].

A C++ implementation using only Open CV is trickier as they don't really expose a function for Cholesky decompostion:

void chol(cv::Mat const& S, cv::Mat& L)
{
    L = cv::Mat::zeros(S.rows, S.rows, cv::DataType<double>::type);
    for (int i = 0; i < S.rows; ++i) {
        for (int j = 0; j <= i ; ++j) {
            double sum = 0;
            for (int k = 0; k < j; ++k)
                sum += L.at<double>(i,k) * L.at<double>(j,k);
            L.at<double>(i,j) = (i == j) ? sqrt(S.at<double>(i,i) - sum) : (S.at<double>(i,j) - sum) / L.at<double>(j,j);
        }
    }
}

void decomposeP(cv::Mat const& P, cv::Mat& K, cv::Mat& R, cv::Mat& t)
{
    cv::Mat M(3, 3, cv::DataType<double>::type);
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            M.at<double>(i, j) = P.at<double>(i ,j);
    cv::Mat Q = cv::Mat::zeros(3, 3, cv::DataType<double>::type);
    Q.at<double>(0, 2) = 1.0;
    Q.at<double>(1, 1) = 1.0;
    Q.at<double>(2, 0) = 1.0;
    cv::Mat O = Q * M * M.t() * Q;
    cv::Mat C;
    chol(O, C);
    cv::Mat B = Q * C * Q;
    K = B / B.at<double>(2,2);
    cv::Mat A = K.inv() * M;
    double l = std::pow((1 / cv::determinant(A)), 1/3);
    R = l * A;
    cv::Mat p(3, 1, cv::DataType<double>::type);
    for (int i = 0; i < 3; ++i)
        p.at<double>(i, 0) = P.at<double>(i ,3);
    t = l * K.inv() * p;
}