Skip to content

Instantly share code, notes, and snippets.

@viveksck
Forked from rougier/fractal-dimension.py
Created March 14, 2017 13:07
Show Gist options
  • Star 18 You must be signed in to star a gist
  • Fork 5 You must be signed in to fork a gist
  • Save viveksck/1110dfca01e4ec2c608515f0d5a5b1d1 to your computer and use it in GitHub Desktop.
Save viveksck/1110dfca01e4ec2c608515f0d5a5b1d1 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)))
@bhatsrinidhi
Copy link

Hi, please could u explain the boxcount method of your code ?? I am confused with "return len(np.where((S > 0) & (S < k*k))[0])"
Thank you in advance

@ALXSIMINSKI
Copy link

ALXSIMINSKI commented Mar 19, 2019

please could u explain the boxcount method of your code ??

Hello. First of all, look at 65 example here https://pythonworld.ru/numpy/100-exercises.html. Here we have binary picture with values of pixels in [ 0, 1 ]. The sum of all pixels in every block is calculating and then we count, how many blocks have sum between 0 (means, that every pixel in this block is black (has value 0)) and kk (means, that each pixel in this block is white (has value 1), so, its sum will be equal to square of block). If block have sum between 0 and kk, then there is a white and black pixels inside it and it means that there is an intersection of object and background in this block.
The same thing in simple form:
len(np.where((S > 0) & (S < k*k))[0])
equals to:

count = 0
       for x_index in range(S.shape[0]):
            for y_index in range(S.shape[1]):
                if S[x_index][y_index] > 0 & S[x_index][y_index] < k * k:
                    count += 1;
return count

@ManuelAlejandrG
Copy link

Hello, Nice contribution.
I have a question.
Why do you use threshold?
Why threshold has 0.9 value?

Thank you in advance.

@dpatikas
Copy link

dpatikas commented Jan 19, 2020

Very clear and helpful contribution.
Please note that scipy.misc.imread has been discontinued after scipy v.1.2.0. It is suggested to use imageio.imread instead.
source

@Yutianzhang888
Copy link

Hello, Nice contribution.
I have a question.
Why do you use threshold?
Why threshold has 0.9 value?

Thank you in advance.

Hi, I also don't know the meaning of the threshold (0.9), and now do you understand the answer?

@Yutianzhang888
Copy link

There is one problem when I input the figure data(array), it always the three-dimension, could you tell me how to obtain the 2d array?

@data-hound
Copy link

data-hound commented Jul 21, 2020

@Yutianzhang888

There is one problem when I input the figure data(array), it always the three-dimension, could you tell me how to obtain the 2d array?

You need to select a single channel of the image.

Hello, Nice contribution.
I have a question.
Why do you use threshold?
Why threshold has 0.9 value?
Thank you in advance.

Hi, I also don't know the meaning of the threshold (0.9), and now do you understand the answer?

The threshold is used to convert the image to binary. Usually images have 8-bit color depths, which map to 256 values for each pixel. Binary image has only 2 shades - black or white.

@Nikita7416
Copy link

Hey guys! May be, I understand something wrong, but my PyCharm writes:

AttributeError: module 'scipy.misc' has no attribute 'imread'

@akshahi17
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'"

@benitez9rh
Copy link

@akashi @Nikita7416, dpatikas above said "Please note that scipy.misc.imread has been discontinued after scipy v.1.2.0. It is suggested to use imageio.imread instead."

@akshahi17
Copy link

I used "imageio.imread" but the problem seems to be not resolved.

Here what I get after running the program

Traceback (most recent call last):
File "C:\Users\Anup Shahi\Documents\pythonProject\venv\lib\site-packages\IPython\core\interactiveshell.py", line 3418, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "", line 1, in
runfile('C:/Users/Anup Shahi/Documents/pythonProject/FD.py', wdir='C:/Users/Anup Shahi/Documents/pythonProject')
File "C:\Program Files\JetBrains\PyCharm 2020.3\plugins\python\helpers\pydev_pydev_bundle\pydev_umd.py", line 197, in runfile
pydev_imports.execfile(filename, global_vars, local_vars) # execute the script
File "C:\Program Files\JetBrains\PyCharm 2020.3\plugins\python\helpers\pydev_pydev_imps_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "C:/Users/Anup Shahi/Documents/pythonProject/FD.py", line 51, in
print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(I))
File "C:/Users/Anup Shahi/Documents/pythonProject/FD.py", line 14, in fractal_dimension
assert(len(Z.shape) == 2)
AssertionError

@EssamWisam
Copy link

I used "imageio.imread" but the problem seems to be not resolved.

Here what I get after running the program

Traceback (most recent call last): File "C:\Users\Anup Shahi\Documents\pythonProject\venv\lib\site-packages\IPython\core\interactiveshell.py", line 3418, in run_code exec(code_obj, self.user_global_ns, self.user_ns) File "", line 1, in runfile('C:/Users/Anup Shahi/Documents/pythonProject/FD.py', wdir='C:/Users/Anup Shahi/Documents/pythonProject') File "C:\Program Files\JetBrains\PyCharm 2020.3\plugins\python\helpers\pydev_pydev_bundle\pydev_umd.py", line 197, in runfile pydev_imports.execfile(filename, global_vars, local_vars) # execute the script File "C:\Program Files\JetBrains\PyCharm 2020.3\plugins\python\helpers\pydev_pydev_imps_pydev_execfile.py", line 18, in execfile exec(compile(contents+"\n", file, 'exec'), glob, loc) File "C:/Users/Anup Shahi/Documents/pythonProject/FD.py", line 51, in print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(I)) File "C:/Users/Anup Shahi/Documents/pythonProject/FD.py", line 14, in fractal_dimension assert(len(Z.shape) == 2) AssertionError

Try

import cv2
img = cv2.imread("sierpinski.png", 0)/256.0
print("Minkowski–Bouligand dimension (computed): ", fractal_dimension(img, 0.9))

@EssamWisam
Copy link

Thank you for sharing this.

I believe we can reduce

    # 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))

to

    # largest n for which 2^n is less than or equal p
    n = int(np.floor(np.log2(p)))

@antoni27
Copy link

antoni27 commented May 13, 2022

sorry, want to ask if this can be used on iris recognition or is there any other version?
btw the image is
forward_look(278)120
forward_look (278)

thanks before and really appreciate it

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