Line data Source code
1 : /**
2 : Copyright (c) 2023-2024 Stappler LLC <admin@stappler.dev>
3 :
4 : Permission is hereby granted, free of charge, to any person obtaining a copy
5 : of this software and associated documentation files (the "Software"), to deal
6 : in the Software without restriction, including without limitation the rights
7 : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 : copies of the Software, and to permit persons to whom the Software is
9 : furnished to do so, subject to the following conditions:
10 :
11 : The above copyright notice and this permission notice shall be included in
12 : all copies or substantial portions of the Software.
13 :
14 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 : THE SOFTWARE.
21 : **/
22 :
23 : #include "SPColorCam16.h"
24 : #include "SPColorHCT.h"
25 :
26 : namespace STAPPLER_VERSIONIZED stappler::geom {
27 :
28 : struct Cam16Vec3 {
29 : Cam16Float a;
30 : Cam16Float b;
31 : Cam16Float c;
32 : };
33 :
34 : constexpr Cam16Float kScaledDiscountFromLinrgb[3][3] = {
35 : {
36 : 0.001200833568784504,
37 : 0.002389694492170889,
38 : 0.0002795742885861124,
39 : },
40 : {
41 : 0.0005891086651375999,
42 : 0.0029785502573438758,
43 : 0.0003270666104008398,
44 : },
45 : {
46 : 0.00010146692491640572,
47 : 0.0005364214359186694,
48 : 0.0032979401770712076,
49 : },
50 : };
51 :
52 : constexpr Cam16Float kLinrgbFromScaledDiscount[3][3] = {
53 : {
54 : 1373.2198709594231,
55 : -1100.4251190754821,
56 : -7.278681089101213,
57 : },
58 : {
59 : -271.815969077903,
60 : 559.6580465940733,
61 : -32.46047482791194,
62 : },
63 : {
64 : 1.9622899599665666,
65 : -57.173814538844006,
66 : 308.7233197812385,
67 : },
68 : };
69 :
70 : constexpr Cam16Float kYFromLinrgb[3] = {0.2126, 0.7152, 0.0722};
71 :
72 : constexpr Cam16Float kCriticalPlanes[255] = {
73 : 0.015176349177441876, 0.045529047532325624, 0.07588174588720938,
74 : 0.10623444424209313, 0.13658714259697685, 0.16693984095186062,
75 : 0.19729253930674434, 0.2276452376616281, 0.2579979360165119,
76 : 0.28835063437139563, 0.3188300904430532, 0.350925934958123,
77 : 0.3848314933096426, 0.42057480301049466, 0.458183274052838,
78 : 0.4976837250274023, 0.5391024159806381, 0.5824650784040898,
79 : 0.6277969426914107, 0.6751227633498623, 0.7244668422128921,
80 : 0.775853049866786, 0.829304845476233, 0.8848452951698498,
81 : 0.942497089126609, 1.0022825574869039, 1.0642236851973577,
82 : 1.1283421258858297, 1.1946592148522128, 1.2631959812511864,
83 : 1.3339731595349034, 1.407011200216447, 1.4823302800086415,
84 : 1.5599503113873272, 1.6398909516233677, 1.7221716113234105,
85 : 1.8068114625156377, 1.8938294463134073, 1.9832442801866852,
86 : 2.075074464868551, 2.1693382909216234, 2.2660538449872063,
87 : 2.36523901573795, 2.4669114995532007, 2.5710888059345764,
88 : 2.6777882626779785, 2.7870270208169257, 2.898822059350997,
89 : 3.0131901897720907, 3.1301480604002863, 3.2497121605402226,
90 : 3.3718988244681087, 3.4967242352587946, 3.624204428461639,
91 : 3.754355295633311, 3.887192587735158, 4.022731918402185,
92 : 4.160988767090289, 4.301978482107941, 4.445716283538092,
93 : 4.592217266055746, 4.741496401646282, 4.893568542229298,
94 : 5.048448422192488, 5.20615066083972, 5.3666897647573375,
95 : 5.5300801301023865, 5.696336044816294, 5.865471690767354,
96 : 6.037501145825082, 6.212438385869475, 6.390297286737924,
97 : 6.571091626112461, 6.7548350853498045, 6.941541251256611,
98 : 7.131223617812143, 7.323895587840543, 7.5195704746346665,
99 : 7.7182615035334345, 7.919981813454504, 8.124744458384042,
100 : 8.332562408825165, 8.543448553206703, 8.757415699253682,
101 : 8.974476575321063, 9.194643831691977, 9.417930041841839,
102 : 9.644347703669503, 9.873909240696694, 10.106627003236781,
103 : 10.342513269534024, 10.58158024687427, 10.8238400726681,
104 : 11.069304815507364, 11.317986476196008, 11.569896988756009,
105 : 11.825048221409341, 12.083451977536606, 12.345119996613247,
106 : 12.610063955123938, 12.878295467455942, 13.149826086772048,
107 : 13.42466730586372, 13.702830557985108, 13.984327217668513,
108 : 14.269168601521828, 14.55736596900856, 14.848930523210871,
109 : 15.143873411576273, 15.44220572664832, 15.743938506781891,
110 : 16.04908273684337, 16.35764934889634, 16.66964922287304,
111 : 16.985093187232053, 17.30399201960269, 17.62635644741625,
112 : 17.95219714852476, 18.281524751807332, 18.614349837764564,
113 : 18.95068293910138, 19.290534541298456, 19.633915083172692,
114 : 19.98083495742689, 20.331304511189067, 20.685334046541502,
115 : 21.042933821039977, 21.404114048223256, 21.76888489811322,
116 : 22.137256497705877, 22.50923893145328, 22.884842241736916,
117 : 23.264076429332462, 23.6469514538663, 24.033477234264016,
118 : 24.42366364919083, 24.817520537484558, 25.21505769858089,
119 : 25.61628489293138, 26.021211842414342, 26.429848230738664,
120 : 26.842203703840827, 27.258287870275353, 27.678110301598522,
121 : 28.10168053274597, 28.529008062403893, 28.96010235337422,
122 : 29.39497283293396, 29.83362889318845, 30.276079891419332,
123 : 30.722335150426627, 31.172403958865512, 31.62629557157785,
124 : 32.08401920991837, 32.54558406207592, 33.010999283389665,
125 : 33.4802739966603, 33.953417292456834, 34.430438229418264,
126 : 34.911345834551085, 35.39614910352207, 35.88485700094671,
127 : 36.37747846067349, 36.87402238606382, 37.37449765026789,
128 : 37.87891309649659, 38.38727753828926, 38.89959975977785,
129 : 39.41588851594697, 39.93615253289054, 40.460400508064545,
130 : 40.98864111053629, 41.520882981230194, 42.05713473317016,
131 : 42.597404951718396, 43.141702194811224, 43.6900349931913,
132 : 44.24241185063697, 44.798841244188324, 45.35933162437017,
133 : 45.92389141541209, 46.49252901546552, 47.065252796817916,
134 : 47.64207110610409, 48.22299226451468, 48.808024568002054,
135 : 49.3971762874833, 49.9904556690408, 50.587870934119984,
136 : 51.189430279724725, 51.79514187861014, 52.40501387947288,
137 : 53.0190544071392, 53.637271562750364, 54.259673423945976,
138 : 54.88626804504493, 55.517063457223934, 56.15206766869424,
139 : 56.79128866487574, 57.43473440856916, 58.08241284012621,
140 : 58.734331877617365, 59.39049941699807, 60.05092333227251,
141 : 60.715611475655585, 61.38457167773311, 62.057811747619894,
142 : 62.7353394731159, 63.417162620860914, 64.10328893648692,
143 : 64.79372614476921, 65.48848194977529, 66.18756403501224,
144 : 66.89098006357258, 67.59873767827808, 68.31084450182222,
145 : 69.02730813691093, 69.74813616640164, 70.47333615344107,
146 : 71.20291564160104, 71.93688215501312, 72.67524319850172,
147 : 73.41800625771542, 74.16517879925733, 74.9167682708136,
148 : 75.67278210128072, 76.43322770089146, 77.1981124613393,
149 : 77.96744375590167, 78.74122893956174, 79.51947534912904,
150 : 80.30219030335869, 81.08938110306934, 81.88105503125999,
151 : 82.67721935322541, 83.4778813166706, 84.28304815182372,
152 : 85.09272707154808, 85.90692527145302, 86.72564993000343,
153 : 87.54890820862819, 88.3767072518277, 89.2090541872801,
154 : 90.04595612594655, 90.88742016217518, 91.73345337380438,
155 : 92.58406282226491, 93.43925555268066, 94.29903859396902,
156 : 95.16341895893969, 96.03240364439274, 96.9059996312159,
157 : 97.78421388448044, 98.6670533535366, 99.55452497210776,
158 : };
159 :
160 1741071 : static float Delinearized(const Cam16Float rgb_component) {
161 1741071 : Cam16Float normalized = rgb_component / 100;
162 : float ret;
163 1741071 : if (normalized <= 0.0031308) {
164 294618 : ret = normalized * 12.92;
165 : } else {
166 1446453 : ret = 1.055 * std::pow(normalized, Cam16Float(1.0 / 2.4)) - 0.055;
167 : }
168 1741071 : return ret;
169 : }
170 :
171 564202 : static Color4F Color4FFromLinrgb(Cam16Vec3 linrgb) {
172 564202 : return Color4F(Delinearized(linrgb.a), Delinearized(linrgb.b), Delinearized(linrgb.c), 1.0f);
173 : }
174 :
175 5889350 : static Cam16Vec3 MatrixMultiply(Cam16Vec3 input, const Cam16Float matrix[3][3]) {
176 5889350 : Cam16Float a = input.a * matrix[0][0] + input.b * matrix[0][1] + input.c * matrix[0][2];
177 5889350 : Cam16Float b = input.a * matrix[1][0] + input.b * matrix[1][1] + input.c * matrix[1][2];
178 5889350 : Cam16Float c = input.a * matrix[2][0] + input.b * matrix[2][1] + input.c * matrix[2][2];
179 5889350 : return {a, b, c};
180 : }
181 :
182 9794604 : static Cam16Float GetAxis(Cam16Vec3 vector, int axis) {
183 9794604 : switch (axis) {
184 3456018 : case 0:
185 3456018 : return vector.a;
186 3587180 : case 1:
187 3587180 : return vector.b;
188 2751406 : case 2:
189 2751406 : return vector.c;
190 0 : default:
191 0 : return -1.0;
192 : }
193 : }
194 :
195 : /**
196 : * Solves the lerp equation.
197 : *
198 : * @param source The starting number.
199 : * @param mid The number in the middle.
200 : * @param target The ending number.
201 : * @return A number t such that lerp(source, target, t) = mid.
202 : */
203 2649392 : Cam16Float Intercept(Cam16Float source, Cam16Float mid, Cam16Float target) {
204 2649392 : return (mid - source) / (target - source);
205 : }
206 :
207 2649392 : Cam16Vec3 LerpPoint(Cam16Vec3 source, Cam16Float t, Cam16Vec3 target) {
208 : return {
209 2649392 : source.a + (target.a - source.a) * t,
210 2649392 : source.b + (target.b - source.b) * t,
211 2649392 : source.c + (target.c - source.c) * t,
212 2649392 : };
213 : }
214 :
215 : /**
216 : * Intersects a segment with a plane.
217 : *
218 : * @param source The coordinates of point A.
219 : * @param coordinate The R-, G-, or B-coordinate of the plane.
220 : * @param target The coordinates of point B.
221 : * @param axis The axis the plane is perpendicular with. (0: R, 1: G, 2: B)
222 : * @return The intersection point of the segment AB with the plane R=coordinate,
223 : * G=coordinate, or B=coordinate
224 : */
225 2649392 : static Cam16Vec3 SetCoordinate(Cam16Vec3 source, Cam16Float coordinate, Cam16Vec3 target, int axis) {
226 2649392 : Cam16Float t = Intercept(GetAxis(source, axis), coordinate, GetAxis(target, axis));
227 2649392 : return LerpPoint(source, t, target);
228 : }
229 :
230 3853560 : static bool IsBounded(Cam16Float x) {
231 3853560 : return 0.0 <= x && x <= 100.0;
232 : }
233 :
234 12654021 : static Cam16Float ChromaticAdaptation(Cam16Float component) {
235 12654021 : Cam16Float af = std::pow(abs(component), Cam16Float(0.42));
236 12654021 : return Cam16::signum(component) * 400.0 * af / (af + 27.13);
237 : }
238 :
239 : /**
240 : * Returns the hue of a linear RGB color in CAM16.
241 : *
242 : * @param linrgb The linear RGB coordinates of a color.
243 : * @return The hue of the color in CAM16, in radians.
244 : */
245 4218007 : static Cam16Float HueOf(Cam16Vec3 linrgb) {
246 4218007 : Cam16Vec3 scaledDiscount = MatrixMultiply(linrgb, kScaledDiscountFromLinrgb);
247 4218007 : Cam16Float r_a = ChromaticAdaptation(scaledDiscount.a);
248 4218007 : Cam16Float g_a = ChromaticAdaptation(scaledDiscount.b);
249 4218007 : Cam16Float b_a = ChromaticAdaptation(scaledDiscount.c);
250 : // redness-greenness
251 4218007 : Cam16Float a = (11.0 * r_a + -12.0 * g_a + b_a) / 11.0;
252 : // yellowness-blueness
253 4218007 : Cam16Float b = (r_a + g_a - 2.0 * b_a) / 9.0;
254 4218007 : return std::atan2(b, a);
255 : }
256 :
257 : /**
258 : * Returns the nth possible vertex of the polygonal intersection.
259 : *
260 : * @param y The Y value of the plane.
261 : * @param n The zero-based index of the point. 0 <= n <= 11.
262 : * @return The nth possible vertex of the polygonal intersection of the y plane
263 : * and the RGB cube, in linear RGB coordinates, if it exists. If this possible
264 : * vertex lies outside of the cube,
265 : * [-1.0, -1.0, -1.0] is returned.
266 : */
267 3853560 : Cam16Vec3 NthVertex(Cam16Float y, int n) {
268 3853560 : const Cam16Float k_r = kYFromLinrgb[0];
269 3853560 : const Cam16Float k_g = kYFromLinrgb[1];
270 3853560 : const Cam16Float k_b = kYFromLinrgb[2];
271 3853560 : const Cam16Float coord_a = n % 4 <= 1 ? 0.0 : 100.0;
272 3853560 : const Cam16Float coord_b = n % 2 == 0 ? 0.0 : 100.0;
273 3853560 : if (n < 4) {
274 1284520 : const Cam16Float g = coord_a;
275 1284520 : const Cam16Float b = coord_b;
276 1284520 : const Cam16Float r = (y - g * k_g - b * k_b) / k_r;
277 1284520 : if (IsBounded(r)) {
278 456363 : return { r, g, b } ;
279 : } else {
280 828157 : return { -1.0, -1.0, -1.0 };
281 : }
282 2569040 : } else if (n < 8) {
283 1284520 : const Cam16Float b = coord_a;
284 1284520 : const Cam16Float r = coord_b;
285 1284520 : const Cam16Float g = (y - r * k_r - b * k_b) / k_g;
286 1284520 : if (IsBounded(g)) {
287 669861 : return { r, g, b };
288 : } else {
289 614659 : return { -1.0, -1.0, -1.0 } ;
290 : }
291 : } else {
292 1284520 : const Cam16Float r = coord_a;
293 1284520 : const Cam16Float g = coord_b;
294 1284520 : const Cam16Float b = (y - r * k_r - g * k_g) / k_b;
295 1284520 : if (IsBounded(b)) {
296 121261 : return {r, g, b};
297 : } else {
298 1163259 : return {-1.0, -1.0, -1.0};
299 : }
300 : }
301 : }
302 :
303 : /**
304 : * Sanitizes a small enough angle in radians.
305 : *
306 : * @param angle An angle in radians; must not deviate too much from 0.
307 : * @return A coterminal angle between 0 and 2pi.
308 : */
309 7976502 : static Cam16Float SanitizeRadians(Cam16Float angle) { return std::fmod(angle + numbers::pi * 8, numbers::pi * 2); }
310 :
311 3988251 : static bool AreInCyclicOrder(Cam16Float a, Cam16Float b, Cam16Float c) {
312 3988251 : Cam16Float delta_a_b = SanitizeRadians(b - a);
313 3988251 : Cam16Float delta_a_c = SanitizeRadians(c - a);
314 3988251 : return delta_a_b < delta_a_c;
315 : }
316 :
317 : /**
318 : * Finds the segment containing the desired color.
319 : *
320 : * @param y The Y value of the color.
321 : * @param target_hue The hue of the color.
322 : * @return A list of two sets of linear RGB coordinates, each corresponding to
323 : * an endpoint of the segment containing the desired color.
324 : */
325 321130 : void BisectToSegment(Cam16Float y, Cam16Float target_hue, Cam16Vec3 out[2]) {
326 321130 : Cam16Vec3 left = { -1.0, -1.0, -1.0 };
327 321130 : Cam16Vec3 right = left;
328 321130 : Cam16Float left_hue = 0.0;
329 321130 : Cam16Float right_hue = 0.0;
330 321130 : bool initialized = false;
331 321130 : bool uncut = true;
332 4174690 : for (int n = 0; n < 12; n++) {
333 3853560 : Cam16Vec3 mid = NthVertex(y, n);
334 3853560 : if (mid.a < 0) {
335 2927205 : continue;
336 : }
337 1247485 : Cam16Float mid_hue = HueOf(mid);
338 1247485 : if (!initialized) {
339 321130 : left = mid;
340 321130 : right = mid;
341 321130 : left_hue = mid_hue;
342 321130 : right_hue = mid_hue;
343 321130 : initialized = true;
344 321130 : continue;
345 : }
346 926355 : if (uncut || AreInCyclicOrder(left_hue, mid_hue, right_hue)) {
347 733634 : uncut = false;
348 733634 : if (AreInCyclicOrder(left_hue, target_hue, mid_hue)) {
349 478709 : right = mid;
350 478709 : right_hue = mid_hue;
351 : } else {
352 254925 : left = mid;
353 254925 : left_hue = mid_hue;
354 : }
355 : }
356 : }
357 321130 : out[0] = left;
358 321130 : out[1] = right;
359 321130 : }
360 :
361 : /**
362 : * Delinearizes an RGB component, returning a floating-point number.
363 : *
364 : * @param rgb_component 0.0 <= rgb_component <= 100.0, represents linear R/G/B
365 : * channel
366 : * @return 0.0 <= output <= 255.0, color channel converted to regular RGB space
367 : */
368 1284520 : static Cam16Float TrueDelinearized(Cam16Float rgb_component) {
369 1284520 : Cam16Float normalized = rgb_component / 100.0;
370 1284520 : Cam16Float delinearized = 0.0;
371 1284520 : if (normalized <= 0.0031308) {
372 200917 : delinearized = normalized * 12.92;
373 : } else {
374 1083603 : delinearized = 1.055 * std::pow(normalized, Cam16Float(1.0 / 2.4)) - 0.055;
375 : }
376 1284520 : return delinearized * 255.0;
377 : }
378 :
379 642260 : static int CriticalPlaneBelow(Cam16Float x) { return (int)std::floor(x - Cam16Float(0.5)); }
380 :
381 642260 : static int CriticalPlaneAbove(Cam16Float x) { return (int)std::ceil(x - Cam16Float(0.5)); }
382 :
383 321130 : static Cam16Vec3 Midpoint(Cam16Vec3 a, Cam16Vec3 b) {
384 : return {
385 321130 : (a.a + b.a) / 2,
386 321130 : (a.b + b.b) / 2,
387 321130 : (a.c + b.c) / 2,
388 321130 : };
389 : }
390 :
391 :
392 : /**
393 : * Finds a color with the given Y and hue on the boundary of the cube.
394 : *
395 : * @param y The Y value of the color.
396 : * @param target_hue The hue of the color.
397 : * @return The desired color, in linear RGB coordinates.
398 : */
399 321130 : static Cam16Vec3 BisectToLimit(Cam16Float y, Cam16Float target_hue) {
400 : Cam16Vec3 segment[2];
401 321130 : BisectToSegment(y, target_hue, segment);
402 321130 : Cam16Vec3 left = segment[0];
403 321130 : Cam16Float left_hue = HueOf(left);
404 321130 : Cam16Vec3 right = segment[1];
405 1284520 : for (int axis = 0; axis < 3; axis++) {
406 963390 : if (GetAxis(left, axis) != GetAxis(right, axis)) {
407 642260 : int l_plane = -1;
408 642260 : int r_plane = 255;
409 642260 : if (GetAxis(left, axis) < GetAxis(right, axis)) {
410 321130 : l_plane = CriticalPlaneBelow(TrueDelinearized(GetAxis(left, axis)));
411 321130 : r_plane = CriticalPlaneAbove(TrueDelinearized(GetAxis(right, axis)));
412 : } else {
413 321130 : l_plane = CriticalPlaneAbove(TrueDelinearized(GetAxis(left, axis)));
414 321130 : r_plane = CriticalPlaneBelow(TrueDelinearized(GetAxis(right, axis)));
415 : }
416 3291652 : for (int i = 0; i < 8; i++) {
417 3218240 : if (abs(r_plane - l_plane) <= 1) {
418 568848 : break;
419 : } else {
420 2649392 : int m_plane = (int) floor((l_plane + r_plane) / 2.0);
421 2649392 : Cam16Float mid_plane_coordinate = kCriticalPlanes[m_plane];
422 2649392 : Cam16Vec3 mid = SetCoordinate(left, mid_plane_coordinate, right, axis);
423 2649392 : Cam16Float mid_hue = HueOf(mid);
424 2649392 : if (AreInCyclicOrder(left_hue, target_hue, mid_hue)) {
425 1514996 : right = mid;
426 1514996 : r_plane = m_plane;
427 : } else {
428 1134396 : left = mid;
429 1134396 : left_hue = mid_hue;
430 1134396 : l_plane = m_plane;
431 : }
432 : }
433 : }
434 : }
435 : }
436 642260 : return Midpoint(left, right);
437 : }
438 :
439 5014029 : static Cam16Float InverseChromaticAdaptation(Cam16Float adapted) {
440 5014029 : Cam16Float adapted_abs = abs(adapted);
441 5014029 : Cam16Float base = std::fmax(Cam16Float(0), Cam16Float(27.13 * adapted_abs / (400.0 - adapted_abs)));
442 5014029 : return Cam16::signum(adapted) * std::pow(base, Cam16Float(1.0 / 0.42));
443 : }
444 :
445 564202 : static Color4F FindResultByJ(Cam16Float hue_radians, Cam16Float chroma, Cam16Float y) {
446 : // Initial estimate of j.
447 564202 : Cam16Float j = std::sqrt(y) * 11.0;
448 : // ===========================================================
449 : // Operations inlined from Cam16 to avoid repeated calculation
450 : // ===========================================================
451 564202 : const Cam16Float t_inner_coeff = 1 / std::pow(Cam16Float(1.64) - std::pow(Cam16Float(0.29),
452 : ViewingConditions::DEFAULT.background_y_to_white_point_y), Cam16Float(0.73));
453 564202 : const Cam16Float e_hue = 0.25 * (std::cos(hue_radians + Cam16Float(2.0)) + 3.8);
454 564202 : const Cam16Float p1 = e_hue * (50000.0 / 13.0) * ViewingConditions::DEFAULT.n_c * ViewingConditions::DEFAULT.ncb;
455 564202 : const Cam16Float h_sin = std::sin(hue_radians);
456 564202 : const Cam16Float h_cos = std::cos(hue_radians);
457 1671343 : for (int iteration_round = 0; iteration_round < 5; ++iteration_round) {
458 : // ===========================================================
459 : // Operations inlined from Cam16 to avoid repeated calculation
460 : // ===========================================================
461 1671343 : Cam16Float j_normalized = j / 100.0;
462 1671343 : Cam16Float alpha = chroma == 0.0 || j == 0.0 ? 0.0 : chroma / sqrt(j_normalized);
463 1671343 : Cam16Float t = std::pow(alpha * t_inner_coeff, Cam16Float(1.0 / 0.9));
464 1671343 : Cam16Float ac = ViewingConditions::DEFAULT.aw * std::pow(j_normalized, Cam16Float(1.0) / ViewingConditions::DEFAULT.c / ViewingConditions::DEFAULT.z);
465 1671343 : Cam16Float p2 = ac / ViewingConditions::DEFAULT.nbb;
466 1671343 : Cam16Float gamma = 23.0 * (p2 + 0.305) * t / (23.0 * p1 + 11 * t * h_cos + 108.0 * t * h_sin);
467 1671343 : Cam16Float a = gamma * h_cos;
468 1671343 : Cam16Float b = gamma * h_sin;
469 1671343 : Cam16Float r_a = (460.0 * p2 + 451.0 * a + 288.0 * b) / 1403.0;
470 1671343 : Cam16Float g_a = (460.0 * p2 - 891.0 * a - 261.0 * b) / 1403.0;
471 1671343 : Cam16Float b_a = (460.0 * p2 - 220.0 * a - 6300.0 * b) / 1403.0;
472 1671343 : Cam16Float r_c_scaled = InverseChromaticAdaptation(r_a);
473 1671343 : Cam16Float g_c_scaled = InverseChromaticAdaptation(g_a);
474 1671343 : Cam16Float b_c_scaled = InverseChromaticAdaptation(b_a);
475 1671343 : Cam16Vec3 scaled { r_c_scaled, g_c_scaled, b_c_scaled };
476 1671343 : Cam16Vec3 linrgb = MatrixMultiply(scaled, kLinrgbFromScaledDiscount);
477 : // ===========================================================
478 : // Operations inlined from Cam16 to avoid repeated calculation
479 : // ===========================================================
480 1671343 : if (linrgb.a < 0 || linrgb.b < 0 || linrgb.c < 0) {
481 564202 : return Color4F::BLACK;
482 : }
483 1402111 : Cam16Float k_r = kYFromLinrgb[0];
484 1402111 : Cam16Float k_g = kYFromLinrgb[1];
485 1402111 : Cam16Float k_b = kYFromLinrgb[2];
486 1402111 : Cam16Float fnj = k_r * linrgb.a + k_g * linrgb.b + k_b * linrgb.c;
487 1402111 : if (fnj <= 0) {
488 0 : return Color4F::BLACK;
489 : }
490 1402111 : if (iteration_round == 4 || abs(fnj - y) < 0.002) {
491 294970 : if (linrgb.a > 100.01 || linrgb.b > 100.01 || linrgb.c > 100.01) {
492 51898 : return Color4F::BLACK;
493 : }
494 243072 : return Color4FFromLinrgb(linrgb);
495 : }
496 : // Iterates with Newton method,
497 : // Using 2 * fn(j) / j as the approximation of fn'(j)
498 1107141 : j = j - (fnj - y) * j / (2 * fnj);
499 : }
500 0 : return Color4F::BLACK;
501 : }
502 :
503 48465 : static Color4F Color4FFromLstar(const Cam16Float lstar) {
504 48465 : auto y = ViewingConditions::YFromLstar(lstar);
505 48465 : auto component = Delinearized(y);
506 48465 : return Color4F(component, component, component, 1.0f);
507 : }
508 :
509 1112 : static float logn(float x, float k) {
510 : float answer;
511 1112 : answer = std::log10(x) / std::log10(k);
512 1112 : return answer;
513 : }
514 :
515 321130 : static void fixTone(Cam16Float &h, Cam16Float &c, Cam16Float &t, Color4F &color) {
516 321130 : const Cam16Float hueOffset = 109.0f;
517 321130 : const Cam16Float hueRange = 30.0f;
518 321130 : const Cam16Float toneOffset = 97.0f;
519 :
520 321130 : if (h >= hueOffset && h <= hueOffset + hueRange && t > toneOffset) {
521 1112 : const float tone = (t - toneOffset) / (100.0f - toneOffset);
522 1112 : const float val = (h - hueOffset) / hueRange * 5.0;
523 1112 : const float log = logn(val, 2.0f);
524 1112 : const float p = std::pow(2.0f, - (log * log));
525 1112 : const float q = (p * 0.95f * std::sqrt(tone));
526 :
527 1112 : color.r = std::pow(color.r, 1.0f - q);
528 1112 : color.g = std::pow(color.g, 1.0f - q);
529 1112 : color.b = std::pow(color.b, 1.0f - q);
530 : }
531 321130 : }
532 :
533 612667 : static Color4F SolveToColor4F(Cam16Float hue_degrees, Cam16Float chroma, Cam16Float lstar) {
534 612667 : if (chroma < 0.0001 || lstar < 0.0001 || lstar > 99.9999) {
535 48465 : return Color4FFromLstar(lstar);
536 : }
537 :
538 564202 : hue_degrees = Cam16::sanitizeDegrees(hue_degrees);
539 : //fixTone(hue_degrees, chroma, lstar);
540 564202 : Cam16Float hue_radians = hue_degrees / 180 * numbers::pi;
541 564202 : Cam16Float y = ViewingConditions::YFromLstar(lstar);
542 564202 : Color4F exact_answer = FindResultByJ(hue_radians, chroma, y);
543 564202 : if (exact_answer != Color4F::BLACK) {
544 243072 : return exact_answer;
545 : }
546 321130 : Cam16Vec3 linrgb = BisectToLimit(y, hue_radians);
547 321130 : auto ret = Color4FFromLinrgb(linrgb);
548 321130 : fixTone(hue_degrees, chroma, lstar, ret);
549 321130 : return ret;
550 : }
551 :
552 11103 : ColorHCT ColorHCT::progress(const ColorHCT &a, const ColorHCT &b, float p) {
553 : return ColorHCT(
554 11103 : (a.data.hue * (1.0f - p) + b.data.hue * p),
555 11103 : (a.data.chroma * (1.0f - p) + b.data.chroma * p),
556 11103 : (a.data.tone * (1.0f - p) + b.data.tone * p),
557 11103 : (a.data.alpha * (1.0f - p) + b.data.alpha * p)
558 11103 : );
559 : }
560 :
561 25 : ColorHCT ColorHCT::solveColorHCT(Cam16Float h, Cam16Float c, Cam16Float t, float a) {
562 25 : auto tmp = SolveToColor4F(h, c, t);
563 25 : tmp.a = a;
564 50 : return ColorHCT(tmp);
565 : }
566 :
567 612642 : Color4F ColorHCT::solveColor4F(Cam16Float h, Cam16Float c, Cam16Float t, float a) {
568 612642 : auto tmp = SolveToColor4F(h, c, t);
569 612642 : tmp.a = a;
570 612642 : return tmp;
571 : }
572 :
573 50 : std::ostream & operator<<(std::ostream & stream, const ColorHCT & obj) {
574 50 : stream << "ColorHCT(h:" << obj.data.hue << " c:" << obj.data.chroma << " t:" << obj.data.tone << " a:" << obj.data.alpha << ");";
575 50 : return stream;
576 : }
577 :
578 : #ifdef __LCC__
579 :
580 : const ViewingConditions ViewingConditions::DEFAULT = {
581 : 11.725676537, 50.000000000, 2.000000000,
582 : false,
583 : 0.184186503, 29.981000900, 1.016919255, 1.016919255, 0.689999998,
584 : 1.000000000, 0.388481468, 0.789482653, 1.909169555,
585 : { 95.047, 100.0, 108.883 },
586 : { 1.021177769, 0.986307740, 0.933960497 }
587 : };
588 :
589 : #endif
590 :
591 : }
|
