# How to segment nuclei from a noisy image

I am trying to segment and count nuclei in a noisy 3-D image stack and it is proving to be somewhat challenging. The strategy I could come up with requires removal of noise, generating seeds (for watershed) and employing watershed segmentation.

(*importing the desired images*)
imgseries = ImageAdjust /@
Import["C:\\Users\\Ali \Hashmi\\Desktop\\emb1_MMStack_Pos0.ome.tif"][[6 ;; -34 ;; 2]];

(* now I subtract background from each 2D slice and fill holes in image  *)
img3dinit = Image3D[ Module[{img = #, background},
background = ImageTake[img, {1, 88}, {1, 204}];
FillingTransform[img - Mean[background]]
] & /@ imgseries];

(* convolving the image with a GaussianMatrix and subtracting extraneous noise *)
img3dconv = ImageConvolve[img3dinit, GaussianMatrix[{{2, 2, 2}}]]
background = ImageTake[img3dconv, All, {1, 58}, {1, 360}];
img3d = DeleteSmallComponents[(img3dconv - Mean[background])];


at this point, I have a stack that looks like the following:

(*binarizing the imaging, deleting small components and filling holes *)
img3dbin = FillingTransform[DeleteSmallComponents[MorphologicalBinarize[img3d, 0.18], 500]];


img3dbin

(* calculating distance transform *)
img3dbin, Padding -> 0], 500];


distance

(* determing seeds from distance transform *)
markers = MaxDetect[distance, 0.1];


markers

(* performing watershed and setting background to zero *)
matr = WatershedComponents[GradientFilter[img3d, 2], markers,Method -> "Rainfall"];
maxind = Last@*First@ComponentMeasurements[matr, "LabelCount"];
backgroundind = Keys /@ MaximalBy[ComponentMeasurements[matr, "Area"], Values, 1];

(* visualizing the segmentation *)
(seg = ArrayComponents[matr, maxind, Thread[backgroundind -> 0]]) // Colorize


below is how the segmentation looks:

having a quick look at the final segmentation I feel it can be improved. Of course I can merge smaller blobs with the nearest neighbours provided the small ones are less than a certain size. All such trivial matters can be taken care of.

However, my main concern is to generate nice seeds and get a more accurate segmentation from the watershed method starting from the noisy stack. Do you consider the use of TotalVariationFilter / WienerFilter or some other form of pre-processing prior to performing DistanceTransform for generating seeds? Furthermore, which Method is best for WatershedComponents function? Any help will be much appreciated.

• The images at that link appear to be all black. And I was unable to import them. Do you have code that will correctly do that? – Daniel Lichtblau Nov 5 '17 at 15:48
• @DanielLichtblau if you click on the link then on the top right you should be able to download the file. By mapping ImageAdjust (as in the code) you should be able to see nuclei. – Ali Hashmi Nov 5 '17 at 18:25
• I did that and when I attempted to open it I got a black square. I am guessing it is a set of tiff files rather than just one, and the .tif extension is making the opener believe it is a single image. But that's a wild guess. Some supporting evidence is that it is 512x512 and purports to occupy around 31 MB. That would be about right for the 60 or so images that are said to be at your link. – Daniel Lichtblau Nov 5 '17 at 18:48
• Update: I can successfully Import a set of 62 images from that link into Mathematica. I can run the first several lines of code above. What I get appears to be darker but maybe it is good enough. Not sure I will make progress on the actual problem but at least I can give it a try. – Daniel Lichtblau Nov 5 '17 at 18:56
• @DanielLichtblau here is the link with the imgseries dropbox.com/sh/41sd3el5zkxvjxx/AAAFg0erau_ySwSr5xGKeddpa?dl=0 – Ali Hashmi Nov 5 '17 at 19:14

Here's my approach. As with Daniel, it's not a complete answer. I got to a point where I needed to define criteria for selecting components which is an arbitrary exercise for me since I don't know what I'm looking at. Anyhow, in the code below, I stored the raw (w/o ImageAdjust) tiff files in the symbol data. First, I want to make sure that all images are being adjusted using the same brightness and contrast. I used a Manipulate to scroll through the data and tweak the brightness/contrast to what I thought was reasonable. Something along the lines of:

Manipulate[ImageAdjust[data[[i]], {c, b}], {{c, 0.45}, -10, 10},
{{b, 68}, 0, 1000}, {{i, 3}, 2, Length@data, 2}]


Then I applied a gradient filter with the 'optimal' brightness/contrast.

img3d = GradientFilter[ImageAdjust[#, {.45, 36}] & /@
data[[2 ;; -1 ;; 2]] // Image3D, {4, 2}] // ImageAdjust


Some tweaking needs to be done, but I think this approach looks promising. Now for MorphologicalComponents:

img3dbin = MorphologicalBinarize[img3d, 0.05];
s = MorphologicalComponents[img3dbin];
Colorize@s


Changing the parameters of the GradientFilter will improve the unique number of components found. With these settings, I find 200, and you can do something like ComponentMeasurements[s, "Circularity"] to get some information about the components (note, I don't know how the properties are calculated for 3D images). Once you know this, you can then pull out the components that meet your criteria and compare them to the original data:

o = SelectComponents[s, "FilledCircularity", 0.8 < # < 1.2 &] // Colorize
ImageAdd[o, img3d] // Image3DSlices


You can see that I've missed a number of objects, but I think you can pick up the process from here.

This is not a complete answer but it might help. I downloaded the 12 tiff images. Somehow I get 3 dimensions (color channels maybe, but ImageDimensions does not show them). I mention this to explain why I use Mean[Mean[...]] in a couple of places. Modify as needed.

(* Import the images *)

images0 =
FileSystemMap[Import[#, "TIFF"] &,

(* Process with the original code but modified to use Mean[Mean[...]] *)

images1 = Map[ImageAdjust, Values[images0]];

img3dinit =
Image3D[Map[
Module[{img = #, background},
background = ImageTake[img, {1, 88}, {1, 204}];
FillingTransform[img - Mean[Mean[background]]]] &, images1]];
img3dconv = ImageConvolve[img3dinit, GaussianMatrix[{{2, 2, 2}}]];
background = ImageTake[img3dconv, All, {1, 58}, {1, 360}];
img3d = DeleteSmallComponents[(img3dconv - Mean[Mean[background]])];

img3dbin =
FillingTransform[
DeleteSmallComponents[MorphologicalBinarize[img3d, 0.18], 500]];

distance =
FillingTransform@
DeleteSmallComponents[


Here I simply sum the image data arrays, make a 2D image, and binarize it.

pixels = ImageData[distance];
pixelSum = Total@pixels;
image2D = Binarize[Image[pixelSum]];
pixels2D = ImageData[image2D];


Here is how it looks, colorized by components.

Colorize[MorphologicalComponents[image2D]]


Now extract the new image data, separate zeros and ones. Make a distance function for the zeros (black background pixels). We will create a new 2D image, modified from this one by making all pixels near the background black. We are in effect shrinking the cells. This has the effect of separating some of the overlapping ones.

pixels2D = ImageData[image2D];
pZero = Position[pixels2D, 0];
pOne = Position[pixels2D, 1];
nf = Nearest[pZero -> "Distance"];
newDistances = Map[Min, nf[pOne]];
newPixels2D = pixels2D;


I chose to make all white pixels within distance 3 of the background into black pixels.

Do[If[newDistances[[j]] <= 3,
newPixels2D[[Apply[Sequence, pOne[[j]]]]] = 0], {j, Length[pOne]}];
newImage2D = Image[newPixels2D];


If we use a value larger than 3 we start to lose cells (we may already lose a few at this level). So you might want to experiment a bit with retaining "small" components, perhaps by not altering any that have a max diameter below some set threshold.

Here is how this one looks.

Colorize[MorphologicalComponents[newImage2D]]


Obviously this is far from a full separation of cells. It might be useful as a start though. From here perhaps the idea of not altering small components can be applied, to do further separation only on the larger ones.