EPI images can be quite distorted due to B0 inhomogeneities, which in the brain are especially pronounced near the sinuses. While signal lost to these inhomogeneities cannot be recovered, signal that has been merely displaced can be returned to its proper location in the image if a map of the B0 field has been acquired. This document discusses two methods used at the center to acquire images that can be used to create field maps, and how to use tools in FSL and SPM to create and apply field maps. This document focuses on the practicalities of how to generate and use the field maps.
Make sure you read this part
Applying distortion correction should make your images look LESS distorted, not more. If after applying distortion correction your images look worse, or about the same, you should look for a mistake somewhere in your processing.
Another really important thing: please don’t take scripts written by someone else in your lab for their project and just assume everything will work. Your acquisition parameters could very well be different from theirs. Things do change.
I already took the data, what kind of field map do I have?
If you are asking this question, the answer is that you don’t have a field map yet. You have data that you will use to calculate a field map. This article covers how to do that for two different methods. Method 1 calculates a field map based on the difference in phase between two different echos in a double echo sequence. Method 2 uses two separate acquisitions with opposite phase encoding directions to calculate a field map based on the difference in distortion between the two acquisitions. If you are an LCNI user, we use both methods. If you don’t know what kind of data you acquired, it should be clear once you look at the results: if you used the stock double-echo gradient echo field map sequence (method 1), you’ll have two series. The first will contain two magnitude images, one for each echo. The second series will contain a single phase difference image, a subtraction of the two phase images from each echo.
If you used the se-epi with opposite phase encode method (method 2), you will also get two series. Each will contain a single magnitude volume, and one will have the opposite phase encoding direction from the other.
Effective echo spacing and total readout time
There are two parameters you will frequently need when calculating and applying field maps: the effective echo spacing and the total readout time for an EPI sequence. We use the “effective” echo spacing, rather than the actual echo spacing, in order to include the effects of parallel imaging, phase oversampling, etc. Siemens has helpfully included this information in the DICOM file with the parameter “Bandwidth Per Pixel Phase Encode”, stored in DICOM tag (0019, 1028). Multiplying this by the size of the reconstructed image in the phase direction (aka the number of reconstructed phase lines) will give you the reciprocal of the effective echo spacing:
Effective Echo Spacing (s) = 1/(BandwidthPerPixelPhaseEncode * MatrixSizePhase)
The size of the image in the phase direction is *usually* the first number in the field (0051, 100b), AcquisitionMatrixText. This is, however, not always the case, so it is better to determine the actual size (rows or columns) of your reconstructed image in the phase encode direction. You can also just read the effective echo spacing from the Series Info screen in MRIConvert if you are running the latest version.
The total readout time (FSL definition) is the time from the center of the first echo to the center of the last:
Total readout time (FSL) = (number of echoes - 1) * echo spacing
In the case of a simple EPI sequence, with no parallel imaging, no partial fourier, no phase oversampling, etc, this is straightforward enough. For the general case, it’s appropriate to consider this an “effective” readout time, calculated using the effective echo spacing and the number of reconstructed phase lines:
Total readout time (FSL) = (MatrixSizePhase - 1) * EffectiveEchoSpacing
For the majority of scans the length of a single echo is so short compared to the readout time that you can use the SPM definition of total readout time. SPM defines it to be from the beginning of the first echo to the end of the last echo. This is just the reciprocal of the Bandwidth Per Pixel Phase Encode:
Total Readout Time (SPM) = 1/(BandwidthPerPixelPhaseEncode)
Since the Bandwidth Per Pixel Phase Encode is in Hz, this will give the readout time in seconds. MRIConvert reports both definitions, and you should find them to be within 2% of each other for most sequences.
Acquiring a field map, method 1: double echo
This method calculates a field map from a double-echo gradient echo sequence. The difference in phase between the two images is proportional to the difference in echo times and the B0 inhomogeneity. The field map is calculated by taking the difference between the two phase images, and dividing that by the echo time difference.
Setting up the sequence
The Skyra conveniently comes with a field map sequence installed. The echo time difference is set to a specific value in order to minimize the effect of the fat-water chemical shift, so don’t change it unless you know what you’re doing. The field map should be in the same orientation as the EPIs you are planning to correct (eg, axial), but does not need to have exactly the same geometry, field of view, or resolution. Make sure that you save both the magnitude images (there will be two, one for each echo) and the phase difference image.
Calculating the field map
FSL has a handy utility for converting a Siemens phase difference image to a field map: fsl_prepare_fieldmap. The syntax is:
fsl_prepare_fieldmap <scanner> <phase_image> <magnitude_image> <out_image> <deltaTE (in ms)>
The magnitude image should be brain extracted, and the brain extraction should be tight. This is discussed in more detail on the FSL wiki if you want more information. You can use either echo from the fieldmap aquisition, or average them together. “deltaTE” will be 2.46 unless you changed the echo times in the fieldmap acquisition. The output image will be in rad/s, ready for use in FEAT or FUGUE.
Acquiring a field map, method 2: opposite phase encoded EPIs
This method acquires two or more echo planar images with opposite phase encode directions. More info on how it uses this information will be added later.
Setting up the sequence
If you are doing this with the multiband sequence, you probably used the SE-EPI. More details on this will be added later, but you should have two series in the end, each with opposite phase encoding. The rest of this document assumes those directions are “AP” and “PA”.
Calculating the field map
You will need to use the FSL tool “topup” to create a field map for use with FUGUE or FEAT. Detailed documentation on topup is available on the FSL wiki, but here is an abbreviated list of necessary steps that should work for most users.
IMPORTANT WARNING: Until FSL release 5.0.9, topup up did not write any orientation information into the outputted nifti header files. The consequences of this omission range from minor to disastrous. Be sure you are using the latest release of FSL, or fix the orientation information in the nifti files after the fact.
1. Use fslmerge to combine the two se-epi images into one 4D file: fslmerge -t [output] [inputAP] [inputPA].
2. Create a “datain” file. This is a text file with four columns. The first three columns indicate the phase encode direction, and the fourth is the total readout time (defined as the time from the center of the first echo to the center of the last) in seconds for the se-epi acquisition. The total readout time is given by the (number of echoes minus one) * echo spacing. See the above discussion on total readout time for how to calculate this value. Remember, this is all from the se-epi sequence, NOT the functional sequence you will be applying the fieldmap to later. When in doubt, ask your local physicist for help. The correct values for the first three columns may change depending on how the files were converted to NIFTI format. You’ll wind up with a file that looks something like this:
0 -1 0 .0351 // This is just an example! 0 1 0 .0351 // Don't just copy these values!
3. Run topup. You can use the default config file b02b0.cnf (you do not need to create this file). The command will look something like this:
topup --imain=se_epi_merged --datain=datain.txt --config=b02b0.cnf --fout=my_fieldmap --iout=se_epi_unwarped
4. The fieldmap will be in units of Hz, but you’ll need it in units of rad/s for use with fugue or feat, so multiply by 2pi:
fslmaths my_fieldmap -mul 6.28 my_fieldmap_rads
5. Finally, if you will be using this with Feat, you’ll need a single magnitude image, and a brain extracted version of that image. You can get that this way:
fslmaths se_epi_unwarped -Tmean my_fieldmap_mag bet2 my_fieldmap_mag my_fieldmap_mag_brain
Make sure the brain extracted image has the same name as the non brain extracted image, but with _brain at the end.
Applying distortion correction: FSL
In Feat, you’ll need to specify the field map image (the one you created with units of rad/s) and the brain extracted magnitude image from the field map acquisition. Feat will look for a non brain extracted magnitude image in the same directory and having the same name, but without _brain at the end, so make sure it’s there. For fieldmaps acquired with the double echo method, it’s very important that the brain extracted image excludes all non-brain regions, as the field map will be quite noisy in those areas. For fieldmaps acquired with the se-epi method this is less vital.
Next you will need the TE from the functional epi series you are processing (not the field map!), and the unwarp direction (-y for AP phase encoding and “FSL NIFTI” output type from MRIConvert). Finally, you’ll need the effective EPI echo spacing in ms (again, from the functional EPI series, not the field map). See the above discussion for how to get the effective echo spacing from the DICOM file.
Calculating and Applying Fieldmaps with SPM
First, let’s use the SPM gui to do one fieldmap, acquired using method 1, for one subject. We’ll cover batch mode later. SPM uses the FieldMap Toolbox to calculate a voxel displacement map that is used to correct the BOLD epi images. The parameters for this can be changed in the GUI, or in a defaults file. There are a bunch of default files already, none of which are likely to be correct for you, but you can save your own. The first parameters you’ll need to change are the short and long echo times.
pm_def.SHORT_ECHO_TIME = 4.37; pm_def.LONG_ECHO_TIME = 6.83;
These are the two echo times in the Siemens double echo stock fieldmap sequence. The values I’ve given here are probably the same ones you used, but double check by looking in the DICOM file (MRIConvert will show you these values in the “series info” screen). Hit “load phase” and select the phase difference image from your field map acquisition. When you are asked “Scale this to radians?” choose yes. Then hit “load mag” and choose one of the two magnitude images from the field map sequence. “Mask brain” can be either “yes” or “no”. Finally, choose “calculate”. A fieldmap should appear in the graphics window.
If your fieldmap was created using method 2, you will have to use FSL & topup to calculate a field map as discussed above. In that case, you will need to load a precalculated fieldmap using the obvious button. The fieldmap needs to already be in Hz, must have a filename starting with “fpm”, and must be in a NIFTI pair file format.
The lower window has information about the EPI image you wish to unwarp. This is where we need to total readout time, which is also in the defaults file like this:
pm_def.TOTAL_EPI_READOUT_TIME = 21.1 // Replace with the right value for your data!
This is the total epi readout time in ms for the BOLD series you are correcting, NOT the total readout time from the field map acquisition. See the above discussion on total readout time for how to get this from the DICOM. “EPI based fieldmap” should be “no”, and “apply jacobian modulation” should also be “no”. Polarity of the blips depends not only on your acquisition but also on how your data was converted to NIFTI, so it could be either + or – ve. The easiest thing to do is to try it both ways. Once these parameters are set, load an EPI image. After a moment, you’ll see the fieldmap, warped EPI, and unwarped EPI in the graphics window.
As I stated way up at the beginning, the unwarped EPI should look objectively “better” than the warped. If you see that a big piece of brain is now missing, it’s because the voxel displacement map is not lined up with the EPI.
You should be able to fix this using the “match VDM to EPI”. Sometimes, however, this doesn’t work, and now I’m afraid I have to get a bit long-winded.
The phase difference between the two echos depends on the local field — that’s how field mapping works. We create a fieldmap that maps out this field, in units of Hz. This includes not only susceptibility effects, but also global drifts in frequency from the nominal center value. If a frequency adjustment is not performed immediately before running the fieldmap sequence, the difference image will include this offset. You may notice that some phase difference images look a little different from the one above:
There’s nothing wrong mathematically, they just aren’t “centered” at the correct “zero”. FSL, in the “fsl_prepare_fieldmaps” script, includes a final demeaning step to account for this. As far as I can tell, SPM does not (or at least does not do it the same way), and it is not able to handle large offsets from the nominal center frequency. It is easy to check this in the fieldmap toolbox gui: after loading an EPI image, switch the polarity of the phase-encode blips. Areas that are not distorted should not change. If, however, the entire brain “jumps” up and down, you have a global offset. If you have a global offset, the voxel displacement map will be wrong and will not unwarp the EPI correctly.
One way of fixing this problem is to calculate your fieldmaps in FSL and load them into SPM. Another is to demean the fieldmaps in SPM by adding the line fm.upm = fm.upm – median(fm.upm(:)); to the file pm_make_fieldmap.m. In SPM12 you’ll want to add it at line 190:
if flags.pad ~= 0 npxs = pxs/pxs(1); ks = floor((flags.pad+1)*[1 1 1]./npxs); if ~mod(ks(1),2) ks(1) = ks(1)+1; end if ~mod(ks(2),2) ks(2) = ks(2)+1; end if ~mod(ks(3),2) ks(3) = ks(3)+1; end kernel = ones(ks); [fm.upm,fm.mask] = pm_pad(fm.upm,fm.mask,kernel); end fm.upm = fm.upm - median(fm.upm(:));
Note that I am still in the process of testing this and no promises are made! Once you have everything all figured out, you’ll probably be using batch mode in SPM instead of the gui. I’m assuming that’s straightforward, but if I run into trouble when I try doing it myself I’ll add that info here.
Please contact Jolinda if you are a LCNI user with more questions about any of this.