---
title: Reading and writing images
output:
  rmarkdown::html_vignette:
    css: albers.css
    includes:
      in_header:
      - albers-header.inc
      - albers-header.html
    toc: yes
    toc_depth: 2.0
resource_files:
- albers.css
- albers-fonts.css
- albers.js
- albers-header.inc
- albers-header.html
- fonts
params:
  family: red
  preset: interaction

vignette: |
  %\VignetteIndexEntry{Reading and writing images}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
source("_common.R", local = TRUE)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

```{r load}
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.

```{r paths}
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 |

```{r readers}
vol <- read_vol(mask_file)
vec <- read_vec(series_file)

dim(vol)
dim(vec)
```

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:

```{r read-index}
dim(read_vol(series_file, index = 3))
```

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

```{r read-image}
class(read_image(series_file))[1]
```

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:

```{r multi}
class(read_vol_list(c(mask_file, mask_file)))[1]
class(read_vec(c(series_file, series_file)))[1]

dim(read_vec(c(series_file, series_file)))
```

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:

```{r header}
hdr <- read_header(mask_file)

dim(hdr)
hdr@data_type
hdr@spacing
```

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

```{r header-trans}
trans(hdr)
```

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

```{r header-slot}
hdr@header$datatype
hdr@header$encoding
hdr@header$vox_offset
```

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

```{r affines-good}
c(qform_code = hdr@header$qform_code, sform_code = hdr@header$sform_code)

all.equal(hdr@header$qform, hdr@header$sform)
```

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

```{r affines-bad}
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]
bad_hdr@header$qform
bad_hdr@header$sform
```

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

```{r affines-consequence}
spacing(read_vol(bad_file))
spacing(read_vol(anat_file))
```

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

```{r affines-divergence}
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)
```

**This is the failure mode to know about**, because no amount of care downstream
will catch it. Check when a file first arrives:

```{r affine-check}
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))
```

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.

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:

```{r repair}
bad <- read_vol(bad_file)
fixed <- NeuroVol(
  as.array(bad),
  NeuroSpace(dim(bad), trans = bad_hdr@header$qform)
)

spacing(fixed)
all.equal(trans(space(fixed)), trans(space(read_vol(anat_file))),
          check.attributes = FALSE)
```

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

```{r write}
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))
```

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:

```{r datatype}
out_byte <- tempfile(fileext = ".nii")
write_vol(vol, out_byte, data_type = "UBYTE")

c(float = file.size(out_nii), ubyte = file.size(out_byte))
read_header(out_byte)@data_type
```

Valid types are `UNKNOWN`, `BINARY`, `UBYTE`, `SHORT`, `INT`, `FLOAT` and
`DOUBLE`. A round trip preserves both data and geometry either way:

```{r roundtrip}
back <- read_vol(out_gz)

all.equal(as.array(back), as.array(vol))
identical(trans(space(back)), trans(space(vol)))
```

4D objects work the same way:

```{r write-vec}
out_vec <- tempfile(fileext = ".nii.gz")
write_vec(sub_vector(vec, 1:2), out_vec)

dim(read_vec(out_vec))
```

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

```{r precision}
sp <- NeuroSpace(c(4L, 4L, 4L), spacing = c(1 / 3, 0.123456789012, pi))

spacing(sp)
```

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

```{r cleanup, include = FALSE}
unlink(c(out_nii, out_gz, out_byte, out_vec))
```

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

```{r afni, eval = FALSE}
vol <- read_vol("subject_anat+orig.HEAD")
names(read_header("subject_anat+orig.HEAD")@header)
```

## Where to go next

- `vignette("spaces-and-coordinates")` — what the affine you just read means
- `vignette("large-data")` — reading files too big to hold in memory
- `?read_vol`, `?read_header`, `?write_vol`