Hi @emilymharriott and welcome to neurostars!
Just to be confirm, are you sure that sub-02_dir-AP_epi.nii is the file that was used for the fmap? That is, is there something like a IntendedFor or B0FieldSource/Identifier in the fmap json that is used to link it to the DWI?
If both your DWI and fMAP are in the same phase encoding direction, that could explain why the ācorrectedā results look bad. There really isnāt a way to use that fmap if it is not truly RPE. However, there are alternatives you can look into.
First (and easiest) is using --use-syn-sdc --force-syn --ignore fieldmaps arguments in QSIPrep, which will use a āfieldmaplessā warp (based on the registration between the b0 and intensity-inverted T1 - https://www.frontiersin.org/articles/10.3389/fninf.2017.00017/full)
Second is by trying SynB0-Disco: GitHub - MASILab/Synb0-DISCO: Distortion correction of diffusion weighted MRI without reverse phase-encoding scans or field-maps
Following the instructions on the repo, you would synthesize the RPE b0, name it as you did in the BIDS directory and make a minimal JSON for it that includes the PhaseEncodingDirection of opposite your DWI. It is going to be incorporated into QSIPrep in subsequent release, but is not there yet. Here is a Python wrapper for it:
from dipy.align.imaffine import (AffineMap,
AffineRegistration,
MutualInformationMetric,
transform_centers_of_mass)
from dipy.align.transforms import (AffineTransform3D,
RigidTransform3D)
from dipy.segment.tissue import TissueClassifierHMRF
from scipy.ndimage import gaussian_filter
import logging
import os
import sys
import nibabel as nib
import numpy as np
import warnings
# Disable tensorflow warnings
with warnings.catch_warnings():
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
warnings.simplefilter("ignore")
from dipy.nn.synb0 import Synb0
def register_image(static, static_grid2world, moving, moving_grid2world,
transformation_type='affine', dwi=None, fine=False):
"""
Register a moving image to a static image using either rigid or affine
transformations. If a DWI (4D) is provided, it applies the transformation
to each volume.
Parameters
----------
static : ndarray
The static image volume to which the moving image will be registered.
static_grid2world : ndarray
The grid-to-world (vox2ras) transformation associated with the static
image.
moving : ndarray
The moving image volume that needs to be registered to the static image.
moving_grid2world : ndarray
The grid-to-world (vox2ras) transformation associated with the moving
image.
transformation_type : str, optional
The type of transformation ('rigid' or 'affine'). Default is 'affine'.
dwi : ndarray, optional
Diffusion-weighted imaging data (if applicable). Default is None.
fine : bool, optional
Whether to use fine or coarse settings for the registration.
Default is False.
Raises
------
ValueError
If the transformation_type is neither 'rigid' nor 'affine'.
Returns
-------
ndarray or tuple
If `dwi` is None, returns transformed moving image and transformation
matrix.
If `dwi` is not None, returns transformed DWI and transformation matrix.
"""
if transformation_type not in ['rigid', 'affine']:
raise ValueError('Transformation type not available in Dipy')
# Set all parameters for registration
nbins = 64 if fine else 32
params0 = None
sampling_prop = None
level_iters = [250, 100, 50, 25] if fine else [50, 25, 5]
sigmas = [8.0, 4.0, 2.0, 1.0] if fine else [8.0, 4.0, 2.0]
factors = [8, 4, 2, 1.0] if fine else [8, 4, 2]
metric = MutualInformationMetric(nbins, sampling_prop)
reg_obj = AffineRegistration(metric=metric, level_iters=level_iters,
sigmas=sigmas, factors=factors, verbosity=0)
# First, align the center of mass of both volume
c_of_mass = transform_centers_of_mass(static, static_grid2world,
moving, moving_grid2world)
# Then, rigid transformation (translation + rotation)
transform = RigidTransform3D()
rigid = reg_obj.optimize(static, moving, transform, params0,
static_grid2world, moving_grid2world,
starting_affine=c_of_mass.affine)
if transformation_type == 'affine':
# Finally, affine transformation (translation + rotation + scaling)
transform = AffineTransform3D()
affine = reg_obj.optimize(static, moving, transform, params0,
static_grid2world, moving_grid2world,
starting_affine=rigid.affine)
mapper = affine
transformation = affine.affine
else:
mapper = rigid
transformation = rigid.affine
if dwi is not None:
trans_dwi = transform_dwi(mapper, static, dwi)
return trans_dwi, transformation
else:
return mapper.transform(moving), transformation
def synb0_synthesis(b0, b0_mask, t1, t1_mask, template, out_b0):
logging.info('The usage of synthetic b0 is not fully tested.'
'Be careful when using it.')
# Load template
template_img = nib.load(template)
template_data = template_img.get_fdata()
# Load b0
b0_img = nib.load(b0)
b0_skull_data = b0_img.get_fdata()
b0_mask_img = nib.load(b0_mask)
b0_mask_data = b0_mask_img.get_fdata()
# Load t1
t1_img = nib.load(t1)
t1_skull_data = t1_img.get_fdata()
t1_mask_img = nib.load(t1_mask)
t1_mask_data = t1_mask_img.get_fdata()
# Apply masks
b0_bet_data = np.zeros(b0_skull_data.shape)
b0_bet_data[b0_mask_data > 0] = b0_skull_data[b0_mask_data > 0]
t1_bet_data = np.zeros(t1_skull_data.shape)
t1_bet_data[t1_mask_data > 0] = t1_skull_data[t1_mask_data > 0]
# Crude estimation of the WM mean intensity and normalization
logging.info('Estimating WM mean intensity')
hmrf = TissueClassifierHMRF()
t1_bet_data = gaussian_filter(t1_bet_data, 2)
_, final_segmentation, _ = hmrf.classify(t1_bet_data, 3, 0.25,
tolerance=1e-4, max_iter=5)
avg_wm = np.mean(t1_skull_data[final_segmentation == 3])
t1_skull_data /= avg_wm
t1_skull_data *= 110
# SyNB0 works only in a standard space, so we need to register the images
logging.info('Registering images')
# Use the BET image for registration
t1_bet_to_b0, t1_bet_to_b0_transform = register_image(b0_bet_data,
b0_img.affine,
t1_bet_data,
t1_img.affine,
fine=True)
affine_map = AffineMap(t1_bet_to_b0_transform,
b0_skull_data.shape, b0_img.affine,
t1_skull_data.shape, t1_img.affine)
t1_skull_to_b0 = affine_map.transform(t1_skull_data.astype(np.float64))
# Then register to MNI (using the BET again)
_, t1_bet_to_b0_to_mni_transform = register_image(template_data,
template_img.affine,
t1_bet_to_b0,
b0_img.affine,
fine=True)
affine_map = AffineMap(t1_bet_to_b0_to_mni_transform,
template_data.shape, template_img.affine,
b0_skull_data.shape, b0_img.affine)
# But for prediction, we want the skull
b0_skull_to_mni = affine_map.transform(b0_skull_data.astype(np.float64))
t1_skull_to_mni = affine_map.transform(t1_skull_to_b0.astype(np.float64))
# SynB0
logging.info('Running SyN-B0')
SyNb0 = Synb0()
rev_b0 = SyNb0.predict(b0_skull_to_mni, t1_skull_to_mni)
# Transform back
rev_b0 = affine_map.transform_inverse(rev_b0.astype(np.float64))
# Save out
dtype = b0_img.get_data_dtype()
nib.save(nib.Nifti1Image(rev_b0.astype(dtype), b0_img.affine),
out_b0)
## RUN IT
template = "/PATH/TO/mni_icbm152_t1_tal_nlin_asym_09c_mask_2_5.nii.gz"
t1 = "/PATH/TO/T1w.nii"
t1_mask = "/PATH/TO/T1_mask.nii"
b0 = "/PATH/TO/b0.nii"
b0_mask = "/PATH/TO/b0_mask.nii"
out_b0 = "/PATH/TO/RPE_b0.nii"
synb0_synthesis(b0, b0_mask, t1, t1_mask, template, out_b0)
Best,
Steven