In the late 1980's, the revival of chaotic systems, fractals and strange atractors was raising high expectations that some complex physical systems like turbulence or financial systems like the stock exchange markets could be described with. This transfer failed however.
But the fractals became in important tool in digital movie and game production to generate virtual landscapes. We implement a simple algorithm using midpoint displacement. It has some similarities with multrgrid methods used in computational fluid dynamics.
The steps to compute
- at the end we have an x-axis with some \(nx\) points, and each x-point has a y-value as mountain hight \(H\), that's all.
- we start with 2 points, each has an x-value and a y-value, \((x,y)\)
- now we compute a new point \((x_m, y_m)\) in the middle (by interpolation) and we add a random deviation to \(y_m\), so that new point is \((x_m, y_m + H \cdot random_m)\), where \(H\) is a scaling height.
- in the next iteration we repeat that for each line segment.
- with each iteration we reduce the scaling height \(H\)
- we stop when we think that we have enough points.
- prolongation is a different word for interpolation only.
- we use fast vectorized forms. So for the interpolation of all mid-points of a line we do not use any for-loop.
- we compute 3 skylines and display them in an overlapping semi-transparent style.
import numpy as np import pandas as pd from scipy.fftpack import fft, ifft import matplotlib.pyplot as plt import matplotlib.colors as mclr
def scale(a): return (a-np.min(a))/(np.max(a)-np.min(a)) def generate_data(NT): t0 = np.linspace(0,1,NT) t1 = t0 * 2 * np.pi t2 = np.pi*np.arange(NT)/2 y = np.sin(t1) return t1,y
def prolong1D(xa,za,nt): [nxa] = xa.shape nxb = nxa + (nxa-1) xb = np.zeros(nxb); zb = np.zeros(nxb) jx = slice(0,nxb,2); xb[jx] = xa; zb[jx] = za jx = slice(1,nxb,2) jxm = slice(0,nxb,2); jxp = slice(2,nxb,2) ixm = slice(0,nxa-1); ixp = slice(1,nxa ) H = (1-nt)**(NT/3.3) xb[jx] = 0.5*(xa[ixm]+xa[ixp]) zb[jx] = 0.5*(za[ixm]+za[ixp]) + (np.random.randn(zb[jx].shape)-0.0)*H return xb,zb
def fractalize(NT=9): x,z = generate_data(2) for nt in np.arange(NT): x,z = prolong1D(x,z,nt/NT) return x,z
def plot_fractals(NL,x,Y): fig, (ax0) = plt.subplots(1, 1, figsize=(25,5)); z_min = np.min(Y) ; z_max = np.max(Y)*1.2 x_min = np.min(x); x_max = np.max(x) plt.fill_between(x,z_min,z_max, color='grey',zorder=0,alpha=0.2) for nl in np.arange(NL): ax0.fill_between(x,Y[nl],z_min,cmap='plasma', zorder=nl+1, alpha=0.5) plt.axis('off'); plt.show() #pic_name='fractalmountain_B_2d_13.png' fig.savefig(pic_name, transparency=True)
def run(NL,NT): x,z = fractalize(NT) nz = z.shape Y = np.zeros((NL,nz)) for nl in np.arange(NL): x,z = fractalize(NT) x = scale(x) z = scale(z) Y[nl] = scale(z) * (1 - 0.2*nl) return x,Y
NL = 3 # number of lines NT = 8 # number of refinement for k in range(6): x,Y = run(NL,NT) plot_fractals(NL,x,Y) print();print();print();print();print();