preprocessFunctional

#!/bin/bash

# wrap bad exit in info message
originalargs="$@"
originalcwd="$(pwd)"
#trap '[ $? -ne 0 ] && echo -e "[$(date +%F\ %H:%M) $originalcwd] ERROR exiting $0\n\t$originalargs"' EXIT
trap 'traperror $? $LINENO $BASH_LINENO "$BASH_COMMAND" $(printf "::%s" ${FUNCNAME[@]:-})' ERR

function printHelp() {

cat <<EndOfHelp
-----------------------------------
preprocessFunctional is a script to preprocess an fMRI EPI scan.

Three files are required for this preprocessing script:
  1) 4d functionals file (-4d command line parameter). Alternative: -dicom
  2) betted mprage scan (-mprage_bet command line parameter)
  3) mprage-to-standard-space warp coefficients (-warpcoef command line parameter)

This script depends on the following tools:
  1) FSL 5.0 or higher
  2) AFNI (preferably 2012 or later)
  3) imagemagick (for -mc_movie)
  4) ffmpeg (for -mc_movie)
  5) python 3.7 or later (for -4d_slice_motion)
        - nipy development version: https://github.com/nipy/nipy/
        - nibabel
        - scipy
  6) awk
  7) perl

Recommended:
  1) ANTS 2.2 or higher (used for epiref -> func registration at the moment)
  2) c3d 1.1 or higher (used for transforming ITK matrices from ANTS)

