Nothing
# TODO test MEGA-PRESS and CSI
# when nechoes = 2 we're probably dealing with MEGA-PRESS data, where edited
# pairs are not interleaved but occupy the first and second half of the data
read_pfile <- function(fname, n_ref_scans = NULL, verbose, extra) {
# check the file size
fbytes <- file.size(fname)
hdr <- read_pfile_header(fname)
# override w ref scans detected in the header
if (!is.null(n_ref_scans)) hdr$rhuser19 <- n_ref_scans
con <- file(fname, "rb")
seek(con, hdr$off_data)
endian <- "little"
# calculate number of data points from file size
Npts <- (fbytes - hdr$off_data) / 4
raw_pts <- readBin(con, "int", n = Npts, size = 4, endian = endian)
close(con)
# calculate coil elements
coils <- 0
for (n in seq(1, 8, 2)) {
if ((hdr$rcv[n] != 0) || (hdr$rcv[n + 1] != 0)) {
coils <- coils + 1 + hdr$rcv[n + 1] - hdr$rcv[n]
}
}
if (coils == 0) coils <- 1
# each echo starts with a frame of zeros followed by w_ref frames followed by
# ws frames
expt_pts <- coils * (hdr$nframes * hdr$nechoes + hdr$nechoes) * hdr$frame_size * 2
if (expt_pts != Npts) {
warning("Unexpected number of data points. Check for data corruption due to autocoil selection bug.")
cat(paste("Expecting :", Npts, "points based on file size.\n"))
cat(paste("Expecting :", expt_pts, "points based on header information.\n"))
cat(paste("Coils :", coils, "\n"))
cat(paste("nframes :", hdr$nframes, "\n"))
cat(paste("nechoes :", hdr$nechoes), "\n")
cat(paste("frame_size :", hdr$frame_size), "\n")
cat(paste("w_frames :", hdr$rhuser19), "\n")
cat(paste("Header rev. :", hdr$hdr_rev), "\n")
cat(paste("Seq. name :", hdr$seq_name), "\n")
cat(paste("Prot. name :", hdr$prot_name), "\n")
cat(paste("CSI dims :", hdr$csi_dims), "\n")
cat(paste("xcsi :", hdr$xcsi), "\n")
cat(paste("ycsi :", hdr$ycsi), "\n")
cat(paste("zcsi :", hdr$zcsi), "\n")
cat(paste("Data offset :", hdr$off_data), "\n")
}
if (verbose) {
cat(paste("Expecting :", Npts, "points based on file size.\n"))
cat(paste("Expecting :", expt_pts, "points based on header information.\n"))
cat(paste("Coils :", coils, "\n"))
cat(paste("nframes :", hdr$nframes, "\n"))
cat(paste("nechoes :", hdr$nechoes), "\n")
cat(paste("frame_size :", hdr$frame_size), "\n")
cat(paste("w_frames :", hdr$rhuser19), "\n")
cat(paste("Header rev. :", hdr$hdr_rev), "\n")
cat(paste("Seq. name :", hdr$seq_name), "\n")
cat(paste("Prot. name :", hdr$prot_name), "\n")
cat(paste("CSI dims :", hdr$csi_dims), "\n")
cat(paste("xcsi :", hdr$xcsi), "\n")
cat(paste("ycsi :", hdr$ycsi), "\n")
cat(paste("zcsi :", hdr$zcsi), "\n\n")
}
data <- raw_pts[c(TRUE, FALSE)] + 1i * raw_pts[c(FALSE, TRUE)]
if (hdr$csi_dims == 0) {
# looks like SVS
dyns <- hdr$nechoes * hdr$nframes + hdr$nechoes
data <- array(data, dim = c(hdr$frame_size, dyns, coils, 1, 1, 1, 1))
data <- aperm(data, c(7, 6, 5, 4, 2, 3, 1))
# remove the empty frame at the start of each coil for each echo
rem <- seq(from = 1, to = dyns, by = dyns / hdr$nechoes)
data <- data[,,,,-rem,,,drop = FALSE]
} else {
# looks like MRSI
# the + 1 is an empty FID at the start of each coil
data <- array(data, dim = c(hdr$frame_size, hdr$xcsi * hdr$ycsi + 1, coils,
1, 1, 1, 1))
data <- aperm(data, c(7, 6, 5, 4, 2, 3, 1))
# remove the empty FID
data <- data[,,,,-1,,,drop = FALSE]
# reshape
new_dim <- dim(data)
new_dim[2] <- hdr$xcsi
new_dim[3] <- hdr$ycsi
new_dim[5] <- 1
dim(data) <- new_dim
# phase cycling correction?
k_sp_corr_x <- array(1, dim(data))
k_sp_corr_x[, c(T, F),,,,,] <- -1
data <- data * k_sp_corr_x
}
res <- c(NA, NA, NA, NA, hdr$tr, NA, 1 / hdr$spec_width)
# freq domain vector vector
freq_domain <- rep(FALSE, 7)
ref <- def_ref()
nuc <- def_nuc()
meta <- list(EchoTime = hdr$te,
SequenceName = hdr$seq_name,
ProtocolName = hdr$prot_name,
Manufacturer = "GE")
if (toupper(hdr$seq_name) == "PROBE-P") {
meta <- append(meta, list(PulseSequenceType = "press"))
}
if (toupper(hdr$seq_name) == "PROBE-S") {
meta <- append(meta, list(PulseSequenceType = "steam"))
}
mrs_data <- mrs_data(data = data, ft = hdr$ps_mps_freq / 10, resolution = res,
ref = ref, nuc = nuc, freq_domain = freq_domain,
affine = NULL, meta = meta, extra = extra)
if (hdr$rhuser19 > 0) {
# split water and metab data for each echo
wref_inds <- rep(FALSE, Ndyns(mrs_data) / hdr$nechoes)
wref_inds[1:hdr$rhuser] <- TRUE
wref_inds <- rep(wref_inds, hdr$nechoes)
ref_mrs <- get_dyns(mrs_data, which(wref_inds))
metab_mrs <- get_dyns(mrs_data, which(!wref_inds))
if (hdr$csi_dims > 0) {
ref_mrs <- img2kspace_xy(ref_mrs)
metab_mrs <- img2kspace_xy(metab_mrs)
}
ref_mrs$meta <- append(ref_mrs$meta,
list(NumberOfTransients = Ndyns(ref_mrs)))
} else {
ref_mrs <- NA
metab_mrs <- mrs_data
if (hdr$csi_dims > 0) {
metab_mrs <- img2kspace_xy(metab_mrs)
}
}
metab_mrs$meta <- append(metab_mrs$meta,
list(NumberOfTransients = Ndyns(metab_mrs)))
out <- list(metab = metab_mrs, ref = ref_mrs)
class(out) <- c("list", "mrs_data")
return(out)
}
read_pfile_header <- function(fname) {
endian <- "little"
vars <- get_pfile_vars()
con <- file(fname, "rb")
vars$hdr_rev <- readBin(con, "numeric", size = 4, endian = endian)
loc <- get_pfile_dict(vars$hdr_rev, con)
seek(con, loc$off_data)
vars$off_data <- readBin(con, "int", size = 4, endian = endian)
seek(con, loc$nechoes)
vars$nechoes <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$nframes)
vars$nframes <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$frame_size)
vars$frame_size <- readBin(con, "int", size = 2, signed = FALSE,
endian = endian)
seek(con, loc$rcv)
vars$rcv <- readBin(con, "int", n = 8, size = 2, endian = endian)
seek(con, loc$rhuser19)
vars$rhuser19 <- readBin(con, "numeric", size = 4, endian = endian)
seek(con, loc$spec_width)
vars$spec_width <- readBin(con, "numeric", size = 4, endian = endian)
seek(con, loc$csi_dims)
vars$csi_dims <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$xcsi)
vars$xcsi <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$ycsi)
vars$ycsi <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$zcsi)
vars$zcsi <- readBin(con, "int", size = 2, endian = endian)
seek(con, loc$ps_mps_freq)
# read as int
ps_mps_freq_bits <- intToBits(readBin(con, "int", size = 4, endian = endian))
# convert to uint
vars$ps_mps_freq <- sum(2^.subset(0:31, as.logical(ps_mps_freq_bits)))
seek(con, loc$te)
vars$te <- readBin(con, "int", size = 4, endian = endian) / 1e6
seek(con, loc$tr)
vars$tr <- readBin(con, "int", size = 4, endian = endian) / 1e6
seek(con, loc$seq_name)
vars$seq_name <- readBin(con, "character", size = 33)
seek(con, loc$prot_name)
vars$prot_name <- readBin(con, "character", size = 65)
close(con)
vars
}
get_pfile_vars <- function() {
vars <- vector(mode = "list", length = 17)
names(vars) <- c("hdr_rev", "off_data", "nechoes", "nframes", "frame_size",
"rcv", "rhuser19", "spec_width", "csi_dims", "xcsi", "ycsi",
"zcsi", "ps_mps_freq", "te", "tr", "seq_name", "prot_name")
vars
}
# following are useful resources for GE p-files
# https://github.com/SIVICLab/sivic/blob/master/libs/src/svkGEPFileReader.cc
# https://github.com/chenkonturek/MRS_MRI_libs/blob/master/MRS_lib/io/mrs_readGEpfile.m
# TARQUIN source is an ok reference for older versions.
# Ralph Noeske also suggests GERead.m from GANNET as a good place to look for
# newer versions.
# Other option is to use a hex editor, and search for known values
# https://hexed.it/
get_pfile_dict <- function(hdr_rev, con) {
loc <- get_pfile_vars()
if (floor(hdr_rev) == 11L) { # version 12
warning("pfile version is untested :", hdr_rev)
loc$hdr_rev <- 0
loc$off_data <- 4
loc$nechoes <- 146
loc$nframes <- 150
loc$frame_size <- 156
loc$rcv <- 264
loc$rhuser19 <- 356
loc$spec_width <- 432
loc$csi_dims <- 436
loc$xcsi <- 438
loc$ycsi <- 440
loc$zcsi <- 442
loc$ps_mps_freq <- 488
loc$te <- 1148
loc$tr <- 199236
loc$seq_name <- 199812
loc$prot_name <- 199812 # TODO
} else if (floor(hdr_rev) == 14L) {
loc$hdr_rev <- 0
loc$off_data <- 1468
loc$nechoes <- 70
loc$nframes <- 74
loc$frame_size <- 80
loc$rcv <- 200
loc$rhuser19 <- 292
loc$spec_width <- 368
loc$csi_dims <- 372
loc$xcsi <- 374
loc$ycsi <- 376
loc$zcsi <- 378
loc$ps_mps_freq <- 424
loc$te <- 144580
loc$tr <- 144572
loc$seq_name <- 145132
loc$prot_name <- 142318
} else if (floor(hdr_rev) == 20L) {
loc$hdr_rev <- 0
loc$off_data <- 1468
loc$nechoes <- 70
loc$nframes <- 74
loc$frame_size <- 80
loc$rcv <- 200
loc$rhuser19 <- 292
loc$spec_width <- 368
loc$csi_dims <- 372
loc$xcsi <- 374
loc$ycsi <- 376
loc$zcsi <- 378
loc$ps_mps_freq <- 424
loc$te <- 1212
loc$tr <- 148396
loc$seq_name <- 148972
loc$prot_name <- 145762
} else if (floor(hdr_rev) == 21L) {
loc$hdr_rev <- 0
loc$off_data <- 1468
loc$nechoes <- 70
loc$nframes <- 74
loc$frame_size <- 80
loc$rcv <- 200
loc$rhuser19 <- 292
loc$spec_width <- 368
loc$csi_dims <- 372
loc$xcsi <- 374
loc$ycsi <- 376
loc$zcsi <- 378
loc$ps_mps_freq <- 424
loc$te <- 1212
loc$tr <- 148944
loc$seq_name <- 149520
loc$prot_name <- 146310
} else if (floor(hdr_rev) == 24L) {
loc$hdr_rev <- 0
loc$off_data <- 1468
loc$nechoes <- 70
loc$nframes <- 74
loc$frame_size <- 80
loc$rcv <- 200
loc$rhuser19 <- 292
loc$spec_width <- 368
loc$csi_dims <- 372
loc$xcsi <- 374
loc$ycsi <- 376
loc$zcsi <- 378
loc$ps_mps_freq <- 424
# loc$te <- 1212
loc$te <- 148404
loc$tr <- 148396
loc$seq_name <- 148972
loc$prot_name <- 145762
} else if (floor(hdr_rev) == 26L) {
loc$hdr_rev <- 0
loc$off_data <- 4
loc$nechoes <- 146
loc$nframes <- 150
loc$frame_size <- 156
loc$rcv <- 264
loc$rhuser19 <- 356
loc$spec_width <- 432
loc$csi_dims <- 436
loc$xcsi <- 438
loc$ycsi <- 440
loc$zcsi <- 442
loc$ps_mps_freq <- 488
loc$te <- 1148
loc$tr <- 199236
loc$seq_name <- 199812
loc$prot_name <- 196602
} else if ((floor(hdr_rev) == 28L) | (floor(hdr_rev) == 30L)) {
# some files from this version are corrupt due to a bug in the GE software
# (relating to auto-coil selection "Air Touch")
# don't spend ages trying to fix files that are inherently borked
loc$hdr_rev <- 0
loc$off_data <- 4
loc$nechoes <- 146
loc$nframes <- 150
loc$frame_size <- 156
loc$rcv <- 264
loc$rhuser19 <- 356
loc$spec_width <- 432
loc$csi_dims <- 436
loc$xcsi <- 438
loc$ycsi <- 440
loc$zcsi <- 442
loc$ps_mps_freq <- 488
loc$te <- 1148
loc$tr <- 207428
loc$seq_name <- 208004
loc$prot_name <- 204794
} else {
warning("pfile version is untested :", hdr_rev)
loc$hdr_rev <- 0
loc$off_data <- 4
loc$nechoes <- 146
loc$nframes <- 150
loc$frame_size <- 156
loc$rcv <- 264
loc$rhuser19 <- 356
loc$spec_width <- 432
loc$csi_dims <- 436
loc$xcsi <- 438
loc$ycsi <- 440
loc$zcsi <- 442
loc$ps_mps_freq <- 488
loc$te <- 1148
loc$tr <- 207428
loc$seq_name <- 208004
loc$prot_name <- 204794
}
return(loc)
}
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.