Extraction of statistics of cracks (craquelure in paintings)

Question

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.

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:

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:

Using @kglr's code I get:

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?

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.stack.imgur.com/xbqhB.png"];

g = MorphologicalGraph@ColorNegate[img]

Needs["IGraphM`"];

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

IGSmoothen[g] // 
 IGVertexMap[vertexAssoc[GraphEmbedding][g], 
  VertexCoordinates -> VertexList]

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

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

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

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

---

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.