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).