Changeset 2736 for libpipi


Ignore:
Timestamp:
Aug 20, 2008, 3:38:39 AM (14 years ago)
Author:
Sam Hocevar
Message:
  • Move the palette reduction algorithm into pipi_reduce().
Location:
libpipi/trunk
Files:
1 added
3 edited
1 copied

Legend:

Unmodified
Added
Removed
  • libpipi/trunk/examples/img2rubik.c

    r2734 r2736  
    99#include <pipi.h>
    1010
    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 
    27 typedef struct
    28 {
    29     double pts[STEPS + 1][MAXCOLORS * (MAXCOLORS - 1) / 2][6];
    30     int hullsize[STEPS + 1];
    31     double bary[STEPS + 1][3];
    32 }
    33 hull_t;
    34 
    35 static double const y[3] = { .299, .587, .114 };
    36 static double u[3], v[3];
    37 static int ylen;
    38 
    39 /*
    40  * Find two base vectors for the chrominance planes.
    41  */
    42 static 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  */
    64 static 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 
    32811int main(int argc, char *argv[])
    32912{
    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 
    33613    static double const mypal[] =
    33714    {
     
    34421    };
    34522
    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 
    35923    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 
     24    pipi_image_t *dst = pipi_reduce(src, 6, mypal);
    48225    pipi_save(dst, argv[2]);
    48326
  • libpipi/trunk/pipi/Makefile.am

    r2718 r2736  
    2828        $(combine_sources) \
    2929        $(filter_sources) \
     30        $(quantize_sources) \
    3031        $(dither_sources) \
    3132        $(NULL)
     
    5657        filter/color.c
    5758
     59quantize_sources = \
     60        quantize/reduce.c
     61
    5862dither_sources = \
    5963        dither/floydsteinberg.c \
  • libpipi/trunk/pipi/pipi.h

    r2725 r2736  
    132132                           int, int, float, float, float, float);
    133133
     134extern pipi_image_t *pipi_reduce(pipi_image_t *, int, double const *);
     135
    134136extern pipi_image_t *pipi_dither_floydsteinberg(pipi_image_t *, pipi_scan_t);
    135137extern pipi_image_t *pipi_dither_jajuni(pipi_image_t *, pipi_scan_t);
  • libpipi/trunk/pipi/quantize/reduce.c

    r2734 r2736  
     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
    119#include "config.h"
    220#include "common.h"
     
    2240
    2341#define MAXCOLORS 16
    24 #define STEPS 256
     42#define STEPS 1024
    2543#define EPSILON (0.000001)
    2644
     
    326344
    327345
    328 int main(int argc, char *argv[])
     346pipi_image_t *pipi_reduce(pipi_image_t *src,
     347                          int ncolors, double const *palette)
    329348{
    330349    static double const rgbpal[] =
     
    334353    };
    335354
    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 
    346355    int i, j;
    347356
     
    349358
    350359    hull_t *rgbhull = compute_hull(8, rgbpal);
    351     hull_t *myhull = compute_hull(6, mypal);
     360    hull_t *myhull = compute_hull(ncolors, palette);
    352361
    353362    /*
     
    357366    debug("\n### PROCESSING IMAGE ###\n\n");
    358367
    359     pipi_image_t *src = pipi_load(argv[1]);
    360368    pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
    361369    float *srcdata = (float *)srcp->pixels;
     
    480488    free(myhull);
    481489
    482     pipi_save(dst, argv[2]);
    483 
    484     return 0;
     490    return dst;
    485491}
    486492
Note: See TracChangeset for help on using the changeset viewer.