nuisance_regression

#!/bin/bash

FSLCLUSTER=fsl-cluster
command -v $FSLCLUSTER >&/dev/null || FSLCLUSTER=cluster

function warp_structmask_to_func {
    structvol="${1}"
    outvol="${2}"

    # need to know the actual file for afni. fsl is cool without the extention
    local templateBrainfile="$templateBrain"
    [ ! -f "$templateBrainfile" ] && templateBrainfile="$templateBrain.nii"
    [ ! -f "$templateBrainfile" ] && templateBrainfile="$templateBrain.nii.gz"
    [ ! -f "$templateBrainfile" ] && rel "ERROR: no file for templateBrain '$templateBrain'" c && exit 1

    # ignore structvol and just use template... assuming that's why we're not warping but still want nuissance regressor
    if [ "x$no_warp" == "x1" -a ! -d unwarp/ ] && all_same_grid "$structvol.nii.gz" "$templateBrainfile"; then
       # all_same_grid from helper_functions
       rel "$FUNCNAME: no warp but want regressor w/o fieldmap and structvol is template" c
       rel "see struct vs template: 3dinfo -same_all_grid '$structvol.nii.gz' '${templateBrain}'" c
       rel "3dcopy $structvol.nii.gz $outvol.nii.gz" # NB. assuming ext is always .nii.gz
    #warp anatomical mask into functional space
    elif [ $( imtest "unwarp/T1_UD_warp_bbr" ) -eq 1 ]; then
	rel "Using inverse of unwarp/T1_UD_warp_bbr to warp $structvol into functional space with distortion" c
	[ ! -f "unwarp/struct_to_func_transform_and_distort${ext}" ] && rel "invwarp -w unwarp/T1_UD_warp_bbr -o unwarp/struct_to_func_transform_and_distort -r unwarp/EF_D_mc_target"
	rel "applywarp --in=\"${structvol}\" --out=\"${outvol}\" --interp=spline --ref=mc_target --warp=unwarp/struct_to_func_transform_and_distort"
    else
       # TODO: test if all_same_grid $structvol as $func
       # otherwise something is wrong
	rel "Using struct_to_func transform to warp ${structvol} mask into functional space" c
	rel "applywarp --in=\"${structvol}\" --out=\"${outvol}\" --interp=spline --ref=mc_target --premat=$funcdir/transforms/struct_to_func.mat"
    fi

    rel "fslmaths ${outvol} -thr 0.5 -bin ${outvol} -odt char" #force to binary mask to remove tiny values due to interpolation
}

#compute motion parameter derivatives for nuisance regression
# a function b/c we call it twice
motionderiv() {
  [ ! -f ".motion_deriv" ] && \
  1d_tool.py -overwrite -infile motion.par -set_nruns 1 -derivative -demean -write .motion_deriv

  return 0
}

COMPUTE_NUISANCE_REGRESSORS_GLOBALS=(nuisance_file nuisance_regressors no_warp \
   postWarp postDespike postSS \
   templateName mprageBet_base ext tr prefix bpLow bpHigh \
   gsr_in_prefix custom_regprefix)

