From f804c9e3465ffa4e5469d7723b2f9ee401099221 Mon Sep 17 00:00:00 2001 From: Romain Guy Date: Thu, 5 Aug 2021 12:44:49 -0700 Subject: [PATCH] Switch white balance to CAT16 and add gamut mapping code (#4457) --- filament/CMakeLists.txt | 2 + filament/src/ColorGrading.cpp | 21 ++- filament/src/ColorSpace.cpp | 266 +++++++++++++++++++++++++++ filament/src/ColorSpace.h | 74 ++++++-- libs/math/include/math/TVecHelpers.h | 7 + libs/viewer/src/SimpleViewer.cpp | 6 +- 6 files changed, 348 insertions(+), 28 deletions(-) create mode 100644 filament/src/ColorSpace.cpp diff --git a/filament/CMakeLists.txt b/filament/CMakeLists.txt index 78480e7a70..cdc7fad59f 100644 --- a/filament/CMakeLists.txt +++ b/filament/CMakeLists.txt @@ -51,6 +51,7 @@ set(SRCS src/Camera.cpp src/Color.cpp src/ColorGrading.cpp + src/ColorSpace.cpp src/Culler.cpp src/DFG.cpp src/DebugRegistry.cpp @@ -100,6 +101,7 @@ set(SRCS ) set(PRIVATE_HDRS + src/ColorSpace.h src/FilamentAPI-impl.h src/FrameHistory.h src/FrameInfo.h diff --git a/filament/src/ColorGrading.cpp b/filament/src/ColorGrading.cpp index 1569e9a967..8fc5732bae 100644 --- a/filament/src/ColorGrading.cpp +++ b/filament/src/ColorGrading.cpp @@ -265,7 +265,7 @@ inline float3 adjustExposure(float3 v, float exposure) { // Return the chromatic adaptation coefficients in LMS space for the given // temperature/tint offsets. The chromatic adaption is perfomed following -// the von Kries method, using the CIECAT02 transform. +// the von Kries method, using the CIECAT16 transform. // See https://en.wikipedia.org/wiki/Chromatic_adaptation // See https://en.wikipedia.org/wiki/CIECAM02#Chromatic_adaptation UTILS_ALWAYS_INLINE @@ -277,13 +277,13 @@ inline float3 adaptationTransform(float2 whiteBalance) { float x = ILLUMINANT_D65_xyY.x - k * (k < 0.0f ? 0.0214f : 0.066f); float y = chromaticityCoordinateIlluminantD(x) + t * 0.066f; - float3 lms = XYZ_to_CIECAT02 * xyY_to_XYZ({x, y, 1.0f}); - return ILLUMINANT_D65_LMS / lms; + float3 lms = XYZ_to_CIECAT16 * xyY_to_XYZ({x, y, 1.0f}); + return ILLUMINANT_D65_LMS_CAT16 / lms; } UTILS_ALWAYS_INLINE inline float3 chromaticAdaptation(float3 v, float3 adaptationTransform) { - return LMS_to_sRGB * (adaptationTransform * (sRGB_to_LMS * v)); + return LMS_CAT16_to_sRGB * (adaptationTransform * (sRGB_to_LMS_CAT16 * v)); } //------------------------------------------------------------------------------ @@ -381,7 +381,7 @@ static float3 luminanceScaling(float3 x, float luminanceIn = dot(x, luminanceWeights); // TODO: We could optimize for the case of single-channel luminance - float luminanceOut = toneMapper(luminanceIn).x; + float luminanceOut = toneMapper(luminanceIn).y; float peak = max(x); float3 chromaRatio = max(x / peak, 0.0f); @@ -512,6 +512,7 @@ FColorGrading::FColorGrading(FEngine& engine, const Builder& builder) { v = LogC_to_linear(v); // TODO: Performed in sRGB, should be in Rec.2020 or AP1 + if (builder->hasAdjustments) { // Exposure v = adjustExposure(v, builder->exposure); @@ -570,7 +571,9 @@ FColorGrading::FColorGrading(FEngine& engine, const Builder& builder) { // TODO: We should convert to the output color space if we use a working // color space that's not sRGB // TODO: Allow the user to customize the output color space + // TODO: Gamut mapping + // We need to clamp for the output transfer function v = saturate(v); // Apply OECF @@ -591,10 +594,10 @@ FColorGrading::FColorGrading(FEngine& engine, const Builder& builder) { #pragma clang loop vectorize_width(8) for (size_t i = 0; i < count; ++i) { float4 v{src[i]}; - uint32_t r = uint32_t(std::floor(v.x * 1023.0f + 0.5f)); - uint32_t g = uint32_t(std::floor(v.y * 1023.0f + 0.5f)); - uint32_t b = uint32_t(std::floor(v.z * 1023.0f + 0.5f)); - dst[i] = (b << 20u) | (g << 10u) | r; + uint32_t pr = uint32_t(std::floor(v.x * 1023.0f + 0.5f)); + uint32_t pg = uint32_t(std::floor(v.y * 1023.0f + 0.5f)); + uint32_t pb = uint32_t(std::floor(v.z * 1023.0f + 0.5f)); + dst[i] = (pb << 20u) | (pg << 10u) | pr; } } diff --git a/filament/src/ColorSpace.cpp b/filament/src/ColorSpace.cpp new file mode 100644 index 0000000000..93d31c662c --- /dev/null +++ b/filament/src/ColorSpace.cpp @@ -0,0 +1,266 @@ +/* + * Copyright (C) 2021 The Android Open Source Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "ColorSpace.h" + +namespace filament { + +using namespace math; + +/* + * The section below is subject to the following license. + * Source: https://bottosson.github.io/posts/gamutclipping/ + * + * Copyright (c) 2021 Björn Ottosson + * + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies + * of the Software, and to permit persons to whom the Software is furnished to do + * so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +// Finds the maximum saturation possible for a given hue that fits in sRGB +// Saturation here is defined as S = C/L +// a and b must be normalized so a^2 + b^2 == 1 +float compute_max_saturation(float a, float b) noexcept { + // Max saturation will be when one of r, g or b goes below zero. + + // Select different coefficients depending on which component goes below zero first + float k0, k1, k2, k3, k4, wl, wm, ws; + + if (-1.88170328f * a - 0.80936493f * b > 1) { + // Red component + 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) { + // Green component + 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 { + // Blue component + k0 = +1.35733652f; + k1 = -0.00915799f; + k2 = -1.15130210f; + k3 = -0.50559606f; + k4 = +0.00692167f; + wl = -0.0041960863f; + wm = -0.7034186147f; + ws = +1.7076147010f; + } + + // Approximate max saturation using a polynomial: + float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b; + + // Do one step Halley's method to get closer + // this gives an error less than 10e6, except for some blue hues where the dS/dh + // is close to infinite + // this should be sufficient for most applications, otherwise do two/three steps + + 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; +} + +// finds L_cusp and C_cusp for a given hue +// a and b must be normalized so a^2 + b^2 == 1 +float2 find_cusp(float a, float b) noexcept { + // First, find the maximum saturation (saturation S = C/L) + float S_cusp = compute_max_saturation(a, b); + + // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1: + float3 rgb_at_max = OkLab_to_sRGB({1.0f, S_cusp * a, S_cusp * b}); + float L_cusp = std::cbrt(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}; +} + +// Finds intersection of the line defined by +// L = L0 * (1 - t) + t * L1; +// C = t * C1; +// a and b must be normalized so a^2 + b^2 == 1 +float find_gamut_intersection(float a, float b, float L1, float C1, float L0) noexcept { + // Find the cusp of the gamut triangle + float2 cusp = find_cusp(a, b); + + // Find the intersection for upper and lower half seprately + float t; + if (((L1 - L0) * cusp.y - (cusp.x - L0) * C1) <= 0.f) { + // Lower half + t = cusp.y * L0 / (C1 * cusp.x + cusp.y * (L0 - L1)); + } else { + // Upper half + + // First intersect with triangle + t = cusp.y * (L0 - 1.f) / (C1 * (cusp.x - 1.f) + cusp.y * (L0 - L1)); + + // Then one step Halley's method + { + 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; + + + // If higher accuracy is required, 2 or 3 iterations of the + // following block can be used: + { + 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 b0 = -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 * b0 * b2); + float t_b = -b0 * u_b; + + t_r = u_r >= 0.f ? t_r : std::numeric_limits::max(); + t_g = u_g >= 0.f ? t_g : std::numeric_limits::max(); + t_b = u_b >= 0.f ? t_b : std::numeric_limits::max(); + + t += min(t_r, min(t_g, t_b)); + } + } + } + + return t; +} + +constexpr float sgn(float x) noexcept { + return (float) (0.f < x) - (float) (x < 0.f); +} + +inline float3 gamut_clip_adaptive_L0_0_5(float3 rgb, float alpha = 0.05f) noexcept { + if (all(lessThanEqual(rgb, float3{1.0f})) && all(greaterThanEqual(rgb, float3{0.0f}))) { + return rgb; + } + + float3 lab = sRGB_to_OkLab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, std::sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float Ld = L - 0.5f; + float e1 = 0.5f + std::abs(Ld) + alpha * C; + float L0 = 0.5f * (1.f + sgn(Ld) * (e1 - std::sqrt(e1 * e1 - 2.f * std::abs(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_sRGB({L_clipped, C_clipped * a_, C_clipped * b_}); +} + +float3 gamutMapping_sRGB(float3 rgb) noexcept { + return gamut_clip_adaptive_L0_0_5(rgb, 0.5f); +} + +/* + * End source: https://bottosson.github.io/posts/gamutclipping/ + */ + +} //namespace filament diff --git a/filament/src/ColorSpace.h b/filament/src/ColorSpace.h index 976a7ed4ae..0d664deaa5 100644 --- a/filament/src/ColorSpace.h +++ b/filament/src/ColorSpace.h @@ -64,16 +64,16 @@ constexpr mat3f XYZ_to_REC2020{ -0.2533663f, 0.0157685f, 0.9421031f }; -constexpr mat3f XYZ_to_CIECAT02{ - 0.7328000f, -0.7036000f, 0.0030000f, - 0.4296000f, 1.6975000f, 0.0136000f, - -0.1624000f, 0.0061000f, 0.9834000f +constexpr mat3f XYZ_to_CIECAT16{ + 0.401288f, -0.250268f, -0.002079f, + 0.650173f, 1.204414f, 0.048952f, + -0.051461f, 0.045854f, 0.953127f }; -constexpr mat3f CIECAT02_to_XYZ{ - 1.0961200f, 0.4543690f, -0.0096276f, - -0.2788690f, 0.4735330f, -0.0056980f, - 0.1827450f, 0.0720978f, 1.0153300f +constexpr mat3f CIECAT16_to_XYZ{ + 1.862068f, 0.387527f, -0.015841f, + -1.011255f, 0.621447f, -0.034123f, + 0.149187f, -0.008974f, 1.049964f }; constexpr mat3f AP1_to_XYZ{ @@ -119,18 +119,42 @@ constexpr mat3f AP0_to_sRGB{ }; constexpr mat3f sRGB_to_AP0{ - 0.4397010f, 0.0897923f, 0.0175440f, - 0.3829780f, 0.8134230f, 0.1115440f, - 0.1773350f, 0.0967616f, 0.8707040f + 0.4397010f, 0.0897923f, 0.0175440f, + 0.3829780f, 0.8134230f, 0.1115440f, + 0.1773350f, 0.0967616f, 0.8707040f +}; + +constexpr mat3f sRGB_to_OkLab_LMS{ + 0.4122214708f, 0.2119034982f, 0.0883024619f, + 0.5363325363f, 0.6806995451f, 0.2817188376f, + 0.0514459929f, 0.1073969566f, 0.6299787005f +}; + +constexpr mat3f OkLab_LMS_to_OkLab{ + 0.2104542553f, 1.9779984951f, 0.0259040371f, + 0.7936177850f, -2.4285922050f, 0.7827717662f, + -0.0040720468f, 0.4505937099f, -0.8086757660f +}; + +constexpr mat3f OkLab_to_OkLab_LMS{ + 1.0000000000f, 1.0000000000f, 1.0000000000f, + 0.3963377774f, -0.1055613458f, -0.0894841775f, + 0.2158037573f, -0.0638541728f, -1.2914855480f +}; + +constexpr mat3f OkLab_LMS_to_sRGB{ + 4.0767416621f, -1.2684380046f, -0.0041960863f, + -3.3077115913f, 2.6097574011f, -0.7034186147f, + 0.2309699292f, -0.3413193965f, 1.7076147010f }; constexpr mat3f sRGB_to_REC2020 = XYZ_to_REC2020 * sRGB_to_XYZ; constexpr mat3f REC2020_to_sRGB = XYZ_to_sRGB * REC2020_to_XYZ; -constexpr mat3f sRGB_to_LMS = XYZ_to_CIECAT02 * sRGB_to_XYZ; +constexpr mat3f sRGB_to_LMS_CAT16 = XYZ_to_CIECAT16 * sRGB_to_XYZ; -constexpr mat3f LMS_to_sRGB = XYZ_to_sRGB * CIECAT02_to_XYZ; +constexpr mat3f LMS_CAT16_to_sRGB = XYZ_to_sRGB * CIECAT16_to_XYZ; constexpr mat3f REC2020_to_AP0 = AP1_to_AP0 * XYZ_to_AP1 * REC2020_to_XYZ; @@ -143,9 +167,9 @@ constexpr mat3f AP1_to_REC2020 = XYZ_to_REC2020 * AP1_to_XYZ; // Standard CIE 1931 2° illuminant D65, in xyY space constexpr float3 ILLUMINANT_D65_xyY{0.31271f, 0.32902f, 1.0f}; -// Standard CIE 1931 2° illuminant D65, in LMS space (CIECAT02) -// Result of: XYZ_to_CIECAT02 * xyY_to_XYZ(ILLUMINANT_D65_xyY); -constexpr float3 ILLUMINANT_D65_LMS{0.949237f, 1.03542f, 1.08728f}; +// Standard CIE 1931 2° illuminant D65, in LMS space (CIECAT16) +// Result of: XYZ_to_CIECAT16 * xyY_to_XYZ(ILLUMINANT_D65_xyY); +constexpr float3 ILLUMINANT_D65_LMS_CAT16{0.975533f, 1.016483f, 1.084837f}; // RGB to luminance coefficients for Rec.2020, from REC2020_to_XYZ constexpr float3 LUMA_REC2020{0.2627002f, 0.6779981f, 0.0593017f}; @@ -187,6 +211,18 @@ inline constexpr xyY XYZ_to_xyY(XYZ v) noexcept { return float3(v.xy / max(v.x + v.y + v.z, 1e-5f), v.y); } +inline float3 sRGB_to_OkLab(float3 x) noexcept { + return OkLab_LMS_to_OkLab * cbrt(sRGB_to_OkLab_LMS * x); +} + +inline constexpr float3 pow3(float3 x) noexcept { + return x * x * x; +} + +inline float3 OkLab_to_sRGB(float3 x) noexcept { + return OkLab_LMS_to_sRGB * pow3(OkLab_to_OkLab_LMS * x); +} + //------------------------------------------------------------------------------ // Conversion functions and encoding/decoding //------------------------------------------------------------------------------ @@ -284,6 +320,12 @@ inline float3 EOCF_PQ(float3 x, float maxPqValue) { return linearColor * maxPqValue; } +//------------------------------------------------------------------------------ +// Gamut mapping +//------------------------------------------------------------------------------ + +float3 gamutMapping_sRGB(float3 rgb) noexcept; + } // namespace filament #endif //TNT_FILAMENT_COLOR_SPACE_H diff --git a/libs/math/include/math/TVecHelpers.h b/libs/math/include/math/TVecHelpers.h index c7b09c3d30..ce66ec1ad3 100644 --- a/libs/math/include/math/TVecHelpers.h +++ b/libs/math/include/math/TVecHelpers.h @@ -463,6 +463,13 @@ private: return v; } + friend inline VECTOR MATH_PURE cbrt(VECTOR v) { + for (size_t i = 0; i < v.size(); i++) { + v[i] = std::cbrt(v[i]); + } + return v; + } + friend inline VECTOR MATH_PURE exp(VECTOR v) { for (size_t i = 0; i < v.size(); i++) { v[i] = std::exp(v[i]); diff --git a/libs/viewer/src/SimpleViewer.cpp b/libs/viewer/src/SimpleViewer.cpp index d67f4fd8c3..18112bdb00 100644 --- a/libs/viewer/src/SimpleViewer.cpp +++ b/libs/viewer/src/SimpleViewer.cpp @@ -379,9 +379,9 @@ SimpleViewer::SimpleViewer(filament::Engine* engine, filament::Scene* scene, fil .intensity(mSettings.lighting.sunlightIntensity) .direction(normalize(mSettings.lighting.sunlightDirection)) .castShadows(true) - .sunAngularRadius(1.9) - .sunHaloSize(10.0) - .sunHaloFalloff(80.0) + .sunAngularRadius(1.0f) + .sunHaloSize(2.0f) + .sunHaloFalloff(80.0f) .build(*engine, mSunlight); if (mSettings.lighting.enableSunlight) { mScene->addEntity(mSunlight);