Get measurements of touching particles in an image

Question

I have a picture of glass spheres. I want to calculate different properties of these spheres such as circularity, size, and count.

I have managed to fill the white centers of the particles to simplify the particle analysis. The problem however is that the particles are touching each other. I managed to end up with an image looking like:

I did this by Erosion to get the particles separated, detection of particles by MorphologicalComponents, and finally by increasing the particle size back to the original size using Dilation:

particles1 = Erosion[nn, DiskMatrix[30]]
  // MorphologicalComponents
  // Colorize
  // DeleteBorderComponents
  // Dilation[#, DiskMatrix[30]] &

My problem is now that I get the particles correctly colored, but somehow Mathematica does not recognize them as separate particles anymore after the Dilation command!

I would like to use ComponentMeasurements to get my desired particle parameters, but the command does not work anymore after the Dilation command, because the particles touch again. I am sure there must be a simple solution to my problem, but can't figure it out. Must it not be possible to use the already obtained color information of the particles??

Thankx a lot in advance!

Answer 1

As I've written in a previous answer, I would start by searching for the white centers:

img = Import["https://i.stack.imgur.com/LaMAg.jpg"];
img = ColorConvert[img, "Grayscale"];

components = 
  ComponentMeasurements[
   MorphologicalBinarize[Closing[img, DiskMatrix[5]]], {"Centroid", 
    "Circularity", "EquivalentDiskRadius", "Mask"}, #2 > 0.65 &];

This is the trickiest part, because the thresholds for Circularity (and possibly other thresholds) have to be adjusted so all centers are found. (Sidenote: For this image, it was relatively easy to find a Circularity-Threshold that found all and only the white circles. If this is more difficult for your other images, this might be a perfect task for Classify, Mathematica's new machine learning functionality. But that would be a topic for a separate question...)

The "Mask" measurement gives a binary mask for each component. If you sum those masks, you get the union of all white centers:

whiteCombined = Image[Total[components[[All, 2, -1]]]];
HighlightImage[img, whiteCombined]

you can then easily fill the white centers by subtracting them from the original image:

bin = MorphologicalBinarize[ImageSubtract[img, whiteCombined]]

(Possible alternative: FillingTransform also takes a markers parameter that you could use: bin=ColorNegate@FillingTransform[ColorNegate[Binarize[img]], whiteCombined] - same result for this image, could be different for others.)

From here, it's standard DistanceTransform + WatershedComponents, where the last argument to WatershedComponents explicitly specifies the basin centers - i.e. WatershedComponents should create one component for each white circle we found in the first step. If you don't do this, WatershedComponents will split some of the elongated components.

centers = components[[All, 2, 1]];
ImageData[ColorNegate@bin]*
  WatershedComponents[ColorNegate@DistanceTransform[ColorNegate@bin], 
   centers] // Colorize

Full code (where img is the source image):

components = 
  ComponentMeasurements[
   MorphologicalBinarize[Closing[img, DiskMatrix[5]]], {"Centroid", 
    "Circularity", "EquivalentDiskRadius", "Mask"}, #2 > 0.65 &];
whiteCombined = Image[Total[components[[All, 2, -1]]]];
bin = MorphologicalBinarize[ImageSubtract[img, whiteCombined]];
centers = components[[All, 2, 1]];
segmentation = 
  ImageData[ColorNegate@bin]*
   WatershedComponents[ColorNegate@DistanceTransform[ColorNegate@bin],
     centers];
Colorize[segmentation]

ADD: Maybe I should have made this clearer: unlike the segmentation in the question, this segmentation gives separate regions, without overlap. No post-processing is required, aside from removing border components and a few isolated pixels:

components = 
  ComponentMeasurements[{DeleteBorderComponents@segmentation, 
    img}, {"Centroid", "EquivalentDiskRadius", "Area", "SemiAxes", 
    "Orientation"}, #2 > 10 &];

ListPlot[Sort[components[[All, 2, 3]]], PlotRange -> All, 
 Filling -> 0]

Show[img,
 Graphics[
  {
   Red,
   Thick,
   Rotate[Circle[#[[1]], #[[4]]], #[[5]]] & /@ components[[All, 2]],
   Text @@@ components[[All, 2, {3, 1}]]
   }]]

Answer 2

First of all, when I run your code to produce particles1 I do not get the same result that you show, but maybe this is a version difference. I'll use this image as the starting point:

Now, the problem, as you correctly identified, is that you have some overlap between particles. However, the overlap is small, and you can easily identify which colors you have "a lot of" (i.e. they correspond to a particle) and which colors you have "little of" (i.e. they correspond to an overlap between colors). To do this, first convert your image to bits

imageBits = ImageData[img, "Byte"];

imageBits is a matrix, where each element is of the form {r,g,b,a} representing the color in RGB-alpha. Now, the next command flattens the matrix to a list, deletes all black pixels, and then sorts them by how many pixels you have of each color:

particleColors=Reverse@SortBy[Tally[DeleteCases[Flatten[imageBits,1],{0,0,0,255}]],Last];

Now you have a list of the form

{{{185, 91, 67, 255},  6311}, 
 {{32, 155, 198, 255}, 6287}, 
 {{189, 229, 135, 255}, 6078},...}

You see that you have 6311 pixels of the color {185, 91, 67, 255}, 6287 pixels of the color {{32, 155, 198, 255}, 6287} and so on. If you plot the number of pixels, you'll see that you have 53 or 54 particles:

ListLinePlot[particleColors[[1 ;; 60, 2]], PlotRange -> All, PlotMarkers -> Automatic]

So now, you can create your label matrix, by replacing each of the particle colors by a positive integer, and replacing all other colors by 0:

labelMatrix = imageBits /. Table[particleColors[[i,1]] -> i, {i, 54}];
labelMatrix = labelMatrix /. {i_Integer, j_Integer, k_Integer, 255} -> 0

The resulting label matrix is:

labelMatrix // MatrixPlot

and you can run, for example, ComponentMeasurements[labelMatrix, "Area"] to get

{1 -> 2329.63, 2 -> 2514.38, 3 -> 2561.88, 4 -> 2594.75, 5 -> 2635.13,
  6 -> 2694.5, 7 -> 2791.25, 8 -> 3100.25, 9 -> 3328.63, 
 10 -> 3330.13, 11 -> 3347.63, 12 -> 3447., 13 -> 3522.25, 
 14 -> 3579., 15 -> 3585., 16 -> 3609.5, 17 -> 3702.63, 18 -> 3718., 
 19 -> 3737.88, 20 -> 3870.63, 21 -> 3891.13, 22 -> 4045.13, 
 23 -> 4051.13, 24 -> 4096.5, 25 -> 4226.38, 26 -> 4529.88, 
 27 -> 4550.13, 28 -> 4566.38, 29 -> 4583.25, 30 -> 4597.25, 
 31 -> 4604.75, 32 -> 4641.25, 33 -> 4655., 34 -> 4768.75, 
 35 -> 4782.25, 36 -> 4783.5, 37 -> 4860.25, 38 -> 4947.25, 
 39 -> 4968.25, 40 -> 4992.13, 41 -> 5034.63, 42 -> 5111.13, 
 43 -> 5114.25, 44 -> 5161.38, 45 -> 5208.63, 46 -> 5266.63, 
 47 -> 5283., 48 -> 5338.88, 49 -> 5523.88, 50 -> 5825.38, 
 51 -> 5967.38, 52 -> 6091.5, 53 -> 6300.38, 54 -> 6323.88}

For convenience, here's only the relevant code without the explanations:

imageBits = ImageData[img, "Byte"];
particleColors = 
  Reverse@SortBy[
    Tally[DeleteCases[Flatten[imageBits, 1], {0, 0, 0, 255}]],
   Last]
labelMatrix = imageBits /. Table[particleColors[[i, 1]] -> i, {i, 55}];
labelMatrix = labelMatrix /. {i_Integer, j_Integer, k_Integer, 255} -> 0;