The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Reading and writing images

library(neuroim2)

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")

Choosing a reader

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  4

Compression 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:

dim(read_vol(series_file, index = 3))
#> [1] 64 64 25

Use read_image() when you do not know what a file holds:

class(read_image(series_file))[1]
#> [1] "DenseNeuroVec"

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  8

Neither returns a list, despite the name — length() on either counts volumes.

Looking before you load

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.7

dim() 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.00

spacing(), 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:

hdr@header$datatype
#> [1] 2
hdr@header$encoding
#> [1] "gzip"
hdr@header$vox_offset
#> [1] 3040

The two affines

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] TRUE

They do not always agree. The anatomical image shipped with this package is a real example — its qform describes 4 mm voxels and its sform describes 1 mm ones:

anat_hdr <- read_header(anat_file)

anat_hdr@header$pixdim[2:4]
#> [1] 4.02083 4.01754 4.02083
anat_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
anat_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    1

Both codes are set, so the sform wins and the image reports 1 mm voxels even though pixdim says otherwise:

spacing(read_vol(anat_file))
#> [1] 1 1 1

Nothing errors. The image loads and plots correctly, and every distance, volume and ROI radius derived from it is wrong by a factor of about four. 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 roughly 140 mm at the far corner:

far <- matrix(dim(read_vol(anat_file)), nrow = 1)

sform_xyz <- grid_to_coord(space(read_vol(anat_file)), far)
qform_xyz <- (anat_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.979

This 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))
#>  mask  anat 
#>  TRUE FALSE

The is.null() guard matters: not every format carries these fields, and returning NA for “cannot tell” is better than erroring on an AFNI file.

If a file is affected and you trust its pixdim, rebuild the geometry:

bad <- read_vol(anat_file)
fixed <- NeuroVol(
  as.array(bad),
  NeuroSpace(dim(bad),
    spacing = anat_hdr@header$pixdim[2:4],
    origin = anat_hdr@header$qform[1:3, 4],
    axes = space(bad)@axes
  )
)

spacing(fixed)
#> [1] 4.02083 4.01754 4.02083

Writing

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    2157

That 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] TRUE

4D 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  2

A precision surprise that is not about files

Affines 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.141590

No file was involved. Because both sides of a round trip get the same rounding, identical() on a re-read affine does succeed, as above. What you cannot rely on is an affine you constructed matching the full-precision numbers you passed in.

AFNI files

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:

vol <- read_vol("subject_anat+orig.HEAD")
names(read_header("subject_anat+orig.HEAD")@header)

Where to go next

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.