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