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

Last change on this file since 2277 was 2277, checked in by Sam Hocevar, 12 years ago
  • Put stuff relevant to the paper in a subdirectory.
File size: 24.3 KB
Line 
1/* test shit */
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
6#include <stdlib.h>
7
8#include <math.h>
9
10#include <SDL_image.h>
11
12#define MAXWIDTH 512
13#define MAXHEIGHT 512
14
15#define true 1
16#define false 0
17
18static int WIDTH, HEIGHT;
19
20static inline float get(float const *src, int x, int y)
21{
22    return src[y * WIDTH + x];
23}
24
25static inline void put(float *src, int x, int y, float p)
26{
27    src[y * WIDTH + x] = p;
28}
29
30static float *new(void)
31{
32    return malloc(WIDTH * HEIGHT * sizeof(float));
33}
34
35static float *copy(float const *src)
36{
37    float *dest = malloc(WIDTH * HEIGHT * sizeof(float));
38    memcpy(dest, src, WIDTH * HEIGHT * sizeof(float));
39    return dest;
40}
41
42#define N 5
43#define NN ((N * 2 + 1))
44
45static void makegauss(float mat[NN][NN], float sigma, float dx, float dy)
46{
47    float t = 0;
48    int i, j;
49
50    sigma = 2. * sigma * sigma;
51
52    for(j = 0; j < NN; j++)
53        for(i = 0; i < NN; i++)
54        {
55            float a = (float)(i - N) + dx;
56            float b = (float)(j - N) + dy;
57            mat[i][j] = pow(M_E, - (a * a + b * b) / sigma);
58            t += mat[i][j];
59        }
60
61    for(j = 0; j < NN; j++)
62        for(i = 0; i < NN; i++)
63            mat[i][j] /= t;
64}
65
66static float *gauss(float const *src, float mat[NN][NN])
67{
68    float *dest = new();
69    int x, y, i, j;
70
71    for(y = N; y < HEIGHT - N; y++)
72        for(x = N; x < WIDTH - N; x++)
73        {
74            float p = 0;
75
76            for(j = 0; j < NN; j++)
77                for(i = 0; i < NN; i++)
78                    p += get(src, x + i - N, y + j - N) * mat[i][j];
79
80            put(dest, x, y, p);
81        }
82
83    return dest;
84}
85
86static float *fullgauss(float const *src, float mat[NN][NN])
87{
88    float *dest = new();
89    int x, y, i, j;
90
91    for(y = 0; y < HEIGHT; y++)
92        for(x = 0; x < WIDTH; x++)
93        {
94            float p = 0;
95
96            for(j = 0; j < NN; j++)
97                for(i = 0; i < NN; i++)
98                    if(x + i >= N && x + i < WIDTH + N
99                       && y + j >= N && y + j < HEIGHT + N)
100                        p += get(src, x + i - N, y + j - N) * mat[i][j];
101
102            put(dest, x, y, p);
103        }
104
105    return dest;
106}
107
108#if 0
109static float fulldist(float const *p1, const float *p2)
110{
111    float error = 0.;
112    int x, y;
113
114    for(y = 0; y < HEIGHT; y++)
115        for(x = 0; x < WIDTH; x++)
116        {
117            float t = get(p1, x, y) - get(p2, x, y);
118            error += t * t;
119        }
120
121    return error / (WIDTH * HEIGHT);
122}
123#endif
124
125static float dist(float const *p1, float const *p2, float max)
126{
127    float error = 0.;
128    int x, y;
129
130    max *= (WIDTH - N) * (HEIGHT - N);
131
132    for(y = N; y < HEIGHT - N; y++)
133    {
134        for(x = N; x < WIDTH - N; x++)
135        {
136            float t = get(p1, x, y) - get(p2, x, y);
137            error += t * t;
138        }
139
140        /* Experience shows that this is almost useless, because the
141         * functions we manipulate are so small that the difference
142         * only happens in the last 1% of the image... */
143        if(error > max)
144            break;
145    }
146
147    return error / ((WIDTH - N) * (HEIGHT - N));
148}
149
150static float *load(char const *name)
151{
152    SDL_Surface *tmp, *surface;
153    uint32_t *pixels;
154    float *floats;
155    int x, y;
156
157    tmp = IMG_Load(name);
158    if(!tmp)
159        return NULL;
160
161    WIDTH = tmp->w > MAXWIDTH ? MAXWIDTH : tmp->w;
162    HEIGHT = tmp->h > MAXHEIGHT ? MAXHEIGHT : tmp->h;
163    floats = malloc(WIDTH * HEIGHT * sizeof(float));
164    if(!floats)
165        return NULL;
166
167    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
168                                   0xff0000, 0xff00, 0xff, 0x0);
169    pixels = (uint32_t *)surface->pixels;
170    SDL_BlitSurface(tmp, NULL, surface, NULL);
171    SDL_FreeSurface(tmp);
172
173    for(y = 0; y < HEIGHT; y++)
174        for(x = 0; x < WIDTH; x++)
175        {
176            int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
177            put(floats, x, y, (float)green / 0xff);
178        }
179
180    return floats;
181}
182
183static void save(float const *src, char const *name)
184{
185    SDL_Surface *surface;
186    uint32_t *pixels;
187    int x, y;
188
189    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
190                                   0xff0000, 0xff00, 0xff, 0x0);
191    pixels = (uint32_t *)surface->pixels;
192
193    for(y = 0; y < HEIGHT; y++)
194    for(x = 0; x < WIDTH; x++)
195    {
196        float p = 255 * get(src, x, y);
197        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
198        pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i;
199    }
200
201    SDL_SaveBMP(surface, name);
202}
203
204static float *ostromoukhov(float const *src)
205{
206    static int const table[][3] =
207    {
208         {13, 0, 5}, {13, 0, 5}, {21, 0, 10}, {7, 0, 4},
209         {8, 0, 5}, {47, 3, 28}, {23, 3, 13}, {15, 3, 8},
210         {22, 6, 11}, {43, 15, 20}, {7, 3, 3}, {501, 224, 211},
211         {249, 116, 103}, {165, 80, 67}, {123, 62, 49}, {489, 256, 191},
212         {81, 44, 31}, {483, 272, 181}, {60, 35, 22}, {53, 32, 19},
213         {237, 148, 83}, {471, 304, 161}, {3, 2, 1}, {481, 314, 185},
214         {354, 226, 155}, {1389, 866, 685}, {227, 138, 125}, {267, 158, 163},
215         {327, 188, 220}, {61, 34, 45}, {627, 338, 505}, {1227, 638, 1075},
216         {20, 10, 19}, {1937, 1000, 1767}, {977, 520, 855}, {657, 360, 551},
217         {71, 40, 57}, {2005, 1160, 1539}, {337, 200, 247}, {2039, 1240, 1425},
218         {257, 160, 171}, {691, 440, 437}, {1045, 680, 627}, {301, 200, 171},
219         {177, 120, 95}, {2141, 1480, 1083}, {1079, 760, 513}, {725, 520, 323},
220         {137, 100, 57}, {2209, 1640, 855}, {53, 40, 19}, {2243, 1720, 741},
221         {565, 440, 171}, {759, 600, 209}, {1147, 920, 285}, {2311, 1880, 513},
222         {97, 80, 19}, {335, 280, 57}, {1181, 1000, 171}, {793, 680, 95},
223         {599, 520, 57}, {2413, 2120, 171}, {405, 360, 19}, {2447, 2200, 57},
224         {11, 10, 0}, {158, 151, 3}, {178, 179, 7}, {1030, 1091, 63},
225         {248, 277, 21}, {318, 375, 35}, {458, 571, 63}, {878, 1159, 147},
226         {5, 7, 1}, {172, 181, 37}, {97, 76, 22}, {72, 41, 17},
227         {119, 47, 29}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
228         {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1},
229         {4, 1, 1}, {4, 1, 1}, {65, 18, 17}, {95, 29, 26},
230         {185, 62, 53}, {30, 11, 9}, {35, 14, 11}, {85, 37, 28},
231         {55, 26, 19}, {80, 41, 29}, {155, 86, 59}, {5, 3, 2},
232         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
233         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
234         {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2},
235         {305, 176, 119}, {155, 86, 59}, {105, 56, 39}, {80, 41, 29},
236         {65, 32, 23}, {55, 26, 19}, {335, 152, 113}, {85, 37, 28},
237         {115, 48, 37}, {35, 14, 11}, {355, 136, 109}, {30, 11, 9},
238         {365, 128, 107}, {185, 62, 53}, {25, 8, 7}, {95, 29, 26},
239         {385, 112, 103}, {65, 18, 17}, {395, 104, 101}, {4, 1, 1}
240    };
241
242    float *dest = new();
243    float *tmp = copy(src);
244    int x, y;
245
246    for(y = 0; y < HEIGHT; y++)
247    {
248        for(x = 0; x < WIDTH; x++)
249        {
250            int x1 = (y & 1) ? WIDTH - 1 - x + 1 : x - 1;
251            int x2 = (y & 1) ? WIDTH - 1 - x : x;
252            int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1;
253
254            float p = get(tmp, x2, y);
255            float q = p > 0.5 ? 1. : 0.;
256            float error = (p - q);
257            int i = p * 255.9999;
258
259            put(dest, x2, y, q);
260
261            if(i > 127)
262                i = 255 - i;
263            if(i < 0)
264                i = 0;
265
266            error /= table[i][0] + table[i][1] + table[i][2];
267
268            if(x < WIDTH - 1)
269                put(tmp, x3, y, get(tmp, x3, y) + error * table[i][0]);
270            if(y < HEIGHT - 1)
271            {
272                if(x > 0)
273                    put(tmp, x1, y + 1,
274                        get(tmp, x1, y + 1) + error * table[i][1]);
275                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * table[i][2]);
276            }
277        }
278    }
279
280    free(tmp);
281
282    return dest;
283}
284
285/* Dither using error diffusion and Floyd-Steinberg-like coefficients:
286         X a b
287     c d e f g
288     h i j k l
289*/
290static float *ed(float const *src, int serpentine,
291                 int a, int b, int c, int d, int e, int f,
292                 int g, int h, int i, int j, int k, int l)
293{
294    float *dest = new();
295    float *tmp = copy(src);
296    int x, y, n;
297
298    n = a + b + c + d + e + f + g + h + i + j + k + l;
299
300    for(y = 0; y < HEIGHT; y++)
301    {
302        int swap = serpentine && (y & 1);
303
304        for(x = 0; x < WIDTH; x++)
305        {
306            int x0 = swap ? WIDTH - 1 - x + 2 : x - 2;
307            int x1 = swap ? WIDTH - 1 - x + 1 : x - 1;
308            int x2 = swap ? WIDTH - 1 - x     : x;
309            int x3 = swap ? WIDTH - 1 - x - 1 : x + 1;
310            int x4 = swap ? WIDTH - 1 - x - 2 : x + 2;
311
312            float p = get(tmp, x2, y);
313            float q = p > 0.5 ? 1. : 0.;
314            float error = (p - q) / n;
315
316            put(dest, x2, y, q);
317
318            if(x < WIDTH - 1)
319                put(tmp, x3, y, get(tmp, x3, y) + error * a);
320            if(x < WIDTH - 2)
321                put(tmp, x4, y, get(tmp, x4, y) + error * b);
322            if(y < HEIGHT - 1)
323            {
324                if(x > 0)
325                {
326                    put(tmp, x1, y + 1,
327                        get(tmp, x1, y + 1) + error * d);
328                    if(x > 1)
329                        put(tmp, x0, y + 1,
330                            get(tmp, x0, y + 1) + error * c);
331                }
332                put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * e);
333                if(x < WIDTH - 1)
334                {
335                    put(tmp, x3, y + 1,
336                        get(tmp, x3, y + 1) + error * f);
337                    if(x < WIDTH - 2)
338                        put(tmp, x4, y + 1,
339                            get(tmp, x4, y + 1) + error * g);
340                }
341            }
342            if(y < HEIGHT - 2)
343            {
344                if(x > 0)
345                {
346                    put(tmp, x1, y + 2,
347                        get(tmp, x1, y + 2) + error * h);
348                    if(x > 1)
349                        put(tmp, x0, y + 2,
350                            get(tmp, x0, y + 2) + error * i);
351                }
352                put(tmp, x2, y + 2, get(tmp, x2, y + 2) + error * j);
353                if(x < WIDTH - 1)
354                {
355                    put(tmp, x3, y + 2,
356                        get(tmp, x3, y + 2) + error * k);
357                    if(x < WIDTH - 2)
358                        put(tmp, x4, y + 2,
359                            get(tmp, x4, y + 2) + error * l);
360                }
361            }
362        }
363    }
364
365    free(tmp);
366
367    return dest;
368}
369
370/* XXX */
371static float *dbs(float const *src, float const *orig,
372                  float sigma, float dx, float dy)
373{
374    float mat[NN][NN];
375    float *dest, *tmp, *tmp2;
376    float error;
377
378    makegauss(mat, sigma, 0., 0.);
379    tmp = fullgauss(src, mat);
380
381    makegauss(mat, sigma, dx, dy);
382    dest = copy(orig);
383    tmp2 = fullgauss(dest, mat);
384
385    error = dist(tmp, tmp2, 1.);
386
387    for(;;)
388    {
389        int changes = 0;
390        int x, y, i, j, n;
391
392        for(y = 0; y < HEIGHT; y++)
393        for(x = 0; x < WIDTH; x++)
394        {
395            float d, d2, e, best = 0.;
396            int opx = -1, opy = -1;
397
398            d = get(dest, x, y);
399
400            /* Compute the effect of a toggle */
401            e = 0.;
402            for(j = -N; j < N + 1; j++)
403            {
404                if(y + j < 0 || y + j >= HEIGHT)
405                    continue;
406
407                for(i = -N; i < N + 1; i++)
408                {
409                    float m, p, q1, q2;
410
411                    if(x + i < 0 || x + i >= WIDTH)
412                        continue;
413                   
414                    m = mat[i + N][j + N];
415                    p = get(tmp, x + i, y + j);
416                    q1 = get(tmp2, x + i, y + j);
417                    q2 = q1 - m * d + m * (1. - d);
418                    e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
419                }
420            }
421            if(e > best)
422            {
423                best = e;
424                opx = opy = 0;
425            }
426
427            /* Compute the effect of swaps */
428            for(n = 0; n < 8; n++)
429            {
430                static int const step[] =
431                    { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 };
432                int idx = step[n * 2], idy = step[n * 2 + 1];
433                if(y + idy < 0 || y + idy >= HEIGHT
434                     || x + idx < 0 || x + idx >= WIDTH)
435                    continue;
436                d2 = get(dest, x + idx, y + idy);
437                if(d2 == d)
438                    continue;
439                e = 0.;
440                for(j = -N; j < N + 1; j++)
441                {
442                    if(y + j < 0 || y + j >= HEIGHT)
443                        continue;
444                    if(j - idy + N < 0 || j - idy + N >= NN)
445                        continue;
446                    for(i = -N; i < N + 1; i++)
447                    {
448                        float ma, mb, p, q1, q2;
449                        if(x + i < 0 || x + i >= WIDTH)
450                            continue;
451                        if(i - idx + N < 0 || i - idx + N >= NN)
452                            continue;
453                        ma = mat[i + N][j + N];
454                        mb = mat[i - idx + N][j - idy + N];
455                        p = get(tmp, x + i, y + j);
456                        q1 = get(tmp2, x + i, y + j);
457                        q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d;
458                        e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p);
459                    }
460                }
461                if(e > best)
462                {
463                    best = e;
464                    opx = idx;
465                    opy = idy;
466                }
467            }
468
469            /* Apply the change if interesting */
470            if(best <= 0.)
471                continue;
472            if(opx || opy)
473            {
474                d2 = get(dest, x + opx, y + opy);
475                put(dest, x + opx, y + opy, d);
476            }
477            else
478                d2 = 1. - d;
479            put(dest, x, y, d2);
480            for(j = -N; j < N + 1; j++)
481            for(i = -N; i < N + 1; i++)
482            {
483                float m = mat[i + N][j + N];
484                if(y + j >= 0 && y + j < HEIGHT
485                    && x + i >= 0 && x + i < WIDTH)
486                {
487                    float t = get(tmp2, x + i, y + j);
488                    put(tmp2, x + i, y + j, t + m * (d2 - d));
489                }
490                if((opx || opy) && y + opy + j >= 0 && y + opy + j < HEIGHT
491                                && x + opx + i >= 0 && x + opx + i < WIDTH)
492                {
493                    float t = get(tmp2, x + opx + i, y + opy + j);
494                    put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2));
495                }
496            }
497
498            changes++;
499        }
500
501        fprintf(stderr, "did %i changes\n", changes);
502
503        if(changes == 0)
504            break;
505    }
506
507    free(tmp);
508    free(tmp2);
509
510    return dest;
511}
512
513static void study(float const *src, float const *dest,
514                  float sigma, float precision)
515{
516#   define Z 3
517    float mat[NN][NN];
518    float *tmp, *tmp2;
519    float e, e0, best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
520    int dx, dy;
521
522    makegauss(mat, sigma, 0., 0.);
523    tmp = gauss(src, mat);
524    tmp2 = gauss(dest, mat);
525
526    e0 = dist(tmp, tmp2, 1.);
527    free(tmp2);
528
529    while(step > precision)
530    {
531        for(dy = 0; dy <= Z; dy++)
532            for(dx = 0; dx <= Z; dx++)
533            {
534                makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z);
535                tmp2 = gauss(dest, mat);
536                e = dist(tmp, tmp2, best);
537                free(tmp2);
538                if(e < best)
539                {
540                    best = e;
541                    bfx = fx + step * dx / Z;
542                    bfy = fy + step * dy / Z;
543                }
544            }
545
546        fx = bfx - step / Z;
547        fy = bfy - step / Z;
548        step = step * 2 / Z;
549    }
550
551    free(tmp);
552
553    printf("E = %g E_min = %g dx = %g dy = %g\n",
554           1000 * e0, 1000 * best, fx, fy);
555    fflush(stdout);
556}
557
558static float *merge(float const *im1, float const *im2, float t)
559{
560    float *dest = new();
561    int x, y;
562
563    for(y = 0; y < HEIGHT; y++)
564    for(x = 0; x < WIDTH; x++)
565        put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y));
566
567    return dest;
568}
569
570static void usage(char *argv[])
571{
572    fprintf(stderr, "Usage: %s <mode> [ARGS...]\n", argv[0]);
573    fprintf(stderr, "Allowed modes:\n");
574    fprintf(stderr, " -1 <src>           raster FS displacement study on src\n");
575    fprintf(stderr, " -2 <src1> <src2>   raster FS displacement study on blends of src1 and src2\n");
576    fprintf(stderr, " -3 <src>           quick (a,b,c,d) ED kernel analysis on src\n");
577    fprintf(stderr, " -4 <src>           exhaustive (a,b,c,d) ED kernel analysis on src\n");
578}
579
580int main(int argc, char *argv[])
581{
582    float *src;
583    int mode, i;
584
585    if(argc < 2)
586    {
587        fprintf(stderr, "%s: too few arguments\n", argv[0]);
588        usage(argv);
589        return EXIT_FAILURE;
590    }
591
592    if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1)))
593    {
594        fprintf(stderr, "%s: invalid mode `%s'\n", argv[0], argv[1]);
595        usage(argv);
596        return EXIT_FAILURE;
597    }
598
599    src = load(argv[2]);
600    if(!src)
601        return 2;
602
603    switch(mode)
604    {
605        case 1:
606        {
607            float *dest = ed(src, false, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
608            study(src, dest, 1.2, 0.001);
609            free(dest);
610            free(src);
611        }
612        break;
613
614        case 2:
615        {
616            float *src2, *dest, *tmp;
617
618            if(argc < 3)
619                return 3;
620            src2 = load(argv[2]);
621            if(!src2)
622                return 4;
623
624            for(i = 0; i <= 100; i++)
625            {
626                tmp = merge(src, src2, (float)i / 100.);
627                dest = ed(tmp, false, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
628                study(tmp, dest, 1.2, 0.001);
629                free(dest);
630                free(tmp);
631            }
632            free(src2);
633            free(src);
634        }
635        break;
636
637        case 3:
638        case 4:
639        {
640            float mat0[NN][NN];
641            float *dest, *tmp, *tmp2;
642            int a, b, c, d, e;
643
644            makegauss(mat0, 1.2, 0, 0);
645            for(a = 0; a <= 16; a++)
646            for(b = 0; b <= 16; b++)
647            for(c = 0; c <= 16; c++)
648            for(d = 0; d <= 16; d++)
649            for(e = 0; e <= 16; e++)
650            {
651                if(a + b + c + d + e != 16)
652                    continue;
653
654                /* Slightly shuffle our coefficients to avoid waiting until
655                 * 75% progress before having an idea of what's going on. */
656                int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e;
657#define SHUFFLE(p,q,n) \
658    if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; }
659                SHUFFLE(c2, d2, 7); SHUFFLE(b2, d2, 5);
660                SHUFFLE(a2, d2, 3); SHUFFLE(b2, c2, 2);
661                SHUFFLE(a2, e2, 3); SHUFFLE(b2, e2, 2);
662                SHUFFLE(a2, c2, 4); SHUFFLE(a2, b2, 6);
663                SHUFFLE(c2, e2, 9); SHUFFLE(d2, e2, 7);
664                SHUFFLE(a2, d2, 9); SHUFFLE(a2, b2, 7);
665
666#if 0
667                if(b2 > c2) continue;
668                if(d2 > c2) continue;
669                if(d2 > a2) continue;
670#endif
671                /* We only want 4-cell kernels for now */
672                if(b2) continue;
673                //printf("K = %d,%d,%d,%d,%d ", a2, b2, c2, d2, e2);
674                printf("K = %d,%d,%d,%d ", a2, c2, d2, e2);
675
676                dest = ed(src, false, a2, 0, b2, c2, d2, e2, 0, 0, 0, 0, 0, 0);
677                if(mode == 3)
678                {
679                    tmp = gauss(src, mat0);
680                    tmp2 = gauss(dest, mat0);
681                    printf("E = %.5g\n", 1000. * dist(tmp, tmp2, 1.));
682                    free(tmp);
683                    free(tmp2);
684                }
685                else
686                    study(src, dest, 1.2, 0.001);
687                fflush(stdout);
688                free(dest);
689            }
690
691            free(src);
692        }
693        break;
694
695#if 0
696    tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
697    //dest = dbs(src, tmp, 0., 0.);
698    dest = dbs(src, tmp, 0.20, 0.30);
699    //dest = dbs(src, tmp, 0.158718, 0.283089);
700    //dest = copy(tmp);
701    free(tmp);
702    study(src, dest, 1.2, 0.00001);
703    save(dest, "output.bmp");
704    free(dest);
705#endif
706
707#if 0
708    //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0);
709    //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
710    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
711    //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
712    //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
713    //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0);
714    //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0);
715    //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0);
716    //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0);
717    dest = ostromoukhov(src);
718    //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
719    //printf("%s: ", argv[1]);
720    //study(src, dest, 1.2, 0.001);
721    save(dest, "output.bmp");
722    free(dest);
723#endif
724
725#if 0
726#   define STEP 32
727    dest = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
728    makegauss(mat, 1.2, 0., 0.);
729    tmp = gauss(src, mat);
730    for(dy = 0; dy < STEP; dy++)
731    {
732        for(dx = 0; dx < STEP; dx++)
733        {
734            float fy = 2. / STEP * (dy - STEP / 2.);
735            float fx = 2. / STEP * (dx - STEP / 2.);
736
737            makegauss(mat, 1.2, fx, fy);
738            tmp2 = gauss(dest, mat);
739            printf("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
740            fflush(stdout);
741            free(tmp2);
742        }
743        printf("\n");
744    }
745
746    save(dest, "output.bmp");
747#endif
748
749#if 0
750    tmp = gauss(src, 0., 0.);
751    for(a = 0; a < 16; a++)
752    for(b = 0; b < 16; b++)
753    for(c = 0; c < 16; c++)
754    for(d = 0; d < 16; d++)
755    {
756        if(a + b + c + d != 16)
757            continue;
758        dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
759        tmp2 = gauss(dest, 0., 0.);
760        printf("%.5g: (%02d %02d %02d %02d)\n",
761               1000. * dist(tmp, tmp2, 1.), a, b, c, d);
762        free(dest);
763        free(tmp2);
764    }
765
766    save(dest, "output.bmp");
767#endif
768
769#if 0
770    printf("%s: ", argv[1]);
771    makegauss(mat0, 1.2, 0, 0);
772    tmp = gauss(src, mat0);
773
774    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
775    tmp2 = gauss(dest, mat0);
776    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
777    free(tmp2);
778    makegauss(mat, 1.2, 0.16, 0.28);
779    tmp2 = gauss(dest, mat);
780    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
781    free(tmp2);
782    free(dest);
783
784    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
785    tmp2 = gauss(dest, mat0);
786    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
787    free(tmp2);
788    makegauss(mat, 1.2, 0.26, 0.76);
789    tmp2 = gauss(dest, mat);
790    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
791    free(tmp2);
792    free(dest);
793
794    dest = ostromoukhov(src);
795    tmp2 = gauss(dest, mat0);
796    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
797    free(tmp2);
798    makegauss(mat, 1.2, 0.0, 0.19);
799    tmp2 = gauss(dest, mat);
800    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
801    free(tmp2);
802    free(dest);
803
804    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
805    tmp2 = gauss(dest, mat0);
806    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
807    free(tmp2);
808    makegauss(mat, 1.2, 0.0, 0.34);
809    tmp2 = gauss(dest, mat);
810    printf("%.5g ", 1000. * dist(tmp, tmp2, 1.));
811    free(tmp2);
812    free(dest);
813
814    printf("\n");
815#endif
816
817#if 0
818    printf("%s\n", argv[1]);
819
820    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
821    study(src, dest, 1.2, 0.01);
822    free(dest);
823
824    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
825    study(src, dest, 1.2, 0.01);
826    free(dest);
827
828    dest = ostromoukhov(src);
829    study(src, dest, 1.2, 0.01);
830    free(dest);
831
832    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
833    study(src, dest, 1.2, 0.01);
834    free(dest);
835
836    printf("\n");
837#endif
838
839#if 0
840    //dest = ostromoukhov(src);
841    //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
842    tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
843    dest = dbs(src, tmp, 0., 0.);
844    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
845    //for(sigma = 0.8; sigma < 2.; sigma *= 1.01)
846    {
847        printf("%g ", sigma);
848        study(src, dest, sigma, 0.01);
849    }
850#endif
851
852    }
853
854#if 0
855    tmp = new();
856    a = 0;
857    for(sigma = 0.8; sigma < 2; sigma *= 1.03)
858    {
859        char buf[1024];
860        printf("%i: %g\n", a, sigma);
861        dest = dbs(src, tmp, sigma, 0., 0.);
862        sprintf(buf, "output-dbs-%i.bmp", a++);
863        save(dest, buf);
864    }
865#endif
866
867    return 0;
868}
869
Note: See TracBrowser for help on using the repository browser.