#!/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)
