lncdprep

#!/usr/bin/env bash

#
# BIDS app
#
# generic formula:
# $0 bids_root output [participant] [options]

# 20200608WF - init and POD docs

:<<=cut
=head1 NAME

lncdprep - run preprocess{Mprage,Functional} on a BIDS directory tree

=head1 SYNOPSIS

lncdprep bidsroot outputdir [--participant_label participant1 [participant2 ...]] [options]  -- [preprocessFunctionalOptions]


=head2 Options

=over 4

=item C<--onlyt1> - just run preprocessMprage

=item C<--onlyt2> - do not run preprocessMprage. use with '-- -no_warp'

=item C<--use_ses> - when finding subjects, look for ses- dir. TODO: autodetect

=item C<--njobs> - number of process to fork. default is 1.

=item C<--fm_cfg> - fieldmap configuration file. default to know. if set will look for phase and mag files

=item C<--fm_patt> - pattern to use when choosing fieldmap files (default to *.nii.gz, might want eg. "*run-1*")

=item C<--ppfmri_args|--> append to default settings. must be given last

=item C<--ppmprage_args> args to send to preprocessMprage (e.g. '-r 1YO_2mm' N.B. all args within quotes)

=item C<--pattern> - only run on pattern for functional images (default '*.nii.gz')

=item C<--help> - detailed help and examples

=item C<--tac> - like tac is to cat, reverse order of inputs

=item C<--no_partial> - skip if we have any thing in the directory we would run in

=item C<--t1> - DEPRICATED 20210519. Folder name within outputdir for preprocessMprage.
default C<MHT1_2mm>.

=item C<--t2> - DEPRICATED 20210519. Folder name within outputdir for preprocessFunctional.

=back

Within this dir., there will be subdirectoires for each C<*.nii.gz>.

=head2 preprocessFunctionalOptions

see C<preprocessFunctional -help>

=cut

# where data will be saved
PATT="*.nii.gz"
FMPATT="*.nii.gz"
FMCFG=""
MPRAGEARGS=""
MAXJOBS=1
ONLYT1=0
ONLYT2=0
USE_SES=0 # default to just sub-*/, instead of sub-*/ses-*
PARTIAL=1 # okay to rerun if something is already in the directory

# pre 20210519. separate roots for T1, T2, and FM
# now stored by protocol name in same sub/ses dir
OLDDERIVFORMAT=0
    T1DNAME="MHT1_2mm"
    T2ROOT="func"
T1TEMPLATE="MNI_2mm" # hard code ppMprage default. this must be changed if changed there
T2TEMPLATE="MNI_2mm" # might want to use 3mm
t1_tmpl_str(){
    perl -lne 'if(m/^(.*)_([0-9.]+mm)$/) {
     $t=$1;$r=$2;
     $t =~ s/_//g;
     print "tmpl-${t}_res-$r"; }' <<< "$1"
} #MNI_2mm tmp-MNI_res-2mm

usage(){ pod2usage "$0"; exit 1; }


## setup dryrun
:<<=cut
=head2 Environmental Variables

C<DRYRUN> - if set, will echo instead of running

=cut
env|grep -q ^DRYRUN=. && DRYRUN=echo || DRYRUN=""
env|grep -q ^DEBUG=. || DEBUG=""
debug() { [ -n "$DEBUG" ] && echo "$@" || :; }


## default globals
:<<=cut
=head1 DESCRIPTION

=head2 Subject list

C<SUBJS> is either provided in as arguments or from traversing the BIDS root at depth 1. matching and removing C<sub->
=cut

SUBJS=""

default_preproc_args(){
  ## default arguements for rest and task
  # single input arg should be 'task' or 'rest'
  DEF_ARGS="\
     -wavelet_despike  \
     -wavelet_threshold 10 \
     -siemens -slice_acquisition interleaved \
     -threshold 98_2 \
     -ica_aroma \
     -rescaling_method 10000_globalmedian \
     -func_struc_dof bbr  \
     -warp_interpolation spline  \
     -constrain_to_template y  \
     -smoothing_kernel default\
     -nuisance_file nuisance_regressors.txt \
     -nuisance_compute 6motion,d6motion,csf,dcsf,wm,dwm"
  
  # rest: stricter motion, bandpass, and regression
  REST_ARGS="\
    -motion_censor fd=0.3 \
    -bandpass_filter 0.009 .08 \
    -nuisance_regression 6motion,d6motion,csf,dcsf,wm,dwm"
  
  # task: more permissive motion. high pass filter 111s ~= .009Hz
  TASK_ARGS="\
     -motion_censor fd=0.8 \
     -hp_filter 111s"
  
  case "$1" in
      rest) echo "$DEF_ARGS $REST_ARGS";;
      task) echo "$DEF_ARGS $TASK_ARGS";;
      *)    echo "$FUNCNAME: unknown type '$1' (want task or rest)" >&2; return 1;;
  esac
}

