## ---------------------------------------------------------------------------
##
## R code to parse and read VTC and BV VOI files and support for filtering
## voxels and ROIs from VTC files.
##
## Hedderik van Rijn & Simone A. Sprenger
##
## 040214 - 040215: HvR: Intial version
## 050104:          HvR: Added some documentation and license. 
## 
## This work is licensed under the Creative Commons Attribution-ShareAlike
## License. To view a copy of this license, visit
## http://creativecommons.org/licenses/by-sa/1.0/nl/ or send a letter to
## Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
## 
## ---------------------------------------------------------------------------
##
## The format of the VTC file is documented here:
## 
## http://www.brainvoyager.com/BV2000OnlineHelp/BrainVoyagerWebHelp/mergedProjects/FileFormats/The_format_of_VTC_files.htm
##
## a backup copy of this file is available here:
##
## http://www.ai.rug.nl/~hedderik/R/VTC/formatVTC.html
##
## ---------------------------------------------------------------------------
##
## This code provides the following functions:
##
## vtc.init(filename,read)
##
## Returns a list with two components. The first component, $info, contains all
## information available in the header of the VTC file specified in the
## argument "filename" (ver.no, FMR.filename, PRT.filename, NumVol, VTC.res,
## xstart, xend, ystart, yend, zstart, zend, hemo.delay, tr, hemo.delta,
## hemo.tau, segment.size, segment.offset, dimX, dimY, dimZ). The second
## element ($data) contains a vector containing all data. This second element
## is only read when the second argument "read" is not F. If read is F, the
## returned $info element contains an additional field "f" that is set to the
## filehandle returned by open(filename). Make sure to close this file using
## vtc.close(...). Rationale for adding this "read" argument is that the amount
## of data can be significant.
##   
## vtc.read(vtc)
## 
## Reads all the data associated with the argument vtc, which is a structure as
## returned by vtc.init(). If the argument "read" of the function vtc.init() is
## not F, this function is called by vtc.init(). 
## 
## vtc.close(vtc)
## 
## Closes the VTC file associated with the argument vtc, which is a structure as
## returned by vtc.init(). If the argument "read" of the function vtc.init() is
## not F, this function is called by vtc.init(). 
## 
## vtc.get.voxel(vtc,x,y,z,coord="system")
## 
## Returns all Volumes for a certain x,y,z voxel given the vtc
## structure. (x,y,z) are by default in system coordinates, if the argument
## coord is not system, Talaraich space is used. 
## 
## vtc.read.BV.voi.file(filename,name)
## 
## Reads and parses a BrainVoyager VOI file and searches in for the region of
## interest named "name". Returns a dataframe containing all voxels associated
## with this region of interest.
## 
## vtc.get.roi(vtc,roi)
## 
## Returns averaged values for the specified ROI. The argument roi should be a
## dataframe with x,y,z columns specifying the voxels of interest (in system
## coordinates).
## 
## ---------------------------------------------------------------------------

vtc.init <- function(filename,read=F) {

  vtc <- list()
  
  vtc$info <- list()
  vtc$data <- NA
  
  vtc$info$f <- file(filename,"rb")
  
  ## Read in version number
  vtc$info$ver.no <- readBin(vtc$info$f,"int",1,2,endian="swap")
  
  ## Read in filename
  vtc$info$FMR.filename <- readBin(vtc$info$f,"char",1,1,endian="swap")
  
  ## and another one:
  vtc$info$PRT.filename <- readBin(vtc$info$f,"char",1,1,endian="swap")

  ## NrOfVolumes
  vtc$info$NumVol <- readBin(vtc$info$f,"int",1,2,endian="swap")

  ## VTC-Resolution 
  vtc$info$VTC.res = readBin(vtc$info$f,"int",1,2,endian="swap");

  vtc$info$xstart = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$xend = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$ystart = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$yend = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$zstart = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$zend = readBin(vtc$info$f,"int",1,2,endian="swap");

  vtc$info$hemo.delay = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$tr = readBin(vtc$info$f,"double",1,4,endian="swap");
  vtc$info$hemo.delta = readBin(vtc$info$f,"double",1,4,endian="swap");
  vtc$info$hemo.tau = readBin(vtc$info$f,"double",1,4,endian="swap");

  vtc$info$segment.size = readBin(vtc$info$f,"int",1,2,endian="swap");
  vtc$info$segment.offset = readBin(vtc$info$f,"int",1,2,endian="swap");

  vtc$info$dimX = (vtc$info$xend - vtc$info$xstart)/vtc$info$VTC.res;
  vtc$info$dimY = (vtc$info$yend - vtc$info$ystart)/vtc$info$VTC.res;
  vtc$info$dimZ = (vtc$info$zend - vtc$info$zstart)/vtc$info$VTC.res;

  if (read) {
    vtc <- vtc.read(vtc)
    vtc <- vtc.close(vtc)
  }
  
  vtc
}

vtc.close <- function(vtc) {
  close(vtc$info$f)
  vtc$info$f <- NULL
  vtc
}

vtc.read <- function(vtc) {
  totValues <- vtc$info$dimX * vtc$info$dimY * vtc$info$dimZ * vtc$info$NumVol
  vtc$data  <- readBin(vtc$info$f,"integer",totValues,2,signed=FALSE,endian="swap")
  vtc
}

## ---------------------------------------------------------------------------

vtc.get.voxel <- function(vtc,x,y,z,coord="system") {

  if(!(x>vtc$info$xstart+vtc$info$VTC.res-1 &
       y>vtc$info$ystart+vtc$info$VTC.res-1 &
       z>vtc$info$zstart+vtc$info$VTC.res-1)) {
    stop("Brainvoyager doesn't alow exporting these values... so neither do we. :-)")
  }

  if(!(x<vtc$info$xend &
       y<vtc$info$yend &
       z<vtc$info$zend)) {
    stop("Brainvoyager doesn't alow exporting these values... so neither do we. :-)")
  }
  
  nX <- vtc$info$dimX
  nY <- vtc$info$dimY
  nZ <- vtc$info$dimZ
  nV <- vtc$info$NumVol

  dX <- nV
  dY <- nX * dX
  dZ <- nY * dY

  if (coord == "system") {
    x <- trunc((x - vtc$info$xstart - 1)/vtc$info$VTC.res)
    y <- trunc((y - vtc$info$ystart - 1)/vtc$info$VTC.res)
    z <- trunc((z - vtc$info$zstart - 1)/vtc$info$VTC.res)
  }
  
  start <- z*dZ+y*dY+x*dX

  vtc$data[(start+1):(start+nV)]
}

vtc.get.roi <- function(vtc,roi) {

  stopifnot(is.data.frame(roi))
  if (sum(c("x","y","z") %in% names(roi))!=3) {
    stop("The dataframe needs to contain x, y, and z columns!")
  }

  nV <- vtc$info$NumVol
  vals <- rep(0,nV)
  
  for (i in 1:nrow(roi)) {
    cz <- roi[i,"z"]
    cx <- roi[i,"x"]
    cy <- roi[i,"y"]
    vals <- vals + vtc.get.voxel(vtc,cx,cy,cz,coord="system")
  }
  
  vals/nrow(roi)
}

vtc.read.BV.voi.file <- function(filename,name) {
  f <- file(filename,"r")
  ws <- "[ \t]+"

  ready <- F
  lineNo <- 0
  while (!ready) {
    lineNo <- lineNo + 1;
    curline <- strsplit(readLines(f,1),": ")

    if (length(curline)==0) {
      stop("Could not file specified VOI");
    }
    
    if (length(curline[[1]]) > 0) {

      ## If not SYStem coordinates, barf
      if (curline[[1]][1] == "CoordsType" &
          curline[[1]][2] != "SYS") {
        stop("Only know how to handle SYStem coordinates!")
      }
        
      ## Determine if this is the VOI we're interested in:
      if (curline[[1]][1] == "NameOfVOI") {
        if (length(curline[[1]]) == 1) {
          stop(paste("Need to specify an name for each VOI! (check line ",lineNo,")",sep=""));
        }
        if (curline[[1]][2] == name) {
          skip <- F
        } else {
          skip <- T
        }
      }
      ## if so (skip==F), and we found the NrOfVoxels field, we're set:
      if (curline[[1]][1] == "NrOfVoxels" &&
          !skip) {
        ready <- T;
        numVox <- as.integer(curline[[1]][2])
      }
    }
  }

  df <- data.frame(x=rep(NA,numVox),
                   y=rep(NA,numVox),
                   z=rep(NA,numVox))

  curline <- strsplit(readLines(f,numVox),ws)
  for (i in 1:numVox) {
    df[i,] <- as.numeric(curline[[i]])+1
  }
  
  close(f)
  df
}

## ---------------------------------------------------------------------------