R/plotting-functions.R

In CeMOS-Mannheim/moleculaR: Spatial Probabilistic Mapping of Metabolite Ensembles in Mass Spectrometry Imaging

#' @include manualSpatstatImport.R
NULL


#' Plot a spectrum focusing on a specific mass-range
#'
#' This function plots a spectrum (`MassPeaks` or `MassSpectrum`) and provides
#' a focus on a mass range around a given m/z of interest taking into
#' consideration the fwhm at that mass.
#'
#' @param x an object of type `MassPeaks` or `MassSpectrum`. Either a single
#' spectrum or a mean spectrum (see `?MALDIquant::mergeMassPeaks`).
#' @param m a numeric (vector), the m/z value(S) at which to do the plotting.
#' @param fwhm a numeric (vector), the fwhm value(s) at `m`, normally generated by `getFwhm`.
#' Must have the same length as `m`.
#' @param sigmalim an integer, how many sigmas (of the calculated Gaussian)
#' around `m` should be considered in the plot. Note that the superimposed
#' Gaussian always spans `m ± sigma`.
#' @param col a character (vector), colors used for the m/z values in `m`.
#'
#' @return
#' a plot.
#'
#' @export
#'

focusPlot <- function(x, m, fwhm, sigmalim = 6, col = NULL){

         m <- sort(m)

         if(length(m) != length(fwhm)){
                  stop("length(m) is not equal to length(fwhm)!\n")
         }


         # find sigma intervals
         sigma <- fwhm/(2*sqrt(2*log(2)))
         s3 <- sigma * 3
         sxlim <- sigma * sigmalim

         # scale the y-axis intensities according to the focus region
         idx <- which(x@mass > min(m-sxlim) & x@mass < max(m+sxlim))

         if(length(idx) > 1){
                  yrange <- range(x@intensity[idx])
         }

         if(length(idx) == 1){
                  yrange <- c(0, x@intensity[idx])
         }

         if(length(idx) == 0){
                  warning("No peaks detected in the given mass range, consider increasing ",
                          "'sigmalim'.\n")
                  yrange <- c(0, median(x@intensity))
         }

         xrange <- c(min(m-sxlim), max(m+sxlim))

         par(mar= c(5,5,5,7))

         MALDIquant::plot(x, xlim = xrange, ylim = yrange,
                          frame = FALSE, axes = FALSE, ylab = "",
                          xlab = "",
                          main = paste0("m/z " ,round(m, 4),
                                        " ± ", round(s3, 4)),
                          lwd = 2, cex.lab = 1.25, cex.axis = 1.25)

         axis(1, at = pretty(xrange), col = "black",
              lwd = 2, cex = 1.25, col.axis = "black",
              cex.axis = 1.25, lty = "solid")

         mtext(substitute(paste(italic("m/z"))), side=1, line=2, cex=1.25, col="black")

         axis(2, at = pretty(yrange),
              col = "black", labels = format(pretty(yrange), scientific = TRUE),
              lwd = 2, cex = 1.25, col.axis = "black",
              cex.axis = 1.25, line = 1.5, lty = "solid")

         mtext("Intensity", side=2, line=4, cex=1.25, col="black")

         if(is.null(col)){
                  col <- rainbow(length(m))
         } else {
                  if(length(col) != length(m)){
                           warning("'col' has to have the same length as 'm'. The input has been ignored.\n")
                           col <- rainbow(length(m))
                  }
         }

         for(ii in 1:length(m)){
                  axis(1, at=m[ii], tick = T, labels="", col=col[ii],
                       lwd = 1.5, cex = 1.25, col.axis = col[ii],
                       gap.axis = 0, cex.axis = 1.25)





                  # gaussian visualization
                  gxx <- seq(m[ii]-s3[ii], m[ii]+s3[ii], 0.0001)
                  gyy <- .gf(x = gxx, mu = m[ii], sigma = sigma[ii])

                  gyy <- .linMap(gyy, yrange[1], yrange[2])

                  polygon(x = c(gxx[1], gxx, gxx[length(gxx)], gxx[1]),
                          y = c(0, gyy, 0, 0),
                          col =  to.transparent(col[ii], 0.25),
                          border = to.transparent(col[ii], 0.25))

                  abline(v = m[ii], lty="dashed", lwd = 0.75, col = col[ii])



         }

         par(new=TRUE)
         plot(1, axes = FALSE, type="n", xlab="", ylab="", xlim=xrange, ylim=c(0, 1), frame = FALSE)
         # axis(4, at=pretty(c(0,1)), col = "coral2",
         #      lwd = 2, cex = 1.25, col.axis = "coral2",
         #      cex.axis = 1.25, line = 1.5)
         axis(4, at=pretty(c(0,1)), col = "black",
              lwd = 2, cex = 1.25, col.axis = "black",
              cex.axis = 1.25, line = 1.5, lty = "dashed")
         mtext("Gaussian Weights", side=4, line=4, cex=1.25, col="black")

         # legend("topright",
         #        legend = c(paste0("fwhm = ", round(fwhm, 4)),
         #                   paste0("sigma = ", round(sigma, 4))),
         #        bty = "n",
         #        col = col)


}

