Message Boards Message Boards

Measuring the fractal dimension of a tree photo

GROUPS:

This is a photo of some tree branches that I took at the park today.

It looks like a fractal, and a nice candidate for measuring the fractal dimension. Here's how we can do that in Mathematica.

Let us start by binarizing the image.

img = Import["~/Downloads/20170304_121322.jpg"];

bin = LocalAdaptiveBinarize[ImageAdjust[img, 0.3], 200, {0.9, 0, 0}];

For a better result, I increased the contrast slightly (ImageAdjust), then used local adaptive binarization. This method chooses a different binarization threshold for each pixel based on its neighbourhood. It helps preserve small branches while avoiding including any part of the clouds (the weather was not very good today).

The easiest way to measure fractal dimension is box counting: overlay a square grid, and see how many of the grid cells contain some of the object in the image. What would we get if the object is a line?

Table[
  Rasterize[Graphics[{Antialiasing -> False, Line[{{0, 0}, {1, 1}}]}],
    ImageSize -> 256, RasterSize -> k],
  {k, 2^Range[2, 8]}
  ] // ListAnimate

The box count would double (i.e. multiply by $2^1$) every time the box size is halved. What would happen if the object is a filled disk?

Table[
  Rasterize[Graphics[{Antialiasing -> False, Disk[]}], 
   ImageSize -> 256, RasterSize -> k],
  {k, 2^Range[2, 8]}
  ] // ListAnimate

The box count would quadruple (i.e. multiply by $2^2$) every time the box size is halved.

In general, for a $d$-dimensional object the box count $n$ is proportional to power $-d$ of the box size $l$:   $n \sim l^{-d}$. For some objetcs, $d$ can be a fractional (non-integer) number. Such objects are called fractals. Let's try it on the tree photo!

First, let us choose box sizes in pixels so that an integer number of boxes fit in the image.

seq = Reverse[Intersection @@ Divisors /@ ImageDimensions[img]]
(* {1008, 504, 336, 252, 168, 144, 126, 112, 84, 72, 63, 56, 48, 42, 36, 28, 24, 21, 18, 16, 14, 12, 9, 8, 7, 6, 4, 3, 2, 1} *)

Luckily, the image width and height had many small prime factors in common, so we have many box sizes to work with. Now let us count boxes.

{width, height} = ImageDimensions[img]
(* {3024, 4032} *)

result = {
  width/#,
  Last@First@ImageLevels@Image[
      BlockMap[Min, ImageData[bin], {#, #}],
      "Bit"]
  } & /@ seq; // AbsoluteTiming

(* {9.79212, Null} *)

BlockMap will partition the image data into blocks (i.e. boxes) and take the minimum of each. Since black pixels are represented by 0 and white pixels by 1, we get 0 only for those boxes that contain at least one black pixel, i.e. they contain a bit of the tree branches. BlockMap was introduced in Mathematica 10.2. Those using older versions can substitute BlockMap[f, matrix, {k,k}] by Map[f, Partition[matrix, {k,k}], {2}].

ImageLevels counts how many times each possible pixel value appears in an image. To make sure that the only possible pixel values are 0 and 1, we create a "Bit" image.

Is the result really a power function of the form $n = \text{(const.)} l^d$? The easy way to test this is to plot it on a log-log plot. Power functions appear as a line when plotted on a logarithmic scale.

ListLogLogPlot[result]

enter image description here

The slope of the line will give the exponent: $\ln n = \ln\text{(const.)} + d \ln l$.

fm = LinearModelFit[Log[result], {1, x}, x]

It is about $\approx 1.85$. Despite the thick tree trunks, the image of the branches behaves as a lower-than-two dimensional structure.

Let us plot the power law fit together with the original data:

plot = Show[
  Plot[fm[x], {x, Log[result[[1, 1]]], Log[result[[-1, 1]]]}, 
   PlotStyle -> Black],
  ListLogLogPlot[result, PlotStyle -> Red],
  AspectRatio -> 1, Axes -> False, Frame -> True
  ]

enter image description here

You can see that the power law is a pretty good fit. The slope of the curve is more or less the same at all size scales. This indicates a certain kind of scale invariance in the image. Indeed, it is clear that if we magnified a small part of the image, it would look similar to the whole: small branches are miniature copies of larger ones. This is a characteristic feature of fractals.

To visualize the box counting method, we can show the series of finer and finer box-grids:

sz = 252;
imgs = Image[BlockMap[Min, ImageData[bin], {#, #}], "Bit"] & /@ seq

Scaling them to the same size and animating them show how they approach the original image better and better:

imgs = If[First@ImageDimensions[#] > sz, 
     ImageResize[Image[#, Real], sz], (* downscale as grayscale *)
     ImageResize[#, sz, Resampling -> "Nearest"]] & /@ imgs; (* upscale as bitmap *)

frm = Table[
  ImageCompose[ImageResize[img, 252], {im, 0.5}], {im, imgs}];

ListAnimate[frm]

enter image description here

We can show the animation alongside the scaling plot:

frames = Table[
   Row[{
     Show[plot, AspectRatio -> 1, ImageSize -> 336, 
      Epilog -> {AbsolutePointSize[10], Red, Point@Log[result[[i]]]}],
      Image[frm[[i]], Magnification -> 1]
     }],
   {i, Length[imgs]}];

With some extra styling it looks like this:

enter image description here

I hope you enjoyed this small demonstration!

POSTED BY: Szabolcs Horvát
Answer
8 months ago

Nice! Thanks for sharing. I was already wondering if you had used BlockMap or Map[...Partition[...]], but well, you thought of both.

POSTED BY: Sander Huisman
Answer
8 months ago

Actually, I used a compiled version of Map + Partition when @Yode in the StackExchange chatroom reminded me of BlockMap. Strangely, their timing scales differently with problem size, but I haven't investigated in detail yet.

POSTED BY: Szabolcs Horvát
Answer
8 months ago

Felicitation, very pedagogic André Dauphiné

POSTED BY: André Dauphiné
Answer
8 months ago

enter image description here - Congratulations! This post is now Staff Pick! Thank you for your wonderful contributions. Please, keep them coming!

POSTED BY: Moderation Team
Answer
8 months ago

This is so cool.

POSTED BY: Eduardo Serna
Answer
8 months ago

Thanks for sharing this sweet post, Szabolcs. It looks like you manage to snap a picture where the parallax is not messing up too much the computation.

Does it make sense to compare structures with similar fractal dimensions? What objects would we find with a dimensions around 1.85 like these trees?

POSTED BY: Matthias Odisio
Answer
8 months ago

This crop of the Amazon river gives 1.96.

enter image description here enter image description here

POSTED BY: Matthias Odisio
Answer
8 months ago

Thanks for the post Matthias! This is a beautiful example of a self-similar structure.

Unfortunately I only had time for a very quick test right now, but if we use more box sizes (a proper way would be to pad the image when necessary), we get a clearer picture about what is going on.

This river network is so dense that if we look at it from far enough, then it appears to fully cover the ground. It looks two-dimensional. But as we zoom in, the self-similarity becomes clear and the fractal nature emerges. Here's the scaling plot I got:

enter image description here

The horizontal axis is proportional to the number of boxes that the image was partitioned into. Thus a large number corresponds to a small box size.

There is a clear break in the middle of the curve. The slope of the curve is ~2 to the left of that break (i.e. for large boxes), which indicates that at this scale the object looks two-dimensional. The slope is ~1.5 to the right of that break (i.e. small boxes), which shows that it looks like a fractal at that scale. The power law is valid separately on both sides of the break.

To plot this, I was sloppy and simply cropped the image to {2580, 1810} to have many common divisors. I think a better way would be to pad the image with white pixels when it cannot be partitioned into an integer number of boxes. This way you can use an arbitrary box size.


If we could zoom in without limit, the exponent would become 1 when the rivers don't branch out anymore. They look like lines at this point. If we could zoom even more, the exponent would become 2 when the width of the rivers becomes visible.

POSTED BY: Szabolcs Horvát
Answer
8 months ago

Group Abstract Group Abstract