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

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

Start writing Visual Studio projects.

File size: 14.9 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("\n### NEW HULL ###\n\n");
92
93    debug("Analysing %i colors\n", ncolors);
94
95    for(i = 0; i < ncolors; i++)
96    {
97        pal[i][R] = palette[i * 3];
98        pal[i][G] = palette[i * 3 + 1];
99        pal[i][B] = palette[i * 3 + 2];
100        debug("  [%i] (%g,%g,%g)\n", i, pal[i][R], pal[i][G], pal[i][B]);
101    }
102
103    /*
104     * 1. Find the darkest and lightest colours
105     */
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    gray[R] = light[R] - dark[R];
122    gray[G] = light[G] - dark[G];
123    gray[B] = light[B] - dark[B];
124
125    debug("  gray axis (%g,%g,%g) - (%g,%g,%g)\n",
126          dark[R], dark[G], dark[B], light[R], light[G], light[B]);
127
128    /*
129     * 3. Browse the grey axis and do stuff
130     */
131    for(n = 0; n <= STEPS; n++)
132    {
133        double pts[MAXCOLORS * (MAXCOLORS - 1) / 2][5];
134        double ptmp[5];
135        double p0[3];
136#define SWAP(p1,p2) do { memcpy(ptmp, p1, sizeof(ptmp)); \
137                         memcpy(p1, p2, sizeof(ptmp)); \
138                         memcpy(p2, ptmp, sizeof(ptmp)); } while(0)
139        double ctx, cty, weight;
140        double t = n * 1.0 / STEPS;
141        int npts = 0, left;
142
143        debug("Slice %i/%i\n", n, STEPS);
144
145        p0[R] = dark[R] + t * gray[R];
146        p0[G] = dark[G] + t * gray[G];
147        p0[B] = dark[B] + t * gray[B];
148
149        debug("  3D gray (%g,%g,%g)\n", p0[R], p0[G], p0[B]);
150
151        /*
152         * 3.1. Find all edges that intersect the t.y + (u,v) plane
153         */
154        for(i = 0; i < ncolors; i++)
155        {
156            double k1[3];
157                        double yk1;
158            k1[R] = pal[i][R] - p0[R];
159            k1[G] = pal[i][G] - p0[G];
160            k1[B] = pal[i][B] - p0[B];
161            tmp = sqrt(k1[R] * k1[R] + k1[G] * k1[G] + k1[B] * k1[B]);
162
163            /* If k1.y > t.y.y, we don't want this point */
164            yk1 = y[R] * k1[R] + y[G] * k1[G] + y[B] * k1[B];
165            if(yk1 > t * ylen * ylen + EPSILON)
166                continue;
167
168            for(j = 0; j < ncolors; j++)
169            {
170                double k2[3];
171                double yk2, s;
172
173                                if(i == j)
174                    continue;
175
176                k2[R] = pal[j][R] - p0[R];
177                k2[G] = pal[j][G] - p0[G];
178                k2[B] = pal[j][B] - p0[B];
179                tmp = sqrt(k2[R] * k2[R] + k2[G] * k2[G] + k2[B] * k2[B]);
180
181                /* If k2.y < t.y.y, we don't want this point */
182                yk2 = y[R] * k2[R] + y[G] * k2[G] + y[B] * k2[B];
183                if(yk2 < t * ylen * ylen - EPSILON)
184                    continue;
185
186                if(yk2 < yk1)
187                    continue;
188
189                s = yk1 == yk2 ? 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        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                        double k1, k2;
253
254            if(i < 2)
255            {
256                i++;
257                continue;
258            }
259
260            k1 = (pts[i - 1][X] - pts[i - 2][X])
261                          * (pts[i][Y] - pts[i - 2][Y]);
262            k2 = (pts[i][X] - pts[i - 2][X])
263                          * (pts[i - 1][Y] - pts[i - 2][Y]);
264            if(k1 <= k2 + EPSILON)
265            {
266                for(j = i - 1; j < npts - 1; j++)
267                    SWAP(pts[j], pts[j + 1]);
268                npts--;
269            }
270            else
271                i++;
272        }
273        /* FIXME: check the last point */
274
275        for(i = 0; i < npts; i++)
276            debug("    2D pt[%i] (%g,%g)\n", i, pts[i][X], pts[i][Y]);
277
278        /* Compute the barycentre coordinates */
279        ctx = 0.;
280        cty = 0.;
281        weight = 0.;
282        for(i = 2; i < npts; i++)
283        {
284            double abx = pts[i - 1][X] - pts[0][X];
285            double aby = pts[i - 1][Y] - pts[0][Y];
286            double acx = pts[i][X] - pts[0][X];
287            double acy = pts[i][Y] - pts[0][Y];
288                        double area;
289            double sqarea = (abx * abx + aby * aby) * (acx * acx + acy * acy)
290                          - (abx * acx + aby * acy) * (abx * acx + aby * acy);
291            if(sqarea <= 0.)
292                continue;
293
294            area = sqrt(sqarea);
295            ctx += area * (abx + acx) / 3;
296            cty += area * (aby + acy) / 3;
297            weight += area;
298        }
299
300        if(weight > EPSILON)
301        {
302            ctx = pts[0][X] + ctx / weight;
303            cty = pts[0][Y] + cty / weight;
304        }
305        else
306        {
307            int right = -1;
308            tmp = -10.;
309            for(i = 0; i < npts; i++)
310                if(pts[i][X] > tmp)
311                {
312                    right = i;
313                    tmp = pts[i][X];
314                }
315            ctx = 0.5 * (pts[0][X] + pts[right][X]);
316            cty = 0.5 * (pts[0][Y] + pts[right][Y]);
317        }
318
319        debug("    2D bary (%g,%g)\n", ctx, cty);
320
321        /*
322         * 3.3. Store the barycentre and convex hull information.
323         */
324
325        ret->bary[n][R] = p0[R] + ctx * u[R] + cty * v[R];
326        ret->bary[n][G] = p0[G] + ctx * u[G] + cty * v[G];
327        ret->bary[n][B] = p0[B] + ctx * u[B] + cty * v[B];
328
329        for(i = 0; i < npts; i++)
330        {
331            ret->pts[n][i][R] = pts[i][R];
332            ret->pts[n][i][G] = pts[i][G];
333            ret->pts[n][i][B] = pts[i][B];
334            ret->pts[n][i][X] = pts[i][X] - ctx;
335            ret->pts[n][i][Y] = pts[i][Y] - cty;
336            ret->pts[n][i][A] = atan2(pts[i][Y] - cty, pts[i][X] - ctx);
337
338            debug("  3D pt[%i] (%g,%g,%g) angle %g\n",
339                  i, pts[i][R], pts[i][G], pts[i][B], ret->pts[n][i][A]);
340        }
341        ret->hullsize[n] = npts;
342
343        debug("  3D bary (%g,%g,%g)\n",
344              ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
345    }
346
347    return ret;
348}
349
350
351pipi_image_t *pipi_reduce(pipi_image_t *src,
352                          int ncolors, double const *palette)
353{
354    static double const rgbpal[] =
355    {
356        0, 0, 0,  0, 0, 1,  0, 1, 0,  0, 1, 1,
357        1, 0, 0,  1, 0, 1,  1, 1, 0,  1, 1, 1,
358    };
359
360    pipi_image_t *dst;
361        pipi_pixels_t *srcp, *dstp;
362        float *srcdata, *dstdata;
363        hull_t *rgbhull, *myhull;
364    int i, j, w, h;
365
366    init_uv();
367
368    rgbhull = compute_hull(8, rgbpal);
369    myhull = compute_hull(ncolors, palette);
370
371    /*
372     * 4. Load image and change its palette.
373     */
374
375    debug("\n### PROCESSING IMAGE ###\n\n");
376
377    srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
378    srcdata = (float *)srcp->pixels;
379
380    w = srcp->w;
381    h = srcp->h;
382
383    dst = pipi_new(w, h);
384    dstp = pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
385    dstdata = (float *)dstp->pixels;
386
387    for(j = 0; j < h; j++)
388        for(i = 0; i < w; i++)
389        {
390            double p[3];
391            double xp, yp, angle, xa, ya, xb, yb, t, s;
392            int slice, n, count;
393
394            /* FIXME: Imlib fucks up the RGB order. */
395            p[B] = srcdata[4 * (j * w + i)];
396            p[G] = srcdata[4 * (j * w + i) + 1];
397            p[R] = srcdata[4 * (j * w + i) + 2];
398
399            debug("Pixel +%i+%i (%g,%g,%g)\n", i, j, p[R], p[G], p[B]);
400
401            slice = (int)(BRIGHT(p) * STEPS + 0.5);
402
403            debug("  slice %i\n", slice);
404
405            /* Convert to 2D. The origin is the slice's barycentre. */
406            xp = (p[R] - rgbhull->bary[slice][R]) * u[R]
407               + (p[G] - rgbhull->bary[slice][G]) * u[G]
408               + (p[B] - rgbhull->bary[slice][B]) * u[B];
409            yp = (p[R] - rgbhull->bary[slice][R]) * v[R]
410               + (p[G] - rgbhull->bary[slice][G]) * v[G]
411               + (p[B] - rgbhull->bary[slice][B]) * v[B];
412
413            debug("    2D pt (%g,%g)\n", xp, yp);
414
415            /* 1. find the excentricity in RGB space. There is an easier
416             * way to do this, which is to find the intersection of our
417             * line with the RGB cube itself, but we'd lose the possibility
418             * of having an original colour space other than RGB. */
419
420            /* First, find the relevant triangle. */
421            count = rgbhull->hullsize[slice];
422            angle = atan2(yp, xp);
423            for(n = 0; n < count; n++)
424            {
425                double a1 = rgbhull->pts[slice][n][A];
426                double a2 = rgbhull->pts[slice][(n + 1) % count][A];
427                if(a1 > a2)
428                {
429                    if(angle >= a1)
430                        a2 += 2 * M_PI;
431                    else
432                        a1 -= 2 * M_PI;
433                }
434                if(angle >= a1 && angle <= a2)
435                    break;
436            }
437
438            /* Now compute the distance to the triangle's edge. If the edge
439             * intersection is M, then t is such as P = t.M (can be zero) */
440            xa = rgbhull->pts[slice][n % count][X];
441            ya = rgbhull->pts[slice][n % count][Y];
442            xb = rgbhull->pts[slice][(n + 1) % count][X];
443            yb = rgbhull->pts[slice][(n + 1) % count][Y];
444            t = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
445
446            if(t > 1.0)
447                t = 1.0;
448
449            debug("    best RGB %g (%g,%g) (%g,%g)\n", t, xa, ya, xb, yb);
450
451            /* 2. apply the excentricity in reduced space. */
452
453            count = myhull->hullsize[slice];
454            for(n = 0; n < count; n++)
455            {
456                double a1 = myhull->pts[slice][n][A];
457                double a2 = myhull->pts[slice][(n + 1) % count][A];
458                if(a1 > a2)
459                {
460                    if(angle >= a1)
461                        a2 += 2 * M_PI;
462                    else
463                        a1 -= 2 * M_PI;
464                }
465                if(angle >= a1 && angle <= a2)
466                    break;
467            }
468
469            /* If the edge intersection is M', s is such as P = s.M'. We
470             * want P' = t.M' = t.P/s  */
471            xa = myhull->pts[slice][n % count][X];
472            ya = myhull->pts[slice][n % count][Y];
473            xb = myhull->pts[slice][(n + 1) % count][X];
474            yb = myhull->pts[slice][(n + 1) % count][Y];
475            s = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
476
477            debug("    best custom %g (%g,%g) (%g,%g)\n", s, xa, ya, xb, yb);
478
479            if(s > 0)
480            {
481                xp *= t / s;
482                yp *= t / s;
483            }
484
485            p[R] = myhull->bary[slice][R] + xp * u[R] + yp * v[R];
486            p[G] = myhull->bary[slice][G] + xp * u[G] + yp * v[G];
487            p[B] = myhull->bary[slice][B] + xp * u[B] + yp * v[B];
488
489            /* Clipping should not be necessary, but the above code
490             * is unfortunately not perfect. */
491            if(p[R] < 0.0) p[R] = 0.0; else if(p[R] > 1.0) p[R] = 1.0;
492            if(p[G] < 0.0) p[G] = 0.0; else if(p[G] > 1.0) p[G] = 1.0;
493            if(p[B] < 0.0) p[B] = 0.0; else if(p[B] > 1.0) p[B] = 1.0;
494
495            dstdata[4 * (j * w + i)] = p[B];
496            dstdata[4 * (j * w + i) + 1] = p[G];
497            dstdata[4 * (j * w + i) + 2] = p[R];
498        }
499
500    free(rgbhull);
501    free(myhull);
502
503    return dst;
504}
505
Note: See TracBrowser for help on using the repository browser.