## parse arguments
ARGS=""
CAT="cat" # default to subjects in order, tac reverses

mkabs() {
   local path="$1";
   [[ "$path" =~ ^/ ]] && echo "$path" && return 0
   local bname="$(dirname "$path")"
   [ ! -d "$bname" ] && echo "path does not exist: '$bname'" && return 1
   echo $(cd "$bname"; pwd)/$(basename "$path")
}

bids_to_old_deriv(){
  # before 20210519 seperate folders for everything
  # bold: $OUTPUTDIR/$T2ROOT/11757/sub-11757_task-SOA_bold 
  # t1w : t1out="$OUTDIR/$T1DNAME/$id" #pre-20210519
  # fm  : distortion_dir=$OUTDIR/FM/$id
  local input="$1"; shift
  case $input in
  *bold.nii*)
    perl -pe "
       s:/func/sub-:/sub-:g;
       s/:/_/g;s:^$BIDSROOT:$OUTDIR/$T2ROOT:;
       s/sub-//;s/.nii(.gz)?$//;" <<< $input;;
  *T1w*)
    local id=$(perl -lne 'print $1 if m:sub-([^/_-]+(/ses-[^/_-]+)?)/:'<<< $input)
    echo "$OUTDIR/$T1DNAME/$id";;
  *)
   echo "$FUNCNAME: unkonwn input not bold or T1w: '$input'" >&2
   return 1;;
 esac
}

