Nevar pievienot vairāk kā 25 tēmas Tēmai ir jāsākas ar burtu vai ciparu, tā var saturēt domu zīmes ('-') un var būt līdz 35 simboliem gara.
 
 
 
 
 
 

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