Computing a fractal dimension with Matlab: 1D, 2D and 3D Box-counting
F. Moisy, 9 july 2008.
Laboratory FAST, University Paris Sud.
Contents
About Fractals and box-counting
A set (e.g. an image) is called "fractal" if it displays self-similarity: it can be split into parts, each of which is (at least approximately) a reduced-size copy of the whole.
A possible characterisation of a fractal set is provided by the "box-counting" method: The number N of boxes of size R needed to cover a fractal set follows a power-law, N = N0 * R^(-DF), with DF<=D (D is the dimension of the space, usually D=1, 2, 3).
DF is known as the Minkowski-Bouligand dimension, or Kolmogorov capacity, or Kolmogorov dimension, or simply box-counting dimension.
To learn more about box-counting, fractals and fractal dimensions:
- http://en.wikipedia.org/wiki/Fractal
- http://en.wikipedia.org/wiki/Box_counting_dimension
- http://mathworld.wolfram.com/Fractal.html
- http://mathworld.wolfram.com/CapacityDimension.html
About the 'boxcount' package for Matlab
The following examples illustrate how to use the Matlab package 'boxcount' to compute the fractal dimension of 1D, 2D or 3D sets, using the 'box-counting' method.
The directory contains the main function 'boxcount', three sample images, and an additional function 'randcantor' to generate 1D, 2D and 3D generalized random Cantor sets.
Type 'help boxcount' or 'help randcantor' for more details.
Box-counting of a 2D image
Let's start with the image 'dla.gif', a 800x800 logical array (i.e., it contains only 0 and 1). It originates from a numerical simulation of a "Diffusion Limited Aggregation" process, in which particles move randomly until they hit a central seed. (see P. Bourke, http://local.wasp.uwa.edu.au/~pbourke/fractals/dla/ )
c = imread('dla.gif'); imagesc(~c) colormap gray axis square
Calling boxcount without output arguments simply displays N (the number of boxes needed to cover the set) as a function of R (the size of the boxes). If the set is a fractal, then a power-law N = N0 * R^(-DF) should appear, with DF the fractal dimension (Kolmogorov capacity).
boxcount(c)
The result of the box count can be obtained using:
[n, r] = boxcount(c) loglog(r, n,'bo-', r, (r/r(end)).^(-2), 'r--') xlabel('r') ylabel('n(r)') legend('actual box-count','space-filling box-count');
n =
Columns 1 through 7
44000 27466 11786 4265 1386 421 121
Columns 8 through 11
37 12 4 1
r =
Columns 1 through 7
1 2 4 8 16 32 64
Columns 8 through 11
128 256 512 1024
The red dotted line shows the scaling N® = R^-2 for comparision, expected for a space-filling 2D image. The discrepancy between the two curves indicates a possible fractal behaviour.
Local scaling exponent
If the set has some fractal properties over a limited range of box size R, this may be appreciated by plotting the local exponent, D = - d ln N / ln R. For this, use the option 'slope':
boxcount(c, 'slope')
Strictly speaking, the local exponent is not constant, but lies in the range [1.6 1.8].
Let's try now with another image, the so-called Apollonian gasket (Wikipedia, http://en.wikipedia.org/wiki/Image:Apollonian_gasket.gif ). The background level is 198 for this image, so this value is used to binarize the image:
c = imread('Apollonian_gasket.gif'); c = (c<198); imagesc(~c) colormap gray axis square figure boxcount(c) figure boxcount(c,'slope')
The local slope shows that the image is indeed approximately fractal, with a fractal dimension DF = 1.4 +/- 0.1 for scales R < 100.
Box-counting of a natural image.
Consider now this RGB (2272x1704) picture of a tree (J.A. Adam, http://epod.usra.edu/archive/images/fractal_tree.jpg ):
c = imread('fractal_tree.jpg'); image(c) axis image
Let's extract a rectangle in the blue (3rd) plane, and binarize the image for levels < 80 (white pixels are logical 'true'):
i = c(1:1200, 120:2150, 3); bi = (i<80); imagesc(bi) colormap gray axis image
[n,r] = boxcount(bi,'slope');
The boxcount shows that the local exponent is approximately constant for less than one decade, in the range 8 < R < 128 (the 'exact' value of Df depends on the threshold, 80 gray levels here):
df = -diff(log(n))./diff(log(r)); disp(['Fractal dimension, Df = ' num2str(mean(df(4:8))) ' +/- ' num2str(std(df(4:8)))]);
Fractal dimension, Df = 1.801 +/- 0.06394
Generalized random Cantor sets
Fractal sets may be obtained from an IFS (iterated function system). For example, the function 'randcantor' provided with the package generates a 1D, 2D or 3D generalized random Cantor set. This set is obtained by iteratively dividing an initial set filled with 1 into 2^D subsets, and setting each subset to 0 with probability P. The result is a fractal set (or "fractal dust") of dimension DF = D + log(P)/log(2) < D.
The following example generates a 2048x2048 image with probability P=0.8, i.e. fractal dimension DF = 1.678.
c = randcantor(0.8, 2048, 2); imagesc(~c) colormap gray axis image
Let's see its box-count and local exponent
boxcount(c)
figure
boxcount(c, 'slope')
For such set generated by an iterated scheme, the local slope shows as expected a well defined plateau, around DF = 1.7.
1D random Cantor set
1D random Cantor sets may also be generated. Here, a 2^18 = 262144 long set with P = 0.9 and expected fractal dimension DF = 1 + log(P)/log(2) = 0.848:
tic c = randcantor(0.9, 2^18, 1, 'show'); figure boxcount(c, 'slope'); toc
Elapsed time is 4.022846 seconds.
3D random Cantor set
Now a 3D random Cantor set of size (2^7)^3 = 128^3 with P = 0.7 and expected fractal dimension DF = 3 + log(P)/log(2) = 2.485. Note that 3D sets cannot be displayed using randcantor.
tic
c = randcantor(0.7, 2^7, 3);
toc
tic
boxcount(c, 'slope');
toc
Elapsed time is 12.197119 seconds. Elapsed time is 0.186419 seconds.
More?
That's all. To learn more about this package, type:
help boxcount.m
BOXCOUNT Box-Counting of a D-dimensional array (with D=1,2,3).
[N, R] = BOXCOUNT(C), where C is a D-dimensional array (with D=1,2,3),
counts the number N of D-dimensional boxes of size R needed to cover
the nonzero elements of C. The box sizes are powers of two, i.e.,
R = 1, 2, 4 ... 2^P, where P is the smallest integer such that
MAX(SIZE(C)) <= 2^P. If the sizes of C over each dimension are smaller
than 2^P, C is padded with zeros to size 2^P over each dimension (e.g.,
a 320-by-200 image is padded to 512-by-512). The output vectors N and R
are of size P+1. For a RGB color image (m-by-n-by-3 array), a summation
over the 3 RGB planes is done first.
The Box-counting method is useful to determine fractal properties of a
1D segment, a 2D image or a 3D array. If C is a fractal set, with
fractal dimension DF < D, then N scales as R^(-DF). DF is known as the
Minkowski-Bouligand dimension, or Kolmogorov capacity, or Kolmogorov
dimension, or simply box-counting dimension.
BOXCOUNT(C,'plot') also shows the log-log plot of N as a function of R
(if no output argument, this option is selected by default).
BOXCOUNT(C,'slope') also shows the semi-log plot of the local slope
DF = - dlnN/dlnR as a function of R. If DF is contant in a certain
range of R, then DF is the fractal dimension of the set C. The
derivative is computed as a 2nd order finite difference (see GRADIENT).
The execution time depends on the sizes of C. It is fastest for powers
of two over each dimension.
Examples:
% Plots the box-count of a vector containing randomly-distributed
% 0 and 1. This set is not fractal: one has N = R^-2 at large R,
% and N = cste at small R.
c = (rand(1,2048)<0.2);
boxcount(c);
% Plots the box-count and the fractal dimension of a 2D fractal set
% of size 512^2 (obtained by RANDCANTOR), with fractal dimension
% DF = 2 + log(P) / log(2) = 1.68 (with P=0.8).
c = randcantor(0.8, 512, 2);
boxcount(c);
figure, boxcount(c, 'slope');
F. Moisy
Revision: 2.10, Date: 2008/07/09