bids_to_deriv_dir(){
 # replace BIDSROOT with our own output directory
 # $BIDSROOT/sub-11757/func/sub-11757_task-SOA_bold.nii.gz
 #  becomes
 # $OUTPUTDIR/sub-11757/task-SOA_bold 

 # still support old format
 # using if-then so bids_to_old_deriv_dir's return 1 is respected
 if [ $OLDDERIVFORMAT -eq 1 ]; then
     bids_to_old_deriv_dir "$1"
     return 0
 fi

 local input="$1"; shift
 local root="$1"; shift
 perl -pe "
   s/:/_/g;
   s:/func/sub-:/sub-:g;
   s:/anat/sub-:/sub-:g;
   s:^$BIDSROOT:$root:;
   # s/sub-//;         # 20210519 keep subject id in deriv folder
   s:sub-[^/_-]+(_ses-[^/_-]+)?_::g; # final folder is just protocol info
   s/.nii(.gz)?$//;" <<< "$input"
}
parse_args() {
   [[ $# -eq 1 && $1 =~ -?-h* ]] && pod2text "$0" && exit 1
   [ $# -lt 2 ] && usage
   BIDSROOT="$(mkabs "$1")";
   OUTDIR="$(mkabs "$2")";
   shift 2;
   while [ $# -gt 0 ]; do
    case "$1" in
      --task) ARGS="$(default_preproc_args task)"; shift 1;;
      --rest) ARGS="$(default_preproc_args rest)"; shift 1;;
      --participant_label) SUBJS="$SUBJS $2"; shift 2;;
      [^-]*) SUBJS="$SUBJS $1"; shift 1;;
      --onlyt1|--t1only) ONLYT1=1;shift;;
      --onlyt2|--t2only) ONLYT2=1;shift;;
      --use_ses) USE_SES=1; shift;;
      --njobs) MAXJOBS="$2"; shift 2;;
      --fm_patt) FMPATT="$2"; shift 2;;
      --fm_cfg) FMCFG="$2"; shift 2;;
      --pattern) PATT="$2"; shift 2;;
      --tac) CAT=tac; shift;;
      --no_partial) PARTIAL=0; shift;;
      # push unknown args onto preprocessFunctional
      --|--ppfmri_args) shift; ARGS="$ARGS $*"; break;;
      --ppmprage_args) MPRAGEARGS="$2"; shift 2;;
      # TODO: dont use a star here? should probably error with nonsense flags
      --t1*) OLDDERIVFORMAT=1; T1DNAME="$2"; shift 2;;
      --t2*) OLDDERIVFORMAT=1; T2ROOT="$2"; shift 2;;
      *) echo "ERROR: unknown arg '$1'" && usage;;
    esac
    debug "now parsing $*"
   done

   if [ -z "$ARGS" ] && [ "$ONLYT1" -ne 1 ]; then
      echo "ERROR: specify --task or --rest. Or use your own args like '$0 $* -- ....';" >&2
      echo "ERROR:    see 'preprocessFunctional -help' for options to fill ...." >&2
      exit 1
   fi
   if [[ ! "$ARGS" =~ "-no_warp" && "$ONLYT2" -eq 1 ]]; then
      echo "ERROR: have --onlyt2 but ppfrmri args do not include '-no_warp'" >&2
      echo "ERROR:    add '.... -- -no_warp' or remove '--onlyt2'" >&2
      exit 1
   fi
   # MPRAGEARGS might be used to change the reference template
   if [[ "$MPRAGEARGS" =~ (-r|-template_brain)\ [\"\']?([^ ]+) ]]; then
       T1TEMPLATE=${BASH_REMATCH[2]} # eg. MNI_2mm, 1YO_2mm
   fi
}

strip_mm(){ perl -lne 's/_[0-9.]+mm$//;print' <<<"$*"; }

t2template(){
   # echos e.g. "-template_brain MNI_3mm"
   # checks that T2TEMPLATE isn't > 1mm from functional
   # checks that T1 and T2 templates are the same space

   # no warp, don't care
   [[ "$ARGS" =~ "-no_warp" ]] && return 0

   local func="$1"; shift
   local t2res=$(3dinfo -adi "$func")
   local defres=2 # could get from T2TEMPLATE, regex _([0-9.]+)mm$

   # explicitly set then use that
   # otherwise check that we are within 1mm 
   if [[ "$ARGS" =~ \-template_brain\ [\"\']?([^ ]+) ]]; then
       T2TEMPLATE=${BASH_REMATCH[1]} # eg. MNI_3mm

   # handy defaults for MNI_2mm
   elif [[ $(bc -l <<< "$t2res >= 3") -eq 1 && $T1TEMPLATE =~ ^MNI_2mm ]]; then
       T2TEMPLATE=MNI_3mm
   elif [[ $(bc -l <<< "$t2res <= 2.0") -eq 1 && $T1TEMPLATE =~ ^MNI_2mm ]]; then
       T2TEMPLATE=MNI_2mm
   elif [[ $(bc -l <<< "$t2res >= 2.0") -eq 1 && $T1TEMPLATE =~ ^MNI_2mm ]]; then
       T2TEMPLATE=MNI_2.3mm

   # for not default. check resolution and force explicit choice if > 1mm diff
   else
       local istoobig=$(bc -l <<< "sqrt($t2res - $defres)^2 >= 1")
       if [ "$istoobig" -eq 1 ]; then
	   echo "ERROR: bold warp template '$T2TEMPLATE' > 1mm different than functional ('$t2res')
if you're okay with this, add: -- '-template_brain $T2TEMPLATE' to lncdprep
" >&2
	   return 1
       fi
   fi
   
   # template space should match across T1 and T2 preproc_functions
   # eg MNI matches MNI
   [[ "$(strip_mm ${T2TEMPLATE})" != "$(strip_mm ${T1TEMPLATE})" ]] &&
     echo "template space mismatch! t1 '$T1TEMPLATE' != t2 '$T2TEMPLATE'" >&2 &&
     return 1

   
   echo "-template_brain $T2TEMPLATE"
}


:<<=cut
=head2 preprocessMprage

We copy the sole T1w nifti to this subjects directory (C</output/t1prefix/>)
and run preprocessMprage with the default settings

=cut
mprage() {
 # uses global MPRAGEARGS. defaults to empty (preprocessMprage defaults)
 local t1out="$1"; shift
 local t1s="$1"
 [ ! -d $t1out ] && mkdir -p $t1out
 cd $t1out
 pwd
 3dcopy $t1s mprage.nii.gz -overwrite
 preprocessMprage -n mprage.nii.gz $MPRAGEARGS
}


:<<=cut
=head2 preprocessDistortion

=cut
fieldmap(){
   [ $# -ne 3 ] && echo "$FUNCNAME: want fmdir, mag, phase. given '$@'." && return 1
   local fmdir="$1"; shift
   local mag="$1"; shift
   local phase="$1"; shift
   [ ! -d "$fmdir" ] && mkdir -p "$fmdir"
   cd $fmdir
   preprocessDistortion -phase "$phase" -mag "$mag" -fm_cfg "$FMCFG" -method "gre.nii.gz"
}

# needs ARGS defined
:<<=cut
=head2 preprocessFunctional

We copy the each nifti to this subjects directory (C</output/t2prefix/nifti_name/>)
and run preprocessFunctional with the a mix of settings not necissaily the default for the sole command.

N.B. repeated command overwrite the previous setting. 
e.g. C<-smoothing_kernel 4 -smoothing_kernel 5> runs at 5mm
=cut
func() {
 local t1out="$1"; shift
 local thisout="$1"; shift
 local func="$1"; shift
 local args="$ARGS $(t2template $func)";

 if [ -n "$DRYRUN" ]; then
    echo "# DRYRUN: fork func '$thisout'"
    sleep 1
    return 0
 fi

 [ ! -d "$thisout" ] && mkdir -p "$thisout"
 cd $thisout


 # in BIDS sbref should be named nearly the same as func if it exists
 sbref=${func/_bold.nii.gz/_sbref.nii.gz}
 if [ -r "$sbref" ]; then
    local sbref_json=${sbref/.nii.gz/.json}
    test -r $sbref_json -a ! -r ref.json  && cp $sbref_json  $_

    ! test -r ref.nii.gz && 3dcopy $sbref $_
    debug "using $sbref reference!"
    args="-func_refimg ref.nii.gz $args"
 else
    debug "no sbref like '$sbref' (in $(pwd))"
 fi

 echo "final args: $args" #todo rmme

 local func_json=${func/.nii.gz/.json}
 test -r $func_json -a ! -r func.json  && cp $func_json  $_

 [ ! -r func.nii.gz ] && 3dcopy $func func.nii.gz

 preprocessFunctional -4d func.nii.gz \
   -mprage_bet $t1out/mprage_bet.nii.gz -warpcoef $t1out/mprage_warpcoef.nii.gz \
   $args

 wait
}

find_func(){
   # returns list of ids either sub-XXX or sub-XXX/ses-YYYY depending on depth ($1)
   debug "# searching $BIDSROOT at depth $1 for dir 'func'" >&2
   find -L $BIDSROOT -maxdepth "$1" -type d -name func |
      perl -lne 'print $1 if m:sub-([^/]*(/ses-[^/]*)?):' |
      $CAT
}

wait_for_jobs() {
  i=1
  SLEEPTIME=30
  while [ $(jobs -p|wc -l) -ge $MAXJOBS ]; do
     echo "[$i@$(date +%FT%H:%m)] $*; $MAXJOBS jobs. waiting $SLEEPTIME";
     let ++i
     sleep $SLEEPTIME 
  done
}

# return false if this is partially completed and we care about that
partial_ok(){
   local thisout="$1"
     [ "$PARTIAL" -eq 0 ] && [ -e "$thisout" ] &&
        [ $(find -L "$thisout" -maxdepth 1 | wc -l) -gt 1 ] &&
        echo "# skipping '$thisout'. not empty. remove --no_partial to continue" &&
        return 1
     return 0
}
find_fm() {
  # find fieldmap files (mag or phase)
  local id="$1"; shift
  local mag_or_phase="$1"; shift;
  case "$mag_or_phase" in 
   mag) mag_or_phase=".*magnitude1?.nii.gz";;
   phase) mag_or_phase=".*phase.nii.gz";;
  esac
  local file="$(find -L "$BIDSROOT/sub-$id/fmap/" -name "$FMPATT" -regex "$mag_or_phase"| sed -n '$p')"
  [ -z "$file" ] &&
	  echo "MISSING FMAP: '$BIDSROOT/sub-$id/fmap/$mag_or_phase' (using FMPATT '$FMPATT')" >&2 &&
	  return 1
  echo $file
}

_lncdprep(){
  local id="$1"; shift
  debug "_lncdprep: $id"
  if [ -n "$FMCFG" ]; then

     [ $OLDDERIVFORMAT -eq 0 ] && distortion_dir=$OUTDIR/sub-$id/FM/  ||
	     distortion_dir=$OUTDIR/FM/$id || # pre 20210519

     ARGS="$ARGS -distortion_dir $distortion_dir/unwarp"
     # dont need to do anything if we already have unwarps
     if [ ! -r $distortion_dir/unwarp/FM_UD_fmap_mag.nii.gz -o \
       ! -r $distortion_dir/unwarp/FM_UD_fmap.nii.gz  ]; then 
	     partial_ok "$distortion_dir" || return 0

	     local mag=$(find_fm $id mag) 
	     local phase=$(find_fm $id phase)
	     local fm_indir="$BIDSROOT/sub-$id/fmap/"
	     [ -z "$mag" -o -z "$phase" ] && echo "MISSING: mag+phase in '$fm_indir' ('$mag', '$phase')" && return 1
	     $DRYRUN fieldmap "$distortion_dir" "$mag" "$phase"
     fi
  fi

  t1out="$id/NO_T1w_FOUND"
  if [ "$ONLYT2" -eq 0 ]; then
     t1s=($(find -L $BIDSROOT/sub-$id -name '*T1w*.nii.gz' -path '*/anat/*'))
     [ ${#t1s[@]} -eq 0 ] && echo "ERROR: $id: no T1 in $BIDSROOT/sub-$id!?" && return
     [ ${#t1s[@]} -gt 1 ] && echo "ERROR: $id: more than one T1 ${t1s[*]}" && return

     t1out=$(bids_to_deriv_dir "${t1s[0]}" "$OUTDIR")
     [ $OLDDERIVFORMAT -eq 0 ] && t1out="${t1out}_$(t1_tmpl_str "$T1TEMPLATE")"
     # previously: t1out="$OUTDIR/$T1DNAME/$id" #pre-20210519
     if [ ! -r $t1out/mprage_warpcoef.nii.gz ]; then
	#test -r $t1out/.preprocessmprage_incomplete && echo "# rm $_ to restart" && return
	partial_ok $t1out || return 0
	$DRYRUN mprage $t1out ${t1s[0]}
     fi
  fi

  # check T1 worked for useful-ish error message
  if test "$ONLYT2" -eq 0 -a ! -r $t1out/mprage_warpcoef.nii.gz; then 
     echo "ERR: T1 failed; missing $_"
     return 1
  fi

  [ $ONLYT1 -ne 0 ] && return

  debug "searching '$BIDSROOT/sub-$id' for '$PATT'"
  find -L "$BIDSROOT/sub-$id" \
     -path '*/func/*'  \
     -iname "$PATT" \
     -name "*nii.gz" \
     -not -iname '*_sbref.nii.gz' |
   while read func; do

     test 1 -ge $(AFNI_NIFTI_TYPE_WARN=NO 3dinfo -nt "$func") &&
        echo "# skipping '$func', too few volumes ($_)" && continue

     thisout=$(bids_to_deriv_dir "$func" "$OUTDIR")
     # $BIDSROOT/sub-11757/func/sub-11757_task-SOA_bold.nii.gz
     #  becomes
     # $OUTPUTDIR/sub-11757/task-SOA_bold 

     debug "func: $func matches '$PATT'"
     partial_ok "$thisout" || continue
     func "$t1out" "$thisout" "$func"
  done
}

prep_all_subjs() {
   if [ -z "$SUBJS" ]; then
      debug "finding subjects"
      [ $USE_SES -eq 0 ] &&
         SUBJS=($(find_func 3)) ||
         SUBJS=($(find_func 4))
      debug "have ${#SUBJS[@]} subjs: ${SUBJS[@]}"
   else
      SUBJS=($SUBJS)
   fi

   total_subjs=${#SUBJS[@]}
   echo "# $total_subjs visits"
   
   for sidx in $(seq 0 $((($total_subjs-1)))); do
     id=$(echo ${SUBJS[$sidx]}| sed s/sub-//)
     debug $id
     if [ $MAXJOBS -eq 1 ]; then
        _lncdprep "$id"
     else
        _lncdprep "$id" &
        sleep .1
        wait_for_jobs "last sent: $id ($((($sidx+1)))/$total_subjs)"
     fi
   done
   wait
}

if [[ $0 =~ lncdprep$ ]]; then
   set -euo pipefail
   trap 'e=$?; [ $e -ne 0 ] && echo "$0 exited in error; wait for cleanup"; wait' EXIT
   parse_args "$@"
   echo "# $(date): Running with $MAXJOBS jobs. Can resume partial? $PARTIAL"
   prep_all_subjs
fi

:<<=cut

=head1 EXAMPLE

lncdprep BIDS/ preproc/ 233/ses-1 --t2dir aus --pattern "*AUS_run*nii.gz" --task -- -smoothing_kernel 6


=over 4

=item look in BIDS for input niftis

=item output to preproc/ directory root

=item run only for sub C<233/ses-1>

=item save to C<aus> directory (i.e. preproc/aus)

=item only run on nifti matching glob "*AUS_run*nii.gz" (i.e. BIDS/sub-233/ses-1/func/*AUS_run*nii.gz; BIDS/**/func/*AUS_run*nii.gz if no sub/ses provided)

=item use task settings

=item but change smoothing_krenl to 6 (default is 4)

=back


=cut