/* test shit */ #include #include #include #include #include #ifdef BYTECODE # include #else # include #endif #define MAXWIDTH 512 #define MAXHEIGHT 512 #define true 1 #define false 0 static int WIDTH, HEIGHT; int main(int, char *[]); #ifdef BYTECODE # define MAXIMAGES 6 static int slots[MAXIMAGES]; static double slotbuf[MAXIMAGES * MAXWIDTH * MAXHEIGHT]; volatile int stop; void sighandler(int signal) { if(sys_write(1, "done", 4) != 4) sys_exit(3); stop = 1; } void bytecode(unsigned char * mem, int heap_size, int stack_size) { char args[10][1024]; char *argv[] = { args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9] }; int i, argc; char c; if(sys_read(0, &c, 1) != 1) sys_exit(-5); argc = (unsigned char)c; for(i = 0; i < argc; i++) { char *p = argv[i]; do if(sys_read(0, p, 1) != 1) sys_exit(-5); while(*p++); } main(argc, argv); c = 0; if(sys_write(1, &c, 1) != 1) sys_exit(-6); } static int myatoi(const char *str) { int i = 0; while(*str) { i *= 10; i += (unsigned char)*str++ - (unsigned char)'0'; } return i; } #define atoi myatoi #endif #define WRITE_INT(i, base) \ do \ { \ char buf[128], *b = buf + 127; \ if(i <= 0) \ sys_write(1, (i = -i) ? "-" : "0", 1); /* XXX: hack here */ \ while(i) \ { \ *b-- = hex2char[i % base]; \ i /= base; \ } \ sys_write(1, b + 1, (int)(buf + 127 - b)); \ } while(0) static void out(FILE *stream, const char *f, va_list args) { #ifdef BYTECODE static char const *hex2char = "0123456789abcdef"; for( ; *f; f++) { if(*f != '%') { sys_write(1, f, 1); continue; } f++; if(!*f) break; if(*f == 'c') { char i = (char)(unsigned char)va_arg(args, int); if(i >= 0x20 && i < 0x7f) sys_write(1, &i, 1); else if(i == '\n') sys_write(1, "\\n", 2); else if(i == '\t') sys_write(1, "\\t", 2); else if(i == '\r') sys_write(1, "\\r", 2); else { sys_write(1, "\\x", 2); sys_write(1, hex2char + ((i & 0xf0) >> 4), 1); sys_write(1, hex2char + (i & 0x0f), 1); } } else if(*f == 'i' || *f == 'd') { int i = va_arg(args, int); WRITE_INT(i, 10); } else if(*f == 'x') { int i = va_arg(args, int); WRITE_INT(i, 16); } else if(f[0] == 'l' && (f[1] == 'i' || f[1] == 'd')) { long int i = va_arg(args, long int); WRITE_INT(i, 10); f++; } else if(f[0] == 'l' && f[1] == 'l' && (f[2] == 'i' || f[1] == 'd')) { long long int i = va_arg(args, long long int); WRITE_INT(i, 10); f += 2; } else if(f[0] == 'g') { double g = va_arg(args, double), h = 0.0000001; int i = g; WRITE_INT(i, 10); for(i = 0; i < 7; i++) { g = (g - (int)g) * 10; h *= 10; if(g < h) break; if(i == 0) sys_write(1, ".", 1); sys_write(1, hex2char + (int)g, 1); } } else if(f[0] == 'p') { uintptr_t i = va_arg(args, uintptr_t); if(!i) sys_write(1, "NULL", 5); else { sys_write(1, "0x", 2); WRITE_INT(i, 16); } } else if(f[0] == 's') { char *s = va_arg(args, char *); if(!s) sys_write(1, "(nil)", 5); else { int l = 0; while(s[l]) l++; sys_write(1, s, l); } } else if(f[0] == '0' && f[1] == '2' && f[2] == 'x') { int i = va_arg(args, int); sys_write(1, hex2char + ((i & 0xf0) >> 4), 1); sys_write(1, hex2char + (i & 0x0f), 1); f += 2; } else { sys_write(1, f - 1, 2); } } #else vfprintf(stream, f, args); fflush(stream); #endif } static void err(const char *f, ...) { va_list args; va_start(args, f); out(stderr, f, args); va_end(args); } static void msg(const char *f, ...) { va_list args; va_start(args, f); out(stdout, f, args); va_end(args); } static inline double get(double const *src, int x, int y) { return src[y * WIDTH + x]; } static inline void put(double *src, int x, int y, double p) { src[y * WIDTH + x] = p; } static double *new(void) { #ifdef BYTECODE int i; for(i = 0; i < MAXIMAGES; i++) if(slots[i] == 0) break; if(i == MAXIMAGES) return NULL; slots[i] = 1; return slotbuf + i * MAXWIDTH * MAXHEIGHT; #else return malloc(WIDTH * HEIGHT * sizeof(double)); #endif } static void del(double *img) { #ifdef BYTECODE int i; for(i = 0; i < MAXIMAGES; i++) if(slotbuf + i * MAXWIDTH * MAXHEIGHT == img) break; if(i == MAXIMAGES) return; slots[i] = 0; #else free(img); #endif } static double *copy(double const *src) { double *dest = new(); memcpy(dest, src, WIDTH * HEIGHT * sizeof(double)); return dest; } #define N 5 #define NN ((N * 2 + 1)) static void makegauss(double mat[NN][NN], double sigma, double dx, double dy) { double t = 0; int i, j; sigma = 2. * sigma * sigma; for(j = 0; j < NN; j++) for(i = 0; i < NN; i++) { double a = (double)(i - N) + dx; double b = (double)(j - N) + dy; mat[i][j] = pow(M_E, - (a * a + b * b) / sigma); t += mat[i][j]; } for(j = 0; j < NN; j++) for(i = 0; i < NN; i++) mat[i][j] /= t; } static double *gauss(double const *src, double mat[NN][NN]) { double *dest = new(); int x, y, i, j; for(y = N; y < HEIGHT - N; y++) for(x = N; x < WIDTH - N; x++) { double p = 0; for(j = 0; j < NN; j++) for(i = 0; i < NN; i++) p += get(src, x + i - N, y + j - N) * mat[i][j]; put(dest, x, y, p); } return dest; } static double *fullgauss(double const *src, double mat[NN][NN]) { double *dest = new(); int x, y, i, j; for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { double p = 0; for(j = 0; j < NN; j++) for(i = 0; i < NN; i++) if(x + i >= N && x + i < WIDTH + N && y + j >= N && y + j < HEIGHT + N) p += get(src, x + i - N, y + j - N) * mat[i][j]; put(dest, x, y, p); } return dest; } #if 0 static double fulldist(double const *p1, const double *p2) { double error = 0.; int x, y; for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { double t = get(p1, x, y) - get(p2, x, y); error += t * t; } return error / (WIDTH * HEIGHT); } #endif static double dist(double const *p1, double const *p2, double max) { double error = 0.; int x, y; max *= (WIDTH - N) * (HEIGHT - N); for(y = N; y < HEIGHT - N; y++) { for(x = N; x < WIDTH - N; x++) { double t = get(p1, x, y) - get(p2, x, y); error += t * t; } /* Experience shows that this is almost useless, because the * functions we manipulate are so small that the difference * only happens in the last 1% of the image... */ if(error > max) break; } return error / ((WIDTH - N) * (HEIGHT - N)); } static double *load(char const *name) { double *floats; int w, h, x, y; #ifdef BYTECODE char c; if(sys_read(0, &c, 1) != 1) sys_exit(-5); w = ((int)(unsigned char)c) << 8; if(sys_read(0, &c, 1) != 1) sys_exit(-5); w |= (int)(unsigned char)c; if(sys_read(0, &c, 1) != 1) sys_exit(-5); h = ((int)(unsigned char)c) << 8; if(sys_read(0, &c, 1) != 1) sys_exit(-5); h |= (int)(unsigned char)c; #else SDL_Surface *tmp, *surface; uint32_t *pixels; tmp = IMG_Load(name); if(!tmp) return NULL; w = tmp->w; h = tmp->h; #endif WIDTH = w > MAXWIDTH ? MAXWIDTH : w; HEIGHT = h > MAXHEIGHT ? MAXHEIGHT : h; floats = new(); if(!floats) return NULL; #ifdef BYTECODE for(y = 0; y < h; y++) for(x = 0; x < w; x++) { if(sys_read(0, &c, 1) != 1) sys_exit(-5); if(x >= WIDTH || y >= HEIGHT) continue; put(floats, x, y, (double)(unsigned char)c / 0xff); } #else surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32, 0xff0000, 0xff00, 0xff, 0x0); pixels = (uint32_t *)surface->pixels; SDL_BlitSurface(tmp, NULL, surface, NULL); SDL_FreeSurface(tmp); for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff; put(floats, x, y, (double)green / 0xff); } #endif return floats; } static void save(double const *src, char const *name) { int x, y; #ifdef BYTECODE for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { double p = 255 * get(src, x, y); uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p; char c = (char)(unsigned char)i; if(sys_write(1, &c, 1) != 1) sys_exit(-6); } #else SDL_Surface *surface; uint32_t *pixels; surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32, 0xff0000, 0xff00, 0xff, 0x0); pixels = (uint32_t *)surface->pixels; for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { double p = 255 * get(src, x, y); uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p; pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i; } SDL_SaveBMP(surface, name); #endif } static double *ostromoukhov(double const *src) { static int const table[][3] = { {13, 0, 5}, {13, 0, 5}, {21, 0, 10}, {7, 0, 4}, {8, 0, 5}, {47, 3, 28}, {23, 3, 13}, {15, 3, 8}, {22, 6, 11}, {43, 15, 20}, {7, 3, 3}, {501, 224, 211}, {249, 116, 103}, {165, 80, 67}, {123, 62, 49}, {489, 256, 191}, {81, 44, 31}, {483, 272, 181}, {60, 35, 22}, {53, 32, 19}, {237, 148, 83}, {471, 304, 161}, {3, 2, 1}, {481, 314, 185}, {354, 226, 155}, {1389, 866, 685}, {227, 138, 125}, {267, 158, 163}, {327, 188, 220}, {61, 34, 45}, {627, 338, 505}, {1227, 638, 1075}, {20, 10, 19}, {1937, 1000, 1767}, {977, 520, 855}, {657, 360, 551}, {71, 40, 57}, {2005, 1160, 1539}, {337, 200, 247}, {2039, 1240, 1425}, {257, 160, 171}, {691, 440, 437}, {1045, 680, 627}, {301, 200, 171}, {177, 120, 95}, {2141, 1480, 1083}, {1079, 760, 513}, {725, 520, 323}, {137, 100, 57}, {2209, 1640, 855}, {53, 40, 19}, {2243, 1720, 741}, {565, 440, 171}, {759, 600, 209}, {1147, 920, 285}, {2311, 1880, 513}, {97, 80, 19}, {335, 280, 57}, {1181, 1000, 171}, {793, 680, 95}, {599, 520, 57}, {2413, 2120, 171}, {405, 360, 19}, {2447, 2200, 57}, {11, 10, 0}, {158, 151, 3}, {178, 179, 7}, {1030, 1091, 63}, {248, 277, 21}, {318, 375, 35}, {458, 571, 63}, {878, 1159, 147}, {5, 7, 1}, {172, 181, 37}, {97, 76, 22}, {72, 41, 17}, {119, 47, 29}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {4, 1, 1}, {65, 18, 17}, {95, 29, 26}, {185, 62, 53}, {30, 11, 9}, {35, 14, 11}, {85, 37, 28}, {55, 26, 19}, {80, 41, 29}, {155, 86, 59}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {5, 3, 2}, {305, 176, 119}, {155, 86, 59}, {105, 56, 39}, {80, 41, 29}, {65, 32, 23}, {55, 26, 19}, {335, 152, 113}, {85, 37, 28}, {115, 48, 37}, {35, 14, 11}, {355, 136, 109}, {30, 11, 9}, {365, 128, 107}, {185, 62, 53}, {25, 8, 7}, {95, 29, 26}, {385, 112, 103}, {65, 18, 17}, {395, 104, 101}, {4, 1, 1} }; double *dest = new(); double *tmp = copy(src); int x, y; for(y = 0; y < HEIGHT; y++) { for(x = 0; x < WIDTH; x++) { int x1 = (y & 1) ? WIDTH - 1 - x + 1 : x - 1; int x2 = (y & 1) ? WIDTH - 1 - x : x; int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1; double p = get(tmp, x2, y); double q = p > 0.5 ? 1. : 0.; double error = (p - q); int i = p * 255.9999; put(dest, x2, y, q); if(i > 127) i = 255 - i; if(i < 0) i = 0; error /= table[i][0] + table[i][1] + table[i][2]; if(x < WIDTH - 1) put(tmp, x3, y, get(tmp, x3, y) + error * table[i][0]); if(y < HEIGHT - 1) { if(x > 0) put(tmp, x1, y + 1, get(tmp, x1, y + 1) + error * table[i][1]); put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * table[i][2]); } } } del(tmp); return dest; } /* Dither using error diffusion and Floyd-Steinberg-like coefficients: X a b c d e f g h i j k l */ static double *ed(double const *src, int serpentine, int a, int b, int c, int d, int e, int f, int g, int h, int i, int j, int k, int l) { double *dest = new(); double *tmp = copy(src); int x, y, n; n = a + b + c + d + e + f + g + h + i + j + k + l; for(y = 0; y < HEIGHT; y++) { int swap = serpentine && (y & 1); for(x = 0; x < WIDTH; x++) { int x0 = swap ? WIDTH - 1 - x + 2 : x - 2; int x1 = swap ? WIDTH - 1 - x + 1 : x - 1; int x2 = swap ? WIDTH - 1 - x : x; int x3 = swap ? WIDTH - 1 - x - 1 : x + 1; int x4 = swap ? WIDTH - 1 - x - 2 : x + 2; double p = get(tmp, x2, y); double q = p > 0.5 ? 1. : 0.; double error = (p - q) / n; put(dest, x2, y, q); if(x < WIDTH - 1) put(tmp, x3, y, get(tmp, x3, y) + error * a); if(x < WIDTH - 2) put(tmp, x4, y, get(tmp, x4, y) + error * b); if(y < HEIGHT - 1) { if(x > 0) { put(tmp, x1, y + 1, get(tmp, x1, y + 1) + error * d); if(x > 1) put(tmp, x0, y + 1, get(tmp, x0, y + 1) + error * c); } put(tmp, x2, y + 1, get(tmp, x2, y + 1) + error * e); if(x < WIDTH - 1) { put(tmp, x3, y + 1, get(tmp, x3, y + 1) + error * f); if(x < WIDTH - 2) put(tmp, x4, y + 1, get(tmp, x4, y + 1) + error * g); } } if(y < HEIGHT - 2) { if(x > 0) { put(tmp, x1, y + 2, get(tmp, x1, y + 2) + error * h); if(x > 1) put(tmp, x0, y + 2, get(tmp, x0, y + 2) + error * i); } put(tmp, x2, y + 2, get(tmp, x2, y + 2) + error * j); if(x < WIDTH - 1) { put(tmp, x3, y + 2, get(tmp, x3, y + 2) + error * k); if(x < WIDTH - 2) put(tmp, x4, y + 2, get(tmp, x4, y + 2) + error * l); } } } } del(tmp); return dest; } /* XXX */ static double *dbs(double const *src, double const *orig, double sigma, double dx, double dy) { double mat[NN][NN]; double *dest, *tmp, *tmp2; double error; makegauss(mat, sigma, 0., 0.); tmp = fullgauss(src, mat); makegauss(mat, sigma, dx, dy); dest = copy(orig); tmp2 = fullgauss(dest, mat); error = dist(tmp, tmp2, 1.); for(;;) { int changes = 0; int x, y, i, j, n; for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) { double d, d2, e, best = 0.; int opx = -1, opy = -1; d = get(dest, x, y); /* Compute the effect of a toggle */ e = 0.; for(j = -N; j < N + 1; j++) { if(y + j < 0 || y + j >= HEIGHT) continue; for(i = -N; i < N + 1; i++) { double m, p, q1, q2; if(x + i < 0 || x + i >= WIDTH) continue; m = mat[i + N][j + N]; p = get(tmp, x + i, y + j); q1 = get(tmp2, x + i, y + j); q2 = q1 - m * d + m * (1. - d); e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p); } } if(e > best) { best = e; opx = opy = 0; } /* Compute the effect of swaps */ for(n = 0; n < 8; n++) { static int const step[] = { 0, 1, 0, -1, -1, 0, 1, 0, -1, -1, -1, 1, 1, -1, 1, 1 }; int idx = step[n * 2], idy = step[n * 2 + 1]; if(y + idy < 0 || y + idy >= HEIGHT || x + idx < 0 || x + idx >= WIDTH) continue; d2 = get(dest, x + idx, y + idy); if(d2 == d) continue; e = 0.; for(j = -N; j < N + 1; j++) { if(y + j < 0 || y + j >= HEIGHT) continue; if(j - idy + N < 0 || j - idy + N >= NN) continue; for(i = -N; i < N + 1; i++) { double ma, mb, p, q1, q2; if(x + i < 0 || x + i >= WIDTH) continue; if(i - idx + N < 0 || i - idx + N >= NN) continue; ma = mat[i + N][j + N]; mb = mat[i - idx + N][j - idy + N]; p = get(tmp, x + i, y + j); q1 = get(tmp2, x + i, y + j); q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d; e += (q1 - p) * (q1 - p) - (q2 - p) * (q2 - p); } } if(e > best) { best = e; opx = idx; opy = idy; } } /* Apply the change if interesting */ if(best <= 0.) continue; if(opx || opy) { d2 = get(dest, x + opx, y + opy); put(dest, x + opx, y + opy, d); } else d2 = 1. - d; put(dest, x, y, d2); for(j = -N; j < N + 1; j++) for(i = -N; i < N + 1; i++) { double m = mat[i + N][j + N]; if(y + j >= 0 && y + j < HEIGHT && x + i >= 0 && x + i < WIDTH) { double t = get(tmp2, x + i, y + j); put(tmp2, x + i, y + j, t + m * (d2 - d)); } if((opx || opy) && y + opy + j >= 0 && y + opy + j < HEIGHT && x + opx + i >= 0 && x + opx + i < WIDTH) { double t = get(tmp2, x + opx + i, y + opy + j); put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2)); } } changes++; } err("did %i changes\n", changes); if(changes == 0) break; } del(tmp); del(tmp2); return dest; } static void study(double const *src, double const *dest, double sigma, double precision, double fdx, double fdy) { # define Z 3 double mat[NN][NN], mat2[NN][NN]; double *tmp, *tmp2; double e, e0, e1; double best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.; int dx, dy; makegauss(mat, sigma, 0., 0.); tmp = gauss(src, mat); tmp2 = gauss(dest, mat); e0 = dist(tmp, tmp2, 1.); del(tmp2); makegauss(mat2, sigma, fdx, fdy); tmp2 = gauss(dest, mat2); e1 = dist(tmp, tmp2, 1.); del(tmp2); while(step > precision) { for(dy = 0; dy <= Z; dy++) for(dx = 0; dx <= Z; dx++) { makegauss(mat, sigma, fx + step * dx / Z, fy + step * dy / Z); tmp2 = gauss(dest, mat); e = dist(tmp, tmp2, best); del(tmp2); if(e < best) { best = e; bfx = fx + step * dx / Z; bfy = fy + step * dy / Z; } } fx = bfx - step / Z; fy = bfy - step / Z; step = step * 2 / Z; } del(tmp); msg("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n", 1000. * e0, 1000. * e1, 1000. * best, fx, fy); } static double *merge(double const *im1, double const *im2, double t) { double *dest = new(); int x, y; for(y = 0; y < HEIGHT; y++) for(x = 0; x < WIDTH; x++) put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y)); return dest; } static void usage(char *argv[]) { msg("Usage: %s [ARGS...]\n", argv[0]); msg("Allowed modes:\n"); msg(" -1 raster FS displacement study on src\n"); msg(" -2 raster FS displacement study on blends of src1 and src2\n"); msg(" -3 quick (a,b,c,d) ED kernel analysis on src\n"); msg(" -4 exhaustive (a,b,c,d) ED kernel analysis on src\n"); msg(" -5 exhaustive displacement study on src\n"); msg(" -6 restrained (a,b,c,d) ED kernel analysis on src\n"); msg(" -7 restrained displacement study on src\n"); msg(" -8 displacement values on src\n"); } int main(int argc, char *argv[]) { double *src; int mode, i; if(argc < 2) { err("%s: too few arguments\n", argv[0]); usage(argv); return EXIT_FAILURE; } if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1))) { err("%s: invalid mode `%s'\n", argv[0], argv[1]); usage(argv); return EXIT_FAILURE; } src = load(argv[2]); if(!src) return 2; msg("### mode %i on `%s' ###\n", mode, argv[2]); switch(mode) { case 1: { double *dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); study(src, dest, 1.2, 0.001, .16, .28); del(dest); del(src); } break; case 2: { double *src2, *dest, *tmp; if(argc < 3) return 3; src2 = load(argv[2]); if(!src2) return 4; for(i = 0; i <= 100; i++) { tmp = merge(src, src2, (double)i / 100.); dest = ed(tmp, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); study(tmp, dest, 1.2, 0.001, .16, .28); del(dest); del(tmp); } del(src2); del(src); } break; case 3: case 4: { double mat[NN][NN]; double *dest, *tmp, *tmp2; int a, b, c, d, e; #define TOTAL 16 makegauss(mat, 1.2, 0, 0); for(a = 0; a <= TOTAL; a++) for(b = 0; b <= TOTAL; b++) for(c = 0; c <= TOTAL; c++) for(d = 0; d <= TOTAL; d++) for(e = 0; e <= TOTAL; e++) { int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e, i; if(a + b + c + d + e != TOTAL) continue; /* Slightly shuffle our coefficients to avoid waiting until * 75% progress before having an idea of what's going on. */ #define SHUFFLE(p,q,n) \ if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; } SHUFFLE(c2, d2, 777); SHUFFLE(b2, d2, 555); SHUFFLE(a2, d2, 333); SHUFFLE(b2, c2, 222); SHUFFLE(a2, e2, 333); SHUFFLE(b2, e2, 222); SHUFFLE(a2, c2, 444); SHUFFLE(a2, b2, 666); SHUFFLE(c2, e2, 999); SHUFFLE(d2, e2, 777); SHUFFLE(a2, d2, 999); SHUFFLE(a2, b2, 777); #if 0 if(b2 > c2) continue; if(d2 > c2) continue; if(d2 > a2) continue; #endif /* We only want 4-cell kernels for now */ if(b2) continue; for(i = 1; i <= 2; i++) { //msg("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2); msg("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2); dest = ed(src, i == 2, a2, 0, b2, c2, d2, e2, 0, 0, 0, 0, 0, 0); if(mode == 4) { /* XXX: E_fast is meaningless, ignore it */ study(src, dest, 1.2, 0.001, 0., 0.); } else { tmp = gauss(src, mat); tmp2 = gauss(dest, mat); msg("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.)); del(tmp); del(tmp2); } fflush(stdout); del(dest); } } del(src); } break; case 5: { double *dest; dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); msg("[1] "); study(src, dest, 1.2, 0.001, .16, .28); del(dest); dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); msg("[2] "); study(src, dest, 1.2, 0.001, .26, .76); del(dest); dest = ostromoukhov(src); msg("[3] "); study(src, dest, 1.2, 0.001, .0, .19); del(dest); dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0); msg("[4] "); study(src, dest, 1.2, 0.001, .0, .34); } break; case 6: case 7: { double mat[NN][NN]; double *dest; int a, ax, b, bx, c, cx, d, dx; if(mode == 6) a = 7, b = 3, c = 5, d = 1; else a = 7, b = 5, c = 4, d = 0; #undef TOTAL #define TOTAL (a+b+c+d) makegauss(mat, 1.2, 0, 0); for(ax = -2; ax <= 2; ax++) for(bx = -2; bx <= 2; bx++) for(cx = -2; cx <= 2; cx++) for(dx = -2; dx <= 3; dx++) { int a2 = a + ax, b2 = b + bx, c2 = c + cx, d2 = d + dx; if(a2 < 0 || a2 > TOTAL || b2 < 0 || b2 > TOTAL || c2 < 0 || c2 > TOTAL || d2 < 0 || d2 > TOTAL) continue; if(a2 + b2 + c2 + d2 != TOTAL) continue; msg("K: %d,%d,%d,%d ", a2, b2, c2, d2); dest = ed(src, mode == 7, a2, 0, 0, b2, c2, d2, 0, 0, 0, 0, 0, 0); /* XXX: E_fast is meaningless, ignore it */ study(src, dest, 1.2, 0.001, 0., 0.); fflush(stdout); del(dest); } del(src); } break; case 8: { const int STEP = 32; double mat[NN][NN]; double *dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); double *tmp, *tmp2; int dx, dy; makegauss(mat, 1.2, 0., 0.); tmp = gauss(src, mat); for(dy = 0; dy <= STEP; dy++) { for(dx = 0; dx <= STEP; dx++) { double fy = 2. / STEP * (dy - STEP / 2.); double fx = 2. / STEP * (dx - STEP / 2.); makegauss(mat, 1.2, fx, fy); tmp2 = gauss(dest, mat); msg("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.)); del(tmp2); } msg("\n"); } del(tmp); del(dest); del(src); } break; #if 0 tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); //dest = dbs(src, tmp, 0., 0.); dest = dbs(src, tmp, 0.20, 0.30); //dest = dbs(src, tmp, 0.158718, 0.283089); //dest = copy(tmp); del(tmp); study(src, dest, 1.2, 0.00001); save(dest, "output.bmp"); del(dest); #endif #if 0 //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0); //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0); //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0); //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0); //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0); //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0); dest = ostromoukhov(src); //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0); //msg("%s: ", argv[1]); //study(src, dest, 1.2, 0.001); save(dest, "output.bmp"); del(dest); #endif #if 0 save(dest, "output.bmp"); #endif #if 0 tmp = gauss(src, 0., 0.); for(a = 0; a < 16; a++) for(b = 0; b < 16; b++) for(c = 0; c < 16; c++) for(d = 0; d < 16; d++) { if(a + b + c + d != 16) continue; dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0); tmp2 = gauss(dest, 0., 0.); msg("%.5g: (%02d %02d %02d %02d)\n", 1000. * dist(tmp, tmp2, 1.), a, b, c, d); del(dest); del(tmp2); } save(dest, "output.bmp"); #endif #if 0 msg("%s\n", argv[1]); dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); study(src, dest, 1.2, 0.01); del(dest); dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); study(src, dest, 1.2, 0.01); del(dest); dest = ostromoukhov(src); study(src, dest, 1.2, 0.01); del(dest); dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0); study(src, dest, 1.2, 0.01); del(dest); msg("\n"); #endif #if 0 //dest = ostromoukhov(src); //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); dest = dbs(src, tmp, 0., 0.); for(sigma = 0.8; sigma < 2; sigma *= 1.03) //for(sigma = 0.8; sigma < 2.; sigma *= 1.01) { msg("%g ", sigma); study(src, dest, sigma, 0.01); } #endif } #if 0 tmp = new(); a = 0; for(sigma = 0.8; sigma < 2; sigma *= 1.03) { char buf[1024]; msg("%i: %g\n", a, sigma); dest = dbs(src, tmp, sigma, 0., 0.); sprintf(buf, "output-dbs-%i.bmp", a++); save(dest, buf); } #endif return 0; }