Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.