Changeset View
Changeset View
Standalone View
Standalone View
src/hsluv.c
- This file was added.
1 | /* | ||||
---|---|---|---|---|---|
2 | * HSLuv-C: Human-friendly HSL | ||||
3 | * <http://github.com/hsluv/hsluv-c> | ||||
4 | * <http://www.hsluv.org/> | ||||
5 | * | ||||
6 | * Copyright (c) 2015 Alexei Boronine (original idea, JavaScript implementation) | ||||
7 | * Copyright (c) 2015 Roger Tallada (Obj-C implementation) | ||||
8 | * Copyright (c) 2017 Martin Mitas (C implementation, based on Obj-C implementation) | ||||
9 | * | ||||
10 | * Permission is hereby granted, free of charge, to any person obtaining a | ||||
11 | * copy of this software and associated documentation files (the "Software"), | ||||
12 | * to deal in the Software without restriction, including without limitation | ||||
13 | * the rights to use, copy, modify, merge, publish, distribute, sublicense, | ||||
14 | * and/or sell copies of the Software, and to permit persons to whom the | ||||
15 | * Software is furnished to do so, subject to the following conditions: | ||||
16 | * | ||||
17 | * The above copyright notice and this permission notice shall be included in | ||||
18 | * all copies or substantial portions of the Software. | ||||
19 | * | ||||
20 | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | ||||
21 | * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | ||||
22 | * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | ||||
23 | * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | ||||
24 | * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | ||||
25 | * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS | ||||
26 | * IN THE SOFTWARE. | ||||
27 | */ | ||||
28 | | ||||
29 | #include "hsluv.h" | ||||
30 | | ||||
31 | #include <float.h> | ||||
32 | #include <math.h> | ||||
33 | | ||||
34 | | ||||
35 | typedef struct Triplet_tag Triplet; | ||||
36 | struct Triplet_tag { | ||||
37 | double a; | ||||
38 | double b; | ||||
39 | double c; | ||||
40 | }; | ||||
41 | | ||||
42 | /* for RGB */ | ||||
43 | static const Triplet m[3] = { | ||||
44 | { 3.24096994190452134377, -1.53738317757009345794, -0.49861076029300328366 }, | ||||
45 | { -0.96924363628087982613, 1.87596750150772066772, 0.04155505740717561247 }, | ||||
46 | { 0.05563007969699360846, -0.20397695888897656435, 1.05697151424287856072 } | ||||
47 | }; | ||||
48 | | ||||
49 | /* for XYZ */ | ||||
50 | static const Triplet m_inv[3] = { | ||||
51 | { 0.41239079926595948129, 0.35758433938387796373, 0.18048078840183428751 }, | ||||
52 | { 0.21263900587151035754, 0.71516867876775592746, 0.07219231536073371500 }, | ||||
53 | { 0.01933081871559185069, 0.11919477979462598791, 0.95053215224966058086 } | ||||
54 | }; | ||||
55 | | ||||
56 | static const double ref_u = 0.19783000664283680764; | ||||
57 | static const double ref_v = 0.46831999493879100370; | ||||
58 | | ||||
59 | static const double kappa = 903.29629629629629629630; | ||||
60 | static const double epsilon = 0.00885645167903563082; | ||||
61 | | ||||
62 | | ||||
63 | typedef struct Bounds_tag Bounds; | ||||
64 | struct Bounds_tag { | ||||
65 | double a; | ||||
66 | double b; | ||||
67 | }; | ||||
68 | | ||||
69 | | ||||
70 | static void | ||||
71 | get_bounds(double l, Bounds bounds[6]) | ||||
72 | { | ||||
73 | double tl = l + 16.0; | ||||
74 | double sub1 = (tl * tl * tl) / 1560896.0; | ||||
75 | double sub2 = (sub1 > epsilon ? sub1 : (l / kappa)); | ||||
76 | int channel; | ||||
77 | int t; | ||||
78 | | ||||
79 | for(channel = 0; channel < 3; channel++) { | ||||
80 | double m1 = m[channel].a; | ||||
81 | double m2 = m[channel].b; | ||||
82 | double m3 = m[channel].c; | ||||
83 | | ||||
84 | for (t = 0; t < 2; t++) { | ||||
85 | double top1 = (284517.0 * m1 - 94839.0 * m3) * sub2; | ||||
86 | double top2 = (838422.0 * m3 + 769860.0 * m2 + 731718.0 * m1) * l * sub2 - 769860.0 * t * l; | ||||
87 | double bottom = (632260.0 * m3 - 126452.0 * m2) * sub2 + 126452.0 * t; | ||||
88 | | ||||
89 | bounds[channel * 2 + t].a = top1 / bottom; | ||||
90 | bounds[channel * 2 + t].b = top2 / bottom; | ||||
91 | } | ||||
92 | } | ||||
93 | } | ||||
94 | | ||||
95 | static double | ||||
96 | intersect_line_line(const Bounds* line1, const Bounds* line2) | ||||
97 | { | ||||
98 | return (line1->b - line2->b) / (line2->a - line1->a); | ||||
99 | } | ||||
100 | | ||||
101 | static double | ||||
102 | dist_from_pole_squared(double x, double y) | ||||
103 | { | ||||
104 | return x * x + y * y; | ||||
105 | } | ||||
106 | | ||||
107 | static double | ||||
108 | ray_length_until_intersect(double theta, const Bounds* line) | ||||
109 | { | ||||
110 | return line->b / (sin(theta) - line->a * cos(theta)); | ||||
111 | } | ||||
112 | | ||||
113 | static double | ||||
114 | max_safe_chroma_for_l(double l) | ||||
115 | { | ||||
116 | double min_len_squared = DBL_MAX; | ||||
117 | Bounds bounds[6]; | ||||
118 | int i; | ||||
119 | | ||||
120 | get_bounds(l, bounds); | ||||
121 | for(i = 0; i < 6; i++) { | ||||
122 | double m1 = bounds[i].a; | ||||
123 | double b1 = bounds[i].b; | ||||
124 | /* x where line intersects with perpendicular running though (0, 0) */ | ||||
125 | Bounds line2 = { -1.0 / m1, 0.0 }; | ||||
126 | double x = intersect_line_line(&bounds[i], &line2); | ||||
127 | double distance = dist_from_pole_squared(x, b1 + x * m1); | ||||
128 | | ||||
129 | if(distance < min_len_squared) | ||||
130 | min_len_squared = distance; | ||||
131 | } | ||||
132 | | ||||
133 | return sqrt(min_len_squared); | ||||
134 | } | ||||
135 | | ||||
136 | static double | ||||
137 | max_chroma_for_lh(double l, double h) | ||||
138 | { | ||||
139 | double min_len = DBL_MAX; | ||||
140 | double hrad = h * 0.01745329251994329577; /* (2 * pi / 360) */ | ||||
141 | Bounds bounds[6]; | ||||
142 | int i; | ||||
143 | | ||||
144 | get_bounds(l, bounds); | ||||
145 | for(i = 0; i < 6; i++) { | ||||
146 | double len = ray_length_until_intersect(hrad, &bounds[i]); | ||||
147 | | ||||
148 | if(len >= 0 && len < min_len) | ||||
149 | min_len = len; | ||||
150 | } | ||||
151 | return min_len; | ||||
152 | } | ||||
153 | | ||||
154 | static double | ||||
155 | dot_product(const Triplet* t1, const Triplet* t2) | ||||
156 | { | ||||
157 | return (t1->a * t2->a + t1->b * t2->b + t1->c * t2->c); | ||||
158 | } | ||||
159 | | ||||
160 | /* Used for rgb conversions */ | ||||
161 | static double | ||||
162 | from_linear(double c) | ||||
163 | { | ||||
164 | if(c <= 0.0031308) | ||||
165 | return 12.92 * c; | ||||
166 | else | ||||
167 | return 1.055 * pow(c, 1.0 / 2.4) - 0.055; | ||||
168 | } | ||||
169 | | ||||
170 | static double | ||||
171 | to_linear(double c) | ||||
172 | { | ||||
173 | if (c > 0.04045) | ||||
174 | return pow((c + 0.055) / 1.055, 2.4); | ||||
175 | else | ||||
176 | return c / 12.92; | ||||
177 | } | ||||
178 | | ||||
179 | static void | ||||
180 | xyz2rgb(Triplet* in_out) | ||||
181 | { | ||||
182 | double r = from_linear(dot_product(&m[0], in_out)); | ||||
183 | double g = from_linear(dot_product(&m[1], in_out)); | ||||
184 | double b = from_linear(dot_product(&m[2], in_out)); | ||||
185 | in_out->a = r; | ||||
186 | in_out->b = g; | ||||
187 | in_out->c = b; | ||||
188 | } | ||||
189 | | ||||
190 | static void | ||||
191 | rgb2xyz(Triplet* in_out) | ||||
192 | { | ||||
193 | Triplet rgbl = { to_linear(in_out->a), to_linear(in_out->b), to_linear(in_out->c) }; | ||||
194 | double x = dot_product(&m_inv[0], &rgbl); | ||||
195 | double y = dot_product(&m_inv[1], &rgbl); | ||||
196 | double z = dot_product(&m_inv[2], &rgbl); | ||||
197 | in_out->a = x; | ||||
198 | in_out->b = y; | ||||
199 | in_out->c = z; | ||||
200 | } | ||||
201 | | ||||
202 | /* http://en.wikipedia.org/wiki/CIELUV | ||||
203 | * In these formulas, Yn refers to the reference white point. We are using | ||||
204 | * illuminant D65, so Yn (see refY in Maxima file) equals 1. The formula is | ||||
205 | * simplified accordingly. | ||||
206 | */ | ||||
207 | static double | ||||
208 | y2l(double y) | ||||
209 | { | ||||
210 | if(y <= epsilon) | ||||
211 | return y * kappa; | ||||
212 | else | ||||
213 | return 116.0 * cbrt(y) - 16.0; | ||||
214 | } | ||||
215 | | ||||
216 | static double | ||||
217 | l2y(double l) | ||||
218 | { | ||||
219 | if(l <= 8.0) { | ||||
220 | return l / kappa; | ||||
221 | } else { | ||||
222 | double x = (l + 16.0) / 116.0; | ||||
223 | return (x * x * x); | ||||
224 | } | ||||
225 | } | ||||
226 | | ||||
227 | static void | ||||
228 | xyz2luv(Triplet* in_out) | ||||
229 | { | ||||
230 | double var_u = (4.0 * in_out->a) / (in_out->a + (15.0 * in_out->b) + (3.0 * in_out->c)); | ||||
231 | double var_v = (9.0 * in_out->b) / (in_out->a + (15.0 * in_out->b) + (3.0 * in_out->c)); | ||||
232 | double l = y2l(in_out->b); | ||||
233 | double u = 13.0 * l * (var_u - ref_u); | ||||
234 | double v = 13.0 * l * (var_v - ref_v); | ||||
235 | | ||||
236 | in_out->a = l; | ||||
237 | if(l < 0.00000001) { | ||||
238 | in_out->b = 0.0; | ||||
239 | in_out->c = 0.0; | ||||
240 | } else { | ||||
241 | in_out->b = u; | ||||
242 | in_out->c = v; | ||||
243 | } | ||||
244 | } | ||||
245 | | ||||
246 | static void | ||||
247 | luv2xyz(Triplet* in_out) | ||||
248 | { | ||||
249 | if(in_out->a <= 0.00000001) { | ||||
250 | /* Black will create a divide-by-zero error. */ | ||||
251 | in_out->a = 0.0; | ||||
252 | in_out->b = 0.0; | ||||
253 | in_out->c = 0.0; | ||||
254 | return; | ||||
255 | } | ||||
256 | | ||||
257 | double var_u = in_out->b / (13.0 * in_out->a) + ref_u; | ||||
258 | double var_v = in_out->c / (13.0 * in_out->a) + ref_v; | ||||
259 | double y = l2y(in_out->a); | ||||
260 | double x = -(9.0 * y * var_u) / ((var_u - 4.0) * var_v - var_u * var_v); | ||||
261 | double z = (9.0 * y - (15.0 * var_v * y) - (var_v * x)) / (3.0 * var_v); | ||||
262 | in_out->a = x; | ||||
263 | in_out->b = y; | ||||
264 | in_out->c = z; | ||||
265 | } | ||||
266 | | ||||
267 | static void | ||||
268 | luv2lch(Triplet* in_out) | ||||
269 | { | ||||
270 | double l = in_out->a; | ||||
271 | double u = in_out->b; | ||||
272 | double v = in_out->c; | ||||
273 | double h; | ||||
274 | double c = sqrt(u * u + v * v); | ||||
275 | | ||||
276 | /* Grays: disambiguate hue */ | ||||
277 | if(c < 0.00000001) { | ||||
278 | h = 0; | ||||
279 | } else { | ||||
280 | h = atan2(v, u) * 57.29577951308232087680; /* (180 / pi) */ | ||||
281 | if(h < 0.0) | ||||
282 | h += 360.0; | ||||
283 | } | ||||
284 | | ||||
285 | in_out->a = l; | ||||
286 | in_out->b = c; | ||||
287 | in_out->c = h; | ||||
288 | } | ||||
289 | | ||||
290 | static void | ||||
291 | lch2luv(Triplet* in_out) | ||||
292 | { | ||||
293 | double hrad = in_out->c * 0.01745329251994329577; /* (pi / 180.0) */ | ||||
294 | double u = cos(hrad) * in_out->b; | ||||
295 | double v = sin(hrad) * in_out->b; | ||||
296 | | ||||
297 | in_out->b = u; | ||||
298 | in_out->c = v; | ||||
299 | } | ||||
300 | | ||||
301 | static void | ||||
302 | hsluv2lch(Triplet* in_out) | ||||
303 | { | ||||
304 | double h = in_out->a; | ||||
305 | double s = in_out->b; | ||||
306 | double l = in_out->c; | ||||
307 | double c; | ||||
308 | | ||||
309 | /* White and black: disambiguate chroma */ | ||||
310 | if(l > 99.9999999 || l < 0.00000001) | ||||
311 | c = 0.0; | ||||
312 | else | ||||
313 | c = max_chroma_for_lh(l, h) / 100.0 * s; | ||||
314 | | ||||
315 | /* Grays: disambiguate hue */ | ||||
316 | if (s < 0.00000001) | ||||
317 | h = 0.0; | ||||
318 | | ||||
319 | in_out->a = l; | ||||
320 | in_out->b = c; | ||||
321 | in_out->c = h; | ||||
322 | } | ||||
323 | | ||||
324 | static void | ||||
325 | lch2hsluv(Triplet* in_out) | ||||
326 | { | ||||
327 | double l = in_out->a; | ||||
328 | double c = in_out->b; | ||||
329 | double h = in_out->c; | ||||
330 | double s; | ||||
331 | | ||||
332 | /* White and black: disambiguate saturation */ | ||||
333 | if(l > 99.9999999 || l < 0.00000001) | ||||
334 | s = 0.0; | ||||
335 | else | ||||
336 | s = c / max_chroma_for_lh(l, h) * 100.0; | ||||
337 | | ||||
338 | /* Grays: disambiguate hue */ | ||||
339 | if (c < 0.00000001) | ||||
340 | h = 0.0; | ||||
341 | | ||||
342 | in_out->a = h; | ||||
343 | in_out->b = s; | ||||
344 | in_out->c = l; | ||||
345 | } | ||||
346 | | ||||
347 | static void | ||||
348 | hpluv2lch(Triplet* in_out) | ||||
349 | { | ||||
350 | double h = in_out->a; | ||||
351 | double s = in_out->b; | ||||
352 | double l = in_out->c; | ||||
353 | double c; | ||||
354 | | ||||
355 | /* White and black: disambiguate chroma */ | ||||
356 | if(l > 99.9999999 || l < 0.00000001) | ||||
357 | c = 0.0; | ||||
358 | else | ||||
359 | c = max_safe_chroma_for_l(l) / 100.0 * s; | ||||
360 | | ||||
361 | /* Grays: disambiguate hue */ | ||||
362 | if (s < 0.00000001) | ||||
363 | h = 0.0; | ||||
364 | | ||||
365 | in_out->a = l; | ||||
366 | in_out->b = c; | ||||
367 | in_out->c = h; | ||||
368 | } | ||||
369 | | ||||
370 | static void | ||||
371 | lch2hpluv(Triplet* in_out) | ||||
372 | { | ||||
373 | double l = in_out->a; | ||||
374 | double c = in_out->b; | ||||
375 | double h = in_out->c; | ||||
376 | double s; | ||||
377 | | ||||
378 | /* White and black: disambiguate saturation */ | ||||
379 | if (l > 99.9999999 || l < 0.00000001) | ||||
380 | s = 0.0; | ||||
381 | else | ||||
382 | s = c / max_safe_chroma_for_l(l) * 100.0; | ||||
383 | | ||||
384 | /* Grays: disambiguate hue */ | ||||
385 | if (c < 0.00000001) | ||||
386 | h = 0.0; | ||||
387 | | ||||
388 | in_out->a = h; | ||||
389 | in_out->b = s; | ||||
390 | in_out->c = l; | ||||
391 | } | ||||
392 | | ||||
393 | | ||||
394 | | ||||
395 | void | ||||
396 | hsluv2rgb(double h, double s, double l, double* pr, double* pg, double* pb) | ||||
397 | { | ||||
398 | Triplet tmp = { h, s, l }; | ||||
399 | | ||||
400 | hsluv2lch(&tmp); | ||||
401 | lch2luv(&tmp); | ||||
402 | luv2xyz(&tmp); | ||||
403 | xyz2rgb(&tmp); | ||||
404 | | ||||
405 | *pr = tmp.a; | ||||
406 | *pg = tmp.b; | ||||
407 | *pb = tmp.c; | ||||
408 | } | ||||
409 | | ||||
410 | void | ||||
411 | hpluv2rgb(double h, double s, double l, double* pr, double* pg, double* pb) | ||||
412 | { | ||||
413 | Triplet tmp = { h, s, l }; | ||||
414 | | ||||
415 | hpluv2lch(&tmp); | ||||
416 | lch2luv(&tmp); | ||||
417 | luv2xyz(&tmp); | ||||
418 | xyz2rgb(&tmp); | ||||
419 | | ||||
420 | *pr = tmp.a; | ||||
421 | *pg = tmp.b; | ||||
422 | *pb = tmp.c; | ||||
423 | } | ||||
424 | | ||||
425 | void | ||||
426 | rgb2hsluv(double r, double g, double b, double* ph, double* ps, double* pl) | ||||
427 | { | ||||
428 | Triplet tmp = { r, g, b }; | ||||
429 | | ||||
430 | rgb2xyz(&tmp); | ||||
431 | xyz2luv(&tmp); | ||||
432 | luv2lch(&tmp); | ||||
433 | lch2hsluv(&tmp); | ||||
434 | | ||||
435 | *ph = tmp.a; | ||||
436 | *ps = tmp.b; | ||||
437 | *pl = tmp.c; | ||||
438 | } | ||||
439 | | ||||
440 | void | ||||
441 | rgb2hpluv(double r, double g, double b, double* ph, double* ps, double* pl) | ||||
442 | { | ||||
443 | Triplet tmp = { r, g, b }; | ||||
444 | | ||||
445 | rgb2xyz(&tmp); | ||||
446 | xyz2luv(&tmp); | ||||
447 | luv2lch(&tmp); | ||||
448 | lch2hpluv(&tmp); | ||||
449 | | ||||
450 | *ph = tmp.a; | ||||
451 | *ps = tmp.b; | ||||
452 | *pl = tmp.c; | ||||
453 | } |