EBImage

We primarily use images via the EBImage package.

# make random EBImage image
f <- system.file("images", "sample.png", package = "EBImage")

# read image with EBImage
img <- readImage(f)

Similar to its counterparts (see writeHDF5Array in HDF5Array), we use the writeImageArray function to write pyramid images to preferred locations. Here, n.levels control the amount of levels that the pyramid image should have.

output_h5 <- tempfile(fileext = ".h5")
imgarray <- writeImageArray(img, output = output_h5, n.levels = 2)
imgarray

## ImageArray Object (x,y) 
## Scales (2): (768,512) (384,256)

Each level of a pyramid can be rasterized at any time, and thus plotted.

bfa.raster <- as.raster(imgarray, level = 2)
plot(bfa.raster)

However, the true functionality of ImageArray is to pick pyramid levels in a memory-efficient manner. Packages, applications and workflows that use ImageArray can rasterize images by defining thresholds for pixel dimensions across all levels.

Here, by using the max.pixel.size, we request ImageArray to return a pyramid level whose both width (X) and height (Y) are lower than, for example, 400. This approach is particularly useful when visualizing:

• 1. the entire image
• 2. a subset of the image with higher resolution

with minimal memory footprint.

# visualize
bfa.raster <- as.raster(imgarray, max.pixel.size = 400)
dim(bfa.raster)

## [1] 256 384

Operations such as rotate, flip, flop or negate will be conducted on all levels (without loading the image in the memory) which then can be rasterized and plotted.

imgarray_rotated <- rotate(imgarray, angle = 90)
imgarray_rotated

## ImageArray Object (x,y) 
## Scales (2): (512,768) (256,384)

# visualize
bfa.raster <- as.raster(imgarray_rotated, max.pixel.size = 400)
plot(bfa.raster)

