x.heatmapPrep

#!/bin/zsh

# Kimberly Hemmerling 2020

function usage()
{
cat << EOF

DESCRIPTION
Perform erosion and mask creation for concentric masks of the spinal cord according \
to the tissue types. The timeseries for each voxel in each mask are outputted. These \
data can be used to create heatmaps in SCheatmap.m. Additional outputs of \
interest are all erosions and masks created, the log file maskDescriptions.txt, \
and vertebralLevels.txt. The vertebral levels can be used in the plot script to \
indicate the levels.

USAGE
  zsh x.heatmapPrep -i ~/data/sub-01/label -f ~/data/sub-01/func.nii.gz -o heatmap_output [-s] [-k x,y,z]

MANDATORY ARGUMENTS
  -f <functional data>    full path to functional data file (same file used for SCT registration)
  -o <output folder>      desired name of output folder
      **Either -i or -m is mandatory, but including both will only use -m**

OPTIONAL ARGUMENTS
  -i <input folder>       full path to label folder - contains the PAM50 template
                          warped to the functional space (from sct_warp_template)
  -m <mask>               full path to binary mask (spinal cord, etc.) to include
  -s                      smooth data with mask of spinal cord using AFNI 3dblurinmask [Default: no smoothing]
  -k <kernel>             define FWHM of isotropic or anisotropic smoothing kernel [Default: 2,2,6]
                          E.g. anisotropic: -k 2,2,6 ; isotropic -k 3
  -c {0,1}                include CSF masks (doesn't work with vertebral levels) [Default: ask user]


EOF
}
# Show usage if no arguments
if [ ! ${#@} -gt 0 ]; then
    usage `basename ${0}`
    exit 1
fi
# Initialize variables
scriptname=${0}
data_path=
outputdir=
smooth=0
include_CSF=2
fwhm=2,2,6 #####default value , make this compatible
while getopts ':hi:f:o:sk:c:m:' OPTION; do
  case $OPTION in
    h)
          usage
          exit 1
          ;;
    i)
          data_path=$OPTARG
          ;;
    f)
          func=$OPTARG
          ;;
    o)
          output_folder=$OPTARG
          ;;
    s)
          smooth=1
          # echo "We will be smoothing"
          ;;
    k)
          fwhm=${OPTARG}
          ;;
    c)
          include_CSF=${OPTARG}
          ;;
    m)
          sc_mask=${OPTARG}
          ;;
    ?)
          echo "\e[91mUnknown flag!! -$OPTARG\e[0m"
          usage
          exit 1
          ;;
  esac
done
# Check inputs exist
if [[ -z ${data_path} && -z ${sc_mask} ]]; then
	 echo "ERROR: Label folder or mask path not specified. Exiting!\n"
     exit 1
fi
if [[ -z ${func} ]]; then
	 echo "ERROR: Functional data file not specified. Exiting!\n"
     exit 1
fi
if [[ -z ${output_folder} ]]; then
	 echo "ERROR: Output folder not specified. Exiting!\n"
     exit 1
fi
# Check for length of fwhm must be 5 (anisotropic) or 1 (isotropic)
STRLENGTH=$(echo -n $fwhm | wc -m)
if [[ ${STRLENGTH} -ne 1 && ${STRLENGTH} -ne 5 ]]; then
  echo "ERROR: Smoothing kernel should either be comma separated for anisotropic\
  (e.g. -k 2,2,6) or single value for isotropic (e.g. -k 3). Exiting!\n"
    exit 1
fi

# If anisotropic smoothing is requested, check orientation for 3dblurinmask
if [[ ${smooth} -eq 1 && ${STRLENGTH} -eq 5 ]]; then
  func_orientation=`3dinfo -short -orient ${func}`
  echo ""
  echo "The allowed order of coordinates are (1) R-L or L-R, (2) A-P or P-A, (3) I-S or S-I"
  echo "The orientation of the input functional data is: ${func_orientation}"
  if [[ ${func_orientation[@]:0:1} =~ ^(R|L)$ && ${func_orientation[@]:1:1} =~ ^(A|P)$ && ${func_orientation[@]:2:1} =~ ^(I|S)$ ]]; then
      echo "This meets requirements for for anisotropic smoothing."
      echo ""
  else
      echo -n "This does not meet requirements for for anisotropic smoothing.\
Would you like to continue without smoothing? Enter y or n: "
      read CONT
      if [[ $CONT =~ ^([yY])$ ]]; then
        smooth=0
        echo "OK, will continue without smoothing!\n"
      elif [[ $CONT =~ ^([nN])$ ]]; then
        echo -e "\e[91mExiting!\e[0m\n"
        exit 1
      else
        echo -e "\e[91mInput not recognized. Exiting!\e[0m\n"
      fi
  fi