# give
# - a nifti get basis functions to remove autocorrelation
# - output
# - optional: mask
# write file to $out with 3dREMLfit to bandpass, remove regressors, remove autocorralation
# USAGE:
# run_3dremlfit dnswdktm_restepi_5.nii.gz Rdnswdktm_restepi_5.nii.gz wktm_restepi_98_2_mask_dil1x_templateTrim.nii.gz
# GLOBAL $tr, $nuisance_file, $bpLow, $bpHigh
# there is lots of boilerplate checks. at the heart this runs 2 commands
# 1) 1dBport to generate bandpass sine/cosine basis functions
# 2) 3dREMLfit to remove these and autocorrelation using ARMA(1,1) model
function run_3dremlfit(){
  img=$1
  out=$2
  mask=$3

  if [ $rmautocorr -eq 1 ]; then
      local outmethod="-Rwherr"
  else
      local outmethod="-Rerrts"
  fi

  [ ! -r $img ]  && echo "$FUNCNAME missing image ($img)" >&2 && return 1
  nvol=$(3dinfo -nt $img)
  [ -z "$nvol" -o $nvol -lt 1 ] && echo "$FUNC_NAME: bad number of timesteps ($nvol) in $img" >&2 && return 1

  #always include baseline (const), linear, and quadratic terms in removal
  Rscript -e "suppressMessages(library(orthopolynom)); \
      unnormalized.p.list <- legendre.polynomials(2, normalized=FALSE); \
      baseline <- polynomial.values(polynomials=unnormalized.p.list, x=seq(-1,1, length.out=$nvol)); \
      baseline[2:length(baseline)] <- lapply(baseline[2:length(baseline)], function(v) { v - mean(v) }); \
      write.table(do.call(cbind, baseline), file='.legendre.1D', col.names=FALSE, row.names=FALSE)"  
  
  bandoption=""
  if [ "$bandpass_filter" -eq 1 ]; then
      if [ -z "$bpLow" -o -z "$bpHigh" ]; then
	  rel "ERROR! $FUNC_NAME has no bandpass low or high parameters, but we said bandpass!" c
	  return 1
      fi
      bandoption="-band $bpLow $bpHigh"
      #MH 30May2017: for reasons I don't know, 1dBport dies when TR is < 1s and is not prefixed by zero (e.g., .894).
      local ztr="${tr/#./0.}"
      rel "1dBport -nodata $nvol $ztr $bandoption -noconst -invert > .basis.1D" #constant handled above in baseline Rscript call
      local filter_file=".basis.1D"
  else
      #MH Jul2018: I am not aware of a reason why bandpassing + ARMA regression should go together
      #rel "WARNING! AFNI gurus say you probably WANT bandpassing. You are NOT bandpassing (by setting no low pass and highpass above nyquist" c      
      #rel "ERROR! I think you need to bandpass to remove autocorrelation. Sorry." c
      #return 1
      local filter_file=""
  fi
  
  reml_in_file=$nuisance_file.3dREMLfit
  
  if [ "$nuisance_regression" -eq 1 ]; then
      rel "paste $nuisance_file .legendre.1D ${filter_file} > $reml_in_file"
  else
      rel "paste .legendre.1D ${filter_file} > $reml_in_file"
  fi

  # could also try below, but get different results!
  # regoption="-stimbase $nuisance_file"

  # are we using a mask?
  maskoption=""
  [ -n "$mask" ] && imtest "$mask" && maskoption="-mask $mask"

  rel "3dREMLfit -input $img $maskoption -matim $reml_in_file $outmethod $out" #-Rwherr returns the pre-whitened residuals; -Rerrts returns conventional residuals Y - \hat{Y}
  [ $rmautocorr -eq 1 ] && rel "date > .remove_autocorr_complete"
  return 0
}

