Welcome, guest | Sign In | My Account | Store | Cart

IFS fractal dimension calculation (Python recipe) by FB36
ActiveState Code (http://code.activestate.com/recipes/578048/)

I found this equation in Wikipedia page for "List of fractals by Hausdorff dimension". The problem is how to calculate the scaling coefficients for any given IFS fractal. You can see the heuristic I used in the code. It gave correct values for all the fractals listed except for Centipede it gave 2.27.
I do not know its correct fractal dimension to compare though.

Python, 148 lines
Download
  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
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
# IFS fractal dimension calculation
# http://en.wikipedia.org/wiki/Iterated_function_system
# http://en.wikipedia.org/wiki/Fractal_dimension
# http://en.wikipedia.org/wiki/List_of_fractals_by_Hausdorff_dimension
# FB - 20120218
import math
import random
from PIL import Image
imgx = 512
imgy = 512 # will be auto-re-adjusted according to aspect ratio of the fractal
maxIt = imgx * imgy * 2

# Circumference of Polygon
# corners must be ordered in clockwise or counter-clockwise direction
def PolygonCircumference(corners):
    n = len(corners) # of corners
    circumference = 0.0
    for i in range(n):
        j = (i + 1) % n
        dx = corners[j][0] - corners[i][0]
        dy = corners[j][1] - corners[i][1]
        circumference += math.hypot(dx, dy)
    return circumference

##fractalName = "Barnsley Fern"
##mat=[[0.0, 0.0, 0.0, 0.16, 0.0, 0.0, 0.01],
##     [0.85, 0.04, -0.04, 0.85, 0.0, 1.6, 0.85],
##     [0.2, -0.26, 0.23, 0.22, 0.0, 1.6, 0.07],
##     [-0.15, 0.28, 0.26, 0.24, 0.0, 0.44, 0.07]]

##fractalName = "Centipede"
##mat = [[0.824074, 0.281482, -0.212346,  0.864198, -1.882290, -0.110607, 0.787473],
##       [0.088272, 0.520988, -0.463889, -0.377778,  0.785360,  8.095795, 0.212527]]

##fractalName = "Levy C curve"
##mat = [[0.5, -0.5, 0.5, 0.5, 0.0, 0.0, 0.5],
##       [0.5, 0.5, -0.5, 0.5, 0.5, 0.5, 0.5]]

##fractalName = "Dragon curve"
##mat = [[0.5, -0.5, 0.5, 0.5, 0.0, 0.0, 0.5],
##       [-0.5, -0.5, 0.5, -0.5, 1.0, 0.0, 0.5]]

fractalName = "Sierpinski Triangle"
mat = [[0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.33],
       [0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.33],
       [0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.33]]

##fractalName = "Sierpinski Carpet"
##mat = [[0.333, 0.0, 0.0, 0.333, 0.000, 0.999, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.333, 0.999, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.666, 0.999, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.000, 0.666, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.666, 0.666, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.000, 0.333, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.333, 0.333, 0.125],
##       [0.333, 0.0, 0.0, 0.333, 0.666, 0.333, 0.125]]

def IFS(x, y, i): # apply ith transformation to given point
    x0 = x * mat[i][0] + y * mat[i][1] + mat[i][4] 
    y  = x * mat[i][2] + y * mat[i][3] + mat[i][5] 
    x = x0
    return (x, y)
    
m = len(mat) # number of IFS transformations
# find xmin, xmax, ymin, ymax of the fractal using IFS algorithm
x = mat[0][4]
y = mat[0][5] 
xa = x
xb = x
ya = y
yb = y
for k in range(maxIt):
    i = random.randint(0, m - 1)
    if random.random() <= mat[i][6]:
        (x, y) = IFS(x, y, i)
        if x < xa:
            xa = x
        if x > xb:
            xb = x
        if y < ya:
            ya = y
        if y > yb:
            yb = y

imgy = int(imgy * (yb - ya) / (xb - xa)) # auto-re-adjust the aspect ratio 
image = Image.new("RGB", (imgx, imgy))
pixels = image.load()

# drawing using IFS algorithm
theColor = (255, 255, 255)
x=0.0
y=0.0 
for k in range(maxIt):
    i = random.randint(0, m - 1)
    if random.random() <= mat[i][6]:
        (x, y) = IFS(x, y, i)
        jx = int((x - xa) / (xb - xa) * (imgx - 1)) 
        jy = (imgy - 1) - int((y - ya) / (yb - ya) * (imgy - 1))
        if jx >= 0 and jx < imgx and jy >= 0 and jy < imgy:
            pixels[jx, jy] = theColor
image.save(fractalName + " fractal.png", "PNG")

# IFS fractal dimension calculation
# calculate scaling coefficients of all transformations
# scaling coefficient =
# circumference of bounding box after transformation applied /
# circumference of bounding box
boundingBoxCircumference = PolygonCircumference([(xa, ya), (xb, ya), (xb, yb), (xa, yb)])
coefficients = [] # list of scaling coefficients
for i in range(m):
    bBox = [IFS(xa, ya, i), IFS(xb, ya, i), IFS(xb, yb, i), IFS(xa, yb, i)]
    coefficients.append(PolygonCircumference(bBox) / boundingBoxCircumference)
    
# Calculate IFS fractal dimension equation value
# if the result is 0 then d is fractal dimension
def IFSfde(d):
    total = 0.0
    for c in coefficients:
        total += c ** d
    return total - 1.0

# Equation Solver using Bisection method
# http://en.wikipedia.org/wiki/Bisection_method
# function values must have opposite signs for f(a) and f(b)
def SolveEq(fun, a, b):
    eps = 1e-7
    ya = fun(a)
    yb = fun(b) 
    y0 = ya
    yc = yb
    while abs(yc - y0) > eps:
        y0 = yc
        c = (a + b) / 2.0
        yc = fun(c)
        if ya * yc < 0:
            b = c
            yb = yc
        elif yb * yc < 0:
            a = c
            ya = yc
        else:
            break
    return c

print "Fractal Name: " + fractalName
print "Scaling Coefficients:"
print coefficients
print "Estimated Fractal Dimension: " + str(SolveEq(IFSfde, 0.0, 3.0))
Tags: fractal, math, mathematics