Results from decomposeProjectionMatrix() when using different scales for P are not consistent #23733
Description
System Information
OpenCV python version: 4.7.0.68
Operating System / Platform: Windows 11
Python version: 3.8.15
Detailed description
When using decomposeProjectionMatrix(), the returned rotation matrix is not consistent for different scales of input matrix P.
Specifically, for small scales of the matrix P, the rotation matrix returned has significant error and does not fulfill the expected orthonormal conditions (i.e. orthogonality of column vectors, RTR = I). Errors are subsequently observed when recomposing the camera matrix P.
In case helpful, I've copied code below showing correct results are obtained when using either RQ decomposition from Numpy (via QR) or Scipy (as provided). In these cases the residuals in RTR - I are of the order of machine precision for all scales of P.
Steps to reproduce
import cv2
import numpy as np
from scipy import linalg
np.set_printoptions(suppress=False)
''' Decomposing camera matrix and analysis of resulting R matrix properties and recomposition '''
print("**R matrix**")
print("*OpenCV*")
# case 1
P = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0]], dtype=np.float64)
print(f"P (case 1): \n{P}")
_, R, _, _, _, _, _ = cv2.decomposeProjectionMatrix(P)
print(f"OpenCV, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _, _, _, _, _ = cv2.decomposeProjectionMatrix(P)
print(f"OpenCV scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
# case 2
P = np.array([[52, -7, 4, 12],
[-6, 49, 12, 8],
[4, 17, 1, 0]], dtype=np.float64)
print(f"P (case 2): \n{P}")
_, R, _, _, _, _, _ = cv2.decomposeProjectionMatrix(P)
print(f"OpenCV, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _, _, _, _, _ = cv2.decomposeProjectionMatrix(P)
print(f"OpenCV scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}\n")
print("*Numpy*")
def rq(M):
Q, R = np.linalg.qr(np.flipud(M).transpose())
R = np.flipud(R.transpose())
R = np.fliplr(R)
Q = Q.transpose()
Q = np.flipud(Q)
R = R * np.linalg.det(Q)
Q = Q * np.linalg.det(Q)
return R, Q
def decompose_proj_numpy(P):
K, R = rq(P[:3, :3])
t = -np.linalg.inv(P[:3, :3]) @ P[:3, 3]
return K, R, t
# case 1
P = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0]], dtype=np.float64)
print(f"P (case 1): \n{P}")
_, R, _ = decompose_proj_numpy(P)
print(f"Numpy, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _ = decompose_proj_numpy(P)
print(f"Numpy scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
# case 2
P = np.array([[52, -7, 4, 12],
[-6, 49, 12, 8],
[4, 17, 1, 0]], dtype=np.float64)
print(f"P (case 2): \n{P}")
_, R, _ = decompose_proj_numpy(P)
print(f"Numpy, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _ = decompose_proj_numpy(P)
print(f"Numpy scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}\n")
print("*Scipy*")
def decompose_proj_scipy(P):
K, R = linalg.rq(P[:3, :3])
t = -np.linalg.inv(P[:3, :3]) @ P[:3, 3]
return K, R, t
# case 1
P = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0]], dtype=np.float64)
print(f"P (case 1): \n{P}")
_, R, _ = decompose_proj_scipy(P)
print(f"Scipy, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _ = decompose_proj_scipy(P)
print(f"Scipy scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
# case 2
P = np.array([[52, -7, 4, 12],
[-6, 49, 12, 8],
[4, 17, 1, 0]], dtype=np.float64)
print(f"P (case 2): \n{P}")
_, R, _ = decompose_proj_scipy(P)
print(f"Scipy, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
P *= 1e-6
_, R, _ = decompose_proj_scipy(P)
print(f"Scipy scaled 1e-6, R.T @ R - I: \n{R.T @ R - np.eye(3)}")
print("\n**Camera matrix recomposition**")
def compose_proj(K, R, t):
P = K @ np.hstack([R, -R @ t[np.newaxis].T])
return P
P = np.array([[52, -7, 4, 12],
[-6, 49, 12, 8],
[4, 17, 1, 0]], dtype=np.float64)
P *= 1e-6
print(f"P: \n{P}\n")
print("*OpenCV*")
K, R, t, _, _, _, _ = cv2.decomposeProjectionMatrix(P)
P_recompose = compose_proj(K, R, (t[:3]/t[3]).flatten())
print(f"OpenCV, P' - P:\n {P_recompose - P}\n")
print("*Numpy*")
K, R, t = decompose_proj_numpy(P)
P_recompose = compose_proj(K, R, t)
print(f"Numpy, P' - P:\n {P_recompose - P}\n")
print("*Scipy*")
K, R, t = decompose_proj_scipy(P)
P_recompose = compose_proj(K, R, t)
print(f"Scipy, P' - P:\n {P_recompose - P}")Issue submission checklist
- I report the issue, it's not a question
- I checked the problem with documentation, FAQ, open issues, forum.opencv.org, Stack Overflow, etc and have not found any solution
- I updated to the latest OpenCV version and the issue is still there
- There is reproducer code and related data files (videos, images, onnx, etc)