From 860d4115d38fb151624e208a95cabecbeb0b0e40 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 1 Sep 2017 10:36:27 -0700 Subject: [PATCH 01/17] ENH: Break mri_coreg out from bbregister --- fmriprep/workflows/bold.py | 2 +- fmriprep/workflows/util.py | 19 ++++++++++++++----- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/fmriprep/workflows/bold.py b/fmriprep/workflows/bold.py index c1d9620fd..71c51b7d1 100755 --- a/fmriprep/workflows/bold.py +++ b/fmriprep/workflows/bold.py @@ -899,7 +899,7 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, ) if freesurfer: - bbr_wf = init_bbreg_wf(bold2t1w_dof) + bbr_wf = init_bbreg_wf(bold2t1w_dof, omp_nthreads=omp_nthreads) else: bbr_wf = init_fsl_bbr_wf(bold2t1w_dof) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 0dbee0917..6e1c94121 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -18,7 +18,7 @@ from niworkflows.nipype.interfaces import utility as niu from niworkflows.interfaces.utils import CopyXForm from niworkflows.nipype.interfaces import fsl, afni, c3, ants, freesurfer as fs -from niworkflows.interfaces.registration import FLIRTRPT, BBRegisterRPT +from niworkflows.interfaces.registration import FLIRTRPT, BBRegisterRPT, MRICoregRPT from niworkflows.interfaces.masks import SimpleShowMaskRPT from ..interfaces.images import extract_wm @@ -174,7 +174,7 @@ def init_skullstrip_bold_wf(name='skullstrip_bold_wf'): return workflow -def init_bbreg_wf(bold2t1w_dof, name='bbreg_wf'): +def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): """ This workflow uses FreeSurfer's ``bbregister`` to register a BOLD image to a T1-weighted structural image. @@ -187,7 +187,7 @@ def init_bbreg_wf(bold2t1w_dof, name='bbreg_wf'): :simple_form: yes from fmriprep.workflows.util import init_bbreg_wf - wf = init_bbreg_wf(bold2t1w_dof=9) + wf = init_bbreg_wf(bold2t1w_dof=9, omp_nthreads=1) Parameters @@ -236,9 +236,14 @@ def init_bbreg_wf(bold2t1w_dof, name='bbreg_wf'): niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report']), name='outputnode') + mri_coreg= pe.Node( + MRICoregRPT(dof=bold2t1w_dof, sep=[4], ftol=0.0001, linmintol=0.01, + num_threads=omp_nthreads, generate_report=True), + name='mri_coreg', n_procs=omp_nthreads) + bbregister = pe.Node( - BBRegisterRPT(dof=bold2t1w_dof, contrast_type='t2', init='coreg', - registered_file=True, out_lta_file=True, generate_report=True), + BBRegisterRPT(dof=bold2t1w_dof, contrast_type='t2', registered_file=True, + out_lta_file=True, generate_report=True), name='bbregister') lta_concat = pe.Node(fs.ConcatenateLTA(out_file='out.lta'), name='lta_concat') @@ -252,9 +257,13 @@ def init_bbreg_wf(bold2t1w_dof, name='bbreg_wf'): name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB) workflow.connect([ + (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id'), + ('in_file', 'source_file')]), (inputnode, bbregister, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), ('in_file', 'source_file')]), + (mri_coreg, bbregister, [('out_lta_file', 'init_reg_file')]), (bbregister, outputnode, [('out_report', 'out_report')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), From 291d2ac49d1d42c91b2a2ca846255af295c91623 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 7 Sep 2017 15:32:19 -0400 Subject: [PATCH 02/17] ENH: Enable first-pass reporting for FLIRT --- fmriprep/workflows/util.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 6e1c94121..ea5e1e715 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -345,11 +345,10 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): name='outputnode') wm_mask = pe.Node(niu.Function(function=extract_wm), name='wm_mask') - flt_bbr_init = pe.Node(FLIRTRPT(dof=6, generate_report=False), name='flt_bbr_init') - flt_bbr = pe.Node(FLIRTRPT(cost_func='bbr', dof=bold2t1w_dof, save_log=True, - generate_report=True), name='flt_bbr') - flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'), - 'etc/flirtsch/bbr.sch') + flt_bbr_init = pe.Node(FLIRTRPT(dof=6, generate_report=True), name='flt_bbr_init') + flt_bbr = pe.Node(FLIRTRPT(cost_func='bbr', dof=bold2t1w_dof, generate_report=True), + name='flt_bbr') + flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'), 'etc/flirtsch/bbr.sch') # make equivalent warp fields invt_bbr = pe.Node(fsl.ConvertXFM(invert_xfm=True), name='invt_bbr', From 15b512d83f6b877041b05cac5f494dbd2f91c919 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 8 Sep 2017 13:44:43 -0400 Subject: [PATCH 03/17] ENH: Compare transforms before/after bbregister --- fmriprep/workflows/util.py | 49 ++++++++++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 5 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index ea5e1e715..6b333d783 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -26,6 +26,27 @@ DEFAULT_MEMORY_MIN_GB = 0.01 +def compare_xforms(fallback_mat, test_mat): + import numpy as np + from transforms3d.affines import decompose44 + from transforms3d.axangles import mat2axangle + # [[R], [A], [S], [1]] = mat1.dot([[x], [y], [z], [1]]) + mat1 = np.loadtxt(fallback_mat) + # [[R'], [A'], [S'], [1]] = mat2.dot([[x], [y], [z], [1]]) + mat2 = np.loadtxt(test_mat) + # [[R'], [A'], [S'], [1]] = mat2.dot(mat1_inv.dot([[R], [A], [S], [1]])) + comp = mat2.dot(np.linalg.pinv(mat1)) + trans, rotation_matrix, scales, shears = decompose44(comp) + + max_trans = np.max(np.abs(trans)) + rot = mat2axangle(rotation_matrix)[1] + max_scale = np.max(np.abs(scales)) + + fallback = any((max_trans > 2, rot > np.pi / 36, max_scale > 1.1)) + + return 0 if fallback else 1 + + def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf', omp_nthreads=1): """ @@ -236,7 +257,7 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report']), name='outputnode') - mri_coreg= pe.Node( + mri_coreg = pe.Node( MRICoregRPT(dof=bold2t1w_dof, sep=[4], ftol=0.0001, linmintol=0.01, num_threads=omp_nthreads, generate_report=True), name='mri_coreg', n_procs=omp_nthreads) @@ -256,6 +277,14 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): fsl2itk_inv = pe.Node(c3.C3dAffineTool(fsl2ras=True, itk_transform=True), name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB) + transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') + reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') + + compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') + + select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') + select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') + workflow.connect([ (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), @@ -264,11 +293,7 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): ('subject_id', 'subject_id'), ('in_file', 'source_file')]), (mri_coreg, bbregister, [('out_lta_file', 'init_reg_file')]), - (bbregister, outputnode, [('out_report', 'out_report')]), - (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), - (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), (inputnode, lta_concat, [('t1_2_fsnative_reverse_transform', 'in_lta2')]), - (bbregister, lta_concat, [('out_lta_file', 'in_lta1')]), (lta_concat, lta2fsl_fwd, [('out_file', 'in_lta')]), (lta_concat, lta2fsl_inv, [('out_file', 'in_lta')]), (inputnode, fsl2itk_fwd, [('t1_brain', 'reference_file'), @@ -277,6 +302,20 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): ('t1_brain', 'source_file')]), (lta2fsl_fwd, fsl2itk_fwd, [('out_fsl', 'transform_file')]), (lta2fsl_inv, fsl2itk_inv, [('out_fsl', 'transform_file')]), + (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), + (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), + (mri_coreg, transforms, [('out_lta_file', 'in1')]), + (bbregister, transforms, [('out_lta_file', 'in2')]), + (mri_coreg, compare_transforms, [('out_lta_file', 'fallback_mat')]), + (bbregister, compare_transforms, [('out_lta_file', 'test_mat')]), + (transforms, select_transform, [('out', 'inlist')]), + (compare_transforms, select_transform, [('out', 'index')]), + (select_transform, lta_concat, [('out', 'in_lta1')]), + (mri_coreg, reports, [('out_report', 'in1')]), + (bbregister, reports, [('out_report', 'in2')]), + (reports, select_report, [('out', 'inlist')]), + (compare_transforms, select_report, [('out', 'index')]), + (select_report, outputnode, [('out', 'out_report')]), ]) return workflow From 3839cde787d09ffabd74d28adea35a1b4970dabe Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 18 Sep 2017 14:41:23 -0400 Subject: [PATCH 04/17] ENH: Compare transforms before/after FLIRT BBR --- fmriprep/workflows/util.py | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 6b333d783..3fdadeec2 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -296,6 +296,7 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (inputnode, lta_concat, [('t1_2_fsnative_reverse_transform', 'in_lta2')]), (lta_concat, lta2fsl_fwd, [('out_file', 'in_lta')]), (lta_concat, lta2fsl_inv, [('out_file', 'in_lta')]), + # Output ITK transforms (inputnode, fsl2itk_fwd, [('t1_brain', 'reference_file'), ('in_file', 'source_file')]), (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), @@ -306,8 +307,10 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), (mri_coreg, transforms, [('out_lta_file', 'in1')]), (bbregister, transforms, [('out_lta_file', 'in2')]), + # Compare transforms (mri_coreg, compare_transforms, [('out_lta_file', 'fallback_mat')]), (bbregister, compare_transforms, [('out_lta_file', 'test_mat')]), + # Select output transform (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, lta_concat, [('out', 'in_lta1')]), @@ -389,7 +392,6 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): name='flt_bbr') flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'), 'etc/flirtsch/bbr.sch') - # make equivalent warp fields invt_bbr = pe.Node(fsl.ConvertXFM(invert_xfm=True), name='invt_bbr', mem_gb=DEFAULT_MEMORY_MIN_GB) @@ -400,6 +402,14 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): fsl2itk_inv = pe.Node(c3.C3dAffineTool(fsl2ras=True, itk_transform=True), name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB) + transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') + reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') + + compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') + + select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') + select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') + workflow.connect([ (inputnode, wm_mask, [('t1_seg', 'in_seg')]), (inputnode, flt_bbr_init, [('in_file', 'in_file'), @@ -412,12 +422,23 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): ('in_file', 'source_file')]), (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), ('t1_brain', 'source_file')]), - (flt_bbr, invt_bbr, [('out_matrix_file', 'in_file')]), - (flt_bbr, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]), + # Compare transforms + (flt_bbr_init, compare_transforms, [('out_matrix_file', 'fallback_mat')]), + (flt_bbr, compare_transforms, [('out_matrix_file', 'test_mat')]), + # Select output transform + (flt_bbr_init, transforms, [('out_matrix_file', 'in1')]), + (flt_bbr, transforms, [('out_matrix_file', 'in2')]), + (transforms, select_transform, [('out', 'inlist')]), + (compare_transforms, select_transform, [('out', 'index')]), + (select_transform, invt_bbr, [('out', 'in_file')]), + (select_transform, fsl2itk_fwd, [('out', 'transform_file')]), (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), - (flt_bbr, outputnode, [('out_report', 'out_report')]), - ]) + (flt_bbr_init, reports, [('out_report', 'in1')]), + (flt_bbr, reports, [('out_report', 'in2')]), + (reports, select_report, [('out', 'inlist')]), + (compare_transforms, select_report, [('out', 'index')]), + (select_report, outputnode, [('out', 'out_report')])]) return workflow From b2979149325b23ba418178cbc3db8bd1cd61b439 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 18 Sep 2017 16:21:13 -0400 Subject: [PATCH 05/17] REPORT: Update report names and summary to reflect registration used --- fmriprep/interfaces/reports.py | 17 ++++++++++----- fmriprep/viz/config.json | 16 ++++++++++++-- fmriprep/workflows/bold.py | 40 +++++++++++++++++++++++----------- fmriprep/workflows/util.py | 40 +++++++++++++++++++--------------- 4 files changed, 75 insertions(+), 38 deletions(-) diff --git a/fmriprep/interfaces/reports.py b/fmriprep/interfaces/reports.py index 19c3b3cc6..8eb9f0266 100644 --- a/fmriprep/interfaces/reports.py +++ b/fmriprep/interfaces/reports.py @@ -143,8 +143,11 @@ class FunctionalSummaryInputSpec(BaseInterfaceInputSpec): mandatory=True) pe_direction = traits.Enum(None, 'i', 'i-', 'j', 'j-', mandatory=True, desc='Phase-encoding direction detected') - registration = traits.Enum('FLIRT', 'bbregister', mandatory=True, + registration = traits.Enum('FSL', 'FreeSurfer', mandatory=True, desc='Functional/anatomical registration method') + fallback = traits.Bool(desc='Boundary-based registration rejected') + registration_dof = traits.Enum(6, 9, 12, desc='Registration degrees of freedom', + mandatory=True) output_spaces = traits.List(desc='Target spaces') confounds = traits.List(desc='Confounds collected') @@ -153,18 +156,22 @@ class FunctionalSummary(SummaryInterface): input_spec = FunctionalSummaryInputSpec def _generate_segment(self): + dof = self.inputs.registration_dof stc = {True: 'Applied', False: 'Not applied', 'TooShort': 'Skipped (too few volumes)'}[self.inputs.slice_timing] - stc = "Applied" if self.inputs.slice_timing else "Not applied" sdc = {'epi': 'Phase-encoding polarity (pepolar)', 'fieldmap': 'Direct fieldmapping', 'phasediff': 'Phase difference', 'SyN': 'Symmetric normalization (SyN) - no fieldmaps', 'None': 'None'}[self.inputs.distortion_correction] - reg = {'FLIRT': 'FLIRT with boundary-based registration (BBR) metric', - 'bbregister': 'FreeSurfer boundary-based registration (bbregister)' - }[self.inputs.registration] + reg = {'FSL': [ + 'FLIRT with boundary-based registration (BBR) metric - %d dof' % dof, + 'FLIRT rigid registration - 6 dof'], + 'FreeSurfer': [ + 'FreeSurfer boundary-based registration (bbregister) - %d dof' % dof, + 'FreeSurfer mri_coreg - %d dof' % dof], + }[self.inputs.registration][self.inputs.fallback] if self.inputs.pe_direction is None: pedir = 'MISSING - Assuming Anterior-Posterior' else: diff --git a/fmriprep/viz/config.json b/fmriprep/viz/config.json index dd2f21136..21e954f57 100644 --- a/fmriprep/viz/config.json +++ b/fmriprep/viz/config.json @@ -104,17 +104,29 @@ "title" : "Skull stripped EPI", "description": "Skull stripped EPI" }, + { + "name": "epi_mean_t1_registration/flirt", + "file_pattern": "func/.*_bold_flirt", + "title": "EPI to T1 registration", + "description": "FLIRT was used to generate transformations from EPI space to T1 Space - BBR refinement rejected" + }, + { + "name": "epi_mean_t1_registration/coreg", + "file_pattern": "func/.*_bold_coreg", + "title": "EPI to T1 registration", + "description": "mri_coreg (FreeSurfer) was used to generate transformations from EPI space to T1 Space - bbregister refinement rejected" + }, { "name": "epi_mean_t1_registration/flt_bbr", "file_pattern": "func/.*_bold_flt_bbr", "title": "EPI to T1 registration", - "description": "EPI was used to generate transformations from EPI space to T1 Space - FAST segmentation used for BBR" + "description": "FLIRT was used to generate transformations from EPI space to T1 Space - FAST segmentation used for BBR" }, { "name": "epi_mean_t1_registration/bbr", "file_pattern": "func/.*_bold_bbr", "title": "EPI to T1 registration", - "description": "EPI was used to generate transformations from EPI space to T1 Space" + "description": "bbregister was used to generate transformations from EPI space to T1 Space" }, { "name": "ica_aroma", diff --git a/fmriprep/workflows/bold.py b/fmriprep/workflows/bold.py index 71c51b7d1..283f8b000 100755 --- a/fmriprep/workflows/bold.py +++ b/fmriprep/workflows/bold.py @@ -258,11 +258,13 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, 'aroma_noise_ics', 'melodic_mix', 'nonaggr_denoised_file']), name='outputnode') - summary = pe.Node(FunctionalSummary(output_spaces=output_spaces, - pe_direction=bold_pe), name='summary', - mem_gb=0.05) - summary.inputs.slice_timing = run_stc - summary.inputs.registration = 'bbregister' if freesurfer else 'FLIRT' + summary = pe.Node( + FunctionalSummary(output_spaces=output_spaces, + slice_timing=run_stc, + registration='FreeSurfer' if freesurfer else 'FSL', + registration_dof=bold2t1w_dof, + pe_direction=bold_pe), + name='summary', mem_gb=DEFAULT_MEMORY_MIN_GB, run_without_submitting=True) func_reports_wf = init_func_reports_wf(reportlets_dir=reportlets_dir, freesurfer=freesurfer, @@ -348,6 +350,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, ('outputnode.validation_report', 'inputnode.validation_report')]), (bold_reg_wf, func_reports_wf, [ ('outputnode.out_report', 'inputnode.bold_reg_report'), + ('outputnode.fallback', 'inputnode.bold_reg_fallback'), ]), (bold_confounds_wf, outputnode, [ ('outputnode.confounds_file', 'confounds'), @@ -362,6 +365,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, ('outputnode.tcompcor_report', 'inputnode.tcompcor_report'), ('outputnode.ica_aroma_report', 'inputnode.ica_aroma_report')]), (bold_confounds_wf, summary, [('outputnode.confounds_list', 'confounds')]), + (bold_reg_wf, summary, [('outputnode.fallback', 'fallback')]), (summary, func_reports_wf, [('out_report', 'inputnode.summary_report')]), ]) @@ -880,6 +884,8 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, BOLD mask in T1 space out_report Reportlet visualizing quality of registration + fallback + Boolean indicating whether BBR was rejected (mri_coreg registration returned) """ workflow = pe.Workflow(name=name) @@ -894,7 +900,7 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, outputnode = pe.Node( niu.IdentityInterface(fields=['itk_bold_to_t1', 'itk_t1_to_bold', 'bold_t1', 'bold_mask_t1', - 'out_report']), + 'out_report', 'fallback']), name='outputnode' ) @@ -913,7 +919,8 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, ('t1_brain', 'inputnode.t1_brain')]), (bbr_wf, outputnode, [('outputnode.itk_bold_to_t1', 'itk_bold_to_t1'), ('outputnode.itk_t1_to_bold', 'itk_t1_to_bold'), - ('outputnode.out_report', 'out_report')]), + ('outputnode.out_report', 'out_report'), + ('outputnode.fallback', 'fallback')]), ]) gen_ref = pe.Node(GenerateSamplingReference(), name='gen_ref', @@ -1495,8 +1502,8 @@ def init_func_reports_wf(reportlets_dir, freesurfer, use_aroma, use_syn, name='f inputnode = pe.Node( niu.IdentityInterface( fields=['source_file', 'summary_report', 'validation_report', 'bold_mask_report', - 'bold_reg_report', 'acompcor_report', 'tcompcor_report', 'syn_sdc_report', - 'ica_aroma_report']), + 'bold_reg_report', 'bold_reg_fallback', 'acompcor_report', 'tcompcor_report', + 'syn_sdc_report', 'ica_aroma_report']), name='inputnode') ds_summary_report = pe.Node( @@ -1523,9 +1530,14 @@ def init_func_reports_wf(reportlets_dir, freesurfer, use_aroma, use_syn, name='f name='ds_syn_sdc_report', run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB) + def _bold_reg_suffix(fallback, freesurfer): + if fallback: + return 'coreg' if freesurfer else 'flirt' + else: + return 'bbr' if freesurfer else 'flt_bbr' + ds_bold_reg_report = pe.Node( - DerivativesDataSink(base_directory=reportlets_dir, - suffix='bbr' if freesurfer else 'flt_bbr'), + DerivativesDataSink(base_directory=reportlets_dir), name='ds_bold_reg_report', run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB) @@ -1554,8 +1566,10 @@ def init_func_reports_wf(reportlets_dir, freesurfer, use_aroma, use_syn, name='f ('validation_report', 'in_file')]), (inputnode, ds_bold_mask_report, [('source_file', 'source_file'), ('bold_mask_report', 'in_file')]), - (inputnode, ds_bold_reg_report, [('source_file', 'source_file'), - ('bold_reg_report', 'in_file')]), + (inputnode, ds_bold_reg_report, [ + ('source_file', 'source_file'), + ('bold_reg_report', 'in_file'), + (('bold_reg_fallback', _bold_reg_suffix, freesurfer), 'suffix')]), (inputnode, ds_acompcor_report, [('source_file', 'source_file'), ('acompcor_report', 'in_file')]), (inputnode, ds_tcompcor_report, [('source_file', 'source_file'), diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 3fdadeec2..e6fef6964 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -26,7 +26,7 @@ DEFAULT_MEMORY_MIN_GB = 0.01 -def compare_xforms(fallback_mat, test_mat): +def compare_xforms(test_mat, fallback_mat): import numpy as np from transforms3d.affines import decompose44 from transforms3d.axangles import mat2axangle @@ -42,9 +42,7 @@ def compare_xforms(fallback_mat, test_mat): rot = mat2axangle(rotation_matrix)[1] max_scale = np.max(np.abs(scales)) - fallback = any((max_trans > 2, rot > np.pi / 36, max_scale > 1.1)) - - return 0 if fallback else 1 + return any((max_trans > 2, rot > np.pi / 36, max_scale > 1.1)) def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf', @@ -242,7 +240,9 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): itk_t1_to_bold Affine transform from T1 space to BOLD space (ITK format) out_report - reportlet for assessing registration quality + Reportlet for assessing registration quality + fallback + Boolean indicating whether BBR was rejected (mri_coreg registration returned) """ workflow = pe.Workflow(name=name) @@ -254,7 +254,7 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): 't1_seg', 't1_brain']), # FLIRT BBR name='inputnode') outputnode = pe.Node( - niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report']), + niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report', 'fallback']), name='outputnode') mri_coreg = pe.Node( @@ -305,17 +305,18 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (lta2fsl_inv, fsl2itk_inv, [('out_fsl', 'transform_file')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), - (mri_coreg, transforms, [('out_lta_file', 'in1')]), - (bbregister, transforms, [('out_lta_file', 'in2')]), + (bbregister, transforms, [('out_lta_file', 'in1')]), + (mri_coreg, transforms, [('out_lta_file', 'in2')]), # Compare transforms - (mri_coreg, compare_transforms, [('out_lta_file', 'fallback_mat')]), (bbregister, compare_transforms, [('out_lta_file', 'test_mat')]), + (mri_coreg, compare_transforms, [('out_lta_file', 'fallback_mat')]), + (compare_transforms, outputnode, [('out', 'fallback')]), # Select output transform (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, lta_concat, [('out', 'in_lta1')]), - (mri_coreg, reports, [('out_report', 'in1')]), - (bbregister, reports, [('out_report', 'in2')]), + (bbregister, reports, [('out_report', 'in1')]), + (mri_coreg, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')]), @@ -371,7 +372,9 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): itk_t1_to_bold Affine transform from T1 space to BOLD space (ITK format) out_report - reportlet for assessing registration quality + Reportlet for assessing registration quality + fallback + Boolean indicating whether BBR was rejected (rigid FLIRT registration returned) """ workflow = pe.Workflow(name=name) @@ -383,7 +386,7 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): 't1_seg', 't1_brain']), # FLIRT BBR name='inputnode') outputnode = pe.Node( - niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report']), + niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'out_report', 'fallback']), name='outputnode') wm_mask = pe.Node(niu.Function(function=extract_wm), name='wm_mask') @@ -423,11 +426,12 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), ('t1_brain', 'source_file')]), # Compare transforms - (flt_bbr_init, compare_transforms, [('out_matrix_file', 'fallback_mat')]), (flt_bbr, compare_transforms, [('out_matrix_file', 'test_mat')]), + (flt_bbr_init, compare_transforms, [('out_matrix_file', 'fallback_mat')]), + (compare_transforms, outputnode, [('out', 'fallback')]), # Select output transform - (flt_bbr_init, transforms, [('out_matrix_file', 'in1')]), - (flt_bbr, transforms, [('out_matrix_file', 'in2')]), + (flt_bbr, transforms, [('out_matrix_file', 'in1')]), + (flt_bbr_init, transforms, [('out_matrix_file', 'in2')]), (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, invt_bbr, [('out', 'in_file')]), @@ -435,8 +439,8 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), - (flt_bbr_init, reports, [('out_report', 'in1')]), - (flt_bbr, reports, [('out_report', 'in2')]), + (flt_bbr, reports, [('out_report', 'in1')]), + (flt_bbr_init, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')])]) From 871679966a3de3f4e1c433c551983212bd037f21 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 19 Sep 2017 14:16:56 -0400 Subject: [PATCH 06/17] RF: Update thresholds --- fmriprep/workflows/util.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index e6fef6964..3785e1bb3 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -38,11 +38,15 @@ def compare_xforms(test_mat, fallback_mat): comp = mat2.dot(np.linalg.pinv(mat1)) trans, rotation_matrix, scales, shears = decompose44(comp) - max_trans = np.max(np.abs(trans)) + shift_thresh = 5 # mm + rot_thresh = np.pi / 36 # 5 degrees + scale_thresh = 1.1 # scale factor + + shift_magnitude = np.sqrt(trans.dot(trans)) # 2-norm rot = mat2axangle(rotation_matrix)[1] max_scale = np.max(np.abs(scales)) - return any((max_trans > 2, rot > np.pi / 36, max_scale > 1.1)) + return shift_magnitude > shift_thresh or rot > rot_thresh or max_scale > scale_thresh def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf', From 0509970fb9ca7c4a4c78b44c95c0771caf50517b Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 19 Sep 2017 16:40:15 -0400 Subject: [PATCH 07/17] DOC: Use simple graph for bold_reg [docs only] --- docs/workflows.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index ca6c461b8..94c9ddd38 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -296,7 +296,7 @@ EPI to T1w registration :mod:`fmriprep.workflows.bold.init_bold_reg_wf` .. workflow:: - :graph2use: colored + :graph2use: orig :simple_form: yes from fmriprep.workflows.bold import init_bold_reg_wf From 3e6cc2f974c7bc440eda9490c339ee3c001c8924 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 19 Sep 2017 16:53:39 -0400 Subject: [PATCH 08/17] LOG: Print transform difference criteria --- fmriprep/workflows/util.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 3785e1bb3..6291b1c6b 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -46,6 +46,9 @@ def compare_xforms(test_mat, fallback_mat): rot = mat2axangle(rotation_matrix)[1] max_scale = np.max(np.abs(scales)) + print("Shift: {:.1g}mm\nRotation: {:.1g}°\nScale: {:.2g}".format( + shift_magnitude, rot * 180 / np.pi, max_scale)) + return shift_magnitude > shift_thresh or rot > rot_thresh or max_scale > scale_thresh From 513b550fe80c664404be9ac9b4e74946e85fca5c Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 19 Sep 2017 17:01:08 -0400 Subject: [PATCH 09/17] CI: Preserve .mat files for out-of-band comparison --- .circle/tests.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.circle/tests.sh b/.circle/tests.sh index be6039466..f5963f2b2 100755 --- a/.circle/tests.sh +++ b/.circle/tests.sh @@ -51,11 +51,11 @@ case ${CIRCLE_NODE_INDEX} in set -e [[ "$RET" -eq "1" ]] # Clean up - find ~/ds054/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -type f -delete + find ~/ds054/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -type f -delete rm -r $HOME/ds054/out/fmriprep/sub-100185/log ;; 1) fmriprep-docker -i poldracklab/fmriprep:latest --config $HOME/nipype.cfg -w $HOME/ds005/scratch $HOME/data/ds005 $HOME/ds005/out participant --debug --write-graph --use-syn-sdc --use-aroma --ignore-aroma-denoising-errors - find ~/ds005/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -type f -delete + find ~/ds005/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -type f -delete ;; esac From ac12f75c387eb2c40e310d177dcc191751d25868 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 21 Sep 2017 11:23:56 -0400 Subject: [PATCH 10/17] ENH: Get LTA transforms for all registration steps --- .circle/tests.sh | 4 ++-- fmriprep/workflows/util.py | 15 ++++++++++++++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/.circle/tests.sh b/.circle/tests.sh index f5963f2b2..7feb1db2a 100755 --- a/.circle/tests.sh +++ b/.circle/tests.sh @@ -51,11 +51,11 @@ case ${CIRCLE_NODE_INDEX} in set -e [[ "$RET" -eq "1" ]] # Clean up - find ~/ds054/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -type f -delete + find ~/ds054/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -not -name "*.lta" -type f -delete rm -r $HOME/ds054/out/fmriprep/sub-100185/log ;; 1) fmriprep-docker -i poldracklab/fmriprep:latest --config $HOME/nipype.cfg -w $HOME/ds005/scratch $HOME/data/ds005 $HOME/ds005/out participant --debug --write-graph --use-syn-sdc --use-aroma --ignore-aroma-denoising-errors - find ~/ds005/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -type f -delete + find ~/ds005/scratch -not -name "*.svg" -not -name "*.html" -not -name "*.rst" -not -name "*.mat" -not -name "*.lta" -type f -delete ;; esac diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 6291b1c6b..d9d4ca626 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -292,6 +292,9 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') + lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_lta'], + name='lta_convert') + workflow.connect([ (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), @@ -327,6 +330,8 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')]), + # Convert VOX2VOX transforms to RAS2RAS transforms + (transforms, lta_convert, [('out', 'in_lta')]), ]) return workflow @@ -420,6 +425,9 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') + lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_fsl'], + name='lta_convert') + workflow.connect([ (inputnode, wm_mask, [('t1_seg', 'in_seg')]), (inputnode, flt_bbr_init, [('in_file', 'in_file'), @@ -450,6 +458,11 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): (flt_bbr_init, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), - (select_report, outputnode, [('out', 'out_report')])]) + (select_report, outputnode, [('out', 'out_report')]), + # Convert FSL transforms to LTA transforms + (transforms, lta_convert, [('out', 'in_fsl')]), + (inputnode, lta_convert, [('in_file', 'source_file'), + ('t1_brain', 'target_file')]), + ]) return workflow From 1e65e259cb2c44f72b4101686a4921b243e11958 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 25 Sep 2017 16:20:35 -0400 Subject: [PATCH 11/17] ENH: Use affine difference norm --- fmriprep/workflows/util.py | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index d9d4ca626..5de6f5bfe 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -28,28 +28,13 @@ def compare_xforms(test_mat, fallback_mat): import numpy as np - from transforms3d.affines import decompose44 - from transforms3d.axangles import mat2axangle - # [[R], [A], [S], [1]] = mat1.dot([[x], [y], [z], [1]]) + from niworkflows.nipype.algorithms.rapidart import _calc_norm_affine mat1 = np.loadtxt(fallback_mat) - # [[R'], [A'], [S'], [1]] = mat2.dot([[x], [y], [z], [1]]) mat2 = np.loadtxt(test_mat) - # [[R'], [A'], [S'], [1]] = mat2.dot(mat1_inv.dot([[R], [A], [S], [1]])) - comp = mat2.dot(np.linalg.pinv(mat1)) - trans, rotation_matrix, scales, shears = decompose44(comp) - shift_thresh = 5 # mm - rot_thresh = np.pi / 36 # 5 degrees - scale_thresh = 1.1 # scale factor + norm, _ = _calc_norm_affine([mat2, mat1], use_differences=True) - shift_magnitude = np.sqrt(trans.dot(trans)) # 2-norm - rot = mat2axangle(rotation_matrix)[1] - max_scale = np.max(np.abs(scales)) - - print("Shift: {:.1g}mm\nRotation: {:.1g}°\nScale: {:.2g}".format( - shift_magnitude, rot * 180 / np.pi, max_scale)) - - return shift_magnitude > shift_thresh or rot > rot_thresh or max_scale > scale_thresh + return norm[1] > 15 def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf', From a6d7870aa255ec7bc6d543ff5f9f05c8b544e3c8 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 26 Sep 2017 09:04:13 -0400 Subject: [PATCH 12/17] CLI: Add --force(-no)-bbr options --- docs/workflows.rst | 3 +++ fmriprep/cli/run.py | 7 +++++++ fmriprep/workflows/base.py | 16 +++++++++++++--- fmriprep/workflows/bold.py | 18 ++++++++++++++---- fmriprep/workflows/tests/test_base.py | 1 + fmriprep/workflows/util.py | 14 ++++++++++---- 6 files changed, 48 insertions(+), 11 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 94c9ddd38..620aab8e3 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -37,6 +37,7 @@ slice-timing information and no fieldmap acquisitions): low_mem=False, anat_only=False, hires=True, + use_bbr=True, bold2t1w_dof=9, fmap_bspline=False, fmap_demean=True, @@ -198,6 +199,7 @@ BOLD preprocessing medial_surface_nan=False, debug=False, low_mem=False, + use_bbr=True, bold2t1w_dof=9, fmap_bspline=True, fmap_demean=True, @@ -303,6 +305,7 @@ EPI to T1w registration wf = init_bold_reg_wf(freesurfer=True, bold_file_size_gb=3, omp_nthreads=1, + use_bbr=True, bold2t1w_dof=9) The reference EPI image of each run is aligned by the ``bbregister`` routine to the diff --git a/fmriprep/cli/run.py b/fmriprep/cli/run.py index 4a47da620..6767199f7 100755 --- a/fmriprep/cli/run.py +++ b/fmriprep/cli/run.py @@ -115,6 +115,12 @@ def get_parser(): ' - fsnative: individual subject surface\n' ' - fsaverage*: FreeSurfer average meshes' ) + g_conf.add_argument( + '--force-bbr', action='store_true', dest='use_bbr', default=None, + help='Always use boundary-based registration (no goodness-of-fit checks)') + g_conf.add_argument( + '--force-no-bbr', action='store_false', dest='use_bbr', default=None, + help='Do not use boundary-based registration (no goodness-of-fit checks)') g_conf.add_argument( '--template', required=False, action='store', choices=['MNI152NLin2009cAsym'], default='MNI152NLin2009cAsym', @@ -310,6 +316,7 @@ def create_workflow(opts): medial_surface_nan=opts.medial_surface_nan, output_grid_ref=opts.output_grid_reference, hires=opts.hires, + use_bbr=opts.use_bbr, bold2t1w_dof=opts.bold2t1w_dof, fmap_bspline=opts.fmap_bspline, fmap_demean=opts.fmap_no_demean, diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index d47cb3ef7..ec14f8554 100755 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -34,7 +34,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, ignore, debug, low_mem, anat_only, longitudinal, omp_nthreads, skull_strip_ants, skull_strip_template, work_dir, output_dir, bids_dir, freesurfer, output_spaces, template, medial_surface_nan, hires, - bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn, + use_bbr, bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn, use_aroma, ignore_aroma_err, output_grid_ref): """ This workflow organizes the execution of FMRIPREP, with a sub-workflow for @@ -68,6 +68,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, template='MNI152NLin2009cAsym', medial_surface_nan=False, hires=True, + use_bbr=True, bold2t1w_dof=9, fmap_bspline=False, fmap_demean=True, @@ -128,6 +129,9 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, Replace medial wall values with NaNs on functional GIFTI files hires : bool Enable sub-millimeter preprocessing in FreeSurfer + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration fmap_bspline : bool @@ -179,6 +183,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, template=template, medial_surface_nan=medial_surface_nan, hires=hires, + use_bbr=use_bbr, bold2t1w_dof=bold2t1w_dof, fmap_bspline=fmap_bspline, fmap_demean=fmap_demean, @@ -206,8 +211,8 @@ def init_single_subject_wf(subject_id, task_id, name, ignore, debug, low_mem, anat_only, longitudinal, omp_nthreads, skull_strip_ants, skull_strip_template, reportlets_dir, output_dir, bids_dir, freesurfer, output_spaces, template, medial_surface_nan, - hires, bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn, - output_grid_ref, use_aroma, ignore_aroma_err): + hires, use_bbr, bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, + force_syn, output_grid_ref, use_aroma, ignore_aroma_err): """ This workflow organizes the preprocessing pipeline for a single subject. It collects and reports information about the subject, and prepares @@ -243,6 +248,7 @@ def init_single_subject_wf(subject_id, task_id, name, low_mem=False, anat_only=False, hires=True, + use_bbr=True, bold2t1w_dof=9, fmap_bspline=False, fmap_demean=True, @@ -302,6 +308,9 @@ def init_single_subject_wf(subject_id, task_id, name, Replace medial wall values with NaNs on functional GIFTI files hires : bool Enable sub-millimeter preprocessing in FreeSurfer + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration fmap_bspline : bool @@ -412,6 +421,7 @@ def init_single_subject_wf(subject_id, task_id, name, layout=layout, ignore=ignore, freesurfer=freesurfer, + use_bbr=use_bbr, bold2t1w_dof=bold2t1w_dof, reportlets_dir=reportlets_dir, output_spaces=output_spaces, diff --git a/fmriprep/workflows/bold.py b/fmriprep/workflows/bold.py index 283f8b000..fc0b7b7a3 100755 --- a/fmriprep/workflows/bold.py +++ b/fmriprep/workflows/bold.py @@ -65,7 +65,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, - bold2t1w_dof, reportlets_dir, + use_bbr, bold2t1w_dof, reportlets_dir, output_spaces, template, output_dir, omp_nthreads, fmap_bspline, fmap_demean, use_syn, force_syn, use_aroma, ignore_aroma_err, medial_surface_nan, @@ -88,6 +88,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, output_spaces=['T1w', 'fsnative', 'template', 'fsaverage5'], debug=False, + use_bbr=True, bold2t1w_dof=9, fmap_bspline=True, fmap_demean=True, @@ -108,6 +109,9 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, freesurfer : bool Enable FreeSurfer functional registration (bbregister) and resampling BOLD series to FreeSurfer surface meshes. + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration reportlets_dir : str @@ -308,6 +312,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer, # mean BOLD registration to T1w bold_reg_wf = init_bold_reg_wf(name='bold_reg_wf', freesurfer=freesurfer, + use_bbr=use_bbr, bold2t1w_dof=bold2t1w_dof, bold_file_size_gb=bold_file_size_gb, omp_nthreads=omp_nthreads, @@ -800,7 +805,7 @@ def init_bold_hmc_wf(bold_file_size_gb, omp_nthreads, name='bold_hmc_wf'): return workflow -def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, +def init_bold_reg_wf(freesurfer, use_bbr, bold2t1w_dof, bold_file_size_gb, omp_nthreads, name='bold_reg_wf', use_compression=True, use_fieldwarp=False): """ @@ -821,12 +826,16 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, wf = init_bold_reg_wf(freesurfer=True, bold_file_size_gb=3, omp_nthreads=1, + use_bbr=True, bold2t1w_dof=9) Parameters freesurfer : bool Enable FreeSurfer functional registration (bbregister) + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration bold_file_size_gb : float @@ -905,9 +914,10 @@ def init_bold_reg_wf(freesurfer, bold2t1w_dof, bold_file_size_gb, omp_nthreads, ) if freesurfer: - bbr_wf = init_bbreg_wf(bold2t1w_dof, omp_nthreads=omp_nthreads) + bbr_wf = init_bbreg_wf(use_bbr=use_bbr, bold2t1w_dof=bold2t1w_dof, + omp_nthreads=omp_nthreads) else: - bbr_wf = init_fsl_bbr_wf(bold2t1w_dof) + bbr_wf = init_fsl_bbr_wf(use_bbr=use_bbr, bold2t1w_dof=bold2t1w_dof) workflow.connect([ (inputnode, bbr_wf, [ diff --git a/fmriprep/workflows/tests/test_base.py b/fmriprep/workflows/tests/test_base.py index 1b12c46c2..003e0d941 100755 --- a/fmriprep/workflows/tests/test_base.py +++ b/fmriprep/workflows/tests/test_base.py @@ -30,6 +30,7 @@ def test_single_subject_wf(self, _): template='MNI152NLin2009cAsym', medial_surface_nan=False, hires=False, + use_bbr=None, bold2t1w_dof=9, fmap_bspline=True, fmap_demean=True, diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 5de6f5bfe..b47eecbb7 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -185,7 +185,7 @@ def init_skullstrip_bold_wf(name='skullstrip_bold_wf'): return workflow -def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): +def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): """ This workflow uses FreeSurfer's ``bbregister`` to register a BOLD image to a T1-weighted structural image. @@ -198,11 +198,14 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): :simple_form: yes from fmriprep.workflows.util import init_bbreg_wf - wf = init_bbreg_wf(bold2t1w_dof=9, omp_nthreads=1) + wf = init_bbreg_wf(use_bbr=True, bold2t1w_dof=9, omp_nthreads=1) Parameters + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration name : str, optional @@ -322,7 +325,7 @@ def init_bbreg_wf(bold2t1w_dof, omp_nthreads, name='bbreg_wf'): return workflow -def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): +def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): """ This workflow uses FSL FLIRT to register a BOLD image to a T1-weighted structural image, using a boundary-based registration (BBR) cost function. @@ -335,11 +338,14 @@ def init_fsl_bbr_wf(bold2t1w_dof, name='fsl_bbr_wf'): :simple_form: yes from fmriprep.workflows.util import init_fsl_bbr_wf - wf = init_fsl_bbr_wf(bold2t1w_dof=9) + wf = init_fsl_bbr_wf(use_bbr=True, bold2t1w_dof=9) Parameters + use_bbr : bool or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. bold2t1w_dof : 6, 9 or 12 Degrees-of-freedom for BOLD-T1w registration name : str, optional From 109dba03ac4621db456679fd56d02ca33877d220 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 26 Sep 2017 09:28:53 -0400 Subject: [PATCH 13/17] RF: Add use_bbr switching to fsl_bbr_wf --- fmriprep/workflows/util.py | 66 +++++++++++++++++++++++++++----------- 1 file changed, 48 insertions(+), 18 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index b47eecbb7..baf6d04c0 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -393,10 +393,7 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): name='outputnode') wm_mask = pe.Node(niu.Function(function=extract_wm), name='wm_mask') - flt_bbr_init = pe.Node(FLIRTRPT(dof=6, generate_report=True), name='flt_bbr_init') - flt_bbr = pe.Node(FLIRTRPT(cost_func='bbr', dof=bold2t1w_dof, generate_report=True), - name='flt_bbr') - flt_bbr.inputs.schedule = op.join(os.getenv('FSLDIR'), 'etc/flirtsch/bbr.sch') + flt_bbr_init = pe.Node(FLIRTRPT(dof=6, generate_report=not use_bbr), name='flt_bbr_init') invt_bbr = pe.Node(fsl.ConvertXFM(invert_xfm=True), name='invt_bbr', mem_gb=DEFAULT_MEMORY_MIN_GB) @@ -408,6 +405,53 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): fsl2itk_inv = pe.Node(c3.C3dAffineTool(fsl2ras=True, itk_transform=True), name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB) + workflow.connect([ + (inputnode, flt_bbr_init, [('in_file', 'in_file'), + ('t1_brain', 'reference')]), + (inputnode, fsl2itk_fwd, [('t1_brain', 'reference_file'), + ('in_file', 'source_file')]), + (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), + ('t1_brain', 'source_file')]), + (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]), + (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), + (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), + ]) + + # Short-circuit workflow building, use rigid registration + if use_bbr is False: + workflow.connect([ + (flt_bbr_init, invt_bbr, [('out_matrix_file', 'in_file')]), + (flt_bbr_init, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]), + (flt_bbr_init, outputnode, [('out_report', 'out_report')]), + ]) + outputnode.inputs.fallback = True + + return workflow + + flt_bbr = pe.Node( + FLIRTRPT(cost_func='bbr', dof=bold2t1w_dof, generate_report=True, + schedule=op.join(os.getenv('FSLDIR'), 'etc/flirtsch/bbr.sch')), + name='flt_bbr') + + workflow.connect([ + (inputnode, wm_mask, [('t1_seg', 'in_seg')]), + (inputnode, flt_bbr, [('in_file', 'in_file'), + ('t1_brain', 'reference')]), + (flt_bbr_init, flt_bbr, [('out_matrix_file', 'in_matrix_file')]), + (wm_mask, flt_bbr, [('out', 'wm_seg')]), + ]) + + # Short-circuit workflow building, use boundary-based registration + if use_bbr is True: + workflow.connect([ + (flt_bbr, invt_bbr, [('out_matrix_file', 'in_file')]), + (flt_bbr, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]), + (flt_bbr, outputnode, [('out_report', 'out_report')]), + ]) + outputnode.inputs.fallback = False + + return workflow + transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') @@ -420,17 +464,6 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): name='lta_convert') workflow.connect([ - (inputnode, wm_mask, [('t1_seg', 'in_seg')]), - (inputnode, flt_bbr_init, [('in_file', 'in_file'), - ('t1_brain', 'reference')]), - (flt_bbr_init, flt_bbr, [('out_matrix_file', 'in_matrix_file')]), - (inputnode, flt_bbr, [('in_file', 'in_file'), - ('t1_brain', 'reference')]), - (wm_mask, flt_bbr, [('out', 'wm_seg')]), - (inputnode, fsl2itk_fwd, [('t1_brain', 'reference_file'), - ('in_file', 'source_file')]), - (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), - ('t1_brain', 'source_file')]), # Compare transforms (flt_bbr, compare_transforms, [('out_matrix_file', 'test_mat')]), (flt_bbr_init, compare_transforms, [('out_matrix_file', 'fallback_mat')]), @@ -442,9 +475,6 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): (compare_transforms, select_transform, [('out', 'index')]), (select_transform, invt_bbr, [('out', 'in_file')]), (select_transform, fsl2itk_fwd, [('out', 'transform_file')]), - (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]), - (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), - (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), (flt_bbr, reports, [('out_report', 'in1')]), (flt_bbr_init, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), From 4497b56263ae003bf02a57e5dc98193d03f691dc Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 26 Sep 2017 09:57:48 -0400 Subject: [PATCH 14/17] RF: Add use_bbr switching to bbreg_wf --- fmriprep/workflows/util.py | 69 ++++++++++++++++++++++++++------------ 1 file changed, 47 insertions(+), 22 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index baf6d04c0..e9f6f381b 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -254,14 +254,9 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): mri_coreg = pe.Node( MRICoregRPT(dof=bold2t1w_dof, sep=[4], ftol=0.0001, linmintol=0.01, - num_threads=omp_nthreads, generate_report=True), + num_threads=omp_nthreads, generate_report=not use_bbr), name='mri_coreg', n_procs=omp_nthreads) - bbregister = pe.Node( - BBRegisterRPT(dof=bold2t1w_dof, contrast_type='t2', registered_file=True, - out_lta_file=True, generate_report=True), - name='bbregister') - lta_concat = pe.Node(fs.ConcatenateLTA(out_file='out.lta'), name='lta_concat') # XXX LTA-FSL-ITK may ultimately be able to be replaced with a straightforward # LTA-ITK transform, but right now the translation parameters are off. @@ -272,29 +267,14 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): fsl2itk_inv = pe.Node(c3.C3dAffineTool(fsl2ras=True, itk_transform=True), name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB) - transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') - reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') - - compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') - - select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') - select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') - - lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_lta'], - name='lta_convert') - workflow.connect([ (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), ('in_file', 'source_file')]), - (inputnode, bbregister, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id'), - ('in_file', 'source_file')]), - (mri_coreg, bbregister, [('out_lta_file', 'init_reg_file')]), + # Output ITK transforms (inputnode, lta_concat, [('t1_2_fsnative_reverse_transform', 'in_lta2')]), (lta_concat, lta2fsl_fwd, [('out_file', 'in_lta')]), (lta_concat, lta2fsl_inv, [('out_file', 'in_lta')]), - # Output ITK transforms (inputnode, fsl2itk_fwd, [('t1_brain', 'reference_file'), ('in_file', 'source_file')]), (inputnode, fsl2itk_inv, [('in_file', 'reference_file'), @@ -303,6 +283,50 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (lta2fsl_inv, fsl2itk_inv, [('out_fsl', 'transform_file')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), + ]) + + # Short-circuit workflow building, use initial registration + if use_bbr is False: + workflow.connect([ + (mri_coreg, outputnode, [('out_report', 'out_report')]), + (mri_coreg, lta_concat, [('out_lta_file', 'in_lta1')])]) + outputnode.inputs.fallback = True + + return workflow + + bbregister = pe.Node( + BBRegisterRPT(dof=bold2t1w_dof, contrast_type='t2', registered_file=True, + out_lta_file=True, generate_report=True), + name='bbregister') + + workflow.connect([ + (inputnode, bbregister, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id'), + ('in_file', 'source_file')]), + (mri_coreg, bbregister, [('out_lta_file', 'init_reg_file')]), + ]) + + # Short-circuit workflow building, use boundary-based registration + if use_bbr is True: + workflow.connect([ + (bbregister, outputnode, [('out_report', 'out_report')]), + (bbregister, lta_concat, [('out_lta_file', 'in_lta1')])]) + outputnode.inputs.fallback = False + + return workflow + + transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') + reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') + + compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') + + select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') + select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') + + lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_lta'], + name='lta_convert') + + workflow.connect([ (bbregister, transforms, [('out_lta_file', 'in1')]), (mri_coreg, transforms, [('out_lta_file', 'in2')]), # Compare transforms @@ -313,6 +337,7 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, lta_concat, [('out', 'in_lta1')]), + # Select output report (bbregister, reports, [('out_report', 'in1')]), (mri_coreg, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), From 8b4ec044dec22b69edb5801e21bd3dacd9b386e3 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 26 Sep 2017 12:17:12 -0400 Subject: [PATCH 15/17] RF: Switch comparisons from FSL to LTA --- fmriprep/workflows/util.py | 52 ++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index e9f6f381b..9ce2b5352 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -26,15 +26,24 @@ DEFAULT_MEMORY_MIN_GB = 0.01 -def compare_xforms(test_mat, fallback_mat): +def compare_xforms(lta_list, norm_threshold=15): import numpy as np from niworkflows.nipype.algorithms.rapidart import _calc_norm_affine - mat1 = np.loadtxt(fallback_mat) - mat2 = np.loadtxt(test_mat) - norm, _ = _calc_norm_affine([mat2, mat1], use_differences=True) + def read_lta(fname): + with open(fname, 'rb') as fobj: + for line in fobj: + if line.strip() == b'1 4 4': + break + lines = fobj.readlines()[:4] + return np.genfromtxt(lines) - return norm[1] > 15 + bbr_affine = read_lta(lta_list[0]) + fallback_affine = read_lta(lta_list[1]) + + norm, _ = _calc_norm_affine([fallback_affine, bbr_affine], use_differences=True) + + return norm[1] > norm_threshold def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf', @@ -318,20 +327,19 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') + lta_ras2ras = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_lta'], + name='lta_ras2ras') compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') - lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_lta'], - name='lta_convert') - workflow.connect([ (bbregister, transforms, [('out_lta_file', 'in1')]), (mri_coreg, transforms, [('out_lta_file', 'in2')]), - # Compare transforms - (bbregister, compare_transforms, [('out_lta_file', 'test_mat')]), - (mri_coreg, compare_transforms, [('out_lta_file', 'fallback_mat')]), + # Normalize LTA transforms to RAS2RAS (inputs are VOX2VOX) and compare + (transforms, lta_ras2ras, [('out', 'in_lta')]), + (lta_ras2ras, compare_transforms, [('out_lta', 'lta_list')]), (compare_transforms, outputnode, [('out', 'fallback')]), # Select output transform (transforms, select_transform, [('out', 'inlist')]), @@ -343,8 +351,6 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')]), - # Convert VOX2VOX transforms to RAS2RAS transforms - (transforms, lta_convert, [('out', 'in_lta')]), ]) return workflow @@ -485,17 +491,19 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') - lta_convert = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_fsl'], - name='lta_convert') + fsl_to_lta = pe.MapNode(fs.utils.LTAConvert(out_lta=True), iterfield=['in_fsl'], + name='fsl_to_lta') workflow.connect([ - # Compare transforms - (flt_bbr, compare_transforms, [('out_matrix_file', 'test_mat')]), - (flt_bbr_init, compare_transforms, [('out_matrix_file', 'fallback_mat')]), - (compare_transforms, outputnode, [('out', 'fallback')]), - # Select output transform (flt_bbr, transforms, [('out_matrix_file', 'in1')]), (flt_bbr_init, transforms, [('out_matrix_file', 'in2')]), + # Convert FSL transforms to LTA (RAS2RAS) transforms and compare + (inputnode, fsl_to_lta, [('in_file', 'source_file'), + ('t1_brain', 'target_file')]), + (transforms, fsl_to_lta, [('out', 'in_fsl')]), + (fsl_to_lta, compare_transforms, [('out_lta', 'lta_list')]), + (compare_transforms, outputnode, [('out', 'fallback')]), + # Select output transform (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, invt_bbr, [('out', 'in_file')]), @@ -505,10 +513,6 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')]), - # Convert FSL transforms to LTA transforms - (transforms, lta_convert, [('out', 'in_fsl')]), - (inputnode, lta_convert, [('in_file', 'source_file'), - ('t1_brain', 'target_file')]), ]) return workflow From bd84448382e6c64923180e71ab119bb958ff3a2b Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 11 Oct 2017 12:41:53 -0400 Subject: [PATCH 16/17] RF: Use load_transform --- fmriprep/workflows/util.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 9ce2b5352..957bd6d19 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -27,19 +27,11 @@ def compare_xforms(lta_list, norm_threshold=15): - import numpy as np + from fmriprep.interfaces.surf import load_transform from niworkflows.nipype.algorithms.rapidart import _calc_norm_affine - def read_lta(fname): - with open(fname, 'rb') as fobj: - for line in fobj: - if line.strip() == b'1 4 4': - break - lines = fobj.readlines()[:4] - return np.genfromtxt(lines) - - bbr_affine = read_lta(lta_list[0]) - fallback_affine = read_lta(lta_list[1]) + bbr_affine = load_transform(lta_list[0]) + fallback_affine = load_transform(lta_list[1]) norm, _ = _calc_norm_affine([fallback_affine, bbr_affine], use_differences=True) From d7f9b29aa0db5a863ba646f0a2f38590ffae399b Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 11 Oct 2017 12:59:42 -0400 Subject: [PATCH 17/17] DOC: Update BBR docs --- fmriprep/workflows/util.py | 50 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/fmriprep/workflows/util.py b/fmriprep/workflows/util.py index 957bd6d19..9f00c7877 100644 --- a/fmriprep/workflows/util.py +++ b/fmriprep/workflows/util.py @@ -27,6 +27,37 @@ def compare_xforms(lta_list, norm_threshold=15): + """ + Computes a normalized displacement between two affine transforms as the + maximum overall displacement of the midpoints of the faces of a cube, when + each transform is applied to the cube. + This combines displacement resulting from scaling, translation and rotation. + + Although the norm is in mm, in a scaling context, it is not necessarily + equivalent to that distance in translation. + + We choose a default threshold of 15mm as a rough heuristic. + Normalized displacement above 20mm showed clear signs of distortion, while + "good" BBR refinements were frequently below 10mm displaced from the rigid + transform. + The 10-20mm range was more ambiguous, and 15mm chosen as a compromise. + This is open to revisiting in either direction. + + See discussion in + `GitHub issue #681`_ `_ + and the `underlying implementation + `_. + + Parameters + ---------- + + lta_list : list or tuple of str + the two given affines in LTA format + norm_threshold : float (default: 15) + the upper bound limit to the normalized displacement caused by the + second transform relative to the first + + """ from fmriprep.interfaces.surf import load_transform from niworkflows.nipype.algorithms.rapidart import _calc_norm_affine @@ -194,6 +225,16 @@ def init_bbreg_wf(use_bbr, bold2t1w_dof, omp_nthreads, name='bbreg_wf'): It is a counterpart to :py:func:`~fmriprep.workflows.util.init_fsl_bbr_wf`, which performs the same task using FSL's FLIRT with a BBR cost function. + The ``use_bbr`` option permits a high degree of control over registration. + If ``False``, standard, affine coregistration will be performed using + FreeSurfer's ``mri_coreg`` tool. + If ``True``, ``bbregister`` will be seeded with the initial transform found + by ``mri_coreg`` (equivalent to running ``bbregister --init-coreg``). + If ``None``, after ``bbregister`` is run, the resulting affine transform + will be compared to the initial transform found by ``mri_coreg``. + Excessive deviation will result in rejecting the BBR refinement and + accepting the original, affine registration. + .. workflow :: :graph2use: orig :simple_form: yes @@ -356,6 +397,15 @@ def init_fsl_bbr_wf(use_bbr, bold2t1w_dof, name='fsl_bbr_wf'): It is a counterpart to :py:func:`~fmriprep.workflows.util.init_bbreg_wf`, which performs the same task using FreeSurfer's ``bbregister``. + The ``use_bbr`` option permits a high degree of control over registration. + If ``False``, standard, rigid coregistration will be performed by FLIRT. + If ``True``, FLIRT-BBR will be seeded with the initial transform found by + the rigid coregistration. + If ``None``, after FLIRT-BBR is run, the resulting affine transform + will be compared to the initial transform found by FLIRT. + Excessive deviation will result in rejecting the BBR refinement and + accepting the original, affine registration. + .. workflow :: :graph2use: orig :simple_form: yes