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

In this vignette, we will look at detrending single images interactively in an R session.

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

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

## An image in need of detrending

The package contains a sample image series which can be found at  
`system.file("extdata", "2ch100frame.tif", package = "detrendr")`. 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", "2ch100frame.tif", package = "detrendr")
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)
```

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(method = "loess", formula = "y ~ x") +
  geom_point() +
  ggtitle("Channel 1") +
  ylim(0.75, 0.95)
pr <- ggplot(filter(plot_tbl, ch == "channel2"), aes(frame, mean)) +
  geom_smooth(method = "loess", formula = "y ~ x") +
  geom_point() +
  ggtitle("Channel 2") +
  ylim(26, 29)
gridExtra::grid.arrange(pl, pr, nrow = 1)
```

However, notice that channel 2 bleaches from `r ch2range[2]` to `r ch2range[1]` which is approximately `r round(abs(diff(ch2range)) / ch2range[2], 2) * 100`%. That's not much, so little correction will be necessary. This is good, you should always endeavour to keep bleaching below 20% and bleaching correction (detrending) should be viewed as a necessary evil.

The mean intensity of channel 1 is `r round(mean(my_img[, , 1, ]), 2)` and the mean intensity of channel 2 is `r round(mean(my_img[, , 2, ]), 2)`. Let's view the mean intensity images of channels 1 and 2. 

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


## 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), so we need to set the background parts to `NA` beforehand to tell the detrending routine that this area should be excluded. `detrendr` has all of the thresholding functionality of the _ImageJ Auto Threshold_ plugin. You can read more about this at https://imagej.net/plugins/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(mean_pillars(autothresholdr::mean_stack_thresh(my_img[, , 1, ], "h")))
display(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.

Thresholding of images is done with the `autothresholdr` package. Thresholding of a stack is preformed with `mean_stack_thresh()` which thresholds away pixels in the stack whose mean intensity is less than a certain threshold value.

```{r autothresholdr-code} 
my_img_threshed <- autothresholdr::mean_stack_thresh(my_img, "Huang")
```

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

## Detrending

The best detrending method is _Robin Hood_, so I'll be using that here. Other detrending methods provided in this package are there for legacy reasons.

To detrend the above thresholded with the _Robin Hood_ algorithm, run

```{r detrend} 
my_detrended_img <- img_detrend_robinhood(my_img_threshed)
```

Let's check out the mean intensity profiles of this detrended image.

```{r plot-detrended-mean-profile, echo=FALSE, fig.height=3, message=FALSE}
plot_tbl <- tibble(
  frame = seq_len(dim(my_detrended_img)[4]),
  channel1 = apply(my_detrended_img[, , 1, ], 3, mean, na.rm = TRUE),
  channel2 = apply(my_detrended_img[, , 2, ], 3, mean, na.rm = TRUE)
) %>%
  gather(ch, mean, -frame)
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)
```

These mean intensity profiles neatly show that the image series has been detrended. Beware, however, that these mean intensity profiles are not a great way to compare detrending routines. There is the temptation to say whichever detrending routine finishes with the flattest mean intensity profile is the best, but this is not true because it is possible to over-detrend (which will over-flatten the mean intensity profile). More rigorous procedures are needed to compare detrending routines, for example using images with simulated photobleaching or other trend-inducing phenomena introduced.