source: libpipi/trunk/pipi/quantize/reduce.c @ 2736

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