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

Reaction Diffusion Simulation (Python recipe) by FB36
ActiveState Code (http://code.activestate.com/recipes/579114/)

Reaction-Diffusion Simulation using Gray-Scott Model.

Python, 93 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
# Reaction-Diffusion Simulation Using Gray-Scott Model
# https://en.wikipedia.org/wiki/Reaction-diffusion_system
# FB - 20151016
import math
import random
from PIL import Image, ImageDraw
imgx = 256; imgy = 256 # image size
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)
pixels = image.load()
steps = 1000
dA = 1.0 # Diffusion Rate A
dB = 0.5 # Diffusion Rate B
f = 0.0545
k = 0.062
h = 0.01 # time step size

# 3x3 Convolution Matrix
weights = [[0.05, 0.2, 0.05], [0.2, -1.0, 0.2], [0.05, 0.2, 0.05]]

##arA = [[[1.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)]
##arB = [[[0.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)]
##arB[0][(imgy - 1) / 2][(imgx - 1) / 2] = 1.0 # seed
arA = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)]
arB = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)]

# simulation
p = 0; print "%" + str(p).zfill(2)
z = 1
for i in range(steps):
    z = 1 - z
    for iy in range(imgy):
        for ix in range(imgx):
            lapA = 0.0
            lapB = 0.0
            for jy in range(3):
                for jx in range(3):
                    kx = ix + (jx - 1)
                    ky = iy + (jy - 1)
                    if kx < 0: kx = imgx - 1
                    if kx > imgx - 1: kx = 0
                    if ky < 0: ky = imgy - 1
                    if ky > imgy - 1: ky = 0
                    lapA += arA[z][ky][kx] * weights[jy][jx]
                    lapB += arB[z][ky][kx] * weights[jy][jx]
            a = arA[z][iy][ix]
            b = arB[z][iy][ix]
            abb = a * b * b
            an = a + (dA * lapA - abb + f * (1.0 - a)) * h
            bn = b + (dB * lapB + abb - (k + f) * b) * h
            # prevent negative chemical amounts
            if an < 0: an = 0.0
            if bn < 0: bn = 0.0
            arA[1 - z][iy][ix] = an
            arB[1 - z][iy][ix] = bn

    pn = 100 * (i + 1) / steps # percent completed
    if pn != p:
        p = pn
        print "%" + str(p).zfill(2)

# paint the final state
z = 1 - z
aMin = arA[z][0][0]; aMax = aMin
bMin = arB[z][0][0]; bMax = bMin
for iy in range(imgy):
    for ix in range(imgx):
        a = arA[z][iy][ix]
        b = arB[z][iy][ix]
        if a < aMin: aMin = a
        if a > aMax: aMax = a
        if b < bMin: bMin = b
        if b > bMax: bMax = b

if aMin != aMax:
    for iy in range(imgy):
        for ix in range(imgx):
            a = arA[z][iy][ix]
            cA = int(255 * (a - aMin) / (aMax - aMin))
            pixels[ix, iy] = (cA, cA, cA)
    # label = "f = " + str(f) + " k = " + str(k)
    # draw.text((0, 0), label, (0, 255, 0))
    image.save("ReactionDiffusionSim_A.png", "PNG")

if bMin != bMax:
    for iy in range(imgy):
        for ix in range(imgx):
            b = arB[z][iy][ix]
            cB = int(255 * (b - bMin) / (bMax - bMin))
            pixels[ix, iy] = (cB, cB, cB)
    # label = "f = " + str(f) + " k = " + str(k)
    # draw.text((0, 0), label, (0, 255, 0))
    image.save("ReactionDiffusionSim_B.png", "PNG")
Tags: chaos, fractal, graphics, image, images, math, mathematics, physics, simulation

1 comment

FB36 (author) 8 years, 3 months ago  # | flag

Here is a much faster NumPy version:

# Reaction-Diffusion Simulation Using Gray-Scott Model
# https://en.wikipedia.org/wiki/Reaction-diffusion_system
# http://www.labri.fr/perso/nrougier/teaching/numpy/numpy.html#
# FB - 20160130
import random
import numpy as np
from PIL import Image, ImageDraw
n = 256
imgx = n; imgy = n # image size
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)
pixels = image.load()
steps = 10000

params = []
params.append((0.16, 0.08, 0.035, 0.065)) # Bacteria 1
params.append((0.14, 0.06, 0.035, 0.065)) # Bacteria 2
params.append((0.16, 0.08, 0.060, 0.062)) # Coral
params.append((0.19, 0.05, 0.060, 0.062)) # Fingerprint
params.append((0.10, 0.10, 0.018, 0.050)) # Spirals
params.append((0.12, 0.08, 0.020, 0.050)) # Spirals Dense
params.append((0.10, 0.16, 0.020, 0.050)) # Spirals Fast
params.append((0.16, 0.08, 0.020, 0.055)) # Unstable
params.append((0.16, 0.08, 0.050, 0.065)) # Worms 1
params.append((0.16, 0.08, 0.054, 0.063)) # Worms 2
params.append((0.16, 0.08, 0.035, 0.060)) # Zebrafish
(Du, Dv, F, k) = random.choice(params)

Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)])
U,V = Z['U'], Z['V']
u,v = U[1:-1,1:-1], V[1:-1,1:-1]

r = 20
u[...] = 1.0
U[n/2-r:n/2+r,n/2-r:n/2+r] = 0.50
V[n/2-r:n/2+r,n/2-r:n/2+r] = 0.25
u += 0.05*np.random.random((n,n))
v += 0.05*np.random.random((n,n))

p = 0
for i in xrange(steps):
    Lu = (                 U[0:-2,1:-1] +
          U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] +
                           U[2:  ,1:-1] )
    Lv = (                 V[0:-2,1:-1] +
          V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] +
                           V[2:  ,1:-1] )
    uvv = u*v*v
    u += (Du*Lu - uvv +  F   *(1-u))
    v += (Dv*Lv + uvv - (F+k)*v    )

    pn = 100 * (i + 1) / steps # percent completed
    if pn != p:
        p = pn
        print "%" + str(p).zfill(2)

# paint the final state
vMin=V.min(); vMax=V.max()
for iy in range(imgy):
    for ix in range(imgx):
        w = V[iy, ix]
        c = int(255 * (w - vMin) / (vMax - vMin))
        pixels[ix, iy] = (c, c, c)
label = "Du=" + str(Du) + " Dv=" + str(Dv) + " F=" + str(F) + " k=" + str(k)
draw.text((0, 0), label, (0, 255, 0))
image.save("ReactionDiffusionSim.png", "PNG")