float compute_max_saturation(float a, float b)
{
float k0, k1, k2, k3, k4, wl, wm, ws;
if (-1.88170328f * a - 0.80936493f * b > 1)
{
k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f;
wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f;
}
else if (1.81444104f * a - 1.19445276f * b > 1)
{
k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f;
wl = -1.2681437731f; wm = +2.6097574011f; ws = -0.3413193965f;
}
else
{
k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f;
wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f;
}
float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;
float k_l = +0.3963377774f * a + 0.2158037573f * b;
float k_m = -0.1055613458f * a - 0.0638541728f * b;
float k_s = -0.0894841775f * a - 1.2914855480f * b;
{
float l_ = 1.f + S * k_l;
float m_ = 1.f + S * k_m;
float s_ = 1.f + S * k_s;
float l = l_ * l_ * l_;
float m = m_ * m_ * m_;
float s = s_ * s_ * s_;
float l_dS = 3.f * k_l * l_ * l_;
float m_dS = 3.f * k_m * m_ * m_;
float s_dS = 3.f * k_s * s_ * s_;
float l_dS2 = 6.f * k_l * k_l * l_;
float m_dS2 = 6.f * k_m * k_m * m_;
float s_dS2 = 6.f * k_s * k_s * s_;
float f = wl * l + wm * m + ws * s;
float f1 = wl * l_dS + wm * m_dS + ws * s_dS;
float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2;
S = S - f * f1 / (f1*f1 - 0.5f * f * f2);
}
return S;
}
struct LC { float L; float C; };
LC find_cusp(float a, float b)
{
float S_cusp = compute_max_saturation(a, b);
RGB rgb_at_max = oklab_to_linear_srgb({ 1, S_cusp * a, S_cusp * b });
float L_cusp = cbrtf(1.f / max(max(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b));
float C_cusp = L_cusp * S_cusp;
return { L_cusp , C_cusp };
}
float find_gamut_intersection(float a, float b, float L1, float C1, float L0)
{
LC cusp = find_cusp(a, b);
float t;
if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0.f)
{
t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1));
}
else
{
t = cusp.C * (L0 - 1.f) / (C1 * (cusp.L - 1.f) + cusp.C * (L0 - L1));
{
float dL = L1 - L0;
float dC = C1;
float k_l = +0.3963377774f * a + 0.2158037573f * b;
float k_m = -0.1055613458f * a - 0.0638541728f * b;
float k_s = -0.0894841775f * a - 1.2914855480f * b;
float l_dt = dL + dC * k_l;
float m_dt = dL + dC * k_m;
float s_dt = dL + dC * k_s;
{
float L = L0 * (1.f - t) + t * L1;
float C = t * C1;
float l_ = L + C * k_l;
float m_ = L + C * k_m;
float s_ = L + C * k_s;
float l = l_ * l_ * l_;
float m = m_ * m_ * m_;
float s = s_ * s_ * s_;
float ldt = 3 * l_dt * l_ * l_;
float mdt = 3 * m_dt * m_ * m_;
float sdt = 3 * s_dt * s_ * s_;
float ldt2 = 6 * l_dt * l_dt * l_;
float mdt2 = 6 * m_dt * m_dt * m_;
float sdt2 = 6 * s_dt * s_dt * s_;
float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1;
float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt;
float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2;
float u_r = r1 / (r1 * r1 - 0.5f * r * r2);
float t_r = -r * u_r;
float g = -1.2681437731f * l + 2.6097574011f * m - 0.3413193965f * s - 1;
float g1 = -1.2681437731f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt;
float g2 = -1.2681437731f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2;
float u_g = g1 / (g1 * g1 - 0.5f * g * g2);
float t_g = -g * u_g;
float b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1;
float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt;
float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2;
float u_b = b1 / (b1 * b1 - 0.5f * b * b2);
float t_b = -b * u_b;
t_r = u_r >= 0.f ? t_r : FLT_MAX;
t_g = u_g >= 0.f ? t_g : FLT_MAX;
t_b = u_b >= 0.f ? t_b : FLT_MAX;
t += min(t_r, min(t_g, t_b));
}
}
}
return t;
}
float clamp(float x, float min, float max)
{
if (x < min)
return min;
if (x > max)
return max;
return x;
}
float sgn(float x)
{
return (float)(0.f < x) - (float)(x < 0.f);
}
RGB gamut_clip_preserve_chroma(RGB rgb)
{
if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
return rgb;
Lab lab = linear_srgb_to_oklab(rgb);
float L = lab.L;
float eps = 0.00001f;
float C = max(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
float a_ = lab.a / C;
float b_ = lab.b / C;
float L0 = clamp(L, 0, 1);
float t = find_gamut_intersection(a_, b_, L, C, L0);
float L_clipped = L0 * (1 - t) + t * L;
float C_clipped = t * C;
return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
}
RGB gamut_clip_project_to_0_5(RGB rgb)
{
if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
return rgb;
Lab lab = linear_srgb_to_oklab(rgb);
float L = lab.L;
float eps = 0.00001f;
float C = max(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
float a_ = lab.a / C;
float b_ = lab.b / C;
float L0 = 0.5;
float t = find_gamut_intersection(a_, b_, L, C, L0);
float L_clipped = L0 * (1 - t) + t * L;
float C_clipped = t * C;
return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
}
RGB gamut_clip_project_to_L_cusp(RGB rgb)
{
if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
return rgb;
Lab lab = linear_srgb_to_oklab(rgb);
float L = lab.L;
float eps = 0.00001f;
float C = max(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
float a_ = lab.a / C;
float b_ = lab.b / C;
LC cusp = find_cusp(a_, b_);
float L0 = cusp.L;
float t = find_gamut_intersection(a_, b_, L, C, L0);
float L_clipped = L0 * (1 - t) + t * L;
float C_clipped = t * C;
return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
}
RGB gamut_clip_adaptive_L0_0_5(RGB rgb, float alpha = 0.05f)
{
if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
return rgb;
Lab lab = linear_srgb_to_oklab(rgb);
float L = lab.L;
float eps = 0.00001f;
float C = max(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
float a_ = lab.a / C;
float b_ = lab.b / C;
float Ld = L - 0.5f;
float e1 = 0.5f + fabs(Ld) + alpha * C;
float L0 = 0.5f*(1.f + sgn(Ld)*(e1 - sqrtf(e1*e1 - 2.f *fabs(Ld))));
float t = find_gamut_intersection(a_, b_, L, C, L0);
float L_clipped = L0 * (1.f - t) + t * L;
float C_clipped = t * C;
return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
}
RGB gamut_clip_adaptive_L0_L_cusp(RGB rgb, float alpha = 0.05f)
{
if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
return rgb;
Lab lab = linear_srgb_to_oklab(rgb);
float L = lab.L;
float eps = 0.00001f;
float C = max(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
float a_ = lab.a / C;
float b_ = lab.b / C;
LC cusp = find_cusp(a_, b_);
float Ld = L - cusp.L;
float k = 2.f * (Ld > 0 ? 1.f - cusp.L : cusp.L);
float e1 = 0.5f*k + fabs(Ld) + alpha * C/k;
float L0 = cusp.L + 0.5f * (sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * k * fabs(Ld))));
float t = find_gamut_intersection(a_, b_, L, C, L0);
float L_clipped = L0 * (1.f - t) + t * L;
float C_clipped = t * C;
return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
}