* More color reduction work: we can now constraint an image's palette to

the convex hull of a given palette, except it currently first converts
    the image to grayscale.

git-svn-id: file:///srv/caca.zoy.org/var/lib/svn/libpipi/trunk@2729 92316355-f0b4-4df1-b90c-862c8a59935f

examples/img2rubik.c

@@ -8,21 +8,25 @@
 
 #include <pipi.h>
 
-#define Y(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
+#define X 3
+#define Y 4
+#define A 5
 
-#define STEPS 16
+#define BRIGHT(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
+
+#define STEPS 1024
 #define EPSILON (0.000001)
 
 int main(int argc, char *argv[])
 {
     double palette[][3] =
     {
-        { 1.0, 0.0, 0.0 }, /* red */
-        { 0.0, 1.0, 0.0 }, /* green */
-        { 0.0, 0.0, 1.0 }, /* blue */
+        { 0.8, 0.0, 0.0 }, /* red */
+        { 0.0, 0.8, 0.0 }, /* green */
+        { 0.0, 0.0, 0.7 }, /* blue */
         { 1.0, 1.0, 1.0 }, /* white */
-        { 1.0, 1.0, 0.0 }, /* yellow */
-        { 1.0, 0.5, 0.0 }, /* orange */
+        { 0.9, 0.9, 0.0 }, /* yellow */
+        { 0.8, 0.4, 0.0 }, /* orange */
     };
 #define NCOLORS ((int)(sizeof(palette)/sizeof(*palette)))
 
@@ -36,7 +40,7 @@ int main(int argc, char *argv[])
     double min = 1.0, max = 0.0;
     for(i = 0; i < NCOLORS; i++)
     {
-        double p = Y(palette[i]);
+        double p = BRIGHT(palette[i]);
         if(p < min)
         {
             dark = palette[i];
@@ -48,8 +52,6 @@ int main(int argc, char *argv[])
             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
@@ -75,24 +77,25 @@ int main(int argc, char *argv[])
     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 bary[STEPS + 1][3];
+    int n;
+    for(n = 0; n <= STEPS; n++)
     {
-        double pts[NCOLORS*NCOLORS/2][3];
-        double p0[3];
+        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 * 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
@@ -135,16 +138,134 @@ int main(int argc, char *argv[])
                 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)
+         * 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] - dark[0]) * u[0]
+                      + (pts[i][1] - dark[1]) * u[1]
+                      + (pts[i][2] - dark[2]) * u[2];
+            pts[i][Y] = (pts[i][0] - dark[0]) * v[0]
+                      + (pts[i][1] - dark[1]) * v[1]
+                      + (pts[i][2] - dark[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]);
+        }
+
+        bary[n][0] = p0[0] + ctx * u[0] + cty * v[0];
+        bary[n][1] = p0[1] + ctx * u[1] + cty * v[1];
+        bary[n][2] = p0[2] + ctx * u[2] + cty * v[2];
     }
 
+
+
+    /*
+     * Last step: load image and change its palette.
+     */
+    pipi_image_t *src = pipi_load(argv[1]);
+    pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_Y_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++)
+        {
+            /* FIXME: Imlib fucks this up */
+            double p = srcdata[j * w + i];
+            dstdata[4 * (j * w + i)] = bary[(int)(p * STEPS)][2];
+            dstdata[4 * (j * w + i) + 1] = bary[(int)(p * STEPS)][1];
+            dstdata[4 * (j * w + i) + 2] = bary[(int)(p * STEPS)][0];
+        }
+
+    pipi_save(dst, argv[2]);
+
     return 0;
 }