Brain Imaging & Analysis Center

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== EPI Correction with polarity images using TOPUP ======

Steps to setup and apply [[https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide/|FSL Topup]] correction on EPI data collected here with the use of your polarity images

==== Create your acquisition parameters ====

There will be two series of single timepoint EPI images.  You can grab relevant info to create the acq_params.txt files from your XML header.

You need to calculate readout time in seconds ( the physical time it takes to get the acquisition matrix of a single slice ) and get the polarity direction ( phase encode direction ).

the readout time in seconds for the parameter file will be:
<code>readout = (echospacing * (acquisitionmatrix[0] * (percentsampling/100))) / 1e6 

echospacing in the BXH header is in microseconds
</code>


the polarity for the entry will have to be determined from the **seriesdescription**, which is typically "field map reverse polarity" or "field map regular"

reverse will be "-1" in the acq_params.txt and regular will be "1"

<code>
 echo '0 1 0 readout' > acq_params.txt
 echo '0 -1 0 readout' >> acq_params.txt

example output:
   0 1 0 0.077312
   0 -1 0 0.077312

</code>

==== Create your Blip Up/Down image ====

Merge the two polarity images together into a 4D file with fsl_merge. Put them in the order of your acq_params.txt file.  In the example above, it would be regular, then reversed.

<code>
fslmerge -t bud ../bia6_00186_006.nii.gz ../bia6_00186_007.nii.gz
    
tesla:00273 cmp12$ fslinfo bud.nii.gz 
    data_type      INT16
    dim1           128
    dim2           128
    dim3           70
    dim4           2
</code>

{{:biac:analysis:epi_polarity.png?800|}}

==== Run topup to calculate the field and coefficients images for correction ====

<code>
#run popup calculate the field/coef
topup --verbose --imain=bud --datain=acq_params.txt --config=$FSLDIR/src/topup/flirtsch/b02b0.cnf --out=rs_topup --fout=topup_field
</code>

==== Now apply the correction to your functional runs ====
**inindex** references the "regular" image index from your acq_params.txt file, which is the same phase encode direction of your functional runs 

<code>
#apply it
applytopup --imain=../bia6_00186_008_01.nii.gz --inindex=1 --method=jac --datain=acq_params.txt --topup=rs_topup --out=run008 --verbose
applytopup --imain=../bia6_00186_009_01.nii.gz --inindex=1 --method=jac --datain=acq_params.txt --topup=rs_topup --out=run009 --verbose
applytopup --imain=../bia6_00186_010_01.nii.gz --inindex=1 --method=jac --datain=acq_params.txt --topup=rs_topup --out=run010 --verbose
</code>

The top two images are an example of the polarity images after correction, which is not necessary.
The bottom two images are the mean functional image of a series uncorrected, then corrected by the topup results.

{{:biac:analysis:epi_corrected.png?800|}}

If you want to re-generate a valid BXH for the corrected data, you can generate a new BXH and merge in the original info

<code>
fslwrapbxh run008
mv run008.bxh tmp.bxh
bxh_merge ../bia6_00186_008_01.bxh tmp.bxh run008.bxh 
rm tmp.bxh
</code>



---------------

====== DWI Correction with RPE B0s using TOPUP ======

Correcting DWI images using reverse phase encoded B0s is almost the same as the above.  Typically you'll have a couple B0s in your normal DWI acquisition, with a separate series of the same number of B0s collected with the phase encoding direction reversed.

==== Create your acq_params.txt file ===

the readout time in seconds for the parameter file will be:
<code>readout = (echospacing * (acquisitionmatrix[0] * (percentsampling/100))) / 1e6 

echospacing in BXH header is in microseconds
</code>

the regular B0s will get the "1" and the reverse phase B0s will get the "-1"

<code>
 echo '0 1 0 [readout]' > acq_params.txt
 echo '0 1 0 [readout]' >> acq_params.txt
 echo '0 -1 0 [readout]' >> acq_params.txt
 echo '0 -1 0 [readout]' >> acq_params.txt

example output:
   0 1 0 0.10656
   0 1 0 0.10656
   0 -1 0 0.10656
   0 -1 0 0.10656
</code>

==== Create your Blip Up/Down image ====

First, you'll need to extract your B0s from the regular DWI acquisition.

If you look at the gradient directions within the DWI acquisition the B0s will be "0 0 0".

Within the BXH the first 2 directions are B0:
<code>
<datapoints label="diffusiondirection">
 <value>0 0 0</value>
 <value>0 0 0</value>
 <value>0.707107 0 0</value> <-- not a B0
</code>

Extract them with bxhselect:
<code>
bxhselect --timeselect 0:1 bia6_00197_012.bxh bu
bxhselect --timeselect 0:1 bia6_00197_013.bxh bd
</code>

This will create bu.nii.gz / bd.nii.gz with the first 2 gradient directions ( the B0s )

Now merge the extracted B0s from the regular acquisition with the images from the reversed acquisition.  Merge these in the same order of the acq_params.txt file.
<code>
fslmerge -t bud bu.nii.gz bd.nii.gz

tesla:HCP cmp12$ fslinfo bud.nii.gz
data_type      INT16
dim1           256
dim2           256
dim3           92
dim4           4
</code>

{{:biac:analysis:dwi_polarity.png?600|}}

==== Run topup to calculate the field and coefficients images for correction ====

<code>
#run popup calculate the field/coef
topup --verbose --imain=bud --datain=acq_params.txt --config=$FSLDIR/src/topup/flirtsch/b02b0.cnf --out=dwi_topup --fout=topup_field
</code>