Command line options:
  -4d: The 4d file containing unprocessed functional volumes for a single subject. Required.
  -4d_slice_motion: use 4d simultaneous slice-timing motion correction algorithm (calls sliceMotion4D)
  -bandpass_filter: Apply bandpass filter to fMRI data using 3dTproject. Specify low and high cutoff frequencies in Hz.
        Example: -bandpass_filter 0.009 .08
  -bids_sidecar: load provided json file. will autocheck for -4d.json. forces 4d_slice_motion unless -no_st.
  -bet_frac: fractional intensity threshold for skull-stripping mean functional (used in func-to-struc warp). Default 0.3.
  -cite: list citation for tools used by preprocessFunctional and preprocessMprage
  -check_dependencies: what tools are avaible on this system. Exit w/ status 0 if everything we could use is found
  -cleanup: Deletes all files used in preprocessing that are unlikely to be needed for analysis, see -keep_motion_files 
  -cleanup_only: Delete all intermediate files, but do not run any processing steps.
  -compute_warp_only: Do not warp data to template space, but do compute the nonlinear transformation to template (and FM unwarp).
        This is useful if one wishes to analyze the data in native space and warp the stats maps to template. This switch also enables -no_warp.
  -constrain_to_template: Constrain brain voxels to the corresponding standard template (e.g., MNI). y or n. Default: y.
  -custom_slice_times: For use with 4d slice timing + motion correction. 
        A csv file containing the times (in seconds summing to TR) of each slice (bottom to top).
        Example: -custom_slice_times "../SPECCMBSliceTimes.csv"
  -delete_dicom: if converting dicom to nifti, whether delete or archive DICOM files. If not
        passed, user will be prompted for action. Options are:
        -delete_dicom no. (leaves DICOM untouched)
        -delete_dicom yes. (deletes DICOM files)
        -delete_dicom archive. (compresses files into archive file): functional_dicom.tar.gz.
  -deoblique_all: Deoblique all functional datasets (3dWarp -deoblique) after all preprocessing complete.
  -despike: Interpolate intensity spikes using 3dDespike. Occurs after motion correction and slice timing correction.
  -despike_thresh: Cut thresholds for 3dDespike (how much to despike). Default to 2.5 and 4.0.
  -dicom: The file pattern to match all dicom files to be converted to nifti. Must be in quotes.
        Example: -dicom "MR*" (This parameter is an alternative to passing the 4d nifti file using -4d)
  -epi_echospacing: effective echo spacing of the EPI data in seconds used in TOPUP distortion correction.
        See SIEMENS Sequence > Echo spacing. Beward partial Fourier and iPAT/GRAPPA.
  -epi_pedir: phase-encoding direction of the EPI data. Should be x/y/z/x-/y-/z-, following FSL FUGUE conventions.
        This is used for TOPUP distortion correction of functional data. For reference, A>>P is -y and P>>A is y.
  -epi_te: echo time of EPI data in milliseconds. Auto-detected in DICOMs. Only needed for fieldmap correction if an -fm_cfg file is not used.
  -flip_topup_pedir: sometimes, the field generated by TOPUP is accurate, but opposite of the EPI phase encoding direction.
        In this case, the files in topup_undistort will look right, but functional -> structural and functional -> template
        files will look doubly distorted. To fix this, pass in this flag, which will flip the expected PE dir for all steps downstream
        that use the TOPUP field. See the 'N.B.' here for details: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide
  -fmap_struct_dof: The number of degrees of freedom to be used when co-registering the fieldmap magnitude image (FM_UD_fmap_mag_brain)
        to the subject structural (T1) scan. Defaults to flirt bbr (boundary based registration) following a 6 dof transformation.
        -fmap_struct_dof bbr. BBR segments the T1 to create a white matter boundary mask, which is used to help align EPI and T1.
        -fmap_struct_dof 6.  Rotation x, y, z; Translation x, y, z.
  -fm_phase: gre field map phase image
  -fm_magnitude: field map magnitude image (mag2.nii.gz)
  -fm_cfg: configuration file for handling fieldmap (TEs, echo spacing, etc.)
  -func_refimg: T2* image (as NIfTI) used as intermediate registration target (i.e., func -> ref -> structural -> template).
        This is useful when the T2* images have poor GM/WM contrast (such as 1.0s MB 5x sequence) and BBR is doing poorly for some subjects.
        Note: this can be a single DICOM image, which will be converted to NIfTI.
        Note: this can also be a directory, in which case the script assumes it is the subject's raw directory and will search for an SBref 
	  image corresponding to the protocol name of the functional data and is either the previous or subsequent series number.
  -func_struct_dof: The number of degrees of freedom to be used when co-registering the mean subject
        functional (EPI) scan to the subject structural (T1) scan. Defaults to bbr (boundary based registration).
        -func_struct_dof bbr. BBR segments the T1 to create a white matter boundary mask, which is used to help align EPI and T1.
        -func_struct_dof 6.  Rotation x, y, z; Translation x, y, z.
        -func_struct_dof 7.  Rotation x, y, z; Translation x, y, z; Global scaling.
        -func_struct_dof 12. Rotation x, y, z; Translation x, y, z; Scaling x, y, z; Skewing x, y, z.
  -gsr: specifically label gsr as a prefix. r becomes gr in filenames. "gs" must be in nuisance_regression/compute
  -keep_motion_files: when doing -cleanup or -cleanup_only, keep mt_*,tm_*, some mc_*, and mc_mats/ files
  -help: Print command help.
  -hp_filter: high-pass filtering frequency (FWHM of filter in volumes). Anything faster than this value will be passed. Default is 40.
        It is preferable to specify this in seconds using suffix "s", such as -hp_filter 100s.
        This allows the script to compute the FWHM of the filter, which = filter_cutoff[seconds]/(sqrt(8*log(2))*TR[seconds]).
  -ica_aroma: automatically remove motion artifacts using ICA-AROMA method of Pruim 2015 NeuroImage. Applied after smoothing
  -ignore_bids: do not read TR or SliceTiming from BIDS .json sidecar.
  -log: Name for log file that documents each command that is run. Default: preprocessFunctional.log
  -mc_first: run motion correction before slice timing correction, regardless of slice acquisition order.
  -mc_movie: create .mpg movies of head movement before and after motion correction (depends on create4DMovie).
  -mc_program: what program to use for motion correction. Options: mcflirt, 3dvolreg. Default: mcflirt.
  -motion_censor: A comma-delimited list of motion/intensity metrics for identifying volumes to censor.
        This calls fsl_motion_outliers for each element in the list and generates both spike regressor matrices (.mat)
        and AFNI-style censor 1D files. Available metrics are dvars, fd, fdrms, relrms, refrms, refmse. If an equals sign is included,
        then this is used as the threshold for flagging volumes to exclude. Otherwise,75th pctile + 1.5*IQR is used.
        Example: -motion_censor fd=0.9,dvars=20 (censor FD at 0.9mm or greater, and DVARS at 20)
  -motion_sinc: (mcflirt only) Complete a 4th stage optimization of motion correction using sinc interpolation (slower but better): y or n. Default y.
  -mprage_bet: The betted structural scan for this participant. Required.
  -no_hp: Do not perform high-pass filtering.
  -no_mc: Do not preform motion correction. (e.g. running after HCP minimum; added 20190916)
  -no_smooth: Do not perform spatial smoothing.
  -no_st: Do not perform slice timing correction.
  -no_warp: Maintain native space and resolution of input image. Do not warp to a template.
  -nuisance_file: Name of file containing all nuisance regressors specified by -nuisance_regression. Default .nuisance_regressors.
  -nuisance_compute: same as -nuisance_regression, but only compute the comma-delimited list of nuisance regressors. This does NOT regress any signals out of the data,
        just creates the file specified by -nuisance_file, which can be used later (e.g., fMRI GLM) to remove nuisance signals.
  -nuisance_regression: regress out signals that are not of interest (most relevant to resting-state processing).
        Comma-delimited list (no spaces!) that can include: 6motion, rx, ry, rz, tx, ty, tz, wm (white matter), csf (cerebrospinal fluid), 
        gs (global signal, use -gsr to label in filename). Derivatives of any of these can be added by using the prefix "d" (e.g., dcsf).
        Motion regressors can also be prefixed with q for quadratic and qd for deriv of quad (e.g. qd6motion)
        If bandpass filtering is specified, nuisance regression is done by 3dTproject (Hallquist et al. 2013, NeuroImage). Otherwise, use 3dDetrend.
        If a number follows csf or wm (e.g., csf5), then an anatomical CompCor approach is taken with the desired number of principal components
	   computed using AFNI 3dmaskSVD.
        Example: -nuisance_regression 6motion,csf,wm,dcsf,dwm
  -output_basename: Base name of file output during DICOM -> NIfTI conversion. Default: functional
  -partial_proc: do not run full pipeline. options are 'motion', 'despike', 'warp'. kills script, and marked as incomplete after desired step (20190729WF)
  -physio_resp: respitory file (1D or .resp [siemens])
  -physio_card: pulse file (1D or .puls [siemens])
  -physio_func_info: dicom directory/bidsjson file, used for timing info to extract resp+puls files
  -prefer_ants: Use ANTS 6dof coregistration when possible (e.g., epiref -> func and fmap_mag -> struct)
  -prefer_flirt: Use flirt 6dof coregistration when possible (e.g., epiref -> func and fmap_mag -> struct)
  -ref_vol: Reference volume to be used in motion correction. Will default to mean volume if not passed.should be a volume (time step) number or either mean or median
  -rescaling_method: The method used to normalize voxel intensities. Options are:
        -rescaling_method 10000_globalmedian (FSL's default, rescale the 4d file by a single quantity: 10000/[global median])
        -rescaling_method 100_voxelmean (Rescale each voxel by 100/[voxel 4d mean]. In principle, approximates % change)
  -resume: Instructs script to resume preprocessing from the last completed step. Should be passed in as the only parameter: preprocessFunctional -resume.
  -rmautocorr: use 3dREMLfit to remove nuisance regressors and filtered frequencies (1dBport) using an ARMA(1,1) model to account for autocorrelation.
  -rmautocorr: use 3dREMLfit to remove nuisance regressors and filtered frequencies (1dBport) using an ARMA(1,1) model to account for autocorrelation (-Rwherr).
  -rmgroup_component: provide a group 1D ts to remove from 4D ts after warping (add c to prefix after warp)
  -se_phasepos: The spin-echo image (or directory) with positive phase encoding direction (typically P>>A or R>>L) to be used in TOPUP distortion correction
  -se_phaseneg: The spin-echo image (or directory) with negative phase encoding direction (typically A>>P or L>>R) to be used in TOPUP distortion correction
  -slice_acquisition: The order in which slices were acquired.
        -slice_acquisition seqasc. Sequential ascending (bottom-up)
        -slice_acquisition seqdesc. Sequential descending (top-down)
        -slice_acquisition interleaved. Interleaved 1,3,5... 2,4,6...
  -smoothing_kernel: The size of gaussian smoothing kernel in mm. Default: 5. 
        -smoothing_kernel 'default' will use 5 if template res >=2.3., and 4 if 2mm.
  -smoother: The smoothing approach to use. Defaults to "susan" (tries to normalize intensities without
        blurring). Alternative is "gaussian", which is regular smoothing.
  -startover: rather than skipping preprocessing stages that have been
        completed, re-run all processing stages, overwriting any existing files.
  -st_first: run slice timing correction first, regardless of slice acquisition order.
        By default, slice timing is run before motion correction for interleaved sequences, but after
        motion correction for sequential acquisition.
  -template_brain: Which standard template brain to use. Default: MNI_3mm. Must match template from mprage normalization. Options:
        -template_brain MNI_2mm (MNI152 brain, 2mm isotropic voxels, template from MNI Fonov et al. 2009). This is the default in preprocessMprage
	-template_brain MNI_2.3mm (MNI152 brain, 2.3mm isotropic voxels, template from MNI Fonov et al. 2009). For 2.3mm MB sequence.
        -template_brain MNI_3mm (MNI152 brain, 3mm isotropic voxels, template from MNI Fonov et al. 2009)
        -template_brain MNI_FSL_3mm (MNI152 brain, 3mm isotropic voxels, template from FSL standard)
        -template_brain MNI_FSL_2.3mm (MNI152 brain, 2.3mm isotropic voxels, template from FSL standard)
        -template_brain MNI_FSL_2mm (MNI152 brain, 2mm isotropic voxels, template from FSL standard)
        -template_brain SPM_2mm (MNI brain from spm8 canonical template)
        -template_brain Tal_3mm (Talairach brain, 3mm isotropic voxels. Brain built by resampling MNI template in Tal space)
                        1YO_2mm, neo_2mm, 1YO_3mm, neo_3mm - UNC baby templates
  -threshold: The method to use for thresholding/masking low intensity voxels. Options are:
        -threshold 98_2 (FSL's default, 2nd %ile + (98th %ile - 2nd %ile)/10
        -threshold 10 (mask any voxels below the 10th percentile)
  -tr: TR duration in seconds.
  -trunc: remove the specified number of volumes (time steps) from the beginning. Some protocols do not discard enough automatically. ( n_rm_firstvols )
  -siemens: only relevant for -4d input. If set, slice timing correction will check for interleaved
        order and odd vs. even number of slices. If even, will create and use --ocustom file 2,4,6...1,3,5...
        This is necessary because of different slice acquisition on Siemens for even versus odd total slices.
  -unifize_funcref: run 3dUnifize on ref before 3dskullstrip. added 20210513
  -verbose: print whats going on (set -xv)
  -warpcoef: The warp coefficients used to transform the structural scan to standard space. Required.
  -warp_interpolation: The intepolation method used to warp functional data to the mprage and to the standard space (e.g., MNI). Default: spline.
        -warp_interpolation spline. Use cubic spline interpolation to apply warp coefficients. Faster than sinc and slightly more blurry, but more robust to mask. (Default)
        -warp_interpolation sinc. Provides higher quality non-blurry resampling relative to trilinear, but is slower. N.B. Only works well when the subject's
               functional mask is precise and contains only brain voxels. Otherwise, sinc is prone to stretching the warp strangely.
        -warp_interpolation trilinear. Faster method, but somewhat blurry.
  -wavelet_despike: Apply wavelet despiking after motion correction and slice timing correction to remove frequency-dependent nonstationary events in fMRI time series (Patel 2014 NeuroImage).
        Note: This is an alternative to 3dDespike and should generally be preferred because it removes intensity artifacts in frequency space, which is useful
        because head movement often produces artifacts that are spread in time (spin history).
  -wavelet_threshold: This is the coefficient threshold value used during maxima and minima coefficient definition. Default: 10. See BWT_Documentation.pdf for details.

ENVIRONMENT
  to control e.g. python:numpy (sliceMotion4d) set environment. defaults to 1

  export OPENBLAS_NUM_THREADS=5
  # or
  export MKL_NUM_THREADS=5

The major steps in the processing pipeline are, in order:
  0)    fslorient (_) and optionally truncate (0)
  1)    Physio/retroicor (optional) - p
  2)    Slice timing correction - t
  3)    Motion correction - m
  4)    Skull stripping and intensity thresholding - k
  5)    Intensity/wavelet despiking (optional) - d
  6)    Field map unwarping (optional) - u
  7)    Co-registration and warping to standard space - w
  8)    Remove group level components (optional) - c  [requires warp]
  9)    Smoothing - s
  10)   ICA-AROMA (optional) - a
  11)   High-pass filtering - f
  12)   Intensity normalization - n
  13)   Nuisance regression (optional) - r
  13.5) Global signal regression (optional, label optional) - g
  14)   Bandpass filtering (optional) - b
  15)  Auto Correlation removal (optional) - A

