I have recently started learning how to use Mathematica's brilliant image processing tools, and the image I've taken as first example is the following (source):

enter image description here

So far I've tried to first use LocalAdaptiveBinarize on the image and then use the MorphologicalGraph for the graph mapping but the results appear quite off since the resulting graph has about $40000$ vertices, whereas we have about $310$ particles in the image. The ideal mapping would be to map each particle to a vertex (and edges between particles in contact) and study the structure of the configuration as a graph.

s2 = MaxDetect@
  LocalAdaptiveBinarize[img, 8, PerformanceGoal -> "Quality"]
gvertex = MorphologicalGraph[s2, VertexCoordinates -> Automatic]

Binarized version:

enter image description here

Trying without the Binarization yields somewhat better results but the resulting graph still has little to do with the image.

  • Is there a way to process the image such that the particles can be more accurately detected? In other words, how should one process such particle based images (where typically like here the particles can be assumed to be spheres) in order to detect the particles positions before invoking MorphologicalGraph?

  • Finally, given we perform the graph mapping, how do assess how close the mapping has been? In other words, other than the basic checks of looking at vertex counts, how can we draw a close comparison between the result and the original image?

Brief update after the wonderful answers:

To compare the two resulting graphs that are obtained by both methods of user LukasLang and NikiEstner, the number of assigned vertices (i.e. detected particles) is $188$ and $273$ respectively, and the degree distributions are shown below (in same order):

enter image description here

enter image description here

I reckon these differences arise from the fact that the starting points were different: in the first mentioned answer a binarized version of the original image was used which meant partially information about the depth of the particles in the image was lost. Generally speaking, it is not immediately clear how to assess whether in determining the neighbourhood of a particle, how the particle's position depth (brightness variation, as neatly demonstrated by LukasLang) we see in the image should be taken into account.

  • $\begingroup$ I understand that the question is about creating a graph from the image. But to decide among the potential approaches, don't you need to know the "true" graph and the "true" number of particles? At least for a single image and preferably for several images. $\endgroup$
    – JimB
    Aug 13, 2019 at 21:49
  • $\begingroup$ @JimB Indeed, perfectly valid point! unfortunately I haven't been given much else to work on. Though I managed to find another example (added) for which the average core size is known, which can be used to assess the accuracy of the detection of either methods. But an example for which the graph properties are properly known and mentioned I cannot seem to find so far... $\endgroup$
    – user52181
    Aug 14, 2019 at 12:21

2 Answers 2


Here is one approach. See the section at the bottom about a few comments on how I chose the most important image processing parameters.

We start with your binarized image:

img = Import["https://i.stack.imgur.com/GAghg.png"]

enter image description here

The basic idea is to use the fact that the borders between particles seem to be nicely separated from the partciles themselves.

Next, we use MorphologicalComponents and SelectComponents to get the background:

bgImg = SelectComponents[MorphologicalComponents[ColorNegate[img], 0.99], Large] //
   Unitize //
  Colorize[#1, ColorRules -> {1 -> White}] &

enter image description here

Next, some cleaning:

procImg = bgImg //
    Dilation[#, 2] & //
   Closing[#, DiskMatrix@6] & //

enter image description here

Now we can apply MorphologicalComponents to get the individual particles, and then we use ArrayFilter with Max to grow them together (Update: I have updated the filter function to only apply Max if the center cell is 0 - this ensures that the individual regions can only grow into the empty space. Additionally, I'm using Nest to apply a filter with a smaller radius multiple times - this should help with growing all particles equally):

comps = procImg //
     ImagePad[#, -2] & //
    MorphologicalComponents[#, 0.5, CornerNeighbors -> False] & //
       If[#[[3, 3]] == 0, Max@#, #[[3, 3]]] &,
       ] &,
     ] &;

enter image description here

The last step is to use ComponentMeasurements with "Neighbours" (to decide which edges to include) and "Centroid" (to position the vertices) to build the graph:

ComponentMeasurements[comps, {"Neighbors", "Centroid"}, "PropertyComponentAssociation"] //
   DeleteDuplicates[Sort /@ Join @@ Thread /@ KeyValueMap[UndirectedEdge]@#Neighbors],
   VertexCoordinates -> Normal@#Centroid,
   VertexSize -> 0.7,
   VertexStyle -> Yellow,
   EdgeStyle -> Directive[Yellow, Thick],
   PlotRange -> Transpose@{{0, 0}, ImageDimensions@img},
   Prolog -> Inset[ImageMultiply[img, 0.7], Automatic, Automatic, Scaled@1]
   ] &

enter image description here

Choosing the parameters

A few notes on how I chose the parameters: The are three key parameters in the process above: The radius for Dilation and Closing, and the nesting parameter used for ArrayFilter. In the following, I will briefly discuss each step. (You will notice that most parameters are not too critical, so making them a bit bigger might help to make the process more robust)


The goal in this step is to make sure the individual particles are cleanly enclosed by the background. We do this by applying Dilation with an appropriate radius. The following shows the effect of a few different values - essentially, as long as the tiny gaps are closed, the parameter is fine.

Row@Table[bgImg // Dilation[#, i] &, {i, 0, 3}]

enter image description here


This step is to remove small gaps in the background that are not real particles. The bigger the radius of the DiskMatrix, the more holes are closed.

Row@Table[bgImg // Dilation[#, 2] & // Closing[#, DiskMatrix@i] &, {i, 2, 8, 2}]

enter image description here


This step is to grow the individual particles together, in order to decide which ones are adjacent. We do this by repeatedly (using Nest) applying Max based ArrayFilter. The more often we apply the filter an the bigger the radius of the filter, the more the particles can be separated and still considered adjacent.

Row@Table[procImg //
      ImagePad[#, -2] & //
     MorphologicalComponents[#, 0.5, CornerNeighbors -> False] & //
    With[{n = i},
       If[#[[n + 1, n + 1]] == 0, Max@#, #[[n + 1, n + 1]]] &,
      ] & // Colorize, {i, 1, 13, 4}]

enter image description here

Note: I chose to use multiple applications of a smaller filter instead of one big one to make sure that all particles are grown more or less equally. Otherwise, the Max part will always choose the particle with the biggest index to grow.

Estimating the z-coordinate of the particles

We can try to estimate the z-position of the particles by looking at the brightness of the particles in the individual image. To do this, we supply the raw image to ComponentMeasurements together with the labeling mask (comps), which allows us to use Mean to get the average brightness of each particle.

rawImg = Import["https://i.stack.imgur.com/rUnvs.jpg"];

     ImageResize[rawImg, ImageDimensions@img],(* make the image the same size *)
     "GrayScale" (* convert to 1-channel image *)
  {"Neighbors", "Centroid", "Mean", "Area"},
  ] //
   Table[Property[i, VertexSize -> Sqrt[#Area[i]/250]], {i, 
     Length@#Neighbors}] (* use the area for the size *),
   DeleteDuplicates[Sort /@ Join @@ Thread /@ KeyValueMap[UndirectedEdge]@#Neighbors],
   VertexCoordinates -> (* use the mean brightness as z-coordinate *)
    Normal@Merge[Apply@Append]@{#Centroid, 500 #Mean},
   EdgeStyle -> Directive[Blue, Thick],
   PlotRange -> Append[All]@Transpose@{{0, 0}, ImageDimensions@img}
   ] &

enter image description here

  • $\begingroup$ very nice answer! can we convert binarized image to mesh in MMA? thanks a lot in advance! $\endgroup$
    Aug 9, 2019 at 20:26
  • $\begingroup$ @ABCDEMMM what exactly do you mean? Which image? And how should the mesh be constructed (i.e. what should the cells/edges correspond to)? $\endgroup$
    – Lukas Lang
    Aug 9, 2019 at 20:40
  • $\begingroup$ e.g. the third picture from your test, and a nice example: static.cambridge.org/resource/id/… $\endgroup$
    Aug 9, 2019 at 21:34
  • 1
    $\begingroup$ Really impressive, what a transformation! thx for including the intermediate steps. If I may ask 1-2 follow-up questions just to understand better: i) during the cleaning, how do you roughly go about choosing those numerical parameters, $2$ for Dilation and $6$ for DiskMatrix? To know how one should adjust those for a different image. ii) For the graph mapping part, if I understood correctly, we model the particles as centroids and but how are edges decided? a distance criterion as in geometric graph? Finally, iii) how in our detection we dealt with fact that particles had different sizes? $\endgroup$
    – user52181
    Aug 10, 2019 at 0:10
  • $\begingroup$ @user929304 i) I've added a section explaining this to the answer, hope that helps. ii) I use ComponentMeasurements on the grown-together particles. This provides us both with the centroid (to position the particles), as well as the neighboring particles (to get the edges). You can look at the documentation for other potentially useful metrics. The advantage of using neighbors instead of a distance threshold is that particle sizes matter less, and we properly respect gaps, which might otherwise lead to extraneous edges. $\endgroup$
    – Lukas Lang
    Aug 10, 2019 at 13:28

@user929304 asked me for a way to solve this question that's not based on his binarization. After playing with the image a little bit, this is the simplest solution I came up with.

The idea is that between the particles, there's a thin dark "ridge" that can be detected with RidgeDetect:

img = Import["https://i.stack.imgur.com/rUnvs.jpg"]    
ridges = RidgeFilter[-img, 5];

(the 5 is an estimate of how thick the dark "ridge" is - but the code isn't very sensitive. I get more or less the same result for filter sizes 2..10.)

I then use a distance transform to get the distance to the nearest ridge for each point:

distRidges = 

and the maxima in this distance image are the centers of the particles we're trying to detect:

distMax = MaxDetect[distRidges, 5];

(5 is the minimum radius of a particle. Again, I get similar results for a range of 2..10.)

and WatershedComponents can find components from these centers (I've written an explanation WatershedComponents of here)

morph = WatershedComponents[ridges, distMax, Method -> "Basins"];

ComponentMeasurements will then find connected components and neighbors for each component:

comp = ComponentMeasurements[{img, morph}, {"Centroid", "Neighbors"}];

in the form

{1 -> {{18.3603, 940.324}, {21, 32}}, 2 -> {{140.395, 943.418}, {16, 21, 24}}, 3 -> {{286.265, 931.95}, {4, 16, 18, 26}}}...

so comp /. (s_ -> {c_, n_}) :> {s -> # & /@ Select[n, # > s &]}] will turn this into a list of graph edges:

graph = Show[img,
  Graph[comp[[All, 1]], 
   Flatten[comp /. (s_ -> {c_, n_}) :> {s -> # & /@ 
        Select[n, # > s &]}], VertexCoordinates -> comp[[All, 2, 1]], 
   EdgeStyle -> Directive[{Red, Thick, Opacity[1]}]]]

and EdgeDetect can be used to find component edges:

edges = Dilation[EdgeDetect[Image[morph], 1, .001], 2];
edgeOverlay = 
 Show[img, SetAlphaChannel[ColorReplace[edges, White -> Red], edges]]

the result then looks like this:

enter image description here

Add: (Response to comment)

does your method differ in how it tackles the fact that the particles in the image are stacked in 3D? Or are we assuming all the particles' centroid to be in the same plane (i.e. purely treated as 2D)? E.g. in the center top, there's a very bright particle which means it is standing on top of the lower stack, does that matter in the above scheme for finding its connected neighbourhood?

If we look at the area you mentioned in 3d, it looks like this:

enter image description here

trim = ImageTrim[img, {{755, 800}}, 150];
Row[{Image[trim, ImageSize -> 400], 
  ListPlot3D[ImageData[trim][[;; , ;; , 1]], PlotTheme -> "ZMesh", 
   ColorFunction -> "SunsetColors", ImageSize -> 500]}]

Now the particles don't have clear "peaks" in the center. That's why looking for local maxima in the brightness image directly doesn't work very well. But they do have "canyons" between them. That's what RidgeDetect looks for. It doesn't assume that the particles are "in the same plane", it just assumes that there's a thin "canyon" between adjacent particles that's "lower" (darker) than both of them


I wanted to ask you about understanding how ComponentMeasurements is in fact finding the neighbours of the particles

The interesting stuff happens in WatershedComponents, not ComponentMeasurements. Imagine the result of RidgeFilter as a 3d landscape:

enter image description here

Now imagine it starts raining on this 3d landscape. Or, alternatively, that someone starts pouring water into each of these valleys. At first, you will have separate pools of water. As the water rises, pools will meet at certain lines. These lines are called watersheds. The components enclosed by these watersheds are the components found by WatershedComponents and then measured by ComponentMeasurements. So the components that share a watershed, where two pools "meet" as the waterlevel rises, are neighbors in the neighborhood graph.

  • $\begingroup$ Dear Niki, I hope all is well! I was wondering if I could ask (admittedly very late) a follow-up question: to find the neighbours of each detected particle, we are using ComponentMeasurements[{img, morph}, {"Centroid", "Neighbors"}], which connects each to all its neighbours. But is there a way to extend this neighbourhood definition such that we can introduce a distance threshold? i.e., centroids/particles more apart than a distance x cannot be connected. Last, if I may, how did you create the colored image after the watershed step? Many thanks in advance! $\endgroup$
    – user52181
    Mar 25, 2020 at 17:58
  • $\begingroup$ @user929304: The color image of the components was just the output Colorize[...]. Finding neighbors with a distance threshold is a bit more tricky. You could use ComponentMeasurements to get the Mask of each component, then use Dilation to enlarge each mask (separately) and use e.g. DeleteDuplicates@Flatten[morph*dilatedMask] to pick out the component indices in the enlarged mask, for each component $\endgroup$ Mar 26, 2020 at 6:26
  • $\begingroup$ Hello! Thank you for these great explanations. I have recently been discussing a similar problem here where I use the image example shown in this post in order to figure out how to compute the particle density from the images. Both the binned approach (counter per bin) and the SmoothDensityHistogram seem to be quite inaccurate for this. I was wondering if this is simply due to how I'm detecting the particles. If you know of a more accurate way of estimating the density, it would definitely be useful to have your input on that post. Thanks $\endgroup$ Aug 14, 2020 at 12:34

Your Answer

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