library("EBImage")

## Parameters used in the segmentation
p = list(
  nuc.athresh.filter = makeBrush(35, shape='box')/(35*35),
  nuc.athresh.t = 0.01,
  nuc.morpho.se = makeBrush(3, shape='diamond'),
  nuc.min.density = 0.15,
  nuc.min.size = 150,
  nuc.max.size = 1000,
  
  cell.thresh.filter = matrix(c(0,1,0,1,2,1,0,1,0)/6,nc=3,nr=3),
  cell.thresh.t = 0.12,
  cell.morpho.se = makeBrush(3, shape='diamond'),
  cell.propagate.mix.power = 0.2,
  cell.propagate.lambda = 0.0001,
  cell.min.density = 0.1,
  cell.min.size = 150,
  cell.max.size = 18000
  )

## Read files
a = readImage('001-02-C03-actin.tif')
t = readImage('001-02-C03-tubulin.tif')
h = readImage('001-02-C03-dna.tif')

## Adjust dynamic range and background level of images.
## Here we use some 'magic numbers', in reality, would 
## determine these in a data-driven way.
a = 2.82*a - 0.17
t = 5.03*t - 0.35
h = 3.00*h - 0.15

## display images
display(a)
display(t)
display(h)

## make composite image 
cal = rgbImage(r=a, g=t, b=h)
display(cal)

######  find the nuclei

## adaptive thresholding
nmask = h > filter2(h, p$nuc.athresh.filter) + p$nuc.athresh.t
display(nmask)

## morphological opening to remove objects smaller than the structuring element
nmask = opening(nmask, p$nuc.morpho.se)
display(nmask)

## fill holes in nuclei
nmask = fillHull(nmask)
display(nmask)

## label each connected set of pixels with a distinct ID
nseg = bwlabel(nmask)

## object sizes
barplot(table(nseg)[-1])

## remove objects that are too small or too large
mom = moments(nseg, h)
irm = which(mom[,'m.pxs']<p$nuc.min.size | mom[,'m.pxs']>p$nuc.max.size)
print(irm)
nseg = rmObjects(nseg, irm)

## display segmented objects with random colors from the RGB cube
##  --  little helper function that does this using index gymnastics --
displayWithRandomColors = function(x) {
  colorIndices = as.vector(x)+1
  n = max(colorIndices)
  rc = function() c(1, runif(n-1))[colorIndices]
  z = Image( c(rc(), rc(), rc()), dim = c(dim(nseg), 3), colormode='Color')
  display(z)
}
displayWithRandomColors(nseg)
  

###### estimate where the boundaries between cells are

## compute composite channel (Euclidean sum)
mix = (a^2+t^2+h^2)^0.5
display(mix)

## compute cell body mask on composite channel
cmask = filter2(mix, p$cell.thresh.filter)>=p$cell.thresh.t
display(cmask)

## clean cell body mask
cmask = closing(cmask, p$cell.morpho.se)
display(cmask)

## display composite view of nucleus mask and cell body mask
display(rgbImage(g=xor(cmask, nmask), b=nmask))

## remove nuclei which are outside of cells
nseg[cmask==0] = 0

## Voronoi segmentation on the composite channel
cseg = propagate(mix^p$cell.propagate.mix.power, nseg, cmask, p$cell.propagate.lambda)

## display the segmented objects (cells) with random colors
displayWithRandomColors(cseg)

## remove objects (cells) that are too small or too large
mom = moments(cseg, h)
irm = which(mom[,'m.pxs']<p$cell.min.size | mom[,'m.pxs']>p$cell.max.size)
print(irm)
cseg = rmObjects(cseg, irm)
nseg = rmObjects(nseg, irm)

####### done !
## display original data together with the segmentations
z = paintObjects(nseg, cal, col=c('#1010ff'))
z = paintObjects(cseg, z, col=c('#ffffff'))
display(z)

## size distribution of cells
##  (see the 'moments' function for many further per-cell
##   statistics that can be computed)
hist( table(cseg)[-1], col='skyblue', breaks=20)