Note that if the slice acquisition order is sequential, then the default is to appply motion correction before slice timing correction, whereas for
interleaved acquisition, the default is slice timing first. The order of these steps can be explicitly specified using -st_first, -mc_first, or -no_st.
Alternatively, consider -4d_slice_motion, which performs both steps at once and is likely more effective than conventional methods.

Characters are prepended as a prefix to the 4d filename, such that the final processed file will have a
prefix such as nfswkmtd, indicating the sequence of processing steps, from right to left.

N.B. For the warping to work properly, you must use preprocessMprage for the mprage scan, which uses the same template set.

Example call:
  preprocessFunctional -4d 10802func.nii.gz -mprage_bet 10802mprage.nii.gz -warpcoef mprage_warpcoef.nii.gz -slice_acquisition seqasc \
      -tr 2 -smoothing_kernel 6 -cleanup

-----------------------------------
EndOfHelp

}

#Written: 2010/5/18
#Changelog:
#
#Last updated: 04/30/2024 - add ignore_bids, slicetime missing message
#Last defaults change: 2018/7/26 - 3dREMLfit is default 'r'
#
#Changelog:
#
#04/30/2024 (WF)
#  - added '-ignore_bids' option and handle slicetiming missing errors
#08/26/2022 (WF)
#  - use no_mc + no_st + no_warp for fmriprep output of task data
#    preprocessFunctional -no_mc -no_st -no_warp -hp_filter 40 -rescaling_method 10000_globalmedian -4d sub*task-*space-MNI*_bold.nii.gz
#05/19/2021 (WF)
#  - added '-smoothing_kernel "default"': kernel size based on template res
#05/13/2021 (WF)
#  - added '-unifize_funcref' 
#01/25/2021 (WF)
#  - if not set, use OPENBLAS_NUM_THREADS=1 / MKL_NUM_THREADS=1
#    set in correct_slice_motion_4d
#    export before running preprocessFunctional to use different settings
#09/16/2019 (WF) 
#  - added -no_mc option for running ontop of another pipeline
#07/29/2019 (WF) 
# - add partial_proc to stop in error at motion or warp 
#   b/c we have >2000 NCANDA scans to preprocess but are not sure what fieldmap should be
#05/02/2019 (WF)
# - merge DEPENd and LNCD
#  $autocorr_with_basis=>$rmautocorr
#  $funcStructFlirtDOF=>$func_struct_dof
#  DEPENd: reml, warp options, psu matlab module | LNCD: preprocessDistortion, explict global signal flag, pyhsio, use_old_mni, rmgroup_component
#11/01/2018 (WF)
# - add old/new template checking in template_funcs (-use_old_mni)
#7/30/2018 (MH)
# - The use of BBR coregistration for fmap -> struct only works when the fieldmap has good WM contrast. Otherwise, it can backfire and yield
#    a very poor result. Thus, expose -fmap_struct_dof to give user control over how to coregister the images.
# - Prefer ANTS rigid registration in general (incl. fmap -> struct), but give user control using -prefer_flirt and -prefer_ants.
# - Use ANTS rigid registration or flirt skull-stripped followed by whole head (-nosearch) for fmap -> struct (this is prior to BBR if requested).
# - Move all qa_image commands to a separate log file to avoid junking up preprocessFunctional.log
#7/26/2018 (MH)
# - Switch to 3dREMLfit as the default for conducting nuisance regression (potentially including bandpass). See Bright & Murphy 2017 for details.
#    If you further specify -rmautocorr, you will receive the ARMA(1,1) whitened residuals after nuisance regression. I wouldn't necessarily recommend this,
#    but it does move toward the econometrics/statistics view that cross-correlations are only valid (non-spurious) if you have white time series.
#    The better/true approach to prewhitening in order to compute region-to-region correlations is implemented in ROI_TempCorr.R at the ROI level.
#7/18/2018 (MH)
# - Support FSL MNI 2.3mm template (also changed the templates so that tissue priors (CSF etc.) are exactly in register.
#7/1/2018 (MH)
# - Implement fsl_regfilt approach (partial regression) manually when ICA spits out so many components that fsl_regfilt crashes. See fslRegfilt.R.
# - Fixes to -rmautocorr to return pre-whitened residuals and not enforce bandpass
#5/11/2018 (WF)
# - add -gsr option to expliclty label gsr in regression step. prefix "rnswktm" becomes "grnswktm"
# - add -rmgroup_component (prefix c) to remove group components after warping
#3/12/2018 (MH)
# - Switch to 3dTproject over 3dBandpass for nuisance regression and bandpass filtering. Implements true simultaneous approach and censoring.
#2/12/2018 (MH)
# - Implement ICA AROMA + no smooth option. Because AROMA benefits from smoothed data to increase SNR, the data are smoothed for AROMA, but then
#    motion components are applied to unsmoothed data, and the pipeline proceeds from there (unsmoothed). Occurs -no_smooth and -ica_aroma are requested.
#7/5/2017 (MH)
# - support TOPUP distortion correction using S-E images (-se_phasepos, -se_phaseneg, -epi_echospacing, -epi_pedir)
# - Switch to 3dSkullStrip for epiref and mc_target due to weak performance of bet in several datasets
# - Produce many labeled overlay images in qa_images directory for checking coregistration, skull-stripping, and unwarping
# - Use BBR coregistration for fmap (magnitude) -> struct alignment. Tests suggest substantial improvement in some subjects.
# - Use ANTS, if available, for func -> epiref registration (excellent performance in tests). Fall back to flirt with skull-strip, then whole head.
# - Move all (linear) transform .mat files to transforms directory to avoid folder bloat
#6/5/2017 (MH)
# - More robust corrections for flawed fieldmap cache (reset to DICOMs if needed)
# - Fix -bptf to be nvol/2.355*TR to achieve FWHM (prior filter was ~2x too slow)
# - Always filter nuisance regressors to match imaging data (even just for compute or regress only)
# - Default to 120s hp filter if not specified (rather than volume-based default, which is TR-dependent)
#5/27/2017 (MH)
# - Support -func_refimg <directory> to have script look for DICOM with the same protocol name, but single-band reference (CMRR), within subject's raw dir
# - Support -ica_aroma to execute ICA-based Automatic Removal of Motion Artifact algorithm (Pruim 2015).
# - Support csf<ncomponents> and wm<ncomponents> syntax in -nuisance_regression to support aCompCor approach by extracting principal components
#4/21/2017 (WF)
# - add auto correlation removal via MJ+MH email exchange: -rmautocorr, prefix="A"
#2/15/2017 (WF)
# - motion nuis. reg.: add quad and quad of deriv 
# - add -trunc option to remove the first n volumes (PNC dataset)
# - ref_vol can be 'median' 
# - add bats tests for the 3 additions
#5/10/2016 (MH)
# - Support -func_refimg as DICOM for CMRR single-band images
# - Prefer standard templates in shared location
# - Allow user to specify basename of NIfTI output: -output_basename
# - Support single DICOM image for -func_refimg
#08/03/2016 (WF)
# - add -keep_motion_files option, mapped to KEEP_MOTION_FILES global
# - add .mean_final_func_complete so we dont redo "important last check to ensure that no blurring has occurred" 
#1/13/2016 (MH)
# - various changes to support running on PSU ACI clusters (module load etc.)
# - move log file setup to parse_args
# - support -custom_slice_times detect to pull slice times from dicom header (useful for CMRR MB sequences where these are recorded already)
# - slight tweaks to dicom header parameter detection to support both Trio and Prisma
#12/23/2015 (WF)
# - allow no_warp for 4dslicemotion+noFM; more verbose
#12/22/2015 (WF)
# - dot complete file for ref bet
#04/30/2015 (MH)
# - initial support for reference/intermediate scan to improve BBR coregistration for low GM/WM T2* contrast (e.g., MB). Use -func_refimg
# - add bias field correction to T2* coregistration targets (mc_target and epiref) to improve coregistration.
#03/18/2015 (WF)
# - added trap exit-on-error message
# - added -verbose switch
# - more involved fieldmap lock checking
#11/17/2014
# -additional cleanup for nuisance regression so that all time series are re-estimated with -startover.
# -explicit TR for 3dBandpass since this was being lost in the header during preprocessing.
#11/14/2014
# -nuisance_regression now computes CSF and WM signals from post-template-warped files using tissue masks from template (to handle poor anatomical segmentation for some scans)
#09/24/2014
# -nuisance_compute: Compute nuisance signals (same syntax as -nuisance_regression), but do not apply to the data.
# -nuisance_file: Provide user ability to name nuisance regressors file
# - Per discussion with Kai, downsample CSF and WM structural masks to functional space to compute nuisance signals (reduce computation burden)
# -compute_warp_only: Initial (slightly clunky) implementation of computing warp to standard space without applying it to preprocessed data.
# -wavelet_m1000: Voxel mean = 1000 normalization prior to wavelet despiking. Per discussion with Ameera Patel, this ensures that the 
#   spike detection threshold is applied similarly throughout the brain.
# -wavelet_threshold: Default threshold is 10, per documentation. This is reasonable for resting-state data, but some folks may with
#   to increase the threshold to decrease denoising.
# - FSL 5.0.7 removes mean component when high-pass filtering. Need to check FSL version and add this back in for normalization.
# - Can specify -hp_filter in seconds, such as -hp_filter 128s.
# - Fix input to onestep_warp for fieldmap only (-no_warp) to grab despiked data, if present.
#08/25/2014
# - Despiking was not being carried forward because warp was using postST data variable. Now that despiking comes after ST and MC, this was
#   not being incorporated in warp_to_template.
# - Also, integrated skull_strip_epi and intensity_threshold into one function. Since these share the common goal of defining a brain mask,
#   it was confusing to maintain them separately.
#08/21/2014
# - Moved despiking after skull strip and intensity threshold because we should only despike brain voxels. This also speeds up wavelet despiking
#   considerably because air voxels are ignored. 
#08/15/2014
# - Generate censor_intersection files in motion_info/ directory (using -motion_censor).
# - Move up check for .preprocessfunctional_complete so that files are untouched.
#07/24/2014
# - Move -wavelet_despike after motion correction and slice timing correction. We want to retain stationary components of fMRI time series
#   from a voxel that represents the same tissue. Running wavelet despiking before motion correction was sorting the large majority of
#   overt movement into the noise part. Further analyses on a single subject with severe movement demonstrated that although the spike
#   percentage was highly consistent for despike before or after motion, the voxelwise time series were not highly correlated (mean r ~ 0.5).
#   Discussing this with SM, BL, and WF, we agreed that one would like the voxel time series to represent the same tissue over time when running
#   the despike, although it remains somewhat magical how uni-voxel despiking is essentially parsing out motion with no knowledge of other voxels.
#   Likewise, some of these arguments apply to time despiking with 3dDespike (e.g., curve fitting a spike using data from a given chunk of tissue
#   seems like it would give the most meaningful intensities in edited time points).
#07/10/2014
# - Allow preprocessFunctional -startover call, which starts the pipeline again based on cached setttings in .preproc_cmd
#07/01/2014
# - Incorporate nuisance regression and bandpass filtering steps
# - Incorporate wavelet despiking (Patel et al. 2014 NI)
# - Better cleanup of intermediate files using -cleanup_only (looks in .preproc_cmd for settings)
#04/27/2014
# - MAJOR OVERHAUL: worker functions for each processing step abstracted to separate files in preproc_functions/
# - Can resume run by just typing preprocessFunctional at the prompt and answering "y" for resume.
# - Use fieldmap in functional -> structural coregistration when BBR is chosen. (recommended)
# - Refactored fieldmap code to be closer to FSL 5.0 (reduce dependence on FM -> EPI direct coregistration)
# - Generate volumewise motion transformation matrices from 3dvolreg using cat_matvec
# - Provide better detection of NIfTI files created by -dicom conversion so that re-running with -dicom does not fail.
# - Provide new -motion_censor option to automatically create .1D and spike regressor .mat files for high motion/signal change volumes.
# - Number of interpolations is minimized by combining motion correction, fieldmap unwarping, and warping to template whenever possible. (onestep_warp)
# - mc_target is now the standard name of the reference volume used for motion correction and functional -> structural.
# - fieldmap conversion to rad/sec now handled internally (prepare_fieldmap), not in fm_cfg files.
# - fieldmap processing works from original DICOM and can detect separate echos when the magnitude directory contains both images.
# - fieldmap processing detects and corrects if phase and magnitude directories were inadvertently swapped.
# - median and 2nd %ile not recomputed in the case of re-running script (creates files .pct2 and .median_intensity)
# - Rename log file if it already exists so that logs are not overwritten.
#
#For previous changes, see preproc_functions/changelog

