ENH: Add number of volumes per *b*-value to QC report

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>Shell distribution: {shell_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,29 @@ def _generate_segment(self):
         )
 
 
+class DiffusionSummaryInputSpec(BaseInterfaceInputSpec):
+    pe_direction = traits.Enum(
+        None, "i", "i-", "j", "j-", "k", "k-", desc="Phase-encoding direction detected"
+    )
+    shell_dist = traits.Dict(mandatory=True, desc="Shell distribution")
+
+
+class DiffusionSummary(SummaryInterface):
+    input_spec = DiffusionSummaryInputSpec
+
+    def _generate_segment(self):
+        pe_direction = self.inputs.pe_direction
+        shell_dist = self.inputs.shell_dist
+        shell_dist_text = ", ".join(
+            f"{shell_dist[key]} directions at b={key} s/mm<sup>2</sup>"
+            for key in shell_dist
+        )
+
+        return DIFFUSION_TEMPLATE.format(
+            pe_direction=pe_direction, shell_dist=shell_dist_text
+        )
+
+
 class AboutSummaryInputSpec(BaseInterfaceInputSpec):
     version = Str(desc="dMRIPrep version")
     command = Str(desc="dMRIPrep command")

dmriprep/interfaces/vectors.py

@@ -20,6 +20,7 @@ class _CheckGradientTableInputSpec(BaseInterfaceInputSpec):
     in_rasb = File(exists=True, xor=["in_bval", "in_bvec"])
     b0_threshold = traits.Float(B0_THRESHOLD, usedefault=True)
     bvec_norm_epsilon = traits.Float(BVEC_NORM_EPSILON, usedefault=True)
+    b_mag = traits.Either(None, traits.Int, usedefault=True)
     b_scale = traits.Bool(True, usedefault=True)
 
 
@@ -29,6 +30,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
+    shell_dist = traits.Dict
     b0_ixs = traits.List(traits.Int)
     b0_mask = traits.List(traits.Bool)
 
@@ -48,6 +51,10 @@ class CheckGradientTable(SimpleInterface):
     (0.0, 0.0, 0.0)
     >>> check.outputs.full_sphere
     True
+    >>> check.outputs.num_shells
+    3
+    >>> check.outputs.shell_dist
+    {0.0: 12, 1200.0: 32, 2500.0: 61}
 
     >>> check = CheckGradientTable(
     ...     dwi_file=str(data_dir / 'dwi.nii.gz'),
@@ -57,6 +64,10 @@ class CheckGradientTable(SimpleInterface):
     (0.0, 0.0, 0.0)
     >>> check.outputs.full_sphere
     True
+    >>> check.outputs.num_shells
+    3
+    >>> check.outputs.shell_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)
@@ -75,14 +86,16 @@ def _run_interface(self, runtime):
             bvecs=_undefined(self.inputs, "in_bvec"),
             bvals=_undefined(self.inputs, "in_bval"),
             rasb_file=rasb_file,
+            b_mag=self.inputs.b_mag,
             b_scale=self.inputs.b_scale,
             bvec_norm_epsilon=self.inputs.bvec_norm_epsilon,
             b0_threshold=self.inputs.b0_threshold,
         )
         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["shell_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
@@ -16,6 +17,7 @@ class DiffusionGradientTable:
     __slots__ = [
         "_affine",
         "_b0_thres",
+        "_b_mag",
         "_b_scale",
         "_bvals",
         "_bvec_norm_epsilon",
@@ -29,6 +31,7 @@ class DiffusionGradientTable:
     def __init__(
         self,
         b0_threshold=B0_THRESHOLD,
+        b_mag=None,
         b_scale=True,
         bvals=None,
         bvec_norm_epsilon=BVEC_NORM_EPSILON,
@@ -45,12 +48,14 @@ def __init__(
         ----------
         b0_threshold : :obj:`float`
             The upper threshold to consider a low-b shell as :math:`b=0`.
+        b_mag : :obj:`int`
+            The order of magnitude to round the b-values.
         b_scale : :obj:`bool`
             Automatically scale the *b*-values with the norm of the corresponding
             *b*-vectors before the latter are normalized.
         bvals : str or os.pathlike or numpy.ndarray
             File path of the b-values.
-        b_vec_norm_epsilon : :obj:`float`
+        bvec_norm_epsilon : :obj:`float`
             The minimum difference in the norm of two *b*-vectors to consider them different.
         bvecs : str or os.pathlike or numpy.ndarray
             File path of the b-vectors.
@@ -92,6 +97,7 @@ def __init__(
         """
         self._affine = None
         self._b0_thres = b0_threshold
+        self._b_mag = b_mag
         self._b_scale = b_scale
         self._bvals = None
         self._bvec_norm_epsilon = bvec_norm_epsilon
@@ -174,6 +180,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."""
@@ -189,6 +200,7 @@ def normalize(self):
             self.bvals,
             b0_threshold=self._b0_thres,
             bvec_norm_epsilon=self._bvec_norm_epsilon,
+            b_mag=self._b_mag,
             b_scale=self._b_scale,
             raise_error=self._raise_inconsistent,
         )
@@ -284,14 +296,15 @@ def normalize_gradients(
     bvals,
     b0_threshold=B0_THRESHOLD,
     bvec_norm_epsilon=BVEC_NORM_EPSILON,
+    b_mag=None,
     b_scale=True,
     raise_error=False,
 ):
     """
     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
@@ -355,17 +368,25 @@ def normalize_gradients(
             raise ValueError(msg)
         config.loggers.cli.warning(msg)
 
-    # Rescale b-vals if requested
+    # Rescale bvals if requested
     if b_scale:
         bvals[~b0s] *= np.linalg.norm(bvecs[~b0s], axis=1) ** 2
 
     # Ensure b0s have (0, 0, 0) vectors
     bvecs[b0s, :3] = np.zeros(3)
 
     # Round bvals
-    bvals = round_bvals(bvals)
+    bvals = round_bvals(bvals, bmag=b_mag)
+
+    # Ensure rounding bvals doesn't change the number of b0s
+    rounded_b0s = bvals == 0
+    if not np.all(b0s == rounded_b0s):
+        msg = f"Inconsistent b0s before ({b0s.sum()}) and after rounding ({rounded_b0s.sum()})."
+        if raise_error:
+            raise ValueError(msg)
+        config.loggers.cli.warning(msg)
 
-    # Rescale b-vecs, skipping b0's, on the appropriate axis to unit-norm length.
+    # Rescale bvecs, skipping b0's, on the appropriate axis to unit-norm length.
     bvecs[~b0s] /= np.linalg.norm(bvecs[~b0s], axis=1)[..., np.newaxis]
     return bvecs, bvals.astype("uint16")
 

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, [("shell_dist", "shell_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