Skip to content

Instantly share code, notes, and snippets.

@rougier
Last active October 18, 2023 12:37
Show Gist options
  • Star 13 You must be signed in to star a gist
  • Fork 13 You must be signed in to fork a gist
  • Save rougier/e5eafc276a4e54f516ed5559df4242c0 to your computer and use it in GitHub Desktop.
Save rougier/e5eafc276a4e54f516ed5559df4242c0 to your computer and use it in GitHub Desktop.
Fractal dimension computing
# -----------------------------------------------------------------------------
# From https://en.wikipedia.org/wiki/Minkowski–Bouligand_dimension:
#
# In fractal geometry, the Minkowski–Bouligand dimension, also known as
# Minkowski dimension or box-counting dimension, is a way of determining the
# fractal dimension of a set S in a Euclidean space Rn, or more generally in a
# metric space (X, d).
# -----------------------------------------------------------------------------
import scipy.misc
import numpy as np
def fractal_dimension(Z, threshold=0.9):
# Only for 2d image
assert(len(Z.shape) == 2)
# From https://github.com/rougier/numpy-100 (#87)
def boxcount(Z, k):
S = np.add.reduceat(
np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),
np.arange(0, Z.shape[1], k), axis=1)
# We count non-empty (0) and non-full boxes (k*k)
return len(np.where((S > 0) & (S < k*k))[0])
# Transform Z into a binary array
Z = (Z < threshold)
# Minimal dimension of image
p = min(Z.shape)
# Greatest power of 2 less than or equal to p
n = 2**np.floor(np.log(p)/np.log(2))
# Extract the exponent
n = int(np.log(n)/np.log(2))
# Build successive box sizes (from 2**n down to 2**1)
sizes = 2**np.arange(n, 1, -1)
# Actual box counting with decreasing size
counts = []
for size in sizes:
counts.append(boxcount(Z, size))
# Fit the successive log(sizes) with log (counts)
coeffs = np.polyfit(np.log(sizes), np.log(counts), 1)
return -coeffs[0]
I = scipy.misc.imread("sierpinski.png")/256.0
print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(I))
print("Haussdorf dimension (theoretical): ", (np.log(3)/np.log(2)))
@rougier
Copy link
Author

rougier commented Mar 17, 2021

I'm not sure but the images seems to be a bit noisy which might impact the computation. Can you try to make black and white (for example I = I > 0.5 and check if the result is the same? If this is the case, this means we need some threshold value when deciding if a box is empty or not.

@jankaWIS
Copy link

jankaWIS commented Mar 17, 2021

Except for the first and last, they are all binary (ie B&W). I have tried and you were right, @rougier, about those, setting the threshold to 0.2 works on them.

What about this note?

Also, I believe that it should go to squares 2x2, that is:
sizes = 2**np.arange(n, 2, -1) should become sizes = 2**np.arange(n, 0, -1).

I am curious because you used to have it that way (at least that is what I see in the history) but you have a discrepancy in the note and in the code. Your code says:

2**np.arange(n, 1, -1) # eg for n=4  you get array([16,  8,  4])

but the comment says
Build successive box sizes (from 2**n down to 2**1) which suggest this code:

2**np.arange(n, 0, -1) # eg for n=4  you get array([16,  8,  4,  2])

And that I believe is the correct code.

Code

import numpy as np
import matplotlib.pyplot as plt

def fractal_dimension(Z, threshold=0.9):

    # Only for 2d image
    assert(len(Z.shape) == 2)

    # From https://github.com/rougier/numpy-100 (#87)
    def boxcount(Z, k):
        S = np.add.reduceat(
            np.add.reduceat(Z, np.arange(0, Z.shape[0], k), axis=0),
                               np.arange(0, Z.shape[1], k), axis=1)

        # We count non-empty (0) and non-full boxes (k*k)
        return len(np.where((S > 0) & (S < k*k))[0])


    # Transform Z into a binary array
    Z = (Z < threshold)
    plt.imshow(Z)
    plt.show()

    # Minimal dimension of image
    p = min(Z.shape)

    # Greatest power of 2 less than or equal to p
    n = 2**np.floor(np.log(p)/np.log(2))

    # Extract the exponent
    n = int(np.log(n)/np.log(2))

    # Build successive box sizes (from 2**n down to 2**1)
    sizes = 2**np.arange(n, 1, -1)

    # Actual box counting with decreasing size
    counts = []
    for size in sizes:
        counts.append(boxcount(Z, size))

    # Fit the successive log(sizes) with log (counts)
    coeffs = np.polyfit(np.log(sizes), np.log(counts), 1)
    return -coeffs[0]

