Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.
 
 
 
 
 
 

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