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