• Learn how to read and write images in and out of R

• Learn how images are represented in R and how to manipulate them

• Understand how to apply filters and transformations to images

• Apply these skills to microscopy images of cells to do segmentation and feature extraction

• Explore spatial distributions of the position of cells

Image processing and analysis toolbox for R

• Reading and writing of image files
• Interactive image viewer
• Image manipulation, transformation and filtering
• Object detection and feature extraction

Since Bioconductor 1.8 (2006)

Original developers: Oleg Sklyar Wolfgang Huber Mike Smith Gregoire Pau

Contributors: Joseph Barry Bernd Fischer Ilia Kats Philip A. Marais

library(EBImage)

f = system.file("images", "sample.png", package="EBImage")
img = readImage(f)

display(img)

Reading images

Images can be read from local files or URLs.

bioc = readImage("http://www.bioconductor.org/images/logo/jpg/bioconductor_logo_rgb.jpg")
display(bioc)

EBImage supports JPEG, PNG and TIFF file formats. For reading proprietary microscopy image data and metadata use RBioFormats.

Displaying images

• interactive JavaScript viewer
• R's build-in plotting device

The default display method can be set by options(EBImage.display).

options(EBImage.display = "raster")

display(img, method = "raster")
text(x = 20, y = 20, label = "Parrots", adj = c(0,1), col = "orange", cex = 2)

filename = "parrots.jpg"
dev.print(jpeg, filename = filename , width = dim(img)[1], height = dim(img)[2])

file.size(filename)

## [1] 37820

Supported file formats: JPEG, PNG and TIFF.

writeImage(img, "sample.jpeg", quality = 85)

writeImage(img, "sample.tiff")
writeImage(img, "sample_compressed.tiff", compression = "deflate")

files = list.files(pattern = "sample*")
data.frame(row.names=files, size=file.size(files))

##                          size
## sample.jpeg             49688
## sample.tiff            393384
## sample_compressed.tiff 293788

str(img)

## Formal class 'Image' [package "EBImage"] with 2 slots
##   ..@ .Data    : num [1:768, 1:512] 0.447 0.451 0.463 0.455 0.463 ...
##   ..@ colormode: int 0

getClassDef("Image")

## Class "Image" [package "EBImage"]
## 
## Slots:
##                           
## Name:      .Data colormode
## Class:     array   integer
## 
## Extends: 
## Class "array", from data part
## Class "structure", by class "array", distance 2
## Class "vector", by class "array", distance 3, with explicit coerce

dim(img)

## [1] 768 512

img

## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 768 512 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.4470588 0.4627451 0.4784314 0.4980392 0.5137255 0.5294118
## [2,] 0.4509804 0.4627451 0.4784314 0.4823529 0.5058824 0.5215686
## [3,] 0.4627451 0.4666667 0.4823529 0.4980392 0.5137255 0.5137255
## [4,] 0.4549020 0.4666667 0.4862745 0.4980392 0.5176471 0.5411765
## [5,] 0.4627451 0.4627451 0.4823529 0.4980392 0.5137255 0.5411765

imageData(img)[1:3, 1:6]

##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.4470588 0.4627451 0.4784314 0.4980392 0.5137255 0.5294118
## [2,] 0.4509804 0.4627451 0.4784314 0.4823529 0.5058824 0.5215686
## [3,] 0.4627451 0.4666667 0.4823529 0.4980392 0.5137255 0.5137255

hist(img)

range(img)

## [1] 0 1

f = system.file("images", "sample-color.png", package="EBImage")
imgcol = readImage(f)
display(imgcol)

print(imgcol, short = TRUE)

## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 768 512 3 
##   frames.total : 3 
##   frames.render: 1

nuc = readImage(system.file("images", "nuclei.tif", package="EBImage"))
print(nuc, short = TRUE)

## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4

display(nuc)

## The image contains more than one frame: only the first one is displayed.
## To display all frames use 'all = TRUE'.

display(nuc, method = "raster", all = TRUE)

Being numeric arrays, images can be conveniently manipulated by any of R's arithmetic operators.

Cropping

img = img[366:749, 58:441]

Inversion

img_neg = max(img) - img
display(img_neg)

Brightness, contrast, and gamma correction

img_comb = combine(
  img,
  img + 0.3,
  img * 2,
  img ^ 0.5
)

