[FIX] Fixes DWI preprocessing using T1 when images are 3D (#1169)

* work on supporting 3D DWI images

* update unit tests

* try filtering dwi images before running the pipeline

* fixes

* try more aggressive filtering

clinica/pipelines/dwi/preprocessing/t1/pipeline.py

@@ -86,6 +86,48 @@ def get_output_fields(self) -> List[str]:
         """
         return ["preproc_dwi", "preproc_bvec", "preproc_bval", "b0_mask"]
 
+    def filter_qc(self) -> tuple[list[str], list[str]]:
+        from clinica.pipelines.dwi.preprocessing.utils import check_dwi_volume
+        from clinica.pipelines.dwi.utils import DWIDataset
+        from clinica.utils.bids import BIDSFileName
+        from clinica.utils.input_files import (
+            DWI_BVAL,
+            DWI_BVEC,
+            DWI_NII,
+        )
+        from clinica.utils.inputs import clinica_list_of_files_reader
+        from clinica.utils.stream import cprint
+
+        subjects = []
+        sessions = []
+        list_bids_files = clinica_list_of_files_reader(
+            self.subjects,
+            self.sessions,
+            self.bids_directory,
+            [DWI_NII, DWI_BVEC, DWI_BVAL],
+            raise_exception=True,
+        )
+        for dwi_image_file, b_vectors_file, b_values_file in zip(*list_bids_files):
+            dwi_image_file_bids = BIDSFileName.from_name(dwi_image_file)
+            try:
+                check_dwi_volume(
+                    DWIDataset(
+                        dwi=dwi_image_file,
+                        b_values=b_values_file,
+                        b_vectors=b_vectors_file,
+                    )
+                )
+                subjects.append(f"sub-{dwi_image_file_bids.subject}")
+                sessions.append(f"ses-{dwi_image_file_bids.session}")
+            except (IOError, ValueError) as e:
+                cprint(
+                    f"Ignoring DWI scan for subject {dwi_image_file_bids.subject} "
+                    f"and session {dwi_image_file_bids.session} for the following reason: {e}",
+                    lvl="warning",
+                )
+
+        return subjects, sessions
+
     def _build_input_node(self):
         """Build and connect an input node to the pipeline."""
         import nipype.interfaces.utility as nutil
@@ -278,6 +320,7 @@ def _build_core_nodes(self):
             use_cuda=self.parameters["use_cuda"],
             initrand=self.parameters["initrand"],
             compute_mask=True,
+            image_id=True,
         )
         # Susceptibility distortion correction using T1w image
         sdc = epi_pipeline(
@@ -337,7 +380,11 @@ def _build_core_nodes(self):
                     eddy_fsl,
                     [
                         (x, f"inputnode.{x}")
-                        for x in ("total_readout_time", "phase_encoding_direction")
+                        for x in (
+                            "total_readout_time",
+                            "phase_encoding_direction",
+                            "image_id",
+                        )
                     ],
                 ),
                 (

clinica/pipelines/dwi/preprocessing/t1/utils.py

@@ -24,7 +24,6 @@
     "configure_working_directory",
     "register_b0",
     "extract_sub_ses_folder_name",
-    "split_dwi_dataset_with_b_values",
 ]
 
 
@@ -333,17 +332,17 @@ def prepare_reference_b0(
 
     dwi_dataset = check_dwi_dataset(dwi_dataset)
     working_directory = configure_working_directory(dwi_dataset.dwi, working_directory)
-    small_b_dataset, large_b_dataset = split_dwi_dataset_with_b_values(
+    small_b_dataset, large_b_dataset = _split_dwi_dataset_with_b_values(
         dwi_dataset,
         b_value_threshold=b_value_threshold,
         working_directory=working_directory,
     )
     reference_b0 = compute_reference_b0(
         small_b_dataset.dwi, dwi_dataset.b_values, b_value_threshold, working_directory
     )
-    reference_dataset = insert_b0_into_dwi(reference_b0, large_b_dataset)
-
-    return reference_b0, reference_dataset
+    if large_b_dataset:
+        return reference_b0, insert_b0_into_dwi(reference_b0, large_b_dataset)
+    return reference_b0, small_b_dataset
 
 
 def insert_b0_into_dwi(b0_filename: PathLike, dwi_dataset: DWIDataset) -> DWIDataset:
@@ -624,11 +623,11 @@ def extract_sub_ses_folder_name(file_path: str) -> str:
     return (Path(Path(file_path).parent).parent).name
 
 
-def split_dwi_dataset_with_b_values(
+def _split_dwi_dataset_with_b_values(
     dwi_dataset: DWIDataset,
     b_value_threshold: float = 5.0,
     working_directory: Optional[Path] = None,
-) -> Tuple[DWIDataset, DWIDataset]:
+) -> Tuple[DWIDataset, Optional[DWIDataset]]:
     """Splits the DWI dataset in two through the B-values.
 
     Splits the DWI volumes into two datasets :
@@ -676,21 +675,24 @@ def split_dwi_dataset_with_b_values(
     large_b_filter = np.array(
         [i for i in range(len(b_values)) if i not in small_b_filter]
     )
-
-    return (
-        _build_dwi_dataset_from_filter(
+    large = None
+    if len(small_b_filter) > 0:
+        small = _build_dwi_dataset_from_filter(
             dwi_dataset,
             "small_b",
             small_b_filter,
             working_directory=working_directory,
-        ),
-        _build_dwi_dataset_from_filter(
+        )
+    else:
+        raise ValueError("No small dataset found.")
+    if len(large_b_filter) > 0:
+        large = _build_dwi_dataset_from_filter(
             dwi_dataset,
             "large_b",
             large_b_filter,
             working_directory=working_directory,
-        ),
-    )
+        )
+    return small, large
 
 
 def _check_b_values_and_b_vectors(
@@ -699,8 +701,8 @@ def _check_b_values_and_b_vectors(
     """Opens the b-values and b-vectors files and transpose the b-vectors if needed."""
     import warnings
 
-    b_values = np.loadtxt(dwi_dataset.b_values)
-    b_vectors = np.loadtxt(dwi_dataset.b_vectors)
+    b_values = np.loadtxt(dwi_dataset.b_values, ndmin=1)
+    b_vectors = np.loadtxt(dwi_dataset.b_vectors, ndmin=2)
     if b_values.shape[0] == b_vectors.shape[0]:
         warnings.warn(
             "Warning: The b-vectors file should be column-wise. The b-vectors will be transposed",
@@ -774,16 +776,21 @@ def _filter_dwi(
     working_directory: Optional[Path] = None,
 ) -> Path:
     """Filters the dwi component of the provided DWI dataset."""
+    import nibabel as nib
+
     from clinica.utils.image import compute_aggregated_volume, get_new_image_like
 
     from ..utils import add_suffix_to_filename
 
     dwi_filename = add_suffix_to_filename(dwi_dataset.dwi, filter_name)
     if working_directory:
         dwi_filename = working_directory / dwi_filename.name
-    data = compute_aggregated_volume(
-        dwi_dataset.dwi, aggregator=None, volumes_to_keep=filter_array
-    )
+    try:
+        data = compute_aggregated_volume(
+            dwi_dataset.dwi, aggregator=None, volumes_to_keep=filter_array
+        )
+    except ValueError:
+        data = nib.load(dwi_dataset.dwi).get_fdata()
     img = get_new_image_like(dwi_dataset.dwi, data)
     img.to_filename(dwi_filename)
 

clinica/pipelines/dwi/preprocessing/utils.py

@@ -57,11 +57,11 @@ def generate_index_file(
     if not b_values_filename.is_file():
         raise FileNotFoundError(f"Unable to find b-values file: {b_values_filename}.")
 
-    b_values = np.loadtxt(b_values_filename)
+    b_values = np.loadtxt(b_values_filename, ndmin=1)
     index_filename = f"{image_id}_index.txt" if image_id else "index.txt"
     output_dir = output_dir or b_values_filename.parent
     index_filename = output_dir / index_filename
-    np.savetxt(index_filename, np.ones(len(b_values)).T)
+    np.savetxt(index_filename, np.ones(len(b_values), dtype=int).T)
 
     return index_filename
 
@@ -296,20 +296,46 @@ def check_dwi_volume(dwi_dataset: DWIDataset) -> None:
     Raises
     ------
     ValueError:
-        if the number of DWI volumes, the number of B-values,
+        If the DWI nifti image is not 4D.
+
+    IOError:
+        If the b-values and/or b-vectors are all zeros.
+
+        If the number of DWI volumes, the number of B-values,
         and the number of B-vectors are not equal.
     """
-    num_b_values = len(np.loadtxt(dwi_dataset.b_values))
-    num_b_vectors = np.loadtxt(dwi_dataset.b_vectors).shape[-1]
-    num_dwi = nib.load(dwi_dataset.dwi).shape[-1]
+    b_values = np.loadtxt(dwi_dataset.b_values, ndmin=1)
+    num_b_values = len(b_values)
+    b_vectors = np.loadtxt(dwi_dataset.b_vectors, ndmin=2)
+    num_b_vectors = b_vectors.shape[-1]
+    num_dwi = _load_nifti_at_least_4d(dwi_dataset.dwi, raise_if_3d=True).shape[-1]
 
+    if np.sum(b_values) == 0 or np.sum(b_vectors) == 0:
+        raise IOError(
+            "The b-values and/or b-vectors have problematic values (all zeros)."
+        )
     if not (num_b_values == num_b_vectors == num_dwi):
         raise IOError(
             f"Number of DWIs, b-vals and b-vecs mismatch "
             f"(# DWI = {num_dwi}, # B-vec = {num_b_vectors}, #B-val = {num_b_values}) "
         )
 
 
+def _load_nifti_at_least_4d(
+    image_path: Path, raise_if_3d: bool = False
+) -> nib.Nifti1Image:
+    image = nib.load(image_path)
+    if len(image.shape) == 3:
+        if not raise_if_3d:
+            data = image.get_fdata()
+            image = nib.Nifti1Image(
+                np.expand_dims(data, axis=-1), image.affine, image.header
+            )
+    if len(image.shape) != 4:
+        raise ValueError(f"DWI image {image_path} has not a valid number of dimension.")
+    return image
+
+
 def check_dwi_dataset(dwi_dataset: DWIDataset) -> DWIDataset:
     """Checks that provided input files of the DWI dataset exist.
 

clinica/pipelines/dwi/preprocessing/workflows.py

@@ -126,7 +126,7 @@ def eddy_fsl_pipeline(
 
     generate_index = pe.Node(
         niu.Function(
-            input_names=["b_values_filename", "output_dir"],
+            input_names=["b_values_filename", "image_id", "output_dir"],
             output_names=["out_file"],
             function=generate_index_file_task,
         ),
@@ -162,7 +162,6 @@ def eddy_fsl_pipeline(
             [("phase_encoding_direction", "fsl_phase_encoding_direction")],
         ),
         (inputnode, generate_index, [("b_values_filename", "b_values_filename")]),
-        (inputnode, generate_index, [("image_id", "image_id")]),
         (inputnode, eddy, [("b_vectors_filename", "in_bvec")]),
         (inputnode, eddy, [("b_values_filename", "in_bval")]),
         (inputnode, eddy, [("dwi_filename", "in_file")]),
@@ -176,6 +175,7 @@ def eddy_fsl_pipeline(
     if image_id:
         connections += [
             (inputnode, generate_acq, [("image_id", "image_id")]),
+            (inputnode, generate_index, [("image_id", "image_id")]),
             (inputnode, eddy, [("image_id", "out_base")]),
         ]
 

clinica/pipelines/engine.py

@@ -477,10 +477,14 @@ def __init__(
             use_session_tsv=False,
             tsv_dir=base_dir,
         )
+        self._subjects, self._sessions = self.filter_qc()
         self._input_node = None
         self._output_node = None
         self._init_nodes()
 
+    def filter_qc(self) -> tuple[list[str], list[str]]:
+        return self._subjects, self._sessions
+
     @property
     def base_dir_was_specified(self) -> bool:
         return self._base_dir_was_specified

test/unittests/pipelines/dwi/preprocessing/test_t1_utils.py

@@ -461,7 +461,7 @@ def dwi_dataset(tmp_path):
 
 def test_split_dwi_dataset_with_b_values_errors(tmp_path, dwi_dataset):
     from clinica.pipelines.dwi.preprocessing.t1.utils import (
-        split_dwi_dataset_with_b_values,
+        _split_dwi_dataset_with_b_values,  # noqa
     )
 
     for filename in ("foo.nii.gz", "foo.bval", "foo.bvec"):
@@ -471,13 +471,13 @@ def test_split_dwi_dataset_with_b_values_errors(tmp_path, dwi_dataset):
         ValueError,
         match="b_value_threshold should be >=0. You provided -1.",
     ):
-        split_dwi_dataset_with_b_values(dwi_dataset, b_value_threshold=-1)
+        _split_dwi_dataset_with_b_values(dwi_dataset, b_value_threshold=-1)
 
 
 @pytest.mark.parametrize("extension", ["nii", "nii.gz"])
 def test_split_dwi_dataset_with_b_values(tmp_path, extension):
     from clinica.pipelines.dwi.preprocessing.t1.utils import (
-        split_dwi_dataset_with_b_values,
+        _split_dwi_dataset_with_b_values,  # noqa
     )
     from clinica.pipelines.dwi.utils import DWIDataset
 
@@ -495,7 +495,7 @@ def test_split_dwi_dataset_with_b_values(tmp_path, extension):
         b_values=tmp_path / "foo.bval",
         b_vectors=tmp_path / "foo.bvec",
     )
-    small_b_dataset, large_b_dataset = split_dwi_dataset_with_b_values(dwi_dataset)
+    small_b_dataset, large_b_dataset = _split_dwi_dataset_with_b_values(dwi_dataset)
     assert small_b_dataset.dwi == tmp_path / f"foo_small_b.{extension}"
     assert small_b_dataset.b_values == tmp_path / "foo_small_b.bval"
     assert small_b_dataset.b_vectors == tmp_path / "foo_small_b.bvec"