%\VignetteIndexEntry{MaxContrastProjection}
%\VignetteKeywords{MaxContrastProjection}
%\VignettePackage{MaxContrastProjection}
%\VignetteEngine{knitr::knitr}

\documentclass[10pt,a4paper]{article}

<<style, echo=FALSE, results='asis'>>=
BiocStyle::latex()
@

\RequirePackage{amsfonts,amsmath,amstext,amssymb,amscd}
\usepackage{graphicx}
\usepackage{color,colortbl}
\usepackage{caption}
\usepackage{subcaption}
%\usepackage{courier}
%\graphicspath{ {../inst/extdata/} }
%\SweaveOpts{keep.source=TRUE,eps=FALSE}

\begin{document}
%\SweaveOpts{concordance=TRUE}

\title{MaxContrastProjection}

\author{Jan Sauer}

\maketitle

\tableofcontents

\section{Getting Started}

\Biocpkg{MaxContrastProjection} is an \R{} package to project 3D image stacks 
including common projections, such as the maximum and minimum intensity 
projections, as well as a novel maximum contrast projection to retain 
information about fine structures that would otherwise be lost in intensity-
based projections. This package is distributed as part of the \Bioconductor{} 
project and can be installed by starting \R{} and running:

<<eval=FALSE>>=
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("MaxContrastProjection")
@

\section{Introduction}

A common problem when recording 3D fluorescent microscopy images is how to 
properly present these results in 2D. Large structures with elements in 
multiple focal planes become very difficult to image properly without the use 
of some form of projection. Objects in the focal plane will be sharp and the 
emitted light focused, whereas objects outside of the focal plane will be 
blurred and the emitted light is spread out over a larger area of the 
detector. Consequently, a maximum intensity projection is commonly used on the 
image stack in order to find the focal plane of an object.

The problem with this approach is that the blurring makes it more difficult to 
tell apart foreground and background of the image. If the exact position of an 
object in a medium is unknown, or if it has a 3D structure, it must be imaged 
at various focal lengths. Fig. \ref{fig:dapizstacks} shows a segment of an 
organoid imaged at different focal lengths.

\begin{figure}[!htb]
\centering
\includegraphics[width=0.95\textwidth]{imgs/DAPI_montage.eps}
\caption{A cluster of cells in an organoid at different focal lengths}
\label{fig:dapizstacks}
\end{figure}

A look at the maximum intensity projection, shown in Fig. 
\ref{fig:dapimaxproj}, reveals that much of the distinction between foreground 
and background is lost in this projection. Consequently, a different type of 
projection is necessary to retain important information. The z-layer in which 
a region of an image is in focus is the z-layer with the least amount of 
blurring, i.e. with the highest contrast of intensity values within that area. 
Fig. \ref{fig:dapicontrastproj} shows this "maximum contrast projection". 
Comparing it to the maximum intensity projection shows that the individual 
cells are much easier to differentiate.

\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.45\textwidth}
    \includegraphics[width=\textwidth]{imgs/DAPI_maxproj.eps}
    \caption{Maximum intensity projection}
    \label{fig:dapimaxproj}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
    \includegraphics[width=\textwidth]{imgs/DAPI_contrastproj.eps}
    \caption{Maximum contrast projection}
    \label{fig:dapicontrastproj}
\end{subfigure}
\caption{A comparison of the maximum intensity projection and the maximum 
contrast projection of the images seen in Fig. \ref{fig:dapizstacks}}
\end{figure}

\section{Contrast Projection}

The contrast $C_{xy}$ at a point $(x,y)$ is defined here as the variance of 
all intensity scores in a defined area around this point. In the simplest 
case, this is is a rectangle centered around $(x,y)$ with side lengths of 
$2 \cdot w_x + 1$ and $2 \cdot w_y + 1$.

\begin{equation}
\label{eq:contrast1}
C_{xy} = \frac{1}{(2 \cdot w_x + 1) \cdot (2 \cdot w_y + 1)} 
        \sum_{k=x-w_x}^{x+w_x} 
        \sum_{l=y-w_y}^{y+w_y} (I_{kl} - \langle I_{xy} \rangle)^2
\end{equation}

\begin{equation}
\label{eq:contrast2}
C_{xy} = \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl}^2 - 2
        \cdot I_{kl} \cdot \langle I_{xy} \rangle + \langle I_{xy} \rangle ^2
\end{equation}

