Get scale, translation and rotation from CATransform3D

Given a CATransform3D transform, I want to extract the scale, translation and rotation as separate transforms. From some digging, I was able to accomplish this for CGAffineTransform in Swift, like so:

extension CGAffineTransform {
    var scaleDelta:CGAffineTransform {
        let xScale = sqrt(a * a + c * c)
        let yScale = sqrt(b * b + d * d)
        return CGAffineTransform(scaleX: xScale, y: yScale)
    }
    var rotationDelta:CGAffineTransform {
        let rotation = CGFloat(atan2f(Float(b), Float(a)))
        return CGAffineTransform(rotationAngle: rotation)
    }
    var translationDelta:CGAffineTransform {
        return CGAffineTransform(translationX: tx, y: ty)
    }
}

How would one do something similar for CATransform3D using math? (I am looking for a solution that doesn't use keypaths.)

(implementation or math-only answers at your discretion)

If you're starting from a proper affine matrix that can be decomposed correctly (if not unambiguously) into a sequence of scale, rotate, translate, this method will perform the decomposition into a tuple of vectors representing the translation, rotation (Euler angles), and scale components:

extension CATransform3D {
    func decomposeTRS() -> (float3, float3, float3) {
        let m0 = float3(Float(self.m11), Float(self.m12), Float(self.m13))
        let m1 = float3(Float(self.m21), Float(self.m22), Float(self.m23))
        let m2 = float3(Float(self.m31), Float(self.m32), Float(self.m33))
        let m3 = float3(Float(self.m41), Float(self.m42), Float(self.m43))

        let t = m3

        let sx = length(m0)
        let sy = length(m1)
        let sz = length(m2)
        let s = float3(sx, sy, sz)

        let rx = m0 / sx
        let ry = m1 / sy
        let rz = m2 / sz

        let pitch = atan2(ry.z, rz.z)
        let yaw = atan2(-rx.z, hypot(ry.z, rz.z))
        let roll = atan2(rx.y, rx.x)
        let r = float3(pitch, yaw, roll)

        return (t, r, s)
    }
}

To show that this routine correctly extracts the various components, construct a transform and ensure that it decomposes as expected:

let rotationX = CATransform3DMakeRotation(.pi / 2, 1, 0, 0)
let rotationY = CATransform3DMakeRotation(.pi / 3, 0, 1, 0)
let rotationZ = CATransform3DMakeRotation(.pi / 4, 0, 0, 1)
let translation = CATransform3DMakeTranslation(1, 2, 3)
let scale = CATransform3DMakeScale(0.1, 0.2, 0.3)
let transform = CATransform3DConcat(CATransform3DConcat(CATransform3DConcat(CATransform3DConcat(scale, rotationX), rotationY), rotationZ), translation)
let (T, R, S) = transform.decomposeTRS()
print("\(T), \(R), \(S))")

This produces:

float3(1.0, 2.0, 3.0), float3(1.5708, 1.0472, 0.785398), float3(0.1, 0.2, 0.3))

Note that this decomposition assumes an Euler multiplication order of XYZ, which is only one of several possible orderings.

Caveat: There are certainly values for which this method is not numerically stable. I haven't tested it extensively enough to know where these pitfalls lie, so caveat emptor.