Wrong registration between fMRI data and fieldmap

[1RuobingLiu]

What happened?

Hi every one, recently I am working on a fMRIPrep for

this example dataset for subeject 33 session Y01: original -containning two imaginary and two real fieldmaps apart from two magnetitude fieldmaps under fmap folder; and processed – I merged imaginary and real fieldmaps into two phase fieldmaps first using Python Nib.arctan function.
I also tried run on my local pc, same results.

The fmriprep-output shows that there is something wrong with registration.

There is something wrong with bold brain mask.
This is the bold mask if we apply fieldmap correction:

This is the bold mask if we do not apply fieldmap:

Can you help me on this, please? How can I fix the registration problem? I have been stuck on this for a long while. Thank you! Appreciate it if you could help!

What command did you use?

#hcc platform
apptainer exec docker://nipreps/fmriprep fmriprep $inPath $outPath participant --participant-label sub-NCANDA00033 --output-spaces MNI152NLin2009cAsym MNI152NLin6Asym:res-2 anat -w $cachePath --fs-license-file $licFile --skip_bids_validation --nthreads 5 --low-mem --cifti-output --ignore slicetiming
#local pc
python -m fmriprep /Users/ruoliu/Downloads/NCANDA-sub95/raw /Users/ruoliu/Downloads/NCANDA-sub95/derivatives/fmriprep participant \
  --participant-label sub-NCANDA00095 \
  --output-spaces MNI152NLin2009cAsym MNI152NLin6Asym:res-2 anat \
  -w /Users/ruoliu/Downloads/NCANDA-sub95/cache/fmriprep \
  --fs-license-file /Users/ruoliu/Downloads/NCANDA-sub95/raw/license.txt \
  --skip_bids_validation \
  --nthreads 5 \
  --low-mem \
  --cifti-output

What version of fMRIPrep are you running?

24.1.1

How are you running fMRIPrep?

Docker

Is your data BIDS valid?

Yes

Are you reusing any previously computed results?

Work directory

Please copy and paste any relevant log output.

No response

Additional information / screenshots

No response

[effigies]

Are the data public? I haven't seen something like this. I could try to reproduce it locally and investigate.

[1RuobingLiu]

Are the data public? I haven't seen something like this. I could try to reproduce it locally and investigate.