#load preprocessing functions
scriptDir=$( dirname "$0" )
source "${scriptDir}/preproc_functions/waitforlock"
source "${scriptDir}/preproc_functions/helper_functions"
source "${scriptDir}/preproc_functions/fast_wmseg"
source "${scriptDir}/preproc_functions/convert_or_use_nii"
source "${scriptDir}/preproc_functions/parse_args"
source "${scriptDir}/preproc_functions/check_requirements"
source "${scriptDir}/preproc_functions/dicom_to_nifti"
source "${scriptDir}/preproc_functions/compute_motion_censor"
source "${scriptDir}/preproc_functions/correct_slice_timing"
source "${scriptDir}/preproc_functions/correct_motion"
source "${scriptDir}/preproc_functions/correct_slice_motion_4d"
source "${scriptDir}/preproc_functions/despike_timeseries"
source "${scriptDir}/preproc_functions/motion_plots"
source "${scriptDir}/preproc_functions/skullstrip_threshold"
source "${scriptDir}/preproc_functions/register_func2struct"
source "${scriptDir}/preproc_functions/remove_first_volumes"
source "${scriptDir}/preproc_functions/prepare_fieldmap"
source "${scriptDir}/preproc_functions/prepare_se_fieldmap"
source "${scriptDir}/preproc_functions/onestep_warp"
source "${scriptDir}/preproc_functions/warp_to_template"
source "${scriptDir}/preproc_functions/spatial_smooth"
source "${scriptDir}/preproc_functions/highpass_filter"
source "${scriptDir}/preproc_functions/intensity_normalization"
source "${scriptDir}/preproc_functions/nuisance_regression"
source "${scriptDir}/preproc_functions/ica_aroma"
source "${scriptDir}/preproc_functions/prepare_mc_target"
source "${scriptDir}/preproc_functions/remove_group_components"
source "${scriptDir}/preproc_functions/reg_physio"
source "${scriptDir}/preproc_functions/template_funcs"