\begin{equation}
\label{eq:contrast3}
C_{xy} = {\langle I_{xy} \rangle}^2 - 2 \cdot \langle I_{xy} \rangle \cdot 
        \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl} + 
        \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl}^2
\end{equation}

where $\langle I_{xy} \rangle$ is the intensity of the image at the position 
$(x,y)$ and $\langle I_{xy} \rangle$ is the mean intensity of the window 
centered around $(x,y)$:

\begin{equation}
\label{eq:mean}
\begin{split}
\langle I_{xy} \rangle 
    &= \frac{1}{(2 \cdot w_x + 1) \cdot (2 \cdot w_y + 1)} 
        \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl} \\
    &= \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl}
\end{split}
\end{equation}

This mean intensity over a window can be quickly determined in R with the 
function \texttt{EBImage::filter2()} by applying a normalized, uniform filter 
to the image. For brevity, the normalizing constant, i.e. the area of the 
window, is abbreviated as $A = (2 \cdot w_x + 1) \cdot (2 \cdot w_y + 1)$.

A comparison of Eqs. \ref{eq:mean} and \ref{eq:contrast3} shows that the 
variance takes on the known form $Var = E[X^2] - (E[X])^2$:

\begin{equation}
\label{eq:contrast4}
\begin{split}
C_{xy} 
    &= {\langle I_{xy} \rangle}^2 + \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} 
        \sum_{l=y-w_y}^{y+w_y} I_{kl}^2 - 2 \cdot 
        {\langle I_{xy} \rangle}^2 \\
    &= \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} \sum_{l=y-w_y}^{y+w_y} I_{kl}^2 - 
        {\langle I_{xy} \rangle}^2 = \langle I^2_{xy} \rangle - 
        {\langle {I}_{xy} \rangle}^2
\end{split}
\end{equation}

where the first term of Eq. \ref{eq:contrast4} is the mean of the squared 
intensities over the same window as used for the mean of the intensities, 
$\langle I^2_{xy} \rangle = \frac{1}{A} \sum_{k=x-w_x}^{x+w_x} 
\sum_{l=y-w_y}^{y+w_y} I_{kl}^2$. This quantity can be calculated by the same 
method as for the intensity using \texttt{EBImage::filter2()}

\section{Application Example}

We will use the image data for the organoids already presented in the 
Introduction of this vignette. In the most general case, the data should be an 
array representing a stack of images, i.e. the first two dimensions correspond 
to the $x-$ and $y-$ dimensions of the image and the third dimension of the 
array corresponds to the number of images in the stack.

The $cells$ dataset represents a stack of 11 images, each with a size of 
261 x 251 pixels (see Fig. \ref{fig:dapizstacks})
<<>>=
library(MaxContrastProjection)
library(EBImage)
data(cells)
dim(cells)
@

The maximum contrast projection returns a single image with the same $x-$ and 
$y-$ dimensions as the input image stack

<<>>=
max_contrast = contrastProjection(imageStack = cells, w_x = 15, w_y = 15, 
                                smoothing = 5, brushShape = "box")
@

$w_x$ and $w_y$ are the radii of the area around each pixel over which to 
calculate the variance. The area of this window should be at least as large 
as the structures of interest in the image, otherwise artefacts may occur in 
the projection. Figs. \ref{fig:windowsizelarge} and \ref{fig:windowsizesmall} 
show the impact of the window size on the projection quality. If the parameter 
$w_y$ is not explicitly set, it is assumed to be equal to $w_x$.

<<>>=
max_contrast_large = contrastProjection(imageStack = cells[30:95, 55:115,], 
                                        w_x = 15, w_y = 15, smoothing = 5, 
                                        brushShape = "box")
max_contrast_small = contrastProjection(imageStack = cells[30:95, 55:115,], 
                                        w_x = 3, w_y = 3, smoothing = 5, 
                                        brushShape = "box")
@

\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.45\textwidth}
<<echo=FALSE>>=
display(max_contrast_large / max(cells), method = "raster")
@
    \caption{Large window with $w_x = w_y = 15$}
    \label{fig:windowsizelarge}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
<<echo=FALSE>>=
display(max_contrast_small / max(cells), method = "raster")
@
    \caption{Small window with $w_x = w_y = 3$}
    \label{fig:windowsizesmall}
