Every analysis starts and ends at a file. This article covers the readers and writers, how to inspect a header before committing to loading anything, and the one header problem that produces wrong results without producing an error.
mask_file <- system.file("extdata", "global_mask2.nii.gz", package = "neuroim2")
series_file <- system.file("extdata", "global_mask_v4.nii", package = "neuroim2")
anat_file <- system.file("extdata", "mni_downsampled.nii.gz", package = "neuroim2")| Function | Returns | Use for |
|---|---|---|
read_vol() |
NeuroVol |
a 3D image, or one volume of a 4D one |
read_vec() |
NeuroVec |
a 4D series |
read_image() |
whichever fits | unknown dimensionality |
read_vol_list() |
NeuroVec |
many 3D files as one series |
read_header() |
FileMetaInfo |
metadata without the data |
vol <- read_vol(mask_file)
vec <- read_vec(series_file)
dim(vol)
#> [1] 64 64 25
dim(vec)
#> [1] 64 64 25 4Compression is handled by extension: .nii and
.nii.gz both work on read and on write.
read_vol() on a 4D file returns its first volume, and
index picks a different one:
Use read_image() when you do not know what a file
holds:
Reading several files gives you a 4D object either way, but not the
same one. read_vol_list() builds a single contiguous
series; read_vec() on a vector of paths keeps the runs as
separable segments:
class(read_vol_list(c(mask_file, mask_file)))[1]
#> [1] "DenseNeuroVec"
class(read_vec(c(series_file, series_file)))[1]
#> [1] "NeuroVecSeq"
dim(read_vec(c(series_file, series_file)))
#> [1] 64 64 25 8Neither returns a list, despite the name — length() on
either counts volumes.
read_header() reads metadata only. On a large file this
costs nothing and tells you what you are about to commit to:
hdr <- read_header(mask_file)
dim(hdr)
#> [1] 64 64 25
hdr@data_type
#> [1] "UBYTE"
hdr@spacing
#> [1] 3.5 3.5 3.7dim() and trans() work on a header, which
is the quickest way to see the affine neuroim2 will actually use without
loading a byte of data:
trans(hdr)
#> [,1] [,2] [,3] [,4]
#> [1,] -3.5 0.0 0.0 112.00
#> [2,] 0.0 3.5 0.0 -108.50
#> [3,] 0.0 0.0 3.7 -46.25
#> [4,] 0.0 0.0 0.0 1.00spacing(), origin() and
space() are not defined for headers — read the slots
(hdr@spacing, hdr@origin) instead. The full
NIfTI header is a plain list in the header slot, so
anything the slots omit is still reachable:
A NIfTI file can carry two spatial transforms: the
qform (a rigid transform stored as a quaternion) and the
sform (a general affine). Each has a code saying whether it
is set. The standard says use the sform when
sform_code > 0, otherwise the qform.
On a well-formed file the two agree:
c(qform_code = hdr@header$qform_code, sform_code = hdr@header$sform_code)
#> qform_code sform_code
#> 1 1
all.equal(hdr@header$qform, hdr@header$sform)
#> [1] TRUEThey do not always agree. The anatomical image shipped with this package was a real example in neuroim2 0.19.1 and earlier: its qform described 4 mm voxels while its sform had lost the scaling and described 1 mm ones. Recreate that fault in a temporary copy:
bad_file <- tempfile(fileext = ".nii.gz")
img <- RNifti::readNifti(anat_file)
broken <- matrix(RNifti::xform(img, useQuaternionFirst = FALSE), 4, 4)
broken[1:3, 1:3] <- diag(3) # identity scaling, translation kept
RNifti::sform(img) <- structure(broken, code = 1L)
RNifti::writeNifti(img, bad_file)
bad_hdr <- read_header(bad_file)
bad_hdr@header$pixdim[2:4]
#> [1] 4.02083 4.01754 4.02083
bad_hdr@header$qform
#> [,1] [,2] [,3] [,4]
#> [1,] 4.02083 0.00000 0.00000 -96
#> [2,] 0.00000 4.01754 0.00000 -132
#> [3,] 0.00000 0.00000 4.02083 -78
#> [4,] 0.00000 0.00000 0.00000 1
bad_hdr@header$sform
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 0 -96
#> [2,] 0 1 0 -132
#> [3,] 0 0 1 -78
#> [4,] 0 0 0 1Both codes are set, so the sform wins and the image reports 1 mm
voxels even though pixdim says otherwise:
spacing(read_vol(bad_file))
#> [1] 1 1 1
spacing(read_vol(anat_file))
#> [1] 4.02083 4.01754 4.02083Nothing errors. The image loads and plots correctly, and every
distance, volume and ROI radius derived from it uses
the wrong physical scale. Lengths differ by a factor of about four and
voxel volumes by about 64. The coordinates themselves are not scaled by
four — both affines share a translation, so they agree exactly at voxel
(1, 1, 1) and diverge by a growing offset that reaches 140
to 170 mm per axis at the far corner:
far <- matrix(dim(read_vol(bad_file)), nrow = 1)
sform_xyz <- grid_to_coord(space(read_vol(bad_file)), far)
qform_xyz <- (bad_hdr@header$qform %*% c(far - 1, 1))[1:3]
rbind(sform = as.vector(sform_xyz), qform = qform_xyz)
#> [,1] [,2] [,3]
#> sform -49.00000 -76.00000 -31.000
#> qform 92.97902 92.98224 110.979This is the failure mode to know about, because no amount of care downstream will catch it. Check when a file first arrives:
affines_agree <- function(file) {
h <- read_header(file)
qc <- h@header$qform_code
sc <- h@header$sform_code
if (is.null(qc) || is.null(sc) || qc <= 0 || sc <= 0) return(NA)
isTRUE(all.equal(h@header$qform, h@header$sform, tolerance = 1e-4))
}
c(mask = affines_agree(mask_file), anat = affines_agree(anat_file),
broken = affines_agree(bad_file))
#> mask anat broken
#> TRUE TRUE FALSEThe is.null() guard matters: not every format carries
these fields, and returning NA for “cannot tell” is better
than erroring on an AFNI file.
Disagreement alone does not tell you which transform is correct; qform and sform may intentionally describe different coordinate spaces. Check the image’s provenance before changing either. If you have established that the qform is the correct geometry, use its full affine so rotation and axis directions survive:
write_vol() and write_vec() take an object
and a path. The extension decides the container and the compression; an
unrecognised or absent extension gives uncompressed NIfTI.
out_nii <- tempfile(fileext = ".nii")
out_gz <- tempfile(fileext = ".nii.gz")
write_vol(vol, out_nii)
write_vol(vol, out_gz)
c(plain = file.size(out_nii), gzipped = file.size(out_gz))
#> plain gzipped
#> 102752 2157That plain file is 400 kB for a binary mask, because
write_vol() promotes to FLOAT by default. Pass
data_type when you know better:
out_byte <- tempfile(fileext = ".nii")
write_vol(vol, out_byte, data_type = "UBYTE")
c(float = file.size(out_nii), ubyte = file.size(out_byte))
#> float ubyte
#> 102752 102752
read_header(out_byte)@data_type
#> [1] "UBYTE"Valid types are UNKNOWN, BINARY,
UBYTE, SHORT, INT,
FLOAT and DOUBLE. A round trip preserves both
data and geometry either way:
back <- read_vol(out_gz)
all.equal(as.array(back), as.array(vol))
#> [1] TRUE
identical(trans(space(back)), trans(space(vol)))
#> [1] TRUE4D objects work the same way:
out_vec <- tempfile(fileext = ".nii.gz")
write_vec(sub_vector(vec, 1:2), out_vec)
dim(read_vec(out_vec))
#> [1] 64 64 25 2Affines round to seven significant figures, and spacing and origin to
six — but this happens in NeuroSpace(), when the object is
built, not when it is written:
sp <- NeuroSpace(c(4L, 4L, 4L), spacing = c(1 / 3, 0.123456789012, pi))
spacing(sp)
#> [1] 0.333333 0.123457 3.141590No file was involved. The mask’s affine above survives this round
trip exactly, but that is not a guarantee for every affine: NIfTI-1 also
stores transform coefficients as 32-bit floats. Use
all.equal() with an appropriate tolerance when checking
general round trips. An affine you construct may also differ from the
full-precision numbers you passed in.
AFNI HEAD/BRIK pairs are read by the same
functions — pass the .HEAD file and the matching
.BRIK is found beside it. read_header()
returns an AFNIMetaInfo whose header slot
holds the AFNI attribute list:
vignette("spaces-and-coordinates") — what the affine
you just read meansvignette("large-data") — reading files too big to hold
in memory?read_vol, ?read_header,
?write_vol