diff --git a/dmriprep/interfaces/reports.py b/dmriprep/interfaces/reports.py index fecf8dd8..4842724e 100644 --- a/dmriprep/interfaces/reports.py +++ b/dmriprep/interfaces/reports.py @@ -35,7 +35,7 @@ \t\tSummary \t\t \t\t """ @@ -133,7 +133,7 @@ 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") + shell_dist = traits.Dict(mandatory=True, desc="Number of distinct shells") class DiffusionSummary(SummaryInterface): @@ -141,10 +141,14 @@ class DiffusionSummary(SummaryInterface): def _generate_segment(self): pe_direction = self.inputs.pe_direction - shells_dist = self.inputs.shells_dist + shell_dist = self.inputs.shell_dist + shell_dist_text = ", ".join( + f"{shell_dist[key]} directions at b={key} s/mm\u00B2" + for key in shell_dist + ) return DIFFUSION_TEMPLATE.format( - pe_direction=pe_direction, shells_dist=shells_dist + pe_direction=pe_direction, shell_dist=shell_dist_text ) diff --git a/dmriprep/interfaces/vectors.py b/dmriprep/interfaces/vectors.py index a3592087..b13d4bcb 100644 --- a/dmriprep/interfaces/vectors.py +++ b/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.Int(None, usedefault=True) b_scale = traits.Bool(True, usedefault=True) @@ -30,7 +31,7 @@ class _CheckGradientTableOutputSpec(TraitedSpec): full_sphere = traits.Bool() pole = traits.Tuple(traits.Float, traits.Float, traits.Float) num_shells = traits.Int - shells_dist = traits.Dict + shell_dist = traits.Dict b0_ixs = traits.List(traits.Int) @@ -51,7 +52,7 @@ class CheckGradientTable(SimpleInterface): True >>> check.outputs.num_shells 3 - >>> check.outputs.shells_dist + >>> check.outputs.shell_dist {0.0: 12, 1200.0: 32, 2500.0: 61} >>> check = CheckGradientTable( @@ -64,7 +65,7 @@ class CheckGradientTable(SimpleInterface): True >>> check.outputs.num_shells 3 - >>> check.outputs.shells_dist + >>> 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) @@ -84,6 +85,7 @@ 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, @@ -92,7 +94,7 @@ def _run_interface(self, runtime): self._results["pole"] = tuple(pole) self._results["full_sphere"] = np.all(pole == 0.0) self._results["num_shells"] = len(table.count_shells) - self._results["shells_dist"] = 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() diff --git a/dmriprep/utils/vectors.py b/dmriprep/utils/vectors.py index 6eeb206c..c86479b1 100644 --- a/dmriprep/utils/vectors.py +++ b/dmriprep/utils/vectors.py @@ -17,6 +17,7 @@ class DiffusionGradientTable: __slots__ = [ "_affine", "_b0_thres", + "_b_mag", "_b_scale", "_bvals", "_bvec_norm_epsilon", @@ -30,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, @@ -46,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. @@ -93,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 @@ -195,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, ) @@ -290,6 +296,7 @@ def normalize_gradients( bvals, b0_threshold=B0_THRESHOLD, bvec_norm_epsilon=BVEC_NORM_EPSILON, + b_mag=None, b_scale=True, raise_error=False, ): @@ -361,7 +368,7 @@ 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 @@ -369,9 +376,17 @@ def normalize_gradients( bvecs[b0s, :3] = np.zeros(3) # Round bvals - bvals = round_bvals(bvals) + bvals = round_bvals(bvals, bmag=b_mag) - # Rescale b-vecs, skipping b0's, on the appropriate axis to unit-norm length. + # 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 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") diff --git a/dmriprep/workflows/dwi/base.py b/dmriprep/workflows/dwi/base.py index 180c2c0f..c0f9f841 100644 --- a/dmriprep/workflows/dwi/base.py +++ b/dmriprep/workflows/dwi/base.py @@ -161,7 +161,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, summary, [("shell_dist", "shell_dist")]), (gradient_table, dwi_reference_wf, [("b0_ixs", "inputnode.b0_ixs")]), (gradient_table, outputnode, [("out_rasb", "gradients_rasb")]), ])