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

\documentclass[10pt,a4paper]{article}

\RequirePackage{amsfonts,amsmath,amstext,amssymb,amscd}
\usepackage{graphicx}
\usepackage{verbatim}
\usepackage{hyperref}
\usepackage{color}
\definecolor{darkblue}{rgb}{0.2,0.0,0.4}

\topmargin -1.5cm
\oddsidemargin -0cm   % read Lamport p.163
\evensidemargin -0cm  % same as oddsidemargin but for left-hand pages
\textwidth 17cm
\textheight 24.5cm

\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}

%% Uncomment the second to show solutions
%\SweaveOpts{echo=FALSE,results=hide}
%\SweaveOpts{echo=FALSE,results=hide,eval=FALSE}
\SweaveOpts{echo=TRUE,results=verbatim,eval=TRUE}

\begin{document}

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

%------------------------------------------------------------
\title{Analysis of biological images with \lib{EBImage}}
%------------------------------------------------------------
\author{Oleg Sklyar, Wolfgang Huber\\\email{osklyar@ebi.ac.uk}}
\maketitle


In this manual the \EBImage{} package is used to analyse a set of images acquired in a large-scale RNAi microscopy screen and to extract numerical descriptors of the cells in these images.  The descriptors define biological phenotypes. The descriptors can be further analysed statistically e.\,g.\ to classify genes by their phenotypic effect.

The images used in this manual are available in the package lab can be downloaded from\\ \web{http://www.ebi.ac.uk/\~{}osklyar/BioC2007/data}.

The images are 16-bit grayscale TIFF's and were measured at two different wavelengths corresponding to two different stainings: DNA content measured in the green channel (suffix G) and one of a cytoplasm protein in red (suffix R). 

<<loadlib,echo=FALSE,results=hide,eval=TRUE>>=
library(EBImage)
@

%%--------------------------------------------
\section{Handling images}
%%--------------------------------------------

\EBImage{} provides two functions to read images, \Rfunction{readImage} and \Rfunction{chooseImage}. They both allow users to specify whether the result should be in grayscale or RGB mode.  

An interactive window for loading images is provided by the \Rfunction{chooseImage} function, which is available if the package was compiled with GTK+ support (always the case on Windows): 

<<example.choose.image,echo=TRUE,results=hide,eval=FALSE>>=
x = chooseImage(TrueColor)
@

The \Rfunction{readImage} function can read local or network files via HTTP or anonymous FTP protocols. Multiple files can be loaded into a single object, an image stack. If the source image has multiple frames they all will be read into a stack. Single, multi-frame or multiple images read are stored in objects of class \code{Image}:
 
<<example.read.image,echo=TRUE,results=verbatim,eval=TRUE>>=
ddir <- paste( system.file(package="EBImage"), "images", sep="/" )
fG = paste(ddir, dir(path=ddir, pattern="_G.tif"), sep="/")
iG = readImage(fG[1], Grayscale)
class(iG)
dim(iG)

fR = paste(ddir, dir(path=ddir, pattern="_R.tif"), sep="/")
iR = readImage(fR[1])
@

Images can be read from remote URL's as in the following example:

<<example.readRemote,echo=TRUE,results=hide,eval=FALSE>>=
baseurl = "http://www.ebi.ac.uk/~osklyar/BioC2007/data/"
a = readImage( paste(baseurl, c("Gene1_R.tif", "Gene2_R.tif"), sep="") )
@

%%--------------------------------------------
\section{Exploring data}
%%--------------------------------------------

Objects of the class \Rclass{Image} can be displayed using either the \Rfunction{display} function or the standard \Rfunction{image} method for the \Rclass{Image} class: 

<<example.imageDisplay,echo=TRUE,results=hide,eval=FALSE>>=
image(iG[,,1])
display(iR)
animate(iR)
@

Beside displaying one can use \Rfunction{hist} and standard \Rfunction{print} methods to explore the data contained in images. For example, for the image stack \code{iR} consisting of 4 images, individual histograms are obtained with the following code (results shown in Figure~1):

<<echo=FALSE,results=hide,eval=TRUE>>=
pdf("01.pdf", width=10, height=10)
@
<<example.histograms,echo=TRUE,results=hide,eval=TRUE>>=
split.screen(c(2,2))
for ( i in 1:4 ) {
  screen(i)
  hist(iR[,,i], xlim=c(0,1))
}
close.screen(all=TRUE)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
dev.off()
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.75\textwidth]{01.pdf}
    \caption{\label{figure:histsbefore} Histograms of individual images in \file{Gene1\_{}R.tif}. The frames have significantly different dynamic ranges and need to be normalized before intensity-based descriptors can be used.}
  \end{center}
