libs/base/src/eu/base/colors.cc
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "eu/base/colors.h" | ||
| 2 | |||
| 3 | #include <algorithm> | ||
| 4 | #include <cmath> | ||
| 5 | #include <limits> | ||
| 6 | |||
| 7 | #include "eu/base/numeric.h" | ||
| 8 | |||
| 9 | namespace eu | ||
| 10 | { | ||
| 11 | |||
| 12 | /// converts from gamma space to linear space | ||
| 13 | ✗ | float linear_from_srgb(float value, float gamma) | |
| 14 | { | ||
| 15 | // aka srgb_transfer_function_inv | ||
| 16 | |||
| 17 | |||
| 18 | // todo(Gustav): is this srgb or a basic gamma2 transformation? what's the difference? verify all code usage! | ||
| 19 | ✗ | return std::pow(value, gamma); | |
| 20 | } | ||
| 21 | |||
| 22 | |||
| 23 | ✗ | Lin_rgb linear_from_srgb(const Rgb& value, float gamma) | |
| 24 | { | ||
| 25 | return | ||
| 26 | { | ||
| 27 | ✗ | .r = linear_from_srgb(value.r, gamma), | |
| 28 | ✗ | .g = linear_from_srgb(value.g, gamma), | |
| 29 | ✗ | .b = linear_from_srgb(value.b, gamma) | |
| 30 | ✗ | }; | |
| 31 | } | ||
| 32 | |||
| 33 | |||
| 34 | ✗ | float linear_from_srgb(float a) | |
| 35 | { | ||
| 36 | ✗ | return .04045f < a ? powf((a + .055f) / 1.055f, 2.4f) : a / 12.92f; | |
| 37 | } | ||
| 38 | |||
| 39 | ✗ | Lin_rgb linear_from_srgb(const Rgb& value) | |
| 40 | { | ||
| 41 | ✗ | return {.r = linear_from_srgb(value.r), .g = linear_from_srgb(value.g), .b = linear_from_srgb(value.b)}; | |
| 42 | } | ||
| 43 | |||
| 44 | |||
| 45 | |||
| 46 | ✗ | float srgb_from_linear(float a) | |
| 47 | { | ||
| 48 | ✗ | return .0031308f >= a ? 12.92f * a : 1.055f * powf(a, .4166666666666667f) - .055f; | |
| 49 | } | ||
| 50 | |||
| 51 | |||
| 52 | ✗ | Rgb srgb_from_linear(const Lin_rgb& value) | |
| 53 | { | ||
| 54 | ✗ | return {srgb_from_linear(value.r), srgb_from_linear(value.g), srgb_from_linear(value.b)}; | |
| 55 | } | ||
| 56 | |||
| 57 | struct MaxSaturationCoefficients | ||
| 58 | { | ||
| 59 | float k0; | ||
| 60 | float k1; | ||
| 61 | float k2; | ||
| 62 | float k3; | ||
| 63 | float k4; | ||
| 64 | float wl; | ||
| 65 | float wm; | ||
| 66 | float ws; | ||
| 67 | }; | ||
| 68 | |||
| 69 | // Finds the maximum saturation possible for a given hue that fits in sRGB | ||
| 70 | // Saturation here is defined as S = C/L | ||
| 71 | // a and b must be normalized so a^2 + b^2 == 1 | ||
| 72 | ✗ | float compute_max_saturation(float a, float b) | |
| 73 | { | ||
| 74 | // Max saturation will be when one of r, g or b goes below zero. | ||
| 75 | |||
| 76 | // Select different coefficients depending on which component goes below zero first | ||
| 77 | ✗ | const auto c = ([&]() -> MaxSaturationCoefficients { | |
| 78 | ✗ | if (-1.88170328f * a - 0.80936493f * b > 1) | |
| 79 | { | ||
| 80 | // Red component | ||
| 81 | return | ||
| 82 | { | ||
| 83 | .k0 = +1.19086277f, | ||
| 84 | .k1 = +1.76576728f, | ||
| 85 | .k2 = +0.59662641f, | ||
| 86 | .k3 = +0.75515197f, | ||
| 87 | .k4 = +0.56771245f, | ||
| 88 | .wl = +4.0767416621f, | ||
| 89 | .wm = -3.3077115913f, | ||
| 90 | .ws = +0.2309699292f, | ||
| 91 | ✗ | }; | |
| 92 | } | ||
| 93 | ✗ | else if (1.81444104f * a - 1.19445276f * b > 1) | |
| 94 | { | ||
| 95 | // Green component | ||
| 96 | return | ||
| 97 | { | ||
| 98 | .k0 = +0.73956515f, | ||
| 99 | .k1 = -0.45954404f, | ||
| 100 | .k2 = +0.08285427f, | ||
| 101 | .k3 = +0.12541070f, | ||
| 102 | .k4 = +0.14503204f, | ||
| 103 | .wl = -1.2684380046f, | ||
| 104 | .wm = +2.6097574011f, | ||
| 105 | .ws = -0.3413193965f, | ||
| 106 | ✗ | }; | |
| 107 | } | ||
| 108 | else | ||
| 109 | { | ||
| 110 | // Blue component | ||
| 111 | return | ||
| 112 | { | ||
| 113 | .k0 = +1.35733652f, | ||
| 114 | .k1 = -0.00915799f, | ||
| 115 | .k2 = -1.15130210f, | ||
| 116 | .k3 = -0.50559606f, | ||
| 117 | .k4 = +0.00692167f, | ||
| 118 | .wl = -0.0041960863f, | ||
| 119 | .wm = -0.7034186147f, | ||
| 120 | .ws = +1.7076147010f, | ||
| 121 | ✗ | }; | |
| 122 | } | ||
| 123 | ✗ | })(); | |
| 124 | |||
| 125 | // Approximate max saturation using a polynomial: | ||
| 126 | ✗ | float big_s = c.k0 + c.k1 * a + c.k2 * b + c.k3 * a * a + c.k4 * a * b; | |
| 127 | |||
| 128 | // Do one step Halley's method to get closer | ||
| 129 | // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite | ||
| 130 | // this should be sufficient for most applications, otherwise do two/three steps | ||
| 131 | |||
| 132 | ✗ | const auto k_l = +0.3963377774f * a + 0.2158037573f * b; | |
| 133 | ✗ | const auto k_m = -0.1055613458f * a - 0.0638541728f * b; | |
| 134 | ✗ | const auto k_s = -0.0894841775f * a - 1.2914855480f * b; | |
| 135 | |||
| 136 | { | ||
| 137 | ✗ | const auto l_prim = 1.f + big_s * k_l; | |
| 138 | ✗ | const auto m_prim = 1.f + big_s * k_m; | |
| 139 | ✗ | const auto s_prim = 1.f + big_s * k_s; | |
| 140 | |||
| 141 | ✗ | const auto l = l_prim * l_prim * l_prim; | |
| 142 | ✗ | const auto m = m_prim * m_prim * m_prim; | |
| 143 | ✗ | const auto s = s_prim * s_prim * s_prim; | |
| 144 | |||
| 145 | ✗ | const auto l_d_big_s = 3.f * k_l * l_prim * l_prim; | |
| 146 | ✗ | const auto m_d_big_s = 3.f * k_m * m_prim * m_prim; | |
| 147 | ✗ | const auto s_d_big_s = 3.f * k_s * s_prim * s_prim; | |
| 148 | |||
| 149 | ✗ | const auto l_d_big_s2 = 6.f * k_l * k_l * l_prim; | |
| 150 | ✗ | const auto m_d_big_s2 = 6.f * k_m * k_m * m_prim; | |
| 151 | ✗ | const auto s_d_big_s2 = 6.f * k_s * k_s * s_prim; | |
| 152 | |||
| 153 | ✗ | const auto f = c.wl * l + c.wm * m + c.ws * s; | |
| 154 | ✗ | const auto f1 = c.wl * l_d_big_s + c.wm * m_d_big_s + c.ws * s_d_big_s; | |
| 155 | ✗ | const auto f2 = c.wl * l_d_big_s2 + c.wm * m_d_big_s2 + c.ws * s_d_big_s2; | |
| 156 | |||
| 157 | ✗ | big_s = big_s - f * f1 / (f1 * f1 - 0.5f * f * f2); | |
| 158 | } | ||
| 159 | |||
| 160 | ✗ | return big_s; | |
| 161 | } | ||
| 162 | |||
| 163 | |||
| 164 | struct Lc | ||
| 165 | { | ||
| 166 | float l; | ||
| 167 | float c; | ||
| 168 | }; | ||
| 169 | |||
| 170 | |||
| 171 | /// finds L_cusp and C_cusp for a given hue | ||
| 172 | /// a and b must be normalized so a^2 + b^2 == 1 | ||
| 173 | ✗ | Lc find_cusp(float a, float b) | |
| 174 | { | ||
| 175 | // First, find the maximum saturation (saturation S = C/L) | ||
| 176 | ✗ | const auto big_s_cusp = compute_max_saturation(a, b); | |
| 177 | |||
| 178 | // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1: | ||
| 179 | ✗ | const auto rgb_at_max = linear_from_oklab({.l = 1, .a = big_s_cusp * a, .b = big_s_cusp * b}); | |
| 180 | ✗ | const auto big_l_cusp = cbrtf(1.f / std::max({rgb_at_max.r, rgb_at_max.g, rgb_at_max.b})); | |
| 181 | ✗ | const auto big_c_cusp = big_l_cusp * big_s_cusp; | |
| 182 | |||
| 183 | ✗ | return {.l = big_l_cusp, .c = big_c_cusp}; | |
| 184 | } | ||
| 185 | |||
| 186 | |||
| 187 | /// Finds intersection of the line defined by | ||
| 188 | /// `L = L0 * (1 - t) + t * L1;` | ||
| 189 | /// `C = t * C1;` | ||
| 190 | /// a and b must be normalized so a^2 + b^2 == 1 | ||
| 191 | ✗ | float find_gamut_intersection(float a, float b, float big_l1, float big_c1, float big_l0) | |
| 192 | { | ||
| 193 | // Find the cusp of the gamut triangle | ||
| 194 | ✗ | const auto cusp = find_cusp(a, b); | |
| 195 | |||
| 196 | // Find the intersection for upper and lower half separately | ||
| 197 | ✗ | if (((big_l1 - big_l0) * cusp.c - (cusp.l - big_l0) * big_c1) <= 0.f) | |
| 198 | { | ||
| 199 | // Lower half | ||
| 200 | ✗ | return cusp.c * big_l0 / (big_c1 * cusp.l + cusp.c * (big_l0 - big_l1)); | |
| 201 | } | ||
| 202 | |||
| 203 | // Upper half | ||
| 204 | |||
| 205 | // First intersect with triangle | ||
| 206 | ✗ | float t = cusp.c * (big_l0 - 1.f) / (big_c1 * (cusp.l - 1.f) + cusp.c * (big_l0 - big_l1)); | |
| 207 | |||
| 208 | // Then one-step Halley's method | ||
| 209 | ✗ | const auto d_big_l = big_l1 - big_l0; | |
| 210 | ✗ | const auto d_big_c = big_c1; | |
| 211 | |||
| 212 | ✗ | const auto k_l = +0.3963377774f * a + 0.2158037573f * b; | |
| 213 | ✗ | const auto k_m = -0.1055613458f * a - 0.0638541728f * b; | |
| 214 | ✗ | const auto k_s = -0.0894841775f * a - 1.2914855480f * b; | |
| 215 | |||
| 216 | ✗ | const auto l_dt = d_big_l + d_big_c * k_l; | |
| 217 | ✗ | const auto m_dt = d_big_l + d_big_c * k_m; | |
| 218 | ✗ | const auto s_dt = d_big_l + d_big_c * k_s; | |
| 219 | |||
| 220 | |||
| 221 | // If higher accuracy is required, 2 or 3 iterations of the following block can be used: | ||
| 222 | { | ||
| 223 | ✗ | const auto big_l = big_l0 * (1.f - t) + t * big_l1; | |
| 224 | ✗ | const auto big_c = t * big_c1; | |
| 225 | |||
| 226 | ✗ | const auto l_prim = big_l + big_c * k_l; | |
| 227 | ✗ | const auto m_prim = big_l + big_c * k_m; | |
| 228 | ✗ | const auto s_prim = big_l + big_c * k_s; | |
| 229 | |||
| 230 | ✗ | const auto l = l_prim * l_prim * l_prim; | |
| 231 | ✗ | const auto m = m_prim * m_prim * m_prim; | |
| 232 | ✗ | const auto s = s_prim * s_prim * s_prim; | |
| 233 | |||
| 234 | ✗ | const auto ldt = 3 * l_dt * l_prim * l_prim; | |
| 235 | ✗ | const auto mdt = 3 * m_dt * m_prim * m_prim; | |
| 236 | ✗ | const auto sdt = 3 * s_dt * s_prim * s_prim; | |
| 237 | |||
| 238 | ✗ | const auto ldt2 = 6 * l_dt * l_dt * l_prim; | |
| 239 | ✗ | const auto mdt2 = 6 * m_dt * m_dt * m_prim; | |
| 240 | ✗ | const auto sdt2 = 6 * s_dt * s_dt * s_prim; | |
| 241 | |||
| 242 | ✗ | const auto r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1; | |
| 243 | ✗ | const auto r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt; | |
| 244 | ✗ | const auto r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2; | |
| 245 | |||
| 246 | ✗ | const auto u_r = r1 / (r1 * r1 - 0.5f * r * r2); | |
| 247 | ✗ | auto t_r = -r * u_r; | |
| 248 | |||
| 249 | ✗ | const auto g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1; | |
| 250 | ✗ | const auto g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt; | |
| 251 | ✗ | const auto g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2; | |
| 252 | |||
| 253 | ✗ | const auto u_g = g1 / (g1 * g1 - 0.5f * g * g2); | |
| 254 | ✗ | auto t_g = -g * u_g; | |
| 255 | |||
| 256 | ✗ | const auto b0 = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1; | |
| 257 | ✗ | const auto b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt; | |
| 258 | ✗ | const auto b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2; | |
| 259 | |||
| 260 | ✗ | const auto u_b = b1 / (b1 * b1 - 0.5f * b0 * b2); | |
| 261 | ✗ | auto t_b = -b0 * u_b; | |
| 262 | |||
| 263 | ✗ | t_r = u_r >= 0.f ? t_r : std::numeric_limits<float>::max(); | |
| 264 | ✗ | t_g = u_g >= 0.f ? t_g : std::numeric_limits<float>::max(); | |
| 265 | ✗ | t_b = u_b >= 0.f ? t_b : std::numeric_limits<float>::max(); | |
| 266 | |||
| 267 | ✗ | t += std::min({t_r, t_g, t_b}); | |
| 268 | } | ||
| 269 | |||
| 270 | ✗ | return t; | |
| 271 | } | ||
| 272 | |||
| 273 | |||
| 274 | ✗ | float clamp(float x, float min, float max) | |
| 275 | { | ||
| 276 | ✗ | if (x < min) | |
| 277 | { | ||
| 278 | ✗ | return min; | |
| 279 | } | ||
| 280 | |||
| 281 | ✗ | if (x > max) | |
| 282 | { | ||
| 283 | ✗ | return max; | |
| 284 | } | ||
| 285 | |||
| 286 | ✗ | return x; | |
| 287 | } | ||
| 288 | |||
| 289 | |||
| 290 | ✗ | float sgn(float x) | |
| 291 | { | ||
| 292 | // todo(Gustav): make this clearer | ||
| 293 | ✗ | return static_cast<float>(0.f < x) - static_cast<float>(x < 0.f); | |
| 294 | } | ||
| 295 | |||
| 296 | |||
| 297 | ✗ | Lin_rgb gamut_clip_preserve_chroma(const Lin_rgb& rgb) | |
| 298 | { | ||
| 299 | ✗ | if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && | |
| 300 | ✗ | rgb.r > 0 && rgb.g > 0 && rgb.b > 0) | |
| 301 | { | ||
| 302 | ✗ | return rgb; | |
| 303 | } | ||
| 304 | |||
| 305 | ✗ | const auto lab = oklab_from_linear(rgb); | |
| 306 | |||
| 307 | ✗ | const auto big_l = lab.l; | |
| 308 | ✗ | const auto eps = 0.00001f; | |
| 309 | ✗ | const auto big_c = std::max(eps, std::sqrt(lab.a * lab.a + lab.b * lab.b)); | |
| 310 | ✗ | const auto a_prim = lab.a / big_c; | |
| 311 | ✗ | const auto b_prim = lab.b / big_c; | |
| 312 | |||
| 313 | ✗ | const auto big_l0 = clamp(big_l, 0, 1); | |
| 314 | |||
| 315 | ✗ | const auto t = find_gamut_intersection(a_prim, b_prim, big_l, big_c, big_l0); | |
| 316 | ✗ | const auto big_l_clipped = big_l0 * (1 - t) + t * big_l; | |
| 317 | ✗ | const auto big_c_clipped = t * big_c; | |
| 318 | |||
| 319 | ✗ | return linear_from_oklab({.l = big_l_clipped, .a = big_c_clipped * a_prim, .b = big_c_clipped * b_prim}); | |
| 320 | } | ||
| 321 | |||
| 322 | |||
| 323 | ✗ | Lin_rgb gamut_clip_project_to_0_5(const Lin_rgb& rgb) | |
| 324 | { | ||
| 325 | ✗ | if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && | |
| 326 | ✗ | rgb.r > 0 && rgb.g > 0 && rgb.b > 0) | |
| 327 | { | ||
| 328 | ✗ | return rgb; | |
| 329 | } | ||
| 330 | |||
| 331 | ✗ | const auto lab = oklab_from_linear(rgb); | |
| 332 | |||
| 333 | ✗ | const auto big_l = lab.l; | |
| 334 | ✗ | const auto eps = 0.00001f; | |
| 335 | ✗ | const auto big_c = std::max(eps, std::sqrt(lab.a * lab.a + lab.b * lab.b)); | |
| 336 | ✗ | const auto a_prim = lab.a / big_c; | |
| 337 | ✗ | const auto b_prim = lab.b / big_c; | |
| 338 | |||
| 339 | ✗ | const auto big_l0 = 0.5f; | |
| 340 | |||
| 341 | ✗ | const auto t = find_gamut_intersection(a_prim, b_prim, big_l, big_c, big_l0); | |
| 342 | ✗ | const auto big_l_clipped = big_l0 * (1 - t) + t * big_l; | |
| 343 | ✗ | const auto big_c_clipped = t * big_c; | |
| 344 | |||
| 345 | ✗ | return linear_from_oklab({.l = big_l_clipped, .a = big_c_clipped * a_prim, .b = big_c_clipped * b_prim}); | |
| 346 | } | ||
| 347 | |||
| 348 | |||
| 349 | ✗ | Lin_rgb gamut_clip_project_to_l_cusp(const Lin_rgb& rgb) | |
| 350 | { | ||
| 351 | ✗ | if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && | |
| 352 | ✗ | rgb.r > 0 && rgb.g > 0 && rgb.b > 0) | |
| 353 | { | ||
| 354 | ✗ | return rgb; | |
| 355 | } | ||
| 356 | |||
| 357 | ✗ | const auto lab = oklab_from_linear(rgb); | |
| 358 | |||
| 359 | ✗ | const auto big_l = lab.l; | |
| 360 | ✗ | const auto eps = 0.00001f; | |
| 361 | ✗ | const auto big_c = std::max(eps, std::sqrt(lab.a * lab.a + lab.b * lab.b)); | |
| 362 | ✗ | const auto a_prim = lab.a / big_c; | |
| 363 | ✗ | const auto b_prim = lab.b / big_c; | |
| 364 | |||
| 365 | // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. | ||
| 366 | ✗ | const auto cusp = find_cusp(a_prim, b_prim); | |
| 367 | |||
| 368 | ✗ | const auto big_l0 = cusp.l; | |
| 369 | |||
| 370 | ✗ | const auto t = find_gamut_intersection(a_prim, b_prim, big_l, big_c, big_l0); | |
| 371 | |||
| 372 | ✗ | const auto big_l_clipped = big_l0 * (1 - t) + t * big_l; | |
| 373 | ✗ | const auto big_c_clipped = t * big_c; | |
| 374 | |||
| 375 | ✗ | return linear_from_oklab({.l = big_l_clipped, .a = big_c_clipped * a_prim, .b = big_c_clipped * b_prim}); | |
| 376 | } | ||
| 377 | |||
| 378 | |||
| 379 | ✗ | Lin_rgb gamut_clip_adaptive_L0_0_5(const Lin_rgb& rgb, float alpha) | |
| 380 | { | ||
| 381 | ✗ | if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && | |
| 382 | ✗ | rgb.r > 0 && rgb.g > 0 && rgb.b > 0) | |
| 383 | { | ||
| 384 | ✗ | return rgb; | |
| 385 | } | ||
| 386 | |||
| 387 | ✗ | const auto lab = oklab_from_linear(rgb); | |
| 388 | |||
| 389 | ✗ | const auto big_l = lab.l; | |
| 390 | ✗ | const auto eps = 0.00001f; | |
| 391 | ✗ | const auto big_c = std::max(eps, std::sqrt(lab.a * lab.a + lab.b * lab.b)); | |
| 392 | ✗ | const auto a_prim = lab.a / big_c; | |
| 393 | ✗ | const auto b_prim = lab.b / big_c; | |
| 394 | |||
| 395 | ✗ | const auto big_ld = big_l - 0.5f; | |
| 396 | ✗ | const auto e1 = 0.5f + std::abs(big_ld) + alpha * big_c; | |
| 397 | ✗ | const auto big_l0 = 0.5f * (1.f + sgn(big_ld) * (e1 - std::sqrt(e1 * e1 - 2.f * std::abs(big_ld)))); | |
| 398 | |||
| 399 | ✗ | const auto t = find_gamut_intersection(a_prim, b_prim, big_l, big_c, big_l0); | |
| 400 | ✗ | const auto big_l_clipped = big_l0 * (1.f - t) + t * big_l; | |
| 401 | ✗ | const auto big_c_clipped = t * big_c; | |
| 402 | |||
| 403 | ✗ | return linear_from_oklab({big_l_clipped, big_c_clipped * a_prim, big_c_clipped * b_prim}); | |
| 404 | } | ||
| 405 | |||
| 406 | |||
| 407 | ✗ | Lin_rgb gamut_clip_adaptive_l0_l_cusp(const Lin_rgb& rgb, float alpha) | |
| 408 | { | ||
| 409 | ✗ | if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && | |
| 410 | ✗ | rgb.r > 0 && rgb.g > 0 && rgb.b > 0) | |
| 411 | { | ||
| 412 | ✗ | return rgb; | |
| 413 | } | ||
| 414 | |||
| 415 | ✗ | const auto lab = oklab_from_linear(rgb); | |
| 416 | |||
| 417 | ✗ | const auto big_l = lab.l; | |
| 418 | ✗ | const auto eps = 0.00001f; | |
| 419 | ✗ | const auto big_c = std::max(eps, std::sqrt(lab.a * lab.a + lab.b * lab.b)); | |
| 420 | ✗ | const auto a_prim = lab.a / big_c; | |
| 421 | ✗ | const auto b_prim = lab.b / big_c; | |
| 422 | |||
| 423 | // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. | ||
| 424 | ✗ | const auto cusp = find_cusp(a_prim, b_prim); | |
| 425 | |||
| 426 | ✗ | const auto big_ld = big_l - cusp.l; | |
| 427 | ✗ | const auto k = 2.f * (big_ld > 0 ? 1.f - cusp.l : cusp.l); | |
| 428 | |||
| 429 | ✗ | const auto e1 = 0.5f * k + std::abs(big_ld) + alpha * big_c / k; | |
| 430 | ✗ | const auto big_l0 = cusp.l + 0.5f * (sgn(big_ld) * (e1 - std::sqrt(e1 * e1 - 2.f * k * std::abs(big_ld)))); | |
| 431 | |||
| 432 | ✗ | const auto t = find_gamut_intersection(a_prim, b_prim, big_l, big_c, big_l0); | |
| 433 | ✗ | const auto big_l_clipped = big_l0 * (1.f - t) + t * big_l; | |
| 434 | ✗ | const auto big_c_clipped = t * big_c; | |
| 435 | |||
| 436 | ✗ | return linear_from_oklab({.l = big_l_clipped, .a = big_c_clipped * a_prim, .b = big_c_clipped * b_prim}); | |
| 437 | } | ||
| 438 | |||
| 439 | |||
| 440 | ✗ | OkLab oklab_from_linear(const Lin_rgb& c) | |
| 441 | { | ||
| 442 | ✗ | const auto cr = c.r; | |
| 443 | ✗ | const auto cg = c.g; | |
| 444 | ✗ | const auto cb = c.b; | |
| 445 | |||
| 446 | ✗ | const auto l = 0.4122214708f * cr + 0.5363325363f * cg + 0.0514459929f * cb; | |
| 447 | ✗ | const auto m = 0.2119034982f * cr + 0.6806995451f * cg + 0.1073969566f * cb; | |
| 448 | ✗ | const auto s = 0.0883024619f * cr + 0.2817188376f * cg + 0.6299787005f * cb; | |
| 449 | |||
| 450 | ✗ | const auto lp = std::cbrtf(l); | |
| 451 | ✗ | const auto mp = std::cbrtf(m); | |
| 452 | ✗ | const auto sp = std::cbrtf(s); | |
| 453 | |||
| 454 | return { | ||
| 455 | ✗ | .l = 0.2104542553f * lp + 0.7936177850f * mp - 0.0040720468f * sp, | |
| 456 | ✗ | .a = 1.9779984951f * lp - 2.4285922050f * mp + 0.4505937099f * sp, | |
| 457 | ✗ | .b = 0.0259040371f * lp + 0.7827717662f * mp - 0.8086757660f * sp, | |
| 458 | ✗ | }; | |
| 459 | } | ||
| 460 | |||
| 461 | ✗ | Lin_rgb linear_from_oklab(const OkLab& c) | |
| 462 | { | ||
| 463 | ✗ | const auto lp = c.l + 0.3963377774f * c.a + 0.2158037573f * c.b; | |
| 464 | ✗ | const auto mp = c.l - 0.1055613458f * c.a - 0.0638541728f * c.b; | |
| 465 | ✗ | const auto sp = c.l - 0.0894841775f * c.a - 1.2914855480f * c.b; | |
| 466 | |||
| 467 | ✗ | const auto l = lp * lp * lp; | |
| 468 | ✗ | const auto m = mp * mp * mp; | |
| 469 | ✗ | const auto s = sp * sp * sp; | |
| 470 | |||
| 471 | return { | ||
| 472 | ✗ | .r = +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s, | |
| 473 | ✗ | .g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s, | |
| 474 | ✗ | .b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s, | |
| 475 | ✗ | }; | |
| 476 | } | ||
| 477 | |||
| 478 | // https://en.wikipedia.org/wiki/Oklab_color_space#Conversion_to_and_from_Oklch | ||
| 479 | ✗ | OkLch oklch_from_oklab(const OkLab& c) | |
| 480 | { | ||
| 481 | ✗ | return {.l = c.l, .c = std::sqrt(c.a * c.a + c.b * c.b), .h = eu::atan2(c.b, c.a)}; | |
| 482 | } | ||
| 483 | |||
| 484 | ✗ | OkLab oklab_from_oklch(const OkLch& c) | |
| 485 | { | ||
| 486 | return | ||
| 487 | { | ||
| 488 | ✗ | .l = c.l, .a = c.c *eu::cos(c.h), .b = c.c *eu::sin(c.h) | |
| 489 | ✗ | }; | |
| 490 | } | ||
| 491 | |||
| 492 | |||
| 493 | // Alternative representation of (L_cusp, C_cusp) | ||
| 494 | // Encoded so S = C_cusp/L_cusp and T = C_cusp/(1-L_cusp) | ||
| 495 | // The maximum value for C in the triangle is then found as std::min(S*L, T*(1-L)), for a given L | ||
| 496 | struct St | ||
| 497 | { | ||
| 498 | float s; | ||
| 499 | float t; | ||
| 500 | }; | ||
| 501 | |||
| 502 | // toe function for L_r | ||
| 503 | ✗ | float toe(float x) | |
| 504 | { | ||
| 505 | ✗ | constexpr float k_1 = 0.206f; | |
| 506 | ✗ | constexpr float k_2 = 0.03f; | |
| 507 | ✗ | constexpr float k_3 = (1.f + k_1) / (1.f + k_2); | |
| 508 | ✗ | return 0.5f * (k_3 * x - k_1 + std::sqrt((k_3 * x - k_1) * (k_3 * x - k_1) + 4 * k_2 * k_3 * x)); | |
| 509 | } | ||
| 510 | |||
| 511 | // inverse toe function for L_r | ||
| 512 | ✗ | float toe_inv(float x) | |
| 513 | { | ||
| 514 | ✗ | constexpr float k_1 = 0.206f; | |
| 515 | ✗ | constexpr float k_2 = 0.03f; | |
| 516 | ✗ | constexpr float k_3 = (1.f + k_1) / (1.f + k_2); | |
| 517 | ✗ | return (x * x + k_1 * x) / (k_3 * (x + k_2)); | |
| 518 | } | ||
| 519 | |||
| 520 | ✗ | St st_from_cusp(const Lc& cusp) | |
| 521 | { | ||
| 522 | ✗ | const auto big_l = cusp.l; | |
| 523 | ✗ | const auto big_c = cusp.c; | |
| 524 | ✗ | return {.s = big_c / big_l, .t = big_c / (1 - big_l)}; | |
| 525 | } | ||
| 526 | |||
| 527 | ✗ | Rgb srgb_from_okhsv(const OkHsv& hsv) | |
| 528 | { | ||
| 529 | ✗ | const auto h = hsv.hue; | |
| 530 | ✗ | const auto s = hsv.saturation; | |
| 531 | ✗ | const auto v = hsv.value; | |
| 532 | |||
| 533 | ✗ | const auto a_prim = eu::cos(h); | |
| 534 | ✗ | const auto b_prim = eu::sin(h); | |
| 535 | |||
| 536 | ✗ | const auto cusp = find_cusp(a_prim, b_prim); | |
| 537 | ✗ | const auto big_st_max = st_from_cusp(cusp); | |
| 538 | ✗ | const auto big_s_max = big_st_max.s; | |
| 539 | ✗ | const auto big_t_max = big_st_max.t; | |
| 540 | ✗ | const auto big_s_0 = 0.5f; | |
| 541 | ✗ | const auto k = 1 - big_s_0 / big_s_max; | |
| 542 | |||
| 543 | // first we compute L and V as if the gamut is a perfect triangle: | ||
| 544 | |||
| 545 | // L, C when v==1: | ||
| 546 | ✗ | const auto big_l_v = 1 - s * big_s_0 / (big_s_0 + big_t_max - big_t_max * k * s); | |
| 547 | ✗ | const auto big_c_v = s * big_t_max * big_s_0 / (big_s_0 + big_t_max - big_t_max * k * s); | |
| 548 | |||
| 549 | ✗ | auto big_l = v * big_l_v; | |
| 550 | ✗ | auto big_c = v * big_c_v; | |
| 551 | |||
| 552 | // then we compensate for both toe and the curved top part of the triangle: | ||
| 553 | ✗ | const auto big_l_vt = toe_inv(big_l_v); | |
| 554 | ✗ | const auto big_c_vt = big_c_v * big_l_vt / big_l_v; | |
| 555 | |||
| 556 | ✗ | const auto big_l_new = toe_inv(big_l); | |
| 557 | ✗ | big_c = big_c * big_l_new / big_l; | |
| 558 | ✗ | big_l = big_l_new; | |
| 559 | |||
| 560 | ✗ | const auto rgb_scale = linear_from_oklab({.l = big_l_vt, .a = a_prim * big_c_vt, .b = b_prim * big_c_vt}); | |
| 561 | ✗ | const auto scale_big_l = cbrtf(1.f / std::max({rgb_scale.r, rgb_scale.g, rgb_scale.b, 0.f})); | |
| 562 | |||
| 563 | ✗ | big_l = big_l * scale_big_l; | |
| 564 | ✗ | big_c = big_c * scale_big_l; | |
| 565 | |||
| 566 | ✗ | const auto rgb = linear_from_oklab({big_l, big_c * a_prim, big_c * b_prim}); | |
| 567 | ✗ | return srgb_from_linear(rgb); | |
| 568 | } | ||
| 569 | |||
| 570 | ✗ | OkHsv okhsv_from_srgb(const Rgb& rgb) | |
| 571 | { | ||
| 572 | ✗ | const auto lab = oklab_from_linear(linear_from_srgb(rgb)); | |
| 573 | |||
| 574 | ✗ | auto big_c = std::sqrt(lab.a * lab.a + lab.b * lab.b); | |
| 575 | ✗ | const auto a_prim = lab.a / big_c; | |
| 576 | ✗ | const auto b_prim = lab.b / big_c; | |
| 577 | |||
| 578 | ✗ | auto big_l = lab.l; | |
| 579 | ✗ | const auto h = eu::atan2(-lab.b, -lab.a); | |
| 580 | |||
| 581 | ✗ | const auto cusp = find_cusp(a_prim, b_prim); | |
| 582 | ✗ | const auto big_st_max = st_from_cusp(cusp); | |
| 583 | ✗ | const auto big_s_max = big_st_max.s; | |
| 584 | ✗ | const auto big_t_max = big_st_max.t; | |
| 585 | ✗ | const auto big_s_0 = 0.5f; | |
| 586 | ✗ | const auto k = 1 - big_s_0 / big_s_max; | |
| 587 | |||
| 588 | // first we find L_v, C_v, L_vt and C_vt | ||
| 589 | |||
| 590 | ✗ | const auto t = big_t_max / (big_c + big_l * big_t_max); | |
| 591 | ✗ | const auto big_l_v = t * big_l; | |
| 592 | ✗ | const auto big_c_v = t * big_c; | |
| 593 | |||
| 594 | ✗ | const auto big_l_vt = toe_inv(big_l_v); | |
| 595 | ✗ | const auto big_c_vt = big_c_v * big_l_vt / big_l_v; | |
| 596 | |||
| 597 | // we can then use these to invert the step that compensates for the toe and the curved top part of the triangle: | ||
| 598 | ✗ | const auto rgb_scale = linear_from_oklab({.l = big_l_vt, .a = a_prim * big_c_vt, .b = b_prim * big_c_vt}); | |
| 599 | ✗ | const auto scale_big_l = cbrtf(1.f / std::max(std::max(rgb_scale.r, rgb_scale.g), std::max(rgb_scale.b, 0.f))); | |
| 600 | |||
| 601 | ✗ | big_l = big_l / scale_big_l; | |
| 602 | // big_c = big_c / scale_big_l; | ||
| 603 | |||
| 604 | // big_c = big_c * toe(big_l) / big_l; | ||
| 605 | ✗ | big_l = toe(big_l); | |
| 606 | |||
| 607 | // we can now compute v and s: | ||
| 608 | |||
| 609 | ✗ | const auto v = big_l / big_l_v; | |
| 610 | ✗ | const auto s = (big_s_0 + big_t_max) * big_c_v / ((big_t_max * big_s_0) + big_t_max * k * big_c_v); | |
| 611 | |||
| 612 | ✗ | return {.hue = h, .saturation = s, .value = v}; | |
| 613 | } | ||
| 614 | |||
| 615 | |||
| 616 | // Returns a smooth approximation of the location of the cusp | ||
| 617 | // This polynomial was created by an optimization process | ||
| 618 | // It has been designed so that S_mid < S_max and T_mid < T_max | ||
| 619 | ✗ | St get_st_mid(float a_prim, float b_prim) | |
| 620 | { | ||
| 621 | ✗ | const auto big_s = 0.11516993f | |
| 622 | + 1.f | ||
| 623 | ✗ | / (+7.44778970f + 4.15901240f * b_prim | |
| 624 | ✗ | + a_prim | |
| 625 | ✗ | * (-2.19557347f + 1.75198401f * b_prim | |
| 626 | ✗ | + a_prim | |
| 627 | ✗ | * (-2.13704948f - 10.02301043f * b_prim | |
| 628 | ✗ | + a_prim * (-4.24894561f + 5.38770819f * b_prim + 4.69891013f * a_prim)))); | |
| 629 | |||
| 630 | ✗ | const auto big_t = 0.11239642f | |
| 631 | + 1.f | ||
| 632 | ✗ | / (+1.61320320f - 0.68124379f * b_prim | |
| 633 | ✗ | + a_prim | |
| 634 | ✗ | * (+0.40370612f + 0.90148123f * b_prim | |
| 635 | ✗ | + a_prim | |
| 636 | ✗ | * (-0.27087943f + 0.61223990f * b_prim | |
| 637 | ✗ | + a_prim * (+0.00299215f - 0.45399568f * b_prim - 0.14661872f * a_prim)))); | |
| 638 | |||
| 639 | ✗ | return {.s = big_s, .t = big_t}; | |
| 640 | } | ||
| 641 | |||
| 642 | struct Cs | ||
| 643 | { | ||
| 644 | float c_0; | ||
| 645 | float c_mid; | ||
| 646 | float c_max; | ||
| 647 | }; | ||
| 648 | |||
| 649 | ✗ | Cs get_cs(float big_l, float a_prim, float b_prim) | |
| 650 | { | ||
| 651 | ✗ | const auto cusp = find_cusp(a_prim, b_prim); | |
| 652 | |||
| 653 | ✗ | const auto big_c_max = find_gamut_intersection(a_prim, b_prim, big_l, 1, big_l); | |
| 654 | ✗ | const auto big_st_max = st_from_cusp(cusp); | |
| 655 | |||
| 656 | // Scale factor to compensate for the curved part of gamut shape: | ||
| 657 | ✗ | const auto k = big_c_max / std::min((big_l * big_st_max.s), (1 - big_l) * big_st_max.t); | |
| 658 | |||
| 659 | ✗ | const auto big_c_mid = [&] | |
| 660 | { | ||
| 661 | ✗ | const auto big_st_mid = get_st_mid(a_prim, b_prim); | |
| 662 | |||
| 663 | // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. | ||
| 664 | ✗ | const auto big_c_a = big_l * big_st_mid.s; | |
| 665 | ✗ | const auto big_c_b = (1.f - big_l) * big_st_mid.t; | |
| 666 | ✗ | return 0.9f * k * std::sqrt(std::sqrt(1.f / (1.f / (big_c_a * big_c_a * big_c_a * big_c_a) + 1.f / (big_c_b * big_c_b * big_c_b * big_c_b)))); | |
| 667 | ✗ | }(); | |
| 668 | |||
| 669 | ✗ | const auto big_c_0 = [&] | |
| 670 | { | ||
| 671 | // for C_0, the shape is independent of hue, so ST are constant. Values picked to roughly be the average values of ST. | ||
| 672 | ✗ | const auto big_c_a = big_l * 0.4f; | |
| 673 | ✗ | const auto big_c_b = (1.f - big_l) * 0.8f; | |
| 674 | |||
| 675 | // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. | ||
| 676 | ✗ | return std::sqrt(1.f / (1.f / (big_c_a * big_c_a) + 1.f / (big_c_b * big_c_b))); | |
| 677 | ✗ | }(); | |
| 678 | |||
| 679 | ✗ | return {.c_0 = big_c_0, .c_mid = big_c_mid, .c_max = big_c_max}; | |
| 680 | } | ||
| 681 | |||
| 682 | ✗ | Rgb srgb_from_okhsl(const OkHsl& hsl) | |
| 683 | { | ||
| 684 | ✗ | const auto h = hsl.hue; | |
| 685 | ✗ | const auto s = hsl.saturation; | |
| 686 | ✗ | const auto l = hsl.lightness; | |
| 687 | |||
| 688 | ✗ | if (l == 1.0f) | |
| 689 | { | ||
| 690 | ✗ | return {1.f, 1.f, 1.f}; | |
| 691 | } | ||
| 692 | |||
| 693 | ✗ | if (l == 0.f) | |
| 694 | { | ||
| 695 | ✗ | return {0.f, 0.f, 0.f}; | |
| 696 | } | ||
| 697 | |||
| 698 | ✗ | const auto ap = eu::cos(h); | |
| 699 | ✗ | const auto bp = eu::sin(h); | |
| 700 | ✗ | const auto big_l = toe_inv(l); | |
| 701 | |||
| 702 | ✗ | const auto cs = get_cs(big_l, ap, bp); | |
| 703 | ✗ | const auto c_0 = cs.c_0; | |
| 704 | ✗ | const auto c_mid = cs.c_mid; | |
| 705 | ✗ | const auto c_max = cs.c_max; | |
| 706 | |||
| 707 | // Interpolate the three values for C so that: | ||
| 708 | // At s=0: dC/ds = C_0, C=0 | ||
| 709 | // At s=0.8: C=C_mid | ||
| 710 | // At s=1.0: C=C_max | ||
| 711 | |||
| 712 | ✗ | const auto mid = 0.8f; | |
| 713 | ✗ | const auto mid_inv = 1.25f; | |
| 714 | |||
| 715 | ✗ | const auto big_c = [&]() { | |
| 716 | ✗ | if (s < mid) | |
| 717 | { | ||
| 718 | ✗ | const auto t = mid_inv * s; | |
| 719 | |||
| 720 | ✗ | const auto k_1 = mid * c_0; | |
| 721 | ✗ | const auto k_2 = (1.f - k_1 / c_mid); | |
| 722 | |||
| 723 | ✗ | return t * k_1 / (1.f - k_2 * t); | |
| 724 | } | ||
| 725 | else | ||
| 726 | { | ||
| 727 | ✗ | const auto t = (s - mid) / (1 - mid); | |
| 728 | |||
| 729 | ✗ | const auto k_0 = c_mid; | |
| 730 | ✗ | const auto k_1 = (1.f - mid) * c_mid * c_mid * mid_inv * mid_inv / c_0; | |
| 731 | ✗ | const auto k_2 = (1.f - (k_1) / (c_max - c_mid)); | |
| 732 | |||
| 733 | ✗ | return k_0 + t * k_1 / (1.f - k_2 * t); | |
| 734 | } | ||
| 735 | ✗ | }(); | |
| 736 | |||
| 737 | ✗ | const auto rgb = linear_from_oklab({.l = big_l, .a = big_c * ap, .b = big_c * bp}); | |
| 738 | ✗ | return srgb_from_linear(rgb); | |
| 739 | } | ||
| 740 | |||
| 741 | ✗ | OkHsl okhsl_from_srgb(const Rgb& rgb) | |
| 742 | { | ||
| 743 | ✗ | const auto lab = oklab_from_linear(linear_from_srgb(rgb)); | |
| 744 | |||
| 745 | ✗ | const auto big_c = std::sqrt(lab.a * lab.a + lab.b * lab.b); | |
| 746 | ✗ | const auto a_prim = lab.a / big_c; | |
| 747 | ✗ | const auto b_prim = lab.b / big_c; | |
| 748 | |||
| 749 | ✗ | const auto big_l = lab.l; | |
| 750 | ✗ | const auto h = eu::atan2(-lab.b, -lab.a); | |
| 751 | |||
| 752 | ✗ | const auto cs = get_cs(big_l, a_prim, b_prim); | |
| 753 | ✗ | const auto big_c_0 = cs.c_0; | |
| 754 | ✗ | const auto big_c_mid = cs.c_mid; | |
| 755 | ✗ | const auto big_c_max = cs.c_max; | |
| 756 | |||
| 757 | // Inverse of the interpolation in okhsl_to_srgb: | ||
| 758 | |||
| 759 | ✗ | const auto mid = 0.8f; | |
| 760 | ✗ | const auto mid_inv = 1.25f; | |
| 761 | |||
| 762 | ✗ | const float s = [&]() { | |
| 763 | ✗ | if (big_c < big_c_mid) | |
| 764 | { | ||
| 765 | ✗ | const auto k_1 = mid * big_c_0; | |
| 766 | ✗ | const auto k_2 = (1.f - k_1 / big_c_mid); | |
| 767 | |||
| 768 | ✗ | const auto t = big_c / (k_1 + k_2 * big_c); | |
| 769 | ✗ | return t * mid; | |
| 770 | } | ||
| 771 | else | ||
| 772 | { | ||
| 773 | ✗ | const auto k_0 = big_c_mid; | |
| 774 | ✗ | const auto k_1 = (1.f - mid) * big_c_mid * big_c_mid * mid_inv * mid_inv / big_c_0; | |
| 775 | ✗ | const auto k_2 = (1.f - (k_1) / (big_c_max - big_c_mid)); | |
| 776 | |||
| 777 | ✗ | const auto t = (big_c - k_0) / (k_1 + k_2 * (big_c - k_0)); | |
| 778 | ✗ | return mid + (1.f - mid) * t; | |
| 779 | ✗ | }}(); | |
| 780 | |||
| 781 | ✗ | const auto l = toe(big_l); | |
| 782 | ✗ | return {.hue = h, .saturation = s, .lightness = l}; | |
| 783 | } | ||
| 784 | |||
| 785 | ✗ | Lin_rgb keep_within(Lin_rgb c) | |
| 786 | { | ||
| 787 | ✗ | const auto ret = Lin_rgb{.r = keep_within01(c.r), .g = keep_within01(c.g), .b = keep_within01(c.b)}; | |
| 788 | ✗ | return ret; | |
| 789 | } | ||
| 790 | |||
| 791 | ✗ | Rgb srgb_from_hsl(const Hsl& hsl) | |
| 792 | { | ||
| 793 | // https://www.rapidtables.com/convert/color/hsl-to-rgb.html | ||
| 794 | |||
| 795 | ✗ | const auto h = hsl.hue.as_degrees(); | |
| 796 | ✗ | const auto s = hsl.saturation; | |
| 797 | ✗ | const auto l = hsl.lightness; | |
| 798 | |||
| 799 | ✗ | const auto c = (1.0f - std::fabs(2.0f * l - 1.0f)) * s; | |
| 800 | ✗ | const auto x = c * (1.0f - std::fabs(std::fmod(h / 60.0f, 2.0f) - 1.0f)); | |
| 801 | |||
| 802 | ✗ | const auto m = l - c / 2.0f; | |
| 803 | |||
| 804 | ✗ | const auto ret = [m](float r1, float g1, float b1) | |
| 805 | { | ||
| 806 | ✗ | const auto r = r1 + m; | |
| 807 | ✗ | const auto g = g1 + m; | |
| 808 | ✗ | const auto b = b1 + m; | |
| 809 | |||
| 810 | ✗ | return Rgb{r, g, b}; | |
| 811 | ✗ | }; | |
| 812 | |||
| 813 | ✗ | if (h < 60.0f) { | |
| 814 | ✗ | return ret(c, x, 0); | |
| 815 | } | ||
| 816 | ✗ | if (h < 120.0f) { | |
| 817 | ✗ | return ret(x, c, 0); | |
| 818 | } | ||
| 819 | ✗ | if (h < 180.0f) { | |
| 820 | ✗ | return ret(0, c, x); | |
| 821 | } | ||
| 822 | ✗ | if (h < 240.0f) { | |
| 823 | ✗ | return ret(0, x, c); | |
| 824 | } | ||
| 825 | ✗ | if (h < 300.0f) { | |
| 826 | ✗ | return ret(x, 0, c); | |
| 827 | } | ||
| 828 | ✗ | return ret(c, 0, x); | |
| 829 | } | ||
| 830 | |||
| 831 | } | ||
| 832 |