libpipi/examples/img2rubik.c

#include "config.h"
#include "common.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <pipi.h>

#define X 3
#define Y 4
#define A 5

#define BRIGHT(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])

#define MAXCOLORS 16
#define STEPS 256
#define EPSILON (0.000001)

typedef struct
{
    double pts[STEPS + 1][MAXCOLORS * (MAXCOLORS - 1) / 2][6];
    int hullsize[STEPS + 1];
    double bary[STEPS + 1][3];
}
hull_t;

static double const y[3] = { .299, .587, .114 };
static double u[3], v[3];
static int ylen;

/*
 * Find two base vectors for the chrominance planes.
 */
static void init_uv(void)
{
    double tmp;

    ylen = sqrt(y[0] * y[0] + y[1] * y[1] + y[2] * y[2]);

    u[0] = y[1];
    u[1] = -y[0];
    u[2] = 0;
    tmp = sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
    u[0] /= tmp; u[1] /= tmp; u[2] /= tmp;

    v[0] = y[1] * u[2] - y[2] * u[1];
    v[1] = y[2] * u[0] - y[0] * u[2];
    v[2] = y[0] * u[1] - y[1] * u[0];
    tmp = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
    v[0] /= tmp; v[1] /= tmp; v[2] /= tmp;
}

/*
 * Compute the convex hull of a given palette.
 */
static hull_t *compute_hull(int ncolors, double const *palette)
{
    hull_t *ret = malloc(sizeof(hull_t));
    double tmp;
    int i, j;

    double pal[ncolors][3];
    for(i = 0; i < ncolors; i++)
    {
        pal[i][0] = palette[i * 3];
        pal[i][1] = palette[i * 3 + 1];
        pal[i][2] = palette[i * 3 + 2];
    }

    /*
     * 1. Find the darkest and lightest colours
     */
    double *dark = NULL, *light = NULL;
    double min = 1.0, max = 0.0;
    for(i = 0; i < ncolors; i++)
    {
        double p = BRIGHT(pal[i]);
        if(p < min)
        {
            dark = pal[i];
            min = p;
        }
        if(p > max)
        {
            light = pal[i];
            max = p;
        }
    }

    double gray[3];

    gray[0] = light[0] - dark[0];
    gray[1] = light[1] - dark[1];
    gray[2] = light[2] - dark[2];

    /*
     * 3. Browse the grey axis and do stuff
     */
    int n;
    for(n = 0; n <= STEPS; n++)
    {
        double pts[ncolors * (ncolors - 1) / 2][6];
        double ptmp[6];
#define SWAP(p1,p2) do { memcpy(ptmp, p1, sizeof(ptmp)); \
                         memcpy(p1, p2, sizeof(ptmp)); \
                         memcpy(p2, ptmp, sizeof(ptmp)); } while(0)
        double t = n * 1.0 / STEPS;
        int npts = 0;

        double p0[3];
        p0[0] = dark[0] + t * gray[0];
        p0[1] = dark[1] + t * gray[1];
        p0[2] = dark[2] + t * gray[2];

        /*
         * 3.1. Find all edges that intersect the t.y + (u,v) plane
         */
        for(i = 0; i < ncolors; i++)
        {
            double k1[3];
            k1[0] = pal[i][0] - p0[0];
            k1[1] = pal[i][1] - p0[1];
            k1[2] = pal[i][2] - p0[2];
            tmp = sqrt(k1[0] * k1[0] + k1[1] * k1[1] + k1[2] * k1[2]);

            /* If k1.y > t.y.y, we don't want this point */
            double yk1 = y[0] * k1[0] + y[1] * k1[1] + y[2] * k1[2];
            if(yk1 > t * ylen * ylen + EPSILON)
                continue;

            for(j = 0; j < ncolors; j++)
            {
                if(i == j)
                    continue;

                double k2[3];
                k2[0] = pal[j][0] - p0[0];
                k2[1] = pal[j][1] - p0[1];
                k2[2] = pal[j][2] - p0[2];
                tmp = sqrt(k2[0] * k2[0] + k2[1] * k2[1] + k2[2] * k2[2]);

                /* If k2.y < t.y.y, we don't want this point */
                double yk2 = y[0] * k2[0] + y[1] * k2[1] + y[2] * k2[2];
                if(yk2 < t * ylen * ylen - EPSILON)
                    continue;

                if(yk2 < yk1)
                    continue;

                double s = yk1 == yk2 ?
                           0.5 : (t * ylen * ylen - yk1) / (yk2 - yk1);

                pts[npts][0] = p0[0] + k1[0] + s * (k2[0] - k1[0]);
                pts[npts][1] = p0[1] + k1[1] + s * (k2[1] - k1[1]);
                pts[npts][2] = p0[2] + k1[2] + s * (k2[2] - k1[2]);
                npts++;
            }
        }

        /*
         * 3.2. Find the barycentre of these points' convex hull. We use
         *      the Graham Scan technique.
         */

        /* Make our problem a 2-D problem. */
        for(i = 0; i < npts; i++)
        {
            pts[i][X] = (pts[i][0] - p0[0]) * u[0]
                           + (pts[i][1] - p0[1]) * u[1]
                           + (pts[i][2] - p0[2]) * u[2];
            pts[i][Y] = (pts[i][0] - p0[0]) * v[0]
                           + (pts[i][1] - p0[1]) * v[1]
                           + (pts[i][2] - p0[2]) * v[2];
        }

        /* Find the leftmost point */
        int left = -1;
        tmp = 10.;
        for(i = 0; i < npts; i++)
            if(pts[i][X] < tmp)
            {
                left = i;
                tmp = pts[i][X];
            }
        SWAP(pts[0], pts[left]);

        /* Sort the remaining points radially around pts[0] */
        for(i = 1; i < npts; i++)
            for(j = 1; j < npts - i; j++)
                if((pts[j][X] - pts[0][X]) * (pts[j + 1][Y] - pts[0][Y])
                 < (pts[j + 1][X] - pts[0][X]) * (pts[j][Y] - pts[0][Y]))
                    SWAP(pts[j], pts[j + 1]);

        /* Remove points not in the convex hull */
        for(i = 2; i < npts; /* */)
        {
            if(i < 2)
            {
                i++;
                continue;
            }

            if((pts[i - 1][X] - pts[i - 2][X]) * (pts[i][Y] - pts[i - 2][Y])
               < (pts[i][X] - pts[i - 2][X]) * (pts[i - 1][Y] - pts[i - 2][Y]))
            {
                for(j = i - 1; j < npts - 1; j++)
                    SWAP(pts[j], pts[j + 1]);
                npts--;
            }
            else
                i++;
        }

        /* Compute the barycentre coordinates */
        double ctx = 0., cty = 0., weight = 0.;
        for(i = 2; i < npts; i++)
        {
            double abx = pts[i - 1][X] - pts[0][X];
            double aby = pts[i - 1][Y] - pts[0][Y];
            double acx = pts[i][X] - pts[0][X];
            double acy = pts[i][Y] - pts[0][Y];
            double sqarea = (abx * abx + aby * aby) * (acx * acx + acy * acy)
                          - (abx * acx + aby * acy) * (abx * acx + aby * acy);
            if(sqarea <= 0.)
                continue;

            double area = sqrt(sqarea);
            ctx += area * (abx + acx) / 3;
            cty += area * (aby + acy) / 3;
            weight += area;
        }

        if(weight > EPSILON)
        {
            ctx = pts[0][X] + ctx / weight;
            cty = pts[0][Y] + cty / weight;
        }
        else
        {
            int right = -1;
            tmp = -10.;
            for(i = 0; i < npts; i++)
                if(pts[i][X] > tmp)
                {
                    right = i;
                    tmp = pts[i][X];
                }
            ctx = 0.5 * (pts[0][X] + pts[right][X]);
            cty = 0.5 * (pts[0][Y] + pts[right][Y]);
        }

        /*
         * 3.3. Store the barycentre and convex hull information.
         */

        ret->bary[n][0] = p0[0] + ctx * u[0] + cty * v[0];
        ret->bary[n][1] = p0[1] + ctx * u[1] + cty * v[1];
        ret->bary[n][2] = p0[2] + ctx * u[2] + cty * v[2];
//if(ncolors == 6) printf("%i %g %g %g\n", n, ret->bary[n][0], ret->bary[n][1], ret->bary[n][2]);

        for(i = 0; i < npts; i++)
        {
            ret->pts[n][i][0] = pts[i][0];
            ret->pts[n][i][1] = pts[i][1];
            ret->pts[n][i][2] = pts[i][2];
        }
        ret->hullsize[n] = npts;
    }

    return ret;
}