# Gaussian function
.gf <- function(x, mu, sigma){
         1/(sqrt(2 * pi) * sigma) * exp(-((x - mu)^2/(2 * sigma^2)))
}

# linear mapping of values into a custom range
.linMap <- function(i, a, b) approxfun(range(i), c(a, b))(i)



#' Plot analytePointPattern
#'
#' A method to plot `anlaytePointPattern` objects. This calls `plot.ppp` with
#' pre-defined defaults. For more control over the plotting please use `plot.ppp`.
#'
#' @param obj S3 object of type `anlaytePointPattern` (and `spp`).
#' @param colourPal the colourmap to be used, see `?viridis::viridis_pal`.
#' @param uniformCol a character specifying a single colour. This will override `colourPal`.
#' @param transpFactor Transparency fraction. Numerical value or vector of values between 0 and 1,
#' giving the opaqueness of a colour. A fully opaque colour has `transpFactor=1`.
#' @param rescale logical, whether to scale the intensities of the output plot
#' to [0,1] interval.
#' @param pch     a positive integer, or a single character. See `?par`.
#' @param size    the size of the symbol: a positive number or zero. see`?symbolmap`.
#' @param main  character, title of the plot. If not given, the m/z value of `obj` is used.
#' @param leg.side position of legend relative to main plot.
#' @param leg.args list of additional arguments passed to control the legend. see
#' `?spatstat.geom::plot.ppp` for more details.
#' @param ... further arguments passed to `spatstat.geom::plot.ppp`.
#'
#' @return nothing, plots only.
#'
#' @export
#'
plotAnalyte <- function(obj, colourPal = "inferno", uniformCol = NULL, transpFactor = 1,
                        rescale = TRUE, pch = 19, size = 0.4, main = NULL,
                        leg.side = "right", leg.args = list(cex = 3, cex.axis = 1.25), ...){

      if(!("analytePointPattern" %in% class(obj))){
            stop("provided obj is not of type 'analytePointPattern'. \n")
      }

      if(obj$n > 30000) {
            cat("plotting ", format(obj$n, big.mark = ","), " points - this takes time! \n")
      }

      if(is.null(main)){
            if(length(obj$metaData$mzVals) > 1){
                  main <- paste0("Collective SPP of ", length(obj$metaData$mzVals)," m/z Values")
            } else{
                  mz <- ifelse(is.numeric(obj$metaData$mzVals),
                               as.character(round(obj$metaData$mzVals, 4)),
                               obj$metaData$mzVals) # to account for simulation cases

                  main <- paste0("m/z ", mz)
            }

      }

      # msi uses reversed y axis
      args <- list(...)

      if("ylim" %in% names(args)){
            yrange <- rev(args[["ylim"]])
      } else {
            yrange <- rev(obj$window$yrange)
      }

      if(obj$n == 0){

            plot.ppp(obj, use.marks = TRUE, which.marks = "intensity",
                     ylim = yrange,
                     main = main)

      } else{

            if(is.null(uniformCol)){

                  if(rescale){

                        obj <- .rescale(obj)
                  }


                  colfun <- colourmap(col = to.transparent((viridis::viridis_pal(option = colourPal)(100)), transpFactor),
                                      range = range(obj$marks$intensity))

                  plot.ppp(obj, use.marks = TRUE, which.marks = "intensity",
                           ylim = yrange,
                           main = main,
                           symap = symbolmap(pch = pch,
                                             cols = colfun,
                                             size = size,
                                             range = range(obj$marks$intensity)),
                           leg.side = leg.side,
                           leg.args = leg.args,
                           ...) # colors according to intensity

            } else{

                  if(rescale){

                        obj <- .rescale(obj)
                  }

                  col <- to.transparent(uniformCol, transpFactor)

                  plot.ppp(obj, use.marks = TRUE, which.marks = "intensity",
                           ylim = yrange,
                           main = main,
                           symap = symbolmap(pch = pch,
                                             cols = col,
                                             size = size),
                           leg.side = leg.side,
                           leg.args = leg.args,
                           ...)

            }


      }

}


