---
title: "Calculating brightness with single images"
author: "Rory Nolan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Calculating brightness with single images}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

In this vignette, we will look at calculating the brightness of single images interactively in an R session.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = "#>", 
                      fig.width = 7, fig.height = 6)
set.seed(1)
pacman::p_load(ggplot2, dplyr, tidyr)
```

First let's load `nandb` and `ijtiff`, which is for reading TIFF files: 
```{r load}
pacman::p_load(nandb, ijtiff)
```


## Ordinary brightness

The package contains a sample image series which can be found at  
`system.file("extdata", "two_ch.tif", package = "nandb")`. It's 2 channels each with 100 frames. Diffusing fluorescent particles are imaged. Protein A is labelled with a red dye and red photons are collected in channel 1. Protein B is labelled in green and green photons are collected in channel 2. 

The image can be read into R with the `read_tif()` command provided by the `ijtiff` package. We'll assign it to a variable called `my_img`.

```{r read}
path <- system.file("extdata", "two_ch.tif", package = "nandb")
print(path)
my_img <- read_tif(path)
```

`my_img` is now a 4-dimensional array. Slots 1 and 2 hold the `y` and `x` pixel positions, slot 3 indexes the channel and slot 4 indexes the frame.

```{r dim}
dim(my_img)
```


### The need for detrending

Plotting the mean intensities of the frames in the two channels, we can see that the second channel has more obvious bleaching.

```{r plot-mean-profile, echo=FALSE, fig.height=3, message=FALSE}
plot_tbl <- tibble(
  frame = seq_len(dim(my_img)[4]),
  channel1 = apply(autothresholdr::mean_stack_thresh(my_img[, , 1, ], "h"), 4, 
                   mean, na.rm = TRUE),
  channel2 = apply(autothresholdr::mean_stack_thresh(my_img[, , 2, ], "h"), 4, 
                   mean, na.rm = TRUE)
  ) %>% 
  gather(ch, mean, -frame)
ch2range <- plot_tbl %>% 
  filter(ch == "channel2") %>% 
  pull(mean) %>% 
  range() %>% 
  round()
pl <- ggplot(filter(plot_tbl, ch == "channel1"), aes(frame, mean)) + 
  geom_smooth() + geom_point() + ggtitle("Channel 1") + ylim(0.75, 0.95)
pr <- ggplot(filter(plot_tbl, ch == "channel2"), aes(frame, mean)) + 
  geom_smooth() + geom_point() + ggtitle("Channel 2") + ylim(26, 29)
gridExtra::grid.arrange(pl, pr, nrow = 1)
```

The presence of bleaching in this image series tells us that we should employ a detrending routine as part of our calculations. Often it's not obvious whether or not detrending is required. `nandb` includes _Robin Hood_ detrending which won't interfere too much with images that don't need detrending, whilst effictively detrending images that do need it. Hence, it's pretty safe to have this detrending option turned on in general, and you should only leave it off if you're sure it's not needed. To turn on detrending, use `detrend = TRUE` in the `brightness()` functions.


### Thresholding

You can see here that this is an image of part of a cell, with the edge of the cell across the top left and hence the top left corner of the image is not cell, just background. It's important to _threshold_ away this background part: the detrending routine assumes that all parts of the image are part of the region of interest (the cell) and later, when you calculate summary statistics like mean/median brightness, you also want to background regions to be excluded. Hence, you need to set the background parts to `NA` before detrending and brightness calculations. `nandb` has all of the thresholding functionality of the _ImageJ Auto Threshold_ plugin. You can read more about this at https://imagej.net/Auto_Threshold. My favourite method is _Huang_. Let's look at both of these channels with _Huang_ thresholding.

```{r display-thresholded-mean, echo=FALSE, fig.height=4}
graphics::par(mfrow = c(1, 2))
display(
  detrendr::mean_pillars(
    autothresholdr::mean_stack_thresh(my_img[, , 1, ], "h")
  )
)
display(
  detrendr::mean_pillars(
    autothresholdr::mean_stack_thresh(my_img[, , 2, ], "h")
  )
)
```

That seems to have worked: the background region in the top left is now greyed out, indicating that it has been set to `NA`. _Huang_ thresholding can be slow, so if you want something similar but faster, try _Triangle_. Always check that your thresholding _looks right_ afterwards.

Note that if all of the image is of interest, detrending is not necessary and should not be done.


### Calculation including thresholding and detrending

To calculate the brightness  of this image with _Huang_ thresholding and _Robin Hood_ detrending, with the _epsilon_ definition of brightness, run

```{r brightness} 
my_brightness_img <- brightness(my_img, def = "e", 
                                thresh = "Huang", detrend = TRUE)
```


#### Timeseries calculation

To calculate a brightness timeseries with 50 frames per set, run

```{r brightness-ts} 
my_brightness_ts_img <- brightness_timeseries(my_img, def = "e", 
                                              frames_per_set = 50,
                                              thresh = "Huang", detrend = TRUE)
```

To calculate an overlapped brightness timeseries with 50 frames per set, run

```{r brightness-ts-overlapped} 
my_brightness_ts_img_overlapped <- brightness_timeseries(my_img, def = "e", 
                                                         frames_per_set = 50,
                                                         overlap = TRUE,
                                                         thresh = "Huang", 
                                                         detrend = TRUE)
```


### Saving Images

We can save the brightness images with `write_tif()`:

```{r write-tif, eval=FALSE}
write_tif(my_brightness_img, "desired/path/of/my-brightness-img")
write_tif(my_brightness_ts_img, "desired/path/of/my-brightness-ts-img")
write_tif(my_brightness_ts_img_overlapped, 
          "desired/path/of/my-brightness-ts-img-overlapped")
```

### Studying the Distribution of Brightnesses

We can take a look at the distribution of brightnesses in channel 1 of the brightness image:
```{r brightness density, message=FALSE, fig.width=7, fig.height=7}
db <- density(my_brightness_img[, , 1, ], na.rm = TRUE)[c("x", "y")] %>% 
  as_tibble()
ggplot(db, aes(x, y)) + geom_line() + 
  labs(x = "brightness", y = "frequency") +
  xlim(-2, 2) + ggtitle("Channel 1 brightness distribution")
```

We can compare that to the distribution of brightnesses from immobile particles:

```{r immobile plot, fig.width=7, fig.height=7}
immobile_brightnesses <- matrix(rpois(50 * 10^6, 50), nrow = 10^6) %>% 
  {matrixStats::rowVars(.) / rowMeans(.) - 1}
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(-2, 2) + 
  ggtitle("Channel 1 brightness distribution compared to immobile entities")
``` 

Note that to create plots like these, you'll need knowledge of the `ggplot2` package.
