25개 이상의 토픽을 선택하실 수 없습니다. Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

493 lines
15 KiB

  1. /*
  2. * libpipi Pathetic image processing interface library
  3. * Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
  4. * All Rights Reserved
  5. *
  6. * $Id$
  7. *
  8. * This library is free software. It comes without any warranty, to
  9. * the extent permitted by applicable law. You can redistribute it
  10. * and/or modify it under the terms of the Do What The Fuck You Want
  11. * To Public License, Version 2, as published by Sam Hocevar. See
  12. * http://sam.zoy.org/wtfpl/COPYING for more details.
  13. */
  14. /*
  15. * reduce.c: palette reduction routines
  16. */
  17. #include "config.h"
  18. #include "common.h"
  19. #include <stdio.h>
  20. #include <stdlib.h>
  21. #include <string.h>
  22. #include <math.h>
  23. #include <pipi.h>
  24. #define R 0
  25. #define G 1
  26. #define B 2
  27. #define X 3
  28. #define Y 4
  29. #define A 5
  30. //#define debug printf
  31. #define debug(...) /* */
  32. #define BRIGHT(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])
  33. #define MAXCOLORS 16
  34. #define STEPS 1024
  35. #define EPSILON (0.000001)
  36. typedef struct
  37. {
  38. double pts[STEPS + 1][MAXCOLORS * (MAXCOLORS - 1) / 2][6];
  39. int hullsize[STEPS + 1];
  40. double bary[STEPS + 1][3];
  41. }
  42. hull_t;
  43. static double const y[3] = { .299, .587, .114 };
  44. static double u[3], v[3];
  45. static int ylen;
  46. /*
  47. * Find two base vectors for the chrominance planes.
  48. */
  49. static void init_uv(void)
  50. {
  51. double tmp;
  52. ylen = sqrt(y[R] * y[R] + y[G] * y[G] + y[B] * y[B]);
  53. u[R] = y[1];
  54. u[G] = -y[0];
  55. u[B] = 0;
  56. tmp = sqrt(u[R] * u[R] + u[G] * u[G] + u[B] * u[B]);
  57. u[R] /= tmp; u[G] /= tmp; u[B] /= tmp;
  58. v[R] = y[G] * u[B] - y[B] * u[G];
  59. v[G] = y[B] * u[R] - y[R] * u[B];
  60. v[B] = y[R] * u[G] - y[G] * u[R];
  61. tmp = sqrt(v[R] * v[R] + v[G] * v[G] + v[B] * v[B]);
  62. v[R] /= tmp; v[G] /= tmp; v[B] /= tmp;
  63. }
  64. /*
  65. * Compute the convex hull of a given palette.
  66. */
  67. static hull_t *compute_hull(int ncolors, double const *palette)
  68. {
  69. hull_t *ret = malloc(sizeof(hull_t));
  70. double tmp;
  71. int i, j;
  72. debug("\n### NEW HULL ###\n\n");
  73. debug("Analysing %i colors\n", ncolors);
  74. double pal[ncolors][3];
  75. for(i = 0; i < ncolors; i++)
  76. {
  77. pal[i][R] = palette[i * 3];
  78. pal[i][G] = palette[i * 3 + 1];
  79. pal[i][B] = palette[i * 3 + 2];
  80. debug(" [%i] (%g,%g,%g)\n", i, pal[i][R], pal[i][G], pal[i][B]);
  81. }
  82. /*
  83. * 1. Find the darkest and lightest colours
  84. */
  85. double *dark = NULL, *light = NULL;
  86. double min = 1.0, max = 0.0;
  87. for(i = 0; i < ncolors; i++)
  88. {
  89. double p = BRIGHT(pal[i]);
  90. if(p < min)
  91. {
  92. dark = pal[i];
  93. min = p;
  94. }
  95. if(p > max)
  96. {
  97. light = pal[i];
  98. max = p;
  99. }
  100. }
  101. double gray[3];
  102. gray[R] = light[R] - dark[R];
  103. gray[G] = light[G] - dark[G];
  104. gray[B] = light[B] - dark[B];
  105. debug(" gray axis (%g,%g,%g) - (%g,%g,%g)\n",
  106. dark[R], dark[G], dark[B], light[R], light[G], light[B]);
  107. /*
  108. * 3. Browse the grey axis and do stuff
  109. */
  110. int n;
  111. for(n = 0; n <= STEPS; n++)
  112. {
  113. double pts[ncolors * (ncolors - 1) / 2][5];
  114. double ptmp[5];
  115. #define SWAP(p1,p2) do { memcpy(ptmp, p1, sizeof(ptmp)); \
  116. memcpy(p1, p2, sizeof(ptmp)); \
  117. memcpy(p2, ptmp, sizeof(ptmp)); } while(0)
  118. double t = n * 1.0 / STEPS;
  119. int npts = 0;
  120. debug("Slice %i/%i\n", n, STEPS);
  121. double p0[3];
  122. p0[R] = dark[R] + t * gray[R];
  123. p0[G] = dark[G] + t * gray[G];
  124. p0[B] = dark[B] + t * gray[B];
  125. debug(" 3D gray (%g,%g,%g)\n", p0[R], p0[G], p0[B]);
  126. /*
  127. * 3.1. Find all edges that intersect the t.y + (u,v) plane
  128. */
  129. for(i = 0; i < ncolors; i++)
  130. {
  131. double k1[3];
  132. k1[R] = pal[i][R] - p0[R];
  133. k1[G] = pal[i][G] - p0[G];
  134. k1[B] = pal[i][B] - p0[B];
  135. tmp = sqrt(k1[R] * k1[R] + k1[G] * k1[G] + k1[B] * k1[B]);
  136. /* If k1.y > t.y.y, we don't want this point */
  137. double yk1 = y[R] * k1[R] + y[G] * k1[G] + y[B] * k1[B];
  138. if(yk1 > t * ylen * ylen + EPSILON)
  139. continue;
  140. for(j = 0; j < ncolors; j++)
  141. {
  142. if(i == j)
  143. continue;
  144. double k2[3];
  145. k2[R] = pal[j][R] - p0[R];
  146. k2[G] = pal[j][G] - p0[G];
  147. k2[B] = pal[j][B] - p0[B];
  148. tmp = sqrt(k2[R] * k2[R] + k2[G] * k2[G] + k2[B] * k2[B]);
  149. /* If k2.y < t.y.y, we don't want this point */
  150. double yk2 = y[R] * k2[R] + y[G] * k2[G] + y[B] * k2[B];
  151. if(yk2 < t * ylen * ylen - EPSILON)
  152. continue;
  153. if(yk2 < yk1)
  154. continue;
  155. double s = yk1 == yk2 ?
  156. 0.5 : (t * ylen * ylen - yk1) / (yk2 - yk1);
  157. pts[npts][R] = p0[R] + k1[R] + s * (k2[R] - k1[R]);
  158. pts[npts][G] = p0[G] + k1[G] + s * (k2[G] - k1[G]);
  159. pts[npts][B] = p0[B] + k1[B] + s * (k2[B] - k1[B]);
  160. npts++;
  161. }
  162. }
  163. /*
  164. * 3.2. Find the barycentre of these points' convex hull. We use
  165. * the Graham Scan technique.
  166. */
  167. /* Make our problem a 2-D problem. */
  168. for(i = 0; i < npts; i++)
  169. {
  170. pts[i][X] = (pts[i][R] - p0[R]) * u[R]
  171. + (pts[i][G] - p0[G]) * u[G]
  172. + (pts[i][B] - p0[B]) * u[B];
  173. pts[i][Y] = (pts[i][R] - p0[R]) * v[R]
  174. + (pts[i][G] - p0[G]) * v[G]
  175. + (pts[i][B] - p0[B]) * v[B];
  176. }
  177. /* Find the leftmost point */
  178. int left = -1;
  179. tmp = 10.;
  180. for(i = 0; i < npts; i++)
  181. if(pts[i][X] < tmp)
  182. {
  183. left = i;
  184. tmp = pts[i][X];
  185. }
  186. SWAP(pts[0], pts[left]);
  187. /* Sort the remaining points radially around pts[0]. Bubble sort
  188. * is okay for small sizes, I don't care. */
  189. for(i = 1; i < npts; i++)
  190. for(j = 1; j < npts - i; j++)
  191. {
  192. double k1 = (pts[j][X] - pts[0][X])
  193. * (pts[j + 1][Y] - pts[0][Y]);
  194. double k2 = (pts[j + 1][X] - pts[0][X])
  195. * (pts[j][Y] - pts[0][Y]);
  196. if(k1 < k2 - EPSILON)
  197. SWAP(pts[j], pts[j + 1]);
  198. else if(k1 < k2 + EPSILON)
  199. {
  200. /* Aligned! keep the farthest point */
  201. double ax = pts[j][X] - pts[0][X];
  202. double ay = pts[j][Y] - pts[0][Y];
  203. double bx = pts[j + 1][X] - pts[0][X];
  204. double by = pts[j + 1][Y] - pts[0][Y];
  205. if(ax * ax + ay * ay > bx * bx + by * by)
  206. SWAP(pts[j], pts[j + 1]);
  207. }
  208. }
  209. /* Remove points not in the convex hull */
  210. for(i = 2; i < npts; /* */)
  211. {
  212. if(i < 2)
  213. {
  214. i++;
  215. continue;
  216. }
  217. double k1 = (pts[i - 1][X] - pts[i - 2][X])
  218. * (pts[i][Y] - pts[i - 2][Y]);
  219. double k2 = (pts[i][X] - pts[i - 2][X])
  220. * (pts[i - 1][Y] - pts[i - 2][Y]);
  221. if(k1 <= k2 + EPSILON)
  222. {
  223. for(j = i - 1; j < npts - 1; j++)
  224. SWAP(pts[j], pts[j + 1]);
  225. npts--;
  226. }
  227. else
  228. i++;
  229. }
  230. /* FIXME: check the last point */
  231. for(i = 0; i < npts; i++)
  232. debug(" 2D pt[%i] (%g,%g)\n", i, pts[i][X], pts[i][Y]);
  233. /* Compute the barycentre coordinates */
  234. double ctx = 0., cty = 0., weight = 0.;
  235. for(i = 2; i < npts; i++)
  236. {
  237. double abx = pts[i - 1][X] - pts[0][X];
  238. double aby = pts[i - 1][Y] - pts[0][Y];
  239. double acx = pts[i][X] - pts[0][X];
  240. double acy = pts[i][Y] - pts[0][Y];
  241. double sqarea = (abx * abx + aby * aby) * (acx * acx + acy * acy)
  242. - (abx * acx + aby * acy) * (abx * acx + aby * acy);
  243. if(sqarea <= 0.)
  244. continue;
  245. double area = sqrt(sqarea);
  246. ctx += area * (abx + acx) / 3;
  247. cty += area * (aby + acy) / 3;
  248. weight += area;
  249. }
  250. if(weight > EPSILON)
  251. {
  252. ctx = pts[0][X] + ctx / weight;
  253. cty = pts[0][Y] + cty / weight;
  254. }
  255. else
  256. {
  257. int right = -1;
  258. tmp = -10.;
  259. for(i = 0; i < npts; i++)
  260. if(pts[i][X] > tmp)
  261. {
  262. right = i;
  263. tmp = pts[i][X];
  264. }
  265. ctx = 0.5 * (pts[0][X] + pts[right][X]);
  266. cty = 0.5 * (pts[0][Y] + pts[right][Y]);
  267. }
  268. debug(" 2D bary (%g,%g)\n", ctx, cty);
  269. /*
  270. * 3.3. Store the barycentre and convex hull information.
  271. */
  272. ret->bary[n][R] = p0[R] + ctx * u[R] + cty * v[R];
  273. ret->bary[n][G] = p0[G] + ctx * u[G] + cty * v[G];
  274. ret->bary[n][B] = p0[B] + ctx * u[B] + cty * v[B];
  275. for(i = 0; i < npts; i++)
  276. {
  277. ret->pts[n][i][R] = pts[i][R];
  278. ret->pts[n][i][G] = pts[i][G];
  279. ret->pts[n][i][B] = pts[i][B];
  280. ret->pts[n][i][X] = pts[i][X] - ctx;
  281. ret->pts[n][i][Y] = pts[i][Y] - cty;
  282. ret->pts[n][i][A] = atan2(pts[i][Y] - cty, pts[i][X] - ctx);
  283. debug(" 3D pt[%i] (%g,%g,%g) angle %g\n",
  284. i, pts[i][R], pts[i][G], pts[i][B], ret->pts[n][i][A]);
  285. }
  286. ret->hullsize[n] = npts;
  287. debug(" 3D bary (%g,%g,%g)\n",
  288. ret->bary[n][R], ret->bary[n][G], ret->bary[n][B]);
  289. }
  290. return ret;
  291. }
  292. pipi_image_t *pipi_reduce(pipi_image_t *src,
  293. int ncolors, double const *palette)
  294. {
  295. static double const rgbpal[] =
  296. {
  297. 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1,
  298. 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1,
  299. };
  300. int i, j;
  301. init_uv();
  302. hull_t *rgbhull = compute_hull(8, rgbpal);
  303. hull_t *myhull = compute_hull(ncolors, palette);
  304. /*
  305. * 4. Load image and change its palette.
  306. */
  307. debug("\n### PROCESSING IMAGE ###\n\n");
  308. pipi_pixels_t *srcp = pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
  309. float *srcdata = (float *)srcp->pixels;
  310. int w = srcp->w, h = srcp->h;
  311. pipi_image_t *dst = pipi_new(w, h);
  312. pipi_pixels_t *dstp = pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
  313. float *dstdata = (float *)dstp->pixels;
  314. for(j = 0; j < h; j++)
  315. for(i = 0; i < w; i++)
  316. {
  317. double p[3];
  318. /* FIXME: Imlib fucks up the RGB order. */
  319. p[B] = srcdata[4 * (j * w + i)];
  320. p[G] = srcdata[4 * (j * w + i) + 1];
  321. p[R] = srcdata[4 * (j * w + i) + 2];
  322. debug("Pixel +%i+%i (%g,%g,%g)\n", i, j, p[R], p[G], p[B]);
  323. int slice = (int)(BRIGHT(p) * STEPS + 0.5);
  324. debug(" slice %i\n", slice);
  325. /* Convert to 2D. The origin is the slice's barycentre. */
  326. double xp = (p[R] - rgbhull->bary[slice][R]) * u[R]
  327. + (p[G] - rgbhull->bary[slice][G]) * u[G]
  328. + (p[B] - rgbhull->bary[slice][B]) * u[B];
  329. double yp = (p[R] - rgbhull->bary[slice][R]) * v[R]
  330. + (p[G] - rgbhull->bary[slice][G]) * v[G]
  331. + (p[B] - rgbhull->bary[slice][B]) * v[B];
  332. debug(" 2D pt (%g,%g)\n", xp, yp);
  333. /* 1. find the excentricity in RGB space. There is an easier
  334. * way to do this, which is to find the intersection of our
  335. * line with the RGB cube itself, but we'd lose the possibility
  336. * of having an original colour space other than RGB. */
  337. /* First, find the relevant triangle. */
  338. int n, count = rgbhull->hullsize[slice];
  339. double angle = atan2(yp, xp);
  340. for(n = 0; n < count; n++)
  341. {
  342. double a1 = rgbhull->pts[slice][n][A];
  343. double a2 = rgbhull->pts[slice][(n + 1) % count][A];
  344. if(a1 > a2)
  345. {
  346. if(angle >= a1)
  347. a2 += 2 * M_PI;
  348. else
  349. a1 -= 2 * M_PI;
  350. }
  351. if(angle >= a1 && angle <= a2)
  352. break;
  353. }
  354. /* Now compute the distance to the triangle's edge. If the edge
  355. * intersection is M, then t is such as P = t.M (can be zero) */
  356. double xa = rgbhull->pts[slice][n % count][X];
  357. double ya = rgbhull->pts[slice][n % count][Y];
  358. double xb = rgbhull->pts[slice][(n + 1) % count][X];
  359. double yb = rgbhull->pts[slice][(n + 1) % count][Y];
  360. double t = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
  361. if(t > 1.0)
  362. t = 1.0;
  363. debug(" best RGB %g (%g,%g) (%g,%g)\n", t, xa, ya, xb, yb);
  364. /* 2. apply the excentricity in reduced space. */
  365. count = myhull->hullsize[slice];
  366. for(n = 0; n < count; n++)
  367. {
  368. double a1 = myhull->pts[slice][n][A];
  369. double a2 = myhull->pts[slice][(n + 1) % count][A];
  370. if(a1 > a2)
  371. {
  372. if(angle >= a1)
  373. a2 += 2 * M_PI;
  374. else
  375. a1 -= 2 * M_PI;
  376. }
  377. if(angle >= a1 && angle <= a2)
  378. break;
  379. }
  380. /* If the edge intersection is M', s is such as P = s.M'. We
  381. * want P' = t.M' = t.P/s */
  382. xa = myhull->pts[slice][n % count][X];
  383. ya = myhull->pts[slice][n % count][Y];
  384. xb = myhull->pts[slice][(n + 1) % count][X];
  385. yb = myhull->pts[slice][(n + 1) % count][Y];
  386. double s = (xp * (yb - ya) - yp * (xb - xa)) / (xa * yb - xb * ya);
  387. debug(" best custom %g (%g,%g) (%g,%g)\n", s, xa, ya, xb, yb);
  388. if(s > 0)
  389. {
  390. xp *= t / s;
  391. yp *= t / s;
  392. }
  393. p[R] = myhull->bary[slice][R] + xp * u[R] + yp * v[R];
  394. p[G] = myhull->bary[slice][G] + xp * u[G] + yp * v[G];
  395. p[B] = myhull->bary[slice][B] + xp * u[B] + yp * v[B];
  396. /* Clipping should not be necessary, but the above code
  397. * is unfortunately not perfect. */
  398. if(p[R] < 0.0) p[R] = 0.0; else if(p[R] > 1.0) p[R] = 1.0;
  399. if(p[G] < 0.0) p[G] = 0.0; else if(p[G] > 1.0) p[G] = 1.0;
  400. if(p[B] < 0.0) p[B] = 0.0; else if(p[B] > 1.0) p[B] = 1.0;
  401. dstdata[4 * (j * w + i)] = p[B];
  402. dstdata[4 * (j * w + i) + 1] = p[G];
  403. dstdata[4 * (j * w + i) + 2] = p[R];
  404. }
  405. free(rgbhull);
  406. free(myhull);
  407. return dst;
  408. }