The below image is of the pattern of cracks in the surface of Johannes Vermeer's Girl with a pearl earring, captured using a method called "grazing incidence illumination" (the light is shined at the painting from a very large angle). There has been some simple image processing to reveal just the cracks. You'll notice the region of her face has very different crack statistics (density, average length of crack, etc.) than outside her face.

Craquelure in "Girl with a pearl earring"

I'd like to develop a Mathematica tool to extract image data so I can compute statistics of such so-called "craquelure" in paintings such as this:

  • histogram of crack length (distances between crack intersections), i.e., lengths of each "path" in a graph
  • orientations of each crack ($0 \to \pi$)
  • overall average density of cracks (number of cracks per unit area)

and so forth. Calculating and plotting and histogramming such data is not the problem; extracting the raw data is.

The first step is to create a mathematical graph of the image data. Define the above image as craquelureImage, and then compute:

paintingGraph= MorphologicalGraph@
 ColorReplace[craquelureImage, {White -> Black, Black -> White}]

This gives a graph with several thousand vertices:

craquelure graph

I suspect there is a problem here, though, because when I highlight the vertices that have degree $2$, there are only a handful, despite the appearance of long "paths."

One can inspect the connectivity (edges), vertices, their location, etc., using common software functions, just as with any traditional graph.

But this graph is not quite appropriate for full analysis. It contains paths (linear sequences of edges) that are really just a single edge. Thus we need to perform what graph theorists call "path contraction"--replace each such linear path with a single edge.

In a separate question on this site, @kglr wrote very clever code that performs such path contraction, and retains the spatial location of the un-contracted vertices. That code works well on graphs generated by traditional code.

However, for some unknown reason it does not work on paintingGraph, defined above. It somehow scrambles the locations of vertices or mis-assigns edges. I have examined the VertexList, EdgeList, and so on of the underlying graph, and cannot determine why it does not work: I get the layout to be completely mixed up:

mixed-up graph layout

Using @kglr's code I get:

enter image description here

Hence my first problem:

Problem 1: How to make the path contraction code work on paintingGraph. I can only presume there is something special about a graph created by MorphologicalGraph that isn't obvious.

Assuming that problem is solved, and the spatial locations of the vertices are proper...

Problem 2: How does one extract from the graph the spatial lengths and orientations of each crack segment?

For reasons of material physics, cracks (in paintings, in dry mud, ...) will always meet at vertices of degree three. That is, they appear like Ts (at some orientation). Again for fundamental physical reasons involving sequential stress relief in the drying paint, the angles of the Ts are almost exactly $90^\circ$. Can we demonstrate that? We can easily extract vertices in the graph of degree $3$. Specifically, at these degree-three vertices....

