# 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")
Diff to Previous Revision
--- revision 1 2015-10-16 19:47:40
+++ revision 2 2015-10-16 19:52:02
@@ -75,7 +75,7 @@
if aMin != aMax:
for iy in range(imgy):
for ix in range(imgx):
- a = arA[1 - z][iy][ix]
+ 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)
@@ -85,7 +85,7 @@
if bMin != bMax:
for iy in range(imgy):
for ix in range(imgx):
- b = arB[1 - z][iy][ix]
+ 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)
