Working with 3D Image Volumes

2026-01-14

Reading a NIFTI formatted image volume

The way to read an volumetric image file is to use read_vol:

    library(neuroim2)
    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 (806.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

read_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] 0

Arithmetic can be performed on images as if they were ordinary arrays:

    
    vol2 <- vol + vol
    sum(vol2) == 2 * sum(vol)
#> [1] TRUE
    
    vol3 <- vol2 - 2*vol
    all(vol3 == 0)
#> [1] TRUE

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

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

A 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] FALSE

Masks 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] 1

We can also create a NeuroVol instance from an array or numeric vector. First we consruct a standard R array:

    x <- array(0, c(64,64,64))

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 slices

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

    z <- ceiling(dim(vol)[3] / 2)
    image(vol[,,z], main = paste("Slice z=", z))
Mid-slice of example volume (grayscale image).

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 remapping

Resample 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 32

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