diff --git a/.docker/fsl-6.0/bin/topup b/.docker/fsl-6.0/bin/topup new file mode 100755 index 00000000..fed01a44 Binary files /dev/null and b/.docker/fsl-6.0/bin/topup differ diff --git a/.docker/fsl-6.0/lib/libgfortran.so.3 b/.docker/fsl-6.0/lib/libgfortran.so.3 new file mode 100755 index 00000000..ae186d1b Binary files /dev/null and b/.docker/fsl-6.0/lib/libgfortran.so.3 differ diff --git a/.docker/fsl-6.0/lib/libopenblas.so.0 b/.docker/fsl-6.0/lib/libopenblas.so.0 new file mode 100755 index 00000000..d898087b Binary files /dev/null and b/.docker/fsl-6.0/lib/libopenblas.so.0 differ diff --git a/Dockerfile b/Dockerfile index 98757650..a885d976 100644 --- a/Dockerfile +++ b/Dockerfile @@ -92,6 +92,9 @@ ENV FSLDIR="/usr/share/fsl/5.0" \ AFNI_PLUGINPATH="/usr/lib/afni/plugins" ENV PATH="/usr/lib/fsl/5.0:/usr/lib/afni/bin:$PATH" +COPY .docker/fsl-6.0/bin/topup /usr/share/fsl/5.0/bin/topup +COPY .docker/fsl-6.0/lib/* /usr/lib/fsl/5.0/ + # Installing ANTs 2.2.0 (NeuroDocker build) ENV ANTSPATH=/usr/lib/ants RUN mkdir -p $ANTSPATH && \ diff --git a/dmriprep/config/reports-spec.yml b/dmriprep/config/reports-spec.yml index cc0cc62c..1a01f242 100644 --- a/dmriprep/config/reports-spec.yml +++ b/dmriprep/config/reports-spec.yml @@ -26,6 +26,19 @@ sections: caption: Surfaces (white and pial) reconstructed with FreeSurfer (recon-all) overlaid on the participant's T1w template. subtitle: Surface reconstruction +- name: Fieldmaps + ordering: session,acquisition,run + reportlets: + - bids: {datatype: dwi, desc: pepolar, suffix: dwi} + caption: Inhomogeneities of the *B0* field introduce (oftentimes severe) spatial distortions + along the phase-encoding direction of the image. Utilizing two or more images with different + phase-encoding polarities (PEPolar) or directions, it is possible to estimate the inhomogeneity + of the field. The plot below shows a reference EPI (echo-planar imaging) volume generated + using two or more EPI images with varying phase-encoding blips. + description: Hover on the panels with the mouse pointer to also visualize the intensity of the + inhomogeneity of the field in Hertz. + static: false + subtitle: "Susceptibility-derived Distortion Correction (SDC): field inhomogeneity estimation" - name: Diffusion ordering: session,acquisition,run reportlets: diff --git a/dmriprep/data/flirtsch/b02b0.cnf b/dmriprep/data/flirtsch/b02b0.cnf new file mode 100644 index 00000000..77c4cd59 --- /dev/null +++ b/dmriprep/data/flirtsch/b02b0.cnf @@ -0,0 +1,26 @@ +# Resolution (knot-spacing) of warps in mm +--warpres=20,16,14,12,10,6,4,4,4 +# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel) +--subsamp=2,2,2,2,2,1,1,1,1 +# FWHM of gaussian smoothing +--fwhm=8,6,4,3,3,2,1,0,0 +# Maximum number of iterations +--miter=5,5,5,5,5,10,10,20,20 +# Relative weight of regularisation +--lambda=0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001 +# If set to 1 lambda is multiplied by the current average squared difference +--ssqlambda=1 +# Regularisation model +--regmod=bending_energy +# If set to 1 movements are estimated along with the field +--estmov=1,1,1,1,1,0,0,0,0 +# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient +--minmet=0,0,0,0,0,1,1,1,1 +# Quadratic or cubic splines +--splineorder=3 +# Precision for calculation and storage of Hessian +--numprec=double +# Linear or spline interpolation +--interp=spline +# If set to 1 the images are individually scaled to a common mean intensity +--scale=1 \ No newline at end of file diff --git a/dmriprep/data/flirtsch/b02b0_1.cnf b/dmriprep/data/flirtsch/b02b0_1.cnf new file mode 100644 index 00000000..07924b41 --- /dev/null +++ b/dmriprep/data/flirtsch/b02b0_1.cnf @@ -0,0 +1,26 @@ +# Resolution (knot-spacing) of warps in mm +--warpres=20,16,14,12,10,6,4,4,4 +# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel) +--subsamp=1,1,1,1,1,1,1,1,1 +# FWHM of gaussian smoothing +--fwhm=8,6,4,3,3,2,1,0,0 +# Maximum number of iterations +--miter=5,5,5,5,5,10,10,20,20 +# Relative weight of regularisation +--lambda=0.0005,0.0001,0.00001,0.0000015,0.0000005,0.0000005,0.00000005,0.0000000005,0.00000000001 +# If set to 1 lambda is multiplied by the current average squared difference +--ssqlambda=1 +# Regularisation model +--regmod=bending_energy +# If set to 1 movements are estimated along with the field +--estmov=1,1,1,1,1,0,0,0,0 +# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient +--minmet=0,0,0,0,0,1,1,1,1 +# Quadratic or cubic splines +--splineorder=3 +# Precision for calculation and storage of Hessian +--numprec=double +# Linear or spline interpolation +--interp=spline +# If set to 1 the images are individually scaled to a common mean intensity +--scale=1 \ No newline at end of file diff --git a/dmriprep/data/flirtsch/b02b0_2.cnf b/dmriprep/data/flirtsch/b02b0_2.cnf new file mode 100644 index 00000000..77c4cd59 --- /dev/null +++ b/dmriprep/data/flirtsch/b02b0_2.cnf @@ -0,0 +1,26 @@ +# Resolution (knot-spacing) of warps in mm +--warpres=20,16,14,12,10,6,4,4,4 +# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel) +--subsamp=2,2,2,2,2,1,1,1,1 +# FWHM of gaussian smoothing +--fwhm=8,6,4,3,3,2,1,0,0 +# Maximum number of iterations +--miter=5,5,5,5,5,10,10,20,20 +# Relative weight of regularisation +--lambda=0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001 +# If set to 1 lambda is multiplied by the current average squared difference +--ssqlambda=1 +# Regularisation model +--regmod=bending_energy +# If set to 1 movements are estimated along with the field +--estmov=1,1,1,1,1,0,0,0,0 +# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient +--minmet=0,0,0,0,0,1,1,1,1 +# Quadratic or cubic splines +--splineorder=3 +# Precision for calculation and storage of Hessian +--numprec=double +# Linear or spline interpolation +--interp=spline +# If set to 1 the images are individually scaled to a common mean intensity +--scale=1 \ No newline at end of file diff --git a/dmriprep/data/flirtsch/b02b0_4.cnf b/dmriprep/data/flirtsch/b02b0_4.cnf new file mode 100644 index 00000000..4e6f4266 --- /dev/null +++ b/dmriprep/data/flirtsch/b02b0_4.cnf @@ -0,0 +1,26 @@ +# Resolution (knot-spacing) of warps in mm +--warpres=20,16,14,12,10,6,4,4,4 +# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel) +--subsamp=4,4,2,2,2,1,1,1,1 +# FWHM of gaussian smoothing +--fwhm=8,6,4,3,3,2,1,0,0 +# Maximum number of iterations +--miter=5,5,5,5,5,10,10,20,20 +# Relative weight of regularisation +--lambda=0.035,0.006,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001 +# If set to 1 lambda is multiplied by the current average squared difference +--ssqlambda=1 +# Regularisation model +--regmod=bending_energy +# If set to 1 movements are estimated along with the field +--estmov=1,1,1,1,1,0,0,0,0 +# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient +--minmet=0,0,0,0,0,1,1,1,1 +# Quadratic or cubic splines +--splineorder=3 +# Precision for calculation and storage of Hessian +--numprec=double +# Linear or spline interpolation +--interp=spline +# If set to 1 the images are individually scaled to a common mean intensity +--scale=1 \ No newline at end of file diff --git a/dmriprep/workflows/base.py b/dmriprep/workflows/base.py index 934bbf51..cb48b851 100755 --- a/dmriprep/workflows/base.py +++ b/dmriprep/workflows/base.py @@ -12,13 +12,14 @@ BIDSInfo, BIDSFreeSurferDir ) from niworkflows.utils.misc import fix_multi_T1w_source_name +from niworkflows.utils.spaces import Reference from smriprep.workflows.anatomical import init_anat_preproc_wf -from niworkflows.utils.spaces import Reference from ..interfaces import DerivativesDataSink, BIDSDataGrabber from ..interfaces.reports import SubjectSummary, AboutSummary from ..utils.bids import collect_data from .dwi.base import init_early_b0ref_wf +from .fmap.base import init_fmap_estimation_wf def init_dmriprep_wf(): @@ -244,21 +245,27 @@ def init_single_subject_wf(subject_id): anat_preproc_wf.__postdesc__ = (anat_preproc_wf.__postdesc__ or "") + f""" Diffusion data preprocessing -: For each of the {len(subject_data["dwi"])} dwi scans found per subject - (across all sessions), the following preprocessing was performed.""" +: For each of the {len(subject_data["dwi"])} DWI scans found per subject + (across all sessions), the gradient table was vetted and converted into the *RASb* +format (i.e., given in RAS+ scanner coordinates, normalized b-vectors and scaled b-values), +and a *b=0* average for reference to the subsequent steps of preprocessing was calculated. +""" layout = config.execution.layout + dwi_data = tuple([ + (dwi, layout.get_metadata(dwi), layout.get_bvec(dwi), layout.get_bval(dwi)) + for dwi in subject_data["dwi"] + ]) inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_data"]), name="inputnode") - inputnode.iterables = [( - "dwi_data", tuple([ - (dwi, layout.get_bvec(dwi), layout.get_bval(dwi), - layout.get_metadata(dwi)["PhaseEncodingDirection"]) - for dwi in subject_data["dwi"] - ]) - )] + inputnode.iterables = [("dwi_data", dwi_data)] + + referencenode = pe.JoinNode(niu.IdentityInterface( + fields=["dwi_file", "metadata", "dwi_reference", "dwi_mask", "gradients_rasb"]), + name="referencenode", joinsource="inputnode", run_without_submitting=True) + split_info = pe.Node(niu.Function( - function=_unpack, output_names=["dwi_file", "bvec", "bval", "pedir"]), + function=_unpack, output_names=["dwi_file", "metadata", "bvec", "bval"]), name="split_info", run_without_submitting=True) early_b0ref_wf = init_early_b0ref_wf() @@ -267,8 +274,21 @@ def init_single_subject_wf(subject_id): (split_info, early_b0ref_wf, [("dwi_file", "inputnode.dwi_file"), ("bvec", "inputnode.in_bvec"), ("bval", "inputnode.in_bval")]), + (split_info, referencenode, [("dwi_file", "dwi_file"), + ("metadata", "metadata")]), + (early_b0ref_wf, referencenode, [ + ("outputnode.dwi_reference", "dwi_reference"), + ("outputnode.dwi_mask", "dwi_mask"), + ("outputnode.gradients_rasb", "gradients_rasb"), + ]) ]) + fmap_estimation_wf = init_fmap_estimation_wf(subject_data["dwi"]) + workflow.connect([ + (referencenode, fmap_estimation_wf, [ + ("dwi_reference", "inputnode.dwi_reference"), + ("dwi_mask", "inputnode.dwi_mask")]), + ]) return workflow diff --git a/dmriprep/workflows/dwi/base.py b/dmriprep/workflows/dwi/base.py index df770097..f24908d3 100644 --- a/dmriprep/workflows/dwi/base.py +++ b/dmriprep/workflows/dwi/base.py @@ -53,11 +53,6 @@ def init_early_b0ref_wf( """ # Build workflow workflow = Workflow(name=name) - workflow.__postdesc__ = """ -For every run and subject, the gradient table was vetted and converted into the *RASb* -format (i.e., given in RAS+ scanner coordinates, normalized b-vectors and scaled b-values), -and a *b=0* average for reference to the subsequent steps of preprocessing was calculated. -""" inputnode = pe.Node(niu.IdentityInterface( fields=['dwi_file', 'in_bvec', 'in_bval']), @@ -84,7 +79,7 @@ def init_early_b0ref_wf( (dwi_reference_wf, outputnode, [ ('outputnode.ref_image', 'dwi_reference'), ('outputnode.dwi_mask', 'dwi_mask')]), - (gradient_table, outputnode, [('out_rasb', 'out_rasb')]) + (gradient_table, outputnode, [('out_rasb', 'gradients_rasb')]) ]) # REPORTING ############################################################ diff --git a/dmriprep/workflows/dwi/util.py b/dmriprep/workflows/dwi/util.py index 1a868dca..30a106b3 100644 --- a/dmriprep/workflows/dwi/util.py +++ b/dmriprep/workflows/dwi/util.py @@ -6,6 +6,7 @@ from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.images import ValidateImage from niworkflows.interfaces.fixes import FixN4BiasFieldCorrection as N4BiasFieldCorrection +from niworkflows.interfaces.nibabel import ApplyMask from niworkflows.interfaces.utils import CopyXForm from ...interfaces.images import ExtractB0, RescaleB0 @@ -212,8 +213,10 @@ def init_enhance_and_skullstrip_dwi_wf( combine_masks = pe.Node(fsl.BinaryMaths(operation='mul'), name='combine_masks') + normalize = pe.Node(niu.Function(function=_normalize), name="normalize") + # Compute masked brain - apply_mask = pe.Node(fsl.ApplyMask(), name='apply_mask') + apply_mask = pe.Node(ApplyMask(), name='apply_mask') workflow.connect([ (inputnode, n4_correct, [('in_file', 'input_image'), @@ -230,11 +233,29 @@ def init_enhance_and_skullstrip_dwi_wf( (skullstrip_first_pass, combine_masks, [('mask_file', 'in_file')]), (skullstrip_second_pass, fixhdr_skullstrip2, [('out_file', 'in_file')]), (fixhdr_skullstrip2, combine_masks, [('out_file', 'operand_file')]), - (fixhdr_unifize, apply_mask, [('out_file', 'in_file')]), - (combine_masks, apply_mask, [('out_file', 'mask_file')]), + (combine_masks, apply_mask, [('out_file', 'in_mask')]), (combine_masks, outputnode, [('out_file', 'mask_file')]), + (n4_correct, normalize, [('output_image', 'in_file')]), + (normalize, apply_mask, [('out', 'in_file')]), + (normalize, outputnode, [('out', 'bias_corrected_file')]), (apply_mask, outputnode, [('out_file', 'skull_stripped_file')]), - (n4_correct, outputnode, [('output_image', 'bias_corrected_file')]), ]) return workflow + + +def _normalize(in_file): + from pathlib import Path + import numpy as np + import nibabel as nb + + nii = nb.load(in_file) + data = nii.get_fdata() + data[data < 0] = 0 + if data.max() >= 2**15 - 1: + data *= 2000 / np.percentile(data.reshape(-1), 98.0) + + hdr = nii.header.copy() + hdr.set_data_dtype('int16') + nii.__class__(data.astype('int16'), nii.affine, hdr).to_filename("normalized.nii.gz") + return Path("normalized.nii.gz").absolute() diff --git a/dmriprep/workflows/fmap/__init__.py b/dmriprep/workflows/fmap/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/dmriprep/workflows/fmap/base.py b/dmriprep/workflows/fmap/base.py new file mode 100644 index 00000000..4de6d5da --- /dev/null +++ b/dmriprep/workflows/fmap/base.py @@ -0,0 +1,136 @@ +"""Deploy a susceptibility-distortion estimation strategy.""" +from ... import config +from pkg_resources import resource_filename as _pkg_fname +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu + +from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from ...interfaces import DerivativesDataSink + + +def init_fmap_estimation_wf( + epi_targets, + generate_report=True, + name="fmap_estimation_wf", +): + """ + Setup a fieldmap estimation strategy and write results to derivatives folder. + + Parameters + ---------- + participant_label : :obj:`str` + The particular subject for which the BIDS layout will be queried. + epi_targets : :obj:`list` of :obj:`os.pathlike` + A list of :abbr:`EPI (echo planar imaging)` scans that will be corrected for + susceptibility distortions with the estimated fieldmaps. + omp_nthreads : :obj:`int` + Number of CPUs available to individual processes for multithreaded execution. + debug : :obj:`bool` + Whether fast (and less accurate) execution parameters should be used whenever available. + name : :obj:`str` + A unique workflow name to build Nipype's workflow hierarchy. + + """ + layout = config.execution.layout + wf = Workflow(name=name) + + inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_reference", "dwi_mask"]), + name="inputnode") + # Create one outputnode with a port for each potential EPI target + outputnode = pe.Node(niu.IdentityInterface(fields=[_fname2outname(p) for p in epi_targets]), + name="outputnode") + + # Return identity transforms for all if fieldmaps are ignored + if "fieldmaps" in config.workflow.ignore: + raise NotImplementedError + + # Set-up PEPOLAR estimators only with EPIs under fmap/ + # fmap_epi = {f: layout.get_metadata(f) + # for f in layout.get( + # subject=participant_label, datatype="fmap", + # suffix="epi", extension=("nii", "nii.gz"))} + + metadata = [layout.get_metadata(p) for p in epi_targets] + pedirs = [m.get("PhaseEncodingDirection", "unknown") for m in metadata] + if len(set(pedirs) - set(("unknown",))) > 1: + if "unknown" in pedirs or len(set(pe[0] for pe in set(pedirs))) > 1: + raise NotImplementedError + + # Get EPI polarities and their metadata + sdc_estimate_wf = init_pepolar_estimate_wf() + sdc_estimate_wf.inputs.inputnode.metadata = metadata + + wf.connect([ + (inputnode, sdc_estimate_wf, [("dwi_reference", "inputnode.in_data")]), + ]) + if generate_report: + from sdcflows.interfaces.reportlets import FieldmapReportlet + pepolar_report = pe.Node(FieldmapReportlet(reference_label="SDC'd B0"), + name="pepolar_report") + ds_report_pepolar = pe.Node(DerivativesDataSink( + base_directory=str(config.execution.work_dir / "reportlets"), keep_dtype=True, + desc="pepolar"), name="ds_report_pepolar") + ds_report_pepolar.inputs.source_file = epi_targets[0] + wf.connect([ + (sdc_estimate_wf, pepolar_report, [ + ("outputnode.fieldmap", "fieldmap"), + ("outputnode.corrected", "reference"), + ("outputnode.corrected_mask", "mask")]), + (pepolar_report, ds_report_pepolar, [("out_report", "in_file")]), + ]) + + return wf + + +def init_pepolar_estimate_wf(generate_report=True, name="pepolar_estimate_wf"): + """Initialize a barebones TOPUP implementation.""" + from nipype.interfaces.afni import Automask + from nipype.interfaces.fsl.epi import TOPUP + from niworkflows.interfaces.nibabel import MergeSeries + from ...interfaces.images import RescaleB0 + wf = Workflow(name=name) + + inputnode = pe.Node(niu.IdentityInterface(fields=["metadata", "in_data"]), + name="inputnode") + outputnode = pe.Node(niu.IdentityInterface(fields=["fieldmap", "corrected", "corrected_mask"]), + name="outputnode") + + concat_blips = pe.Node(MergeSeries(), name="concat_blips") + + topup = pe.Node(TOPUP(config=_pkg_fname("dmriprep", "data/flirtsch/b02b0.cnf")), + name="topup") + + pre_mask = pe.Node(Automask(dilate=1, outputtype="NIFTI_GZ"), + name="pre_mask") + rescale_corrected = pe.Node(RescaleB0(), name="rescale_corrected") + post_mask = pe.Node(Automask(outputtype="NIFTI_GZ"), + name="post_mask") + wf.connect([ + (inputnode, concat_blips, [("in_data", "in_files")]), + (inputnode, topup, [(("metadata", _get_ro), "readout_times"), + (("metadata", _get_pedir), "encoding_direction")]), + (concat_blips, topup, [("out_file", "in_file")]), + (topup, pre_mask, [("out_corrected", "in_file")]), + (pre_mask, rescale_corrected, [("out_file", "mask_file")]), + (topup, rescale_corrected, [("out_corrected", "in_file")]), + (topup, outputnode, [("out_field", "fieldmap")]), + (rescale_corrected, post_mask, [("out_ref", "in_file")]), + (rescale_corrected, outputnode, [("out_ref", "corrected")]), + (post_mask, outputnode, [("out_file", "corrected_mask")]), + ]) + + return wf + + +def _get_ro(metadata): + return [m["TotalReadoutTime"] for m in metadata] + + +def _get_pedir(metadata): + return [m["PhaseEncodingDirection"].replace("j", "y").replace("i", "x").replace("k", "z") + for m in metadata] + + +def _fname2outname(in_file): + from pathlib import Path + return Path(in_file).name.rstrip(".gz").rstrip(".nii").replace("-", "_") diff --git a/setup.cfg b/setup.cfg index 429b151b..cdcae2cf 100644 --- a/setup.cfg +++ b/setup.cfg @@ -24,7 +24,7 @@ install_requires = indexed_gzip >=0.8.8 nibabel ~= 3.0 nipype ~= 1.4 - niworkflows >=1.2.0rc4,<1.3 + niworkflows @ git+https://github.com/nipreps/niworkflows@master numpy pybids >=0.9.4 pyyaml @@ -82,6 +82,7 @@ dmriprep = VERSION config/reports-spec.yml data/boilerplate.bib + data/flirtsch/*.cnf data/tests/config.toml data/tests/THP/* data/tests/THP/sub-THP0005/anat/*