include number of distinct shells in reportlet

Author: josephmje

dmriprep/config/reports-spec.yml

@@ -72,6 +72,7 @@ sections:
 - name: Diffusion
   ordering: session,acquisition,run
   reportlets:
+  - bids: {datatype: figures, desc: summary, suffix: dwi}
   - bids: {datatype: figures, desc: validation, suffix: dwi}
   - bids: {datatype: figures, desc: brain, suffix: mask}
     caption: Average <em>b=0</em> that serves for reference in early preprocessing steps.

dmriprep/interfaces/reports.py

@@ -30,6 +30,16 @@
 \t</ul>
 """
 
+DIFFUSION_TEMPLATE = """\
+\t\t<details open>
+\t\t<summary>Summary</summary>
+\t\t<ul class="elem-desc">
+\t\t\t<li>Phase-encoding (PE) direction: {pe_direction}</li>
+\t\t\t<li>Distinct shells: {shells_dist}</li>
+\t\t</ul>
+\t\t</details>
+"""
+
 ABOUT_TEMPLATE = """\t<ul>
 \t\t<li>dMRIPrep version: {version}</li>
 \t\t<li>dMRIPrep command: <code>{command}</code></li>
@@ -119,6 +129,25 @@ def _generate_segment(self):
         )
 
 
+class DiffusionSummaryInputSpec(BaseInterfaceInputSpec):
+    pe_direction = traits.Enum(
+        None, "i", "i-", "j", "j-", "k", "k-", desc="Phase-encoding direction detected"
+    )
+    shells_dist = traits.Dict(mandatory=True, desc="Number of distinct shells")
+
+
+class DiffusionSummary(SummaryInterface):
+    input_spec = DiffusionSummaryInputSpec
+
+    def _generate_segment(self):
+        pe_direction = self.inputs.pe_direction
+        shells_dist = self.inputs.shells_dist
+
+        return DIFFUSION_TEMPLATE.format(
+            pe_direction=pe_direction, shells_dist=shells_dist
+        )
+
+
 class AboutSummaryInputSpec(BaseInterfaceInputSpec):
     version = Str(desc="dMRIPrep version")
     command = Str(desc="dMRIPrep command")

dmriprep/interfaces/vectors.py

@@ -29,6 +29,8 @@ class _CheckGradientTableOutputSpec(TraitedSpec):
     out_bvec = File(exists=True)
     full_sphere = traits.Bool()
     pole = traits.Tuple(traits.Float, traits.Float, traits.Float)
+    num_shells = traits.Int
+    shells_dist = traits.Dict
     b0_ixs = traits.List(traits.Int)
     b0_mask = traits.List(traits.Bool)
 
@@ -48,6 +50,10 @@ class CheckGradientTable(SimpleInterface):
     (0.0, 0.0, 0.0)
     >>> check.outputs.full_sphere
     True
+    >>> check.outputs.num_shells
+    3
+    >>> check.outputs.shells_dist
+    {0.0: 12, 1200.0: 32, 2500.0: 61}
 
     >>> check = CheckGradientTable(
     ...     dwi_file=str(data_dir / 'dwi.nii.gz'),
@@ -57,6 +63,10 @@ class CheckGradientTable(SimpleInterface):
     (0.0, 0.0, 0.0)
     >>> check.outputs.full_sphere
     True
+    >>> check.outputs.num_shells
+    3
+    >>> check.outputs.shells_dist
+    {0: 12, 1200: 32, 2500: 61}
     >>> newrasb = np.loadtxt(check.outputs.out_rasb, skiprows=1)
     >>> oldrasb = np.loadtxt(str(data_dir / 'dwi.tsv'), skiprows=1)
     >>> np.allclose(newrasb, oldrasb, rtol=1.e-3)
@@ -82,7 +92,8 @@ def _run_interface(self, runtime):
         pole = table.pole
         self._results["pole"] = tuple(pole)
         self._results["full_sphere"] = np.all(pole == 0.0)
-        self._results["b0_mask"] = table.b0mask.tolist()
+        self._results["num_shells"] = len(table.count_shells)
+        self._results["shells_dist"] = table.count_shells
         self._results["b0_ixs"] = np.where(table.b0mask)[0].tolist()
 
         cwd = Path(runtime.cwd).absolute()

dmriprep/utils/vectors.py

@@ -1,5 +1,6 @@
 """Utilities to operate on diffusion gradients."""
 from .. import config
+from collections import Counter
 from pathlib import Path
 from itertools import permutations
 import nibabel as nb
@@ -174,6 +175,11 @@ def bvals(self, value):
             raise ValueError("The number of b-vectors and b-values do not match")
         self._bvals = np.array(value)
 
+    @property
+    def count_shells(self):
+        """Count the number of volumes per b-value."""
+        return Counter(sorted(self._bvals))
+
     @property
     def b0mask(self):
         """Get a mask of low-b frames."""
@@ -291,7 +297,7 @@ def normalize_gradients(
     Normalize b-vectors and b-values.
 
     The resulting b-vectors will be of unit length for the non-zero b-values.
-    The resultinb b-values will be normalized by the square of the
+    The resulting b-values will be normalized by the square of the
     corresponding vector amplitude.
 
     Parameters

dmriprep/workflows/dwi/base.py

@@ -6,6 +6,7 @@
 
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 from ...interfaces import DerivativesDataSink
+from ...interfaces.reports import DiffusionSummary
 
 
 def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
@@ -41,6 +42,8 @@ def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
         File path of the b-vectors
     in_bval
         File path of the b-values
+    metadata
+        dwi metadata
     fmap
         File path of the fieldmap
     fmap_ref
@@ -82,6 +85,7 @@ def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
     config.loggers.workflow.debug(
         f"Creating DWI preprocessing workflow for <{dwi_file.name}>"
     )
+    metadata = layout.get_metadata(str(dwi_file))
 
     if has_fieldmap:
         import re
@@ -107,6 +111,7 @@ def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
                 "dwi_file",
                 "in_bvec",
                 "in_bval",
+                "metadata",
                 # From SDCFlows
                 "fmap",
                 "fmap_ref",
@@ -134,12 +139,19 @@ def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
     inputnode.inputs.dwi_file = str(dwi_file.absolute())
     inputnode.inputs.in_bvec = str(layout.get_bvec(dwi_file))
     inputnode.inputs.in_bval = str(layout.get_bval(dwi_file))
+    inputnode.metadata = metadata
 
     outputnode = pe.Node(
         niu.IdentityInterface(fields=["dwi_reference", "dwi_mask", "gradients_rasb"]),
         name="outputnode",
     )
 
+    summary = pe.Node(
+        DiffusionSummary(pe_direction=metadata.get("PhaseEncodingDirection")),
+        name="dwi_summary",
+        run_without_submitting=True,
+    )
+
     gradient_table = pe.Node(CheckGradientTable(), name="gradient_table")
 
     dwi_reference_wf = init_dwi_reference_wf(
@@ -153,6 +165,7 @@ def init_dwi_preproc_wf(dwi_file, has_fieldmap=False):
                                      ("in_bvec", "in_bvec"),
                                      ("in_bval", "in_bval")]),
         (inputnode, dwi_reference_wf, [("dwi_file", "inputnode.dwi_file")]),
+        (gradient_table, summary, [("shells_dist", "shells_dist")]),
         (gradient_table, dwi_reference_wf, [("b0_ixs", "inputnode.b0_ixs")]),
         (gradient_table, outputnode, [("out_rasb", "gradients_rasb")]),
     ])
