Reading a NIFTI formatted image volume
The way to read an volumetric image file is to use
read_vol:
library(neuroim2)
#> Loading required package: Matrix
#>
#> Attaching package: 'neuroim2'
#> The following object is masked from 'package:base':
#>
#> scale
file_name <- system.file("extdata", "global_mask2.nii.gz", package="neuroim2")
vol <- read_vol(file_name)Working with image volumes
Information about the geometry of the image volume is shown here:
print(vol)
#>
#> === NeuroVol Object ===
#>
#> * Basic Information
#> Type: DenseNeuroVol
#> Dimensions: 64 x 64 x 25 (406.6 Kb)
#> Total Voxels: 102,400
#>
#> * Data Properties
#> Value Range: [0.00, 1.00]
#>
#> * Spatial Properties
#> Spacing: 3.50 x 3.50 x 3.70 mm
#> Origin: 112.0, -108.5, -46.2 mm
#> Axes: Right-to-Left x Posterior-to-Anterior x Inferior-to-Superior
#>
#> ======================================
#>
#> Access Methods:
#> . Get Slice: slice(object, zlevel=10)
#> . Get Value: object[i, j, k]
#> . Plot: plot(object) # shows multiple slicesread_vol returns an object of class
NeuroVol object which extends an R array and
has 3 dimensions (x,y,z).
class(vol)
#> [1] "DenseNeuroVol"
#> attr(,"package")
#> [1] "neuroim2"
is.array(vol)
#> [1] TRUE
dim(vol)
#> [1] 64 64 25
vol[1,1,1]
#> [1] 0
vol[64,64,24]
#> [1] 0Arithmetic can be performed on images as if they were ordinary
arrays:
Inspecting geometry and spatial metadata
Each NeuroVol has an associated NeuroSpace
describing its geometry (dimensions, spacing, origin,
axes/orientation).
sp <- space(vol)
sp # human-readable summary
#>
#> NeuroSpace Object
#>
#> >> Dimensions
#> Grid Size: 64 x 64 x 25
#> Memory: 5.9 KB
#>
#> >> Spatial Properties
#> Spacing: 3.50 x 3.50 x 3.70 mm
#> Origin: 112.00 x -108.50 x -46.25 mm
#>
#> >> Anatomical Orientation
#> X: Right-to-Left | Y: Posterior-to-Anterior | Z: Inferior-to-Superior
#>
#> >> World Transformation
#> Forward (Voxel to World):
#> -3.500 -0.000 -0.000 112.000
#> 0.000 3.500 -0.000 -108.500
#> -0.000 0.000 3.700 -46.250
#> 0.000 0.000 0.000 1.000
#> Inverse (World to Voxel):
#> -0.286 -0.000 -0.000 32.000
#> 0.000 0.286 0.000 31.000
#> 0.000 0.000 0.270 12.500
#> 0.000 0.000 0.000 1.000
#>
#> >> Bounding Box
#> Min Corner: -108.5, -108.5, -46.2 mm
#> Max Corner: 112.0, 112.0, 42.6 mm
#>
#> ==================================================
dim(vol) # spatial dimensions (x, y, z)
#> [1] 64 64 25
spacing(vol) # voxel size in mm
#> [1] 3.5 3.5 3.7
origin(vol) # image origin
#> [1] 112.00 -108.50 -46.25You can convert between indices, voxel grid coordinates, and real-world coordinates:
idx <- 1:5
g <- index_to_grid(vol, idx) # 1D index -> (i,j,k)
w <- index_to_coord(vol, idx) # 1D index -> world coords
idx2 <- coord_to_index(vol, w) # back to indices
all.equal(idx, idx2)
#> [1] TRUEA numeric image volume can be converted to a binary image as follows:
vol2 <- as.logical(vol)
class(vol2)
#> [1] "LogicalNeuroVol"
#> attr(,"package")
#> [1] "neuroim2"
print(vol2[1,1,1])
#> [1] FALSEMasks and LogicalNeuroVol
Create a mask from a threshold or an explicit set of indices. Masks
are LogicalNeuroVol and align with the 3D space.
# Threshold-based mask
mask1 <- as.mask(vol > 0.5)
mask1
#>
#> === NeuroVol Object ===
#>
#> * Basic Information
#> Type: LogicalNeuroVol
#> Dimensions: 64 x 64 x 25 (406.6 Kb)
#> Total Voxels: 102,400
#>
#> * Data Properties
#> Value Range: [0.00, 1.00]
#>
#> * Spatial Properties
#> Spacing: 3.50 x 3.50 x 3.70 mm
#> Origin: 112.0, -108.5, -46.2 mm
#> Axes: Right-to-Left x Posterior-to-Anterior x Inferior-to-Superior
#>
#> ======================================
#>
#> Access Methods:
#> . Get Slice: slice(object, zlevel=10)
#> . Get Value: object[i, j, k]
#> . Plot: plot(object) # shows multiple slices
# Index-based mask
idx_hi <- which(vol > 0.8)
mask2 <- as.mask(vol, idx_hi)
sum(mask2) == length(idx_hi)
#> [1] TRUE
# Use a mask to compute a summary
mean_in_mask <- mean(vol[mask1@.Data])
mean_in_mask
#> [1] 1We can also create a NeuroVol instance from an
array or numeric vector. First we consruct a
standard R array:
Now we reate a NeuroSpace instance that describes the
geometry of the image including, at minimum, its dimensions and voxel
spacing.
bspace <- NeuroSpace(dim=c(64,64,64), spacing=c(1,1,1))
vol <- NeuroVol(x, bspace)
vol
#>
#> === NeuroVol Object ===
#>
#> * Basic Information
#> Type: DenseNeuroVol
#> Dimensions: 64 x 64 x 64 (2 Mb)
#> Total Voxels: 262,144
#>
#> * Data Properties
#> Value Range: [0.00, 0.00]
#>
#> * Spatial Properties
#> Spacing: 1.00 x 1.00 x 1.00 mm
#> Origin: 0.0, 0.0, 0.0 mm
#> Axes: Left-to-Right x Posterior-to-Anterior x Inferior-to-Superior
#>
#> ======================================
#>
#> Access Methods:
#> . Get Slice: slice(object, zlevel=10)
#> . Get Value: object[i, j, k]
#> . Plot: plot(object) # shows multiple slicesWe do not usually have to create NeuroSpace objects,
because geometric information about an image is automatically determined
from information stored in the image file header. Thus,
NeuroSpace objects are usually copied from existing images
using the space extractor function when needed:
vol2 <- NeuroVol((vol+1)*25, space(vol))
max(vol2)
#> [1] 25
space(vol2)
#>
#> NeuroSpace Object
#>
#> >> Dimensions
#> Grid Size: 64 x 64 x 64
#> Memory: 5.9 KB
#>
#> >> Spatial Properties
#> Spacing: 1.00 x 1.00 x 1.00 mm
#> Origin: 0.00 x 0.00 x 0.00 mm
#>
#> >> Anatomical Orientation
#> X: Left-to-Right | Y: Posterior-to-Anterior | Z: Inferior-to-Superior
#>
#> >> World Transformation
#> Forward (Voxel to World):
#> 1.000 0.000 0.000 0.000
#> 0.000 1.000 0.000 0.000
#> 0.000 0.000 1.000 0.000
#> 0.000 0.000 0.000 1.000
#> Inverse (World to Voxel):
#> 1.000 0.000 0.000 0.000
#> 0.000 1.000 0.000 0.000
#> 0.000 0.000 1.000 0.000
#> 0.000 0.000 0.000 1.000
#>
#> >> Bounding Box
#> Min Corner: 0.0, 0.0, 0.0 mm
#> Max Corner: 63.0, 63.0, 63.0 mm
#>
#> ==================================================Slicing and quick visualization
Extract a single 2D slice for display using standard array indexing:
Mid-slice of example volume
Reorienting and resampling
You can change an image’s orientation and voxel spacing. Use
reorient() to remap axes (e.g., to RAS) and
resample_to() to match a target space.
# Reorient the space (LPI -> RAS) and compare coordinate mappings
sp_lpi <- space(vol)
sp_ras <- reorient(sp_lpi, c("R","A","S"))
g <- t(matrix(c(10, 10, 10)))
world_lpi <- grid_to_coord(sp_lpi, g)
world_ras <- grid_to_coord(sp_ras, g)
# world_lpi and world_ras differ due to axis remappingResample to a new spacing or match a target
NeuroSpace:
# Create a target space with 2x finer resolution
sp <- space(vol)
sp2 <- NeuroSpace(sp@dim * c(2,2,2), sp@spacing/2, origin=sp@origin, trans=trans(vol))
# Resample (trilinear)
vol_resamp <- resample_to(vol, sp2, method = "linear")
dim(vol_resamp)Downsampling
Reduce spatial resolution to speed up downstream operations.
# Downsample by target spacing
vol_ds1 <- downsample(vol, spacing = spacing(vol)[1:3] * 2)
dim(vol_ds1)
#> [1] 32 32 32
# Or by factor
vol_ds2 <- downsample(vol, factor = 0.5)
dim(vol_ds2)
#> [1] 32 32 32Writing a NIFTI formatted image volume
When we’re ready to write an image volume to disk, we use
write_vol
write_vol(vol2, "output.nii")
## adding a '.gz' extension results ina gzipped file.
write_vol(vol2, "output.nii.gz")You can also write to a temporary file during workflows:
tmp <- tempfile(fileext = ".nii.gz")
write_vol(vol2, tmp)
file.exists(tmp)
#> [1] TRUE
unlink(tmp)