// Ported from http://cs.brown.edu/people/pfelzens/dt/
sclass DistanceTransform {
  static float INF = 1E20f;

  // dt of 1d function using squared distance
  float[] dt(float[] f, int n) {
    float[] d = new float[n], z = new float[n+1];
    int[] v = new int[n];
    int k = 0;
    z[0] = -INF;
    z[1] = INF;
    for (int q = 1; q <= n-1; q++) {
      float s = ((f[q]+sqr(q))-(f[v[k]]+sqr(v[k])))/(2*q-2*v[k]);
      while (s <= z[k]) {
        k--;
        s = ((f[q]+sqr(q))-(f[v[k]]+sqr(v[k])))/(2*q-2*v[k]);
      }
      k++;
      v[k] = q;
      z[k] = s;
      z[k+1] = INF;
    }
  
    k = 0;
    for (int q = 0; q <= n-1; q++) {
      while (z[k+1] < q)
        k++;
      d[q] = sqr(q-v[k]) + f[v[k]];
    }
  
    ret d;
  }
  
  // dt of 2d function using squared distance
  void dt(FloatBWImage im) {
    int width = im.getWidth();
    int height = im.getHeight();
    float[] f = new[max(width,height)];
  
    // transform along columns
    for x to width: {
      for y to height:
        f[y] = im.getPixel(x, y);
      float[] d = dt(f, height);
      for y to height:
        im.setPixel(x, y, d[y]);
    }
  
    // transform along rows
    for y to height: {
      for x to width:
        f[x] = im.getPixel(x, y);
      float[] d = dt(f, width);
      for x to width:
        im.setPixel(x, y, d[x]);
    }
  }
  
  /* dt of binary image using squared distance */
  FloatBWImage dt(BWImage im, float threshold default .5f, bool brightIsOn) {
    int width = im.getWidth();
    int height = im.getHeight();
  
    var out = new FloatBWImage(width, height);
    for y to height:
      for x to width:
        out.setPixel(x, y,
          (im.getFloatPixel(x, y) >= threshold == brightIsOn)
            ? 0 : INF);

    dt(out);
    ret out;
  }
}