R/analysis.R

In jefferis/scanimage: read and analyse 2 photon images from ScanImage

#' Construct a mean image, blur and threshold to make an ROI for analysis
#'
#' @details This is currently targeted at single plane imaging data.
#'
#' @param x A hyperstack or path to one or more TIFF files (see
#'   \code{\link{make_hyperstack}})
#' @param blur The gaussian sigma to use
#' @param thr A thresholding specification (passed to
#'   \code{imager::\link[imager]{threshold}})
#' @param ... Additional arguments passed to \code{\link{make_hyperstack}}
#'
#' @return A 2D mask in the form of an imager::cimg object
#' @seealso \code{\link[imager]{cimg}}, \code{\link{t_profile}},
#'   \code{\link{make_hyperstack}}
#' @export
#'
#' @importFrom imager threshold as.cimg isoblur
make_roi_from_mean <- function(x, blur=5, thr="95%", ...) {
  if(is.character(x)) {
    x <- if(length(x)==1) make_hyperstack(x, ...) else make_hyperstack_multifile(x, ...)
  }
  # make mean image
  mx=apply(x, 1:2, mean)

  cx=as.cimg(mx)
  bcx=isoblur(cx, sigma=blur)
  imager::threshold(bcx, thr)
}

#' Construct a temporal profile of voxels within an ROI
#'
#' @details When the \code{baseline} frames are specified the result will be
#'   normalised to have a mean value of 1.0 in this period.
#'
#'   When \code{freq=TRUE} or is passed a numeric value, the return type will be
#'   a time series \code{\link[stats]{ts}} object. This can then be used in
#'   plots or subjected to filtering etc.
#'
#' @param x A hyperstack with imaging data
#' @param mask An ROI mask e.g. from \code{\link{make_roi_from_mean}}
#' @param baseline Optionally specifies which samples are the baseline (see
#'   details).
#' @param freq Logical specifying whether to use frequency information in
#'   attributes of \code{x} (default \code{TRUE}) or numeric value specifying
#'   frequency directly. See details.
#'
#' @return A matrix with the third and fourth dimensions of \code{x}, optionally
#' @export
#' @importFrom stats ts
#'
#' @seealso \code{\link{make_roi_from_mean}}, \code{\link[stats]{ts}},
#'   \code{\link[stats]{plot.ts}}
t_profile <- function(x, mask=NULL, baseline=NULL, freq=TRUE) {
  if(!is.null(mask)) {
    x2=apply(x, 3:4, "[", mask)
  } else {
    x2=apply(x, 3:4, c)
  }
  cmx=colMeans(x2)
  if(!is.null(baseline)) {
    cmx=scale(cmx, center = F, scale = colMeans(cmx[baseline,,drop=F]))
  }

  if(isTRUE(freq)) freq=scanimageinfo(x)$freq
  else if(isTRUE(!freq)) freq=NULL

  if(!is.null(freq)){
    ts(cmx, frequency = freq, start = 0)
  } else {
    cmx
  }
}

#' Make a tidy data frame from time profile output e.g. for use with ggplot
#'
#' @details The idea is that you can take results generated by
#'   \code{\link{t_profile}} (usually as a time series object) and turn them
#'   into a \bold{tidy} data.frame suitable for \code{tidyverse} style analysis.
#'
#'   If you have a \code{list} of t_profile results generated for multiple
#'   separate experiments, you can gather those together as well to make a
#'   single grand summary data.frame.
#'
#' @param x A time profile as generated by \code{\link{t_profile}}.
#' @param ... Additional argument passed to methods and eventually to
#'   \code{\link[tidyr]{gather}}.
#' @export
#' @seealso \code{\link[tidyr]{gather}}, \code{\link{t_profile}}.
gather_t_profile <- function(x, ...) UseMethod("gather_t_profile")

#' @export
#' @importFrom stats time
gather_t_profile.ts <- function(x, ...) {
  df=as.data.frame(x)
  names(df)=make.unique(names(df))
  res=cbind(df, T=c(time(x[,1])))
  gather(res, "odour", "F", -T, ...)
}

#' @export
#' @importFrom tidyr gather
gather_t_profile.default <- function(x, ...) {
  df=as.data.frame(x)
  names(df)=make.unique(names(df))
  res=cbind(df, i=1:nrow(df))
  gather(res, "odour", "F", -i, ...)
}

#' @export
gather_t_profile.list <- function(x, ...) {
  lgp=lapply(x, gather_t_profile, ...)
  resdf=do.call(rbind, lgp)
  if(is.null(names(x))) names(x)=seq_along(x)
  resdf$expt=rep(names(x), sapply(lgp, nrow))
  resdf
}

#' @export
gather_t_profile.by <- gather_t_profile.list

#' @importFrom graphics rect par
#' @importFrom grDevices rgb
stimrect <- function(x, col=rgb(1,0,0,0.3), border=NA, ...) {
  rect(x[1], par("usr")[3], x[2], par("usr")[4], col = col, border=border, ...)
}