21
$\begingroup$

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

enter image description here

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:

enter image description here

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!

$\endgroup$
2
  • $\begingroup$ So, if you combine nikie's and yohbs's answers, you get the entire deal. $\endgroup$ Jul 13, 2015 at 12:45
  • $\begingroup$ @Guesswhoitis ...which means I'll have a tough time deciding to whom I'll award the bounty. Perhaps I'll just award two of them. $\endgroup$
    – march
    Jul 13, 2015 at 17:28

2 Answers 2

20
+100
$\begingroup$

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]

enter image description here

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

bin = MorphologicalBinarize[ImageSubtract[img, whiteCombined]]

enter image description here

(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

enter image description here

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]

enter image description here

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

enter image description here

$\endgroup$
10
+50
$\begingroup$

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:

enter image description here

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]

enter image description here

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

enter image description here

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;
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.