fi
################################################################################
# End of checks, beginning analysis
################################################################################
# # Change directory to func folder and create output directory - need to do this mask thing...
# cd ${data_path}
# echo "\nLooking at functional data in... ${data_path} \n"
# mkdir ../${output_folder}
# echo '\nCreated new output directory called: ' ${output_folder} '\n'
# cd ../${output_folder}


# p=abc/def/ghi/jkl.nii.gz # func
func_dir=`dirname ${func}`
echo "${func_dir}"
mkdir ${func_dir}/${output_folder}
cd ${func_dir}/${output_folder}
outputdir=`pwd`

# Do checks then
if [[ -f "${outputdir}/maskDescriptions.txt" ]]; then
    echo -n "Looks like you've already done this analysis, continuing will cause files to be overwritten. Would you like to continue? \e[42mEnter y or n:\e[0m "
    read CONT
    if [[ $CONT =~ ^([yY])$ ]]; then
      echo "OK, let's keep going!\n"
    elif [[ $CONT =~ ^([nN])$ ]]; then
      echo -e "\e[91mExiting!\e[0m\n"
      exit 1
    else
      echo -e "\e[91mInput not recognized. Exiting!\e[0m\n"
    fi
fi
# Initialize mask descriptions log file
date > ${outputdir}/maskDescriptions.txt
# Check if mask input exists to do ANALYSIS IN MASK
if [[ -n ${sc_mask} ]]; then
   echo "Will output tissue timeseries within provided mask.\n"
   # Add to mask descriptions log file
   echo "Path to functional data input: ${func}" >> ${outputdir}/maskDescriptions.txt
   echo "Path to output directory: ${outputdir}" >> ${outputdir}/maskDescriptions.txt
   echo "" >> ${outputdir}/maskDescriptions.txt

   # Re-Save input mask
   m=00
   cp ${sc_mask} ${outputdir}/mask${m}.nii.gz
   echo "Input mask is mask${m}.nii.gz"
   # Calculate timeseries for mask:
   fslmeants -i ${func} -m ${sc_mask} --showall > ${outputdir}/mask${m}ts.txt
   echo "INPUT MASK" >> ${outputdir}/maskDescriptions.txt
   echo "mask${m}.nii.gz" >> ${outputdir}/maskDescriptions.txt
   # Calculate DVARS trace using PAM50_cord as the mask
   3dTto1D -input ${func} -mask ${sc_mask} -method dvars -prefix ${outputdir}/dvars.txt

     exit 1
fi

