Here's the source I used, written by myself. The choice of 16384x16384 might have been overkill; I cannot spot any significant differences between this and using 1024x1024 taps all over the image. The generation took about an hour on one core of an Intel Q9450.
#include <stdio.h> #include <stdlib.h> #include <math.h> /* * Do a finer sampling the higher up we get in the picture, * since there's much more detail there. */ #define SAMPLES_FROM_Y(y) ((16384 * 512) / (y*y + 512)) #define MAX_SAMPLES SAMPLES_FROM_Y(0) float weights[MAX_SAMPLES][MAX_SAMPLES]; static inline int color(double x, double y) { double t, z; int i, j, k; x = x * (1.0 / 128.0) - 0.5; y = y * (1.0 / 2048.0); t = 1.0 / (y + 0.001); z = t * x; i = floor(t); j = floor(z); k = i + j; return ((k % 2) != 0); } static double sinc(double x) { static const double cutoff = 1.220703668e-4; /* sqrt(sqrt(eps)) */ if (abs(x) < cutoff) { /* For small |x|, use Taylor series instead */ const double x2 = x * x; const double x4 = x2 * x2; return 1.0 - x2 / 6.0 + x4 / 120.0; } else { return sin(x) / x; } } static double lanczos_tap(double x) { if (x < -3.0 || x > 3.0) return 0.0; if (x < 0.0) return sinc(-x*M_PI) * sinc(-x*M_PI / 3.0); else return sinc(x*M_PI) * sinc(x*M_PI / 3.0); } static void precalculate_weights(unsigned num_samples) { double total = 0.0; for (int yi = 0; yi < num_samples; ++yi) { double sy = 6.0 * ((double)yi / num_samples - 0.5); double wy = lanczos_tap(sy); for (int xi = 0; xi < num_samples; ++xi) { double sx = 6.0 * ((double)xi / num_samples - 0.5); weights[yi][xi] = wy * lanczos_tap(sx); total += weights[yi][xi]; } } double inv_total = 1.0 / total; for (int yi = 0; yi < num_samples; ++yi) { for (int xi = 0; xi < num_samples; ++xi) { weights[yi][xi] *= inv_total; } } } static double antialiased_color(double x, double y, unsigned num_samples) { double acc = 0.0; double delta = 6.0 / num_samples; int xi, yi; double sx, sy; for (yi = 0, sy = y - 3.0; yi < num_samples; ++yi, sy += delta) { for (xi = 0, sx = x - 3.0; xi < num_samples; ++xi, sx += delta) { if (color(sx, sy)) { acc += weights[yi][xi]; } } } return acc; }