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

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