%\VignetteIndexEntry{Introduction to EBImage}
%\VignetteKeywords{image processing, visualization}
%\VignettePackage{EBImage}

\documentclass[10pt,a4paper]{article}

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

\RequirePackage{amsfonts,amsmath,amstext,amssymb,amscd}
\RequirePackage{graphicx}
\RequirePackage{verbatim}
\RequirePackage{hyperref}

\newcommand{\code}[1]{\Rcode{#1}} % mapping to BiocStyle

\hypersetup{
 baseurl={http://www.bioconductor.org/packages/release/bioc/html/EBImage.html},%
 pdfsubject={EBImage},%
 pdfkeywords={image processing}%
}

\SweaveOpts{keep.source=TRUE,eps=FALSE}


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

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

\title{Introduction to \Rpackage{EBImage},\\
an image processing and analysis toolkit for R}
\renewcommand{\thetitle}{Introduction to \Rpackage{EBImage}} % short title used in headers
\author{Gregoire Pau, Oleg Sklyar, Wolfgang Huber\\\email{gpau@ebi.ac.uk}}

\maketitle

\tableofcontents


\section{Reading/displaying/writing images}

The package \Biocpkg{EBImage} is loaded by the following command.
<<library,results=hide>>=
library("EBImage")
@
<<display-hack,echo=FALSE>>=
display = function(...) if (interactive()) EBImage::display(...)
@

The function \Rfunction{readImage} is able to read images from files
or URLs. Current supported image formats are JPEG, PNG and TIFF.

<<readImage1>>=
f = system.file("images", "lena.png", package="EBImage")
lena = readImage(f)
@

Images can be displayed using the function \Rfunction{display}. Pixel
intensities should range from 0 (black) to 1 (white).
<<display>>=
display(lena)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{lena.jpeg}
\includegraphics[width=0.48\textwidth]{lenac.jpeg}
\caption{\code{lena}, \code{lenac}}
\end{center}
\end{figure}

Color images or images with multiple frames can also be read with
\Rfunction{readImage}.
<<readImage2>>=
lenac = readImage(system.file("images", "lena-color.png", package="EBImage"))
display(lenac)
nuc = readImage(system.file('images', 'nuclei.tif', package='EBImage'))
display(nuc)
@
<<readImage2h,echo=FALSE>>=
writeImage(nuc, 'nuc.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{nuc-0.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-1.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-2.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-3.jpeg}
\caption{\code{nuc}}
\end{center}
\end{figure}

Images can be written with \Rfunction{writeImage}. The file format
is deduced from the file name extension. This is useful to convert
image formats, here from PNG format to JPEG format.

<<writeImage>>=
writeImage(lena,  'lena.jpeg', quality=85)
writeImage(lenac, 'lenac.jpeg', quality=85)
@

\pagebreak
\newpage
\section{Image objects and matrices}

The package \Rpackage{EBImage} uses the class \code{Image} to store and process
images. Images are stored as multi-dimensional arrays containing the
pixel intensities. All EBImage functions are also able to work with
matrices and arrays.

<<print>>=
print(lena)
@

As matrices, images can be manipulated with all R mathematical operators.
This includes \code{+} to control the brightness of an image, \code{*} to
control the contrast of an image or \code{\^} to control the gamma
correction parameter.

<<math1>>=
lena1 = lena+0.5
lena2 = 3*lena
lena3 = (0.2+lena)^3
@
<<math1h, echo=FALSE>>=
writeImage(lena1, 'lena1.jpeg', quality=85)
writeImage(lena2, 'lena2.jpeg', quality=85)
writeImage(lena3, 'lena3.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{lena.jpeg}
\includegraphics[width=0.24\textwidth]{lena1.jpeg}
\includegraphics[width=0.24\textwidth]{lena2.jpeg}
\includegraphics[width=0.24\textwidth]{lena3.jpeg}
\caption{\code{lena}, \code{lena1}, \code{lena2}, \code{lena3}}
\end{center}
\end{figure}

Others operators include \code{[} to crop images, \code{<} to threshold
images or \code{t} to transpose images.

<<math2>>=
lena4 = lena[299:376, 224:301]
lena5 = lena>0.5
lena6 = t(lena)
print(median(lena))
@
<<math2h, echo=FALSE>>=
writeImage(lena4, 'lena4.jpeg', quality=85)
writeImage(lena5, 'lena5.png')
writeImage(lena6, 'lena6.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{lena.jpeg}
\includegraphics[width=0.24\textwidth]{lena4.jpeg}
\includegraphics[width=0.24\textwidth]{lena5.png}
\includegraphics[width=0.24\textwidth]{lena6.jpeg}
\caption{\code{lena}, \code{lena4}, \code{lena5}, \code{lena6}}
\end{center}
\end{figure}

Images with multiple frames are created using \code{combine} which merges images.
<<combine>>=
lenacomb = combine(lena, lena*2, lena*3, lena*4)
display(lenacomb)
@
<<combineh, echo=FALSE>>=
writeImage(lenacomb, 'lenacomb.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{lenacomb-0.jpeg}
\includegraphics[width=0.24\textwidth]{lenacomb-1.jpeg}
\includegraphics[width=0.24\textwidth]{lenacomb-2.jpeg}
\includegraphics[width=0.24\textwidth]{lenacomb-3.jpeg}
\caption{\code{lenacomb}}
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Spatial transformations}

Specific spatial image transformations are done with the functions
\code{resize}, \code{rotate}, \code{translate} and the functions
\code{flip} and \code{flop} to reflect images.

<<spatial>>=
lena7 = rotate(lena, 30)
lena8 = translate(lena, c(40, 70))
lena9 = flip(lena)
@
<<spatialh, echo=FALSE>>=
writeImage(lena7, 'lena7.jpeg', quality=85)
writeImage(lena8, 'lena8.jpeg', quality=85)
writeImage(lena9, 'lena9.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{lena.jpeg}
\includegraphics[width=0.24\textwidth]{lena7.jpeg}
\includegraphics[width=0.24\textwidth]{lena8.jpeg}
\includegraphics[width=0.24\textwidth]{lena9.jpeg}
\caption{\code{lena}, \code{lena7}, \code{lena8}, \code{lena9}}
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Color management}

The class \code{Image} extends the base class \code{array} and uses
the \code{colormode} slot to store how the color information
of the multi-dimensional data should be handled.

As an example, the color image \code{lenac} is a 512x512x3 array, with a
\code{colormode} slot equals to \code{Color}. The object is understood as
a color image by \Rpackage{EBImage} functions.
<<print>>=
print(lenac)
@

The function \code{colorMode} can access and change the value of the slot
\code{colormode}, modifying the rendering mode of an image. In the next
example, the \code{Color} image \code{lenac} with one frame is changed
into a \code{Grayscale} image with 3 frames, corresponding to the red,
green and blue channels. The function \code{colorMode} does not change
the content of the image but changes only the way the image is rendered
by \Rpackage{EBImage}.

<<colorMode>>=
colorMode(lenac) = Grayscale
display(lenac)
@
<<colorModeh,echo=FALSE>>=
writeImage(lenac, 'lenac.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{lenac.jpeg}
\includegraphics[width=0.24\textwidth]{lenac-0.jpeg}
\includegraphics[width=0.24\textwidth]{lenac-1.jpeg}
\includegraphics[width=0.24\textwidth]{lenac-2.jpeg}
\caption{\code{lenac}, rendered as a \code{Color} image and as a \code{Grayscale} image with 3 frames (red channel, green channel, blue channel)}
\end{center}
\end{figure}

The color mode of image \code{lenac} is reverted back to \code{Color}.
<<colorMode2>>=
colorMode(lenac) = Color
@

The function \code{channel} performs colorspace conversion and can convert
\code{Grayscale} images into \code{Color} ones both ways and can extract color
channels from \code{Color} images. Unlike \code{colorMode}, \code{channel}
changes the pixel intensity values of the image. The function \code{rgbImage}
is able to combine 3 \code{Grayscale} images into a \code{Color} one.

<<channel>>=
lenak = channel(lena, 'rgb')
lenak[236:276, 106:146, 1] = 1
lenak[236:276, 156:196, 2] = 1
lenak[236:276, 206:246, 3] = 1
lenab = rgbImage(red=lena, green=flip(lena), blue=flop(lena))
@
<<channelh,echo=FALSE>>=
writeImage(lenak, 'lenak.jpeg', quality=85)
writeImage(lenab, 'lenab.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{lenak.jpeg}
\includegraphics[width=0.48\textwidth]{lenab.jpeg}
\caption{\code{lenak}, \code{lenab}}
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Image filtering}

Images can be linearly filtered using \code{filter2}. \code{filter2}
convolves the image with a matrix filter. Linear filtering is useful
to perform low-pass filtering (to blur images, remove noise, ...) and
high-pass filtering (to detect edges, sharpen images, ...). Various filter
shapes can be generated using \code{makeBrush}.

<<filter>>=
flo = makeBrush(21, shape='disc', step=FALSE)^2
flo = flo/sum(flo)
lenaflo = filter2(lenac, flo)

fhi =  matrix(1, nc=3, nr=3)
fhi[2,2] = -8
lenafhi = filter2(lenac, fhi)
@
<<filterh,echo=FALSE>>=
writeImage(lenaflo, 'lenaflo.jpeg', quality=85)
writeImage(lenafhi, 'lenafhi.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{lenaflo.jpeg}
\includegraphics[width=0.48\textwidth]{lenafhi.jpeg}
\caption{Low-pass filtered \code{lenaflo} and high-pass filtered \code{lenafhi}}
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Morphological operations}
Binary images are images where the pixels of value 0 constitute the
background and the other ones constitute the foreground. These images
are subject to several non-linear  mathematical operators called
morphological operators, able to \code{erode} and \code{dilate} an
image.
<<morpho>>=
ei = readImage(system.file('images', 'shapes.png', package='EBImage'))
ei = ei[110:512,1:130]
display(ei)

kern = makeBrush(5, shape='diamond')
eierode = erode(ei, kern)
eidilat = dilate(ei, kern)
@
<<morphoh,echo=FALSE>>=
writeImage(ei, 'ei.png')
writeImage(kern, 'kern.png')
writeImage(eierode, 'eierode.png')
writeImage(eidilat, 'eidilat.png')
@
\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{ei.png}\\{\tiny $ $}\\
\includegraphics[width=0.48\textwidth]{eierode.png}\\{\tiny $ $}\\
\includegraphics[width=0.48\textwidth]{eidilat.png}
\caption{\code{ei} ; \code{eierode} ; \code{eidilat}}
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Segmentation}

Segmentation consists in extracting objects from an image. The function
\code{bwlabel} is a simple function able to extract every connected sets
of pixels from an image and relabel these sets with a unique increasing
integer. \code{bwlabel} can be used on binary images and is useful after
thresholding.

<<segmentation>>=
eilabel = bwlabel(ei)
cat('Number of objects=', max(eilabel),'\n')

nuct = nuc[,,1]>0.2
nuclabel = bwlabel(nuct)
cat('Number of nuclei=', max(nuclabel),'\n')
@
<<segmentationh,echo=FALSE>>=
writeImage(eilabel/max(eilabel), 'eilabel.png')
writeImage(nuclabel/max(nuclabel), 'nuclabel.png')
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{ei.png}
\includegraphics[width=0.48\textwidth]{eilabel.png}
\end{center}
\caption{\code{ei}, \code{eilabel/max(eilabel)}}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{nuc-0.jpeg}
\includegraphics[width=0.48\textwidth]{nuclabel.png}
\caption{\code{nuc[ , ,1]}, \code{nuclabel/max(nuclabel)}}
\end{center}
\end{figure}

Since the images \code{eilabel} and \code{nuclabel} range from 0 to the
number of object they contain (given by \code{max(eilabel)} and
\code{max(nucabel)}), they have to be divided by these number before
displaying, in order to fit the [0,1] range needed by \code{display}.

The grayscale top-bottom gradient observable in \code{eilabel} and \code{nuclabel}
is due to the way \code{bwlabel} labels the connected sets, from top-left to
bottom-right.

Adaptive thresholding consists in comparing the intensity of pixels with
their neighbors, where the neighborhood is specified by a filter matrix.
The function \code{thresh} performs a fast adaptive thresholding of an image
with a rectangular window while the combination of \code{filter2} and \code{<}
allows a finer control. Adaptive thresholding allows a better segmentation
when objects are close together.

<<segmentation2>>=
nuct2 =  thresh(nuc[,,1], w=10, h=10, offset=0.05)
kern = makeBrush(5, shape='disc')
nuct2 = dilate(erode(nuct2, kern), kern)
nuclabel2 = bwlabel(nuct2)
cat('Number of nuclei=', max(nuclabel2),'\n')
@
<<segmentation2h,echo=FALSE>>=
writeImage(nuclabel2/max(nuclabel2), 'nuclabel2.png')
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{nuc-0.jpeg}
\includegraphics[width=0.48\textwidth]{nuclabel2.png}
\caption{\code{nuc[ , ,1]}, \code{nuclabel2/max(nuclabel)} }
\end{center}
\end{figure}

\pagebreak
\newpage
\section{Object manipulation}

Objects, defined as sets of pixels with the same unique integer value
can be outlined and painted using \code{paintObjects}. Some holes are
present in objects of \code{nuclabel2} which can be filled using
\code{fillHull}.

<< manip >>=
nucgray = channel(nuc[,,1], 'rgb')
nuch1 = paintObjects(nuclabel2, nucgray, col='#ff00ff')
nuclabel3 = fillHull(nuclabel2)
nuch2 = paintObjects(nuclabel3, nucgray, col='#ff00ff')
@
<<maniph,echo=FALSE>>=
writeImage(nuch1, 'nuch1.jpeg', quality=85)
writeImage(nuch2, 'nuch2.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{nuch1.jpeg}
\includegraphics[width=0.48\textwidth]{nuch2.jpeg}
\caption{\code{nuch1}, \code{nuch2} }
\end{center}
\end{figure}

A broad variety of objects features (basic, image moments, shape, Haralick features) can
be computed using \code{computeFeatures}. In particular, object
coordinates are computed with the function \code{computeFeatures.moment}.
<< manip2 >>=
xy = computeFeatures.moment(nuclabel3)[, c("m.cx", "m.cy")]
xy[1:4,]
@

\pagebreak
\newpage
\section{Cell segmentation example}

This is a complete example of segmentation of cells (nucleus + cell bodies) using
the functions described before and the function \code{propagate}, able to
perform Voronoi-based region segmentation.

Images of nuclei and cell bodies are first loaded:
<< cs1 >>=
nuc = readImage(system.file('images', 'nuclei.tif', package='EBImage'))
cel = readImage(system.file('images', 'cells.tif', package='EBImage'))
img = rgbImage(green=1.5*cel, blue=nuc)
@
<<cs1h, echo=FALSE>>=
writeImage(cel, 'cel.jpeg', quality=85)
writeImage(img, 'img.jpeg', quality=85)
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{nuc-0.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-1.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-2.jpeg}
\includegraphics[width=0.24\textwidth]{nuc-3.jpeg}
\caption{\code{nuc}}
\end{center}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{cel-0.jpeg}
\includegraphics[width=0.24\textwidth]{cel-1.jpeg}
\includegraphics[width=0.24\textwidth]{cel-2.jpeg}
\includegraphics[width=0.24\textwidth]{cel-3.jpeg}
\caption{\code{cel}}
\end{center}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{img-0.jpeg}
\includegraphics[width=0.24\textwidth]{img-1.jpeg}
\includegraphics[width=0.24\textwidth]{img-2.jpeg}
\includegraphics[width=0.24\textwidth]{img-3.jpeg}
\caption{\code{img}}
\end{center}
\end{figure}


Nuclei are first segmented using \code{thresh}, \code{fillHull}, \code{bwlabel}
and \code{opening}, which is an \code{erosion} followed by a \code{dilatation}.
<< cs2 >>=
nmask = thresh(nuc, w=10, h=10, offset=0.05)
nmask = opening(nmask, makeBrush(5, shape='disc'))
nmask = fillHull(nmask)
nmask = bwlabel(nmask)
@
<<cs2h, echo=FALSE>>=
writeImage(nmask/max(nmask), 'nmask.png')
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{nmask-0.png}
\includegraphics[width=0.24\textwidth]{nmask-1.png}
\includegraphics[width=0.24\textwidth]{nmask-2.png}
\includegraphics[width=0.24\textwidth]{nmask-3.png}
\caption{\code{nmask/max(nmask)}}
\end{center}
\end{figure}

Cell bodies are segmented using \code{propagate}.
<< cs3 >>=
ctmask = opening(cel>0.1, makeBrush(5, shape='disc'))
cmask = propagate(cel, seeds=nmask, mask=ctmask)
@
<<cs3h, echo=FALSE>>=
writeImage(cmask/max(cmask), 'cmask.png')
@

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.24\textwidth]{cmask-0.png}
\includegraphics[width=0.24\textwidth]{cmask-1.png}
\includegraphics[width=0.24\textwidth]{cmask-2.png}
\includegraphics[width=0.24\textwidth]{cmask-3.png}
\caption{\code{cmask/max(cmask)}}
\end{center}
\end{figure}

Cells are outlined using \code{paintObjects}.

<< cs4 >>=
res = paintObjects(cmask, img, col='#ff00ff')
res = paintObjects(nmask, res, col='#ffff00')
@
<<cs4h, echo=FALSE>>=
writeImage(res, 'res.jpeg', quality=85)
@


\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.48\textwidth]{res-0.jpeg}
\includegraphics[width=0.48\textwidth]{res-1.jpeg}\\{\tiny$ $}\\
\includegraphics[width=0.48\textwidth]{res-2.jpeg}
\includegraphics[width=0.48\textwidth]{res-3.jpeg}
\caption{Final segmentation \code{res}}
\end{center}
\end{figure}

\end{document}