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

Last change on this file since 2902 was 2902, checked in by Sam Hocevar, 11 years ago

Support C99 types on Win32 through the same hacks as in libcaca.

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