Measuring fractal dimension of natural objects from digital images

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:

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[FDL]],
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.

-
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) – Dauphine Oct 16 '12 at 9:32
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? – Sjoerd C. de Vries Oct 16 '12 at 12:59
@SjoerdC.deVries Updated, it is a single question now. – Vitaliy Kaurov Oct 16 '12 at 15:06
could you perhaps describe or link some examples of use of fractal dimension? – magma Oct 17 '12 at 10:27
@magma: Wikipedia lists some examples here: en.wikipedia.org/wiki/… – shrx Jun 14 '13 at 20:22

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.

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.

-
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) – Dr. belisarius Dec 23 '13 at 21:03
@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. – carlosayam Dec 23 '13 at 22:40
sorry, I missed that link. Thanks – Dr. belisarius Dec 24 '13 at 0:05
+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. – Vitaliy Kaurov Feb 14 '14 at 20:02
@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. – carlosayam Feb 16 '14 at 12:40