* img2rubik.c:

+ fix RGB-space UV projection.
    + free allocated memory at the end.
    + clean up the debug messages and make them optional.

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

examples/img2rubik.c

@@ -15,10 +15,13 @@
 #define Y 4
 #define A 5
 
+//#define debug printf
+#define debug(...) /* */
+
 #define BRIGHT(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
 
 #define MAXCOLORS 16
-#define STEPS 32
+#define STEPS 25
 #define EPSILON (0.000001)
 
 typedef struct
@@ -64,12 +67,17 @@ static hull_t *compute_hull(int ncolors, double const *palette)
     double tmp;
     int i, j;
 
+    debug("\n### NEW HULL ###\n\n");
+
+    debug("Analysing %i colors\n", ncolors);
+
     double pal[ncolors][3];
     for(i = 0; i < ncolors; i++)
     {
         pal[i][R] = palette[i * 3];
         pal[i][G] = palette[i * 3 + 1];
         pal[i][B] = palette[i * 3 + 2];
+        debug("  [%i] (%g,%g,%g)\n", i, pal[i][R], pal[i][G], pal[i][B]);
     }
 
     /*
@@ -98,13 +106,15 @@ static hull_t *compute_hull(int ncolors, double const *palette)
     gray[G] = light[G] - dark[G];
     gray[B] = light[B] - dark[B];
 
+    debug("  gray axis (%g,%g,%g) - (%g,%g,%g)\n",
+          dark[R], dark[G], dark[B], light[R], light[G], light[B]);
+
     /*
      * 3. Browse the grey axis and do stuff
      */
     int n;
     for(n = 0; n <= STEPS; n++)
     {
-printf("slice %i\n", n);
         double pts[ncolors * (ncolors - 1) / 2][5];
         double ptmp[5];
 #define SWAP(p1,p2) do { memcpy(ptmp, p1, sizeof(ptmp)); \
@@ -113,11 +123,15 @@ printf("slice %i\n", n);
         double t = n * 1.0 / STEPS;
         int npts = 0;
 
+        debug("Slice %i/%i\n", n, STEPS);
+
         double p0[3];
         p0[R] = dark[R] + t * gray[R];
         p0[G] = dark[G] + t * gray[G];
         p0[B] = dark[B] + t * gray[B];
 
+        debug("  3D gray (%g,%g,%g)\n", p0[R], p0[G], p0[B]);
+
         /*
          * 3.1. Find all edges that intersect the t.y + (u,v) plane
          */
@@ -214,15 +228,11 @@ printf("slice %i\n", n);
                 }
             }
 
-for(i = 0; i < npts; i++)
-   printf("  : %g %g\n", pts[i][X], pts[i][Y]);
         /* Remove points not in the convex hull */
         for(i = 2; i < npts; /* */)
         {
-printf("  # %g %g", pts[i][X], pts[i][Y]);
             if(i < 2)
             {
-printf("  skipping\n");
                 i++;
                 continue;
             }
@@ -231,22 +241,18 @@ printf("  skipping\n");
                           * (pts[i][Y] - pts[i - 2][Y]);
             double k2 = (pts[i][X] - pts[i - 2][X])
                           * (pts[i - 1][Y] - pts[i - 2][Y]);
-printf("  %g < %g ? ", k1, k2);
             if(k1 <= k2 + EPSILON)
             {
-printf("yes, removing\n");
                 for(j = i - 1; j < npts - 1; j++)
                     SWAP(pts[j], pts[j + 1]);
                 npts--;
             }
             else
-{
                 i++;
-printf("no\n");
-}
         }
-for(i = 0; i < npts; i++)
-   printf("  : %g %g\n", pts[i][X], pts[i][Y]);
+
+        for(i = 0; i < npts; i++)
+            debug("    2D pt[%i] (%g,%g)\n", i, pts[i][X], pts[i][Y]);
 
         /* Compute the barycentre coordinates */
         double ctx = 0., cty = 0., weight = 0.;
@@ -286,6 +292,8 @@ for(i = 0; i < npts; i++)
             cty = 0.5 * (pts[0][Y] + pts[right][Y]);
         }
 
+        debug("    2D bary (%g,%g)\n", ctx, cty);
+
         /*
          * 3.3. Store the barycentre and convex hull information.
          */
@@ -293,7 +301,6 @@ for(i = 0; i < npts; i++)
         ret->bary[n][R] = p0[R] + ctx * u[R] + cty * v[R];
         ret->bary[n][G] = p0[G] + ctx * u[G] + cty * v[G];
         ret->bary[n][B] = p0[B] + ctx * u[B] + cty * v[B];
