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 | |
---|
22 | static int WIDTH, HEIGHT; |
---|
23 | |
---|
24 | int main(int, char *[]); |
---|
25 | |
---|
26 | #ifdef BYTECODE |
---|
27 | # define MAXIMAGES 6 |
---|
28 | static int slots[MAXIMAGES]; |
---|
29 | static double slotbuf[MAXIMAGES * MAXWIDTH * MAXHEIGHT]; |
---|
30 | |
---|
31 | volatile int stop; |
---|
32 | |
---|
33 | void sighandler(int signal) |
---|
34 | { |
---|
35 | if(sys_write(1, "done", 4) != 4) |
---|
36 | sys_exit(3); |
---|
37 | stop = 1; |
---|
38 | } |
---|
39 | |
---|
40 | void 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 | |
---|
68 | static 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 | |
---|
95 | static 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 | |
---|
210 | static 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 | |
---|
218 | static 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 | |
---|
226 | static inline double get(double const *src, int x, int y) |
---|
227 | { |
---|
228 | return src[y * WIDTH + x]; |
---|
229 | } |
---|
230 | |
---|
231 | static inline void put(double *src, int x, int y, double p) |
---|
232 | { |
---|
233 | src[y * WIDTH + x] = p; |
---|
234 | } |
---|
235 | |
---|
236 | static 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 | |
---|
252 | static 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 | |
---|
267 | static 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 | |
---|
277 | static 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 | |
---|
298 | static 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 | |
---|
318 | static 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 |
---|
341 | static 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 | |
---|
357 | static 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 | |
---|
382 | static 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 | |
---|
450 | static 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 | |
---|
483 | static 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 | */ |
---|
569 | static 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 */ |
---|
650 | static 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 | |
---|
792 | static 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 | |
---|
842 | static 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 | |
---|
854 | static 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 | |
---|
868 | int 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 | |
---|