#TODO: read dirname from the input and cd to that directory
funcdir=$(pwd) #define working directory for functional processing

##load necessary modules if running on cluster (specific to ACI at the moment)
if command -v module >/dev/null && uname -a | grep -q aci.ics.psu.edu ; then
    #setup DEPENd lab environment and programs
    source /gpfs/group/mnh5174/default/lab_resources/ni_path.bash
fi

set -e #exit if any error occurs (stop processing)
set -o errtrace # inherits trap on ERR in function and subshell

datefmt='+%F+%I:%M'
#grab command, arguments, script version (date), and start time
_datever=$(perl -ne 'print $1 if /^#Last updated: (.*)/' $0)
_gitver=$(preproc_git_ver)  # from helper_functions
thiscommandinfo="$0 $@ \n# preprocessFunctional v$_datever commit $_gitver\n# Run started $(date $datefmt) in $(pwd)"
dotfile="${funcdir}/.preprocessfunctional"

#exit script if processing complete
if [[ -f "${dotfile}_complete" && "$@" != *-startover* && "$@" != *-cleanup_only* && "$@" != *-help* ]]; then
    echo -e "\n--- preprocessFunctional already complete.\n"
    exit 0
fi

#setup default values and process all command line inputs
parse_args "$@" # 20260326: ME_INPUTS ME_ECHOTIMES

