Resampling an Image

All code for this document is located at here.

Goal

In this short tutorial, we will show various ways of resampling an image to different voxel sizes and image dimensions. This procedure may be necessary in some analyses and can also be used to reduce the size of the image for demonstration purposes.

Reading in one file

We will read in one T1-weighted NIfTI image to see the dimensions of the image. We can read in the data a number of ways, including the neurobase::readnii function, which relies on the oro.nifti::readNIfTI function, the RNifti::readNifti, and ANTsRCore::antsImageRead functions. Note, all of these functions give back different types of objects.

fname = kirby21.t1::get_t1_filenames()[1]
if (file.exists(fname)) {
  oro_img = readnii(fname)
  print(oro_img)
  rnif_img = readNifti(fname)
  print(rnif_img)
  ants_img = antsImageRead(fname)
  print(ants_img)
}
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 170 x 256 x 256
  Pixel Dimension : 1.2 x 1 x 1
  Voxel Units     : mm
  Time Units      : sec
Image array of mode "double" (85 Mb)
- 170 x 256 x 256 voxels
- 1.2 x 1 x 1 mm per voxel
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 170x256x256 
  Voxel Spacing       : 1.20000004768372x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 
  Filename           : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/kirby21.t1/visit_1/113/113-01-T1.nii.gz 

Conversion between objects

Here we will show some ways to convert from one format to another. Note, you should always check cases of conversion, such as the data type, especsially when there are missing values or NaN values.

oro.nifti nifti RNifti niftiImage format

The RNifti function asNifti allows for a direct conversion of formats.

asNifti(oro_img)
Image array of mode "double" (85 Mb)
- 170 x 256 x 256 voxels
- 1.2 x 1 x 1 mm per voxel
all(asNifti(oro_img) == rnif_img)
[1] TRUE

RNifti niftiImage to oro.nifti nifti format

The oro.nifti function nii2oro will convert from niftiImage object to nifti objects.

nii2oro(rnif_img)
NIfTI-1 format
  Type            : nifti
  Data Type       : 64 (FLOAT64)
  Bits per Pixel  : 64
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 170 x 256 x 256
  Pixel Dimension : 1.2 x 1 x 1
  Voxel Units     : mm
  Time Units      : sec
all(oro_img == nii2oro(rnif_img))
[1] TRUE

Note, the slots for the header are directly copied from the output of niftiHeader, even if they are inconsistent in what readnii would give. For example, the scl_slope argument may not be the same:

scl_slope(oro_img)
[1] 1
scl_slope(nii2oro(rnif_img))
[1] 0

ANTsRCore antsImage to oro.nifti nifti format

The extrantsr functions ants2oro will convert from ANTsRCore antsImage to oro.nifti nifti objects:

ants2oro(ants_img)
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 170 x 256 x 256
  Pixel Dimension : 1.2 x 1 x 1
  Voxel Units     : mm
  Time Units      : Unknown

By default, this will write the image to disk and then read it in using readnii, but you can pass in a reference image that header information will be copied from (but not checked for consistency):

ants2oro(ants_img, reference = oro_img)
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 170 x 256 x 256
  Pixel Dimension : 1.2 x 1 x 1
  Voxel Units     : mm
  Time Units      : sec

which we see is much faster:

system.time(ants2oro(ants_img))
   user  system elapsed 
  3.336   0.101   3.441 
system.time(ants2oro(ants_img, reference = oro_img))
   user  system elapsed 
  0.578   0.027   0.606 

oro.nifti nifti to ANTsRCore antsImage format

The extrantsr functions oro2ants will convert from oro.nifti nifti to ANTsRCore antsImage objects. Note, this essentially writes the output to disk, then reads it back in using antsImageRead, but is a convenience image for this:

oro2ants(oro_img)
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 170x256x256 
  Voxel Spacing       : 1.20000004768372x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 
  Filename           : /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmp253ES1/file367e658de1c6.nii.gz 

The as.antsImage function should work as well (and is much faster), but doesn’t carry the appropriate header information:

as.antsImage(oro_img)
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 170x256x256 
  Voxel Spacing       : 1x1x1 
  Origin              : 0 0 0 
  Direction           : 1 0 0 0 1 0 0 0 1 

And can error as the header information is different:

as.antsImage(oro_img) == ants_img
Error in as.antsImage(oro_img) == ants_img: Images do not occupy the same physical space

You can specify this with a reference, or the arguments of spacing, origin, and direction (which may be difficult), however:

all(as.array(as.antsImage(oro_img, reference = ants_img) == ants_img) == 1)
[1] TRUE

Also, niftiImage objects work as well using oro2ants:

oro2ants(oro_img)
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 170x256x256 
  Voxel Spacing       : 1.20000004768372x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 
  Filename           : /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmp253ES1/file367e3720266.nii.gz 

but does not work directly with as.antsImage

as.antsImage(rnif_img)
Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'as.antsImage' for signature '"niftiImage"'

Resampling Images

Using FSL and fslr

The fslr package calls FSL to perform operations on images that are on disk. For functions that pass in nifti objects in the fslr package, these images are written to disk using in a temporary location and then the FSL function is applied. If a resulting image is returned, then it can be read back in as a nifti object.

Let’s just display the pixel dimensions of the original image:

oro_img@pixdim[2:4]
[1] 1.2 1.0 1.0

The fsl_resample function allows you to resample an image on disk. Thus, if you have the character path of the image or nifti object, these images can be resampled using fsl_resample. This code calls flirt from FSL and uses the applyisoxfm to perform isotoropic sampling (all the same size). Thus, fsl_resample can only resample images to cubes with all the same side length. Some of the examples below have more flexibility.

res = fsl_resample(file = fname, voxel_size = 1)
flirt -in "/Library/Frameworks/R.framework/Versions/4.0/Resources/library/kirby21.t1/visit_1/113/113-01-T1.nii.gz"  -ref /Library/Frameworks/R.framework/Versions/4.0/Resources/library/kirby21.t1/visit_1/113/113-01-T1.nii.gz -applyisoxfm 1 -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmp253ES1/file367e3883b55d";
res
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 1 (Scanner_Anat)
  Dimension       : 204 x 256 x 256
  Pixel Dimension : 1 x 1 x 1
  Voxel Units     : mm
  Time Units      : sec
oro_res = fsl_resample(oro_img, voxel_size = 1)
flirt -in "/private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmp253ES1/file367e4ef809a7.nii.gz"  -ref /private/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T/Rtmp253ES1/file367e4ef809a7.nii.gz -applyisoxfm 1 -out "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//Rtmp253ES1/file367e7e41db71";
all.equal(res, oro_res)
[1] TRUE
res@pixdim[2:4]
[1] 1 1 1

Using RNiftyReg rescale function

To resample a niftiImage object, you need to use the RNiftyReg package. In this case, the scales argument is the “Scale factors along each axis”, which is not the voxel sizes. The output dimensions will be the current dimensions multiplied by the scale factors. Alternatively, you can think of the output voxel dimensions as the current voxel dimensions divided by the scale factors. Thus, if we run the following code:

niftiHeader(rnif_img)$pixdim[2:4]
[1] 1.2 1.0 1.0
RNiftyReg::rescale(rnif_img, c(2,1.5,1))
Image array of mode "double" (255 Mb)
- 340 x 384 x 256 voxels
- 0.6 x 0.6667 x 1 mm per voxel

The image will be scaled by 2 in the x-direction, 1.5 in the y-direction and stay the same for the z-direction. Thus, if we want an image to be 1x1x1, then we need to scale the dimensions as follows:

RNiftyReg::rescale(rnif_img, c(1.2,1,1))
Image array of mode "double" (102 Mb)
- 204 x 256 x 256 voxels
- 1 x 1 x 1 mm per voxel

If we want to make the output image 2.5 mm\(^3\), then we can run:

pdim = niftiHeader(rnif_img)$pixdim[2:4]
RNiftyReg::rescale(rnif_img, pdim / c(2.5, 2.5, 2.5))
Image array of mode "double" (6.4 Mb)
- 81 x 102 x 102 voxels
- 2.5 x 2.5 x 2.5 mm per voxel

Or we can refactor the arguments such that if we want an image to be 2.5x2.3x3 we can run:

rescale2 = function(img, voxel_size = c(1,1,1), ...) {
  pdim = niftiHeader(rnif_img)$pixdim[2:4]
  scales = pdim / voxel_size
  RNiftyReg::rescale(img, scales = scales, ...)
}
rescale2(rnif_img, c(2.5, 2.3, 3))
Image array of mode "double" (5.8 Mb)
- 81 x 111 x 85 voxels
- 2.5 x 2.3 x 3 mm per voxel

NB: you can pass in nifti objects as well, and the output is converted to a niftiImage

rescale2(oro_img, c(2.5, 2.3, 3))
Image array of mode "double" (5.8 Mb)
- 81 x 111 x 85 voxels
- 2.5 x 2.3 x 3 mm per voxel

which you can wrap in nii2oro:

nii2oro(rescale2(oro_img, c(2.5, 2.3, 3)))
NIfTI-1 format
  Type            : nifti
  Data Type       : 64 (FLOAT64)
  Bits per Pixel  : 64
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 81 x 111 x 85
  Pixel Dimension : 2.5 x 2.3 x 3
  Voxel Units     : mm
  Time Units      : sec

Additional arguments - interpolation

The interpolation argument in rescale, which is passed to applyTransform in RNiftyReg can also specify the interpolation done, using 0 (nearest neighbour), 1 (trilinear) or 3 (cubic spline), where cubic spline is the default:

interp0 = nii2oro(rescale2(oro_img, c(2.5, 2.3, 3), interpolation = 0L))
interp3 = nii2oro(rescale2(oro_img, c(2.5, 2.3, 3), interpolation = 3L))
double_ortho(interp0, interp3)

Using ANTsRCore resampleImage and extrantsr resample_image

The resampleImage provides a flexible function to resample antsImage objects. The useVoxels argument determines what the resample parameters are specified in, either millimeters (default) or output voxel spacing

resampleImage(ants_img, resampleParams = c(1, 1, 1))
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 204x256x256 
  Voxel Spacing       : 1x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 
resampleImage(ants_img, resampleParams = c(1, 1, 1), useVoxels = TRUE)
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 1x1x1 
  Voxel Spacing       : 204.000008106232x256x256 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 
resampleImage(ants_img, resampleParams = c(150, 200, 200), useVoxels = TRUE)
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 150x200x200 
  Voxel Spacing       : 1.36000005404154x1.28x1.28 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 

The extrantsr::resample_image wrapper function refactors the arguments so that you can pass in the target voxel sizes in mm or dimensions.

resample_image(oro_img, parameters = c(1, 1, 1))
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 1 (Scanner_Anat)
  Sform Code      : 0 (Unknown)
  Dimension       : 204 x 256 x 256
  Pixel Dimension : 1 x 1 x 1
  Voxel Units     : mm
  Time Units      : Unknown
resample_image(ants_img, parameters = c(1, 1, 1))
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 204x256x256 
  Voxel Spacing       : 1x1x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 

where we see the output depends on the input. Also, we can perform operations where the output is not isotropic:

resample_image(ants_img, parameters = c(2.5, 2.3, 1))
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 82x111x256 
  Voxel Spacing       : 2.5x2.3x1 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 

We can also specify the resampling by image dimensions:

resample_image(ants_img, parameters = c(150, 200, 200), parameter_type = "voxels")
antsImage
  Pixel Type          : float 
  Components Per Pixel: 1 
  Dimensions          : 150x200x200 
  Voxel Spacing       : 1.36000005404154x1.28x1.28 
  Origin              : 202.8 0 0 
  Direction           : -1 0 0 0 1 0 0 0 1 

Additional arguments - interpolation

As with rescale, we can perform different interpolations:

formals(resample_image)$interpolator
c("nearestneighbor", "linear", "gaussian", "windowedsinc", "bspline")
formals(resampleImage)$interpType
[1] "nearestneighbor"

where we see the nearest neighbor interpolation is the default:

interp0 = resample_image(ants_img, parameters = c(2.5, 2.3, 3), interpolator = "nearestneighbor")
interp3 = resample_image(ants_img, parameters = c(2.5, 2.3, 3), interpolator = "bspline")
double_ortho(interp0, interp3)

ortho2(interp0, interp3 - interp0, col.y = scales::alpha(hotmetal(), 0.5))