diff --git a/examples/img2rubik.c b/examples/img2rubik.c
index 815cbeb..2b0742e 100644
--- a/examples/img2rubik.c
+++ b/examples/img2rubik.c
@@ -4,26 +4,146 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 
 #include <pipi.h>
 
+#define Y(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
+
+#define STEPS 16
+#define EPSILON (0.000001)
+
 int main(int argc, char *argv[])
 {
-    char *srcname = NULL, *dstname = NULL;
-    pipi_image_t *img;
+    double palette[][3] =
+    {
+        { 1.0, 0.0, 0.0 }, /* red */
+        { 0.0, 1.0, 0.0 }, /* green */
+        { 0.0, 0.0, 1.0 }, /* blue */
+        { 1.0, 1.0, 1.0 }, /* white */
+        { 1.0, 1.0, 0.0 }, /* yellow */
+        { 1.0, 0.5, 0.0 }, /* orange */
+    };
+#define NCOLORS ((int)(sizeof(palette)/sizeof(*palette)))
 
-    if(argc < 3)
+    double tmp;
+    int i, j;
+
+    /*
+     * 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++)
     {
-        fprintf(stderr, "%s: too few arguments\n", argv[0]);
-        return EXIT_FAILURE;
+        double p = Y(palette[i]);
+        if(p < min)
+        {
+            dark = palette[i];
+            min = p;
+        }
+        if(p > max)
+        {
+            light = palette[i];
+            max = p;
+        }
     }
+    printf("lightest (%g,%g,%g)\n", light[0], light[1], light[2]);
+    printf("darkest (%g,%g,%g)\n", dark[0], dark[1], dark[2]);
+
+    /*
+     * 2. Find two base vectors for the chrominance planes
+     * FIXME: this doesn't work in all cases because u can be null
+     */
+    double y[3], u[3], v[3];
+    double ylen;
+
+    y[0] = light[0] - dark[0];
+    y[1] = light[1] - dark[1];
+    y[2] = light[2] - dark[2];
+    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;
+
+    printf("gray axis (%g,%g,%g) - length %g\n", y[0], y[1], y[2], ylen);
+    printf("u (%g,%g,%g)\n", u[0], u[1], u[2]);
+    printf("v (%g,%g,%g)\n", v[0], v[1], v[2]);
+
+    /*
+     * 3. Browse the grey axis and do stuff
+     */
+    double t;
+    for(t = 0.; t <= 1.0; t += 1. / STEPS)
+    {
+        double pts[NCOLORS*NCOLORS/2][3];
+        double p0[3];
+        int npts = 0;
+
+        p0[0] = dark[0] + t * y[0];
+        p0[1] = dark[1] + t * y[1];
+        p0[2] = dark[2] + t * y[2];
+        printf("%g,%g,%g\n", p0[0], p0[1], p0[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] = palette[i][0] - dark[0];
+            k1[1] = palette[i][1] - dark[1];
+            k1[2] = palette[i][2] - dark[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;
 
-    srcname = argv[1];
-    dstname = argv[2];
+                double k2[3];
+                k2[0] = palette[j][0] - dark[0];
+                k2[1] = palette[j][1] - dark[1];
+                k2[2] = palette[j][2] - dark[2];
+                tmp = sqrt(k2[0] * k2[0] + k2[1] * k2[1] + k2[2] * k2[2]);
 
-    img = pipi_load(srcname);
-    pipi_save(img, dstname);
-    pipi_free(img);
+                /* 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] = dark[0] + k1[0] + s * (k2[0] - k1[0]);
+                pts[npts][1] = dark[1] + k1[1] + s * (k2[1] - k1[1]);
+                pts[npts][2] = dark[2] + k1[2] + s * (k2[2] - k1[2]);
+                printf("  %i,%i: s = (%g,%g,%g)\n", i, j, pts[npts][0], pts[npts][1], pts[npts][2]);
+                npts++;
+            }
+        }
+
+        /*
+         * 3.2. Find the barycentre of these points' convex hull. (TODO)
+         */
+    }
 
     return 0;
 }