research/2008-displacement/main.c

/* test shit */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>

#include <math.h>

#ifdef BYTECODE
#   include <seccomp-loader.h>
#else
#   include <SDL_image.h>
#endif

#define MAXWIDTH 512
#define MAXHEIGHT 512

#define true 1
#define false 0

static int WIDTH, HEIGHT;

int main(int, char *[]);

#ifdef BYTECODE
#   define MAXIMAGES 6
static int slots[MAXIMAGES];
static double slotbuf[MAXIMAGES * MAXWIDTH * MAXHEIGHT];

volatile int stop;

void sighandler(int signal)
{
    if(sys_write(1, "done", 4) != 4)
        sys_exit(3);
    stop = 1;
}

void bytecode(unsigned char * mem, int heap_size, int stack_size)
{
    char args[10][1024];
    char *argv[] = { args[0], args[1], args[2], args[3], args[4],
                     args[5], args[6], args[7], args[8], args[9] };
    int i, argc;
    char c;

    if(sys_read(0, &c, 1) != 1)
        sys_exit(-5);
    argc = (unsigned char)c;

    for(i = 0; i < argc; i++)
    {
        char *p = argv[i];
        do
            if(sys_read(0, p, 1) != 1)
                sys_exit(-5);
        while(*p++);
    }

    main(argc, argv);

    c = 0;
    if(sys_write(1, &c, 1) != 1)
        sys_exit(-6);
}

static int myatoi(const char *str)
{
    int i = 0;
    while(*str)
    {
        i *= 10;
        i += (unsigned char)*str++ - (unsigned char)'0';
    }
    return i;
}
#define atoi myatoi
#endif

#define WRITE_INT(i, base) \
    do \
    { \
        char buf[128], *b = buf + 127; \
        if(i <= 0) \
            sys_write(1, (i = -i) ? "-" : "0", 1); /* XXX: hack here */ \
        while(i) \
        { \
            *b-- = hex2char[i % base]; \
            i /= base; \
        } \
        sys_write(1, b + 1, (int)(buf + 127 - b)); \
    } while(0)

static void out(FILE *stream, const char *f, va_list args)
{
#ifdef BYTECODE
    static char const *hex2char = "0123456789abcdef";

    for( ; *f; f++)
    {
        if(*f != '%')
        {
            sys_write(1, f, 1);
            continue;
        }

        f++;
        if(!*f)
            break;

        if(*f == 'c')
        {
            char i = (char)(unsigned char)va_arg(args, int);
            if(i >= 0x20 && i < 0x7f)
                sys_write(1, &i, 1);
            else if(i == '\n')
                sys_write(1, "\\n", 2);
            else if(i == '\t')
                sys_write(1, "\\t", 2);
            else if(i == '\r')
                sys_write(1, "\\r", 2);
            else
            {
                sys_write(1, "\\x", 2);
                sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
                sys_write(1, hex2char + (i & 0x0f), 1);
            }
        }
        else if(*f == 'i' || *f == 'd')
        {
            int i = va_arg(args, int);
            WRITE_INT(i, 10);
        }
        else if(*f == 'x')
        {
            int i = va_arg(args, int);
            WRITE_INT(i, 16);
        }
        else if(f[0] == 'l' && (f[1] == 'i' || f[1] == 'd'))
        {
            long int i = va_arg(args, long int);
            WRITE_INT(i, 10);
            f++;
        }
        else if(f[0] == 'l' && f[1] == 'l' && (f[2] == 'i' || f[1] == 'd'))
        {
            long long int i = va_arg(args, long long int);
            WRITE_INT(i, 10);
            f += 2;
        }
        else if(f[0] == 'g')
        {
            double g = va_arg(args, double), h = 0.0000001;
            int i = g;
            WRITE_INT(i, 10);
            for(i = 0; i < 7; i++)
            {
                g = (g - (int)g) * 10;
                h *= 10;
                if(g < h)
                    break;
                if(i == 0)
                    sys_write(1, ".", 1);
                sys_write(1, hex2char + (int)g, 1);
            }
        }
        else if(f[0] == 'p')
        {
            uintptr_t i = va_arg(args, uintptr_t);
            if(!i)
                sys_write(1, "NULL", 5);
            else
            {
                sys_write(1, "0x", 2);
                WRITE_INT(i, 16);
            }
        }
        else if(f[0] == 's')
        {
            char *s = va_arg(args, char *);
            if(!s)
                sys_write(1, "(nil)", 5);
            else
            {
                int l = 0;
                while(s[l])
                    l++;
                sys_write(1, s, l);
            }
        }
        else if(f[0] == '0' && f[1] == '2' && f[2] == 'x')
        {
            int i = va_arg(args, int);
            sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
            sys_write(1, hex2char + (i & 0x0f), 1);
            f += 2;
        }
        else
        {
            sys_write(1, f - 1, 2);
        }
    }
#else
    vfprintf(stream, f, args);
    fflush(stream);
#endif
}

static void err(const char *f, ...)
{
    va_list args;
    va_start(args, f);
    out(stderr, f, args);
    va_end(args);
}

static void msg(const char *f, ...)
{
    va_list args;
    va_start(args, f);
    out(stdout, f, args);
    va_end(args);
}

static inline double get(double const *src, int x, int y)
{
    return src[y * WIDTH + x];
}

static inline void put(double *src, int x, int y, double p)
{
    src[y * WIDTH + x] = p;
}

static double *new(void)
{
#ifdef BYTECODE
    int i;
    for(i = 0; i < MAXIMAGES; i++)
        if(slots[i] == 0)
            break;
    if(i == MAXIMAGES)
        return NULL;
    slots[i] = 1;
    return slotbuf + i * MAXWIDTH * MAXHEIGHT;
#else
    return malloc(WIDTH * HEIGHT * sizeof(double));
#endif
}

static void del(double *img)
{
#ifdef BYTECODE
    int i;
    for(i = 0; i < MAXIMAGES; i++)
        if(slotbuf + i * MAXWIDTH * MAXHEIGHT == img)
            break;
    if(i == MAXIMAGES)
        return;
    slots[i] = 0;
#else
    free(img);
#endif
}

static double *copy(double const *src)
{
    double *dest = new();
    memcpy(dest, src, WIDTH * HEIGHT * sizeof(double)); 
    return dest;
}

#define N 5
#define NN ((N * 2 + 1))

static void makegauss(double mat[NN][NN], double sigma, double dx, double dy)
{
    double t = 0;
    int i, j;

    sigma = 2. * sigma * sigma;

    for(j = 0; j < NN; j++)
        for(i = 0; i < NN; i++)
        {
            double a = (double)(i - N) + dx;
            double b = (double)(j - N) + dy;
            mat[i][j] = pow(M_E, - (a * a + b * b) / sigma);
            t += mat[i][j];
        }

    for(j = 0; j < NN; j++)
        for(i = 0; i < NN; i++)
            mat[i][j] /= t;
}

static double *gauss(double const *src, double mat[NN][NN])
{
    double *dest = new();
    int x, y, i, j;

    for(y = N; y < HEIGHT - N; y++)
        for(x = N; x < WIDTH - N; x++)
        {
            double p = 0;

            for(j = 0; j < NN; j++)
                for(i = 0; i < NN; i++)
                    p += get(src, x + i - N, y + j - N) * mat[i][j];

            put(dest, x, y, p);
        }

    return dest;
}

static double *fullgauss(double const *src, double mat[NN][NN])
{
    double *dest = new();
    int x, y, i, j;

    for(y = 0; y < HEIGHT; y++)
        for(x = 0; x < WIDTH; x++)
        {
            double p = 0;

            for(j = 0; j < NN; j++)
                for(i = 0; i < NN; i++)
                    if(x + i >= N && x + i < WIDTH + N
                       && y + j >= N && y + j < HEIGHT + N)
                        p += get(src, x + i - N, y + j - N) * mat[i][j];

            put(dest, x, y, p);
        }

    return dest;
}

#if 0
static double fulldist(double const *p1, const double *p2)
{
    double error = 0.;
    int x, y;

    for(y = 0; y < HEIGHT; y++)
        for(x = 0; x < WIDTH; x++)
        {
            double t = get(p1, x, y) - get(p2, x, y);
            error += t * t;
        }

    return error / (WIDTH * HEIGHT);
}
#endif

static double dist(double const *p1, double const *p2, double max)
{
    double error = 0.;
    int x, y;

    max *= (WIDTH - N) * (HEIGHT - N);

    for(y = N; y < HEIGHT - N; y++)
    {
        for(x = N; x < WIDTH - N; x++)
        {
            double t = get(p1, x, y) - get(p2, x, y);
            error += t * t;
        }

        /* Experience shows that this is almost useless, because the
         * functions we manipulate are so small that the difference
         * only happens in the last 1% of the image... */
        if(error > max)
            break;
    }

    return error / ((WIDTH - N) * (HEIGHT - N));
}

static double *load(char const *name)
{
    double *floats;
    int w, h, x, y;

#ifdef BYTECODE
    char c;

    if(sys_read(0, &c, 1) != 1)
        sys_exit(-5);
    w = ((int)(unsigned char)c) << 8;
    if(sys_read(0, &c, 1) != 1)
        sys_exit(-5);
    w |= (int)(unsigned char)c;

    if(sys_read(0, &c, 1) != 1)
        sys_exit(-5);
    h = ((int)(unsigned char)c) << 8;
    if(sys_read(0, &c, 1) != 1)
        sys_exit(-5);
    h |= (int)(unsigned char)c;
#else
    SDL_Surface *tmp, *surface;
    uint32_t *pixels;

    tmp = IMG_Load(name);
    if(!tmp)
        return NULL;

    w = tmp->w;
    h = tmp->h;
#endif

    WIDTH = w > MAXWIDTH ? MAXWIDTH : w;
    HEIGHT = h > MAXHEIGHT ? MAXHEIGHT : h;

    floats = new();
    if(!floats)
        return NULL;

#ifdef BYTECODE
    for(y = 0; y < h; y++)
    for(x = 0; x < w; x++)
    {
        if(sys_read(0, &c, 1) != 1)
            sys_exit(-5);
        if(x >= WIDTH || y >= HEIGHT)
            continue;
        put(floats, x, y, (double)(unsigned char)c / 0xff);
    }
#else
    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
                                   0xff0000, 0xff00, 0xff, 0x0);
    pixels = (uint32_t *)surface->pixels;
    SDL_BlitSurface(tmp, NULL, surface, NULL);
    SDL_FreeSurface(tmp);

    for(y = 0; y < HEIGHT; y++)
    for(x = 0; x < WIDTH; x++)
    {
        int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
        put(floats, x, y, (double)green / 0xff);
    }
#endif

    return floats;
}

static void save(double const *src, char const *name)
{
    int x, y;
#ifdef BYTECODE
    for(y = 0; y < HEIGHT; y++)
    for(x = 0; x < WIDTH; x++)
    {
        double p = 255 * get(src, x, y);
        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
        char c = (char)(unsigned char)i;
        if(sys_write(1, &c, 1) != 1)
            sys_exit(-6);
    }
#else
    SDL_Surface *surface;
    uint32_t *pixels;

    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
                                   0xff0000, 0xff00, 0xff, 0x0);
    pixels = (uint32_t *)surface->pixels;

    for(y = 0; y < HEIGHT; y++)
    for(x = 0; x < WIDTH; x++)
    {
        double p = 255 * get(src, x, y);
        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
        pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i;
    }

    SDL_SaveBMP(surface, name);
#endif
}

static double *ostromoukhov(double const *src)
{
    static int const table[][3] =
    {
         {13, 0, 5}, {13, 0, 5}, {21, 0, 10}, {7, 0, 4},
         {8, 0, 5}, {47, 3, 28}, {23, 3, 13}, {15, 3, 8},
         {22, 6, 11}, {43, 15, 20}, {7, 3, 3}, {501, 224, 211},
         {249, 116, 103}, {165, 80, 67}, {123, 62, 49}, {489, 256, 191},
         {81, 44, 31}, {483, 272, 181}, {60, 35, 22}, {53, 32, 19},
         {237, 148, 83}, {471, 304, 161}, {3, 2, 1}, {481, 314, 185},
         {354, 226, 155}, {1389, 866, 685}, {227, 138, 125}, {267, 158, 163},
         {327, 188, 220}, {61, 34, 45}, {627, 338, 505}, {1227, 638, 1075},
         {20, 10, 19}, {1937, 1000, 1767}, {977, 520, 855}, {657, 360, 551},
         {71, 40, 57}, {2005, 1160, 1539}, {337, 200, 247}, {2039, 1240, 1425},
         {257, 160, 171}, {691, 440, 437}, {1045, 680, 627}, {301, 200, 171},
         {177, 120, 95}, {2141, 1480, 1083}, {1079, 760, 513}, {725, 520, 323},
         {137, 100, 57}, {2209, 1640, 855}, {53, 40, 19}, {2243, 1720, 741},
         {565, 440, 171}, {759, 600, 209}, {1147, 920, 285}, {2311, 1880, 513},
         {97, 80, 19}, {335, 280, 57}, {1181, 1000, 171}, {793, 680, 95},
         {599, 520, 57}, {2413, 2120, 171}, {405, 360, 19}, {2447, 2200, 57},
         {11, 10, 0}, {158, 151, 3}, {178, 179, 7}, {1030, 1091, 63},
         {248, 277, 21}, {318, 375, 35}, {458, 571, 63}, {878, 1159, 147},
         {5, 7, 1}, {172, 181, 37}, {97, 76, 22}, {72, 41, 17},
         {119, 47, 29}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
         {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
         {4, 1, 1}, {4, 1, 1}, {65, 18, 17}, {95, 29, 26},
         {185, 62, 53}, {30, 11, 9}, {35, 14, 11}, {85, 37, 28},
         {55, 26, 19}, {80, 41, 29}, {155, 86, 59}, {5, 3, 2},
         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
         {305, 176, 119}, {155, 86, 59}, {105, 56, 39}, {80, 41, 29},
         {65, 32, 23}, {55, 26, 19}, {335, 152, 113}, {85, 37, 28},
         {115, 48, 37}, {35, 14, 11}, {355, 136, 109}, {30, 11, 9},
         {365, 128, 107}, {185, 62, 53}, {25, 8, 7}, {95, 29, 26},
         {385, 112, 103}, {65, 18, 17}, {395, 104, 101}, {4, 1, 1}
    };

    double *dest = new();
    double *tmp = copy(src);
    int x, y;

    for(y = 0; y < HEIGHT; y++)
    {
        for(x = 0; x < WIDTH; x++)
        {
            int x1 = (y & 1) ? WIDTH - 1 - x + 1 : x - 1;
            int x2 = (y & 1) ? WIDTH - 1 - x : x;
            int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1;

            double p = get(tmp, x2, y);
            double q = p > 0.5 ? 1. : 0.;
            double error = (p - q);
            int i = p * 255.9999;

            put(dest, x2, y, q);

            if(i > 127)
                i = 255 - i;
            if(i < 0)
                i = 0;

            error /= table[i][0] + table[i][1] + table[i][2];

            if(x < WIDTH - 1)
                put(tmp, x3, y, get(tmp, x3, y) + error * table[i][0]);
            if(y < HEIGHT - 1)
            {
                if(x > 0)
                    put(tmp, x1, y + 1,
                        get(tmp, x1, y + 1) + error * table[i][1]);
                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * table[i][2]);
            }
        }
    }

    del(tmp);

    return dest;
}

