Macuyiko · January 8, 2017 16:33 · Jan 8, 2017
diff --git a/smoothlife.pde b/smoothlife.pde
@@ -0,0 +1,315 @@
+//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;
+}