This tutorial demonstrates how to run a single registration for each subject and then use this in your first-level FEAT analyses. The benefits of this approach are:
FEAT first level analyses finish faster as they don't run registration, especially if using FNIRT non-linear registration to standard space.
The FEAT behavior whereby higher-level analyses are forced to use a standard space can be avoided.
Registrations will be perfectly consistent across each subject's first-level analysis directories.
The procedure for this approach is the following:
Register the standard brain to the high-res anatomical image.
Register the anatomical image to the functional image.
Place all registration files in a single directory.
Run FEAT first-level analyses without performing registration or post-stats.
Make symbolic links to our registration directory inside each lower-level FEAT directory.
If desired, run post-stats on the lower-level FEAT directories.
Run FEAT higher-level analyses as normal.
This code is in python and assumes that you have already installed FSL and FreeSurfer tools and they are accessible via the command line. It will require minimal modification if you have arranged your data into BIDS format and performed surface reconstruction using Freesurfer, but can be easily adapted to suit your own project file structure. Note that not all these steps may be necessary, but I have tried to include all the reg outputs that FEAT automatically makes, plus transformations between all spaces in both directions.
First, import the necessary modules.
import os
import os.path as op
import glob
import shutil
import subprocessNow state paths to the data you should already have.
PROJ_DIR = '/path/to/your/project/directory'
subject = 'M015'
# native functional space reference image
ref_func = f'/path/to/func_space_reference.nii.gz'
# subject's FreeSurfer directory and anatomical space reference images
fs_dir = f'{os.environ["SUBJECTS_DIR"]}/sub-{subject}'
ref_anat = f'{fs_dir}/mri/orig/001.nii'
ref_anat_brain = f'{fs_dir}/mri/orig/001_brain.nii.gz'
# standard space reference images
ref_std = f'{os.environ["FSLDIR"]}/data/standard/MNI152_T1_2mm.nii.gz'
ref_std_brain = f'{ref_std[:-7]}_brain.nii.gz'
ref_std_mask = f'{ref_std[:-7]}_brain_mask_dil.nii.gz'Next, we will perform the registration between the subject's high-res anatomical space and the standard MNI152 2mm space. Outputs will be saved in the subject's freesurfer directory so that they can be used in other projects. The contents of this directory will all end up in a 'reg' dir that FSL FEAT uses, so we need to use their naming scheme.
# make directory for outputs of this transformation
fnirt_dir = f'{fs_dir}/mri/transforms/fnirt'
os.makedirs(fnirt_dir, exist_ok=True)
# make links to the reference images
for in_path, label in zip(
[ref_anat, ref_std, ref_std_brain],
['highres', 'standard_head', 'standard']):
out_path = f'{fnirt_dir}/{label}.nii.gz'
if not op.exists(out_path):
os.system(f'ln -s {op.abspath(in_path)} {out_path}')
# use FLIRT to obtain a linear transformation as a starting point for FNIRT
highres2standard = f'{fnirt_dir}/highres2standard.mat'
if not op.isfile(highres2standard):
os.system(f'flirt '
f'-in {ref_anat} '
f'-ref {ref_std} '
f'-omat {highres2standard} '
f'-cost corratio '
f'-dof 12 '
f'-searchrx -90 90 '
f'-searchry -90 90 '
f'-searchrz -90 90 '
f'-interp trilinear')
standard2highres = f'{fnirt_dir}/standard2highres.mat'
# invert the transformation
if not op.isfile(standard2highres):
os.system(f'convert_xfm '
f'-inverse '
f'-omat {standard2highres} '
f'{highres2standard}')
# perform FNIRT using the standard T1_2_MNI152_2mm config
highres2standard_warp = f'{fnirt_dir}/highres2standard_warp.nii.gz'
if not op.isfile(highres2standard_warp):
os.system(f'fnirt '
f'--in={ref_anat} '
f'--ref={ref_std} '
f'--refmask={ref_std_mask} '
f'--config=T1_2_MNI152_2mm '
f'--aff={highres2standard} '
f'--cout={highres2standard_warp} '
f'--iout={fnirt_dir}/highres2standard_head '
f'--jout={fnirt_dir}/highres2highres_jac '
f'--warpres=10,10,10')
standard2highres_warp = f'{fnirt_dir}/standard2highres_warp.nii.gz'
# invert the transformation
standard2highres_warp = f'{fnirt_dir}/standard2highres_warp.nii.gz'
if not op.isfile(standard2highres_warp):
os.system(f'invwarp '
f'-w {highres2standard_warp} '
f'-o {standard2highres_warp} '
f'-r {ref_anat}')
# Apply the inverted transformation to the anatomical image
highres2standard_img = f'{fnirt_dir}/highres2standard.nii.gz'
if not op.isfile(highres2standard_img):
os.system(f'applywarp '
f'-i {ref_anat_brain} '
f'-r {ref_std_brain} '
f'-o {highres2standard_img} '
f'-w {highres2standard_warp}')The quality of the transformations can be inspected by creating the same plots that FEAT makes in its registration report. An example is shown at the top of this page.
# make reg images
d = fnirt_dir
if not op.isfile(f'{d}/highres2standard.png'):
slicer_str = (
f'-x 0.35 {d}/sla.png -x 0.45 {d}/slb.png '
f'-x 0.55 {d}/slc.png -x 0.65 {d}/sld.png '
f'-y 0.35 {d}/sle.png -y 0.45 {d}/slf.png '
f'-y 0.55 {d}/slg.png -y 0.65 {d}/slh.png '
f'-z 0.35 {d}/sli.png -z 0.45 {d}/slj.png '
f'-z 0.55 {d}/slk.png -z 0.65 {d}/sll.png')
append_str = ' + '.join([f'{d}/sl{i}.png' for i in 'abcdefghijkl'])
cmd = f'slicer {d}/highres2standard {d}/standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/highres2standard1.png'
subprocess.Popen(cmd.split()).wait()
cmd = f'slicer {d}/standard {d}/highres2standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/highres2standard2.png'
subprocess.Popen(cmd.split()).wait()
cmd = (f'pngappend {d}/highres2standard1.png - '
f'{d}/highres2standard2.png {d}/highres2standard.png')
subprocess.Popen(cmd.split()).wait()
# clean up
for search in ['*1.png', '*2.png', 'sl?.png']:
imgs = glob.glob(f'{d}/{search}')
for img in imgs:
os.remove(img)
In our local project directory, we will need a folder that will serve as the 'reg' directory for FSL first-level analyses. This will contain all the files we just made plus the registration to functional space which is yet to be made. To do this, we can make a new folder containing symbolic links to the registration files:
reg_dir = f'{PROJ_DIR}/registration/sub-{subject}'
os.makedirs(reg_dir, exist_ok=True)
# put links to everything in local project reg dir
in_paths = glob.glob(f'{fnirt_dir}/*')
for path in in_paths:
outpath = f'{reg_dir}/{op.basename(path)}'
if not op.isfile(outpath):
os.system(f'ln -s {path} {outpath}')
# link to example_func
out_path = f'{reg_dir}/example_func.nii.gz'
if not op.exists(out_path):
os.system(f'ln -s {op.abspath(ref_func)} {out_path}')The next transformation is between the native anatomical and functional spaces. These are stored in the local project dir as they are unlikely to be useful for other experiments that this subject participated in. There are different options for this transformation. I use FSL tools by default, and Freesurfer tools when this fails, e.g. if the field of view in the func image is small, which FSL doesn't like.
method = 'FSL'
if method == 'FSL':
example_func2highres = f'{reg_dir}/example_func2highres.mat'
if not op.isfile(example_func2highres):
# BBR method
os.system(
f'epi_reg '
f'--epi={ref_func} '
f'--t1={ref_anat} '
f'--t1brain={ref_anat_brain} '
f'--out={example_func2highres[:-4]}')
# linear search method
# os.system(f'flirt -in {ref_anat_brain} -ref {ref_func}
# -omat {highres2example_func}')
# invert the transformation
highres2example_func = f'{reg_dir}/highres2example_func.mat'
if not op.isfile(highres2example_func):
os.system(f'convert_xfm '
f'-omat {highres2example_func} '
f'-inverse '
f'{example_func2highres}')
elif method == 'FreeSurfer':
# make lta file using bbregister
lta = f'{reg_dir}/example_func2highres.lta'
if not op.isfile(lta):
os.system(f'bbregister '
f'--s sub-{subject} '
f'--mov {ref_func} '
f'--init-fsl '
f'--lta {lta} '
f'--bold')
# convert to an FSL-friendly format
example_func2highres = f'{reg_dir}/example_func2highres.mat'
if not op.isfile(example_func2highres):
os.system(f'lta_convert '
f'--inlta {lta} '
f'--outfsl {example_func2highres} '
f'--src {ref_func} '
f'--trg {ref_anat}')
# invert the transformation
lta_inv = f'{reg_dir}/highres2example_func.lta'
if not op.isfile(lta_inv):
os.system(f'lta_convert '
f'--inlta {lta} '
f'--outlta {lta_inv} '
f'--invert')
# convert the inverted transform to an FSL-friendly format
highres2example_func = f'{reg_dir}/highres2example_func.mat'
if not op.isfile(highres2example_func):
os.system(f'lta_convert '
f'--inlta {lta_inv} '
f'--outfsl {highres2example_func} '
f'--src {ref_anat} '
f'--trg {ref_func}')Now that we can combine the two transforms we have made to convert from standard space directly to functional space. These can then be checked by making the same plots as before.
example_func2standard_warp =f'{reg_dir}/example_func2standard_warp.nii.gz'
if not op.isfile(example_func2standard_warp):
os.system(f'convertwarp '
f'--ref={ref_std_brain} '
f'--premat={example_func2highres} '
f'--warp1={highres2standard_warp} '
f'--out={example_func2standard_warp}')
standard2example_func_warp =f'{reg_dir}/standard2example_func_warp.nii.gz'
# invert the transformation
if not op.isfile(standard2example_func_warp):
os.system(f'invwarp '
f'-w {example_func2standard_warp} '
f'-o {standard2example_func_warp} '
f'-r {ref_func}')
# apply transform to reference functional image
example_func2standard_img = f'{reg_dir}/example_func2standard.nii.gz'
if not op.isfile(example_func2standard_img):
os.system(f'applywarp -i {ref_func} -r {ref_std_brain} '
f'-o {reg_dir}/example_func2standard '
f'-w {reg_dir}/example_func2standard_warp')
# make reg images
d = reg_dir
if not op.isfile(f'{d}/example_func2standard.png'):
slicer_str = (
f'-x 0.35 {d}/sla.png -x 0.45 {d}/slb.png '
f'-x 0.55 {d}/slc.png -x 0.65 {d}/sld.png '
f'-y 0.35 {d}/sle.png -y 0.45 {d}/slf.png '
f'-y 0.55 {d}/slg.png -y 0.65 {d}/slh.png '
f'-z 0.35 {d}/sli.png -z 0.45 {d}/slj.png '
f'-z 0.55 {d}/slk.png -z 0.65 {d}/sll.png')
append_str = ' + '.join([f'{d}/sl{i}.png' for i in 'abcdefghijkl'])
cmd = f'slicer {d}/example_func2standard {d}/standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/example_func2standard1.png'
subprocess.Popen(cmd.split()).wait()
cmd = f'slicer {d}/standard {d}/example_func2standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/example_func2standard2.png'
subprocess.Popen(cmd.split()).wait()
cmd = (f'pngappend {d}/example_func2standard1.png - '
f'{d}/example_func2standard2.png {d}/example_func2standard.png')
subprocess.Popen(cmd.split()).wait()
# clean up
for search in ['*1.png', '*2.png', 'sl?.png']:
imgs = glob.glob(f'{d}/{search}')
for img in imgs:
os.remove(img)
# make a fake reg dir that tricks fsl into doing higher-level analyses in
# native func space
no_reg_dir = f'{PROJ_DIR}/derivatives/registration/sub-{subject}_no-reg'
if not op.isdir(no_reg_dir):
shutil.copytree(reg_dir, no_reg_dir)
for reg_mat in glob.glob(f'{no_reg_dir}/*.mat'):
os.remove(reg_mat)
shutil.copy(f'{os.environ["FSLDIR"]}/etc/flirtsch/ident.mat',
reg_mat)
for reg_warp in glob.glob(f'{no_reg_dir}/*warp.nii.gz'):
os.remove(reg_warp)
for space in ['standard', 'highres']:
os.remove(f'{no_reg_dir}/{space}.nii.gz')
shutil.copy(f'{no_reg_dir}/example_func.nii.gz',
f'{no_reg_dir}/{space}.nii.gz')Again, we can make the registration plots to check the final transformation.
d = reg_dir
if not op.isfile(f'{d}/example_func2standard.png'):
slicer_str = (
f'-x 0.35 {d}/sla.png -x 0.45 {d}/slb.png '
f'-x 0.55 {d}/slc.png -x 0.65 {d}/sld.png '
f'-y 0.35 {d}/sle.png -y 0.45 {d}/slf.png '
f'-y 0.55 {d}/slg.png -y 0.65 {d}/slh.png '
f'-z 0.35 {d}/sli.png -z 0.45 {d}/slj.png '
f'-z 0.55 {d}/slk.png -z 0.65 {d}/sll.png')
append_str = ' + '.join([f'{d}/sl{i}.png' for i in 'abcdefghijkl'])
cmd = f'slicer {d}/example_func2standard {d}/standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/example_func2standard1.png'
subprocess.Popen(cmd.split()).wait()
cmd = f'slicer {d}/standard {d}/example_func2standard -s 2 {slicer_str}'
subprocess.Popen(cmd.split()).wait()
cmd = f'pngappend {append_str} {d}/example_func2standard2.png'
subprocess.Popen(cmd.split()).wait()
cmd = (f'pngappend {d}/example_func2standard1.png - '
f'{d}/example_func2standard2.png {d}/example_func2standard.png')
subprocess.Popen(cmd.split()).wait()
# clean up
for search in ['*1.png', '*2.png', 'sl?.png']:
imgs = glob.glob(f'{d}/{search}')
for img in imgs:
os.remove(img)To use this directory in FSL FEAT, you need to configure the first-level analyses to skip registration and post-stats. Run these analyses, then make a link to the registration directory inside each first-level output directory.
os.system(f'ln -s {op.abspath(reg_dir)} /path/to/first-level.feat')Higher-level analyses should now run as normal, with outputs in standard space. If you were interested in first-level post-stats, these will need to be run at this point. Finally, if you wanted to run higher-level analyses in native functional space, you can make a 'fake' reg directory, which tricks FEAT into thinking it is converting everything into standard space when really it all stays in native functional space.
no_reg_dir = f'{PROJ_DIR}/derivatives/registration/sub-{subject}_no-reg'
if not op.isdir(no_reg_dir):
shutil.copytree(reg_dir, no_reg_dir)
for reg_mat in glob.glob(f'{no_reg_dir}/*.mat'):
os.remove(reg_mat)
shutil.copy(f'{os.environ["FSLDIR"]}/etc/flirtsch/ident.mat',
reg_mat)
for reg_warp in glob.glob(f'{no_reg_dir}/*warp.nii.gz'):
os.remove(reg_warp)
for space in ['standard', 'highres']:
os.remove(f'{no_reg_dir}/{space}.nii.gz')
shutil.copy(f'{no_reg_dir}/example_func.nii.gz',
f'{no_reg_dir}/{space}.nii.gz')
os.system(f'ln -s {op.abspath(no_reg_dir)} /path/to/first-level.feat')Note that higher-level analyses across different subjects will not work if they haven't been registered to the same standard space. That's it!
Dave
Comentarios