-//if(ncolors == 6) printf("%i %g %g %g\n", n, ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
 
         for(i = 0; i < npts; i++)
         {
@@ -303,10 +310,14 @@ for(i = 0; i < npts; i++)
             ret->pts[n][i][X] = pts[i][X] - ctx;
             ret->pts[n][i][Y] = pts[i][Y] - cty;
             ret->pts[n][i][A] = atan2(pts[i][Y] - cty, pts[i][X] - ctx);
-printf("%g %g -> %g\n", ret->pts[n][i][X], ret->pts[n][i][Y], ret->pts[n][i][A]);
+
+            debug("  3D pt[%i] (%g,%g,%g) angle %g\n",
+                  i, pts[i][R], pts[i][G], pts[i][B], ret->pts[n][i][A]);
         }
-printf("\n");
         ret->hullsize[n] = npts;
+
+        debug("  3D bary (%g,%g,%g)\n",
+              ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
     }
 
     return ret;
@@ -323,12 +334,12 @@ int main(int argc, char *argv[])
 
     static double const mypal[] =
     {
-        0.8, 0.001, 0.001, /* red */
-        0.001, 0.8, 0.001, /* green */
-        0.05, 0.001, 0.4, /* blue */
-        0.9, 0.9, 0.9, /* white */
-        0.9, 0.9, 0.001, /* yellow */
-        0.8, 0.4, 0.001, /* orange */
+        0.900, 0.001, 0.001, /* red */
+        0.001, 0.800, 0.001, /* green */
+        0.005, 0.001, 0.500, /* blue */
+        0.900, 0.900, 0.900, /* white */
+        0.900, 0.900, 0.001, /* yellow */
+        0.800, 0.400, 0.001, /* orange */
     };
 
     int i, j;
@@ -342,6 +353,8 @@ int main(int argc, char *argv[])
      * 4. Load image and change its palette.
      */
 
+    debug("\n### PROCESSING IMAGE ###\n\n");
+
     pipi_image_t *src = pipi_load(argv[1]);
     pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
     float *srcdata = (float *)srcp->pixels;
@@ -360,18 +373,22 @@ int main(int argc, char *argv[])
             p[B] = srcdata[4 * (j * w + i)];
             p[G] = srcdata[4 * (j * w + i) + 1];
             p[R] = srcdata[4 * (j * w + i) + 2];
-            double gray = BRIGHT(p);
 
-            /* Convert to 2-D */
-            double xp = (p[R] - gray) * u[R]
-                      + (p[G] - gray) * u[G]
-                      + (p[B] - gray) * u[B];
-            double yp = (p[R] - gray) * v[R]
-                      + (p[G] - gray) * v[G]
-                      + (p[B] - gray) * v[B];
-printf("PIX(%g,%g,%g): gray %g + P(%g,%g)\n", p[R], p[G], p[B], gray, xp, yp);
+            debug("Pixel +%i+%i (%g,%g,%g)\n", i, j, p[R], p[G], p[B]);
+
+            int slice = (int)(BRIGHT(p) * STEPS + 0.5);
 
-            int pt = (int)(gray * STEPS + 0.5);
+            debug("  slice %i\n", slice);
+
+            /* Convert to 2D. The origin is the slice's barycentre. */
+            double xp = (p[R] - rgbhull->bary[slice][R]) * u[R]
+                      + (p[G] - rgbhull->bary[slice][G]) * u[G]
+                      + (p[B] - rgbhull->bary[slice][B]) * u[B];
+            double yp = (p[R] - rgbhull->bary[slice][R]) * v[R]
+                      + (p[G] - rgbhull->bary[slice][G]) * v[G]
+                      + (p[B] - rgbhull->bary[slice][B]) * v[B];
+
+            debug("    2D pt (%g,%g)\n", xp, yp);
 
             /* 1. find the excentricity in RGB space. There is an easier
              * way to do this, which is to find the intersection of our
@@ -379,13 +396,12 @@ printf("PIX(%g,%g,%g): gray %g + P(%g,%g)\n", p[R], p[G], p[B], gray, xp, yp);
              * of having an original colour space other than RGB. */
 
             /* First, find the relevant triangle. */
-            int n, count = rgbhull->hullsize[pt];
+            int n, count = rgbhull->hullsize[slice];
             double angle = atan2(yp, xp);
