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

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