dmriprep
dmriprep copied to clipboard
Framewise displacement calculation
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.
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