It is not public, but may I share one subject to your google drive ([email protected]) please? I have tried to debug it, it seems to be random walker from skimage.segmentation problem in bold_fit_wf (unwarp_wf (brainextraction_wf (masker (brainmask). Could you help me please? Thank you!

[effigies]

This username at gmail works.

[1RuobingLiu]

This username at gmail works.

Already sent. Thank you

[effigies]

Reproduced:

[1RuobingLiu]

Reproduced:

Yes, that’s what I got. It’s a bit weird right. Think it’s bold mask generation problem

[1RuobingLiu]

Reproduced:

I found the issue could be related to bold unwarp brainmasker https://github.com/nipreps/sdcflows/blob/master/sdcflows/utils/tools.py#L154, and I tried different ball and padding number, but the masks all seem not good. Could you help me please? Thank you.

Here is the code in the bold mask generation:

def brain_masker(in_file, out_file=None, padding=5):
    """Use grayscale morphological operations to obtain a quick mask of EPI data with visualization."""
    from pathlib import Path
    import re
    import nibabel as nb
    import numpy as np
    from scipy import ndimage
    from skimage.morphology import ball
    from skimage.filters import threshold_otsu
    from skimage.segmentation import random_walker
    import matplotlib.pyplot as plt
    def visualize(step_title, data, slice_idx=None, window_id=None):
        """Helper function to visualize a single slice of 3D data in a new window."""
        if slice_idx is None:
            slice_idx = data.shape[2] // 2 # Use the middle slice by default
        plt.figure(window_id if window_id else step_title, figsize=(8, 6)) # Ensure a new window
        plt.imshow(data[..., slice_idx], cmap="gray")
        plt.title(step_title)
        plt.axis("off")
        plt.show()
    def save_intermediate(step_num, step_title, data, prefix="intermediate"):
        """Save intermediate 3D data as a NIfTI file and append the step number to the filename."""
        step_filename = f"{prefix}_step{step_num}_{step_title.replace(' ', '_').lower()}.nii.gz"
        img = nb.Nifti1Image(data, np.eye(4))
        img.to_filename(step_filename)
    # Load data
    img = nb.load(in_file)
    data = np.pad(img.get_fdata(dtype="float32"), padding)
    hdr = img.header.copy()
    visualize("Original Image (Padded)", data, window_id="Step 1")
    save_intermediate(1, "original_image_padded", data)
    # Cleanup background and invert intensity
    data[data < np.percentile(data[data > 0], 15)] = 0
    visualize("After Background Cleanup", data, window_id="Step 2")
    save_intermediate(2, "background_cleanup", data)
    data[data > 0] -= data[data > 0].min()
    datainv = -data.copy()
    datainv -= datainv.min()
    datainv /= datainv.max()
    visualize("Inverted Intensity", datainv, window_id="Step 3")
    save_intermediate(3, "inverted_intensity", datainv)
    # Grayscale closing to enhance CSF layer surrounding the brain
    closed = ndimage.grey_closing(datainv, structure=ball(1))
    visualize("After Grayscale Closing", closed, window_id="Step 4")
    save_intermediate(4, "grayscale_closing", closed)
    denoised = ndimage.median_filter(closed, footprint=ball(3))
    visualize("After Median Filtering", denoised, window_id="Step 5")
    save_intermediate(5, "median_filtering", denoised)
    th = threshold_otsu(denoised)
    # Rough binary mask
    closedbin = np.zeros_like(closed)
    closedbin[closed < th] = 1
    visualize("Rough Binary Mask", closedbin, window_id="Step 6")
    save_intermediate(6, "rough_binary_mask", closedbin)
    closedbin = ndimage.binary_opening(closedbin, ball(3)).astype("uint8")
    visualize("After Binary Opening", closedbin, window_id="Step 7")
    save_intermediate(7, "binary_opening", closedbin)
    # Extract largest connected component
    label_im, nb_labels = ndimage.label(closedbin)
    sizes = ndimage.sum(closedbin, label_im, range(nb_labels + 1))
    mask = sizes == sizes.max()
    closedbin = mask[label_im]
    closedbin = ndimage.binary_closing(closedbin, ball(5)).astype("uint8")
    visualize("Largest Connected Component", closedbin, window_id="Step 8")
    save_intermediate(8, "largest_connected_component", closedbin)
    # Prepare markers
    markers = np.ones_like(closed, dtype="int8") * 2
    markers[1:-1, 1:-1, 1:-1] = 0
    closedbin_dil = ndimage.binary_dilation(closedbin, ball(3)) # change ball(5) to ball(1)
    markers[closedbin_dil] = 0
    closed_eroded = ndimage.binary_erosion(closedbin, structure=ball(3)) # change ball(5) to ball(1)
    markers[closed_eroded] = 1
    visualize("Markers for Random Walker", markers, window_id="Step 9")
    save_intermediate(9, "markers_random_walker", markers)
    # Run random walker
    closed[closed > 0.0] -= closed[closed > 0.0].min()
    segtarget = (2 * closed / closed.max()) - 1.0
    labels = random_walker(
    segtarget, markers, spacing=img.header.get_zooms()[:3], return_full_prob=True
    )[..., padding:-padding, padding:-padding, padding:-padding]
    visualize("Random Walker Result (Probability Map)", labels[0, ...], window_id="Step 10")
    save_intermediate(10, "random_walker_result", labels[0, ...])
    out_mask = Path(out_file or "brain_mask.nii.gz").absolute()
    hdr.set_data_dtype("uint8")
    img.__class__((labels[0, ...] >= 0.5).astype("uint8"), img.affine, hdr).to_filename(
    out_mask
    )
    out_probseg = re.sub(
    r"\.nii(\.gz)$", r"_probseg.nii\1", str(out_mask).replace("_mask.", ".")
    )
    hdr.set_data_dtype("float32")
    img.__class__((labels[0, ...]), img.affine, hdr).to_filename(out_probseg)
    out_brain = re.sub(
    r"\.nii(\.gz)$", r"_brainmasked.nii\1", str(out_mask).replace("_mask.", ".")
    )
    data = np.asanyarray(img.dataobj)
    data[labels[0, ...] < 0.5] = 0
    img.__class__(data, img.affine, img.header).to_filename(out_brain)
    visualize("Final Brain Masked Image", data, window_id="Step 11")
    save_intermediate(11, "final_brain_masked", data)
    return str(out_brain), str(out_probseg), str(out_mask)

Here is a sample mask that it generated:

[effigies]

With #3415:

[1RuobingLiu]

With #3415:

Thank you!