\end{figure}

Try additionally the following commands to investigate the structure of the data:

<<example.print.methods,echo=TRUE,results=hide,eval=FALSE>>=
str(iR)
print(iR)
dim(iR)
range(iR)
@
%
Try also to do the same for \code{iG}, the images of the DNA content.

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

The \Rfunction{normalize} function allows you to perform normalization of image data to a given range (the default is $[0,1]$). Images in a stack can be normalized either separately from each other or simultaneously as parts of a large dataset. In our case, we want to normalize each image separately. Try the following two normalizations and compare the outputs:
 
<<example.normalizations,echo=TRUE,results=verbatim,eval=TRUE>>=
iRn.wrong = normalize(iR, separate=FALSE)
apply(iRn.wrong, 3, range)
iRn = normalize(iR, separate=TRUE)
apply(iRn, 3, range)
iGn = normalize(iG, separate=TRUE)
@

\EBImage{} provides a set of functions to manipulate the quality of images like \Rfunction{enhance}, \Rfunction{contrast}, \Rfunction{despeckle} or \Rfunction{denoise}. Try them out on the images in \code{iRn}. Note that such manipulations can be useful for visualisation, but they are not appropriate for quantitative analyses, such as when intensity-derived features are used as phenotypic descriptors.
 
We can also apply non-linear transformations to improve image quality. Another reason for some of these transformations is to achieve a higher object contrast against the background, which turns useful in consequent image segmentation. Please bear in mind the same note: images transformed in such a way cannot (should not) be used to extract intensity-based descriptors and are assumed for visualisation or segmentation purposes.

Let us define several functions for transforming the range $[0,1]$ into itself.

<<example.modifFunctions,echo=TRUE,results=verbatim,eval=TRUE>>=
modif1 = function(x) sin((x-1.2)^3)+1
modif2 = function(x, s) (exp(-s*x) - 1) / (exp(-s) - 1)
modif3 = function(x) x^1.5
@

The graphs of these functions are shown in Figure~\ref{figure:modif} on page~\pageref{figure:modif}. Non-linear transformations of the range $[0,1]$ into itself are often used in image processing to prepare an image for subsequent operations. I remark as a sidenote that this is also one of the most useful filters in Photoshop or GIMP when it comes to processing holiday photos. 

<<modifPlots,echo=FALSE,results=hide,eval=TRUE>>=
pdf("02.pdf", width=6,height=6)
  x = (0:100)/100
  plot(x, x, type="l", xlim=c(0,1), ylim=c(0,1), col="gray", lwd=1, lty=2,
     xlab="source intensity", ylab="target intensity")
  lines(x, modif2(x,4), col="#10508B", lwd=2)
  text(0.3, 0.83, "modif2(s=4)", col="#10508B")
  lines(x, modif2(x,-4), col="#4E82B5", lwd=2)
  text(0.7, 0.16, "modif2(s=-4)", col="#4E82B5")
  lines(x, sqrt(x), col="#24A072", lwd=2)
  text(0.2, 0.37, "sqrt")
  lines(x, x^1.5, col="#421C80", lwd=2)
  text(0.52, 0.3, "modif3", col="#421C80")
  lines(x, modif1(x), col="red", lwd=2)
  text(0.52, 0.6, "modif1", col="red")