/* Dither using error diffusion and Floyd-Steinberg-like coefficients:
         X a b
     c d e f g
     h i j k l
*/
static double *ed(double const *src, int serpentine,
                 int a, int b, int c, int d, int e, int f,
                 int g, int h, int i, int j, int k, int l)
{
    double *dest = new();
    double *tmp = copy(src);
    int x, y, n;

    n = a + b + c + d + e + f + g + h + i + j + k + l;

    for(y = 0; y < HEIGHT; y++)
    {
        int swap = serpentine && (y & 1);

        for(x = 0; x < WIDTH; x++)
        {
            int x0 = swap ? WIDTH - 1 - x + 2 : x - 2;
            int x1 = swap ? WIDTH - 1 - x + 1 : x - 1;
            int x2 = swap ? WIDTH - 1 - x     : x;
            int x3 = swap ? WIDTH - 1 - x - 1 : x + 1;
            int x4 = swap ? WIDTH - 1 - x - 2 : x + 2;

            double p = get(tmp, x2, y);
            double q = p > 0.5 ? 1. : 0.;
            double error = (p - q) / n;

            put(dest, x2, y, q);

            if(x < WIDTH - 1)
                put(tmp, x3, y, get(tmp, x3, y) + error * a);
            if(x < WIDTH - 2)
                put(tmp, x4, y, get(tmp, x4, y) + error * b);
            if(y < HEIGHT - 1)
            {
                if(x > 0)
                {
                    put(tmp, x1, y + 1,
                        get(tmp, x1, y + 1) + error * d);
                    if(x > 1)
                        put(tmp, x0, y + 1,
                            get(tmp, x0, y + 1) + error * c);
                }
                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * e);
                if(x < WIDTH - 1)
                {
                    put(tmp, x3, y + 1,
                        get(tmp, x3, y + 1) + error * f);
                    if(x < WIDTH - 2)
                        put(tmp, x4, y + 1,
                            get(tmp, x4, y + 1) + error * g);
                }
            }
            if(y < HEIGHT - 2)
            {
                if(x > 0)
                {
                    put(tmp, x1, y + 2,
                        get(tmp, x1, y + 2) + error * h);
                    if(x > 1)
                        put(tmp, x0, y + 2,
                            get(tmp, x0, y + 2) + error * i);
                }
                put(tmp, x2, y + 2, get(tmp, x2, y + 2) + error * j);
                if(x < WIDTH - 1)
                {
                    put(tmp, x3, y + 2,
                        get(tmp, x3, y + 2) + error * k);
                    if(x < WIDTH - 2)
                        put(tmp, x4, y + 2,
                            get(tmp, x4, y + 2) + error * l);
                }
            }
        }
    }

    del(tmp);

    return dest;
}

/* XXX */
static double *dbs(double const *src, double const *orig,
                  double sigma, double dx, double dy)
{
    double mat[NN][NN];
    double *dest, *tmp, *tmp2;
    double error;

    makegauss(mat, sigma, 0., 0.);
    tmp = fullgauss(src, mat);

    makegauss(mat, sigma, dx, dy);
    dest = copy(orig);
    tmp2 = fullgauss(dest, mat);

    error = dist(tmp, tmp2, 1.);

    for(;;)
    {
        int changes = 0;
        int x, y, i, j, n;

        for(y = 0; y < HEIGHT; y++)
        for(x = 0; x < WIDTH; x++)
        {
            double d, d2, e, best = 0.;
            int opx = -1, opy = -1;

            d = get(dest, x, y);

            /* Compute the effect of a toggle */
            e = 0.;
            for(j = -N; j < N + 1; j++)
            {
                if(y + j < 0 || y + j >= HEIGHT)
                    continue;

                for(i = -N; i < N + 1; i++)
                {
                    double m, p, q1, q2;

                    if(x + i < 0 || x + i >= WIDTH)
                        continue;
                    
                    m = mat[i + N][j + N];
                    p = get(tmp, x + i, y + j);
                    q1 = get(tmp2, x + i, y + j);
                    q2 = q1 - m * d + m * (1. - d);
                    e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
                }
            }
            if(e > best)
            {
                best = e;
                opx = opy = 0;
            }

            /* Compute the effect of swaps */
            for(n = 0; n < 8; n++)
            {
                static int const step[] =
                    { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
                int idx = step[n * 2], idy = step[n * 2 + 1];
                if(y + idy < 0 || y + idy >= HEIGHT
                     || x + idx < 0 || x + idx >= WIDTH)
                    continue;
                d2 = get(dest, x + idx, y + idy);
                if(d2 == d)
                    continue;
                e = 0.;
                for(j = -N; j < N + 1; j++)
                {
                    if(y + j < 0 || y + j >= HEIGHT)
                        continue;
                    if(j - idy + N < 0 || j - idy + N >= NN)
                        continue;
                    for(i = -N; i < N + 1; i++)
                    {
                        double ma, mb, p, q1, q2;
                        if(x + i < 0 || x + i >= WIDTH)
                            continue;
                        if(i - idx + N < 0 || i - idx + N >= NN)
                            continue;
                        ma = mat[i + N][j + N];
                        mb = mat[i - idx + N][j - idy + N];
                        p = get(tmp, x + i, y + j);
                        q1 = get(tmp2, x + i, y + j);
                        q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
                    }
                }
                if(e > best)
                {
                    best = e;
                    opx = idx;
                    opy = idy;
                }
            }

            /* Apply the change if interesting */
            if(best <= 0.)
                continue;
            if(opx || opy)
            {
                d2 = get(dest, x + opx, y + opy);
                put(dest, x + opx, y + opy, d);
            }
            else
                d2 = 1. - d;
            put(dest, x, y, d2);
            for(j = -N; j < N + 1; j++)
            for(i = -N; i < N + 1; i++)
            {
                double m = mat[i + N][j + N];
                if(y + j >= 0 && y + j < HEIGHT
                    && x + i >= 0 && x + i < WIDTH)
                {
                    double t = get(tmp2, x + i, y + j);
                    put(tmp2, x + i, y + j, t + m * (d2 - d));
                }
                if((opx || opy) && y + opy + j >= 0 && y + opy + j < HEIGHT
                                && x + opx + i >= 0 && x + opx + i < WIDTH)
                {
                    double t = get(tmp2, x + opx + i, y + opy + j);
                    put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2));
                }
            }

            changes++;
        }

        err("did %i changes\n", changes);

        if(changes == 0)
            break;
    }

    del(tmp);
    del(tmp2);

    return dest;
}

static void study(double const *src, double const *dest,
                  double sigma, double precision, double fdx, double fdy)
{
#   define Z 3
    double mat[NN][NN], mat2[NN][NN];
    double *tmp, *tmp2;
    double e, e0, e1;
    double best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
    int dx, dy;

    makegauss(mat, sigma, 0., 0.);
    tmp = gauss(src, mat);

    tmp2 = gauss(dest, mat);
    e0 = dist(tmp, tmp2, 1.);
    del(tmp2);

    makegauss(mat2, sigma, fdx, fdy);
    tmp2 = gauss(dest, mat2);
    e1 = dist(tmp, tmp2, 1.);
    del(tmp2);

    while(step > precision)
    {
        for(dy = 0; dy <= Z; dy++)
            for(dx = 0; dx <= Z; dx++)
            {
                makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z);
                tmp2 = gauss(dest, mat);
                e = dist(tmp, tmp2, best);
                del(tmp2);
                if(e < best)
                {
                    best = e;
                    bfx = fx + step * dx / Z;
                    bfy = fy + step * dy / Z;
                }
            }

        fx = bfx - step / Z;
        fy = bfy - step / Z;
        step = step * 2 / Z;
    }

    del(tmp);

    msg("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n",
        1000. * e0, 1000. * e1, 1000. * best, fx, fy);
}

static double *merge(double const *im1, double const *im2, double t)
{
    double *dest = new();
    int x, y;

    for(y = 0; y < HEIGHT; y++)
    for(x = 0; x < WIDTH; x++)
        put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y));

    return dest;
}

static void usage(char *argv[])
{
    msg("Usage: %s <mode> [ARGS...]\n", argv[0]);
    msg("Allowed modes:\n");
    msg(" -1 <src>           raster FS displacement study on src\n");
    msg(" -2 <src1> <src2>   raster FS displacement study on blends of src1 and src2\n");
    msg(" -3 <src>           quick (a,b,c,d) ED kernel analysis on src\n");
    msg(" -4 <src>           exhaustive (a,b,c,d) ED kernel analysis on src\n");
    msg(" -5 <src>           exhaustive displacement study on src\n");
    msg(" -6 <src>           restrained (a,b,c,d) ED kernel analysis on src\n");
    msg(" -7 <src>           restrained displacement study on src\n");
    msg(" -8 <src>           displacement values on src\n");
}

