for plotImage, add argument to control contrast (cl), including doc examples, unit tests etc.

R/plotImage.R

@@ -12,9 +12,12 @@
 #'   800 x 800px; use Inf to plot the lowest resolution available.
 #' @param ch image channel(s) to be used for plotting (defaults to 
 #'   the first channel(s) available); use \code{channels()} to see 
-#'   which channels are available for a given \code{ImageArray} 
-#'   
-#' @param c plotting aesthetics; color 
+#'   which channels are available for a given \code{ImageArray}
+#' @param c character vector; colors to use for each channel. 
+#' @param cl list of length-2 numeric vectors (non-negative, increasing); 
+#'   specifies channel-wise contrast limits - defaults to [0, 1] for all 
+#'   (ignored when \code{image(x, i)} is an RGB image; 
+#'   for convenience, any NULL = [0, 1], and n = [0, n]).
 #'
 #' @return ggplot
 #'
@@ -28,6 +31,14 @@
 #'   plotImage(x, i=2, k=.))
 #' patchwork::wrap_plots(ms)
 #' 
+#' # custom colors
+#' cmy <- c("cyan", "magenta", "yellow")
+#' plotSpatialData() + plotImage(x, c=cmy)
+#' 
+#' # contrast limits
+#' plotSpatialData() + plotImage(x, c=cmy, 
+#'   cl=list(c(0.2,1), c(0,0.8), c(0,1)))
+#' 
 #' @import SpatialData
 NULL
 
@@ -36,22 +47,53 @@ NULL
 #' @export
 plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme 
 
+.check_cl <- \(cl, d) {
+    if (is.null(cl)) {
+        # default to [0, 1] for all channels
+        cl <- replicate(d, c(0, 1), FALSE)
+    } else {
+        # should be a list with as many elements as channels
+        if (!is.list(cl)) stop("'cl' should be a list")
+        if (length(cl) != d) stop("'cl' should be of length ", d)
+        for (. in seq_len(d)) {
+            # replace NULL by [0, 1] & n by [0, n]
+            if (is.null(cl[[.]])) cl[[.]] <- c(0, 1)
+            if (length(cl[[.]]) == 1) {
+                if (cl[[.]] < 0) stop("scalar 'cl' can't be < 0")
+                cl[[.]] <- c(0, cl[[.]])
+            }
+        }
+        # elements should be length-2, numeric, non-negative, increasing
+        .f <- \(.) length(.) == 2 && is.numeric(.) && all(. >= 0) && .[2] > .[1]
+        if (!all(vapply(cl, .f, logical(1))))
+            stop("elements of 'cl' should be length-2,",
+                " non-negative, increasing numeric vectors")
+    }
+    return(cl)
+}
+
 # merge/manage image channels
 # if no colors and channels defined, return the first channel
 #' @importFrom grDevices col2rgb
 #' @noRd
-.manage_channels <- \(a, ch, c=NULL){
+.chs2rgb <- \(a, ch, c=NULL, cl=NULL) {
+    cl <- .check_cl(cl, d <- dim(a)[1])
     if (length(ch) > (n <- length(.DEFAULT_COLORS)) && is.null(c))
         stop("Only ", n, " default colors available, but",
             length(ch), " are needed; please specify 'c'")
     if (!is.null(c) || (is.null(c) && length(ch) > 1)) {
         if (is.null(c)) c <- .DEFAULT_COLORS[seq_along(ch)] 
         c <- col2rgb(c)/255
         b <- array(0, dim=c(3, dim(a)[-1]))
-        for (i in seq_len(dim(a)[1])) {
-            b[1,,] <- b[1,,,drop=FALSE] + a[i,,,drop=FALSE]*c[1,i]
-            b[2,,] <- b[2,,,drop=FALSE] + a[i,,,drop=FALSE]*c[2,i]
-            b[3,,] <- b[3,,,drop=FALSE] + a[i,,,drop=FALSE]*c[3,i]
+        for (i in seq_len(d)) {
+            for (j in seq_len(3)) {
+                rgb <- a[i,,,drop=FALSE]*c[j,i]
+                # apply upper contrast lim.
+                rgb <- rgb*(1/cl[[i]][2]) 
+                b[j,,] <- b[j,,,drop=FALSE] + rgb
+                # apply lower contrast lim.
+                b[j,,][b[j,,] < cl[[i]][1]] <- 0
+            }
         }
         a <- pmin(b, 1)
     } else {
@@ -64,21 +106,22 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 #' @importFrom SpatialData getZarrArrayPath
 #' @importFrom Rarr zarr_overview
 #' @noRd
-.get_image_dtype <- \(a) {
+.get_img_dt <- \(a) {
     pa <- getZarrArrayPath(a)
     df <- zarr_overview(pa, as_data_frame=TRUE)
     if (!is.null(dt <- df$data_type)) return(dt)
 }
 
 # normalize the image data given its data type
 #' @noRd
-.normalize_image_array <- \(a, dt){
-  if (dt %in% names(.DTYPE_MAX_VALUES)) {
-      a <- a/.DTYPE_MAX_VALUES[dt]
-  } else if (max(a) > 1) {
-      for (i in seq_len(dim(a)[1])) 
-          a[i,,] <- a[i,,]/max(a[i,,])
-  }
+.norm_ia <- \(a, dt) {
+    d <- dim(a)[1]
+    if (dt %in% names(.DTYPE_MAX_VALUES)) {
+        a <- a / .DTYPE_MAX_VALUES[dt]
+    } else if (max(a) > 1) {
+        for (i in seq_len(d))
+            a[i,,] <- a[i,,] / max(a[i,,])
+    }
   return(a)
 }
 
@@ -124,24 +167,26 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
     which.min(d)
 }
 
-.get_plot_data <- \(x, k=NULL, w=800, h=800) {
+.get_img_data <- \(x, k=NULL, w=800, h=800) {
     if (!is.null(k)) return(data(x, k))
     data(x, .guess_scale(x, w, h))
 }
 
 #' @importFrom methods as
 #' @importFrom grDevices rgb
 #' @importFrom DelayedArray realize
-.df_i <- \(x, k=NULL, ch=NULL, c=NULL) {
-    a <- .get_plot_data(x, k)
-    ch_i <- .ch_idx(x, ch)
-    if (!.is_rgb(x))
-        a <- a[ch_i, , , drop=FALSE]
-    dt <- .get_image_dtype(a)
-    a <- realize(as(a, "DelayedArray"))
-    a <- .normalize_image_array(a, dt)
+.df_i <- \(x, k=NULL, ch=NULL, c=NULL, cl=NULL) {
+    a <- .get_img_data(x, k)
+    ch <- .ch_idx(x, ch)
     if (!.is_rgb(x))
-        a <- .manage_channels(a, ch_i, c)
+        a <- a[ch, , , drop=FALSE]
+    dt <- .get_img_dt(a)
+    a <- as(a, "DelayedArray")
+    a <- .norm_ia(realize(a), dt)
+    # enter when image isn't RGB already, either
+    # custom colors or contrasts are specified
+    if (!.is_rgb(x) || !is.null(c) || !is.null(cl))
+        a <- .chs2rgb(a, ch, c, cl)
     apply(a, c(2, 3), \(.) do.call(rgb, as.list(.))) 
 }
 
@@ -163,13 +208,13 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 
 #' @rdname plotImage
 #' @export
-setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=NULL) {
+setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=NULL, cl=NULL) {
     if (is.numeric(i))
         i <- imageNames(x)[i]
     y <- image(x, i)
     if (is.numeric(j))
         j <- CTname(y)[j]
-    df <- .df_i(y, k, ch, c)
+    df <- .df_i(y, k, ch, c, cl)
     wh <- .get_wh(x, i, j)
     .gg_i(df, wh$w, wh$h)
 })

inst/NEWS

@@ -1,6 +1,7 @@
 changes in version 0.99.2
 
-- in 'plotImage', added support to visualize channels of choice
+- in 'plotImage', added support to visualize channels of choice,
+  as well as 'cl' argument to control constrast limits
 - updated vignette to include corresponding examples
 
 changes in version 0.99.1

tests/testthat/test-plotArray.R

@@ -22,11 +22,11 @@ test_that(".guess_scale", {
             dim <- lapply(c(6, 3), \(.) c(3, rep(., 2))), \(.) 
             array(sample(seq_len(255), prod(.), replace=TRUE), dim=.)))
     # manual scale
-    expect_identical(.get_plot_data(img, k=1), lys[[1]]) 
-    expect_identical(.get_plot_data(img, k=2), lys[[2]])
+    expect_identical(.get_img_data(img, k=1), lys[[1]]) 
+    expect_identical(.get_img_data(img, k=2), lys[[2]])
     # automatic scale
-    expect_identical(.get_plot_data(img, k=NULL, w=5, h=7), lys[[1]]) 
-    expect_identical(.get_plot_data(img, k=NULL, w=2, h=2), lys[[2]])
+    expect_identical(.get_img_data(img, k=NULL, w=5, h=7), lys[[1]]) 
+    expect_identical(.get_img_data(img, k=NULL, w=2, h=2), lys[[2]])
 })
 
 test_that("plotImage()", {

tests/testthat/test-plotImage.R

@@ -24,7 +24,84 @@ test_that(".ch_idx()", {
   expect_warning(expect_equal(.ch_idx(image(x,1), ch=99), 1)) 
 })
 
-# TODO: any tests for image array normalization ? 
-test_that(".normalize_image_array", {
-    skip()
+test_that(".check_cl", {
+    # valid
+    n <- sample(seq(3, 9), 1)
+    v <- replicate(n, sort(runif(2)), FALSE)
+    expect_identical(.check_cl(v, n), v)
+    # single NULL
+    n <- sample(seq(3, 9), 1)
+    l <- .check_cl(NULL, n)
+    expect_is(l, "list")
+    expect_identical(l, replicate(n, c(0, 1), FALSE))
+    # one NULL, rest scalar
+    n <- sample(seq(3, 9), 1)
+    i <- sample(n, 1)
+    . <- replicate(n, NULL, FALSE)
+    .[[i]] <- v <- c(0.2, 0.8)
+    l <- .check_cl(., n)
+    expect_is(l, "list")
+    expect_identical(l[[i]], v)
+    expect_identical(l[-i], replicate(n-1, c(0, 1), FALSE))
+    # invalid
+    expect_error(.check_cl(c(0.2, 0.4, 0.6), 3)) # non-list
+    expect_error(.check_cl(as.list(seq_len(4)), 3)) # wrong length
+    expect_error(.check_cl(list(NULL, NULL, c(-1, 1)), 3)) # negative entry
+    expect_error(.check_cl(as.list(letters[seq_len(3)]), 3)) # non-numeric
+    expect_error(.check_cl(list(NULL, NULL, c(1, 0)), 3)) # decreasing
+    expect_error(.check_cl(list(NULL, NULL, -1), 3)) # negative scalar
+    expect_error(.check_cl(list(NULL, NULL, 0), 3)) # zero scalar
 })
+
+dir.create(td <- tempfile())
+pa <- unzip_spd_demo(zipname="steinbock_io.zip", dest=td, source="biocOSN")
+x <- readSpatialData(pa, images=1, labels=FALSE, points=FALSE, shapes=FALSE, tables=FALSE)
+a <- data(image(x)[seq_len(3), seq_len(100), seq_len(100)], 1)
+
+test_that(".get_img_dt", {
+    da <- (za <- data(image(x)))[1,,]
+    dt_za <- .get_img_dt(za) # from 'ZarrArray'
+    dt_da <- .get_img_dt(da) # from 'DelayedArray'
+    expect_is(dt_za, "character")
+    expect_identical(dt_za, dt_da)
+})
+
+test_that(".norm_ia", {
+    # valid data type
+    dt <- .get_img_dt(a)
+    b <- .norm_ia(realize(a), dt)
+    expect_equal(
+        apply(b, 1, range), 
+        replicate(3, c(0, 1)))
+    # invalid data type
+    b <- .norm_ia(realize(a), "")
+    expect_equal(
+        apply(b, 1, range), 
+        replicate(3, c(0, 1)))
+})
+
+test_that(".chs2rgb", {
+    dt <- .get_img_dt(a)
+    ch <- seq_len(d <- dim(a)[1])
+    a <- .norm_ia(realize(a), dt)
+    # no colors, no contrasts
+    b <- .chs2rgb(a, ch)
+    expect_identical(a, b)
+    # colors
+    cmy <- c("cyan", "magenta", "yellow")
+    b <- .chs2rgb(a, ch, c=cmy)
+    expect_equal(dim(a), dim(b))
+    expect_equal(
+        apply(b, 1, range), 
+        replicate(d, c(0, 1)))
+    # lower contrast lim.
+    lim <- list(c(0.5, 1), NULL, NULL)
+    b <- .chs2rgb(a, ch, cl=lim)
+    expect_identical(b[-1,,], a[-1,,])
+    expect_true(sum(b[1,,] == 0) > sum(a[1,,] == 0))
+    # upper contrast lim.
+    lim <- list(c(0, 0.5), NULL, NULL)
+    b <- .chs2rgb(a, ch, cl=lim)
+    fac <- mean(b[1,,]/a[1,,], na.rm=TRUE)
+    expect_equal(fac, 2, tolerance=0.05)
+})

vignettes/SpatialData.plot.Rmd

@@ -276,15 +276,30 @@ pa <- unzip_spd_demo(
 x <- readSpatialData(pa, anndataR=FALSE)
 ```
 
-Plotting with multiple image channels.
+### channels
 
-```{r steinbock-plot}
-plotSpatialData() + plotImage(x, 
+```{r steinbock-ch}
+plotSpatialData() + plotImage(x,
     i="Patient3_003_image",
     ch=c(6, 22, 39),
     c=c("blue", "cyan", "yellow"))
 ```
 
+### contrasts
+
+```{r steinbock-cl, fig.width=9, fig.height=3}
+i <- image(x, "Patient3_003_image")
+image(x, "crop") <- i[, 200:400, 200:400]
+lapply(list(c(0.2, 1), c(0, 0.8), c(0, 1.2)), \(.) {
+    plotSpatialData() + plotImage(x, 
+        i="crop", 
+        ch=c(6, 22, 39), 
+        cl=list(1, 1, .),
+        c=c("blue", "cyan", "yellow")) +
+        ggtitle(sprintf("[%s, %s]", .[1], .[2]))
+}) |> wrap_plots(nrow=1)
+```
+
 # Masking
 
 Back to blobs...