int main(int argc, char *argv[])
{
    static double const rgbpal[] =
    {
        0, 0, 0,  0, 0, 1,  0, 1, 0,  0, 1, 1,
        1, 0, 0,  1, 0, 1,  1, 1, 0,  1, 1, 1,
    };

    static double const mypal[] =
    {
        0.8, 0.0, 0.0, /* red */
        0.0, 0.8, 0.0, /* green */
        0.0, 0.0, 0.6, /* blue */
        1.0, 1.0, 1.0, /* white */
        0.9, 0.9, 0.0, /* yellow */
        0.8, 0.4, 0.0, /* orange */
    };

    int i, j;

    init_uv();

    hull_t *rgbhull = compute_hull(8, rgbpal);
    hull_t *myhull = compute_hull(6, mypal);

    /*
     * 4. Load image and change its palette.
     */

    pipi_image_t *src = pipi_load(argv[1]);
    pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
    float *srcdata = (float *)srcp->pixels;

    int w = srcp->w, h = srcp->h;

    pipi_image_t *dst = pipi_new(w, h);
    pipi_pixels_t *dstp = pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
    float *dstdata = (float *)dstp->pixels;

    for(j = 0; j < h; j++)
        for(i = 0; i < w; i++)
        {
            double p[3];
            /* FIXME: Imlib fucks up the RGB order. */
            p[2] = srcdata[4 * (j * w + i)];
            p[1] = srcdata[4 * (j * w + i) + 1];
            p[0] = srcdata[4 * (j * w + i) + 2];
            double gray = BRIGHT(p);

            double dx = (p[0] - gray * y[0]) * u[0]
                      + (p[1] - gray * y[1]) * u[1]
                      + (p[2] - gray * y[2]) * u[2];
            double dy = (p[0] - gray * y[0]) * v[0]
                      + (p[1] - gray * y[1]) * v[1]
                      + (p[2] - gray * y[2]) * v[2];

            dstdata[4 * (j * w + i)] = myhull->bary[(int)(gray * STEPS)][2]
                                     + dx * u[2] * .3 + dy * v[2] * .3;
            dstdata[4 * (j * w + i) + 1] = myhull->bary[(int)(gray * STEPS)][1]
                                         + dx * u[1] * .3 + dy * v[1] * .3;
            dstdata[4 * (j * w + i) + 2] = myhull->bary[(int)(gray * STEPS)][0]
                                         + dx * u[0] * .3 + dy * v[0] * .3;
        }

    pipi_save(dst, argv[2]);

    return 0;
}