source: libpipi/trunk/examples/img2rubik.c @ 2734

Last change on this file since 2734 was 2734, checked in by sam, 6 years ago
  • img2rubik.c: improve stability by clipping t.
File size: 14.6 KB
Line 
1#include "config.h"
2#include "common.h"
3
4#include <stdio.h>
5#include <stdlib.h>
6#include <string.h>
7#include <math.h>
8
9#include <pipi.h>
10
11#define R 0
12#define G 1
13#define B 2
14#define X 3
15#define Y 4
16#define A 5
17
18//#define debug printf
19#define debug(...) /* */
20
21#define BRIGHT(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
22
23#define MAXCOLORS 16
24#define STEPS 256
25#define EPSILON (0.000001)
26
27typedef struct
28{
29    double pts[STEPS + 1][MAXCOLORS * (MAXCOLORS - 1) / 2][6];
30    int hullsize[STEPS + 1];
31    double bary[STEPS + 1][3];
32}
33hull_t;
34
35static double const y[3] = { .299, .587, .114 };
36static double u[3], v[3];
37static int ylen;
38
39/*
40 * Find two base vectors for the chrominance planes.
41 */
42static void init_uv(void)
43{
44    double tmp;
45
46    ylen = sqrt(y[R] * y[R] + y[G] * y[G] + y[B] * y[B]);
47
48    u[R] = y[1];
49    u[G] = -y[0];
50    u[B] = 0;
51    tmp = sqrt(u[R] * u[R] + u[G] * u[G] + u[B] * u[B]);
52    u[R] /= tmp; u[G] /= tmp; u[B] /= tmp;
53
54    v[R] = y[G] * u[B] - y[B] * u[G];
55    v[G] = y[B] * u[R] - y[R] * u[B];
56    v[B] = y[R] * u[G] - y[G] * u[R];
57    tmp = sqrt(v[R] * v[R] + v[G] * v[G] + v[B] * v[B]);
58    v[R] /= tmp; v[G] /= tmp; v[B] /= tmp;
59}
60
61/*
62 * Compute the convex hull of a given palette.
63 */
64static hull_t *compute_hull(int ncolors, double const *palette)
65{
66    hull_t *ret = malloc(sizeof(hull_t));
67    double tmp;
68    int i, j;
69
70    debug("\n### NEW HULL ###\n\n");
71
72    debug("Analysing %i colors\n", ncolors);
73
74    double pal[ncolors][3];
75    for(i = 0; i < ncolors; i++)
76    {
77        pal[i][R] = palette[i * 3];
78        pal[i][G] = palette[i * 3 + 1];
79        pal[i][B] = palette[i * 3 + 2];
80        debug("  [%i] (%g,%g,%g)\n", i, pal[i][R], pal[i][G], pal[i][B]);
81    }
82
83    /*
84     * 1. Find the darkest and lightest colours
85     */
86    double *dark = NULL, *light = NULL;
87    double min = 1.0, max = 0.0;
88    for(i = 0; i < ncolors; i++)
89    {
90        double p = BRIGHT(pal[i]);
91        if(p < min)
92        {
93            dark = pal[i];
94            min = p;
95        }
96        if(p > max)
97        {
98            light = pal[i];
99            max = p;
100        }
101    }
102
103    double gray[3];
104
105    gray[R] = light[R] - dark[R];
106    gray[G] = light[G] - dark[G];
107    gray[B] = light[B] - dark[B];
108
109    debug("  gray axis (%g,%g,%g) - (%g,%g,%g)\n",
110          dark[R], dark[G], dark[B], light[R], light[G], light[B]);
111
112    /*
113     * 3. Browse the grey axis and do stuff
114     */
115    int n;
116    for(n = 0; n <= STEPS; n++)
117    {
118        double pts[ncolors * (ncolors - 1) / 2][5];
119        double ptmp[5];
120#define SWAP(p1,p2) do { memcpy(ptmp, p1, sizeof(ptmp)); \
121                         memcpy(p1, p2, sizeof(ptmp)); \
122                         memcpy(p2, ptmp, sizeof(ptmp)); } while(0)
123        double t = n * 1.0 / STEPS;
124        int npts = 0;
125
126        debug("Slice %i/%i\n", n, STEPS);
127
128        double p0[3];
129        p0[R] = dark[R] + t * gray[R];
130        p0[G] = dark[G] + t * gray[G];
131        p0[B] = dark[B] + t * gray[B];
132
133        debug("  3D gray (%g,%g,%g)\n", p0[R], p0[G], p0[B]);
134
135        /*
136         * 3.1. Find all edges that intersect the t.y + (u,v) plane
137         */
138        for(i = 0; i < ncolors; i++)
139        {
140            double k1[3];
141            k1[R] = pal[i][R] - p0[R];
142            k1[G] = pal[i][G] - p0[G];
143            k1[B] = pal[i][B] - p0[B];
144            tmp = sqrt(k1[R] * k1[R] + k1[G] * k1[G] + k1[B] * k1[B]);
145
146            /* If k1.y > t.y.y, we don't want this point */
147            double yk1 = y[R] * k1[R] + y[G] * k1[G] + y[B] * k1[B];
148            if(yk1 > t * ylen * ylen + EPSILON)
149                continue;
150
151            for(j = 0; j < ncolors; j++)
152            {
153                if(i == j)
154                    continue;
155
156                double k2[3];
157                k2[R] = pal[j][R] - p0[R];
158                k2[G] = pal[j][G] - p0[G];
159                k2[B] = pal[j][B] - p0[B];
160                tmp = sqrt(k2[R] * k2[R] + k2[G] * k2[G] + k2[B] * k2[B]);
161
162                /* If k2.y < t.y.y, we don't want this point */
163                double yk2 = y[R] * k2[R] + y[G] * k2[G] + y[B] * k2[B];
164                if(yk2 < t * ylen * ylen - EPSILON)
165                    continue;
166
167                if(yk2 < yk1)
168                    continue;
169
170                double s = yk1 == yk2 ?
171                           0.5 : (t * ylen * ylen - yk1) / (yk2 - yk1);
172
173                pts[npts][R] = p0[R] + k1[R] + s * (k2[R] - k1[R]);
174                pts[npts][G] = p0[G] + k1[G] + s * (k2[G] - k1[G]);
175                pts[npts][B] = p0[B] + k1[B] + s * (k2[B] - k1[B]);
176                npts++;
177            }
178        }
179
180        /*
181         * 3.2. Find the barycentre of these points' convex hull. We use
182         *      the Graham Scan technique.
183         */
184
185        /* Make our problem a 2-D problem. */
186        for(i = 0; i < npts; i++)
187        {
188            pts[i][X] = (pts[i][R] - p0[R]) * u[R]
189                           + (pts[i][G] - p0[G]) * u[G]
190                           + (pts[i][B] - p0[B]) * u[B];
191            pts[i][Y] = (pts[i][R] - p0[R]) * v[R]
192                           + (pts[i][G] - p0[G]) * v[G]
193                           + (pts[i][B] - p0[B]) * v[B];
194        }
195
196        /* Find the leftmost point */
197        int left = -1;
198        tmp = 10.;
199        for(i = 0; i < npts; i++)
200            if(pts[i][X] < tmp)
201            {
202                left = i;
203                tmp = pts[i][X];
204            }
205        SWAP(pts[0], pts[left]);
206
207        /* Sort the remaining points radially around pts[0]. Bubble sort
208         * is okay for small sizes, I don't care. */
209        for(i = 1; i < npts; i++)
210            for(j = 1; j < npts - i; j++)
211            {
212                double k1 = (pts[j][X] - pts[0][X])
213                              * (pts[j + 1][Y] - pts[0][Y]);
214                double k2 = (pts[j + 1][X] - pts[0][X])
215                              * (pts[j][Y] - pts[0][Y]);
216                if(k1 < k2 - EPSILON)
217                    SWAP(pts[j], pts[j + 1]);
218                else if(k1 < k2 + EPSILON)
219                {
220                    /* Aligned! keep the farthest point */
221                    double ax = pts[j][X] - pts[0][X];
222                    double ay = pts[j][Y] - pts[0][Y];
223                    double bx = pts[j + 1][X] - pts[0][X];
224                    double by = pts[j + 1][Y] - pts[0][Y];
225
226                    if(ax * ax + ay * ay > bx * bx + by * by)
227                        SWAP(pts[j], pts[j + 1]);
228                }
229            }
230
231        /* Remove points not in the convex hull */
232        for(i = 2; i < npts; /* */)
233        {
234            if(i < 2)
235            {
236                i++;
237                continue;
238            }
239
240            double k1 = (pts[i - 1][X] - pts[i - 2][X])
241                          * (pts[i][Y] - pts[i - 2][Y]);
242            double k2 = (pts[i][X] - pts[i - 2][X])
243                          * (pts[i - 1][Y] - pts[i - 2][Y]);
244            if(k1 <= k2 + EPSILON)
245            {
246                for(j = i - 1; j < npts - 1; j++)
247                    SWAP(pts[j], pts[j + 1]);
248                npts--;
249            }
250            else
251                i++;
252        }
253        /* FIXME: check the last point */
254
255        for(i = 0; i < npts; i++)
256            debug("    2D pt[%i] (%g,%g)\n", i, pts[i][X], pts[i][Y]);
257
258        /* Compute the barycentre coordinates */
259        double ctx = 0., cty = 0., weight = 0.;
260        for(i = 2; i < npts; i++)
261        {
262            double abx = pts[i - 1][X] - pts[0][X];
263            double aby = pts[i - 1][Y] - pts[0][Y];
264            double acx = pts[i][X] - pts[0][X];
265            double acy = pts[i][Y] - pts[0][Y];
266            double sqarea = (abx * abx + aby * aby) * (acx * acx + acy * acy)
267                          - (abx * acx + aby * acy) * (abx * acx + aby * acy);
268            if(sqarea <= 0.)
269                continue;
270
271            double area = sqrt(sqarea);
272            ctx += area * (abx + acx) / 3;
273            cty += area * (aby + acy) / 3;
274            weight += area;
275        }
276
277        if(weight > EPSILON)
278        {
279            ctx = pts[0][X] + ctx / weight;
280            cty = pts[0][Y] + cty / weight;
281        }
282        else
283        {
284            int right = -1;
285            tmp = -10.;
286            for(i = 0; i < npts; i++)
287                if(pts[i][X] > tmp)
288                {
289                    right = i;
290                    tmp = pts[i][X];
291                }
292            ctx = 0.5 * (pts[0][X] + pts[right][X]);
293            cty = 0.5 * (pts[0][Y] + pts[right][Y]);
294        }
295
296        debug("    2D bary (%g,%g)\n", ctx, cty);
297
298        /*
299         * 3.3. Store the barycentre and convex hull information.
300         */
301
302        ret->bary[n][R] = p0[R] + ctx * u[R] + cty * v[R];
303        ret->bary[n][G] = p0[G] + ctx * u[G] + cty * v[G];
304        ret->bary[n][B] = p0[B] + ctx * u[B] + cty * v[B];
305
306        for(i = 0; i < npts; i++)
307        {
308            ret->pts[n][i][R] = pts[i][R];
309            ret->pts[n][i][G] = pts[i][G];
310            ret->pts[n][i][B] = pts[i][B];
311            ret->pts[n][i][X] = pts[i][X] - ctx;
312            ret->pts[n][i][Y] = pts[i][Y] - cty;
313            ret->pts[n][i][A] = atan2(pts[i][Y] - cty, pts[i][X] - ctx);
314
315            debug("  3D pt[%i] (%g,%g,%g) angle %g\n",
316                  i, pts[i][R], pts[i][G], pts[i][B], ret->pts[n][i][A]);
317        }
318        ret->hullsize[n] = npts;
319
320        debug("  3D bary (%g,%g,%g)\n",
321              ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
322    }
323
324    return ret;
325}
326
327
328int main(int argc, char *argv[])
329{
330    static double const rgbpal[] =
331    {
332        0, 0, 0,  0, 0, 1,  0, 1, 0,  0, 1, 1,
333        1, 0, 0,  1, 0, 1,  1, 1, 0,  1, 1, 1,
334    };
335
336    static double const mypal[] =
337    {
338        0.900, 0.001, 0.001, /* red */
339        0.001, 0.800, 0.001, /* green */
340        0.005, 0.001, 0.500, /* blue */
341        0.900, 0.900, 0.900, /* white */
342        0.900, 0.900, 0.001, /* yellow */
343        0.800, 0.400, 0.001, /* orange */
344    };
345
346    int i, j;
347
348    init_uv();
349
350    hull_t *rgbhull = compute_hull(8, rgbpal);
351    hull_t *myhull = compute_hull(6, mypal);
352
353    /*
354     * 4. Load image and change its palette.
355     */
356
357    debug("\n### PROCESSING IMAGE ###\n\n");
358
359    pipi_image_t *src = pipi_load(argv[1]);
360    pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
361    float *srcdata = (float *)srcp->pixels;
362
363    int w = srcp->w, h = srcp->h;
364
365    pipi_image_t *dst = pipi_new(w, h);
366    pipi_pixels_t *dstp = pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
367    float *dstdata = (float *)dstp->pixels;
368
369    for(j = 0; j < h; j++)
370        for(i = 0; i < w; i++)
371        {
372            double p[3];
373            /* FIXME: Imlib fucks up the RGB order. */
374            p[B] = srcdata[4 * (j * w + i)];
375            p[G] = srcdata[4 * (j * w + i) + 1];
376            p[R] = srcdata[4 * (j * w + i) + 2];
377
378            debug("Pixel +%i+%i (%g,%g,%g)\n", i, j, p[R], p[G], p[B]);
379
380            int slice = (int)(BRIGHT(p) * STEPS + 0.5);
381
382            debug("  slice %i\n", slice);
383
384            /* Convert to 2D. The origin is the slice's barycentre. */
385            double xp = (p[R] - rgbhull->bary[slice][R]) * u[R]
386                      + (p[G] - rgbhull->bary[slice][G]) * u[G]
387                      + (p[B] - rgbhull->bary[slice][B]) * u[B];
388            double yp = (p[R] - rgbhull->bary[slice][R]) * v[R]
389                      + (p[G] - rgbhull->bary[slice][G]) * v[G]
390                      + (p[B] - rgbhull->bary[slice][B]) * v[B];
391
392            debug("    2D pt (%g,%g)\n", xp, yp);
393
394            /* 1. find the excentricity in RGB space. There is an easier
395             * way to do this, which is to find the intersection of our
396             * line with the RGB cube itself, but we'd lose the possibility
397             * of having an original colour space other than RGB. */
398
399            /* First, find the relevant triangle. */
400            int n, count = rgbhull->hullsize[slice];
401            double angle = atan2(yp, xp);
402            for(n = 0; n < count; n++)
403            {
404                double a1 = rgbhull->pts[slice][n][A];
405                double a2 = rgbhull->pts[slice][(n + 1) % count][A];
406                if(a1 > a2)
407                {
408                    if(angle >= a1)
409                        a2 += 2 * M_PI;
410                    else
411                        a1 -= 2 * M_PI;
412                }
413                if(angle >= a1 && angle <= a2)
414                    break;
415            }
416
417            /* Now compute the distance to the triangle's edge. If the edge
418             * intersection is M, then t is such as P = t.M (can be zero) */
419            double xa = rgbhull->pts[slice][n % count][X];
420            double ya = rgbhull->pts[slice][n % count][Y];
421            double xb = rgbhull->pts[slice][(n + 1) % count][X];
422            double yb = rgbhull->pts[slice][(n + 1) % count][Y];
423            double t = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
424
425            if(t > 1.0)
426                t = 1.0;
427
428            debug("    best RGB %g (%g,%g) (%g,%g)\n", t, xa, ya, xb, yb);
429
430            /* 2. apply the excentricity in reduced space. */
431
432            count = myhull->hullsize[slice];
433            for(n = 0; n < count; n++)
434            {
435                double a1 = myhull->pts[slice][n][A];
436                double a2 = myhull->pts[slice][(n + 1) % count][A];
437                if(a1 > a2)
438                {
439                    if(angle >= a1)
440                        a2 += 2 * M_PI;
441                    else
442                        a1 -= 2 * M_PI;
443                }
444                if(angle >= a1 && angle <= a2)
445                    break;
446            }
447
448            /* If the edge intersection is M', s is such as P = s.M'. We
449             * want P' = t.M' = t.P/s  */
450            xa = myhull->pts[slice][n % count][X];
451            ya = myhull->pts[slice][n % count][Y];
452            xb = myhull->pts[slice][(n + 1) % count][X];
453            yb = myhull->pts[slice][(n + 1) % count][Y];
454            double s = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
455
456            debug("    best custom %g (%g,%g) (%g,%g)\n", s, xa, ya, xb, yb);
457
458            if(s > 0)
459            {
460                xp *= t / s;
461                yp *= t / s;
462            }
463
464            p[R] = myhull->bary[slice][R] + xp * u[R] + yp * v[R];
465            p[G] = myhull->bary[slice][G] + xp * u[G] + yp * v[G];
466            p[B] = myhull->bary[slice][B] + xp * u[B] + yp * v[B];
467
468            /* Clipping should not be necessary, but the above code
469             * is unfortunately not perfect. */
470            if(p[R] < 0.0) p[R] = 0.0; else if(p[R] > 1.0) p[R] = 1.0;
471            if(p[G] < 0.0) p[G] = 0.0; else if(p[G] > 1.0) p[G] = 1.0;
472            if(p[B] < 0.0) p[B] = 0.0; else if(p[B] > 1.0) p[B] = 1.0;
473
474            dstdata[4 * (j * w + i)] = p[B];
475            dstdata[4 * (j * w + i) + 1] = p[G];
476            dstdata[4 * (j * w + i) + 2] = p[R];
477        }
478
479    free(rgbhull);
480    free(myhull);
481
482    pipi_save(dst, argv[2]);
483
484    return 0;
485}
486
Note: See TracBrowser for help on using the repository browser.