%\VignetteIndexEntry{EBImage: image processing and analysis toolkit for R}
%\VignetteDepends{}
%\VignetteKeywords{image processing, visualization}
%\VignettePackage{EBImage}

\documentclass[11pt]{article}
%% times,10pt,twocolumn
\usepackage{graphicx}
\usepackage{verbatim}
\topmargin -1.5cm        % read Lamport p.163
\oddsidemargin -0.04cm   % read Lamport p.163
\evensidemargin -0.04cm  % same as oddsidemargin but for left-hand pages
\textwidth 16.59cm
\textheight 21.94cm

\usepackage{color}
\usepackage{hyperref}
\definecolor{darkblue}{rgb}{0.2,0.0,0.4}

\newcommand{\lib}[1]{{\mbox{\normalfont\textsf{#1}}}}
\newcommand{\file}[1]{{\mbox{\normalfont\textsf{'#1'}}}}
\newcommand{\R}{{\mbox{\normalfont\textsf{R}}}}
\newcommand{\EBImage}{{\mbox{\normalfont\textsf{EBImage}}}}
\newcommand{\Rfunction}[1]{{\mbox{\normalfont\texttt{#1}}}}
\newcommand{\Robject}[1]{{\mbox{\normalfont\texttt{#1}}}}
\newcommand{\Rpackage}[1]{{\mbox{\normalfont\textsf{#1}}}}
\newcommand{\Rclass}[1]{{\mbox{\normalfont\textit{#1}}}}
\newcommand{\code}[1]{{\mbox{\normalfont\texttt{#1}}}}

\newcommand{\email}[1]{\mbox{\href{mailto:#1}{\textcolor{darkblue}{\normalfont{#1}}}}}
\newcommand{\web}[2]{\mbox{\href{#2}{\textcolor{darkblue}{\normalfont{#1}}}}}

\usepackage{theorem}
\theoremstyle{break} \newtheorem{Ex}{Exercise}
\newenvironment{Sol}{\tiny}{}

%\usepackage[baseurl={http://www.ebi.ac.uk/~osklyar/EBImage/},pdftitle={EBImage - Image Processing Toolkit For R},pdfauthor={Oleg Sklyar},pdfsubject={EBImage},pdfkeywords={image processing},pagebackref,bookmarks,colorlinks,linkcolor=darkblue,citecolor=darkblue,pagecolor=darkblue,raiselinks,plainpages,pdftex]{hyperref}

\bibliographystyle{abbrvnat}
\begin{document}

\begin{figure}
%\vspace*{.1in}
\begin{center}
\scalebox{0.2}{\includegraphics{logo.png}}
\end{center}
\end{figure}

\title{Analysing Microscopy Screens with \Rpackage{EBImage}}
\author{Oleg Sklyar and Wolfgang Huber\\
    \emph{European Bioinformatics Institute - EMBL}\\
    \emph{Wellcome Trust Genome Campus}\\
    \emph{Cambridge CB10 1SD, England}\\
    \email{osklyar@ebi.ac.uk} }

\maketitle


The \R{}-package \Rpackage{EBImage} provides functionality to perform {\em image processing} and {\em image analysis} on large sets of images in a programmatic fashion using the \R{} language.

We use the term {\em image analysis} to describe the extraction of numeric features ({\em image descriptors}) from images and image collections~\cite{R:Russ:2002}. Image descriptors can then be used for statistical analysis, such as classification, clustering and hypothesis testing, using the resources of \R{} and its contributed packages.

Image analysis is not an easy task, and the definition of image descriptors depends on the problem. Analysis algorithms need to be adapted correspondingly. We find it desirable to develop and optimize such algorithms in conjunction with the subsequent statistical analysis, rather than as separate tasks. This is one of our motivations for writing \Rpackage{EBImage} as an \R{} package.

We use the term {\em image processing} for operations that turn images into images with the goals of enhancing, manipulating, sharpening, denoising or similar~\cite{R:Russ:2002}. While some image processing is often needed as a preliminary step for image analysis, image processing is not the primary goal of the package. We focus on methods that do not require interactive user input, such as selecting image regions with a pointer etc. Whereas interactive methods can be very effective for small sets of images, they tend to have limited throughput and reproducibility.

\Rpackage{EBImage} uses the \cite{R:ImageMagick:2006} library API's to implement much of its functionality in image processing and input/output operations.

%--------------------------------------------------
\section*{Cell-based assays}
%--------------------------------------------------

Advances in automated microscopy have made it possible to conduct large scale cell-based assays with image-type phenotypic readouts. In such an assay, cells are grown in the wells of a microtitre plate (often a 96- or 384-well format is used) under a condition or stimulus of interest. Each well is treated with one of the reagents from the screening library and the cells' response is monitored, for which in many cases certain proteins of interest are antibody-stained or labeled with a GFP-tag~\cite{CarpenterSabatini2004, Wiemann2004, MoffatSabatini2006, R:Neumann:2006}.

The resulting imaging data can be in the form of two-dimensional (2D) still images, three-dimensional (3D) image stacks or image-based time courses. Such assays can be used to screen compound libraries for the effect of potential drugs on the cellular system of interest. Similarly, RNA interference (RNAi) libraries can be used to screen a set of genes (in many cases the whole genome) for the effect of their loss of function in a certain biological process~\cite{R:Boutros:2004}.

%--------------------------------------------------
\section*{Importing and handling images}
%--------------------------------------------------

Images in \Rpackage{EBImage} are stored in objects of \lib{S4} class \Rclass{Image} that extends the \R{} \Rclass{array} class. The class encapsulates the most important image properies like size, resolution and colour mode, which are accessible by means of \emph{accessor methods}. It is recommended not to access slot data directly as their structure might change in the future releases and only use accessor methods. The \code{get} and \code{set} accessor methods for the \Rclass{Image}class are: \Rfunction{colorMode}, \Rfunction{imageData}, \Rfunction{fileName}, \Rfunction{compression}, \Rfunction{resolution}, and the \code{get}-only \Rfunction{features}. The colour mode can be either \code{Grayscale} or \code{TrueColor} with the corresponding integer constants defined in the package. See help pages for the corresponding functions for more detail.

New images can be created with the standard \R{} constructor \Rfunction{new}, or using the wrapper function \Rfunction{Image}. The following example code generates a \normalfont\textsf{100x100} pixel grayscale image with black and white vertical stripes \footnote{All examples in this vignette are based on version 1.9.12 (RC) of the package. Older versions have similar functions, that can have slightly different names.}:

\begin{verbatim}
> im <- Image(0.0, c(100,100), Grayscale)
> im[ c(1:20, 40:60, 80:100), , ] = 1
\end{verbatim}

As mentioned above, for input/output operations \Rpackage{EBImage} uses \lib{ImageMagick} API, which support reading and writing of more than 95 image formats including JPEG, TIFF and PNG. The package can read and write multi-page images (image stacks or 3D images) or process multiple files simultaneously. For example, the following code demonstrates how to read all PNG files in the working directory into a single object of class \code{Image}, convert them to grayscales and save the output as a single multi-page TIFF file:

\begin{verbatim}
> files <- dir(pattern=".png")
> im <- read.image(files, TrueColor)
> img <- channel(im, "gray")
> write.image(img, "single_multipage.tif")
\end{verbatim}

\noindent The package also provides an interactive way to select and load (multiple) images using the \Rfunction{choose.image} function and an alternative way of specifying the file name to write as in the following example:

\begin{verbatim}
> img <- choose.image()
> fileName(img) <- "single_multipage.tif"
> write.image(img)
\end{verbatim}

If the output format does not support multiple frames (stacks), \lib{ImageMagick} will automatically extend the provided filename with counter suffixes \code{-0, -1, ...} placed either before or after the file extension depending on your installation of \lib{ImageMagick}. Alternatively, one can specify the same number of files as image frames in call to \Rfunction{write.image} or \Rfunction{fileName}. File names supplied to \Rfunction{write.image} will take precedence over

Besides operations on local image files, anonymous HTTP and FTP protocols are supported. The package can read from both and it can write to FTP only. These protocols are supported internally by \lib{ImageMagick} and do not use \R{} connections.

%--------------------------------------------------
\section*{Displaying images}
%--------------------------------------------------

\Rpackage{EBImage} implements the method \Rfunction{display} to show images on screen interactively. By default (and if the package was compiled with \lib{GTK+} support) \Rpackage{EBImage} uses the \lib{GTK+} implementation of this function. It allows to display as many images simultaneously as the user needs and its interface is easily expandable programmatically. The alternative implementation uses the \lib{ImageMagick} API's \Rfunction{display} function. This implementation is only intended for systems where \lib{GTK+} interface is not available and can only show one image at a time. The \code{GTK+} interface is known to be working on all supported architectures, including \lib{MS Windows}.

In some cases users may want to use the generic \R{} \Rfunction{image} function, the package provides a simple wrapper, which keeps the image aspect ration and allows to specify which frame has to be displayed as in the following example:

\begin{verbatim}
> display(abc)  # 3 images, interactively
> image(abc, 2) # displays the 2nd frame
\end{verbatim}

%--------------------------------------------------
\section*{Image processing}
%--------------------------------------------------

The \lib{ImageMagick} library provides a number of image processing routines, so called {\em filters}. Many of those are ported to \R{} in \Rpackage{EBImage}. The missing ones can be added as opportunity arises. Additional image processing routines implementing distance map and watershed transforms and others have been added to the package to power our work on cell-based assays.

Filters are implemented as \lib{S4} methods for class \Rclass{Image}. The return value is usually an object of class \Rclass{Image} as well. One can divide all the filters into four categories: image enhancement, segmentation and transformation, and colour correction. Some examples are listed below.

\begin{description}
    \setlength{\itemsep}{0in}
    \item[channel] enables colour mode and colour channel conversions
    \item[sharpen, umask] generate sharpened versions of an image
    \item[gblur] applies the Gaussian blur to an image softening sharp edges and noise
    \item[thresh] segments a grayscale image into a binary black-and-white image by the adaptive thresholding algorithm
    \item[opening, closing] sequentially use erosion and dilation operators to enhance edges of objects in binary images and to reduce noise
    \item[distmap] performs a Euclidean distance transform of a binary image, also known as {\em distance map}. On a distance map, values of pixels indicate how far are they away from the nearest background
    \item[normalize] shifts and scales colours of images to a specified range, normally $[0,1]$ for grayscales
    \item[resample, resize] changes image size, proportionally or arbitrary
    \item[watershed] performs a \emph{watershed}-like image segmentation (used also for object detection)
\end{description}

The storage mode of grayscale images is \Rclass{double} (or \Rclass{numeric}), and \R{} functions that work with arrays can be directly applied to grayscale images. This includes the arithmetic functions, subsetting, histograms, Fast Fourier transformation, (local) regression etc. For example, the sharpened image in Figure~\ref{figure:sharpening}{\bf c} can be obtained by subtracting the slightly blurred, scaled in colour version of the original image (Figure~\ref{figure:sharpening}{\bf b}) from its source in Figure~\ref{figure:sharpening}{\bf a}. All pixels that become negative after subtraction are then re-set to background. The source image is obtained by subsetting the original microscopic image. Hereafter, variables in the code are given the same literal names as the corresponding image labels (e.g. data of variable \code{a} are shown in Figure ~\ref{figure:sharpening}~{\bf a}, \code{b} -- in {\bf b}, and \code{C} -- in {\bf c}, etc).

\begin{verbatim}
> orig <- read.image("ch2.png")
> a <- orig[150:550, 120:520,]
> b <- blur(0.5 * a, 80, 5)
> C <- a - b
> C[C < 0] = 0
> C <- normalize(C)
\end{verbatim}

One can think of this code as of a naive but fast and effective version of the {\em unsharp mask} filter; a more sophisticated implementation from the \lib{ImageMagick} library is provided by the \Rfunction{umask} function in the package.

\begin{figure}
\vspace*{.1in}
\begin{center}
\scalebox{0.75}{\includegraphics{umask.jpg}}
\end{center}
\caption{\label{figure:sharpening}
    Implementation of a simple {\em unsharp mask} filter: (a) source image, (b) blurred colour-scaled image, (c) sharpened image after normalization}
\end{figure}

Some of the image analysis routines in \Rpackage{EBImage} assume grayscale data in the interval $[0,1]$, but there are no formal restrictions on the range.

The storage mode of true colour images is \Rclass{integer}, and we use the three lowest bytes to store the red (R), green (G) and blue (B) values, each in the \code{unsigned char} (of \lib{C}) range of $[0, 255]$. Because of this, application of arithmetic and other functions on true colour images is generally meaningless; although they can be useful in some special cases, as shown in the example code in the following section. The support for true colour images is included mainly to enhance the display of the analysis results and to support further extensions of the package by third parties. Most analysis routines require grayscale data though.

The following code demonstrates how grayscale images recorded using three different microscope filters (Figure~\ref{figure:channels}~{\bf a}, {\bf b} and {\bf c}) can be put together into a single {\em false-colour} representation (Figure~\ref{figure:channels}~{\bf d}), and conversely, how a single false-colour image can be decomposed into its individual channels.

\begin{figure}
\vspace*{.1in}
\begin{center}
\scalebox{0.75}{\includegraphics{channels.jpg}}
\end{center}
\caption{\label{figure:channels}
    Composing a false-colour image (d) from a set of grayscale microscopy images for three different luminescent compounds: (a) -- DAPI, (b) -- tubulin and (c) -- phalloidin
}
\end{figure}

\begin{verbatim}
> files <- c("ch1.png","ch2.png","ch3.png")
> orig <- read.image(files, Grayscale)
> abc <- orig[150:550, 120:520,]
> a <- channel(abc[,,1], "asgreen")  # RGB
> b <- channel(abc[,,2], "asred")    # RGB
> d <- a + b + channel(abc[,,1], "asblue")
> C <- channel(d, "blue")            # gray
\end{verbatim}

%------------------------------------------------------------
\section*{Analysing an RNAi screen}
%------------------------------------------------------------

Consider an experiment in which images like in Figure~\ref{figure:channels} were recorded for each of $\approx$ 20,000 genes, using a whole-genome RNAi library to test the effect of gene knock-down on cell viability and appearance. Among the image descriptors of interest are the number, position, size, shape and the fluorescent intensities of cells and nuclei.

As a side effect of the watershed transform, the function \Rfunction{watershed}, \Rpackage{EBImage} identifies separable objects in an image. For this particular purpose it is advisable to use the result of the distance map transform (\Rfunction{distmap}) as input. Object descriptors can be obtained by calling \Rfunction{getObjects} on the result. One can see that an image has data of detected object by printing it (typing its name and hitting \code{Enter}): this information will be included in the object description. The features are a list of matrices, one per image frame, each containing different objects in its rows and different object descriptors in the columns. If a reference grayscale image was supplied to the watershed routine, object gray value intensities will be calculated as well.

For every gene, the image analysis workflow looks, therefore, as follows: load and normalize images, perform image segmentation, enhance the segmented images by morphological opening and closing, generate distance maps and use them to identify cells and nuclei and to extract image descriptors, and, finally, generate image previews with the identified objects marked.

Object descriptors can then be analysed statistically to cluster genes by their phenotypic effect, generate a list of genes that should be studied further in more detail (hit list), e.~g. genes that have a specific phenotypic effect of interest, etc. The image previews can be used to verify and audit the performance of the algorithm through visual inspection.

A schematic implementation is illustrated in the following example code and in Figure~\ref{figure:analysis}. Here we omit the step of nuclei detection (object \Robject{x1}), from where the matrix of nuclei coordinates (object \Robject{seeds}) is retrieved to serve as starting points for the cell detection. The nuclei detection is done analogously to the cell detection without specifying starting points, seeds.

\begin{figure}
\vspace*{.1in}
\begin{center}
\scalebox{0.75}{\includegraphics{analysis.jpg}}
\end{center}
\caption{\label{figure:analysis}
Illustration of the object detection algorithm: (a) -- {\em sqrt}-brightened combined image of nuclei (DAPI from Figure~\ref{figure:channels}{\bf a}) and cells (phalloidin from Figure~\ref{figure:channels}{\bf c}); (b) -- image {\bf a} after {\em blur} and {\em adaptive thresholding}; (c) -- image {\bf b} after morphological {\em opening} followed by {\em closing}; (d) -- normalized {\em distance map} generated from image {\bf c}; (e) -- outlines of cells detected using \Rfunction{watershed} drawn on top of the RGB image from Figure~\ref{figure:channels}{\bf d}; (f) -- colour-mapped cells and nuclei as detected with \Rfunction{watershed} (one unique colour per object)
}
\end{figure}

\begin{verbatim}
> for (X in genes) {
+  files <- dir(pattern=X)
+  orig <- read.image(files)
+  abc <- normalize(orig)
+  i1 <- abc[,,1]
+  i2 <- abc[,,2]
+  i3 <- abc[,,3]
+  a <- sqrt(normalize(i1 + i3))
+  b <- thresh( blur(a, 4, 2), 50, 50)
+  C <- opening(b, morphKern(7))
+  C <- closing(C, morphKern(7))
+  d <- distmap(C)  # shown normalized
+ ## x1 <- watershed(...) # nuclei detection
+  features(x1) <- getObjects(x1, i1)
+  x2 <- propagate(x1, a, d)
+  features(x2) <- getObjects(x2, i3)
+  rgb <- channel(i1, "asgreen") + channel(i2, "asred")
+  rgb <- rgb + channel(i3, "asblue")
+  e <- paintObjects(x2, rgb)
+  f <- paintObjects(x1, paintObjects(x2, i3) )
+ }
\end{verbatim}
%$

Note that here we adopted the \emph{record-at-a-time} approach: image data, which can be huge, are stored on a mass-storage device and are loaded into RAM in portions of just a few images at a time.

\section*{Summary}

\Rpackage{EBImage} brings image processing and image analysis capabilities to \R. Its focus is the programmatic (non-interactive) analysis of large sets of similar images, such as those that are obtained in cell-based assays for gene function via RNAi knock-down. Image descriptors extracted as the result of analysis can be analysed further using existing \R{}-functionality in machine learning (clustering, classification) and hypothesis testing.

Our future developments in image analysis will focus primarily on more accurate object detection and on algorithms for feature/descriptor extraction, and eventually on image registration, alignment and object tracking. Algorithms for the statistical analysis of image descriptors will be developed as part of a separate package that uses \Rpackage{EBImage} for image processing and analysis. In addition, one can imagine many other useful features, for example, support for more \lib{ImageMagick} functions. Contributions or collaborations on these or other topics are welcome.

%--------------------------------------------------
\section*{Acknowledgements and Availability}
%--------------------------------------------------

We thank F.~Fuchs and M.~Boutros for providing their miscroscopy data and for many stimulating discussions about the technology and the biology, and the European Bioinformatics Institute (EBI), Cambridge, UK, for financial support.

\Rpackage{EBImage} is available from \lib{BioConductor} (\url{www.bioconductor.org}). The package was tested to work on 32 and 64 bit UNIX and Linux systems, on both PowerPC and Intel based \lib{MacOS} and on \lib{MS Windows}. The package requires external dependencies, \lib{ImageMagick} and \lib{GTK+}, therefore, please read the included file \file{INSTALL} in the package tarball for installation instructions on different systems.

Parts of this text for older versions of \Rpackage{EBImage} were used in \R{}-News~\cite{Sklyar2006}.


\begin{thebibliography}{1}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{{\tt #1}}\fi

\bibitem[Sklyar2006]{Sklyar2006}
  O.~Sklyar, W.~Huber.
  \newblock Image Analysis for microscopy Screens.
  \newblock {\em R-News}, 6/5:12--16, 2006.

\bibitem[Boutros2004]{R:Boutros:2004}
  M.~Boutros, A.~Kiger, S.~Armknecht, \textit{et al.}
  \newblock Genome-wide RNAi analysis of cell growth and viability in Drosophila.
  \newblock {\em Science}, 303:832--835, 2004.

%\bibitem[Fuchs and Boutros(2006)]{R:Fuchs+Boutros:2006}
%    F.~Fuchs and M.~Boutros.
%    \newblock Cellular phenotyping by RNAi.
%    \newblock {\em Brief Funct Genomic Proteomic.}, 5:52--56, 2006.

\bibitem[Carpenter2004]{CarpenterSabatini2004}
  A.~E.~Carpenter and D.M.~Sabatini.
  \newblock Systematic genome-wide screens of gene function.
  \newblock {\em Nature Reviews Genetics}, 5:11--22, 2004.

%\bibitem[Carpenter et al.(in preparation)]{R:Carpenter:prep}
%    A.~E.~Carpenter, T.~R.~Jones, M.~Lamprecht, \textit{et al.}
%    \newblock CellProfiler: image analysis for high throughput microscopy.
%    \newblock {\em In preparation}.
%    \newblock URL \url{http://www.cellprofiler.org/}

\bibitem[ImageMagick]{R:ImageMagick:2006}
    ImageMagick: software to convert, edit, and compose images.
    \newblock {\em Copyright: ImageMagick Studio LLC}, 1999-2006.
    \newblock URL \url{http://www.imagemagick.org/}

%\bibitem[Jones et al.(2005)]{R:Jones:2005}
%    T.~R.~Jones and A.~E.~Carpenter and P.~Golland.
%    \newblock Voronoi-based segmentation of cells on image manifolds.
%    \newblock {\em ICCV Workshop on Computer Vision for Biomedical
%    Image Applications}, 535--543, 2005 (conference proceedings).

%\bibitem[Lotufo and Zampirolli(2001)]{R:Lotufo+Zampirolli:2001}
%    R.~Lotufo and F.~Zampirolli.
%    \newblock Fast multidimensional parallel Euclidean distance transform based on mathematical %morphology.
%    \newblock {\em SIBGRAPI-2001/Brazil}, 100--105, 2001.

\bibitem[Moffat2006]{MoffatSabatini2006}
  J.~Moffat and D.M.~Sabatini.
  \newblock Building mammalian signalling pathways with RNAi screens.
  \newblock {\em Nature Reviews Mol. Cell Biol.}, 7:177--187, 2006.

\bibitem[Neumann2006]{R:Neumann:2006}
    B.~Neumann, M.~Held, U.~Liebel, \textit{et al.}
    \newblock High-throughput RNAi screening by time-lapse imaging of live human cells.
    \newblock {\em Nature Mathods}, 3\penalty0 (5):\penalty0 385--390, 2006.

\bibitem[Russ2002]{R:Russ:2002}
    J.~C.~Russ.
   \newblock The image processing handbook -- 4th ed.
   \newblock CRC Press, {\em Boca Raton}. 732 p., 2002

\bibitem[Wiemann2004]{Wiemann2004}
   S.~Wiemann, D.~Arlt, W.~Huber,  \textit{et al.}
   \newblock From ORFeome to biology: a functional genomics pipeline.
   \newblock {\em Genome Res.} 14(10B):2136--2144, 2004.

\end{thebibliography}

\begin{center}
\scalebox{0.2}{\includegraphics{logo.png}}
\end{center}

\end{document}
