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