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

Last change on this file since 3161 was 3161, checked in by Sam Hocevar, 12 years ago

libpipi: remove extra newlines in debug messages.

File size: 15.0 KB
Line 
1/*
2 *  libpipi       Pathetic image processing interface 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
21#include <stdio.h>
22#include <stdlib.h>
23#include <string.h>
24#include <math.h>
25#ifndef M_PI
26#   define M_PI 3.14159265358979323846
27#endif
28
29#include "pipi.h"
30#include "pipi_internals.h"
31
32#define R 0
33#define G 1
34#define B 2
35#define X 3
36#define Y 4
37#define A 5
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    double pal[MAXCOLORS][3];
85    double gray[3];
86    double *dark = NULL, *light = NULL;
87    double tmp, min = 1.0, max = 0.0;
88    hull_t *ret = malloc(sizeof(hull_t));
89    int i, j, n;
90
91    debug("");
92    debug("### NEW HULL ###");
93    debug("");
94
95    debug("Analysing %i colors", ncolors);
96
97    for(i = 0; i < ncolors; i++)
98    {
99        pal[i][R] = palette[i * 3];
100        pal[i][G] = palette[i * 3 + 1];
101        pal[i][B] = palette[i * 3 + 2];
102        debug("  [%i] (%g,%g,%g)", i, pal[i][R], pal[i][G], pal[i][B]);
103    }
104
105    /*
106     * 1. Find the darkest and lightest colours
107     */
108    for(i = 0; i < ncolors; i++)
109    {
110        double p = BRIGHT(pal[i]);
111        if(p < min)
112        {
113            dark = pal[i];
114            min = p;
115        }
116        if(p > max)
117        {
118            light = pal[i];
119            max = p;
120        }
121    }
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)",
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    for(n = 0; n <= STEPS; n++)
134    {
135        double pts[MAXCOLORS * (MAXCOLORS - 1) / 2][5];
136        double ptmp[5];
137        double p0[3];
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 ctx, cty, weight;
142        double t = n * 1.0 / STEPS;
143        int npts = 0, left;
144
145        debug("Slice %i/%i", n, STEPS);
146
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)", 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            double yk1;
160            k1[R] = pal[i][R] - p0[R];
161            k1[G] = pal[i][G] - p0[G];
162            k1[B] = pal[i][B] - p0[B];
163            tmp = sqrt(k1[R] * k1[R] + k1[G] * k1[G] + k1[B] * k1[B]);
164
165            /* If k1.y > t.y.y, we don't want this point */
166            yk1 = y[R] * k1[R] + y[G] * k1[G] + y[B] * k1[B];
167            if(yk1 > t * ylen * ylen + EPSILON)
168                continue;
169
170            for(j = 0; j < ncolors; j++)
171            {
172                double k2[3];
173                double yk2, s;
174
175                if(i == j)
176                    continue;
177
178                k2[R] = pal[j][R] - p0[R];
179                k2[G] = pal[j][G] - p0[G];
180                k2[B] = pal[j][B] - p0[B];
181                tmp = sqrt(k2[R] * k2[R] + k2[G] * k2[G] + k2[B] * k2[B]);
182
183                /* If k2.y < t.y.y, we don't want this point */
184                yk2 = y[R] * k2[R] + y[G] * k2[G] + y[B] * k2[B];
185                if(yk2 < t * ylen * ylen - EPSILON)
186                    continue;
187
188                if(yk2 < yk1)
189                    continue;
190
191                s = yk1 == yk2 ? 0.5 : (t * ylen * ylen - yk1) / (yk2 - yk1);
192
193                pts[npts][R] = p0[R] + k1[R] + s * (k2[R] - k1[R]);
194                pts[npts][G] = p0[G] + k1[G] + s * (k2[G] - k1[G]);
195                pts[npts][B] = p0[B] + k1[B] + s * (k2[B] - k1[B]);
196                npts++;
197            }
198        }
199
200        /*
201         * 3.2. Find the barycentre of these points' convex hull. We use
202         *      the Graham Scan technique.
203         */
204
205        /* Make our problem a 2-D problem. */
206        for(i = 0; i < npts; i++)
207        {
208            pts[i][X] = (pts[i][R] - p0[R]) * u[R]
209                           + (pts[i][G] - p0[G]) * u[G]
210                           + (pts[i][B] - p0[B]) * u[B];
211            pts[i][Y] = (pts[i][R] - p0[R]) * v[R]
212                           + (pts[i][G] - p0[G]) * v[G]
213                           + (pts[i][B] - p0[B]) * v[B];
214        }
215
216        /* Find the leftmost point */
217        left = -1;
218        tmp = 10.;
219        for(i = 0; i < npts; i++)
220            if(pts[i][X] < tmp)
221            {
222                left = i;
223                tmp = pts[i][X];
224            }
225        SWAP(pts[0], pts[left]);
226
227        /* Sort the remaining points radially around pts[0]. Bubble sort
228         * is okay for small sizes, I don't care. */
229        for(i = 1; i < npts; i++)
230            for(j = 1; j < npts - i; j++)
231            {
232                double k1 = (pts[j][X] - pts[0][X])
233                              * (pts[j + 1][Y] - pts[0][Y]);
234                double k2 = (pts[j + 1][X] - pts[0][X])
235                              * (pts[j][Y] - pts[0][Y]);
236                if(k1 < k2 - EPSILON)
237                    SWAP(pts[j], pts[j + 1]);
238                else if(k1 < k2 + EPSILON)
239                {
240                    /* Aligned! keep the farthest point */
241                    double ax = pts[j][X] - pts[0][X];
242                    double ay = pts[j][Y] - pts[0][Y];
243                    double bx = pts[j + 1][X] - pts[0][X];
244                    double by = pts[j + 1][Y] - pts[0][Y];
245
246                    if(ax * ax + ay * ay > bx * bx + by * by)
247                        SWAP(pts[j], pts[j + 1]);
248                }
249            }
250
251        /* Remove points not in the convex hull */
252        for(i = 2; i < npts; /* */)
253        {
254            double k1, k2;
255
256            if(i < 2)
257            {
258                i++;
259                continue;
260            }
261
262            k1 = (pts[i - 1][X] - pts[i - 2][X])
263                          * (pts[i][Y] - pts[i - 2][Y]);
264            k2 = (pts[i][X] - pts[i - 2][X])
265                          * (pts[i - 1][Y] - pts[i - 2][Y]);
266            if(k1 <= k2 + EPSILON)
267            {
268                for(j = i - 1; j < npts - 1; j++)
269                    SWAP(pts[j], pts[j + 1]);
270                npts--;
271            }
272            else
273                i++;
274        }
275        /* FIXME: check the last point */
276
277        for(i = 0; i < npts; i++)
278            debug("    2D pt[%i] (%g,%g)", i, pts[i][X], pts[i][Y]);
279
280        /* Compute the barycentre coordinates */
281        ctx = 0.;
282        cty = 0.;
283        weight = 0.;
284        for(i = 2; i < npts; i++)
285        {
286            double abx = pts[i - 1][X] - pts[0][X];
287            double aby = pts[i - 1][Y] - pts[0][Y];
288            double acx = pts[i][X] - pts[0][X];
289            double acy = pts[i][Y] - pts[0][Y];
290            double area;
291            double sqarea = (abx * abx + aby * aby) * (acx * acx + acy * acy)
292                          - (abx * acx + aby * acy) * (abx * acx + aby * acy);
293            if(sqarea <= 0.)
294                continue;
295
296            area = sqrt(sqarea);
297            ctx += area * (abx + acx) / 3;
298            cty += area * (aby + acy) / 3;
299            weight += area;
300        }
301
302        if(weight > EPSILON)
303        {
304            ctx = pts[0][X] + ctx / weight;
305            cty = pts[0][Y] + cty / weight;
306        }
307        else
308        {
309            int right = -1;
310            tmp = -10.;
311            for(i = 0; i < npts; i++)
312                if(pts[i][X] > tmp)
313                {
314                    right = i;
315                    tmp = pts[i][X];
316                }
317            ctx = 0.5 * (pts[0][X] + pts[right][X]);
318            cty = 0.5 * (pts[0][Y] + pts[right][Y]);
319        }
320
321        debug("    2D bary (%g,%g)", ctx, cty);
322
323        /*
324         * 3.3. Store the barycentre and convex hull information.
325         */
326
327        ret->bary[n][R] = p0[R] + ctx * u[R] + cty * v[R];
328        ret->bary[n][G] = p0[G] + ctx * u[G] + cty * v[G];
329        ret->bary[n][B] = p0[B] + ctx * u[B] + cty * v[B];
330
331        for(i = 0; i < npts; i++)
332        {
333            ret->pts[n][i][R] = pts[i][R];
334            ret->pts[n][i][G] = pts[i][G];
335            ret->pts[n][i][B] = pts[i][B];
336            ret->pts[n][i][X] = pts[i][X] - ctx;
337            ret->pts[n][i][Y] = pts[i][Y] - cty;
338            ret->pts[n][i][A] = atan2(pts[i][Y] - cty, pts[i][X] - ctx);
339
340            debug("  3D pt[%i] (%g,%g,%g) angle %g",
341                  i, pts[i][R], pts[i][G], pts[i][B], ret->pts[n][i][A]);
342        }
343        ret->hullsize[n] = npts;
344
345        debug("  3D bary (%g,%g,%g)",
346              ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
347    }
348
349    return ret;
350}
351
352
353pipi_image_t *pipi_reduce(pipi_image_t *src,
354                          int ncolors, double const *palette)
355{
356    static double const rgbpal[] =
357    {
358        0, 0, 0,  0, 0, 1,  0, 1, 0,  0, 1, 1,
359        1, 0, 0,  1, 0, 1,  1, 1, 0,  1, 1, 1,
360    };
361
362    pipi_image_t *dst;
363    pipi_pixels_t *srcp, *dstp;
364    float *srcdata, *dstdata;
365    hull_t *rgbhull, *myhull;
366    int i, j, w, h;
367
368    init_uv();
369
370    rgbhull = compute_hull(8, rgbpal);
371    myhull = compute_hull(ncolors, palette);
372
373    /*
374     * 4. Load image and change its palette.
375     */
376
377    debug("");
378    debug("### PROCESSING IMAGE ###");
379    debug("");
380
381    srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
382    srcdata = (float *)srcp->pixels;
383
384    w = srcp->w;
385    h = srcp->h;
386
387    dst = pipi_new(w, h);
388    dstp = pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
389    dstdata = (float *)dstp->pixels;
390
391    for(j = 0; j < h; j++)
392        for(i = 0; i < w; i++)
393        {
394            double p[3];
395            double xp, yp, angle, xa, ya, xb, yb, t, s;
396            int slice, n, count;
397
398            /* FIXME: Imlib fucks up the RGB order. */
399            p[B] = srcdata[4 * (j * w + i)];
400            p[G] = srcdata[4 * (j * w + i) + 1];
401            p[R] = srcdata[4 * (j * w + i) + 2];
402
403            debug("Pixel +%i+%i (%g,%g,%g)", i, j, p[R], p[G], p[B]);
404
405            slice = (int)(BRIGHT(p) * STEPS + 0.5);
406
407            debug("  slice %i", slice);
408
409            /* Convert to 2D. The origin is the slice's barycentre. */
410            xp = (p[R] - rgbhull->bary[slice][R]) * u[R]
411               + (p[G] - rgbhull->bary[slice][G]) * u[G]
412               + (p[B] - rgbhull->bary[slice][B]) * u[B];
413            yp = (p[R] - rgbhull->bary[slice][R]) * v[R]
414               + (p[G] - rgbhull->bary[slice][G]) * v[G]
415               + (p[B] - rgbhull->bary[slice][B]) * v[B];
416
417            debug("    2D pt (%g,%g)", xp, yp);
418
419            /* 1. find the excentricity in RGB space. There is an easier
420             * way to do this, which is to find the intersection of our
421             * line with the RGB cube itself, but we'd lose the possibility
422             * of having an original colour space other than RGB. */
423
424            /* First, find the relevant triangle. */
425            count = rgbhull->hullsize[slice];
426            angle = atan2(yp, xp);
427            for(n = 0; n < count; n++)
428            {
429                double a1 = rgbhull->pts[slice][n][A];
430                double a2 = rgbhull->pts[slice][(n + 1) % count][A];
431                if(a1 > a2)
432                {
433                    if(angle >= a1)
434                        a2 += 2 * M_PI;
435                    else
436                        a1 -= 2 * M_PI;
437                }
438                if(angle >= a1 && angle <= a2)
439                    break;
440            }
441
442            /* Now compute the distance to the triangle's edge. If the edge
443             * intersection is M, then t is such as P = t.M (can be zero) */
444            xa = rgbhull->pts[slice][n % count][X];
445            ya = rgbhull->pts[slice][n % count][Y];
446            xb = rgbhull->pts[slice][(n + 1) % count][X];
447            yb = rgbhull->pts[slice][(n + 1) % count][Y];
448            t = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
449
450            if(t > 1.0)
451                t = 1.0;
452
453            debug("    best RGB %g (%g,%g) (%g,%g)", t, xa, ya, xb, yb);
454
455            /* 2. apply the excentricity in reduced space. */
456
457            count = myhull->hullsize[slice];
458            for(n = 0; n < count; n++)
459            {
460                double a1 = myhull->pts[slice][n][A];
461                double a2 = myhull->pts[slice][(n + 1) % count][A];
462                if(a1 > a2)
463                {
464                    if(angle >= a1)
465                        a2 += 2 * M_PI;
466                    else
467                        a1 -= 2 * M_PI;
468                }
469                if(angle >= a1 && angle <= a2)
470                    break;
471            }
472
473            /* If the edge intersection is M', s is such as P = s.M'. We
474             * want P' = t.M' = t.P/s  */
475            xa = myhull->pts[slice][n % count][X];
476            ya = myhull->pts[slice][n % count][Y];
477            xb = myhull->pts[slice][(n + 1) % count][X];
478            yb = myhull->pts[slice][(n + 1) % count][Y];
479            s = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
480
481            debug("    best custom %g (%g,%g) (%g,%g)", s, xa, ya, xb, yb);
482
483            if(s > 0)
484            {
485                xp *= t / s;
486                yp *= t / s;
487            }
488
489            p[R] = myhull->bary[slice][R] + xp * u[R] + yp * v[R];
490            p[G] = myhull->bary[slice][G] + xp * u[G] + yp * v[G];
491            p[B] = myhull->bary[slice][B] + xp * u[B] + yp * v[B];
492
493            /* Clipping should not be necessary, but the above code
494             * is unfortunately not perfect. */
495            if(p[R] < 0.0) p[R] = 0.0; else if(p[R] > 1.0) p[R] = 1.0;
496            if(p[G] < 0.0) p[G] = 0.0; else if(p[G] > 1.0) p[G] = 1.0;
497            if(p[B] < 0.0) p[B] = 0.0; else if(p[B] > 1.0) p[B] = 1.0;
498
499            dstdata[4 * (j * w + i)] = p[B];
500            dstdata[4 * (j * w + i) + 1] = p[G];
501            dstdata[4 * (j * w + i) + 2] = p[R];
502        }
503
504    free(rgbhull);
505    free(myhull);
506
507    return dst;
508}
509
Note: See TracBrowser for help on using the repository browser.