How to segment nuclei from a noisy image

Question

(data can be downloaded from this link.)

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 *)
distance = FillingTransform@DeleteSmallComponents[ImageAdjust@DistanceTransform[
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.

Answer 1

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.

Answer 2

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"] &, 
   "C:\\Users\\Daniel\\Downloads\\nuclei"];

(* 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[
        ImageAdjust@DistanceTransform[img3dbin, Padding -> 0], 500];

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.