Saturday, June 21, 2008

Python: Isotropic fractal surface generator

#! /usr/bin/env python
from __future__ import division
import Image
from scipy import *

class FractalSurface(object):
'''Generate isotropic fractal surface image using
spectral synthesis method [1, p.]
References:
1. Yuval Fisher, Michael McGuire,
The Science of Fractal Images, 1988
'''

def __init__(self,fd=2.5, N=256):
self.N=N
self.H=1-(fd-2);
self.X=zeros((self.N,self.N),complex)
self.A=zeros((self.N,self.N),complex)
self.img=Image.Image()

def genSurface(self):
'''Spectral synthesis method
'''
N=self.N; A=self.A
powerr=-(self.H+1.0)/2.0

for i in range(int(N/2)+1):
for j in range(int(N/2)+1):

phase=2*pi*rand()

if i is not 0 or j is not 0:
rad=(i*i+j*j)**powerr*random.normal()
else:
rad=0.0

self.A[i,j]=complex(rad*cos(phase),rad*sin(phase))

if i is 0:
i0=0.0
else:
i0=N-i
if j is 0:
j0=0.0
else:
j0=N-j

self.A[i0,j0]=complex(rad*cos(phase),-rad*sin(phase))


self.A.imag[N/2][0]=0.0
self.A.imag[0,N/2]=0.0
self.A.imag[N/2][N/2]=0.0

for i in range(1,int(N/2)):
for j in range(1,int(N/2)):
phase=2*pi*rand()
rad=(i*i+j*j)**powerr*random.normal()
self.A[i,N-j]=complex(rad*cos(phase),rad*sin(phase))
self.A[N-i,j]=complex(rad*cos(phase),-rad*sin(phase))

itemp=fftpack.ifft2(self.A)
itemp=itemp-itemp.min()
self.X=itemp


def genImage(self):
#Aa=abs(Aa)
Aa=self.X
im=Aa.real/Aa.real.max()*255.0
self.img=Image.fromarray(uint8(im))
#img2=Image.fromstring("L",(N,N),uint8(im).tostring())

def showImg(self):
self.img.show()

def saveImg(self,fname="fs.tiff"):
self.img.save(fname)

def getFSimg(self):
return self.img

def main():
fs=FractalSurface()
fs.genSurface()
fs.genImage()
fs.saveImg()
fs.showImg()

if __name__ == '__main__':
main()

Example


FD=2.1




FD=2.5




FD=2.9