display( tile(img_comb, numberOfFrames(img_comb), lwd = 20, fg.col = "white") )

Binary images

Images which contain only two sets of pixels, which represent the background and the foreground pixels.

img_thresh = img > .5
display(img_thresh)

Binary images

Images which contain only two sets of pixels, which represent the background and the foreground pixels.

img_thresh

## Image 
##   colorMode    : Grayscale 
##   storage.mode : logical 
##   dim          : 384 384 
##   frames.total : 1 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6]
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] FALSE FALSE FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE FALSE

Translation

img_translate = translate(img, v = c(100, -50))
display(img_translate)

Rotation

img_rotate = rotate(img, angle = 30, bg.col = "white")
display(img_rotate)

Scaling

img_resize = resize(img, w = 512, h = 256)
display(img_resize)

img_resize = resize(img, 256)
display(img_resize)

Vertical and horizontal reflection

display( flip(img) )

display( flop(img) )

Described by a 3x2 transformation matrix \(\textbf{m}\)

Example: horizontal sheer mapping

angle = pi/6
m = matrix(c(1, -sin(angle), sin(angle)*dim(img)[2]/2, 0, 1, 0), nrow = 3, ncol = 2)
m

##      [,1] [,2]
## [1,]  1.0    0
## [2,] -0.5    1
## [3,] 96.0    0

img_affine = affine(img, m)
display(img_affine)

imgcol_t = transpose(imgcol)
display(imgcol_t)

Removing local artifacts or noise – average in a rectangular window

Generalization: weighted average using a weight function \(w\)

The weight function can be generated using makeBrush.

w = makeBrush(size = 31, shape = 'gaussian', sigma = 5)
plot(w[(nrow(w)+1)/2, ], ylab = "w", xlab = "")

Other available brush shapes: "box" (default), "disc", "diamond", "line"

2D linear convolution is implemented by filter2 (uses FFT)

Example: Smoothing an image using a Gaussian filter of width 5

img_flo = filter2(img, w)
display(img_flo)

Detecting edges and sharpening images

fhi = matrix(1, nrow = 3, ncol = 3)
fhi[2, 2] = -8
fhi

##      [,1] [,2] [,3]
## [1,]    1    1    1
## [2,]    1   -8    1
## [3,]    1    1    1

img_fhi = filter2(img, fhi)
display(img_fhi)

Compensate for spatial dependencies of the background signal

disc = makeBrush(21, "disc")
disc = disc / sum(disc)
nuc = getFrame(nuc, 1)
nuc_bg = filter2(nuc, disc)
offset = 0.02
nuc_thresh = (nuc - nuc_bg) > offset

img_comb = combine(nuc, nuc_bg, nuc_thresh)
display( tile(img_comb, numberOfFrames(img_comb), lwd = 20, fg.col = "white") )

Non-linear noise removal technique

• particularly effective in the case of speckle noise
• preserves edges

l = length(img)
n = l/10
img_noisy = img
img_noisy[sample(l, n)] = runif(n, min = 0, max = 1)
img_median = medianFilter(img_noisy, size = 1)

display( combine(img_noisy, img_median), all=TRUE)

Fluorescent microscopy images from two channels of HeLa cells.

dna = readImage(system.file("images", "nuclei.tif", package="EBImage"))
print(dna, short=TRUE)

## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4

tub = readImage(system.file("images", "cells.tif", package="EBImage"))
print(tub, short=TRUE)

## Image 
##   colorMode    : Grayscale 
##   storage.mode : double 
##   dim          : 510 510 4 
##   frames.total : 4 
##   frames.render: 4

display( combine(getFrame(dna, 3), getFrame(tub, 3)), all=TRUE )

rgb = rgbImage(green = 1.5 * tub, blue = dna)
displayTiles(rgb)

• Adaptive thresholding of the DNA channel
• Morphological opening and filling of holes
• Distance map computation and watershed transformation

Adaptive thresholding and morphological filtering

nmaskt = thresh(dna, w = 15, h = 15, offset = 0.05)
nmaskf = fillHull( opening(nmaskt, makeBrush(5, shape='disc')) )

display( combine(getFrame(nmaskt, 3), getFrame(nmaskf, 3)), all=TRUE )

Distance map computation

