`#! /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