# Multifractal analysis of an image using wavelets

I want to measure the fractal character of images using a wavelet approach (if possible) in python. (E.g., maybe things like roughness, anisotropy, or just the fractal dimension - perhaps just separating images as fractal vs. not-fractal)

The goal is to produce a few measurements for each image (there are thousands) which I can use for spatial analysis or clustering/classification purposes.

I am stuck with the current code (see error below) but I suspect there is a better way to do this..

Here is an example image (whick I pre-processed to binary):

import numpy as np
from math import log
import matplotlib.pyplot as plt

from skimage import io

import pywt
import pywt.data

# NOTE: I have this as a numpy array - but showing here as a png Using code from Pywavelets I decomposed the image into coefficients (LL, LH, HL, HH):

# Wavelet transform of image, and plot approximation and details
titles = ['Approximation', ' Horizontal detail',
'Vertical detail', 'Diagonal detail']

coeffs = pywt.dwt2(mx, 'haar') # using haar b/c of binary input array

LL, (LH, HL, HH) = coeffs
fig = plt.figure(figsize=(12, 3))

for i, a in enumerate([LL, LH, HL, HH]):
ax = fig.add_subplot(1, 4, i + 1)
ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
ax.set_title(titles[i], fontsize=10)
ax.set_xticks([])
ax.set_yticks([])

fig.tight_layout()
plt.show() Using modified code from an ancient newsgroup (converting calls to the 'numerical' library to numpy):

# Note: I have hashed out the older code, and replaced the line below
# Included here in case I have messed something up!

#M = np.array([[0,-1], [1,0]]) * np.sqrt(3)/2
M = LH # the 'Low-High' pass wavelet coefficient
#------------------------------------------------------------------

def N(points, scale):
"""Return (num of coverage, scale) of points by boxes of size 1/scale"""
unique = {}
for point in points:
#box = tuple((point * scale).astype(Int))
box = tuple((point * scale).astype(int))
unique[box] = 1
return float(len(unique)), scale

def dim(points):
""" Calculate dimensions of points at various scale = 2**level"""
f0 = 1, 1
#for level in xrange(1, 12):
for level in range(1, 12):
f1 = N(points, 2.0**level)
dim = log(f1/f0) / log(f1/f0)
print ("%2d:  %.4g" % (level, dim))
f0 = f1

def makekoch(x1, x2, level, flist=[]):
""" Return a list of points on Koch curve, subdivided levels """
if level == 0: return flist
xd = (x2 - x1) / 3.
xa = x1 + xd
xb = x2 - xd
#xm = (x1 + x2) / 2.0 - matrixmultiply(M, xd)
xm = (x1 + x2) / 2.0 - np.dot(M, xd)
flist.append(xa)
makekoch(x1, xa, level-1, flist)
makekoch(xa, xm, level-1, flist)
flist.append(xm)
makekoch(xm, xb, level-1, flist)
makekoch(xb, x2, level-1, flist)
flist.append(xb)
return flist

#------------------------------------------------------------------
if __name__ == "__main__":
print ("line dimension")
#x = (np.arange(100)/100.)[:,NewAxis]
x = (np.arange(100)/100.)[:,np.newaxis]
dim(x)

print ("plane dimension")
list2d = []
#for x in xrange(100):
for x in range(100):
#for y in xrange(100):
for y in range(100):
list2d.append(np.array((x, y))/100.0)
dim(list2d)

print ("space dimension")
list3d = []
#for x in xrange(30):
for x in range(30):
#for y in xrange(30):
for y in range(30):
#for z in xrange(30):
for z in range(30):
list3d.append(np.array((x, y, z))/30.0)
dim(list3d)

print ("koch dimension")
listkoch = makekoch(np.array((0,0)), np.array((1,0)), 7)
dim(listkoch)


Unfortunately this results in an error when running the koch dimension function:

ValueError: shapes (752,1002) and (2,) not aligned: 1002 (dim 1) != 2 (dim 0)


And I am not clever enough to know how to modify xd to be adaptable to variously shaped 2D arrays!