special = [
    'https://upload.wikimedia.org/wikipedia/commons/thumb/7/7c/Peano_Sierpinski_carpet_4.svg/440px-Peano_Sierpinski_carpet_4.svg.png',
    'https://upload.wikimedia.org/wikipedia/commons/thumb/3/3b/Quadratic_Koch_2.svg/2880px-Quadratic_Koch_2.svg.png',
    'https://upload.wikimedia.org/wikipedia/commons/thumb/5/55/Quadratic_Koch.svg/300px-Quadratic_Koch.svg.png'
]
images = [
    'https://upload.wikimedia.org/wikipedia/commons/thumb/7/7c/Peano_Sierpinski_carpet_4.svg/440px-Peano_Sierpinski_carpet_4.svg.png',
    'https://upload.wikimedia.org/wikipedia/commons/thumb/b/b7/Sierpinski_carpet_3.svg/244px-Sierpinski_carpet_3.svg.png',
    'https://upload.wikimedia.org/wikipedia/commons/f/f0/Flocke.PNG',
    'https://upload.wikimedia.org/wikipedia/commons/thumb/3/3b/Quadratic_Koch_2.svg/2880px-Quadratic_Koch_2.svg.png',
    'https://upload.wikimedia.org/wikipedia/commons/thumb/5/55/Quadratic_Koch.svg/300px-Quadratic_Koch.svg.png'
]

dimensions = [(np.log(3)/np.log(2)), (np.log(3)/np.log(2)), (np.log(4)/np.log(3)), 1.46, 1.5]


for dim, im in zip(dimensions, images):
    # load im
    x = plt.imread(im)
    if im in special:
        print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(x[:,:,3], threshold=0.2))
        print("theoretical:        ", dim)
        print(f"min {x[:,:,3].min()} and max {x[:,:,3].max()} in the array/image")
        
    else:
        print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(x[:,:,1],threshold=0.2))
        print("theoretical:        ", dim)
        print(f"min {x[:,:,1].min()} and max {x[:,:,1].max()} in the array/image")

image

@jankaWIS
Copy link

I am getting this error message when executing this programme in PyCharm.
1)"AttributeError: module 'scipy.misc' has no attribute 'imread'"
2) "AttributeError: module 'scipy' has no attribute 'misc'"

@akshahi17 check this answer https://stackoverflow.com/questions/15345790/scipy-misc-module-has-no-attribute-imread
the function has been deprecated so you can either downgrade (not recommended) or use what @rougier recommended (imageio) or what I'm using, imread in matplotlib.pyplot.

@barbaroo
Copy link

Hi! I found this code very useful and I would like to cite it in my paper. Is there a correct way to cite this? Please help! I have to submit soon!

@rougier
Copy link
Author

rougier commented Jan 31, 2022

Glad the code is useful but I'm afraid there there's no proper way to cite a it. Maybe as a misc bibitem with the relevant url/author.

@antoni27
Copy link

antoni27 commented Apr 9, 2022

hello

sorry, can anyone tell me how to change the data that I want to get from the link to my laptop data?
I want to try to extract features of an iris

@jankaWIS
Copy link

jankaWIS commented Apr 17, 2022

Hi @antoni27, I am not really sure what data you have and how you want to use it but simply changing the path in line 52 of the original code to link to another image should be enough.

@antoni27
Copy link

antoni27 commented May 9, 2022

thanks, @jankaWIS
and sorry if I ask again is the 52 line is
images = [
'https://upload.wikimedia.org/wikipedia/commons/thumb/7/7c/Peano_Sierpinski_carpet_4.svg/440px-Peano_Sierpinski_carpet_4.svg.png',
'https://upload.wikimedia.org/wikipedia/commons/thumb/b/b7/Sierpinski_carpet_3.svg/244px-Sierpinski_carpet_3.svg.png',
'https://upload.wikimedia.org/wikipedia/commons/f/f0/Flocke.PNG',
'https://upload.wikimedia.org/wikipedia/commons/thumb/3/3b/Quadratic_Koch_2.svg/2880px-Quadratic_Koch_2.svg.png',
'https://upload.wikimedia.org/wikipedia/commons/thumb/5/55/Quadratic_Koch.svg/300px-Quadratic_Koch.svg.png'
]
?

sorry I using the jupyter notebook so not to know about the number

@zyzhoux
Copy link

zyzhoux commented Jun 21, 2022

Why are boxes of 2x2 excluded?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment