def rotation_matrix_from_point(point, ret_inv=False):
    """Compute rotation matrix to go from [0,0,1] to `point`.

    First, the matrix rotates to in the polar direction. Then,
    a rotation about the y-axis is performed to match the
    azimuthal angle in the x-z-plane.

    This rotation matrix is required for the correct 3D orientation
    of the backpropagated projections.

    Parameters
    ----------
    points: list-like, length 3
        The coordinates of the point in 3D.
    ret_inv: bool
        Also return the inverse of the rotation matrix. The inverse
        is required for :func:`scipy.ndimage.interpolation.affine_transform`
        which maps the output coordinates to the input coordinates.

    Returns
    -------
    Rmat [, Rmat_inv]: 3x3 ndarrays
        The rotation matrix that rotates [0,0,1] to `point` and
        optionally its inverse.
    """
    x, y, z = point
    # azimuthal angle
    phi = np.arctan2(x, z)
    # angle in polar direction (negative)
    theta = -np.arctan2(y, np.sqrt(x**2+z**2))

    # Rotation in polar direction
    Rtheta = np.array([
        [1,             0,              0],
        [0, np.cos(theta), -np.sin(theta)],
        [0, np.sin(theta),  np.cos(theta)],
    ])

    # rotation in x-z-plane
    Rphi = np.array([
        [np.cos(phi), 0, -np.sin(phi)],
        [0, 1,            0],
        [np.sin(phi), 0, np.cos(phi)],
    ])

    D = np.dot(Rphi, Rtheta)
    # The inverse of D
    Dinv = np.dot(Rtheta.T, Rphi.T)

    if ret_inv:
        return D, Dinv
    else:
        return D