🦠 Conway's game of life(ish) in 3D, with OpenVDB output for every generation

conway_3d.py

#!/usr/bin/env python3

import pyopenvdb as vdb
import numpy as np
import scipy.signal

# PARAMETERS

kernel = '''
XXX|XXX|XXX
XXX|X X|XXX
XXX|XXX|XXX
'''

rule = lambda state, neighbors: \
    (neighbors == 5) | (~state & (neighbors == 4))

generations = 20
filename = lambda n: f'test-{n:03}.vdb'

initial = '''
     |     |     |     |     
     |     |  X  |     |     
     |  X  | XXX |  X  |     
     |     |  X  |     |     
     |     |     |     |     
'''

# VOXEL PARSING

def parse_3d(desc: str) -> np.array:
    if desc.startswith('\n'): desc = desc[1:]

    # split segments
    lines = desc.splitlines()
    if not lines: return np.zeros((0, 0, 0))
    segs = [ line.split('|') for line in lines ]

    # verify conformance
    markers = set( tuple(len(x) for x in line[:-1]) for line in segs )
    assert len(markers) == 1
    markers = set( next(iter(markers)) )
    assert len(markers) == 1
    dim = next(iter(markers))
    ldim = max( len(line[-1].rstrip()) for line in segs )
    assert ldim <= dim

    # parse segments into array
    parse_seg = lambda seg: [ { ' ': False, 'X': True }[x] for x in seg[:dim].ljust(dim) ]
    parse_line = lambda line: [ parse_seg(seg) for seg in line ]
    return np.array([ parse_line(line) for line in segs ])

# SIMULATION

kernel = parse_3d(kernel).astype('uint8')

def autocrop(arr: np.array) -> tuple[np.array, tuple[int, ...]]:
    slices = ()
    for axis in range(arr.ndim):
        bound = arr
        for other_axis in reversed(range(arr.ndim)):
            if axis != other_axis:
                bound = np.any(bound, other_axis)
        idxs = np.where(bound)[0]
        if not len(idxs): return None, None
        slices += ( slice(*(idxs[[0, -1]] + (0, 1))), )
    return arr[slices], tuple( s.start for s in slices )

state = parse_3d(initial)
pos = -np.array(state.shape) // 2

for i in range(generations):
    # crop and pad state
    state, dpos = autocrop(state)
    if state is None: break
    pos += dpos
    state = np.pad(state, 1)
    pos -= 1

    # compute next generation
    print(f'Generation {i+1:3} / {generations:3}... {state.shape}')
    neighbors = scipy.signal.convolve(state.astype('uint8'), kernel, 'same')
    oldstate, state = state, rule(state, neighbors)

    # output OpenVDB
    state_grid = vdb.BoolGrid()
    state_grid.name = 'state'
    state_grid.copyFromArray(oldstate, ijk=tuple(map(int, pos)))
    vdb.write(filename(i + 1), grids=[state_grid])