dmap = distmap(nmaskf)
range(dmap)

## [1]  0.00000 14.31782

display(normalize(dmap), frame = 3)

Watershed transformation

nmask = watershed(dmap, tolerance = 2)

display( combine(
  toRGB( getFrame(dna, 3) ), 
  colorLabels( getFrame(nmask, 3) )
), all=TRUE )

bwlabel - a simpler and faster function for segmenting connected objects

Identification of cell bodies by Voronoi partitioning using the segmented nuclei as seeds.

cmaskt = closing( gblur(tub, 1) > 0.105, makeBrush(5, shape='disc') )

cmask  = propagate(tub, seeds=nmask, mask=cmaskt, lambda = 0.001)

display( combine(
  toRGB( getFrame(cmaskt, 3) ), 
  colorLabels( getFrame(cmask, 3) )
), all=TRUE )

Voronoi diagram is a partitioning of a plane into regions based on distance to some specific points of that plane (seeds).

Voronoi tesselation on a Riemann manifold

• arbitrary shape (provided in mask)
• arbitrary metric (dependent on image gradient \(z\))

display( paintObjects(nmask,
            paintObjects(cmask, rgb, col = "magenta", thick = TRUE),
         col = "yellow", thick = TRUE), all = TRUE)

st = stackObjects(cmask, rgb)

display(st, all = TRUE)

Quantitative cellular descriptors

• basic statistics on pixel intensities computeFeatures.basic
• shape characteristics (area, perimeter, radius) computeFeatures.shape
• moments (center of mass, eccentricity, …) computeFeatures.moment
• Haralick textural features computeFeatures.haralick

head( computeFeatures.shape(cmask[,,1], tub[,,1]) )

##   s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1   4216         408      43.52530   14.254634    16.133188     73.48543
## 2   3079         263      34.60167    9.061043    19.146184     53.49633
## 3   3092         211      31.00933    4.636759    16.798016     38.91896
## 4   1091         161      21.16105    8.163347     7.496897     35.43516
## 5    355          78      10.65161    2.772033     4.568584     15.93175
## 6   2377         209      29.37053    8.030861    16.278131     44.85961

Using multivariate analysis methods we can:

• detect cell subpopulations (clustering)
• classify cells into pre-defined cell types or phenotypes (classification)
• based on the frequencies of the subpopulations compare different biological conditions

Lymph nodes

• immunologic filter composed of B cells, T cells, dendric cells (DC) and macrophages
• clinical significance in cancer staging (the sentinel lymph node)

Breast cancer lymph node biopsy (Setiadi et al. 2010). Cell typing was done by staining the cells with protein antibodies that provide specific signatures.

data(brcalymphnode, package = "BioC2015Oles")
head(brcalymphnode)

##     class    x     y
## 1 T_cells 6355 10382
## 2 T_cells 6356 10850
## 3 T_cells 6357 11070
## 4 T_cells 6357 11082
## 5 T_cells 6358 10600
## 6 T_cells 6361 10301

nrow(brcalymphnode)

## [1] 209462

table(brcalymphnode$class)

## 
##     T_cells       Tumor         DCs other_cells 
##      103681       27822         878       77081

Plot of the x and y positions of the T cells (left) and tumor cells (right)

par(mfrow = c(1, 2), mar = c(4.5, 4.5, 0.5, 0.5))
xlim = range(brcalymphnode$x)
ylim = range(brcalymphnode$y)
cols = c(`T_cells` = "dodgerblue4", `Tumor` = "darkmagenta")
for(i in seq_along(cols))
  plot(subset(brcalymphnode, class==names(cols)[i])[, c("x", "y")], 
       pch = ".", asp = 1, xlim = xlim, ylim = ylim, col = cols[i])

Marked spatial point process:

• set of isolated points located in a mathematical space (xlim, ylim)
• points marked with certain properties (factor class)

library("spatstat")

ln = with(brcalymphnode, 
  ppp(x = x, y = y, marks = class, xrange = xlim, yrange = ylim)
)

ln

## Marked planar point pattern: 209462 points
## Multitype, with levels = T_cells, Tumor, DCs, other_cells 
## window: rectangle = [3839, 17276] x [6713, 23006] units

chull = convexhull(ln)