function compute_nuisance_regressors {
    # 20170421WF: skipping might be a problem if we change what regressors we want
    if [ -r .nuisance_compute_complete${custom_regprefix:-} ]; then
       rel "already computed nuisance regression" c
       return 0
    fi

    declare -a nuisancecols #nuisance regressors to be concatenated

    #if despiking was applied, use this file for extracting CSF, WM, GS time series. Otherwise, use the file after skull strip, motion, and slice timing
    [ -n "${postDespike}" ] && local prewarp_func_ts="${postDespike}" || local prewarp_func_ts="${postSS}"
	
    if [ $no_warp -eq 0 ]; then
	#extracting time series from post-warp files to leverage tissue probability maps for templates
	local func_ts="${postWarp}"	
    else
	local func_ts="${prewarp_func_ts}" #no warped file to refer to
    fi

    # 2017114 - c/o BTC - dont want to put aroma removed noise back in 
    if [ "$ica_aroma" -eq 1 ]; then
      local func_ts="$postAroma"
    fi

    #parse comma-delimited regressors of interest into array
    OLDIFS="${IFS}"
    IFS=',' read -ra reg <<< "$nuisance_regressors"
    IFS="${OLDIFS}"
    #reg=$( printf "%s\n" "${reg[@]}" | sort -u ) #eliminate any duplicate regressors
    #switched to perl because it does not reorder regressors
    reg=$( printf "%s\n" "${reg[@]}" | perl -ne '$a{$_}++; END{print join"\n",(sort keys %a)}' ) #eliminate any duplicate regressors

    for r in $reg; do
	if [[ "$r" =~ ^(6motion|rx|ry|rz|tx|ty|tz)$ && ! -f ".motion_demean" ]]; then 
	    #compute de-meaned motion parameters for nuisance regression
	    rel "1d_tool.py -overwrite -infile motion.par -set_nruns 1 -demean -write .motion_demean"
	fi
	
	if [[ "$r" =~ ^(d6motion|drx|dry|drz|dtx|dty|dtz)$ && ! -f ".motion_deriv" ]]; then 
	    #compute motion parameter derivatives for nuisance regression
            motionderiv 
	fi
	
	# quad (q6motion, qrx,qrt...)
	if [[ "$r" =~ ^q(6m|r|t) && ! -f ".motion_quad" ]]; then 
      	    [ ! -L .motion.1D ] && ln -s motion.par .motion.1D
      	    3dcalc -a .motion.1D\' -expr 'a**2' -prefix - |
		1d_tool.py -transpose -infile - -write .motion_quad
	fi
	
	# quad of deriv (qd6motion, qdrx,qdrt...)
	if [[ "$r" =~ ^qd(6m|r|t) && ! -f ".motion_deriv_quad" ]]; then 
            # need motion deriv to do quad on it
            motionderiv 
	    
      	    [ ! -L .motion_dev.1D ] && ln -s .motion_deriv .motion_dev.1D
      	    3dcalc -a .motion_dev.1D\' -expr 'a**2' -prefix - |
      		1d_tool.py -transpose -infile - -write .motion_deriv_quad 
	fi
	
	if [[ "$r" =~ ^(csf|dcsf|csf[0-9]+|dcsf[0-9]+)$ && ! -f ".csf_ts" ]]; then
	    #technically should probably change this to allow for use of templates for -warp_compute
	    if [ $no_warp -eq 0 ]; then
		#if we have a warp to standard space then use the inverse warp to transform the csf mask from the probabilistic atlas to subject space
		if [ $templateName = MNI_3mm ]; then
		    #thr 0.95 without erosian is best: 95 voxels
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_csf_tal_nlin_asym_09c_3mm\" -thr 0.95 -bin .template_csf_prob"
		elif [ $templateName = MNI_2.3mm ]; then
		    #thr 0.98 without erosion works well here: 274 voxels
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_csf_tal_nlin_asym_09c_2.3mm\" -thr 0.98 -bin .template_csf_prob"
		elif [ $templateName = MNI_2mm ]; then
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_csf_tal_nlin_asym_09c_2mm\" -thr 0.98 -bin .template_csf_prob"
		elif [ $templateName = SPM_2mm ]; then
		    #thr 0.6 works best here: 193 voxels
		    rel "fslmaths \"$stddir/spm8_mni/csf.nii\" -thr 0.6 -bin -eroF .template_csf_prob"	    
		elif [ $templateName = MNI_FSL_3mm ]; then
		    #for the FSL MNI, thr 0.7 works well visually: 87 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_csf_3mm.nii\" -thr 0.7 -bin .template_csf_prob"	    
		elif [ $templateName = MNI_FSL_2.3mm ]; then
		    #for the FSL MNI 2.3mm, thr 0.73 works well visually: 189 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_csf_2.3mm.nii\" -thr 0.73 -bin .template_csf_prob"
		elif [ $templateName = MNI_FSL_2mm ]; then
		    #for the FSL MNI 2mm, thr 0.75 works well visually: 193 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_csf_2mm.nii\" -thr 0.75 -bin .template_csf_prob"
		elif [ $templateName = 1YO_2mm ]; then
		    # numbered 1 to N for likelyhood?
		    rel "fslmaths \"$stddir/UNCInfant/2mm_1YO/infant-1yr-csf.nii.gz\" -thr 0.99 -bin -eroF .template_csf_prob"
		else
		    rel "Unsupported template brain for CSF extraction: $templateName" c
		    exit 1
		fi

		local csfmask=.template_csf_prob${ext}    
	    else
		#use the subject's segmented anatomical scan to identify the csf voxels
		
		#fast segmentation should have been run during preprocessMprage
		#use pve_0 for csf
		#re-run FAST here if needed
		[ $( imtest "${mprageBet_base}_fast_seg_0" ) -eq 0 ] && rel "fast -g -o \"${mprageBet_base}_fast\" \"${mprageBet_base}\""

		warp_structmask_to_func "${mprageBet_base}_fast_pve_0" .csf_pve

		#threshold at 0.8 to ensure that only voxels likely to be csf are retained	    
		rel "fslmaths .csf_pve -thr 0.8 -bin .csf_pve_thr0p8"

		#erode csf mask once to reduce risk of partial volume
		rel "3dmask_tool -overwrite -input .csf_pve_thr0p8${ext} -dilate_result -1 -prefix .csf_pve_thr0p8_ero1${ext}"

		#cluster csf mask and only retain clusters of 5 or more contiguous voxels
		rel "$FSLCLUSTER -i .csf_pve_thr0p8_ero1 -t 1.0 --osize=.csf_pve_clustsize --no_table"
		rel "fslmaths .csf_pve_clustsize -thr 5.01 -bin -mas ${postSS}_tmean_mask .csf_mask -odt char" #multiply by tight brain mask to eliminate non-brain voxels
		
		#if fewer than 20 CSF voxels (in functional space) are present, revert to a non-dilated csf mask
		#too few voxels may result in undue influece of a given CSF voxel on the mask average
		local nvox=$( 3dBrickStat -non-zero -count .csf_mask${ext} )

		if [ $nvox -lt 20 ]; then
		    rel "Warning: fewer than 20 CSF voxels were found in .csf_mask. Reverting to a non-dilated mask" c	
		    rel "$FSLCLUSTER -i .csf_pve_thr0p8 -t 1.0 --osize=.csf_pve_clustsize_noero --no_table"
		    rel "fslmaths .csf_pve_clustsize_noero -thr 15.01 -bin -mas ${postSS}_tmean_mask .csf_mask_noero -odt char"

		    local csfmask=.csf_mask_noero${ext}
		else
		    local csfmask=.csf_mask${ext}
		fi
		
		#cleanup some intermediate files
		rel "imrm .csf_pve_thr0p8 .csf_pve_thr0p8_ero1 .csf_pve_clustsize .csf_pve_clustsize_noero"
	    fi
	    
	    #extract CSF time series. If a number follows the CSF specification, this is an aCompCor approach with the specified number of components
	    if [[ "$r" =~ ^(csf[0-9]+|dcsf[0-9]+)$ ]]; then
		rel "3dmaskSVD -vnorm -mask $csfmask -sval $(( ${r/?(d)csf/} - 1 )) -polort 2 ${func_ts}${ext} > .csf_ts" || ( echo "3dmaskSVD failed. Unable to extract CSF voxels. Check your mask! $csfmask" && exit 1 )
	    else
		rel "3dmaskave -mask $csfmask -q ${func_ts}${ext} > .csf_ts" || ( echo "3dmaskave failed. Unable to extract CSF voxels. Check your mask! $csfmask" && exit 1 )
	    fi

	    rel "1d_tool.py -overwrite -infile .csf_ts -demean -write .csf_ts" #demean	    
	    rel "1d_tool.py -overwrite -infile .csf_ts -derivative -demean -write .csf_ts_deriv" #derivative
	fi

	if [[ "$r" =~ ^(wm|dwm|wm[0-9]+|dwm[0-9]+)$ && ! -f ".wm_ts" ]]; then
	    if [ $no_warp -eq 0 ]; then
		#if we have a warp to standard space then use the inverse warp to transform the wm mask from the probabilistic atlas to subject space
		if [ $templateName = MNI_3mm ]; then
		    #0.9 with erosion gives nice mask: 793 voxels
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_wm_tal_nlin_asym_09c_3mm\" -thr 0.9 -bin -eroF .template_wm_prob"
		elif [ $templateName = MNI_2.3mm ]; then
		    #0.95 with erosion: 2681 voxels
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_wm_tal_nlin_asym_09c_2.3mm\" -thr 0.95 -bin -eroF .template_wm_prob"
		elif [ $templateName = MNI_2mm ]; then		    
		    rel "fslmaths \"$stddir/mni_icbm152_nlin_asym_09c/mni_icbm152_wm_tal_nlin_asym_09c_2mm\" -thr 0.95 -bin -eroF .template_wm_prob"
		elif [ $templateName = SPM_2mm ]; then
		    #0.85 with erosion: 2981 voxels
		    rel "fslmaths \"$stddir/spm8_mni/white.nii\" -thr 0.85 -bin -eroF .template_wm_prob"
		elif [ $templateName = MNI_FSL_3mm ]; then
		    #0.8 with erosion: 659 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_white_3mm.nii\" -thr 0.8 -bin -eroF .template_wm_prob"
		elif [ $templateName = MNI_FSL_2.3mm ]; then
		    #0.85 with erosion: 1546 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_white_2.3mm.nii\" -thr 0.85 -bin -eroF .template_wm_prob"
		elif [ $templateName = MNI_FSL_2mm ]; then
		    #0.85 with erosion: 3113 voxels
		    rel "fslmaths \"$stddir/fsl_mni152/avg152T1_white_2mm.nii\" -thr 0.85 -bin -eroF .template_wm_prob"	    
		elif [ $templateName = 1YO_2mm ]; then
		    # integers 0 to 246. 75% is 184
		    rel "fslmaths \"$stddir/UNCInfant/2mm_1YO/infant-1yr-wm.nii.gz\" -thr 184 -bin -eroF .template_wm_prob"	    
		else
		    rel "Unsupported template brain for WM extraction: $templateName" c
		    exit 1
		fi

		local wmmask=.template_wm_prob${ext}
	    else

		#fast segmentation should have been run during preprocessMprage
		#use seg_2 for wm
		#re-run FAST here if needed
		[ $( imtest "${mprageBet_base}_fast_seg_2" ) -eq 0 ] && rel "fast -g -o \"${mprageBet_base}_fast\" \"${mprageBet_base}\""

		warp_structmask_to_func "${mprageBet_base}_fast_pve_2" .wm_pve

		#threshold at 0.8 to ensure that only voxels likely to be wm are retained	    
		rel "fslmaths .wm_pve -thr 0.8 -bin .wm_pve_thr0p8"

		#erode wm mask twice to reduce risk of partial WM voxels
		#more contiguous wm voxels than csf voxels -- hence 2x dilation works well in general as a start to retain deep cerebral white matter
		rel "3dmask_tool -overwrite -input .wm_pve_thr0p8${ext} -dilate_result -2 -prefix .wm_pve_thr0p8_ero2${ext}"

		#cluster wm mask and only retain clusters of 15 or more contiguous voxels
		rel "$FSLCLUSTER -i .wm_pve_thr0p8_ero2 -t 1.0 --osize=.wm_pve_clustsize --no_table"
		rel "fslmaths .wm_pve_clustsize -thr 15.01 -bin -mas ${postSS}_tmean_mask .wm_mask -odt char"
		
		#check that at least 50 voxels are present in mask to get a reasonable sampling distribution of WM timeseries
		local nvox=$( 3dBrickStat -non-zero -count .wm_mask${ext} )

		if [ $nvox -lt 50 ]; then
		    #go back to a 1x erosion, but up voxel clusters to 50 minimum
		    rel "Warning: fewer than 50 WM voxels were found in .wm_mask. Reverting to a 1x eroded WM mask" c
		    rel "3dmask_tool -overwrite -input .wm_pve_thr0p8${ext} -dilate_result -1 -prefix .wm_pve_thr0p8_ero1${ext}"
		    
		    rel "$FSLCLUSTER -i .wm_pve_thr0p8_ero1 -t 1.0 --osize=.wm_pve_clustsize_ero1 --no_table"
		    rel "fslmaths .wm_pve_clustsize_ero1 -thr 50.01 -bin -mas ${postSS}_tmean_mask .wm_mask_ero1 -odt char"

		    local wmmask=.wm_mask_ero1${ext}
		else
		    local wmmask=.wm_mask${ext}
		fi    

		#cleanup some intermediate files
		rel "imrm .wm_pve_thr0p8 .wm_pve_thr0p8_ero2 .wm_pve_thr0p8_ero1 .wm_pve_clustsize .wm_pve_clustsize_noero_ero1"
	    fi

	    #extract WM time series. If a number follows the WM specification, this is an aCompCor approach with the specified number of components
	    if [[ "$r" =~ ^(wm[0-9]+|wm[0-9]+)$ ]]; then
		rel "3dmaskSVD -vnorm -mask $wmmask -sval $(( ${r/?(d)wm/} - 1 )) -polort 2 ${func_ts}${ext} > .wm_ts" || ( echo "3dmaskSVD failed. Unable to extract WM voxels. Check your mask! $wmmask" && exit 1 )
	    else
		rel "3dmaskave -mask $wmmask -q ${func_ts}${ext} > .wm_ts" || ( echo "3dmaskave failed. Unable to extract WM voxels. Check your mask! $wmmask" && exit 1 )
	    fi

	    rel "1d_tool.py -overwrite -infile .wm_ts -demean -write .wm_ts" #demean	    
	    rel "1d_tool.py -overwrite -infile .wm_ts -derivative -demean -write .wm_ts_deriv" #compute derivative

	fi

	# NOTE: for now, I believe these will blow up if specified before the basic regressor (e.g., wm)
	# quadratic of csf (qcsf*)
	if [[ "$r" =~ ^(qcsf|qcsf[0-9]+)$ && ! -f ".csf_ts_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.csf_ts')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.csf_ts_quad')\""
	fi

	# quadratic of dcsf (qdcsf*)
	if [[ "$r" =~ ^(qdcsf|qdcsf[0-9]+)$ && ! -f ".csf_ts_deriv_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.csf_ts_deriv')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.csf_ts_deriv_quad')\""
	fi
	
	# quadratic of wm (qwm*)
	if [[ "$r" =~ ^(qwm|qwm[0-9]+)$ && ! -f ".wm_ts_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.wm_ts')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.wm_ts_quad')\""
	fi

	# quadratic of dwm (qdwm*)
	if [[ "$r" =~ ^(qdwm|qdwm[0-9]+)$ && ! -f ".wm_ts_deriv_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.wm_ts_deriv')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.wm_ts_deriv_quad')\""
	fi
	    
	#global signal
	if [[ "$r" =~ ^(gs|dgs)$ && ! -f ".gs_ts" ]]; then
	    #erode the tight skull-stripped brain mask from the structural scan to identify brain voxels
	    #rel "fslmaths ${mprageBet_base} -bin -eroF -eroF .brainmask_ero2x_anat -odt char" #eroding the binarized mprage_bet

	    #warp_structmask_to_func .brainmask_ero2x_anat .brainmask_ero2x

	    #rel "fslmaths .brainmask_ero2x -mas ${postSS}_tmean_mask .brainmask_ero2x -odt char" #mask by functional brain mask

	    #MH May2017: On further consideration, it's probably best to simply erode the mask from betting the functionals since
	    #this does not depend on the quality of the functional -> structural coregistration (esp. concerning in high distortion areas)
	    rel "fslmaths ${postSS}_tmean_mask -eroF -eroF .brainmask_ero2x -odt char"
	    
	    #GS uses pre-warp data since we are using a brain mask in functional space
	    rel "3dmaskave -mask .brainmask_ero2x${ext} -q ${prewarp_func_ts}${ext} > .gs_ts"
	    rel "1d_tool.py -overwrite -infile .gs_ts -demean -write .gs_ts" #demean	    
	    rel "1d_tool.py -overwrite -infile .gs_ts -derivative -demean -write .gs_ts_deriv"
	fi

	# quadratic of gs
	if [[ "$r" == "qgs" && ! -f ".gs_ts_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.gs_ts')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.gs_ts_quad')\""
	fi

	# quadratic of dgs
	if [[ "$r" == "qdgs" && ! -f ".gs_ts_deriv_quad" ]]; then
	    rel "Rscript -e \"x <- read.table('.gs_ts_deriv')$V1; x2 <- x^2; cat(x2 - mean(x2), sep='\n', file='.gs_ts_deriv_quad')\""
	fi
	
	case "$r" in 
	    6motion)
		nuisancecols+=(".motion_demean[0]")
		nuisancecols+=(".motion_demean[1]")
		nuisancecols+=(".motion_demean[2]")
		nuisancecols+=(".motion_demean[3]")
		nuisancecols+=(".motion_demean[4]")
		nuisancecols+=(".motion_demean[5]")
		;;
	    rx)
		nuisancecols+=(".motion_demean[0]");;
	    ry)
		nuisancecols+=(".motion_demean[1]");;
	    rz)
		nuisancecols+=(".motion_demean[2]");;
	    tx)
		nuisancecols+=(".motion_demean[3]");;
	    ty)
		nuisancecols+=(".motion_demean[4]");;
	    tz)
		nuisancecols+=(".motion_demean[5]");;
	    q6motion)
	      nuisancecols+=(".motion_quad[0]")
	      nuisancecols+=(".motion_quad[1]")
	      nuisancecols+=(".motion_quad[2]")
	      nuisancecols+=(".motion_quad[3]")
	      nuisancecols+=(".motion_quad[4]")
	      nuisancecols+=(".motion_quad[5]")
	      ;;
	    qrx)
	      nuisancecols+=(".motion_quad[0]");;
	    qry)
	      nuisancecols+=(".motion_quad[1]");;
	    qrz)
	      nuisancecols+=(".motion_quad[2]");;
	    qtx)
	      nuisancecols+=(".motion_quad[3]");;
	    qty)
	      nuisancecols+=(".motion_quad[4]");;
	    qtz)
	      nuisancecols+=(".motion_quad[5]");;
	    qd6motion)
		nuisancecols+=(".motion_deriv_quad[0]")
		nuisancecols+=(".motion_deriv_quad[1]")
		nuisancecols+=(".motion_deriv_quad[2]")
		nuisancecols+=(".motion_deriv_quad[3]")
		nuisancecols+=(".motion_deriv_quad[4]")
		nuisancecols+=(".motion_deriv_quad[5]")
		;;
	    qdrx)
		nuisancecols+=(".motion_deriv_quad[0]");;
	    qdry)
		nuisancecols+=(".motion_deriv_quad[1]");;
	    qdrz)
		nuisancecols+=(".motion_deriv_quad[2]");;
	    qdtx)
		nuisancecols+=(".motion_deriv_quad[3]");;
	    qdty)
		nuisancecols+=(".motion_deriv_quad[4]");;
	    qdtz)
		nuisancecols+=(".motion_deriv_quad[5]");;
	    d6motion)
		nuisancecols+=(".motion_deriv[0]")
		nuisancecols+=(".motion_deriv[1]")
		nuisancecols+=(".motion_deriv[2]")
		nuisancecols+=(".motion_deriv[3]")
		nuisancecols+=(".motion_deriv[4]")
		nuisancecols+=(".motion_deriv[5]")
		;;
	    drx)
		nuisancecols+=(".motion_deriv[0]");;
	    dry)
		nuisancecols+=(".motion_deriv[1]");;
	    drz)
		nuisancecols+=(".motion_deriv[2]");;
	    dtx)
		nuisancecols+=(".motion_deriv[3]");;
	    dty)
		nuisancecols+=(".motion_deriv[4]");;
	    dtz)
		nuisancecols+=(".motion_deriv[5]");;
	    wm*)
		nuisancecols+=(".wm_ts");;
	    qwm*)
		nuisancecols+=(".wm_ts_quad");;
	    dwm*)
		nuisancecols+=(".wm_ts_deriv");;
	    qdwm*)
		nuisancecols+=(".wm_ts_deriv_quad");;
	    csf*)
		nuisancecols+=(".csf_ts");;
	    qcsf*)
		nuisancecols+=(".csf_ts_quad");;
	    dcsf*)
		nuisancecols+=(".csf_ts_deriv");;
	    qdcsf*)
		nuisancecols+=(".csf_ts_deriv_quad");;		    
	    gs)
		nuisancecols+=(".gs_ts");;
	    qgs)
		nuisancecols+=(".gs_ts_quad");;
	    dgs)
		nuisancecols+=(".gs_ts_deriv");;
	    qdgs)
		nuisancecols+=(".gs_ts_deriv_quad");;
	esac

    done

    rel "Creating ${nuisance_file} containing regressors ($(echo $reg|tr '\n' ' ')) as columns: ${nuisancecols[*]}" c
    #rel "1dcat -overwrite ${nuisancecols[@]} > ${nuisance_file}"
    # 20151218 WF -- "star" instead of "at" so single shell-world (one argument) sent to rel
    rel "1dcat ${nuisancecols[*]} > unfiltered_${nuisance_file}"

    #MH Jun2017: Need to apply the same filtering to nuisance regressors as the data
    #3dBandpass will take care of this internally, but 3dDetrend, or the user later at some point,
    #may not, which could cause frequency mismatch

    if [ "$bandpass_filter" -eq 1 ]; then
	rel "1dBandpass -dt $tr $bpLow $bpHigh unfiltered_${nuisance_file} > ${nuisance_file}"
    elif [ "$no_hp" -eq 0 ]; then
	#hpFilter has already been converted in parse_args to FWHM volumes. Pass forward verbatim to 1dbptf"
	rel "1dbptf -matrix unfiltered_${nuisance_file} -tr $tr -time_along_rows -out_file ${nuisance_file} -hp_volumes $hpFilter"
    fi    
    
    # TODO: do multiple passes to build nuisancecols so we can check that would should be redoing?
    rel "echo '${nuisancecols[*]}' > .regressors_in_use${custom_regprefix:-}"

    rel "date > .nuisance_compute_complete${custom_regprefix:-}"

    return 0
}