# legacy check for old (bad, pre 201810) mni template
old_template_check
warp_template_check || exit 1

#handle cleanup only scenario
#TODO: tends to leave behind some mask files since logic for which mask is subject mask is below.
if [ $cleanup_only -eq 1 ]; then
    echo "Cleaning up intermediate files, then exiting."
    cleanup_preprocessFunctional
    exit 0
fi

#setup output directories for transform matrices and qa images
[ ! -d "${funcdir}/qa_images" ] && mkdir "${funcdir}/qa_images"
[ ! -d "${funcdir}/transforms" ] && mkdir "${funcdir}/transforms"
export qa_imgdir="${funcdir}/qa_images" #use export so that TopupPreprocessingAll picks up the directory
export qa_imglog="${qa_imgdir}/qa_images.log" #use export so that TopupPreprocessingAll picks up the path

# write to an incomplete file (later to be changed to complete if process finishes)
echo -e "$thiscommandinfo " > ${dotfile}_incomplete

#handle conversion of dicom to nifti
#this step will be skipped when -4d is passed in
dicom_to_nifti # preproc_functions/dicom_to_nifti

#check for required files and settings
check_requirements

# 20260326 - add support for multiecho
#            TODO: physio also needs an array of files
#                  check_partial will only run on single echo
preproc_dir=$PWD
n_echos=${#ME_INPUTS[@]}
for echo_num in $(seq 0 $n_echos); do
   # per echo directory if we have more than one echo
   if [ $n_echos -gt 1 ]; then
      # start with 00e- so we can find with a glob later
      # will anyone ever have more than 10 echos? 02d is probably overkill
      thisechodir=$preproc_dir/$(printf "%02d" ${echo_num})e-$(basename ${ME_INPUTS[$echo_num]} .nii.gz)
      mkdir $thisechodir
      ln -s "$(abspath ${ME_INPUTS[$echo_num]})" $thisechodir
      cd $thisechodir
      ME_INPUTS[$echo_num]=$(basename ${ME_INPUTS[$echo_num]})

   fi
   funcFile=${ME_INPUTS[$echo_num]}
   #reorient images to match LPI/RPI orientation, which is useful to prevent huge rotations/translations in warp
   prefix="_" #tracks file prefix indicating which processing steps have been performed


   if [[ ! -f .reorient2std_complete || $( imtest "${prefix}${funcFile}" ) -eq 0 ]]; then
       rel "fslreorient2std ${funcFile} ${prefix}${funcFile}" #add underscore to $funcFile to make it easy later to prepend letters corresponding to each processing step
       #rel "fslmaths ${funcFile} ${prefix}${funcFile}" #alternate to skip this step.
       date > .reorient2std_complete
   fi

   # truncate volumes
   [ $n_rm_firstvols -gt 0 ] && remove_first_volumes


   # TODO: remove middleVol calc. Not used anywhere?
   imtestln $prefix${funcFile} || echo "ERROR: $prefix${funcFile} DNE, middleVol calulation will be off!"
   numVols=$( fslhd $prefix${funcFile} | grep '^dim4' | perl -pe 's/dim4\s+(\d+)/\1/' )
   middleVol=$( echo "$numVols/2" | bc )

   ######
   # motion correction
   # slice timing correction
   # fieldmap unwarping
   # registration of functional to structural
   #
   # These steps can be combined in a number of ways depending on the user specification,
   # and the code below determines the order.

   # NOTE: Jesper Andersson <jesper@fmrib.ox.ac.uk> (FSL mailing list)
   # In short: The correct thing is to do a 4D interpolation, but that is practically cumbersome. 
   # When doing them sequentially you start with the one that introduces the biggest variance/error,
   # and that is definitely movement.
   # My suggestion is 1) Movement+fieldmap 2) Slice timing.

   if [[ -z $st_first && -z $mc_first && $sliceMotion4D -eq 0 && $no_st -eq 0 ]]; then
       #no specification of order for slice timing and motion correction
       #Default:
       #  if sequential: 1) motion, 2) slice timing
       #  if interleaved: 1) slice timing, 2) motion
       if [ $sliceAcquisition == interleaved ]; then
          rel "Defaulting to slice timing correction first because of interleaved slice acquistion order." c
          st_first=1
          mc_first=0
       else
          rel "Defaulting to motion correction first because of sequential slice acquistion order." c
          mc_first=1
          st_first=0
       fi
   fi

   # regress out physio (init: 20180620WF)
   # afni_proc.py puts the 'ricor' before 'tshift' (and before 'volreg'), but after 'despike'.
   reg_physio

   # add fm config settings to gobal environment (incase we already ran prepare_fieldmap, but need it for register_func2struct
   [ -n "$fm_cfg" ] && find_and_source_fmconfig $fm_cfg

   if [ $no_mc -eq 1 ]; then
      rel "WARNING: no motion correcton. assume you're just applying this pipeline ontop of another one. This is not well tested!" c
       preMC="${prefix}${funcNifti}"
       postMC="${prefix}${funcNifti}"   

       # mc_target is used in a few places. easier to make it here
       ! imtestln mc_target && rel "3dTstat -overwrite -mean -prefix mc_target.nii.gz $preMC"
       ! imtestln mc_target_brain_restore && rel "bet  mc_target.nii.gz mc_target_brain_restore.nii.gz"
       prepare_mc_target

       [ $no_st -eq 0 ] && correct_slice_timing
       prepare_fieldmap
       func_struct_dof=6
       register_func2struct
   elif [ $sliceMotion4D -eq 1 ]; then
       #Option 1: 4-D slice timing + motion correction
       correct_slice_motion_4d
       prepare_mc_target
       prepare_fieldmap #prepare fieldmap must always come after motion correction so that mc_target exists (coregistration target)
       register_func2struct
   elif [ $no_st -eq 1 ]; then
       #just motion and FM unwarp (skip slice timing)
       correct_motion
       prepare_mc_target
       prepare_fieldmap
       register_func2struct
   elif [ $st_first -eq 1 ]; then	
       correct_slice_timing
       correct_motion
       prepare_mc_target
       prepare_fieldmap
       register_func2struct
   elif [ $mc_first -eq 1 ]; then
       correct_motion
       prepare_mc_target
       prepare_fieldmap
       register_func2struct
       onestep_warp mc_target #apply MC + FM unwarping before slice timing correction
       correct_slice_timing
   fi



   #make plots (and optional movie) of head motion
   [ $no_mc -eq 0 ] && motion_plots

   # exit in error if we want to stop after motion correction
   check_partialproc "motion" "$partial_proc"

   #Skull strip and intensity threshold EPI
   skullstrip_threshold           # preproc_functions/skullstrip_threshold

   #wavelet or time despiking after motion and slice timing correction, as well as intensity thresholding
   #want to zero non-brain voxels prior to despiking so that small fluctuations in air voxels are not flagged
   despike_timeseries

   #compute motion/intensity outliers to be excluded from analysis
   compute_motion_censor

   check_partialproc "despike" "$partial_proc"
done

if [ $n_echos -gt 1 ]; then
   preproc_echos=([0-9][0-9]e-*/$prefix*.nii.gz)
   [ ${#preproc_echos[@]} -ne $n_echos ] && rel "ERROR: expect $n_echos files instead have ${preproc_echos[*]}" c && exit 1
   skip-exist tedana/ \
      niinote __SKIPFILE \
      tedana -d "${prerpoc_echos[@]}" -e "${ME_ECHOTIMES[@]}" --out-dir tedana
   ln -s tedana/ c$prefix${funcFile}.nii.gz
   prefix=c$prefix
fi
# Warp functional to template
# NOTE: This step is skipped for -no_warp.
# NOTE: When relevant, this step also incorporates motion correction and fieldmap correction.
# These steps are combined so that there is a single interpolation into template space.
warp_to_template               # preproc_functions/warp_to_template

# 20180511 (WF) method for FC/DM to remove top group PCAs
# if we have a group componet to remove, now's the time
remove_group_components 

# handle request for computing, but not applying, warp to standard space.
# This is useful if one wishes to warp files to template space later (e.g., after GLM),
# but to complete preprocessing pipeline in native space.
if [[ "$no_warp" -eq 0 && "$compute_warp" -eq 1 && ! -f .compute_warp_complete ]]; then

    # Note: these warp coefficients only represent EPI -> template and do not include fieldmap unwarping or
    # motion correction. These steps are applied upstream by warp_to_template
    [ -z "$templateBrain" -o -z "$warpCoef" ] &&
       rel "ERROR! templateBrain ($templateBrain) or warp coef ($warpCoef) is not set!" c &&
       exit 1
    # 20220222 - mat file moved to transforms dir a while ago. never caught in here?
    func_to_struct=transforms/func_to_struct.mat
    test ! -r $func_to_struct &&
       rel "ERROR: cannot find affine matrix '$_'" c && exit 1

    rel "Computing warp to ${templateName} via $func_to_struct and $warpCoef." c
    rel "convertwarp --ref=${templateBrain} --premat=$func_to_struct --warp1=$warpCoef \
          --relout --out=func_to_${templateName}_warpfield"
    rel "Remember to use --rel with applywarp when applying this file to an image." c
    date > .compute_warp_complete
fi

# exit in error if we want to stop after warping
check_partialproc "warp" "$partial_proc"

# spatial smoothing
spatial_smooth                 # preproc_functions/spatial_smooth


check_partialproc "smooth" "$partial_proc"

# ICA-AROMA (denoising of motion artifacts)
ica_aroma                      # preproc_functions/ica_aroma

# high pass filtering
highpass_filter                # preproc_functions/highpass_filter

# normalize intensities (to allow for comparability across subjects)
intensity_normalization

# handle nuisance regression and/or bandpass filtering
nuisance_regression


#####
#Ensure that the resulting 4D NIFTI file has the correct TR. (history of FSL stripping this information out)

nuisance_reg_out_vol=$(imglob -extension ${prefix}${funcFile}${smoothing_suffix})
if [ -z "$nuisance_reg_out_vol" -o ! -r "$nuisance_reg_out_vol" ]; then
   rel "WARNING: cannot read expected nuisance_regession output: ${prefix}${funcFile}${smoothing_suffix} (matches '$nuisance_reg_out_vol')" c 
   rel "         last 2 modified files are $(ls -at |head -n2| paste -sd, )" c 
else
   #tr_hdr=$( fslhd ${prefix}${funcFile}${smoothing_suffix} | grep pixdim4 | perl -pe 's/^.*pixdim4\s+([\d+\.]+).*$/\1/' )
   tr_hdr=$(3dinfo -tr $nuisance_reg_out_vol)
   tr_hdr_nomatch=$( echo "${tr_hdr} != $tr" | bc ) 
   if [ -z "$tr_hdr" -o "$tr_hdr_nomatch" -eq 1 ]; then
       rel "Setting TR in ${prefix}${funcFile}${smoothing_suffix}.nii.gz to: $tr" c
       rel "3drefit -TR $tr \"${prefix}${funcFile}${smoothing_suffix}.nii.gz\""
   else 
       rel "nuisance_regression set tr to '$tr_hdr' matches '$tr', not adjusting '${prefix}${funcFile}${smoothing_suffix}.nii.gz' " c
   fi
fi
#####
#Compute mean functional of final 4d file and mask to be used in GLM
#First mask the final preprocessed file using the warped mask from above.

#probably not totally necessary, but seems like an important last check to ensure that no blurring has occurred.
#extents mask above, in principle, should help with this.
if [ ! -r .mean_final_func_complete ]; then
   rel "fslmaths \"${prefix}${funcFile}${smoothing_suffix}\" -mas ${subjMask} \"${prefix}${funcFile}${smoothing_suffix}\""
   rel "fslmaths \"${prefix}${funcFile}${smoothing_suffix}\" -Tmean \"${prefix}mean_func${smoothing_suffix}\""
   date > .mean_final_func_complete
fi

#if mask file exists, and subject_mask.nii is not a file (but could be a symlink), then create subject_mask symlink
if [[ -f "${subjMask}.nii"  && ! -f ./subject_mask.nii ]]; then
    rel "ln -sfn \"${subjMask}.nii\" ./subject_mask.nii"
elif [[ -f "${subjMask}.nii.gz" && ! -f ./subject_mask.nii.gz ]]; then
    rel "ln -sfn \"${subjMask}.nii.gz\" ./subject_mask.nii.gz"
fi

# 9/27/2011: note that with sinc interpolation, there tend to be some small negative values ringing the brain post-smoothing
# with an outer ring of even tinier positive values. When binarizing this, it results in a bizarre mask with an asteroid belt of sorts.
# The sensible solution is to use the mask from the thresholding step, which represents the brain voxels retained after skull stripping and
# intensity thresholding. In the warping step, I use that to mask the functional output. Need to use the same mask again here
# to mask the final file that has been smoothed. Computing a mask after all steps have finished, as I had earlier, makes little sense.
# we don't want smoothing, temporal filtering, or rescaling to influence what constitutes a brain voxel.


if [ ${deoblique_all} -eq 1 ]; then
    allNii=$( ls | grep ".*\.nii\.gz" )

    for nii in ${allNii}; do
	mv "${nii}" "toDeoblique_${nii}"

        # Run AFNI's deoblique so that it doesn't complain when looking at oblique datasets
        # Default is linear interpolation (blurry), so upgrade to quintic.
        # Note that this step necessarily resamples and interpolates the data.
	# Thus, in general, I would not recommend it since its main purpose is to quiet AFNI's warning
	3dWarp -deoblique -prefix "${nii}" -quintic "toDeoblique_${nii}"
    
	rm -f  "toDeoblique_${nii}"
    done
fi

#Cleanup files, if requested
if [ $cleanup -eq 1 ]; then
    cleanup_preprocessFunctional
fi

#move registration-related images to qa_images
if [ -n "$qa_imgdir" ]; then
    rel "Moving registration-related NIfTIs to $qa_imgdir" c
    [ $(imtest epiref_to_template) -eq 1 ] && rel "immv epiref_to_template $qa_imgdir"
    [ $(imtest epiref_to_struct) -eq 1 ] && rel "immv epiref_to_struct $qa_imgdir"
    [ $(imtest epiref_to_struct_init) -eq 1 ] && rel "immv epiref_to_struct_init $qa_imgdir"
    [ $(imtest func_to_epiref) -eq 1 ] && rel "immv func_to_epiref $qa_imgdir"
    [ $(imtest func_to_struct) -eq 1 ] && rel "immv func_to_struct $qa_imgdir"
    [ $(imtest unwarp/fmap_to_struct) -eq 1 ] && rel "immv unwarp/fmap_to_struct $qa_imgdir"
fi

# write finish time, mv incomplete to complete
echo "# finished $(date $datefmt)" >> ${dotfile}_incomplete
mv ${dotfile}_incomplete ${dotfile}_complete

#ensure that any previous crash is cleared since we've completed processing
[ -r "${dotfile}_crash" ] && rm -f "${dotfile}_crash"

exit 0

# vim: set tabstop=7: