R/io.R

In fmri: Analysis of fMRI Experiments

read.ANALYZE.header <- function(filename) {
  if (is.na(file.info(filename)$size) | (file.info(filename)$size != 348))
    stop("Hmmm! This does not seem to be an ANALYZE header! Wrong size or does not exist!");

  con <- file(filename,"rb")

  header <- list()

  # read 4 bytes and get the endianess by comparing filesize with 348
  endian <- if ((sizeofhdr <- readBin(con,"int",1,4,endian="little")) == 348) "little" else "big"
  header$sizeofhdr <- 348
  header$endian <- endian

  header$datatype1 <- readChar(con,10, TRUE)
  header$dbname <- readChar(con,18, TRUE)
  header$extents <- readBin(con,"int",1,4,endian=endian)
  header$sessionerror <- readBin(con,"int",1,2,endian=endian)
  header$regular <- readChar(con,1, TRUE)
  header$hkey <- readChar(con,1, TRUE)
  header$dimension <- readBin(con,"int",8,2,endian=endian)
  header$unused <- readBin(con,"int",7,2,endian=endian)
  header$datatype <- readBin(con,"int",1,2,endian=endian)
  header$bitpix <- readBin(con,"int",1,2,endian=endian)
  header$dimun0 <- readBin(con,"int",1,2,endian=endian)
  header$pixdim <- readBin(con,"double",8,4,endian=endian)
  header$voxoffset <- readBin(con,"double",1,4,endian=endian)
  header$funused <- readBin(con,"double",3,4,endian=endian)
  header$calmax <- readBin(con,"double",1,4,endian=endian)
  header$calmin <- readBin(con,"double",1,4,endian=endian)
  header$compressed <- readBin(con,"double",1,4,endian=endian)
  header$verified <- readBin(con,"double",1,4,endian=endian)
  header$glmax <- readBin(con,"int",1,4,endian=endian)
  header$glmin <- readBin(con,"int",1,4,endian=endian)
  header$describ <- readChar(con,80, TRUE)
  header$auxfile<- readChar(con,24, TRUE)
  header$orient <- readChar(con,1, TRUE) # is this really a character?
  header$originator <- readBin(con,"int",5,2,endian=endian) # documented as 10 byte character!!
  header$generated <- readChar(con,10, TRUE)
  header$scannum <- readChar(con,10, TRUE)
  header$patientid <- readChar(con,10, TRUE)
  header$expdate <- readChar(con,10, TRUE)
  header$exptime <- readChar(con,10, TRUE)
  header$histun0 <- readChar(con,3, TRUE)
  header$views <- readBin(con,"int",1,4,endian=endian)
  header$voladded<- readBin(con,"int",1,4,endian=endian)
  header$startfield <- readBin(con,"int",1,4,endian=endian)
  header$fieldskip <- readBin(con,"int",1,4,endian=endian)
  header$omax <- readBin(con,"int",1,4,endian=endian)
  header$omin <- readBin(con,"int",1,4,endian=endian)
  header$smax <- readBin(con,"int",1,4,endian=endian)
  header$smin <- readBin(con,"int",1,4,endian=endian)

  close(con)

  header
}

read.ANALYZE.volume <- function(filename) {
  file.name <- substring(filename, 1, nchar(filename) - 4)
  file.hdr <- paste(file.name, ".hdr", sep = "")
  file.img <- paste(file.name, ".img", sep = "")

  header <- read.ANALYZE.header(file.hdr)
  (dx <- header$dimension[2]) || (dx <- 1)
  (dy <- header$dimension[3]) || (dy <- 1)
  (dz <- header$dimension[4]) || (dz <- 1)
  (dt <- header$dimension[5]) || (dt <- 1)
  endian <- header$endian
  if (header$datatype == 1) { # logical
    what <- "raw"
    signed <- TRUE
    size <- 1
  } else if (header$datatype == 2) { # unsigned int 1-byte
    what <- "int"
    signed <- FALSE
    size <- if (header$bitpix) header$bitpix/8 else 1
  } else if (header$datatype == 4) { # signed short
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 2
  } else if (header$datatype == 8) { # signed integer
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 4
  } else if (header$datatype == 16) { # float
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 4
  } else if (header$datatype == 32) { # complex
    what <- "complex"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 8
  } else if (header$datatype == 64) { # double
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 8
  } else { # all other
    what <- "raw"
    signed <- TRUE
    size <- 1
  }

  con <- file(filename,"rb")
  ttt <- readBin(con, what, n=dx*dy*dz*dt*size, size=size, signed=signed, endian=endian)
  close(con)

  dim(ttt) <- c(dx,dy,dz,dt)

  invisible(list(ttt=ttt,header=header))
}

write.ANALYZE.header <- function(header,filename) {
  con <- file(paste(filename, ".hdr", sep=""), "wb")

  writeBin(as.integer(348), con, 4)
  writeChar(header$datatype1, con, 10, NULL)
  writeChar(header$dbname, con, 18, NULL)
  writeBin(as.integer(header$extents), con, 4)
  writeBin(as.integer(header$sessionerror), con, 2)
  writeChar(header$regular, con, 1, NULL)
  writeChar(header$hkey, con, 1, NULL)
  writeBin(as.integer(header$dimension), con, 2)
  writeBin(as.integer(header$unused), con, 2)
  writeBin(as.integer(header$datatype), con, 2)
  writeBin(as.integer(header$bitpix), con, 2)
  writeBin(as.integer(header$dimun0), con, 2)
  writeBin(header$pixdim, con, 4)
  writeBin(header$voxoffset, con, 4)
  writeBin(header$funused, con, 4)
  writeBin(header$calmax, con, 4)
  writeBin(header$calmin, con, 4)
  writeBin(header$compressed, con, 4)
  writeBin(header$verified, con, 4)
  writeBin(as.integer(header$glmax), con, 4)
  writeBin(as.integer(header$glmin), con, 4)
  writeChar(header$describ, con, 80, NULL)
  writeChar(header$auxfile, con, 24, NULL)
  writeChar(header$orient, con, 1, NULL) # is this really a character?
  writeBin(as.integer(header$originator), con, 2) # documented as 10 byte character!!
  writeChar(header$generated, con, 10, NULL)
  writeChar(header$scannum, con, 10, NULL)
  writeChar(header$patientid, con, 10, NULL)
  writeChar(header$expdate, con, 10, NULL)
  writeChar(header$exptime, con, 10, NULL)
  writeChar(header$histun0, con, 3, NULL)
  writeBin(as.integer(header$views), con, 4)
  writeBin(as.integer(header$voladded), con, 4)
  writeBin(as.integer(header$startfield), con, 4)
  writeBin(as.integer(header$fieldskip), con, 4)
  writeBin(as.integer(header$omax), con, 4)
  writeBin(as.integer(header$omin), con, 4)
  writeBin(as.integer(header$smax), con, 4)
  writeBin(as.integer(header$smin), con, 4)

  close(con)
}

write.ANALYZE.volume <- function(ttt,filename,size) {
  con <- file(paste(filename, ".img", sep=""), "wb")
  dim(ttt) <- NULL
  writeBin(as.integer(ttt), con, size)
  close(con)
}

read.ANALYZE <- function(prefix = c(""), numbered = FALSE, postfix = "",
              picstart = 1, numbpic = 1, level=0.75, mask=NULL, setmask=TRUE) {
  counter <- c(paste("00", 1:9, sep=""), paste("0", 10:99, sep=""),paste(100:999,sep=""));

  prefix <- strsplit(prefix,".img")
  filename <- character(length(prefix))
  if (numbered) {
    for (i in picstart:(picstart+numbpic-1)) {
      filename[i-picstart+1] <- paste(prefix[[1]][1], counter[i], postfix, ".img", sep="")
    }
  } else {
    for (i in 1:length(prefix)) {
      if (length(prefix[[i]]) > 1)
        warning("filename",paste(prefix[[i]],collapse=""),"probably misinterpreted due to extension-like .img in name!")
      filename[i] <- paste(prefix[[i]][1], ".img", sep="")
    }
  }

  if (!is.na(file.info(filename[1])$size)) {
    analyze <- read.ANALYZE.volume(filename[1]);
    ttt <- analyze$ttt
    ddim <- dim(ttt)
    cat(".")
    header <- analyze$header;

    if ((numbpic > 1) && numbered) {
      for (i in 2:numbpic) {
        a <- read.ANALYZE.volume(filename[i])$ttt
        if (sum() != 0)
          cat("Error: wrong spatial dimension in picture",i)
        ttt <- c(ttt,a);
        ddim[4] <- ddim[4] + dim(a)[4]
        cat(".")
      }
    }

    dim(ttt) <- ddim

    if (min(abs(header$pixdim[2:4])) != 0) {
      weights <-
        abs(header$pixdim[2:4]/min(abs(header$pixdim[2:4])))
    } else {
      weights <- NULL
    }
    delta <- header$pixdim[2:4]
    mask0 <- array(TRUE,ddim[1:3])
    if (setmask) {
      mask0[ttt[,,,1] < quantile(ttt[,,,1],level,na.rm = TRUE)] <- FALSE
      dim(ttt) <- c(prod(dim(ttt)[1:3]),dim(ttt)[4])
      na <- ttt %*% rep(1,dim(ttt)[2])
      mask0[is.na(na)] <- FALSE
      ttt[is.na(na),] <- 0
      dim(mask0) <- ddim[1:3]
      mask0 <- connect.mask(mask0)
    }
    z <- list(ttt=writeBin(as.numeric(ttt),raw(),4),
              format="ANALYZE",
              delta=delta,
              origin=NULL,
              orient=NULL,
              dim=ddim,
              roixa=1,
              roixe=ddim[1],
              roiya=1,
              roiye=ddim[2],
              roiza=1,
              roize=ddim[3],
              roit=1:ddim[4],
              weights=weights,
              header=header,
              mask=mask0)
  } else {
    warning(paste("Error: file",filename[1],"does not exist!"))
    z <- list(ttt=NULL,
              format="ANALYZE",
              delta=NULL,
              origin=NULL,
              orient=NULL,
              dim=NULL,
              roixa=NULL,
              roixe=NULL,
              roiya=NULL,
              roiye=NULL,
              roiza=NULL,
              roize=NULL,
              roit=NULL,
              weights=NULL,
              header=NULL,
              mask=NULL)
  }

  z <- complete.fmriHeader(z)
  class(z) <- "fmridata"

  if (length(filename) > 1) {
    attr(z,"file") <- paste(filename[1], "...",
                            filename[length(filename)],
                            sep="")
  } else {
    attr(z,"file") <- paste(filename[1], sep="")
  }
  if(!is.null(mask)) z <- setmask(z,mask)
  invisible(z)
}



write.ANALYZE <- function(ttt, header=NULL, filename) {
  if (is.null(header)) header <- list()

  if (!("datatype1" %in% names(header))) header$datatype1 <- paste(rep(" ",10),collapse="")
  if (!("dbname" %in% names(header))) header$dbname <- paste(rep(" ",18),collapse="")
  if (!("extents" %in% names(header))) header$extents <- c(0)
  if (!("sessionerror" %in% names(header))) header$sessionerror <- c(0)
  if (!("regular" %in% names(header))) header$regular <- "r"
  if (!("hkey" %in% names(header))) header$hkey <- " "
  if (!("dimension" %in% names(header))) {header$dimension <- rep(0,8); header$dimension <- c(length(dim(ttt)),dim(ttt))}
  if (length(header$dimension) < 8) header$dimension[(length(header$dimension)+1):8] <- 0
  if (!("unused" %in% names(header))) header$unused <- rep(0,7)
  if (length(header$unused) < 7) header$unused[(length(header$unused)+1):7] <- 0
  if (!("datatype" %in% names(header))) header$datatype <- c(4)
  if (!("bitpix" %in% names(header))) header$bitpix <- c(0)
  if (!("dimun0" %in% names(header))) header$dimun0 <- c(0)
  if (!("pixdim" %in% names(header))) header$pixdim <- c(0,4,4,4,rep(0,4))
  if (length(header$pixdim) < 8) header$pixdim[(length(header$pixdim)+1):8] <- 0
  if (!("voxoffset" %in% names(header))) header$voxoffset <- c(0)
  if (!("funused" %in% names(header))) header$funused <- rep(0,3)
  if (length(header$funused) < 3) header$funused[(length(header$funused)+1):3] <- 0
  if (!("calmax" %in% names(header))) header$calmax <- c(0)
  if (!("calmin" %in% names(header))) header$calmin <- c(0)
  if (!("compressed" %in% names(header))) header$compressed <- c(0)
  if (!("verified" %in% names(header))) header$verified <- c(0)
  if (!("glmax" %in% names(header))) header$glmax <- c(0)
  if (!("glmin" %in% names(header))) header$glmin <- c(0)
  if (!("describ" %in% names(header))) header$describ <- paste(rep(" ",80),collapse="")
  if (!("auxfile" %in% names(header))) header$auxfile <- paste(rep(" ",24),collapse="")
  if (!("orient" %in% names(header))) header$orient <- " "
  if (!("originator" %in% names(header))) header$originator <- rep(0,5)
  if (length(header$originator) < 5) header$originator[(length(header$originator)+1):8] <- 0
  if (!("generated" %in% names(header))) header$generated <- paste(rep(" ",10),collapse="")
  if (!("scannum" %in% names(header))) header$scannum <- paste(rep(" ",10),collapse="")
  if (!("patientid" %in% names(header))) header$patientid <- paste(rep(" ",10),collapse="")
  if (!("expdate" %in% names(header))) header$expdate <- paste(rep(" ",10),collapse="")
  if (!("exptime" %in% names(header))) header$exptime <- paste(rep(" ",10),collapse="")
  if (!("histun0" %in% names(header))) header$histun0 <- paste(rep(" ",3),collapse="")
  if (!("views" %in% names(header))) header$views <- c(0)
  if (!("voladded" %in% names(header))) header$voladded <- c(0)
  if (!("startfield" %in% names(header))) header$startfield <- c(0)
  if (!("fieldskip" %in% names(header))) header$fieldskip <- c(0)
  if (!("omax" %in% names(header))) header$omax <- c(0)
  if (!("omin" %in% names(header))) header$omin <- c(0)
  if (!("smax" %in% names(header))) header$smax <- c(0)
  if (!("smin" %in% names(header))) header$smin <- c(0)

  if (header$datatype == 1) { # logical
    what <- "raw"
    signed <- TRUE
    size <- 1
  } else if (header$datatype == 2) { # unsigned char????
    what <- "int"
    signed <- FALSE
    size <- if (header$bitpix) header$bitpix/8 else 2
  } else if (header$datatype == 4) { # signed short
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 2
  } else if (header$datatype == 8) { # signed integer
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 4
  } else if (header$datatype == 16) { # float
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 4
  } else if (header$datatype == 32) { # complex
    what <- "complex"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 8
  } else if (header$datatype == 64) { # double
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8 else 8
  } else { # all other
    what <- "raw"
    signed <- TRUE
    size <- 1
  }

  write.ANALYZE.header(header,filename)

  write.ANALYZE.volume(ttt, filename,size)
}



read.AFNI <- function(filename, vol=NULL, level=0.75, mask=NULL, setmask=TRUE){
  fileparts <- strsplit(filename,"\\.")[[1]]

  if (length(fileparts) == 1) {
    filename.head <- paste(fileparts[1],"HEAD",sep=".")
    filename.brik <- paste(fileparts[1],"BRIK",sep=".")
  } else {
    ext <- tolower(fileparts[length(fileparts)])
    if (ext == "head") {
      filename.head <- filename
      filename.brik <- paste(c(fileparts[-length(fileparts)],"BRIK"),collapse=".")
    } else if (ext == "brik") {
      filename.head <- paste(c(fileparts[-length(fileparts)],"HEAD"),collapse=".")
      filename.brik <- filename
    } else if (ext == "") {
      filename.head <- paste(c(fileparts,"HEAD"),collapse=".")
      filename.brik <- paste(c(fileparts,"BRIK"),collapse=".")
    }
  }

  conhead <- file(filename.head,"r")
  header <- readLines(conhead)
  close(conhead)

  types <- NULL
  args <- NULL
  counts <- NULL
  values <- NULL

  for (i in 1:length(header)) {
    if (regexpr("^type *= *", header[i]) != -1) {
      tmptype <- strsplit(header[i]," *= *")[[1]][2]
      types <- c(types,tmptype)
      args <- c(args,strsplit(header[i+1]," *= *")[[1]][2])
      tmpcounts <- as.numeric(strsplit(header[i+2]," *= *")[[1]][2])
      counts <- c(counts,tmpcounts)
      i <- i+3
      tmpvalue <- ""
      while ((regexpr("^$", header[i]) == -1) && (i <= length(header))) {
        tmpvalue <- paste(tmpvalue,header[i])
        i <- i+1
      }
      tmpvalue <- sub("^ +","",tmpvalue)
      if ((tmptype == "integer-attribute") || (tmptype == "float-attribute")) {
        tmpvalue <- as.numeric(strsplit(tmpvalue," +")[[1]])
      } else {
        tmpvalue <- sub("~$","",sub("^\'","",tmpvalue))
      }
      values <- c(values,list(value=tmpvalue))
    }
  }

  names(values) <- args

  dx <- values$DATASET_DIMENSIONS[1]
  dy <- values$DATASET_DIMENSIONS[2]
  dz <- values$DATASET_DIMENSIONS[3]
  ddim <- c(dx,dy,dz)
  dt <- values$DATASET_RANK[2]
  scale <- values$BRICK_FLOAT_FACS
  size <- file.info(filename.brik)$size/(dx*dy*dz*dt)
  bricktypes <- values$BRICK_TYPES[1]
#  orientation <- values$ORIENT_SPECIFIC
#  if (any(sort((orientation)%/%2) != 0:2)) stop("invalid orientation",orientation,"found! \n")
  delta <- values$DELTA

  if (is.null(bricktypes)) {
    if (size == 2) { what <- "int" }
    else if (size == 4) { what <- "double" }
    else if (size == 16) { what <- "complex" }
    else { what <- "int" }
  } else {
    if (bricktypes == 1) { what <- "int" }
    else if (bricktypes == 3) { what <- "double" }
    else if (bricktypes == 5) { what <- "complex" }
    else { what <- "int" }
  }

  if (regexpr("MSB",values$BYTEORDER_STRING[1]) != -1) {
    endian <- "big"
  } else {
    endian <- "little"
  }

  if (min(abs(delta)) != 0) {
    weights <-
      abs(delta/min(abs(delta)))
  } else {
    weights <- NULL
  }

  if (is.null(vol)) vol <- 1:dt

  if ((as.integer(size) == size) && (length(vol) > 0)) {
    myttt <- array(0,dim=c(ddim,length(vol)))
    kk <- 1
    conbrik <- file(filename.brik,"rb")
    for (k in 1:dt) {
      if (k %in% vol) {
        temp <- readBin(conbrik, what, n=prod(ddim), size=size,
                        signed=TRUE, endian=endian)
        dim(temp) <- c(dx,dy,dz)
        myttt[,,,kk] <- temp
        rm(temp)
        if (scale[k] != 0) {
          cat("scale",k,"with",scale[k],"\n")
          cat(range(myttt[,,,kk]),"\n")
          myttt[,,,kk] <- scale[k] * myttt[,,,kk]
          cat(range(myttt[,,,kk]),"\n")
        }
        kk <- kk + 1
      } else {
        readBin(conbrik, what, n=dx*dy*dz, size=size, signed=TRUE, endian=endian)
      }
    }
    close(conbrik)

#
#   set correct orientation
#   DO IT IN plot.fmridata()
#
#    xyz <- (orientation)%/%2+1
#    swap <- orientation-2*(orientation%/%2)
#    if(any(xyz!=1:3)) {
#      abc <- 1:3
#      abc[xyz] <- abc
#      myttt <- aperm(myttt,c(abc,4))
#      swap[xyz] <- swap
#      delta[xyz] <- delta
#      weights[xyz] <- weights
#      ddim[xyz]<- ddim
#    }
#    if(swap[1]==1) {
#      myttt <- myttt[ddim[1]:1,,,,drop=FALSE]
#    }
#    if(swap[2]==1) {
#      myttt <- myttt[,ddim[2]:1,,,drop=FALSE]
#    }
#    if(swap[3]==0) {
#      myttt <- myttt[,,ddim[3]:1,,drop=FALSE]
#    }

    mask0 <- array(TRUE,ddim)
    if (setmask) {
      mask0[myttt[,,,1] < quantile(myttt[,,,1],level,na.rm = TRUE)] <- FALSE
      dim(myttt) <- c(prod(ddim),dim(myttt)[4])
      na <- myttt %*% rep(1,dim(myttt)[2])
      mask0[is.na(na)] <- FALSE
      myttt[is.na(na),] <- 0
      dim(mask0) <- ddim
      mask0 <- connect.mask(mask0)
    }
    z <-
      list(ttt=writeBin(as.numeric(myttt),raw(),4),
           format="HEAD/BRIK",
           delta=delta,
           origin=values$ORIGIN,
           orient=c(0,2,5),
           dim=c(ddim,length(vol)),
           roixa=1,
           roixe=ddim[1],
           roiya=1,
           roiye=ddim[2],
           roiza=1,
           roize=ddim[3],
           roit=vol,
           weights=weights,
           header=values,
           mask=mask0)
  } else {
    warning("Error reading file: Could not detect size per voxel\n")
    z <- list(ttt=NULL,
              format="HEAD/BRIK",
              delta=NULL,
              origin=NULL,
              orient=NULL,
              dim=NULL,
              roixa=NULL,
              roixe=NULL,
              roiya=NULL,
              roiye=NULL,
              roiza=NULL,
              roize=NULL,
              roit=NULL,
              weights=NULL,
              header=values,
              mask=NULL)
  }

  z <- complete.fmriHeader(z)
  class(z) <- "fmridata"
  attr(z,"file") <- paste(filename,".HEAD/BRIK",sep="")
  if(!is.null(mask)) z <- setmask(z,mask)
  invisible(z)
}

write.AFNI <- function(filename, ttt, label=NULL, note=NULL, origin=NULL, delta=NULL, idcode=NULL, header=NULL, taxis = FALSE) {

  afni.header <- list()
  # afni.header$NAME <- c("NAME", "type", category, count, reserved-count, multiple)
  # "NAME"         : Attribute name (string)
  # "type"         : Attribute type (one of "string", "integer", "float")
  # category       : according to AFNI Doc. 0 is mandatory, 1 only for time series
  # count          : number of parameters (used), 0 is unknown, or not implemented
  # reserved-count : number of parameters (reserved), 0 is unknown, or not implemented
  # multiple       : is this a multiple argument like NOTE_NUMBER_001

  # Mandatory Attributes
  afni.header$DATASET_RANK <- c("DATASET_RANK", "integer", 0, 2, 8, FALSE)
  afni.header$DATASET_DIMENSIONS <- c("DATASET_DIMENSION", "integer", 0, 3, 5, FALSE)
  afni.header$TYPESTRING <- c("TYPESTRING", "string", 0, 1, 1, FALSE)
  afni.header$SCENE_DATA <- c("SCENE_DATA", "integer", 0, 3, 8, FALSE)
  afni.header$ORIENT_SPECIFIC <- c("ORIENT_SPECIFIC", "integer", 0, 3, 3, FALSE)
  afni.header$ORIGIN <- c("ORIGIN", "float", 0, 3, 3, FALSE)
  afni.header$DELTA <- c("DELTA", "float", 0, 3, 3, FALSE)

  # Time-dependent Dataset Attributes
  afni.header$TAXIS_NUMS <- c("TAXIS_NUMS", "integer", 1, 3, 8, FALSE)
  afni.header$TAXIS_FLOATS <- c("TAXIS_FLOATS", "float", 1, 5, 8, FALSE)
  afni.header$TAXIS_OFFSETS <- c("TAXIS_OFFSETS", "float", 1, 0, 0, FALSE)

  # Almost Mandatory Attributes
  afni.header$IDCODE_STRING <- c("IDCODE_STRING", "string", 2, 1, 1, FALSE)
  afni.header$IDCODE_DATE <- c("IDCODE_DATE", "string", 2, 1, 1, FALSE)
  afni.header$BYTEORDER_STRING <- c("BYTEORDER_STRING", "string", 2, 1, 1, FALSE)
  afni.header$BRICK_STATS <- c("BRICK_STATS", "float", 2, 0, 0, FALSE)
  afni.header$BRICK_TYPES <- c("BRICK_TYPES", "integer", 2, 0, 0, FALSE)
  afni.header$BRICK_FLOAT_FACS <- c("BRICK_FLOAT_FACS", "float", 2, 0, 0, FALSE)
  afni.header$BRICK_LABS <- c("BRICK_LABS", "string", 2, 0, 0, FALSE)
  afni.header$BRICK_STATAUX <- c("BRICK_STATAUX", "float", 2, 0, 0, FALSE)
  afni.header$STAT_AUX <- c("STAT_AUX", "float", 2, 0, 0, FALSE) # depreciated

  # Note Attributes
  afni.header$HISTORY_NOTE <- c("HISTORY_NOTE", "string", 3, 1, 1, FALSE)
  afni.header$NOTES_COUNT <- c("NOTES_COUNT", "integer", 3, 1, 1, FALSE)
  afni.header$NOTE_NUMBER <- c("NOTE_NUMBER", "string", 3, 0, 0, TRUE) # more than one!!!

  # Registration Attributes
  afni.header$TAGALIGN_MATVEC <- c("TAGALIGN_MATVEC", "float", 4, 12, 12, FALSE)
  afni.header$VOLREG_MATVEC <- c("VOLREG_MATVEC", "float", 4, 12, 12, TRUE) # more than one!!!
  afni.header$VOLREG_ROTCOM <- c("VOLREG_ROTCOM", "string", 4, 1, 1, TRUE) # more than one!!!
  afni.header$VOLREG_CENTER_OLD <- c("VOLREG_CENTER_OLD", "float", 4, 3, 3, FALSE)
  afni.header$VOLREG_CENTER_BASE <- c("VOLREG_CENTER_BASE", "float", 4, 3, 3, FALSE)
  afni.header$VOLREG_ROTPARENT_IDCODE <- c("VOLREG_ROTPARENT_IDCODE", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_ROTPARENT_NAME <- c("VOLREG_ROTPARENT_NAME", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_GRIDPARENT_IDCODE <- c("VOLREG_GRIDPARENT_IDCODE", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_GRIDPARENT_NAME <- c("VOLREG_GRIDPARENT_NAME", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_INPUT_IDCODE <- c("VOLREG_INPUT_IDCODE", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_INPUT_NAME <- c("VOLREG_INPUT_NAME", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_BASE_IDCODE <- c("VOLREG_BASE_IDCODE", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_BASE_NAME <- c("VOLREG_BASE_NAME", "string", 4, 1, 1, FALSE)
  afni.header$VOLREG_ROTCOM_NUM <- c("VOLREG_ROTCOM_NUM", "integer", 4, 1, 1, FALSE)

  # Miscellaneous Attributes
  afni.header$IDCOE_ANAT_PARENT <- c("IDCOE_ANAT_PARENT", "string", 5, 1, 1, FALSE)
  afni.header$TO3D_ZPAD <- c("TO3D_ZPAD", "integer", 5, 3, 3, FALSE)

  # Warping Attributes
  afni.header$IDCOE_WARP_PARENT <- c("IDCOE_WARP_PARENT", "string", 6, 1, 1, FALSE)
  afni.header$WARP_TYPE <- c("WARP_TYPE", "integer", 6, 2, 2, FALSE)
  afni.header$WARP_DATA <- c("WARP_DATA", "float", 6, 0, 0, FALSE)

  # Talairach Markers Attributes
  afni.header$MARKS_XYZ <- c("MARKS_XYZ", "float", 7, 30, 30, FALSE)
  afni.header$MARKS_LAB <- c("MARKS_LAB", "string", 7, 1, 1, FALSE)
  afni.header$MARKS_HELP <- c("MARKS_HELP", "string", 7, 1, 1, FALSE)
  afni.header$MARKS_FLAGS <- c("MARKS_FLAGS", "integer", 7, 2, 2, FALSE)

  # Attributes for User-Defined Tags
  afni.header$TAGSET_NUM <- c("TAGSET_NUM", "integer", 8, 2, 2, FALSE)
  afni.header$TAGSET_FLOATS <- c("TAGSET_FLOATS", "floats", 8, 0, 0, FALSE)
  afni.header$TAGSET_LABELS <- c("TAGSET_LABELS", "string", 8, 1, 1, FALSE)

  # Nearly Useless Attributes
  afni.header$LABEL_1 <- c("LABEL_1", "string", 9, 1, 1, FALSE)
  afni.header$LABEL_2 <- c("LABEL_2", "string", 9, 1, 1, FALSE)
  afni.header$DATASET_NAME <- c("DATASET_NAME", "string", 9, 1, 1, FALSE)
  afni.header$DATASET_KEYWORDS <- c("DATASET_KEYWORDS", "string", 9, 1, 1, FALSE)
  afni.header$BRICK_KEYWORDS <- c("BRICK_KEYWORDS", "string", 9, 1, 1, FALSE)

  AFNIheaderpart <- function(name, value, conhead, check=NULL) {
    if (regexpr("_[0-9]*$", name)) if (!is.null(afni.header[[sub("_[0-9]*$","", name)]])) {
      if (as.logical(afni.header[[sub("_[0-9]*$","", name)]][6])) {
        header.entry <- afni.header[[sub("_[0-9]*$","", name)]]
      } else {
        header.entry <- afni.header[[name]]
      }
    } else {
      header.entry <- afni.header[[name]]
    }
    if (!(is.null(header.entry[2]))) {
      # now we know, that entry "name" is valid attribute
      a <- "\n"
      type <- header.entry[2]
      a <- paste(a, "type = ", type, "-attribute\n", sep="")
      a <- paste(a, "name = ", name, "\n", sep="")
      if (is.null(value)) {
        if (header.entry[3] == 0) {
          stop("not found mandatory attribute ",name," in header", call.=FALSE)
        } else if ((header.entry[3] == 1) && taxis) {
          stop("not found mandatory attribute ",name," in header", call.=FALSE)
        } else {
          stop("attempt to write attribute ",name,", which was not found", call.=FALSE)
        }
      }

      if (regexpr("string",type) == 1) {
        if (!is.null(check)) if (!(value %in% check))
          warning("value ", value, " for attribute ", name," does not seem to be valid! Please check!", call.=FALSE)
        value <- paste("'", value, "~", sep="")                         # make syntax highlightening work:'
        a <- paste(a, "count = ", nchar(value) - 1, "\n", sep ="")
        a <- paste(a, value, "\n", sep="")
      } else {
        if (header.entry[4] != 0) {
          if (length(value) < header.entry[4]) {
            warning("too few values for ",name,", expecting ", header.entry[4] ," values, filling with zeros!\n", call.=FALSE)
          }
          length(value) <- as.integer(header.entry[5])
          value[is.na(value)] <- 0
          if (!is.null(check)) {
            length(check) <- as.integer(header.entry[5])
            check[is.na(check)] <- 0
            if (!as.logical(prod((value[1:as.integer(header.entry[4])] == check[1:as.integer(header.entry[4])]))))
              warning("value ", value, " for attribute ", name," does not seem to be valid or match dataset! Please check!", call.=FALSE)
          }
        }
        a <- paste(a, "count = ", length(value), "\n", sep ="")
        j <- 0
        while (j<length(value)) {
          left <- length(value) - j
          if (left>4) left <- 5
          a <- paste(a, paste(value[(j+1):(j+left)],collapse="  "), "\n", sep="  ")
          j <- j+5
        }
      }
      writeChar(a,conhead,eos=NULL)
    } else {
      warning("attempt to write unknown attribute ",name," ...skipping!", call.=FALSE)
    }
  }

  if (!is.null(c(label, note, origin, delta, idcode))) warning("The use of any of the arguments label, note, origin, delta, idcode is depreciated. Please use header and taxis argument, see documentation. THEY WILL VANISH IN SOME FUTURE RELAESE OF THIS SOFTWARE!", call.=FALSE)

  # checking dataset!!!
  if ((length(dim(ttt)) < 3) || (length(dim(ttt)) > 4) || any(dim(ttt)[1:3] == 1)) stop("bad dimension",dim(ttt) ,"for dataset!")
  if (length(dim(ttt)) == 3) dim(ttt) <- c(dim(ttt),1)

  if (is.null(header)) header <- list()

  conhead <- file(paste(filename, ".HEAD", sep=""), "w")
  # write mandatory attributes
  if (is.null(header$DATASET_RANK)) header$DATASET_RANK <- c(3,dim(ttt)[4])
  AFNIheaderpart("DATASET_RANK",header$DATASET_RANK, conhead,c(3,dim(ttt)[4]))
  header$DATASET_RANK <- NULL
  if (is.null(header$DATASET_DIMENSIONS)) header$DATASET_DIMENSIONS <- c(dim(ttt)[1:3])
  AFNIheaderpart("DATASET_DIMENSIONS",header$DATASET_DIMENSIONS, conhead,c(dim(ttt)[1:3]))
  header$DATASET_DIMENSIONS <- NULL
  if (is.null(header$TYPESTRING)) {header$TYPESTRING <- "3DIM_HEAD_FUNC";
          warning("TYPESTRING not given, setting to 3DIM_HEAD_FUNC for backward compatibility", call.=FALSE) }
  AFNIheaderpart("TYPESTRING",header$TYPESTRING, conhead,c("3DIM_HEAD_ANAT","3DIM_HEAD_FUNC","3DIM_GEN_ANAT","3DIM_GEN_ANAT"))
  header$TYPESTRING <- NULL
  if (is.null(header$SCENE_DATA)) {header$SCENE_DATA <- c(0,11,1,-999,-999,-999,-999,-999);
           warning("SCENE_DATA not given, setting to c(0,11,1,-999,-999,-999,-999,-999) for backward compatibility", call.=FALSE) }
  AFNIheaderpart("SCENE_DATA",header$SCENE_DATA, conhead)
  header$SCENE_DATA <- NULL
  if (is.null(header$ORIENT_SPECIFIC)) {header$ORIENT_SPECIFIC <- c(0,3,4);
             warning("ORIENT_SPECIFIC not given, setting to c(0,3,4) for backward compatibility", call.=FALSE) }
  AFNIheaderpart("ORIENT_SPECIFIC",header$ORIENT_SPECIFIC, conhead)
  header$ORIENT_SPECIFIC <- NULL
  if (!is.null(origin)) header$ORIGIN <- origin
  if (is.null(header$ORIGIN)) { header$ORIGIN <- c(0,0,0);
             warning("ORIGIN not given, setting to c(0,0,0) for backward compatibility", call.=FALSE) }
  AFNIheaderpart("ORIGIN",header$ORIGIN, conhead)
  header$ORIGIN <- NULL
  if (!is.null(delta)) header$DELTA <- delta
  if (is.null(header$DELTA)) { header$DELTA <- c(4,4,4);
          warning("DELTA not given, setting to c(4,4,4) for backward compatibility", call.=FALSE) }
  AFNIheaderpart("DELTA",header$DELTA, conhead)
  header$DELTA <- NULL

  # write mandatory attributes for time series
  if (taxis) {
    if (is.null(header$TAXIS_NUMS)) stop("TAXIS_NUMS not given")
    AFNIheaderpart("TAXIS_NUMS",header$TAXIS_NUMS, conhead)
    if (is.null(header$TAXIS_FLOATS)) stop("TAXIS_FLOATS not given")
    AFNIheaderpart("TAXIS_FLOATS",header$TAXIS_FLOATS, conhead)
    header$TAXIS_FLOATS <- NULL
    if (!is.null(header$TAXIS_NUMS[2])) if (header$TAXIS_NUMS[2] != 0) if (is.null(header$TAXIS_OFFSETS)) {
      stop("TAXIS_OFFSETS not given", call.=FALSE) } else { AFNIheaderpart("TAXIS_OFFSETS",header$TAXIS_OFFSETS, conhead); header$TAXIS_OFFSETS <- NULL }
    header$TAXIS_NUMS <- NULL
  }

  # almost mandatory attributes
  if (!is.null(idcode)) header$IDCODE_STRING <- idcode
  if (is.null(header$IDCODE_STRING)) {header$IDCODE_STRING <- "WIAS_noid"; warning("IDCODE_STRING not given, setting to WIAS_noid", call.=FALSE) }
  AFNIheaderpart("IDCODE_STRING",header$IDCODE_STRING, conhead)
  header$IDCODE_STRING <- NULL
  header$IDCODE_DATE <- date()
  AFNIheaderpart("IDCODE_DATE",header$IDCODE_DATE, conhead)
  header$IDCODE_DATE <- NULL
  if (is.null(header$BYTEORDER_STRING)) {
    endian <- .Platform$endian
    header$BYTEORDER_STRING <- switch(endian, "little" = "LSB_FIRST",
                                       "big" = "MSB_FIRST")
  } else {
    if (!(header$BYTEORDER_STRING%in% c("MSB_FIRST","LSB_FIRST"))) header$BYTEORDER <- "MSB_FIRST"
    endian <- switch(header$BYTEORDER_STRING, "MSB_FIRST" = "big",
                                       "LSB_FIRST" = "little")
  }
  AFNIheaderpart("BYTEORDER_STRING",header$BYTEORDER_STRING, conhead,c("MSB_FIRST","LSB_FIRST"))
  header$BYTEORDER_STRING <- NULL
  if (is.null(header$BRICK_STATS)) header$BRICK_STATS <- apply(ttt,4,range)
  dim(header$BRICK_STATS) <- NULL
  AFNIheaderpart("BRICK_STATS",header$BRICK_STATS, conhead)
  if (is.null(header$BRICK_TYPES)) {
    warning("no BRICK_TYPES given, assuming short", call.=FALSE)
    header$BRICK_TYPES <- rep(1,dim(ttt)[4])
  }
  bricktypes <- header$BRICK_TYPES[1]
  if ((bricktypes == 1) && (max(abs(header$BRICK_STATS)) > 32767)) {
    for (k in 1:dim(ttt)[4]) {
      header$BRICK_FLOAT_FACS[k] <- max(abs(header$BRICK_STATS[2*k-1]),abs(header$BRICK_STATS[2*k]))/32767
      ttt[,,,k] <- ttt[,,,k] / header$BRICK_FLOAT_FACS[k]
    }
  }
  header$BRICK_STATS <- NULL
  AFNIheaderpart("BRICK_TYPES",header$BRICK_TYPES, conhead)
  header$BRICK_TYPES <- NULL
  if (is.null(header$BRICK_FLOAT_FACS)) header$BRICK_FLOAT_FACS <- rep(0,dim(ttt)[4])
  AFNIheaderpart("BRICK_FLOAT_FACS",header$BRICK_FLOAT_FACS, conhead)
  header$BRICK_FLOAT_FACS <- NULL
  if (!is.null(label)) header$BRICK_LABS <- paste(label,collapse="~")
  if (!is.null(header$BRICK_LABS)) AFNIheaderpart("BRICK_LABS",header$BRICK_LABS, conhead)
  header$BRICK_LABS <- NULL

  # note attributes
  if (!(is.null(note))) header$HISTORY_NOTE <- paste(header$HISTORY_NOTE,note)
  if (is.null(header$HISTORY_NOTE)) header$HISTORY_NOTE <- ""
  AFNIheaderpart("HISTORY_NOTE",header$HISTORY_NOTE, conhead)
  header$HISTORY_NOTE <- NULL

  for (name in names(header)) {
    AFNIheaderpart(name,header[[name]], conhead)
    header[[name]] <- NULL
  }
  close(conhead)

  if (!(bricktypes %in% c(1,3,5))) stop("Sorry, cannot write this BRICK_TYPES.", call.=FALSE)
  conbrik <- file(paste(filename, ".BRIK", sep=""), "wb")
  dim(ttt) <- NULL
  switch(as.character(bricktypes), "1" = writeBin(as.integer(ttt), conbrik, size=2, endian=endian),
                                   "3" = writeBin(as.numeric(ttt), conbrik, size=4, endian=endian),
                                   "5" = writeBin(as.complex(ttt), conbrik, size=16, endian=endian))
  close(conbrik)
}


read.DICOM <- function(filename, includedata=TRUE) {
  read.DICOM.groupelement <- function(con,endian="little") {
    if (endian == "little") {
      paste(paste(rev(readBin(con,"raw",2,1)),collapse=""),paste(rev(readBin(con,"raw",2,1)),collapse=""),sep=",")
    } else {
      paste(paste(readBin(con,"raw",2,1),collapse=""),paste(readBin(con,"raw",2,1),collapse=""),sep=",")
    }
  }

  con <- file(filename,"rb")

  endian <- "little"

  headerdetails <- list()
  bytes <- 0

  empty <- paste(readBin(con,"character",128,1),collapse="")
  bytes <- bytes + 128
  if (empty != "") {
    warning("This does not seem to be a DICOM file\n")
  }
  dicom <- readChar(con, 4, TRUE)
  bytes <- bytes + 4
  if (dicom != "DICM") {
    warning("This does not seem to be a DICOM file\n")
  }

  groupelement <- 0
  while (TRUE) {
    groupelement <- read.DICOM.groupelement(con)
    bytes <- bytes + 4
    if (groupelement == "7fe0,0010") break
    vr <- readChar(con, 2, TRUE)
    bytes <- bytes + 2

    if (vr %in% c("OB","OW","OF","SQ","UT","UN")) {
      reserved <- readBin(con,"raw",2)
      bytes <- bytes + 2
      length <- readBin(con,"integer",1,4,signed=FALSE,endian=endian)
      bytes <- bytes + 4
      if (length == -1) {
	while (TRUE) {
          groupelement <- read.DICOM.groupelement(con)
          bytes <- bytes + 4
          if (groupelement == "fffe,e0dd") break
          length <- readBin(con,"integer",1,4,signed=FALSE,endian=endian)
          bytes <- bytes + 4
          if (length == -1) {
            sqv <- readBin(con,"raw",4,1)
            bytes <- bytes + 4
            while (TRUE) {
              if (endian == "little") {
                testelement <- paste(paste(rev(sqv[(length(sqv)-3):(length(sqv)-2)]),collapse=""),
                                     paste(rev(sqv[(length(sqv)-1):length(sqv)]),collapse=""),sep=",")
              } else {
                testelement <- paste(paste(sqv[(length(sqv)-3):(length(sqv)-2)],collapse=""),
                                     paste(sqv[(length(sqv)-1):length(sqv)],collapse=""),sep=",")
              }
              if (testelement == "fffe,e00d") {
                length <- readBin(con,"integer",1,4,signed=FALSE,endian=endian)
                bytes <- bytes + 4
                if (length != 0) warning("no valid delimiter in data sequence")
                break
              }
              sqv <- c(sqv,readBin(con,"raw",1,1))
              bytes <- bytes + 1
            }
          } else {
            value <- readBin(con,"raw",length,1)
            bytes <- bytes + length
          }
        }
        length <- readBin(con,"integer",1,4,signed=FALSE,endian=endian)
        bytes <- bytes + 4
        value <- "undecoded sequence"
      }
    } else {
      length <- readBin(con,"integer",1,2,signed=FALSE,endian=endian)
      bytes <- bytes + 2
    }

    if (length != 0) {
      if (vr %in% c("UI","DS","SH","IS")) {
        value <- readChar(con, length, TRUE)
        bytes <- bytes + length
      } else if (vr %in% c("US")) {
        total <- 0
        value <- ""
        while (total < length) {
          value <- paste(value,readBin(con,"integer",1,2,signed=FALSE,endian=endian),sep="")
          bytes <- bytes + 2
          total <- total + 2
        }
      } else {
        value <- readBin(con,"raw",length,1)
        bytes <- bytes + length
      }
    }
#    cat(groupelement,vr,length,paste(value,collapse=""),"\n")
    if (groupelement %in% c("0002,0010",
                            "0018,0050",
                            "0020,0010",
                            "0020,0011",
                            "0020,0012",
                            "0020,0013",
                            "0020,0032",
                            "0020,1041",
                            "0028,0002",
                            "0028,0010",
                            "0028,0011",
                            "0028,0030",
                            "0028,0100",
                            "0028,1050",
                            "0028,1051",
                            "0028,1052",
                            "0028,1053")) {
      headerdetails[[groupelement]] <- value
    }
  }

  # belongs to last groupelement "7fe0,0010"
  vr <- readChar(con, 2, TRUE)
  bytes <- bytes + 2
  if (vr %in% c("OB","OW","OF","SQ","UT","UN")) {
    reserved <- readBin(con,"raw",2)
    bytes <- bytes + 2
    length <- readBin(con,"integer",1,4,signed=FALSE,endian=endian)
    bytes <- bytes + 4
  } else {
    length <- readBin(con,"integer",1,2,signed=FALSE,endian=endian)
    bytes <- bytes + 2
  }
  headerdetails[[groupelement]] <- length
#  cat("Bytes for Header:",bytes,"\n")

  close(con)

  if (!is.null(headerdetails[["0028,0010"]])) {
    xdim <- as.integer(headerdetails[["0028,0010"]])
  } else {
    xdim <- NULL
  }
  if (!is.null(headerdetails[["0028,0011"]])) {
    ydim <- as.integer(headerdetails[["0028,0011"]])
  } else {
    ydim <- NULL
  }
  if (!is.null(headerdetails[["0028,0100"]])) {
    if (headerdetails[["0028,0100"]] == 16) {
      depth <- 2
    } else {
      depth <- 1
    }
  } else {
    depth <- 1
  }
  if (!is.null(headerdetails[["0028,0030"]])) {
    delta <- as.numeric(strsplit(headerdetails[["0028,0030"]],"\\",fixed=TRUE)[[1]])
    if (!is.null(headerdetails[["0018,0050"]])) {
      delta <- c(delta,as.numeric(headerdetails[["0018,0050"]]))
    }
  } else {
    delta <- NULL
  }
  if (!is.null(headerdetails[["0020,0011"]])) {
    series <- as.integer(headerdetails[["0020,0011"]])
  } else {
    series <- NULL
  }
  if (!is.null(headerdetails[["0020,0013"]])) {
    image <- as.integer(headerdetails[["0020,0013"]])
  } else {
    image <- NULL
  }

  # read again
  con <- file(filename,"rb")

  header <- readBin(con,"raw",bytes)

  if (includedata) {
    if (is.null(xdim) || is.null(ydim)) {
      ttt <- readBin(con,"integer",length/depth,depth,signed=FALSE,endian=endian)
      warning("Cannot assign dimension to image because information was not found!")
    } else {
      ttt <- array(readBin(con,"integer",length/depth,depth,signed=FALSE,endian=endian),c(xdim,ydim))
    }

    z <- list(header=headerdetails,ttt=writeBin(as.numeric(ttt),raw(),4),format="DICOM",delta=delta,series=series,image=image,dim=c(xdim,ydim))
  } else {
    z <- list(header=headerdetails,format="DICOM",delta=delta,series=series,image=image,dim=c(xdim,ydim))
  }
  close(con)
#  class(z) <- "fmridata"
#  z <- complete.fmriHeader(z)

  attr(z,"file") <- filename
  invisible(z)
}

read.NIFTI.header <- function(con) {
  header <- list()

  # read 4 bytes and get the endianess by comparing filesize with 348
  endian <- if ((sizeofhdr <- readBin(con,"int",1,4,endian="little")) == 348) "little" else "big"
  header$sizeofhdr <- 348
  header$endian <- endian

  header$datatype1 <- readChar(con,10, TRUE)
  header$dbname <- readChar(con,18, TRUE)
  header$extents <- readBin(con,"int",1,4,endian=endian)
  header$sessionerror <- readBin(con,"int",1,2,endian=endian)
  header$regular <- readChar(con,1, TRUE)
  header$diminfo <- readChar(con,1, TRUE)
  header$dimension <- readBin(con,"int",8,2,endian=endian)
  header$intentp1 <- readBin(con,"double",1,4,endian=endian)
  header$intentp2 <- readBin(con,"double",1,4,endian=endian)
  header$intentp3 <- readBin(con,"double",1,4,endian=endian)
  header$intentcode <- readBin(con,"int",1,2,endian=endian)
  header$datatype <- readBin(con,"int",1,2,endian=endian)
  header$bitpix <- readBin(con,"int",1,2,endian=endian)
  header$slicestart <- readBin(con,"int",1,2,endian=endian)
  header$pixdim <- readBin(con,"double",8,4,endian=endian)
  header$voxoffset <- readBin(con,"double",1,4,endian=endian)
  header$sclslope <- readBin(con,"double",1,4,endian=endian)
  header$sclinter <- readBin(con,"double",1,4,endian=endian)
  header$sliceend <- readBin(con,"int",1,2,endian=endian)
  header$slicecode <- readChar(con,1, TRUE)
  header$xyztunits <- readChar(con,1, TRUE)
  header$calmax <- readBin(con,"double",1,4,endian=endian)
  header$calmin <- readBin(con,"double",1,4,endian=endian)
  header$sliceduration <- readBin(con,"double",1,4,endian=endian)
  header$toffset <- readBin(con,"double",1,4,endian=endian)
  header$glmax <- readBin(con,"int",1,4,endian=endian)
  header$glmin <- readBin(con,"int",1,4,endian=endian)
  header$describ <- readChar(con,80, TRUE)
  header$auxfile <- readChar(con,24, TRUE)
  header$qform <- readBin(con,"int",1,2,endian=endian)
  header$sform <- readBin(con,"int",1,2,endian=endian)
  header$quaternb <- readBin(con,"double",1,4,endian=endian)
  header$quaternc <- readBin(con,"double",1,4,endian=endian)
  header$quaternd <- readBin(con,"double",1,4,endian=endian)
  header$qoffsetx <- readBin(con,"double",1,4,endian=endian)
  header$qoffsety <- readBin(con,"double",1,4,endian=endian)
  header$qoffsetz <- readBin(con,"double",1,4,endian=endian)
  header$srowx <- readBin(con,"double",4,4,endian=endian)
  header$srowy <- readBin(con,"double",4,4,endian=endian)
  header$srowz <- readBin(con,"double",4,4,endian=endian)
  header$intentname <- readChar(con,16, TRUE)
  header$magic <- readChar(con,4, TRUE)

  header
}

read.NIFTI <- function(filename, level=0.75, mask=NULL, setmask=TRUE) {
  fileparts <- strsplit(filename,"\\.")[[1]]
  ext <- tolower(fileparts[length(fileparts)])

  if (ext == "nii") {
    filename.nii <- filename
    filename.hdr <- paste(c(fileparts[-length(fileparts)],"hdr"),collapse=".")
    filename.img <- paste(c(fileparts[-length(fileparts)],"img"),collapse=".")
  } else if (ext == "hdr") {
    filename.hdr <- filename
    filename.img <- paste(c(fileparts[-length(fileparts)],"img"),collapse=".")
  } else if (ext == "img") {
    filename.hdr <- paste(c(fileparts[-length(fileparts)],"hdr"),collapse=".")
    filename.img <- filename
  } else {
    filename.nii <- paste(filename,".nii",sep="")
    filename.hdr <- paste(filename,".hdr",sep="")
    filename.img <- paste(filename,".img",sep="")
  }

  if ((ext != "hdr") && (ext != "img") && (!is.na(file.info(filename.nii)$size))) {
    con <- file(filename.nii,"rb")
    header <- read.NIFTI.header(con)
    if (!(header$magic == "n+1") && !(header$magic == "ni1"))
      warning("Hmmm! Dont see the magic NIFTI string! Try to proceed, but maybe some weird results will occur!");
    bytes <- header$voxoffset - 348
    header$extension <- readBin(con,"raw",bytes)
  } else {
    if (is.na(file.info(filename.hdr)$size) | (file.info(filename.hdr)$size < 348))
      stop("Hmmm! This does not seem to be a NIFTI header (hdr/img-pair)! Wrong size or does not exist!");
    con <- file(filename.hdr,"rb")
    header <- read.NIFTI.header(con)
    header$extension <- NULL
    close(con)
    if (is.na(file.info(filename.img)$size))
      stop("Hmmm! This does not seem to be a NIFTI header (hdr/img-pair)! img-file not found!");
    con <- file(filename.img,"rb")
  }

  dx <- header$dimension[2]
  dy <- header$dimension[3]
  dz <- header$dimension[4]
  dt <- header$dimension[5]
  dd <- if (header$dimension[1] == 5) header$dimension[6] else 1

  endian <- header$endian
  if (header$datatype == 1) { # logical
    what <- "raw"
    signed <- TRUE
    size <- 1
  } else if (header$datatype == 2) { # unsigned char????
    what <- "int"
    signed <- FALSE
    size <- if (header$bitpix) header$bitpix/8/dd else 2
  } else if (header$datatype == 4) { # signed short
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 2
  } else if (header$datatype == 8) { # signed integer
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 4
  } else if (header$datatype == 16) { # float
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 4
  } else if (header$datatype == 32) { # complex
    what <- "complex"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 8
  } else if (header$datatype == 64) { # double
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 8
  } else { # all other
    what <- "raw"
    signed <- TRUE
    size <- 1
  }
  ttt <- readBin(con, what, n=dx*dy*dz*dt*dd, size=size, signed=signed, endian=endian)
  close(con)

  if (min(abs(header$pixdim[2:4])) != 0) {
    weights <-
      abs(header$pixdim[2:4]/min(abs(header$pixdim[2:4])))
  } else {
    weights <- NULL
  }
  dim(ttt) <- if (dd == 1) c(dx,dy,dz,dt) else if (dt == 1) c(dx,dy,dz,dd) else c(dx,dy,dz,dt,dd)

  if (dd == 1) {
    mask0 <- array(TRUE,c(dx,dy,dz))
    if (setmask) {
      mask0[ttt[,,,1] < quantile(ttt[,,,1],level,na.rm = TRUE)] <- FALSE
      dim(ttt) <- c(prod(dim(ttt)[1:3]),dim(ttt)[4])
      na <- ttt %*% rep(1,dim(ttt)[2])
      mask0[is.na(na)] <- FALSE
      ttt[is.na(na),] <- 0
      dim(mask0) <- c(dx,dy,dz)
      mask0 <- connect.mask(mask0)
    }
    z <- list(ttt=writeBin(as.numeric(ttt),raw(),4),
              format="NIFTI",
              delta=header$pixdim[2:4],
              origin=c(header$qoffsetx,header$qoffsety,header$qoffsetz),
              orient=NULL,
              dim=header$dimension[2:5],
              roixa=1,
              roixe=dx,
              roiya=1,
              roiye=dy,
              roiza=1,
              roize=dz,
              roit=1:dt,
              weights=weights,
              header=header,
              mask=mask0)

    class(z) <- "fmridata"
  } else {
    z <- list(ttt=writeBin(as.numeric(ttt),raw(),4),
              format="NIFTI",
              delta=header$pixdim[2:4],
              origin=c(header$qoffsetx,header$qoffsety,header$qoffsetz),
              orient=NULL,
              dim=c(dx,dy,dz,dd),
              roixa=1,
              roixe=dx,
              roiya=1,
              roiye=dy,
              roiza=1,
              roize=dz,
              roit=1:dd,
              weights=weights,
              header=header)
  }

  attr(z,"file") <- paste(filename, sep="")
  z <- complete.fmriHeader(z)
  if(!is.null(mask)) z <- setmask(z,mask)
  invisible(z)
}

complete.fmriHeader <- function(fmriobj){
   header <- fmriobj$header
     if (is.null(header)) header <- list()

  if (!("datatype1" %in% names(header))) header$datatype1 <- paste(rep(" ",10),collapse="")
  if (!("dbname" %in% names(header))) header$dbname <- paste(rep(" ",18),collapse="")
  if (!("extents" %in% names(header))) header$extents <- c(0)
  if (!("sessionerror" %in% names(header))) header$sessionerror <- c(0)
  if (!("regular" %in% names(header))) header$regular <- "r"
  if (!("diminfo" %in% names(header))) header$diminfo <- " "
  if (!("dimension" %in% names(header))) {
    if (is.null(header$dim0))
      header$dimension <- rep(0,8)
    else
      header$dimension <- header$dim0
  }
  if (length(header$dimension) < 8) header$dimension[(length(header$dimension)+1):8] <- 0
  if (!("intentp1" %in% names(header))) header$intentp1 <- 0
  if (!("intentp2" %in% names(header))) header$intentp2 <- 0
  if (!("intentp3" %in% names(header))) header$intentp3 <- 0
  if (!("intentcode" %in% names(header))) header$intentcode <- 0
  if (!("datatype" %in% names(header))) header$datatype <- c(4)
  if (!("bitpix" %in% names(header))) header$bitpix <- c(0)
  if (!("slicestart" %in% names(header))) header$slicestart <- c(0)
  if (!("pixdim" %in% names(header))) header$pixdim <- c(0,4,4,4,rep(0,4))
  if (length(header$pixdim) < 8) header$pixdim[(length(header$pixdim)+1):8] <- 0
  if (!("voxoffset" %in% names(header))||header$voxoffset<348) header$voxoffset <- c(352)
### check if test should be for 348 or 352  ###
  if (!("sclslope" %in% names(header))) header$sclslope <- 0
  if (!("sclinter" %in% names(header))) header$sclinter <- 0
  if (!("sliceend" %in% names(header))) header$sliceend <- 0
  if (!("slicecode" %in% names(header))) header$slicecode <- " "
  if (!("xyztunits" %in% names(header))) header$xyztunits <- " "
  if (!("calmax" %in% names(header))) header$calmax <- c(0)
  if (!("calmin" %in% names(header))) header$calmin <- c(0)
  if (!("sliceduration" %in% names(header))) header$sliceduration <- c(0)
  if (!("toffset" %in% names(header))) header$toffset <- c(0)
  if (!("glmax" %in% names(header))) header$glmax <- c(0)
  if (!("glmin" %in% names(header))) header$glmin <- c(0)
  if (!("describ" %in% names(header))) header$describ <- paste(rep(" ",80),collapse="")
  if (!("auxfile" %in% names(header))) header$auxfile <- paste(rep(" ",24),collapse="")
  if (!("qform" %in% names(header))) header$qform <- 0
  if (!("sform" %in% names(header))) header$sform <- 0
  if (!("quaternb" %in% names(header))) header$quaternb <- 0
  if (!("quaternc" %in% names(header))) header$quaternc <- 0
  if (!("quaternd" %in% names(header))) header$quaternd <- 0
  if (!("qoffsetx" %in% names(header))) header$qoffsetx <- 0
  if (!("qoffsety" %in% names(header))) header$qoffsety <- 0
  if (!("qoffsetz" %in% names(header))) header$qoffsetz <- 0
  if (!("srowx" %in% names(header))) header$srowx <- c(0,0,0,0)
  if (!("srowy" %in% names(header))) header$srowy <- c(0,0,0,0)
  if (!("srowz" %in% names(header))) header$srowz <- c(0,0,0,0)
  if (!("intentname" %in% names(header))) header$intentname <- paste(rep(" ",16),collapse="")
  header$magic <- "n+1"
  if (!("extension" %in% names(header))) header$extension <- as.raw(rep(0,header$voxoffset-348))
  fmriobj$header <- header
  fmriobj
}
write.NIFTI <- function(ttt, header=NULL, filename) {
  if (is.null(header)) header <- list()

  if (!("datatype1" %in% names(header))) header$datatype1 <- paste(rep(" ",10),collapse="")
  if (!("dbname" %in% names(header))) header$dbname <- paste(rep(" ",18),collapse="")
  if (!("extents" %in% names(header))) header$extents <- c(0)
  if (!("sessionerror" %in% names(header))) header$sessionerror <- c(0)
  if (!("regular" %in% names(header))) header$regular <- "r"
  if (!("diminfo" %in% names(header))) header$diminfo <- " "
  if (!("dimension" %in% names(header))) {header$dimension <- rep(0,8); header$dimension <- c(length(dim(ttt)),dim(ttt))}
  if (length(header$dimension) < 8) header$dimension[(length(header$dimension)+1):8] <- 0
  if (!("intentp1" %in% names(header))) header$intentp1 <- 0
  if (!("intentp2" %in% names(header))) header$intentp2 <- 0
  if (!("intentp3" %in% names(header))) header$intentp3 <- 0
  if (!("intentcode" %in% names(header))) header$intentcode <- 0
  if (!("datatype" %in% names(header))) header$datatype <- c(4)
  if (!("bitpix" %in% names(header))) header$bitpix <- c(0)
  if (!("slicestart" %in% names(header))) header$slicestart <- c(0)
  if (!("pixdim" %in% names(header))) header$pixdim <- c(0,4,4,4,rep(0,4))
  if (length(header$pixdim) < 8) header$pixdim[(length(header$pixdim)+1):8] <- 0
  if (!("voxoffset" %in% names(header))) header$voxoffset <- c(352)
  if (!("sclslope" %in% names(header))) header$sclslope <- 0
  if (!("sclinter" %in% names(header))) header$sclinter <- 0
  if (!("sliceend" %in% names(header))) header$sliceend <- 0
  if (!("slicecode" %in% names(header))) header$slicecode <- " "
  if (!("xyztunits" %in% names(header))) header$xyztunits <- " "
  if (!("calmax" %in% names(header))) header$calmax <- c(0)
  if (!("calmin" %in% names(header))) header$calmin <- c(0)
  if (!("sliceduration" %in% names(header))) header$sliceduration <- c(0)
  if (!("toffset" %in% names(header))) header$toffset <- c(0)
  if (!("glmax" %in% names(header))) header$glmax <- c(0)
  if (!("glmin" %in% names(header))) header$glmin <- c(0)
  if (!("describ" %in% names(header))) header$describ <- paste(rep(" ",80),collapse="")
  if (!("auxfile" %in% names(header))) header$auxfile <- paste(rep(" ",24),collapse="")
  if (!("qform" %in% names(header))) header$qform <- 0
  if (!("sform" %in% names(header))) header$sform <- 0
  if (!("quaternb" %in% names(header))) header$quaternb <- 0
  if (!("quaternc" %in% names(header))) header$quaternc <- 0
  if (!("quaternd" %in% names(header))) header$quaternd <- 0
  if (!("qoffsetx" %in% names(header))) header$qoffsetx <- 0
  if (!("qoffsety" %in% names(header))) header$qoffsety <- 0
  if (!("qoffsetz" %in% names(header))) header$qoffsetz <- 0
  if (!("srowx" %in% names(header))) header$srowx <- c(0,0,0,0)
  if (!("srowy" %in% names(header))) header$srowy <- c(0,0,0,0)
  if (!("srowz" %in% names(header))) header$srowz <- c(0,0,0,0)
  if (!("intentname" %in% names(header))) header$intentname <- paste(rep(" ",16),collapse="")
  header$magic <- "n+1"
  if (!("extension" %in% names(header))) header$extension <- as.raw(rep(0,header$voxoffset-348))

  dd <- if (header$dimension[1] == 5) header$dimension[6] else 1
  if (header$datatype == 1) { # logical
    what <- "raw"
    signed <- TRUE
    size <- 1
  } else if (header$datatype == 2) { # unsigned char????
    what <- "int"
    signed <- FALSE
    size <- if (header$bitpix) header$bitpix/8/dd else 2
  } else if (header$datatype == 4) { # signed short
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 2
  } else if (header$datatype == 8) { # signed integer
    what <- "int"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 4
  } else if (header$datatype == 16) { # float
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 4
  } else if (header$datatype == 32) { # complex
    what <- "complex"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 8
  } else if (header$datatype == 64) { # double
    what <- "double"
    signed <- TRUE
    size <- if (header$bitpix) header$bitpix/8/dd else 8
  } else { # all other
    what <- "raw"
    signed <- TRUE
    size <- 1
  }

  fileparts <- strsplit(filename,"\\.")[[1]]

  if (length(fileparts) == 1) {
    filename <- paste(fileparts[1],"nii",sep=".")
  } else {
    ext <- tolower(fileparts[length(fileparts)])
    if (ext != "nii") filename <- paste(c(fileparts,"nii"),collapse=".")
  }
  con <- file(filename, "wb")

  writeBin(as.integer(348), con, 4)
  writeChar(header$datatype1, con, 10, NULL)
  writeChar(header$dbname, con, 18, NULL)
  writeBin(as.integer(header$extents), con, 4)
  writeBin(as.integer(header$sessionerror), con, 2)
  writeChar(header$regular, con, 1, NULL)
  writeChar(header$diminfo, con, 1, NULL)
  writeBin(as.integer(header$dimension), con, 2)
  writeBin(header$intentp1, con, 4)
  writeBin(header$intentp2, con, 4)
  writeBin(header$intentp3, con, 4)
  writeBin(as.integer(header$intentcode), con, 2)
  writeBin(as.integer(header$datatype), con, 2)
  writeBin(as.integer(header$bitpix), con, 2)
  writeBin(as.integer(header$slicestart), con, 2)
  writeBin(header$pixdim, con, 4)
  writeBin(header$voxoffset, con, 4)
  writeBin(header$sclslope, con, 4)
  writeBin(header$sclinter, con, 4)
  writeBin(as.integer(header$sliceend), con, 2)
  writeChar(header$slicecode, con, 1, NULL)
  writeChar(header$xyztunits, con, 1, NULL)
  writeBin(header$calmax, con, 4)
  writeBin(header$calmin, con, 4)
  writeBin(header$sliceduration, con, 4)
  writeBin(header$toffset, con, 4)
  writeBin(as.integer(header$glmax), con, 4)
  writeBin(as.integer(header$glmin), con, 4)
  writeChar(header$describ, con, 80, NULL)
  writeChar(header$auxfile, con, 24, NULL)
  writeBin(as.integer(header$qform), con, 2)
  writeBin(as.integer(header$sform), con, 2)
  writeBin(header$quaternb, con, 4)
  writeBin(header$quaternc, con, 4)
  writeBin(header$quaternd, con, 4)
  writeBin(header$qoffsetx, con, 4)
  writeBin(header$qoffsety, con, 4)
  writeBin(header$qoffsetz, con, 4)
  writeBin(header$srowx, con, 4)
  writeBin(header$srowy, con, 4)
  writeBin(header$srowz, con, 4)
  writeChar(header$intentname, con, 16, NULL)
  writeChar(header$magic, con, 3)
  bytes <- header$voxoffset - 348
  if (bytes != length(header$extension)) warning("header$extension is not of expected size ",bytes," as given by header$voxoffset. cutting!")
  writeBin(header$extension[1:bytes], con)

  dim(ttt) <- NULL
  writeBin(ttt, con, size)
  close(con)
}

extractData <- function(z, what="data", maskOnly=FALSE){
  #
  #  extract data or residuals in either compact (mask only) or expanded (full images) form
  #
    if(is.null(z$maskOnly)) z$maskOnly <- FALSE
  # z$maskOnly was not set, i.e, we have complete data cubes
    expand <- z$maskOnly & !maskOnly
    condense <- !z$maskOnly & maskOnly
    if (z$maskOnly){
       mask <- as.logical(z$mask)
       nvoxel <- sum(mask)
    } else {
       nvoxel <- prod(z$dim[1:3])
       mask <- array(TRUE,z$dim[1:3])
    }
    nt <- z$dim[4]
    n <- nt*nvoxel
    if (what == "residuals") {
        if (!is.null(z$resscale)) {
            ttt <- readBin(z$residuals, "integer", n, 2) * z$resscale
            dim(ttt) <- c(nt,nvoxel)
        }
        else {
            warning("extractData: No residuals available, returning NULL")
            ttt <- NULL
        }
    }
    else {
        if (!is.null(z$ttt)) {
            if (is.null(z$datascale)) {
                ttt <- readBin(z$ttt, "numeric", n, 4)
            }
            else {
                ttt <- readBin(z$ttt, "integer", n, 2) * z$datascale
            }
            dim(ttt) <- c(nvoxel,nt)
        }
        else {
            warning("extractData: No residuals available, returning NULL")
            ttt <- NULL
        }
    }
## now expand to full arrays if needed
    if(!is.null(ttt)&expand){
      ttt0 <- ttt
      if(what=="data"){
         ttt <- matrix(0,n/nt,nt)
         ttt[mask,] <- ttt0
      } else {
        ttt <- matrix(0,nt,n/nt)
        ttt[,mask] <- ttt0
      }
    }
## now condense to mask if required
    if(!is.null(ttt)&condense){
      if(what=="data"){
        ttt <- ttt[mask,]
      } else {
        ttt <- ttt[,mask]
      }
    }
    if(!is.null(ttt)&!maskOnly){
      if(what=="data") ind <- 1:4 else ind <- c(4,1:3)
      dim(ttt) <- z$dim[ind]
    }
    invisible(ttt)
}

expandfMRI <- function(z){
# expand fmridataobj to contain full images
if(!is.null(z$maskOnly)&z$maskOnly){
  mask <- z$mask
  nvoxel <- sum(mask)
  nt <- z$dim[4]
  n <- nt*nvoxel
  if(!is.null(z$residuals)){
    ttt0 <- readBin(z$residuals, "integer", n, 2)
    ttt <- matrix(0L, nt, nvoxel)
    ttt[,mask] <- ttt
    z$residuals <- writeBin(as.integer(ttt),raw(),2)
  }
  if(!is.null(z$ttt)){
#    object contains data
     if (is.null(z$datascale)) {
         ttt0 <- readBin(z$ttt, "numeric", n, 4)
     }
     else {
         ttt0 <- readBin(z$ttt, "integer", n, 2) * z$datascale
     }
     datascale <- max(abs(range(ttt)))/32767
     ttt <- matrix(0L, nvoxel, nt)
     ttt[mask,] <- ttt
     z$ttt <- writeBin(as.integer(ttt/datascale),raw(),2)
     z$datascale <- datascale
  }
}
z$maskOnly <- FALSE
invisible(z)
}

condensefMRI <- function(z, mask=NULL){
#
#  reduce size of fmri object to voxel in mask
#  can also be used to set a more restrictive mask
#
  if(is.null(z$maskOnly)) z$maskOnly <- FALSE
  if(is.null(z$mask)) z$mask <- array(TRUE, z$dim[1:3])
  if(z$maskOnly&is.null(mask)){
    warning("condensefMRI: nothing to do\n returning unchanged object")
    invisible(z)
  }
  if(z$maskOnly) z <- expandfMRI(z)
  if(!is.null(z$mask)&!is.null(mask)){
    if(any(mask&!z$mask))
    warning("condensefMRI: new mask is not a subset of old mask, using intesection")
    mask[!z$mask] <- FALSE
  }
  if(is.null(mask)) mask <- z$mask
  if(is.null(mask)){
    warning("condensefMRI: no mask available\n returning unchanged object")
    invisible(z)
  }
  z$mask <- mask
  nt <- z$dim[4]
  ncube <- prod(z$dim[1:3])
  n <- nt*ncube
  nvoxel <- sum(mask)
  if(!is.null(z$residuals)){
    ttt <- matrix(readBin(z$residuals, "integer", n, 2), nt, ncube)
    z$residuals <- writeBin(as.integer(ttt[,mask]),raw(),2)
  }
  if(!is.null(z$ttt)){
#    object contains data
     if (is.null(z$datascale)) {
         ttt <- readBin(z$ttt, "numeric", n, 4)
     }
     else {
         ttt <- readBin(z$ttt, "integer", n, 2) * z$datascale
     }
     datascale <- max(abs(range(ttt)))/32767
     dim(ttt) <- c(ncube, nt)
     z$ttt <- writeBin(as.integer(ttt[mask,]/datascale),raw(),2)
     z$datascale <- datascale
  }
  z$maskOnly <- TRUE
  invisible(z)
}