---
title: "`ImageArray`"
date: "`r format(Sys.Date(), '%B %d, %Y')`"
package: "`r BiocStyle::pkg_ver('ImageArray')`"
author: 
  - name: Artur Manukyan
output:
  BiocStyle::html_document
vignette: |
  %\VignetteIndexEntry{ImageArray}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

The `r Biocpkg("ImageArray")` package provides a **unified, memory-efficient** 
framework for working with both **pyramidal** and **non-pyramidal** images in R 
by leveraging the `r Biocpkg("DelayedArray")` package. With 
`r Biocpkg("ImageArray")` you can store large images on disk (as HDF5 or Zarr), 
treat them as array-like objects, perform delayed/lazy operations (e.g., 
rotate, flip, crop) across all pyramid levels, and seamlessly integrate with 
other R package and workflows that use large bioimages. 

# Installation

You can install `r Biocpkg("ImageArray")` from Bioconductor with:

```{r install-BiocManager, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("ImageArray")
```

# Why pyramidal images?

Image pyramids are multi-resolution representations that start with a full 
resolution image followed by a series of down-sampled levels (e.g., half 
resolutions). These are common in microscopy, digital pathology and large-scale 
imaging because:

- They allow fast zooming and visualization (you don’t always need full 
resolution to inspect large images).  
- They support scale-aware analysis (you may operate on coarse levels when 
doing overview tasks, fine levels for detail).  
- Storing data in a pyramidal stack often reduces memory footprint and 
enables on-disk processing.

In our context, an `r Biocpkg("ImageArray")` object typically contain 
**multiple levels** (each a resolution level) internally, and most operations 
are applied lazily (or delayed) across all levels.

# Usage

```{r load-libs, message=FALSE, warning=FALSE}
library(ImageArray)
library(BiocFileCache)
library(RBioFormats)
library(EBImage)
library(magick)
library(ggplot2)
library(shiny)
```

## EBImage

We primarily use images via the `r Biocpkg("EBImage")` package. 

```{r read}
# 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 `r Biocpkg("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. 

```{r write}
output_h5 <- tempfile(fileext = ".h5")
imgarray <- writeImageArray(img, output = output_h5, n.levels = 2)
imgarray
```

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

```{r visualize_read}
bfa.raster <- as.raster(imgarray, level = 2)
plot(bfa.raster)
```

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

Here, by using the `max.pixel.size`, we request `r Biocpkg("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:

- (i) the entire image 
- (ii) a subset of the image with higher resolution 

with minimal memory footprint.

```{r visualize_read2}
# visualize
bfa.raster <- as.raster(imgarray, max.pixel.size = 400)
dim(bfa.raster)
```

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.

```{r rotate}
imgarray_rotated <- rotate(imgarray, angle = 90)
imgarray_rotated

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

```{r crop}
imgarray_cropped <- imgarray[300:500, 200:300]
imgarray_cropped

# 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
```

Cropping can also be performed with named indices, thus allowing users
to define indices for only the target dimensions. For now, 
`r Biocpkg("ImageArray")` only includes `x` and `y` dimensions and 'c' for
image channels. 

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

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)`).

```{r crop_named2}
# 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

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

## magick images

`r Biocpkg("ImageArray")` is also compatible with the `r CRANpkg("magick")` 
package where we retain the dimensionality of the image similar to 
`r Biocpkg("EBImage")`.

```{r read3}
# 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

# 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 `r Biocpkg("RBioFormats")`) with already defined layers. 

```{r ometiff}
ome.tiff.file <- system.file("extdata", "xy_12bit__plant.ome.tiff",
                             package = "ImageArray")
read.metadata(ome.tiff.file)

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

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

```{r ometiffvis2, out.width="50%"}
bfa.raster <- as.raster(imgarray, max.pixel.size = 300)
plot(bfa.raster)
dim(bfa.raster)

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

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

# Use cases

The delayed pyramid scheme introduced by `ImageArray` objects can also 
be used to generate scalable plots and interactive Shiny application 
for visualizing large images without loading them into the memory. 

## 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).

```{r he}
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)
```

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.

```{r he_imgarray}
imgarray <- createImageArray(image_file,
                             series = 1,
                             resolution = 1:6)
imgarray
```

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.

```{r speedy_vis, out.width="70%"}
# 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)

# 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`).

```{r speedy_vis_shiny}
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)
}
```

# Session info

```{r sessionInfo, echo=FALSE}
sessionInfo()
```