Problem 3: How do we calculate the relative angle of the intersection at each T?

  • 1
    $\begingroup$ Jeepers, I'm still struggling with the last puzzle you set us ... $\endgroup$ Commented Mar 26, 2020 at 15:48
  • $\begingroup$ @HighPerformanceMark: Tell me about it... I've spent many hours on all of this!... and am deeply grateful for the help already offered. I assure you, that many art scholars will want this software tool, once I perfect it. $\endgroup$ Commented Mar 26, 2020 at 15:52
  • 2
    $\begingroup$ "I suspect there is a problem here, though, because when I highlight the vertices that have degree 2, there are only a handful" There is no problem. MorphologicalGraph isn't really meant to return degree-2 vertices. I expect these are artefacts. The vertices should all be either branching points or endpoints. $\endgroup$
    – Szabolcs
    Commented Mar 26, 2020 at 20:16
  • $\begingroup$ Re Problem 1: can you try defining vCoords = AssociationThread[VertexList[#], GraphEmbedding[#]] &; and changing VertexCoordinates -> {v_ :> GraphEmbedding[g][[v]]}to VertexCoordinates -> {v_ :> vCoords[g]@v}? $\endgroup$
    – kglr
    Commented Mar 26, 2020 at 20:51

1 Answer 1


David, you say that my solution to path contraction does not work on this graph either. However, I see no problem at all.

img = Import["https://i.sstatic.net/xbqhB.png"];

g = MorphologicalGraph@ColorNegate[img]

enter image description here


vertexAssoc[fun_][g_] := AssociationThread[VertexList[g], fun[g]]

IGSmoothen[g] // 
  VertexCoordinates -> VertexList]

enter image description here

Tested in M11.3, M12.0, M12.1.

Can you please show the code that you tried and that failed? The solution I gave to the path contraction is only two lines—there's not much room for variation.

Segment lengths and angles

The result of MorphologicalGraph has self-loops. First, we must get rid of them.

The following line requires Mathematica 12.0+ to preserve vertex coordinates and edge weights:

g = SimpleGraph[g];

If you have a Mathematica version earlier than 12.0, then do not use the above line. It would remove edge weights and vertex coordinates from g. Instead, use:

g = IGWeightedSimpleGraph[g, SelfLoops -> False, VertexCoordinates -> GraphEmbedding[g]];

Now we can extract the line segments and endpoint coordinates.

points = GraphEmbedding[g];
segments = IGIndexEdgeList[g];

A very fast way to visualize the segments is now

Graphics@GraphicsComplex[points, Line[segments]]

We can get the lengths and angles:

lengths = EuclideanDistance @@ points[[#]] & /@ segments;

angles = Mod[ArcTan @@ Subtract @@ points[[#]] & /@ segments, Pi];

Some of the segments are curved. lengths contains not their true length but the distance between their endpoints. The true lengths of the skeletonized image (the number of pixels in each segment) are stored in the edge weights and can be extracted as lengths2 = IGEdgeProp[EdgeWeight][g].

I believe MorphologicalGraph actually works with something like Thinning@Binarize@ColorNegate[img]. It counts pixels in a skeleton like this one. You may want to do this step yourself, to control the outcome better, then apply MorphologicalGraph to this skeleton (not to the original image). Pruning is also useful to get rid of very short segments, some of which may even be artefacts from Thinning.

Now you can do e.g. Histogram[angles], but first you may want to filter short segments. I recommend my BoolEval package for this. E.g., to show the angle histogram of segments longer than 3, use

Histogram[BoolPick[angles, lengths > 3], {0, Pi, 2 Pi/30}]

enter image description here

To find a good length threshold, you can use something like:

  GraphicsComplex[points, Line@BoolPick[segments, lengths > t]],
 {t, 0, 20}

enter image description here

To get better and more detailed results than what Mathematica is capable of, I recommend Fiji. The rough steps are: convert to 8-bit, threshold, skeletonize, analyse skeleton (but I'm not very fluent with Fiji). Unfortunately, I know of no convenient way to use Fiji from Mathematica. Fiji is very GUI-oriented.

  • $\begingroup$ The core line was: Fold[VertexContract, g, contractComponent[g, #] & /@ components2[g]]. But look at the lower-left in your solution: there appear to be several paths of length greater than 2. Why are they not contracted? I think (now) that the core problem is with image processing and MorphologicalGraph, which erroneously breaks paths. As you say, there are very few vertices of degree 2... which doesn't make sense. $\endgroup$ Commented Mar 26, 2020 at 20:12
  • $\begingroup$ @DavidG.Stork That is not my code. I understand that the code from someone else's solution may not work, but have you tried my simple 2-line solution (which is well-tested, as I've been using it for years, and faster than all others)? You said in your comment that my solution did not work either. $\endgroup$
    – Szabolcs
    Commented Mar 26, 2020 at 20:18
  • $\begingroup$ Ooops.. my mistake. I tried: Needs["IGraphM`"]; mg = MorphologicalGraph @ ColorNegate[VermeerCraquelure]; qq = IGSmoothen[Graph[mg]]; Graph[qq, GraphEmbedding[mg]] and got the rectangular "error." $\endgroup$ Commented Mar 26, 2020 at 20:23
  • 1
    $\begingroup$ @DavidG.Stork Are you using 11.3? Warning: SimpleGraph does not preserve vertex coordinates or edge weights in that version! $\endgroup$
    – Szabolcs
    Commented Mar 26, 2020 at 20:54
  • 3
    $\begingroup$ @DavidG.Stork Mathematica 12.0 brought a major change: now many graph functions preserve properties. E.g., if you take a subgraph, edge weights are no longer discarded. In earlier Mathematica, everything is discarded. For some reason, the developers refuse to document this change ... I asked many times during the 12.0 prerelease period, yet they never did it. I asked many times during the 12.1 prerelease period, still they never did it ... If you are using kglr's code and it didn't work as you expected, it may also be because you are using 11.3. $\endgroup$
    – Szabolcs
    Commented Mar 27, 2020 at 9:00

Your Answer

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

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