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 | #define TOTAL_POINTS 138 |
---|
15 | #define POINTS_PER_CELL 2 |
---|
16 | |
---|
17 | #define RANGE_X 16 |
---|
18 | #define RANGE_Y 16 |
---|
19 | #define RANGE_R 5 |
---|
20 | #define RANGE_G 5 |
---|
21 | #define RANGE_B 5 |
---|
22 | #define RANGE_S 1 |
---|
23 | |
---|
24 | #define RANGE_SY (RANGE_S*RANGE_Y) |
---|
25 | #define RANGE_SYX (RANGE_S*RANGE_Y*RANGE_X) |
---|
26 | #define RANGE_SYXR (RANGE_S*RANGE_Y*RANGE_X*RANGE_R) |
---|
27 | #define RANGE_SYXRG (RANGE_S*RANGE_Y*RANGE_X*RANGE_R*RANGE_G) |
---|
28 | #define RANGE_SYXRGB (RANGE_S*RANGE_Y*RANGE_X*RANGE_R*RANGE_G*RANGE_B) |
---|
29 | |
---|
30 | struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; |
---|
31 | typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation; |
---|
32 | typedef std::vector<std::pair<K::Point_2, K::FT> > Point_coordinate_vector; |
---|
33 | |
---|
34 | /* Global aspect ratio */ |
---|
35 | static int dw, dh; |
---|
36 | |
---|
37 | /* Global point encoding */ |
---|
38 | static uint32_t points[1024]; |
---|
39 | static int npoints = 0; |
---|
40 | |
---|
41 | /* Global triangulation */ |
---|
42 | static Delaunay_triangulation dt; |
---|
43 | |
---|
44 | static unsigned int det_rand(unsigned int mod) |
---|
45 | { |
---|
46 | static unsigned long next = 1; |
---|
47 | next = next * 1103515245 + 12345; |
---|
48 | return ((unsigned)(next / 65536) % 32768) % mod; |
---|
49 | } |
---|
50 | |
---|
51 | static inline int float2int(float val, int range) |
---|
52 | { |
---|
53 | int ret = (int)(val * ((float)range - 0.0001)); |
---|
54 | return ret < 0 ? 0 : ret > range - 1? range - 1 : ret; |
---|
55 | } |
---|
56 | |
---|
57 | static inline float int2float(int val, int range) |
---|
58 | { |
---|
59 | return (float)(1 + 2 * val) / (float)(2 * range); |
---|
60 | } |
---|
61 | |
---|
62 | static inline uint32_t set_point(int index, float x, float y, float r, |
---|
63 | float g, float b, float s) |
---|
64 | { |
---|
65 | int dx = (index / POINTS_PER_CELL) % dw; |
---|
66 | int dy = (index / POINTS_PER_CELL) / dw; |
---|
67 | |
---|
68 | float fx = (x - dx * RANGE_X) / RANGE_X; |
---|
69 | float fy = (y - dy * RANGE_Y) / RANGE_Y; |
---|
70 | |
---|
71 | int is = float2int(s, RANGE_S); |
---|
72 | |
---|
73 | int ix = float2int(fx, RANGE_X); |
---|
74 | int iy = float2int(fy, RANGE_Y); |
---|
75 | |
---|
76 | int ir = float2int(r, RANGE_R); |
---|
77 | int ig = float2int(g, RANGE_G); |
---|
78 | int ib = float2int(b, RANGE_B); |
---|
79 | |
---|
80 | points[index] = is + RANGE_S * (iy + RANGE_Y * (ix + RANGE_X * |
---|
81 | (ib + RANGE_B * (ig + (RANGE_R * ir))))); |
---|
82 | } |
---|
83 | |
---|
84 | static inline void get_point(int index, float *x, float *y, float *r, |
---|
85 | float *g, float *b, float *s) |
---|
86 | { |
---|
87 | uint32_t pt = points[index]; |
---|
88 | |
---|
89 | unsigned int dx = (index / POINTS_PER_CELL) % dw; |
---|
90 | unsigned int dy = (index / POINTS_PER_CELL) / dw; |
---|
91 | |
---|
92 | float fs = int2float(pt % RANGE_S, RANGE_S); pt /= RANGE_S; |
---|
93 | |
---|
94 | *s = fs < 0.5 ? 0.0 : 1.0; |
---|
95 | |
---|
96 | float fy = int2float(pt % RANGE_Y, RANGE_Y); pt /= RANGE_Y; |
---|
97 | float fx = int2float(pt % RANGE_X, RANGE_X); pt /= RANGE_X; |
---|
98 | |
---|
99 | *x = (fx + dx) * RANGE_X; |
---|
100 | *y = (fy + dy) * RANGE_Y; |
---|
101 | |
---|
102 | *b = int2float(pt % RANGE_R, RANGE_R); pt /= RANGE_R; |
---|
103 | *g = int2float(pt % RANGE_G, RANGE_G); pt /= RANGE_G; |
---|
104 | *r = int2float(pt % RANGE_B, RANGE_B); pt /= RANGE_B; |
---|
105 | } |
---|
106 | |
---|
107 | static inline float clip(float x, int modulo) |
---|
108 | { |
---|
109 | float mul = (float)modulo + 0.9999; |
---|
110 | int round = (int)(x * mul); |
---|
111 | return (float)round / (float)modulo; |
---|
112 | } |
---|
113 | |
---|
114 | static void add_point(float x, float y, float r, float g, float b, float s) |
---|
115 | { |
---|
116 | set_point(npoints, x, y, r, g, b, s); |
---|
117 | get_point(npoints, &x, &y, &r, &g, &b, &s); |
---|
118 | npoints++; |
---|
119 | } |
---|
120 | |
---|
121 | #define MAX_OPS 15 |
---|
122 | |
---|
123 | static uint32_t apply_op(uint8_t op, uint32_t val) |
---|
124 | { |
---|
125 | uint32_t rem, ext; |
---|
126 | |
---|
127 | switch(op) |
---|
128 | { |
---|
129 | case 0: /* Flip strength value */ |
---|
130 | return val ^ 1; |
---|
131 | case 1: /* Move up; if impossible, down */ |
---|
132 | rem = val % RANGE_S; |
---|
133 | ext = (val / RANGE_S) % RANGE_Y; |
---|
134 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
135 | return (val / RANGE_SY * RANGE_Y + ext) * RANGE_S + rem; |
---|
136 | case 2: /* Move down; if impossible, up */ |
---|
137 | rem = val % RANGE_S; |
---|
138 | ext = (val / RANGE_S) % RANGE_Y; |
---|
139 | ext = ext < RANGE_Y - 1 ? ext + 1 : ext - 1; |
---|
140 | return (val / RANGE_SY * RANGE_Y + ext) * RANGE_S + rem; |
---|
141 | case 3: /* Move left; if impossible, right */ |
---|
142 | rem = val % RANGE_SY; |
---|
143 | ext = (val / RANGE_SY) % RANGE_X; |
---|
144 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
145 | return (val / RANGE_SYX * RANGE_X + ext) * RANGE_SY + rem; |
---|
146 | case 4: /* Move left; if impossible, right */ |
---|
147 | rem = val % RANGE_SY; |
---|
148 | ext = (val / RANGE_SY) % RANGE_X; |
---|
149 | ext = ext < RANGE_X - 1 ? ext + 1 : ext - 1; |
---|
150 | return (val / RANGE_SYX * RANGE_X + ext) * RANGE_SY + rem; |
---|
151 | case 5: /* Corner 1 */ |
---|
152 | return apply_op(1, apply_op(3, val)); |
---|
153 | case 6: /* Corner 2 */ |
---|
154 | return apply_op(1, apply_op(4, val)); |
---|
155 | case 7: /* Corner 3 */ |
---|
156 | return apply_op(2, apply_op(4, val)); |
---|
157 | case 8: /* Corner 4 */ |
---|
158 | return apply_op(2, apply_op(3, val)); |
---|
159 | case 9: /* R-- (or R++) */ |
---|
160 | rem = val % RANGE_SYX; |
---|
161 | ext = (val / RANGE_SYX) % RANGE_R; |
---|
162 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
163 | return (val / RANGE_SYXR * RANGE_R + ext) * RANGE_SYX + rem; |
---|
164 | case 10: /* R++ (or R--) */ |
---|
165 | rem = val % RANGE_SYX; |
---|
166 | ext = (val / RANGE_SYX) % RANGE_R; |
---|
167 | ext = ext < RANGE_R - 1 ? ext + 1 : ext - 1; |
---|
168 | return (val / RANGE_SYXR * RANGE_R + ext) * RANGE_SYX + rem; |
---|
169 | case 11: /* G-- (or G++) */ |
---|
170 | rem = val % RANGE_SYXR; |
---|
171 | ext = (val / RANGE_SYXR) % RANGE_G; |
---|
172 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
173 | return (val / RANGE_SYXRG * RANGE_G + ext) * RANGE_SYXR + rem; |
---|
174 | case 12: /* G++ (or G--) */ |
---|
175 | rem = val % RANGE_SYXR; |
---|
176 | ext = (val / RANGE_SYXR) % RANGE_G; |
---|
177 | ext = ext < RANGE_G - 1 ? ext + 1 : ext - 1; |
---|
178 | return (val / RANGE_SYXRG * RANGE_G + ext) * RANGE_SYXR + rem; |
---|
179 | case 13: /* B-- (or B++) */ |
---|
180 | rem = val % RANGE_SYXRG; |
---|
181 | ext = (val / RANGE_SYXRG) % RANGE_B; |
---|
182 | ext = ext > 0 ? ext - 1 : ext + 1; |
---|
183 | return ext * RANGE_SYXRG + rem; |
---|
184 | case 14: /* B++ (or B--) */ |
---|
185 | rem = val % RANGE_SYXRG; |
---|
186 | ext = (val / RANGE_SYXRG) % RANGE_B; |
---|
187 | ext = ext < RANGE_B - 1 ? ext + 1 : ext - 1; |
---|
188 | return ext * RANGE_SYXRG + rem; |
---|
189 | #if 0 |
---|
190 | case 15: /* Brightness-- */ |
---|
191 | return apply_op(9, apply_op(11, apply_op(13, val))); |
---|
192 | case 16: /* Brightness++ */ |
---|
193 | return apply_op(10, apply_op(12, apply_op(14, val))); |
---|
194 | case 17: /* RG-- */ |
---|
195 | return apply_op(9, apply_op(11, val)); |
---|
196 | case 18: /* RG++ */ |
---|
197 | return apply_op(10, apply_op(12, val)); |
---|
198 | case 19: /* GB-- */ |
---|
199 | return apply_op(11, apply_op(13, val)); |
---|
200 | case 20: /* GB++ */ |
---|
201 | return apply_op(12, apply_op(14, val)); |
---|
202 | case 21: /* RB-- */ |
---|
203 | return apply_op(9, apply_op(13, val)); |
---|
204 | case 22: /* RB++ */ |
---|
205 | return apply_op(10, apply_op(14, val)); |
---|
206 | #endif |
---|
207 | default: |
---|
208 | return val; |
---|
209 | } |
---|
210 | } |
---|
211 | |
---|
212 | static void render(pipi_image_t *dst, int rx, int ry, int rw, int rh) |
---|
213 | { |
---|
214 | uint8_t lookup[TOTAL_POINTS / POINTS_PER_CELL * RANGE_X * RANGE_Y]; |
---|
215 | pipi_pixels_t *p = pipi_get_pixels(dst, PIPI_PIXELS_RGBA_F32); |
---|
216 | float *data = (float *)p->pixels; |
---|
217 | float fx, fy, fr, fg, fb, fs; |
---|
218 | int i, x, y; |
---|
219 | |
---|
220 | memset(lookup, 0, sizeof(lookup)); |
---|
221 | dt.clear(); |
---|
222 | for(i = 0; i < npoints; i++) |
---|
223 | { |
---|
224 | get_point(i, &fx, &fy, &fr, &fg, &fb, &fs); |
---|
225 | lookup[(int)fx + dw * RANGE_X * (int)fy] = i; /* Keep link to point */ |
---|
226 | dt.insert(K::Point_2(fx, fy)); |
---|
227 | } |
---|
228 | |
---|
229 | /* Add fake points to close the triangulation */ |
---|
230 | dt.insert(K::Point_2(-p->w, -p->h)); |
---|
231 | dt.insert(K::Point_2(2 * p->w, -p->h)); |
---|
232 | dt.insert(K::Point_2(-p->w, 2 * p->h)); |
---|
233 | dt.insert(K::Point_2(2 * p->w, 2 * p->h)); |
---|
234 | |
---|
235 | for(y = ry; y < ry + rh; y++) |
---|
236 | { |
---|
237 | for(x = rx; x < rx + rw; x++) |
---|
238 | { |
---|
239 | int i1 = 0, i2 = 0, i3 = 0; |
---|
240 | float d1 = 1000000, d2 = 0, d3 = 0; |
---|
241 | |
---|
242 | K::Point_2 m(x, y); |
---|
243 | Point_coordinate_vector coords; |
---|
244 | CGAL::Triple< |
---|
245 | std::back_insert_iterator<Point_coordinate_vector>, |
---|
246 | K::FT, bool> result = |
---|
247 | CGAL::natural_neighbor_coordinates_2(dt, m, |
---|
248 | std::back_inserter(coords)); |
---|
249 | |
---|
250 | float r = 0.0f, g = 0.0f, b = 0.0f, norm = 0.0f; |
---|
251 | |
---|
252 | Point_coordinate_vector::iterator it; |
---|
253 | for(it = coords.begin(); it != coords.end(); ++it) |
---|
254 | { |
---|
255 | float fx = (*it).first.x(); |
---|
256 | float fy = (*it).first.y(); |
---|
257 | |
---|
258 | if(fx < 0 || fy < 0 || fx > p->w - 1 || fy > p->h - 1) |
---|
259 | continue; |
---|
260 | |
---|
261 | int index = lookup[(int)fx + dw * RANGE_X * (int)fy]; |
---|
262 | |
---|
263 | get_point(index, &fx, &fy, &fr, &fg, &fb, &fs); |
---|
264 | |
---|
265 | float k = (*it).second; |
---|
266 | if(fs > 0.5) k = pow(k, 1.2); |
---|
267 | else k = pow(k, 0.8); |
---|
268 | //float k = (*it).second * (1.0f + fs); |
---|
269 | |
---|
270 | r += k * fr; |
---|
271 | g += k * fg; |
---|
272 | b += k * fb; |
---|
273 | norm += k; |
---|
274 | } |
---|
275 | |
---|
276 | data[4 * (x + y * p->w) + 0] = r / norm; |
---|
277 | data[4 * (x + y * p->w) + 1] = g / norm; |
---|
278 | data[4 * (x + y * p->w) + 2] = b / norm; |
---|
279 | data[4 * (x + y * p->w) + 3] = 0.0; |
---|
280 | } |
---|
281 | } |
---|
282 | |
---|
283 | pipi_release_pixels(dst, p); |
---|
284 | } |
---|
285 | |
---|
286 | static void analyse(pipi_image_t *src) |
---|
287 | { |
---|
288 | pipi_pixels_t *p = pipi_get_pixels(src, PIPI_PIXELS_RGBA_F32); |
---|
289 | float *data = (float *)p->pixels; |
---|
290 | int i; |
---|
291 | |
---|
292 | for(int dy = 0; dy < dh; dy++) |
---|
293 | for(int dx = 0; dx < dw; dx++) |
---|
294 | { |
---|
295 | float min = 1.1f, max = -0.1f; |
---|
296 | float total = 0.0; |
---|
297 | int xmin = 0, xmax = 0, ymin = 0, ymax = 0; |
---|
298 | int npixels = 0; |
---|
299 | |
---|
300 | for(int iy = RANGE_Y * dy; iy < RANGE_Y * (dy + 1); iy++) |
---|
301 | for(int ix = RANGE_X * dx; ix < RANGE_X * (dx + 1); ix++) |
---|
302 | { |
---|
303 | float lum = 0.0f; |
---|
304 | |
---|
305 | lum += data[4 * (ix + iy * p->w) + 0]; |
---|
306 | lum += data[4 * (ix + iy * p->w) + 1]; |
---|
307 | lum += data[4 * (ix + iy * p->w) + 2]; |
---|
308 | |
---|
309 | if(lum < min) |
---|
310 | { |
---|
311 | min = lum; |
---|
312 | xmin = ix; |
---|
313 | ymin = iy; |
---|
314 | } |
---|
315 | |
---|
316 | if(lum > max) |
---|
317 | { |
---|
318 | max = lum; |
---|
319 | xmax = ix; |
---|
320 | ymax = iy; |
---|
321 | } |
---|
322 | |
---|
323 | total += lum; |
---|
324 | npixels++; |
---|
325 | } |
---|
326 | |
---|
327 | total /= npixels; |
---|
328 | |
---|
329 | float wmin, wmax; |
---|
330 | |
---|
331 | if(total < min + (max - min) / 4) |
---|
332 | wmin = 1.0, wmax = 0.0; |
---|
333 | else if(total < min + (max - min) / 4 * 2) |
---|
334 | wmin = 0.0, wmax = 0.0; |
---|
335 | else if(total < min + (max - min) / 4 * 3) |
---|
336 | wmin = 0.0, wmax = 0.0; |
---|
337 | else |
---|
338 | wmin = 0.0, wmax = 1.0; |
---|
339 | |
---|
340 | //wmin = wmax = 1.0; |
---|
341 | //if(total < min + (max - min) /3 ) |
---|
342 | //if((dx + dy) & 1) |
---|
343 | { |
---|
344 | add_point(xmin, ymin, |
---|
345 | data[4 * (xmin + ymin * p->w) + 0], |
---|
346 | data[4 * (xmin + ymin * p->w) + 1], |
---|
347 | data[4 * (xmin + ymin * p->w) + 2], |
---|
348 | wmin); |
---|
349 | } |
---|
350 | //else |
---|
351 | { |
---|
352 | add_point(xmax, ymax, |
---|
353 | data[4 * (xmax + ymax * p->w) + 0], |
---|
354 | data[4 * (xmax + ymax * p->w) + 1], |
---|
355 | data[4 * (xmax + ymax * p->w) + 2], |
---|
356 | wmax); |
---|
357 | } |
---|
358 | } |
---|
359 | } |
---|
360 | |
---|
361 | int main(int argc, char *argv[]) |
---|
362 | { |
---|
363 | int opstats[2 * MAX_OPS]; |
---|
364 | pipi_image_t *src, *tmp, *dst; |
---|
365 | double error = 1.0; |
---|
366 | int width, height, ret = 0; |
---|
367 | |
---|
368 | /* Load image */ |
---|
369 | pipi_set_gamma(1.0); |
---|
370 | src = pipi_load(argv[1]); |
---|
371 | width = pipi_get_image_width(src); |
---|
372 | height = pipi_get_image_height(src); |
---|
373 | |
---|
374 | /* Compute best w/h ratio */ |
---|
375 | dw = 1; |
---|
376 | for(int i = 1; i <= TOTAL_POINTS / POINTS_PER_CELL; i++) |
---|
377 | { |
---|
378 | float r = (float)width / (float)height; |
---|
379 | float ir = (float)i / (float)(TOTAL_POINTS / POINTS_PER_CELL / i); |
---|
380 | float dwr = (float)dw / (float)(TOTAL_POINTS / POINTS_PER_CELL / dw); |
---|
381 | |
---|
382 | if(fabs(logf(r / ir)) < fabs(logf(r / dwr))) |
---|
383 | dw = i; |
---|
384 | } |
---|
385 | dh = TOTAL_POINTS / POINTS_PER_CELL / dw; |
---|
386 | fprintf(stderr, "Chosen image ratio: %i:%i\n", dw, dh); |
---|
387 | |
---|
388 | /* Resize and filter image to better state */ |
---|
389 | tmp = pipi_resize(src, dw * RANGE_X, dh * RANGE_Y); |
---|
390 | pipi_free(src); |
---|
391 | src = pipi_median_ext(tmp, 2, 2); |
---|
392 | |
---|
393 | /* Analyse image */ |
---|
394 | analyse(src); |
---|
395 | |
---|
396 | /* Render what we just computed */ |
---|
397 | tmp = pipi_new(dw * RANGE_X, dh * RANGE_Y); |
---|
398 | render(tmp, 0, 0, dw * RANGE_X, dh * RANGE_Y); |
---|
399 | error = pipi_measure_rmsd(src, tmp); |
---|
400 | |
---|
401 | fprintf(stderr, "Distance: %2.10g\n", error); |
---|
402 | |
---|
403 | memset(opstats, 0, sizeof(opstats)); |
---|
404 | for(int iter = 0, failures = 0; /*failures < 200 &&*/ iter < 3000; iter++) |
---|
405 | { |
---|
406 | uint32_t oldval; |
---|
407 | pipi_image_t *scrap = pipi_copy(tmp); |
---|
408 | |
---|
409 | /* Choose a point at random */ |
---|
410 | int pt = det_rand(npoints); |
---|
411 | oldval = points[pt]; |
---|
412 | |
---|
413 | /* Apply a random operation and measure its effect */ |
---|
414 | uint8_t op = det_rand(MAX_OPS); |
---|
415 | points[pt] = apply_op(op, oldval); |
---|
416 | render(scrap, 0, 0, dw * RANGE_X, dh * RANGE_Y); |
---|
417 | |
---|
418 | opstats[op * 2]++; |
---|
419 | |
---|
420 | double newerr = pipi_measure_rmsd(src, scrap); |
---|
421 | if(newerr < error) |
---|
422 | { |
---|
423 | pipi_free(tmp); |
---|
424 | tmp = scrap; |
---|
425 | error = newerr; |
---|
426 | fprintf(stderr, "%06i %2.010g after op%i(%i)\n", |
---|
427 | iter, error, op, pt); |
---|
428 | opstats[op * 2 + 1]++; |
---|
429 | failures = 0; |
---|
430 | } |
---|
431 | else |
---|
432 | { |
---|
433 | pipi_free(scrap); |
---|
434 | points[pt] = oldval; |
---|
435 | failures++; |
---|
436 | } |
---|
437 | } |
---|
438 | |
---|
439 | fprintf(stderr, "operation: "); |
---|
440 | for(int i = 0; i < MAX_OPS; i++) |
---|
441 | fprintf(stderr, "%3i ", i); |
---|
442 | fprintf(stderr, "\nattempts: "); |
---|
443 | for(int i = 0; i < MAX_OPS; i++) |
---|
444 | fprintf(stderr, "%3i ", opstats[i * 2]); |
---|
445 | fprintf(stderr, "\nsuccesses: "); |
---|
446 | for(int i = 0; i < MAX_OPS; i++) |
---|
447 | fprintf(stderr, "%3i ", opstats[i * 2 + 1]); |
---|
448 | fprintf(stderr, "\n"); |
---|
449 | |
---|
450 | fprintf(stderr, "Distance: %2.10g\n", error); |
---|
451 | |
---|
452 | dst = pipi_resize(tmp, width, height); |
---|
453 | pipi_free(tmp); |
---|
454 | |
---|
455 | /* Save image and bail out */ |
---|
456 | pipi_save(dst, "lol.bmp"); |
---|
457 | pipi_free(dst); |
---|
458 | |
---|
459 | return ret; |
---|
460 | } |
---|
461 | |
---|