Converting a rotation matrix to euler angles and back - special case?

Question

I'm disassembling a rotation matrix to Euler angles (Tait-Bryan angles more specifically in the order x-y-z, that is rotation around x axis first) and back to a rotation matrix. I used the transforms3d python library (https://github.com/matthew-brett/transforms3d) and also followed this tutorial www.gregslabaugh.net/publications/euler.pdf Both give the same result.

The problem is that the reassambled rotation matrix doesn't match the one that I started with.

The matrix I'm working with was created by the "decomposeHomographyMat" function from openCV, so I expect it to be a valid rotation matrix. Maybe it is a special case? The matrix is

The three angles are [-1.8710997 , 0.04623301, -0.03679793]. If I convert them back to a rotation matrix I get

where R_23 cannot be a rounding error.

Following the paper above, rotation around the y axis (beta) can be calculated by asin(-R_31). Another valid angle would be pi-asin(-R_31). The angle around the x axis (alpha) can be calculated by atan2(R_32,R_33). I could also get alpha by asin(R_32/cos(beta)) or by acos(R_33/cos(beta)). If I use the latter two equations I only get the same result for alpha if I use beta=pi-arcsin(-R_31), which implies that there is only one valid solution for beta. atan2(R_32,R_33) gives a different result from both.

Anyway something seems to be wrong with my matrix or I cannot figure out why the disassambly doesn't work.

import numpy as np

def rot2eul(R):
    beta = -np.arcsin(R[2,0])
    alpha = np.arctan2(R[2,1]/np.cos(beta),R[2,2]/np.cos(beta))
    gamma = np.arctan2(R[1,0]/np.cos(beta),R[0,0]/np.cos(beta))
    return np.array((alpha, beta, gamma))

def eul2rot(theta) :

    R = np.array([[np.cos(theta[1])*np.cos(theta[2]),       np.sin(theta[0])*np.sin(theta[1])*np.cos(theta[2]) - np.sin(theta[2])*np.cos(theta[0]),      np.sin(theta[1])*np.cos(theta[0])*np.cos(theta[2]) + np.sin(theta[0])*np.sin(theta[2])],
                  [np.sin(theta[2])*np.cos(theta[1]),       np.sin(theta[0])*np.sin(theta[1])*np.sin(theta[2]) + np.cos(theta[0])*np.cos(theta[2]),      np.sin(theta[1])*np.sin(theta[2])*np.cos(theta[0]) - np.sin(theta[0])*np.cos(theta[2])],
                  [-np.sin(theta[1]),                        np.sin(theta[0])*np.cos(theta[1]),                                                           np.cos(theta[0])*np.cos(theta[1])]])

    return R


R = np.array([[ 0.9982552 , -0.03323557, -0.04880523],
       [-0.03675031,  0.29723396, -0.95409716],
       [-0.04621654, -0.95422606, -0.29549393]])

ang = rot2eul(R)
eul2rot(ang)

import transforms3d.euler as eul
ang = eul.mat2euler(R, axes='sxyz')
eul.euler2mat(ang[0], ang[1], ang[2], axes='sxyz')

Answer 1

It turns out the rotation matrix has a negative determinant, which makes it an improper rotation matrix. The openCV function "decomposeHomographyMat" has a bug: https://github.com/opencv/opencv/issues/4978