diff --git a/R/aaa-class-GCIMSDataset.R b/R/aaa-class-GCIMSDataset.R
index d367742..8cc5834 100644
--- a/R/aaa-class-GCIMSDataset.R
+++ b/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
     },
diff --git a/R/getTIS_getRIC-GCIMSDataset.R b/R/getTIS_getRIC-GCIMSDataset.R
index 1cfa79b..0d7c372 100644
--- a/R/getTIS_getRIC-GCIMSDataset.R
+++ b/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
diff --git a/R/plot-GCIMSDataset.R b/R/plot-GCIMSDataset.R
new file mode 100644
index 0000000..9e8f0d3
--- /dev/null
+++ b/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)
+  }
+)
diff --git a/R/plot-GCIMSSample.R b/R/plot-GCIMSSample.R
index a3c9e3e..42c6ec5 100644
--- a/R/plot-GCIMSSample.R
+++ b/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
diff --git a/R/utils-plot.R b/R/utils-plot.R
index 41559e5..326b456 100644
--- a/R/utils-plot.R
+++ b/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`
diff --git a/man/GCIMSDataset.Rd b/man/GCIMSDataset.Rd
index bbe2408..5e2ca40 100644
--- a/man/GCIMSDataset.Rd
+++ b/man/GCIMSDataset.Rd
@@ -96,6 +96,11 @@ dummy_dataset <- GCIMSDataset$new(
 
 \item{\code{rt_ref}}{A numeric retention time of reference}
 
+\item{\code{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 \code{plot(GCIMSDataset)}'s
+shared color scale across samples.}
+
 \item{\code{userData}}{A list to store arbitrary data in the dataset}
 }
 \if{html}{\out{</div>}}
diff --git a/man/plot-GCIMSDataset-ANY-method.Rd b/man/plot-GCIMSDataset-ANY-method.Rd
new file mode 100644
index 0000000..6685ebe
--- /dev/null
+++ b/man/plot-GCIMSDataset-ANY-method.Rd
@@ -0,0 +1,77 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot-GCIMSDataset.R
+\name{plot,GCIMSDataset,ANY-method}
+\alias{plot,GCIMSDataset,ANY-method}
+\title{Topographical plot of several samples of a GCIMSDataset, on a shared color scale}
+\usage{
+\S4method{plot}{GCIMSDataset,ANY}(
+  x,
+  sample = NULL,
+  dt_range = NULL,
+  rt_range = NULL,
+  ...,
+  remove_baseline = FALSE,
+  trans = "cubic_root",
+  intensity_range = "global",
+  nrow = NULL,
+  ncol = NULL,
+  page = 1
+)
+}
+\arguments{
+\item{x}{A \link{GCIMSDataset} object}
+
+\item{sample}{A number, a string, or a vector of numbers/strings with the
+sample index(es) or name(s) to plot. If \code{NULL} (the default), all samples
+are plotted (across as many pages as needed).}
+
+\item{dt_range}{The minimum and maximum drift times to extract (length 2 vector)}
+
+\item{rt_range}{The minimum and maximum retention times to extract (length 2 vector)}
+
+\item{...}{Ignored}
+
+\item{remove_baseline}{Set to \code{TRUE} to subtract each sample's estimated
+baseline first}
+
+\item{trans}{The transformation to the intensity values, see \verb{plot,GCIMSSample-method}}
+
+\item{intensity_range}{Controls the shared color scale across all plotted
+samples (all pages, not just the current one). One of:
+\itemize{
+\item \code{"global"} (the default): each sample's full, uncropped intensity range,
+regardless of \code{sample}/\code{dt_range}/\code{rt_range}. With \code{remove_baseline = FALSE}, this comes from a cache computed once, alongside TIS/RIC, so
+it's free. With \code{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).
+\item \code{"ranged"}: the range of exactly what's plotted (respecting \code{dt_range},
+\code{rt_range} and \code{remove_baseline}), across \emph{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).
+\item A numeric vector of length 2, \code{c(min, max)}: fixed limits.
+\item A length-2 vector/list whose elements are independently a number,
+\code{"global"} or \code{"ranged"}, e.g. \code{list(min = 0, max = "global")}.
+}}
+
+\item{nrow, ncol}{Page grid shape, forwarded to \code{\link[cowplot:plot_grid]{cowplot::plot_grid()}}. Page
+capacity is \code{nrow * ncol}. If both are \code{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
+\code{nrow}/\code{ncol} is given, the other becomes \code{3} when there are more than 9
+selected samples, or just enough to fit them all on one page otherwise.}
+
+\item{page}{Which page to plot, 1-based. Errors if out of bounds.}
+}
+\value{
+A combined plot with one panel per sample, for the requested page
+}
+\description{
+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.
+}
diff --git a/man/plot-GCIMSSample-ANY-method.Rd b/man/plot-GCIMSSample-ANY-method.Rd
index 49194da..c3fd336 100644
--- a/man/plot-GCIMSSample-ANY-method.Rd
+++ b/man/plot-GCIMSSample-ANY-method.Rd
@@ -10,7 +10,8 @@
   rt_range = NULL,
   ...,
   remove_baseline = FALSE,
-  trans = "cubic_root"
+  trans = "cubic_root",
+  intensity_range = "ranged"
 )
 }
 \arguments{
@@ -26,6 +27,20 @@
 
 \item{trans}{The transformation to the intensity values. "cubic_root" is the default. "intensity" is also valid.
 See the \code{trans} argument in \code{\link[ggplot2:continuous_scale]{ggplot2::continuous_scale()}} for other possibilities.}
+
+\item{intensity_range}{Controls the color scale limits. One of:
+\itemize{
+\item \code{"ranged"} (the default): the range of exactly what's plotted (respecting
+\code{dt_range}, \code{rt_range} and \code{remove_baseline}).
+\item \code{"global"}: this sample's own full, uncropped intensity range, regardless
+of \code{dt_range}/\code{rt_range} -- useful to see how strong a cropped region is
+relative to the whole sample. With \code{remove_baseline = TRUE}, this is the
+full, uncropped, baseline-removed range.
+\item A numeric vector of length 2, \code{c(min, max)}: fixed limits. Set this
+explicitly to make several plots comparable on the same color scale.
+\item A length-2 vector/list whose elements are independently a number,
+\code{"global"} or \code{"ranged"}, e.g. \code{list(min = 0, max = "global")}.
+}}
 }
 \value{
 A plot of the GCIMSSample
diff --git a/tests/testthat/test-plot-GCIMSDataset.R b/tests/testthat/test-plot-GCIMSDataset.R
new file mode 100644
index 0000000..5e16915
--- /dev/null
+++ b/tests/testthat/test-plot-GCIMSDataset.R
@@ -0,0 +1,202 @@
+# --- resolve_intensity_range() -----------------------------------------
+
+make_counting_range <- function(value) {
+  n_calls <- 0L
+  fn <- function() {
+    n_calls <<- n_calls + 1L
+    value
+  }
+  attr(fn, "n_calls") <- function() n_calls
+  fn
+}
+
+test_that("resolve_intensity_range('global'/'ranged') calls only the matching closure, once", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range("global", global_range, ranged_range), c(0, 10))
+  expect_equal(attr(global_range, "n_calls")(), 1L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+})
+
+test_that("resolve_intensity_range() with a fixed numeric vector never calls either closure", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range(c(1, 9), global_range, ranged_range), c(1, 9))
+  expect_equal(attr(global_range, "n_calls")(), 0L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+})
+
+test_that("resolve_intensity_range() resolves mixed list(min=, max=) elements independently", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_equal(resolve_intensity_range(list(min = 1, max = "global"), global_range, ranged_range), c(1, 10))
+  expect_equal(attr(global_range, "n_calls")(), 1L)
+  expect_equal(attr(ranged_range, "n_calls")(), 0L)
+
+  expect_equal(resolve_intensity_range(list("ranged", 99), global_range, ranged_range), c(2, 99))
+  expect_equal(attr(ranged_range, "n_calls")(), 1L)
+})
+
+test_that("resolve_intensity_range() rejects invalid specs", {
+  global_range <- make_counting_range(c(0, 10))
+  ranged_range <- make_counting_range(c(2, 8))
+
+  expect_error(resolve_intensity_range("nope", global_range, ranged_range), "intensity_range should be")
+  expect_error(resolve_intensity_range(c(1, 2, 3), global_range, ranged_range), "intensity_range should be")
+  expect_error(resolve_intensity_range(list("nope", 1), global_range, ranged_range), "should be a number")
+})
+
+# --- resolve_page_grid() --------------------------------------------------
+
+test_that("resolve_page_grid() with neither nrow nor ncol given follows the exact-fit-then-3x3-cap table", {
+  expected <- list(
+    `1` = c(1, 1), `2` = c(1, 2), `3` = c(1, 3), `4` = c(2, 2),
+    `5` = c(2, 3), `6` = c(2, 3), `7` = c(3, 3), `9` = c(3, 3), `25` = c(3, 3)
+  )
+  for (n in names(expected)) {
+    g <- resolve_page_grid(NULL, NULL, as.integer(n))
+    expect_equal(c(g$nrow, g$ncol), expected[[n]], info = n)
+  }
+})
+
+test_that("resolve_page_grid() with exactly one of nrow/ncol given fills the other", {
+  # num_samples <= 9: fit everyone on one page
+  expect_equal(resolve_page_grid(NULL, 2, 5)$nrow, 3)
+  expect_equal(resolve_page_grid(4, NULL, 9)$ncol, 3)
+
+  # num_samples > 9: missing dimension is fixed at 3, regardless of the given one
+  expect_equal(resolve_page_grid(NULL, 5, 20)$nrow, 3)
+  expect_equal(resolve_page_grid(7, NULL, 50)$ncol, 3)
+})
+
+test_that("resolve_page_grid() with both given uses them as-is", {
+  g <- resolve_page_grid(4, 5, 100)
+  expect_equal(c(g$nrow, g$ncol), c(4, 5))
+})
+
+# --- intensity_range caching ----------------------------------------------
+
+make_range_dataset <- function() {
+  s1 <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
+  s2 <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(100:124, nrow = 5))
+  GCIMSDataset$new_from_list(samples = list(s1 = s1, s2 = s2), on_ram = TRUE, scratch_dir = NULL)
+}
+
+test_that("realizing a dataset populates intensity_range alongside TIS/RIC", {
+  old_bpparam <- BiocParallel::bpparam()
+  BiocParallel::register(BiocParallel::SerialParam())
+  on.exit(BiocParallel::register(old_bpparam))
+
+  ds <- make_range_dataset()
+  ds$realize()
+
+  expect_equal(unname(ds$intensity_range["s1", ]), c(1, 25))
+  expect_equal(unname(ds$intensity_range["s2", ]), c(100, 124))
+})
+
+test_that("subsetting a dataset resets the cached intensity_range", {
+  old_bpparam <- BiocParallel::bpparam()
+  BiocParallel::register(BiocParallel::SerialParam())
+  on.exit(BiocParallel::register(old_bpparam))
+
+  ds <- make_range_dataset()
+  ds$realize()
+  expect_false(is.null(ds$intensity_range))
+
+  ds$subset("s1", inplace = TRUE)
+
+  expect_null(ds$intensity_range)
+})
+
+# --- plot(GCIMSDataset) ----------------------------------------------------
+
+test_that("plot(GCIMSDataset) returns a combined plot using the cached global range by default", {
+  ds <- make_range_dataset()
+
+  p <- plot(ds)
+
+  expect_s3_class(p, "ggplot")
+})
+
+test_that("plot(GCIMSDataset, sample =) restricts which samples are plotted", {
+  ds <- make_range_dataset()
+
+  expect_no_error(plot(ds, sample = "s1"))
+  expect_no_error(plot(ds, sample = 1))
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'global') computes a baseline-removed range with remove_baseline = TRUE", {
+  ds <- make_range_dataset()
+  ds <- estimateBaseline(ds, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+  ds$realize()
+
+  expect_no_error(plot(ds, remove_baseline = TRUE))
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'ranged') works with remove_baseline = TRUE", {
+  ds <- make_range_dataset()
+  ds <- estimateBaseline(ds, dt_peak_fwhm_ms = 1, dt_region_multiplier = 2, rt_length_s = 2)
+  ds$realize()
+
+  expect_no_error(plot(ds, remove_baseline = TRUE, intensity_range = "ranged"))
+})
+
+test_that("plot(GCIMSDataset, intensity_range = c(min, max)) accepts fixed limits without touching the dataset", {
+  ds <- make_range_dataset()
+
+  expect_no_error(plot(ds, intensity_range = c(0, 200)))
+})
+
+# --- pagination -------------------------------------------------------------
+
+make_paged_dataset <- function(n = 11) {
+  samples <- stats::setNames(
+    purrr::map(seq_len(n), function(i) {
+      GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(rep(i, 25), nrow = 5))
+    }),
+    paste0("s", seq_len(n))
+  )
+  GCIMSDataset$new_from_list(samples = samples, on_ram = TRUE, scratch_dir = NULL)
+}
+
+test_that("plot(GCIMSDataset) with more than 9 samples defaults to a 3x3 first page", {
+  ds <- make_paged_dataset(11)
+
+  expect_no_error(plot(ds))
+  expect_no_error(plot(ds, page = 1))
+  expect_no_error(plot(ds, page = 2))
+})
+
+test_that("plot(GCIMSDataset, page =) errors clearly when out of bounds", {
+  ds <- make_paged_dataset(11)
+
+  expect_error(plot(ds, page = 3), "page.*3.*out of bounds")
+  expect_error(plot(ds, page = 0), "page.*0.*out of bounds")
+})
+
+test_that("plot(GCIMSDataset, intensity_range = 'ranged') is computed over every selected sample, not just the current page", {
+  # 11 samples, values 1..11 -> page 1 (samples 1-9) would only see 1..9 if
+  # ranged were page-scoped, but sample 11 (value 11) is on page 2.
+  ds <- make_paged_dataset(11)
+
+  captured <- list()
+  testthat::local_mocked_bindings(
+    mat_to_gplot = function(intmat, ..., intensity_range = NULL) {
+      captured[[length(captured) + 1]] <<- intensity_range
+      ggplot2::ggplot()
+    },
+    .package = "GCIMS"
+  )
+
+  plot(ds, intensity_range = "ranged", page = 1)
+  range_page1 <- captured[[1]]
+  captured <<- list()
+  plot(ds, intensity_range = "ranged", page = 2)
+  range_page2 <- captured[[1]]
+
+  expect_equal(range_page1, c(1, 11))
+  expect_equal(range_page2, c(1, 11))
+})
diff --git a/tests/testthat/test-plot-GCIMSSample.R b/tests/testthat/test-plot-GCIMSSample.R
index 4687f3f..fd3459f 100644
--- a/tests/testthat/test-plot-GCIMSSample.R
+++ b/tests/testthat/test-plot-GCIMSSample.R
@@ -81,6 +81,55 @@ test_that("plot(remove_baseline = TRUE) on a GCIMSSample does not error once a b
   expect_s3_class(p, "ggplot")
 })
 
+test_that("plot(intensity_range = ) overrides the auto-computed color scale limits", {
+  s <- GCIMSSample(drift_time = 1:5, retention_time = 1:5, data = matrix(1:25, nrow = 5))
+
+  p_auto <- plot(s)
+  p_fixed <- plot(s, intensity_range = c(0, 200))
+
+  tr <- cubic_root_trans()
+  expect_equal(p_auto$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 25)))
+  expect_equal(p_fixed$scales$get_scales("fill")$get_limits(), tr$transform(c(0, 200)))
+})
+
+test_that("plot(intensity_range = 'ranged') is the default and scales to the cropped view", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p_cropped <- plot(s, dt_range = c(1, 1))
+
+  tr <- cubic_root_trans()
+  expect_equal(p_cropped$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 91)))
+})
+
+test_that("plot(intensity_range = 'global') scales to the sample's full range, ignoring dt_range/rt_range", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p_global <- plot(s, dt_range = c(1, 1), intensity_range = "global")
+
+  tr <- cubic_root_trans()
+  expect_equal(p_global$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 100)))
+})
+
+test_that("plot(intensity_range = 'global') computes the full, uncropped, baseline-removed range with remove_baseline = TRUE", {
+  s <- make_sample()
+  s <- estimateBaseline(s, dt_peak_fwhm_ms = 0.3, dt_region_multiplier = 6, rt_length_s = 10, remove = FALSE)
+
+  p <- plot(s, dt_range = c(0, 1), remove_baseline = TRUE, intensity_range = "global")
+
+  expected_range <- range(intensity(s) - baseline(s))
+  tr <- cubic_root_trans()
+  expect_equal(p$scales$get_scales("fill")$get_limits(), tr$transform(expected_range))
+})
+
+test_that("plot(intensity_range = list(min=, max=)) resolves each endpoint independently", {
+  s <- GCIMSSample(drift_time = 1:10, retention_time = 1:10, data = matrix(1:100, nrow = 10))
+
+  p <- plot(s, dt_range = c(1, 1), intensity_range = list(min = "ranged", max = "global"))
+
+  tr <- cubic_root_trans()
+  expect_equal(p$scales$get_scales("fill")$get_limits(), tr$transform(c(1, 100)))
+})
+
 test_that("plot() on a GCIMSSample accepts trans as a string, or as a transform object directly", {
   s <- make_sample()