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/data/flirtsch/b02b0_quick.cnf b/dmriprep/data/flirtsch/b02b0_quick.cnf new file mode 100644 index 00000000..5eef46f8 --- /dev/null +++ b/dmriprep/data/flirtsch/b02b0_quick.cnf @@ -0,0 +1,26 @@ +# Resolution (knot-spacing) of warps in mm +--warpres=20,14,10 +# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel) +--subsamp=2,2,2 +# FWHM of gaussian smoothing +--fwhm=8,2,1 +# Maximum number of iterations +--miter=5,5,5 +# Relative weight of regularisation +--lambda=0.005,0.000005,0.0000005 +# 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=0,0,0 +# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient +--minmet=0,0,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..2b36ae40 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,13 @@ def init_single_subject_wf(subject_id): ]), ]) + fmap_estimation_wf = init_fmap_estimation_wf( + subject_data["dwi"], debug=config.execution.debug) + 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..0c306ab2 --- /dev/null +++ b/dmriprep/workflows/fmap/base.py @@ -0,0 +1,145 @@ +"""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, + debug=False, + 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") + wf.add_nodes([inputnode]) # TODO: remove when fully implemented + # 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: + return wf + + # 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] + if any("TotalReadoutTime" not in m for m in metadata): + return wf + + 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(debug=debug) + 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(debug=False, 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 sdcflows.interfaces.fmap import get_trt + 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") + readout_time = pe.MapNode(niu.Function( + input_names=["in_meta", "in_file"], function=get_trt), name="readout_time", + iterfield=["in_meta", "in_file"], run_without_submitting=True + ) + + topup = pe.Node(TOPUP(config=_pkg_fname( + "dmriprep", f"data/flirtsch/b02b0{'_quick' * debug}.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, readout_time, [("in_data", "in_file"), + ("metadata", "in_meta")]), + (inputnode, topup, [(("metadata", _get_pedir), "encoding_direction")]), + (readout_time, topup, [("out", "readout_times")]), + (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_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/*