b0 harmonization error if there is only one b0 image in a series

Summary

Incorrect b0 harmonization when there is only one b=0 image in a series. The reported mean of b0 is the mean of ALL volumes in a DWI series if there is only one b=0 in a series.

Additional details

QSIPrep version: 19.1 and current

What were you trying to do?

Run QSIPrep on a diffusion acquisition with two diffusion sequences. The first sequence contains one b=0 and multiple b=1000 volumes with phase encoding PA. The second sequence contains four b=0 volumes and multiple b-values up to 5000 with phase encoding AP. The sequences are matched on voxel size, TE, etc.

What did you expect to happen?

I expect the b0 harmonization to scale two series' intensities such that the b=0 intensities match. I expect this scaling to be very subtle since the raw b=0 volumes are inherently very similar in intensity.

What actually happened?

The harmonized dwi data contains b=0 volumes which differ greatly in intensity (a scaling factor of about 3). The result is incorrect FA maps, etc because the b=0 intensities (and DWI intensities) are not consistent across acquisition series.

Reproducing the bug

Please share any steps you performed that revealed the bug--> Please include any code snippets. Enclose them in triple back-ticks (```)

This error occurs if you have a diffusion series with one b=0 volume. The code below seems to have an error where if the number of b=0 volumes is 1, the b0_mean is the average of all dwi volumes (not just the b0 volumes). The error is most severe if you have mismatched diffusion series (like I do) where one series has one b=0 and the other has multiple b=0. With mismatched series, one series gets the mean of the b=0 volumes and the other series sets the mean of all the DWI volumes, producing a large scaling factor.

The error is suppressed by setting qsiprep to not perform b0 harmonization.

https://github.com/PennLINC/qsiprep/blob/36b93fe470272294b88069742f9fa1295bd69133/qsiprep/interfaces/dwi_merge.py#L690C1-L705C33

# Load the dwi data and get the mean values from the b=0 images
num_dwis = len(dwi_files)
dwi_niis = []
b0_means = []
for dwi_file, bval_file in zip(dwi_files, bvals):
    dwi_nii = load_img(dwi_file)
    _bvals = np.loadtxt(bval_file)
    b0_indices = np.flatnonzero(_bvals < b0_threshold)
    if b0_indices.size == 0:
        b0_mean = np.nan
    else:
        if len(b0_indices) > 1:
            b0_mean = index_img(dwi_nii, b0_indices).get_fdata().mean()
        else:
            b0_mean = dwi_nii.get_fdata().mean()
    b0_means.append(b0_mean)
    dwi_niis.append(dwi_nii)

PennLINC / qsiprep