# Box counting method for fractal dimension

First I saw here that the box counting fractal dimension defined by

$D = \lim_{\epsilon \rightarrow 0}{ {\log N( \epsilon)} \over {\log { {1}\over{ \epsilon }}}}$

which makes sense for me. Then I saw here a mathematica code for calculating it. This method does a fitting of the data $\{ x= \log { {1}\over{ \epsilon }}, y= \log N( \epsilon) \}$ to a line, giving two parameter $a$ and $b$ such that

$y = a x + b$

and $a$ is interpreted as the box counting fractal dimension, which I don't understand because I don't think there is any reason to assume linear behaviour.

However, this method is "correct", because

${{y} \over {x}} = a + {{b}\over{x}} \quad$ and $\quad \lim_{\epsilon \rightarrow 0} = \lim_{x \rightarrow \infty}$

but by the same reason there are infinite other "correct" methods assuming

$y = a x + b + {{c} \over {x}} + {{d} \over {x^2}} + ...$

or other more complicated things.

Why do they choose a linear behaviour? (Instead of just plotting ${ {\log N( \epsilon)} \over {\log { {1}\over{ \epsilon }}}}$ against $\epsilon$ and taking the $y$-intercept)

• The limit is really only defined for mathematically idealized objects. For real objects, we should look at a sufficiently wide range of scales, from very zoomed out up to very zoomed in, and compute the dimension at each one. A shape is considered to be fractal only when the measured dimension stays approximately constant across multiple scales. The slope computation accomplishes just that. This youtube video describes the situation nicely. Nov 10 '17 at 15:39
• Precisely in that video it is noticed how slow the fractal dimension converges when calculated that way. The value $a$ in my question (the order of the polynomial in the video), only approaches the fractal dimension when $x \rightarrow 0$ (when the scaling factor goes to infinity in the video). So why not just take the computational or approximate limit? Nov 12 '17 at 9:48

Like Euclidean objects, fractals are idealized abstractions of reality. A tree is not a fractal any more than its trunk is a line segment. As a result, we cannot compute an actual limit to find the box-counting dimension of an object. Even an attempt to approximate the limit by using a small value of $\varepsilon$ in $$\frac{\log(N({\varepsilon}))}{\log(1/\varepsilon)}$$ is hazardous because we don't know a priori what scale of $\varepsilon$ works well with the object. The process of fitting a line to the data $$\{ x= \log { {1}/ \varepsilon_k}, y= \log N( \epsilon_k) \}$$ for some terms chosen from a sequence $\varepsilon_k$ that tends down to zero allows us to estimate the dimension in a way that accounts for behavior of the set over a large scale of values of $\varepsilon$.

We can illustrate this with an example which, since this question was originally posted on mathematica.stackexchange, we present with Mathematica. Consider the following image of the Sierpinski carpet:

Import[
"https://www.marksmath.org/FractalGeometryPackages/FractalGeometry/IteratedFunctionSystems.m"
];
A = {{1, 0}, {0, 1}}/3;
IFS = {A, #} & /@ {
{-1, 1}, {0, 1},  {1, 1},
{-1, 0},          {1, 0},
{-1, -1}, {0, -1},{1, -1}
};
sierpPic = ShowIFS[IFS, 5, PlotRange -> All,
Initiator -> Polygon[{{-3, -3}, {3, -3}, {3, 3}, {-3, 3}}/2],


Here's an easy way to box-count the image:

pic = Binarize[ImportString[ExportString[sierpPic, "PNG",
ImageSize -> 2^10]]];
data = ImageData[pic];
test[matrix_] := Boole[Length[Cases[matrix, 0, {2}, 1]] > 0];
count[n_] := Total[Flatten[Map[test, Partition[data, {n, n}], {2}]]];
radii = Table[2^k, {k, 0, 8}];

(* Out:
{626688, 167476, 45276, 12124, 3287, 864, 236, 60, 16}
*)


Here's the "approximate limit" computation:

Log[First[counts]]/Log[2^10] // N

(* Out: 1.92574 *)


Here's a slope fit computation:

Fit[Transpose[Log[{1/radii[[2 ;; -3]], counts[[2 ;; -3]]}]], {1, x}, x]

(* Out: 13.3465 + 1.89636 x *)


The actual dimension is $$\frac{\log(8)}{\log(3)} \approx 1.89279.$$ Thus, this choice of slope fitting does a better job.

Now, I freely admit that I fiddled a bit with the choice of points to use in my fit of the slope. The point, though, is that slope fitting can do a better job. Also, it's easy to see why the "approximate limit" approach might fail. If $\varepsilon$ is too small, the small blocks in the approximate Sierpinski carpet are solid and contribute too much to the box-count. Thus, we get too large a value. That's exactly the point when it comes to a real world object. An apparent, approximate self-similar structure is only valid down to some scale and you have no way of knowing how far down you should look.

Ultimately, though, box-counting with natural objects is necessarily approximate, error prone, and not necessarily even well defined.

If you prefer an example with Python, you find one over on Stackoverflow.