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"]

Lake Nasser boundary

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]]}

outline of Lake Nasser

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]

plot of fractal dimension line

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.

  • 4
    $\begingroup$ You have a lot of programs in mathematica to measure fractal dimensions and multi fractal spectrum of an image in the book: Fractal Geography, Andre Dauphine, Wiley, 2012 See the book on the Wolfram Mathematica | Books or Amazon ![enter image description here](wolfram.com/books/profile.cgi?id=8108) $\endgroup$
    – Dauphine
    Oct 16, 2012 at 9:32
  • $\begingroup$ Vitaliy, with all due respect, the ending bullets make the question (actually there isn't a question mark anywhere, right?) rather wide ranging. Perhaps you could focus it somewhat and make it sound more like a real question? $\endgroup$ Oct 16, 2012 at 12:59
  • 1
    $\begingroup$ @SjoerdC.deVries Updated, it is a single question now. $\endgroup$ Oct 16, 2012 at 15:06
  • 1
    $\begingroup$ could you perhaps describe or link some examples of use of fractal dimension? $\endgroup$
    – magma
    Oct 17, 2012 at 10:27
  • $\begingroup$ @magma: Wikipedia lists some examples here: en.wikipedia.org/wiki/… $\endgroup$
    – shrx
    Jun 14, 2013 at 20:22

2 Answers 2


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.


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.

  • $\begingroup$ There isn't a "snowflake" image. All three images come from the same object (and I guess the fractal dimension should be the same for all of them) $\endgroup$ Dec 23, 2013 at 21:03
  • $\begingroup$ @belisarius, the PO published a Koch snowflake image as Benchmark (see last part of post) and that's the one I used. I called it 'snowflake' for short. On the other hand, this is still a numerical procedure, so the numbers will be different depending on approximation. The Update2 refers to the eqypt2.jpg image, also made available by the PO, as he asked for an improvement suitable for real-life images. $\endgroup$
    – carlosayam
    Dec 23, 2013 at 22:40
  • $\begingroup$ sorry, I missed that link. Thanks $\endgroup$ Dec 24, 2013 at 0:05
  • $\begingroup$ +1 This is very neat @caya, thank you! I will wait "a bit" hoping that someone can implement things like wavelet multi-fractals or similar. $\endgroup$ Feb 14, 2014 at 20:02
  • $\begingroup$ @VitaliyKaurov, glad you liked it. Note that Viola & Jones rightly pointed out a link between integral image and Haar wavelet basis in their paper; although this wasn't explored further. It is unclear to me the link between the paper you mentioned and fractal dimension, but certainly worth reading as I am interested in these kind of problems. Cheers. $\endgroup$
    – carlosayam
    Feb 16, 2014 at 12:40

"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. enter image description here

Then I run the Egypt2.jpg file and as you can see all $D_q$s are not exactly equal. enter image description here 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.


  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.

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