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

Revision 2288, 32.4 KB checked in by sam, 5 years ago (diff)
  • Add mode 8, for the very first graphic in the paper.
  • Add results for mode 8 in part0/lena-values.txt.
  • Make the bytecode version use less memory.
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 6
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 w, h, x, y;
386
387#ifdef BYTECODE
388    char c;
389
390    if(sys_read(0, &c, 1) != 1)
391        sys_exit(-5);
392    w = ((int)(unsigned char)c) << 8;
393    if(sys_read(0, &c, 1) != 1)
394        sys_exit(-5);
395    w |= (int)(unsigned char)c;
396
397    if(sys_read(0, &c, 1) != 1)
398        sys_exit(-5);
399    h = ((int)(unsigned char)c) << 8;
400    if(sys_read(0, &c, 1) != 1)
401        sys_exit(-5);
402    h |= (int)(unsigned char)c;
403#else
404    SDL_Surface *tmp, *surface;
405    uint32_t *pixels;
406
407    tmp = IMG_Load(name);
408    if(!tmp)
409        return NULL;
410
411    w = tmp->w;
412    h = tmp->h;
413#endif
414
415    WIDTH = w > MAXWIDTH ? MAXWIDTH : w;
416    HEIGHT = h > MAXHEIGHT ? MAXHEIGHT : h;
417
418    floats = new();
419    if(!floats)
420        return NULL;
421
422#ifdef BYTECODE
423    for(y = 0; y < h; y++)
424    for(x = 0; x < w; x++)
425    {
426        if(sys_read(0, &c, 1) != 1)
427            sys_exit(-5);
428        if(x >= WIDTH || y >= HEIGHT)
429            continue;
430        put(floats, x, y, (double)(unsigned char)c / 0xff);
431    }
432#else
433    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
434                                   0xff0000, 0xff00, 0xff, 0x0);
435    pixels = (uint32_t *)surface->pixels;
436    SDL_BlitSurface(tmp, NULL, surface, NULL);
437    SDL_FreeSurface(tmp);
438
439    for(y = 0; y < HEIGHT; y++)
440    for(x = 0; x < WIDTH; x++)
441    {
442        int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
443        put(floats, x, y, (double)green / 0xff);
444    }
445#endif
446
447    return floats;
448}
449
450static void save(double const *src, char const *name)
451{
452    int x, y;
453#ifdef BYTECODE
454    for(y = 0; y < HEIGHT; y++)
455    for(x = 0; x < WIDTH; x++)
456    {
457        double p = 255 * get(src, x, y);
458        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
459        char c = (char)(unsigned char)i;
460        if(sys_write(1, &c, 1) != 1)
461            sys_exit(-6);
462    }
463#else
464    SDL_Surface *surface;
465    uint32_t *pixels;
466
467    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
468                                   0xff0000, 0xff00, 0xff, 0x0);
469    pixels = (uint32_t *)surface->pixels;
470
471    for(y = 0; y < HEIGHT; y++)
472    for(x = 0; x < WIDTH; x++)
473    {
474        double p = 255 * get(src, x, y);
475        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
476        pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i;
477    }
478
479    SDL_SaveBMP(surface, name);
480#endif
481}
482
483static double *ostromoukhov(double const *src)
484{
485    static int const table[][3] =
486    {
487         {13, 0, 5}, {13, 0, 5}, {21, 0, 10}, {7, 0, 4},
488         {8, 0, 5}, {47, 3, 28}, {23, 3, 13}, {15, 3, 8},
489         {22, 6, 11}, {43, 15, 20}, {7, 3, 3}, {501, 224, 211},
490         {249, 116, 103}, {165, 80, 67}, {123, 62, 49}, {489, 256, 191},
491         {81, 44, 31}, {483, 272, 181}, {60, 35, 22}, {53, 32, 19},
492         {237, 148, 83}, {471, 304, 161}, {3, 2, 1}, {481, 314, 185},
493         {354, 226, 155}, {1389, 866, 685}, {227, 138, 125}, {267, 158, 163},
494         {327, 188, 220}, {61, 34, 45}, {627, 338, 505}, {1227, 638, 1075},
495         {20, 10, 19}, {1937, 1000, 1767}, {977, 520, 855}, {657, 360, 551},
496         {71, 40, 57}, {2005, 1160, 1539}, {337, 200, 247}, {2039, 1240, 1425},
497         {257, 160, 171}, {691, 440, 437}, {1045, 680, 627}, {301, 200, 171},
498         {177, 120, 95}, {2141, 1480, 1083}, {1079, 760, 513}, {725, 520, 323},
499         {137, 100, 57}, {2209, 1640, 855}, {53, 40, 19}, {2243, 1720, 741},
500         {565, 440, 171}, {759, 600, 209}, {1147, 920, 285}, {2311, 1880, 513},
501         {97, 80, 19}, {335, 280, 57}, {1181, 1000, 171}, {793, 680, 95},
502         {599, 520, 57}, {2413, 2120, 171}, {405, 360, 19}, {2447, 2200, 57},
503         {11, 10, 0}, {158, 151, 3}, {178, 179, 7}, {1030, 1091, 63},
504         {248, 277, 21}, {318, 375, 35}, {458, 571, 63}, {878, 1159, 147},
505         {5, 7, 1}, {172, 181, 37}, {97, 76, 22}, {72, 41, 17},
506         {119, 47, 29}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
507         {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
508         {4, 1, 1}, {4, 1, 1}, {65, 18, 17}, {95, 29, 26},
509         {185, 62, 53}, {30, 11, 9}, {35, 14, 11}, {85, 37, 28},
510         {55, 26, 19}, {80, 41, 29}, {155, 86, 59}, {5, 3, 2},
511         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
512         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
513         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
514         {305, 176, 119}, {155, 86, 59}, {105, 56, 39}, {80, 41, 29},
515         {65, 32, 23}, {55, 26, 19}, {335, 152, 113}, {85, 37, 28},
516         {115, 48, 37}, {35, 14, 11}, {355, 136, 109}, {30, 11, 9},
517         {365, 128, 107}, {185, 62, 53}, {25, 8, 7}, {95, 29, 26},
518         {385, 112, 103}, {65, 18, 17}, {395, 104, 101}, {4, 1, 1}
519    };
520
521    double *dest = new();
522    double *tmp = copy(src);
523    int x, y;
524
525    for(y = 0; y < HEIGHT; y++)
526    {
527        for(x = 0; x < WIDTH; x++)
528        {
529            int x1 = (y & 1) ? WIDTH - 1 - x + 1 : x - 1;
530            int x2 = (y & 1) ? WIDTH - 1 - x : x;
531            int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1;
532
533            double p = get(tmp, x2, y);
534            double q = p > 0.5 ? 1. : 0.;
535            double error = (p - q);
536            int i = p * 255.9999;
537
538            put(dest, x2, y, q);
539
540            if(i > 127)
541                i = 255 - i;
542            if(i < 0)
543                i = 0;
544
545            error /= table[i][0] + table[i][1] + table[i][2];
546
547            if(x < WIDTH - 1)
548                put(tmp, x3, y, get(tmp, x3, y) + error * table[i][0]);
549            if(y < HEIGHT - 1)
550            {
551                if(x > 0)
552                    put(tmp, x1, y + 1,
553                        get(tmp, x1, y + 1) + error * table[i][1]);
554                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * table[i][2]);
555            }
556        }
557    }
558
559    del(tmp);
560
561    return dest;
562}
563
564/* Dither using error diffusion and Floyd-Steinberg-like coefficients:
565         X a b
566     c d e f g
567     h i j k l
568*/
569static double *ed(double const *src, int serpentine,
570                 int a, int b, int c, int d, int e, int f,
571                 int g, int h, int i, int j, int k, int l)
572{
573    double *dest = new();
574    double *tmp = copy(src);
575    int x, y, n;
576
577    n = a + b + c + d + e + f + g + h + i + j + k + l;
578
579    for(y = 0; y < HEIGHT; y++)
580    {
581        int swap = serpentine && (y & 1);
582
583        for(x = 0; x < WIDTH; x++)
584        {
585            int x0 = swap ? WIDTH - 1 - x + 2 : x - 2;
586            int x1 = swap ? WIDTH - 1 - x + 1 : x - 1;
587            int x2 = swap ? WIDTH - 1 - x     : x;
588            int x3 = swap ? WIDTH - 1 - x - 1 : x + 1;
589            int x4 = swap ? WIDTH - 1 - x - 2 : x + 2;
590
591            double p = get(tmp, x2, y);
592            double q = p > 0.5 ? 1. : 0.;
593            double error = (p - q) / n;
594
595            put(dest, x2, y, q);
596
597            if(x < WIDTH - 1)
598                put(tmp, x3, y, get(tmp, x3, y) + error * a);
599            if(x < WIDTH - 2)
600                put(tmp, x4, y, get(tmp, x4, y) + error * b);
601            if(y < HEIGHT - 1)
602            {
603                if(x > 0)
604                {
605                    put(tmp, x1, y + 1,
606                        get(tmp, x1, y + 1) + error * d);
607                    if(x > 1)
608                        put(tmp, x0, y + 1,
609                            get(tmp, x0, y + 1) + error * c);
610                }
611                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * e);
612                if(x < WIDTH - 1)
613                {
614                    put(tmp, x3, y + 1,
615                        get(tmp, x3, y + 1) + error * f);
616                    if(x < WIDTH - 2)
617                        put(tmp, x4, y + 1,
618                            get(tmp, x4, y + 1) + error * g);
619                }
620            }
621            if(y < HEIGHT - 2)
622            {
623                if(x > 0)
624                {
625                    put(tmp, x1, y + 2,
626                        get(tmp, x1, y + 2) + error * h);
627                    if(x > 1)
628                        put(tmp, x0, y + 2,
629                            get(tmp, x0, y + 2) + error * i);
630                }
631                put(tmp, x2, y + 2, get(tmp, x2, y + 2) + error * j);
632                if(x < WIDTH - 1)
633                {
634                    put(tmp, x3, y + 2,
635                        get(tmp, x3, y + 2) + error * k);
636                    if(x < WIDTH - 2)
637                        put(tmp, x4, y + 2,
638                            get(tmp, x4, y + 2) + error * l);
639                }
640            }
641        }
642    }
643
644    del(tmp);
645
646    return dest;
647}
648
649/* XXX */
650static double *dbs(double const *src, double const *orig,
651                  double sigma, double dx, double dy)
652{
653    double mat[NN][NN];
654    double *dest, *tmp, *tmp2;
655    double error;
656
657    makegauss(mat, sigma, 0., 0.);
658    tmp = fullgauss(src, mat);
659
660    makegauss(mat, sigma, dx, dy);
661    dest = copy(orig);
662    tmp2 = fullgauss(dest, mat);
663
664    error = dist(tmp, tmp2, 1.);
665
666    for(;;)
667    {
668        int changes = 0;
669        int x, y, i, j, n;
670
671        for(y = 0; y < HEIGHT; y++)
672        for(x = 0; x < WIDTH; x++)
673        {
674            double d, d2, e, best = 0.;
675            int opx = -1, opy = -1;
676
677            d = get(dest, x, y);
678
679            /* Compute the effect of a toggle */
680            e = 0.;
681            for(j = -N; j < N + 1; j++)
682            {
683                if(y + j < 0 || y + j >= HEIGHT)
684                    continue;
685
686                for(i = -N; i < N + 1; i++)
687                {
688                    double m, p, q1, q2;
689
690                    if(x + i < 0 || x + i >= WIDTH)
691                        continue;
692                   
693                    m = mat[i + N][j + N];
694                    p = get(tmp, x + i, y + j);
695                    q1 = get(tmp2, x + i, y + j);
696                    q2 = q1 - m * d + m * (1. - d);
697                    e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
698                }
699            }
700            if(e > best)
701            {
702                best = e;
703                opx = opy = 0;
704            }
705
706            /* Compute the effect of swaps */
707            for(n = 0; n < 8; n++)
708            {
709                static int const step[] =
710                    { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
711                int idx = step[n * 2], idy = step[n * 2 + 1];
712                if(y + idy < 0 || y + idy >= HEIGHT
713                     || x + idx < 0 || x + idx >= WIDTH)
714                    continue;
715                d2 = get(dest, x + idx, y + idy);
716                if(d2 == d)
717                    continue;
718                e = 0.;
719                for(j = -N; j < N + 1; j++)
720                {
721                    if(y + j < 0 || y + j >= HEIGHT)
722                        continue;
723                    if(j - idy + N < 0 || j - idy + N >= NN)
724                        continue;
725                    for(i = -N; i < N + 1; i++)
726                    {
727                        double ma, mb, p, q1, q2;
728                        if(x + i < 0 || x + i >= WIDTH)
729                            continue;
730                        if(i - idx + N < 0 || i - idx + N >= NN)
731                            continue;
732                        ma = mat[i + N][j + N];
733                        mb = mat[i - idx + N][j - idy + N];
734                        p = get(tmp, x + i, y + j);
735                        q1 = get(tmp2, x + i, y + j);
736                        q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
737                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
738                    }
739                }
740                if(e > best)
741                {
742                    best = e;
743                    opx = idx;
744                    opy = idy;
745                }
746            }
747
748            /* Apply the change if interesting */
749            if(best <= 0.)
750                continue;
751            if(opx || opy)
752            {
753                d2 = get(dest, x + opx, y + opy);
754                put(dest, x + opx, y + opy, d);
755            }
756            else
757                d2 = 1. - d;
758            put(dest, x, y, d2);
759            for(j = -N; j < N + 1; j++)
760            for(i = -N; i < N + 1; i++)
761            {
762                double m = mat[i + N][j + N];
763                if(y + j >= 0 && y + j < HEIGHT
764                    && x + i >= 0 && x + i < WIDTH)
765                {
766                    double t = get(tmp2, x + i, y + j);
767                    put(tmp2, x + i, y + j, t + m * (d2 - d));
768                }
769                if((opx || opy) && y + opy + j >= 0 && y + opy + j < HEIGHT
770                                && x + opx + i >= 0 && x + opx + i < WIDTH)
771                {
772                    double t = get(tmp2, x + opx + i, y + opy + j);
773                    put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2));
774                }
775            }
776
777            changes++;
778        }
779
780        err("did %i changes\n", changes);
781
782        if(changes == 0)
783            break;
784    }
785
786    del(tmp);
787    del(tmp2);
788
789    return dest;
790}
791
792static void study(double const *src, double const *dest,
793                  double sigma, double precision, double fdx, double fdy)
794{
795#   define Z 3
796    double mat[NN][NN], mat2[NN][NN];
797    double *tmp, *tmp2;
798    double e, e0, e1;
799    double best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
800    int dx, dy;
801
802    makegauss(mat, sigma, 0., 0.);
803    tmp = gauss(src, mat);
804
805    tmp2 = gauss(dest, mat);
806    e0 = dist(tmp, tmp2, 1.);
807    del(tmp2);
808
809    makegauss(mat2, sigma, fdx, fdy);
810    tmp2 = gauss(dest, mat2);
811    e1 = dist(tmp, tmp2, 1.);
812    del(tmp2);
813
814    while(step > precision)
815    {
816        for(dy = 0; dy <= Z; dy++)
817            for(dx = 0; dx <= Z; dx++)
818            {
819                makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z);
820                tmp2 = gauss(dest, mat);
821                e = dist(tmp, tmp2, best);
822                del(tmp2);
823                if(e < best)
824                {
825                    best = e;
826                    bfx = fx + step * dx / Z;
827                    bfy = fy + step * dy / Z;
828                }
829            }
830
831        fx = bfx - step / Z;
832        fy = bfy - step / Z;
833        step = step * 2 / Z;
834    }
835
836    del(tmp);
837
838    msg("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n",
839        1000. * e0, 1000. * e1, 1000. * best, fx, fy);
840}
841
842static double *merge(double const *im1, double const *im2, double t)
843{
844    double *dest = new();
845    int x, y;
846
847    for(y = 0; y < HEIGHT; y++)
848    for(x = 0; x < WIDTH; x++)
849        put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y));
850
851    return dest;
852}
853
854static void usage(char *argv[])
855{
856    msg("Usage: %s <mode> [ARGS...]\n", argv[0]);
857    msg("Allowed modes:\n");
858    msg(" -1 <src>           raster FS displacement study on src\n");
859    msg(" -2 <src1> <src2>   raster FS displacement study on blends of src1 and src2\n");
860    msg(" -3 <src>           quick (a,b,c,d) ED kernel analysis on src\n");
861    msg(" -4 <src>           exhaustive (a,b,c,d) ED kernel analysis on src\n");
862    msg(" -5 <src>           exhaustive displacement study on src\n");
863    msg(" -6 <src>           restrained (a,b,c,d) ED kernel analysis on src\n");
864    msg(" -7 <src>           restrained displacement study on src\n");
865    msg(" -8 <src>           displacement values on src\n");
866}
867
868int main(int argc, char *argv[])
869{
870    double *src;
871    int mode, i;
872
873    if(argc < 2)
874    {
875        err("%s: too few arguments\n", argv[0]);
876        usage(argv);
877        return EXIT_FAILURE;
878    }
879
880    if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1)))
881    {
882        err("%s: invalid mode `%s'\n", argv[0], argv[1]);
883        usage(argv);
884        return EXIT_FAILURE;
885    }
886
887    src = load(argv[2]);
888    if(!src)
889        return 2;
890
891    msg("### mode %i on `%s' ###\n", mode, argv[2]);
892
893    switch(mode)
894    {
895        case 1:
896        {
897            double *dest = ed(src, false, 7, 0,
898                                 0, 3, 5, 1, 0,
899                                 0, 0, 0, 0, 0);
900            study(src, dest, 1.2, 0.001, .16, .28);
901            del(dest);
902            del(src);
903        }
904        break;
905
906        case 2:
907        {
908            double *src2, *dest, *tmp;
909
910            if(argc < 3)
911                return 3;
912            src2 = load(argv[2]);
913            if(!src2)
914                return 4;
915
916            for(i = 0; i <= 100; i++)
917            {
918                tmp = merge(src, src2, (double)i / 100.);
919                dest = ed(tmp, false, 7, 0,
920                             0, 3, 5, 1, 0,
921                             0, 0, 0, 0, 0);
922                study(tmp, dest, 1.2, 0.001, .16, .28);
923                del(dest);
924                del(tmp);
925            }
926            del(src2);
927            del(src);
928        }
929        break;
930
931        case 3:
932        case 4:
933        {
934            double mat[NN][NN];
935            double *dest, *tmp, *tmp2;
936            int a, b, c, d, e;
937
938#define TOTAL 16
939            makegauss(mat, 1.2, 0, 0);
940            for(a = 0; a <= TOTAL; a++)
941            for(b = 0; b <= TOTAL; b++)
942            for(c = 0; c <= TOTAL; c++)
943            for(d = 0; d <= TOTAL; d++)
944            for(e = 0; e <= TOTAL; e++)
945            {
946                int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e, i;
947
948                if(a + b + c + d + e != TOTAL)
949                    continue;
950
951                /* Slightly shuffle our coefficients to avoid waiting until
952                 * 75% progress before having an idea of what's going on. */
953#define SHUFFLE(p,q,n) \
954    if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; }
955                SHUFFLE(c2, d2, 777); SHUFFLE(b2, d2, 555);
956                SHUFFLE(a2, d2, 333); SHUFFLE(b2, c2, 222);
957                SHUFFLE(a2, e2, 333); SHUFFLE(b2, e2, 222);
958                SHUFFLE(a2, c2, 444); SHUFFLE(a2, b2, 666);
959                SHUFFLE(c2, e2, 999); SHUFFLE(d2, e2, 777);
960                SHUFFLE(a2, d2, 999); SHUFFLE(a2, b2, 777);
961
962#if 0
963                if(b2 > c2) continue;
964                if(d2 > c2) continue;
965                if(d2 > a2) continue;
966#endif
967                /* We only want 4-cell kernels for now */
968                if(b2) continue;
969
970                for(i = 1; i <= 2; i++)
971                {
972                    //msg("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2);
973                    msg("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2);
974
975                    dest = ed(src, i == 2, a2, 0,
976                               b2, c2, d2, e2, 0,
977                                0,  0,  0,  0, 0);
978                    if(mode == 4)
979                    {
980                        /* XXX: E_fast is meaningless, ignore it */
981                        study(src, dest, 1.2, 0.001, 0., 0.);
982                    }
983                    else
984                    {
985                        tmp = gauss(src, mat);
986                        tmp2 = gauss(dest, mat);
987                        msg("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.));
988                        del(tmp);
989                        del(tmp2);
990                    }
991                    fflush(stdout);
992                    del(dest);
993                }
994            }
995
996            del(src);
997        }
998        break;
999
1000        case 5:
1001        {
1002            double *dest;
1003
1004            dest = ed(src, false, 7, 0,
1005                         0, 3, 5, 1, 0,
1006                         0, 0, 0, 0, 0);
1007            msg("[1] ");
1008            study(src, dest, 1.2, 0.001, .16, .28);
1009            del(dest);
1010
1011            dest = ed(src, false, 7, 5,
1012                         3, 5, 7, 5, 3,
1013                         1, 3, 5, 3, 1);
1014            msg("[2] ");
1015            study(src, dest, 1.2, 0.001, .26, .76);
1016            del(dest);
1017
1018            dest = ostromoukhov(src);
1019            msg("[3] ");
1020            study(src, dest, 1.2, 0.001, .0, .19);
1021            del(dest);
1022
1023            dest = ed(src, true, 2911, 0,
1024                  0, 1373, 3457, 2258, 0,
1025                  0,    0,    0,    0, 0);
1026            msg("[4] ");
1027            study(src, dest, 1.2, 0.001, .0, .34);
1028        }
1029        break;
1030
1031        case 6:
1032        case 7:
1033        {
1034            double mat[NN][NN];
1035            double *dest;
1036            int a, ax, b, bx, c, cx, d, dx;
1037
1038            if(mode == 6)
1039                a = 7, b = 3, c = 5, d = 1;
1040            else
1041                a = 7, b = 5, c = 4, d = 0;
1042
1043#undef TOTAL
1044#define TOTAL (a+b+c+d)
1045            makegauss(mat, 1.2, 0, 0);
1046            for(ax = -2; ax <= 2; ax++)
1047            for(bx = -2; bx <= 2; bx++)
1048            for(cx = -2; cx <= 2; cx++)
1049            for(dx = -2; dx <= 3; dx++)
1050            {
1051                int a2 = a + ax, b2 = b + bx, c2 = c + cx, d2 = d + dx;
1052
1053                if(a2 < 0 || a2 > TOTAL || b2 < 0 || b2 > TOTAL ||
1054                     c2 < 0 || c2 > TOTAL || d2 < 0 || d2 > TOTAL)
1055                    continue;
1056
1057                if(a2 + b2 + c2 + d2 != TOTAL)
1058                    continue;
1059
1060                msg("K: %d,%d,%d,%d ", a2, b2, c2, d2);
1061
1062                dest = ed(src, mode == 7, a2, 0,
1063                               0, b2, c2, d2, 0,
1064                               0,  0,  0,  0, 0);
1065                /* XXX: E_fast is meaningless, ignore it */
1066                study(src, dest, 1.2, 0.001, 0., 0.);
1067                fflush(stdout);
1068                del(dest);
1069            }
1070
1071            del(src);
1072        }
1073        break;
1074
1075        case 8:
1076        {
1077            const int STEP = 32;
1078            double mat[NN][NN];
1079            double *dest = ed(src, false, 7, 0,
1080                                 0, 3, 5, 1, 0,
1081                                 0, 0, 0, 0, 0);
1082            double *tmp, *tmp2;
1083            int dx, dy;
1084
1085            makegauss(mat, 1.2, 0., 0.);
1086            tmp = gauss(src, mat);
1087            for(dy = 0; dy <= STEP; dy++)
1088            {
1089                for(dx = 0; dx <= STEP; dx++)
1090                {
1091                    double fy = 2. / STEP * (dy - STEP / 2.);
1092                    double fx = 2. / STEP * (dx - STEP / 2.);
1093
1094                    makegauss(mat, 1.2, fx, fy);
1095                    tmp2 = gauss(dest, mat);
1096                    msg("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
1097                    del(tmp2);
1098                }
1099                msg("\n");
1100            }
1101            del(tmp);
1102            del(dest);
1103            del(src);
1104        }
1105        break;
1106
1107#if 0
1108    tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1109    //dest = dbs(src, tmp, 0., 0.);
1110    dest = dbs(src, tmp, 0.20, 0.30);
1111    //dest = dbs(src, tmp, 0.158718, 0.283089);
1112    //dest = copy(tmp);
1113    del(tmp);
1114    study(src, dest, 1.2, 0.00001);
1115    save(dest, "output.bmp");
1116    del(dest);
1117#endif
1118
1119#if 0
1120    //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0);
1121    //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1122    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1123    //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
1124    //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
1125    //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
1126    //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0);
1127    //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0);
1128    //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0);
1129    dest = ostromoukhov(src);
1130    //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
1131    //msg("%s: ", argv[1]);
1132    //study(src, dest, 1.2, 0.001);
1133    save(dest, "output.bmp");
1134    del(dest);
1135#endif
1136
1137#if 0
1138
1139    save(dest, "output.bmp");
1140#endif
1141
1142#if 0
1143    tmp = gauss(src, 0., 0.);
1144    for(a = 0; a < 16; a++)
1145    for(b = 0; b < 16; b++)
1146    for(c = 0; c < 16; c++)
1147    for(d = 0; d < 16; d++)
1148    {
1149        if(a + b + c + d != 16)
1150            continue;
1151        dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
1152        tmp2 = gauss(dest, 0., 0.);
1153        msg("%.5g: (%02d %02d %02d %02d)\n",
1154            1000. * dist(tmp, tmp2, 1.), a, b, c, d);
1155        del(dest);
1156        del(tmp2);
1157    }
1158
1159    save(dest, "output.bmp");
1160#endif
1161
1162#if 0
1163    msg("%s\n", argv[1]);
1164
1165    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1166    study(src, dest, 1.2, 0.01);
1167    del(dest);
1168
1169    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
1170    study(src, dest, 1.2, 0.01);
1171    del(dest);
1172
1173    dest = ostromoukhov(src);
1174    study(src, dest, 1.2, 0.01);
1175    del(dest);
1176
1177    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
1178    study(src, dest, 1.2, 0.01);
1179    del(dest);
1180
1181    msg("\n");
1182#endif
1183
1184#if 0
1185    //dest = ostromoukhov(src);
1186    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1187    tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
1188    dest = dbs(src, tmp, 0., 0.);
1189    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
1190    //for(sigma = 0.8; sigma < 2.; sigma *= 1.01)
1191    {
1192        msg("%g ", sigma);
1193        study(src, dest, sigma, 0.01);
1194    }
1195#endif
1196
1197    }
1198
1199#if 0
1200    tmp = new();
1201    a = 0;
1202    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
1203    {
1204        char buf[1024];
1205        msg("%i: %g\n", a, sigma);
1206        dest = dbs(src, tmp, sigma, 0., 0.);
1207        sprintf(buf, "output-dbs-%i.bmp", a++);
1208        save(dest, buf);
1209    }
1210#endif
1211
1212    return 0;
1213}
1214
Note: See TracBrowser for help on using the repository browser.