fd_calc

#!/usr/bin/env Rscript

# calculate framewise displacement from motion paramters on stdin
# 20201016WF - init
# TODO: sanity checks
#   columns are reasonable magnitude
#   deg/rad is correct pick
# 20201115 - BUG! never projected rotatin angles to displacement
# also see FIACH::fd() - https://rdrr.io/cran/FIACH/src/R/fd.R
#
# and fsl_motion_outliers
#    # multiply rots (radians) by 50mm and add up with abs translations to get FD in Power et al, 2011
#    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_rot -abs -mul 50 ${outdir}_mc/res_mse_par_rot
#    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_rot -Tmean -mul 3 ${outdir}_mc/res_mse_par_rotsum
#    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_trans -abs -Tmean -mul 3 ${outdir}_mc/res_mse_par_transsum
#    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_transsum -add ${outdir}_mc/res_mse_par_rotsum ${outdir}_mc/res_mse_diffZ
#    ${FSLDIR}/bin/fsl2ascii ${outdir}_mc/res_mse_diffZ ${outdir}_mc/res_mse_diff.txt
ANGLE_DISP_SPHERE <- 50 # mm
args <- commandArgs(trailingOnly=TRUE)

if(length(args) < 3 || isatty(stdin())) {
  cat('USAGE: fd_calc x:y rotx:rotz rad|deg [thres] < motion.par

 # get fd (preprocessFunctional motion.par. order: 3rads,3trans)
 fd_calc 4:6 1:3 rad    < motion.par       > fd.txt
 # get censor (ABCD motion par file)
 fd_calc 1:3 4:6 deg .5 < motion.par       > censor_fd.5.txt
 # remove header
 sed 1d motion.par | fd_calc 1:3 4:6 deg   > fd.txt

for more complicated situations. like censoring previous see
1d_tool.py and 1deval
  echo -e "1\\n0\\n1\\n" | 1d_tool.py -censor_prev_TR -infile - -write -
  echo -e "1\\n0\\n1\\n1\\n0" | 1deval -a stdin: -expr "step(a*z)"
')
  quit("no", status=1)
}

# currently forcing at least 4 to be provided
# not necessary
opt <- list()
opt$tran_cols <- ifelse(length(args)>=2,args[1], '1:3')
opt$rot_cols  <- ifelse(length(args)>=2,args[2], '4:6')
opt$rot_meas  <- ifelse(length(args)>=3,args[3], 'deg')
opt$thres     <- ifelse(length(args)>=4,args[4], '0')

mkidx <- function(s)  eval(parse(text=paste0('c(', s ,')')))
opt$tran_cols <- mkidx(opt$tran_cols)
opt$rot_cols  <- mkidx(opt$rot_cols)
opt$thres     <- as.numeric(opt$thres)

d <- read.table(file('stdin', 'r'))
trans <- d[,opt$tran_cols]
rot_cov <- ifelse(opt$rot_meas=='deg', pi/180, 1)
rots <- rot_cov * d[,opt$rot_cols]

# check rot
rots_mean <-mean(apply(rots,1, function(x) mean(abs(x))))
if(rots_mean > .01 || rots_mean < 1e-4){
   message("#WARNING fd_calc input: rot mean rot is ", rots_mean,
       "expect 1e-4 < x < .01. check deg vs rad\n")
}

rot_disp <- sin(rots) * ANGLE_DISP_SPHERE
x <- cbind(trans, rot_disp)
diffcols <- function(x) abs(apply(x,2,diff))
#fd <- c(0, sqrt(rowSums(diffcols(x)^2)))
fd <- c(0, rowSums(diffcols(x)))

if(any(fd<0)) stop("negative fd! something is wrong!")

if(opt$thres > 0) {
   output <- ifelse(fd > opt$thres, 0, 1)
} else {
   output <- fd
}
outstr <- paste(collapse="\n", output)
cat(outstr)
cat("\n")
on.exit(close(stdin()))