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 a3bad5c8..8beea1db 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,run + reportlets: + - bids: {datatype: figures, desc: pepolar, suffix: fieldmap} + 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 b2892478..22fff3c0 100755 --- a/dmriprep/workflows/base.py +++ b/dmriprep/workflows/base.py @@ -19,6 +19,7 @@ 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(): @@ -292,6 +293,12 @@ def init_single_subject_wf(subject_id): ]), ]) + 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/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..4e97f147 --- /dev/null +++ b/dmriprep/workflows/fmap/base.py @@ -0,0 +1,137 @@ +"""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.output_dir), datatype="figures", + suffix="fieldmap", desc="pepolar", dismiss_entities=("acquisition", "dir")), + 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 8f58ddc0..61b04423 100644 --- a/setup.cfg +++ b/setup.cfg @@ -50,6 +50,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/*