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