#' Plot raster images of type `im`
#'
#' A method to plot `im` objects of the `spatstat` package. This calls `plot.im` with
#' pre-defined defaults. For more control over the plotting please use `plot.im`.
#'
#' @param obj S3 object of type `im`.
#' @param colourPal the colourmap to be used, see `?viridis::viridis_pal`.
#' @param uniformCol a character specifying a single colour. This will override `colourPal`.
#' @param rescale logical, whether to scale the intensities of the output plot
#' to (0,1] interval.
#' @param smooth logical, whether to apply image smoothing via an isotropic
#' Gaussian kernel.
#' @param transpFactor Transparency fraction. Numerical value or vector of values between 0 and 1,
#' giving the opaqueness of a colour. A fully opaque colour has `transpFactor=1`.
#' @param irange a numeric of length 2, a custome intensity range for plotting. Incompatible with `rescale`;
#' if provided, `rescale` will be set to `FALSE`.
#' @param ribargs a list of additional arguments to control the display of the ribbon. see
#' `?spatstat.geom::plot.im` for details.
#' @param ribsep a numeric controlling the space between the ribbon and the image.
#' @param ... further arguments passed to `plot.im`.
#'
#' @return nothing, plots only.
#'
#' @export
#'
#'
plotImg <- function(obj, colourPal = "inferno", uniformCol = NULL,
                    rescale = TRUE, smooth = FALSE, transpFactor = 1,
                    irange = NULL, ribargs = list(las = 2, cex.axis =1.25),
                    ribsep = 0.05, ...){

      if(class(obj) != "im"){
            stop("provided obj must be of type 'im'.\n")
      }

      # msi uses reversed y axis
      imageArgs <- list(...)

      if("ylim" %in% names(imageArgs)){
            yrange <- rev(imageArgs[["ylim"]])
      } else {
            yrange <- rev(obj$yrange)
      }


      if(is.null(uniformCol)){

            if(smooth){
                  obj <- Smooth.im(obj, sigma = 0.5,
                                   bleed = FALSE, normalise = TRUE)
            }

            if(rescale & is.null(irange)){

                  obj <- .rescale(obj)
            }

            if(is.null(irange)){
                  irange <- range(obj$v, na.rm = T)
            }



            plot.im(x = obj,
                    col = colourmap(to.transparent((viridis::viridis_pal(option = colourPal)(100)), transpFactor),
                                    range = irange),
                    ylim = yrange,
                    box = FALSE,
                    ribargs = ribargs,
                    ribsep = ribsep,
                    ...)


      } else{



            col <- to.transparent(uniformCol, transpFactor)


            plot.im(x = obj,
                    col = col,
                    ylim = yrange,
                    box = FALSE,
                    ribargs = ribargs,
                    ribsep = ribsep,
                    ...)

      }



}

#' Plot molecular Probability Maps
#'
#' A method to plot MPMs with different combinations.
#'
#' @param obj S3 object of type `molProbMap`.
#' @param what What to plot, c("detailed", "analytePointPattern", "csrPointPattern", "analyteDensityImage",
#' "csrDensityImage", "kdeIntensitiesDistr", "MPM"). The "detailed" option plots all options.
#' @param densityAsBase a logical, if set to `TRUE`, the base image, ontop of which
#' hotspot/coldspot contours are plotted, is taken from the estimated density image
#' of m/z MOI. Otherwise (default), the Gaussian-weighted ion image is used.
#' @param sppArgs a named list of arguments to be passed to `plotAnalyte` or plotting
#' "analytePointPattern" and "csrPointPattern". There are sensible defaults.
#' @param imageArgs a named list of arguments to be passed to `spatstat.geom::plot.im` for plotting
#' "analyteDensityImage", "csrDensityImage" and "MPM". There are sensible defaults.
#' @param fArgs a named list of arguments to be passed to `base::plot` for plotting
#' "kdeIntensitiesDistr". There are sensible defaults.
#' @param rescale logical, whether to scale the intensities of the output plot
#' to [0,1] interval.
#' @param smooth logical, whether to apply image smoothing via an isotropic
#' Gaussian kernel.
#' @param zero.rm logical, whether to remove zero-valued pixels from the intensity
#' image.
#' @param signifArea a character indicating which significance area to plot i.e. c("both", "hotspot", "coldspot").
#' @param ionImage an optional rastered image of type `im` of the corresponding "regular"
#' ion image, used for comparison. Could be generated via `moleculaR::searchAnalyte(..., wMethod = "sum")`
#' and subsequently using `pixellate`.
#' @param figLegend logical, whether to show the hotspot and coldspot legends.
#' @param lwd.signifArea numeric, the width of the hotspot and coldspot contours.
#' @param col.hotspot a character specifying the color of the hotspot contour.
#' @param col.coldspot a character specifying the color of the coldspot contour.
#'
#' @return nothing, plots only.
#'
#' @method plot molProbMap
#'
#' @export
#'
plot.molProbMap <- function(obj,
                            what = "detailed",
                            densityAsBase = FALSE,
                            sppArgs = list(),
                            imageArgs = list(),
                            fArgs = list(),
                            rescale = TRUE,
                            smooth = FALSE,
                            zero.rm = FALSE,
                            signifArea = "both",
                            ionImage = NA,
                            figLegend = TRUE,
                            lwd.signifArea = 5,
                            col.hotspot = "red",
                            col.coldspot = "blue"){

      # checks
      if(!(signifArea %in% c("hotspot", "coldspot", "both"))){
            stop("'signifArea' must be one of c('hotspot', 'coldspot', 'both').\n")
      }



      switch (what,
              "analytePointPattern" = {

                    if(rescale){
                          obj$sppMoi <- .rescale(obj$sppMoi)
                    }

                    # workout the plotting arguments
                    # sppDefaults    <- list(colourPal = "inferno", uniformCol = NULL, transpFactor = 1,
                    #                        pch = 19, size = 0.4, main = NULL)
                    # sppArgs        <- .mergeArgs(sppArgs, sppDefaults)

                    # call the plotting function
                    do.call(plotAnalyte, c(list(obj = obj$sppMoi), sppArgs))


              },
              "csrPointPattern" = {

                    if(rescale){
                          obj$csrMoi <- .rescale(obj$csrMoi)
                    }

                    # workout the plotting arguments
                    # sppDefaults    <- list(colourPal = "inferno", uniformCol = NULL, transpFactor = 1,
                    #                        pch = 19, size = 0.4, main = NULL)
                    # sppArgs        <- .mergeArgs(sppArgs, sppDefaults)



                    # call the plotting function
                    do.call(plotAnalyte, c(list(obj = obj$csrMoi), sppArgs))


              },
              "analyteDensityImage" = {

                    if(rescale){
                          obj$rhoMoi <- .rescale(obj$rhoMoi)
                    }

                    #workout the plotting arguments
                    imageDefaults  <- list(main = expression(paste(rho["MOI"], "(x,y)")))

                    imageArgs      <- .mergeArgs(imageArgs, imageDefaults)


                    # call the plotting function
                    do.call(plotImg, c(list(obj = obj$rhoMoi), imageArgs))


              },
              "csrDensityImage" = {

                    if(rescale){
                          obj$rhoCsr <- .rescale(obj$rhoCsr)
                    }

                    # workout the plotting arguments
                    imageDefaults  <- list(main = expression(paste(rho["CSR"], "(x,y)")))

                    imageArgs      <- .mergeArgs(imageArgs, imageDefaults)




                    # call the plotting function
                    do.call(plotImg, c(list(obj = obj$rhoCsr), imageArgs))


              },
              "kdeIntensitiesDistr" = {

                    if(rescale){
                          obj$rhoMoi <- .rescale(obj$rhoMoi)
                          obj$rhoCsr <- .rescale(obj$rhoCsr)
                    }

                    fmoi <- density(c(obj$rhoMoi$v), na.rm = TRUE)
                    fcsr <- density(c(obj$rhoCsr$v), na.rm = TRUE)
                    muCsr <- mean(obj$rhoCsr$v, na.rm = TRUE)
                    sigmaCsr <- sd(obj$rhoCsr$v, na.rm = TRUE)
                    lowerCutoff <- qnorm(0.05, muCsr, sigmaCsr)
                    upperCutoff <- qnorm(0.05, muCsr, sigmaCsr, lower.tail = FALSE)

                    ymax <- max(max(fmoi$y), max(fcsr$y))
                    ymin <- min(min(fmoi$y), min(fcsr$y))
                    xmax <- max(max(fmoi$x), max(fcsr$x))
                    xmin <- min(min(fmoi$x), min(fcsr$x))

                    # Create two panels side by side
                    #layout(t(1:1), widths=c(5,1))

                    # Set margins and turn all axis labels horizontally (with `las=1`)
                    #par(mar=rep(3, 4), oma=rep(4, 4), las=1)

                    # workout the plotting arguments
                    fDefaults  <- list(main = bquote(f[CSR]*(k) ~ and ~ f[MOI]*(k) ~ at ~ bw==.(obj$bw)),
                                       col = "black",
                                       ylim = c(ymin, ymax),
                                       xlim = c(xmin, xmax),
                                       xlab = "Intensities",
                                       ylab = "Density",
                                       bty = "n",
                                       lwd = 3)

                    fArgs      <- .mergeArgs(fArgs, fDefaults)

                    # call the plotting function
                    do.call(plot, c(list(fcsr), fArgs))


                    lines(fmoi, col =  "forestgreen", lwd = fArgs$lwd)

                    polygon(x = c(upperCutoff, xmax, xmax, upperCutoff),
                            y = c(ymin, ymin, ymax, ymax), col = to.transparent("coral2", 0.5),
                            border = NA)

                    polygon(x = c(xmin, lowerCutoff, lowerCutoff, xmin),
                            y = c(ymin, ymin, ymax, ymax), col = to.transparent("cornflowerblue", 0.5),
                            border = NA)

                    polygon(x = c(fcsr$x[fcsr$x >= upperCutoff],  upperCutoff),
                            y = c(fcsr$y[fcsr$x >= upperCutoff],  0),
                            col = to.transparent("red", 0.5),
                            border = NA)

                    polygon(x = c(fcsr$x[fcsr$x <= lowerCutoff], lowerCutoff),
                            y = c(fcsr$y[fcsr$x <= lowerCutoff],  0),
                            col = to.transparent("blue", 0.5),
                            border = NA)

                    if(figLegend){
                          legend("topright", legend = c(expression(paste(italic(f["CSR"]), italic("(k)"))),
                                                        expression(paste(italic(f["MOI"]), italic("(k)")))),
                                 lty = c("solid"), lwd = 2,
                                 col = c("black", "forestgreen"),
                                 bty = "n", horiz = FALSE)
                    }


              },
              "MPM" = {

                    #spwin <- obj$sppMoi$window

                    if(densityAsBase){

                          imgMpm <- obj$rhoMoi

                    } else {

                          if(smooth){
                                imgMpm <- spp2im(obj$sppMoi, zero.rm = zero.rm)
                                imgMpm <- Smooth.im(imgMpm, sigma = 0.5,
                                                    bleed = FALSE, normalise = TRUE)
                          }else{
                                imgMpm <- spp2im(obj$sppMoi, zero.rm = zero.rm)
                          }
                    }

                    if(rescale){
                          imgMpm <- .rescale(imgMpm)
                    }

                    # workout the plotting arguments
                    imageDefaults  <- list(main = NULL)

                    imageArgs      <- .mergeArgs(imageArgs, imageDefaults)

                    # fix title
                    if(is.null(imageArgs$main)){
                          if(length(obj$sppMoi$metaData$mzVals) > 1){
                                imageArgs$main <- paste0("CPPM of ", length(obj$sppMoi$metaData$mzVals)," m/z Values")
                          } else{
                                mz <- ifelse(is.numeric(obj$sppMoi$metaData$mzVals),
                                             as.character(round(obj$sppMoi$metaData$mzVals, 4)),
                                             obj$sppMoi$metaData$mzVals) # to account for simulation cases

                                imageArgs$main <- paste0("m/z ", mz, " MPM")
                          }
                    }



                    # call the plotting function
                    do.call(plotImg, c(list(obj = imgMpm), imageArgs))


                    if(signifArea == "both" | signifArea == "hotspot"){
                          plot.owin(obj$hotspotpp$window, col = rgb(1,1,1,0.0), border = "white", lwd = lwd.signifArea,  add = TRUE)
                          plot.owin(obj$hotspotpp$window, col = rgb(1,1,1,0.0), border = col.hotspot, lwd = lwd.signifArea/2, lty = "dashed",add = TRUE)
                    }

                    if(signifArea == "both" | signifArea == "coldspot"){
                          plot.owin(obj$coldspotpp$window, col = rgb(1,1,1,0.0), border = "white", lwd = lwd.signifArea,  add = TRUE)
                          plot.owin(obj$coldspotpp$window, col = rgb(1,1,1,0.0), border = col.coldspot, lwd = lwd.signifArea/2, lty = "dashed",add = TRUE)
                    }

                    if(figLegend){
                          legend("bottom", legend = c("Analyte Hotspot", "Analyte Coldspot"), lty = c("dashed"),
                                 col = c("red", "blue"), bty = "n", horiz = TRUE, inset = c(-0.5,0))

                    }


              },
              "detailed" = {

                    if(obj$sppMoi$n > 30000) {
                          cat("plotting ", format(obj$sppMoi$n, big.mark = ","), " points - this takes time! \n")
                    }

                    par(mfrow = c(3, 2))

                    plot(obj = obj, what = "csrPointPattern",
                         densityAsBase = densityAsBase,
                         sppArgs = sppArgs,
                         imageArgs = imageArgs,
                         fArgs = fArgs,
                         rescale = rescale,
                         signifArea = signifArea,
                         ionImage = ionImage,
                         figLegend = figLegend,
                         lwd.signifArea = lwd.signifArea,
                         col.hotspot = col.hotspot,
                         col.coldspot = col.coldspot)

                    plot(obj = obj, what = "analytePointPattern",
                         densityAsBase = densityAsBase,
                         sppArgs = sppArgs,
                         imageArgs = imageArgs,
                         fArgs = fArgs,
                         rescale = rescale,
                         signifArea = signifArea,
                         ionImage = ionImage,
                         figLegend = figLegend,
                         lwd.signifArea = lwd.signifArea,
                         col.hotspot = col.hotspot,
                         col.coldspot = col.coldspot)

                    plot(obj = obj, what = "csrDensityImage",
                         densityAsBase = densityAsBase,
                         sppArgs = sppArgs,
                         imageArgs = imageArgs,
                         fArgs = fArgs,
                         rescale = rescale,
                         signifArea = signifArea,
                         ionImage = ionImage,
                         figLegend = figLegend,
                         lwd.signifArea = lwd.signifArea,
                         col.hotspot = col.hotspot,
                         col.coldspot = col.coldspot)

                    plot(obj = obj, what = "analyteDensityImage",
                         densityAsBase = densityAsBase,
                         sppArgs = sppArgs,
                         imageArgs = imageArgs,
                         fArgs = fArgs,
                         rescale = rescale,
                         signifArea = signifArea,
                         ionImage = ionImage,
                         figLegend = figLegend,
                         lwd.signifArea = lwd.signifArea,
                         col.hotspot = col.hotspot,
                         col.coldspot = col.coldspot)

                    if(identical(ionImage, NA)){
                          plot(obj = obj, what = "kdeIntensitiesDistr",
                               densityAsBase = densityAsBase,
                               sppArgs = sppArgs,
                               imageArgs = imageArgs,
                               fArgs = fArgs,
                               rescale = rescale,
                               signifArea = signifArea,
                               ionImage = ionImage,
                               figLegend = figLegend,
                               lwd.signifArea = lwd.signifArea,
                               col.hotspot = col.hotspot,
                               col.coldspot = col.coldspot)
                    }

                    plot(obj = obj, what = "MPM",
                         densityAsBase = densityAsBase,
                         sppArgs = sppArgs,
                         imageArgs = imageArgs,
                         fArgs = fArgs,
                         rescale = rescale,
                         zero.rm = zero.rm,
                         smooth = smooth,
                         signifArea = signifArea,
                         ionImage = ionImage,
                         figLegend = figLegend,
                         lwd.signifArea = lwd.signifArea,
                         col.hotspot = col.hotspot,
                         col.coldspot = col.coldspot)

                    if(!(identical(ionImage, NA))){

                          if(rescale){
                                ionImage <- .rescale(ionImage)
                          }

                          # workout the plotting arguments
                          imageDefaults  <- list(main = "Corresponding Ion Image")

                          imageArgs      <- .mergeArgs(imageArgs, imageDefaults)



                          # call the plotting function
                          do.call(plotImg, c(list(obj = ionImage), imageArgs))

                    }


              },
              stop("argument 'what' must be one of c('detailed', 'analytePointPattern', 'csrPointPattern', ",
                   "'analyteDensityImage','csrDensityImage', 'kdeIntensitiesDistr' or 'MPM').\n")
      )

}

# merge default arguments with provided ones
# providedArgs and defualt args are named lists
.mergeArgs <- function(providedArgs, defaultArgs){

      if(length(providedArgs) == 0){
            providedArgs <- defaultArgs
      } else {
            providedArgs <- c(providedArgs, defaultArgs[!(names(defaultArgs) %in% names(providedArgs))])
      }

      return(providedArgs)
}