diff --git a/.Rbuildignore b/.Rbuildignore
index c503c4f..d325fbf 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1 +1,3 @@
 ^\.github$
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..52c1c38
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+*.Rproj
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index 5960769..e35dbac 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -41,7 +41,8 @@ Imports:
     rlang,
     sf,
     S4Vectors,
-    SingleCellExperiment
+    SingleCellExperiment,
+    Rarr
 Suggests:
     BiocStyle,
     ggnewscale,
diff --git a/NAMESPACE b/NAMESPACE
index c8ecce5..e5188d0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,9 +10,11 @@ exportMethods(scale)
 exportMethods(translation)
 import(SpatialData)
 importFrom(DelayedArray,realize)
+importFrom(Rarr,zarr_overview)
 importFrom(S4Vectors,metadata)
 importFrom(SingleCellExperiment,int_colData)
 importFrom(SingleCellExperiment,int_metadata)
+importFrom(SpatialData,getZarrArrayPath)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(ggforce,geom_circle)
diff --git a/R/plotImage.R b/R/plotImage.R
index 77d3b54..2f7bd3d 100644
--- a/R/plotImage.R
+++ b/R/plotImage.R
@@ -10,6 +10,8 @@
 #' @param k index of the scale of an image; by default (NULL), will auto-select 
 #'   scale in order to minimize memory-usage and blurring for a target size of 
 #'   800 x 800px; use Inf to plot the lowest resolution available.
+#' @param ch the image channels to be used for plotting (default: first channel)
+#' @param c plotting aesthetics; color 
 #'
 #' @return ggplot
 #'
@@ -31,6 +33,89 @@ NULL
 #' @export
 plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme 
 
+# merge/manage image channels
+# if no colors and channels defined, return the first channel
+#' @noRd
+.manage_channels <- \(a, ch, c=NULL){
+  if(length(ch) > length(.DEFAULT_COLORS) && is.null(c))
+    stop("You can only choose at most seven default colors!")
+  if(!is.null(c) || (is.null(c) && length(ch) > 1)) {
+    if(is.null(c))
+      c <- .DEFAULT_COLORS[1:length(ch)] 
+    c <- col2rgb(c)/255
+    a_new <- array(0, dim = c(3,dim(a)[-1]))
+    for(i in 1:dim(a)[1]){
+      a_new[1,,] <-  a_new[1,,,drop = FALSE] + a[i,,,drop = FALSE]*c[1,i]
+      a_new[2,,] <-  a_new[2,,,drop = FALSE] + a[i,,,drop = FALSE]*c[2,i]
+      a_new[3,,] <-  a_new[3,,,drop = FALSE] + a[i,,,drop = FALSE]*c[3,i]
+    }
+    a <- pmin(a_new,1)
+  } else {
+    a <- a[rep(1,3),,]
+  }
+  a
+}
+
+# check if an image is rgb or not
+#' @importFrom SpatialData getZarrArrayPath
+#' @importFrom Rarr zarr_overview
+#' @noRd
+.get_image_dtype <- \(a){
+  zarray_spec <- Rarr::zarr_overview(getZarrArrayPath(a), 
+                                     as_data_frame = TRUE)
+  if("data_type" %in% names(zarray_spec))
+    return(zarray_spec$data_type)
+  return(NULL)
+}
+
+# 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 1:dim(a)[1])
+      a[i,,] <- a[i,,]/max(a[i,,])
+  }
+  a 
+}
+
+# check if an image is rgb or not
+# NOTE: some rgb channels are named as 0:2
+#' @noRd
+.is.rgb <- \(x){
+  if(!is.null(md <- x@meta))
+    labels <- md[[2]]$channels$label
+  if(length(labels) == 3)
+    if(all(labels %in% c("r", "g", "b")) || all(labels %in% seq(0,2))) {
+      return(TRUE)
+    }
+  return(FALSE)
+}
+
+channelNames <- function(x){
+  if(!is.null(md <- attr(x, "meta")))
+    return(md[[2]]$channels$label)
+  return(NULL)
+}
+  
+# check if channels are indices or channel names
+#' @noRd
+.ch_ind <- \(x, ch){
+  if(is.null(ch))
+    return(1)
+  lbs <- channelNames(x)
+  if(all(ch %in% lbs)){
+    return(match(ch,lbs))
+  } else if(!any(ch %in% lbs)){
+    warning("Some channels are not found, picking first one!")
+    return(1)
+  } else {
+    warning("Channels are not found, picking first one!")
+    return(1)
+  }
+  return(NULL)
+}
+
 .guess_scale <- \(x, w, h) {
     n <- length(dim(x))
     i <- ifelse(n == 3, -1, TRUE)
@@ -47,14 +132,17 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 #' @importFrom methods as
 #' @importFrom grDevices rgb
 #' @importFrom DelayedArray realize
-.df_i <- \(x, k=NULL) {
-    a <- .get_plot_data(x, k)
-    a <- if (dim(a)[1] == 1) a[rep(1,3),,] else a
-    a <- realize(as(a, "DelayedArray"))
-    img <- rgb(
-        maxColorValue=max(a),
-        c(a[1,,]), c(a[2,,]), c(a[3,,]))
-    array(img, dim(a)[-1])
+.df_i <- \(x, k=NULL, ch=NULL, c=NULL) {
+  a <- .get_plot_data(x, k)
+  ch_i <- .ch_ind(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)
+  if(!.is.rgb(x))
+    a <- .manage_channels(a, ch_i, c)
+  apply(a, c(2, 3), \(.) do.call(rgb, as.list(.))) 
 }
 
 .get_wh <- \(x, i, j) {
@@ -75,13 +163,13 @@ plotSpatialData <- \() ggplot() + scale_y_reverse() + .theme
 
 #' @rdname plotImage
 #' @export
-setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=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)
-    wh <- .get_wh(x, i, j)
-    .gg_i(df, wh$w, wh$h)
+setMethod("plotImage", "SpatialData", \(x, i=1, j=1, k=NULL, ch=NULL, c=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)
+  wh <- .get_wh(x, i, j)
+  .gg_i(df, wh$w, wh$h)
 })
\ No newline at end of file
diff --git a/R/utils.R b/R/utils.R
index efdbfef..6eca93e 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -33,4 +33,15 @@
         plot.title=element_text(hjust=0.5),
         axis.text=element_text(color="grey"),
         axis.ticks=element_line(color="grey"))
