Accessing Information from a nifti

To access additional information, called a slot, you can use the @ operator. We do not recommend this, as there should be a function implemented to “access” this slot. These are hence called accessor functions (they access things!). For example, if you want to get the cal_max slot of a nifti object, you should use the cal_max function. If an accessor function is not implemented, you should still use the slot(object, name) syntax over @.

Here’s an example where we make an array of random normal data, and put that array into a nifti object with the nifti function:

nim@cal_max

## [1] 2.706524

cal_max(nim)

## [1] 2.706524

slot(nim, "cal_max")

## [1] 2.706524

Accessing the “data”

If you want to access the “data” of the image, you can access that using:

data = slot(nim, ".Data")
class(data)

## [1] "array"

With newer versions of oro.nifti (especially that on GitHub and in Neuroconductor), there is a img_data function to access the data:

data = oro.nifti::img_data(nim)
class(data)

## [1] "array"

dim(data)

## [1] 10 10 10

This array is 3-dimensional and can be subset using normal square-bracket notations ([row, column, slice]). Thus, if we want the 3rd “slice” of the image, we can use:

slice = data[,,3]
class(slice)

## [1] "matrix" "array"

Thus we see we get a matrix of values from the 3rd “slice”. We should note that we generally reference an image by x, y, and z planes (in that order). Most of the time, the x direction refers to going left/right on an image, y refers to front/back (or anterior/posterior), and the z direction refers to up/down (superior/inferior). The actual direction depends on the header information of the NIfTI image.

slice = data[,,3, drop = FALSE]
class(slice)

## [1] "array"

Show all slots

You can see which slots exist for a nifti object by using slotNames

slotNames(nim)

##  [1] ".Data"          "sizeof_hdr"     "data_type"      "db_name"       
##  [5] "extents"        "session_error"  "regular"        "dim_info"      
##  [9] "dim_"           "intent_p1"      "intent_p2"      "intent_p3"     
## [13] "intent_code"    "datatype"       "bitpix"         "slice_start"   
## [17] "pixdim"         "vox_offset"     "scl_slope"      "scl_inter"     
## [21] "slice_end"      "slice_code"     "xyzt_units"     "cal_max"       
## [25] "cal_min"        "slice_duration" "toffset"        "glmax"         
## [29] "glmin"          "descrip"        "aux_file"       "qform_code"    
## [33] "sform_code"     "quatern_b"      "quatern_c"      "quatern_d"     
## [37] "qoffset_x"      "qoffset_y"      "qoffset_z"      "srow_x"        
## [41] "srow_y"         "srow_z"         "intent_name"    "magic"         
## [45] "extender"       "reoriented"

If you would like to see information about each one of these slots, please see this blog post about the NIfTI header.

Differences between orthographic and ortho2

The above code does not fully illustrate the differences between orthographic and ortho2. One marked difference is when you would like to “overlay” an image on top of another in an orthographic view. Here we will highlight voxels greater than the 90th quantile of the image:

orthographic(eve, y = eve > quantile(eve, 0.9))

We see that the white matter is represented here, along with the burned in orientations, but we would like to see areas of the brain that are not over this quantile to be shown as the image. Let us contrast this with:

ortho2(eve, y = eve > quantile(eve, 0.9))

We see the image where the mask is 0 shows the original image. This is due to the NA.y argument in ortho2. The ortho2 (and orthograhic) function is based on the graphics::image function in R, as well as many other functions we will discuss below. When graphics::image sees an NA value, it does not plot anything there. The NA.y argument in ortho2 makes it so any values in the y argument (in this case the mask) that are equal to zero are turned to NA.

Bright values

If you have artifacts or simply large values of an image, it can “dampen” the viewing of an image. Let’s make one value of the eve template very large. We will set the voxel with the largest value to be that value times 5 :

eve2 = eve
eve2[which.max(eve)] = eve2[which.max(eve)] * 5

Let’s plot this image again:

ortho2(eve2)

We see a faint outline of the image, but this single large value affects how we view the image. The function robust_window calculates quantiles of an image, by default the 0 (min) and 99.9th quantile, and sets values outside of this range to that quantile. If you are familiar with the process of Winsorizing, this is the exact same procedure. Many times we use this function to plotting, but could be thought of an outlier dampening procedure. Let’s plot this windowed image:

ortho2(robust_window(eve2))

Changing the probs argument in robust_window, which is passed to quantile, can also be used to limit artifacts with remarkably low values. The zlim option can also denote which range of intensities that can be plotted:

ortho2(eve2, zlim = quantile(eve2, probs = c(0, 0.999)))

