Thresholding Image Stacks

if (utils::packageVersion("knitr") >= "1.20.15") {
  knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7, fig.height = 6,
    tidy = "styler"
  )
} else {
  knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7, fig.height = 6
  )
}
library(magrittr)
apply_on_pillars <- function(arr3d, FUN) {
  if (length(dim(arr3d)) == 4 && dim(arr3d)[3] == 1) arr3d %<>% {.[, , 1, ]}
  apply(arr3d, c(1, 2), FUN) %>% {
    if (length(dim(.)) == 3) {
      aperm(., c(2, 3, 1))
    } else {
      .
    }
  }
}

Stacks of Images

50.tif is a TIFF file which is a stack of 50 images of a bit of a cell taken over a short space of time (i.e. it's a video).

img <- ijtiff::read_tif(system.file("extdata", "50.tif", 
                                    package = "autothresholdr"))
dim(img)

Let's display the first 3 frames:

first3 <- matrix(max(img), ncol = 3 * dim(img)[2] + 2, nrow = dim(img)[1])
for (i in 1:3) {
  first3[, (i - 1) * dim(img)[2] + i + seq_len(dim(img)[2]) - 1] <- 
    img[, , 1, i]
}
ijtiff::display(first3)

and the last 3 frames:

last3 <- matrix(max(img), ncol = 3 * dim(img)[2] + 2, nrow = dim(img)[1])
for (i in 1:3) {
  last3[, (i - 1) * dim(img)[2] + i + seq_len(dim(img)[2]) - 1] <- 
    img[, , 1, dim(img)[4] - (i - 1)]
}
ijtiff::display(last3)

You'll notice that these images are almost identical. That's because they're images of the same area taken very quickly one after another. There are two ways to threshold an image stack like this. One is the naiive way: find a threshold based on every pixel in the stack, then pixels below the threshold are excluded (set to NA). Let's try that:

library(autothresholdr)
naiively_threshed_img <- apply_mask(img, "tri")
attr(naiively_threshed_img, "thresh")  # The threshold chosen by "Triangle" is 4

Now let's display the first 3 frames and the last 3 frames:

first3 <- matrix(max(naiively_threshed_img), ncol = 3 * dim(naiively_threshed_img)[2] + 2, nrow = dim(naiively_threshed_img)[1])
for (i in 1:3) {
  first3[, (i - 1) * dim(naiively_threshed_img)[2] + i + seq_len(dim(naiively_threshed_img)[2]) - 1] <- 
    naiively_threshed_img[, , 1, i]
}
ijtiff::display(first3)
last3 <- matrix(max(naiively_threshed_img), ncol = 3 * dim(naiively_threshed_img)[2] + 2, nrow = dim(naiively_threshed_img)[1])
for (i in 1:3) {
  last3[, (i - 1) * dim(naiively_threshed_img)[2] + i + seq_len(dim(naiively_threshed_img)[2]) - 1] <- 
    naiively_threshed_img[, , 1, dim(naiively_threshed_img)[4] - (i - 1)]
}
ijtiff::display(last3)

If you look closely, you can see that the threshold mask is (slightly) different for different frames. Let's highlight the pixels which are sometimes thresholded away, sometimes not:

naiively_threshed_img %>% 
  apply_on_pillars(function(x) {
    (sum(is.na(x)) > 0) && (sum(is.na(x)) < length(x))}) %>% 
  ijtiff::display()

So (unsurprisingly), it seems that around the edges of the cell (where the signal from the cell is more feint), the pixels are sometimes thresholded away, sometimes not. There are also some seemingly random pixels within the cell which are sometimes thresholded away, sometimes not.

Now, given that you know that the cell is more or less stationary and you want the threshold to get rid of the non-cell bits and keep the cell bits, its reasonable to assert that the mask should be the same for every frame. It's possible to apply the same mask to every frame, and to compute this mask, it makes sense to incorporate information from all of the frames. This is what mean_stack_thresh() and med_stack_thresh() do. mean_stack_thresh() computes the mask based on the mean of all of the frames (gotten by calculating the mean intensity of the stack at each pixel position). med_stack_thresh() uses the median instead of the mean. They're both very similar. If you don't know which one to use, just use mean_stack_thresh().

Let's visualize them both:

ijtiff::display(mean_stack_thresh(img, "tri"))
ijtiff::display(med_stack_thresh(img, "tri"))

You can see that the results of mean_stack_thresh() and med_stack_thresh() are similar but not identical. Both do a fine job.

Note that if the cell (or whatever is recorded over the course of several frames) is not stationary (or almost stationary), then mean_stack_thresh() and med_stack_thresh() are not appropriate.



Try the autothresholdr package in your browser

Any scripts or data that you put into this service are public.

autothresholdr documentation built on Feb. 16, 2023, 6:12 p.m.