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

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