\end{subfigure}
\caption{A comparison of window sizes for the maximum contrast projection on a 
segment of the \texttt{cells} data set. The small window size 
($w_x = w_y = 3$) shows some artefacts, i.e. several cells seem to be 
incorrectly projected from different z-stacks. The large window size 
($w_x = w_y = 15$) ensures that any single cell is projected from the same 
z-stack.}
\end{figure}

Similarily, the \texttt{smoothing} parameter applies a median filter onto the 
projection index map to avoid harsh jumps between z-stacks. It is unlikely 
that any given object would rapidly fluctuate between different focal plains 
and the median filter ensures the smoothness of the projection.

The \texttt{brushShape} parameter indicates the shape of the window to be used 
when determining the contrast at each pixel. By default ("box"), this is a 
rectangle with the side lengths $2 \cdot w_x + 1$ and $2 \cdot w_y + 1$. 
Currently supported is also a circle ("disc"). Depending on the geometry of 
the structures of interest in the image, different window shapes may have an 
impact on the quality of the projection. For radially symmetrical shapes, 
such as "disc", the parameter $w_y$ is assumed to be equal to $w_x$ regardless 
of whether it was explicitly defined.

It is also possible to retrieve the actual contrast values as well as the 
projection index map. The contrast stack is simply the contrast, i.e. 
variance, for every single pixel of every image in the input image stack. 
The projection index map then effectively indicates the z-layer with the 
maximum contrast for every pixel in the $x-y$-plane (after optional 
smoothing).

<<>>=
contrastStack = getContrastStack(imageStack = cells, w_x = 15, w_y = 15, 
                                brushShape = "box")
indexMap = getIndexMap(contrastStack = contrastStack, smoothing = 5)
max_contrast_fromMap = projection_fromMap(imageStack = cells, 
                                        indexMap = indexMap)
@

The maximum contrast projection can then be executed with just the image 
stack as well as the index map using \texttt{projection{\_}fromMap}. In 
fact, any projection index map can be used here so long as all its values 
lie between $1$ and the number of z-stacks in the input image stack.

A problem of the contrast projection is that if an object lies in multiple 
focal planes, the contrast projection generates artifical boundaries between 
these individual areas. This is an unavoidable consequence of imaging an 
object in discrete layers. Fig. \ref{fig:interpolation} shows what this looks 
like. In principle, the projection index map is blurred, so that boundaries 
between focal regions are softened. The projection is then a weighted linear 
combination of two values. For example, if the projection index map at a 
given pixel $(x, y)$ has a value of $7.4$, then the value of that pixel in the 
projection is $I_{xy} = 0.6 \cdot I_{xy}^{(7)} + 0.4 \cdot I_{xy}^{(8)}$, 
where $I_{xy}^{(j)}$ is the intensity of the pixel in the $j$-th layer of 
image stack.

\begin{figure}[!htb]
\centering
\begin{subfigure}[b]{0.45\textwidth}
    \includegraphics[width=\textwidth]{imgs/no_interpolation.eps}
    \caption{No interpolation}
    \label{fig:interpolation_no}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
    \includegraphics[width=\textwidth]{imgs/with_interpolation.eps}
    \caption{With interpolation}
    \label{fig:interpolation_yes}
\end{subfigure}
\caption{A comparison of the effect of linear interpolation on the maximum 
contrast projection. Note that the central area in Fig. 
\ref{fig:interpolation_no} is distinctly separated from the outside of the 
organoid. This separation is undesirable as it may lead to issues with image 
segmentation later on.}
\label{fig:interpolation}
\end{figure}

Lastly, this package also includes the possibility for other projection types. 
Currently, these are the maximum and minimum intensity projection, the mean 
and median intensity projection, the standard deviation projection, which 
calculates the standard deviation of the z-stack at each pixel in the 
$x-y$-plane, and the sum projection, which simply adds together all the values 
of the z-stack at each pixel in the $x-y$-plane.

<<>>=
max_intensity_proj = intensityProjection(imageStack = cells, projType = "max")
min_intensity_proj = intensityProjection(imageStack = cells, projType = "min")
mean_intensity_proj = intensityProjection(imageStack = cells, 
                                        projType = "mean")
median_intensity_proj = intensityProjection(imageStack = cells, 
                                            projType = "median")
sd_intensity_proj = intensityProjection(imageStack = cells, projType = "sd")
sum_intensity_proj = intensityProjection(imageStack = cells, projType = "sum")
@

\section{Session Info}
<<sessioninfo, results='asis'>>=
toLatex(sessionInfo())
@

\end{document}