You can even crop images, either by using crop or [ methods.

imgarray_cropped <- imgarray[300:500, 200:300]
imgarray_cropped

## ImageArray Object (x,y) 
## Scales (2): (201,101) (101,51)

# visualize
bfa.raster <- as.raster(imgarray_cropped, max.pixel.size = 120)
plot(bfa.raster)

# cropping with indices
imgarray_cropped <- crop(imgarray, index = list(300:500, 200:300))
imgarray_cropped

## ImageArray Object (x,y) 
## Scales (2): (201,101) (101,51)

Cropping can also be performed with named indices, thus allowing users to define indices for only the target dimensions. For now, ImageArray only includes x and y dimensions and ‘c’ for image channels.

# cropping with named indices
imgarray_cropped <- crop(imgarray, index = list(x=300:500, y=200:300))
imgarray_cropped

## ImageArray Object (x,y) 
## Scales (2): (201,101) (101,51)

The crop function serves as a more flexible alternative to the traditional subsetting method [. The named list of index can be of any order while [ requires the order matching with the order of axis (as revealed in the show method, i.e. (x,y)).

# index in yx order, results in same slice as xy order
imgarray_cropped <- crop(imgarray, index = list(y=200:300, x=300:500))
imgarray_cropped

## ImageArray Object (x,y) 
## Scales (2): (201,101) (101,51)

# subsetting with `[` provides a different slice
imgarray_cropped <- imgarray[200:300, 300:500]
imgarray_cropped

## ImageArray Object (x,y) 
## Scales (2): (101,201) (51,101)

magick images

ImageArray is also compatible with the magick package where we retain the dimensionality of the image similar to EBImage.

# make random EBImage image
f <- system.file("images", "sample.png", package = "EBImage")

# read image with magick
img <- image_read(f)

# create image array
output_h5 <- tempfile(fileext = ".h5")
imgarray <- writeImageArray(img, output = output_h5)
imgarray

## ImageArray Object (c,x,y) 
## Scales (2): (1,768,512) (1,384,256)

# plot
bfa.raster <- as.raster(imgarray, level = 2)
plot(bfa.raster)

OME.TIFF (Bio-Formats) images

You can create ImageArray objects from OME-TIFF files (or any file compatible with Bio-formats, see RBioFormats) with already defined layers.

ome.tiff.file <- system.file("extdata", "xy_12bit__plant.ome.tiff",
                             package = "ImageArray")
read.metadata(ome.tiff.file)

## ImageMetadata list of length 2 
## 
##  series res sizeX sizeY sizeC sizeZ sizeT total
##  1      1   512   512   1     1     1     1    
##  1      2   256   256   1     1     1     1    
## 
## globalMetadata:List of 1058
##  $ xy_12bit__plant.oir inTrigger triggerNo #5                                : chr "0"
##  $ xy_12bit__plant.oir inTrigger triggerNo #4                                : chr "1"
##  $ xy_12bit__plant.oir laser name #7                                         : chr "LD640"
##  $ xy_12bit__plant.oir inTrigger triggerNo #3                                : chr "0"
##  $ xy_12bit__plant.oir laser name #6                                         : chr "LD594"
##   [list output truncated]

# define ImageArray object
imgarray <- ImageArray(ome.tiff.file, series = 1, resolution = 1:2)
imgarray

## ImageArray Object (x,y) 
## Scales (2): (512,512) (256,256)

You can again use either max.pixel.size or min.pixel.size to control the size of the image that loaded into the memory.

bfa.raster <- as.raster(imgarray, max.pixel.size = 300)
plot(bfa.raster)

dim(bfa.raster)

## [1] 256 256

bfa.raster <- as.raster(imgarray, min.pixel.size = 300)
plot(bfa.raster)

dim(bfa.raster)

## [1] 512 512

bfa.raster <- as.raster(imgarray, level = 1)
plot(bfa.raster)
dim(bfa.raster)

## [1] 512 512

ImageArray can also read images with all five dimensions (XYZCT) available in the OME (Open Microscopy Environment) specifications. See OME specifications webpage for more information here.

ome.tiff.file <- system.file("extdata", "4D-series.ome.tiff",
                             package = "ImageArray")

# read metadata
read.metadata(ome.tiff.file, proprietary.metadata = TRUE, 
              filter.metadata = TRUE)

## ImageMetadata
##  $ coreMetadata:List of 18
##   ..$ sizeX    : int 439
##   ..$ sizeY    : int 167
##   ..$ sizeZ    : int 5
##   ..$ sizeC    : int 1
##   ..$ sizeT    : int 7
##   ..$ pixelType: chr "int8"
##   .. [list output truncated]

# define ImageArray object
imgarray <- ImageArray(ome.tiff.file, series = 1, resolution = 1)
imgarray

## ImageArray Object (x,y,z,t) 
## Scales (1): (439,167,5,7)

ImageArray is also capable of reading other BioFormats compatible image formats (e.g. qptiff, svs). For more information, visit the Supported Formats webpage of BioFormats.

translation

Translation transformation will not alter the pixel dimensionality of the image, but will change the physical space that the image occupies. You can view the physical space that the image spans using the extent function.

imgarray_trans <- translation(imgarray,
                              shift = c(100,20))
imgarray_trans

## ImageArray Object (x,y) 
## Scales (2): (512,512) (256,256)

extent(imgarray_trans)

## $x
## [1] 100 612
## 
## $y
## [1]  20 532

scale

Scale transformations are performed providing the dimensionality of the new scaled ImageArray object, i.e. output.dim for x and y dimensions. Here, we use the filter argument to determine the method of interpolation of pixel intensity values after transformations. Regular images typically use bilinear whereas none invokes a nearest neighbor interpolation approach more suitable for image masks.

imgarray_scale <- scale(imgarray, 
                        output.dim = c(200,300),
                        filter = "bilinear")
imgarray_scale

## ImageArray Object (x,y) 
## Scales (2): (200,300) (100,150)

rotate

Rotation requires specifying the angle that rotates the image in clock-wise fashion and filter for interpolation method.

imgarray_rotate <- rotate(imgarray, 
                          angle = 45, 
                          filter = "bilinear")
imgarray_rotate

## ImageArray Object (x,y) 
## Scales (2): (724,724) (362,362)

plot(as.raster(imgarray_rotate))

affine

Affine transformations are also available where users have to provide a 2x3 matrix that applies an affine transformation on x and y axes.

m <- matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
imgarray_affine <- affine(imgarray, 
                          m = m,
                          filter = "bilinear")
imgarray_affine

## ImageArray Object (x,y) 
## Scales (2): (768,512) (384,256)

plot(as.raster(imgarray_affine))

sequence

Because the transformations are delayed, they can also be performed in sequence. This is useful when an image must first be then rotated, scaled and finally shifted onto a new coordinate system.

m <- matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
imgarray <- affine(imgarray, 
                   m = m, 
                   filter = "bilinear")
imgarray <- scale(imgarray,
                  output.dim = c(200,300),
                  filter = "bilinear")
imgarray <- translation(imgarray,
                        shift = c(100,20))
plot(as.raster(imgarray))

Visualizing H&E images

For this example, we will use a large H&E image used by 10x Genomics, generated after a Xenium in Situ platform run (i.e. postXenium H&E).

image_file <- paste(
  "https://cf.10xgenomics.com/samples/xenium/1.0.1",
  "Xenium_FFPE_Human_Breast_Cancer_Rep1",
  "Xenium_FFPE_Human_Breast_Cancer_Rep1_he_image.ome.tif",
  sep = "/")
library(BiocFileCache)
bfc <- BiocFileCache()
image_file <- bfcrpath(bfc, image_file)

## adding rname 'https://cf.10xgenomics.com/samples/xenium/1.0.1/Xenium_FFPE_Human_Breast_Cancer_Rep1/Xenium_FFPE_Human_Breast_Cancer_Rep1_he_image.ome.tif'

Let us create an ImageArray object first with the image above. The OME-TIFF file includes multiple resolutions that we can ask ImageArray object to include.

imgarray <- ImageArray(image_file,
                       series = 1,
                       resolution = 1:6)
                             
imgarray

## ImageArray Object (x,y,c) 
## Scales (6): (30786,24241,3) (15393,12120,3) ... (1924,1515,3)
##   (962,757,3)

You can visualize subsets of large images really quickly using crop and as.raster. Again, here we use max.pixel.size to optimize the level parsed from the pyramid; that is, as defined below, only resolutions whose width and height are both under size 800 would be returned, whether it is the entire image, or a subset.

# crop image
imgarray_sub <- crop(imgarray, index = list(16000:19000, 7000:10000, NULL))

# convert to raster
img_raster <- as.raster(imgarray_sub, max.pixel.size = 800)
dim(img_raster)

## [1] 752 752

# plot with ggplot
ggplot(data.frame(x = 0, y = 0), aes(x, y)) +
  coord_fixed(expand = FALSE,
              xlim = c(0, dim(img_raster)[2]),
              ylim = c(0, dim(img_raster)[1])) +
  annotation_raster(img_raster, 
                    0, dim(img_raster)[2],
                    dim(img_raster)[1], 0, interpolate = FALSE)

Interactive Shiny Applications

Now let us create a shiny application where we can interactively subset and visualize the OME-TIFF image quickly.

The Shiny application lets the user to define image subsets (or slices) before visualizing without loading large images in memory.

This is because as.raster determines the layer with the pixel dimensions of he specified subset not exceeding the defined threshold (max.pixel.size).

if (interactive()) {

  # variables
  dimimg <- dim(imgarray)
  max.pixel.size <- 800

  # Define UI
  ui <- fluidPage(
    sliderInput("x_slider", label = "X", min = 1,
                max = dimimg[1], value = c(16000, 19000)),
    sliderInput("y_slider", label = "Y", min = 1,
                max = dimimg[2], value = c(7000, 10000)),
    mainPanel(
      plotOutput(outputId = "scatterPlot")
    )
  )

  # Define server logic
  server <- function(input, output) {

    output$scatterPlot <- renderPlot({

      # crop image
      indx <- seq(input$x_slider[1], input$x_slider[2])
      indy <- seq(input$y_slider[1], input$y_slider[2])
      imgarray_sub <- crop(imgarray, index = list(indx, indy, NULL))

      # convert to raster
      img_raster <- as.raster(imgarray_sub, max.pixel.size = max.pixel.size)

      # plot with ggplot
      imgggplot <- ggplot(data.frame(x = 0, y = 0), aes(x, y)) +
        coord_fixed(expand = FALSE,
                    xlim = c(0, dim(img_raster)[2]),
                    ylim = c(0, dim(img_raster)[1])) +
        annotation_raster(img_raster, 
                          0, dim(img_raster)[2],
                          dim(img_raster)[1], 0, interpolate = FALSE)
      imgggplot
    })
  }

  # Run the app
  shinyApp(ui = ui, server = server)
}