2D Interpolation, Part 1: The Digital Differential Analyzer

To my Planet OCaml readers, I apologize for veering into Java. I've been playing with Processing of lately, because it's far easier to prototype silly, simple animations (fractals! Fractals everywhere!) in it than by using OCamlSDL, say. Last time I showed a basic 2D interpolator; now it's time to start massaging it into shape. Let's recall the original code:

void interp0() {
  final float dy = (float) (ny - 3) / height;
  final float dx = (float) (nx - 3) / width;
  int offset = 0;
  for (int i = 0; i < height; i++) {
    float yr = 1 + dy * i;
    int yi = (int) yr;
    yr -= yi;
    for (int j = 0; j < width; j++) {
      float xr = 1 + dx * j;
      int xi = (int) xr;
      xr -= xi;
      float c = 0;
      switch (method) {
      case NEAREST:
        c = data[yi + (int) (yr + 0.5)][xi + (int) (xr + 0.5)];
        break;
      case BILINEAR:
        c = linear(yr,
            linear(xr, data[yi+0][xi+0], data[yi+0][xi+1]),
            linear(xr, data[yi+1][xi+0], data[yi+1][xi+1]) );
        break;
      case BICUBIC:
        c = cubic(yr,
            cubic(xr, data[yi-1][xi-1], data[yi-1][xi+0], data[yi-1][xi+1], data[yi-1][xi+2]),
            cubic(xr, data[yi+0][xi-1], data[yi+0][xi+0], data[yi+0][xi+1], data[yi+0][xi+2]),
            cubic(xr, data[yi+1][xi-1], data[yi+1][xi+0], data[yi+1][xi+1], data[yi+1][xi+2]),
            cubic(xr, data[yi+2][xi-1], data[yi+2][xi+0], data[yi+2][xi+1], data[yi+2][xi+2]) );
        break;
      }
      pixels[offset++] = hue(c);
    }
  }
}

The coordinate transformations from destination array to source matrix involve a linear transformation s = s₀ + k·(d - d₀), where k is a floating-point approximation to a rational number. This can be done using the all-integer Digital Differential Analyzer or DDA. First of all, note that there's no reason for yr, xr to start at 1; only yi, xi need to be offset:

void interp0() {
  final float dy = (float) (ny - 3) / height;
  final float dx = (float) (nx - 3) / width;
  int offset = 0;
  for (int i = 0; i < height; i++) {
    float yr = dy * i;
    int yi = (int) yr;
    yr -= yi;
    yi += 1;
    for (int j = 0; j < width; j++) {
      float xr = dx * j;
      int xi = (int) xr;
      xr -= xi;
      xi += 1;
      float c = 0;
      switch (method) {
      case NEAREST:
        c = data[yi + (int) (yr + 0.5)][xi + (int) (xr + 0.5)];
        break;
      case BILINEAR:
        c = linear(yr,
            linear(xr, data[yi+0][xi+0], data[yi+0][xi+1]),
            linear(xr, data[yi+1][xi+0], data[yi+1][xi+1]) );
        break;
      case BICUBIC:
        c = cubic(yr,
            cubic(xr, data[yi-1][xi-1], data[yi-1][xi+0], data[yi-1][xi+1], data[yi-1][xi+2]),
            cubic(xr, data[yi+0][xi-1], data[yi+0][xi+0], data[yi+0][xi+1], data[yi+0][xi+2]),
            cubic(xr, data[yi+1][xi-1], data[yi+1][xi+0], data[yi+1][xi+1], data[yi+1][xi+2]),
            cubic(xr, data[yi+2][xi-1], data[yi+2][xi+0], data[yi+2][xi+1], data[yi+2][xi+2]) );
        break;
      }
      pixels[offset++] = hue(c);
    }
  }
}

At this point you will notice error pixels along some contours. Those are the pixels for which the floating-point quotient (ny−3)/height, resp. (nx−3)/width is not exactly representable, and shifting the floating point values by one makes the error manifest.

Variables yr, xr could start at 0 and be incremented by dy, dx respectively, at the end of the corresponding loop. This is the essence of the DDA. Unfortunately this would accumulate the error inherent in computing dy and dx very quickly.

Introduce a variable yn = yr·height so that yr = yn/height. Every time yr is incremented by dy = (ny−3)/height, yn is incremented by ny−3 which is an integer. Now yi is the integer part of yr, which is to say the integer part of yn/height. This means that yi will increase whenever yn ≥ height; we can keep track of the quotient of i·(ny−3) by height in yi and let yn be the remainder, adjusting it at each step.

Exactly the same reasoning applies to xr and xi by introducing an integer variable xn representing the remainder of j·(nx−3) when divided by width. Making both changes results in yr, xr being simple constants:

void interp0() {
  int offset = 0;
  int yn = 0;
  int yi = 1;
  for (int i = 0; i < height; i++) {
    final float yr = (float) yn / (float) height;
    int xn = 0;
    int xi = 1;
    for (int j = 0; j < width; j++) {
      final float xr = (float) xn / (float) width;
      float c = 0;
      switch (method) {
      case NEAREST:
        c = data[yi + (int) (yr + 0.5)][xi + (int) (xr + 0.5)];
        break;
      case BILINEAR:
        c = linear(yr,
            linear(xr, data[yi+0][xi+0], data[yi+0][xi+1]),
            linear(xr, data[yi+1][xi+0], data[yi+1][xi+1]) );
        break;
      case BICUBIC:
        c = cubic(yr,
            cubic(xr, data[yi-1][xi-1], data[yi-1][xi+0], data[yi-1][xi+1], data[yi-1][xi+2]),
            cubic(xr, data[yi+0][xi-1], data[yi+0][xi+0], data[yi+0][xi+1], data[yi+0][xi+2]),
            cubic(xr, data[yi+1][xi-1], data[yi+1][xi+0], data[yi+1][xi+1], data[yi+1][xi+2]),
            cubic(xr, data[yi+2][xi-1], data[yi+2][xi+0], data[yi+2][xi+1], data[yi+2][xi+2]) );
        break;
      }
      pixels[offset++] = hue(c);
      xn += nx - 3;
      if (xn >= width) {
        xn -= width;
        xi += 1;
      }
    }
    yn += ny - 3;
    if (yn >= height) {
      yn -= height;
      yi += 1;
    }
  }
}

Next time I'll show how to put the array accesses on a diet.