cd ${data_path}
# All inputs are data in functional space
# Calculate functional mean image for visualizations later
fslmaths ${func} -Tmean func_mean.nii.gz
# The data_path should be the "label" folder containing these PAM50 files
cord=template/PAM50_cord.nii.gz # cord mask
levels=template/PAM50_levels.nii.gz # vertebral levels (NOT CORD)
# Add to mask descriptions file
echo "Path to analyzed data: ${data_path}" >> ${outputdir}/maskDescriptions.txt
echo "Path to functional data input: ${func}" >> ${outputdir}/maskDescriptions.txt
echo "Path to output directory: ${outputdir}" >> ${outputdir}/maskDescriptions.txt
echo "" >> ${outputdir}/maskDescriptions.txt
# Output vertebral levels text file
fslmaths ${levels} -thr 0.5 -bin ${outputdir}/tempLevelsToOne.nii.gz
fslmeants -i ${levels} -m ${outputdir}/tempLevelsToOne.nii.gz --showall > ${outputdir}/vertebralLevels.txt
rm ${outputdir}/tempLevelsToOne.nii.gz
# Need to binarize these tissue specific (probabilistic) masks
fslmaths template/PAM50_csf.nii.gz -thr 0.5 -bin ${outputdir}/full_csf_mask.nii.gz
fslmaths template/PAM50_gm.nii.gz -thr 0.5 -bin ${outputdir}/full_gm_mask.nii.gz
fslmaths template/PAM50_wm.nii.gz -thr 0.5 -bin ${outputdir}/full_wm_mask.nii.gz
csf=${outputdir}/full_csf_mask.nii.gz
gm=${outputdir}/full_gm_mask.nii.gz
wm=${outputdir}/full_wm_mask.nii.gz
###############################################################################
# OPTIONAL SECTION FOR CSF EROSION AND MASK CREATION:
# Option is prompted in command window for user input. More detailed comments are
# provided in the WM section below.
###############################################################################
o=00 #output
i=00 #input
m=00 #mask
if [[ ${include_CSF} -eq 2 ]]; then
    echo -n "Creating tissue type masks. GM/WM masks will be created. Do you want to include CSF? \e[42mEnter y or n:\e[0m "
    read YN
    if [[ $YN =~ ^([yY])$ ]]; then
      include_CSF=1
    elif [[ $YN =~ ^([nN])$ ]]; then
      include_CSF=0
      echo "\n OK, no CSF. \n"
    else
      echo -e "\e[91mInput not recognized. Exiting!\e[0m\n"
      exit 1
    fi
fi
if [[ ${include_CSF} -eq 1 ]]; then
  echo "\n OK, will create CSF masks..."
  echo -e "\e[1;35m    Edit CSF mask. Fill in top/bottom slices inside CSF \e[0m"
  echo -e "\e[1;35m    Save to same file name. Close window when done. \n \e[0m"
  cp ${outputdir}/full_csf_mask.nii.gz ${outputdir}/csfAndCord_mask.nii.gz
  csfCord=${outputdir}/csfAndCord_mask.nii.gz
  fsleyes ${outputdir}/csfAndCord_mask.nii.gz -cm red-yellow
  fslmaths ${outputdir}/csfAndCord_mask.nii.gz -fillh ${outputdir}/csfAndCord_mask.nii.gz
  fsleyes ${outputdir}/csfAndCord_mask.nii.gz -cm red-yellow &
  # After this, ${csfCord} is really csf+gm+wm
  # Erode and make masks (descriptions for each of these steps in WM section)
  # Use cord for subtraction (input)
  # In sct 4.2.2 sct_maths used : -erode 3,3,1. Updated to sct v. 5.2.0 : -erode 1 -shape disk -dim 2
  sct_maths -i ${csfCord} -o ${outputdir}/erosion${o}.nii.gz -erode 1 -shape disk -dim 2
  sct_maths -i ${csfCord} -o ${outputdir}/mask_TEMP.nii.gz -sub ${outputdir}/erosion${o}.nii.gz
  sct_maths -i ${outputdir}/mask_TEMP.nii.gz -o ${outputdir}/sub_cord_TEMP.nii.gz -sub ${cord}
  fslmaths ${outputdir}/sub_cord_TEMP.nii.gz -thr 0 -bin ${outputdir}/mask${o}.nii.gz
  fslstats ${outputdir}/mask${o}.nii.gz -m >> ${outputdir}/mask_averages.txt
  mask_avg=$( tail -1 ${outputdir}/mask_averages.txt )
  rm ${outputdir}/sub_cord_TEMP.nii.gz ${outputdir}/mask_TEMP.nii.gz
  # Prepare file to store descriptions of masks
  echo "CEREBROSPINAL FLUID" >> ${outputdir}/maskDescriptions.txt
  echo "mask${o}.nii.gz" >> ${outputdir}/maskDescriptions.txt
  while (( ${mask_avg} > 0 ))
  do
    o=$(( o+1 ))
    if [[ ${o} -le 9 ]]; then
      o="0"${o}
    fi
    sct_maths -i ${outputdir}/erosion${i}.nii.gz -o ${outputdir}/erosion${o}.nii.gz -erode 1 -shape disk -dim 2
    sct_maths -i ${outputdir}/erosion${i}.nii.gz -o ${outputdir}/mask_TEMP.nii.gz -sub ${outputdir}/erosion${o}.nii.gz
    sct_maths -i ${outputdir}/mask_TEMP.nii.gz -o ${outputdir}/sub_cord_TEMP.nii.gz -sub ${cord}
    fslmaths ${outputdir}/sub_cord_TEMP.nii.gz -thr 0 -bin ${outputdir}/mask${o}.nii.gz
    fslstats ${outputdir}/mask${o}.nii.gz -m >> ${outputdir}/mask_averages.txt
    mask_avg=$( tail -1 ${outputdir}/mask_averages.txt )
    rm ${outputdir}/sub_cord_TEMP.nii.gz ${outputdir}/mask_TEMP.nii.gz
    i=$(( i+1 ))
    if [[ ${i} -le 9 ]]; then
      i="0"${i}
    fi
    echo "CSF mask ${i} with average: " ${mask_avg}
    echo "mask${i}.nii.gz" >> ${outputdir}/maskDescriptions.txt
  done
  # TBD delte next 3 lines
  echo -e "\e[1;35mFYI: this mask should be empty. Check file: checkIfEmpty_csf.nii.gz!!! \e[0m"
  echo -e "\e[1;35mThis file is left here to ensure no data is missed./debug. \e[0m"
  cp ${outputdir}/mask${i}.nii.gz ${outputdir}/checkIfEmpty_csf.nii.gz
  rm ${outputdir}/mask${i}.nii.gz ${outputdir}/erosion${i}.nii.gz
  # Remove last line of maskDescriptions.txt file (since it should be empty)
  sed "s/mask${i}.nii.gz/------------/" ${outputdir}/maskDescriptions.txt > ${outputdir}/tempTemp.txt
  mv ${outputdir}/tempTemp.txt ${outputdir}/maskDescriptions.txt
  # Make actual full CSF mask for later optional smoothing within mask
  # rm ${outputdir}/full_csf_mask.nii.gz
  sct_maths -i ${csfCord} -o ${outputdir}/temp_csf_mask1.nii.gz -sub ${wm}
  sct_maths -i ${outputdir}/temp_csf_mask1.nii.gz -o ${outputdir}/temp_csf_mask2.nii.gz -sub ${gm}
  fslmaths ${outputdir}/temp_csf_mask2.nii.gz -thr 0 -bin ${outputdir}/full_csf_mask.nii.gz
  rm ${outputdir}/temp_csf_mask1.nii.gz ${outputdir}/temp_csf_mask2.nii.gz
  csf=${outputdir}/full_csf_mask.nii.gz
elif [[ include_CSF -eq 0 ]]; then
  echo "\n OK, no CSF. \n"
fi
###############################################################################
# WM EROSION AND MASK CREATION:
###############################################################################
# Create first mask using full cord mask:
# Erode voxel layer from outer edges
sct_maths -i ${cord} -o ${outputdir}/erosion${o}.nii.gz -erode 1 -shape disk -dim 2
# Subtract erosion from full mask to get only outer edge mask
sct_maths -i ${cord} -o ${outputdir}/mask_TEMP.nii.gz -sub ${outputdir}/erosion${o}.nii.gz
# Subtract gm mask from mask to have no gm/wm overlap
sct_maths -i ${outputdir}/mask_TEMP.nii.gz -o ${outputdir}/sub_gm_TEMP.nii.gz -sub ${gm}
# Threshold out negative values (from previous step)
fslmaths ${outputdir}/sub_gm_TEMP.nii.gz -thr 0.1 -bin ${outputdir}/mask${o}.nii.gz
# Save average value in mask to check if it's empty
fslstats ${outputdir}/mask${o}.nii.gz -m >> ${outputdir}/mask_averages.txt
# Keep current mask avg as variable
mask_avg=$( tail -1 ${outputdir}/mask_averages.txt )
# Remove temp files
rm ${outputdir}/sub_gm_TEMP.nii.gz ${outputdir}/mask_TEMP.nii.gz
# Prepare file to store descriptions of masks
echo "WHITE MATTER" >> ${outputdir}/maskDescriptions.txt
echo "mask${o}.nii.gz" >> ${outputdir}/maskDescriptions.txt
# Continue to create masks by eroding and subtracting the eroded mask iteratively
# until the center is reached. The average within a mask signifies whether the
# mask is empty yet or not, i.e. whether we're done with this section.
while (( ${mask_avg} > 0 ))
do
  o=$(( o+1 ))
  if [[ ${o} -le 9 ]]; then
    o="0"${o}
  fi
  sct_maths -i ${outputdir}/erosion${i}.nii.gz -o ${outputdir}/erosion${o}.nii.gz -erode 1 -shape disk -dim 2
  sct_maths -i ${outputdir}/erosion${i}.nii.gz -o ${outputdir}/mask_TEMP.nii.gz -sub ${outputdir}/erosion${o}.nii.gz
  sct_maths -i ${outputdir}/mask_TEMP.nii.gz -o ${outputdir}/sub_gm_TEMP.nii.gz -sub ${gm}
  fslmaths ${outputdir}/sub_gm_TEMP.nii.gz -thr 0.1 -bin ${outputdir}/mask${o}.nii.gz
  fslstats ${outputdir}/mask${o}.nii.gz -m >> ${outputdir}/mask_averages.txt
  mask_avg=$( tail -1 ${outputdir}/mask_averages.txt )
  rm ${outputdir}/sub_gm_TEMP.nii.gz ${outputdir}/mask_TEMP.nii.gz
  i=$(( i+1 ))
  if [[ ${i} -le 9 ]]; then
    i="0"${i}
  fi
  echo "WM mask ${i} with average: " ${mask_avg}
  echo "mask${i}.nii.gz" >> ${outputdir}/maskDescriptions.txt
done
echo -e "\e[1;35mFYI: this mask should be empty. Check file: checkIfEmpty_wm.nii.gz!!! \e[0m"
echo -e "\e[1;35mIf not empty... need to troubleshoot/debug. \e[0m"
cp ${outputdir}/mask${i}.nii.gz ${outputdir}/checkIfEmpty_wm.nii.gz
# Remove last line of maskDescriptions.txt file (since it should be empty)
sed "s/mask${i}.nii.gz/------------/" ${outputdir}/maskDescriptions.txt > ${outputdir}/tempTemp.txt
mv ${outputdir}/tempTemp.txt ${outputdir}/maskDescriptions.txt
# Calculate timeseries for each mask:
m=00
for mask in ${outputdir}/mask*.nii.gz
do
	echo "Processing mask $m"
  fslmeants -i ${func} -m ${mask} --showall > ${outputdir}/mask${m}ts.txt
  m=$(( m+1 ))
  if [[ ${m} -le 9 ]]; then
    m="0"${m}
  fi
done
###############################################################################
# Gray matter:
###############################################################################
# No erosion needed for GM because it is quite small already
m=$(( m-1 ))
if [[ ${m} -le 9 ]]; then
  m="0"${m}
fi
# Save gray matter mask as last mask #
cp ${gm} ${outputdir}/mask${m}.nii.gz
echo "GM mask is mask${m}.nii.gz"
# Calculate timeseries for mask:
fslmeants -i ${func} -m ${gm} --showall > ${outputdir}/mask${m}ts.txt
echo "GRAY MATTER" >> ${outputdir}/maskDescriptions.txt
echo "mask${m}.nii.gz" >> ${outputdir}/maskDescriptions.txt
###############################################################################
# Create mask visualization file allMaskVisualization.nii.gz
v=00
cp ${outputdir}/mask00.nii.gz ${outputdir}/allMaskVisualization.nii.gz
for mask in ${outputdir}/mask*.nii.gz
do
  if (( $v == 00 )); then
    # Skip mask0 since it was already included as the base
    v=$(( v+1 ))
    if [[ ${v} -le 9 ]]; then
      v="0"${v}
    fi
  else
    v=$(( v+1 ))
    if [[ ${v} -le 9 ]]; then
      v="0"${v}
    fi
    echo "${v}"
    echo "\n ${mask}"
    fslmaths ${mask} -mul $v -add ${outputdir}/allMaskVisualization.nii.gz ${outputdir}/temp.nii.gz
    cp ${outputdir}/temp.nii.gz ${outputdir}/allMaskVisualization.nii.gz
  fi
