def quaternion_from_matrix(matrix, isprecise=False):
M = numpy.array(matrix, dtype=numpy.float64, copy=False)[:<NUM_LIT:4>, :<NUM_LIT:4>]<EOL>if isprecise:<EOL><INDENT>q = numpy.empty((<NUM_LIT:4>, ))<EOL>t = numpy.trace(M)<EOL>if t > M[<NUM_LIT:3>, <NUM_LIT:3>]:<EOL><INDENT>q[<NUM_LIT:0>] = t<EOL>q[<NUM_LIT:3>] = M[<NUM_LIT:1>, <NUM_LIT:0>] - M[<NUM_LIT:0>, <NUM_LIT:1>]<EOL>q[<NUM_LIT:2>] = M[<NUM_LIT:0>, <NUM_LIT:2>] - M[<NUM_LIT:2>, <NUM_LIT:0>]<EOL>q[<NUM_LIT:1>] = M[<NUM_LIT:2>, <NUM_LIT:1>] - M[<NUM_LIT:1>, <NUM_LIT:2>]<EOL><DEDENT>else:<EOL><INDENT>i, j, k = <NUM_LIT:1>, <NUM_LIT:2>, <NUM_LIT:3><EOL>if M[<NUM_LIT:1>, <NUM_LIT:1>] > M[<NUM_LIT:0>, <NUM_LIT:0>]:<EOL><INDENT>i, j, k = <NUM_LIT:2>, <NUM_LIT:3>, <NUM_LIT:1><EOL><DEDENT>if M[<NUM_LIT:2>, <NUM_LIT:2>] > M[i, i]:<EOL><INDENT>i, j, k = <NUM_LIT:3>, <NUM_LIT:1>, <NUM_LIT:2><EOL><DEDENT>t = M[i, i] - (M[j, j] + M[k, k]) + M[<NUM_LIT:3>, <NUM_LIT:3>]<EOL>q[i] = t<EOL>q[j] = M[i, j] + M[j, i]<EOL>q[k] = M[k, i] + M[i, k]<EOL>q[<NUM_LIT:3>] = M[k, j] - M[j, k]<EOL><DEDENT>q *= <NUM_LIT:0.5> / math.sqrt(t * M[<NUM_LIT:3>, <NUM_LIT:3>])<EOL><DEDENT>else:<EOL><INDENT>m00 = M[<NUM_LIT:0>, <NUM_LIT:0>]<EOL>m01 = M[<NUM_LIT:0>, <NUM_LIT:1>]<EOL>m02 = M[<NUM_LIT:0>, <NUM_LIT:2>]<EOL>m10 = M[<NUM_LIT:1>, <NUM_LIT:0>]<EOL>m11 = M[<NUM_LIT:1>, <NUM_LIT:1>]<EOL>m12 = M[<NUM_LIT:1>, <NUM_LIT:2>]<EOL>m20 = M[<NUM_LIT:2>, <NUM_LIT:0>]<EOL>m21 = M[<NUM_LIT:2>, <NUM_LIT:1>]<EOL>m22 = M[<NUM_LIT:2>, <NUM_LIT:2>]<EOL>K = numpy.array([[m00-m11-m22, <NUM_LIT:0.0>,         <NUM_LIT:0.0>,         <NUM_LIT:0.0>],<EOL>[m01+m10,     m11-m00-m22, <NUM_LIT:0.0>,         <NUM_LIT:0.0>],<EOL>[m02+m20,     m12+m21,     m22-m00-m11, <NUM_LIT:0.0>],<EOL>[m21-m12,     m02-m20,     m10-m01,     m00+m11+m22]])<EOL>K /= <NUM_LIT><EOL>w, V = numpy.linalg.eigh(K)<EOL>q = V[[<NUM_LIT:3>, <NUM_LIT:0>, <NUM_LIT:1>, <NUM_LIT:2>], numpy.argmax(w)]<EOL><DEDENT>if q[<NUM_LIT:0>] < <NUM_LIT:0.0>:<EOL><INDENT>numpy.negative(q, q)<EOL><DEDENT>return q<EOL>
Return quaternion from rotation matrix.

    If isprecise is True, the input matrix is assumed to be a precise rotation
    matrix and a faster algorithm is used.

    >>> q = quaternion_from_matrix(numpy.identity(4), True)
    >>> numpy.allclose(q, [1, 0, 0, 0])
    True
    >>> q = quaternion_from_matrix(numpy.diag([1, -1, -1, 1]))
    >>> numpy.allclose(q, [0, 1, 0, 0]) or numpy.allclose(q, [0, -1, 0, 0])
    True
    >>> R = rotation_matrix(0.123, (1, 2, 3))
    >>> q = quaternion_from_matrix(R, True)
    >>> numpy.allclose(q, [0.9981095, 0.0164262, 0.0328524, 0.0492786])
    True
    >>> R = [[-0.545, 0.797, 0.260, 0], [0.733, 0.603, -0.313, 0],
    ...      [-0.407, 0.021, -0.913, 0], [0, 0, 0, 1]]
    >>> q = quaternion_from_matrix(R)
    >>> numpy.allclose(q, [0.19069, 0.43736, 0.87485, -0.083611])
    True
    >>> R = [[0.395, 0.362, 0.843, 0], [-0.626, 0.796, -0.056, 0],
    ...      [-0.677, -0.498, 0.529, 0], [0, 0, 0, 1]]
    >>> q = quaternion_from_matrix(R)
    >>> numpy.allclose(q, [0.82336615, -0.13610694, 0.46344705, -0.29792603])
    True
    >>> R = random_rotation_matrix()
    >>> q = quaternion_from_matrix(R)
    >>> is_same_transform(R, quaternion_matrix(q))
    True
    >>> R = euler_matrix(0.0, 0.0, numpy.pi/2.0)
    >>> numpy.allclose(quaternion_from_matrix(R, isprecise=False),
    ...                quaternion_from_matrix(R, isprecise=True))
    True