In [17]:
from __future__ import division
import os
import sys
import datetime
import tempfile
import shutil
import subprocess
import logging
import nibabel
import nilearn.image
import numpy as np
from docopt import docopt
import ciftify
from ciftify.utils import WorkFlowSettings, run, get_stdout, section_header, FWHM2Sigma
from ciftify.filenames import *
logger = logging.getLogger('ciftify')
logger.setLevel(logging.DEBUG)
DRYRUN = False
ch = logging.StreamHandler()
formatter = logging.Formatter('%(message)s')
ch.setFormatter(formatter)
logger.addHandler(ch)
In [3]:
tmpdir = '/scratch/edickie/tmp_20180421'
ciftify.utils.make_dir(tmpdir)
In [4]:
arguments = {'--DilateBelowPct': None,
'--FLIRT-to-T1w': True,
'--OutputSurfDiagnostics': True,
'--SmoothingFWHM': '12',
'--already-in-MNI': False,
'--ciftify-conf': None,
'--ciftify-work-dir': None,
'--debug': True,
'--dry-run': False,
'--func-ref': 'first_vol',
'--hcp-data-dir': '/scratch/edickie/ciftify_intergration_tests/run_2018-04-16/hcp',
'--help': False,
'--surf-reg': 'FS',
'--verbose': False,
'<func.nii.gz>': '/scratch/edickie/ciftify_intergration_tests/src_data/ds000030_R1.0.4/derivatives/fmriprep/sub-50005/func/sub-50005_task-rest_bold_space-native_preproc.nii.gz',
'<subject>': 'sub-50005',
'<task_label>': 'rest_test2'}
In [5]:
from ciftify.bin.ciftify_subject_fmri import Settings
settings = Settings(arguments)
In [23]:
import yaml
print(yaml.dump(settings.__dict__))
In [6]:
def define_func_3D(settings, tmpdir):
"""
create the func_3D file as per subject instrutions
"""
logger.info(section_header("Getting fMRI reference volume"))
if settings.func_ref.mode == "first_vol":
'''take the fist image (default)'''
native_func_3D = os.path.join(tmpdir, "native_func_first.nii.gz")
run(['wb_command', '-volume-math "(x)"', native_func_3D,
'-var','x', settings.func_4D, '-subvolume', '1'])
elif settings.func_ref.mode == "median":
'''take the median over time, if indicated'''
native_func_3D = os.path.join(tmpdir, "native_func_median.nii.gz")
run(['wb_command', '-volume-reduce',
settings.func_4D, 'MEDIAN', native_func_3D])
elif settings.func_ref.mode == "path":
'''use the file indicated by the user..after checking the dimension'''
cifify.meants.verify_nifti_dimensions_match(settings.ref_vol.path,
settings.func_4D)
native_func_3D = settings.func_ref.path
else:
sys.exit('Failed to define the ref volume')
return native_func_3D
In [13]:
def calc_sform_differences(native_func_3D, settings, tmpdir):
"""
adjust for any differences between the func sform and the T1w image
using nilearn resample to find and intermidiate
"""
logger.info('---Adjusting for differences between underlying sform---')
resampled_ref = os.path.join(tmpdir, 'vol_ref_rT1w.nii.gz')
logger.info("Using nilearn to create resampled image {}".format(resampled_ref))
resampled_ref_vol = nilearn.image.resample_to_img(
source_img = settings.ref_vol.path,
target_img = os.path.join(
settings.vol_reg['src_dir'],
settings.vol_reg['T1wImage']))
resampled_ref_vol.to_filename(resampled_ref)
logger.info("Calculating linear transform between resampled reference and reference vols")
func2T1w_mat = os.path.join(settings.results_dir, 'native','mat_EPI_to_T1.mat')
run(['mkdir','-p',os.path.join(settings.results_dir,'native')])
run(['flirt',
'-in', native_func_3D,
'-ref', resampled_ref,
'-omat', func2T1w_mat,
'-dof', "6",
'-cost', "corratio", '-searchcost', "corratio"])
return func2T1w_mat
In [36]:
def run_flirt_to_T1w(native_func_3D, settings,
cost_function = "corratio", degrees_of_freedom = "12"):
"""
Use FSL's FLIRT to calc a transform to the T1w Image.. not ideal transform
"""
logger.info('Running FLIRT transform to T1w with costfunction {} and dof {}'.format(cost_function, degrees_of_freedom))
func2T1w_mat = os.path.join(settings.results_dir, 'native','mat_EPI_to_T1.mat')
## calculate the fMRI to native T1w transform
run(['mkdir','-p',os.path.join(settings.results_dir,'native')])
run(['flirt',
'-in', native_func_3D,
'-ref', os.path.join(
settings.vol_reg['src_dir'],
settings.vol_reg['T1wBrain']),
'-omat', func2T1w_mat,
'-o', os.path.join(tmpdir,"epi_in_T1w.nii.gz"),
'-dof', str(degrees_of_freedom),
'-cost', cost_function, '-searchcost', cost_function,
'-searchrx', '-180', '180',
'-searchry', '-180', '180',
'-searchrz', '-180', '180'])
return func2T1w_mat
In [30]:
def transform_to_MNI(func2T1w_mat, atlas_fMRI_4D, settings):
'''
transform the fMRI image to MNI space 2x2x2mm using FSL
RegTemplate An optional 3D MRI Image from the functional to use for registration
'''
## make the directory to hold the transforms if it doesn't exit
run(['mkdir','-p',os.path.join(settings.results_dir,'native')])
## concatenate the transforms
func2MNI_mat = os.path.join(settings.results_dir,'native','mat_EPI_to_TAL.mat')
run(['convert_xfm','-omat', func2MNI_mat, '-concat',
os.path.join(
settings.vol_reg['xfms_dir'],
settings.vol_reg['AtlasTransform_Linear']),
func2T1w_mat])
## now apply the warp!!
logger.info("Transforming to MNI space and resampling to 2x2x2mm")
run(['applywarp',
'--ref={}'.format(os.path.join(
settings.vol_reg['dest_dir'],
settings.vol_reg['T1wImage'])),
'--in={}'.format(settings.func_4D),
'--warp={}'.format(os.path.join(
settings.vol_reg['xfms_dir'],
settings.vol_reg['AtlasTransform_NonLinear'])),
'--premat={}'.format(func2MNI_mat),
'--interp=spline',
'--out={}'.format(atlas_fMRI_4D)])
In [29]:
# define some mesh paths as a dictionary
meshes = define_meshes(settings.subject.path, tmpdir,
low_res_meshes = settings.low_res)
logger.info('The functional data will first be projected to '
'the surfaces in {}'.format(meshes['AtlasSpaceNative']['Folder']))
logger.info('The data is then resampled to the {} surfaces in {} folder using the {} sphere'.format(
meshes['AtlasSpaceNative']['meshname'],
meshes['AtlasSpaceNative']['Folder'],
settings.surf_reg))
if not settings.diagnostics.path:
settings.diagnostics.path = tmpdir
else:
logger.info('cifti files for surface mapping QA will'
' be written to {}:'.format(settings.diagnostics.path))
## either transform or copy the input_fMRI
native_func_3D = define_func_3D(settings, tmpdir)
In [37]:
atlas_fMRI_4D = os.path.join(settings.results_dir,
'{}.nii.gz'.format(settings.fmri_label))
if settings.already_atlas_transformed:
run(['cp', settings.func_4D, atlas_fMRI_4D])
else:
logger.info(section_header('MNI Transform'))
if not settings.run_flirt:
func2T1w_mat = calc_sform_differences(native_func_3D, settings, tmpdir)
else:
func2T1w_mat = run_flirt_to_T1w(native_func_3D, settings)
transform_to_MNI(func2T1w_mat, atlas_fMRI_4D, settings)
In [38]:
def make_cortical_ribbon(ref_vol, ribbon_vol, settings, mesh_settings):
''' make left and right cortical ribbons and combine '''
logger.info(section_header('Making fMRI Ribbon'))
with ciftify.utils.TempDir() as ribbon_tmp:
for Hemisphere in ['L', 'R']:
hemisphere_cortical_ribbon(Hemisphere, settings.subject.id,
ref_vol, mesh_settings,
os.path.join(ribbon_tmp,'{}.ribbon.nii.gz'.format(Hemisphere)),
ribbon_tmp)
# combine the left and right ribbons into one mask
run(['fslmaths', os.path.join(ribbon_tmp,'L.ribbon.nii.gz'),
'-add', os.path.join(ribbon_tmp,'R.ribbon.nii.gz'),
ribbon_vol])
In [40]:
def hemisphere_cortical_ribbon(hemisphere, subject, ref_vol, mesh_settings,
ribbon_out, ribbon_tmp, GreyRibbonValue = 1):
'''
builds a cortical ribbon mask for that hemisphere
'''
## create a volume of distances from the surface
tmp_white_vol = os.path.join(ribbon_tmp,'{}.white.native.nii.gz'.format(hemisphere))
tmp_pial_vol = os.path.join(ribbon_tmp,'{}.pial.native.nii.gz'.format(hemisphere))
run(['wb_command', '-create-signed-distance-volume',
surf_file(subject, 'white', hemisphere, mesh_settings),
ref_vol, tmp_white_vol])
run(['wb_command', '-create-signed-distance-volume',
surf_file(subject, 'pial', hemisphere, mesh_settings),
ref_vol, tmp_pial_vol])
## threshold and binarise these distance files
tmp_white_vol_thr = os.path.join(ribbon_tmp,'{}.white_thr0.native.nii.gz'.format(hemisphere))
tmp_pial_vol_thr = os.path.join(ribbon_tmp,'{}.pial_uthr0.native.nii.gz'.format(hemisphere))
run(['fslmaths', tmp_white_vol, '-thr', '0', '-bin', '-mul', '255',
tmp_white_vol_thr])
run(['fslmaths', tmp_white_vol_thr, '-bin', tmp_white_vol_thr])
run(['fslmaths', tmp_pial_vol, '-uthr', '0', '-abs', '-bin', '-mul', '255',
tmp_pial_vol_thr])
run(['fslmaths', tmp_pial_vol_thr, '-bin', tmp_pial_vol_thr])
## combine the pial and white to get the ribbon
run(['fslmaths', tmp_pial_vol_thr,
'-mas', tmp_white_vol_thr,
'-mul', '255',
ribbon_out])
run(['fslmaths', ribbon_out, '-bin', '-mul', str(GreyRibbonValue), ribbon_out])
In [41]:
logger.info(section_header('Making fMRI Ribbon'))
ribbon_vol=os.path.join(settings.diagnostics.path,'ribbon_only.nii.gz')
make_cortical_ribbon(ref_vol = os.path.join(settings.vol_reg['dest_dir'],
settings.vol_reg['T1wImage']),
ribbon_vol = ribbon_vol,
settings = settings,
mesh_settings = meshes['AtlasSpaceNative'])
In [42]:
from ciftify.bin.ciftify_subject_fmri import define_good_voxels
In [43]:
logger.info(section_header('Determining Noisy fMRI voxels'))
goodvoxels_vol = os.path.join(settings.diagnostics.path, 'goodvoxels.nii.gz')
tmean_vol, cov_vol = define_good_voxels(
atlas_fMRI_4D, ribbon_vol, goodvoxels_vol, tmpdir)