# look into regression options and build todo string
# A b g (r || custom_prefix)
function regression_todo {
    local todo=""
    [ $rmautocorr -eq 1 ] && todo="A"
    [ $bandpass_filter -eq 1 ] && todo="${todo}b"
    [ $gsr_in_prefix       -eq 1 ] && todo="${todo}g"
    [ $nuisance_regression -eq 1 ] && todo="${todo}${custom_regprefix:-r}"
    echo "$todo"
}
function reg_already_done {
    local todo="${1:?$FUNC_NAME needs regression operations by prefix, eg Abr}"
    local reg_prefix=${custom_regprefix:-r}
    local complete_file=.nuisance_regression${custom_regprefix:-}_complete # default .nuisance_regression_complete

    [[ -f $complete_file && $todo=$reg_prefix ]] && return 0
    [[ -f .bandpass_filter_complete && $todo=b ]] && return 0
    [[ -f .bandpass_filter_complete && .remove_autocorr_complete && $todo=Ab ]] && return 0
    [[ -f $complete_file && -f .bandpass_filter_complete && $todo=b$reg_prefix ]] && return 0
    [[ -f $complete_file && -f .remove_autocorr_complete && $todo=A$reg_prefix ]] && return 0
    [[ -f $complete_file && -f .bandpass_filter_complete && -f .remove_autocorr_complete && $todo=Ab$reg_prefix ]] && return 0

    # not already done
    return 1
}
function reg_todo_to_desc {
    [[ $# -eq 0 || -z $1 ]] && echo "none" && return

    local todo="${1:?$FUNC_NAME first arg is regression operations by prefix/todo. eg. Abgr }"
    #bandpass + prewhitening only
    if [ "$todo" = "Ab" ]; then
       echo "bandpass+prewhite"
    #bandpass only
    elif [[ "$todo" = "b" ]]; then
       echo "bandpass"
    # regress only (withour or without 'A'utocorrelation rm'ed (w/reml)
    elif [[ $todo =~ "A?${custom_regprefix:-r}" ]]; then
       echo "regressonly"
    elif [[ "$todo" =~ A?bg?${custom_regprefix:-r} ]]; then
       echo "bandpass+regress"
    else
       echo "none"
    fi
}

function nuisance_regression {

    #handle nuisance regression and/or bandpass filtering
    if [[ \
         "$nuisance_compute"    -eq 0 &&  \
         "$nuisance_regression" -eq 0 &&  \
         "$bandpass_filter"     -eq 0 &&  \
         "$rmautocorr" -eq 0     \
    ]]; then
      #rel "(Not trying to deal with nuisance regressors, bandpassing, or autocorrelation)" c
      return 0 #no nuisance, bandpass, or basis correct. exit
    fi

    # 20230111 - wish had these for regenerating single missing in 7T mgsencmem
    print_vars_in COMPUTE_NUISANCE_REGRESSORS_GLOBALS || :

    preNRBP="${prefix}${funcFile}${smoothing_suffix}${ext}"

    todo=$(regression_todo)
    # todo is 
    #  only: A b r 
    #  combined: Ab Abr br
    #  'A' alone is not allowed, will throw error -- checked by parseargs or preprocessFunctional
    #  gsr_in_prefix (-gsr) must have gs in $nuisance_regressors, check by parseargs
    
    #check whether requisite steps are complete
    reg_already_done "$todo" && return 0
    
    prefix="${todo}${prefix}"
    postNRBP="${prefix}${funcFile}${smoothing_suffix}${ext}"

    if [[ ! -f "${nuisance_file}" && ( $nuisance_regression == 1 || $nuisance_compute == 1 ) ]]; then
	compute_nuisance_regressors
    fi

    #bandpass + prewhitening only
    case $(reg_todo_to_desc "$todo") in
       bandpass+prewhite)
          run_3dremlfit $preNRBP $postNRBP ${subjMask}${ext} #function should pick up this circumstance
          rel "date > .bandpass_filter_complete"
          ;;
       bandpass)
          rel "3dBandpass -input \"$preNRBP\" -mask \"${subjMask}${ext}\" -dt $tr \
                 -prefix \"$postNRBP\" $bpLow $bpHigh"
          rel "date > .bandpass_filter_complete"	
          ;;
       regressonly)
          if [ $arma_nuisance_regression -eq 0 ]; then
              rel "3dDetrend -overwrite -verb -polort 2 -vector ${nuisance_file} \
                 -prefix \"$postNRBP\" \"$preNRBP\""
              #need to mask detrend by brain mask as with 3dBandpass (3dDetrend doesn't support -mask)
              rel "fslmaths \"$postNRBP\" -mas \"$subjMask\" \"$postNRBP\""
          else
              #use 3dREMLfit for regression (default)
              run_3dremlfit $preNRBP $postNRBP ${subjMask}${ext}	    
          fi
          rel "date > .nuisance_regression${custom_regprefix:-}_complete"
          ;;
       bandpass+regress) 
	   if [ $arma_nuisance_regression -eq 0 ]; then
	       rel "3dTproject -overwrite -input \"$preNRBP\" -mask \"${subjMask}${ext}\" -dt $tr \
	       	-prefix \"$postNRBP\" -polort 2 -ort ${nuisance_file} -passband $bpLow $bpHigh"

	       #Deprecated -- 3dTproject performs true simultaneous approach, not filter all -> bandpass
               #rel "3dBandpass -overwrite -input \"$preNRBP\" -mask \"${subjMask}${ext}\" -dt $tr \
	       #     	-prefix \"$postNRBP\" -ort ${nuisance_file} $bpLow $bpHigh"
	   else
	       #use 3dREMLfit for bandpass + regression (default)
	       run_3dremlfit $preNRBP $postNRBP ${subjMask}${ext}
	   fi
	   rel "date > .bandpass_filter_complete"
	   rel "date > .nuisance_regression${custom_regprefix:-}_complete"
          ;;
       none) rel "(Regressors computed but no bandpass or nuisance regression applied)" c ;;
       *)
          rel "ERROR: don't know what to do with $todo = '$(reg_todo_to_desc "$todo")'" c
          return 1
        ;;
    esac
       
    
    #explicit return code needed to avoid implicit status of prior command
    return 0
    
}

# vim: set tabstop=7: