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 | |
---|
18 | int WIDTH, HEIGHT; |
---|
19 | |
---|
20 | static inline float get(float *src, int x, int y) |
---|
21 | { |
---|
22 | return src[y * WIDTH + x]; |
---|
23 | } |
---|
24 | |
---|
25 | static inline void put(float *src, int x, int y, float p) |
---|
26 | { |
---|
27 | src[y * WIDTH + x] = p; |
---|
28 | } |
---|
29 | |
---|
30 | static float *new(void) |
---|
31 | { |
---|
32 | return malloc(WIDTH * HEIGHT * sizeof(float)); |
---|
33 | } |
---|
34 | |
---|
35 | static 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 | |
---|
45 | static 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 | |
---|
66 | static 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 | |
---|
86 | static 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 | |
---|
108 | static 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 | |
---|
123 | static 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 | |
---|
138 | static 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 | |
---|
171 | static 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 | |
---|
192 | static 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 | */ |
---|
278 | static 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 */ |
---|
359 | static 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 | |
---|
500 | static 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 = 0., 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("E = %g E_min = %g dx = %g dy = %g\n", |
---|
544 | 1000 * e0, 1000 * best, fx, fy); |
---|
545 | fflush(stdout); |
---|
546 | } |
---|
547 | |
---|
548 | static float *merge(float *im1, float *im2, float t) |
---|
549 | { |
---|
550 | float *dest = new(); |
---|
551 | int x, y; |
---|
552 | |
---|
553 | for(y = 0; y < HEIGHT; y++) |
---|
554 | for(x = 0; x < WIDTH; x++) |
---|
555 | put(dest, x, y, t * get(im1, x, y) + (1. - t) * get(im2, x, y)); |
---|
556 | |
---|
557 | return dest; |
---|
558 | } |
---|
559 | |
---|
560 | static void usage(char *argv[]) |
---|
561 | { |
---|
562 | fprintf(stderr, "Usage: %s <mode> [ARGS...]\n", argv[0]); |
---|
563 | fprintf(stderr, "Allowed modes:\n"); |
---|
564 | fprintf(stderr, " -1 <src1> <src2> raster FS displacement study on blends of src1 and src2\n"); |
---|
565 | fprintf(stderr, " -2 <src> quick (a,b,c,d) ED kernel analysis on src\n"); |
---|
566 | fprintf(stderr, " -3 <src> exhaustive (a,b,c,d) ED kernel analysis on src\n"); |
---|
567 | } |
---|
568 | |
---|
569 | int main(int argc, char *argv[]) |
---|
570 | { |
---|
571 | float mat0[NN][NN]; |
---|
572 | float mat[NN][NN]; |
---|
573 | float *src; |
---|
574 | float sigma; |
---|
575 | int mode, dx, dy, i; |
---|
576 | |
---|
577 | if(argc < 2) |
---|
578 | { |
---|
579 | fprintf(stderr, "%s: too few arguments\n", argv[0]); |
---|
580 | usage(argv); |
---|
581 | return EXIT_FAILURE; |
---|
582 | } |
---|
583 | |
---|
584 | if(argv[1][0] != '-' || !(mode = atoi(argv[1] + 1))) |
---|
585 | { |
---|
586 | fprintf(stderr, "%s: invalid mode `%s'\n", argv[0], argv[1]); |
---|
587 | usage(argv); |
---|
588 | return EXIT_FAILURE; |
---|
589 | } |
---|
590 | |
---|
591 | src = load(argv[2]); |
---|
592 | if(!src) |
---|
593 | return 2; |
---|
594 | |
---|
595 | switch(mode) |
---|
596 | { |
---|
597 | case 1: |
---|
598 | { |
---|
599 | float *src2, *dest, *tmp; |
---|
600 | |
---|
601 | if(argc < 3) |
---|
602 | return 3; |
---|
603 | src2 = load(argv[2]); |
---|
604 | if(!src2) |
---|
605 | return 4; |
---|
606 | |
---|
607 | for(i = 0; i <= 100; i++) |
---|
608 | { |
---|
609 | tmp = merge(src, src2, (float)i / 100.); |
---|
610 | dest = ed(tmp, false, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0); |
---|
611 | study(tmp, dest, 1.2, 0.001); |
---|
612 | free(dest); |
---|
613 | free(tmp); |
---|
614 | } |
---|
615 | free(src2); |
---|
616 | free(src); |
---|
617 | } |
---|
618 | break; |
---|
619 | |
---|
620 | case 2: |
---|
621 | case 3: |
---|
622 | { |
---|
623 | float *dest, *tmp, *tmp2; |
---|
624 | int a, b, c, d, e; |
---|
625 | |
---|
626 | makegauss(mat0, 1.2, 0, 0); |
---|
627 | for(a = 0; a <= 16; a++) |
---|
628 | for(b = 0; b <= 16; b++) |
---|
629 | for(c = 0; c <= 16; c++) |
---|
630 | for(d = 0; d <= 16; d++) |
---|
631 | for(e = 0; e <= 16; e++) |
---|
632 | { |
---|
633 | if(a + b + c + d + e != 16) |
---|
634 | continue; |
---|
635 | |
---|
636 | /* Slightly shuffle our coefficients to avoid waiting until |
---|
637 | * 75% progress before having an idea of what's going on. */ |
---|
638 | int a2 = a, b2 = b, c2 = c, d2 = d, e2 = e; |
---|
639 | #define SHUFFLE(p,q,n) \ |
---|
640 | if(p+q) { int tmp = p+q; p = (p+n) % (tmp+1); q = tmp-p; } |
---|
641 | SHUFFLE(c2, d2, 7); SHUFFLE(b2, d2, 5); |
---|
642 | SHUFFLE(a2, d2, 3); SHUFFLE(b2, c2, 2); |
---|
643 | SHUFFLE(a2, e2, 3); SHUFFLE(b2, e2, 2); |
---|
644 | SHUFFLE(a2, c2, 4); SHUFFLE(a2, b2, 6); |
---|
645 | SHUFFLE(c2, e2, 9); SHUFFLE(d2, e2, 7); |
---|
646 | SHUFFLE(a2, d2, 9); SHUFFLE(a2, b2, 7); |
---|
647 | |
---|
648 | #if 0 |
---|
649 | if(b2 > c2) continue; |
---|
650 | if(d2 > c2) continue; |
---|
651 | if(d2 > a2) continue; |
---|
652 | #endif |
---|
653 | /* We only want 4-cell kernels for now */ |
---|
654 | if(b2) continue; |
---|
655 | //printf("K = %d,%d,%d,%d,%d ", a2, b2, c2, d2, e2); |
---|
656 | printf("K = %d,%d,%d,%d ", a2, c2, d2, e2); |
---|
657 | |
---|
658 | dest = ed(src, false, a2, 0, b2, c2, d2, e2, 0, 0, 0, 0, 0, 0); |
---|
659 | if(mode == 2) |
---|
660 | { |
---|
661 | tmp = gauss(src, mat0); |
---|
662 | tmp2 = gauss(dest, mat0); |
---|
663 | printf("E = %.5g\n", 1000. * dist(tmp, tmp2)); |
---|
664 | free(tmp); |
---|
665 | free(tmp2); |
---|
666 | } |
---|
667 | else |
---|
668 | study(src, dest, 1.2, 0.01); |
---|
669 | fflush(stdout); |
---|
670 | free(dest); |
---|
671 | } |
---|
672 | |
---|
673 | free(src); |
---|
674 | } |
---|
675 | break; |
---|
676 | |
---|
677 | #if 0 |
---|
678 | tmp = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
679 | //dest = dbs(src, tmp, 0., 0.); |
---|
680 | dest = dbs(src, tmp, 0.20, 0.30); |
---|
681 | //dest = dbs(src, tmp, 0.158718, 0.283089); |
---|
682 | //dest = copy(tmp); |
---|
683 | free(tmp); |
---|
684 | study(src, dest, 1.2, 0.00001); |
---|
685 | save(dest, "output.bmp"); |
---|
686 | free(dest); |
---|
687 | #endif |
---|
688 | |
---|
689 | #if 0 |
---|
690 | //dest = ed(src, 5, 0, 0, 3, 5, 3, 0, 0, 0, 0, 0, 0); |
---|
691 | //dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
692 | //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
693 | //dest = ed(src, true, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0); |
---|
694 | //dest = ed(src, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); |
---|
695 | //dest = ed(src, 7, 0, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0); |
---|
696 | //dest = ed(src, 7, 0, 1, 4, 4, 0, 0, 0, 0, 0, 0, 0); |
---|
697 | //dest = ed(src, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0); |
---|
698 | //dest = ed(src, 2911, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0, 0); |
---|
699 | dest = ostromoukhov(src); |
---|
700 | //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0); |
---|
701 | //printf("%s: ", argv[1]); |
---|
702 | //study(src, dest, 1.2, 0.001); |
---|
703 | save(dest, "output.bmp"); |
---|
704 | free(dest); |
---|
705 | #endif |
---|
706 | |
---|
707 | #if 0 |
---|
708 | # define STEP 32 |
---|
709 | dest = ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
710 | makegauss(mat, 1.2, 0., 0.); |
---|
711 | tmp = gauss(src, mat); |
---|
712 | for(dy = 0; dy < STEP; dy++) |
---|
713 | { |
---|
714 | for(dx = 0; dx < STEP; dx++) |
---|
715 | { |
---|
716 | float fy = 2. / STEP * (dy - STEP / 2.); |
---|
717 | float fx = 2. / STEP * (dx - STEP / 2.); |
---|
718 | |
---|
719 | makegauss(mat, 1.2, fx, fy); |
---|
720 | tmp2 = gauss(dest, mat); |
---|
721 | printf("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2)); |
---|
722 | fflush(stdout); |
---|
723 | free(tmp2); |
---|
724 | } |
---|
725 | printf("\n"); |
---|
726 | } |
---|
727 | |
---|
728 | save(dest, "output.bmp"); |
---|
729 | #endif |
---|
730 | |
---|
731 | #if 0 |
---|
732 | tmp = gauss(src, 0., 0.); |
---|
733 | for(a = 0; a < 16; a++) |
---|
734 | for(b = 0; b < 16; b++) |
---|
735 | for(c = 0; c < 16; c++) |
---|
736 | for(d = 0; d < 16; d++) |
---|
737 | { |
---|
738 | if(a + b + c + d != 16) |
---|
739 | continue; |
---|
740 | dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0); |
---|
741 | tmp2 = gauss(dest, 0., 0.); |
---|
742 | printf("%.5g: (%02d %02d %02d %02d)\n", |
---|
743 | 1000. * dist(tmp, tmp2), a, b, c, d); |
---|
744 | free(dest); |
---|
745 | free(tmp2); |
---|
746 | } |
---|
747 | |
---|
748 | save(dest, "output.bmp"); |
---|
749 | #endif |
---|
750 | |
---|
751 | #if 0 |
---|
752 | printf("%s: ", argv[1]); |
---|
753 | makegauss(mat0, 1.2, 0, 0); |
---|
754 | tmp = gauss(src, mat0); |
---|
755 | |
---|
756 | dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
757 | tmp2 = gauss(dest, mat0); |
---|
758 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
759 | free(tmp2); |
---|
760 | makegauss(mat, 1.2, 0.16, 0.28); |
---|
761 | tmp2 = gauss(dest, mat); |
---|
762 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
763 | free(tmp2); |
---|
764 | free(dest); |
---|
765 | |
---|
766 | dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); |
---|
767 | tmp2 = gauss(dest, mat0); |
---|
768 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
769 | free(tmp2); |
---|
770 | makegauss(mat, 1.2, 0.26, 0.76); |
---|
771 | tmp2 = gauss(dest, mat); |
---|
772 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
773 | free(tmp2); |
---|
774 | free(dest); |
---|
775 | |
---|
776 | dest = ostromoukhov(src); |
---|
777 | tmp2 = gauss(dest, mat0); |
---|
778 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
779 | free(tmp2); |
---|
780 | makegauss(mat, 1.2, 0.0, 0.19); |
---|
781 | tmp2 = gauss(dest, mat); |
---|
782 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
783 | free(tmp2); |
---|
784 | free(dest); |
---|
785 | |
---|
786 | dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0); |
---|
787 | tmp2 = gauss(dest, mat0); |
---|
788 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
789 | free(tmp2); |
---|
790 | makegauss(mat, 1.2, 0.0, 0.34); |
---|
791 | tmp2 = gauss(dest, mat); |
---|
792 | printf("%.5g ", 1000. * dist(tmp, tmp2)); |
---|
793 | free(tmp2); |
---|
794 | free(dest); |
---|
795 | |
---|
796 | printf("\n"); |
---|
797 | #endif |
---|
798 | |
---|
799 | #if 0 |
---|
800 | printf("%s\n", argv[1]); |
---|
801 | |
---|
802 | dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
803 | study(src, dest, 1.2, 0.01); |
---|
804 | free(dest); |
---|
805 | |
---|
806 | dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1); |
---|
807 | study(src, dest, 1.2, 0.01); |
---|
808 | free(dest); |
---|
809 | |
---|
810 | dest = ostromoukhov(src); |
---|
811 | study(src, dest, 1.2, 0.01); |
---|
812 | free(dest); |
---|
813 | |
---|
814 | dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0); |
---|
815 | study(src, dest, 1.2, 0.01); |
---|
816 | free(dest); |
---|
817 | |
---|
818 | printf("\n"); |
---|
819 | #endif |
---|
820 | |
---|
821 | #if 0 |
---|
822 | //dest = ostromoukhov(src); |
---|
823 | //dest = ed(src, true, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
824 | tmp = new();//ed(src, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0); |
---|
825 | dest = dbs(src, tmp, 0., 0.); |
---|
826 | for(sigma = 0.8; sigma < 2; sigma *= 1.03) |
---|
827 | //for(sigma = 0.8; sigma < 2.; sigma *= 1.01) |
---|
828 | { |
---|
829 | printf("%g ", sigma); |
---|
830 | study(src, dest, sigma, 0.01); |
---|
831 | } |
---|
832 | #endif |
---|
833 | |
---|
834 | } |
---|
835 | |
---|
836 | #if 0 |
---|
837 | tmp = new(); |
---|
838 | a = 0; |
---|
839 | for(sigma = 0.8; sigma < 2; sigma *= 1.03) |
---|
840 | { |
---|
841 | char buf[1024]; |
---|
842 | printf("%i: %g\n", a, sigma); |
---|
843 | dest = dbs(src, tmp, sigma, 0., 0.); |
---|
844 | sprintf(buf, "output-dbs-%i.bmp", a++); |
---|
845 | save(dest, buf); |
---|
846 | } |
---|
847 | #endif |
---|
848 | |
---|
849 | return 0; |
---|
850 | } |
---|
851 | |
---|