-)
\ No newline at end of file
+)
+
+# default colors (from ImageJ/Fiji)
+.DEFAULT_COLORS <- c("red", "green", "blue", "gray", "cyan", "magenta", "yellow")
+
+# image data type factors (max values)
+# TODO: add more cases from other data types
+# https://doc.embedded-wizard.de/uint-type
+.DTYPE_MAX_VALUES <- list("uint8" = 255,
+                          "uint16" = 65535,
+                          "uint32" = 4294967295,
+                          "uint64" = 2^64 - 1)
\ No newline at end of file
diff --git a/man/plotImage.Rd b/man/plotImage.Rd
index 96f65be..daca083 100644
--- a/man/plotImage.Rd
+++ b/man/plotImage.Rd
@@ -8,7 +8,7 @@
 \usage{
 plotSpatialData()
 
-\S4method{plotImage}{SpatialData}(x, i = 1, j = 1, k = NULL)
+\S4method{plotImage}{SpatialData}(x, i = 1, j = 1, k = NULL, ch = NULL, c = NULL)
 }
 \arguments{
 \item{x}{\code{\link{SpatialData}} object.}
@@ -20,6 +20,10 @@ plotSpatialData()
 \item{k}{index of the scale of an image; by default (NULL), will auto-select 
 scale in order to minimize memory-usage and blurring for a target size of 
 800 x 800px; use Inf to plot the lowest resolution available.}
+
+\item{ch}{the image channels to be used for plotting (default: first channel)}
+
+\item{c}{plotting aesthetics; color}
 }
 \value{
 ggplot
diff --git a/tests/testthat/test-plotImage.R b/tests/testthat/test-plotImage.R
new file mode 100644
index 0000000..aa8e473
--- /dev/null
+++ b/tests/testthat/test-plotImage.R
@@ -0,0 +1,22 @@
+require(SpatialData, quietly=TRUE)
+x <- file.path("extdata", "blobs.zarr")
+x <- system.file(x, package="SpatialData")
+x <- readSpatialData(x, tables=FALSE)
+
+test_that("get/check channel names", {
+  
+  # get channel names
+  expect_equal(channelNames(image(x,1)), c(0,1,2))
+
+  # get indices of channels
+  expect_equal(.ch_ind(image(x,1), ch = c(2,0,1)), c(3,1,2))
+  expect_warning(expect_equal(.ch_ind(image(x,1), ch = 45), 1)) # return first if no matching channel
+  
+  # .is.rgb
+  expect_true(.is.rgb(image(x,1)))
+})
+
+# TODO: any tests for image array normalization ? 
+test_that(".normalize_image_array", {
+  skip()
+})
diff --git a/vignettes/SpatialData.plot.Rmd b/vignettes/SpatialData.plot.Rmd
index 4a770d1..bff6cce 100644
--- a/vignettes/SpatialData.plot.Rmd
+++ b/vignettes/SpatialData.plot.Rmd
@@ -238,6 +238,49 @@ wrap_plots(nrow=1, lapply(seq(3), \(.)
     plot_layout(guides="collect")
 ```
 
+## CyCIF (MCMICRO)
+
+Small lung adenocarcinoma, 250 MB; 1 image, 2 labels, 2 tables.
+
+```{r mcmicro-read}
+dir.create(td <- tempfile())
+pa <- unzip_spd_demo(
+    zipname="mcmicro_io.zip", 
+    dest=td, source="biocOSN")
+(x <- readSpatialData(pa, anndataR=FALSE))
+```
+
+Getting channel names for the image
+
+```{r mcmicro-channels}
+# TODO: add this to SpatialData as an ImageArray method
+# channelNames(image(x,1))
+```
+
+Plotting with multiple image channels.
+
+```{r mcmicro-plot}
+plotSpatialData() + plotImage(x, 1, ch = c("DNA_6", "ELANE", "CD57"), c = c("blue", "cyan", "yellow"))
+```
+
+## IMC (Steinbock)
+
+4 different cancers (SCCHN, BCC, NSCLC, CRC), 820 MB; 14 images, 14 labels, 1 table.
+
+```{r steinbock-read}
+dir.create(td <- tempfile())
+pa <- unzip_spd_demo(
+    zipname="steinbock_io.zip", 
+    dest=td, source="biocOSN")
+(x <- readSpatialData(pa, anndataR=FALSE))
+```
+
+Plotting with multiple image channels.
+
+```{r steinbock-plot}
+plotSpatialData() + plotImage(x, 1, ch = c(0,1,2), c = c("blue", "cyan", "yellow"))
+```
+
 # Masking
 
 Back to blobs...