Source code for cmind.utils.fsl_mcflirt_opt_ASLBOLD

#!/usr/bin/env python
"""
**mcflirt motion correction for simultaneous ASL/BOLD data**
"""
from __future__ import division, print_function, absolute_import
import subprocess
import os
import glob
import shutil
import multiprocessing  #can run ASL and BOLD timeseries processing in parallel
from time import time

import numpy as np

from cmind.utils.file_utils import imexist, split_multiple_ext
from cmind.utils.cmind_utils import find_optimal_FSL_refvol, run_mcflirt, run_mcflirt_parallel
from cmind.utils.logging_utils import log_cmd

#scriptpath = '/home/lee8rx/src_repositories/svn_stuff/cmind-matlab/trunk/lib/'  #os.getcwd(); #"../Test/MyModule.py"
#sys.path.append(os.path.abspath(scriptpath))

VERBOSE=True  #print extra timing info

[docs]def fsl_mcflirt_opt_ASLBOLD(ASLinput_name,BOLDinput_name,cost_type='normcorr',out_name_ASL="",out_name_BOLD="",FORCE_SINGLEPROCESS=False,output_dir=None, verbose=False, logger=None): """Utility that tries to improve the robustness of mcflirt by running in 3 stages: This version runs ASL and BOLD registrations in parallel, forcing the same reference frame to be used Parameters ---------- ASLinput_name : str ASL 4D NIFTI or NIFTI-GZ volume BOLDinput_name : str BOLD 4D NIFTI or NIFTI-GZ volume cost_type : {'normcorr','mutualinfo','woods','corratio','normmi','leastsquares'}, optional cost function used by mcflirt (default = 'normcorr') out_name_ASL : str, optional output filename to be used for the ASL motion-corrected timeseries out_name_BOLD : str, optional output filename to be used for the BOLD motion-corrected timeseries output_dir : str, optional if None, output will be stored in current working directory via os.getcwd() FORCE_SINGLEPROCESS : bool, optional if False will try to run ASL & BOLD motion correction in parallel via multiprocessing verbose : bool, optional print additional output (to terminal and log) logger : logging.Logger or str, optional logging.Logger object (or string of a filename to log to) Returns ------- out_name_ASL : str output filename of the ASL motion-corrected timeseries out_name_BOLD : str output filename of the BOLD motion-corrected timeseries Notes ----- The mcflirt stages run are: 1.) Fast initial registration to initial frame without sinc interpolation to find the frame closest to mean position 2.) Rerun mcflirt using this "optimal" reference frame and with sinc interpolation to reduce blurring 3.) Take a timeseries average of the reference frames and then run mcflirt a 3rd time using this high SNR timeseries average as the reference mcflirt will duplicate the top and bottom slice during the calculations to allow some corrections to be made to the edge slices cost_type options are: normcorr (default), mutualinfo, woods, corratio, normmi, leastsquares """ if not FORCE_SINGLEPROCESS: cpus=multiprocessing.cpu_count() #of CPUs on the system (actually seems to be number of threads for CPUs with hyperthreading) if cpus>1: FORCE_SINGLEPROCESS=False else: FORCE_SINGLEPROCESS=True if cost_type=="": cost_type='normcorr' if VERBOSE: starttime = time() (exists,fname_full_ASL,fname_root_ASL,ASLpth,ASLext)=imexist(ASLinput_name,'nifti'); if not exists: raise IOError("Volume {} doesn't exist".format(ASLinput_name)) (exists,fname_full_BOLD,fname_root_BOLD,BOLDpth,BOLDext)=imexist(BOLDinput_name,'nifti'); if not exists: raise IOError("Volume {} doesn't exist".format(BOLDinput_name)) #Note: FSL doesn't care if inputs have no file extension or .nii or .nii.gz if not output_dir: #avoid os.getcwd() if output_dir is passed as input output_dir=os.getcwd() if not out_name_ASL: #reuse input ASL filename, but with _mcf appended out_pth_ASL=output_dir fname_out_root_ASL=os.path.join(out_pth_ASL,fname_root_ASL + '_mcf') out_name_ASL=fname_out_root_ASL+ASLext #fname_full #os.path.join(os.getcwd(),fname_root) else: fname_out_root_ASL=split_multiple_ext(out_name_ASL)[0] #.replace('.nii','').replace('.gz','') out_pth_ASL=os.path.dirname(out_name_ASL) if not out_pth_ASL: out_pth_ASL=output_dir out_name_ASL=os.path.join(out_pth_ASL,out_name_ASL) if not out_name_BOLD: #reuse input BOLD filename, but with _mcf appended out_pth_BOLD=output_dir fname_out_root_BOLD=os.path.join(out_pth_BOLD,fname_root_BOLD + '_mcf') out_name_BOLD=fname_out_root_BOLD+BOLDext #fname_full #os.path.join(os.getcwd(),fname_root) else: fname_out_root_BOLD=split_multiple_ext(out_name_BOLD)[0] #.replace('.nii','').replace('.gz','') out_pth_BOLD=os.path.dirname(out_name_BOLD) if not out_pth_BOLD: out_pth_BOLD=output_dir out_name_BOLD=os.path.join(out_pth_BOLD,out_name_BOLD) if FORCE_SINGLEPROCESS: run_mcflirt(fname_full_ASL,out_name_ASL,cost_type) run_mcflirt(fname_full_BOLD,out_name_BOLD,cost_type) else: #print((fname_full_ASL,out_name_ASL,cost_type,0,'',verbose,None)) #print((fname_full_BOLD,out_name_BOLD,cost_type,0,'',verbose,None)) run_mcflirt_parallel([(fname_full_ASL,out_name_ASL,cost_type,0,'',verbose,None),(fname_full_BOLD,out_name_BOLD,cost_type,0,'',verbose,None)]) if VERBOSE: print('MCFLIRT 1 of 3, time elapsed:', time() - starttime) (my_refvol_ASL,dist_ASL) = find_optimal_FSL_refvol(out_name_ASL + '.par')[0:2] (my_refvol_BOLD,dist_BOLD) = find_optimal_FSL_refvol(out_name_BOLD + '.par')[0:2] my_refvol=np.argmin(dist_ASL+dist_BOLD) if VERBOSE: starttime = time() #redo the registration using this "optimal" minimal-distance reference volume if FORCE_SINGLEPROCESS: run_mcflirt(fname_full_ASL,out_name_ASL,cost_type,my_refvol, verbose=verbose, logger=logger) run_mcflirt(fname_full_BOLD,out_name_BOLD,cost_type,my_refvol, verbose=verbose, logger=logger) else: run_mcflirt_parallel([(fname_full_ASL,out_name_ASL,cost_type,my_refvol,'',verbose, None),(fname_full_BOLD,out_name_BOLD,cost_type,my_refvol,'', verbose, None)]) if VERBOSE: print('MCFLIRT 2 of 3, time elapsed:', time() - starttime) #Now, may want to rerun the registration to this higher SNR, timeseries averaged frame #probably of minimal benefit, but doesn't take too long to run if VERBOSE: starttime = time() #Create a time-series average of the "optimally" referenced timecourse if FORCE_SINGLEPROCESS: out = log_cmd('$FSLDIR/bin/fslmaths "%s" -Tmean "%s_avg"' % (out_name_ASL,out_name_ASL), verbose=verbose, logger=logger) out = log_cmd('$FSLDIR/bin/fslmaths "%s" -Tmean "%s_avg"' % (out_name_BOLD,out_name_BOLD), verbose=verbose, logger=logger) else: print("computing Tmeans") p1 = subprocess.Popen('$FSLDIR/bin/fslmaths "%s" -Tmean "%s_avg"' % (out_name_ASL,out_name_ASL), shell=True) p2 = subprocess.Popen('$FSLDIR/bin/fslmaths "%s" -Tmean "%s_avg"' % (out_name_BOLD,out_name_BOLD), shell=True) for p in [p1, p2]: p.wait() if VERBOSE: print('Tmean, time elapsed:', time() - starttime) starttime = time() #finally, rerun the registration, but to the average image from the previous registration if FORCE_SINGLEPROCESS: run_mcflirt(fname_full_ASL,out_name_ASL,cost_type,refvol=out_name_ASL+"_avg",extra_args="-sinc_final",verbose=verbose,logger=logger) run_mcflirt(fname_full_BOLD,out_name_BOLD,cost_type,refvol=out_name_BOLD+"_avg",extra_args="-sinc_final",verbose=verbose,logger=logger) else: run_mcflirt_parallel([(fname_full_ASL,out_name_ASL,cost_type,out_name_ASL+"_avg","-sinc_final",verbose,None), (fname_full_BOLD,out_name_BOLD,cost_type,out_name_BOLD+"_avg","-sinc_final",verbose,None)]) if VERBOSE: print('MCFLIRT 3 of 3, time elapsed:', time() - starttime) #out = log_cmd('$FSLDIR/bin/mcflirt -in "%s" -out "%s" -plots -rmsrel -rmsabs -reffile "%s_avg_mcf" -sinc_final -cost %s' % (fname_full_ASL,out_name_ASL,out_name_ASL,cost_type), verbose=verbose, logger=logger) #out = log_cmd('$FSLDIR/bin/mcflirt -in "%s" -out "%s" -plots -rmsrel -rmsabs -reffile "%s_avg_mcf" -sinc_final -cost %s' % (fname_full_BOLD,out_name_BOLD,out_name_BOLD,cost_type), verbose=verbose, logger=logger) if(fname_out_root_ASL.find(' ')==-1): #for safety don't allow spaces in the argument to rm mcf_fnames=glob.glob(os.path.join(out_pth_ASL,fname_out_root_ASL,'*_avg.*')); [os.remove(f) for f in mcf_fnames] temp_dirs=[os.path.join(out_pth_ASL,o) for o in os.listdir(out_pth_ASL) if (os.path.isdir(os.path.join(out_pth_ASL,o)) and o.find('.mat')!=-1)] #find all subdirectories with .mat in the directory name [shutil.rmtree(d) for d in temp_dirs] if(fname_out_root_BOLD.find(' ')==-1): #for safety don't allow spaces in the argument to rm mcf_fnames=glob.glob(os.path.join(out_pth_BOLD,fname_out_root_BOLD,'*_avg.*')); [os.remove(f) for f in mcf_fnames] temp_dirs=[os.path.join(out_pth_BOLD,o) for o in os.listdir(out_pth_BOLD) if (os.path.isdir(os.path.join(out_pth_BOLD,o)) and o.find('.mat')!=-1)] #find all subdirectories with .mat in the directory name [shutil.rmtree(d) for d in temp_dirs] #os.rmdir(d) #note: doesn't work with non-empty directories return (out_name_ASL,out_name_BOLD)