Measuring fractal dimension of natural objects from digital images

Question

This is a useful topic. A college physics lab, medical diagnostics, urban growth, etc. - there is a lot of applications. On this site by Paul Bourke about Google Earth fractals we can get a high resolution images (in this post they are low res - import from source for experiments). For example, around Lake Nasser in Egypt:

img = Import["http://paulbourke.net/fractals/googleearth/egypt2.jpg"]

The simplest method I know is Box Counting Method which has a lot of shortcomings. We start from extracting the boundary - which is the fractal object:

{Binarize[img], iEdge = EdgeDetect[Binarize[img]]}

Now we could partition image into boxes and see how many boxes have at least 1 white pixel. This is a very rudimentary implementation:

MinS = Floor[Min[ImageDimensions[iEdge]]/2];
data = ParallelTable[{1/size, Total[Sign /@ (Total[#, 2] & /@ (ImageData /@ 
         Flatten[ImagePartition[iEdge, size]]))]}, {size, 10, MinS/2, 10}];

From this the slope is 1.69415 which is a fractal dimension that makes sense

line = Fit[Log[data], {1, x}, x]

13.0276 + 1.69415 x

Plot[line, {x, -6, -2}, Epilog -> Point[Log[data]], 
  PlotStyle -> Red, Frame -> True, Axes -> False]

Benchmark: if I run this on high res of Koch snowflake i get something like ~ 1.3 with more exact number being 4/log 3 ≈ 1.26186.

Question: can we improve or go beyond the above box counting method?

All approaches are acceptable if they find fractal dimension from any image of natural fractal.

Answer 1

You can still use box count, but doing it smarter :)

Counting boxes with at least 1 white pixel from ImagePartition can be done more efficiently using Integral Image, a technique used by Viola-Jones (2004) in their now popular face recognition framework. For a mathematical motivation (and proof), Viola and Jones point to this source.

Actually, someone already asked about a Mathematica implementation here.

What Integral Image allows you to do is to compute efficiently the total mass of any rectangle in an image. So, you can define the following:

IntegralImage[d_] := Map[Accumulate, d, {0, 1}];
data = ImageData[ColorConvert[img, "Grayscale"]]; (* img: your snowflake image *)
ii = IntegralImage[data];

Then, the mass (white content) of a region, is

(* PixelCount: total mass in region delimited by two corner points, 
   given ii, the IntegralImage *)
PixelCount[ii_, {p0x_, p0y_}, {p1x_, p1y_}] := 
   ii[[p1x, p1y]] + ii[[p0x, p0y]] - ii[[p1x, p0y]] - ii[[p0x, p1y]];

So, instead of partitioning the image using ImagePartition, you can get a list of all the boxes of a given size by

PartitionBoxes[{rows_, cols_}, size_] := 
   Transpose /@ Tuples[{Partition[Range[1, rows, size], 2, 1], 
                        Partition[Range[1, cols, size], 2, 1]}];

If you apply PixelCount to above, as in your algorithm, you should have the same data but calculated faster.

PixelCountsAtSize[{rows_, cols_}, ii_, size_] :=
    ((PixelCount [ii, #1, #2] &) @@ # &) /@ PartitionBoxes[{rows, cols}, size];

Following your approach here, we should then do

fractalDimensionData = 
    Table[{1/size, Total[Sign /@ PixelCountsAtSize[Dimensions[ii], ii, size]]}, 
          {size, 3, Floor[Min[Dimensions[ii]]/10]}];
line = Fit[Log[fractalDimensionData], {1, x}, x]

Out:= 10.4414 + 1.27104 x

which is very close to the actual fractal dimension of the snowflake (which I used as input).

Two things to note. Because this is faster, I dared to generate the table at box size 3. Also, unlike ImagePartition, my partition boxes are all of the same size and therefore, it excludes uneven boxes at the edges. So, instead of doing minSize/2 as you did, I put minSize/10 - excluding bigger and misleading values for big boxes.

Hope this helps.

Update

Just ran the algorithm starting with 2 and got this 10.4371 + 1.27008 x. And starting with 1 is 10.4332 + 1.26919 x, much better. Of course, it takes longer but still under or around 1 min for your snowflake image.

Update 2

And finally, for your image from Google Earth (eqypt2.jpg) the output is (starting at 1-pixel boxes)

12.1578 + 1.47597 x

It ran in 43.5 secs in my laptop. Using ParallelTable is faster: around 28 secs.

Answer 2

"Can we improve the box-counting method?" Most certainly.

But first of all, it is important to point out that if you are only checking if a box is empty or not you are effectively measuring the $D_0$ (the capacity or in other words the box-counting dimension). By referring to $D_0$ as the fractal dimension you are assuming that your object is a monofractal, i.e. that all its generalized dimensions $D_q$ (ref.1) $$ D_q=\lim_{\epsilon \rightarrow \infty}\frac{\frac{1}{1-q}\log{\sum_i^Np_i^q}}{\log{\frac{1}{\epsilon}}} $$

are equal; $D_0$=$D_1$(entropy dimension)=$D_2$(correlation dimension) and so on...

Normally you should refrain from such assumptions unless you know the generating process. However, even if you assume that the object is monofractal, review literature of fractal dimension estimation concludes that the box-counting method is a bad choice. As Thelier puts it in "box-counting algorithm is a poor way to estimate the box-counting dimension" (ref. 2).

So what can we do? The basic principle in fractal dimension estimation is to study the object at different scales. The question is how to sample this continuum of scales. There are two different class of methods. The first one is the fixed-sized methods. These count the number of points in a fixed size region (boxes/spheres). They perform poorly because you never know how to grow your boxes (should I double, triple?). The second school is the fixed-mass methods. These grow their extend by reaching over to their mth nearest neighbor. Thus they always assure that they are covering the same mass (hence the name). even though fixed-mass methods are more robust than fixed-size methods, both suffer from finite-size and edge effect. Thus you need some extra measures to avoid the edges of your object, otherwise at large scales (of mass/size) you start to get this leveling off effect (that you also have in your last figure).

A few years ago we developed a complicated method (non-overlapping barycentric fixed-mass method) to do exactly that (ref. 3). Here is a summary of the fractal dimension estimations with the BFM method:

I first gave it a go with your Koch snowflake and I get all $D_q$s to be more or less equal to $1.25$ as you would expect from a monofractal.

Then I run the Egypt2.jpg file and as you can see all $D_q$s are not exactly equal. This can be used to argue that the coastline might be multifractal (even if slightly). Another thing to notice is that the slope is not the same throughout the whole mass range. If you consider only the range $[10 - 631]$ you would get a different $D_q$ curve. This implies that the object is not self-similar throughout the entire scale. This could also also explains why the previous answer finds $D_0$ something around $1.47$ while you find $1.69$. In other words there is a scaling break.

References:

1. Renyi, A. (1970). Probability Theory. Amsterdam: North-Holland.
2. Theiler, J. (1990). Estimating fractal dimension. Journal of the Optical Society of America A, 7(6), 1055.
3. Kamer, Y., Ouillon, G., & Sornette, D. (2013). Barycentric fixed-mass method for multifractal analysis. Physical Review E, 88(2), 022922.