# based on CPAC generate_motion_statistics and scrubbing
# check whether the motion parameter file matches the dimensions of the image.
def check(motparfile, input_image):
import numpy as np
import nibabel as nb
# load motion pars
motpardata = np.loadtxt(motparfile)
# load input image
img = nb.load(input_image)
if not motpardata.shape[0] == img.get_shape()[3]:
print motpardata.shape[0], img.get_shape()[3]
print 'The number of motion parameters does not correspond to the number of timepoints in the image.'
return False
else:
return True
# calculate FD as in Power et al. 2012
def calc_FD(motparfile, rotcols):
"""
Method to calculate Framewise Displacement (FD)
as outlined in Power et al., 2012, NeuroImage.
Parameters
----------
in_file : string
movement parameters vector file path
rotcols : list
list of indices specify which columns have the rotations, first column = 1
Returns
-------
FD : array
numpy array holding the FD values for each timepoint
"""
import numpy as np
# read in file
data = np.loadtxt(motparfile)
rots = []
# determine what the rotations are
for x in range(0, data.shape[1]):
# range started at 0, so adjust to accomodate what user entered
if x + 1 in rotcols:
if len(rots) == 0:
rots = [ y for y in data[:,x]]
else:
rots = np.column_stack([rots,data[:,x]])
# convert rot to translations on 50mm sphere
rots2trans = (rots/360)*2*50*np.pi
# determine the translations
trans = []
for x in range(0, data.shape[1]):
if not x + 1 in rotcols:
if len(trans) == 0:
trans = [ y for y in data[:,x] ]
else:
trans = np.column_stack([trans,data[:,x]])
#put back together as translations [0,1,2] and converted rotations [3,4,5]
transdata = np.hstack([trans,rots2trans])
# obtain framewise displacement
# N+1 - N
FDtransdata = transdata[1:,:] - transdata[:transdata.shape[0]-1,:]
# calculate the FD measure specified in Power et al.
# FD = sum of absolute 'translations' per timepoint
FD = [ np.sum(np.abs(transdata[x,:])) for x in range(0, transdata.shape[0]) ]
# convert to array for easy indexing
FD = np.array(FD)
return FD
# select which frames to keep and which ones should be deleted.
def prepare_scrubbing(FD, threshold=0.2, frames_before=1, frames_after=2):
"""
Based on the Framewise Displacement (FD)
as outlined in Power et al., 2012, NeuroImage, select which timepoints (frames)
to keep and which ones to delete.
Parameters
----------
FD : array
numpy array holding one FD value per timepoint
threshold : a float
scrubbing threshold value
frames_before : an integer
number of frames preceding the offending time frame
default value is 1
frames_after : an integer
number of frames following the offending time frame
default value is 2
Returns
-------
FD : array
numpy array holding the FD values for each timepoint
"""
import numpy as np
# initial exclusion
indices = [i[0] for i in (np.argwhere(FD >= 3)).tolist()]
# remove remaining indices
extra_indices = []
for i in indices:
#remove preceding frames
if i > 0 :
count = 1
while count <= frames_before:
extra_indices.append(i-count)
count+=1
#remove following frames
count = 1
while count <= frames_after:
extra_indices.append(i+count)
count+=1
# take union of indices and extra indices
frames_ex = list(set(indices) | set(extra_indices))
frames_ex.sort()
# frames to keep
frames_in = range(len(FD))
temp = list(frames_in)
# go through all frames, start from the back that way index deletion works
# remove indices from temp if they are in frames_ex
for i in sorted(frames_in, reverse = True):
if i in frames_ex:
del temp[i]
frames_in = temp
output = {'marked_frames': indices, 'frames_in': frames_in, 'frames_ex': frames_ex}
return output
# adjust motion parameters according to frames_in
def adjust_motpar(motparfile, prep_scrub_output):
from numpy import loadtxt
# read in file
data = loadtxt(motparfile)
# adjust
data = data[prep_scrub_output['frames_in'],:]
# output
print data.shape
return data
# save frames and adjusted motparfile to files for reference
def save_frames(outfile, prep_scrub_output, adjusted_motpar):
import os
from numpy import savetxt
#write frames_ex to a file
with open(os.path.join(os.path.dirname(outfile),'scrubbing_excluded_frames.1D'),'w') as myfile:
mystring = [ str(x) for x in prep_scrub_output['frames_ex'] ]
myfile.write(','.join(mystring))
#write frames_in to a file
with open(os.path.join(os.path.dirname(outfile),'scrubbing_included_frames.1D'),'w') as myfile:
mystring = [ str(x) for x in prep_scrub_output['frames_in'] ]
myfile.write(','.join(mystring))
#write marked_frames to a file
with open(os.path.join(os.path.dirname(outfile),'scrubbing_marked_frames.1D'),'w') as myfile:
mystring = [ str(x) for x in prep_scrub_output['marked_frames'] ]
myfile.write(','.join(mystring))
#write adjusted motpars
savetxt(outfile.split('.nii.gz')[0] + '_motpar.1D', adjusted_motpar, fmt='%.8f', delimiter='\t', newline='\n')
# do the actual scrubbing
def scrub(input_file, output_file, keepers):
import nibabel as nb
# load image and data
img = nb.load(input_file)
data = img.get_data()
# only keep timepoints that are in the keepers list
data = data[:,:,:,keepers]
# save image
temp_image = nb.Nifti1Image(data, img.get_affine())
nb.save(temp_image, output_file)
return data
# actions when the script is called from the command line.
# example: python2.7 /home/mrstats/maamen/DCCN/Scripts/NeuroImage/processing_ndes/motionn_scrubbing.py filtered_func_data.nii.gz scrub_test.nii.gz mc/prefiltered_func_data_mcf.par [1,2,3]
if __name__ == "__main__":
import sys
infile = sys.argv[1]
outfile = sys.argv[2]
motparfile = sys.argv[3]
rotcols = sys.argv[4]
rotcols = [ int(x) for x in rotcols.strip(']').strip('[').split(',') ]
# check consistency of files
c = check(motparfile, infile)
if c == False:
sys.exit()
# calculate FD
FD = calc_FD(motparfile, rotcols)
# determine which timepoints to keep
to_scrub = prepare_scrubbing(FD)
# scrub the data
scrub(infile, outfile, to_scrub['frames_in'])
# adjust motion parameters
adjusted_motpar = adjust_motpar(motparfile, to_scrub)
# save motion parameter
save_frames(outfile, to_scrub, adjusted_motpar)
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent