Brain Imaging & Analysis Center

User Tools

Site Tools

You are here: biac » biac » analysis » topup_correction
biac:analysis:topup_correction

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision		Previous revision
Next revisionBoth sides next revision
biac:analysis:topup_correction [2019/02/19 16:05] – created cmp12biac:analysis:topup_correction [2019/02/22 14:31] – cmp12
Line 1: Line 1:
 ====== EPI Correction with polarity images using TOPUP ====== ====== 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) a short acquisition with RPE B0s
 +2) an entire acquisition 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_007_LAS.bxh bia6_00260_008_LAS.bxh both
 +
 +#extract the gradients
 +extractdiffdirs --fsl both 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 ALL 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?600 |}}
  
biac/analysis/topup_correction.txt · Last modified: 2024/06/21 15:44 by 127.0.0.1