done
###############################################################################
# Calculate DVARS trace using PAM50_cord as the mask
3dTto1D -input ${func} -mask ${cord} -method dvars -prefix ${outputdir}/dvars.txt
###############################################################################
# Extract 8 largest tracts by volume, binarize then extract ts (WIP)
# Binarize masks
mkdir ${outputdir}/tracts
fslmaths atlas/PAM50_atlas_04.nii.gz -thr 0.5 -bin ${outputdir}/tracts/LLC_mask.nii.gz # Left lateral corticospinal
fslmaths atlas/PAM50_atlas_05.nii.gz -thr 0.5 -bin ${outputdir}/tracts/RLC_mask.nii.gz # Right lateral corticospinal
fslmaths atlas/PAM50_atlas_22.nii.gz -thr 0.5 -bin ${outputdir}/tracts/LVC_mask.nii.gz # Left ventral corticospinal
fslmaths atlas/PAM50_atlas_23.nii.gz -thr 0.5 -bin ${outputdir}/tracts/RVC_mask.nii.gz # Right ventral corticospinal
fslmaths atlas/PAM50_atlas_02.nii.gz -thr 0.5 -bin ${outputdir}/tracts/LFC_mask.nii.gz # Left fasciculus cuneatus
fslmaths atlas/PAM50_atlas_03.nii.gz -thr 0.5 -bin ${outputdir}/tracts/RFC_mask.nii.gz # Right fasciculus cuneatus
fslmaths atlas/PAM50_atlas_00.nii.gz -thr 0.5 -bin ${outputdir}/tracts/LFG_mask.nii.gz # Left fasciculus gracilis
fslmaths atlas/PAM50_atlas_01.nii.gz -thr 0.5 -bin ${outputdir}/tracts/RFG_mask.nii.gz # Right fasciculus gracilis
# Extract timeseries
fslmeants -i ${func} -m ${outputdir}/tracts/LLC_mask.nii.gz --showall > ${outputdir}/tracts/LLCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/RLC_mask.nii.gz --showall > ${outputdir}/tracts/RLCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/LVC_mask.nii.gz --showall > ${outputdir}/tracts/LVCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/RVC_mask.nii.gz --showall > ${outputdir}/tracts/RVCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/LFC_mask.nii.gz --showall > ${outputdir}/tracts/LFCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/RFC_mask.nii.gz --showall > ${outputdir}/tracts/RFCts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/LFG_mask.nii.gz --showall > ${outputdir}/tracts/LFGts.txt
fslmeants -i ${func} -m ${outputdir}/tracts/RFG_mask.nii.gz --showall > ${outputdir}/tracts/RFGts.txt
###############################################################################
# SMOOTHING USING AFNI 3dblurinmask
###############################################################################
if [[ ${smooth} -eq 1 ]]; then
  echo "\nSmoothing was requested. Beginning smoothing... \n"
  # Blur in cord mask (Find a way to automate this (maybe fslroi?)):
  echo "\nDelete top and bottom slice in mask (for smoothing).\nOverwrite file when saving (PAM50_cord_forSmoothing.nii.gz).\n"
  cp ${cord} template/PAM50_cord_forSmoothing.nii.gz
  blurMask=template/PAM50_cord_forSmoothing.nii.gz
  fsleyes ${blurMask}
  # Loop through -k cord to retrive x,y,z FWHM kernel values
  # MAKE THIS SO THAT IT ONLY GOES THROUGH THIS LOOP IF THE INPUT IS LENGH STRLENGTH=5
  if [[ ${smooth} -eq 1 && ${STRLENGTH} -eq 5 ]]; then
    # If anisotropic:
    for i in $(echo $fwhm | sed "s/,/ /g")
    do
       args+=("$i")
    done
    # @ means all, 0 is offset, 1 is how many desired (accessible for bash & zsh)
    FWHMx=${args[@]:0:1}
    FWHMy=${args[@]:1:1}
    FWHMz=${args[@]:2:1}
    echo "FWHM split: X:${FWHMx}, Y:${FWHMy}, Z:${FWHMz}"
    echo "\nBlurring func data in mask (anisotropic)...\n"
    3dblurinmask -input ${func} -FWHMxyz ${FWHMx} ${FWHMy} ${FWHMz} -Mmask  ${blurMask} \
    -prefix ${outputdir}/func_blur_cord.nii.gz -preserve -quiet
    echo "Anisotropic smoothing complete."
  elif [[ ${smooth} -eq 1 && ${STRLENGTH} -eq 1 ]]; then
    # If isotropic:
    FWHMxyz=${fwhm}
    echo "FWHM: ${FWHMxyz}"
    echo "\nBlurring func data in mask (isotropic)...\n"
    3dblurinmask -input ${func} -FWHM ${FWHMxyz} -Mmask  ${blurMask} \
    -prefix ${outputdir}/func_blur_cord.nii.gz -preserve -quiet
    echo "Isotropic smoothing complete."
  fi
  # Calculate blurred timeseries for each mask (use for restricted or undrestricted methods)
  m=00
  for mask in ${outputdir}/mask*.nii.gz
  do
  	echo "Processing mask $m in func_blur_cord.nii.gz"
    fslmeants -i ${outputdir}/func_blur_cord.nii.gz -m ${mask} --showall > ${outputdir}/blur_mask${m}ts.txt
    m=$(( m+1 ))
    if [[ ${m} -le 9 ]]; then
      m="0"${m}
    fi
  done