par(mar = c(0, 0, 0, 0))
plot(Window(ln), main = NULL, asp = 1)
plot(chull, lty = 2, col = "lightblue", add = TRUE)

Window(ln) = chull
ln

## Marked planar point pattern: 209462 points
## Multitype, with levels = T_cells, Tumor, DCs, other_cells 
## window: polygonal boundary
## enclosing rectangle: [3839, 17276] x [6713, 23006] units

• Number of points lying in a circle of area \(a\) around a point \(p=(x,y)\) \[N(p, a)\]

• Local intensity \[\lambda(p) = \lim_{a\rightarrow 0} \frac{E[ N(p, a)]}{a}\]

• For a stationary process, i.e., homogeneous all over the region \[\lambda(p) = \text{const.}\] Then the intensity in an area \(A\) is proportional to the area \[E(N(\cdot, A)) = \lambda A\]

densities = solist( 
  `Diggle's edge correction` = density(subset(ln, marks=="Tumor"), diggle = TRUE),
  `No edge correction`       = density(subset(ln, marks=="Tumor"), edge = FALSE)
)
plot(densities, equal.ribbon = TRUE, col = topo.colors, main = "")

Estimate of the spatially-varying probability of a particular event type

rr = relrisk(ln, sigma = 250)

plot(rr, equal.ribbon = TRUE, col = topo.colors, nrows = 1, main = "")

Spatial clustering (or anticlustering) of points: whether and to what extent neighboring points are more/less dense than expected by chance.

Poisson process

The nearest neighbour (NN) distance distribution function \(G\)

Cumulative distribution function of the distance \(W\) from a typical random point to its nearest neighbor. For a Poisson process

\[G(r) = P(W \leq r) = 1-e^{-\lambda \pi r^2}\]

gln = Gest(ln)
gln

## Function value object (class 'fv')
## for the function r -> G(r)
## .....................................................................
##         Math.label      Description                                  
## r       r               distance argument r                          
## theo    G[pois](r)      theoretical Poisson G(r)                     
## han     hat(G)[han](r)  Hanisch estimate of G(r)                     
## rs      hat(G)[bord](r) border corrected estimate of G(r)            
## km      hat(G)[km](r)   Kaplan-Meier estimate of G(r)                
## hazard  hat(h)[km](r)   Kaplan-Meier estimate of hazard function h(r)
## theohaz h[pois](r)      theoretical Poisson hazard function h(r)     
## .....................................................................
## Default plot formula:  .~r
## where "." stands for 'km', 'rs', 'han', 'theo'
## Recommended range of argument r: [0, 20.998]
## Available range of argument r: [0, 52.443]

library(RColorBrewer)
par(mar = c(4.5, 4.5, 0.5, 0.5))

plot(gln, lty = 1, col = brewer.pal(4, "Set1"), main = "")

For a stationary point process \(\lambda K(r)\) is the expected number of additional points within a distance \(r\) of a given, randomly picked point.

Generalization to inhomogenous point processes

where \(d(p_i, p_j)\) is the distance between points \(p_i\) and \(p_j\), and \(e(p_i, p_j, r)\) is an edge correction factor.

For a Poisson point pattern, the theoretical value is \(L(r) = r\).

Lln = Linhom( subset(ln, marks=="T_cells") )
Lln

## Function value object (class 'fv')
## for the function r -> L[inhom](r)
## ...........................................................................
##            Math.label                
## r          r                         
## theo       L[pois](r)                
## border     {hat(L)[inhom]^{bord}}(r) 
## bord.modif {hat(L)[inhom]^{bordm}}(r)
##            Description                                      
## r          distance argument r                              
## theo       theoretical Poisson L[inhom](r)                  
## border     border-corrected estimate of L[inhom](r)         
## bord.modif modified border-corrected estimate of L[inhom](r)
## ...........................................................................
## Default plot formula:  .~r
## where "." stands for 'bord.modif', 'border', 'theo'
## Recommended range of argument r: [0, 694.7]
## Available range of argument r: [0, 694.7]

par(mar = c(4.5, 4.5, 0.5, 0.5))

plot(Lln, lty = 1, col = brewer.pal(3, "Set1"), main = "")

For a stationary Poisson process \(g(r) = 1\).

pcfln = pcf( Kinhom(subset(ln, marks=="T_cells")) )

plot(pcfln, lty = 1, log = "x")