-printf("  slice %i angle %g\n", pt, angle);
             for(n = 0; n < count; n++)
             {
-                double a1 = rgbhull->pts[pt][n][A];
-                double a2 = rgbhull->pts[pt][(n + 1) % count][A];
+                double a1 = rgbhull->pts[slice][n][A];
+                double a2 = rgbhull->pts[slice][(n + 1) % count][A];
                 if(a1 > a2)
                 {
                     if(angle >= a1)
@@ -399,20 +415,21 @@ printf("  slice %i angle %g\n", pt, angle);
 
             /* Now compute the distance to the triangle's edge. If the edge
              * intersection is M, then t is such as P = t.M (can be zero) */
-            double xa = rgbhull->pts[pt][n % count][X];
-            double ya = rgbhull->pts[pt][n % count][Y];
-            double xb = rgbhull->pts[pt][(n + 1) % count][X];
-            double yb = rgbhull->pts[pt][(n + 1) % count][Y];
+            double xa = rgbhull->pts[slice][n % count][X];
+            double ya = rgbhull->pts[slice][n % count][Y];
+            double xb = rgbhull->pts[slice][(n + 1) % count][X];
+            double yb = rgbhull->pts[slice][(n + 1) % count][Y];
             double t = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
-printf("  best is %g %g (%g) -- %g %g (%g)\n", xa, ya, rgbhull->pts[pt][n % count][A], xb, yb, rgbhull->pts[pt][(n + 1) % count][A]);
+
+            debug("    best RGB %g (%g,%g) (%g,%g)\n", t, xa, ya, xb, yb);
 
             /* 2. apply the excentricity in reduced space. */
 
-            count = myhull->hullsize[pt];
+            count = myhull->hullsize[slice];
             for(n = 0; n < count; n++)
             {
-                double a1 = myhull->pts[pt][n][A];
-                double a2 = myhull->pts[pt][(n + 1) % count][A];
+                double a1 = myhull->pts[slice][n][A];
+                double a2 = myhull->pts[slice][(n + 1) % count][A];
                 if(a1 > a2)
                 {
                     if(angle >= a1)
@@ -426,34 +443,38 @@ printf("  best is %g %g (%g) -- %g %g (%g)\n", xa, ya, rgbhull->pts[pt][n % coun
 
             /* If the edge intersection is M', s is such as P = s.M'. We
              * want P' = t.M' = t.P/s  */
-            xa = myhull->pts[pt][n % count][X];
-            ya = myhull->pts[pt][n % count][Y];
-            xb = myhull->pts[pt][(n + 1) % count][X];
-            yb = myhull->pts[pt][(n + 1) % count][Y];
+            xa = myhull->pts[slice][n % count][X];
+            ya = myhull->pts[slice][n % count][Y];
+            xb = myhull->pts[slice][(n + 1) % count][X];
+            yb = myhull->pts[slice][(n + 1) % count][Y];
             double s = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
-printf("  apply to %g %g (%g) -- %g %g (%g)\n", xa, ya, myhull->pts[pt][n % count][A], xb, yb, myhull->pts[pt][(n + 1) % count][A]);
+
+            debug("    best custom %g (%g,%g) (%g,%g)\n", s, xa, ya, xb, yb);
 
             if(s > 0)
             {
                 xp *= t / s;
                 yp *= t / s;
             }
-printf("    t = %g   s = %g  -> (%g,%g)\n", t, s, xp, yp);
-
-            dstdata[4 * (j * w + i)] = myhull->bary[(int)(gray * STEPS)][B]
-                                     + xp * u[B] + yp * v[B];
-            dstdata[4 * (j * w + i) + 1] = myhull->bary[(int)(gray * STEPS)][G]
-                                         + xp * u[G] + yp * v[G];
-            dstdata[4 * (j * w + i) + 2] = myhull->bary[(int)(gray * STEPS)][R]
-                                         + xp * u[R] + yp * v[R];
-if(dstdata[4 * (j * w + i)] < 0.) dstdata[4 * (j * w + i)] = 0.;
-if(dstdata[4 * (j * w + i)] > 1.) dstdata[4 * (j * w + i)] = 1.;
-if(dstdata[4 * (j * w + i) + 1] < 0.) dstdata[4 * (j * w + i) + 1] = 0.;
-if(dstdata[4 * (j * w + i) + 1] > 1.) dstdata[4 * (j * w + i) + 1] = 1.;
-if(dstdata[4 * (j * w + i) + 2] < 0.) dstdata[4 * (j * w + i) + 2] = 0.;
-if(dstdata[4 * (j * w + i) + 2] > 1.) dstdata[4 * (j * w + i) + 2] = 1.;
+
+            p[R] = myhull->bary[slice][R] + xp * u[R] + yp * v[R];
+            p[G] = myhull->bary[slice][G] + xp * u[G] + yp * v[G];
+            p[B] = myhull->bary[slice][B] + xp * u[B] + yp * v[B];
+
+            /* Clipping should not be necessary, but the above code
+             * is unfortunately not perfect. */
+            if(p[R] < 0.0) p[R] = 0.0; else if(p[R] > 1.0) p[R] = 1.0;
+            if(p[G] < 0.0) p[G] = 0.0; else if(p[G] > 1.0) p[G] = 1.0;
+            if(p[B] < 0.0) p[B] = 0.0; else if(p[B] > 1.0) p[B] = 1.0;
+
+            dstdata[4 * (j * w + i)] = p[B];
+            dstdata[4 * (j * w + i) + 1] = p[G];
+            dstdata[4 * (j * w + i) + 2] = p[R];
         }
 
+    free(rgbhull);
+    free(myhull);
+
     pipi_save(dst, argv[2]);
 
     return 0;