int main(int argc, char *argv[])
{
    double *src;
    int mode, i;

    if(argc < 2)
    {
        err("%s: too few arguments\n", argv[0]);
        usage(argv);
        return EXIT_FAILURE;
    }

    if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1)))
    {
        err("%s: invalid mode `%s'\n", argv[0], argv[1]);
        usage(argv);
        return EXIT_FAILURE;
    }

    src = load(argv[2]);
    if(!src)
        return 2;

    msg("### mode %i on `%s' ###\n", mode, argv[2]);

    switch(mode)
    {
        case 1:
        {
            double *dest = ed(src, false, 7, 0,
                                 0, 3, 5, 1, 0,
                                 0, 0, 0, 0, 0);
            study(src, dest, 1.2, 0.001, .16, .28);
            del(dest);
            del(src);
        }
        break;

        case 2:
        {
            double *src2, *dest, *tmp;

            if(argc < 3)
                return 3;
            src2 = load(argv[2]);
            if(!src2)
                return 4;

            for(i = 0; i <= 100; i++)
            {
                tmp = merge(src, src2, (double)i / 100.);
                dest = ed(tmp, false, 7, 0,
                             0, 3, 5, 1, 0,
                             0, 0, 0, 0, 0);
                study(tmp, dest, 1.2, 0.001, .16, .28);
                del(dest);
                del(tmp);
            }
            del(src2);
            del(src);
        }
        break;

        case 3:
        case 4:
        {
            double mat[NN][NN];
            double *dest, *tmp, *tmp2;
            int a, b, c, d, e;

#define TOTAL 16
            makegauss(mat, 1.2, 0, 0);
            for(a = 0; a <= TOTAL; a++)
            for(b = 0; b <= TOTAL; b++)
            for(c = 0; c <= TOTAL; c++)
            for(d = 0; d <= TOTAL; d++)
            for(e = 0; e <= TOTAL; e++)
            {
                int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e, i;

                if(a + b + c + d + e != TOTAL)
                    continue;

                /* Slightly shuffle our coefficients to avoid waiting until
                 * 75% progress before having an idea of what's going on. */
#define SHUFFLE(p,q,n) \
    if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; }
                SHUFFLE(c2, d2, 777); SHUFFLE(b2, d2, 555);
                SHUFFLE(a2, d2, 333); SHUFFLE(b2, c2, 222);
                SHUFFLE(a2, e2, 333); SHUFFLE(b2, e2, 222);
                SHUFFLE(a2, c2, 444); SHUFFLE(a2, b2, 666);
                SHUFFLE(c2, e2, 999); SHUFFLE(d2, e2, 777);
                SHUFFLE(a2, d2, 999); SHUFFLE(a2, b2, 777);

#if 0
                if(b2 > c2) continue;
                if(d2 > c2) continue;
                if(d2 > a2) continue;
#endif
                /* We only want 4-cell kernels for now */
                if(b2) continue;

                for(i = 1; i <= 2; i++)
                {
                    //msg("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2);
                    msg("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2);

                    dest = ed(src, i == 2, a2, 0,
                               b2, c2, d2, e2, 0,
                                0,  0,  0,  0, 0);
                    if(mode == 4)
                    {
                        /* XXX: E_fast is meaningless, ignore it */
                        study(src, dest, 1.2, 0.001, 0., 0.);
                    }
                    else
                    {
                        tmp = gauss(src, mat);
                        tmp2 = gauss(dest, mat);
                        msg("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.));
                        del(tmp);
                        del(tmp2);
                    }
                    fflush(stdout);
                    del(dest);
                }
            }

            del(src);
        }
        break;

        case 5:
        {
            double *dest;

            dest = ed(src, false, 7, 0,
                         0, 3, 5, 1, 0,
                         0, 0, 0, 0, 0);
            msg("[1] ");
            study(src, dest, 1.2, 0.001, .16, .28);
            del(dest);

            dest = ed(src, false, 7, 5,
                         3, 5, 7, 5, 3,
                         1, 3, 5, 3, 1);
            msg("[2] ");
            study(src, dest, 1.2, 0.001, .26, .76);
            del(dest);

            dest = ostromoukhov(src);
            msg("[3] ");
            study(src, dest, 1.2, 0.001, .0, .19);
            del(dest);

            dest = ed(src, true, 2911, 0,
                  0, 1373, 3457, 2258, 0,
                  0,    0,    0,    0, 0);
            msg("[4] ");
            study(src, dest, 1.2, 0.001, .0, .34);
        }
        break;

        case 6:
        case 7:
        {
            double mat[NN][NN];
            double *dest;
            int a, ax, b, bx, c, cx, d, dx;

            if(mode == 6)
                a = 7, b = 3, c = 5, d = 1;
            else
                a = 7, b = 5, c = 4, d = 0;

#undef TOTAL
#define TOTAL (a+b+c+d)
            makegauss(mat, 1.2, 0, 0);
            for(ax = -2; ax <= 2; ax++)
            for(bx = -2; bx <= 2; bx++)
            for(cx = -2; cx <= 2; cx++)
            for(dx = -2; dx <= 3; dx++)
            {
                int a2 = a + ax, b2 = b + bx, c2 = c + cx, d2 = d + dx;

                if(a2 < 0 || a2 > TOTAL || b2 < 0 || b2 > TOTAL ||
                     c2 < 0 || c2 > TOTAL || d2 < 0 || d2 > TOTAL)
                    continue;

                if(a2 + b2 + c2 + d2 != TOTAL)
                    continue;

                msg("K: %d,%d,%d,%d ", a2, b2, c2, d2);

                dest = ed(src, mode == 7, a2, 0,
                               0, b2, c2, d2, 0,
                               0,  0,  0,  0, 0);
                /* XXX: E_fast is meaningless, ignore it */
                study(src, dest, 1.2, 0.001, 0., 0.);
                fflush(stdout);
                del(dest);
            }

            del(src);
        }
        break;

        case 8:
        {
            const int STEP = 32;
            double mat[NN][NN];
            double *dest = ed(src, false, 7, 0,
                                 0, 3, 5, 1, 0,
                                 0, 0, 0, 0, 0);
            double *tmp, *tmp2;
            int dx, dy;

            makegauss(mat, 1.2, 0., 0.);
            tmp = gauss(src, mat);
            for(dy = 0; dy <= STEP; dy++)
            {
                for(dx = 0; dx <= STEP; dx++)
                {
                    double fy = 2. / STEP * (dy - STEP / 2.);
                    double fx = 2. / STEP * (dx - STEP / 2.);

                    makegauss(mat, 1.2, fx, fy);
                    tmp2 = gauss(dest, mat);
                    msg("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
                    del(tmp2);
                }
                msg("\n");
            }
            del(tmp);
            del(dest);
            del(src);
        }
        break;

#if 0
    tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    //dest = dbs(src, tmp, 0., 0.);
    dest = dbs(src, tmp, 0.20, 0.30);
    //dest = dbs(src, tmp, 0.158718, 0.283089);
    //dest = copy(tmp);
    del(tmp);
    study(src, dest, 1.2, 0.00001);
    save(dest, "output.bmp");
    del(dest);
#endif

#if 0
    //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
    //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0);
    //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0);
    dest = ostromoukhov(src);
    //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
    //msg("%s: ", argv[1]);
    //study(src, dest, 1.2, 0.001);
    save(dest, "output.bmp");
    del(dest);
#endif

#if 0

    save(dest, "output.bmp");
#endif

#if 0
    tmp = gauss(src, 0., 0.);
    for(a = 0; a < 16; a++)
    for(b = 0; b < 16; b++)
    for(c = 0; c < 16; c++)
    for(d = 0; d < 16; d++)
    {
        if(a + b + c + d != 16)
            continue;
        dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
        tmp2 = gauss(dest, 0., 0.);
        msg("%.5g: (%02d %02d %02d %02d)\n",
            1000. * dist(tmp, tmp2, 1.), a, b, c, d);
        del(dest);
        del(tmp2);
    }

    save(dest, "output.bmp");
#endif

#if 0
    msg("%s\n", argv[1]);

    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    study(src, dest, 1.2, 0.01);
    del(dest);

    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
    study(src, dest, 1.2, 0.01);
    del(dest);

    dest = ostromoukhov(src);
    study(src, dest, 1.2, 0.01);
    del(dest);

    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
    study(src, dest, 1.2, 0.01);
    del(dest);

    msg("\n");
#endif

#if 0
    //dest = ostromoukhov(src);
    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    dest = dbs(src, tmp, 0., 0.);
    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
    //for(sigma = 0.8; sigma < 2.; sigma *= 1.01)
    {
        msg("%g ", sigma);
        study(src, dest, sigma, 0.01);
    }
#endif

    }

#if 0
    tmp = new();
    a = 0;
    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
    {
        char buf[1024];
        msg("%i: %g\n", a, sigma);
        dest = dbs(src, tmp, sigma, 0., 0.);
        sprintf(buf, "output-dbs-%i.bmp", a++);
        save(dest, buf);
    }
#endif

    return 0;
}