Number and Brightness Image Analysis

Rory Nolan

2017-07-11

Loading the libraries

required <- c("nandb", "EBImage", "dplyr", "ggplot2", "matrixStats")
sapply(required, library, character.only = TRUE)
#> Warning: package 'dplyr' was built under R version 3.4.1

Reading and Displaying Images

Now let’s read in an image:

img <- system.file("extdata", "50.tif", package = "nandb") %>% ReadImageData
dim(img)  # get img's dimensions
#> [1] 50 50 50

So we can see that img is a 3-dimensional array of 70 slices each of which is a \(50 \times 50\) image. We can view the first slice:

display(normalize(img[, , 1]), method = "raster")

Preprocessing, Brightness Calculation and Post-Filtering

Now we can calculate the brightness image based on this image series, with an exponential filtering detrend with a time constant of (say) 10 frames and using the Huang thresholding method via

img_brightness <- Brightness(img, tau = "auto", mst = "Huang")
MatrixRasterPlot(img_brightness, log.trans = TRUE)

and we can see what happens when we median filter.

Brightness(img, tau = "auto", mst = "Huang", filt = "median") %>% 
  MatrixRasterPlot(scale.name = "brightness", log.trans = TRUE)

Saving Text Images

We can save the brightness image (let’s do it without median filtering this time). It saves as a text csv file (not as a tiff image, since the brightness values can be any non-negative real number, but tiff images can only store pixel values as integers).

WriteImageTxt(img_brightness, "BrightnessImage")
list.files(pattern = "\\.csv$")  # check that the csv file was written to the folder
#> [1] "BrightnessImage.csv"

Studying the Distribution of Brightnesses

We can take a look at the distribution of brightnesses in that image:

db <- density(img_brightness, na.rm = TRUE) %>% 
  {data.frame(x = .$x, y = .$y)}
ggplot(db, aes(x, y)) + geom_line() + 
  labs(x = "brightness", y = "frequency") +
  xlim(0.5, 2)
#> Warning: Removed 268 rows containing missing values (geom_path).

and we can compare it to the distribution of brightnesses from immobile particles:

immobile_brightnesses <- matrix(rpois(50 * 10^6, 50), nrow = 10^6) %>% 
  {rowVars(.) / rowMeans(.)}
di <- density(immobile_brightnesses) %>% {data.frame(x = .$x, y = .$y)}
rbind(mutate(db, mobility = "mobile"), mutate(di, mobility = "immobile")) %>%      mutate(mobility = factor(mobility)) %>% 
  ggplot(aes(x, y)) + geom_line(aes(colour = mobility)) + 
  labs(x = "brightness", y = "frequency") +
  xlim(0.5, 2)
#> Warning: Removed 406 rows containing missing values (geom_path).

Calculating the Brightness of a Monomer

Say we know that our monomeric brightness is 1.02. Then we can create a plot of all kmers present:

KmerPlot(img_brightness, 1.04)