###################################################
### chunk number 1: loadlib
###################################################
library("EBImage")


###################################################
### chunk number 3: reada
###################################################
a = readImage('images/001-01-C03-actin.jpeg')
display(a)


###################################################
### chunk number 5: readd
###################################################
d = readImage('images/001-01-C03-dna.jpeg')
display(d)


###################################################
### chunk number 7: matrix
###################################################
display(a*2)
display(a+d)
display(d+0.3)


###################################################
### chunk number 8: transf
###################################################
display(resize(a, w=256))
display(rotate(a, 30))


###################################################
### chunk number 9: rgbImage
###################################################
t = readImage('images/001-01-C03-tubulin.jpeg')
x = rgbImage(red=a, green=t, blue=d)
display(x)

###################################################
### chunk number 11: writeImage
###################################################
writeImage(x, 'images/001-01-C03-rgb.jpeg')


###################################################
### chunk number 12: batch-rgbImage
###################################################
ids = unique(substr(dir('images'), 1, 10))
for (id in ids) {
  cat(id, '\n')
  a = readImage(file.path('images', paste(id, '-actin.jpeg', sep='')))
  t = readImage(file.path('images', paste(id, '-tubulin.jpeg', sep='')))
  d = readImage(file.path('images', paste(id, '-dna.jpeg', sep='')))
  x = rgbImage(red=a, green=t, blue=d)
  writeImage(x, file.path('images', paste(id, '-rgb.jpeg', sep='')))
}

###################################################
### chunk number 14: batch-resize
###################################################
dir.create('thumbs')
img = dir('images')
for (i in img) {
  cat(i, '\n')
  y = readImage(file.path('images', i)) 
  writeImage(resize(y, w=256), file.path('thumbs', i))
}

###################################################
### chunk number 16: hwriter
###################################################
library(hwriter)
m1 = matrix(file.path('thumbs', img), nrow=length(ids), byrow=TRUE)
m2 = matrix(file.path('images', img), nrow=length(ids), byrow=TRUE)
m1 = hwriteImage(m1, link=m2, table=FALSE)
rownames(m1) = ids
hwrite(m1, 'gallery.html')
browseURL('gallery.html')


###################################################
### chunk number 17: globthreshold
###################################################
display(d)
nmask1 = d>0.15
display(nmask1)


###################################################
### chunk number 18: adapt
###################################################
nmask2 = thresh(d, 15, 15, 0.005)
display(nmask2)
display(rgbImage(red=nmask1, green=nmask2))


###################################################
### chunk number 19: clean
###################################################
mk3 = makeBrush(3, shape='diamond')
nmask2 = opening(nmask2, mk3)
nmask2 = fillHull(nmask2)
display(nmask2)


###################################################
### chunk number 20: nseg
###################################################
nseg = bwlabel(nmask2)
nn = max(nseg)
display(nseg/nn)
print(nn)


###################################################
### chunk number 21: nf
###################################################
nf = hullFeatures(nseg)
print(nf)


###################################################
### chunk number 22: remove small nuc
###################################################
nr = which(nf[,'g.s'] < 150)
nseg = rmObjects(nseg, nr)
nn = max(nseg)
print(nn)
display(nseg/nn)


###################################################
### chunk number 24: cmask
###################################################
mix = sqrt(a^2+t^2+d^2)
display(mix)
cmask = filter2(mix, mk3/5)>=0.12
cmask = closing(cmask, mk3)
nseg[!(cmask>0)] = 0
display(rgbImage(red=cmask, green=nseg))
display(x)


###################################################
### chunk number 26: cseg
###################################################
cseg = propagate(mix^0.2, nseg, cmask, 0.1, 0)
cseg = fillHull(cseg)
display(cseg/max(cseg))


###################################################
### chunk number 27: cf
###################################################
cf = hullFeatures(cseg)
cr = which(cf[,'g.s'] < 150)
cseg = rmObjects(cseg, cr)
nseg = rmObjects(nseg, cr)
acf = list()
acf[[id]] = cf


###################################################
### chunk number 28: seg
###################################################
print(max(cseg))
seg = paintObjects(nseg, x, col=c('#ffff00'))
seg = paintObjects(cseg, seg, col=c('#ff00ff'))
display(seg)
writeImage(seg, file.path('images', paste(id, '-seg.jpeg', sep='')))


###################################################
### chunk number 30: batch seg
###################################################
for (id in ids) {
  cat(id, '\n')
  a = readImage(file.path('images', paste(id, '-actin.jpeg', sep='')))
  t = readImage(file.path('images', paste(id, '-tubulin.jpeg', sep='')))
  d = readImage(file.path('images', paste(id, '-dna.jpeg', sep='')))
  x = rgbImage(red=a, green=t, blue=d)
  nmask2 = thresh(d, 15, 15, 0.005)
  
  mk3 = makeBrush(3, shape='diamond')
  nmask2 = opening(nmask2, mk3)
  nmask2 = fillHull(nmask2)
  nseg = bwlabel(nmask2)
  nf = hullFeatures(nseg)
  nr = which(nf[,'g.s'] < 150)
  nseg = rmObjects(nseg, nr)
  mix = (a^2+t^2+d^2)^0.5
  cmask = filter2(mix, mk3/5)>=0.12
  cmask = closing(cmask, mk3)
  nseg[!(cmask>0)] = 0
  cseg = propagate(mix^0.2, nseg, cmask, 0.1, 0)
  cseg = fillHull(cseg)
  cf = hullFeatures(cseg)
  cr = which(cf[,'g.s'] < 150)
  cseg = rmObjects(cseg, cr)
  nseg = rmObjects(nseg, cr)
  
  ## done
  seg = paintObjects(nseg, x, col=c('#ffff00'))
  seg = paintObjects(cseg, seg, col=c('#ff00ff'))
  writeImage(seg, file.path('images', paste(id, '-seg.jpeg', sep='')))
  
  acf[[id]] = cf
}


###################################################
### chunk number 31: seg results
###################################################
img = dir('images')
img = img[sort(c(grep('seg', img), grep('rgb', img)))]
i2 = file.path('images', img)
m2 = matrix(i2, nrow=length(ids), byrow=TRUE)
m2 = hwriteImage(m2, table=FALSE)
rownames(m2) = ids
hwrite(m2, 'gallery2.html')
browseURL('gallery2.html')


###################################################
### chunk number 32: multidensity
###################################################
library(geneplotter)
cp = sapply(acf, function(ac) ac[,'g.p'])
multidensity(cp[c('001-01-C03', '001-01-D03', '001-01-M04', '003-01-H03')], xlim=c(0,1000), bw=100,
             xlab='Cell perimeter (in pixels)', main='Distribution of cell perimeters')


###################################################
### chunk number 34: wilcox.test
###################################################
wilcox.test(cp[['001-01-C03']],cp[['003-01-H03']])
wilcox.test(cp[['001-01-C03']],cp[['001-01-M04']])