This is a bit more like trimming, however.

Double orthographic view

Sometimes you would like to represent 2 images side by side, of the same dimensions and orientation of course. The double_ortho function allows you to do this. Let’s read in the flipped image:

flipped_fname = system.file(
  file.path("nifti", "mniLR.nii.gz"),
  package = "oro.nifti")
flipped_eve = readnii(flipped_fname)

We can view the flipped T1 alongside the original image:

double_ortho(flipped_eve, eve)

Different Planes

We can specify z the same way but change the plane to be different to get a different slice of the brain:

image(eve, z = 52, plot.type = "single", plane = "sagittal")

We can similarly do the same for “coronal” slices.

Dropping with multiple images

You can pass in other images in the other.imgs function to applying this dropping procedure. For example, let’s say you have 3 images, a T1, T2, and FLAIR image that are all registered to the T1 image and have a brain mask from the T1 image. You can pass in the mask image and pass in the other images into other.imgs so they all drop the same slices and are the same dimensions after the procedure (as they were the same prior), whereas if you performed the operation you may not be ensured to drop exactly the same slices due to some modalities allowing values of zero (albeit highly unlikely).

Adding back the dropped dimensions

To reverse this procedure, the replace_dropped_dimensions function will add back dimensions to the image correctly using the indices from drop_empty_dim. Here is an example:

dd = dropEmptyImageDimensions(eve, keep_ind = TRUE)
reduced = dd$outimg
reversed = replace_dropped_dimensions(img = reduced, 
                                      inds = dd$inds,
                                      orig.dim = dd$orig.dim)
all(reversed == eve)

## [1] TRUE

Copying header information

If you have an array, you can simply write:

dims = dim(eve)
arr = array(rnorm(prod(dims)), dim = dims)
nim = nifti(arr)
nim

## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 2 (UINT8)
##   Bits per Pixel  : 8
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 0 (Unknown)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 91 x 109 x 91
##   Pixel Dimension : 1 x 1 x 1
##   Voxel Units     : Unknown
##   Time Units      : Unknown

but the header information of the nifti output nim does not match that of eve. The copyNIfTIHeader function allows you to…copy the NIfTI header:

nim = copyNIfTIHeader(img = eve, arr = nim)
nim

## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 4 (INT16)
##   Bits per Pixel  : 16
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 0 (Unknown)
##   Sform Code      : 4 (MNI_152)
##   Dimension       : 91 x 109 x 91
##   Pixel Dimension : 2 x 2 x 2
##   Voxel Units     : mm
##   Time Units      : sec

Creating a “shell”

The niftiarr function does much of the same functionality of copyNIfTIHeader, but copyNIfTIHeader you have an array or nifti in the arr argument. If you wanted a nifti object with the same header as eve, but all zeroes, you can use niftiarr:

zeroes = niftiarr(eve, 0)

The main difference is the line in niftiarr:

if (!is(arr, "array")){
  arr = array(arr, dim=dim(img))
}

which implies that if you pass in a vector instead of an array, it will create an array on the fly. Sometimes you want an operation to error if arr is not an array (as in copyNIfTIHeader) or be able to pass in a vector and get the correct output niftiarr.

Technical note: This code is legacy and somewhat old and probably can (and may be) replaced by making copyNIfTIHeader a generic and having different versions for when arr is an array or vector.

Subsetting via masks

In a previous example, we calculated the mean of eve over the entire image. Many times we want to calculate values over a mask. For example, let’s get the mean of all voxels in the mask, where the mask is any value of eve greater than zero. We can do this by subsetting the values in the mask and then calculating the mean:

vals = eve[ eve > 0 ]
mean(vals)

## [1] 82.90267

mean(eve)

## [1] 82.89716

We see that the mean of the voxels in the mask versus all voxels is very different, because the voxels not in the mask are simply adding zeroes to the calculation. We could simply do the same by making zero voxels NA and adding the na.rm argument to TRUE for the mean.

mask = eve > 0
mask[ mask == 0 ] = NA
eve_masked = mask_img(eve, mask)
mean(eve_masked)

## [1] NA

mean(eve_masked, na.rm = TRUE)

## [1] 82.90267

Again, as nifti objects inherits properties from array objects, we can subset using logical indices, as above in mask == 0 and could use indices such as which(mask == 0).

We can do a marginal density of the values only in the mask:

vals = eve[ eve > 0 ]

plot(density(vals))

For most plotting functions, we can also add in a mask argument to plot values only for that masked data.

plot(density(eve, mask = eve > 0))