From c07bc573889dc95d25a352d87818e05af694687f Mon Sep 17 00:00:00 2001 From: Timothy Brown Date: Tue, 18 Nov 2014 15:30:02 -0600 Subject: [PATCH] Integrated code to support General Electric Gradient Echo field map --- ...enericfMRIVolumeProcessingPipelineBatch.sh | 71 ++-- .../Scripts/PreFreeSurferPipelineBatch.sh | 209 +++++++++--- PreFreeSurfer/PreFreeSurferPipeline.sh | 257 ++++++++------ .../T2wToT1wDistortionCorrectAndReg.sh | 226 ++++++++---- .../GenericfMRIVolumeProcessingPipeline.sh | 40 ++- ...IToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh | 323 +++++++++++------- ...GeneralElectricFieldMapPreprocessingAll.sh | 127 +++++++ ....sh => SiemensFieldMapPreprocessingAll.sh} | 5 +- 8 files changed, 871 insertions(+), 387 deletions(-) create mode 100755 global/scripts/GeneralElectricFieldMapPreprocessingAll.sh rename global/scripts/{FieldMapPreprocessingAll.sh => SiemensFieldMapPreprocessingAll.sh} (97%) diff --git a/Examples/Scripts/GenericfMRIVolumeProcessingPipelineBatch.sh b/Examples/Scripts/GenericfMRIVolumeProcessingPipelineBatch.sh index dff2253c5..ce1cc8601 100755 --- a/Examples/Scripts/GenericfMRIVolumeProcessingPipelineBatch.sh +++ b/Examples/Scripts/GenericfMRIVolumeProcessingPipelineBatch.sh @@ -65,33 +65,50 @@ PRINTCOM="" ########################################## INPUTS ########################################## -#Scripts called by this script do NOT assume anything about the form of the input names or paths. -#This batch script assumes the HCP raw data naming convention, e.g. for tfMRI_EMOTION_LR and tfMRI_EMOTION_RL: - +# Scripts called by this script do NOT assume anything about the form of the input names or paths. +# This batch script assumes the HCP raw data naming convention. +# +# For example, if phase encoding directions are LR and RL, for tfMRI_EMOTION_LR and tfMRI_EMOTION_RL: +# # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_LR/${Subject}_3T_tfMRI_EMOTION_LR.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_LR/${Subject}_3T_tfMRI_EMOTION_LR_SBRef.nii.gz - +# # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_RL/${Subject}_3T_tfMRI_EMOTION_RL.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_RL/${Subject}_3T_tfMRI_EMOTION_RL_SBRef.nii.gz - +# # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_LR/${Subject}_3T_SpinEchoFieldMap_LR.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_LR/${Subject}_3T_SpinEchoFieldMap_RL.nii.gz - +# # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_RL/${Subject}_3T_SpinEchoFieldMap_LR.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_RL/${Subject}_3T_SpinEchoFieldMap_RL.nii.gz - -#Change Scan Settings: Dwelltime, FieldMap Delta TE (if using), and $PhaseEncodinglist to match your images -#These are set to match the HCP Protocol by default - -#If using gradient distortion correction, use the coefficents from your scanner -#The HCP gradient distortion coefficents are only available through Siemens -#Gradient distortion in standard scanners like the Trio is much less than for the HCP Skyra. - -#To get accurate EPI distortion correction with TOPUP, the flags in PhaseEncodinglist must match the phase encoding -#direction of the EPI scan, and you must have used the correct images in SpinEchoPhaseEncodeNegative and Positive -#variables. If the distortion is twice as bad as in the original images, flip either the order of the spin echo -#images or reverse the phase encoding list flag. The pipeline expects you to have used the same phase encoding -#axis in the fMRI data as in the spin echo field map data (x/-x or y/-y). +# +# If phase encoding directions are PA and AP: +# +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_PA/${Subject}_3T_tfMRI_EMOTION_PA.nii.gz +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_PA/${Subject}_3T_tfMRI_EMOTION_PA_SBRef.nii.gz +# +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_AP/${Subject}_3T_tfMRI_EMOTION_AP.nii.gz +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_AP/${Subject}_3T_tfMRI_EMOTION_AP_SBRef.nii.gz +# +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_PA/${Subject}_3T_SpinEchoFieldMap_PA.nii.gz +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_PA/${Subject}_3T_SpinEchoFieldMap_AP.nii.gz +# +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_AP/${Subject}_3T_SpinEchoFieldMap_PA.nii.gz +# ${StudyFolder}/${Subject}/unprocessed/3T/tfMRI_EMOTION_AP/${Subject}_3T_SpinEchoFieldMap_AP.nii.gz +# +# +# Change Scan Settings: Dwelltime, FieldMap Delta TE (if using), and $PhaseEncodinglist to match your images +# These are set to match the HCP Protocol by default +# +# If using gradient distortion correction, use the coefficents from your scanner +# The HCP gradient distortion coefficents are only available through Siemens +# Gradient distortion in standard scanners like the Trio is much less than for the HCP Skyra. +# +# To get accurate EPI distortion correction with TOPUP, the flags in PhaseEncodinglist must match the phase encoding +# direction of the EPI scan, and you must have used the correct images in SpinEchoPhaseEncodeNegative and Positive +# variables. If the distortion is twice as bad as in the original images, flip either the order of the spin echo +# images or reverse the phase encoding list flag. The pipeline expects you to have used the same phase encoding +# axis in the fMRI data as in the spin echo field map data (x/-x or y/-y). ######################################### DO WORK ########################################## @@ -108,15 +125,25 @@ for Subject in $Subjlist ; do fMRITimeSeries="${StudyFolder}/${Subject}/unprocessed/3T/${fMRIName}/${Subject}_3T_${fMRIName}.nii.gz" fMRISBRef="${StudyFolder}/${Subject}/unprocessed/3T/${fMRIName}/${Subject}_3T_${fMRIName}_SBRef.nii.gz" #A single band reference image (SBRef) is recommended if using multiband, set to NONE if you want to use the first volume of the timeseries for motion correction DwellTime="0.00058" #Echo Spacing or Dwelltime of fMRI image, set to NONE if not used. Dwelltime = 1/(BandwidthPerPixelPhaseEncode * # of phase encoding samples): DICOM field (0019,1028) = BandwidthPerPixelPhaseEncode, DICOM field (0051,100b) AcquisitionMatrixText first value (# of phase encoding samples). On Siemens, iPAT/GRAPPA factors have already been accounted for. - DistortionCorrection="TOPUP" #FIELDMAP or TOPUP, distortion correction is required for accurate processing + DistortionCorrection="TOPUP" # FIELDMAP, SiemensFieldMap, GeneralElectricFieldMap, or TOPUP: distortion correction is required for accurate processing SpinEchoPhaseEncodeNegative="${StudyFolder}/${Subject}/unprocessed/3T/${fMRIName}/${Subject}_3T_SpinEchoFieldMap_LR.nii.gz" #For the spin echo field map volume with a negative phase encoding direction (LR in HCP data, AP in 7T HCP data), set to NONE if using regular FIELDMAP SpinEchoPhaseEncodePositive="${StudyFolder}/${Subject}/unprocessed/3T/${fMRIName}/${Subject}_3T_SpinEchoFieldMap_RL.nii.gz" #For the spin echo field map volume with a positive phase encoding direction (RL in HCP data, PA in 7T HCP data), set to NONE if using regular FIELDMAP MagnitudeInputName="NONE" #Expects 4D Magnitude volume with two 3D timepoints, set to NONE if using TOPUP PhaseInputName="NONE" #Expects a 3D Phase volume, set to NONE if using TOPUP + + # Path to General Electric style B0 fieldmap with two volumes + # 1. field map in degrees + # 2. magnitude + # Set to "NONE" if not using "GeneralElectricFieldMap" as the value for the DistortionCorrection variable + # + # Example Value: + # GEB0InputName="${StudyFolder}/${Subject}/unprocessed/3T/${fMRIName}/${Subject}_3T_GradientEchoFieldMap.nii.gz" + GEB0InputName="NONE" + DeltaTE="NONE" #2.46ms for 3T, 1.02ms for 7T, set to NONE if using TOPUP FinalFMRIResolution="2" #Target final resolution of fMRI data. 2mm is recommended for 3T HCP data, 1.6mm for 7T HCP data (i.e. should match acquired resolution). Use 2.0 or 1.0 to avoid standard FSL templates # GradientDistortionCoeffs="${HCPPIPEDIR_Config}/coeff_SC72C_Skyra.grad" #Gradient distortion correction coefficents, set to NONE to turn off - GradientDistortionCoeffs="NONE" # SEt to NONE to skip gradient distortion correction + GradientDistortionCoeffs="NONE" # Set to NONE to skip gradient distortion correction TopUpConfig="${HCPPIPEDIR_Config}/b02b0.cnf" #Topup config if using TOPUP, set to NONE if using regular FIELDMAP if [ -n "${command_line_specified_run_local}" ] ; then @@ -137,6 +164,7 @@ for Subject in $Subjlist ; do --SEPhasePos=$SpinEchoPhaseEncodePositive \ --fmapmag=$MagnitudeInputName \ --fmapphase=$PhaseInputName \ + --fmapgeneralelectric=$GEB0InputName \ --echospacing=$DwellTime \ --echodiff=$DeltaTE \ --unwarpdir=$UnwarpDir \ @@ -157,6 +185,7 @@ for Subject in $Subjlist ; do --SEPhasePos=$SpinEchoPhaseEncodePositive \ --fmapmag=$MagnitudeInputName \ --fmapphase=$PhaseInputName \ + --fmapgeneralelectric=$GEB0InputName \ --echospacing=$DwellTime \ --echodiff=$DeltaTE \ --unwarpdir=$UnwarpDir \ diff --git a/Examples/Scripts/PreFreeSurferPipelineBatch.sh b/Examples/Scripts/PreFreeSurferPipelineBatch.sh index fb7ab7502..e94ebabe8 100755 --- a/Examples/Scripts/PreFreeSurferPipelineBatch.sh +++ b/Examples/Scripts/PreFreeSurferPipelineBatch.sh @@ -49,7 +49,7 @@ fi # installed versions of: FSL (version 5.0.6), FreeSurfer (version 5.3.0-HCP), gradunwarp (HCP version 1.0.2) if doing gradient distortion correction # environment: FSLDIR , FREESURFER_HOME , HCPPIPEDIR , CARET7DIR , PATH (for gradient_unwarp.py) -#Set up pipeline environment variables and software +# Set up pipeline environment variables and software . ${EnvironmentScript} # Log the originating call @@ -66,29 +66,40 @@ PRINTCOM="" ########################################## INPUTS ########################################## -#Scripts called by this script do NOT assume anything about the form of the input names or paths. -#This batch script assumes the HCP raw data naming convention, e.g. - +# Scripts called by this script do NOT assume anything about the form of the input names or paths. +# This batch script assumes the HCP raw data naming convention, e.g. +# # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_T1w_MPR1.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR2/${Subject}_3T_T1w_MPR2.nii.gz - +# # ${StudyFolder}/${Subject}/unprocessed/3T/T2w_SPC1/${Subject}_3T_T2w_SPC1.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/T2w_SPC2/${Subject}_3T_T2w_SPC2.nii.gz - +# # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Magnitude.nii.gz # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Phase.nii.gz -#Change Scan Settings: Sample Spacings, and $UnwarpDir to match your images -#These are set to match the HCP Protocol by default - -#You have the option of using either gradient echo field maps or spin echo field maps to -#correct your structural images for readout distortion, or not to do this correction at all -#Change either the gradient echo field map or spin echo field map scan settings to match your data -#The default is to use gradient echo field maps using the HCP Protocol - -#If using gradient distortion correction, use the coefficents from your scanner -#The HCP gradient distortion coefficents are only available through Siemens -#Gradient distortion in standard scanners like the Trio is much less than for the HCP Skyra. +# Change Scan Settings: Sample Spacings, and $UnwarpDir to match your images +# These are set to match the HCP Protocol by default +# +# Readout Distortion Correction: +# +# You have the option of using either gradient echo field maps or spin echo field maps to +# perform readout distorction correction on your structural images, or not to do readout +# distortion correction at all. +# +# The HCP Pipeline Scripts currently support the use of gradient echo field maps or spin echo +# field maps as they are produced by the Siemens Connectom Scanner. They also support the +# use of gradient echo field maps as generated by General Electric scanners. +# +# Change either the gradient echo field map or spin echo field map scan settings to match +# your data. This script is setup to use gradient echo field maps from the Siemens Connectom +# Scanner using the HCP Protocol. +# +# Gradient Distortion Correction: +# +# If using gradient distortion correction, use the coefficents from your scanner. +# The HCP gradient distortion coefficents are only available through Siemens +# Gradient distortion in standard scanners like the Trio is much less than for the HCP Skyra. ######################################### DO WORK ########################################## @@ -97,8 +108,8 @@ PRINTCOM="" for Subject in $Subjlist ; do echo $Subject - #Input Images - #Detect Number of T1w Images + # Input Images + # Detect Number of T1w Images numT1ws=`ls ${StudyFolder}/${Subject}/unprocessed/3T | grep T1w_MPR | wc -l` echo "Found ${numT1ws} T1w Images for subject ${Subject}" T1wInputImages="" @@ -108,7 +119,7 @@ for Subject in $Subjlist ; do i=$(($i+1)) done - #Detect Number of T2w Images + # Detect Number of T2w Images numT2ws=`ls ${StudyFolder}/${Subject}/unprocessed/3T | grep T2w_SPC | wc -l` echo "Found ${numT2ws} T2w Images for subject ${Subject}" T2wInputImages="" @@ -118,22 +129,132 @@ for Subject in $Subjlist ; do i=$(($i+1)) done - #Readout Distortion Correction: - AvgrdcSTRING="FIELDMAP" #Averaging and readout distortion correction methods: "NONE" = average any repeats with no readout correction "FIELDMAP" = average any repeats and use field map for readout correction "TOPUP" = use spin echo field map + # Readout Distortion Correction: + # + # Currently supported Averaging and readout distortion correction methods: + # (i.e. supported values for the AvgrdcSTRING variable in this script and the + # --avgrdcmethod= command line option for the PreFreeSurferPipeline.sh script.) + # + # "NONE" + # Average any repeats but do no readout distortion correction + # + # "FIELDMAP" + # This value is equivalent to the "SiemensFieldMap" value described below. + # Use of the "SiemensFieldMap" value is prefered, but "FIELDMAP" is + # included for backward compatibility with the versions of these + # scripts that only supported use of Siemens-specific Gradient Echo + # Field Maps and did not support Gradient Echo Field Maps from any + # other scanner vendor. + # + # "TOPUP" + # Average any repeats and use Spin Echo Field Maps for readout distortion + # correction + # + # "GeneralElectricFieldMap" + # Average any repeats and use General Electric specific Gradient Echo + # Field Map for readout distortion correction + # + # "SiemensFieldMap" + # Average any repeats and use Siemens specific Gradient Echo Field Maps + # for readout distortion correction + + # + # Current Setup is for Siemens specific Gradient Echo Field Maps + # + # The following settings for AvgrdcSTRING, MagnitudeInputName, PhaseInputName, + # and TE are for using the Siemens specific Gradient Echo Field Maps that are + # collected and used in the standard HCP protocol. + # + # Note: The AvgrdcSTRING variable could also be set to the value "FIELDMAP" + # which is equivalent to "SiemensFieldMap". + AvgrdcSTRING="SiemensFieldMap" + + # ------------------------------------------------------------------------ + # Variables related to using Siemens specific Gradient Echo Field Maps + # ------------------------------------------------------------------------ + + # The MagnitudeInputName variable should be set to a 4D magitude volume + # with two 3D timepoints or "NONE" if not used + MagnitudeInputName="${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Magnitude.nii.gz" + + # The PhaseInputName variable should be set to a 3D phase difference volume + # or "NONE" if not used + PhaseInputName="${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Phase.nii.gz" + + # The TE variable should be set to 2.46ms for 3T scanner, 1.02ms for 7T + # scanner or "NONE" if not using + TE="2.46" + + # --------------------------------------------------- + # Variables related to using Spin Echo Field Maps + # --------------------------------------------------- + + # The following variables would be set to values other than "NONE" for + # using Spin Echo Field Maps (i.e. when AvgrdcSTRING="TOPUP") + + # The SpinEchoPhaseEncodeNegative variable should be set to the + # spin echo field map volume with a negative phase encoding direction + # (LR in 3T HCP data, AP in 7T HCP data), and set to "NONE" if not + # using Spin Echo Field Maps (i.e. if AvgrdcSTRING is not equal to "TOPUP") + # + # Example values for when using Spin Echo Field Maps: + # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_SpinEchoFieldMap_LR.nii.gz + # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_SpinEchoFieldMap_AP.nii.gz + SpinEchoPhaseEncodeNegative="NONE" + + # The SpinEchoPhaseEncodePositive variable should be set to the + # spin echo field map volume with positive phase encoding direction + # (RL in 3T HCP data, PA in 7T HCP data), and set to "NONE" if not + # using Spin Echo Field Maps (i.e. if AvgrdcSTRING is not equal to "TOPUP") + # + # Example values for when using Spin Echo Field Maps: + # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_SpinEchoFieldMap_RL.nii.gz + # ${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_SpinEchoFieldMap_PA.nii.gz + SpinEchoPhaseEncodePositive="NONE" + + # Echo Spacing or Dwelltime of spin echo EPI MRI image. Specified in seconds. + # Set to "NONE" if not used. + # + # Dwelltime = 1/(BandwidthPerPixelPhaseEncode * # of phase encoding samples) + # DICOM field (0019,1028) = BandwidthPerPixelPhaseEncode + # DICOM field (0051,100b) = AcquisitionMatrixText first value (# of phase encoding samples). + # On Siemens, iPAT/GRAPPA factors have already been accounted for. + # + # Example value for when using Spin Echo Field Maps: + # 0.000580002668012 + DwellTime="NONE" + + # Spin Echo Unwarping Direction + # x or y (minus or not does not matter) + # "NONE" if not used + # + # Example values for when using Spin Echo Field Maps: x, -x, y, -y + # Note: +x or +y are not supported. For positive values, do not include the + sign + SEUnwarpDir="NONE" + + # Topup Configuration file + # "NONE" if not used + TopupConfig="NONE" + + # --------------------------------------------------------------------------------- + # Variables related to using General Electric specific Gradient Echo Field Maps + # --------------------------------------------------------------------------------- + + # The following variables would be set to values other than "NONE" for + # using General Electric specific Gradient Echo Field Maps (i.e. when + # AvgrdcSTRING="GeneralElectricFieldMap") - #Using Regular Gradient Echo Field Maps (same as for fMRIVolume pipeline) - MagnitudeInputName="${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Magnitude.nii.gz" #Expects 4D magitude volume with two 3D timepoints or "NONE" if not used - PhaseInputName="${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_FieldMap_Phase.nii.gz" #Expects 3D phase difference volume or "NONE" if not used - TE="2.46" #2.46ms for 3T, 1.02ms for 7T, set to NONE if not using - - #Using Spin Echo Field Maps (same as for fMRIVolume pipeline) - SpinEchoPhaseEncodeNegative="NONE" #For the spin echo field map volume with a negative phase encoding direction (LR in HCP data, AP in 7T HCP data), set to NONE if using regular FIELDMAP - SpinEchoPhaseEncodePositive="NONE" #For the spin echo field map volume with a positive phase encoding direction (RL in HCP data, PA in 7T HCP data), set to NONE if using regular FIELDMAP - DwellTime="NONE" #Echo Spacing or Dwelltime of spin echo EPI MRI image, set to NONE if not used. Dwelltime = 1/(BandwidthPerPixelPhaseEncode * # of phase encoding samples): DICOM field (0019,1028) = BandwidthPerPixelPhaseEncode, DICOM field (0051,100b) AcquisitionMatrixText first value (# of phase encoding samples). On Siemens, iPAT/GRAPPA factors have already been accounted for. - SEUnwarpDir="NONE" #x or y (minus or not does not matter) "NONE" if not used - TopupConfig="NONE" #Config for topup or "NONE" if not used - - #Templates + # Example value for when using General Electric Gradient Echo Field Map + # + # GEB0InputName should be a General Electric style B0 fielmap with two volumes + # 1) fieldmap in deg and + # 2) magnitude, + # set to NONE if using TOPUP or FIELDMAP/SiemensFieldMap + # + # GEB0InputName="${StudyFolder}/${Subject}/unprocessed/3T/T1w_MPR1/${Subject}_3T_GradientEchoFieldMap.nii.gz" + GEB0InputName="NONE" + + # Templates T1wTemplate="${HCPPIPEDIR_Templates}/MNI152_T1_0.7mm.nii.gz" #Hires T1w MNI template T1wTemplateBrain="${HCPPIPEDIR_Templates}/MNI152_T1_0.7mm_brain.nii.gz" #Hires brain extracted MNI template T1wTemplate2mm="${HCPPIPEDIR_Templates}/MNI152_T1_2mm.nii.gz" #Lowres T1w MNI template @@ -143,14 +264,22 @@ for Subject in $Subjlist ; do TemplateMask="${HCPPIPEDIR_Templates}/MNI152_T1_0.7mm_brain_mask.nii.gz" #Hires MNI brain mask template Template2mmMask="${HCPPIPEDIR_Templates}/MNI152_T1_2mm_brain_mask_dil.nii.gz" #Lowres MNI brain mask template - #Structural Scan Settings (set all to NONE if not doing readout distortion correction) + # Structural Scan Settings (set all to NONE if not doing readout distortion correction) + # + # Sample values for when using General Electric Gradient Echo Field Maps + # T1wSampleSpacing="0.000011999" # For General Electric scanners, 1/((0018,0095)*(0028,0010)) + # T2wSampleSpacing="0.000008000" # For General Electric scanners, 1/((0018,0095)*(0028,0010)) + # UnwarpDir="y" + # + # The values set below are for the HCP Protocol using the Siemens Connectom Scanner T1wSampleSpacing="0.0000074" #DICOM field (0019,1018) in s or "NONE" if not used T2wSampleSpacing="0.0000021" #DICOM field (0019,1018) in s or "NONE" if not used - UnwarpDir="z" #z appears to be best or "NONE" if not used + UnwarpDir="z" # z appears to be best for Siemens Gradient Echo Field Maps or "NONE" if not used - #Other Config Settings + # Other Config Settings BrainSize="150" #BrainSize in mm, 150 for humans FNIRTConfig="${HCPPIPEDIR_Config}/T1_2_MNI152_2mm.cnf" #FNIRT 2mm T1w Config + # GradientDistortionCoeffs="${HCPPIPEDIR_Config}/coeff_SC72C_Skyra.grad" #Location of Coeffs file or "NONE" to skip GradientDistortionCoeffs="NONE" # Set to NONE to skip gradient distortion correction @@ -179,6 +308,7 @@ for Subject in $Subjlist ; do --fnirtconfig="$FNIRTConfig" \ --fmapmag="$MagnitudeInputName" \ --fmapphase="$PhaseInputName" \ + --fmapgeneralelectric="$GEB0InputName" \ --echodiff="$TE" \ --SEPhaseNeg="$SpinEchoPhaseEncodeNegative" \ --SEPhasePos="$SpinEchoPhaseEncodePositive" \ @@ -210,11 +340,12 @@ for Subject in $Subjlist ; do --fnirtconfig=${FNIRTConfig} \ --fmapmag=${MagnitudeInputName} \ --fmapphase=${PhaseInputName} \ + --fmapgeneralelectric=${GEB0InputName} \ --echodiff=${TE} \ --SEPhaseNeg=${SpinEchoPhaseEncodeNegative} \ --SEPhasePos=${SpinEchoPhaseEncodePositive} \ --echospacing=${DwellTime} \ - --seunwarpdir=${SEUnwarpDir} \ + --seunwarpdir=${SEUnwarpDir} \ --t1samplespacing=${T1wSampleSpacing} \ --t2samplespacing=${T2wSampleSpacing} \ --unwarpdir=${UnwarpDir} \ diff --git a/PreFreeSurfer/PreFreeSurferPipeline.sh b/PreFreeSurfer/PreFreeSurferPipeline.sh index 00a05e291..f18c8d587 100755 --- a/PreFreeSurfer/PreFreeSurferPipeline.sh +++ b/PreFreeSurfer/PreFreeSurferPipeline.sh @@ -17,7 +17,9 @@ # * Matthew F. Glasser, Department of Anatomy and Neurobiology, Washington University in St. Louis # * Mark Jenkinson, FMRIB Centre, Oxford University # * Timothy B. Brown, Neuroinformatics Research Group, Washington University in St. Louis -# +# * Modifications to support General Electric Gradient Echo field maps for readout distortion correction +# are based on example code provided by Gaurav Patel, Columbia University +# # ## Product # # [Human Connectome Project][HCP] (HCP) Pipelines @@ -151,6 +153,16 @@ # script itself exits and does not attempt any further processing. set -e +# ----------------------------------------------------------------------------------- +# Constants for specification of Averaging and Readout Distortion Correction Method +# ----------------------------------------------------------------------------------- + +NONE_METHOD_OPT="NONE" +FIELDMAP_METHOD_OPT="FIELDMAP" +SIEMENS_METHOD_OPT="SiemensFieldMap" +SPIN_ECHO_METHOD_OPT="TOPUP" +GENERAL_ELECTRIC_METHOD_OPT="GeneralElectricFieldMap" + # ------------------------------------------------------------------------------ # Load Function Libraries # ------------------------------------------------------------------------------ @@ -179,48 +191,67 @@ Usage: PreeFreeSurferPipeline.sh [options] (T1w) structural images for the subject (required) --t2= An @ symbol separated list of full paths to T2-weighted (T2w) structural images for the subject (required) - --t1template= MNI T1w template - --t1templatebrain= Brain extracted MNI T1wTemplate - --t1template2mm= MNI 2mm T1wTemplate - --t2template= MNI T2w template - --t2templatebrain= Brain extracted MNI T2wTemplate - --t2template2mm= MNI 2mm T2wTemplate - --templatemask= Brain mask MNI Template - --template2mmmask= Brain mask MNI 2mm Template - --brainsize= Brain size estimate in mm, 150 for humans - --fnirtconfig= FNIRT 2mm T1w Configuration file - --fmapmag= Fieldmap magnitude file - --fmapphase= Fieldmap phase file - --echodiff= Delta TE in ms for field map or "NONE" if - not used - --SEPhaseNeg={LR, RL, NONE} For spin echo field map volume with a negative - phase encoding direction (LR in HCP data), set - to "NONE" if using regular FIELDMAP - --SEPhasePos={LR, RL, NONE} For the spin echo field map volume with a - positive phase encoding direction (RL in HCP - data), set to "NONE" if using regular FIELDMAP - --echospacing= Echo Spacing or Dwelltime of Spin Echo Field - Map or "NONE" if not used - --seunwarpdir={x, y, NONE} Phase encoding direction of the spin echo - field map. (Only applies when using a spin echo - field map.) - --t1samplespacing= T1 image sample spacing, "NONE" if not used - --t2samplespacing= T2 image sample spacing, "NONE" if not used - --unwarpdir={x, y, z} Readout direction of the T1w and T2w images - (Used with either a regular field map or a spin - echo field map) - --gdcoeffs= File containing gradient distortion - coefficients, Set to "NONE" to turn off - --avgrdcmethod={NONE, FIELDMAP, TOPUP} - Averaging and readout distortion correction - method. - "NONE" = average any repeats with no readout - correction - "FIELDMAP" = average any repeats and use field - map for readout correction - "TOPUP" = average any repeats and use spin - echo field map for readout - correction + --t1template= MNI T1w template + --t1templatebrain= Brain extracted MNI T1wTemplate + --t1template2mm= MNI 2mm T1wTemplate + --t2template= MNI T2w template + --t2templatebrain= Brain extracted MNI T2wTemplate + --t2template2mm= MNI 2mm T2wTemplate + --templatemask= Brain mask MNI Template + --template2mmmask= Brain mask MNI 2mm Template + --brainsize= Brain size estimate in mm, 150 for humans + --fnirtconfig= FNIRT 2mm T1w Configuration file + --fmapmag= Siemens Gradient Echo Fieldmap magnitude file + --fmapphase= Siemens Gradient Echo Fieldmap phase file + --fmapgeneralelectric= General Electric Gradient Echo Field Map file + Two volumes in one file + 1. field map in deg + 2. magnitude + --echodiff= Delta TE in ms for field map or "NONE" if + not used + --SEPhaseNeg={, NONE} For spin echo field map, path to volume with + a negative phase encoding direction (LR in + HCP data), set to "NONE" if not using Spin + Echo Field Maps + --SEPhasePos={, NONE} For spin echo field map, path to volume with + a positive phase encoding direction (RL in + HCP data), set to "NONE" if not using Spin + Echo Field Maps + --echospacing= Echo Spacing or Dwelltime of Spin Echo Field + Map or "NONE" if not used + --seunwarpdir={x, y, NONE} Phase encoding direction of the spin echo + field map. (Only applies when using a spin echo + field map.) + --t1samplespacing= T1 image sample spacing, "NONE" if not used + --t2samplespacing= T2 image sample spacing, "NONE" if not used + --unwarpdir={x, y, z} Readout direction of the T1w and T2w images + (Used with either a gradient echo field map + or a spin echo field map) + --gdcoeffs= File containing gradient distortion + coefficients, Set to "NONE" to turn off + --avgrdcmethod= Averaging and readout distortion correction + method. See below for supported values. + + "${NONE_METHOD_OPT}" + average any repeats with no readout distortion correction + + "${FIELDMAP_METHOD_OPT}" + equivalent to "${SIEMENS_METHOD_OPT}" (see below) + SiemensFieldMap is preferred. This option value is maintained for + backward compatibility. + + "${SPIN_ECHO_METHOD_OPT}" + average any repeats and use Spin Echo Field Maps for readout + distortion correction + + "${GENERAL_ELECTRIC_METHOD_OPT}" + average any repeats and use General Electric specific Gradient + Echo Field Maps for readout distortion correction + + "${SIEMENS_METHOD_OPT}" + average any repeats and use Siemens specific Gradient Echo + Field Maps for readout distortion correction + --topupconfig= Configuration file for topup or "NONE" if not used --bfsigma= Bias Field Smoothing Sigma (optional) @@ -264,6 +295,7 @@ BrainSize=`opts_GetOpt1 "--brainsize" $@` FNIRTConfig=`opts_GetOpt1 "--fnirtconfig" $@` MagnitudeInputName=`opts_GetOpt1 "--fmapmag" $@` PhaseInputName=`opts_GetOpt1 "--fmapphase" $@` +GEB0InputName=`opts_GetOpt1 "--fmapgeneralelectric" $@` TE=`opts_GetOpt1 "--echodiff" $@` SpinEchoPhaseEncodeNegative=`opts_GetOpt1 "--SEPhaseNeg" $@` SpinEchoPhaseEncodePositive=`opts_GetOpt1 "--SEPhasePos" $@` @@ -302,6 +334,7 @@ log_Msg "BrainSize: ${BrainSize}" log_Msg "FNIRTConfig: ${FNIRTConfig}" log_Msg "MagnitudeInputName: ${MagnitudeInputName}" log_Msg "PhaseInputName: ${PhaseInputName}" +log_Msg "GEB0InputName: ${GEB0InputName}" log_Msg "TE: ${TE}" log_Msg "SpinEchoPhaseEncodeNegative: ${SpinEchoPhaseEncodeNegative}" log_Msg "SpinEchoPhaseEncodePositive: ${SpinEchoPhaseEncodePositive}" @@ -360,7 +393,7 @@ log_Msg "POSIXLY_CORRECT="${POSIXLY_CORRECT} # Loop over the processing for T1w and T2w (just with different names). # For each modality, perform # - Gradient Nonlinearity Correction (Unless no gradient distortion -# coefficients are available +# coefficients are available) # - Average same modality images (if more than one is available) # - Rigidly align images to 0.7mm MNI Template to create native volume space # - Perform Brain Extraction(FNIRT-based Masking) @@ -428,9 +461,8 @@ for TXw in ${Modalities} ; do log_Msg "Averaging ${TXw} Images" log_Msg "mkdir -p ${TXwFolder}/Average${TXw}Images" mkdir -p ${TXwFolder}/Average${TXw}Images - log_Msg "PERFORMING SIMPLE AVERAGING" + log_Msg "PERFORMING SIMPLE AVERAGING" ${RUN} ${HCPPIPEDIR_PreFS}/AnatomicalAverage.sh -o ${TXwFolder}/${TXwImage} -s ${TXwTemplate} -m ${TemplateMask} -n -w ${TXwFolder}/Average${TXw}Images --noclean -v -b $BrainSize $OutputTXwImageSTRING - # fi else log_Msg "Not Averaging ${TXw} Images" log_Msg "ONLY ONE AVERAGE FOUND: COPYING" @@ -472,68 +504,77 @@ done # T2w to T1w Registration and Optional Readout Distortion Correction # ------------------------------------------------------------------------------ -if [[ ${AvgrdcSTRING} = "FIELDMAP" || ${AvgrdcSTRING} = "TOPUP" ]] ; then - log_Msg "Performing ${AvgrdcSTRING} Readout Distortion Correction" - wdir=${T2wFolder}/T2wToT1wDistortionCorrectAndReg - if [ -d ${wdir} ] ; then - # DO NOT change the following line to "rm -r ${wdir}" because the - # chances of something going wrong with that are much higher, and - # rm -r always needs to be treated with the utmost caution - rm -r ${T2wFolder}/T2wToT1wDistortionCorrectAndReg - fi - - log_Msg "mkdir -p ${wdir}" - mkdir -p ${wdir} - - ${RUN} ${HCPPIPEDIR_PreFS}/T2wToT1wDistortionCorrectAndReg.sh \ - --workingdir=${wdir} \ - --t1=${T1wFolder}/${T1wImage}_acpc \ - --t1brain=${T1wFolder}/${T1wImage}_acpc_brain \ - --t2=${T2wFolder}/${T2wImage}_acpc \ - --t2brain=${T2wFolder}/${T2wImage}_acpc_brain \ - --fmapmag=${MagnitudeInputName} \ - --fmapphase=${PhaseInputName} \ - --echodiff=${TE} \ - --SEPhaseNeg=${SpinEchoPhaseEncodeNegative} \ - --SEPhasePos=${SpinEchoPhaseEncodePositive} \ - --echospacing=${DwellTime} \ - --seunwarpdir=${SEUnwarpDir} \ - --t1sampspacing=${T1wSampleSpacing} \ - --t2sampspacing=${T2wSampleSpacing} \ - --unwarpdir=${UnwarpDir} \ - --ot1=${T1wFolder}/${T1wImage}_acpc_dc \ - --ot1brain=${T1wFolder}/${T1wImage}_acpc_dc_brain \ - --ot1warp=${T1wFolder}/xfms/${T1wImage}_dc \ - --ot2=${T1wFolder}/${T2wImage}_acpc_dc \ - --ot2warp=${T1wFolder}/xfms/${T2wImage}_reg_dc \ - --method=${AvgrdcSTRING} \ - --topupconfig=${TopupConfig} \ - --gdcoeffs=${GradientDistortionCoeffs} -else - log_Msg "NOT PERFORMING READOUT DISTORTION CORRECTION" - wdir=${T2wFolder}/T2wToT1wReg - if [ -e ${wdir} ] ; then - # DO NOT change the following line to "rm -r ${wdir}" because the - # chances of something going wrong with that are much higher, and - # rm -r always needs to be treated with the utmost caution - rm -r ${T2wFolder}/T2wToT1wReg - fi - - log_Msg "mkdir -p ${wdir}" - mkdir -p ${wdir} - - ${RUN} ${HCPPIPEDIR_PreFS}/T2wToT1wReg.sh \ - ${wdir} \ - ${T1wFolder}/${T1wImage}_acpc \ - ${T1wFolder}/${T1wImage}_acpc_brain \ - ${T2wFolder}/${T2wImage}_acpc \ - ${T2wFolder}/${T2wImage}_acpc_brain \ - ${T1wFolder}/${T1wImage}_acpc_dc \ - ${T1wFolder}/${T1wImage}_acpc_dc_brain \ - ${T1wFolder}/xfms/${T1wImage}_dc \ - ${T1wFolder}/${T2wImage}_acpc_dc \ - ${T1wFolder}/xfms/${T2wImage}_reg_dc -fi +case $AvgrdcSTRING in + + ${FIELDMAP_METHOD_OPT} | ${SPIN_ECHO_METHOD_OPT} | ${GENERAL_ELECTRIC_METHOD_OPT} | ${SIEMENS_METHOD_OPT}) + + log_Msg "Performing ${AvgrdcSTRING} Readout Distortion Correction" + wdir=${T2wFolder}/T2wToT1wDistortionCorrectAndReg + if [ -d ${wdir} ] ; then + # DO NOT change the following line to "rm -r ${wdir}" because the + # chances of something going wrong with that are much higher, and + # rm -r always needs to be treated with the utmost caution + rm -r ${T2wFolder}/T2wToT1wDistortionCorrectAndReg + fi + + log_Msg "mkdir -p ${wdir}" + mkdir -p ${wdir} + + ${RUN} ${HCPPIPEDIR_PreFS}/T2wToT1wDistortionCorrectAndReg.sh \ + --workingdir=${wdir} \ + --t1=${T1wFolder}/${T1wImage}_acpc \ + --t1brain=${T1wFolder}/${T1wImage}_acpc_brain \ + --t2=${T2wFolder}/${T2wImage}_acpc \ + --t2brain=${T2wFolder}/${T2wImage}_acpc_brain \ + --fmapmag=${MagnitudeInputName} \ + --fmapphase=${PhaseInputName} \ + --fmapgeneralelectric=${GEB0InputName} \ + --echodiff=${TE} \ + --SEPhaseNeg=${SpinEchoPhaseEncodeNegative} \ + --SEPhasePos=${SpinEchoPhaseEncodePositive} \ + --echospacing=${DwellTime} \ + --seunwarpdir=${SEUnwarpDir} \ + --t1sampspacing=${T1wSampleSpacing} \ + --t2sampspacing=${T2wSampleSpacing} \ + --unwarpdir=${UnwarpDir} \ + --ot1=${T1wFolder}/${T1wImage}_acpc_dc \ + --ot1brain=${T1wFolder}/${T1wImage}_acpc_dc_brain \ + --ot1warp=${T1wFolder}/xfms/${T1wImage}_dc \ + --ot2=${T1wFolder}/${T2wImage}_acpc_dc \ + --ot2warp=${T1wFolder}/xfms/${T2wImage}_reg_dc \ + --method=${AvgrdcSTRING} \ + --topupconfig=${TopupConfig} \ + --gdcoeffs=${GradientDistortionCoeffs} + + ;; + + *) + + log_Msg "NOT PERFORMING READOUT DISTORTION CORRECTION" + wdir=${T2wFolder}/T2wToT1wReg + if [ -e ${wdir} ] ; then + # DO NOT change the following line to "rm -r ${wdir}" because the + # chances of something going wrong with that are much higher, and + # rm -r always needs to be treated with the utmost caution + rm -r ${T2wFolder}/T2wToT1wReg + fi + + log_Msg "mkdir -p ${wdir}" + mkdir -p ${wdir} + + ${RUN} ${HCPPIPEDIR_PreFS}/T2wToT1wReg.sh \ + ${wdir} \ + ${T1wFolder}/${T1wImage}_acpc \ + ${T1wFolder}/${T1wImage}_acpc_brain \ + ${T2wFolder}/${T2wImage}_acpc \ + ${T2wFolder}/${T2wImage}_acpc_brain \ + ${T1wFolder}/${T1wImage}_acpc_dc \ + ${T1wFolder}/${T1wImage}_acpc_dc_brain \ + ${T1wFolder}/xfms/${T1wImage}_dc \ + ${T1wFolder}/${T2wImage}_acpc_dc \ + ${T1wFolder}/xfms/${T2wImage}_reg_dc + +esac # ------------------------------------------------------------------------------ # Bias Field Correction: Calculate bias field using square root of the product diff --git a/PreFreeSurfer/scripts/T2wToT1wDistortionCorrectAndReg.sh b/PreFreeSurfer/scripts/T2wToT1wDistortionCorrectAndReg.sh index 58ae4a3c7..0341a7c34 100755 --- a/PreFreeSurfer/scripts/T2wToT1wDistortionCorrectAndReg.sh +++ b/PreFreeSurfer/scripts/T2wToT1wDistortionCorrectAndReg.sh @@ -1,9 +1,20 @@ -#!/bin/bash +!/bin/bash set -e # Requirements for this script # installed versions of: FSL (version 5.0.6), HCP-gradunwarp (HCP version 1.0.2) # environment: FSLDIR and PATH for gradient_unwarp.py +SCRIPT_NAME="T2WToT1wDistortionCorrectAndReg.sh" + +# ----------------------------------------------------------------------------------- +# Constants for specification of Averaging and Readout Distortion Correction Method +# ----------------------------------------------------------------------------------- + +SIEMENS_METHOD_OPT="SiemensFieldMap" +SPIN_ECHO_METHOD_OPT="TOPUP" +GENERAL_ELECTRIC_METHOD_OPT="GeneralElectricFieldMap" +FIELDMAP_METHOD_OPT="FIELDMAP" + ################################################ SUPPORT FUNCTIONS ################################################## Usage() { @@ -16,6 +27,7 @@ Usage() { echo " --t2brain=" echo " [--fmapmag=]" echo " [--fmapphase=]" + echo " [--fmapgeneralelectric=]" echo " [--echodiff=]" echo " [--SEPhaseNeg=]" echo " [--SEPhasePos=]" @@ -30,7 +42,21 @@ Usage() { echo " --ot2=" echo " --ot2brain=" echo " --ot2warp=" - echo " --method=" + echo " --method=" + echo "" + echo " ${FIELDMAP_METHOD_OPT}" + echo " equivalent to ${SIEMENS_METHOD_OPT} (see below)" + echo " ${SIEMENS_METHOD_OPT} is preferred. This option is maintained for" + echo " backward compatibility." + echo " ${SPIN_ECHO_METHOD_OPT}" + echo " use Spin Echo Field Maps for readout distortion correction" + echo " ${GENERAL_ELECTRIC_METHOD_OPT}" + echo " use General Electric specific Gradient Echo Field Maps for" + echo " readout distortion correction" + echo " ${SIEMENS_METHOD_OPT}" + echo " use Siemens specific Gradient Echo Field Maps for readout" + echo " distortion correction" + echo "" echo " [--topupconfig=]" echo " [--gdcoeffs=]" } @@ -82,29 +108,30 @@ if [ $# -eq 0 ] ; then Usage; exit 0; fi if [ $# -lt 17 ] ; then Usage; exit 1; fi # parse arguments -WD=`getopt1 "--workingdir" $@` # "$1" -T1wImage=`getopt1 "--t1" $@` # "$2" -T1wImageBrain=`getopt1 "--t1brain" $@` # "$3" -T2wImage=`getopt1 "--t2" $@` # "$4" -T2wImageBrain=`getopt1 "--t2brain" $@` # "$5" -MagnitudeInputName=`getopt1 "--fmapmag" $@` # "$6" -PhaseInputName=`getopt1 "--fmapphase" $@` # "$7" -TE=`getopt1 "--echodiff" $@` # "$8" -SpinEchoPhaseEncodeNegative=`getopt1 "--SEPhaseNeg" $@` # "$7" -SpinEchoPhaseEncodePositive=`getopt1 "--SEPhasePos" $@` # "$5" -DwellTime=`getopt1 "--echospacing" $@` # "$9" -SEUnwarpDir=`getopt1 "--seunwarpdir" $@` # "${11}" -T1wSampleSpacing=`getopt1 "--t1sampspacing" $@` # "$9" -T2wSampleSpacing=`getopt1 "--t2sampspacing" $@` # "${10}" -UnwarpDir=`getopt1 "--unwarpdir" $@` # "${11}" -OutputT1wImage=`getopt1 "--ot1" $@` # "${12}" -OutputT1wImageBrain=`getopt1 "--ot1brain" $@` # "${13}" -OutputT1wTransform=`getopt1 "--ot1warp" $@` # "${14}" -OutputT2wImage=`getopt1 "--ot2" $@` # "${15}" -OutputT2wTransform=`getopt1 "--ot2warp" $@` # "${16}" -DistortionCorrection=`getopt1 "--method" $@` # "${21}" -TopupConfig=`getopt1 "--topupconfig" $@` # "${22}" -GradientDistortionCoeffs=`getopt1 "--gdcoeffs" $@` # "${18}" +WD=`getopt1 "--workingdir" $@` +T1wImage=`getopt1 "--t1" $@` +T1wImageBrain=`getopt1 "--t1brain" $@` +T2wImage=`getopt1 "--t2" $@` +T2wImageBrain=`getopt1 "--t2brain" $@` +MagnitudeInputName=`getopt1 "--fmapmag" $@` +PhaseInputName=`getopt1 "--fmapphase" $@` +GEB0InputName=`getopt1 "--fmapgeneralelectric" $@` +TE=`getopt1 "--echodiff" $@` +SpinEchoPhaseEncodeNegative=`getopt1 "--SEPhaseNeg" $@` +SpinEchoPhaseEncodePositive=`getopt1 "--SEPhasePos" $@` +DwellTime=`getopt1 "--echospacing" $@` +SEUnwarpDir=`getopt1 "--seunwarpdir" $@` +T1wSampleSpacing=`getopt1 "--t1sampspacing" $@` +T2wSampleSpacing=`getopt1 "--t2sampspacing" $@` +UnwarpDir=`getopt1 "--unwarpdir" $@` +OutputT1wImage=`getopt1 "--ot1" $@` +OutputT1wImageBrain=`getopt1 "--ot1brain" $@` +OutputT1wTransform=`getopt1 "--ot1warp" $@` +OutputT2wImage=`getopt1 "--ot2" $@` +OutputT2wTransform=`getopt1 "--ot2warp" $@` +DistortionCorrection=`getopt1 "--method" $@` +TopupConfig=`getopt1 "--topupconfig" $@` +GradientDistortionCoeffs=`getopt1 "--gdcoeffs" $@` # default parameters WD=`defaultopt $WD .` @@ -122,7 +149,7 @@ T2wImageBasename=`basename "$T2wImage"` Modalities="T1w T2w" echo " " -echo " START: T2wToT1wDistortionCorrectionAndReg" +echo " START: ${SCRIPT_NAME}" mkdir -p $WD mkdir -p ${WD}/FieldMap @@ -136,46 +163,87 @@ echo " " >> $WD/log.txt ########################################## DO WORK ########################################## -###### FIELDMAP VERSION (GE FIELDMAPS) ###### -if [ $DistortionCorrection = "FIELDMAP" ] ; then - ### Create fieldmaps (and apply gradient non-linearity distortion correction) - echo " " - echo " " - echo " " - #echo ${HCPPIPEDIR_Global}/FieldMapPreprocessingAll.sh ${WD}/FieldMap ${MagnitudeInputName} ${PhaseInputName} ${TE} ${WD}/Magnitude ${WD}/Magnitude_brain ${WD}/Phase ${WD}/FieldMap ${GradientDistortionCoeffs} ${GlobalScripts} - - ${HCPPIPEDIR_Global}/FieldMapPreprocessingAll.sh \ - --workingdir=${WD}/FieldMap \ - --fmapmag=${MagnitudeInputName} \ - --fmapphase=${PhaseInputName} \ - --echodiff=${TE} \ - --ofmapmag=${WD}/Magnitude \ - --ofmapmagbrain=${WD}/Magnitude_brain \ - --ofmap=${WD}/FieldMap \ - --gdcoeffs=${GradientDistortionCoeffs} -###### TOPUP VERSION (SE FIELDMAPS) ###### -elif [ $DistortionCorrection = "TOPUP" ] ; then - if [[ ${SEUnwarpDir} = "x" || ${SEUnwarpDir} = "y" ]] ; then - ScoutInputName="${SpinEchoPhaseEncodePositive}" - elif [[ ${SEUnwarpDir} = "-x" || ${SEUnwarpDir} = "-y" || ${SEUnwarpDir} = "x-" || ${SEUnwarpDir} = "y-" ]] ; then - ScoutInputName="${SpinEchoPhaseEncodeNegative}" - fi - # Use topup to distortion correct the scout scans - # using a blip-reversed SE pair "fieldmap" sequence - ${HCPPIPEDIR_Global}/TopupPreprocessingAll.sh \ - --workingdir=${WD}/FieldMap \ - --phaseone=${SpinEchoPhaseEncodeNegative} \ - --phasetwo=${SpinEchoPhaseEncodePositive} \ - --scoutin=${ScoutInputName} \ - --echospacing=${DwellTime} \ - --unwarpdir=${SEUnwarpDir} \ - --ofmapmag=${WD}/Magnitude \ - --ofmapmagbrain=${WD}/Magnitude_brain \ - --ofmap=${WD}/FieldMap \ - --ojacobian=${WD}/Jacobian \ - --gdcoeffs=${GradientDistortionCoeffs} \ - --topupconfig=${TopupConfig} -fi +case $DistortionCorrection in + + ${FIELDMAP_METHOD_OPT} | ${SIEMENS_METHOD_OPT}) + + # -------------------------------------- + # -- Siemens Gradient Echo Field Maps -- + # -------------------------------------- + + ### Create fieldmaps (and apply gradient non-linearity distortion correction) + echo " " + echo " " + echo " " + + ${HCPPIPEDIR_Global}/SiemensFieldMapPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --fmapmag=${MagnitudeInputName} \ + --fmapphase=${PhaseInputName} \ + --echodiff=${TE} \ + --ofmapmag=${WD}/Magnitude \ + --ofmapmagbrain=${WD}/Magnitude_brain \ + --ofmap=${WD}/FieldMap \ + --gdcoeffs=${GradientDistortionCoeffs} + + ;; + + ${GENERAL_ELECTRIC_METHOD_OPT}) + + # ----------------------------------------------- + # -- General Electric Gradient Echo Field Maps -- + # ----------------------------------------------- + + ### Create fieldmaps (and apply gradient non-linearity distortion correction) + echo " " + echo " " + echo " " + + ${HCPPIPEDIR_Global}/GeneralElectricFieldMapPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --fmap=${GEB0InputName} \ + --ofmapmag=${WD}/Magnitude \ + --ofmapmagbrain=${WD}/Magnitude_brain \ + --ofmap=${WD}/FieldMap \ + --gdcoeffs=${GradientDistortionCoeffs} + + ;; + + ${SPIN_ECHO_METHOD_OPT}) + + # -------------------------- + # -- Spin Echo Field Maps -- + # -------------------------- + + if [[ ${SEUnwarpDir} = "x" || ${SEUnwarpDir} = "y" ]] ; then + ScoutInputName="${SpinEchoPhaseEncodePositive}" + elif [[ ${SEUnwarpDir} = "-x" || ${SEUnwarpDir} = "-y" || ${SEUnwarpDir} = "x-" || ${SEUnwarpDir} = "y-" ]] ; then + ScoutInputName="${SpinEchoPhaseEncodeNegative}" + fi + + # Use topup to distortion correct the scout scans + # using a blip-reversed SE pair "fieldmap" sequence + ${HCPPIPEDIR_Global}/TopupPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --phaseone=${SpinEchoPhaseEncodeNegative} \ + --phasetwo=${SpinEchoPhaseEncodePositive} \ + --scoutin=${ScoutInputName} \ + --echospacing=${DwellTime} \ + --unwarpdir=${SEUnwarpDir} \ + --ofmapmag=${WD}/Magnitude \ + --ofmapmagbrain=${WD}/Magnitude_brain \ + --ofmap=${WD}/FieldMap \ + --ojacobian=${WD}/Jacobian \ + --gdcoeffs=${GradientDistortionCoeffs} \ + --topupconfig=${TopupConfig} + + ;; + + *) + echo "${SCRIPT_NAME} - ERROR - Unable to create FSL-suitable readout distortion correction field map" + echo "${SCRIPT_NAME} Unrecognized distortion correction method: ${DistortionCorrection}" + exit 1 +esac ### LOOP over available modalities ### @@ -199,13 +267,23 @@ for TXw in $Modalities ; do ${FSLDIR}/bin/fugue --loadfmap=${WD}/FieldMap --dwell=${TXwSampleSpacing} --saveshift=${WD}/FieldMap_ShiftMap${TXw}.nii.gz ${FSLDIR}/bin/convertwarp --relout --rel --ref=${WD}/Magnitude --shiftmap=${WD}/FieldMap_ShiftMap${TXw}.nii.gz --shiftdir=${UnwarpDir} --out=${WD}/FieldMap_Warp${TXw}.nii.gz - if [ $DistortionCorrection = "FIELDMAP" ] ; then - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Magnitude -r ${WD}/Magnitude -w ${WD}/FieldMap_Warp${TXw}.nii.gz -o ${WD}/Magnitude_warpped${TXw} - ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_warpped${TXw} -ref ${TXwImage} -out ${WD}/Magnitude_warpped${TXw}2${TXwImageBasename} -omat ${WD}/Fieldmap2${TXwImageBasename}.mat -searchrx -30 30 -searchry -30 30 -searchrz -30 30 - elif [ $DistortionCorrection = "TOPUP" ] ; then - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Magnitude_brain -r ${WD}/Magnitude_brain -w ${WD}/FieldMap_Warp${TXw}.nii.gz -o ${WD}/Magnitude_brain_warpped${TXw} - ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_brain_warpped${TXw} -ref ${TXwImageBrain} -out ${WD}/Magnitude_brain_warpped${TXw}2${TXwImageBasename} -omat ${WD}/Fieldmap2${TXwImageBasename}.mat -searchrx -30 30 -searchry -30 30 -searchrz -30 30 - fi + case $DistortionCorrection in + + ${FIELDMAP_METHOD_OPT} | ${SIEMENS_METHOD_OPT} | ${GENERAL_ELECTRIC_METHOD_OPT}) + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Magnitude -r ${WD}/Magnitude -w ${WD}/FieldMap_Warp${TXw}.nii.gz -o ${WD}/Magnitude_warpped${TXw} + ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_warpped${TXw} -ref ${TXwImage} -out ${WD}/Magnitude_warpped${TXw}2${TXwImageBasename} -omat ${WD}/Fieldmap2${TXwImageBasename}.mat -searchrx -30 30 -searchry -30 30 -searchrz -30 30 + ;; + + ${SPIN_ECHO_METHOD_OPT}) + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Magnitude_brain -r ${WD}/Magnitude_brain -w ${WD}/FieldMap_Warp${TXw}.nii.gz -o ${WD}/Magnitude_brain_warpped${TXw} + ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_brain_warpped${TXw} -ref ${TXwImageBrain} -out ${WD}/Magnitude_brain_warpped${TXw}2${TXwImageBasename} -omat ${WD}/Fieldmap2${TXwImageBasename}.mat -searchrx -30 30 -searchry -30 30 -searchrz -30 30 + ;; + + *) + echo "${SCRIPT_NAME} - ERROR - Unable to apply readout distortion correction" + echo "${SCRIPT_NAME} Unrecognized distortion correction method: ${DistortionCorrection}" + exit 1 + esac ${FSLDIR}/bin/flirt -in ${WD}/FieldMap.nii.gz -ref ${TXwImage} -applyxfm -init ${WD}/Fieldmap2${TXwImageBasename}.mat -out ${WD}/FieldMap2${TXwImageBasename} @@ -252,7 +330,7 @@ ${FSLDIR}/bin/imcp ${WD}/T2w2T1w/T2w_dc_reg ${OutputT2wTransform} ${FSLDIR}/bin/imcp ${WD}/T2w2T1w/T2w_reg ${OutputT2wImage} echo " " -echo " END: T2wToT1wDistortionCorrectionAndReg" +echo " END: ${SCRIPT_NAME}" echo " END: `date`" >> $WD/log.txt ########################################## QA STUFF ########################################## diff --git a/fMRIVolume/GenericfMRIVolumeProcessingPipeline.sh b/fMRIVolume/GenericfMRIVolumeProcessingPipeline.sh index c6d24498e..40d126258 100755 --- a/fMRIVolume/GenericfMRIVolumeProcessingPipeline.sh +++ b/fMRIVolume/GenericfMRIVolumeProcessingPipeline.sh @@ -47,22 +47,27 @@ fi log_Msg "Parsing Command Line Options" # parse arguments -Path=`opts_GetOpt1 "--path" $@` # "$1" -Subject=`opts_GetOpt1 "--subject" $@` # "$2" -NameOffMRI=`opts_GetOpt1 "--fmriname" $@` # "$6" -fMRITimeSeries=`opts_GetOpt1 "--fmritcs" $@` # "$3" -fMRIScout=`opts_GetOpt1 "--fmriscout" $@` # "$4" -SpinEchoPhaseEncodeNegative=`opts_GetOpt1 "--SEPhaseNeg" $@` # "$7" -SpinEchoPhaseEncodePositive=`opts_GetOpt1 "--SEPhasePos" $@` # "$5" -MagnitudeInputName=`opts_GetOpt1 "--fmapmag" $@` # "$8" #Expects 4D volume with two 3D timepoints -PhaseInputName=`opts_GetOpt1 "--fmapphase" $@` # "$9" -DwellTime=`opts_GetOpt1 "--echospacing" $@` # "${11}" -deltaTE=`opts_GetOpt1 "--echodiff" $@` # "${12}" -UnwarpDir=`opts_GetOpt1 "--unwarpdir" $@` # "${13}" -FinalfMRIResolution=`opts_GetOpt1 "--fmrires" $@` # "${14}" -DistortionCorrection=`opts_GetOpt1 "--dcmethod" $@` # "${17}" #FIELDMAP or TOPUP -GradientDistortionCoeffs=`opts_GetOpt1 "--gdcoeffs" $@` # "${18}" -TopupConfig=`opts_GetOpt1 "--topupconfig" $@` # "${20}" #NONE if Topup is not being used +Path=`opts_GetOpt1 "--path" $@` +Subject=`opts_GetOpt1 "--subject" $@` +NameOffMRI=`opts_GetOpt1 "--fmriname" $@` +fMRITimeSeries=`opts_GetOpt1 "--fmritcs" $@` +fMRIScout=`opts_GetOpt1 "--fmriscout" $@` +SpinEchoPhaseEncodeNegative=`opts_GetOpt1 "--SEPhaseNeg" $@` +SpinEchoPhaseEncodePositive=`opts_GetOpt1 "--SEPhasePos" $@` +MagnitudeInputName=`opts_GetOpt1 "--fmapmag" $@` # Expects 4D volume with two 3D timepoints +PhaseInputName=`opts_GetOpt1 "--fmapphase" $@` +GEB0InputName=`opts_GetOpt1 "--fmapgeneralelectric" $@` +DwellTime=`opts_GetOpt1 "--echospacing" $@` +deltaTE=`opts_GetOpt1 "--echodiff" $@` +UnwarpDir=`opts_GetOpt1 "--unwarpdir" $@` +FinalfMRIResolution=`opts_GetOpt1 "--fmrires" $@` + +# FIELDMAP, SiemensFieldMap, GeneralElectricFieldMap, or TOPUP +# Note: FIELDMAP and SiemensFieldMap are equivalent +DistortionCorrection=`opts_GetOpt1 "--dcmethod" $@` + +GradientDistortionCoeffs=`opts_GetOpt1 "--gdcoeffs" $@` +TopupConfig=`opts_GetOpt1 "--topupconfig" $@` # NONE if Topup is not being used RUN=`opts_GetOpt1 "--printcom" $@` # use ="echo" for just printing everything and not running the commands (default is to run) # Setup PATHS @@ -157,7 +162,7 @@ ${RUN} "$PipelineScripts"/MotionCorrection_FLIRTbased.sh \ "$fMRIFolder"/"$MotionMatrixFolder" \ "$MotionMatrixPrefix" -#EPI Distortion Correction and EPI to T1w Registration +# EPI Distortion Correction and EPI to T1w Registration log_Msg "EPI Distortion Correction and EPI to T1w Registration" if [ -e ${fMRIFolder}/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased ] ; then rm -r ${fMRIFolder}/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased @@ -173,6 +178,7 @@ ${RUN} ${PipelineScripts}/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurf --t1brain=${T1wFolder}/${T1wRestoreImageBrain} \ --fmapmag=${MagnitudeInputName} \ --fmapphase=${PhaseInputName} \ + --fmapgeneralelectric=${GEB0InputName} \ --echodiff=${deltaTE} \ --SEPhaseNeg=${SpinEchoPhaseEncodeNegative} \ --SEPhasePos=${SpinEchoPhaseEncodePositive} \ diff --git a/fMRIVolume/scripts/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh b/fMRIVolume/scripts/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh index 6664f8a87..b95c92d50 100755 --- a/fMRIVolume/scripts/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh +++ b/fMRIVolume/scripts/DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh @@ -5,6 +5,17 @@ set -e # installed versions of: FSL (version 5.0.6) and FreeSurfer (version 5.3.0-HCP) # environment: FSLDIR, FREESURFER_HOME + others +SCRIPT_NAME="DistortionCorrectionAndEPIToT1wReg_FLIRTBBRAndFreeSurferBBRbased.sh" + +# --------------------------------------------------------------------- +# Constants for specification of Readout Distortion Correction Method +# --------------------------------------------------------------------- + +FIELDMAP_METHOD_OPT="FIELDMAP" +SIEMENS_METHOD_OPT="SiemensFieldMap" +GENERAL_ELECTRIC_METHOD_OPT="GeneralElectricFieldMap" +SPIN_ECHO_METHOD_OPT="TOPUP" + ################################################ SUPPORT FUNCTIONS ################################################## Usage() { @@ -15,8 +26,9 @@ Usage() { echo " --t1=" echo " --t1restore=" echo " --t1brain=" - echo " --fmapmag=" - echo " --fmapphase=" + echo " --fmapmag=" + echo " --fmapphase=" + echo " --fmapgeneralelectric=" echo " --echodiff=" echo " --SEPhaseNeg=" echo " --SEPhasePos=" @@ -29,7 +41,23 @@ Usage() { echo " --freesurfersubjectid=" echo " --gdcoeffs=" echo " [--qaimage=]" - echo " --method=" + echo "" + echo " --method=" + echo "" + echo " \"${FIELDMAP_METHOD_OPT}\"" + echo " equivalent to ${SIEMENS_METHOD_OPT} (see below)" + echo "" + echo " \"${SIEMENS_METHOD_OPT}\"" + echo " use Siemens specific Gradient Echo Field Maps for" + echo " readout distortion correction" + echo "" + echo " \"${SPIN_ECHO_METHOD_OPT}\"" + echo " use Spin Echo Field Maps for readout distortion correction" + echo "" + echo " \"${GENERAL_ELECTRIC_METHOD_OPT}\"" + echo " use General Electric specific Gradient Echo Field Maps" + echo " for readout distortion correction" + echo "" echo " [--topupconfig=]" echo " --ojacobian=" @@ -55,7 +83,7 @@ defaultopt() { # Outputs (in $WD): # -# FIELDMAP section only: +# FIELDMAP, SiemensFieldMap, and GeneralElectricFieldMap: # Magnitude Magnitude_brain FieldMap # # FIELDMAP and TOPUP sections: @@ -72,38 +100,38 @@ defaultopt() { # # ${RegOutput} ${OutputTransform} ${JacobianOut} ${QAImage} - - ################################################## OPTION PARSING ##################################################### + # Just give usage if no arguments specified if [ $# -eq 0 ] ; then Usage; exit 0; fi # check for correct options if [ $# -lt 21 ] ; then Usage; exit 1; fi # parse arguments -WD=`getopt1 "--workingdir" $@` # "$1" -ScoutInputName=`getopt1 "--scoutin" $@` # "$2" -T1wImage=`getopt1 "--t1" $@` # "$3" -T1wRestoreImage=`getopt1 "--t1restore" $@` # "$4" -T1wBrainImage=`getopt1 "--t1brain" $@` # "$5" -SpinEchoPhaseEncodeNegative=`getopt1 "--SEPhaseNeg" $@` # "$7" -SpinEchoPhaseEncodePositive=`getopt1 "--SEPhasePos" $@` # "$5" -DwellTime=`getopt1 "--echospacing" $@` # "$9" -MagnitudeInputName=`getopt1 "--fmapmag" $@` # "$6" -PhaseInputName=`getopt1 "--fmapphase" $@` # "$7" -deltaTE=`getopt1 "--echodiff" $@` # "$8" -UnwarpDir=`getopt1 "--unwarpdir" $@` # "${10}" -OutputTransform=`getopt1 "--owarp" $@` # "${11}" -BiasField=`getopt1 "--biasfield" $@` # "${12}" -RegOutput=`getopt1 "--oregim" $@` # "${13}" -FreeSurferSubjectFolder=`getopt1 "--freesurferfolder" $@` # "${14}" -FreeSurferSubjectID=`getopt1 "--freesurfersubjectid" $@` # "${15}" -GradientDistortionCoeffs=`getopt1 "--gdcoeffs" $@` # "${17}" -QAImage=`getopt1 "--qaimage" $@` # "${20}" -DistortionCorrection=`getopt1 "--method" $@` # "${21}" -TopupConfig=`getopt1 "--topupconfig" $@` # "${22}" -JacobianOut=`getopt1 "--ojacobian" $@` # "${23}" +WD=`getopt1 "--workingdir" $@` +ScoutInputName=`getopt1 "--scoutin" $@` +T1wImage=`getopt1 "--t1" $@` +T1wRestoreImage=`getopt1 "--t1restore" $@` +T1wBrainImage=`getopt1 "--t1brain" $@` +SpinEchoPhaseEncodeNegative=`getopt1 "--SEPhaseNeg" $@` +SpinEchoPhaseEncodePositive=`getopt1 "--SEPhasePos" $@` +DwellTime=`getopt1 "--echospacing" $@` +MagnitudeInputName=`getopt1 "--fmapmag" $@` +PhaseInputName=`getopt1 "--fmapphase" $@` +GEB0InputName=`getopt1 "--fmapgeneralelectric" $@` +deltaTE=`getopt1 "--echodiff" $@` +UnwarpDir=`getopt1 "--unwarpdir" $@` +OutputTransform=`getopt1 "--owarp" $@` +BiasField=`getopt1 "--biasfield" $@` +RegOutput=`getopt1 "--oregim" $@` +FreeSurferSubjectFolder=`getopt1 "--freesurferfolder" $@` +FreeSurferSubjectID=`getopt1 "--freesurfersubjectid" $@` +GradientDistortionCoeffs=`getopt1 "--gdcoeffs" $@` +QAImage=`getopt1 "--qaimage" $@` +DistortionCorrection=`getopt1 "--method" $@` +TopupConfig=`getopt1 "--topupconfig" $@` +JacobianOut=`getopt1 "--ojacobian" $@` ScoutInputFile=`basename $ScoutInputName` T1wBrainImageFile=`basename $T1wBrainImage` @@ -117,7 +145,7 @@ TopupConfig=`defaultopt $TopupConfig ${HCPPIPEDIR_Config}/b02b0.cnf` UseJacobian=false echo " " -echo " START: DistortionCorrectionEpiToT1wReg_FLIRTBBRAndFreeSurferBBRBased" +echo " START: ${SCRIPT_NAME}" mkdir -p $WD @@ -135,106 +163,149 @@ fi cp ${T1wBrainImage}.nii.gz ${WD}/${T1wBrainImageFile}.nii.gz -###### FIELDMAP VERSION (GE FIELDMAPS) ###### -if [ $DistortionCorrection = "FIELDMAP" ] ; then - # process fieldmap with gradient non-linearity distortion correction - ${GlobalScripts}/FieldMapPreprocessingAll.sh \ - --workingdir=${WD}/FieldMap \ - --fmapmag=${MagnitudeInputName} \ - --fmapphase=${PhaseInputName} \ - --echodiff=${deltaTE} \ - --ofmapmag=${WD}/Magnitude \ - --ofmapmagbrain=${WD}/Magnitude_brain \ - --ofmap=${WD}/FieldMap \ - --gdcoeffs=${GradientDistortionCoeffs} - cp ${ScoutInputName}.nii.gz ${WD}/Scout.nii.gz - #Test if Magnitude Brain and T1w Brain Are Similar in Size, if not, assume Magnitude Brain Extraction Failed and Must Be Retried After Removing Bias Field - MagnitudeBrainSize=`${FSLDIR}/bin/fslstats ${WD}/Magnitude_brain -V | cut -d " " -f 2` - T1wBrainSize=`${FSLDIR}/bin/fslstats ${WD}/${T1wBrainImageFile} -V | cut -d " " -f 2` - - if [[ X`echo "if ( (${MagnitudeBrainSize} / ${T1wBrainSize}) > 1.25 ) {1}" | bc -l` = X1 || X`echo "if ( (${MagnitudeBrainSize} / ${T1wBrainSize}) < 0.75 ) {1}" | bc -l` = X1 ]] ; then - ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude.nii.gz -ref ${T1wImage} -omat "$WD"/Mag2T1w.mat -out ${WD}/Magnitude2T1w.nii.gz -searchrx -30 30 -searchry -30 30 -searchrz -30 30 - ${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Mag.mat -inverse "$WD"/Mag2T1w.mat - #${FSLDIR}/bin/applywarp --interp=spline -i ${BiasField} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_Bias.nii.gz - #mv ${WD}/Magnitude.nii.gz ${WD}/MagnitudeOrig.nii.gz - #fslmaths ${WD}/MagnitudeOrig.nii.gz -div ${WD}/Magnitude_Bias.nii.gz ${WD}/Magnitude.nii.gz - ${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_brain_mask.nii.gz - #fslmaths ${WD}/Magnitude_brain_mask.nii.gz -bin -dilD ${WD}/Magnitude_brain_mask.nii.gz - #fslmaths ${WD}/Magnitude.nii.gz -mas ${WD}/Magnitude_brain_mask.nii.gz ${WD}/Magnitude_brain.nii.gz - #fslmaths ${WD}/${T1wBrainImageFile} -bin ${WD}/${T1wBrainImageFile}_mask - #fslmaths ${T1wImage} -mas ${WD}/${T1wBrainImageFile}_mask ${WD}/${T1wBrainImageFile}_norestore_brain - #${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_brain.nii.gz -ref ${WD}/${T1wBrainImageFile}_norestore_brain -init "$WD"/Mag2T1w.mat -omat "$WD"/Mag2T1w.mat -out ${WD}/Magnitude2T1w.nii.gz -nosearch - #${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Mag.mat -inverse "$WD"/Mag2T1w.mat - #${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_brain_mask.nii.gz - #${FSLDIR}/bin/bet ${WD}/Magnitude.nii.gz ${WD}/Magnitude_brain.nii.gz -f 0.35 -m #Brain extract the magnitude image - fslmaths ${WD}/Magnitude_brain_mask.nii.gz -bin ${WD}/Magnitude_brain_mask.nii.gz - fslmaths ${WD}/Magnitude.nii.gz -mas ${WD}/Magnitude_brain_mask.nii.gz ${WD}/Magnitude_brain.nii.gz - - ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Scout.nii.gz -ref ${T1wImage} -omat "$WD"/Scout2T1w.mat -out ${WD}/Scout2T1w.nii.gz -searchrx -30 30 -searchry -30 30 -searchrz -30 30 - ${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Scout.mat -inverse "$WD"/Scout2T1w.mat - #${FSLDIR}/bin/applywarp --interp=spline -i ${BiasField} -r ${WD}/Scout.nii.gz --premat="$WD"/T1w2Scout.mat -o ${WD}/Scout_Bias.nii.gz - #mv ${WD}/Scout.nii.gz ${WD}/ScoutOrig.nii.gz - #fslmaths ${WD}/ScoutOrig.nii.gz -div ${WD}/Scout_Bias.nii.gz ${WD}/Scout.nii.gz - ${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Scout.nii.gz --premat="$WD"/T1w2Scout.mat -o ${WD}/Scout_brain_mask.nii.gz - #${FSLDIR}/bin/bet ${WD}/Scout.nii.gz ${WD}/Scout_brain.nii.gz -f 0.35 -m #Brain extract the magnitude image - fslmaths ${WD}/Scout_brain_mask.nii.gz -bin ${WD}/Scout_brain_mask.nii.gz - fslmaths ${WD}/Scout.nii.gz -mas ${WD}/Scout_brain_mask.nii.gz ${WD}/Scout_brain.nii.gz +case $DistortionCorrection in + + ${FIELDMAP_METHOD_OPT} | ${SIEMENS_METHOD_OPT} | ${GENERAL_ELECTRIC_METHOD_OPT}) + + if [ $DistortionCorrection = "${FIELDMAP_METHOD_OPT}" ] || [ $DistortionCorrection = "${SIEMENS_METHOD_OPT}" ] ; then + # -------------------------------------- + # -- Siemens Gradient Echo Field Maps -- + # -------------------------------------- + + # process fieldmap with gradient non-linearity distortion correction + ${GlobalScripts}/SiemensFieldMapPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --fmapmag=${MagnitudeInputName} \ + --fmapphase=${PhaseInputName} \ + --echodiff=${deltaTE} \ + --ofmapmag=${WD}/Magnitude \ + --ofmapmagbrain=${WD}/Magnitude_brain \ + --ofmap=${WD}/FieldMap \ + --gdcoeffs=${GradientDistortionCoeffs} + + elif [ $DistortionCorrection = "${GENERAL_ELECTRIC_METHOD_OPT}" ] ; then + # ----------------------------------------------- + # -- General Electric Gradient Echo Field Maps -- + # ----------------------------------------------- + + # process fieldmap with gradient non-linearity distortion correction + ${GlobalScripts}/GeneralElectricFieldMapPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --fmap=${GEB0InputName} \ + --ofmapmag=${WD}/Magnitude \ + --ofmapmagbrain=${WD}/Magnitude_brain \ + --ofmap=${WD}/FieldMap \ + --gdcoeffs=${GradientDistortionCoeffs} + + else + echo "${SCRIPT_NAME}: Script programming error. Unhandled Distortion Correction Method: ${DistortionCorrection}" + exit 1 + fi + + cp ${ScoutInputName}.nii.gz ${WD}/Scout.nii.gz + + # Test if Magnitude Brain and T1w Brain Are Similar in Size, if not, assume Magnitude Brain Extraction + # Failed and Must Be Retried After Removing Bias Field + MagnitudeBrainSize=`${FSLDIR}/bin/fslstats ${WD}/Magnitude_brain -V | cut -d " " -f 2` + T1wBrainSize=`${FSLDIR}/bin/fslstats ${WD}/${T1wBrainImageFile} -V | cut -d " " -f 2` + + if [[ X`echo "if ( (${MagnitudeBrainSize} / ${T1wBrainSize}) > 1.25 ) {1}" | bc -l` = X1 || X`echo "if ( (${MagnitudeBrainSize} / ${T1wBrainSize}) < 0.75 ) {1}" | bc -l` = X1 ]] ; then + ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude.nii.gz -ref ${T1wImage} -omat "$WD"/Mag2T1w.mat -out ${WD}/Magnitude2T1w.nii.gz -searchrx -30 30 -searchry -30 30 -searchrz -30 30 + ${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Mag.mat -inverse "$WD"/Mag2T1w.mat + #${FSLDIR}/bin/applywarp --interp=spline -i ${BiasField} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_Bias.nii.gz + #mv ${WD}/Magnitude.nii.gz ${WD}/MagnitudeOrig.nii.gz + #fslmaths ${WD}/MagnitudeOrig.nii.gz -div ${WD}/Magnitude_Bias.nii.gz ${WD}/Magnitude.nii.gz + ${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_brain_mask.nii.gz + #fslmaths ${WD}/Magnitude_brain_mask.nii.gz -bin -dilD ${WD}/Magnitude_brain_mask.nii.gz + #fslmaths ${WD}/Magnitude.nii.gz -mas ${WD}/Magnitude_brain_mask.nii.gz ${WD}/Magnitude_brain.nii.gz + #fslmaths ${WD}/${T1wBrainImageFile} -bin ${WD}/${T1wBrainImageFile}_mask + #fslmaths ${T1wImage} -mas ${WD}/${T1wBrainImageFile}_mask ${WD}/${T1wBrainImageFile}_norestore_brain + #${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Magnitude_brain.nii.gz -ref ${WD}/${T1wBrainImageFile}_norestore_brain -init "$WD"/Mag2T1w.mat -omat "$WD"/Mag2T1w.mat -out ${WD}/Magnitude2T1w.nii.gz -nosearch + #${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Mag.mat -inverse "$WD"/Mag2T1w.mat + #${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Magnitude.nii.gz --premat="$WD"/T1w2Mag.mat -o ${WD}/Magnitude_brain_mask.nii.gz + #${FSLDIR}/bin/bet ${WD}/Magnitude.nii.gz ${WD}/Magnitude_brain.nii.gz -f 0.35 -m #Brain extract the magnitude image + fslmaths ${WD}/Magnitude_brain_mask.nii.gz -bin ${WD}/Magnitude_brain_mask.nii.gz + fslmaths ${WD}/Magnitude.nii.gz -mas ${WD}/Magnitude_brain_mask.nii.gz ${WD}/Magnitude_brain.nii.gz + + ${FSLDIR}/bin/flirt -interp spline -dof 6 -in ${WD}/Scout.nii.gz -ref ${T1wImage} -omat "$WD"/Scout2T1w.mat -out ${WD}/Scout2T1w.nii.gz -searchrx -30 30 -searchry -30 30 -searchrz -30 30 + ${FSLDIR}/bin/convert_xfm -omat "$WD"/T1w2Scout.mat -inverse "$WD"/Scout2T1w.mat + #${FSLDIR}/bin/applywarp --interp=spline -i ${BiasField} -r ${WD}/Scout.nii.gz --premat="$WD"/T1w2Scout.mat -o ${WD}/Scout_Bias.nii.gz + #mv ${WD}/Scout.nii.gz ${WD}/ScoutOrig.nii.gz + #fslmaths ${WD}/ScoutOrig.nii.gz -div ${WD}/Scout_Bias.nii.gz ${WD}/Scout.nii.gz + ${FSLDIR}/bin/applywarp --interp=nn -i ${WD}/${T1wBrainImageFile} -r ${WD}/Scout.nii.gz --premat="$WD"/T1w2Scout.mat -o ${WD}/Scout_brain_mask.nii.gz + #${FSLDIR}/bin/bet ${WD}/Scout.nii.gz ${WD}/Scout_brain.nii.gz -f 0.35 -m #Brain extract the magnitude image + fslmaths ${WD}/Scout_brain_mask.nii.gz -bin ${WD}/Scout_brain_mask.nii.gz + fslmaths ${WD}/Scout.nii.gz -mas ${WD}/Scout_brain_mask.nii.gz ${WD}/Scout_brain.nii.gz - # register scout to T1w image using fieldmap - ${FSLDIR}/bin/epi_reg --epi=${WD}/Scout_brain.nii.gz --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted --fmap=${WD}/FieldMap.nii.gz --fmapmag=${WD}/Magnitude.nii.gz --fmapmagbrain=${WD}/Magnitude_brain.nii.gz --echospacing=${DwellTime} --pedir=${UnwarpDir} - else - # register scout to T1w image using fieldmap - ${FSLDIR}/bin/epi_reg --epi=${WD}/Scout.nii.gz --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted --fmap=${WD}/FieldMap.nii.gz --fmapmag=${WD}/Magnitude.nii.gz --fmapmagbrain=${WD}/Magnitude_brain.nii.gz --echospacing=${DwellTime} --pedir=${UnwarpDir} - fi - # convert epi_reg warpfield from abs to rel convention (NB: this is the current convention for epi_reg but it may change in the future, or take an option) - #${FSLDIR}/bin/immv ${WD}/${ScoutInputFile}_undistorted_warp ${WD}/${ScoutInputFile}_undistorted_warp_abs - #${FSLDIR}/bin/convertwarp --relout --abs -r ${WD}/${ScoutInputFile}_undistorted_warp_abs -w ${WD}/${ScoutInputFile}_undistorted_warp_abs -o ${WD}/${ScoutInputFile}_undistorted_warp - # create spline interpolated output for scout to T1w + apply bias field correction - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${T1wImage} -w ${WD}/${ScoutInputFile}_undistorted_warp.nii.gz -o ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz - ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz -div ${BiasField} ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz - ${FSLDIR}/bin/immv ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz - ###Jacobian Volume FAKED for Regular Fieldmaps (all ones) ### - ${FSLDIR}/bin/fslmaths ${T1wImage} -abs -add 1 -bin ${WD}/Jacobian2T1w.nii.gz + # register scout to T1w image using fieldmap + ${FSLDIR}/bin/epi_reg --epi=${WD}/Scout_brain.nii.gz --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted --fmap=${WD}/FieldMap.nii.gz --fmapmag=${WD}/Magnitude.nii.gz --fmapmagbrain=${WD}/Magnitude_brain.nii.gz --echospacing=${DwellTime} --pedir=${UnwarpDir} + else + # register scout to T1w image using fieldmap + ${FSLDIR}/bin/epi_reg --epi=${WD}/Scout.nii.gz --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted --fmap=${WD}/FieldMap.nii.gz --fmapmag=${WD}/Magnitude.nii.gz --fmapmagbrain=${WD}/Magnitude_brain.nii.gz --echospacing=${DwellTime} --pedir=${UnwarpDir} + fi + + # convert epi_reg warpfield from abs to rel convention (NB: this is the current convention for epi_reg but it may change in the future, or take an option) + #${FSLDIR}/bin/immv ${WD}/${ScoutInputFile}_undistorted_warp ${WD}/${ScoutInputFile}_undistorted_warp_abs + #${FSLDIR}/bin/convertwarp --relout --abs -r ${WD}/${ScoutInputFile}_undistorted_warp_abs -w ${WD}/${ScoutInputFile}_undistorted_warp_abs -o ${WD}/${ScoutInputFile}_undistorted_warp + # create spline interpolated output for scout to T1w + apply bias field correction + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${T1wImage} -w ${WD}/${ScoutInputFile}_undistorted_warp.nii.gz -o ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz + ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz -div ${BiasField} ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz + ${FSLDIR}/bin/immv ${WD}/${ScoutInputFile}_undistorted_1vol.nii.gz ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz + ###Jacobian Volume FAKED for Regular Fieldmaps (all ones) ### + ${FSLDIR}/bin/fslmaths ${T1wImage} -abs -add 1 -bin ${WD}/Jacobian2T1w.nii.gz -###### TOPUP VERSION (SE FIELDMAPS) ###### -elif [ $DistortionCorrection = "TOPUP" ] ; then - # Use topup to distortion correct the scout scans - # using a blip-reversed SE pair "fieldmap" sequence - ${GlobalScripts}/TopupPreprocessingAll.sh \ - --workingdir=${WD}/FieldMap \ - --phaseone=${SpinEchoPhaseEncodeNegative} \ - --phasetwo=${SpinEchoPhaseEncodePositive} \ - --scoutin=${ScoutInputName} \ - --echospacing=${DwellTime} \ - --unwarpdir=${UnwarpDir} \ - --owarp=${WD}/WarpField \ - --ojacobian=${WD}/Jacobian \ - --gdcoeffs=${GradientDistortionCoeffs} \ - --topupconfig=${TopupConfig} - - # create a spline interpolated image of scout (distortion corrected in same space) - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${ScoutInputName} -w ${WD}/WarpField.nii.gz -o ${WD}/${ScoutInputFile}_undistorted - # apply Jacobian correction to scout image (optional) - if [ $UseJacobian = true ] ; then - ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -mul ${WD}/Jacobian.nii.gz ${WD}/${ScoutInputFile}_undistorted - fi - # register undistorted scout image to T1w - ${FSLDIR}/bin/epi_reg --epi=${WD}/${ScoutInputFile}_undistorted --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted - # generate combined warpfields and spline interpolated images + apply bias field correction - ${FSLDIR}/bin/convertwarp --relout --rel -r ${T1wImage} --warp1=${WD}/WarpField.nii.gz --postmat=${WD}/${ScoutInputFile}_undistorted.mat -o ${WD}/${ScoutInputFile}_undistorted_warp - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Jacobian.nii.gz -r ${T1wImage} --premat=${WD}/${ScoutInputFile}_undistorted.mat -o ${WD}/Jacobian2T1w.nii.gz - ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${T1wImage} -w ${WD}/${ScoutInputFile}_undistorted_warp -o ${WD}/${ScoutInputFile}_undistorted - # apply Jacobian correction to scout image (optional) - if [ $UseJacobian = true ] ; then - ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -div ${BiasField} -mul ${WD}/Jacobian2T1w.nii.gz ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz - else - ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -div ${BiasField} ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz - fi -else - echo "UNKNOWN DISTORTION CORRECTION METHOD" - exit -fi + ;; + ${SPIN_ECHO_METHOD_OPT}) + + # -------------------------- + # -- Spin Echo Field Maps -- + # -------------------------- + + # Use topup to distortion correct the scout scans + # using a blip-reversed SE pair "fieldmap" sequence + ${GlobalScripts}/TopupPreprocessingAll.sh \ + --workingdir=${WD}/FieldMap \ + --phaseone=${SpinEchoPhaseEncodeNegative} \ + --phasetwo=${SpinEchoPhaseEncodePositive} \ + --scoutin=${ScoutInputName} \ + --echospacing=${DwellTime} \ + --unwarpdir=${UnwarpDir} \ + --owarp=${WD}/WarpField \ + --ojacobian=${WD}/Jacobian \ + --gdcoeffs=${GradientDistortionCoeffs} \ + --topupconfig=${TopupConfig} + + # create a spline interpolated image of scout (distortion corrected in same space) + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${ScoutInputName} -w ${WD}/WarpField.nii.gz -o ${WD}/${ScoutInputFile}_undistorted + + # apply Jacobian correction to scout image (optional) + if [ $UseJacobian = true ] ; then + ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -mul ${WD}/Jacobian.nii.gz ${WD}/${ScoutInputFile}_undistorted + fi + + # register undistorted scout image to T1w + ${FSLDIR}/bin/epi_reg --epi=${WD}/${ScoutInputFile}_undistorted --t1=${T1wImage} --t1brain=${WD}/${T1wBrainImageFile} --out=${WD}/${ScoutInputFile}_undistorted + # generate combined warpfields and spline interpolated images + apply bias field correction + ${FSLDIR}/bin/convertwarp --relout --rel -r ${T1wImage} --warp1=${WD}/WarpField.nii.gz --postmat=${WD}/${ScoutInputFile}_undistorted.mat -o ${WD}/${ScoutInputFile}_undistorted_warp + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/Jacobian.nii.gz -r ${T1wImage} --premat=${WD}/${ScoutInputFile}_undistorted.mat -o ${WD}/Jacobian2T1w.nii.gz + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${ScoutInputName} -r ${T1wImage} -w ${WD}/${ScoutInputFile}_undistorted_warp -o ${WD}/${ScoutInputFile}_undistorted + + # apply Jacobian correction to scout image (optional) + if [ $UseJacobian = true ] ; then + ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -div ${BiasField} -mul ${WD}/Jacobian2T1w.nii.gz ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz + else + ${FSLDIR}/bin/fslmaths ${WD}/${ScoutInputFile}_undistorted -div ${BiasField} ${WD}/${ScoutInputFile}_undistorted2T1w_init.nii.gz + fi + + ;; + + *) + + echo "${SCRIPT_NAME}: UNKNOWN DISTORTION CORRECTION METHOD: ${DistortionCorrection}" + exit 1 + +esac ### FREESURFER BBR - found to be an improvement, probably due to better GM/WM boundary SUBJECTS_DIR=${FreeSurferSubjectFolder} diff --git a/global/scripts/GeneralElectricFieldMapPreprocessingAll.sh b/global/scripts/GeneralElectricFieldMapPreprocessingAll.sh new file mode 100755 index 000000000..4792ee7b8 --- /dev/null +++ b/global/scripts/GeneralElectricFieldMapPreprocessingAll.sh @@ -0,0 +1,127 @@ +#!/bin/bash + +set -e + +# Requirements for this script +# installed versions of: FSL5.0.1 or higher, gradunwarp python package (from MGH) +# environment: as in SetUpHCPPipeline.sh (or individually: FSLDIR, HCPPIPEDIR_Global and PATH for gradient_unwarp.py) + +################################################ SUPPORT FUNCTIONS ################################################## + +Usage() { + echo "`basename $0`: Script for generating a fieldmap suitable for FSL from General Electric Gradient Echo field map," + echo " and also do gradient non-linearity distortion correction of these" + echo " " + echo "Usage: `basename $0` " + echo " [--workingdir=]" + echo " --fmap=" + echo " --ofmapmag=" + echo " --ofmapmagbrain=" + echo " --ofmap=" + echo " [--gdcoeffs=]" +} + +# function for parsing options +getopt1() { + sopt="$1" + shift 1 + for fn in $@ ; do + if [ `echo $fn | grep -- "^${sopt}=" | wc -w` -gt 0 ] ; then + echo $fn | sed "s/^${sopt}=//" + return 0 + fi + done +} + +defaultopt() { + echo $1 +} + + +################################################### OUTPUT FILES ##################################################### + +# Output images (in $WD): Magnitude Magnitude_brain Magnitude_brain_mask FieldMap +# Plus the following if gradient distortion correction is run: +# Magnitude_gdc Magnitude_gdc_warp Magnitude_brain_gdc FieldMap_gdc +# Output images (not in $WD): ${MagnitudeOutput} ${MagnitudeBrainOutput} ${FieldMapOutput} + +################################################## OPTION PARSING ##################################################### + +# Just give usage if no arguments specified +if [ $# -eq 0 ] ; then Usage; exit 0; fi +# check for correct options +if [ $# -lt 5 ] ; then Usage; exit 1; fi + +# parse arguments +WD=`getopt1 "--workingdir" $@` +GEB0InputName=`getopt1 "--fmap" $@` +MagnitudeOutput=`getopt1 "--ofmapmag" $@` +MagnitudeBrainOutput=`getopt1 "--ofmapmagbrain" $@` +FieldMapOutput=`getopt1 "--ofmap" $@` +GradientDistortionCoeffs=`getopt1 "--gdcoeffs" $@` + +# default parameters +GlobalScripts=${HCPPIPEDIR_Global} +WD=`defaultopt $WD .` +GradientDistortionCoeffs=`defaultopt $GradientDistortionCoeffs "NONE"` + +echo " " +echo " START: Field Map Preprocessing and Gradient Unwarping" + +mkdir -p $WD + +# Record the input options in a log file +echo "$0 $@" >> $WD/log.txt +echo "PWD = `pwd`" >> $WD/log.txt +echo "date: `date`" >> $WD/log.txt +echo " " >> $WD/log.txt + +########################################## DO WORK ########################################## + +${FSLDIR}/bin/fslsplit ${GEB0InputName} # split image into vol0000=fieldmap and vol0001=magnitude +mv vol0000.nii.gz ${WD}/FieldMap_deg.nii.gz +mv vol0001.nii.gz ${WD}/Magnitude.nii.gz +${FSLDIR}/bin/bet ${WD}/Magnitude ${WD}/Magnitude_brain -f 0.35 -m #Brain extract the magnitude image +${FSLDIR}/bin/fslmaths ${WD}/FieldMap_deg.nii.gz -mul 6.28 ${WD}/FieldMap.nii.gz + +echo "DONE: preparing General Electric fieldmap" + +if [ ! $GradientDistortionCoeffs = "NONE" ] ; then + ${GlobalScripts}/GradientDistortionUnwarp.sh \ + --workingdir=${WD} \ + --coeffs=${GradientDistortionCoeffs} \ + --in=${WD}/Magnitude \ + --out=${WD}/Magnitude_gdc \ + --owarp=${WD}/Magnitude_gdc_warp + + ${FSLDIR}/bin/applywarp --rel --interp=nn -i ${WD}/Magnitude_brain -r ${WD}/Magnitude_brain -w ${WD}/Magnitude_gdc_warp -o ${WD}/Magnitude_brain_gdc + ${FSLDIR}/bin/fslmaths ${WD}/Magnitude_gdc -mas ${WD}/Magnitude_brain_gdc ${WD}/Magnitude_brain_gdc + ${FSLDIR}/bin/fslmaths ${WD}/Magnitude_brain_gdc -bin -ero -ero ${WD}/Magnitude_brain_gdc_ero + ${FSLDIR}/bin/fslmaths ${WD}/FieldMap -mas ${WD}/Magnitude_brain_gdc_ero -dilM -dilM -dilM -dilM -dilM ${WD}/FieldMap_dil + ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${WD}/FieldMap_dil -r ${WD}/FieldMap_dil -w ${WD}/Magnitude_gdc_warp -o ${WD}/FieldMap_gdc + ${FSLDIR}/bin/fslmaths ${WD}/FieldMap_gdc -mas ${WD}/Magnitude_brain_gdc ${WD}/FieldMap_gdc + + ${FSLDIR}/bin/imcp ${WD}/Magnitude_gdc ${MagnitudeOutput} + ${FSLDIR}/bin/imcp ${WD}/Magnitude_brain_gdc ${MagnitudeBrainOutput} + cp ${WD}/FieldMap_gdc.nii.gz ${FieldMapOutput}.nii.gz +else + ${FSLDIR}/bin/imcp ${WD}/Magnitude ${MagnitudeOutput} + ${FSLDIR}/bin/imcp ${WD}/Magnitude_brain ${MagnitudeBrainOutput} + cp ${WD}/FieldMap.nii.gz ${FieldMapOutput}.nii.gz +fi + +echo " " +echo " END: Field Map Preprocessing and Gradient Unwarping" +echo " END: `date`" >> $WD/log.txt + +########################################## QA STUFF ########################################## + +if [ -e $WD/qa.txt ] ; then rm -f $WD/qa.txt ; fi +echo "cd `pwd`" >> $WD/qa.txt +echo "# Check the brain extraction and distortion correction of the fieldmap magnitude image" >> $WD/qa.txt +echo "fslview ${WD}/Magnitude ${MagnitudeOutput} ${MagnitudeBrainOutput} -l Red -t 0.5" >> $WD/qa.txt +echo "# Check the range (largish values around 600 rad/s) and general smoothness/look of fieldmap (should be large in inferior/temporal areas mainly)" >> $WD/qa.txt +echo "fslview ${FieldMapOutput}" >> $WD/qa.txt + +############################################################################################## + diff --git a/global/scripts/FieldMapPreprocessingAll.sh b/global/scripts/SiemensFieldMapPreprocessingAll.sh similarity index 97% rename from global/scripts/FieldMapPreprocessingAll.sh rename to global/scripts/SiemensFieldMapPreprocessingAll.sh index a6fe398f9..b42c3a5c0 100755 --- a/global/scripts/FieldMapPreprocessingAll.sh +++ b/global/scripts/SiemensFieldMapPreprocessingAll.sh @@ -8,7 +8,8 @@ set -e ################################################ SUPPORT FUNCTIONS ################################################## Usage() { - echo "`basename $0`: Script for generating a fieldmap suitable for FSL, and also do gradient non-linearity distortion correction of these" + echo "`basename $0`: Script for generating a fieldmap suitable for FSL from Siemens Gradient Echo field map," + echo " and also do gradient non-linearity distortion correction of these" echo " " echo "Usage: `basename $0` [--workingdir=]" echo " --fmapmag=" @@ -17,7 +18,7 @@ Usage() { echo " --ofmapmag=" echo " --ofmapmagbrain=" echo " --ofmap=" - echo " [--gdcoeffs=]" + echo " [--gdcoeffs=]" } # function for parsing options