fi

# # Blur unrestricted (delete later):
# 3dblurinmask -input ${func} -FWHM 2.9 \
# -prefix ${outputdir}/func_blur.nii.gz -preserve -quiet

# # Blur with SCT function (smooth is only supported at 3D)
# # Make new folder for splitting the functional data
# mkdir ${outputdir}/split
# cd ${outputdir}/split
# # Split data into volumes for smoothing
# sct_image -i ${func} -split t -o func_vol.nii.gz
# v=0
# concatInput=
# for FILE in *
# do
#   echo "Working on ${FILE}..."
#   echo ${concatInput} #Delete this lline later
#   # Get centerline of each 3D image and ouptut numbered volumes
#   printf -v vol "%04d" ${v}
#   sct_get_centerline -i ${FILE} -c t2s -o center_T${vol} -v 0
#   # Make string for input to concatenation
#   concatInput="${concatInput} func_vol_T${vol}_smooth.nii.gz" # NEED TO SWITCH THIS TO SMOOTHING files
#   # Smooth each 3D image (SCT 5.2.0 has -o option)
#   sct_smooth_spinalcord -i ${FILE} -s center_T${vol}.nii.gz -smooth 2,2,6 -v 0
#   v=$(( v+1 ))
# done
#
# # Concatenate 3D volumes to 4D timeseries (THIS STEP DOES NOT WORK)
# sct_image -i ${value} -concat t -o func_SCTsmoothed.nii.gz
# mv func_SCTsmoothed.nii.gz ${outputdir}
# cd ${data_path}

# # Calculate blurred timeseries for each mask (use for SCT method)
# m=0
# for mask in ${outputdir}/mask*.nii.gz
# do
# 	echo "Processing mask $m in func_SCTsmoothed.nii.gz"
#   fslmeants -i ${outputdir}/func_SCTsmoothed.nii.gz -m ${mask} --showall > ${outputdir}/blur_mask${m}ts.txt
#   m=$(( m+1 ))
# done
###############################################################################
# Tissue type mask visualization
fsleyes func_mean.nii.gz ${outputdir}/allMaskVisualization.nii.gz -cm brain_colours_greengray &
###############################################################################
# Vertebral levels mask visualization
fsleyes func_mean.nii.gz ${levels} -cm blue-lightblue &
###############################################################################
echo ""
echo -e "\e[1;35m... done! Runtime: $SECONDS seconds.\e[0m"
echo ""
echo "Note: the checkIfEmpty files can be deleted if empty.\n"