Source code for cmind.pipeline.cmind_fieldmap_process

#!/usr/bin/env python
from __future__ import division, print_function, absolute_import

[docs]def cmind_fieldmap_process(output_dir,fieldmap_nii_rads,fieldmap_nii_mag,T1_head_vol,T1_brain_vol, deltaTE=0.0023, unwarpdir='y', ForceUpdate=False, verbose=False, logger=None): """coregister functional images to a structural image, optionally using a fieldmap to correct distortions along the phase encoding direction Parameters ---------- output_dir : str directory in which to store the output fieldmap_nii_rads : str fieldmap volume in rad/s fieldmap_nii_mag : str magnitude image corresponding to `fieldmap_nii_rads` T1_head_vol : str filename of the "fixed" head image (before brain extraction) T1_brain_vol : str filename of the "fixed" brain image deltaTE : str, optional echo time difference used to generate `fieldmap_nii_rads` unwarpdir : {'x','y','z','x-','y-','z-'}, optional phase encoding direction during the readout ForceUpdate : bool,optional if True, rerun and overwrite any previously existing results 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 ------- fieldmap_vol : str processed fieldmap fieldmap_unmasked : str unmasked processed fieldmap fieldmap2struct_affine : str affine transform from fieldmap to structural space See Also -------- cmind_fieldmap_regBBR Notes ----- uses fugue and fslmaths to preprocess the fieldmap image .. figure:: ../img/IRC04H_06M008_P_1_fieldmaprads2struct.png :align: center :scale: 100% Processed fieldmap volume """ import os, warnings import numpy as np import nibabel as nib from cmind.utils.utils import input2bool from cmind.utils.file_utils import imexist from cmind.utils.nifti_utils import cmind_save_nii_mod_header from cmind.utils.cmind_utils import cmind_reg_report_img from cmind.utils.logging_utils import cmind_init_logging, log_cmd, cmind_func_info #TODO: incorporate anything from $FSLDIR/bin/fsl_prepare_fieldmap?: # Demean to avoid gross shifting & clean up edge #convert various strings to boolean ForceUpdate, verbose = input2bool([ForceUpdate, verbose]) if verbose: cmind_func_info(logger=logger) module_logger=cmind_init_logging(verbose=verbose, logger=logger) module_logger.info("Starting Fieldmap preprocessing") if not fieldmap_nii_rads: warnings.warn("No fieldmap specified. Skipping fieldmap preprocessing") return (None, None, None) if not os.path.isdir(output_dir): os.makedirs(output_dir); if(not imexist(fieldmap_nii_rads,'nifti-1')[0]): raise IOError("specified fieldmap_nii_rads, {} doesn't exist!".format(fieldmap_nii_rads)) # if(not imexist(fieldmap_nii_mag,'nifti-1')[0]): raise IOError("specified fieldmap_nii_mag, {} doesn't exist!".format(fieldmap_nii_mag)) if not deltaTE: raise ValueError("must specify deltaTE") o=output_dir; deltaTE=float(deltaTE) fieldmap_vol=os.path.join(output_dir, 'Fieldmap_medianfilt.nii.gz') fieldmap_unmasked=os.path.join(output_dir, 'fieldmaprads_unmasked.nii.gz') fieldmap2struct_affine=os.path.join(output_dir, 'fieldmap2struct.mat') Fieldmap={} if (not imexist(fieldmap_unmasked,'nifti-1')[0]) or (not os.path.exists(fieldmap2struct_affine)) or ForceUpdate: fieldmap_nii=nib.load(fieldmap_nii_mag); Fieldmap['mag']=fieldmap_nii.get_data() fieldmap_rads_nii=nib.load(fieldmap_nii_rads); Fieldmap['fmap_rads']=fieldmap_rads_nii.get_data() Fieldmap['phdiff_rad']=Fieldmap['fmap_rads']*deltaTE; nii_type=16; #float cmind_save_nii_mod_header(fieldmap_rads_nii,Fieldmap['fmap_rads'],os.path.join(output_dir,'Fieldmap.nii.gz'),nii_type) cmind_save_nii_mod_header(fieldmap_nii,Fieldmap['mag'],os.path.join(output_dir,'Fieldmap_mag.nii.gz'),nii_type) log_cmd('$FSLDIR/bin/bet2 "%s/Fieldmap_mag" "%s/Fieldmap_mag_brain" -f 0.5 -g 0 -o -m' % (o,o), verbose=verbose, logger=module_logger) #Despike and median filter the fieldmap log_cmd('$FSLDIR/bin/fugue --loadfmap="%s/Fieldmap" --despike -m --savefmap="%s/Fieldmap_medianfilt"' % (o,o), verbose=verbose, logger=module_logger) #log_cmd('$FSLDIR/bin/fugue --loadfmap=Fieldmap --despike --poly=12 --savefmap=Fieldmap_polyfit Fieldmap_vol=os.path.join(output_dir, 'Fieldmap_medianfilt') FieldmapMag_vol=os.path.join(output_dir, 'Fieldmap_mag') FieldmapMag_brain_vol=os.path.join(output_dir, 'Fieldmap_mag_brain') #PROCESS THE FIELDMAP #register fmap to structural image log_cmd('$FSLDIR/bin/flirt -in "%s" -ref "%s" -dof 6 -omat "%s/fieldmap2struct_init.mat"' % (FieldmapMag_brain_vol, T1_brain_vol,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/flirt -in "%s" -ref "%s" -dof 6 -init "%s/fieldmap2struct_init.mat" -omat "%s/fieldmap2struct.mat" -out "%s/fieldmap2struct" -nosearch' % (FieldmapMag_vol, T1_head_vol,o,o,o), verbose=verbose, logger=module_logger) cmind_reg_report_img(os.path.join(output_dir,'fieldmap2struct.nii.gz'), T1_head_vol, os.path.join(output_dir,'fieldmap_mag2struct.png'),[],1) #unmask the fieldmap (necessary to avoid edge effects) log_cmd('$FSLDIR/bin/fslmaths "%s" -abs -bin "%s/fieldmaprads_mask"' % (FieldmapMag_brain_vol,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/fslmaths "%s" -abs -bin -mul "%s/fieldmaprads_mask" "%s/fieldmaprads_mask"' % (Fieldmap_vol,o,o), verbose=verbose, logger=module_logger) if False: #custom approach to interpolating fieldmap values outside brain #addpath('/media/Data1/matlab/other_toolboxes/moving_average_v3.1') from scipy import ndimage mask_nii=nib.load(os.path.join(output_dir,'fieldmaprads_mask.nii.gz')); mask=mask_nii.get_data()>0; fmap_mask=nib.load(os.path.join(output_dir,'Fieldmap_medianfilt.nii.gz')); fmap_mask=fmap_mask.get_data(); #fmap_mask[~mask]=np.NaN; fmap_mask2=fmap_mask.copy(); if True: #nearest neighbor interpolation ind = ndimage.distance_transform_edt(~mask, return_distances=False, return_indices=True) fmap_mask2=fmap_mask2[tuple(ind)] else: #gray value dilation via max filter fmap_mask_filt=ndimage.grey_dilation(fmap_mask,size=(5,5,5)); fmap_mask2[~mask]=fmap_mask_filt[~mask]; cmind_save_nii_mod_header(fieldmap_nii,fmap_mask2,os.path.join(output_dir,'Fieldmap_customfilt'),nii_type) log_cmd('$FSLDIR/bin/fugue --loadfmap="%s/Fieldmap_customfilt" --despike -m --savefmap="%s/Fieldmap_customfilt2"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/slicer "%s/Fieldmap_customfilt2" -a "%s/Fieldmap_customfilt2.png"' % (o,o), verbose=verbose, logger=module_logger) #Fieldmap_vol_unmasked=os.path.join(output_dir, 'Fieldmap_customfilt2') #log_cmd('$FSLDIR/bin/applywarp -i Fieldmap_customfilt2 -r %s --premat=fieldmap2struct.mat -o fieldmaprads2struct_pad0_customfilt2' % (T1_head_vol), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/applywarp -i "%s/Fieldmap_customfilt2" -r "%s" --premat="%s/fieldmap2struct.mat" -o "%s/fieldmaprads2struct_pad0"' % (o,T1_head_vol,o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/immv "%s/Fieldmap_customfilt2" "%s/fieldmaprads_unmasked"' % (o,o), verbose=verbose, logger=module_logger) else: #FSL approach #TODO: come up with something better? unwarpdir=unwarpdir.lower() if unwarpdir in ['x','x-','y','y-','z','z-']: log_cmd('$FSLDIR/bin/fugue --loadfmap="%s" --mask="%s/fieldmaprads_mask" --unmaskfmap --savefmap="%s/fieldmaprads_unmasked" --unwarpdir="%s"' % (Fieldmap_vol,o,o,unwarpdir), verbose=verbose, logger=module_logger) else: raise ValueError("Invalid unwarpdir") log_cmd('$FSLDIR/bin/slicer "%s/fieldmaprads_unmasked" -a "%s/fieldmaprads_unmasked.png"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/applywarp -i "%s/fieldmaprads_unmasked" -r "%s" --premat="%s/fieldmap2struct.mat" -o "%s/fieldmaprads2struct_pad0"' % (o,T1_head_vol,o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/slicer "%s/Fieldmap" -a "%s/Fieldmap.png"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/slicer "%s/Fieldmap_medianfilt" -a "%s/Fieldmap_medianfilt.png"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/slicer "%s/fieldmaprads_mask" -a "%s/fieldmaprads_mask.png"' % (o,o), verbose=verbose, logger=module_logger) # the following is a NEW HACK to fix extrapolation when fieldmap is too small log_cmd('$FSLDIR/bin/fslmaths "%s/fieldmaprads2struct_pad0" -abs -bin "%s/fieldmaprads2struct_innermask"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/fugue --loadfmap="%s/fieldmaprads2struct_pad0" --mask="%s/fieldmaprads2struct_innermask" --unmaskfmap --unwarpdir=%s --savefmap="%s"/fieldmaprads2struct_dilated' % (o,o,unwarpdir,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/fslmaths "%s/fieldmaprads2struct_dilated" "%s/fieldmaprads2struct"' % (o,o), verbose=verbose, logger=module_logger) log_cmd('$FSLDIR/bin/slicer "%s/fieldmaprads2struct" -a "%s/fieldmaprads2struct.png"' % (o,o), verbose=verbose, logger=module_logger) cmind_reg_report_img(os.path.join(o,'fieldmaprads2struct.nii.gz'), T1_head_vol, os.path.join(o,'fieldmap2struct.png'),[],1) log_cmd('$FSLDIR/bin/imrm "%s/fieldmaprads2struct_pad0" "%s/fieldmaprads2struct_innermask" "%s/fieldmaprads2struct_dilated"' % (o,o,o), verbose=verbose, logger=module_logger) #print any filename outputs for capture by LONI print("fieldmap_vol:{}".format(fieldmap_vol)) print("fieldmap_unmasked:{}".format(fieldmap_unmasked)) print("fieldmap2struct_affine:{}".format(fieldmap2struct_affine)) return (fieldmap_vol, fieldmap_unmasked, fieldmap2struct_affine)
def _testme(): import os, tempfile from cmind.utils.logging_utils import cmind_logger from cmind.globals import cmind_example_dir,cmind_example_output_dir deltaTE=0.0023 unwarpdir="y" ForceUpdate=False; verbose=True logger=cmind_logger(log_level_console='INFO',logfile=os.path.join(tempfile.gettempdir(),'cmind_log.txt'),log_level_file='DEBUG',file_mode='w') if verbose: func=cmind_timer(cmind_fieldmap_process,logger=logger) else: func=cmind_fieldmap_process test_cases=['P','F'] for test_case in test_cases: if test_case=='P': output_dir=os.path.join(cmind_example_output_dir,'P','PhilipsGRE') fieldmap_nii_rads=os.path.join(cmind_example_dir,'IRC04H_06M008_P_1_WIP_Philips_GRE_Map_SENSE_6_1_fmap_rads.nii.gz') fieldmap_nii_mag=os.path.join(cmind_example_dir,'IRC04H_06M008_P_1_WIP_Philips_GRE_Map_SENSE_6_1_fmap_abs.nii.gz') elif test_case=='F': output_dir=os.path.join(cmind_example_output_dir,'F','PhilipsGRE') fieldmap_nii_rads=os.path.join(cmind_example_dir,'IRC04H_06M008_F_1_WIP_Philips_GRE_Map_SENSE_6_1_fmap_rads.nii.gz') fieldmap_nii_mag=os.path.join(cmind_example_dir,'IRC04H_06M008_F_1_WIP_Philips_GRE_Map_SENSE_6_1_fmap_abs.nii.gz') T1_head_vol=os.path.join(output_dir,'..','T1_proc','T1W_Structural_N4.nii.gz') T1_brain_vol=os.path.join(output_dir,'..','T1_proc','T1W_Structural_N4_brain.nii.gz') func(output_dir,fieldmap_nii_rads,fieldmap_nii_mag,T1_head_vol,T1_brain_vol, deltaTE=deltaTE, unwarpdir=unwarpdir, ForceUpdate=ForceUpdate,verbose=verbose,logger=logger) if __name__ == '__main__': import sys, argparse from cmind.utils.utils import _parser_to_loni from cmind.utils.decorators import cmind_timer if len(sys.argv) == 2 and (sys.argv[1]=='test' or sys.argv[1]=='testF'): _testme() else: parser = argparse.ArgumentParser(description="coregister functional images to a structural image, optionally using a fieldmap to correct distortions along the phase encoding direction", epilog="") #-h, --help exist by default #example of a positional argument parser.add_argument("-o","--output_dir",required=True, help="directory in which to store the output") parser.add_argument("-rads","--fieldmap_nii_rads",required=True, help="fieldmap volume in rad/s") parser.add_argument("-mag","--fieldmap_nii_mag",required=True, help="magnitude image corresponding to `fieldmap_nii_rads`") parser.add_argument("-head","-t1","--T1_head_vol",required=True, help="filename of the 'fixed' head image (before brain extraction)") parser.add_argument("-brain","-t1brain","--T1_brain_vol",required=True, help="filename of the 'fixed' brain image") parser.add_argument("-dte","--deltaTE", type=float, default=0.0023, help="echo time difference used to generate `fieldmap_nii_rads`") parser.add_argument("-unw","--unwarpdir", type=str, default="y", help= "{'x','y','z','x-','y-','z-'}, optional phase encoding direction during the readout") parser.add_argument("-u","--ForceUpdate",type=str,default="False", help="if True, rerun and overwrite any previously existing results") parser.add_argument("-v","-debug","--verbose", help="print additional output (to terminal and log)", action="store_true") parser.add_argument("-log","--logfile", dest='logger', type=str, default="", help="logging.Logger object (or string of a filename to log to)") #if first command line argument is 'loni_bash', print bash template instead of proceeding _parser_to_loni(parser, sys.argv, __file__) #kludge to auto-generate case/esac and usage statements for the corresponding LONI shell scripts args, unknown = parser.parse_known_args() #allow unrecognized arguments now for backward compatibility with a prior version output_dir=args.output_dir; fieldmap_nii_rads=args.fieldmap_nii_rads; fieldmap_nii_mag=args.fieldmap_nii_mag; T1_head_vol=args.T1_head_vol; T1_brain_vol=args.T1_brain_vol; deltaTE=args.deltaTE; unwarpdir=args.unwarpdir; ForceUpdate=args.ForceUpdate; verbose=args.verbose; logger=args.logger if verbose: cmind_fieldmap_process=cmind_timer(cmind_fieldmap_process,logger=logger) cmind_fieldmap_process(output_dir,fieldmap_nii_rads,fieldmap_nii_mag,T1_head_vol,T1_brain_vol, deltaTE=deltaTE, unwarpdir=unwarpdir, ForceUpdate=ForceUpdate,verbose=verbose,logger=logger)