Slime Simulations - Create Stunning Animations from Simple Math

Introduction

After watching Sebastian Lague's most recent coding adventure, I was inspired to try my hand at the slime simulations in the 2nd half of the video. This page was another great resource and source of inspiration. I've messed around with generative art before (maybe a future post..) and I think it's one of the coolest things to work on.

The Basics

The previous links explain it pretty well, but just to sum it up:

a collection of "slime" particles move around and leave a trail everywhere they go
other particles are attracted to this trail
the slime trail spreads out and dissipates over time
the results are quite stunning for such a simple algorithm

The Code

I spent the last few hours writing the code and playing around with settings. It's still pretty basic and probably needs a lot of work, but that's what most of my projects look like. 😉
I tried to leave useful comments (denoted by #$%) and a few options to play around with. It doesn't run very fast because it's just one thread on the CPU, so maybe I'll adapt it with Numba's CUDA functionality later. Anyway, here it is:

import math
import random
import numpy as np
from numba import jit, stencil
from PIL import Image

seed = 0 #$% random seed
random.seed(seed)
step_size = 2 #$% forward motion per timestep
sensor_dist = 10 #$% sensor distance from particle
sensor_angle = math.pi / 8 #$% right and left sensor angle
rotation_angle = math.pi / 8 #$% particle turn angle
deposit_amount = 32 #$% additional intensity per particle per pixel per timestep
decay_factor = 0.75 #$% decay multiplier
step_fudge = 0.05 #$% randomize particle step
angle_fudge = math.pi / 32 #$% randomize particle angle
width = 400 #$% image width
height = 400 #$% image height
particle_count = 2 * width * height

def init_data():
    data = np.zeros((particle_count, 3), dtype=np.float32)
    for i in range(particle_count):
        #$% uniform distribution
        data[i] = [
            random.uniform(0, width-0), #$% x
            random.uniform(0, height-0), #$% y
            random.uniform(0, 2*math.pi) #$% angle
        ]
        #$% optional adjustable distribution 
        # data[i] = [
            # random.uniform(4 * width / 10, 6 * width / 10),
            # random.uniform(4 * height / 10, 6 * height / 10),
            # random.uniform(0, 2*math.pi)
        # ]
    return data

@jit
def sense(d, trail):
    x, y, phi = d
    
    #$% middle sensor
    dx1 = sensor_dist * math.cos(phi)
    dy1 = sensor_dist * math.sin(phi)
    x1 = round((x + dx1) % width)
    y1 = round((y + dy1) % height)
    s1 = trail[y1, x1]
    
    #$% left sensor
    dx2 = sensor_dist * math.cos(phi - sensor_angle)
    dy2 = sensor_dist * math.sin(phi - sensor_angle)
    x2 = round((x + dx2) % width)
    y2 = round((y + dy2) % height)
    s2 = trail[y2, x2]
    
    #$% right sensor
    dx3 = sensor_dist * math.cos(phi + sensor_angle)
    dy3 = sensor_dist * math.sin(phi + sensor_angle)
    x3 = round((x + dx3) % width)
    y3 = round((y + dy3) % height)
    s3 = trail[y3, x3]
    return (s1, s2, s3)

@jit
def move(data, trail):
    for d in data:
        s1, s2, s3 = sense(d, trail)
        if (s1 > s2) and (s1 > s3):
            #$% keep going straight
            pass
        elif (s2 > s1) and (s3 > s1):
            #$% random turn
            d[2] += random.uniform(-rotation_angle, rotation_angle)
        elif (s2 > s3) and (s2 > s1):
            #$% left turn
            d[2] -= rotation_angle
        elif (s3 > s2) and (s3 > s1):
            #$% right turn
            d[2] += rotation_angle
        
        #$% step forward & randomize step
        d[0] += step_size * math.cos(d[2]) + random.uniform(0, step_fudge)
        d[1] += step_size * math.sin(d[2]) + random.uniform(0, step_fudge)
        d[2] += random.uniform(-angle_fudge, angle_fudge) #$% randomize angle
        
        #$% wrap around
        d[0] %= width
        d[1] %= height

@jit
def deposit(data, trail):
    for d in data:
        x = min(max(round(d[0]), 0), width - 1)
        y = min(max(round(d[1]), 0), height - 1)
        #$% increase trail intensity up to max 255 after decay
        trail[y, x] = min(trail[y, x] + deposit_amount, int(255 / decay_factor))

@stencil
def kernel(a):
    #$% 3x3 average
    return (1 / 9) * (
        a[0, 0] + a[-1, 0] + a[0, -1] + 
        a[1, 0] + a[0, 1] + a[-1, -1] + 
        a[-1, 1] + a[1, -1] + a[1, 1]
    )

@jit
def diffuse(a):
    #$% apply stencil to the whole array
    return kernel(a).astype(np.int16)

@jit
def decay(a):
    #$% multiplicative decay
    return (a * decay_factor).astype(np.int16)

def do_image(trail):
    #$% create image from np array
    img = Image.fromarray(trail.astype(np.uint8))
    palette = []
    for i in range(256):
        palette.append(i) #$% red
        palette.append(i) #$% green
        palette.append(i) #$% blue
    # palette.reverse() #$% optional inversion
    img.putpalette(palette)
    return img

def make_gif(images, name):
    #$% make and save animated gif
    images[0].save(f'gifs/{name}.gif', save_all=True, append_images=images[1:], optimize=False, duration=40, loop=0)

def main():
    images = []
    data = init_data()
    trail = np.zeros((height, width), dtype=np.int16)
    for i in range(500):
        move(data, trail)
        deposit(data, trail)
        trail = diffuse(trail)
        trail = decay(trail)
        images.append(do_image(trail))
    name = f'{seed}_{step_size}_{sensor_dist}_{sensor_angle:.2f}_{rotation_angle:.2f}_{deposit_amount}_{decay_factor}_{width}_{height}_{particle_count}.png'
    make_gif(images, name)
    images[-1].save(f'pics/{name}') #$% save last image

if __name__ == '__main__':
    main()

The Results

Here are a few of the resulting images. There are infinite variations, and I certainly haven't found the best settings yet.

The real magic is in the animations, which are unfortunately hard to show in a post like this. Here's a small sample though. It's so creepily organic. I'll try to render a much larger video and post it in a comment when I can.