Framewise displacement calculation

[oesteban]

Based on the estimation of #94, extract some score of the relative rigid-body misalignment of each orientation.

We will identify volumes that are outliers in terms of head motion, or other severe artifacts that make them likely candidates for exclusion from further analysis.

[dPys]

Here's some of my scratchwork we can work with-- it will need to be modified to use rasb-ts. get_params may be a part of niworkflows already?

def id_outliers_fn(outlier_report, threshold, dwi_file):
    """Get list of scans that exceed threshold for number of outliers
    Parameters
    ----------
    outlier_report: string
        Path to the fsl_eddy outlier report
    threshold: int or float
        If threshold is an int, it is treated as number of allowed outlier
        slices. If threshold is a float between 0 and 1 (exclusive), it is
        treated the fraction of allowed outlier slices before we drop the
        whole volume.
    dwi_file: string
        Path to nii dwi file to determine total number of slices
    Returns
    -------
    drop_scans: numpy.ndarray
        List of scan indices to drop
    """
    import nibabel as nib
    import numpy as np
    import os.path as op
    import parse

    with open(op.abspath(outlier_report), "r") as fp:
        lines = fp.readlines()

    p = parse.compile(
        "Slice {slice:d} in scan {scan:d} is an outlier with "
        "mean {mean_sd:f} standard deviations off, and mean "
        "squared {mean_sq_sd:f} standard deviations off."
    )

    outliers = [p.parse(l).named for l in lines]
    scans = {d["scan"] for d in outliers}

    def num_outliers(scan, outliers):
        return len([d for d in outliers if d["scan"] == scan])

    if 0 < threshold < 1:
        img = nib.load(dwi_file)
        try:
            threshold *= img.header.get_n_slices()
        except nib.spatialimages.HeaderDataError:
            print(
                "WARNING. We are not sure which dimension has the "
                "slices in this image. So we are using the 3rd dim.",
                img.shape,
            )
            threshold *= img.shape[2]

    drop_scans = np.array([s for s in scans if num_outliers(s, outliers) > threshold])

    outpath = op.abspath("dropped_scans.txt")
    np.savetxt(outpath, drop_scans, fmt="%d")

    return drop_scans, outpath


def drop_outliers_fn(in_file, in_bval, in_bvec, drop_scans, in_sigma=None, perc_missing=0.15):
    """Drop outlier volumes from dwi file
    Parameters
    ----------
    in_file: string
        Path to nii dwi file to drop outlier volumes from
    in_bval: string
        Path to bval file to drop outlier volumes from
    in_bvec: string
        Path to bvec file to drop outlier volumes from
    drop_scans: numpy.ndarray or str
        List of scan indices to drop. If str, then assume path to text file
        listing outlier volumes.

    Returns
    -------
    out_file: string
        Path to "thinned" output dwi file
    out_bval: string
        Path to "thinned" output bval file
    out_bvec: string
        Path to "thinned" output bvec file
    """
    import nibabel as nib
    import numpy as np
    import os.path as op
    from nipype.utils.filemanip import fname_presuffix

    if isinstance(drop_scans, str):
        try:
            drop_scans = np.genfromtxt(drop_scans).tolist()
            if not isinstance(drop_scans, list):
                drop_scans = [drop_scans]

        except TypeError:
            print(
                "Unrecognized file format. Unable to load vector from drop_scans file."
            )

    print("%s%s" % ('Dropping outlier volumes:\n', drop_scans))

    img = nib.load(op.abspath(in_file))
    drop_perc = (len(drop_scans))/float(img.shape[-1])
    if drop_perc > perc_missing:
        raise ValueError('Missing > ' + str(np.round(100*drop_perc, 2)) + '% of volumes after outlier removal. '
                                                                          'This dataset is unuseable based on the '
                                                                          'given rejection threshold.')
    # drop 4d outliers from dwi
    img_data = img.get_data()
    img_data_thinned = np.delete(img_data, drop_scans, axis=3)
    if isinstance(img, nib.nifti1.Nifti1Image):
        img_thinned = nib.Nifti1Image(
            img_data_thinned.astype(np.float64), img.affine, header=img.header
        )
    elif isinstance(img, nib.nifti2.Nifti2Image):
        img_thinned = nib.Nifti2Image(
            img_data_thinned.astype(np.float64), img.affine, header=img.header
        )
    else:
        raise TypeError("in_file does not contain Nifti image datatype.")

    out_file = fname_presuffix(in_file, suffix="_thinned", newpath=op.abspath("."))
    nib.save(img_thinned, op.abspath(out_file))

    # drop outliers from sigma (if 4d)
    if in_sigma is not None:
        sigma = np.load(op.abspath(in_sigma))
        if len(sigma.shape) == 4:
            sigma_thinned = np.delete(sigma, drop_scans, axis=3)
            out_sigma = fname_presuffix(in_sigma, suffix="_thinned", newpath=op.abspath("."))
            np.save(sigma_thinned, op.abspath(out_sigma))
        else:
            out_sigma = in_sigma
    else:
        out_sigma = None

    bval = np.loadtxt(in_bval)
    bval_thinned = np.delete(bval, drop_scans, axis=0)
    out_bval = fname_presuffix(in_bval, suffix="_thinned", newpath=op.abspath("."))
    np.savetxt(out_bval, bval_thinned.astype('int'), fmt='%i')

    bvec = np.loadtxt(in_bvec)
    if bvec.shape[0] == 3:
        bvec_thinned = np.delete(bvec, drop_scans, axis=1)
    else:
        bvec_thinned = np.delete(bvec, drop_scans, axis=0)
    out_bvec = fname_presuffix(in_bvec, suffix="_thinned", newpath=op.abspath("."))
    np.savetxt(out_bvec, bvec_thinned.astype('float'), fmt='%10f')

    return out_file, out_bval, out_bvec, out_sigma


def drop_scans_from_4d(in_file, drop_scans):
    import nibabel as nib
    import numpy as np
    import os.path as op
    from nipype.utils.filemanip import fname_presuffix
    img = nib.load(op.abspath(in_file))
    img_data = img.get_data()
    img_data_thinned = np.delete(img_data, drop_scans, axis=3)
    if isinstance(img, nib.nifti1.Nifti1Image):
        img_thinned = nib.Nifti1Image(
            img_data_thinned.astype(np.float64), img.affine, header=img.header
        )
    elif isinstance(img, nib.nifti2.Nifti2Image):
        img_thinned = nib.Nifti2Image(
            img_data_thinned.astype(np.float64), img.affine, header=img.header
        )
    else:
        raise TypeError("in_file does not contain Nifti image datatype.")

    out_file = fname_presuffix(in_file, suffix="_thinned", newpath=op.abspath("."))
    nib.save(img_thinned, op.abspath(out_file))
    return out_file


def get_params(A):
    """This is a copy of spm's spm_imatrix where
    we already know the rotations and translations matrix,
    shears and zooms (as outputs from fsl FLIRT/avscale)
    Let A = the 4x4 rotation and translation matrix
    R = [          c5*c6,           c5*s6, s5]
        [-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
        [-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]
    """

    def rang(b):
        a = min(max(b, -1), 1)
        return a

    Ry = np.arcsin(A[0, 2])
    # Rx = np.arcsin(A[1, 2] / np.cos(Ry))
    # Rz = np.arccos(A[0, 1] / np.sin(Ry))

    if (abs(Ry) - np.pi / 2) ** 2 < 1e-9:
        Rx = 0
        Rz = np.arctan2(-rang(A[1, 0]), rang(-A[2, 0] / A[0, 2]))
    else:
        c = np.cos(Ry)
        Rx = np.arctan2(rang(A[1, 2] / c), rang(A[2, 2] / c))
        Rz = np.arctan2(rang(A[0, 1] / c), rang(A[0, 0] / c))

    rotations = [Rx, Ry, Rz]
    translations = [A[0, 3], A[1, 3], A[2, 3]]

    return rotations, translations


def get_flirt_motion_parameters(flirt_mats):
    import os.path as op
    from nipype.interfaces.fsl.utils import AvScale
    from dmriprep.utils.core import get_params

    motion_params = open(op.abspath("motion_parameters.par"), "w")
    for mat in flirt_mats:
        res = AvScale(mat_file=mat).run()
        A = np.asarray(res.outputs.rotation_translation_matrix)
        rotations, translations = get_params(A)
        for i in rotations + translations:
            motion_params.write("%f " % i)
        motion_params.write("\n")
    motion_params.close()
    motion_params = op.abspath("motion_parameters.par")
    return motion_params