source: libpipi/trunk/examples/img2twit.cpp @ 3518

Last change on this file since 3518 was 3518, checked in by Sam Hocevar, 14 years ago

First attempt at a super-compressor for the purpose of sending images to
Twitter, rendering this service slightly more useful. It's still full of
crap, of course.
  • Property svn:keywords set to Id
File size: 13.9 KB
Line 
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
30struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
31typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation;
32typedef std::vector<std::pair<K::Point_2, K::FT> > Point_coordinate_vector;
33
34/* Global aspect ratio */
35static int dw, dh;
36
37/* Global point encoding */
38static uint32_t points[1024];
39static int npoints = 0;
40
41/* Global triangulation */
42static Delaunay_triangulation dt;
43
44static 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
51static 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
57static inline float int2float(int val, int range)
58{
59    return (float)(1 + 2 * val) / (float)(2 * range);
60}
61
62static 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
84static 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
107static 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
114static 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
123static 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
212static 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
286static 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
361int 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
Note: See TracBrowser for help on using the repository browser.