18
$\begingroup$

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.

$\endgroup$
  • $\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 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$ – user929304 Aug 14 at 12:21
21
$\begingroup$

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] & //
  ColorNegate

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] & //
   Nest[
     ArrayFilter[
       If[#[[3, 3]] == 0, Max@#, #[[3, 3]]] &,
       #,
       2
       ] &,
     #,
     2
     ] &;
Colorize@comps

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"] //
 Graph[
   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)

Dilation:

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

Closing:

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

ArrayFilter:

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},
      ArrayFilter[
       If[#[[n + 1, n + 1]] == 0, Max@#, #[[n + 1, n + 1]]] &,
       #,
       n
       ]
      ] & // 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"];

ComponentMeasurements[
  {
   ImagePad[
    ColorConvert[
     ImageResize[rawImg, ImageDimensions@img],(* make the image the same size *)
     "GrayScale" (* convert to 1-channel image *)
     ], 
    -2
    ],
   comps
   },
  {"Neighbors", "Centroid", "Mean", "Area"},
  "PropertyComponentAssociation"
  ] //
 Graph3D[
   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

$\endgroup$
  • $\begingroup$ very nice answer! can we convert binarized image to mesh in MMA? thanks a lot in advance! $\endgroup$ – ABCDEMMM Aug 9 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 at 20:40
  • $\begingroup$ e.g. the third picture from your test, and a nice example: static.cambridge.org/resource/id/… $\endgroup$ – ABCDEMMM Aug 9 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$ – user929304 Aug 10 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 at 13:28
13
$\begingroup$

@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 = 
  DistanceTransform@ColorNegate@MorphologicalBinarize[ridges];

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


ADD 2

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.

$\endgroup$
  • $\begingroup$ Many thanks for taking the time and coming up with this new approach, really beautiful & interesting! If I may ask: compared to the old binarization I and Lukas had adopted, 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? Thanks again! $\endgroup$ – user929304 Aug 13 at 18:38
  • $\begingroup$ Dear Niki, many many thanks for the 2nd addition, beautifully explained! If I may ask 2 follow-up questions: i) From graph point of view, do the adjacencies obtained from the watershed correspond to nearest neighbours or contact neighbours (particles overlapping)? or possibly another definition of connectivity? In other words, for the obtained graph, how do we describe what connects two nodes, i.e. what kind of connectivity does the watershed algorithmic approach yield? $\endgroup$ – user929304 Aug 23 at 12:34
  • $\begingroup$ ii) we can see that some of the completely dark regions have also been assigned a node (see e.g. near top left region), although they might be really empty space between particles, so could we somehow tune how deep of a particle we still assign a node to? or on the safe-side, in the limit of not including any of those background particles, we'd have the very upper fully visible and bright layer. $\endgroup$ – user929304 Aug 23 at 12:36
  • $\begingroup$ Hi, sorry for all the follow-up questions, I really appreciate your patience and throughout getting back to me on them, really kind! Regarding the landscape made out of the ridges (to demonstrate the watershed approach), how do you determine the height of the barrier walls (corresponding to ridge lines separating particles)? Would the found neighbourhood be sensitive to the chosen heights, i.e. lead to different pools of water connectivity? Thanks again. $\endgroup$ – user929304 Aug 28 at 0:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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