The fractal dimension of an object is a single scalar number that allows us to quantify how compact an object is , *i.e. *how wiggly a curve is, how wrinkled a surface is, how porous a complex volume is.

However estimating the fractal dimension of an object is not an easy task, and many methods exist. The simplest method is** **box counting: the idea is to fully cover the object with many boxes of a given size, count how many boxes are needed to cover the object and repeat for many box sizes. The scaling of the number of boxes covering the object with the size of the boxes gives an estimate for the fractal dimension of the object.

The algorithm has many limitation, but in its simplest form it can be easily implemented in Python. The idea is to simply bin the object in an histogram of variable bin sizes. This can be easily generalised to any dimensions, thanks to Numpy’s `histrogramdd ` function.

In the following code, I test the idea with a known fractal, Sierpinski’s triangle, which has an exact (Hausdorff) fractal dimension of log(3)/log(2).

import numpy as np import pylab as pl def rgb2gray(rgb): r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2] gray = 0.2989 * r + 0.5870 * g + 0.1140 * b return gray image=rgb2gray(pl.imread(&quot;Sierpinski.png&quot;)) # finding all the non-zero pixels pixels=[] for i in range(image.shape[0]): for j in range(image.shape[1]): if image[i,j]&gt;0: pixels.append((i,j)) Lx=image.shape[1] Ly=image.shape[0] print Lx, Ly pixels=pl.array(pixels) # pl.plot(pixels[:,1], pixels[:,0], '.', ms=0.01) # pl.show() print pixels.shape # computing the fractal dimension #considering only scales in a logarithmic list scales=np.logspace(1, 8, num=20, endpoint=False, base=2) Ns=[] # looping over several scales for scale in scales: print &quot;======= Scale :&quot;,scale # computing the histogram H, edges=np.histogramdd(pixels, bins=(np.arange(0,Lx,scale),np.arange(0,Ly,scale))) Ns.append(np.sum(H&gt;0)) # linear fit, polynomial of degree 1 coeffs=np.polyfit(np.log(scales), np.log(Ns), 1) pl.plot(np.log(scales),np.log(Ns), 'o', mfc='none') pl.plot(np.log(scales), np.polyval(coeffs,np.log(scales))) pl.xlabel('log $\epsilon$') pl.ylabel('log N') pl.savefig('sierpinski_dimension.pdf') print &quot;The Hausdorff dimension is&quot;, -coeffs[0] #the fractal dimension is the OPPOSITE of the fitting coefficient np.savetxt(&quot;scaling.txt&quot;, zip(scales,Ns))