forked from rougier/from-python-to-numpy
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathfractal_dimension.py
More file actions
60 lines (51 loc) · 2.11 KB
/
fractal_dimension.py
File metadata and controls
60 lines (51 loc) · 2.11 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# -----------------------------------------------------------------------------
# From Numpy to Python
# Copyright (2017) Nicolas P. Rougier - BSD license
# More information at https://github.com/rougier/numpy-book
# -----------------------------------------------------------------------------
import numpy as np
def fractal_dimension(Z, threshold=0.9):
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)
return len(np.where((S > 0) & (S < k*k))[0])
Z = (Z < threshold)
p = min(Z.shape)
n = 2**np.floor(np.log(p)/np.log(2))
n = int(np.log(n)/np.log(2))
sizes = 2**np.arange(n, 1, -1)
counts = []
for size in sizes:
counts.append(boxcount(Z, size))
coeffs = np.polyfit(np.log(sizes), np.log(counts), 1)
return -coeffs[0]
if __name__ == '__main__':
from scipy import misc
import matplotlib.pyplot as plt
import matplotlib.patches as patches
Z = 1.0 - misc.imread("../data/Great-Britain.png")/255
print(fractal_dimension(Z, threshold=0.25))
sizes = 128, 64, 32
xmin, xmax = 0, Z.shape[1]
ymin, ymax = 0, Z.shape[0]
fig = plt.figure(figsize=(10, 5))
for i, size in enumerate(sizes):
ax = plt.subplot(1, len(sizes), i+1, frameon=False)
ax.imshow(1-Z, plt.cm.gray, interpolation="bicubic", vmin=0, vmax=1,
extent=[xmin, xmax, ymin, ymax], origin="upper")
ax.set_xticks([])
ax.set_yticks([])
for y in range(Z.shape[0]//size+1):
for x in range(Z.shape[1]//size+1):
s = (Z[y*size:(y+1)*size, x*size:(x+1)*size] > 0.25).sum()
if s > 0 and s < size*size:
rect = patches.Rectangle(
(x*size, Z.shape[0]-1-(y+1)*size),
width=size, height=size,
linewidth=.5, edgecolor='.25',
facecolor='.75', alpha=.5)
ax.add_patch(rect)
plt.tight_layout()
plt.savefig("fractal-dimension.png")
plt.show()