1  #!/usr/bin/env Rscript
 2  
 3  # calculate framewise displacement from motion paramters on stdin
 4  # 20201016WF - init
 5  # TODO: sanity checks
 6  #   columns are reasonable magnitude
 7  #   deg/rad is correct pick
 8  # 20201115 - BUG! never projected rotatin angles to displacement
 9  # also see FIACH::fd() - https://rdrr.io/cran/FIACH/src/R/fd.R
10  #
11  # and fsl_motion_outliers
12  #    # multiply rots (radians) by 50mm and add up with abs translations to get FD in Power et al, 2011
13  #    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_rot -abs -mul 50 ${outdir}_mc/res_mse_par_rot
14  #    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_rot -Tmean -mul 3 ${outdir}_mc/res_mse_par_rotsum
15  #    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_trans -abs -Tmean -mul 3 ${outdir}_mc/res_mse_par_transsum
16  #    ${FSLDIR}/bin/fslmaths ${outdir}_mc/res_mse_par_transsum -add ${outdir}_mc/res_mse_par_rotsum ${outdir}_mc/res_mse_diffZ
17  #    ${FSLDIR}/bin/fsl2ascii ${outdir}_mc/res_mse_diffZ ${outdir}_mc/res_mse_diff.txt
18  ANGLE_DISP_SPHERE <- 50 # mm
19  args <- commandArgs(trailingOnly=TRUE)
20  
21  if(length(args) < 3 || isatty(stdin())) {
22    cat('USAGE: fd_calc x:y rotx:rotz rad|deg [thres] < motion.par
23  
24   # get fd (preprocessFunctional motion.par. order: 3rads,3trans)
25   fd_calc 4:6 1:3 rad    < motion.par       > fd.txt
26   # get censor (ABCD motion par file)
27   fd_calc 1:3 4:6 deg .5 < motion.par       > censor_fd.5.txt
28   # remove header
29   sed 1d motion.par | fd_calc 1:3 4:6 deg   > fd.txt
30  
31  for more complicated situations. like censoring previous see
32  1d_tool.py and 1deval
33    echo -e "1\\n0\\n1\\n" | 1d_tool.py -censor_prev_TR -infile - -write -
34    echo -e "1\\n0\\n1\\n1\\n0" | 1deval -a stdin: -expr "step(a*z)"
35  ')
36    quit("no", status=1)
37  }
38  
39  # currently forcing at least 4 to be provided
40  # not necessary
41  opt <- list()
42  opt$tran_cols <- ifelse(length(args)>=2,args[1], '1:3')
43  opt$rot_cols  <- ifelse(length(args)>=2,args[2], '4:6')
44  opt$rot_meas  <- ifelse(length(args)>=3,args[3], 'deg')
45  opt$thres     <- ifelse(length(args)>=4,args[4], '0')
46  
47  mkidx <- function(s)  eval(parse(text=paste0('c(', s ,')')))
48  opt$tran_cols <- mkidx(opt$tran_cols)
49  opt$rot_cols  <- mkidx(opt$rot_cols)
50  opt$thres     <- as.numeric(opt$thres)
51  
52  d <- read.table(file('stdin', 'r'))
53  trans <- d[,opt$tran_cols]
54  rot_cov <- ifelse(opt$rot_meas=='deg', pi/180, 1)
55  rots <- rot_cov * d[,opt$rot_cols]
56  
57  # check rot
58  rots_mean <-mean(apply(rots,1, function(x) mean(abs(x))))
59  if(rots_mean > .01 || rots_mean < 1e-4){
60     message("#WARNING fd_calc input: rot mean rot is ", rots_mean,
61         "expect 1e-4 < x < .01. check deg vs rad\n")
62  }
63  
64  rot_disp <- sin(rots) * ANGLE_DISP_SPHERE
65  x <- cbind(trans, rot_disp)
66  diffcols <- function(x) abs(apply(x,2,diff))
67  #fd <- c(0, sqrt(rowSums(diffcols(x)^2)))
68  fd <- c(0, rowSums(diffcols(x)))
69  
70  if(any(fd<0)) stop("negative fd! something is wrong!")
71  
72  if(opt$thres > 0) {
73     output <- ifelse(fd > opt$thres, 0, 1)
74  } else {
75     output <- fd
76  }
77  outstr <- paste(collapse="\n", output)
78  cat(outstr)
79  cat("\n")
80  on.exit(close(stdin()))