1  #!/usr/bin/env bash
  2  
  3  deconreml(){
  4     local prefix="$1"; shift
  5     local ts4d="$1";shift
  6     local seed="$1"; shift
  7     local motionfile="$1"; shift
  8     # only use a mask if provided
  9     local maskfile="";
 10     [ $# -gt 0 ] && maskfile="$1"
 11  
 12  	# all files are readable right?
 13     for infile in ts4d seed motionfile; do
 14        ! test -r "${!infile}" -a -s "${!infile}"  && echo "Cannot open $infile '$_'" && return 1
 15     done
 16  
 17     # use mask only we hav eit to use
 18     # speeds things up potentially
 19     # but doesn't change values (TODO: is that true for reml?)
 20     [ -n "$maskfile" -a -r "$maskfile" ] &&
 21        maskswitch="-mask '$maskfile'" ||
 22        maskswitch="" 
 23  
 24  
 25     test -r "${prefix}.xmat.1D" ||
 26      3dDeconvolve \
 27        -input "$ts4d" -polort 3 -jobs "${DECON_NJOBS:-10}" \
 28        -num_stimts 1 -stim_file 1 "${seed}" -stim_label 1 "seed" \
 29        $maskswitch \
 30        -rout \
 31        -bucket "${prefix}_seeddecon.nii.gz" \
 32        -x1D "${prefix}.xmat.1D" \
 33        -censor "${motionfile}" \
 34        -overwrite
 35     
 36     test -r "${prefix}_REMLvar.nii.gz" ||
 37      OMP_NUM_THREADS=${DECON_NJOBS:-10} 3dREMLfit -matrix "${prefix}.xmat.1D"\
 38        -input "$ts4d" \
 39        $maskswitch -rout -tout \
 40        -Rbuck "${prefix}_REML.nii.gz" \
 41        -Rvar "${prefix}_REMLvar.nii.gz" \
 42        -verb \
 43        -overwrite
 44     
 45     3dcalc -r "${prefix}_REML.nii.gz[seed_R^2]" \
 46            -c "${prefix}_REML.nii.gz[seed#0_Coef]" \
 47            -expr 'ispositive(c)*sqrt(r)-isnegative(c)*sqrt(r)' \
 48            -prefix "${prefix}.nii.gz" \
 49            -overwrite
 50  }
 51  cleanup(){
 52     local prefix="$1"
 53     local bname=$(basename "$prefix" .nii.gz)
 54     mapfile -t TORM < <(find "$(dirname "$prefix")" -maxdepth 1 \
 55        -name "$bname*" \
 56        -not -name "$bname.nii.gz") #-exec rm {} \+
 57  
 58     [ ${#TORM[@]} -eq 5 ] &&
 59        rm "${TORM[@]}"  ||
 60        warn "cannot cleanup '$bname*', have ${#TORM[@]} != 5 matches"
 61        
 62  }
 63  
 64  tcor1d(){
 65     local prefix="$1"; shift
 66     local ts4d="$1";shift
 67     local seed="$1"; shift
 68     local cen="$1"; shift
 69     gidx=$(awk '($1==1){print NR-1}' "$cen" |paste -sd,)
 70     # shellcheck disable=SC1087 # afni subbrik selectors, not bash array
 71     3dTcorr1D -prefix "$prefix" "$ts4d[$gidx]" "$seed{$gidx}"
 72  }
 73  
 74  isnii(){
 75     # remove subbrik and value selector
 76     img=$(sed 's/[<[].*//' <<< "$1");
 77     [[ $CEN =~ nii(.gz)?$ ]] && return 0
 78     [[ $CEN =~ (HEAD|BRIK(.gz)?)$ ]] && return 0
 79     return 1
 80  }
 81  
 82  
 83  usage() { echo "USAGE:
 84     $0 -prefix xxx -ts yyy.nii.gz -seed t.1d -cen cen.1d -mask subj.nii.gz [-jobs 10]
 85  
 86  generate single 3d volume with per-voxel correlation between 4d timeseries and 1d seed timeseries
 87  with high motion timepoints remove (censoring)
 88  
 89  REQUIRED
 90    -prefix file.nii.gz   output. will create temporary files like file_REML*
 91    -ts     ts.nii.gz     input 4d time series 
 92    -seed   ts.1d         input 1d seed time series 
 93    -cen    cen.1d        1 include/0 excluded nrows = seed nrows = ts nvols
 94  EITHER
 95    -reml                 3dDeconvolve + 3dREMLfit
 96    -tcor                 run 3dTcorr1D instead of decon+REML (still censors, faster)
 97  
 98  OPTIONAL
 99    -noclean              dont rm temporary decon and reml files
100    -mask   subj.nii.gz   optional subject mask
101    -jobs   10            set number of cores to use for decon and reml
102  
103  
104     ";}
105  dr_main(){
106     local MASK="" TS="" SEED="" CEN=""
107     local INPUT="$*" tcor=0 reml=0 clean=1
108     #NB DECON_NJOBS not local
109     [ $# -eq 0 ] && usage && exit 1
110     while [ $# -gt 0 ]; do
111        case "$1" in
112           -prefix)   PREFIX="$2"; shift 2;;
113           -ts)       TS="$2"; shift 2;;
114           -seed)     SEED="$2"; shift 2;;
115           -cen)      CEN="$2"; shift 2;;
116           -mask)     MASK="$2"; shift 2;;
117           -tcor)     tcor=1; shift 1;;
118           -reml)     reml=1; shift 1;;
119           -noclean)  clean=0; shift 1;;
120           -jobs)     DECON_NJOBS="$2"; shift 2;;
121           -help) usage; exit;;
122        esac
123     done
124  
125     [ -z "$PREFIX" ] && echo "$0: must provide -prefix" && exit 1
126     ! [[ "$PREFIX" =~ nii.gz$ ]] && echo "-prefix should end in nii.gz" && exit 1
127     [ -r "$PREFIX" ] && warn "# '$PREFIX' exists. skipping" && exit 0
128     prefix=$(sed s/.nii.gz$// <<< "$PREFIX")
129  
130     if isnii "$SEED"; then
131        echo "# seed is an image. making ts 1D"
132        3dmaskave -queit -mask "$SEED" "$TS" > "$PREFIX.roits.1d"
133        SEED="$prefix.roits.1d"
134     fi
135  
136     # TODO:
137     # if cen is mot.tsv awk NF>0, get censor file
138  
139     # fast or slow?
140     if [ $tcor -eq 1 ]; then
141       tcor1d "$PREFIX" "$TS" "$SEED" "$CEN"
142     elif [ $reml -eq 1 ]; then
143       deconreml "$prefix" "$TS" "$SEED" "$CEN" "$MASK"
144       # remove intermediate steps
145       [ $clean -eq 1 ] && cleanup "$PREFIX"
146     else
147        warn "# must specify either -tcor or -reml!"
148        exit 1
149     fi
150  
151     # add not to final output
152     3dNotes -h "$0 $INPUT" "$PREFIX"
153  
154  }
155  
156  
157  eval "$(iffmain "dr_main")"