==== Now apply the correction to your DWI ====

**inindex** references the "regular" dwi B0s index from your acq_para
ms.
txt file

<code>
#apply it
applytopup --imain=bia6_00197_012.nii.gz --inindex=1 --method=jac --datain=acq_params.txt --topup=dwi_topup --out=dwi_corrected --verbose
</code>

{{:biac:analysis:dwi_corrected.png?600|}}

If you want to re-generate a valid BXH for the corrected data, you can generate a new BXH and merge in the original info

<code>
fslwrapbxh dwi_corrected.nii.gz
mv dwi_corrected.bxh tmp.bxh
bxh_addgradients --fsl tmp.bxh bvecs bvals dwi_corrected.bxh
rm tmp.bxh
</code>


---------------

====== DWI Correction with MRTRIX3 and RPE Data ======

[[https://mrtrix.readthedocs.io/en/latest/|mrtrix3]] has a nice wrapper for FSL's eddy/topup/apply_topup that can do the above corrections and along with eddy corrections

If you are running dwidenoise, do it BEFORE dwipreproc. 

==== Calculate your readout time ====
<code>readout = (echospacing * (acquisitionmatrix[0] * (percentsampling/100))) / 1e6 

echospacing in BXH header is in microseconds
</code>

==== Prepare your input datasets ====

there are 2 scenarios that typically apply here
1) DWI and a short acquisition with RPE B0s
2) 2 entire acquisitions with the same gradient table and reverse phase encoding 

In both scenarios it is important to create your data with normal phase encoding direction first, followed by reversed.

===== Scenario 1 =====

==== Prepare your input datasets ====

Specify that a set of images (typically b=0 volumes) will be provided for use in inhomogeneity field estimation only

For scenario 1, create your blip up / blip down B0 data the same way as above.
<code>
bxhselect --timeselect 0 bi
a6_00197_012.bxh bu 
bxhselect --timeselect 0 bia6_00197_013.bxh bd
fslmerge -t bud bu.nii.gz bd.nii.gz
rm bu.* bd.*
</code>

Prepare your dwi series for input to mrtrix3.  If you use their mif data format, then the gradient tables will be automatically passed throughout their software with additional flags.

<code>
#pull the gradient table from BXH in FSL format
extractdiffdirs --fsl bia6_00414_013.bxh bvecs bvals

#convert nii.gz data into mif while inserting the gradients
mrconvert -fslgrad bvecs bvals bia6_00414_013.nii.gz dwi.mif 

#export the table back out in their format 
mrinfo dwi.mif -export_grad_mrtrix grads.b

#put THOSE back into the mif header
mrconvert -grad grads.b dwi.mif dwi.mif -force

</code>

The last step seems redundant, but there is an issue with how mrtrix3 tools are handling the multi-shell GE data produced on our scanners.  If you data is NOT multi-shell, then you don't need it.  If your data IS multi-shell, then putting it back in can workaround other issues until their RC4 code is released.  [[https://github.com/MRtrix3/mrtrix3/issues/1487|discussion here]]

==== Run dwipreproc -rpe_pair ====

Run [[https://mrtrix.readthedocs.io/en/latest/reference/scripts/dwipreproc.html|dwipreproc]] with the "-rpe_pair" option

<code>
dwipre
proc dwi.mif dwi_corr.mif -rpe_pair -se_epi bud.nii.gz -pe_dir AP -readout_time 0.10656  -debug

-rpe_pair specifies you're providing a pair of B0s ( regular, reversed )
-se_epi is you bud image
-pe_dir is your dwi's phase encode direction
-readout_time is from calculation above
</code>
 
dwipreproc will run topup, eddy and apply topup.  eddy will output update diffusion directions and insert them into your dwi_corr.mif output header.

You can generate a BXH header for a mif file with:
<code>
mif2bxh dwi_corr.mif dwi_corr.bxh
</code>

===== Scenario 2 =====

You have an entire acquisition with the same gradient table as your DWI sequence, but with reverse phase encoding.  This information will be used to perform a recombination of image volumes (each pair of volumes with the same b-vector but different phase encoding directions will be combined together into a single volume).

Calculate the readout time in the same way as above

==== Prepare your input datasets ====

This time you'll concatenate the 2 DWI acquisitions.  regular, then reversed 

<code>
#concat the 2 series ( regular, reversed ) 
bxh_concat bia6_00260_00
7_LAS.bxh bia6_00260_008_LAS.bxh both

#extract the gradients
extractdiffdirs --fsl both.bxh bvecs bvals

#convert to mif
mrconvert -fslgrad bvecs bvals both.nii.gz dwi.mif

#re-extract table
mrinfo dwi.mif -export_grad_mrtrix grads.b

#put corrected table back in
mrconvert -grad grads.b dwi.mif dwi.mif -force

#remove the temporary files ( bvecs, bvals, both.* )
</code>


==== Run dwipreproc -rpe_all ====

Run dwipreproc with the "-rpe_all" option

<code>
dwipreproc dwi.mif dwi_corr.mif -rpe_all -pe_dir AP -readout_time 0.10656 -debug

-rpe_all signals that you've replicated 
A
LL the directions with a rpe acquisition
-pe_dir is the phase encode direction of your regular acquisition
-readout_time from above  
</code>


{{ :biac:analysis:trace_raw_corr.png?800 |}}
Brain Imaging & Analysis Center

User Tools

Site Tools

Page Tools
