## ----echo=FALSE, message=FALSE-------------------------------------------
library(BioC2015Oles)
library(knitr)

opts_chunk$set(error=FALSE)
set.seed(7)
.dpi = 100

## ---- message=FALSE, fig.width=768/.dpi, fig.height=512/.dpi, dpi=.dpi/2----
library(EBImage)

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

display(img)

## ---- fig.width=480/.dpi, fig.height=138/.dpi, dpi=.dpi, eval=FALSE------
#  bioc = readImage("http://www.bioconductor.org/images/logo/jpg/bioconductor_logo_rgb.jpg")
#  display(bioc)

## ---- eval=FALSE---------------------------------------------------------
#  options(EBImage.display = "raster")

## ---- fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2, results='hide'----
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)

## ------------------------------------------------------------------------
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))

## ------------------------------------------------------------------------
str(img)
getClassDef("Image")
dim(img)

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

## ---- fig.width=4, fig.height=4------------------------------------------
hist(img)
range(img)

## ---- fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
f = system.file("images", "sample-color.png", package="EBImage")
imgcol = readImage(f)
display(imgcol)
print(imgcol, short = TRUE)

## ---- echo=FALSE---------------------------------------------------------
nuc = readImage(system.file("images", "nuclei.tif", package="EBImage"))

## ---- fig.width=dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2----
nuc = readImage(system.file("images", "nuclei.tif", package="EBImage"))
print(nuc, short = TRUE)
display(nuc)

## ---- fig.width=dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi----
display(nuc, method = "raster", all = TRUE)

## ---- fig.width=384L/.dpi, fig.height=384L/.dpi, dpi=.dpi/2--------------
img = img[366:749, 58:441]

## ----negative, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
img_neg = max(img) - img
display(img_neg)

## ----arithmetic, fig.width=(4*dim(img)[1L]+100)/.dpi, fig.height=(dim(img)[2L]+40)/.dpi, dpi=.dpi/2----
img_comb = combine(
  img,
  img + 0.3,
  img * 2,
  img ^ 0.5
)

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

## ---- fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
img_thresh = img > .5
display(img_thresh)

## ------------------------------------------------------------------------
img_thresh

## ----translate, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
img_translate = translate(img, v = c(100, -50))
display(img_translate)

## ----rotate-pre, echo=FALSE----------------------------------------------
img_rotate = rotate(img, 30)

## ----rotate, fig.width=dim(img_rotate)[1L]/.dpi, fig.height=dim(img_rotate)[2L]/.dpi, dpi=.dpi/2----
img_rotate = rotate(img, angle = 30, bg.col = "white")
display(img_rotate)

## ----resize, fig.width=512/.dpi, fig.height=256/.dpi, dpi=.dpi/2---------
img_resize = resize(img, w = 512, h = 256)
display(img_resize)

## ----resize2, fig.width=256/.dpi, fig.height=256/.dpi, dpi=.dpi/2--------
img_resize = resize(img, 256)
display(img_resize)

## ----flipflop, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
display( flip(img) )
display( flop(img) )

## ---- fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
angle = pi/6
m = matrix(c(1, -sin(angle), sin(angle)*dim(img)[2]/2, 0, 1, 0), nrow = 3, ncol = 2)
m
img_affine = affine(img, m)
display(img_affine)

## ----transpose, fig.width=dim(imgcol)[2L]/.dpi, fig.height=dim(imgcol)[1L]/.dpi, dpi=.dpi/2----
imgcol_t = transpose(imgcol)
display(imgcol_t)

## ---- echo=FALSE---------------------------------------------------------
opts_knit$set(global.par=TRUE)

## ---- echo=FALSE---------------------------------------------------------
.par=par(mai = c(.45, .75, 0.05, 0.05))

## ----makeBrush, fig.width=3.8, fig.height=3.5, dev="svg"-----------------
w = makeBrush(size = 31, shape = 'gaussian', sigma = 5)
plot(w[(nrow(w)+1)/2, ], ylab = "w", xlab = "")

## ---- echo=FALSE---------------------------------------------------------
opts_knit$set(global.par=FALSE)

## ----lopass, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
img_flo = filter2(img, w)
display(img_flo)

## ----highpass, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2----
fhi = matrix(1, nrow = 3, ncol = 3)
fhi[2, 2] = -8
fhi
img_fhi = filter2(img, fhi)
display(img_fhi)

## ---- fig.width=(3*dim(nuc)[1L]+80)/.dpi, fig.height=(dim(nuc)[2L]+40)/.dpi, dpi=.dpi/2----
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") )

## ----medianFilter, fig.width=2*dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2, eval=TRUE----
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)

## ----echo=FALSE----------------------------------------------------------
leaf = readImage( system.file('images', 'leaf.png', package='BioC2015Oles') )
.leafDim = BioC2015Oles:::tilesDim(dim(leaf), 3)

## ---- fig.width=.leafDim[1]/.dpi, fig.height=.leafDim[2]/.dpi, dpi=.dpi, out.width=2*.leafDim[1], results='hide'----
leaf = readImage( system.file('images', 'leaf.png', package='BioC2015Oles') )
kern = makeBrush(size = 3)
morph = combine(
  leaf,
  erode(leaf, kern),
  dilate(leaf, kern)
)

displayTiles(morph)

## ---- fig.width=.leafDim[1]/.dpi, fig.height=.leafDim[2]/.dpi, dpi=.dpi, out.width=2*.leafDim[1]----
morph = combine(
  leaf,
  opening(leaf, kern),
  closing(leaf, kern)
)

displayTiles(morph) 

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

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/3----
display( combine(getFrame(dna, 3), getFrame(tub, 3)), all=TRUE )

## ---- fig.width=BioC2015Oles:::tilesDim(dim(dna), numberOfFrames(dna))[1L]/.dpi, fig.height=BioC2015Oles:::tilesDim(dim(dna), numberOfFrames(dna))[2L]/.dpi, dpi=.dpi/3----
rgb = rgbImage(green = 1.5 * tub, blue = dna)
displayTiles(rgb)

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2----
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 )

## ---- fig.width=dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2----
dmap = distmap(nmaskf)
range(dmap)

display(normalize(dmap), frame = 3)

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2----
nmask = watershed(dmap, tolerance = 2)

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

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2----
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 )

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=2*dim(nuc)[2L]/.dpi, dpi=.dpi/2----
display( paintObjects(nmask,
            paintObjects(cmask, rgb, col = "magenta", thick = TRUE),
         col = "yellow", thick = TRUE), all = TRUE)

## ---- fig.width=2*dim(nuc)[1L]/.dpi, fig.height=2*dim(nuc)[2L]/.dpi, dpi=.dpi/2, eval=TRUE----
st = stackObjects(cmask, rgb)

display(st, all = TRUE)

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

## ---- echo=FALSE---------------------------------------------------------
.lymphnode = readImage( system.file("images", "testscan_1_2_RGB99-4525D.jpg", package = "BioC2015Oles") )

## ---- echo=FALSE, fig.width=dim(.lymphnode)[1L]/(4*.dpi), fig.height=dim(.lymphnode)[2L]/(4*.dpi), dpi=.dpi----
display(.lymphnode)

## ------------------------------------------------------------------------
data(brcalymphnode, package = "BioC2015Oles")
head(brcalymphnode)
nrow(brcalymphnode)
table(brcalymphnode$class)

## ---- fig.width=7, fig.height=3.5----------------------------------------
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])

## ---- message=FALSE------------------------------------------------------
library("spatstat")

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

ln

## ---- fig.width=2, fig.height=2------------------------------------------
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

## ----density.ppp1, fig.width=5.25, fig.height=2.65-----------------------
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 = "")

## ---- fig.width=9.5, fig.height=2.8--------------------------------------
rr = relrisk(ln, sigma = 250)

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

## ----Gest----------------------------------------------------------------
gln = Gest(ln)
gln

## ---- fig.width=4, fig.height=4, message=FALSE---------------------------
library(RColorBrewer)
par(mar = c(4.5, 4.5, 0.5, 0.5))

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

## ---- eval=FALSE---------------------------------------------------------
#  Lln = Linhom( subset(ln, marks=="T_cells") )
#  Lln

## ---- echo=FALSE---------------------------------------------------------
data(Lln, package = "BioC2015Oles")
Lln

## ---- fig.width=4.5, fig.height=4.5--------------------------------------
par(mar = c(4.5, 4.5, 0.5, 0.5))

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

## ---- eval=FALSE---------------------------------------------------------
#  pcfln = pcf( Kinhom(subset(ln, marks=="T_cells")) )
#  
#  plot(pcfln, lty = 1, log = "x")

## ---- fig.width=4.5, fig.height=4.5, echo=FALSE--------------------------
par(mar = c(4.5, 4.5, 0.5, 0.5))

data(pcfln, package = "BioC2015Oles")
plot(pcfln, lty = 1, log = "x", main = "")