My version of the code from that site below.
It gives quite good results, but I don't think they are 100% correct.
You will at least have to experiment with the signs for the rotations to get it correct.
I'm also not sure what to do with the determinant, so if you can help out it will be much appreciated.
I have the same question here:
/Community/posts.php?topic=82644#932186
Method Decompose()
_Decompose()
ry = heading '(rotate about y) 'roll
rz = attitude '(rotate about z) 'pitch
rx = bank '(rotate about x) yaw
End Method
Method _Decompose()
' get translation
Tx = a4
Ty = b4
Tz = c4
Local rm:aiMatrix3x3 = New aiMatrix3x3
'transpose to row format instead of colums if needed (not shown here)
rm.a1 = a1 ; rm.a2 = a2 ; rm.a3 = a3
rm.b1 = b1 ; rm.b2 = b2 ; rm.b3 = b3
rm.c1 = c1 ; rm.c2 = c2 ; rm.c3 = c3
Sx = Sqr( rm.a1*rm.a1 + rm.b1*rm.b1 + rm.c1*rm.c1 )
Sy = Sqr( rm.a2*rm.a2 + rm.b2*rm.b2 + rm.c2*rm.c2 )
Sz = Sqr( rm.a3*rm.a3 + rm.b3*rm.b3 + rm.c3*rm.c3 )
Local det:Float = rm.a1*rm.b2*rm.c3 + rm.a2*rm.b3*rm.c1 + rm.a3*rm.b1*rm.c2 - rm.a1*rm.b3*rm.c2 - rm.a2*rm.b1*rm.c3 - rm.a3*rm.b2*rm.c1
If sx Then
rm.a1:/sx
rm.b1:/sx
rm.c1:/sx
EndIf
If sy Then
rm.a2:/sy
rm.b2:/sy
rm.c2:/sy
EndIf
If sz Then
rm.a3:/sz
rm.b3:/sz
rm.c3:/sz
EndIf
' Sx:* Sgn( det)
' Sy:* Sgn( det)
' Sz:* Sgn( det)
If (rm.b1 > 0.998) ' singularity at north pole
heading = ATan2(rm.a3,rm.c3)
attitude = 90 'Pi/2
bank = 0
DebugLog "' singularity at north pole ******************************************************"
Return
EndIf
If (rm.b1 < -0.998)' singularity at south pole
heading = ATan2(rm.a3,rm.c3)
attitude = - 90 '-Pi/2
bank = 0
DebugLog "' singularity at south pole ******************************************************"
Return
EndIf
heading = ATan2(-rm.c1,rm.a1)
bank = ATan2(-rm.b3,rm.b2)
attitude = ASin(rm.b1)
Return