source: research/2008-displacement/main.c @ 2286

Last change on this file since 2286 was 2286, checked in by Sam Hocevar, 14 years ago
  • More scripts for part 3.
  • Implemented main.c as a seccomp bytecode for CPUShare.
File size: 31.5 KB
Line 
1/* test shit */
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
6#include <stdarg.h>
7
8#include <math.h>
9
10#ifdef BYTECODE
11#   include <seccomp-loader.h>
12#else
13#   include <SDL_image.h>
14#endif
15
16#define MAXWIDTH 512
17#define MAXHEIGHT 512
18
19#define true 1
20#define false 0
21
22static int WIDTH, HEIGHT;
23
24int main(int, char *[]);
25
26#ifdef BYTECODE
27#   define MAXIMAGES 128
28static int slots[MAXIMAGES];
29static double slotbuf[MAXIMAGES * MAXWIDTH * MAXHEIGHT];
30
31volatile int stop;
32
33void sighandler(int signal)
34{
35    if(sys_write(1, "done", 4) != 4)
36        sys_exit(3);
37    stop = 1;
38}
39
40void bytecode(unsigned char * mem, int heap_size, int stack_size)
41{
42    char mode[] = "-0";
43    char *argv[] = { "program", mode, "arg" };
44    char c;
45
46    if(sys_read(0, mode + 1, 1) != 1)
47        sys_exit(-5);
48
49    main(sizeof(argv)/sizeof(*argv), argv);
50
51    c = 0;
52    if(sys_write(1, &c, 1) != 1)
53        sys_exit(-6);
54}
55
56static int myatoi(const char *str)
57{
58    int i = 0;
59    while(*str)
60    {
61        i *= 10;
62        i += (unsigned char)*str++ - (unsigned char)'0';
63    }
64    return i;
65}
66#define atoi myatoi
67#endif
68
69#define WRITE_INT(i, base) \
70    do \
71    { \
72        char buf[128], *b = buf + 127; \
73        if(i <= 0) \
74            sys_write(1, (i = -i) ? "-" : "0", 1); /* XXX: hack here */ \
75        while(i) \
76        { \
77            *b-- = hex2char[i % base]; \
78            i /= base; \
79        } \
80        sys_write(1, b + 1, (int)(buf + 127 - b)); \
81    } while(0)
82
83static void msg(const char *f, ...)
84{
85    va_list args;
86    va_start(args, f);
87#ifdef BYTECODE
88    static char const *hex2char = "0123456789abcdef";
89
90    for( ; *f; f++)
91    {
92        if(*f != '%')
93        {
94            sys_write(1, f, 1);
95            continue;
96        }
97
98        f++;
99        if(!*f)
100            break;
101
102        if(*f == 'c')
103        {
104            char i = (char)(unsigned char)va_arg(args, int);
105            if(i >= 0x20 && i < 0x7f)
106                sys_write(1, &i, 1);
107            else if(i == '\n')
108                sys_write(1, "\\n", 2);
109            else if(i == '\t')
110                sys_write(1, "\\t", 2);
111            else if(i == '\r')
112                sys_write(1, "\\r", 2);
113            else
114            {
115                sys_write(1, "\\x", 2);
116                sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
117                sys_write(1, hex2char + (i & 0x0f), 1);
118            }
119        }
120        else if(*f == 'i' || *f == 'd')
121        {
122            int i = va_arg(args, int);
123            WRITE_INT(i, 10);
124        }
125        else if(*f == 'x')
126        {
127            int i = va_arg(args, int);
128            WRITE_INT(i, 16);
129        }
130        else if(f[0] == 'l' && (f[1] == 'i' || f[1] == 'd'))
131        {
132            long int i = va_arg(args, long int);
133            WRITE_INT(i, 10);
134            f++;
135        }
136        else if(f[0] == 'l' && f[1] == 'l' && (f[2] == 'i' || f[1] == 'd'))
137        {
138            long long int i = va_arg(args, long long int);
139            WRITE_INT(i, 10);
140            f += 2;
141        }
142        else if(f[0] == 'g')
143        {
144            double g = va_arg(args, double), h = 9.9f;
145            int i;
146            if(g < 0.f)
147            {
148                sys_write(1, "-", 1);
149                g = -g;
150            }
151            for(i = 0; i < 7; i++)
152            {
153                sys_write(1, hex2char + (int)g, 1);
154                if(i == 0)
155                    sys_write(1, ".", 1);
156                g = (g - (int)g) * 10;
157                h = h / 10.f;
158                if(g < h)
159                    break;
160            }
161        }
162        else if(f[0] == 'p')
163        {
164            uintptr_t i = va_arg(args, uintptr_t);
165            if(!i)
166                sys_write(1, "NULL", 5);
167            else
168            {
169                sys_write(1, "0x", 2);
170                WRITE_INT(i, 16);
171            }
172        }
173        else if(f[0] == 's')
174        {
175            char *s = va_arg(args, char *);
176            if(!s)
177                sys_write(1, "(nil)", 5);
178            else
179            {
180                int l = 0;
181                while(s[l])
182                    l++;
183                sys_write(1, s, l);
184            }
185        }
186        else if(f[0] == '0' && f[1] == '2' && f[2] == 'x')
187        {
188            int i = va_arg(args, int);
189            sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
190            sys_write(1, hex2char + (i & 0x0f), 1);
191            f += 2;
192        }
193        else
194        {
195            sys_write(1, f - 1, 2);
196        }
197    }
198#else
199    vprintf(f, args);
200    fflush(stdout);
201#endif
202    va_end(args);
203}
204
205static inline double get(double const *src, int x, int y)
206{
207    return src[y * WIDTH + x];
208}
209
210static inline void put(double *src, int x, int y, double p)
211{
212    src[y * WIDTH + x] = p;
213}
214
215static double *new(void)
216{
217#ifdef BYTECODE
218    int i;
219    for(i = 0; i < MAXIMAGES; i++)
220        if(slots[i] == 0)
221            break;
222    if(i == MAXIMAGES)
223        return NULL;
224    slots[i] = 1;
225    return slotbuf + i * MAXWIDTH * MAXHEIGHT;
226#else
227    return malloc(WIDTH * HEIGHT * sizeof(double));
228#endif
229}
230
231static void del(double *img)
232{
233#ifdef BYTECODE
234    int i;
235    for(i = 0; i < MAXIMAGES; i++)
236        if(slotbuf + i * MAXWIDTH * MAXHEIGHT == img)
237            break;
238    if(i == MAXIMAGES)
239        return;
240    slots[i] = 0;
241#else
242    free(img);
243#endif
244}
245
246static double *copy(double const *src)
247{
248    double *dest = new();
249    memcpy(dest, src, WIDTH * HEIGHT * sizeof(double));
250    return dest;
251}
252
253#define N 5
254#define NN ((N * 2 + 1))
255
256static void makegauss(double mat[NN][NN], double sigma, double dx, double dy)
257{
258    double t = 0;
259    int i, j;
260
261    sigma = 2. * sigma * sigma;
262
263    for(j = 0; j < NN; j++)
264        for(i = 0; i < NN; i++)
265        {
266            double a = (double)(i - N) + dx;
267            double b = (double)(j - N) + dy;
268            mat[i][j] = pow(M_E, - (a * a + b * b) / sigma);
269            t += mat[i][j];
270        }
271
272    for(j = 0; j < NN; j++)
273        for(i = 0; i < NN; i++)
274            mat[i][j] /= t;
275}
276
277static double *gauss(double const *src, double mat[NN][NN])
278{
279    double *dest = new();
280    int x, y, i, j;
281
282    for(y = N; y < HEIGHT - N; y++)
283        for(x = N; x < WIDTH - N; x++)
284        {
285            double p = 0;
286
287            for(j = 0; j < NN; j++)
288                for(i = 0; i < NN; i++)
289                    p += get(src, x + i - N, y + j - N) * mat[i][j];
290
291            put(dest, x, y, p);
292        }
293
294    return dest;
295}
296
297static double *fullgauss(double const *src, double mat[NN][NN])
298{
299    double *dest = new();
300    int x, y, i, j;
301
302    for(y = 0; y < HEIGHT; y++)
303        for(x = 0; x < WIDTH; x++)
304        {
305            double p = 0;
306
307            for(j = 0; j < NN; j++)
308                for(i = 0; i < NN; i++)
309                    if(x + i >= N && x + i < WIDTH + N
310                       && y + j >= N && y + j < HEIGHT + N)
311                        p += get(src, x + i - N, y + j - N) * mat[i][j];
312
313            put(dest, x, y, p);
314        }
315
316    return dest;
317}
318
319#if 0
320static double fulldist(double const *p1, const double *p2)
321{
322    double error = 0.;
323    int x, y;
324
325    for(y = 0; y < HEIGHT; y++)
326        for(x = 0; x < WIDTH; x++)
327        {
328            double t = get(p1, x, y) - get(p2, x, y);
329            error += t * t;
330        }
331
332    return error / (WIDTH * HEIGHT);
333}
334#endif
335
336static double dist(double const *p1, double const *p2, double max)
337{
338    double error = 0.;
339    int x, y;
340
341    max *= (WIDTH - N) * (HEIGHT - N);
342
343    for(y = N; y < HEIGHT - N; y++)
344    {
345        for(x = N; x < WIDTH - N; x++)
346        {
347            double t = get(p1, x, y) - get(p2, x, y);
348            error += t * t;
349        }
350
351        /* Experience shows that this is almost useless, because the
352         * functions we manipulate are so small that the difference
353         * only happens in the last 1% of the image... */
354        if(error > max)
355            break;
356    }
357
358    return error / ((WIDTH - N) * (HEIGHT - N));
359}
360
361static double *load(char const *name)
362{
363    double *floats;
364    int x, y;
365
366#ifdef BYTECODE
367    char c;
368
369    if(sys_read(0, &c, 1) != 1)
370        sys_exit(-5);
371    WIDTH = ((int)(unsigned char)c) << 8;
372    if(sys_read(0, &c, 1) != 1)
373        sys_exit(-5);
374    WIDTH |= (int)(unsigned char)c;
375
376    if(sys_read(0, &c, 1) != 1)
377        sys_exit(-5);
378    HEIGHT = ((int)(unsigned char)c) << 8;
379    if(sys_read(0, &c, 1) != 1)
380        sys_exit(-5);
381    HEIGHT |= (int)(unsigned char)c;
382
383    floats = new();
384    if(!floats)
385        return NULL;
386
387    for(y = 0; y < HEIGHT; y++)
388    for(x = 0; x < WIDTH; x++)
389    {
390        if(sys_read(0, &c, 1) != 1)
391            sys_exit(-5);
392        put(floats, x, y, (double)(unsigned char)c / 0xff);
393    }   
394#else
395    SDL_Surface *tmp, *surface;
396    uint32_t *pixels;
397
398    tmp = IMG_Load(name);
399    if(!tmp)
400        return NULL;
401
402    WIDTH = tmp->w > MAXWIDTH ? MAXWIDTH : tmp->w;
403    HEIGHT = tmp->h > MAXHEIGHT ? MAXHEIGHT : tmp->h;
404    floats = new();
405    if(!floats)
406        return NULL;
407
408    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
409                                   0xff0000, 0xff00, 0xff, 0x0);
410    pixels = (uint32_t *)surface->pixels;
411    SDL_BlitSurface(tmp, NULL, surface, NULL);
412    SDL_FreeSurface(tmp);
413
414    for(y = 0; y < HEIGHT; y++)
415    for(x = 0; x < WIDTH; x++)
416    {
417        int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
418        put(floats, x, y, (double)green / 0xff);
419    }
420#endif
421
422    return floats;
423}
424
425static void save(double const *src, char const *name)
426{
427    int x, y;
428#ifdef BYTECODE
429    for(y = 0; y < HEIGHT; y++)
430    for(x = 0; x < WIDTH; x++)
431    {
432        double p = 255 * get(src, x, y);
433        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
434        char c = (char)(unsigned char)i;
435        if(sys_write(1, &c, 1) != 1)
436            sys_exit(-6);
437    }
438#else
439    SDL_Surface *surface;
440    uint32_t *pixels;
441
442    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
443                                   0xff0000, 0xff00, 0xff, 0x0);
444    pixels = (uint32_t *)surface->pixels;
445
446    for(y = 0; y < HEIGHT; y++)
447    for(x = 0; x < WIDTH; x++)
448    {
449        double p = 255 * get(src, x, y);
450        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
451        pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i;
452    }
453
454    SDL_SaveBMP(surface, name);
455#endif
456}
457
458static double *ostromoukhov(double const *src)
459{
460    static int const table[][3] =
461    {
462         {13, 0, 5}, {13, 0, 5}, {21, 0, 10}, {7, 0, 4},
463         {8, 0, 5}, {47, 3, 28}, {23, 3, 13}, {15, 3, 8},
464         {22, 6, 11}, {43, 15, 20}, {7, 3, 3}, {501, 224, 211},
465         {249, 116, 103}, {165, 80, 67}, {123, 62, 49}, {489, 256, 191},
466         {81, 44, 31}, {483, 272, 181}, {60, 35, 22}, {53, 32, 19},
467         {237, 148, 83}, {471, 304, 161}, {3, 2, 1}, {481, 314, 185},
468         {354, 226, 155}, {1389, 866, 685}, {227, 138, 125}, {267, 158, 163},
469         {327, 188, 220}, {61, 34, 45}, {627, 338, 505}, {1227, 638, 1075},
470         {20, 10, 19}, {1937, 1000, 1767}, {977, 520, 855}, {657, 360, 551},
471         {71, 40, 57}, {2005, 1160, 1539}, {337, 200, 247}, {2039, 1240, 1425},
472         {257, 160, 171}, {691, 440, 437}, {1045, 680, 627}, {301, 200, 171},
473         {177, 120, 95}, {2141, 1480, 1083}, {1079, 760, 513}, {725, 520, 323},
474         {137, 100, 57}, {2209, 1640, 855}, {53, 40, 19}, {2243, 1720, 741},
475         {565, 440, 171}, {759, 600, 209}, {1147, 920, 285}, {2311, 1880, 513},
476         {97, 80, 19}, {335, 280, 57}, {1181, 1000, 171}, {793, 680, 95},
477         {599, 520, 57}, {2413, 2120, 171}, {405, 360, 19}, {2447, 2200, 57},
478         {11, 10, 0}, {158, 151, 3}, {178, 179, 7}, {1030, 1091, 63},
479         {248, 277, 21}, {318, 375, 35}, {458, 571, 63}, {878, 1159, 147},
480         {5, 7, 1}, {172, 181, 37}, {97, 76, 22}, {72, 41, 17},
481         {119, 47, 29}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
482         {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
483         {4, 1, 1}, {4, 1, 1}, {65, 18, 17}, {95, 29, 26},
484         {185, 62, 53}, {30, 11, 9}, {35, 14, 11}, {85, 37, 28},
485         {55, 26, 19}, {80, 41, 29}, {155, 86, 59}, {5, 3, 2},
486         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
487         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
488         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
489         {305, 176, 119}, {155, 86, 59}, {105, 56, 39}, {80, 41, 29},
490         {65, 32, 23}, {55, 26, 19}, {335, 152, 113}, {85, 37, 28},
491         {115, 48, 37}, {35, 14, 11}, {355, 136, 109}, {30, 11, 9},
492         {365, 128, 107}, {185, 62, 53}, {25, 8, 7}, {95, 29, 26},
493         {385, 112, 103}, {65, 18, 17}, {395, 104, 101}, {4, 1, 1}
494    };
495
496    double *dest = new();
497    double *tmp = copy(src);
498    int x, y;
499
500    for(y = 0; y < HEIGHT; y++)
501    {
502        for(x = 0; x < WIDTH; x++)
503        {
504            int x1 = (y & 1) ? WIDTH - 1 - x + 1 : x - 1;
505            int x2 = (y & 1) ? WIDTH - 1 - x : x;
506            int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1;
507
508            double p = get(tmp, x2, y);
509            double q = p > 0.5 ? 1. : 0.;
510            double error = (p - q);
511            int i = p * 255.9999;
512
513            put(dest, x2, y, q);
514
515            if(i > 127)
516                i = 255 - i;
517            if(i < 0)
518                i = 0;
519
520            error /= table[i][0] + table[i][1] + table[i][2];
521
522            if(x < WIDTH - 1)
523                put(tmp, x3, y, get(tmp, x3, y) + error * table[i][0]);
524            if(y < HEIGHT - 1)
525            {
526                if(x > 0)
527                    put(tmp, x1, y + 1,
528                        get(tmp, x1, y + 1) + error * table[i][1]);
529                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * table[i][2]);
530            }
531        }
532    }
533
534    del(tmp);
535
536    return dest;
537}
538
539/* Dither using error diffusion and Floyd-Steinberg-like coefficients:
540         X a b
541     c d e f g
542     h i j k l
543*/
544static double *ed(double const *src, int serpentine,
545                 int a, int b, int c, int d, int e, int f,
546                 int g, int h, int i, int j, int k, int l)
547{
548    double *dest = new();
549    double *tmp = copy(src);
550    int x, y, n;
551
552    n = a + b + c + d + e + f + g + h + i + j + k + l;
553
554    for(y = 0; y < HEIGHT; y++)
555    {
556        int swap = serpentine && (y & 1);
557
558        for(x = 0; x < WIDTH; x++)
559        {
560            int x0 = swap ? WIDTH - 1 - x + 2 : x - 2;
561            int x1 = swap ? WIDTH - 1 - x + 1 : x - 1;
562            int x2 = swap ? WIDTH - 1 - x     : x;
563            int x3 = swap ? WIDTH - 1 - x - 1 : x + 1;
564            int x4 = swap ? WIDTH - 1 - x - 2 : x + 2;
565
566            double p = get(tmp, x2, y);
567            double q = p > 0.5 ? 1. : 0.;
568            double error = (p - q) / n;
569
570            put(dest, x2, y, q);
571
572            if(x < WIDTH - 1)
573                put(tmp, x3, y, get(tmp, x3, y) + error * a);
574            if(x < WIDTH - 2)
575                put(tmp, x4, y, get(tmp, x4, y) + error * b);
576            if(y < HEIGHT - 1)
577            {
578                if(x > 0)
579                {
580                    put(tmp, x1, y + 1,
581                        get(tmp, x1, y + 1) + error * d);
582                    if(x > 1)
583                        put(tmp, x0, y + 1,
584                            get(tmp, x0, y + 1) + error * c);
585                }
586                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * e);
587                if(x < WIDTH - 1)
588                {
589                    put(tmp, x3, y + 1,
590                        get(tmp, x3, y + 1) + error * f);
591                    if(x < WIDTH - 2)
592                        put(tmp, x4, y + 1,
593                            get(tmp, x4, y + 1) + error * g);
594                }
595            }
596            if(y < HEIGHT - 2)
597            {
598                if(x > 0)
599                {
600                    put(tmp, x1, y + 2,
601                        get(tmp, x1, y + 2) + error * h);
602                    if(x > 1)
603                        put(tmp, x0, y + 2,
604                            get(tmp, x0, y + 2) + error * i);
605                }
606                put(tmp, x2, y + 2, get(tmp, x2, y + 2) + error * j);
607                if(x < WIDTH - 1)
608                {
609                    put(tmp, x3, y + 2,
610                        get(tmp, x3, y + 2) + error * k);
611                    if(x < WIDTH - 2)
612                        put(tmp, x4, y + 2,
613                            get(tmp, x4, y + 2) + error * l);
614                }
615            }
616        }
617    }
618
619    del(tmp);
620
621    return dest;
622}
623
624/* XXX */
625static double *dbs(double const *src, double const *orig,
626                  double sigma, double dx, double dy)
627{
628    double mat[NN][NN];
629    double *dest, *tmp, *tmp2;
630    double error;
631
632    makegauss(mat, sigma, 0., 0.);
633    tmp = fullgauss(src, mat);
634
635    makegauss(mat, sigma, dx, dy);
636    dest = copy(orig);
637    tmp2 = fullgauss(dest, mat);
638
639    error = dist(tmp, tmp2, 1.);
640
641    for(;;)
642    {
643        int changes = 0;
644        int x, y, i, j, n;
645
646        for(y = 0; y < HEIGHT; y++)
647        for(x = 0; x < WIDTH; x++)
648        {
649            double d, d2, e, best = 0.;
650            int opx = -1, opy = -1;
651
652            d = get(dest, x, y);
653
654            /* Compute the effect of a toggle */
655            e = 0.;
656            for(j = -N; j < N + 1; j++)
657            {
658                if(y + j < 0 || y + j >= HEIGHT)
659                    continue;
660
661                for(i = -N; i < N + 1; i++)
662                {
663                    double m, p, q1, q2;
664
665                    if(x + i < 0 || x + i >= WIDTH)
666                        continue;
667                   
668                    m = mat[i + N][j + N];
669                    p = get(tmp, x + i, y + j);
670                    q1 = get(tmp2, x + i, y + j);
671                    q2 = q1 - m * d + m * (1. - d);
672                    e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
673                }
674            }
675            if(e > best)
676            {
677                best = e;
678                opx = opy = 0;
679            }
680
681            /* Compute the effect of swaps */
682            for(n = 0; n < 8; n++)
683            {
684                static int const step[] =
685                    { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
686                int idx = step[n * 2], idy = step[n * 2 + 1];
687                if(y + idy < 0 || y + idy >= HEIGHT
688                     || x + idx < 0 || x + idx >= WIDTH)
689                    continue;
690                d2 = get(dest, x + idx, y + idy);
691                if(d2 == d)
692                    continue;
693                e = 0.;
694                for(j = -N; j < N + 1; j++)
695                {
696                    if(y + j < 0 || y + j >= HEIGHT)
697                        continue;
698                    if(j - idy + N < 0 || j - idy + N >= NN)
699                        continue;
700                    for(i = -N; i < N + 1; i++)
701                    {
702                        double ma, mb, p, q1, q2;
703                        if(x + i < 0 || x + i >= WIDTH)
704                            continue;
705                        if(i - idx + N < 0 || i - idx + N >= NN)
706                            continue;
707                        ma = mat[i + N][j + N];
708                        mb = mat[i - idx + N][j - idy + N];
709                        p = get(tmp, x + i, y + j);
710                        q1 = get(tmp2, x + i, y + j);
711                        q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
712                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
713                    }
714                }
715                if(e > best)
716                {
717                    best = e;
718                    opx = idx;
719                    opy = idy;
720                }
721            }
722
723            /* Apply the change if interesting */
724            if(best <= 0.)
725                continue;
726            if(opx || opy)
727            {
728                d2 = get(dest, x + opx, y + opy);
729                put(dest, x + opx, y + opy, d);
730            }
731            else
732                d2 = 1. - d;
733            put(dest, x, y, d2);
734            for(j = -N; j < N + 1; j++)
735            for(i = -N; i < N + 1; i++)
736            {
737                double m = mat[i + N][j + N];
738                if(y + j >= 0 && y + j < HEIGHT
739                    && x + i >= 0 && x + i < WIDTH)
740                {
741                    double t = get(tmp2, x + i, y + j);
742                    put(tmp2, x + i, y + j, t + m * (d2 - d));
743                }
744                if((opx || opy) && y + opy + j >= 0 && y + opy + j < HEIGHT
745                                && x + opx + i >= 0 && x + opx + i < WIDTH)
746                {
747                    double t = get(tmp2, x + opx + i, y + opy + j);
748                    put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2));
749                }
750            }
751
752            changes++;
753        }
754
755        fprintf(stderr, "did %i changes\n", changes);
756
757        if(changes == 0)
758            break;
759    }
760
761    del(tmp);
762    del(tmp2);
763
764    return dest;
765}
766
767static void study(double const *src, double const *dest,
768                  double sigma, double precision, double fdx, double fdy)
769{
770#   define Z 3
771    double mat[NN][NN], mat2[NN][NN];
772    double *tmp, *tmp2;
773    double e, e0, e1;
774    double best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
775    int dx, dy;
776
777    makegauss(mat, sigma, 0., 0.);
778    tmp = gauss(src, mat);
779
780    tmp2 = gauss(dest, mat);
781    e0 = dist(tmp, tmp2, 1.);
782    del(tmp2);
783
784    makegauss(mat2, sigma, fdx, fdy);
785    tmp2 = gauss(dest, mat2);
786    e1 = dist(tmp, tmp2, 1.);
787    del(tmp2);
788
789    while(step > precision)
790    {
791        for(dy = 0; dy <= Z; dy++)
792            for(dx = 0; dx <= Z; dx++)
793            {
794                makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z);
795                tmp2 = gauss(dest, mat);
796                e = dist(tmp, tmp2, best);
797                del(tmp2);
798                if(e < best)
799                {
800                    best = e;
801                    bfx = fx + step * dx / Z;
802                    bfy = fy + step * dy / Z;
803                }
804            }
805
806        fx = bfx - step / Z;
807        fy = bfy - step / Z;
808        step = step * 2 / Z;
809    }
810
811    del(tmp);
812
813    msg("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n",
814        1000. * e0, 1000. * e1, 1000. * best, fx, fy);
815}
816
817static double *merge(double const *im1, double const *im2, double t)
818{
819    double *dest = new();
820    int x, y;
821
822    for(y = 0; y < HEIGHT; y++)
823    for(x = 0; x < WIDTH; x++)
824        put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y));
825
826    return dest;
827}
828
829static void usage(char *argv[])
830{
831    fprintf(stderr, "Usage: %s <mode> [ARGS...]\n", argv[0]);
832    fprintf(stderr, "Allowed modes:\n");
833    fprintf(stderr, " -1 <src>           raster FS displacement study on src\n");
834    fprintf(stderr, " -2 <src1> <src2>   raster FS displacement study on blends of src1 and src2\n");
835    fprintf(stderr, " -3 <src>           quick (a,b,c,d) ED kernel analysis on src\n");
836    fprintf(stderr, " -4 <src>           exhaustive (a,b,c,d) ED kernel analysis on src\n");
837    fprintf(stderr, " -5 <src>           exhaustive displacement study on src\n");
838}
839
840int main(int argc, char *argv[])
841{
842    double *src;
843    int mode, i;
844
845    if(argc < 2)
846    {
847        fprintf(stderr, "%s: too few arguments\n", argv[0]);
848        usage(argv);
849        return EXIT_FAILURE;
850    }
851
852    if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1)))
853    {
854        fprintf(stderr, "%s: invalid mode `%s'\n", argv[0], argv[1]);
855        usage(argv);
856        return EXIT_FAILURE;
857    }
858
859    src = load(argv[2]);
860    if(!src)
861        return 2;
862
863    msg("### mode %i on `%s' ###\n", mode, argv[2]);
864
865    switch(mode)
866    {
867        case 1:
868        {
869            double *dest = ed(src, false, 7, 0,
870                                0, 3, 5, 1, 0,
871                                0, 0, 0, 0, 0);
872            study(src, dest, 1.2, 0.001, .16, .28);
873            del(dest);
874            del(src);
875        }
876        break;
877
878        case 2:
879        {
880            double *src2, *dest, *tmp;
881
882            if(argc < 3)
883                return 3;
884            src2 = load(argv[2]);
885            if(!src2)
886                return 4;
887
888            for(i = 0; i <= 100; i++)
889            {
890                tmp = merge(src, src2, (double)i / 100.);
891                dest = ed(tmp, false, 7, 0,
892                             0, 3, 5, 1, 0,
893                             0, 0, 0, 0, 0);
894                study(tmp, dest, 1.2, 0.001, .16, .28);
895                del(dest);
896                del(tmp);
897            }
898            del(src2);
899            del(src);
900        }
901        break;
902
903        case 3:
904        case 4:
905        {
906            double mat[NN][NN];
907            double *dest, *tmp, *tmp2;
908            int a, b, c, d, e;
909
910#define TOTAL 16
911            makegauss(mat, 1.2, 0, 0);
912            for(a = 0; a <= TOTAL; a++)
913            for(b = 0; b <= TOTAL; b++)
914            for(c = 0; c <= TOTAL; c++)
915            for(d = 0; d <= TOTAL; d++)
916            for(e = 0; e <= TOTAL; e++)
917            {
918                int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e, i;
919
920                if(a + b + c + d + e != TOTAL)
921                    continue;
922
923                /* Slightly shuffle our coefficients to avoid waiting until
924                 * 75% progress before having an idea of what's going on. */
925#define SHUFFLE(p,q,n) \
926    if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; }
927                SHUFFLE(c2, d2, 777); SHUFFLE(b2, d2, 555);
928                SHUFFLE(a2, d2, 333); SHUFFLE(b2, c2, 222);
929                SHUFFLE(a2, e2, 333); SHUFFLE(b2, e2, 222);
930                SHUFFLE(a2, c2, 444); SHUFFLE(a2, b2, 666);
931                SHUFFLE(c2, e2, 999); SHUFFLE(d2, e2, 777);
932                SHUFFLE(a2, d2, 999); SHUFFLE(a2, b2, 777);
933
934#if 0
935                if(b2 > c2) continue;
936                if(d2 > c2) continue;
937                if(d2 > a2) continue;
938#endif
939                /* We only want 4-cell kernels for now */
940                if(b2) continue;
941
942                for(i = 1; i <= 2; i++)
943                {
944                    //msg("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2);
945                    msg("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2);
946
947                    dest = ed(src, i == 2, a2, 0,
948                               b2, c2, d2, e2, 0,
949                                0,  0,  0,  0, 0);
950                    if(mode == 4)
951                    {
952                        /* XXX: E_fast is meaningless, ignore it */
953                        study(src, dest, 1.2, 0.001, 0., 0.);
954                    }
955                    else
956                    {
957                        tmp = gauss(src, mat);
958                        tmp2 = gauss(dest, mat);
959                        msg("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.));
960                        del(tmp);
961                        del(tmp2);
962                    }
963                    fflush(stdout);
964                    del(dest);
965                }
966            }
967
968            del(src);
969        }
970        break;
971
972        case 5:
973        {
974            double *dest;
975
976            dest = ed(src, false, 7, 0,
977                         0, 3, 5, 1, 0,
978                         0, 0, 0, 0, 0);
979            msg("[1] ");
980            study(src, dest, 1.2, 0.001, .16, .28);
981            del(dest);
982
983            dest = ed(src, false, 7, 5,
984                         3, 5, 7, 5, 3,
985                         1, 3, 5, 3, 1);
986            msg("[2] ");
987            study(src, dest, 1.2, 0.001, .26, .76);
988            del(dest);
989
990            dest = ostromoukhov(src);
991            msg("[3] ");
992            study(src, dest, 1.2, 0.001, .0, .19);
993            del(dest);
994
995            dest = ed(src, true, 2911, 0,
996                  0, 1373, 3457, 2258, 0,
997                  0,    0,    0,    0, 0);
998            msg("[4] ");
999            study(src, dest, 1.2, 0.001, .0, .34);
1000        }
1001        break;
1002
1003        case 6:
1004        case 7:
1005        {
1006            double mat[NN][NN];
1007            double *dest;
1008            int a, ax, b, bx, c, cx, d, dx;
1009
1010            if(mode == 6)
1011                a = 7, b = 3, c = 5, d = 1;
1012            else
1013                a = 7, b = 5, c = 4, d = 0;
1014
1015#undef TOTAL
1016#define TOTAL (a+b+c+d)
1017            makegauss(mat, 1.2, 0, 0);
1018            for(ax = -2; ax <= 2; ax++)
1019            for(bx = -2; bx <= 2; bx++)
1020            for(cx = -2; cx <= 2; cx++)
1021            for(dx = -2; dx <= 3; dx++)
1022            {
1023                int a2 = a + ax, b2 = b + bx, c2 = c + cx, d2 = d + dx;
1024
1025                if(a2 < 0 || a2 > TOTAL || b2 < 0 || b2 > TOTAL ||
1026                     c2 < 0 || c2 > TOTAL || d2 < 0 || d2 > TOTAL)
1027                    continue;
1028
1029                if(a2 + b2 + c2 + d2 != TOTAL)
1030                    continue;
1031
1032                msg("K: %d,%d,%d,%d ", a2, b2, c2, d2);
1033
1034                dest = ed(src, mode == 7, a2, 0,
1035                               0, b2, c2, d2, 0,
1036                               0,  0,  0,  0, 0);
1037                /* XXX: E_fast is meaningless, ignore it */
1038                study(src, dest, 1.2, 0.001, 0., 0.);
1039                fflush(stdout);
1040                del(dest);
1041            }
1042
1043            del(src);
1044        }
1045        break;
1046
1047#if 0
1048    tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1049    //dest = dbs(src, tmp, 0., 0.);
1050    dest = dbs(src, tmp, 0.20, 0.30);
1051    //dest = dbs(src, tmp, 0.158718, 0.283089);
1052    //dest = copy(tmp);
1053    del(tmp);
1054    study(src, dest, 1.2, 0.00001);
1055    save(dest, "output.bmp");
1056    del(dest);
1057#endif
1058
1059#if 0
1060    //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0);
1061    //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1062    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1063    //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
1064    //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
1065    //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
1066    //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0);
1067    //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0);
1068    //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0);
1069    dest = ostromoukhov(src);
1070    //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
1071    //msg("%s: ", argv[1]);
1072    //study(src, dest, 1.2, 0.001);
1073    save(dest, "output.bmp");
1074    del(dest);
1075#endif
1076
1077#if 0
1078#   define STEP 32
1079    dest = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1080    makegauss(mat, 1.2, 0., 0.);
1081    tmp = gauss(src, mat);
1082    for(dy = 0; dy < STEP; dy++)
1083    {
1084        for(dx = 0; dx < STEP; dx++)
1085        {
1086            double fy = 2. / STEP * (dy - STEP / 2.);
1087            double fx = 2. / STEP * (dx - STEP / 2.);
1088
1089            makegauss(mat, 1.2, fx, fy);
1090            tmp2 = gauss(dest, mat);
1091            msg("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
1092            fflush(stdout);
1093            del(tmp2);
1094        }
1095        msg("\n");
1096    }
1097
1098    save(dest, "output.bmp");
1099#endif
1100
1101#if 0
1102    tmp = gauss(src, 0., 0.);
1103    for(a = 0; a < 16; a++)
1104    for(b = 0; b < 16; b++)
1105    for(c = 0; c < 16; c++)
1106    for(d = 0; d < 16; d++)
1107    {
1108        if(a + b + c + d != 16)
1109            continue;
1110        dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
1111        tmp2 = gauss(dest, 0., 0.);
1112        msg("%.5g: (%02d %02d %02d %02d)\n",
1113            1000. * dist(tmp, tmp2, 1.), a, b, c, d);
1114        del(dest);
1115        del(tmp2);
1116    }
1117
1118    save(dest, "output.bmp");
1119#endif
1120
1121#if 0
1122    msg("%s\n", argv[1]);
1123
1124    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1125    study(src, dest, 1.2, 0.01);
1126    del(dest);
1127
1128    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
1129    study(src, dest, 1.2, 0.01);
1130    del(dest);
1131
1132    dest = ostromoukhov(src);
1133    study(src, dest, 1.2, 0.01);
1134    del(dest);
1135
1136    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
1137    study(src, dest, 1.2, 0.01);
1138    del(dest);
1139
1140    msg("\n");
1141#endif
1142
1143#if 0
1144    //dest = ostromoukhov(src);
1145    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1146    tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1147    dest = dbs(src, tmp, 0., 0.);
1148    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
1149    //for(sigma = 0.8; sigma < 2.; sigma *= 1.01)
1150    {
1151        msg("%g ", sigma);
1152        study(src, dest, sigma, 0.01);
1153    }
1154#endif
1155
1156    }
1157
1158#if 0
1159    tmp = new();
1160    a = 0;
1161    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
1162    {
1163        char buf[1024];
1164        msg("%i: %g\n", a, sigma);
1165        dest = dbs(src, tmp, sigma, 0., 0.);
1166        sprintf(buf, "output-dbs-%i.bmp", a++);
1167        save(dest, buf);
1168    }
1169#endif
1170
1171    return 0;
1172}
1173
Note: See TracBrowser for help on using the repository browser.