ENH: Add a shell data property to DWI data class

Add a shell data property to `DWI` data class that returns a list of
pairs consisting of the estimated b-value and the associated DWI data.

Author: jhlegarreta

src/nifreeze/data/dmri.py

@@ -241,6 +241,46 @@ def to_nifti(self, filename: Path | str, insert_b0: bool = False) -> None:
         np.savetxt(bvecs_file, self.gradients[:3, ...].T, fmt="%.6f")
         np.savetxt(bvals_file, self.gradients[:3, ...], fmt="%.6f")
 
+    @property
+    def shells(
+        self,
+        num_bins: int = DEFAULT_NUM_BINS,
+        multishell_nonempty_bin_count_thr: int = DEFAULT_MULTISHELL_BIN_COUNT_THR,
+        bval_cap: int = DEFAULT_HIGHB_THRESHOLD,
+    ):
+        """Get the shell data according to the b-value groups.
+
+        Bin the shell data according to the b-value groups found by `~find_shelling_scheme`.
+
+        Parameters
+        ----------
+        num_bins : :obj:`int`, optional
+            Number of bins.
+        multishell_nonempty_bin_count_thr : :obj:`int`, optional
+            Bin count to consider a multi-shell scheme.
+        bval_cap : :obj:`int`, optional
+            Maximum b-value to be considered in a multi-shell scheme.
+
+        Returns
+        -------
+        :obj:`list`
+            Tuples of binned b-values and corresponding shell data.
+        """
+
+        _, bval_groups, bval_estimated = find_shelling_scheme(
+            self.gradients[-1, ...],
+            num_bins=num_bins,
+            multishell_nonempty_bin_count_thr=multishell_nonempty_bin_count_thr,
+            bval_cap=bval_cap,
+        )
+        indices = [
+            np.hstack(np.where(np.isin(self.gradients[-1, ...], bvals))) for bvals in bval_groups
+        ]
+        return [
+            (bval_estimated[idx], self.dataobj[indices, ...])
+            for idx, indices in enumerate(indices)
+        ]
+
 
 def load(
     filename: Path | str,

test/test_data_dmri.py

@@ -182,6 +182,29 @@ def test_equality_operator(tmp_path):
     assert round_trip_dwi_obj == dwi_obj
 
 
+def test_shells(datadir):
+
+    dwi_h5 = load(datadir / "dwi.h5")
+    num_bins = 3
+
+    _, expected_bval_groups, expected_bval_est = find_shelling_scheme(
+        dwi_h5.gradients[-1, ...], num_bins=num_bins
+    )
+
+    indices = [
+        np.hstack(np.where(np.isin(dwi_h5.gradients[-1, ...], bvals)))
+        for bvals in expected_bval_groups
+    ]
+    expected_shell_data = [dwi_h5.dataobj[indices, ...] for indices in indices]
+
+    shells = dwi_h5.shells(num_bins=num_bins)
+    obtained_bval_est, obtaine_shell_data = zip(*shells, strict=True)
+
+    assert len(shells) == num_bins
+    assert obtained_bval_est == expected_bval_est
+    assert np.allclose(obtaine_shell_data, expected_shell_data)
+
+
 @pytest.mark.parametrize(
     ("bvals", "exp_scheme", "exp_bval_groups", "exp_bval_estimated"),
     [