diff --git a/examples/img2rubik.c b/examples/img2rubik.c index ee2391b..de7a053 100644 --- a/examples/img2rubik.c +++ b/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;