1 | #include "config.h" |
---|
2 | |
---|
3 | #include <stdio.h> |
---|
4 | #include <stdlib.h> |
---|
5 | #include <string.h> |
---|
6 | #include <math.h> |
---|
7 | |
---|
8 | #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> |
---|
9 | #include <CGAL/Delaunay_triangulation_2.h> |
---|
10 | #include <CGAL/natural_neighbor_coordinates_2.h> |
---|
11 | |
---|
12 | #include <pipi.h> |
---|
13 | |
---|
14 | /* |
---|
15 | * User-definable settings. |
---|
16 | */ |
---|
17 | |
---|
18 | /* The maximum message length */ |
---|
19 | #define MAX_MSG_LEN 140 |
---|
20 | |
---|
21 | /* The number of characters at disposal */ |
---|
22 | //#define NUM_CHARACTERS 0x7fffffff // The sky's the limit |
---|
23 | //#define NUM_CHARACTERS 1111998 // Full valid Unicode set |
---|
24 | //#define NUM_CHARACTERS 100507 // Full graphic Unicode |
---|
25 | #define NUM_CHARACTERS 32768 // Chinese characters |
---|
26 | //#define NUM_CHARACTERS 127 // ASCII |
---|
27 | |
---|
28 | /* The maximum image size we want to support */ |
---|
29 | #define MAX_W 4000 |
---|
30 | #define MAX_H 4000 |
---|
31 | |
---|
32 | /* How does the algorithm work: one point per cell, or two */ |
---|
33 | #define POINTS_PER_CELL 2 |
---|
34 | |
---|
35 | /* The range value for point parameters: X Y, red/green/blue, "strength" |
---|
36 | * Tested values (on Mona Lisa) are: |
---|
37 | * 16 16 5 5 5 2 -> 0.06511725914 |
---|
38 | * 16 16 6 7 6 1 -> 0.05731491348 * |
---|
39 | * 16 16 7 6 6 1 -> 0.06450513783 |
---|
40 | * 14 14 7 7 6 1 -> 0.0637207893 |
---|
41 | * 19 19 6 6 5 1 -> 0.06801999094 */ |
---|
42 | static unsigned int RANGE_X = 16; |
---|
43 | static unsigned int RANGE_Y = 16; |
---|
44 | static unsigned int RANGE_R = 6; |
---|
45 | static unsigned int RANGE_G = 6; |
---|
46 | static unsigned int RANGE_B = 6; |
---|
47 | static unsigned int RANGE_S = 1; |
---|
48 | |
---|
49 | /* |
---|
50 | * These values are computed at runtime |
---|
51 | */ |
---|
52 | |
---|
53 | static float TOTAL_BITS; |
---|
54 | static float HEADER_BITS; |
---|
55 | static float DATA_BITS; |
---|
56 | static float POINT_BITS; |
---|
57 | |
---|
58 | static unsigned int TOTAL_CELLS; |
---|
59 | |
---|
60 | #define RANGE_SY (RANGE_S*RANGE_Y) |
---|
61 | #define RANGE_SYX (RANGE_S*RANGE_Y*RANGE_X) |
---|
62 | #define RANGE_SYXR (RANGE_S*RANGE_Y*RANGE_X*RANGE_R) |
---|
63 | #define RANGE_SYXRG (RANGE_S*RANGE_Y*RANGE_X*RANGE_R*RANGE_G) |
---|
64 | #define RANGE_SYXRGB (RANGE_S*RANGE_Y*RANGE_X*RANGE_R*RANGE_G*RANGE_B) |
---|
65 | |
---|
66 | struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; |
---|
67 | typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; |
---|
68 | typedef std::vector<std::pair<K::Point_2, K::FT> > Point_coordinate_vector; |
---|
69 | |
---|
70 | /* Global aspect ratio */ |
---|
71 | static unsigned int dw, dh; |
---|
72 | |
---|
73 | /* Global point encoding */ |
---|
74 | static uint32_t points[1024]; |
---|
75 | static int npoints = 0; |
---|
76 | |
---|
77 | /* Global triangulation */ |
---|
78 | static Delaunay_triangulation dt; |
---|
79 | |
---|
80 | static unsigned int det_rand(unsigned int mod) |
---|
81 | { |
---|
82 | static unsigned long next = 1; |
---|
83 | next = next * 1103515245 + 12345; |
---|
84 | return ((unsigned)(next / 65536) % 32768) % mod; |
---|
85 | } |
---|
86 | |
---|
87 | static inline int range2int(float val, int range) |
---|
88 | { |
---|
89 | int ret = (int)(val * ((float)range - 0.0001)); |
---|
90 | return ret < 0 ? 0 : ret > range - 1 ? range - 1 : ret; |
---|
91 | } |
---|
92 | |
---|
93 | static inline float int2midrange(int val, int range) |
---|
94 | { |
---|
95 | return (float)(1 + 2 * val) / (float)(2 * range); |
---|
96 | } |
---|
97 | |
---|
98 | static inline float int2fullrange(int val, int range) |
---|
99 | { |
---|
100 | return range > 1 ? (float)val / (float)(range - 1) : 0.0; |
---|
101 | } |
---|
102 | |
---|
103 | static inline void set_point(int index, float x, float y, float r, |
---|
104 | float g, float b, float s) |
---|
105 | { |
---|
106 | int dx = (index / POINTS_PER_CELL) % dw; |
---|
107 | int dy = (index / POINTS_PER_CELL) / dw; |
---|
108 | |
---|
109 | float fx = (x - dx * RANGE_X) / RANGE_X; |
---|
110 | float fy = (y - dy * RANGE_Y) / RANGE_Y; |
---|
111 | |
---|
112 | int is = range2int(s, RANGE_S); |
---|
113 | |
---|
114 | int ix = range2int(fx, RANGE_X); |
---|
115 | int iy = range2int(fy, RANGE_Y); |
---|
116 | |
---|
117 | int ir = range2int(r, RANGE_R); |
---|
118 | int ig = range2int(g, RANGE_G); |
---|
119 | int ib = range2int(b, RANGE_B); |
---|
120 | |
---|
121 | points[index] = is + RANGE_S * (iy + RANGE_Y * (ix + RANGE_X * |
---|
122 | (ib + RANGE_B * (ig + (RANGE_R * ir))))); |
---|
123 | } |
---|
124 | |
---|
125 | static inline void get_point(int index, float *x, float *y, float *r, |
---|
126 | float *g, float *b, float *s) |
---|
127 | { |
---|
128 | uint32_t pt = points[index]; |
---|
129 | |
---|
130 | unsigned int dx = (index / POINTS_PER_CELL) % dw; |
---|
131 | unsigned int dy = (index / POINTS_PER_CELL) / dw; |
---|
132 | |
---|
133 | *s = int2fullrange(pt % RANGE_S, RANGE_S); pt /= RANGE_S; |
---|
134 | |
---|
135 | float fy = int2midrange(pt % RANGE_Y, RANGE_Y); pt /= RANGE_Y; |
---|
136 | float fx = int2midrange(pt % RANGE_X, RANGE_X); pt /= RANGE_X; |
---|
137 | |
---|
138 | *x = (fx + dx) * RANGE_X; |
---|
139 | *y = (fy + dy) * RANGE_Y; |
---|
140 | |
---|
141 | *b = int2midrange(pt % RANGE_R, RANGE_R); pt /= RANGE_R; |
---|
142 | *g = int2midrange(pt % RANGE_G, RANGE_G); pt /= RANGE_G; |
---|
143 | *r = int2midrange(pt % RANGE_B, RANGE_B); pt /= RANGE_B; |
---|
144 | } |
---|
145 | |
---|
146 | static inline float clip(float x, int modulo) |
---|
147 | { |
---|
148 | float mul = (float)modulo + 0.9999; |
---|
149 | int round = (int)(x * mul); |
---|
150 | return (float)round / (float)modulo; |
---|
151 | } |
---|
152 | |
---|
153 | static void add_point(float x, float y, float r, float g, float b, float s) |
---|
154 | { |
---|
155 | set_point(npoints, x, y, r, g, b, s); |
---|
156 | npoints++; |
---|
157 | } |
---|
158 | |
---|
159 | static void add_random_point() |
---|
160 | { |
---|
161 | points[npoints] = det_rand(RANGE_SYXRGB); |
---|
162 | npoints++; |
---|
163 | } |
---|
164 | |
---|
165 | #define NB_OPS 20 |
---|
166 | |
---|
167 | static uint8_t rand_op(void) |
---|
168 | { |
---|
169 | uint8_t x = det_rand(NB_OPS); |
---|
170 | |
---|
171 | /* Randomly ignore statistically less efficient ops */ |
---|
172 | if(x == 0) |
---|
173 | return rand_op(); |
---|
174 | if(x == 1 && (RANGE_S == 1 || det_rand(2))) |
---|
175 | return rand_op(); |
---|
176 | if(x <= 5 && det_rand(2)) |
---|
177 | return rand_op(); |
---|
178 | //if((x < 10 || x > 15) && !det_rand(4)) /* Favour colour changes */ |
---|
179 | // return rand_op(); |
---|
180 | |
---|
181 | return x; |
---|
182 | } |
---|
183 | |
---|
184 | static uint32_t apply_op(uint8_t op, uint32_t val) |
---|
185 | { |
---|
186 | uint32_t rem, ext; |
---|
187 | |
---|
188 | switch(op) |
---|
189 | { |
---|
190 | case 0: /* Flip strength value */ |
---|
191 | case 1: |
---|
192 | /* Statistics show that this helps often, but does not reduce |
---|
193 | * the error significantly. */ |
---|
194 | return val ^ 1; |
---|
195 | case 2: /* Move up; if impossible, down */ |
---|
196 | rem = val % RANGE_S; |
---|
197 | ext = (val / RANGE_S) % RANGE_Y; |
---|
198 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
199 | return (val / RANGE_SY * RANGE_Y + ext) * RANGE_S + rem; |
---|
200 | case 3: /* Move down; if impossible, up */ |
---|
201 | rem = val % RANGE_S; |
---|
202 | ext = (val / RANGE_S) % RANGE_Y; |
---|
203 | ext = ext < RANGE_Y - 1 ? ext + 1 : ext - 1; |
---|
204 | return (val / RANGE_SY * RANGE_Y + ext) * RANGE_S + rem; |
---|
205 | case 4: /* Move left; if impossible, right */ |
---|
206 | rem = val % RANGE_SY; |
---|
207 | ext = (val / RANGE_SY) % RANGE_X; |
---|
208 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
209 | return (val / RANGE_SYX * RANGE_X + ext) * RANGE_SY + rem; |
---|
210 | case 5: /* Move left; if impossible, right */ |
---|
211 | rem = val % RANGE_SY; |
---|
212 | ext = (val / RANGE_SY) % RANGE_X; |
---|
213 | ext = ext < RANGE_X - 1 ? ext + 1 : ext - 1; |
---|
214 | return (val / RANGE_SYX * RANGE_X + ext) * RANGE_SY + rem; |
---|
215 | case 6: /* Corner 1 */ |
---|
216 | return apply_op(2, apply_op(4, val)); |
---|
217 | case 7: /* Corner 2 */ |
---|
218 | return apply_op(2, apply_op(5, val)); |
---|
219 | case 8: /* Corner 3 */ |
---|
220 | return apply_op(3, apply_op(5, val)); |
---|
221 | case 9: /* Corner 4 */ |
---|
222 | return apply_op(3, apply_op(4, val)); |
---|
223 | case 16: /* Double up */ |
---|
224 | return apply_op(2, apply_op(2, val)); |
---|
225 | case 17: /* Double down */ |
---|
226 | return apply_op(3, apply_op(3, val)); |
---|
227 | case 18: /* Double left */ |
---|
228 | return apply_op(4, apply_op(4, val)); |
---|
229 | case 19: /* Double right */ |
---|
230 | return apply_op(5, apply_op(5, val)); |
---|
231 | case 10: /* R-- (or R++) */ |
---|
232 | rem = val % RANGE_SYX; |
---|
233 | ext = (val / RANGE_SYX) % RANGE_R; |
---|
234 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
235 | return (val / RANGE_SYXR * RANGE_R + ext) * RANGE_SYX + rem; |
---|
236 | case 11: /* R++ (or R--) */ |
---|
237 | rem = val % RANGE_SYX; |
---|
238 | ext = (val / RANGE_SYX) % RANGE_R; |
---|
239 | ext = ext < RANGE_R - 1 ? ext + 1 : ext - 1; |
---|
240 | return (val / RANGE_SYXR * RANGE_R + ext) * RANGE_SYX + rem; |
---|
241 | case 12: /* G-- (or G++) */ |
---|
242 | rem = val % RANGE_SYXR; |
---|
243 | ext = (val / RANGE_SYXR) % RANGE_G; |
---|
244 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
245 | return (val / RANGE_SYXRG * RANGE_G + ext) * RANGE_SYXR + rem; |
---|
246 | case 13: /* G++ (or G--) */ |
---|
247 | rem = val % RANGE_SYXR; |
---|
248 | ext = (val / RANGE_SYXR) % RANGE_G; |
---|
249 | ext = ext < RANGE_G - 1 ? ext + 1 : ext - 1; |
---|
250 | return (val / RANGE_SYXRG * RANGE_G + ext) * RANGE_SYXR + rem; |
---|
251 | case 14: /* B-- (or B++) */ |
---|
252 | rem = val % RANGE_SYXRG; |
---|
253 | ext = (val / RANGE_SYXRG) % RANGE_B; |
---|
254 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
255 | return ext * RANGE_SYXRG + rem; |
---|
256 | case 15: /* B++ (or B--) */ |
---|
257 | rem = val % RANGE_SYXRG; |
---|
258 | ext = (val / RANGE_SYXRG) % RANGE_B; |
---|
259 | ext = ext < RANGE_B - 1 ? ext + 1 : ext - 1; |
---|
260 | return ext * RANGE_SYXRG + rem; |
---|
261 | #if 0 |
---|
262 | case 15: /* Brightness-- */ |
---|
263 | return apply_op(9, apply_op(11, apply_op(13, val))); |
---|
264 | case 16: /* Brightness++ */ |
---|
265 | return apply_op(10, apply_op(12, apply_op(14, val))); |
---|
266 | case 17: /* RG-- */ |
---|
267 | return apply_op(9, apply_op(11, val)); |
---|
268 | case 18: /* RG++ */ |
---|
269 | return apply_op(10, apply_op(12, val)); |
---|
270 | case 19: /* GB-- */ |
---|
271 | return apply_op(11, apply_op(13, val)); |
---|
272 | case 20: /* GB++ */ |
---|
273 | return apply_op(12, apply_op(14, val)); |
---|
274 | case 21: /* RB-- */ |
---|
275 | return apply_op(9, apply_op(13, val)); |
---|
276 | case 22: /* RB++ */ |
---|
277 | return apply_op(10, apply_op(14, val)); |
---|
278 | #endif |
---|
279 | default: |
---|
280 | return val; |
---|
281 | } |
---|
282 | } |
---|
283 | |
---|
284 | static void render(pipi_image_t *dst, int rx, int ry, int rw, int rh) |
---|
285 | { |
---|
286 | uint8_t lookup[TOTAL_CELLS * RANGE_X * RANGE_Y]; |
---|
287 | pipi_pixels_t *p = pipi_get_pixels(dst, PIPI_PIXELS_RGBA_F32); |
---|
288 | float *data = (float *)p->pixels; |
---|
289 | int i, x, y; |
---|
290 | |
---|
291 | memset(lookup, 0, sizeof(lookup)); |
---|
292 | dt.clear(); |
---|
293 | for(i = 0; i < npoints; i++) |
---|
294 | { |
---|
295 | float fx, fy, fr, fg, fb, fs; |
---|
296 | get_point(i, &fx, &fy, &fr, &fg, &fb, &fs); |
---|
297 | lookup[(int)fx + dw * RANGE_X * (int)fy] = i; /* Keep link to point */ |
---|
298 | dt.insert(K::Point_2(fx, fy)); |
---|
299 | } |
---|
300 | |
---|
301 | /* Add fake points to close the triangulation */ |
---|
302 | dt.insert(K::Point_2(-p->w, -p->h)); |
---|
303 | dt.insert(K::Point_2(2 * p->w, -p->h)); |
---|
304 | dt.insert(K::Point_2(-p->w, 2 * p->h)); |
---|
305 | dt.insert(K::Point_2(2 * p->w, 2 * p->h)); |
---|
306 | |
---|
307 | for(y = ry; y < ry + rh; y++) |
---|
308 | { |
---|
309 | for(x = rx; x < rx + rw; x++) |
---|
310 | { |
---|
311 | K::Point_2 m(x, y); |
---|
312 | Point_coordinate_vector coords; |
---|
313 | CGAL::Triple< |
---|
314 | std::back_insert_iterator<Point_coordinate_vector>, |
---|
315 | K::FT, bool> result = |
---|
316 | CGAL::natural_neighbor_coordinates_2(dt, m, |
---|
317 | std::back_inserter(coords)); |
---|
318 | |
---|
319 | float r = 0.0f, g = 0.0f, b = 0.0f, norm = 0.0f; |
---|
320 | |
---|
321 | Point_coordinate_vector::iterator it; |
---|
322 | for(it = coords.begin(); it != coords.end(); ++it) |
---|
323 | { |
---|
324 | float fx, fy, fr, fg, fb, fs; |
---|
325 | |
---|
326 | fx = (*it).first.x(); |
---|
327 | fy = (*it).first.y(); |
---|
328 | |
---|
329 | if(fx < 0 || fy < 0 || fx > p->w - 1 || fy > p->h - 1) |
---|
330 | continue; |
---|
331 | |
---|
332 | int index = lookup[(int)fx + dw * RANGE_X * (int)fy]; |
---|
333 | |
---|
334 | get_point(index, &fx, &fy, &fr, &fg, &fb, &fs); |
---|
335 | |
---|
336 | //float k = pow((*it).second * (1.0 + fs), 1.2); |
---|
337 | float k = (*it).second * (1.00f + fs); |
---|
338 | //float k = (*it).second * (0.60f + fs); |
---|
339 | //float k = pow((*it).second, (1.0f + fs)); |
---|
340 | |
---|
341 | r += k * fr; |
---|
342 | g += k * fg; |
---|
343 | b += k * fb; |
---|
344 | norm += k; |
---|
345 | } |
---|
346 | |
---|
347 | data[4 * (x + y * p->w) + 0] = r / norm; |
---|
348 | data[4 * (x + y * p->w) + 1] = g / norm; |
---|
349 | data[4 * (x + y * p->w) + 2] = b / norm; |
---|
350 | data[4 * (x + y * p->w) + 3] = 0.0; |
---|
351 | } |
---|
352 | } |
---|
353 | |
---|
354 | pipi_release_pixels(dst, p); |
---|
355 | } |
---|
356 | |
---|
357 | static void analyse(pipi_image_t *src) |
---|
358 | { |
---|
359 | pipi_pixels_t *p = pipi_get_pixels(src, PIPI_PIXELS_RGBA_F32); |
---|
360 | float *data = (float *)p->pixels; |
---|
361 | |
---|
362 | for(unsigned int dy = 0; dy < dh; dy++) |
---|
363 | for(unsigned int dx = 0; dx < dw; dx++) |
---|
364 | { |
---|
365 | float min = 1.1f, max = -0.1f; |
---|
366 | float total = 0.0; |
---|
367 | int xmin = 0, xmax = 0, ymin = 0, ymax = 0; |
---|
368 | int npixels = 0; |
---|
369 | |
---|
370 | for(unsigned int iy = RANGE_Y * dy; iy < RANGE_Y * (dy + 1); iy++) |
---|
371 | for(unsigned int ix = RANGE_X * dx; ix < RANGE_X * (dx + 1); ix++) |
---|
372 | { |
---|
373 | float lum = 0.0f; |
---|
374 | |
---|
375 | lum += data[4 * (ix + iy * p->w) + 0]; |
---|
376 | lum += data[4 * (ix + iy * p->w) + 1]; |
---|
377 | lum += data[4 * (ix + iy * p->w) + 2]; |
---|
378 | |
---|
379 | if(lum < min) |
---|
380 | { |
---|
381 | min = lum; |
---|
382 | xmin = ix; |
---|
383 | ymin = iy; |
---|
384 | } |
---|
385 | |
---|
386 | if(lum > max) |
---|
387 | { |
---|
388 | max = lum; |
---|
389 | xmax = ix; |
---|
390 | ymax = iy; |
---|
391 | } |
---|
392 | |
---|
393 | total += lum; |
---|
394 | npixels++; |
---|
395 | } |
---|
396 | |
---|
397 | total /= npixels; |
---|
398 | |
---|
399 | float wmin, wmax; |
---|
400 | |
---|
401 | if(total < min + (max - min) / 4) |
---|
402 | wmin = 1.0, wmax = 0.0; |
---|
403 | else if(total < min + (max - min) / 4 * 3) |
---|
404 | wmin = 0.0, wmax = 0.0; |
---|
405 | else |
---|
406 | wmin = 0.0, wmax = 1.0; |
---|
407 | |
---|
408 | #if 0 |
---|
409 | add_random_point(); |
---|
410 | add_random_point(); |
---|
411 | #else |
---|
412 | #if POINTS_PER_CELL == 1 |
---|
413 | if(total < min + (max - min) / 2) |
---|
414 | { |
---|
415 | #endif |
---|
416 | add_point(xmin, ymin, |
---|
417 | data[4 * (xmin + ymin * p->w) + 0], |
---|
418 | data[4 * (xmin + ymin * p->w) + 1], |
---|
419 | data[4 * (xmin + ymin * p->w) + 2], |
---|
420 | wmin); |
---|
421 | #if POINTS_PER_CELL == 1 |
---|
422 | } |
---|
423 | else |
---|
424 | { |
---|
425 | #endif |
---|
426 | add_point(xmax, ymax, |
---|
427 | data[4 * (xmax + ymax * p->w) + 0], |
---|
428 | data[4 * (xmax + ymax * p->w) + 1], |
---|
429 | data[4 * (xmax + ymax * p->w) + 2], |
---|
430 | wmax); |
---|
431 | #if POINTS_PER_CELL == 1 |
---|
432 | } |
---|
433 | #endif |
---|
434 | #endif |
---|
435 | } |
---|
436 | } |
---|
437 | |
---|
438 | int main(int argc, char *argv[]) |
---|
439 | { |
---|
440 | int opstats[2 * NB_OPS]; |
---|
441 | pipi_image_t *src, *tmp, *dst; |
---|
442 | double error = 1.0; |
---|
443 | int width, height, ret = 0; |
---|
444 | |
---|
445 | /* Compute bit allocation */ |
---|
446 | fprintf(stderr, "Available characters: %i\n", NUM_CHARACTERS); |
---|
447 | fprintf(stderr, "Maximum message size: %i\n", MAX_MSG_LEN); |
---|
448 | TOTAL_BITS = MAX_MSG_LEN * logf(NUM_CHARACTERS) / logf(2); |
---|
449 | fprintf(stderr, "Available bits: %f\n", TOTAL_BITS); |
---|
450 | fprintf(stderr, "Maximum image resolution: %ix%i\n", MAX_W, MAX_H); |
---|
451 | HEADER_BITS = logf(MAX_W * MAX_H) / logf(2); |
---|
452 | fprintf(stderr, "Header bits: %f\n", HEADER_BITS); |
---|
453 | DATA_BITS = TOTAL_BITS - HEADER_BITS; |
---|
454 | fprintf(stderr, "Bits available for data: %f\n", DATA_BITS); |
---|
455 | #if POINTS_PER_CELL == 1 |
---|
456 | POINT_BITS = logf(RANGE_SYXRGB) / logf(2); |
---|
457 | #else |
---|
458 | float coord_bits = logf((RANGE_Y * RANGE_X) * (RANGE_Y * RANGE_X + 1) / 2); |
---|
459 | float other_bits = logf(RANGE_R * RANGE_G * RANGE_B * RANGE_S); |
---|
460 | POINT_BITS = (coord_bits + 2 * other_bits) / logf(2); |
---|
461 | #endif |
---|
462 | fprintf(stderr, "Cell bits: %f\n", POINT_BITS); |
---|
463 | TOTAL_CELLS = (int)(DATA_BITS / POINT_BITS); |
---|
464 | fprintf(stderr, "Available cells: %i\n", TOTAL_CELLS); |
---|
465 | fprintf(stderr, "Wasted bits: %f\n", DATA_BITS - POINT_BITS * TOTAL_CELLS); |
---|
466 | |
---|
467 | /* Load image */ |
---|
468 | pipi_set_gamma(1.0); |
---|
469 | src = pipi_load(argv[1]); |
---|
470 | width = pipi_get_image_width(src); |
---|
471 | height = pipi_get_image_height(src); |
---|
472 | |
---|
473 | /* Compute best w/h ratio */ |
---|
474 | dw = 1; dh = TOTAL_CELLS; |
---|
475 | for(unsigned int i = 1; i <= TOTAL_CELLS; i++) |
---|
476 | { |
---|
477 | int j = TOTAL_CELLS / i; |
---|
478 | |
---|
479 | float r = (float)width / (float)height; |
---|
480 | float ir = (float)i / (float)j; |
---|
481 | float dwr = (float)dw / (float)dh; |
---|
482 | |
---|
483 | if(fabs(logf(r / ir)) < fabs(logf(r / dwr))) |
---|
484 | { |
---|
485 | dw = i; |
---|
486 | dh = TOTAL_CELLS / dw; |
---|
487 | } |
---|
488 | } |
---|
489 | while((dh + 1) * dw <= TOTAL_CELLS) dh++; |
---|
490 | while(dw * (dh + 1) <= TOTAL_CELLS) dw++; |
---|
491 | fprintf(stderr, "Chosen image ratio: %i:%i (wasting %i point cells)\n", |
---|
492 | dw, dh, TOTAL_CELLS - dw * dh); |
---|
493 | fprintf(stderr, "Total wasted bits: %f\n", |
---|
494 | DATA_BITS - POINT_BITS * dw * dh); |
---|
495 | |
---|
496 | /* Resize and filter image to better state */ |
---|
497 | tmp = pipi_resize(src, dw * RANGE_X, dh * RANGE_Y); |
---|
498 | pipi_free(src); |
---|
499 | src = pipi_median_ext(tmp, 1, 1); |
---|
500 | pipi_free(tmp); |
---|
501 | |
---|
502 | /* Analyse image */ |
---|
503 | analyse(src); |
---|
504 | |
---|
505 | /* Render what we just computed */ |
---|
506 | tmp = pipi_new(dw * RANGE_X, dh * RANGE_Y); |
---|
507 | render(tmp, 0, 0, dw * RANGE_X, dh * RANGE_Y); |
---|
508 | error = pipi_measure_rmsd(src, tmp); |
---|
509 | |
---|
510 | fprintf(stderr, "Distance: %2.10g\n", error); |
---|
511 | |
---|
512 | memset(opstats, 0, sizeof(opstats)); |
---|
513 | for(int iter = 0, stuck = 0, failures = 0, success = 0; |
---|
514 | /*stuck < 5 && */iter < 10000; |
---|
515 | iter++) |
---|
516 | { |
---|
517 | if(failures > 500) |
---|
518 | { |
---|
519 | stuck++; |
---|
520 | failures = 0; |
---|
521 | } |
---|
522 | |
---|
523 | pipi_image_t *scrap = pipi_copy(tmp); |
---|
524 | |
---|
525 | /* Choose a point at random */ |
---|
526 | int pt = det_rand(npoints); |
---|
527 | uint32_t oldval = points[pt]; |
---|
528 | |
---|
529 | /* Compute the affected image zone */ |
---|
530 | float fx, fy, fr, fg, fb, fs; |
---|
531 | get_point(pt, &fx, &fy, &fr, &fg, &fb, &fs); |
---|
532 | int zonex = (int)fx / RANGE_X - 1; |
---|
533 | int zoney = (int)fy / RANGE_Y - 1; |
---|
534 | int zonew = 3; |
---|
535 | int zoneh = 3; |
---|
536 | if(zonex < 0) { zonex = 0; zonew--; } |
---|
537 | if(zoney < 0) { zoney = 0; zoneh--; } |
---|
538 | if(zonex + zonew >= (int)dw) { zonew--; } |
---|
539 | if(zoney + zoneh >= (int)dh) { zoneh--; } |
---|
540 | |
---|
541 | /* Choose random operations and measure their effect */ |
---|
542 | uint8_t op1 = rand_op(); |
---|
543 | //uint8_t op2 = rand_op(); |
---|
544 | |
---|
545 | uint32_t candidates[3]; |
---|
546 | double besterr = error + 1.0; |
---|
547 | int bestop = -1; |
---|
548 | candidates[0] = apply_op(op1, oldval); |
---|
549 | //candidates[1] = apply_op(op2, oldval); |
---|
550 | //candidates[2] = apply_op(op1, apply_op(op2, oldval)); |
---|
551 | |
---|
552 | for(int i = 0; i < 1; i++) |
---|
553 | //for(int i = 0; i < 3; i++) |
---|
554 | { |
---|
555 | if(oldval == candidates[i]) |
---|
556 | continue; |
---|
557 | |
---|
558 | points[pt] = candidates[i]; |
---|
559 | |
---|
560 | render(scrap, zonex * RANGE_X, zoney * RANGE_Y, |
---|
561 | zonew * RANGE_X, zoneh * RANGE_Y); |
---|
562 | |
---|
563 | double newerr = pipi_measure_rmsd(src, scrap); |
---|
564 | if(newerr < besterr) |
---|
565 | { |
---|
566 | besterr = newerr; |
---|
567 | bestop = i; |
---|
568 | } |
---|
569 | } |
---|
570 | |
---|
571 | opstats[op1 * 2]++; |
---|
572 | //opstats[op2 * 2]++; |
---|
573 | |
---|
574 | if(besterr < error) |
---|
575 | { |
---|
576 | points[pt] = candidates[bestop]; |
---|
577 | /* Redraw image if the last check wasn't the best one */ |
---|
578 | if(bestop != 2) |
---|
579 | render(scrap, zonex * RANGE_X, zoney * RANGE_Y, |
---|
580 | zonew * RANGE_X, zoneh * RANGE_Y); |
---|
581 | |
---|
582 | pipi_free(tmp); |
---|
583 | tmp = scrap; |
---|
584 | //fprintf(stderr, "%08i %2.010g %2.010g after op%i(%i)\n", |
---|
585 | // iter, besterr - error, error, op1, pt); |
---|
586 | fprintf(stderr, "%08i -.%08i %2.010g after op%i(%i)\n", iter, |
---|
587 | (int)((error - besterr) * 100000000), error, op1, pt); |
---|
588 | error = besterr; |
---|
589 | opstats[op1 * 2 + 1]++; |
---|
590 | //opstats[op2 * 2 + 1]++; |
---|
591 | failures = 0; |
---|
592 | success++; |
---|
593 | |
---|
594 | /* Save image! */ |
---|
595 | //char buf[128]; |
---|
596 | //sprintf(buf, "twit%08i.bmp", success); |
---|
597 | //if((success % 10) == 0) |
---|
598 | // pipi_save(tmp, buf); |
---|
599 | } |
---|
600 | else |
---|
601 | { |
---|
602 | pipi_free(scrap); |
---|
603 | points[pt] = oldval; |
---|
604 | failures++; |
---|
605 | } |
---|
606 | } |
---|
607 | |
---|
608 | for(int j = 0; j < 2; j++) |
---|
609 | { |
---|
610 | fprintf(stderr, "operation: "); |
---|
611 | for(int i = NB_OPS / 2 * j; i < NB_OPS / 2 * (j + 1); i++) |
---|
612 | fprintf(stderr, "%4i ", i); |
---|
613 | fprintf(stderr, "\nattempts: "); |
---|
614 | for(int i = NB_OPS / 2 * j; i < NB_OPS / 2 * (j + 1); i++) |
---|
615 | fprintf(stderr, "%4i ", opstats[i * 2]); |
---|
616 | fprintf(stderr, "\nsuccesses: "); |
---|
617 | for(int i = NB_OPS / 2 * j; i < NB_OPS / 2 * (j + 1); i++) |
---|
618 | fprintf(stderr, "%4i ", opstats[i * 2 + 1]); |
---|
619 | fprintf(stderr, "\n"); |
---|
620 | } |
---|
621 | |
---|
622 | fprintf(stderr, "Distance: %2.10g\n", error); |
---|
623 | |
---|
624 | dst = pipi_resize(tmp, width, height); |
---|
625 | pipi_free(tmp); |
---|
626 | |
---|
627 | /* Save image and bail out */ |
---|
628 | pipi_save(dst, "lol.bmp"); |
---|
629 | pipi_free(dst); |
---|
630 | |
---|
631 | return ret; |
---|
632 | } |
---|
633 | |
---|