int w, h, d = 3, e = 24, n;      float[][][] X;      float[][] Y;
int[] r, c, q;                   PImage I, J;        String f = "susan";

import java.util.*;

void setup() { size(1440, 1440); w = width; h = height; colorMode(RGB, 256);
  X = new float[w][h][d]; Y = new float[w][h];
  I = loadImage(f + ".jpeg"); J = createImage(w, h, RGB);
  for (int j = 0; j < h; j++) { for (int i = 0; i < w; i++) {
    color lab = rgb2lab(I.get(i,j)); Y[i][j] = red(lab); J.set(i,j,color(Y[i][j]));
    X[i][j][0] = Y[i][j] -  ceil(e/2.0); X[i][j][1] = Y[i][j];
    X[i][j][2] = Y[i][j] + floor(e/2.0); } }
  J.updatePixels(); image(J, 0, 0);
  n = (w-2)*(h-2); r = new int[n]; c = new int[n]; q = new int[n];
  for (int j = 1; j < h-1; j++) { for (int i = 1; i < w-1; i++) {
    int k = i-1 + (w-2)*(j-1); c[k] = i; r[k] = j; q[k] = k; } } }

void draw() { float t = 0.0, s = 0.025; Collections.shuffle(Arrays.asList(q));
for (int k = 0; k < n-1; k++) { int i = c[q[k]]; int j = r[q[k]];
  float a = X[i][j][0]; float b = X[i][j][2]; float[] y = new float[4];
  y[0] = (min(max(X[i-1][j  ][1], a), b) + min(max(X[i+1][j  ][1], a), b))/2.0;
  y[1] = (min(max(X[i  ][j-1][1], a), b) + min(max(X[i  ][j+1][1], a), b))/2.0;
  y[2] = (min(max(X[i-1][j-1][1], a), b) + min(max(X[i+1][j+1][1], a), b))/2.0;
  y[3] = (min(max(X[i-1][j+1][1], a), b) + min(max(X[i+1][j-1][1], a), b))/2.0;
  float x = pow(2, 10); for (int q = 0; q < 4; q++) { 
    float p = y[q] - X[i][j][1]; if (abs(x) > abs(p)) { x = p; } }
  t = max(t, abs(x)); X[i][j][1] += x/2.0; }
println(t); if (t < s) { for (int j = 0; j < h; j++) { for (int i = 0; i < w; i++) {
    J.set(i, j, color((min(max(X[i][j][1], 0.0), 255.0) + Y[i][j])/2.0)); } }
  J.updatePixels(); J.save(f + ".tif"); exit(); } }

color rgb2lab(color rgb) { float[] X = new float[3];
  float[] C = { red(rgb)/255.0, green(rgb)/255.0, blue(rgb)/255.0 };
  for (int i = 0; i < 3; i++) { if (C[i] <= 0.0405) { C[i] = C[i]/12.92;
    } else { C[i] = pow((C[i] + 0.055)/1.055, 2.4); } }
  X[0] = (0.4124*C[0] + 0.3576*C[1] + 0.1805*C[2]) / 0.95047;
  X[1] =  0.2126*C[0] + 0.7152*C[1] + 0.0722*C[2];
  X[2] = (0.0193*C[0] + 0.1192*C[1] + 0.9505*C[2]) / 1.08883;
  for (int i = 0; i < 3; i++) { if (X[i] > 0.008856) { X[i] = pow(X[i], 1.0/3.0);
    } else { X[i] = 7.787*X[i] + 16.0/116.0; } }
  return color(round(min(255.0, max(0.0, 2.55*(116.0*X[1] - 16.0)))),
               round(min(255.0, max(0.0, 2.55*(500.0*(X[0] - X[1]))))),
               round(min(255.0, max(0.0, 2.55*(200.0*(X[1] - X[2])))))); }