`data.dmri` has mismatches between reading/writting/representation

[jhlegarreta]

What happened?

The DWI reading/saving/representation logic has some bugs.

• The test_data_mri::test_load is reading the bvals/bvecs files
  

  nifreeze/test/test_data_dmri.py

  Lines 126 to 129 in 7322e93

  	dwi_from_nifti2 = load(
  	dwi_nifti_path,
  	bvec_file=bvecs_path,
  	bval_file=bvals_path,

  
  that are being written by the test function
  

  nifreeze/test/test_data_dmri.py

  Lines 112 to 113 in 7322e93

  	np.savetxt(str(bvecs_path), grad_table[:3])
  	np.savetxt(str(bvals_path), grad_table[-1])

  
  based on the grad_table variable, which is effectively inserting a b0 value to the bvals/bvecs file (having shapes (1, 103), (3, 103)), instead of writing the bval/bvec file pair that is written by the DWI.to_nifti method, which is what we should be testing.
• When insert_b0 is True, as it is the case in the test, the DWI.to_nifti method:
  

  nifreeze/src/nifreeze/data/dmri.py

  Lines 210 to 237 in 7322e93

  	if not insert_b0:
  	# Parent's to_nifti to handle the primary NIfTI export.
  	super().to_nifti(filename)
  	else:
  	data = np.concatenate((self.bzero[..., np.newaxis], self.dataobj), axis=-1)
  	nii = nb.Nifti1Image(data, self.affine, self.datahdr)
  	if self.datahdr is None:
  	nii.header.set_xyzt_units("mm")
  	nii.to_filename(filename)
  	# Convert filename to a Path object.
  	out_root = Path(filename).absolute()
  	# Remove .gz if present, then remove .nii if present.
  	# This yields the base stem for writing .bvec / .bval.
  	if out_root.suffix == ".gz":
  	out_root = out_root.with_suffix("") # remove '.gz'
  	if out_root.suffix == ".nii":
  	out_root = out_root.with_suffix("") # remove '.nii'
  	# Construct sidecar file paths.
  	bvecs_file = out_root.with_suffix(".bvec")
  	bvals_file = out_root.with_suffix(".bval")
  	# Save bvecs and bvals to text files
  	# Each row of bvecs is one direction (3 rows, N columns).
  	np.savetxt(bvecs_file, self.gradients[:3, ...].T, fmt="%.6f")
  	np.savetxt(bvals_file, self.gradients[:3, ...], fmt="%.6f")

  
  does not insert any null gradient to the bval/bvecs files (having shapes (1, 102), (3, 102)). Thus, if we had read those files in the test, the test would be failing due to a shape mismatch error in

  nifreeze/src/nifreeze/data/dmri.py

  Line 336 in 7322e93

  dataobj=fulldata[..., gradmsk],

as the data contained by the NIfTI files will have 103 DWI volumes (102 corresponding to non-b0 data, and the one inserted by virtue of insert_b0 being True), and the gradients read will have 102 values. Note that the represented DWI data has non-null (not None), built-in bzero data as well (besides the b0 data that the writer is asked to write following the insert_b0 flag).

However, there is a deeper bug here:

• If we change the DWI.to_nifti method to add a null gradient to the bval/bvecs files:

       bvecs = self.bvecs
       bvals = self.bvals

        if not insert_b0:
            # Parent's to_nifti to handle the primary NIfTI export.
            super().to_nifti(filename)
        else:
            data = np.concatenate((self.bzero[..., np.newaxis], self.dataobj), axis=-1)
            bvecs = np.concatenate((np.zeros(3)[:, np.newaxis], bvecs), axis=-1)
            bvals = np.concatenate((np.zeros(1), bvals))
            nii = nb.Nifti1Image(data, self.affine, self.datahdr)
            if self.datahdr is None:
                nii.header.set_xyzt_units("mm")
            nii.to_filename(filename)

        # Convert filename to a Path object.
        out_root = Path(filename).absolute()

        # Get the base stem for writing .bvec / .bval.
        out_root = out_root.parent / out_root.name.replace("".join(out_root.suffixes), "")

        # Construct sidecar file paths.
        bvecs_file = out_root.with_suffix(".bvec")
        bvals_file = out_root.with_suffix(".bval")

        # Save bvecs and bvals to text files
        # Each row of bvecs is one direction (3 rows, N columns).
        np.savetxt(bvecs_file, bvecs, fmt="%.6f")
        np.savetxt(bvals_file, bvals[np.newaxis, :], fmt="%.6f")

and if insert_b0 is True, the test will pass. However, it will not pass if insert_b0 is False. None of the cases on lines

nifreeze/src/nifreeze/data/dmri.py

Lines 345 to 356 in 7322e93

	if b0_file:
	b0img = nb.load(str(b0_file))
	b0vol = np.asanyarray(b0img.dataobj)
	# We'll assume your DWI class has a bzero: np.ndarray | None attribute
	dwi_obj.bzero = b0vol
	# Otherwise, if any volumes remain outside gradmsk, compute a median B0:
	elif np.any(~gradmsk):
	# The b=0 volumes are those that did NOT pass b0_thres
	b0_volumes = fulldata[..., ~gradmsk]
	# A simple approach is to take the median across that last dimension
	# Note that axis=3 is valid only if your data is 4D (x, y, z, volumes).
	dwi_obj.bzero = np.median(b0_volumes, axis=3)

are True and thus the read dwi_obj will have a None value in its bzero property, whereas the test object dwi_h5 does have a bzero volume, thus assert np.allclose(dwi_h5.bzero, dwi_from_nifti.bzero) will fail.

So the logic behind this needs to be reworked. Some thoughts:

• The above shows a clear mismatch of what the object holds, and what is told to write: if we tell to insert a bzero, the written file becomes out of sync with its original representation.
• It is confusing to have a bzero value (as the original instance read form the h5 file has), but not to have a 0-valued bval/bvec in the gradients.
• Commonly, we expect the non-weighted values to be present in the bval/bvec data.
• It is equally confusing the fact that the load method accepts an optional b0_file, but no b0_file is being written. Not writing an additional b0 file seems natural, as one expects the b0 data to be saved together with the non-b0 data, but then allowing to read one does not seem natural.
• It is confusing to ask to insert a b0 when writing to some data that already has b0 data: how should this be interpreted and written?
• What happens when a DWI volume has multiple b0 volumes?
• test_load should be parameterized with insert_b0 to [False, True] to test all of the above properly.

What command did you use?

pytest test_data_mri.py

What version of the software are you running?

HEAD (commit 7322e93) with the necessary fixes in PR #72

How are you running this software?

Local installation ("bare-metal")

Is your data BIDS valid?

Yes

Are you reusing any previously computed results?

No

Please copy and paste any relevant log output.

For the issue (i.e. not writing the bval/bvec files that should be read, and not writing the null gradient to the bval/bvec files):

>       dwi_from_nifti1 = load(
            dwi_nifti_path,
            bvec_file=bvecs_path,
            bval_file=bvals_path,
        )
(...)

        # The shape checking is somewhat flexible: (4, N) or (N, 4)
        dwi_obj = DWI(
>           dataobj=fulldata[..., gradmsk],
            affine=affine,
            # We'll assign the filtered gradients below.
        )
E       IndexError: boolean index did not match indexed array along dimension 3; dimension is 103 but corresponding boolean dimension is 102

For the second issue (having made such that the above passes add a null gradient to the written bval and bvec files but using False for insert_b0):

>       assert np.allclose(dwi_h5.bzero, dwi_from_nifti1.bzero)
(...)
E       TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Additional information / screenshots

X-ref #72 (comment).

[jhlegarreta]

After discussing with Oscar:

Commonly, we expect the non-weighted values to be present in the bval/bvec data.

The rationale behind excluding those is that we want NiFreeze to have the possibility to load "good" b0 data computed using some clever strategy and to free NiFreeze from the burden of dealing with b0s, as the motion/ECD correction needs to be done on non-null diffusion-weighted data.

It is equally confusing the fact that the load method accepts an optional b0_file, but no b0_file is being written. Not writing an additional b0 file seems natural, as one expects the b0 data to be saved together with the non-b0 data, but then allowing to read one does not seem natural.

If insert_b0 is True, the b0 data is written both to the NIfTI (together with the rest of the DWI data) and bval/bvec files. Re the first point about the rationale behind offering to load b0 data from an external file.

It is confusing to ask to insert a b0 when writing to some data that already has b0 data: how should this be interpreted and written?

When the DWI instance is build, a non-null gradient mask is created, and the bzero attribute is assumed to have some data: no new data is inserted.

What happens when a DWI volume has multiple b0 volumes?

It is dealt with here:

nifreeze/src/nifreeze/data/dmri.py

Lines 373 to 377 in e964b96

	dwi_obj = DWI(
	dataobj=fulldata[..., gradmsk],
	affine=affine,
	# We'll assign the filtered gradients below.
	)