dev.off()
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{02.pdf}
    \caption{\label{figure:modif} Non-linear transformations for data in the range $[0,1]$ onto $[0,1]$. The \emph{curves} filter in Photoshop or The GIMP.}
  \end{center}
\end{figure}

We will use \code{modif1} to improve the contrast and proceed as follows. The modification used (if any) is problem specific and is determined by the quality of original images:

<<example.do.modif1,echo=TRUE,results=hide,eval=TRUE>>=
iGn = modif1(iGn)
iRn = modif1(iRn)
@

The histograms after the normalization and the above transformation are given in Figure~\ref{figure:histsafter}:

<<echo=FALSE,results=hide,eval=TRUE>>=
pdf("03.pdf", width=10,height=10)
@
<<example.histograms2,echo=TRUE,results=hide,eval=TRUE>>=
split.screen(c(2,2))
for ( i in 1:4 ) {
  screen(i)
  hist(iRn[,,i], xlim=c(0,1))
}
close.screen(all=TRUE)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
dev.off()
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.75\textwidth]{03.pdf}
    \caption{\label{figure:histsafter} Histograms of individual images in \file{Gene1\_{}R.tif} after the normalization.}
  \end{center}
\end{figure}

%%--------------------------------------------
\section{Colour modes}
%%--------------------------------------------

For visual representations, images can be converted between the grayscale and true colour modes; a grayscale image can also be converted into one of the RGB channels if required and channels can be added together as in the following example:

<<example.grayscaleToRGB,echo=TRUE,results=hide,eval=TRUE>>=
iRG = channel(iRn, "asred") + channel(iGn, "asgreen")
@
<<example.displayConversions,echo=TRUE,results=hide,eval=FALSE>>=
display( iRG )
display( channel(iRG, "gray") )
display( channel(iRG, "asred") )
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(iRG[,,1], "04.png")
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{04.png}
    \caption{\label{figure:screenshot} False-colour representation of
      the data from two independent channels, red and green.}
  \end{center}
\end{figure}

The \Rfunction{channel} function can also be used to convert vectors containing colour data from one format to another. Grayscale to RGB integers: 

<<example.channel1,echo=TRUE,results=verbatim,eval=TRUE>>=
ch=channel(c(0.2,0.5,0.7), "rgb") 
ch
sprintf("%X", ch)
@ 
\noindent Grayscale to X11 hexadecimal color strings:
<<example.channel2,echo=TRUE,results=verbatim,eval=TRUE>>=
channel(c(0.2,0.5,0.7), "x11") 
@ 
\noindent Color strings to RGB:
<<example.channel3,echo=TRUE,results=verbatim,eval=TRUE>>=
channel(c("red","green","#0000FF"), "rgb") 
@ 
\noindent RGB to grayscale:
<<example.channel4,echo=TRUE,results=verbatim,eval=TRUE>>=
channel(as.integer(c(3355443, 8355711, 11711154)), "gray") 
# channel(c(3355443L, 8355711L, 11711154L), "gray") # in newer R
@

Images can be stored in files as in the following example for \code{iRG}:

<<example.write.image,echo=TRUE,results=verbatim,eval=FALSE>>=
writeImage(iRG, "composite.tif")
compression(iRG) = "JPEG"
writeImage(iRG, "composite.jpg", quality=98)
@

%%--------------------------------------------
\section{Creating images and further data manipulation}
%%--------------------------------------------

Images can be created either using the default constructor \Rfunction{new} for class \code{Image} or using a wrapper function, \Rfunction{Image}:

<<echo=FALSE,results=verbatim,eval=TRUE>>=
pdf("05.pdf", width=6,height=4)
par(mar=c(1,1,1,1))
@
<<example.runif.Image,echo=TRUE,results=verbatim,eval=TRUE>>=
a = Image(runif(200*100), c(200,100))
image(a)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
dev.off()
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{05.pdf}
    \caption{\label{figure:imx} An image of uniform random numbers.}
  \end{center}
\end{figure}

\noindent One can also use the data of other images to create new ones:

<<example.newFromExisting,echo=TRUE,results=verbatim,eval=TRUE>>=
a = Image(iRn, dim(iRn), colormode=colorMode(iRn))
@
<<echo=TRUE,results=verbatim,eval=FALSE>>=
display(a)
@

\noindent In the above example, it is not necessary to supply the dimensions as those will be automatically detected from the data if that has the \code{dim} method defined.

Transformations can be performed using \Rfunction{resize}, \Rfunction{rotate} etc. Try the following two to increase the image by a factor of 1.3 or to rotate it by 15 degrees:

<<example.resize,echo=TRUE,results=verbatim,eval=TRUE>>=
a = resize(iRn, dim(iRn)[1]*1.3)
@
<<echo=TRUE,results=verbatim,eval=FALSE>>=
display(a)
@
<<example.rotate,echo=TRUE,results=verbatim,eval=TRUE>>=
a = rotate(iRn, 15)
@
<<echo=TRUE,results=verbatim,eval=FALSE>>=
display(a)
@

Simple data manipulations can be performed by subsetting. For example, the following lines represent simple thresholding:

<<example.masking,echo=TRUE,results=verbatim,eval=TRUE>>=
a = iRn
a[ a > 0.6 ] = 1.0
a[ a <= 0.6] = 0.0
@

\noindent On a grayscale image the values of $0$ and $1$ in the above example create a black-and-white mask. If now we want to mark the background e.g. in blue and foreground in red, this can be achieved as follows:

<<example.colorMask,echo=TRUE,results=verbatim,eval=TRUE>>=
b = channel(a, "rgb")
b[ a >= 0.1] = channel("red", "rgb")
b[ a < 0.1 ] = channel("#114D90", "rgb")
@
<<echo=TRUE,results=verbatim,eval=FALSE>>=
display(b)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(b[,,1], "06.png")
rm(a,b,x)
gc()
@

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{06.png}
    \caption{\label{figure:marked} Colour-marked binary mask.}
  \end{center}
\end{figure}

%%--------------------------------------------
\section{Image segmentation and image analysis}
%%--------------------------------------------

The purpose of segmentation is to mask the objects of interest from the background prior to identifying them. The quality of the segmentation will generally define the quality of the subsequent object indexing and feature extraction. We need something better than the mask in Figure~\ref{figure:marked}.
 
For further object indexing We will make use of the fact that we have two images corresponding to the same location of the microscope -- one of the nuclei staining and one for the cytoplasm protein. Assuming that every cell has a nucleus we will use indexed nuclei (after segmentation, indexing and feature extraction) to index cells. Therefore, we start start with segmenting nuclei, images \code{iG}.
 
The function \Rfunction{thresh} provides an implementation of an adaptive threshold filter that takes into account inequalities in background intensity across the image. For \code{iG} the segmented image can be obtained as follows:

<<example.threshG,echo=TRUE,results=verbatim,eval=TRUE>>=
mask = thresh(iGn, 15, 15, 0.002)
@

The parameters \Robject{w}, \Robject{h} of the function are related to the size of the objects we expect to find in the image: objects of different size would require adjustment of these parameters. The \Robject{offset} is determined by the local intensity differences. Try using different parameters and compare the segmentation. The quality of segmentation is vital for good quality of object indexing and feature extraction, therefore it worth to spend time tuning the parameters. For comparable results, images across the experiment should be segmented using the same parameter set. This might lead to artifacts in segmentation in some cases, but will ensure that same types of cells look similar in different images! The result of the above segmentation is shown in Figure~\ref{figure:nt1}.
 
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(mask[,,1], "07.png")
@

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{07.png}
    \caption{\label{figure:nt1} Preliminary nuclei segmentation.}
  \end{center}
\end{figure}

Some further smoothing of the mask is necessary. A useful set of instruments for this is provided by \emph{mathematical morphology} implemented in the morphological operators \Rfunction{dilate}, \Rfunction{erode}, \Rfunction{opening} and \Rfunction{closing}:

<<example.dilate,echo=TRUE,results=verbatim,eval=TRUE>>=
mk3 = morphKern(3)
mk5 = morphKern(5)
mask = dilate(erode(closing(mask, mk5), mk3), mk5)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(mask[,,1], "08.png")
@

\noindent Here, several operators were used sequentially. You can observe the results of each of these operators separately by looking at the intermediate images. You can also try different kernels, i.\,e.\ different parameters for the function \Rfunction{morphKern}. The current result is shown in Figure~\ref{figure:nt2}.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{08.png}
    \caption{\label{figure:nt2} Nuclei segmentation after smoothing and noise removal.}
  \end{center}
\end{figure}

As the next step, one needs to index regions in the segmented image that correspond to the different objects.  A classic algorithm for this is computing the distance map transform followed by the \code{watershed} transform (see Figure~\ref{figure:wsn1}):

<<example.watershedOnDistmap,echo=TRUE,results=verbatim,eval=TRUE>>=
sG = watershed( distmap(mask), 1.5, 1)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(normalize(sG[,,1]), "09.png")
@

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{09.png}
    \caption{\label{figure:wsn1} Nuclei segmentation by watershed 
  (before artifact removal).}
  \end{center}
\end{figure}

Finally, when we are happy with the result of the watershed transform, we can remove nuclei that are either too small or too dark or fall on the edge of the images, etc (see Figure~\ref{figure:wsn2}): 

<<example.removeBadG,echo=TRUE,results=verbatim,eval=TRUE>>=
ft = hullFeatures(sG)  ## need these for edge and perimeter
mf = moments(sG, iGn)   ## need these for intensity and size
for ( i in seq_along(ft) ) ft[[i]] = cbind(ft[[i]], mf[[i]])
sG = rmObjects(sG, 
        lapply(ft, 
          function(x) 
            which(x[,"h.s"] < 150 | x[,"h.s"] > 10000 | x[,"int"] < 30 |
                  0.4 * x[,"h.p"] < x[,"h.edge"] )
      ))
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(normalize(sG[,,1]), "10.png")
@

\noindent here \code{h.s} and \code{h.p} stand for the hull size and perimeter, \code{h.edge} for the number of pixels at the image edge and \code{int} for the intensity of the region (as returned by \Rfunction{hullFeatures} and \Rfunction{moments}). Investigate the structure of \code{ft} and \code{mf} and explain what kind of objects were removed.

What we have finally obtained is an \code{IndexedImage} for the nuclei, where each nucleus is given an index from $1$ to \code{max(sG)}. One can now directly use functions like \Rfunction{getFeatures} or \Rfunction{moments} etc.\ to obtain numerical descriptors of each nucleus.
 
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{10.png}
    \caption{\label{figure:wsn2} Nuclei segmentation by watershed 
    (after artefact removal).}
  \end{center}
\end{figure}

In principle, the same distance-map/watershed algorithm could be used to segment the cells, however we often find that neighbouring cells are touching and lead to segmentation errors. We can use the already identified nuclei as seed points to detect corresponding cells -- assuming that each cell has exactly one nucleus. This method falls short of detecting multi-nuclear cells, but it improves the quality of detection for all other cells tremendously. We start similarly to the nuclei segmentation, however instead of using \code{watershed}, we use \code{propagate}, supplying it with an \code{IndexedImage} of the seed points (nuclei). The function implements an elegant algorithm that produces a Voronoi segmentation using a metric that combines Euclidean distance and intensity differences between different pixels in the image:

<<example.propagate,echo=TRUE,results=verbatim,eval=TRUE>>=
mask = thresh(blur(iRn,4,1.5), 25, 25, 0.005)
mask = erode( erode( dilate(mask,mk5), mk5), mk3 )
sR = propagate( iRn, sG, mask, 1e-5, 1.6)
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(normalize(sR[,,1]), "11.png")
@

\noindent A weighting factor is used in propagate to either give more weight to the Euclidean distance or otherwise to the intensity-driven one. We use a very low value of $1e-5$ basically minimizing the effect of the Euclidean. Also please note that we used the \Rfunction{blur} filter to obtain the original mask. In case of cells we use seed points and we know where the cells are, therefore the mask we use is larger and smoother to accommodate more tiny settled changed in the overall image of every individual cell. The result is shown in Figure~\ref{figure:cells}.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{11.png}
    \caption{\label{figure:cells} Cell segmentation by the
      \Rfunction{propagate} function (before artifact removal).}
  \end{center}
\end{figure}

Again, some artifacts need to be removed. In all consequent computations, the number of objects in an indexed image (as obtained from \Rfunction{watershed} or \Rfunction{propagate}) is determined by finding the maximum value. Consider that this value is $N$ for \code{sR}. If the image contains pixels indexed with $N$, but is missing pixels with some other indexes, smaller than $N$, the corresponding objects will be identified with $0$ data, first of all $0$ size. $N$ can be smaller than the original number of nuclei as it could happen that for some nuclei no cells were identified. There can be many reasons for this: cells masked out or too small or too dark etc. In order to preserve the 1-t-1 match of nuclei to cells, \code{max(sG)} must be equal $N$, so we mask out all nuclei with indexes larger than $N$:

<<example.matching,echo=TRUE,results=verbatim,eval=TRUE>>=
for ( i in 1:dim(sR)[3] ) {
  x = sG[,,i]
  x[ x > max(sR[,,i]) ] = 0.0
  sG[,,i] = x
}
@

Now as we ensured the 1-to-1 match of nuclei to cells, we can remove cells that are too small or too large to be plausible, are on the edge of the image, are too dark, etc. We also remove the corresponding nuclei (the result is given in Figure~\ref{figure:wsn3}):
 
<<example.removeBadR,echo=TRUE,results=verbatim,eval=TRUE>>=
ft = hullFeatures(sR)
mf = moments(sR, iRn)
for ( i in seq_along(ft) ) ft[[i]] = cbind(ft[[i]], mf[[i]])
index = lapply(ft, 
           function(x) 
             which( x[,"h.s"] < 150 | x[,"h.s"] > 15000 | 
                    x[,"int"]/x[,"h.s"] < 0.1 | 
                    0.3 * x[,"h.p"] < x[,"h.edge"]
         ))
sR = rmObjects(sR, index)
sG = rmObjects(sG, index)
@

\noindent See above for the notations of the column names in \code{x}.

<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(normalize(sR[,,1]), "12.png")
@
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{12.png}
    \caption{\label{figure:wsn3} Cell segmentation by propagate (after artefact removal).}
  \end{center}
\end{figure}

Finally, having the indexed images for cells and nuclei, the full set of descriptors can be extracted using the \Rfunction{getFeatures} function:
 
<<example.getFeatures,echo=TRUE,results=verbatim,eval=TRUE>>=
sG = getFeatures(sG, iGn)
sR = getFeatures(sR, iRn)
nucl = do.call("rbind", features(sG))
cells = do.call("rbind", features(sR))
stopifnot(identical(dim(nucl), dim(cells)))
@
The resulting matrices have \Sexpr{nrow(nucl)} rows (one for each of cell/nucleus) and \Sexpr{ncol(nucl)} columns (one for each object descriptor).  

You can now try out the following visualisations, with the first one shown in Figure~\ref{figure:done}:

<<example.visualize,echo=TRUE,results=verbatim,eval=TRUE>>=
rgb = paintObjects(sR, iRG)
rgb = paintObjects(sG, rgb)

ct = tile( combine(stackObjects(sR, iRn)) )
nt = tile( combine(stackObjects(sG, iGn)) )
@
<<echo=FALSE,results=hide,eval=TRUE>>=
writeImage(rgb[,,1], "13.png")
@
The result is shown in Figure~\ref{figure:done}.
\begin{figure}[htbp]
  \begin{center}
    \includegraphics[width=0.5\textwidth]{13.png}
    \caption{\label{figure:done} Results of detection.}
  \end{center}
\end{figure}



\end{document}


