Add plot(GCIMSDataset) with a shared, comparable color scale

R/aaa-class-GCIMSDataset.R

@@ -79,6 +79,11 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
     dt_ref = NULL, # numeric or NULL
     #' @field rt_ref A numeric retention time of reference
     rt_ref = NULL, # numeric or NULL
+    #' @field intensity_range A matrix of n_samples x c(min, max) with each
+    #' sample's raw intensity range, extracted alongside TIS/RIC. Internal
+    #' only: no S4 generic exposes this, it only backs `plot(GCIMSDataset)`'s
+    #' shared color scale across samples.
+    intensity_range = NULL, # matrix or NULL
     #' @field userData A list to store arbitrary data in the dataset
     userData = list(), # list
     # Methods:
@@ -221,6 +226,7 @@ GCIMSDataset <- R6::R6Class("GCIMSDataset",
       self$RIC <- NULL
       self$dt_ref <- NULL
       self$rt_ref <- NULL
+      self$intensity_range <- NULL
       private$delayed_dataset$subset(samples)
       self
     },

R/getTIS_getRIC-GCIMSDataset.R

@@ -114,7 +114,11 @@ setMethod(
   ric <- intmat[ric_pos, ]
   ric <- max(ric) - ric
   ric <- ric/sum(ric)
-  list(ric = ric, tis = tis, rt = rt, dt = dt)
+  # intmat is already loaded here, so this rides along for free -- it backs
+  # plot(GCIMSDataset)'s shared intensity scale across samples (see
+  # dataset_intensity_range()), without a dedicated extraction pass.
+  intensity_range <- range(intmat)
+  list(ric = ric, tis = tis, rt = rt, dt = dt, intensity_range = intensity_range)
 }
 
 .extract_RIC_and_TIS_fun_aggregate <- function(ds, objs) {
@@ -123,6 +127,7 @@ setMethod(
   tiss <- purrr::map(objs, "tis")
   dtimes <- purrr::map(objs, "dt")
   rtimes <- purrr::map(objs, "rt")
+  intensity_ranges <- purrr::map(objs, "intensity_range")
 
   dt_ref <- ds$dt_ref
   rt_ref <- ds$rt_ref
@@ -144,9 +149,23 @@ setMethod(
   }
   stopifnot(nrow(ds$TIS) == num_samples)
   stopifnot(nrow(ds$RIC) == num_samples)
+  ds$intensity_range <- do.call(rbind, intensity_ranges)
+  colnames(ds$intensity_range) <- c("min", "max")
   ds
 }
 
+#' The dataset's overall raw intensity range, realizing/refreshing the cache if needed
+#' @param object A [GCIMSDataset] object
+#' @return A length-2 numeric vector `c(min, max)`
+#' @noRd
+dataset_intensity_range <- function(object) {
+  if (object$hasDelayedOps() || is.null(object$intensity_range)) {
+    object$extract_RIC_and_TIS()
+    object$realize()
+  }
+  range(object$intensity_range)
+}
+
 #' Extract the Reverse Ion Chromatogram and Total Ion Spectrum from the samples
 #'
 #' @param object A GCIMSDataset object

R/plot-GCIMSDataset.R

@@ -0,0 +1,138 @@
+#' Topographical plot of several samples of a GCIMSDataset, on a shared color scale
+#'
+#' Several samples are plotted side by side, on the same intensity color
+#' scale, so they are visually comparable. Paginated: only one page's worth
+#' of samples is ever loaded into memory at a time for rendering.
+#'
+#' @param x A [GCIMSDataset] object
+#' @param sample A number, a string, or a vector of numbers/strings with the
+#' sample index(es) or name(s) to plot. If `NULL` (the default), all samples
+#' are plotted (across as many pages as needed).
+#' @inheritParams dt_rt_range_normalization
+#' @param ... Ignored
+#' @param remove_baseline Set to `TRUE` to subtract each sample's estimated
+#' baseline first
+#' @param trans The transformation to the intensity values, see `plot,GCIMSSample-method`
+#' @param intensity_range Controls the shared color scale across all plotted
+#' samples (all pages, not just the current one). One of:
+#' - `"global"` (the default): each sample's full, uncropped intensity range,
+#'   regardless of `sample`/`dt_range`/`rt_range`. With `remove_baseline =
+#'   FALSE`, this comes from a cache computed once, alongside TIS/RIC, so
+#'   it's free. With `remove_baseline = TRUE`, the cache (raw intensities)
+#'   doesn't apply, so each selected sample's full, uncropped, baseline-removed
+#'   range is computed instead (loading every selected sample's data an extra
+#'   time).
+#' - `"ranged"`: the range of exactly what's plotted (respecting `dt_range`,
+#'   `rt_range` and `remove_baseline`), across *every* selected sample, not
+#'   just the current page -- so the scale is still comparable page to page.
+#'   This is deliberately not page-scoped: a cheaper, page-local range would
+#'   make pages incomparable to each other, defeating the purpose of paging
+#'   through a dataset to compare samples. Costs loading every selected
+#'   sample's data twice (once to compute the range, once to plot).
+#' - A numeric vector of length 2, `c(min, max)`: fixed limits.
+#' - A length-2 vector/list whose elements are independently a number,
+#'   `"global"` or `"ranged"`, e.g. `list(min = 0, max = "global")`.
+#' @param nrow,ncol Page grid shape, forwarded to [cowplot::plot_grid()]. Page
+#' capacity is `nrow * ncol`. If both are `NULL` (the default), they're picked
+#' from the number of selected samples: an exact fit for 1-6 samples (e.g. 5
+#' or 6 samples -> 2x3), otherwise 3x3 (so page capacity maxes out at 9, and
+#' more than 9 selected samples span multiple pages). If only one of
+#' `nrow`/`ncol` is given, the other becomes `3` when there are more than 9
+#' selected samples, or just enough to fit them all on one page otherwise.
+#' @param page Which page to plot, 1-based. Errors if out of bounds.
+#' @return A combined plot with one panel per sample, for the requested page
+#' @export
+setMethod(
+  "plot",
+  "GCIMSDataset",
+  function(x, sample = NULL, dt_range = NULL, rt_range = NULL, ...,
+           remove_baseline = FALSE, trans = "cubic_root",
+           intensity_range = "global", nrow = NULL, ncol = NULL, page = 1) {
+    require_pkgs("cowplot")
+    sample_names_all <- sampleNames(x)
+    if (is.null(sample)) {
+      sample <- sample_names_all
+    }
+    selected_names <- sample_name_or_number_to_both(sample, sample_names_all)$name
+    num_samples <- length(selected_names)
+    if (num_samples == 0) {
+      cli_abort("No samples to plot")
+    }
+
+    grid_dims <- resolve_page_grid(nrow, ncol, num_samples)
+    nrow <- grid_dims[["nrow"]]
+    ncol <- grid_dims[["ncol"]]
+    page_size <- nrow * ncol
+    n_pages <- ceiling(num_samples / page_size)
+    if (page < 1 || page > n_pages) {
+      cli_abort(
+        c(
+          "{.arg page} = {page} is out of bounds",
+          "i" = "With {num_samples} sample(s) laid out on a {nrow}x{ncol} grid, there are {n_pages} page(s)"
+        )
+      )
+    }
+    page_start <- (page - 1) * page_size + 1
+    page_end <- min(page * page_size, num_samples)
+    page_names <- selected_names[page_start:page_end]
+
+    cropped_intensity <- function(s) {
+      dt <- dtime(s)
+      rt <- rtime(s)
+      idx <- dt_rt_range_normalization(dt, rt, dt_range, rt_range)
+      intmat <- intensity(s, idx)
+      if (isTRUE(remove_baseline)) {
+        basel <- baseline(s)[idx[["dt_logical"]], idx[["rt_logical"]]]
+        intmat <- intmat - basel
+      }
+      intmat
+    }
+
+    ranged_range_cache <- NULL
+    get_ranged_range <- function() {
+      if (is.null(ranged_range_cache)) {
+        # Scoped to every selected sample (all pages), not just the current
+        # one, so the scale stays comparable across pages. Each sample is
+        # loaded, ranged and discarded one at a time to avoid holding them
+        # all in RAM at once.
+        ranges <- purrr::map(selected_names, function(nm) {
+          range(cropped_intensity(x$getSample(nm)))
+        })
+        ranged_range_cache <<- range(unlist(ranges))
+      }
+      ranged_range_cache
+    }
+    global_range_cache <- NULL
+    get_global_range <- function() {
+      if (is.null(global_range_cache)) {
+        global_range_cache <<- if (isTRUE(remove_baseline)) {
+          ranges <- purrr::map(selected_names, function(nm) {
+            s <- x$getSample(nm)
+            range(intensity(s) - baseline(s))
+          })
+          range(unlist(ranges))
+        } else {
+          dataset_intensity_range(x)
+        }
+      }
+      global_range_cache
+    }
+
+    limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
+
+    # Only the current page's samples are loaded for rendering, one at a
+    # time, so peak memory never holds more than one raw sample plus the
+    # (much smaller, native-raster-encoded) panels built so far.
+    panels <- purrr::map(page_names, function(nm) {
+      s <- x$getSample(nm)
+      plot(
+        s,
+        dt_range = dt_range, rt_range = rt_range,
+        remove_baseline = remove_baseline, trans = trans,
+        intensity_range = limits
+      ) + ggplot2::labs(title = nm)
+    })
+
+    cowplot::plot_grid(plotlist = panels, nrow = nrow, ncol = ncol)
+  }
+)

R/plot-GCIMSSample.R

@@ -7,6 +7,17 @@
 #' @param remove_baseline Set to `TRUE` to subtract the estimated baseline first
 #' @param trans The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 #' See the `trans` argument in [ggplot2::continuous_scale()] for other possibilities.
+#' @param intensity_range Controls the color scale limits. One of:
+#' - `"ranged"` (the default): the range of exactly what's plotted (respecting
+#'   `dt_range`, `rt_range` and `remove_baseline`).
+#' - `"global"`: this sample's own full, uncropped intensity range, regardless
+#'   of `dt_range`/`rt_range` -- useful to see how strong a cropped region is
+#'   relative to the whole sample. With `remove_baseline = TRUE`, this is the
+#'   full, uncropped, baseline-removed range.
+#' - A numeric vector of length 2, `c(min, max)`: fixed limits. Set this
+#'   explicitly to make several plots comparable on the same color scale.
+#' - A length-2 vector/list whose elements are independently a number,
+#'   `"global"` or `"ranged"`, e.g. `list(min = 0, max = "global")`.
 #' @return A plot of the GCIMSSample
 #' @examples
 #' dummy_obj <-GCIMSSample(
@@ -22,7 +33,7 @@
 setMethod(
   "plot",
   "GCIMSSample",
-  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root") {
+  function(x, dt_range = NULL, rt_range = NULL, ..., remove_baseline = FALSE, trans = "cubic_root", intensity_range = "ranged") {
   dt <- dtime(x)
   rt <- rtime(x)
   idx <- dt_rt_range_normalization(dt, rt, dt_range, rt_range)
@@ -32,18 +43,31 @@ setMethod(
     basel <- baseline(x)[idx$dt_logical, idx$rt_logical]
     intmat <- intmat - basel
   }
+
+  get_global_range <- function() {
+    if (isTRUE(remove_baseline)) {
+      range(intensity(x) - baseline(x))
+    } else {
+      range(intensity(x))
+    }
+  }
+  get_ranged_range <- function() range(intmat)
+
+  limits <- resolve_intensity_range(intensity_range, get_global_range, get_ranged_range)
+
   mat_to_gplot(
     intmat,
     dt_min = idx$dt_ms_min,
     dt_max = idx$dt_ms_max,
     rt_min = idx$rt_s_min,
     rt_max = idx$rt_s_max,
-    trans = trans
+    trans = trans,
+    intensity_range = limits
   )
 })
 
 
-mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root") {
+mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt_max = NULL, trans = "cubic_root", intensity_range = NULL) {
   require_pkgs(c("farver", "viridisLite"))
   if (is.null(dt_min)) {
     dt_min <- as.numeric(rownames(intmat)[1L])
@@ -57,7 +81,7 @@ mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt
   if (is.null(rt_max)) {
     rt_max <- as.numeric(colnames(intmat)[ncol(intmat)])
   }
-  minmax <- range(intmat)
+  minmax <- if (is.null(intensity_range)) range(intmat) else intensity_range
 
   if (is.character(trans)) {
     trans_func <- paste0(trans, "_trans")
@@ -76,7 +100,7 @@ mat_to_gplot <- function(intmat, dt_min = NULL, dt_max = NULL, rt_min = NULL, rt
   colormap <- farver::encode_native(
     viridisLite::viridis(256L, direction = -1, option = "G")
   )
-  nr <- mat_to_nativeRaster(intmat_trans, colormap)
+  nr <- mat_to_nativeRaster(intmat_trans, colormap, rangex = trans$transform(minmax))
 
   # The geom_rect is fake and it is only used to force the fill legend to appear
   # The geom_rect  limits are used to help set the plot limits

R/utils-plot.R

@@ -35,6 +35,87 @@ mat_to_nativeRaster <- function(x, colormap, rangex = NULL)  {
 }
 
 
+#' Resolve the (`nrow`, `ncol`) page grid for `plot,GCIMSDataset-method`,
+#' auto-filling whichever of `nrow`/`ncol` isn't given
+#'
+#' @param nrow Number of rows, or `NULL`
+#' @param ncol Number of columns, or `NULL`
+#' @param num_samples Number of samples to lay out (before pagination)
+#' @return A list `list(nrow =, ncol =)`
+#' @noRd
+resolve_page_grid <- function(nrow, ncol, num_samples) {
+  if (is.null(nrow) && is.null(ncol)) {
+    dims <- if (num_samples <= 1) {
+      c(1L, 1L)
+    } else if (num_samples == 2) {
+      c(1L, 2L)
+    } else if (num_samples == 3) {
+      c(1L, 3L)
+    } else if (num_samples == 4) {
+      c(2L, 2L)
+    } else if (num_samples <= 6) {
+      c(2L, 3L)
+    } else {
+      c(3L, 3L)
+    }
+    return(list(nrow = dims[1], ncol = dims[2]))
+  }
+  if (is.null(nrow)) {
+    nrow <- if (num_samples > 9) 3L else ceiling(num_samples / ncol)
+    return(list(nrow = nrow, ncol = ncol))
+  }
+  if (is.null(ncol)) {
+    ncol <- if (num_samples > 9) 3L else ceiling(num_samples / nrow)
+    return(list(nrow = nrow, ncol = ncol))
+  }
+  list(nrow = nrow, ncol = ncol)
+}
+
+#' Resolve an `intensity_range` argument (shared by `plot,GCIMSSample-method`
+#' and `plot,GCIMSDataset-method`) into `c(min, max)`
+#'
+#' @param intensity_range One of `"global"`, `"ranged"`, a numeric vector of
+#' length 2, or a length-2 vector/list whose elements are independently a
+#' number, `"global"` or `"ranged"`.
+#' @param global_range A zero-argument function returning the full/uncropped
+#' `c(min, max)`. Only called if actually referenced.
+#' @param ranged_range A zero-argument function returning the `c(min, max)`
+#' of exactly what will be plotted. Only called if actually referenced.
+#' @return A numeric vector `c(min, max)`
+#' @noRd
+resolve_intensity_range <- function(intensity_range, global_range, ranged_range) {
+  resolve_endpoint <- function(value, position) {
+    if (is.numeric(value) && length(value) == 1) {
+      return(value)
+    }
+    if (identical(value, "global")) {
+      return(global_range()[position])
+    }
+    if (identical(value, "ranged")) {
+      return(ranged_range()[position])
+    }
+    cli_abort('Each element of intensity_range should be a number, "global" or "ranged"')
+  }
+
+  if (is.character(intensity_range) && length(intensity_range) == 1) {
+    if (intensity_range == "global") {
+      return(global_range())
+    } else if (intensity_range == "ranged") {
+      return(ranged_range())
+    }
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  if (length(intensity_range) != 2) {
+    cli_abort('intensity_range should be "global", "ranged", or a length-2 vector/list of numbers, "global" or "ranged"')
+  }
+
+  c(
+    resolve_endpoint(intensity_range[[1]], 1L),
+    resolve_endpoint(intensity_range[[2]], 2L)
+  )
+}
+
 #' Make a plot interactive
 #'
 #' Wraps the `plt` with [plotly::ggplotly()] and sets the `xaxis` and `yaxis`