@@ -254,6 +267,7 @@ def _bold_reg_suffix(fallback):
     # fmt: off
     workflow.connect([
         (inputnode, reportlets_wf, [("dwi_file", "inputnode.source_file")]),
+        (summary, reportlets_wf, [("out_report", "inputnode.summary_report")]),
         (dwi_reference_wf, reportlets_wf, [
             ("outputnode.validation_report", "inputnode.validation_report"),
         ]),
@@ -285,7 +299,6 @@ def _bold_reg_suffix(fallback):
     unwarp_wf = init_unwarp_wf(
         debug=config.execution.debug, omp_nthreads=config.nipype.omp_nthreads
     )
-    unwarp_wf.inputs.inputnode.metadata = layout.get_metadata(str(dwi_file))
 
     output_select = pe.Node(
         KeySelect(fields=["fmap", "fmap_ref", "fmap_coeff", "fmap_mask"]),
@@ -322,6 +335,7 @@ def _bold_reg_suffix(fallback):
         (dwi_reference_wf, coeff2epi_wf, [
             ("outputnode.ref_image", "inputnode.target_ref"),
             ("outputnode.dwi_mask", "inputnode.target_mask")]),
+        (inputnode, unwarp_wf, [("metadata", "inputnode.metadata")]),
         (dwi_reference_wf, unwarp_wf, [("outputnode.ref_image", "inputnode.distorted")]),
         (coeff2epi_wf, unwarp_wf, [
             ("outputnode.fmap_coeff", "inputnode.fmap_coeff")]),

dmriprep/workflows/dwi/outputs.py

@@ -15,6 +15,7 @@ def init_reportlets_wf(output_dir, sdc_report=False, name="reportlets_wf"):
         niu.IdentityInterface(
             fields=[
                 "source_file",
+                "summary_report",
                 "dwi_ref",
                 "dwi_mask",
                 "validation_report",
@@ -23,6 +24,15 @@ def init_reportlets_wf(output_dir, sdc_report=False, name="reportlets_wf"):
         ),
         name="inputnode",
     )
+
+    ds_report_summary = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir, desc="summary", datatype="figures"
+        ),
+        name="ds_report_summary",
+        run_without_submitting=True,
+    )
+
     mask_reportlet = pe.Node(SimpleShowMaskRPT(), name="mask_reportlet")
 
     ds_report_mask = pe.Node(
@@ -42,11 +52,13 @@ def init_reportlets_wf(output_dir, sdc_report=False, name="reportlets_wf"):
 
     # fmt:off
     workflow.connect([
+        (inputnode, ds_report_summary, [("source_file", "source_file"),
+                                        ("summary_report", "in_file")]),
         (inputnode, mask_reportlet, [("dwi_ref", "background_file"),
                                      ("dwi_mask", "mask_file")]),
-        (inputnode, ds_report_validation, [("source_file", "source_file")]),
+        (inputnode, ds_report_validation, [("source_file", "source_file"),
+                                           ("validation_report", "in_file")]),
         (inputnode, ds_report_mask, [("source_file", "source_file")]),
-        (inputnode, ds_report_validation, [("validation_report", "in_file")]),
         (mask_reportlet, ds_report_mask, [("out_report", "in_file")]),
     ])
     # fmt:on