distortion correction using field maps in AFNI

AFNI Message Board
1

Hello,

I need to re-analyze a multiband multi-echo dataset (previously analyzed with FSL) using AFNI. We only have field maps for this dataset, so I have been trying to find a way to do fieldmap-based distortion correction with AFNI and integrate it with our current afni_proc pipeline. I know this topic has been brought up several times on the message board, but it seems so far there has not been any official solution integrated with afni_proc.

I found this repository on GitHub (https://github.com/nih-fmrif/bids-b0-tools) which implements B0-correction using AFNI tools (‘distortionFix.py’). Has anyone used this repository? Has anyone integrated ‘distortionFix.py’ with afni_proc.py? If so, how?

Thank you,
Golbarg

2

Hello,

I am responsible for most of the content of this repository, and got guidance implementing the distortion correction tools from various members of the AFNI team.

The B[sub]0[/sub] correction module using AFNI tools runs from approximately lines 250 to 300 in the “distortionFix.py” file.

The input to this module are a BIDS ID for 1 subject’s data for a given session, the magnitude image from a B[sub]0[/sub] field map acquisition, the frequency map from the B[sub]0[/sub] field map acquisition (essentially a scaled phase map - more details below), a mask of the B[sub]0[/sub] (usually generated by running 3dAutomask on the magnitude B[sub]0[/sub] data), the effective echo spacing between lines of EPI acquisition (in seconds), and the EPI phase FOV (in mm).

The first step here is to take the information in the frequency map, scale it “appropriately”, according to the following formula:

Frequency Map * EPI Echo Spacing (in seconds) * EPI Phase FOV (in mm)

The field map is then smoothed using gaussian blurring.

The scaled and smoothed field map is then used to “self distort” - so that its distortions now match the EPI distortions. This is the “fmapInHz-smoothed-warped2epi_” data set.

Finally, this distorted field map is used to unwarp the EPI data.

In addition to the timing and geometric data put in above, you would also need to know your phase encoding direction. As far as I know, the vendors don’t have a standard way to denote +y / -y / +AP / -AP. I would say try both orientations on a handful of subjects - it should be pretty apparent which is correct from the visual changes in distortion. Once you’ve figured this out, and the acquisition for your study remained consistent, this should “just work”.

As to integration with the AFNI processing stream, I’ve not done this work. The distorted field map generated above is essentially another warp field that should be able to be concatenated with other warps, but I think other members of the AFNI team can better answer that particular part of your question.

Glad you found this repository - and hope you can use the methods in there!

Vinai

3

Hi Vinai,

Thanks a lot for the helpful reply.

I noticed that the scaled frequency map is in units of Hz rather than rad/s, is this because you convert the units of the phase image to Hz to get a frequency map? We collected our data using a Siemens Prisma scanner. Is “distortionFix.py” written for GE data?

Does scaling and smoothing of the frequency map replicate what FSL function “fsl_prepare_fieldmap” does to calibrate the phase image?

What is the rational for setting “1blur_sigma” to 9 mm?

How do you calculate the unwarping shift map? (i.e., an image similar to feat output “FM_UD_fmap2epi_shift” found in reg/unwarp)

Do you know whether any output in reg/unwarp could be used directly with 3dNwarpApply to unwarp the EPI data? I appreciate your thoughts and suggestions.

Thank you,
Golbarg

4

Hello,

The unit of Hz/s was chosen as it would give an easy way to scale the field offset, without having to pass along the additional parameter of the echo-time difference used to generate the B[sub]0[/sub] map. This was not for “GE data.” This was a custom sequence written by me to generate B[sub]0[/sub] map data. Frequency to radians (per second) is a trivial scaling (lines 352 - 357 in this file). I am not sure what Siemens’ B[sub]0[/sub] sequence produces - I thought it was the phase difference between the 2 echoes. In that case, you would need to use ‘3dcalc’ to scale that data appropriately to frequency offset, so the afniB0 module in distortionFix.py can be used directly, or to RPS …

The choices made for smoothing and generally dealing with the B[sub]0[/sub] maps were the results of several weeks of iterations and trying to match a set of parameters that would be given to FSL’s distortion correction pipeline, so that field maps would be regularized / smoothed to similar extents by both pipelines. You can take a look at the “fslB0” function in the same file to see what parameter choices were passed along to fugue. If you are accustomed to using different field map smoothings with FSL/fugue, then you would need to figure out the appropriate mappings of your favorite parameter sets to AFNI commands.

I am not sure about using reg/unwarp outputs, primarily b/c I don’t know what these are. If it is a warp field, then it is technically possible that they might be usable in AFNI by 3dNWarpApply. But you would need to ensure the outputs of FSL are consistent (from geometry, orientations, scaling, etc) would be consistent with what AFNI’s commands expect.

Cheers!

5

Hello,
Thank you for the helpful reply.
In the past, we applied distortion correction by Feat pipeline rather than calling Fugue directly, which would generate the folder I mentioned in my previous post (e.g., run1.feat/reg/unwarp).

I re-wrote afniB0 module for a quick test on our data. I’d appreciate it if you could take a look at my steps since the warped field map image does not look right. I just skipped smoothing since no smoothing was done when Fugue was called via Feat pipeline. The phase encoding direction is A >> P, and we got undistorted EPIs by setting the unwarping direction to “-y” in FSL.
Thank you,
Golbarg

Convert field map phase image (in units of rad) to rad/s (echo time difference on Siemens= 2.46ms)

3dcalc -a gre_fieldmap_phase -expr ‘a/(0.00246)’ -datum float -prefix fmap_radsec

Convert rad/s to Hz

3dcalc -a fmap_radsec -expr ‘a * 1/(2*PI)’ -datum float -prefix fmapInHz

Scale the frequency map

3dcalc -a fmapInHz -b mag_brain_mask_ero -expr ‘a * 0.00049 * 205 * b’ -datum float -prefix fmapInHz_scaled

Distort scaled field map in phase encoding direction (A >> P)

3dNwarpApply -warp AP:1:fmapInHz_scaled -prefix fmapInHz_scaled_warped2epi -source fmapInHz_scaled

Unwarp EPI using warped field map (unwarp direction: -y)

3dNwarpApply -warp AP:-1:fmapInHz_scaled_warped2epi -prefix echo2_fixed -source echo2

6

This looks good. From everything here, the only thing I’ve seen missing is some kind of field map smoothing / regularization. This was accomplished in our GitHub repository using the ‘3dmerge’ command in the ‘afniB0’ function in the ‘distortionFix.py’ file.

If your input phase maps above (i.e. ‘gre_fieldmap_phase’) are already smoothed/regularized, this is probably not as critical. Where it might show problems is at the edges of the brain mask. The blurring option in ‘3dmerge’ will blur the field map at the edges of the mask, and should improve the behavior of the warping at those points. With the sharp edges in image masks, warping can generate … ‘interesting’ artifacts …

Just one other comment about your echo spacing above (490 µs). If your EPI acquisition is accelerated in-plane (with either SENSE or GRAPPA), then you would have to divide this by your in-plane acceleration factor (R). Mutli-band acquisition does not affect this. So if you used R = 2, then your effective echo spacing would be 490 µs / 2 = 245 µs. If you have not used any kind of in-plane acceleration (such as the Human Connectome or ABCD protocols), then the way you’ve specified the above scaling is correct.

7

Thank you for the feedback.
I smoothed the field map using 3dmerge, but the distorted field map image still looks very wrong. Could you take a look at the images? I’ll send you those shortly.

Thank you,
Golbarg

8

Hello

I thought a bit more about the images you sent, and then it struck me about how the warped B[sub]0[/sub] field maps overlaid on the original B[sub]0[/sub] field map.

I believe that the output field map from Siemens is not in radians. It is in scaled to span the full range of image units in Siemens DICOM (I believe either 0 - 4096, or -2048 to 2047. Someone else should confirm this).

If you look at the “siemens_process” function in “fsl_prepare_fieldmap” script – you can see that to convert the Siemens generated B[sub]0[/sub] map to radians, there is a line calling “fslmaths” that does:

$FSLDIR/bin/fslmaths ${newphaseroot} -div 2048 -sub 1 -mul 3.14159 -mas ${maskim} ${tmpnm}_tmp_ph_radians -odt float

So if I am reading this correctly, to convert the Siemens B[sub]0[/sub] field map to Hz, you would run, on that field map:

3dcalc -a siemens_b0_map.nii -expr ((a / 2048 – 1)) / ∆TE

It seems that Siemens computes their output field map by dividing by 2 * PI, then multiplying by 4096 – to span the image range. So to get it back in radians, you would need to reverse this. And since you want to convert from radians/s to Hz, there’s no need to multiply and then divide by PI – canceling those factors out.

You should be able to combine multiple 3dcalc commands to avoid keeping track of multiple intermediate files, so combining this with the commands you previously wrote, you should have something like:

3dcalc -a siemens_b0_map.nii -b mag_brain_mask_ero -expr ‘(a / 2048 - 1) * 0.00049 * 205 * b / ∆TE’ -datum float -prefix fmap_masked_and_scaled

I am not so familiar with working with Siemens’ B[sub]0[/sub] field maps, so I did not catch this. But try this, and see if the processed and warped field maps now look closer to the shape of the brain in the EPI data. I think this will get you closer to where you need to be.

9

Hello,

Thank you for the helpful post.
You are correct. Based on this webpage, Siemens field map data seems to be in the range -4096 to +4096
https://www.mccauslandcenter.sc.edu/crnl/tools/fieldmap

I tried the following
1- converted the prepared field map image (in units of rad/s, generated by “fsl_prepare_fieldmap” command) to Hz, then applied scaling, smoothing and distortion

2- ran 3dcalc on the original field map phase image as you suggested, then applied smoothing and distortion
3dcalc -a gre_fieldmap_phase -b mag_brain_mask_ero -expr ‘(a / 4096 - 1) * 0.00049 * 205 * b / ∆TE’ -datum float -prefix fmap_scaled

The new distorted field map images make more sense; but they look somewhat different. Also, FOV is rotated and there is banding in the images. I’d appreciate it if you could take a look. I’ll send you the images shortly.

Thank you,
Golbarg