Sun 08 January 2017, by Seppe "Macuyiko" vanden Broucke
SmoothLife in Processing

Conway’s Game of Life is one of the most popular examples of cellular automaton. However, did you know that Stephan Rafler proposed a version of Game of Life in 2011 which works over a contious domain, called SmoothLife?

Cool, isn’t it? The basic idea is that the grid of cells is replaced here with an effective grid where each cell occupies a continuous coordinate with a very small finite size and with the neighbors calculated based on a radius around that cell.

This blog post summarised the idea well and also provides a JavaScript implementation.

Since I’m toying around with Processing again these days, it seemed like a fun quick project to convert the source to Java:

//Quick adaption from:
//- http://jsfiddle.net/mikola/aj2vq/
//- https://0fps.net/2012/11/19/conways-game-of-life-for-curved-surfaces-part-1/
import java.util.*;
int INNER_RADIUS = 7;
int OUTER_RADIUS = 3 * INNER_RADIUS;
double B1 = 0.278;
double B2 = 0.365;
double D1 = 0.267;
double D2 = 0.445;
double ALPHA_N = 0.028;
double ALPHA_M = 0.147;
int LOG_RES = 9;
int TRUE_RES = (1<<LOG_RES);
int[] COLOR_SHIFT = new int[]{0, 0, 0};
int[] COLOR_SCALE = new int[]{256, 256, 256};
int NR_OF_FIELDS = 2;
int current_field;
List<double[][]> fields = new ArrayList<double[][]>();
double[][] imaginary_field, M_re_buffer, M_im_buffer, N_re_buffer, N_im_buffer;
double[][] M_re, M_im, N_re, N_im;
void setup() {
size(800, 800);
current_field = 0;
for (int i = 0; i < NR_OF_FIELDS; i++) {
fields.add(zeromatrix());
}
imaginary_field = zeromatrix();
M_re_buffer = zeromatrix();
M_im_buffer = zeromatrix();
N_re_buffer = zeromatrix();
N_im_buffer = zeromatrix();
//Precalculate multipliers for m,n
BesselJ inner_bessel = BesselJ(INNER_RADIUS);
BesselJ outer_bessel = BesselJ(OUTER_RADIUS);
double inner_w = 1.0 / inner_bessel.w;
double outer_w = 1.0 / (outer_bessel.w - inner_bessel.w);
M_re = inner_bessel.re;
M_im = inner_bessel.im;
N_re = outer_bessel.re;
N_im = outer_bessel.im;
for (int i=0; i < TRUE_RES; ++i) {
for (int j=0; j < TRUE_RES; ++j) {
N_re[i][j] = outer_w * (N_re[i][j] - M_re[i][j]);
N_im[i][j] = outer_w * (N_im[i][j] - M_im[i][j]);
M_re[i][j] *= inner_w;
M_im[i][j] *= inner_w;
}
}
background(0);
frameRate(30);
init_field(0);
speckle_field(2000, 1);
}
void draw() {
perform_calculation_step();
double[][] cur_field = fields.get(current_field);
if (mousePressed){
int v = (int) (((double) mouseX / (double) width) * (double) TRUE_RES);
int u = (int) (((double) mouseY / (double) height) * (double) TRUE_RES);
for (int x = -OUTER_RADIUS/2; x < OUTER_RADIUS/2; x++) {
for (int y = -OUTER_RADIUS/2; y < OUTER_RADIUS/2; y++) {
cur_field[u+x][v+y] = 1;
}
}
for (int x = -INNER_RADIUS/2; x < INNER_RADIUS/2; x++) {
for (int y = -INNER_RADIUS/2; y < INNER_RADIUS/2; y++) {
cur_field[u+x][v+y] = 0;
}
}
}
int image_ptr = 0;
PImage img = createImage(TRUE_RES, TRUE_RES, RGB);
img.loadPixels();
for (int i = 0; i < TRUE_RES; i++) {
for (int j = 0; j < TRUE_RES; j++) {
double s = cur_field[i][j];
int[] clr = new int[]{255, 255, 255};
for (int k = 0; k < 3; k++) {
clr[k] = (int) Math.max(0, Math.min(255, Math.floor(COLOR_SHIFT[k] + COLOR_SCALE[k]*s)));
}
img.pixels[image_ptr] = color(clr[0], clr[1], clr[2]);
image_ptr++;
}
}
img.updatePixels();
//Blit to screen
img.resize(width, height);
image(img, 0, 0);
}
void init_field(double x) {
double[][] cur_field = fields.get(current_field);
for (int i = 0; i < TRUE_RES; i++) {
for (int j = 0; j < TRUE_RES; j++) {
cur_field[i][j] = x;
}
}
}
void speckle_field(int count, double intensity) {
double[][] cur_field = fields.get(current_field);
for (int i = 0; i < count; i++) {
int u = (int) Math.floor(Math.random() * (TRUE_RES - INNER_RADIUS));
int v = (int) Math.floor(Math.random() * (TRUE_RES - INNER_RADIUS));
for (int x = 0; x < INNER_RADIUS; x++) {
for (int y = 0; y < INNER_RADIUS; y++) {
cur_field[u+x][v+y] = intensity;
}
}
}
}
void perform_calculation_step() {
//Read in fields
double[][] cur_field = fields.get(current_field);
current_field = (current_field + 1) % fields.size();
double[][] next_field = fields.get(current_field);
//Clear extra imaginary field
for (int i = 0; i < TRUE_RES; i++) {
for (int j = 0; j < TRUE_RES; j++) {
imaginary_field[i][j] = 0.0;
}
}
//Compute m,n fields
fft2(true, LOG_RES, cur_field, imaginary_field);
field_multiply(cur_field, imaginary_field, M_re, M_im, M_re_buffer, M_im_buffer);
fft2(false, LOG_RES, M_re_buffer, M_im_buffer);
field_multiply(cur_field, imaginary_field, N_re, N_im, N_re_buffer, N_im_buffer);
fft2(false, LOG_RES, N_re_buffer, N_im_buffer);
//Step s
for (int i=0; i<next_field.length; ++i) {
for (int j=0; j<next_field.length; ++j) {
next_field[i][j] = S(N_re_buffer[i][j], M_re_buffer[i][j]);
}
}
}
class BesselJ {
double[][] re;
double[][] im;
double w;
BesselJ(double[][] re, double[][] im, double w) {
this.re = re;
this.im = im;
this.w = w;
}
}
BesselJ BesselJ(double radius) {
double[][] field = zeromatrix();
double weight = 0.0;
for (int i = 0; i < field.length; i++) {
for (int j = 0; j < field.length; j++) {
double ii = ((i + field.length/2) % field.length) - field.length/2;
double jj = ((j + field.length/2) % field.length) - field.length/2;
double r = Math.sqrt(ii*ii + jj*jj) - radius;
double v = 1.0 / (1.0 + Math.exp(LOG_RES * r));
weight += v;
field[i][j] = v;
}
}
double[][] imag_field = zeromatrix();
fft2(true, LOG_RES, field, imag_field);
return new BesselJ(field, imag_field, weight);
}
double sigma(double x, double a, double alpha) {
return 1.0 / (1.0 + Math.exp(-4.0/alpha * (x - a)));
}
double sigma_2(double x, double a, double b) {
return sigma(x, a, ALPHA_N) * (1.0 - sigma(x, b, ALPHA_N));
}
double lerp(double a, double b, double t) {
return (1.0-t)*a + t*b;
}
double S(double n, double m) {
double alive = sigma(m, 0.5, ALPHA_M);
return sigma_2(n, lerp(B1, D1, alive), lerp(B2, D2, alive));
}
void field_multiply(double[][] a_r, double[][] a_i,
double[][] b_r, double[][] b_i,
double[][] c_r, double[][] c_i) {
for (int i = 0; i < TRUE_RES; i++) {
double[] Ar = a_r[i];
double[] Ai = a_i[i];
double[] Br = b_r[i];
double[] Bi = b_i[i];
double[] Cr = c_r[i];
double[] Ci = c_i[i];
for (int j = 0; j < TRUE_RES; j++) {
double a = Ar[j];
double b = Ai[j];
double c = Br[j];
double d = Bi[j];
double t = a * (c + d);
Cr[j] = t - d*(a+b);
Ci[j] = t + c*(b-a);
}
}
}
void fft(boolean dir, int m, double[] x, double[] y) {
int nn, i, i1, j, k, i2, l, l1, l2;
double c1, c2, tx, ty, t1, t2, u1, u2, z;
nn = x.length;
/* Do the bit reversal */
i2 = nn >> 1;
j = 0;
for (i=0; i<nn-1; i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l = 0; l < m; l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0; j<l1; j++) {
for (i=j; i<nn; i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = Math.sqrt((1.0 - c1) / 2.0);
if (dir)
c2 = -c2;
c1 = Math.sqrt((1.0 + c1) / 2.0);
}
/* Scaling for forward transform */
if (!dir) {
double scale_f = 1.0 / nn;
for (i=0; i<nn; i++) {
x[i] *= scale_f;
y[i] *= scale_f;
}
}
}
void fft2(boolean dir, int m, double[][] x, double[][] y) {
/* In place 2D FFT */
for (int i = 0; i < x.length; ++i) {
fft(dir, m, x[i], y[i]);
}
for (int i=0; i<x.length; ++i) {
for (int j=0; j<i; ++j) {
double t = x[i][j];
x[i][j] = x[j][i];
x[j][i] = t;
}
}
for (int i=0; i<y.length; ++i) {
for (int j=0; j<i; ++j) {
double t = y[i][j];
y[i][j] = y[j][i];
y[j][i] = t;
}
}
for (int i=0; i < x.length; ++i) {
fft(dir, m, x[i], y[i]);
}
}
double[][] zeromatrix() {
return matrix(TRUE_RES, TRUE_RES, 0.0);
}
static double[][] matrix(int rows, int cols, double value) {
double[][] matrix = new double[rows][cols];
for (int row = 0; row < rows; row++)
for (int col = 0; col < cols; col++)
matrix[row][col] = value;
return matrix;
}
view raw smoothlife.pde hosted with ❤ by GitHub

Here’s what the results look like. Higher resolutions are possible but can’t be rendered in real-time due to the slower FFT implementation (it would be fun to import FFTW and see how that fares).

Here’s what the result looks like: