diff options
author | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-19 00:47:55 +0000 |
---|---|---|
committer | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-19 00:47:55 +0000 |
commit | 26a029d407be480d791972afb5975cf62c9360a6 (patch) | |
tree | f435a8308119effd964b339f76abb83a57c29483 /third_party/jpeg-xl/lib/jxl/cms | |
parent | Initial commit. (diff) | |
download | firefox-26a029d407be480d791972afb5975cf62c9360a6.tar.xz firefox-26a029d407be480d791972afb5975cf62c9360a6.zip |
Adding upstream version 124.0.1.upstream/124.0.1
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'third_party/jpeg-xl/lib/jxl/cms')
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/color_encoding_cms.h | 623 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/jxl_cms.cc | 1343 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/jxl_cms_internal.h | 1083 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/opsin_params.h | 160 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/tone_mapping-inl.h | 191 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/tone_mapping.h | 179 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/tone_mapping_test.cc | 147 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/transfer_functions-inl.h | 334 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/transfer_functions.h | 131 | ||||
-rw-r--r-- | third_party/jpeg-xl/lib/jxl/cms/transfer_functions_test.cc | 94 |
10 files changed, 4285 insertions, 0 deletions
diff --git a/third_party/jpeg-xl/lib/jxl/cms/color_encoding_cms.h b/third_party/jpeg-xl/lib/jxl/cms/color_encoding_cms.h new file mode 100644 index 0000000000..db61f820ca --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/color_encoding_cms.h @@ -0,0 +1,623 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#ifndef LIB_JXL_CMS_COLOR_ENCODING_CMS_H_ +#define LIB_JXL_CMS_COLOR_ENCODING_CMS_H_ + +#include <jxl/cms_interface.h> +#include <jxl/color_encoding.h> +#include <jxl/types.h> + +#include <cmath> +#include <cstdint> +#include <cstring> +#include <utility> +#include <vector> + +#include "lib/jxl/base/status.h" + +namespace jxl { +namespace cms { + +using IccBytes = std::vector<uint8_t>; + +// Returns whether the two inputs are approximately equal. +static inline bool ApproxEq(const double a, const double b, + double max_l1 = 1E-3) { + // Threshold should be sufficient for ICC's 15-bit fixed-point numbers. + // We have seen differences of 7.1E-5 with lcms2 and 1E-3 with skcms. + return std::abs(a - b) <= max_l1; +} + +// (All CIE units are for the standard 1931 2 degree observer) + +// Color space the color pixel data is encoded in. The color pixel data is +// 3-channel in all cases except in case of kGray, where it uses only 1 channel. +// This also determines the amount of channels used in modular encoding. +enum class ColorSpace : uint32_t { + // Trichromatic color data. This also includes CMYK if a kBlack + // ExtraChannelInfo is present. This implies, if there is an ICC profile, that + // the ICC profile uses a 3-channel color space if no kBlack extra channel is + // present, or uses color space 'CMYK' if a kBlack extra channel is present. + kRGB, + // Single-channel data. This implies, if there is an ICC profile, that the ICC + // profile also represents single-channel data and has the appropriate color + // space ('GRAY'). + kGray, + // Like kRGB, but implies fixed values for primaries etc. + kXYB, + // For non-RGB/gray data, e.g. from non-electro-optical sensors. Otherwise + // the same conditions as kRGB apply. + kUnknown + // NB: don't forget to update EnumBits! +}; + +// Values from CICP ColourPrimaries. +enum class WhitePoint : uint32_t { + kD65 = 1, // sRGB/BT.709/Display P3/BT.2020 + kCustom = 2, // Actual values encoded in separate fields + kE = 10, // XYZ + kDCI = 11, // DCI-P3 + // NB: don't forget to update EnumBits! +}; + +// Values from CICP ColourPrimaries +enum class Primaries : uint32_t { + kSRGB = 1, // Same as BT.709 + kCustom = 2, // Actual values encoded in separate fields + k2100 = 9, // Same as BT.2020 + kP3 = 11, + // NB: don't forget to update EnumBits! +}; + +// Values from CICP TransferCharacteristics +enum class TransferFunction : uint32_t { + k709 = 1, + kUnknown = 2, + kLinear = 8, + kSRGB = 13, + kPQ = 16, // from BT.2100 + kDCI = 17, // from SMPTE RP 431-2 reference projector + kHLG = 18, // from BT.2100 + // NB: don't forget to update EnumBits! +}; + +enum class RenderingIntent : uint32_t { + // Values match ICC sRGB encodings. + kPerceptual = 0, // good for photos, requires a profile with LUT. + kRelative, // good for logos. + kSaturation, // perhaps useful for CG with fully saturated colors. + kAbsolute, // leaves white point unchanged; good for proofing. + // NB: don't forget to update EnumBits! +}; + +// Chromaticity (Y is omitted because it is 1 for white points and implicit for +// primaries) +struct CIExy { + double x = 0.0; + double y = 0.0; +}; + +struct PrimariesCIExy { + CIExy r; + CIExy g; + CIExy b; +}; + +// Serializable form of CIExy. +struct Customxy { + static constexpr uint32_t kMul = 1000000; + static constexpr double kRoughLimit = 4.0; + static constexpr int32_t kMin = -0x200000; + static constexpr int32_t kMax = 0x1FFFFF; + + int32_t x = 0; + int32_t y = 0; + + CIExy GetValue() const { + CIExy xy; + xy.x = x * (1.0 / kMul); + xy.y = y * (1.0 / kMul); + return xy; + } + + Status SetValue(const CIExy& xy) { + bool ok = (std::abs(xy.x) < kRoughLimit) && (std::abs(xy.y) < kRoughLimit); + if (!ok) return JXL_FAILURE("X or Y is out of bounds"); + x = static_cast<int32_t>(roundf(xy.x * kMul)); + if (x < kMin || x > kMax) return JXL_FAILURE("X is out of bounds"); + y = static_cast<int32_t>(roundf(xy.y * kMul)); + if (y < kMin || y > kMax) return JXL_FAILURE("Y is out of bounds"); + return true; + } + + bool IsSame(const Customxy& other) const { + return (x == other.x) && (y == other.y); + } +}; + +static inline Status WhitePointFromExternal(const JxlWhitePoint external, + WhitePoint* out) { + switch (external) { + case JXL_WHITE_POINT_D65: + *out = WhitePoint::kD65; + return true; + case JXL_WHITE_POINT_CUSTOM: + *out = WhitePoint::kCustom; + return true; + case JXL_WHITE_POINT_E: + *out = WhitePoint::kE; + return true; + case JXL_WHITE_POINT_DCI: + *out = WhitePoint::kDCI; + return true; + } + return JXL_FAILURE("Invalid WhitePoint enum value %d", + static_cast<int>(external)); +} + +static inline Status PrimariesFromExternal(const JxlPrimaries external, + Primaries* out) { + switch (external) { + case JXL_PRIMARIES_SRGB: + *out = Primaries::kSRGB; + return true; + case JXL_PRIMARIES_CUSTOM: + *out = Primaries::kCustom; + return true; + case JXL_PRIMARIES_2100: + *out = Primaries::k2100; + return true; + case JXL_PRIMARIES_P3: + *out = Primaries::kP3; + return true; + } + return JXL_FAILURE("Invalid Primaries enum value"); +} + +static inline Status RenderingIntentFromExternal( + const JxlRenderingIntent external, RenderingIntent* out) { + switch (external) { + case JXL_RENDERING_INTENT_PERCEPTUAL: + *out = RenderingIntent::kPerceptual; + return true; + case JXL_RENDERING_INTENT_RELATIVE: + *out = RenderingIntent::kRelative; + return true; + case JXL_RENDERING_INTENT_SATURATION: + *out = RenderingIntent::kSaturation; + return true; + case JXL_RENDERING_INTENT_ABSOLUTE: + *out = RenderingIntent::kAbsolute; + return true; + } + return JXL_FAILURE("Invalid RenderingIntent enum value"); +} + +struct CustomTransferFunction { + // Highest reasonable value for the gamma of a transfer curve. + static constexpr uint32_t kMaxGamma = 8192; + static constexpr uint32_t kGammaMul = 10000000; + + bool have_gamma = false; + + // OETF exponent to go from linear to gamma-compressed. + uint32_t gamma = 0; // Only used if have_gamma_. + + // Can be kUnknown. + TransferFunction transfer_function = + TransferFunction::kSRGB; // Only used if !have_gamma_. + + TransferFunction GetTransferFunction() const { + JXL_ASSERT(!have_gamma); + return transfer_function; + } + void SetTransferFunction(const TransferFunction tf) { + have_gamma = false; + transfer_function = tf; + } + + bool IsUnknown() const { + return !have_gamma && (transfer_function == TransferFunction::kUnknown); + } + bool IsSRGB() const { + return !have_gamma && (transfer_function == TransferFunction::kSRGB); + } + bool IsLinear() const { + return !have_gamma && (transfer_function == TransferFunction::kLinear); + } + bool IsPQ() const { + return !have_gamma && (transfer_function == TransferFunction::kPQ); + } + bool IsHLG() const { + return !have_gamma && (transfer_function == TransferFunction::kHLG); + } + bool Is709() const { + return !have_gamma && (transfer_function == TransferFunction::k709); + } + bool IsDCI() const { + return !have_gamma && (transfer_function == TransferFunction::kDCI); + } + + double GetGamma() const { + JXL_ASSERT(have_gamma); + return gamma * (1.0 / kGammaMul); // (0, 1) + } + Status SetGamma(double new_gamma) { + if (new_gamma < (1.0 / kMaxGamma) || new_gamma > 1.0) { + return JXL_FAILURE("Invalid gamma %f", new_gamma); + } + + have_gamma = false; + if (ApproxEq(new_gamma, 1.0)) { + transfer_function = TransferFunction::kLinear; + return true; + } + if (ApproxEq(new_gamma, 1.0 / 2.6)) { + transfer_function = TransferFunction::kDCI; + return true; + } + // Don't translate 0.45.. to kSRGB nor k709 - that might change pixel + // values because those curves also have a linear part. + + have_gamma = true; + gamma = roundf(new_gamma * kGammaMul); + transfer_function = TransferFunction::kUnknown; + return true; + } + + bool IsSame(const CustomTransferFunction& other) const { + if (have_gamma != other.have_gamma) { + return false; + } + if (have_gamma) { + if (gamma != other.gamma) { + return false; + } + } else { + if (transfer_function != other.transfer_function) { + return false; + } + } + return true; + } +}; + +static inline Status ConvertExternalToInternalTransferFunction( + const JxlTransferFunction external, TransferFunction* internal) { + switch (external) { + case JXL_TRANSFER_FUNCTION_709: + *internal = TransferFunction::k709; + return true; + case JXL_TRANSFER_FUNCTION_UNKNOWN: + *internal = TransferFunction::kUnknown; + return true; + case JXL_TRANSFER_FUNCTION_LINEAR: + *internal = TransferFunction::kLinear; + return true; + case JXL_TRANSFER_FUNCTION_SRGB: + *internal = TransferFunction::kSRGB; + return true; + case JXL_TRANSFER_FUNCTION_PQ: + *internal = TransferFunction::kPQ; + return true; + case JXL_TRANSFER_FUNCTION_DCI: + *internal = TransferFunction::kDCI; + return true; + case JXL_TRANSFER_FUNCTION_HLG: + *internal = TransferFunction::kHLG; + return true; + case JXL_TRANSFER_FUNCTION_GAMMA: + return JXL_FAILURE("Gamma should be handled separately"); + } + return JXL_FAILURE("Invalid TransferFunction enum value"); +} + +// Compact encoding of data required to interpret and translate pixels to a +// known color space. Stored in Metadata. Thread-compatible. +struct ColorEncoding { + // Only valid if HaveFields() + WhitePoint white_point = WhitePoint::kD65; + Primaries primaries = Primaries::kSRGB; // Only valid if HasPrimaries() + RenderingIntent rendering_intent = RenderingIntent::kRelative; + + // When false, fields such as white_point and tf are invalid and must not be + // used. This occurs after setting a raw bytes-only ICC profile, only the + // ICC bytes may be used. The color_space_ field is still valid. + bool have_fields = true; + + IccBytes icc; // Valid ICC profile + + ColorSpace color_space = ColorSpace::kRGB; // Can be kUnknown + bool cmyk = false; + + // "late sync" fields + CustomTransferFunction tf; + Customxy white; // Only used if white_point == kCustom + Customxy red; // Only used if primaries == kCustom + Customxy green; // Only used if primaries == kCustom + Customxy blue; // Only used if primaries == kCustom + + // Returns false if the field is invalid and unusable. + bool HasPrimaries() const { + return (color_space != ColorSpace::kGray) && + (color_space != ColorSpace::kXYB); + } + + size_t Channels() const { return (color_space == ColorSpace::kGray) ? 1 : 3; } + + PrimariesCIExy GetPrimaries() const { + JXL_DASSERT(have_fields); + JXL_ASSERT(HasPrimaries()); + PrimariesCIExy xy; + switch (primaries) { + case Primaries::kCustom: + xy.r = red.GetValue(); + xy.g = green.GetValue(); + xy.b = blue.GetValue(); + return xy; + + case Primaries::kSRGB: + xy.r.x = 0.639998686; + xy.r.y = 0.330010138; + xy.g.x = 0.300003784; + xy.g.y = 0.600003357; + xy.b.x = 0.150002046; + xy.b.y = 0.059997204; + return xy; + + case Primaries::k2100: + xy.r.x = 0.708; + xy.r.y = 0.292; + xy.g.x = 0.170; + xy.g.y = 0.797; + xy.b.x = 0.131; + xy.b.y = 0.046; + return xy; + + case Primaries::kP3: + xy.r.x = 0.680; + xy.r.y = 0.320; + xy.g.x = 0.265; + xy.g.y = 0.690; + xy.b.x = 0.150; + xy.b.y = 0.060; + return xy; + } + JXL_UNREACHABLE("Invalid Primaries %u", static_cast<uint32_t>(primaries)); + } + + Status SetPrimaries(const PrimariesCIExy& xy) { + JXL_DASSERT(have_fields); + JXL_ASSERT(HasPrimaries()); + if (xy.r.x == 0.0 || xy.r.y == 0.0 || xy.g.x == 0.0 || xy.g.y == 0.0 || + xy.b.x == 0.0 || xy.b.y == 0.0) { + return JXL_FAILURE("Invalid primaries %f %f %f %f %f %f", xy.r.x, xy.r.y, + xy.g.x, xy.g.y, xy.b.x, xy.b.y); + } + + if (ApproxEq(xy.r.x, 0.64) && ApproxEq(xy.r.y, 0.33) && + ApproxEq(xy.g.x, 0.30) && ApproxEq(xy.g.y, 0.60) && + ApproxEq(xy.b.x, 0.15) && ApproxEq(xy.b.y, 0.06)) { + primaries = Primaries::kSRGB; + return true; + } + + if (ApproxEq(xy.r.x, 0.708) && ApproxEq(xy.r.y, 0.292) && + ApproxEq(xy.g.x, 0.170) && ApproxEq(xy.g.y, 0.797) && + ApproxEq(xy.b.x, 0.131) && ApproxEq(xy.b.y, 0.046)) { + primaries = Primaries::k2100; + return true; + } + if (ApproxEq(xy.r.x, 0.680) && ApproxEq(xy.r.y, 0.320) && + ApproxEq(xy.g.x, 0.265) && ApproxEq(xy.g.y, 0.690) && + ApproxEq(xy.b.x, 0.150) && ApproxEq(xy.b.y, 0.060)) { + primaries = Primaries::kP3; + return true; + } + + primaries = Primaries::kCustom; + JXL_RETURN_IF_ERROR(red.SetValue(xy.r)); + JXL_RETURN_IF_ERROR(green.SetValue(xy.g)); + JXL_RETURN_IF_ERROR(blue.SetValue(xy.b)); + return true; + } + + CIExy GetWhitePoint() const { + JXL_DASSERT(have_fields); + CIExy xy; + switch (white_point) { + case WhitePoint::kCustom: + return white.GetValue(); + + case WhitePoint::kD65: + xy.x = 0.3127; + xy.y = 0.3290; + return xy; + + case WhitePoint::kDCI: + // From https://ieeexplore.ieee.org/document/7290729 C.2 page 11 + xy.x = 0.314; + xy.y = 0.351; + return xy; + + case WhitePoint::kE: + xy.x = xy.y = 1.0 / 3; + return xy; + } + JXL_UNREACHABLE("Invalid WhitePoint %u", + static_cast<uint32_t>(white_point)); + } + + Status SetWhitePoint(const CIExy& xy) { + JXL_DASSERT(have_fields); + if (xy.x == 0.0 || xy.y == 0.0) { + return JXL_FAILURE("Invalid white point %f %f", xy.x, xy.y); + } + if (ApproxEq(xy.x, 0.3127) && ApproxEq(xy.y, 0.3290)) { + white_point = WhitePoint::kD65; + return true; + } + if (ApproxEq(xy.x, 1.0 / 3) && ApproxEq(xy.y, 1.0 / 3)) { + white_point = WhitePoint::kE; + return true; + } + if (ApproxEq(xy.x, 0.314) && ApproxEq(xy.y, 0.351)) { + white_point = WhitePoint::kDCI; + return true; + } + white_point = WhitePoint::kCustom; + return white.SetValue(xy); + } + + // Checks if the color spaces (including white point / primaries) are the + // same, but ignores the transfer function, rendering intent and ICC bytes. + bool SameColorSpace(const ColorEncoding& other) const { + if (color_space != other.color_space) return false; + + if (white_point != other.white_point) return false; + if (white_point == WhitePoint::kCustom) { + if (!white.IsSame(other.white)) { + return false; + } + } + + if (HasPrimaries() != other.HasPrimaries()) return false; + if (HasPrimaries()) { + if (primaries != other.primaries) return false; + if (primaries == Primaries::kCustom) { + if (!red.IsSame(other.red)) return false; + if (!green.IsSame(other.green)) return false; + if (!blue.IsSame(other.blue)) return false; + } + } + return true; + } + + // Checks if the color space and transfer function are the same, ignoring + // rendering intent and ICC bytes + bool SameColorEncoding(const ColorEncoding& other) const { + return SameColorSpace(other) && tf.IsSame(other.tf); + } + + // Returns true if all fields have been initialized (possibly to kUnknown). + // Returns false if the ICC profile is invalid or decoding it fails. + Status SetFieldsFromICC(IccBytes&& new_icc, const JxlCmsInterface& cms) { + // In case parsing fails, mark the ColorEncoding as invalid. + JXL_ASSERT(!new_icc.empty()); + color_space = ColorSpace::kUnknown; + tf.transfer_function = TransferFunction::kUnknown; + icc.clear(); + + JxlColorEncoding external; + JXL_BOOL new_cmyk; + JXL_RETURN_IF_ERROR(cms.set_fields_from_icc(cms.set_fields_data, + new_icc.data(), new_icc.size(), + &external, &new_cmyk)); + cmyk = new_cmyk; + JXL_RETURN_IF_ERROR(FromExternal(external)); + icc = std::move(new_icc); + return true; + } + + JxlColorEncoding ToExternal() const { + JxlColorEncoding external = {}; + if (!have_fields) { + external.color_space = JXL_COLOR_SPACE_UNKNOWN; + external.primaries = JXL_PRIMARIES_CUSTOM; + external.rendering_intent = JXL_RENDERING_INTENT_PERCEPTUAL; //? + external.transfer_function = JXL_TRANSFER_FUNCTION_UNKNOWN; + external.white_point = JXL_WHITE_POINT_CUSTOM; + return external; + } + external.color_space = static_cast<JxlColorSpace>(color_space); + + external.white_point = static_cast<JxlWhitePoint>(white_point); + + CIExy wp = GetWhitePoint(); + external.white_point_xy[0] = wp.x; + external.white_point_xy[1] = wp.y; + + if (external.color_space == JXL_COLOR_SPACE_RGB || + external.color_space == JXL_COLOR_SPACE_UNKNOWN) { + external.primaries = static_cast<JxlPrimaries>(primaries); + PrimariesCIExy p = GetPrimaries(); + external.primaries_red_xy[0] = p.r.x; + external.primaries_red_xy[1] = p.r.y; + external.primaries_green_xy[0] = p.g.x; + external.primaries_green_xy[1] = p.g.y; + external.primaries_blue_xy[0] = p.b.x; + external.primaries_blue_xy[1] = p.b.y; + } + + if (tf.have_gamma) { + external.transfer_function = JXL_TRANSFER_FUNCTION_GAMMA; + external.gamma = tf.GetGamma(); + } else { + external.transfer_function = + static_cast<JxlTransferFunction>(tf.GetTransferFunction()); + external.gamma = 0; + } + + external.rendering_intent = + static_cast<JxlRenderingIntent>(rendering_intent); + return external; + } + + // NB: does not create ICC. + Status FromExternal(const JxlColorEncoding& external) { + // TODO(eustas): update non-serializable on call-site + color_space = static_cast<ColorSpace>(external.color_space); + + JXL_RETURN_IF_ERROR( + WhitePointFromExternal(external.white_point, &white_point)); + if (external.white_point == JXL_WHITE_POINT_CUSTOM) { + CIExy wp; + wp.x = external.white_point_xy[0]; + wp.y = external.white_point_xy[1]; + JXL_RETURN_IF_ERROR(SetWhitePoint(wp)); + } + + if (external.color_space == JXL_COLOR_SPACE_RGB || + external.color_space == JXL_COLOR_SPACE_UNKNOWN) { + JXL_RETURN_IF_ERROR( + PrimariesFromExternal(external.primaries, &primaries)); + if (external.primaries == JXL_PRIMARIES_CUSTOM) { + PrimariesCIExy primaries; + primaries.r.x = external.primaries_red_xy[0]; + primaries.r.y = external.primaries_red_xy[1]; + primaries.g.x = external.primaries_green_xy[0]; + primaries.g.y = external.primaries_green_xy[1]; + primaries.b.x = external.primaries_blue_xy[0]; + primaries.b.y = external.primaries_blue_xy[1]; + JXL_RETURN_IF_ERROR(SetPrimaries(primaries)); + } + } + CustomTransferFunction tf; + if (external.transfer_function == JXL_TRANSFER_FUNCTION_GAMMA) { + JXL_RETURN_IF_ERROR(tf.SetGamma(external.gamma)); + } else { + TransferFunction tf_enum; + // JXL_TRANSFER_FUNCTION_GAMMA is not handled by this function since + // there's no internal enum value for it. + JXL_RETURN_IF_ERROR(ConvertExternalToInternalTransferFunction( + external.transfer_function, &tf_enum)); + tf.SetTransferFunction(tf_enum); + } + this->tf = tf; + + JXL_RETURN_IF_ERROR(RenderingIntentFromExternal(external.rendering_intent, + &rendering_intent)); + + icc.clear(); + + return true; + } +}; + +} // namespace cms +} // namespace jxl + +#endif // LIB_JXL_CMS_COLOR_ENCODING_CMS_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/jxl_cms.cc b/third_party/jpeg-xl/lib/jxl/cms/jxl_cms.cc new file mode 100644 index 0000000000..dd00b8b81f --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/jxl_cms.cc @@ -0,0 +1,1343 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#include <jxl/cms.h> + +#ifndef JPEGXL_ENABLE_SKCMS +#define JPEGXL_ENABLE_SKCMS 0 +#endif + +#include <jxl/cms_interface.h> + +#include <algorithm> +#include <array> +#include <cmath> +#include <cstddef> +#include <cstdint> +#include <cstring> +#include <memory> + +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "lib/jxl/cms/jxl_cms.cc" +#include <hwy/foreach_target.h> +#include <hwy/highway.h> + +#include "lib/jxl/base/compiler_specific.h" +#include "lib/jxl/base/matrix_ops.h" +#include "lib/jxl/base/printf_macros.h" +#include "lib/jxl/base/span.h" +#include "lib/jxl/base/status.h" +#include "lib/jxl/cms/jxl_cms_internal.h" +#include "lib/jxl/cms/transfer_functions-inl.h" +#include "lib/jxl/color_encoding_internal.h" +#if JPEGXL_ENABLE_SKCMS +#include "skcms.h" +#else // JPEGXL_ENABLE_SKCMS +#include "lcms2.h" +#include "lcms2_plugin.h" +#endif // JPEGXL_ENABLE_SKCMS + +#define JXL_CMS_VERBOSE 0 + +// Define these only once. We can't use HWY_ONCE here because it is defined as +// 1 only on the last pass. +#ifndef LIB_JXL_JXL_CMS_CC +#define LIB_JXL_JXL_CMS_CC + +namespace jxl { +namespace { + +using ::jxl::cms::ColorEncoding; + +struct JxlCms { +#if JPEGXL_ENABLE_SKCMS + IccBytes icc_src, icc_dst; + skcms_ICCProfile profile_src, profile_dst; +#else + void* lcms_transform; +#endif + + // These fields are used when the HLG OOTF or inverse OOTF must be applied. + bool apply_hlg_ootf; + size_t hlg_ootf_num_channels; + // Y component of the primaries. + std::array<float, 3> hlg_ootf_luminances; + + size_t channels_src; + size_t channels_dst; + + std::vector<float> src_storage; + std::vector<float*> buf_src; + std::vector<float> dst_storage; + std::vector<float*> buf_dst; + + float intensity_target; + bool skip_lcms = false; + ExtraTF preprocess = ExtraTF::kNone; + ExtraTF postprocess = ExtraTF::kNone; +}; + +Status ApplyHlgOotf(JxlCms* t, float* JXL_RESTRICT buf, size_t xsize, + bool forward); +} // namespace +} // namespace jxl + +#endif // LIB_JXL_JXL_CMS_CC + +HWY_BEFORE_NAMESPACE(); +namespace jxl { +namespace HWY_NAMESPACE { + +#if JXL_CMS_VERBOSE >= 2 +const size_t kX = 0; // pixel index, multiplied by 3 for RGB +#endif + +// xform_src = UndoGammaCompression(buf_src). +Status BeforeTransform(JxlCms* t, const float* buf_src, float* xform_src, + size_t buf_size) { + switch (t->preprocess) { + case ExtraTF::kNone: + JXL_DASSERT(false); // unreachable + break; + + case ExtraTF::kPQ: { + HWY_FULL(float) df; + TF_PQ tf_pq(t->intensity_target); + for (size_t i = 0; i < buf_size; i += Lanes(df)) { + const auto val = Load(df, buf_src + i); + const auto result = tf_pq.DisplayFromEncoded(df, val); + Store(result, df, xform_src + i); + } +#if JXL_CMS_VERBOSE >= 2 + printf("pre in %.4f %.4f %.4f undoPQ %.4f %.4f %.4f\n", buf_src[3 * kX], + buf_src[3 * kX + 1], buf_src[3 * kX + 2], xform_src[3 * kX], + xform_src[3 * kX + 1], xform_src[3 * kX + 2]); +#endif + break; + } + + case ExtraTF::kHLG: + for (size_t i = 0; i < buf_size; ++i) { + xform_src[i] = static_cast<float>( + TF_HLG_Base::DisplayFromEncoded(static_cast<double>(buf_src[i]))); + } + if (t->apply_hlg_ootf) { + JXL_RETURN_IF_ERROR( + ApplyHlgOotf(t, xform_src, buf_size, /*forward=*/true)); + } +#if JXL_CMS_VERBOSE >= 2 + printf("pre in %.4f %.4f %.4f undoHLG %.4f %.4f %.4f\n", buf_src[3 * kX], + buf_src[3 * kX + 1], buf_src[3 * kX + 2], xform_src[3 * kX], + xform_src[3 * kX + 1], xform_src[3 * kX + 2]); +#endif + break; + + case ExtraTF::kSRGB: + HWY_FULL(float) df; + for (size_t i = 0; i < buf_size; i += Lanes(df)) { + const auto val = Load(df, buf_src + i); + const auto result = TF_SRGB().DisplayFromEncoded(val); + Store(result, df, xform_src + i); + } +#if JXL_CMS_VERBOSE >= 2 + printf("pre in %.4f %.4f %.4f undoSRGB %.4f %.4f %.4f\n", buf_src[3 * kX], + buf_src[3 * kX + 1], buf_src[3 * kX + 2], xform_src[3 * kX], + xform_src[3 * kX + 1], xform_src[3 * kX + 2]); +#endif + break; + } + return true; +} + +// Applies gamma compression in-place. +Status AfterTransform(JxlCms* t, float* JXL_RESTRICT buf_dst, size_t buf_size) { + switch (t->postprocess) { + case ExtraTF::kNone: + JXL_DASSERT(false); // unreachable + break; + case ExtraTF::kPQ: { + HWY_FULL(float) df; + TF_PQ tf_pq(t->intensity_target); + for (size_t i = 0; i < buf_size; i += Lanes(df)) { + const auto val = Load(df, buf_dst + i); + const auto result = tf_pq.EncodedFromDisplay(df, val); + Store(result, df, buf_dst + i); + } +#if JXL_CMS_VERBOSE >= 2 + printf("after PQ enc %.4f %.4f %.4f\n", buf_dst[3 * kX], + buf_dst[3 * kX + 1], buf_dst[3 * kX + 2]); +#endif + break; + } + case ExtraTF::kHLG: + if (t->apply_hlg_ootf) { + JXL_RETURN_IF_ERROR( + ApplyHlgOotf(t, buf_dst, buf_size, /*forward=*/false)); + } + for (size_t i = 0; i < buf_size; ++i) { + buf_dst[i] = static_cast<float>( + TF_HLG_Base::EncodedFromDisplay(static_cast<double>(buf_dst[i]))); + } +#if JXL_CMS_VERBOSE >= 2 + printf("after HLG enc %.4f %.4f %.4f\n", buf_dst[3 * kX], + buf_dst[3 * kX + 1], buf_dst[3 * kX + 2]); +#endif + break; + case ExtraTF::kSRGB: + HWY_FULL(float) df; + for (size_t i = 0; i < buf_size; i += Lanes(df)) { + const auto val = Load(df, buf_dst + i); + const auto result = TF_SRGB().EncodedFromDisplay(df, val); + Store(result, df, buf_dst + i); + } +#if JXL_CMS_VERBOSE >= 2 + printf("after SRGB enc %.4f %.4f %.4f\n", buf_dst[3 * kX], + buf_dst[3 * kX + 1], buf_dst[3 * kX + 2]); +#endif + break; + } + return true; +} + +Status DoColorSpaceTransform(void* cms_data, const size_t thread, + const float* buf_src, float* buf_dst, + size_t xsize) { + // No lock needed. + JxlCms* t = reinterpret_cast<JxlCms*>(cms_data); + + const float* xform_src = buf_src; // Read-only. + if (t->preprocess != ExtraTF::kNone) { + float* mutable_xform_src = t->buf_src[thread]; // Writable buffer. + JXL_RETURN_IF_ERROR(BeforeTransform(t, buf_src, mutable_xform_src, + xsize * t->channels_src)); + xform_src = mutable_xform_src; + } + +#if JPEGXL_ENABLE_SKCMS + if (t->channels_src == 1 && !t->skip_lcms) { + // Expand from 1 to 3 channels, starting from the end in case + // xform_src == t->buf_src[thread]. + float* mutable_xform_src = t->buf_src[thread]; + for (size_t i = 0; i < xsize; ++i) { + const size_t x = xsize - i - 1; + mutable_xform_src[x * 3] = mutable_xform_src[x * 3 + 1] = + mutable_xform_src[x * 3 + 2] = xform_src[x]; + } + xform_src = mutable_xform_src; + } +#else + if (t->channels_src == 4 && !t->skip_lcms) { + // LCMS does CMYK in a weird way: 0 = white, 100 = max ink + float* mutable_xform_src = t->buf_src[thread]; + for (size_t x = 0; x < xsize * 4; ++x) { + mutable_xform_src[x] = 100.f - 100.f * mutable_xform_src[x]; + } + xform_src = mutable_xform_src; + } +#endif + +#if JXL_CMS_VERBOSE >= 2 + // Save inputs for printing before in-place transforms overwrite them. + const float in0 = xform_src[3 * kX + 0]; + const float in1 = xform_src[3 * kX + 1]; + const float in2 = xform_src[3 * kX + 2]; +#endif + + if (t->skip_lcms) { + if (buf_dst != xform_src) { + memcpy(buf_dst, xform_src, xsize * t->channels_src * sizeof(*buf_dst)); + } // else: in-place, no need to copy + } else { +#if JPEGXL_ENABLE_SKCMS + JXL_CHECK( + skcms_Transform(xform_src, + (t->channels_src == 4 ? skcms_PixelFormat_RGBA_ffff + : skcms_PixelFormat_RGB_fff), + skcms_AlphaFormat_Opaque, &t->profile_src, buf_dst, + skcms_PixelFormat_RGB_fff, skcms_AlphaFormat_Opaque, + &t->profile_dst, xsize)); +#else // JPEGXL_ENABLE_SKCMS + cmsDoTransform(t->lcms_transform, xform_src, buf_dst, + static_cast<cmsUInt32Number>(xsize)); +#endif // JPEGXL_ENABLE_SKCMS + } +#if JXL_CMS_VERBOSE >= 2 + printf("xform skip%d: %.4f %.4f %.4f (%p) -> (%p) %.4f %.4f %.4f\n", + t->skip_lcms, in0, in1, in2, xform_src, buf_dst, buf_dst[3 * kX], + buf_dst[3 * kX + 1], buf_dst[3 * kX + 2]); +#endif + +#if JPEGXL_ENABLE_SKCMS + if (t->channels_dst == 1 && !t->skip_lcms) { + // Contract back from 3 to 1 channel, this time forward. + float* grayscale_buf_dst = t->buf_dst[thread]; + for (size_t x = 0; x < xsize; ++x) { + grayscale_buf_dst[x] = buf_dst[x * 3]; + } + buf_dst = grayscale_buf_dst; + } +#endif + + if (t->postprocess != ExtraTF::kNone) { + JXL_RETURN_IF_ERROR(AfterTransform(t, buf_dst, xsize * t->channels_dst)); + } + return true; +} + +// NOLINTNEXTLINE(google-readability-namespace-comments) +} // namespace HWY_NAMESPACE +} // namespace jxl +HWY_AFTER_NAMESPACE(); + +#if HWY_ONCE +namespace jxl { +namespace { + +HWY_EXPORT(DoColorSpaceTransform); +int DoColorSpaceTransform(void* t, size_t thread, const float* buf_src, + float* buf_dst, size_t xsize) { + return HWY_DYNAMIC_DISPATCH(DoColorSpaceTransform)(t, thread, buf_src, + buf_dst, xsize); +} + +// Define to 1 on OS X as a workaround for older LCMS lacking MD5. +#define JXL_CMS_OLD_VERSION 0 + +#if JPEGXL_ENABLE_SKCMS + +JXL_MUST_USE_RESULT CIExy CIExyFromXYZ(const float XYZ[3]) { + const float factor = 1.f / (XYZ[0] + XYZ[1] + XYZ[2]); + CIExy xy; + xy.x = XYZ[0] * factor; + xy.y = XYZ[1] * factor; + return xy; +} + +#else // JPEGXL_ENABLE_SKCMS +// (LCMS interface requires xyY but we omit the Y for white points/primaries.) + +JXL_MUST_USE_RESULT CIExy CIExyFromxyY(const cmsCIExyY& xyY) { + CIExy xy; + xy.x = xyY.x; + xy.y = xyY.y; + return xy; +} + +JXL_MUST_USE_RESULT CIExy CIExyFromXYZ(const cmsCIEXYZ& XYZ) { + cmsCIExyY xyY; + cmsXYZ2xyY(/*Dest=*/&xyY, /*Source=*/&XYZ); + return CIExyFromxyY(xyY); +} + +JXL_MUST_USE_RESULT cmsCIEXYZ D50_XYZ() { + // Quantized D50 as stored in ICC profiles. + return {0.96420288, 1.0, 0.82490540}; +} + +// RAII + +struct ProfileDeleter { + void operator()(void* p) { cmsCloseProfile(p); } +}; +using Profile = std::unique_ptr<void, ProfileDeleter>; + +struct TransformDeleter { + void operator()(void* p) { cmsDeleteTransform(p); } +}; +using Transform = std::unique_ptr<void, TransformDeleter>; + +struct CurveDeleter { + void operator()(cmsToneCurve* p) { cmsFreeToneCurve(p); } +}; +using Curve = std::unique_ptr<cmsToneCurve, CurveDeleter>; + +Status CreateProfileXYZ(const cmsContext context, + Profile* JXL_RESTRICT profile) { + profile->reset(cmsCreateXYZProfileTHR(context)); + if (profile->get() == nullptr) return JXL_FAILURE("Failed to create XYZ"); + return true; +} + +#endif // !JPEGXL_ENABLE_SKCMS + +#if JPEGXL_ENABLE_SKCMS +// IMPORTANT: icc must outlive profile. +Status DecodeProfile(const uint8_t* icc, size_t size, + skcms_ICCProfile* const profile) { + if (!skcms_Parse(icc, size, profile)) { + return JXL_FAILURE("Failed to parse ICC profile with %" PRIuS " bytes", + size); + } + return true; +} +#else // JPEGXL_ENABLE_SKCMS +Status DecodeProfile(const cmsContext context, Span<const uint8_t> icc, + Profile* profile) { + profile->reset(cmsOpenProfileFromMemTHR(context, icc.data(), icc.size())); + if (profile->get() == nullptr) { + return JXL_FAILURE("Failed to decode profile"); + } + + // WARNING: due to the LCMS MD5 issue mentioned above, many existing + // profiles have incorrect MD5, so do not even bother checking them nor + // generating warning clutter. + + return true; +} +#endif // JPEGXL_ENABLE_SKCMS + +#if JPEGXL_ENABLE_SKCMS + +ColorSpace ColorSpaceFromProfile(const skcms_ICCProfile& profile) { + switch (profile.data_color_space) { + case skcms_Signature_RGB: + case skcms_Signature_CMYK: + // spec says CMYK is encoded as RGB (the kBlack extra channel signals that + // it is actually CMYK) + return ColorSpace::kRGB; + case skcms_Signature_Gray: + return ColorSpace::kGray; + default: + return ColorSpace::kUnknown; + } +} + +// vector_out := matmul(matrix, vector_in) +void MatrixProduct(const skcms_Matrix3x3& matrix, const float vector_in[3], + float vector_out[3]) { + for (int i = 0; i < 3; ++i) { + vector_out[i] = 0; + for (int j = 0; j < 3; ++j) { + vector_out[i] += matrix.vals[i][j] * vector_in[j]; + } + } +} + +// Returns white point that was specified when creating the profile. +JXL_MUST_USE_RESULT Status UnadaptedWhitePoint(const skcms_ICCProfile& profile, + CIExy* out) { + float media_white_point_XYZ[3]; + if (!skcms_GetWTPT(&profile, media_white_point_XYZ)) { + return JXL_FAILURE("ICC profile does not contain WhitePoint tag"); + } + skcms_Matrix3x3 CHAD; + if (!skcms_GetCHAD(&profile, &CHAD)) { + // If there is no chromatic adaptation matrix, it means that the white point + // is already unadapted. + *out = CIExyFromXYZ(media_white_point_XYZ); + return true; + } + // Otherwise, it has been adapted to the PCS white point using said matrix, + // and the adaptation needs to be undone. + skcms_Matrix3x3 inverse_CHAD; + if (!skcms_Matrix3x3_invert(&CHAD, &inverse_CHAD)) { + return JXL_FAILURE("Non-invertible ChromaticAdaptation matrix"); + } + float unadapted_white_point_XYZ[3]; + MatrixProduct(inverse_CHAD, media_white_point_XYZ, unadapted_white_point_XYZ); + *out = CIExyFromXYZ(unadapted_white_point_XYZ); + return true; +} + +Status IdentifyPrimaries(const skcms_ICCProfile& profile, + const CIExy& wp_unadapted, ColorEncoding* c) { + if (!c->HasPrimaries()) return true; + + skcms_Matrix3x3 CHAD, inverse_CHAD; + if (skcms_GetCHAD(&profile, &CHAD)) { + JXL_RETURN_IF_ERROR(skcms_Matrix3x3_invert(&CHAD, &inverse_CHAD)); + } else { + static constexpr skcms_Matrix3x3 kLMSFromXYZ = { + {{0.8951, 0.2664, -0.1614}, + {-0.7502, 1.7135, 0.0367}, + {0.0389, -0.0685, 1.0296}}}; + static constexpr skcms_Matrix3x3 kXYZFromLMS = { + {{0.9869929, -0.1470543, 0.1599627}, + {0.4323053, 0.5183603, 0.0492912}, + {-0.0085287, 0.0400428, 0.9684867}}}; + static constexpr float kWpD50XYZ[3] = {0.96420288, 1.0, 0.82490540}; + float wp_unadapted_XYZ[3]; + JXL_RETURN_IF_ERROR( + CIEXYZFromWhiteCIExy(wp_unadapted.x, wp_unadapted.y, wp_unadapted_XYZ)); + float wp_D50_LMS[3], wp_unadapted_LMS[3]; + MatrixProduct(kLMSFromXYZ, kWpD50XYZ, wp_D50_LMS); + MatrixProduct(kLMSFromXYZ, wp_unadapted_XYZ, wp_unadapted_LMS); + inverse_CHAD = {{{wp_unadapted_LMS[0] / wp_D50_LMS[0], 0, 0}, + {0, wp_unadapted_LMS[1] / wp_D50_LMS[1], 0}, + {0, 0, wp_unadapted_LMS[2] / wp_D50_LMS[2]}}}; + inverse_CHAD = skcms_Matrix3x3_concat(&kXYZFromLMS, &inverse_CHAD); + inverse_CHAD = skcms_Matrix3x3_concat(&inverse_CHAD, &kLMSFromXYZ); + } + + float XYZ[3]; + PrimariesCIExy primaries; + CIExy* const chromaticities[] = {&primaries.r, &primaries.g, &primaries.b}; + for (int i = 0; i < 3; ++i) { + float RGB[3] = {}; + RGB[i] = 1; + skcms_Transform(RGB, skcms_PixelFormat_RGB_fff, skcms_AlphaFormat_Opaque, + &profile, XYZ, skcms_PixelFormat_RGB_fff, + skcms_AlphaFormat_Opaque, skcms_XYZD50_profile(), 1); + float unadapted_XYZ[3]; + MatrixProduct(inverse_CHAD, XYZ, unadapted_XYZ); + *chromaticities[i] = CIExyFromXYZ(unadapted_XYZ); + } + return c->SetPrimaries(primaries); +} + +bool IsApproximatelyEqual(const skcms_ICCProfile& profile, + const ColorEncoding& JXL_RESTRICT c) { + IccBytes bytes; + if (!MaybeCreateProfile(c.ToExternal(), &bytes)) { + return false; + } + + skcms_ICCProfile profile_test; + if (!DecodeProfile(bytes.data(), bytes.size(), &profile_test)) { + return false; + } + + if (!skcms_ApproximatelyEqualProfiles(&profile_test, &profile)) { + return false; + } + + return true; +} + +void DetectTransferFunction(const skcms_ICCProfile& profile, + ColorEncoding* JXL_RESTRICT c) { + JXL_CHECK(c->color_space != ColorSpace::kXYB); + + float gamma[3] = {}; + if (profile.has_trc) { + const auto IsGamma = [](const skcms_TransferFunction& tf) { + return tf.a == 1 && tf.b == 0 && + /* if b and d are zero, it is fine for c not to be */ tf.d == 0 && + tf.e == 0 && tf.f == 0; + }; + for (int i = 0; i < 3; ++i) { + if (profile.trc[i].table_entries == 0 && + IsGamma(profile.trc->parametric)) { + gamma[i] = 1.f / profile.trc->parametric.g; + } else { + skcms_TransferFunction approximate_tf; + float max_error; + if (skcms_ApproximateCurve(&profile.trc[i], &approximate_tf, + &max_error)) { + if (IsGamma(approximate_tf)) { + gamma[i] = 1.f / approximate_tf.g; + } + } + } + } + } + if (gamma[0] != 0 && std::abs(gamma[0] - gamma[1]) < 1e-4f && + std::abs(gamma[1] - gamma[2]) < 1e-4f) { + if (c->tf.SetGamma(gamma[0])) { + if (IsApproximatelyEqual(profile, *c)) return; + } + } + + for (TransferFunction tf : Values<TransferFunction>()) { + // Can only create profile from known transfer function. + if (tf == TransferFunction::kUnknown) continue; + c->tf.SetTransferFunction(tf); + if (IsApproximatelyEqual(profile, *c)) return; + } + + c->tf.SetTransferFunction(TransferFunction::kUnknown); +} + +#else // JPEGXL_ENABLE_SKCMS + +uint32_t Type32(const ColorEncoding& c, bool cmyk) { + if (cmyk) return TYPE_CMYK_FLT; + if (c.color_space == ColorSpace::kGray) return TYPE_GRAY_FLT; + return TYPE_RGB_FLT; +} + +uint32_t Type64(const ColorEncoding& c) { + if (c.color_space == ColorSpace::kGray) return TYPE_GRAY_DBL; + return TYPE_RGB_DBL; +} + +ColorSpace ColorSpaceFromProfile(const Profile& profile) { + switch (cmsGetColorSpace(profile.get())) { + case cmsSigRgbData: + case cmsSigCmykData: + return ColorSpace::kRGB; + case cmsSigGrayData: + return ColorSpace::kGray; + default: + return ColorSpace::kUnknown; + } +} + +// "profile1" is pre-decoded to save time in DetectTransferFunction. +Status ProfileEquivalentToICC(const cmsContext context, const Profile& profile1, + const IccBytes& icc, const ColorEncoding& c) { + const uint32_t type_src = Type64(c); + + Profile profile2; + JXL_RETURN_IF_ERROR(DecodeProfile(context, Bytes(icc), &profile2)); + + Profile profile_xyz; + JXL_RETURN_IF_ERROR(CreateProfileXYZ(context, &profile_xyz)); + + const uint32_t intent = INTENT_RELATIVE_COLORIMETRIC; + const uint32_t flags = cmsFLAGS_NOOPTIMIZE | cmsFLAGS_BLACKPOINTCOMPENSATION | + cmsFLAGS_HIGHRESPRECALC; + Transform xform1(cmsCreateTransformTHR(context, profile1.get(), type_src, + profile_xyz.get(), TYPE_XYZ_DBL, + intent, flags)); + Transform xform2(cmsCreateTransformTHR(context, profile2.get(), type_src, + profile_xyz.get(), TYPE_XYZ_DBL, + intent, flags)); + if (xform1 == nullptr || xform2 == nullptr) { + return JXL_FAILURE("Failed to create transform"); + } + + double in[3]; + double out1[3]; + double out2[3]; + + // Uniformly spaced samples from very dark to almost fully bright. + const double init = 1E-3; + const double step = 0.2; + + if (c.color_space == ColorSpace::kGray) { + // Finer sampling and replicate each component. + for (in[0] = init; in[0] < 1.0; in[0] += step / 8) { + cmsDoTransform(xform1.get(), in, out1, 1); + cmsDoTransform(xform2.get(), in, out2, 1); + if (!cms::ApproxEq(out1[0], out2[0], 2E-4)) { + return false; + } + } + } else { + for (in[0] = init; in[0] < 1.0; in[0] += step) { + for (in[1] = init; in[1] < 1.0; in[1] += step) { + for (in[2] = init; in[2] < 1.0; in[2] += step) { + cmsDoTransform(xform1.get(), in, out1, 1); + cmsDoTransform(xform2.get(), in, out2, 1); + for (size_t i = 0; i < 3; ++i) { + if (!cms::ApproxEq(out1[i], out2[i], 2E-4)) { + return false; + } + } + } + } + } + } + + return true; +} + +// Returns white point that was specified when creating the profile. +// NOTE: we can't just use cmsSigMediaWhitePointTag because its interpretation +// differs between ICC versions. +JXL_MUST_USE_RESULT cmsCIEXYZ UnadaptedWhitePoint(const cmsContext context, + const Profile& profile, + const ColorEncoding& c) { + const cmsCIEXYZ* white_point = static_cast<const cmsCIEXYZ*>( + cmsReadTag(profile.get(), cmsSigMediaWhitePointTag)); + if (white_point != nullptr && + cmsReadTag(profile.get(), cmsSigChromaticAdaptationTag) == nullptr) { + // No chromatic adaptation matrix: the white point is already unadapted. + return *white_point; + } + + cmsCIEXYZ XYZ = {1.0, 1.0, 1.0}; + Profile profile_xyz; + if (!CreateProfileXYZ(context, &profile_xyz)) return XYZ; + // Array arguments are one per profile. + cmsHPROFILE profiles[2] = {profile.get(), profile_xyz.get()}; + // Leave white point unchanged - that is what we're trying to extract. + cmsUInt32Number intents[2] = {INTENT_ABSOLUTE_COLORIMETRIC, + INTENT_ABSOLUTE_COLORIMETRIC}; + cmsBool black_compensation[2] = {0, 0}; + cmsFloat64Number adaption[2] = {0.0, 0.0}; + // Only transforming a single pixel, so skip expensive optimizations. + cmsUInt32Number flags = cmsFLAGS_NOOPTIMIZE | cmsFLAGS_HIGHRESPRECALC; + Transform xform(cmsCreateExtendedTransform( + context, 2, profiles, black_compensation, intents, adaption, nullptr, 0, + Type64(c), TYPE_XYZ_DBL, flags)); + if (!xform) return XYZ; // TODO(lode): return error + + // xy are relative, so magnitude does not matter if we ignore output Y. + const cmsFloat64Number in[3] = {1.0, 1.0, 1.0}; + cmsDoTransform(xform.get(), in, &XYZ.X, 1); + return XYZ; +} + +Status IdentifyPrimaries(const cmsContext context, const Profile& profile, + const cmsCIEXYZ& wp_unadapted, ColorEncoding* c) { + if (!c->HasPrimaries()) return true; + if (ColorSpaceFromProfile(profile) == ColorSpace::kUnknown) return true; + + // These were adapted to the profile illuminant before storing in the profile. + const cmsCIEXYZ* adapted_r = static_cast<const cmsCIEXYZ*>( + cmsReadTag(profile.get(), cmsSigRedColorantTag)); + const cmsCIEXYZ* adapted_g = static_cast<const cmsCIEXYZ*>( + cmsReadTag(profile.get(), cmsSigGreenColorantTag)); + const cmsCIEXYZ* adapted_b = static_cast<const cmsCIEXYZ*>( + cmsReadTag(profile.get(), cmsSigBlueColorantTag)); + + cmsCIEXYZ converted_rgb[3]; + if (adapted_r == nullptr || adapted_g == nullptr || adapted_b == nullptr) { + // No colorant tag, determine the XYZ coordinates of the primaries by + // converting from the colorspace. + Profile profile_xyz; + if (!CreateProfileXYZ(context, &profile_xyz)) { + return JXL_FAILURE("Failed to retrieve colorants"); + } + // Array arguments are one per profile. + cmsHPROFILE profiles[2] = {profile.get(), profile_xyz.get()}; + cmsUInt32Number intents[2] = {INTENT_RELATIVE_COLORIMETRIC, + INTENT_RELATIVE_COLORIMETRIC}; + cmsBool black_compensation[2] = {0, 0}; + cmsFloat64Number adaption[2] = {0.0, 0.0}; + // Only transforming three pixels, so skip expensive optimizations. + cmsUInt32Number flags = cmsFLAGS_NOOPTIMIZE | cmsFLAGS_HIGHRESPRECALC; + Transform xform(cmsCreateExtendedTransform( + context, 2, profiles, black_compensation, intents, adaption, nullptr, 0, + Type64(*c), TYPE_XYZ_DBL, flags)); + if (!xform) return JXL_FAILURE("Failed to retrieve colorants"); + + const cmsFloat64Number in[9] = {1.0, 0.0, 0.0, 0.0, 1.0, + 0.0, 0.0, 0.0, 1.0}; + cmsDoTransform(xform.get(), in, &converted_rgb->X, 3); + adapted_r = &converted_rgb[0]; + adapted_g = &converted_rgb[1]; + adapted_b = &converted_rgb[2]; + } + + // TODO(janwas): no longer assume Bradford and D50. + // Undo the chromatic adaptation. + const cmsCIEXYZ d50 = D50_XYZ(); + + cmsCIEXYZ r, g, b; + cmsAdaptToIlluminant(&r, &d50, &wp_unadapted, adapted_r); + cmsAdaptToIlluminant(&g, &d50, &wp_unadapted, adapted_g); + cmsAdaptToIlluminant(&b, &d50, &wp_unadapted, adapted_b); + + const PrimariesCIExy rgb = {CIExyFromXYZ(r), CIExyFromXYZ(g), + CIExyFromXYZ(b)}; + return c->SetPrimaries(rgb); +} + +void DetectTransferFunction(const cmsContext context, const Profile& profile, + ColorEncoding* JXL_RESTRICT c) { + JXL_CHECK(c->color_space != ColorSpace::kXYB); + + float gamma = 0; + if (const auto* gray_trc = reinterpret_cast<const cmsToneCurve*>( + cmsReadTag(profile.get(), cmsSigGrayTRCTag))) { + const double estimated_gamma = + cmsEstimateGamma(gray_trc, /*precision=*/1e-4); + if (estimated_gamma > 0) { + gamma = 1. / estimated_gamma; + } + } else { + float rgb_gamma[3] = {}; + int i = 0; + for (const auto tag : + {cmsSigRedTRCTag, cmsSigGreenTRCTag, cmsSigBlueTRCTag}) { + if (const auto* trc = reinterpret_cast<const cmsToneCurve*>( + cmsReadTag(profile.get(), tag))) { + const double estimated_gamma = + cmsEstimateGamma(trc, /*precision=*/1e-4); + if (estimated_gamma > 0) { + rgb_gamma[i] = 1. / estimated_gamma; + } + } + ++i; + } + if (rgb_gamma[0] != 0 && std::abs(rgb_gamma[0] - rgb_gamma[1]) < 1e-4f && + std::abs(rgb_gamma[1] - rgb_gamma[2]) < 1e-4f) { + gamma = rgb_gamma[0]; + } + } + + if (gamma != 0 && c->tf.SetGamma(gamma)) { + IccBytes icc_test; + if (MaybeCreateProfile(c->ToExternal(), &icc_test) && + ProfileEquivalentToICC(context, profile, icc_test, *c)) { + return; + } + } + + for (TransferFunction tf : Values<TransferFunction>()) { + // Can only create profile from known transfer function. + if (tf == TransferFunction::kUnknown) continue; + + c->tf.SetTransferFunction(tf); + + IccBytes icc_test; + if (MaybeCreateProfile(c->ToExternal(), &icc_test) && + ProfileEquivalentToICC(context, profile, icc_test, *c)) { + return; + } + } + + c->tf.SetTransferFunction(TransferFunction::kUnknown); +} + +void ErrorHandler(cmsContext context, cmsUInt32Number code, const char* text) { + JXL_WARNING("LCMS error %u: %s", code, text); +} + +// Returns a context for the current thread, creating it if necessary. +cmsContext GetContext() { + static thread_local void* context_; + if (context_ == nullptr) { + context_ = cmsCreateContext(nullptr, nullptr); + JXL_ASSERT(context_ != nullptr); + + cmsSetLogErrorHandlerTHR(static_cast<cmsContext>(context_), &ErrorHandler); + } + return static_cast<cmsContext>(context_); +} + +#endif // JPEGXL_ENABLE_SKCMS + +Status GetPrimariesLuminances(const ColorEncoding& encoding, + float luminances[3]) { + // Explanation: + // We know that the three primaries must sum to white: + // + // [Xr, Xg, Xb; [1; [Xw; + // Yr, Yg, Yb; × 1; = Yw; + // Zr, Zg, Zb] 1] Zw] + // + // By noting that X = x·(X+Y+Z), Y = y·(X+Y+Z) and Z = z·(X+Y+Z) (note the + // lower case indicating chromaticity), and factoring the totals (X+Y+Z) out + // of the left matrix and into the all-ones vector, we get: + // + // [xr, xg, xb; [Xr + Yr + Zr; [Xw; + // yr, yg, yb; × Xg + Yg + Zg; = Yw; + // zr, zg, zb] Xb + Yb + Zb] Zw] + // + // Which makes it apparent that we can compute those totals as: + // + // [Xr + Yr + Zr; inv([xr, xg, xb; [Xw; + // Xg + Yg + Zg; = yr, yg, yb; × Yw; + // Xb + Yb + Zb] zr, zg, zb]) Zw] + // + // From there, by multiplying each total by its corresponding y, we get Y for + // that primary. + + float white_XYZ[3]; + CIExy wp = encoding.GetWhitePoint(); + JXL_RETURN_IF_ERROR(CIEXYZFromWhiteCIExy(wp.x, wp.y, white_XYZ)); + + const PrimariesCIExy primaries = encoding.GetPrimaries(); + double chromaticities[3][3] = { + {primaries.r.x, primaries.g.x, primaries.b.x}, + {primaries.r.y, primaries.g.y, primaries.b.y}, + {1 - primaries.r.x - primaries.r.y, 1 - primaries.g.x - primaries.g.y, + 1 - primaries.b.x - primaries.b.y}}; + JXL_RETURN_IF_ERROR(Inv3x3Matrix(&chromaticities[0][0])); + const double ys[3] = {primaries.r.y, primaries.g.y, primaries.b.y}; + for (size_t i = 0; i < 3; ++i) { + luminances[i] = ys[i] * (chromaticities[i][0] * white_XYZ[0] + + chromaticities[i][1] * white_XYZ[1] + + chromaticities[i][2] * white_XYZ[2]); + } + return true; +} + +Status ApplyHlgOotf(JxlCms* t, float* JXL_RESTRICT buf, size_t xsize, + bool forward) { + if (295 <= t->intensity_target && t->intensity_target <= 305) { + // The gamma is approximately 1 so this can essentially be skipped. + return true; + } + float gamma = 1.2f * std::pow(1.111f, std::log2(t->intensity_target * 1e-3f)); + if (!forward) gamma = 1.f / gamma; + + switch (t->hlg_ootf_num_channels) { + case 1: + for (size_t x = 0; x < xsize; ++x) { + buf[x] = std::pow(buf[x], gamma); + } + break; + + case 3: + for (size_t x = 0; x < xsize; x += 3) { + const float luminance = buf[x] * t->hlg_ootf_luminances[0] + + buf[x + 1] * t->hlg_ootf_luminances[1] + + buf[x + 2] * t->hlg_ootf_luminances[2]; + const float ratio = std::pow(luminance, gamma - 1); + if (std::isfinite(ratio)) { + buf[x] *= ratio; + buf[x + 1] *= ratio; + buf[x + 2] *= ratio; + if (forward && gamma < 1) { + // If gamma < 1, the ratio above will be > 1 which can push bright + // saturated highlights out of gamut. There are several possible + // ways to bring them back in-gamut; this one preserves hue and + // saturation at the slight expense of luminance. If !forward, the + // previously-applied forward OOTF with gamma > 1 already pushed + // those highlights down and we are simply putting them back where + // they were so this is not necessary. + const float maximum = + std::max(buf[x], std::max(buf[x + 1], buf[x + 2])); + if (maximum > 1) { + const float normalizer = 1.f / maximum; + buf[x] *= normalizer; + buf[x + 1] *= normalizer; + buf[x + 2] *= normalizer; + } + } + } + } + break; + + default: + return JXL_FAILURE("HLG OOTF not implemented for %" PRIuS " channels", + t->hlg_ootf_num_channels); + } + return true; +} + +bool IsKnownTransferFunction(jxl::cms::TransferFunction tf) { + using TF = jxl::cms::TransferFunction; + // All but kUnknown + return tf == TF::k709 || tf == TF::kLinear || tf == TF::kSRGB || + tf == TF::kPQ || tf == TF::kDCI || tf == TF::kHLG; +} + +constexpr uint8_t kColorPrimariesP3_D65 = 12; + +bool IsKnownColorPrimaries(uint8_t color_primaries) { + using P = jxl::cms::Primaries; + // All but kCustom + if (color_primaries == kColorPrimariesP3_D65) return true; + const auto p = static_cast<Primaries>(color_primaries); + return p == P::kSRGB || p == P::k2100 || p == P::kP3; +} + +bool ApplyCICP(const uint8_t color_primaries, + const uint8_t transfer_characteristics, + const uint8_t matrix_coefficients, const uint8_t full_range, + ColorEncoding* JXL_RESTRICT c) { + if (matrix_coefficients != 0) return false; + if (full_range != 1) return false; + + const auto primaries = static_cast<Primaries>(color_primaries); + const auto tf = static_cast<TransferFunction>(transfer_characteristics); + if (!IsKnownTransferFunction(tf)) return false; + if (!IsKnownColorPrimaries(color_primaries)) return false; + c->color_space = ColorSpace::kRGB; + c->tf.SetTransferFunction(tf); + if (primaries == Primaries::kP3) { + c->white_point = WhitePoint::kDCI; + c->primaries = Primaries::kP3; + } else if (color_primaries == kColorPrimariesP3_D65) { + c->white_point = WhitePoint::kD65; + c->primaries = Primaries::kP3; + } else { + c->white_point = WhitePoint::kD65; + c->primaries = primaries; + } + return true; +} + +JXL_BOOL JxlCmsSetFieldsFromICC(void* user_data, const uint8_t* icc_data, + size_t icc_size, JxlColorEncoding* c, + JXL_BOOL* cmyk) { + if (c == nullptr) return JXL_FALSE; + if (cmyk == nullptr) return JXL_FALSE; + + *cmyk = JXL_FALSE; + + // In case parsing fails, mark the ColorEncoding as invalid. + c->color_space = JXL_COLOR_SPACE_UNKNOWN; + c->transfer_function = JXL_TRANSFER_FUNCTION_UNKNOWN; + + if (icc_size == 0) return JXL_FAILURE("Empty ICC profile"); + + ColorEncoding c_enc; + +#if JPEGXL_ENABLE_SKCMS + if (icc_size < 128) { + return JXL_FAILURE("ICC file too small"); + } + + skcms_ICCProfile profile; + JXL_RETURN_IF_ERROR(skcms_Parse(icc_data, icc_size, &profile)); + + // skcms does not return the rendering intent, so get it from the file. It + // is encoded as big-endian 32-bit integer in bytes 60..63. + uint32_t rendering_intent32 = icc_data[67]; + if (rendering_intent32 > 3 || icc_data[64] != 0 || icc_data[65] != 0 || + icc_data[66] != 0) { + return JXL_FAILURE("Invalid rendering intent %u\n", rendering_intent32); + } + // ICC and RenderingIntent have the same values (0..3). + c_enc.rendering_intent = static_cast<RenderingIntent>(rendering_intent32); + + if (profile.has_CICP && + ApplyCICP(profile.CICP.color_primaries, + profile.CICP.transfer_characteristics, + profile.CICP.matrix_coefficients, + profile.CICP.video_full_range_flag, &c_enc)) { + *c = c_enc.ToExternal(); + return true; + } + + c_enc.color_space = ColorSpaceFromProfile(profile); + *cmyk = (profile.data_color_space == skcms_Signature_CMYK); + + CIExy wp_unadapted; + JXL_RETURN_IF_ERROR(UnadaptedWhitePoint(profile, &wp_unadapted)); + JXL_RETURN_IF_ERROR(c_enc.SetWhitePoint(wp_unadapted)); + + // Relies on color_space. + JXL_RETURN_IF_ERROR(IdentifyPrimaries(profile, wp_unadapted, &c_enc)); + + // Relies on color_space/white point/primaries being set already. + DetectTransferFunction(profile, &c_enc); +#else // JPEGXL_ENABLE_SKCMS + + const cmsContext context = GetContext(); + + Profile profile; + JXL_RETURN_IF_ERROR( + DecodeProfile(context, Bytes(icc_data, icc_size), &profile)); + + const cmsUInt32Number rendering_intent32 = + cmsGetHeaderRenderingIntent(profile.get()); + if (rendering_intent32 > 3) { + return JXL_FAILURE("Invalid rendering intent %u\n", rendering_intent32); + } + // ICC and RenderingIntent have the same values (0..3). + c_enc.rendering_intent = static_cast<RenderingIntent>(rendering_intent32); + + static constexpr size_t kCICPSize = 12; + static constexpr auto kCICPSignature = + static_cast<cmsTagSignature>(0x63696370); + uint8_t cicp_buffer[kCICPSize]; + if (cmsReadRawTag(profile.get(), kCICPSignature, cicp_buffer, kCICPSize) == + kCICPSize && + ApplyCICP(cicp_buffer[8], cicp_buffer[9], cicp_buffer[10], + cicp_buffer[11], &c_enc)) { + *c = c_enc.ToExternal(); + return true; + } + + c_enc.color_space = ColorSpaceFromProfile(profile); + if (cmsGetColorSpace(profile.get()) == cmsSigCmykData) { + *cmyk = JXL_TRUE; + *c = c_enc.ToExternal(); + return true; + } + + const cmsCIEXYZ wp_unadapted = UnadaptedWhitePoint(context, profile, c_enc); + JXL_RETURN_IF_ERROR(c_enc.SetWhitePoint(CIExyFromXYZ(wp_unadapted))); + + // Relies on color_space. + JXL_RETURN_IF_ERROR( + IdentifyPrimaries(context, profile, wp_unadapted, &c_enc)); + + // Relies on color_space/white point/primaries being set already. + DetectTransferFunction(context, profile, &c_enc); + +#endif // JPEGXL_ENABLE_SKCMS + + *c = c_enc.ToExternal(); + return true; +} + +} // namespace + +namespace { + +void JxlCmsDestroy(void* cms_data) { + if (cms_data == nullptr) return; + JxlCms* t = reinterpret_cast<JxlCms*>(cms_data); +#if !JPEGXL_ENABLE_SKCMS + TransformDeleter()(t->lcms_transform); +#endif + delete t; +} + +void AllocateBuffer(size_t length, size_t num_threads, + std::vector<float>* storage, std::vector<float*>* view) { + constexpr size_t kAlign = 128 / sizeof(float); + size_t stride = RoundUpTo(length, kAlign); + storage->resize(stride * num_threads + kAlign); + intptr_t addr = reinterpret_cast<intptr_t>(storage->data()); + size_t offset = + (RoundUpTo(addr, kAlign * sizeof(float)) - addr) / sizeof(float); + view->clear(); + view->reserve(num_threads); + for (size_t i = 0; i < num_threads; ++i) { + view->emplace_back(storage->data() + offset + i * stride); + } +} + +void* JxlCmsInit(void* init_data, size_t num_threads, size_t xsize, + const JxlColorProfile* input, const JxlColorProfile* output, + float intensity_target) { + JXL_ASSERT(init_data != nullptr); + auto cms = static_cast<const JxlCmsInterface*>(init_data); + auto t = jxl::make_unique<JxlCms>(); + IccBytes icc_src, icc_dst; + if (input->icc.size == 0) { + JXL_NOTIFY_ERROR("JxlCmsInit: empty input ICC"); + return nullptr; + } + if (output->icc.size == 0) { + JXL_NOTIFY_ERROR("JxlCmsInit: empty OUTPUT ICC"); + return nullptr; + } + icc_src.assign(input->icc.data, input->icc.data + input->icc.size); + ColorEncoding c_src; + if (!c_src.SetFieldsFromICC(std::move(icc_src), *cms)) { + JXL_NOTIFY_ERROR("JxlCmsInit: failed to parse input ICC"); + return nullptr; + } + icc_dst.assign(output->icc.data, output->icc.data + output->icc.size); + ColorEncoding c_dst; + if (!c_dst.SetFieldsFromICC(std::move(icc_dst), *cms)) { + JXL_NOTIFY_ERROR("JxlCmsInit: failed to parse output ICC"); + return nullptr; + } +#if JXL_CMS_VERBOSE + printf("%s -> %s\n", Description(c_src).c_str(), Description(c_dst).c_str()); +#endif + +#if JPEGXL_ENABLE_SKCMS + if (!DecodeProfile(input->icc.data, input->icc.size, &t->profile_src)) { + JXL_NOTIFY_ERROR("JxlCmsInit: skcms failed to parse input ICC"); + return nullptr; + } + if (!DecodeProfile(output->icc.data, output->icc.size, &t->profile_dst)) { + JXL_NOTIFY_ERROR("JxlCmsInit: skcms failed to parse output ICC"); + return nullptr; + } +#else // JPEGXL_ENABLE_SKCMS + const cmsContext context = GetContext(); + Profile profile_src, profile_dst; + if (!DecodeProfile(context, Bytes(c_src.icc), &profile_src)) { + JXL_NOTIFY_ERROR("JxlCmsInit: lcms failed to parse input ICC"); + return nullptr; + } + if (!DecodeProfile(context, Bytes(c_dst.icc), &profile_dst)) { + JXL_NOTIFY_ERROR("JxlCmsInit: lcms failed to parse output ICC"); + return nullptr; + } +#endif // JPEGXL_ENABLE_SKCMS + + t->skip_lcms = false; + if (c_src.SameColorEncoding(c_dst)) { + t->skip_lcms = true; +#if JXL_CMS_VERBOSE + printf("Skip CMS\n"); +#endif + } + + t->apply_hlg_ootf = c_src.tf.IsHLG() != c_dst.tf.IsHLG(); + if (t->apply_hlg_ootf) { + const ColorEncoding* c_hlg = c_src.tf.IsHLG() ? &c_src : &c_dst; + t->hlg_ootf_num_channels = c_hlg->Channels(); + if (t->hlg_ootf_num_channels == 3 && + !GetPrimariesLuminances(*c_hlg, t->hlg_ootf_luminances.data())) { + JXL_NOTIFY_ERROR( + "JxlCmsInit: failed to compute the luminances of primaries"); + return nullptr; + } + } + + // Special-case SRGB <=> linear if the primaries / white point are the same, + // or any conversion where PQ or HLG is involved: + bool src_linear = c_src.tf.IsLinear(); + const bool dst_linear = c_dst.tf.IsLinear(); + + if (c_src.tf.IsPQ() || c_src.tf.IsHLG() || + (c_src.tf.IsSRGB() && dst_linear && c_src.SameColorSpace(c_dst))) { + // Construct new profile as if the data were already/still linear. + ColorEncoding c_linear_src = c_src; + c_linear_src.tf.SetTransferFunction(TransferFunction::kLinear); +#if JPEGXL_ENABLE_SKCMS + skcms_ICCProfile new_src; +#else // JPEGXL_ENABLE_SKCMS + Profile new_src; +#endif // JPEGXL_ENABLE_SKCMS + // Only enable ExtraTF if profile creation succeeded. + if (MaybeCreateProfile(c_linear_src.ToExternal(), &icc_src) && +#if JPEGXL_ENABLE_SKCMS + DecodeProfile(icc_src.data(), icc_src.size(), &new_src)) { +#else // JPEGXL_ENABLE_SKCMS + DecodeProfile(context, Bytes(icc_src), &new_src)) { +#endif // JPEGXL_ENABLE_SKCMS +#if JXL_CMS_VERBOSE + printf("Special HLG/PQ/sRGB -> linear\n"); +#endif +#if JPEGXL_ENABLE_SKCMS + t->icc_src = std::move(icc_src); + t->profile_src = new_src; +#else // JPEGXL_ENABLE_SKCMS + profile_src.swap(new_src); +#endif // JPEGXL_ENABLE_SKCMS + t->preprocess = c_src.tf.IsSRGB() + ? ExtraTF::kSRGB + : (c_src.tf.IsPQ() ? ExtraTF::kPQ : ExtraTF::kHLG); + c_src = c_linear_src; + src_linear = true; + } else { + if (t->apply_hlg_ootf) { + JXL_NOTIFY_ERROR( + "Failed to create extra linear source profile, and HLG OOTF " + "required"); + return nullptr; + } + JXL_WARNING("Failed to create extra linear destination profile"); + } + } + + if (c_dst.tf.IsPQ() || c_dst.tf.IsHLG() || + (c_dst.tf.IsSRGB() && src_linear && c_src.SameColorSpace(c_dst))) { + ColorEncoding c_linear_dst = c_dst; + c_linear_dst.tf.SetTransferFunction(TransferFunction::kLinear); +#if JPEGXL_ENABLE_SKCMS + skcms_ICCProfile new_dst; +#else // JPEGXL_ENABLE_SKCMS + Profile new_dst; +#endif // JPEGXL_ENABLE_SKCMS + // Only enable ExtraTF if profile creation succeeded. + if (MaybeCreateProfile(c_linear_dst.ToExternal(), &icc_dst) && +#if JPEGXL_ENABLE_SKCMS + DecodeProfile(icc_dst.data(), icc_dst.size(), &new_dst)) { +#else // JPEGXL_ENABLE_SKCMS + DecodeProfile(context, Bytes(icc_dst), &new_dst)) { +#endif // JPEGXL_ENABLE_SKCMS +#if JXL_CMS_VERBOSE + printf("Special linear -> HLG/PQ/sRGB\n"); +#endif +#if JPEGXL_ENABLE_SKCMS + t->icc_dst = std::move(icc_dst); + t->profile_dst = new_dst; +#else // JPEGXL_ENABLE_SKCMS + profile_dst.swap(new_dst); +#endif // JPEGXL_ENABLE_SKCMS + t->postprocess = c_dst.tf.IsSRGB() + ? ExtraTF::kSRGB + : (c_dst.tf.IsPQ() ? ExtraTF::kPQ : ExtraTF::kHLG); + c_dst = c_linear_dst; + } else { + if (t->apply_hlg_ootf) { + JXL_NOTIFY_ERROR( + "Failed to create extra linear destination profile, and inverse " + "HLG OOTF required"); + return nullptr; + } + JXL_WARNING("Failed to create extra linear destination profile"); + } + } + + if (c_src.SameColorEncoding(c_dst)) { +#if JXL_CMS_VERBOSE + printf("Same intermediary linear profiles, skipping CMS\n"); +#endif + t->skip_lcms = true; + } + +#if JPEGXL_ENABLE_SKCMS + if (!skcms_MakeUsableAsDestination(&t->profile_dst)) { + JXL_NOTIFY_ERROR( + "Failed to make %s usable as a color transform destination", + ColorEncodingDescription(c_dst.ToExternal()).c_str()); + return nullptr; + } +#endif // JPEGXL_ENABLE_SKCMS + + // Not including alpha channel (copied separately). + const size_t channels_src = (c_src.cmyk ? 4 : c_src.Channels()); + const size_t channels_dst = c_dst.Channels(); + JXL_CHECK(channels_src == channels_dst || + (channels_src == 4 && channels_dst == 3)); +#if JXL_CMS_VERBOSE + printf("Channels: %" PRIuS "; Threads: %" PRIuS "\n", channels_src, + num_threads); +#endif + +#if !JPEGXL_ENABLE_SKCMS + // Type includes color space (XYZ vs RGB), so can be different. + const uint32_t type_src = Type32(c_src, channels_src == 4); + const uint32_t type_dst = Type32(c_dst, false); + const uint32_t intent = static_cast<uint32_t>(c_dst.rendering_intent); + // Use cmsFLAGS_NOCACHE to disable the 1-pixel cache and make calling + // cmsDoTransform() thread-safe. + const uint32_t flags = cmsFLAGS_NOCACHE | cmsFLAGS_BLACKPOINTCOMPENSATION | + cmsFLAGS_HIGHRESPRECALC; + t->lcms_transform = + cmsCreateTransformTHR(context, profile_src.get(), type_src, + profile_dst.get(), type_dst, intent, flags); + if (t->lcms_transform == nullptr) { + JXL_NOTIFY_ERROR("Failed to create transform"); + return nullptr; + } +#endif // !JPEGXL_ENABLE_SKCMS + + // Ideally LCMS would convert directly from External to Image3. However, + // cmsDoTransformLineStride only accepts 32-bit BytesPerPlaneIn, whereas our + // planes can be more than 4 GiB apart. Hence, transform inputs/outputs must + // be interleaved. Calling cmsDoTransform for each pixel is expensive + // (indirect call). We therefore transform rows, which requires per-thread + // buffers. To avoid separate allocations, we use the rows of an image. + // Because LCMS apparently also cannot handle <= 16 bit inputs and 32-bit + // outputs (or vice versa), we use floating point input/output. + t->channels_src = channels_src; + t->channels_dst = channels_dst; + size_t actual_channels_src = channels_src; + size_t actual_channels_dst = channels_dst; +#if JPEGXL_ENABLE_SKCMS + // SkiaCMS doesn't support grayscale float buffers, so we create space for RGB + // float buffers anyway. + actual_channels_src = (channels_src == 4 ? 4 : 3); + actual_channels_dst = 3; +#endif + AllocateBuffer(xsize * actual_channels_src, num_threads, &t->src_storage, + &t->buf_src); + AllocateBuffer(xsize * actual_channels_dst, num_threads, &t->dst_storage, + &t->buf_dst); + t->intensity_target = intensity_target; + return t.release(); +} + +float* JxlCmsGetSrcBuf(void* cms_data, size_t thread) { + JxlCms* t = reinterpret_cast<JxlCms*>(cms_data); + return t->buf_src[thread]; +} + +float* JxlCmsGetDstBuf(void* cms_data, size_t thread) { + JxlCms* t = reinterpret_cast<JxlCms*>(cms_data); + return t->buf_dst[thread]; +} + +} // namespace + +extern "C" { + +JXL_CMS_EXPORT const JxlCmsInterface* JxlGetDefaultCms() { + static constexpr JxlCmsInterface kInterface = { + /*set_fields_data=*/nullptr, + /*set_fields_from_icc=*/&JxlCmsSetFieldsFromICC, + /*init_data=*/const_cast<void*>(static_cast<const void*>(&kInterface)), + /*init=*/&JxlCmsInit, + /*get_src_buf=*/&JxlCmsGetSrcBuf, + /*get_dst_buf=*/&JxlCmsGetDstBuf, + /*run=*/&DoColorSpaceTransform, + /*destroy=*/&JxlCmsDestroy}; + return &kInterface; +} + +} // extern "C" + +} // namespace jxl +#endif // HWY_ONCE diff --git a/third_party/jpeg-xl/lib/jxl/cms/jxl_cms_internal.h b/third_party/jpeg-xl/lib/jxl/cms/jxl_cms_internal.h new file mode 100644 index 0000000000..c00fe82d8c --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/jxl_cms_internal.h @@ -0,0 +1,1083 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#ifndef LIB_JXL_CMS_JXL_CMS_INTERNAL_H_ +#define LIB_JXL_CMS_JXL_CMS_INTERNAL_H_ + +// ICC profiles and color space conversions. + +#include <jxl/color_encoding.h> + +#include <algorithm> +#include <cmath> +#include <cstddef> +#include <cstdint> +#include <cstring> +#include <string> +#include <vector> + +#include "lib/jxl/base/common.h" +#include "lib/jxl/base/compiler_specific.h" +#include "lib/jxl/base/matrix_ops.h" +#include "lib/jxl/base/span.h" // Bytes +#include "lib/jxl/base/status.h" +#include "lib/jxl/cms/opsin_params.h" +#include "lib/jxl/cms/tone_mapping.h" +#include "lib/jxl/cms/transfer_functions.h" + +#ifndef JXL_ENABLE_3D_ICC_TONEMAPPING +#define JXL_ENABLE_3D_ICC_TONEMAPPING 1 +#endif + +namespace jxl { + +enum class ExtraTF { + kNone, + kPQ, + kHLG, + kSRGB, +}; + +static Status PrimariesToXYZ(float rx, float ry, float gx, float gy, float bx, + float by, float wx, float wy, float matrix[9]) { + bool ok = (wx >= 0) && (wx <= 1) && (wy > 0) && (wy <= 1); + if (!ok) { + return JXL_FAILURE("Invalid white point"); + } + // TODO(lode): also require rx, ry, gx, gy, bx, to be in range 0-1? ICC + // profiles in theory forbid negative XYZ values, but in practice the ACES P0 + // color space uses a negative y for the blue primary. + float primaries[9] = { + rx, gx, bx, ry, gy, by, 1.0f - rx - ry, 1.0f - gx - gy, 1.0f - bx - by}; + float primaries_inv[9]; + memcpy(primaries_inv, primaries, sizeof(float) * 9); + JXL_RETURN_IF_ERROR(Inv3x3Matrix(primaries_inv)); + + float w[3] = {wx / wy, 1.0f, (1.0f - wx - wy) / wy}; + // 1 / tiny float can still overflow + JXL_RETURN_IF_ERROR(std::isfinite(w[0]) && std::isfinite(w[2])); + float xyz[3]; + Mul3x3Vector(primaries_inv, w, xyz); + + float a[9] = { + xyz[0], 0, 0, 0, xyz[1], 0, 0, 0, xyz[2], + }; + + Mul3x3Matrix(primaries, a, matrix); + return true; +} + +/* Chromatic adaptation matrices*/ +constexpr float kBradford[9] = { + 0.8951f, 0.2664f, -0.1614f, -0.7502f, 1.7135f, + 0.0367f, 0.0389f, -0.0685f, 1.0296f, +}; +constexpr float kBradfordInv[9] = { + 0.9869929f, -0.1470543f, 0.1599627f, 0.4323053f, 0.5183603f, + 0.0492912f, -0.0085287f, 0.0400428f, 0.9684867f, +}; + +// Adapts whitepoint x, y to D50 +static Status AdaptToXYZD50(float wx, float wy, float matrix[9]) { + bool ok = (wx >= 0) && (wx <= 1) && (wy > 0) && (wy <= 1); + if (!ok) { + // Out of range values can cause division through zero + // further down with the bradford adaptation too. + return JXL_FAILURE("Invalid white point"); + } + float w[3] = {wx / wy, 1.0f, (1.0f - wx - wy) / wy}; + // 1 / tiny float can still overflow + JXL_RETURN_IF_ERROR(std::isfinite(w[0]) && std::isfinite(w[2])); + float w50[3] = {0.96422f, 1.0f, 0.82521f}; + + float lms[3]; + float lms50[3]; + + Mul3x3Vector(kBradford, w, lms); + Mul3x3Vector(kBradford, w50, lms50); + + if (lms[0] == 0 || lms[1] == 0 || lms[2] == 0) { + return JXL_FAILURE("Invalid white point"); + } + float a[9] = { + // /----> 0, 1, 2, 3, /----> 4, 5, 6, 7, /----> 8, + lms50[0] / lms[0], 0, 0, 0, lms50[1] / lms[1], 0, 0, 0, lms50[2] / lms[2], + }; + if (!std::isfinite(a[0]) || !std::isfinite(a[4]) || !std::isfinite(a[8])) { + return JXL_FAILURE("Invalid white point"); + } + + float b[9]; + Mul3x3Matrix(a, kBradford, b); + Mul3x3Matrix(kBradfordInv, b, matrix); + + return true; +} + +static Status PrimariesToXYZD50(float rx, float ry, float gx, float gy, + float bx, float by, float wx, float wy, + float matrix[9]) { + float toXYZ[9]; + JXL_RETURN_IF_ERROR(PrimariesToXYZ(rx, ry, gx, gy, bx, by, wx, wy, toXYZ)); + float d50[9]; + JXL_RETURN_IF_ERROR(AdaptToXYZD50(wx, wy, d50)); + + Mul3x3Matrix(d50, toXYZ, matrix); + return true; +} + +static Status ToneMapPixel(const JxlColorEncoding& c, const float in[3], + uint8_t pcslab_out[3]) { + float primaries_XYZ[9]; + JXL_RETURN_IF_ERROR(PrimariesToXYZ( + c.primaries_red_xy[0], c.primaries_red_xy[1], c.primaries_green_xy[0], + c.primaries_green_xy[1], c.primaries_blue_xy[0], c.primaries_blue_xy[1], + c.white_point_xy[0], c.white_point_xy[1], primaries_XYZ)); + const float luminances[3] = {primaries_XYZ[3], primaries_XYZ[4], + primaries_XYZ[5]}; + float linear[3]; + JxlTransferFunction tf = c.transfer_function; + if (tf == JXL_TRANSFER_FUNCTION_PQ) { + for (size_t i = 0; i < 3; ++i) { + linear[i] = TF_PQ_Base::DisplayFromEncoded( + /*display_intensity_target=*/10000.0, in[i]); + } + } else { + for (size_t i = 0; i < 3; ++i) { + linear[i] = TF_HLG_Base::DisplayFromEncoded(in[i]); + } + } + if (tf == JXL_TRANSFER_FUNCTION_PQ) { + Rec2408ToneMapperBase tone_mapper({0, 10000}, {0, 250}, luminances); + tone_mapper.ToneMap(&linear[0], &linear[1], &linear[2]); + } else { + HlgOOTF_Base ootf(/*source_luminance=*/300, /*target_luminance=*/80, + luminances); + ootf.Apply(&linear[0], &linear[1], &linear[2]); + } + GamutMapScalar(&linear[0], &linear[1], &linear[2], luminances, + /*preserve_saturation=*/0.3f); + + float chad[9]; + JXL_RETURN_IF_ERROR( + AdaptToXYZD50(c.white_point_xy[0], c.white_point_xy[1], chad)); + float to_xyzd50[9]; + Mul3x3Matrix(chad, primaries_XYZ, to_xyzd50); + + float xyz[3] = {0, 0, 0}; + for (size_t xyz_c = 0; xyz_c < 3; ++xyz_c) { + for (size_t rgb_c = 0; rgb_c < 3; ++rgb_c) { + xyz[xyz_c] += linear[rgb_c] * to_xyzd50[3 * xyz_c + rgb_c]; + } + } + + const auto lab_f = [](const float x) { + static constexpr float kDelta = 6. / 29; + return x <= kDelta * kDelta * kDelta + ? x * (1 / (3 * kDelta * kDelta)) + 4.f / 29 + : std::cbrt(x); + }; + static constexpr float kXn = 0.964212; + static constexpr float kYn = 1; + static constexpr float kZn = 0.825188; + + const float f_x = lab_f(xyz[0] / kXn); + const float f_y = lab_f(xyz[1] / kYn); + const float f_z = lab_f(xyz[2] / kZn); + + pcslab_out[0] = + static_cast<uint8_t>(.5f + 255.f * Clamp1(1.16f * f_y - .16f, 0.f, 1.f)); + pcslab_out[1] = static_cast<uint8_t>( + .5f + 128.f + Clamp1(500 * (f_x - f_y), -128.f, 127.f)); + pcslab_out[2] = static_cast<uint8_t>( + .5f + 128.f + Clamp1(200 * (f_y - f_z), -128.f, 127.f)); + + return true; +} + +static std::vector<uint16_t> CreateTableCurve(uint32_t N, const ExtraTF tf, + bool tone_map) { + // The generated PQ curve will make room for highlights up to this luminance. + // TODO(sboukortt): make this variable? + static constexpr float kPQIntensityTarget = 10000; + + JXL_ASSERT(N <= 4096); // ICC MFT2 only allows 4K entries + JXL_ASSERT(tf == ExtraTF::kPQ || tf == ExtraTF::kHLG); + + static constexpr float kLuminances[] = {1.f / 3, 1.f / 3, 1.f / 3}; + Rec2408ToneMapperBase tone_mapper({0, kPQIntensityTarget}, + {0, kDefaultIntensityTarget}, kLuminances); + // No point using float - LCMS converts to 16-bit for A2B/MFT. + std::vector<uint16_t> table(N); + for (uint32_t i = 0; i < N; ++i) { + const float x = static_cast<float>(i) / (N - 1); // 1.0 at index N - 1. + const double dx = static_cast<double>(x); + // LCMS requires EOTF (e.g. 2.4 exponent). + double y = (tf == ExtraTF::kHLG) + ? TF_HLG_Base::DisplayFromEncoded(dx) + : TF_PQ_Base::DisplayFromEncoded(kPQIntensityTarget, dx); + if (tone_map && tf == ExtraTF::kPQ && + kPQIntensityTarget > kDefaultIntensityTarget) { + float r = y * 10000 / kPQIntensityTarget, g = r, b = r; + tone_mapper.ToneMap(&r, &g, &b); + y = r; + } + JXL_ASSERT(y >= 0.0); + // Clamp to table range - necessary for HLG. + if (y > 1.0) y = 1.0; + // 1.0 corresponds to table value 0xFFFF. + table[i] = static_cast<uint16_t>(roundf(y * 65535.0)); + } + return table; +} + +static Status CIEXYZFromWhiteCIExy(double wx, double wy, float XYZ[3]) { + // Target Y = 1. + if (std::abs(wy) < 1e-12) return JXL_FAILURE("Y value is too small"); + const float factor = 1 / wy; + XYZ[0] = wx * factor; + XYZ[1] = 1; + XYZ[2] = (1 - wx - wy) * factor; + return true; +} + +namespace detail { + +constexpr bool kEnable3DToneMapping = JXL_ENABLE_3D_ICC_TONEMAPPING; + +static bool CanToneMap(const JxlColorEncoding& encoding) { + // If the color space cannot be represented by a CICP tag in the ICC profile + // then the rest of the profile must unambiguously identify it; we have less + // freedom to do use it for tone mapping. + JxlTransferFunction tf = encoding.transfer_function; + JxlPrimaries p = encoding.primaries; + JxlWhitePoint wp = encoding.white_point; + return encoding.color_space == JXL_COLOR_SPACE_RGB && + (tf == JXL_TRANSFER_FUNCTION_PQ || tf == JXL_TRANSFER_FUNCTION_HLG) && + ((p == JXL_PRIMARIES_P3 && + (wp == JXL_WHITE_POINT_D65 || wp == JXL_WHITE_POINT_DCI)) || + (p != JXL_PRIMARIES_CUSTOM && wp == JXL_WHITE_POINT_D65)); +} + +static void ICCComputeMD5(const std::vector<uint8_t>& data, uint8_t sum[16]) + JXL_NO_SANITIZE("unsigned-integer-overflow") { + std::vector<uint8_t> data64 = data; + data64.push_back(128); + // Add bytes such that ((size + 8) & 63) == 0. + size_t extra = ((64 - ((data64.size() + 8) & 63)) & 63); + data64.resize(data64.size() + extra, 0); + for (uint64_t i = 0; i < 64; i += 8) { + data64.push_back(static_cast<uint64_t>(data.size() << 3u) >> i); + } + + static const uint32_t sineparts[64] = { + 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, + 0xa8304613, 0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, + 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, 0xf61e2562, 0xc040b340, + 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8, + 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, + 0x676f02d9, 0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, + 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70, 0x289b7ec6, 0xeaa127fa, + 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665, + 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, + 0xffeff47d, 0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, + 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391, + }; + static const uint32_t shift[64] = { + 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, + 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, + 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, + 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, + }; + + uint32_t a0 = 0x67452301, b0 = 0xefcdab89, c0 = 0x98badcfe, d0 = 0x10325476; + + for (size_t i = 0; i < data64.size(); i += 64) { + uint32_t a = a0, b = b0, c = c0, d = d0, f, g; + for (size_t j = 0; j < 64; j++) { + if (j < 16) { + f = (b & c) | ((~b) & d); + g = j; + } else if (j < 32) { + f = (d & b) | ((~d) & c); + g = (5 * j + 1) & 0xf; + } else if (j < 48) { + f = b ^ c ^ d; + g = (3 * j + 5) & 0xf; + } else { + f = c ^ (b | (~d)); + g = (7 * j) & 0xf; + } + uint32_t dg0 = data64[i + g * 4 + 0], dg1 = data64[i + g * 4 + 1], + dg2 = data64[i + g * 4 + 2], dg3 = data64[i + g * 4 + 3]; + uint32_t u = dg0 | (dg1 << 8u) | (dg2 << 16u) | (dg3 << 24u); + f += a + sineparts[j] + u; + a = d; + d = c; + c = b; + b += (f << shift[j]) | (f >> (32u - shift[j])); + } + a0 += a; + b0 += b; + c0 += c; + d0 += d; + } + sum[0] = a0; + sum[1] = a0 >> 8u; + sum[2] = a0 >> 16u; + sum[3] = a0 >> 24u; + sum[4] = b0; + sum[5] = b0 >> 8u; + sum[6] = b0 >> 16u; + sum[7] = b0 >> 24u; + sum[8] = c0; + sum[9] = c0 >> 8u; + sum[10] = c0 >> 16u; + sum[11] = c0 >> 24u; + sum[12] = d0; + sum[13] = d0 >> 8u; + sum[14] = d0 >> 16u; + sum[15] = d0 >> 24u; +} + +static Status CreateICCChadMatrix(double wx, double wy, float result[9]) { + float m[9]; + if (wy == 0) { // WhitePoint can not be pitch-black. + return JXL_FAILURE("Invalid WhitePoint"); + } + JXL_RETURN_IF_ERROR(AdaptToXYZD50(wx, wy, m)); + memcpy(result, m, sizeof(float) * 9); + return true; +} + +// Creates RGB to XYZ matrix given RGB primaries and whitepoint in xy. +static Status CreateICCRGBMatrix(double rx, double ry, double gx, double gy, + double bx, double by, double wx, double wy, + float result[9]) { + float m[9]; + JXL_RETURN_IF_ERROR(PrimariesToXYZD50(rx, ry, gx, gy, bx, by, wx, wy, m)); + memcpy(result, m, sizeof(float) * 9); + return true; +} + +static void WriteICCUint32(uint32_t value, size_t pos, + std::vector<uint8_t>* icc) { + if (icc->size() < pos + 4) icc->resize(pos + 4); + (*icc)[pos + 0] = (value >> 24u) & 255; + (*icc)[pos + 1] = (value >> 16u) & 255; + (*icc)[pos + 2] = (value >> 8u) & 255; + (*icc)[pos + 3] = value & 255; +} + +static void WriteICCUint16(uint16_t value, size_t pos, + std::vector<uint8_t>* icc) { + if (icc->size() < pos + 2) icc->resize(pos + 2); + (*icc)[pos + 0] = (value >> 8u) & 255; + (*icc)[pos + 1] = value & 255; +} + +static void WriteICCUint8(uint8_t value, size_t pos, + std::vector<uint8_t>* icc) { + if (icc->size() < pos + 1) icc->resize(pos + 1); + (*icc)[pos] = value; +} + +// Writes a 4-character tag +static void WriteICCTag(const char* value, size_t pos, + std::vector<uint8_t>* icc) { + if (icc->size() < pos + 4) icc->resize(pos + 4); + memcpy(icc->data() + pos, value, 4); +} + +static Status WriteICCS15Fixed16(float value, size_t pos, + std::vector<uint8_t>* icc) { + // "nextafterf" for 32768.0f towards zero are: + // 32767.998046875, 32767.99609375, 32767.994140625 + // Even the first value works well,... + bool ok = (-32767.995f <= value) && (value <= 32767.995f); + if (!ok) return JXL_FAILURE("ICC value is out of range / NaN"); + int32_t i = value * 65536.0f + 0.5f; + // Use two's complement + uint32_t u = static_cast<uint32_t>(i); + WriteICCUint32(u, pos, icc); + return true; +} + +static Status CreateICCHeader(const JxlColorEncoding& c, + std::vector<uint8_t>* header) { + // TODO(lode): choose color management engine name, e.g. "skia" if + // integrated in skia. + static const char* kCmm = "jxl "; + + header->resize(128, 0); + + WriteICCUint32(0, 0, header); // size, correct value filled in at end + WriteICCTag(kCmm, 4, header); + WriteICCUint32(0x04400000u, 8, header); + const char* profile_type = + c.color_space == JXL_COLOR_SPACE_XYB ? "scnr" : "mntr"; + WriteICCTag(profile_type, 12, header); + WriteICCTag(c.color_space == JXL_COLOR_SPACE_GRAY ? "GRAY" : "RGB ", 16, + header); + if (kEnable3DToneMapping && CanToneMap(c)) { + // We are going to use a 3D LUT for tone mapping, which will be more compact + // with an 8-bit LUT to CIELAB than with a 16-bit LUT to XYZ. 8-bit XYZ + // would not be viable due to XYZ being linear, whereas it is fine with + // CIELAB's ~cube root. + WriteICCTag("Lab ", 20, header); + } else { + WriteICCTag("XYZ ", 20, header); + } + + // Three uint32_t's date/time encoding. + // TODO(lode): encode actual date and time, this is a placeholder + uint32_t year = 2019, month = 12, day = 1; + uint32_t hour = 0, minute = 0, second = 0; + WriteICCUint16(year, 24, header); + WriteICCUint16(month, 26, header); + WriteICCUint16(day, 28, header); + WriteICCUint16(hour, 30, header); + WriteICCUint16(minute, 32, header); + WriteICCUint16(second, 34, header); + + WriteICCTag("acsp", 36, header); + WriteICCTag("APPL", 40, header); + WriteICCUint32(0, 44, header); // flags + WriteICCUint32(0, 48, header); // device manufacturer + WriteICCUint32(0, 52, header); // device model + WriteICCUint32(0, 56, header); // device attributes + WriteICCUint32(0, 60, header); // device attributes + WriteICCUint32(static_cast<uint32_t>(c.rendering_intent), 64, header); + + // Mandatory D50 white point of profile connection space + WriteICCUint32(0x0000f6d6, 68, header); + WriteICCUint32(0x00010000, 72, header); + WriteICCUint32(0x0000d32d, 76, header); + + WriteICCTag(kCmm, 80, header); + + return true; +} + +static void AddToICCTagTable(const char* tag, size_t offset, size_t size, + std::vector<uint8_t>* tagtable, + std::vector<size_t>* offsets) { + WriteICCTag(tag, tagtable->size(), tagtable); + // writing true offset deferred to later + WriteICCUint32(0, tagtable->size(), tagtable); + offsets->push_back(offset); + WriteICCUint32(size, tagtable->size(), tagtable); +} + +static void FinalizeICCTag(std::vector<uint8_t>* tags, size_t* offset, + size_t* size) { + while ((tags->size() & 3) != 0) { + tags->push_back(0); + } + *offset += *size; + *size = tags->size() - *offset; +} + +// The input text must be ASCII, writing other characters to UTF-16 is not +// implemented. +static void CreateICCMlucTag(const std::string& text, + std::vector<uint8_t>* tags) { + WriteICCTag("mluc", tags->size(), tags); + WriteICCUint32(0, tags->size(), tags); + WriteICCUint32(1, tags->size(), tags); + WriteICCUint32(12, tags->size(), tags); + WriteICCTag("enUS", tags->size(), tags); + WriteICCUint32(text.size() * 2, tags->size(), tags); + WriteICCUint32(28, tags->size(), tags); + for (size_t i = 0; i < text.size(); i++) { + tags->push_back(0); // prepend 0 for UTF-16 + tags->push_back(text[i]); + } +} + +static Status CreateICCXYZTag(float xyz[3], std::vector<uint8_t>* tags) { + WriteICCTag("XYZ ", tags->size(), tags); + WriteICCUint32(0, tags->size(), tags); + for (size_t i = 0; i < 3; ++i) { + JXL_RETURN_IF_ERROR(WriteICCS15Fixed16(xyz[i], tags->size(), tags)); + } + return true; +} + +static Status CreateICCChadTag(float chad[9], std::vector<uint8_t>* tags) { + WriteICCTag("sf32", tags->size(), tags); + WriteICCUint32(0, tags->size(), tags); + for (size_t i = 0; i < 9; i++) { + JXL_RETURN_IF_ERROR(WriteICCS15Fixed16(chad[i], tags->size(), tags)); + } + return true; +} + +static void MaybeCreateICCCICPTag(const JxlColorEncoding& c, + std::vector<uint8_t>* tags, size_t* offset, + size_t* size, std::vector<uint8_t>* tagtable, + std::vector<size_t>* offsets) { + if (c.color_space != JXL_COLOR_SPACE_RGB) { + return; + } + uint8_t primaries = 0; + if (c.primaries == JXL_PRIMARIES_P3) { + if (c.white_point == JXL_WHITE_POINT_D65) { + primaries = 12; + } else if (c.white_point == JXL_WHITE_POINT_DCI) { + primaries = 11; + } else { + return; + } + } else if (c.primaries != JXL_PRIMARIES_CUSTOM && + c.white_point == JXL_WHITE_POINT_D65) { + primaries = static_cast<uint8_t>(c.primaries); + } else { + return; + } + JxlTransferFunction tf = c.transfer_function; + if (tf == JXL_TRANSFER_FUNCTION_UNKNOWN || + tf == JXL_TRANSFER_FUNCTION_GAMMA) { + return; + } + WriteICCTag("cicp", tags->size(), tags); + WriteICCUint32(0, tags->size(), tags); + WriteICCUint8(primaries, tags->size(), tags); + WriteICCUint8(static_cast<uint8_t>(tf), tags->size(), tags); + // Matrix + WriteICCUint8(0, tags->size(), tags); + // Full range + WriteICCUint8(1, tags->size(), tags); + FinalizeICCTag(tags, offset, size); + AddToICCTagTable("cicp", *offset, *size, tagtable, offsets); +} + +static void CreateICCCurvCurvTag(const std::vector<uint16_t>& curve, + std::vector<uint8_t>* tags) { + size_t pos = tags->size(); + tags->resize(tags->size() + 12 + curve.size() * 2, 0); + WriteICCTag("curv", pos, tags); + WriteICCUint32(0, pos + 4, tags); + WriteICCUint32(curve.size(), pos + 8, tags); + for (size_t i = 0; i < curve.size(); i++) { + WriteICCUint16(curve[i], pos + 12 + i * 2, tags); + } +} + +// Writes 12 + 4*params.size() bytes +static Status CreateICCCurvParaTag(std::vector<float> params, size_t curve_type, + std::vector<uint8_t>* tags) { + WriteICCTag("para", tags->size(), tags); + WriteICCUint32(0, tags->size(), tags); + WriteICCUint16(curve_type, tags->size(), tags); + WriteICCUint16(0, tags->size(), tags); + for (size_t i = 0; i < params.size(); i++) { + JXL_RETURN_IF_ERROR(WriteICCS15Fixed16(params[i], tags->size(), tags)); + } + return true; +} + +static Status CreateICCLutAtoBTagForXYB(std::vector<uint8_t>* tags) { + WriteICCTag("mAB ", tags->size(), tags); + // 4 reserved bytes set to 0 + WriteICCUint32(0, tags->size(), tags); + // number of input channels + WriteICCUint8(3, tags->size(), tags); + // number of output channels + WriteICCUint8(3, tags->size(), tags); + // 2 reserved bytes for padding + WriteICCUint16(0, tags->size(), tags); + // offset to first B curve + WriteICCUint32(32, tags->size(), tags); + // offset to matrix + WriteICCUint32(244, tags->size(), tags); + // offset to first M curve + WriteICCUint32(148, tags->size(), tags); + // offset to CLUT + WriteICCUint32(80, tags->size(), tags); + // offset to first A curve + // (reuse linear B curves) + WriteICCUint32(32, tags->size(), tags); + + // offset = 32 + // no-op curves + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + // offset = 80 + // number of grid points for each input channel + for (int i = 0; i < 16; ++i) { + WriteICCUint8(i < 3 ? 2 : 0, tags->size(), tags); + } + // precision = 2 + WriteICCUint8(2, tags->size(), tags); + // 3 bytes of padding + WriteICCUint8(0, tags->size(), tags); + WriteICCUint16(0, tags->size(), tags); + // 2*2*2*3 entries of 2 bytes each = 48 bytes + const jxl::cms::ColorCube3D& cube = jxl::cms::UnscaledA2BCube(); + for (size_t ix = 0; ix < 2; ++ix) { + for (size_t iy = 0; iy < 2; ++iy) { + for (size_t ib = 0; ib < 2; ++ib) { + const jxl::cms::ColorCube0D& out_f = cube[ix][iy][ib]; + for (int i = 0; i < 3; ++i) { + int32_t val = static_cast<int32_t>(0.5f + 65535 * out_f[i]); + JXL_DASSERT(val >= 0 && val <= 65535); + WriteICCUint16(val, tags->size(), tags); + } + } + } + } + // offset = 148 + // 3 curves with 5 parameters = 3 * (12 + 5 * 4) = 96 bytes + for (size_t i = 0; i < 3; ++i) { + const float b = -jxl::cms::kXYBOffset[i] - + std::cbrt(jxl::cms::kNegOpsinAbsorbanceBiasRGB[i]); + std::vector<float> params = { + 3, + 1.0f / jxl::cms::kXYBScale[i], + b, + 0, // unused + std::max(0.f, -b * jxl::cms::kXYBScale[i]), // make skcms happy + }; + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag(params, 3, tags)); + } + // offset = 244 + const double matrix[] = {1.5170095, -1.1065225, 0.071623, + -0.050022, 0.5683655, -0.018344, + -1.387676, 1.1145555, 0.6857255}; + // 12 * 4 = 48 bytes + for (size_t i = 0; i < 9; ++i) { + JXL_RETURN_IF_ERROR(WriteICCS15Fixed16(matrix[i], tags->size(), tags)); + } + for (size_t i = 0; i < 3; ++i) { + float intercept = 0; + for (size_t j = 0; j < 3; ++j) { + intercept += matrix[i * 3 + j] * jxl::cms::kNegOpsinAbsorbanceBiasRGB[j]; + } + JXL_RETURN_IF_ERROR(WriteICCS15Fixed16(intercept, tags->size(), tags)); + } + return true; +} + +static Status CreateICCLutAtoBTagForHDR(JxlColorEncoding c, + std::vector<uint8_t>* tags) { + static constexpr size_t k3DLutDim = 9; + WriteICCTag("mft1", tags->size(), tags); + // 4 reserved bytes set to 0 + WriteICCUint32(0, tags->size(), tags); + // number of input channels + WriteICCUint8(3, tags->size(), tags); + // number of output channels + WriteICCUint8(3, tags->size(), tags); + // number of CLUT grid points + WriteICCUint8(k3DLutDim, tags->size(), tags); + // 1 reserved bytes for padding + WriteICCUint8(0, tags->size(), tags); + + // Matrix (per specification, must be identity if input is not XYZ) + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + JXL_RETURN_IF_ERROR( + WriteICCS15Fixed16(i == j ? 1.f : 0.f, tags->size(), tags)); + } + } + + // Input tables + for (size_t c = 0; c < 3; ++c) { + for (size_t i = 0; i < 256; ++i) { + WriteICCUint8(i, tags->size(), tags); + } + } + + for (size_t ix = 0; ix < k3DLutDim; ++ix) { + for (size_t iy = 0; iy < k3DLutDim; ++iy) { + for (size_t ib = 0; ib < k3DLutDim; ++ib) { + float f[3] = {ix * (1.0f / (k3DLutDim - 1)), + iy * (1.0f / (k3DLutDim - 1)), + ib * (1.0f / (k3DLutDim - 1))}; + uint8_t pcslab_out[3]; + JXL_RETURN_IF_ERROR(ToneMapPixel(c, f, pcslab_out)); + for (uint8_t val : pcslab_out) { + WriteICCUint8(val, tags->size(), tags); + } + } + } + } + + // Output tables + for (size_t c = 0; c < 3; ++c) { + for (size_t i = 0; i < 256; ++i) { + WriteICCUint8(i, tags->size(), tags); + } + } + + return true; +} + +// Some software (Apple Safari, Preview) requires this. +static Status CreateICCNoOpBToATag(std::vector<uint8_t>* tags) { + WriteICCTag("mBA ", tags->size(), tags); + // 4 reserved bytes set to 0 + WriteICCUint32(0, tags->size(), tags); + // number of input channels + WriteICCUint8(3, tags->size(), tags); + // number of output channels + WriteICCUint8(3, tags->size(), tags); + // 2 reserved bytes for padding + WriteICCUint16(0, tags->size(), tags); + // offset to first B curve + WriteICCUint32(32, tags->size(), tags); + // offset to matrix + WriteICCUint32(0, tags->size(), tags); + // offset to first M curve + WriteICCUint32(0, tags->size(), tags); + // offset to CLUT + WriteICCUint32(0, tags->size(), tags); + // offset to first A curve + WriteICCUint32(0, tags->size(), tags); + + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({1.0f}, 0, tags)); + + return true; +} + +// These strings are baked into Description - do not change. + +static std::string ToString(JxlColorSpace color_space) { + switch (color_space) { + case JXL_COLOR_SPACE_RGB: + return "RGB"; + case JXL_COLOR_SPACE_GRAY: + return "Gra"; + case JXL_COLOR_SPACE_XYB: + return "XYB"; + case JXL_COLOR_SPACE_UNKNOWN: + return "CS?"; + } + // Should not happen - visitor fails if enum is invalid. + JXL_UNREACHABLE("Invalid ColorSpace %u", static_cast<uint32_t>(color_space)); +} + +static std::string ToString(JxlWhitePoint white_point) { + switch (white_point) { + case JXL_WHITE_POINT_D65: + return "D65"; + case JXL_WHITE_POINT_CUSTOM: + return "Cst"; + case JXL_WHITE_POINT_E: + return "EER"; + case JXL_WHITE_POINT_DCI: + return "DCI"; + } + // Should not happen - visitor fails if enum is invalid. + JXL_UNREACHABLE("Invalid WhitePoint %u", static_cast<uint32_t>(white_point)); +} + +static std::string ToString(JxlPrimaries primaries) { + switch (primaries) { + case JXL_PRIMARIES_SRGB: + return "SRG"; + case JXL_PRIMARIES_2100: + return "202"; + case JXL_PRIMARIES_P3: + return "DCI"; + case JXL_PRIMARIES_CUSTOM: + return "Cst"; + } + // Should not happen - visitor fails if enum is invalid. + JXL_UNREACHABLE("Invalid Primaries %u", static_cast<uint32_t>(primaries)); +} + +static std::string ToString(JxlTransferFunction transfer_function) { + switch (transfer_function) { + case JXL_TRANSFER_FUNCTION_SRGB: + return "SRG"; + case JXL_TRANSFER_FUNCTION_LINEAR: + return "Lin"; + case JXL_TRANSFER_FUNCTION_709: + return "709"; + case JXL_TRANSFER_FUNCTION_PQ: + return "PeQ"; + case JXL_TRANSFER_FUNCTION_HLG: + return "HLG"; + case JXL_TRANSFER_FUNCTION_DCI: + return "DCI"; + case JXL_TRANSFER_FUNCTION_UNKNOWN: + return "TF?"; + case JXL_TRANSFER_FUNCTION_GAMMA: + JXL_UNREACHABLE("Invalid TransferFunction: gamma"); + } + // Should not happen - visitor fails if enum is invalid. + JXL_UNREACHABLE("Invalid TransferFunction %u", + static_cast<uint32_t>(transfer_function)); +} + +static std::string ToString(JxlRenderingIntent rendering_intent) { + switch (rendering_intent) { + case JXL_RENDERING_INTENT_PERCEPTUAL: + return "Per"; + case JXL_RENDERING_INTENT_RELATIVE: + return "Rel"; + case JXL_RENDERING_INTENT_SATURATION: + return "Sat"; + case JXL_RENDERING_INTENT_ABSOLUTE: + return "Abs"; + } + // Should not happen - visitor fails if enum is invalid. + JXL_UNREACHABLE("Invalid RenderingIntent %u", + static_cast<uint32_t>(rendering_intent)); +} + +static std::string ColorEncodingDescriptionImpl(const JxlColorEncoding& c) { + std::string d = ToString(c.color_space); + + bool explicit_wp_tf = (c.color_space != JXL_COLOR_SPACE_XYB); + if (explicit_wp_tf) { + d += '_'; + if (c.white_point == JXL_WHITE_POINT_CUSTOM) { + d += jxl::ToString(c.white_point_xy[0]) + ';'; + d += jxl::ToString(c.white_point_xy[1]); + } else { + d += ToString(c.white_point); + } + } + + if ((c.color_space != JXL_COLOR_SPACE_GRAY) && + (c.color_space != JXL_COLOR_SPACE_XYB)) { + d += '_'; + if (c.primaries == JXL_PRIMARIES_CUSTOM) { + d += jxl::ToString(c.primaries_red_xy[0]) + ';'; + d += jxl::ToString(c.primaries_red_xy[1]) + ';'; + d += jxl::ToString(c.primaries_green_xy[0]) + ';'; + d += jxl::ToString(c.primaries_green_xy[1]) + ';'; + d += jxl::ToString(c.primaries_blue_xy[0]) + ';'; + d += jxl::ToString(c.primaries_blue_xy[1]); + } else { + d += ToString(c.primaries); + } + } + + d += '_'; + d += ToString(c.rendering_intent); + + if (explicit_wp_tf) { + JxlTransferFunction tf = c.transfer_function; + d += '_'; + if (tf == JXL_TRANSFER_FUNCTION_GAMMA) { + d += 'g'; + d += jxl::ToString(c.gamma); + } else { + d += ToString(tf); + } + } + return d; +} + +static Status MaybeCreateProfileImpl(const JxlColorEncoding& c, + std::vector<uint8_t>* icc) { + std::vector<uint8_t> header, tagtable, tags; + JxlTransferFunction tf = c.transfer_function; + if (c.color_space == JXL_COLOR_SPACE_UNKNOWN || + tf == JXL_TRANSFER_FUNCTION_UNKNOWN) { + return false; // Not an error + } + + switch (c.color_space) { + case JXL_COLOR_SPACE_RGB: + case JXL_COLOR_SPACE_GRAY: + case JXL_COLOR_SPACE_XYB: + break; // OK + default: + return JXL_FAILURE("Invalid CS %u", + static_cast<unsigned int>(c.color_space)); + } + + if (c.color_space == JXL_COLOR_SPACE_XYB && + c.rendering_intent != JXL_RENDERING_INTENT_PERCEPTUAL) { + return JXL_FAILURE( + "Only perceptual rendering intent implemented for XYB " + "ICC profile."); + } + + JXL_RETURN_IF_ERROR(CreateICCHeader(c, &header)); + + std::vector<size_t> offsets; + // tag count, deferred to later + WriteICCUint32(0, tagtable.size(), &tagtable); + + size_t tag_offset = 0, tag_size = 0; + + CreateICCMlucTag(ColorEncodingDescriptionImpl(c), &tags); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("desc", tag_offset, tag_size, &tagtable, &offsets); + + const std::string copyright = "CC0"; + CreateICCMlucTag(copyright, &tags); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("cprt", tag_offset, tag_size, &tagtable, &offsets); + + // TODO(eustas): isn't it the other way round: gray image has d50 WhitePoint? + if (c.color_space == JXL_COLOR_SPACE_GRAY) { + float wtpt[3]; + JXL_RETURN_IF_ERROR( + CIEXYZFromWhiteCIExy(c.white_point_xy[0], c.white_point_xy[1], wtpt)); + JXL_RETURN_IF_ERROR(CreateICCXYZTag(wtpt, &tags)); + } else { + float d50[3] = {0.964203, 1.0, 0.824905}; + JXL_RETURN_IF_ERROR(CreateICCXYZTag(d50, &tags)); + } + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("wtpt", tag_offset, tag_size, &tagtable, &offsets); + + if (c.color_space != JXL_COLOR_SPACE_GRAY) { + // Chromatic adaptation matrix + float chad[9]; + JXL_RETURN_IF_ERROR( + CreateICCChadMatrix(c.white_point_xy[0], c.white_point_xy[1], chad)); + + JXL_RETURN_IF_ERROR(CreateICCChadTag(chad, &tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("chad", tag_offset, tag_size, &tagtable, &offsets); + } + + if (c.color_space == JXL_COLOR_SPACE_RGB) { + MaybeCreateICCCICPTag(c, &tags, &tag_offset, &tag_size, &tagtable, + &offsets); + + float m[9]; + JXL_RETURN_IF_ERROR(CreateICCRGBMatrix( + c.primaries_red_xy[0], c.primaries_red_xy[1], c.primaries_green_xy[0], + c.primaries_green_xy[1], c.primaries_blue_xy[0], c.primaries_blue_xy[1], + c.white_point_xy[0], c.white_point_xy[1], m)); + float r[3] = {m[0], m[3], m[6]}; + float g[3] = {m[1], m[4], m[7]}; + float b[3] = {m[2], m[5], m[8]}; + + JXL_RETURN_IF_ERROR(CreateICCXYZTag(r, &tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("rXYZ", tag_offset, tag_size, &tagtable, &offsets); + + JXL_RETURN_IF_ERROR(CreateICCXYZTag(g, &tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("gXYZ", tag_offset, tag_size, &tagtable, &offsets); + + JXL_RETURN_IF_ERROR(CreateICCXYZTag(b, &tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("bXYZ", tag_offset, tag_size, &tagtable, &offsets); + } + + if (c.color_space == JXL_COLOR_SPACE_XYB) { + JXL_RETURN_IF_ERROR(CreateICCLutAtoBTagForXYB(&tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("A2B0", tag_offset, tag_size, &tagtable, &offsets); + JXL_RETURN_IF_ERROR(CreateICCNoOpBToATag(&tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("B2A0", tag_offset, tag_size, &tagtable, &offsets); + } else if (kEnable3DToneMapping && CanToneMap(c)) { + JXL_RETURN_IF_ERROR(CreateICCLutAtoBTagForHDR(c, &tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("A2B0", tag_offset, tag_size, &tagtable, &offsets); + JXL_RETURN_IF_ERROR(CreateICCNoOpBToATag(&tags)); + FinalizeICCTag(&tags, &tag_offset, &tag_size); + AddToICCTagTable("B2A0", tag_offset, tag_size, &tagtable, &offsets); + } else { + if (tf == JXL_TRANSFER_FUNCTION_GAMMA) { + float gamma = 1.0 / c.gamma; + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag({gamma}, 0, &tags)); + } else if (c.color_space != JXL_COLOR_SPACE_XYB) { + switch (tf) { + case JXL_TRANSFER_FUNCTION_HLG: + CreateICCCurvCurvTag( + CreateTableCurve(64, ExtraTF::kHLG, CanToneMap(c)), &tags); + break; + case JXL_TRANSFER_FUNCTION_PQ: + CreateICCCurvCurvTag( + CreateTableCurve(64, ExtraTF::kPQ, CanToneMap(c)), &tags); + break; + case JXL_TRANSFER_FUNCTION_SRGB: + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag( + {2.4, 1.0 / 1.055, 0.055 / 1.055, 1.0 / 12.92, 0.04045}, 3, + &tags)); + break; + case JXL_TRANSFER_FUNCTION_709: + JXL_RETURN_IF_ERROR(CreateICCCurvParaTag( + {1.0 / 0.45, 1.0 / 1.099, 0.099 / 1.099, 1.0 / 4.5, 0.081}, 3, + &tags)); + break; + case JXL_TRANSFER_FUNCTION_LINEAR: + JXL_RETURN_IF_ERROR( + CreateICCCurvParaTag({1.0, 1.0, 0.0, 1.0, 0.0}, 3, &tags)); + break; + case JXL_TRANSFER_FUNCTION_DCI: + JXL_RETURN_IF_ERROR( + CreateICCCurvParaTag({2.6, 1.0, 0.0, 1.0, 0.0}, 3, &tags)); + break; + default: + JXL_UNREACHABLE("Unknown TF %u", static_cast<unsigned int>(tf)); + } + } + FinalizeICCTag(&tags, &tag_offset, &tag_size); + if (c.color_space == JXL_COLOR_SPACE_GRAY) { + AddToICCTagTable("kTRC", tag_offset, tag_size, &tagtable, &offsets); + } else { + AddToICCTagTable("rTRC", tag_offset, tag_size, &tagtable, &offsets); + AddToICCTagTable("gTRC", tag_offset, tag_size, &tagtable, &offsets); + AddToICCTagTable("bTRC", tag_offset, tag_size, &tagtable, &offsets); + } + } + + // Tag count + WriteICCUint32(offsets.size(), 0, &tagtable); + for (size_t i = 0; i < offsets.size(); i++) { + WriteICCUint32(offsets[i] + header.size() + tagtable.size(), 4 + 12 * i + 4, + &tagtable); + } + + // ICC profile size + WriteICCUint32(header.size() + tagtable.size() + tags.size(), 0, &header); + + *icc = header; + Bytes(tagtable).AppendTo(icc); + Bytes(tags).AppendTo(icc); + + // The MD5 checksum must be computed on the profile with profile flags, + // rendering intent, and region of the checksum itself, set to 0. + // TODO(lode): manually verify with a reliable tool that this creates correct + // signature (profile id) for ICC profiles. + std::vector<uint8_t> icc_sum = *icc; + if (icc_sum.size() >= 64 + 4) { + memset(icc_sum.data() + 44, 0, 4); + memset(icc_sum.data() + 64, 0, 4); + } + uint8_t checksum[16]; + detail::ICCComputeMD5(icc_sum, checksum); + + memcpy(icc->data() + 84, checksum, sizeof(checksum)); + + return true; +} + +} // namespace detail + +// Returns a representation of the ColorEncoding fields (not icc). +// Example description: "RGB_D65_SRG_Rel_Lin" +static JXL_MAYBE_UNUSED std::string ColorEncodingDescription( + const JxlColorEncoding& c) { + return detail::ColorEncodingDescriptionImpl(c); +} + +// NOTE: for XYB colorspace, the created profile can be used to transform a +// *scaled* XYB image (created by ScaleXYB()) to another colorspace. +static JXL_MAYBE_UNUSED Status MaybeCreateProfile(const JxlColorEncoding& c, + std::vector<uint8_t>* icc) { + return detail::MaybeCreateProfileImpl(c, icc); +} + +} // namespace jxl + +#endif // LIB_JXL_CMS_JXL_CMS_INTERNAL_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/opsin_params.h b/third_party/jpeg-xl/lib/jxl/cms/opsin_params.h new file mode 100644 index 0000000000..48e8e254f7 --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/opsin_params.h @@ -0,0 +1,160 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#ifndef LIB_JXL_CMS_OPSIN_PARAMS_H_ +#define LIB_JXL_CMS_OPSIN_PARAMS_H_ + +#include <array> + +// Constants that define the XYB color space. + +namespace jxl { +namespace cms { + +// Parameters for opsin absorbance. +constexpr float kM02 = 0.078f; +constexpr float kM00 = 0.30f; +constexpr float kM01 = 1.0f - kM02 - kM00; + +constexpr float kM12 = 0.078f; +constexpr float kM10 = 0.23f; +constexpr float kM11 = 1.0f - kM12 - kM10; + +constexpr float kM20 = 0.24342268924547819f; +constexpr float kM21 = 0.20476744424496821f; +constexpr float kM22 = 1.0f - kM20 - kM21; + +constexpr float kBScale = 1.0f; +constexpr float kYToBRatio = 1.0f; // works better with 0.50017729543783418 +constexpr float kBToYRatio = 1.0f / kYToBRatio; + +constexpr float kOpsinAbsorbanceBias0 = 0.0037930732552754493f; +constexpr float kOpsinAbsorbanceBias1 = kOpsinAbsorbanceBias0; +constexpr float kOpsinAbsorbanceBias2 = kOpsinAbsorbanceBias0; + +// Opsin absorbance matrix is now frozen. +constexpr std::array<float, 9> kOpsinAbsorbanceMatrix = { + kM00, kM01, kM02, kM10, kM11, kM12, kM20, kM21, kM22, +}; + +constexpr std::array<float, 9> kDefaultInverseOpsinAbsorbanceMatrix = { + 11.031566901960783f, -9.866943921568629f, -0.16462299647058826f, + -3.254147380392157f, 4.418770392156863f, -0.16462299647058826f, + -3.6588512862745097f, 2.7129230470588235f, 1.9459282392156863f}; + +// Must be the inverse matrix of kOpsinAbsorbanceMatrix and match the spec. +static inline const float* DefaultInverseOpsinAbsorbanceMatrix() { + return kDefaultInverseOpsinAbsorbanceMatrix.data(); +} + +constexpr std::array<float, 3> kOpsinAbsorbanceBias = { + kOpsinAbsorbanceBias0, + kOpsinAbsorbanceBias1, + kOpsinAbsorbanceBias2, +}; + +constexpr std::array<float, 4> kNegOpsinAbsorbanceBiasRGB = { + -kOpsinAbsorbanceBias0, -kOpsinAbsorbanceBias1, -kOpsinAbsorbanceBias2, + 1.0f}; + +constexpr float kScaledXYBOffset0 = 0.015386134f; +constexpr float kScaledXYBOffset1 = 0.0f; +constexpr float kScaledXYBOffset2 = 0.27770459f; + +constexpr std::array<float, 3> kScaledXYBOffset = { + kScaledXYBOffset0, kScaledXYBOffset1, kScaledXYBOffset2}; + +constexpr float kScaledXYBScale0 = 22.995788804f; +constexpr float kScaledXYBScale1 = 1.183000077f; +constexpr float kScaledXYBScale2 = 1.502141333f; + +constexpr std::array<float, 3> kScaledXYBScale = { + kScaledXYBScale0, + kScaledXYBScale1, + kScaledXYBScale2, +}; + +// NB(eustas): following function/variable names are just "namos". + +// More precise calculation of 1 / ((1 / r1) + (1 / r2)) +constexpr float ReciprocialSum(float r1, float r2) { + return (r1 * r2) / (r1 + r2); +} + +constexpr float kXYBOffset0 = kScaledXYBOffset0 + kScaledXYBOffset1; +constexpr float kXYBOffset1 = + kScaledXYBOffset1 - kScaledXYBOffset0 + (1.0f / kScaledXYBScale0); +constexpr float kXYBOffset2 = kScaledXYBOffset1 + kScaledXYBOffset2; + +constexpr std::array<float, 3> kXYBOffset = {kXYBOffset0, kXYBOffset1, + kXYBOffset2}; + +constexpr float kXYBScale0 = ReciprocialSum(kScaledXYBScale0, kScaledXYBScale1); +constexpr float kXYBScale1 = ReciprocialSum(kScaledXYBScale0, kScaledXYBScale1); +constexpr float kXYBScale2 = ReciprocialSum(kScaledXYBScale1, kScaledXYBScale2); + +constexpr std::array<float, 3> kXYBScale = {kXYBScale0, kXYBScale1, kXYBScale2}; + +template <size_t idx> +constexpr float ScaledXYBScale() { + return (idx == 0) ? kScaledXYBScale0 + : (idx == 1) ? kScaledXYBScale1 + : kScaledXYBScale2; +} + +template <size_t idx> +constexpr float ScaledXYBOffset() { + return (idx == 0) ? kScaledXYBOffset0 + : (idx == 1) ? kScaledXYBOffset1 + : kScaledXYBOffset2; +} + +template <size_t x, size_t y, size_t b, size_t idx> +constexpr float XYBCorner() { + return (((idx == 0) ? x + : (idx == 1) ? y + : b) / + ScaledXYBScale<idx>()) - + ScaledXYBOffset<idx>(); +} + +template <size_t x, size_t y, size_t b, size_t idx> +constexpr float ScaledA2BCorner() { + return (idx == 0) ? (XYBCorner<x, y, b, 1>() + XYBCorner<x, y, b, 0>()) + : (idx == 1) ? (XYBCorner<x, y, b, 1>() - XYBCorner<x, y, b, 0>()) + : (XYBCorner<x, y, b, 2>() + XYBCorner<x, y, b, 1>()); +} + +typedef std::array<float, 3> ColorCube0D; +template <size_t x, size_t y, size_t b> +constexpr ColorCube0D UnscaledA2BCorner() { + return {(ScaledA2BCorner<x, y, b, 0>() + kXYBOffset0) * kXYBScale0, + (ScaledA2BCorner<x, y, b, 1>() + kXYBOffset1) * kXYBScale1, + (ScaledA2BCorner<x, y, b, 2>() + kXYBOffset2) * kXYBScale2}; +} + +typedef std::array<ColorCube0D, 2> ColorCube1D; +template <size_t x, size_t y> +constexpr ColorCube1D UnscaledA2BCubeXY() { + return {UnscaledA2BCorner<x, y, 0>(), UnscaledA2BCorner<x, y, 1>()}; +} + +typedef std::array<ColorCube1D, 2> ColorCube2D; +template <size_t x> +constexpr ColorCube2D UnscaledA2BCubeX() { + return {UnscaledA2BCubeXY<x, 0>(), UnscaledA2BCubeXY<x, 1>()}; +} + +typedef std::array<ColorCube2D, 2> ColorCube3D; +constexpr ColorCube3D UnscaledA2BCube() { + return {UnscaledA2BCubeX<0>(), UnscaledA2BCubeX<1>()}; +} + +constexpr ColorCube3D kUnscaledA2BCube = UnscaledA2BCube(); + +} // namespace cms +} // namespace jxl + +#endif // LIB_JXL_CMS_OPSIN_PARAMS_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/tone_mapping-inl.h b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping-inl.h new file mode 100644 index 0000000000..3d94ccea12 --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping-inl.h @@ -0,0 +1,191 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#if defined(LIB_JXL_CMS_TONE_MAPPING_INL_H_) == defined(HWY_TARGET_TOGGLE) +#ifdef LIB_JXL_CMS_TONE_MAPPING_INL_H_ +#undef LIB_JXL_CMS_TONE_MAPPING_INL_H_ +#else +#define LIB_JXL_CMS_TONE_MAPPING_INL_H_ +#endif + +#include <hwy/highway.h> + +#include "lib/jxl/cms/tone_mapping.h" +#include "lib/jxl/cms/transfer_functions-inl.h" + +HWY_BEFORE_NAMESPACE(); +namespace jxl { +namespace HWY_NAMESPACE { +namespace { + +// These templates are not found via ADL. +using hwy::HWY_NAMESPACE::Clamp; +using hwy::HWY_NAMESPACE::Max; +using hwy::HWY_NAMESPACE::ZeroIfNegative; + +template <typename D> +class Rec2408ToneMapper : Rec2408ToneMapperBase { + private: + using V = hwy::HWY_NAMESPACE::Vec<D>; + + public: + using Rec2408ToneMapperBase::Rec2408ToneMapperBase; + + void ToneMap(V* red, V* green, V* blue) const { + const V luminance = Mul(Set(df_, source_range_.second), + (MulAdd(Set(df_, red_Y_), *red, + MulAdd(Set(df_, green_Y_), *green, + Mul(Set(df_, blue_Y_), *blue))))); + const V pq_mastering_min = Set(df_, pq_mastering_min_); + const V inv_pq_mastering_range = Set(df_, inv_pq_mastering_range_); + const V normalized_pq = Min( + Set(df_, 1.f), + Mul(Sub(InvEOTF(luminance), pq_mastering_min), inv_pq_mastering_range)); + const V ks = Set(df_, ks_); + const V e2 = + IfThenElse(Lt(normalized_pq, ks), normalized_pq, P(normalized_pq)); + const V one_minus_e2 = Sub(Set(df_, 1), e2); + const V one_minus_e2_2 = Mul(one_minus_e2, one_minus_e2); + const V one_minus_e2_4 = Mul(one_minus_e2_2, one_minus_e2_2); + const V b = Set(df_, min_lum_); + const V e3 = MulAdd(b, one_minus_e2_4, e2); + const V pq_mastering_range = Set(df_, pq_mastering_range_); + const V e4 = MulAdd(e3, pq_mastering_range, pq_mastering_min); + const V new_luminance = + Min(Set(df_, target_range_.second), + ZeroIfNegative(tf_pq_.DisplayFromEncoded(df_, e4))); + const V min_luminance = Set(df_, 1e-6f); + const auto use_cap = Le(luminance, min_luminance); + const V ratio = Div(new_luminance, Max(luminance, min_luminance)); + const V cap = Mul(new_luminance, Set(df_, inv_target_peak_)); + const V normalizer = Set(df_, normalizer_); + const V multiplier = Mul(ratio, normalizer); + for (V* const val : {red, green, blue}) { + *val = IfThenElse(use_cap, cap, Mul(*val, multiplier)); + } + } + + private: + V InvEOTF(const V luminance) const { + return tf_pq_.EncodedFromDisplay(df_, luminance); + } + V T(const V a) const { + const V ks = Set(df_, ks_); + const V inv_one_minus_ks = Set(df_, inv_one_minus_ks_); + return Mul(Sub(a, ks), inv_one_minus_ks); + } + V P(const V b) const { + const V t_b = T(b); + const V t_b_2 = Mul(t_b, t_b); + const V t_b_3 = Mul(t_b_2, t_b); + const V ks = Set(df_, ks_); + const V max_lum = Set(df_, max_lum_); + return MulAdd( + MulAdd(Set(df_, 2), t_b_3, MulAdd(Set(df_, -3), t_b_2, Set(df_, 1))), + ks, + MulAdd(Add(t_b_3, MulAdd(Set(df_, -2), t_b_2, t_b)), + Sub(Set(df_, 1), ks), + Mul(MulAdd(Set(df_, -2), t_b_3, Mul(Set(df_, 3), t_b_2)), + max_lum))); + } + + D df_; + const TF_PQ tf_pq_ = TF_PQ(/*display_intensity_target=*/1.0); +}; + +class HlgOOTF : HlgOOTF_Base { + public: + using HlgOOTF_Base::HlgOOTF_Base; + + static HlgOOTF FromSceneLight(float display_luminance, + const float primaries_luminances[3]) { + return HlgOOTF(/*gamma=*/1.2f * + std::pow(1.111f, std::log2(display_luminance / 1000.f)), + primaries_luminances); + } + + static HlgOOTF ToSceneLight(float display_luminance, + const float primaries_luminances[3]) { + return HlgOOTF( + /*gamma=*/(1 / 1.2f) * + std::pow(1.111f, -std::log2(display_luminance / 1000.f)), + primaries_luminances); + } + + template <typename V> + void Apply(V* red, V* green, V* blue) const { + hwy::HWY_NAMESPACE::DFromV<V> df; + if (!apply_ootf_) return; + const V luminance = + MulAdd(Set(df, red_Y_), *red, + MulAdd(Set(df, green_Y_), *green, Mul(Set(df, blue_Y_), *blue))); + const V ratio = + Min(FastPowf(df, luminance, Set(df, exponent_)), Set(df, 1e9)); + *red = Mul(*red, ratio); + *green = Mul(*green, ratio); + *blue = Mul(*blue, ratio); + } + + bool WarrantsGamutMapping() const { return apply_ootf_ && exponent_ < 0; } +}; + +template <typename V> +void GamutMap(V* red, V* green, V* blue, const float primaries_luminances[3], + float preserve_saturation = 0.1f) { + hwy::HWY_NAMESPACE::DFromV<V> df; + const V luminance = + MulAdd(Set(df, primaries_luminances[0]), *red, + MulAdd(Set(df, primaries_luminances[1]), *green, + Mul(Set(df, primaries_luminances[2]), *blue))); + + // Desaturate out-of-gamut pixels. This is done by mixing each pixel + // with just enough gray of the target luminance to make all + // components non-negative. + // - For saturation preservation, if a component is still larger than + // 1 then the pixel is normalized to have a maximum component of 1. + // That will reduce its luminance. + // - For luminance preservation, getting all components below 1 is + // done by mixing in yet more gray. That will desaturate it further. + const V zero = Zero(df); + const V one = Set(df, 1); + V gray_mix_saturation = zero; + V gray_mix_luminance = zero; + for (const V* ch : {red, green, blue}) { + const V& val = *ch; + const V val_minus_gray = Sub(val, luminance); + const V inv_val_minus_gray = + Div(one, IfThenElse(Eq(val_minus_gray, zero), one, val_minus_gray)); + const V val_over_val_minus_gray = Mul(val, inv_val_minus_gray); + gray_mix_saturation = + IfThenElse(Ge(val_minus_gray, zero), gray_mix_saturation, + Max(gray_mix_saturation, val_over_val_minus_gray)); + gray_mix_luminance = + Max(gray_mix_luminance, + IfThenElse(Le(val_minus_gray, zero), gray_mix_saturation, + Sub(val_over_val_minus_gray, inv_val_minus_gray))); + } + const V gray_mix = Clamp( + MulAdd(Set(df, preserve_saturation), + Sub(gray_mix_saturation, gray_mix_luminance), gray_mix_luminance), + zero, one); + for (V* const ch : {red, green, blue}) { + V& val = *ch; + val = MulAdd(gray_mix, Sub(luminance, val), val); + } + const V max_clr = Max(Max(one, *red), Max(*green, *blue)); + const V normalizer = Div(one, max_clr); + for (V* const ch : {red, green, blue}) { + V& val = *ch; + val = Mul(val, normalizer); + } +} + +} // namespace +// NOLINTNEXTLINE(google-readability-namespace-comments) +} // namespace HWY_NAMESPACE +} // namespace jxl +HWY_AFTER_NAMESPACE(); + +#endif // LIB_JXL_CMS_TONE_MAPPING_INL_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/tone_mapping.h b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping.h new file mode 100644 index 0000000000..a114109ea6 --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping.h @@ -0,0 +1,179 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#ifndef LIB_JXL_CMS_TONE_MAPPING_H_ +#define LIB_JXL_CMS_TONE_MAPPING_H_ + +#include <algorithm> +#include <cmath> +#include <utility> + +#include "lib/jxl/base/common.h" +#include "lib/jxl/base/compiler_specific.h" +#include "lib/jxl/cms/transfer_functions.h" + +namespace jxl { + +class Rec2408ToneMapperBase { + public: + explicit Rec2408ToneMapperBase(std::pair<float, float> source_range, + std::pair<float, float> target_range, + const float primaries_luminances[3]) + : source_range_(source_range), + target_range_(target_range), + red_Y_(primaries_luminances[0]), + green_Y_(primaries_luminances[1]), + blue_Y_(primaries_luminances[2]) {} + + // TODO(eustas): test me + void ToneMap(float* red, float* green, float* blue) const { + const float luminance = + source_range_.second * + (red_Y_ * *red + green_Y_ * *green + blue_Y_ * *blue); + const float normalized_pq = + std::min(1.f, (InvEOTF(luminance) - pq_mastering_min_) * + inv_pq_mastering_range_); + const float e2 = (normalized_pq < ks_) ? normalized_pq : P(normalized_pq); + const float one_minus_e2 = 1 - e2; + const float one_minus_e2_2 = one_minus_e2 * one_minus_e2; + const float one_minus_e2_4 = one_minus_e2_2 * one_minus_e2_2; + const float e3 = min_lum_ * one_minus_e2_4 + e2; + const float e4 = e3 * pq_mastering_range_ + pq_mastering_min_; + const float d4 = + TF_PQ_Base::DisplayFromEncoded(/*display_intensity_target=*/1.0, e4); + const float new_luminance = Clamp1(d4, 0.f, target_range_.second); + const float min_luminance = 1e-6f; + const bool use_cap = (luminance <= min_luminance); + const float ratio = new_luminance / std::max(luminance, min_luminance); + const float cap = new_luminance * inv_target_peak_; + const float multiplier = ratio * normalizer_; + for (float* const val : {red, green, blue}) { + *val = use_cap ? cap : *val * multiplier; + } + } + + protected: + float InvEOTF(const float luminance) const { + return TF_PQ_Base::EncodedFromDisplay(/*display_intensity_target=*/1.0, + luminance); + } + float T(const float a) const { return (a - ks_) * inv_one_minus_ks_; } + float P(const float b) const { + const float t_b = T(b); + const float t_b_2 = t_b * t_b; + const float t_b_3 = t_b_2 * t_b; + return (2 * t_b_3 - 3 * t_b_2 + 1) * ks_ + + (t_b_3 - 2 * t_b_2 + t_b) * (1 - ks_) + + (-2 * t_b_3 + 3 * t_b_2) * max_lum_; + } + + const std::pair<float, float> source_range_; + const std::pair<float, float> target_range_; + const float red_Y_; + const float green_Y_; + const float blue_Y_; + + const float pq_mastering_min_ = InvEOTF(source_range_.first); + const float pq_mastering_max_ = InvEOTF(source_range_.second); + const float pq_mastering_range_ = pq_mastering_max_ - pq_mastering_min_; + const float inv_pq_mastering_range_ = 1.0f / pq_mastering_range_; + // TODO(eustas): divide instead of inverse-multiply? + const float min_lum_ = (InvEOTF(target_range_.first) - pq_mastering_min_) * + inv_pq_mastering_range_; + // TODO(eustas): divide instead of inverse-multiply? + const float max_lum_ = (InvEOTF(target_range_.second) - pq_mastering_min_) * + inv_pq_mastering_range_; + const float ks_ = 1.5f * max_lum_ - 0.5f; + + const float inv_one_minus_ks_ = 1.0f / std::max(1e-6f, 1.0f - ks_); + + const float normalizer_ = source_range_.second / target_range_.second; + const float inv_target_peak_ = 1.f / target_range_.second; +}; + +class HlgOOTF_Base { + public: + explicit HlgOOTF_Base(float source_luminance, float target_luminance, + const float primaries_luminances[3]) + : HlgOOTF_Base(/*gamma=*/std::pow(1.111f, std::log2(target_luminance / + source_luminance)), + primaries_luminances) {} + + // TODO(eustas): test me + void Apply(float* red, float* green, float* blue) const { + if (!apply_ootf_) return; + const float luminance = red_Y_ * *red + green_Y_ * *green + blue_Y_ * *blue; + const float ratio = std::min<float>(powf(luminance, exponent_), 1e9); + *red *= ratio; + *green *= ratio; + *blue *= ratio; + } + + protected: + explicit HlgOOTF_Base(float gamma, const float luminances[3]) + : exponent_(gamma - 1), + red_Y_(luminances[0]), + green_Y_(luminances[1]), + blue_Y_(luminances[2]) {} + const float exponent_; + const bool apply_ootf_ = exponent_ < -0.01f || 0.01f < exponent_; + const float red_Y_; + const float green_Y_; + const float blue_Y_; +}; + +static JXL_MAYBE_UNUSED void GamutMapScalar(float* red, float* green, + float* blue, + const float primaries_luminances[3], + float preserve_saturation = 0.1f) { + const float luminance = primaries_luminances[0] * *red + + primaries_luminances[1] * *green + + primaries_luminances[2] * *blue; + + // Desaturate out-of-gamut pixels. This is done by mixing each pixel + // with just enough gray of the target luminance to make all + // components non-negative. + // - For saturation preservation, if a component is still larger than + // 1 then the pixel is normalized to have a maximum component of 1. + // That will reduce its luminance. + // - For luminance preservation, getting all components below 1 is + // done by mixing in yet more gray. That will desaturate it further. + float gray_mix_saturation = 0.0f; + float gray_mix_luminance = 0.0f; + for (const float* ch : {red, green, blue}) { + const float& val = *ch; + const float val_minus_gray = val - luminance; + const float inv_val_minus_gray = + 1.0f / ((val_minus_gray == 0.0f) ? 1.0f : val_minus_gray); + const float val_over_val_minus_gray = val * inv_val_minus_gray; + gray_mix_saturation = + (val_minus_gray >= 0.0f) + ? gray_mix_saturation + : std::max(gray_mix_saturation, val_over_val_minus_gray); + gray_mix_luminance = + std::max(gray_mix_luminance, + (val_minus_gray <= 0.0f) + ? gray_mix_saturation + : (val_over_val_minus_gray - inv_val_minus_gray)); + } + const float gray_mix = + Clamp1((preserve_saturation * (gray_mix_saturation - gray_mix_luminance) + + gray_mix_luminance), + 0.0f, 1.0f); + for (float* const ch : {red, green, blue}) { + float& val = *ch; + val = gray_mix * (luminance - val) + val; + } + const float max_clr = std::max({1.0f, *red, *green, *blue}); + const float normalizer = 1.0f / max_clr; + for (float* const ch : {red, green, blue}) { + float& val = *ch; + val *= normalizer; + } +} + +} // namespace jxl + +#endif // LIB_JXL_CMS_TONE_MAPPING_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/tone_mapping_test.cc b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping_test.cc new file mode 100644 index 0000000000..dda2bbb0aa --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/tone_mapping_test.cc @@ -0,0 +1,147 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "lib/jxl/cms/tone_mapping_test.cc" +#include "lib/jxl/cms/tone_mapping.h" + +#include <cstdio> +#include <hwy/foreach_target.h> + +#include "lib/jxl/base/random.h" +#include "lib/jxl/cms/tone_mapping-inl.h" +#include "lib/jxl/testing.h" + +// Test utils +#include <hwy/highway.h> +#include <hwy/tests/hwy_gtest.h> +HWY_BEFORE_NAMESPACE(); +namespace jxl { +namespace HWY_NAMESPACE { +namespace { + +HWY_NOINLINE void TestRec2408ToneMap() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + float src = 11000.0 + rng.UniformF(-150.0f, 150.0f); + float tgt = 250 + rng.UniformF(-5.0f, 5.0f); + float luminances[3] = {rng.UniformF(0.2f, 0.4f), rng.UniformF(0.2f, 0.4f), + rng.UniformF(0.2f, 0.4f)}; + float rgb[3] = {rng.UniformF(0.0f, 1.0f), rng.UniformF(0.0f, 1.0f), + rng.UniformF(0.0f, 1.0f)}; + Rec2408ToneMapper<decltype(d)> tone_mapper({0, src}, {0, tgt}, luminances); + auto r = Set(d, rgb[0]); + auto g = Set(d, rgb[1]); + auto b = Set(d, rgb[2]); + tone_mapper.ToneMap(&r, &g, &b); + Rec2408ToneMapperBase tone_mapper_base({0, src}, {0, tgt}, luminances); + tone_mapper_base.ToneMap(&rgb[0], &rgb[1], &rgb[2]); + const float actual_r = GetLane(r); + const float expected_r = rgb[0]; + const float abs_err_r = std::abs(expected_r - actual_r); + EXPECT_LT(abs_err_r, 2.75e-5); + const float actual_g = GetLane(g); + const float expected_g = rgb[1]; + const float abs_err_g = std::abs(expected_g - actual_g); + EXPECT_LT(abs_err_g, 2.75e-5); + const float actual_b = GetLane(b); + const float expected_b = rgb[2]; + const float abs_err_b = std::abs(expected_b - actual_b); + EXPECT_LT(abs_err_b, 2.75e-5); + max_abs_err = std::max({max_abs_err, abs_err_r, abs_err_g, abs_err_b}); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +HWY_NOINLINE void TestHlgOotfApply() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + float src = 300.0 + rng.UniformF(-50.0f, 50.0f); + float tgt = 80 + rng.UniformF(-5.0f, 5.0f); + float luminances[3] = {rng.UniformF(0.2f, 0.4f), rng.UniformF(0.2f, 0.4f), + rng.UniformF(0.2f, 0.4f)}; + float rgb[3] = {rng.UniformF(0.0f, 1.0f), rng.UniformF(0.0f, 1.0f), + rng.UniformF(0.0f, 1.0f)}; + HlgOOTF ootf(src, tgt, luminances); + auto r = Set(d, rgb[0]); + auto g = Set(d, rgb[1]); + auto b = Set(d, rgb[2]); + ootf.Apply(&r, &g, &b); + HlgOOTF_Base ootf_base(src, tgt, luminances); + ootf_base.Apply(&rgb[0], &rgb[1], &rgb[2]); + const float actual_r = GetLane(r); + const float expected_r = rgb[0]; + const float abs_err_r = std::abs(expected_r - actual_r); + EXPECT_LT(abs_err_r, 7.2e-7); + const float actual_g = GetLane(g); + const float expected_g = rgb[1]; + const float abs_err_g = std::abs(expected_g - actual_g); + EXPECT_LT(abs_err_g, 7.2e-7); + const float actual_b = GetLane(b); + const float expected_b = rgb[2]; + const float abs_err_b = std::abs(expected_b - actual_b); + EXPECT_LT(abs_err_b, 7.2e-7); + max_abs_err = std::max({max_abs_err, abs_err_r, abs_err_g, abs_err_b}); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +HWY_NOINLINE void TestGamutMap() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + float preserve_saturation = rng.UniformF(0.2f, 0.4f); + float luminances[3] = {rng.UniformF(0.2f, 0.4f), rng.UniformF(0.2f, 0.4f), + rng.UniformF(0.2f, 0.4f)}; + float rgb[3] = {rng.UniformF(0.0f, 1.0f), rng.UniformF(0.0f, 1.0f), + rng.UniformF(0.0f, 1.0f)}; + auto r = Set(d, rgb[0]); + auto g = Set(d, rgb[1]); + auto b = Set(d, rgb[2]); + GamutMap(&r, &g, &b, luminances, preserve_saturation); + GamutMapScalar(&rgb[0], &rgb[1], &rgb[2], luminances, preserve_saturation); + const float actual_r = GetLane(r); + const float expected_r = rgb[0]; + const float abs_err_r = std::abs(expected_r - actual_r); + EXPECT_LT(abs_err_r, 1e-10); + const float actual_g = GetLane(g); + const float expected_g = rgb[1]; + const float abs_err_g = std::abs(expected_g - actual_g); + EXPECT_LT(abs_err_g, 1e-10); + const float actual_b = GetLane(b); + const float expected_b = rgb[2]; + const float abs_err_b = std::abs(expected_b - actual_b); + EXPECT_LT(abs_err_b, 1e-10); + max_abs_err = std::max({max_abs_err, abs_err_r, abs_err_g, abs_err_b}); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +} // namespace +// NOLINTNEXTLINE(google-readability-namespace-comments) +} // namespace HWY_NAMESPACE +} // namespace jxl +HWY_AFTER_NAMESPACE(); + +#if HWY_ONCE +namespace jxl { + +class ToneMappingTargetTest : public hwy::TestWithParamTarget {}; +HWY_TARGET_INSTANTIATE_TEST_SUITE_P(ToneMappingTargetTest); + +HWY_EXPORT_AND_TEST_P(ToneMappingTargetTest, TestRec2408ToneMap); +HWY_EXPORT_AND_TEST_P(ToneMappingTargetTest, TestHlgOotfApply); +HWY_EXPORT_AND_TEST_P(ToneMappingTargetTest, TestGamutMap); + +} // namespace jxl +#endif // HWY_ONCE diff --git a/third_party/jpeg-xl/lib/jxl/cms/transfer_functions-inl.h b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions-inl.h new file mode 100644 index 0000000000..84bcbb45ed --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions-inl.h @@ -0,0 +1,334 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// Transfer functions for color encodings. + +#if defined(LIB_JXL_CMS_TRANSFER_FUNCTIONS_INL_H_) == defined(HWY_TARGET_TOGGLE) +#ifdef LIB_JXL_CMS_TRANSFER_FUNCTIONS_INL_H_ +#undef LIB_JXL_CMS_TRANSFER_FUNCTIONS_INL_H_ +#else +#define LIB_JXL_CMS_TRANSFER_FUNCTIONS_INL_H_ +#endif + +#include <algorithm> +#include <cmath> +#include <hwy/highway.h> + +#include "lib/jxl/base/compiler_specific.h" +#include "lib/jxl/base/fast_math-inl.h" +#include "lib/jxl/base/rational_polynomial-inl.h" +#include "lib/jxl/base/status.h" +#include "lib/jxl/cms/transfer_functions.h" + +HWY_BEFORE_NAMESPACE(); +namespace jxl { +namespace HWY_NAMESPACE { + +// These templates are not found via ADL. +using hwy::HWY_NAMESPACE::And; +using hwy::HWY_NAMESPACE::AndNot; +using hwy::HWY_NAMESPACE::Gt; +using hwy::HWY_NAMESPACE::IfThenElse; +using hwy::HWY_NAMESPACE::Lt; +using hwy::HWY_NAMESPACE::Or; +using hwy::HWY_NAMESPACE::Sqrt; +using hwy::HWY_NAMESPACE::TableLookupBytes; + +// Definitions for BT.2100-2 transfer functions (used inside/outside SIMD): +// "display" is linear light (nits) normalized to [0, 1]. +// "encoded" is a nonlinear encoding (e.g. PQ) in [0, 1]. +// "scene" is a linear function of photon counts, normalized to [0, 1]. + +// Despite the stated ranges, we need unbounded transfer functions: see +// http://www.littlecms.com/CIC18_UnboundedCMM.pdf. Inputs can be negative or +// above 1 due to chromatic adaptation. To avoid severe round-trip errors caused +// by clamping, we mirror negative inputs via copysign (f(-x) = -f(x), see +// https://developer.apple.com/documentation/coregraphics/cgcolorspace/1644735-extendedsrgb) +// and extend the function domains above 1. + +// Hybrid Log-Gamma. +class TF_HLG : TF_HLG_Base { + public: + // Maximum error 5e-7. + template <class D, class V> + JXL_INLINE V EncodedFromDisplay(D d, V x) const { + const hwy::HWY_NAMESPACE::Rebind<uint32_t, D> du; + const V kSign = BitCast(d, Set(du, 0x80000000u)); + const V original_sign = And(x, kSign); + x = AndNot(kSign, x); // abs + const V below_div12 = Sqrt(Mul(Set(d, 3.0f), x)); + const V e = + MulAdd(Set(d, kA * 0.693147181f), + FastLog2f(d, MulAdd(Set(d, 12), x, Set(d, -kB))), Set(d, kC)); + const V magnitude = IfThenElse(Le(x, Set(d, kDiv12)), below_div12, e); + return Or(AndNot(kSign, magnitude), original_sign); + } +}; + +class TF_709 { + public: + JXL_INLINE double EncodedFromDisplay(const double d) const { + if (d < kThresh) return kMulLow * d; + return kMulHi * std::pow(d, kPowHi) + kSub; + } + + // Maximum error 1e-6. + template <class D, class V> + JXL_INLINE V EncodedFromDisplay(D d, V x) const { + auto low = Mul(Set(d, kMulLow), x); + auto hi = + MulAdd(Set(d, kMulHi), FastPowf(d, x, Set(d, kPowHi)), Set(d, kSub)); + return IfThenElse(Le(x, Set(d, kThresh)), low, hi); + } + + template <class D, class V> + JXL_INLINE V DisplayFromEncoded(D d, V x) const { + auto low = Mul(Set(d, kInvMulLow), x); + auto hi = FastPowf(d, MulAdd(x, Set(d, kInvMulHi), Set(d, kInvAdd)), + Set(d, kInvPowHi)); + return IfThenElse(Lt(x, Set(d, kInvThresh)), low, hi); + } + + private: + static constexpr double kThresh = 0.018; + static constexpr double kMulLow = 4.5; + static constexpr double kMulHi = 1.099; + static constexpr double kPowHi = 0.45; + static constexpr double kSub = -0.099; + + static constexpr double kInvThresh = 0.081; + static constexpr double kInvMulLow = 1 / 4.5; + static constexpr double kInvMulHi = 1 / 1.099; + static constexpr double kInvPowHi = 1 / 0.45; + static constexpr double kInvAdd = 0.099 * kInvMulHi; +}; + +// Perceptual Quantization +class TF_PQ : TF_PQ_Base { + public: + explicit TF_PQ(float display_intensity_target = kDefaultIntensityTarget) + : display_scaling_factor_to_10000_nits_(display_intensity_target * + (1.0f / 10000.0f)), + display_scaling_factor_from_10000_nits_(10000.0f / + display_intensity_target) {} + + // Maximum error 3e-6 + template <class D, class V> + JXL_INLINE V DisplayFromEncoded(D d, V x) const { + const hwy::HWY_NAMESPACE::Rebind<uint32_t, D> du; + const V kSign = BitCast(d, Set(du, 0x80000000u)); + const V original_sign = And(x, kSign); + x = AndNot(kSign, x); // abs + // 4-over-4-degree rational polynomial approximation on x+x*x. This improves + // the maximum error by about 5x over a rational polynomial for x. + auto xpxx = MulAdd(x, x, x); + HWY_ALIGN constexpr float p[(4 + 1) * 4] = { + HWY_REP4(2.62975656e-04f), HWY_REP4(-6.23553089e-03f), + HWY_REP4(7.38602301e-01f), HWY_REP4(2.64553172e+00f), + HWY_REP4(5.50034862e-01f), + }; + HWY_ALIGN constexpr float q[(4 + 1) * 4] = { + HWY_REP4(4.21350107e+02f), HWY_REP4(-4.28736818e+02f), + HWY_REP4(1.74364667e+02f), HWY_REP4(-3.39078883e+01f), + HWY_REP4(2.67718770e+00f), + }; + auto magnitude = EvalRationalPolynomial(d, xpxx, p, q); + return Or( + AndNot(kSign, + Mul(magnitude, Set(d, display_scaling_factor_from_10000_nits_))), + original_sign); + } + + // Maximum error 7e-7. + template <class D, class V> + JXL_INLINE V EncodedFromDisplay(D d, V x) const { + const hwy::HWY_NAMESPACE::Rebind<uint32_t, D> du; + const V kSign = BitCast(d, Set(du, 0x80000000u)); + const V original_sign = And(x, kSign); + x = AndNot(kSign, x); // abs + // 4-over-4-degree rational polynomial approximation on x**0.25, with two + // different polynomials above and below 1e-4. + auto xto025 = + Sqrt(Sqrt(Mul(x, Set(d, display_scaling_factor_to_10000_nits_)))); + HWY_ALIGN constexpr float p[(4 + 1) * 4] = { + HWY_REP4(1.351392e-02f), HWY_REP4(-1.095778e+00f), + HWY_REP4(5.522776e+01f), HWY_REP4(1.492516e+02f), + HWY_REP4(4.838434e+01f), + }; + HWY_ALIGN constexpr float q[(4 + 1) * 4] = { + HWY_REP4(1.012416e+00f), HWY_REP4(2.016708e+01f), + HWY_REP4(9.263710e+01f), HWY_REP4(1.120607e+02f), + HWY_REP4(2.590418e+01f), + }; + + HWY_ALIGN constexpr float plo[(4 + 1) * 4] = { + HWY_REP4(9.863406e-06f), HWY_REP4(3.881234e-01f), + HWY_REP4(1.352821e+02f), HWY_REP4(6.889862e+04f), + HWY_REP4(-2.864824e+05f), + }; + HWY_ALIGN constexpr float qlo[(4 + 1) * 4] = { + HWY_REP4(3.371868e+01f), HWY_REP4(1.477719e+03f), + HWY_REP4(1.608477e+04f), HWY_REP4(-4.389884e+04f), + HWY_REP4(-2.072546e+05f), + }; + + auto magnitude = IfThenElse(Lt(x, Set(d, 1e-4f)), + EvalRationalPolynomial(d, xto025, plo, qlo), + EvalRationalPolynomial(d, xto025, p, q)); + return Or(AndNot(kSign, magnitude), original_sign); + } + + private: + const float display_scaling_factor_to_10000_nits_; + const float display_scaling_factor_from_10000_nits_; +}; + +// sRGB +class TF_SRGB { + public: + template <typename V> + JXL_INLINE V DisplayFromEncoded(V x) const { + const HWY_FULL(float) d; + const HWY_FULL(uint32_t) du; + const V kSign = BitCast(d, Set(du, 0x80000000u)); + const V original_sign = And(x, kSign); + x = AndNot(kSign, x); // abs + + // TODO(janwas): range reduction + // Computed via af_cheb_rational (k=100); replicated 4x. + HWY_ALIGN constexpr float p[(4 + 1) * 4] = { + 2.200248328e-04f, 2.200248328e-04f, 2.200248328e-04f, 2.200248328e-04f, + 1.043637593e-02f, 1.043637593e-02f, 1.043637593e-02f, 1.043637593e-02f, + 1.624820318e-01f, 1.624820318e-01f, 1.624820318e-01f, 1.624820318e-01f, + 7.961564959e-01f, 7.961564959e-01f, 7.961564959e-01f, 7.961564959e-01f, + 8.210152774e-01f, 8.210152774e-01f, 8.210152774e-01f, 8.210152774e-01f, + }; + HWY_ALIGN constexpr float q[(4 + 1) * 4] = { + 2.631846970e-01f, 2.631846970e-01f, 2.631846970e-01f, + 2.631846970e-01f, 1.076976492e+00f, 1.076976492e+00f, + 1.076976492e+00f, 1.076976492e+00f, 4.987528350e-01f, + 4.987528350e-01f, 4.987528350e-01f, 4.987528350e-01f, + -5.512498495e-02f, -5.512498495e-02f, -5.512498495e-02f, + -5.512498495e-02f, 6.521209011e-03f, 6.521209011e-03f, + 6.521209011e-03f, 6.521209011e-03f, + }; + const V linear = Mul(x, Set(d, kLowDivInv)); + const V poly = EvalRationalPolynomial(d, x, p, q); + const V magnitude = + IfThenElse(Gt(x, Set(d, kThreshSRGBToLinear)), poly, linear); + return Or(AndNot(kSign, magnitude), original_sign); + } + + // Error ~5e-07 + template <class D, class V> + JXL_INLINE V EncodedFromDisplay(D d, V x) const { + const hwy::HWY_NAMESPACE::Rebind<uint32_t, D> du; + const V kSign = BitCast(d, Set(du, 0x80000000u)); + const V original_sign = And(x, kSign); + x = AndNot(kSign, x); // abs + + // Computed via af_cheb_rational (k=100); replicated 4x. + HWY_ALIGN constexpr float p[(4 + 1) * 4] = { + -5.135152395e-04f, -5.135152395e-04f, -5.135152395e-04f, + -5.135152395e-04f, 5.287254571e-03f, 5.287254571e-03f, + 5.287254571e-03f, 5.287254571e-03f, 3.903842876e-01f, + 3.903842876e-01f, 3.903842876e-01f, 3.903842876e-01f, + 1.474205315e+00f, 1.474205315e+00f, 1.474205315e+00f, + 1.474205315e+00f, 7.352629620e-01f, 7.352629620e-01f, + 7.352629620e-01f, 7.352629620e-01f, + }; + HWY_ALIGN constexpr float q[(4 + 1) * 4] = { + 1.004519624e-02f, 1.004519624e-02f, 1.004519624e-02f, 1.004519624e-02f, + 3.036675394e-01f, 3.036675394e-01f, 3.036675394e-01f, 3.036675394e-01f, + 1.340816930e+00f, 1.340816930e+00f, 1.340816930e+00f, 1.340816930e+00f, + 9.258482155e-01f, 9.258482155e-01f, 9.258482155e-01f, 9.258482155e-01f, + 2.424867759e-02f, 2.424867759e-02f, 2.424867759e-02f, 2.424867759e-02f, + }; + const V linear = Mul(x, Set(d, kLowDiv)); + const V poly = EvalRationalPolynomial(d, Sqrt(x), p, q); + const V magnitude = + IfThenElse(Gt(x, Set(d, kThreshLinearToSRGB)), poly, linear); + return Or(AndNot(kSign, magnitude), original_sign); + } + + private: + static constexpr float kThreshSRGBToLinear = 0.04045f; + static constexpr float kThreshLinearToSRGB = 0.0031308f; + static constexpr float kLowDiv = 12.92f; + static constexpr float kLowDivInv = 1.0f / kLowDiv; +}; + +// Linear to sRGB conversion with error of at most 1.2e-4. +template <typename D, typename V> +V FastLinearToSRGB(D d, V v) { + const hwy::HWY_NAMESPACE::Rebind<uint32_t, D> du; + const hwy::HWY_NAMESPACE::Rebind<int32_t, D> di; + // Convert to 0.25 - 0.5 range. + auto v025_05 = BitCast( + d, And(Or(BitCast(du, v), Set(du, 0x3e800000)), Set(du, 0x3effffff))); + // third degree polynomial approximation between 0.25 and 0.5 + // of 1.055/2^(7/2.4) * x^(1/2.4) * 0.5. A degree 4 polynomial only improves + // accuracy by about 3x. + auto d1 = MulAdd(v025_05, Set(d, 0.059914046f), Set(d, -0.108894556f)); + auto d2 = MulAdd(d1, v025_05, Set(d, 0.107963754f)); + auto pow = MulAdd(d2, v025_05, Set(d, 0.018092343f)); + // Compute extra multiplier depending on exponent. Valid exponent range for + // [0.0031308f, 1.0) is 0...8 after subtracting 118. + // The next three constants contain a representation of the powers of + // 2**(1/2.4) = 2**(5/12) times two; in particular, bits from 26 to 31 are + // always the same and in k2to512powers_basebits, and the two arrays contain + // the next groups of 8 bits. This ends up being a 22-bit representation (with + // a mantissa of 13 bits). The choice of polynomial to approximate is such + // that the multiplication factor has the highest 5 bits constant, and that + // the factor for the lowest possible exponent is a power of two (thus making + // the additional bits 0, which is used to correctly merge back together the + // floats). + constexpr uint32_t k2to512powers_basebits = 0x40000000; + HWY_ALIGN constexpr uint8_t k2to512powers_25to18bits[16] = { + 0x0, 0xa, 0x19, 0x26, 0x32, 0x41, 0x4d, 0x5c, + 0x68, 0x75, 0x83, 0x8f, 0xa0, 0xaa, 0xb9, 0xc6, + }; + HWY_ALIGN constexpr uint8_t k2to512powers_17to10bits[16] = { + 0x0, 0xb7, 0x4, 0xd, 0xcb, 0xe7, 0x41, 0x68, + 0x51, 0xd1, 0xeb, 0xf2, 0x0, 0xb7, 0x4, 0xd, + }; + // Note that vld1q_s8_x2 on ARM seems to actually be slower. +#if HWY_TARGET != HWY_SCALAR + using hwy::HWY_NAMESPACE::ShiftLeft; + using hwy::HWY_NAMESPACE::ShiftRight; + // Every lane of exp is now (if cast to byte) {0, 0, 0, <index for lookup>}. + auto exp = Sub(ShiftRight<23>(BitCast(di, v)), Set(di, 118)); + auto pow25to18bits = TableLookupBytes( + LoadDup128(di, + reinterpret_cast<const int32_t*>(k2to512powers_25to18bits)), + exp); + auto pow17to10bits = TableLookupBytes( + LoadDup128(di, + reinterpret_cast<const int32_t*>(k2to512powers_17to10bits)), + exp); + // Now, pow* contain {0, 0, 0, <part of float repr of multiplier>}. Here + // we take advantage of the fact that each table has its position 0 equal to + // 0. + // We can now just reassemble the float. + auto mul = BitCast( + d, Or(Or(ShiftLeft<18>(pow25to18bits), ShiftLeft<10>(pow17to10bits)), + Set(di, k2to512powers_basebits))); +#else + // Fallback for scalar. + uint32_t exp = ((BitCast(di, v).raw >> 23) - 118) & 0xf; + auto mul = BitCast(d, Set(di, (k2to512powers_25to18bits[exp] << 18) | + (k2to512powers_17to10bits[exp] << 10) | + k2to512powers_basebits)); +#endif + return IfThenElse(Lt(v, Set(d, 0.0031308f)), Mul(v, Set(d, 12.92f)), + MulAdd(pow, mul, Set(d, -0.055))); +} + +// NOLINTNEXTLINE(google-readability-namespace-comments) +} // namespace HWY_NAMESPACE +} // namespace jxl +HWY_AFTER_NAMESPACE(); + +#endif // LIB_JXL_CMS_TRANSFER_FUNCTIONS_INL_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/transfer_functions.h b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions.h new file mode 100644 index 0000000000..4e5273d5d3 --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions.h @@ -0,0 +1,131 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// Transfer functions for color encodings. + +#ifndef LIB_JXL_CMS_TRANSFER_FUNCTIONS_H_ +#define LIB_JXL_CMS_TRANSFER_FUNCTIONS_H_ + +#include <algorithm> +#include <cmath> + +#include "lib/jxl/base/status.h" + +namespace jxl { + +// Definitions for BT.2100-2 transfer functions (used inside/outside SIMD): +// "display" is linear light (nits) normalized to [0, 1]. +// "encoded" is a nonlinear encoding (e.g. PQ) in [0, 1]. +// "scene" is a linear function of photon counts, normalized to [0, 1]. + +// Despite the stated ranges, we need unbounded transfer functions: see +// http://www.littlecms.com/CIC18_UnboundedCMM.pdf. Inputs can be negative or +// above 1 due to chromatic adaptation. To avoid severe round-trip errors caused +// by clamping, we mirror negative inputs via copysign (f(-x) = -f(x), see +// https://developer.apple.com/documentation/coregraphics/cgcolorspace/1644735-extendedsrgb) +// and extend the function domains above 1. + +// Hybrid Log-Gamma. +class TF_HLG_Base { + public: + // EOTF. e = encoded. + static double DisplayFromEncoded(const double e) { return OOTF(InvOETF(e)); } + + // Inverse EOTF. d = display. + static double EncodedFromDisplay(const double d) { return OETF(InvOOTF(d)); } + + private: + // OETF (defines the HLG approach). s = scene, returns encoded. + static double OETF(double s) { + if (s == 0.0) return 0.0; + const double original_sign = s; + s = std::abs(s); + + if (s <= kDiv12) return copysignf(std::sqrt(3.0 * s), original_sign); + + const double e = kA * std::log(12 * s - kB) + kC; + JXL_ASSERT(e > 0.0); + return copysignf(e, original_sign); + } + + // e = encoded, returns scene. + static double InvOETF(double e) { + if (e == 0.0) return 0.0; + const double original_sign = e; + e = std::abs(e); + + if (e <= 0.5) return copysignf(e * e * (1.0 / 3), original_sign); + + const double s = (std::exp((e - kC) * kRA) + kB) * kDiv12; + JXL_ASSERT(s >= 0); + return copysignf(s, original_sign); + } + + // s = scene, returns display. + static double OOTF(const double s) { + // The actual (red channel) OOTF is RD = alpha * YS^(gamma-1) * RS, where + // YS = 0.2627 * RS + 0.6780 * GS + 0.0593 * BS. Let alpha = 1 so we return + // "display" (normalized [0, 1]) instead of nits. Our transfer function + // interface does not allow a dependency on YS. Fortunately, the system + // gamma at 334 nits is 1.0, so this reduces to RD = RS. + return s; + } + + // d = display, returns scene. + static double InvOOTF(const double d) { + return d; // see OOTF(). + } + + protected: + static constexpr double kA = 0.17883277; + static constexpr double kRA = 1.0 / kA; + static constexpr double kB = 1 - 4 * kA; + static constexpr double kC = 0.5599107295; + static constexpr double kDiv12 = 1.0 / 12; +}; + +// Perceptual Quantization +class TF_PQ_Base { + public: + static double DisplayFromEncoded(float display_intensity_target, double e) { + if (e == 0.0) return 0.0; + const double original_sign = e; + e = std::abs(e); + + const double xp = std::pow(e, 1.0 / kM2); + const double num = std::max(xp - kC1, 0.0); + const double den = kC2 - kC3 * xp; + JXL_DASSERT(den != 0.0); + const double d = std::pow(num / den, 1.0 / kM1); + JXL_DASSERT(d >= 0.0); // Equal for e ~= 1E-9 + return copysignf(d * (10000.0f / display_intensity_target), original_sign); + } + + // Inverse EOTF. d = display. + static double EncodedFromDisplay(float display_intensity_target, double d) { + if (d == 0.0) return 0.0; + const double original_sign = d; + d = std::abs(d); + + const double xp = + std::pow(d * (display_intensity_target * (1.0f / 10000.0f)), kM1); + const double num = kC1 + xp * kC2; + const double den = 1.0 + xp * kC3; + const double e = std::pow(num / den, kM2); + JXL_DASSERT(e > 0.0); + return copysignf(e, original_sign); + } + + protected: + static constexpr double kM1 = 2610.0 / 16384; + static constexpr double kM2 = (2523.0 / 4096) * 128; + static constexpr double kC1 = 3424.0 / 4096; + static constexpr double kC2 = (2413.0 / 4096) * 32; + static constexpr double kC3 = (2392.0 / 4096) * 32; +}; + +} // namespace jxl + +#endif // LIB_JXL_CMS_TRANSFER_FUNCTIONS_H_ diff --git a/third_party/jpeg-xl/lib/jxl/cms/transfer_functions_test.cc b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions_test.cc new file mode 100644 index 0000000000..26de409a4e --- /dev/null +++ b/third_party/jpeg-xl/lib/jxl/cms/transfer_functions_test.cc @@ -0,0 +1,94 @@ +// Copyright (c) the JPEG XL Project Authors. All rights reserved. +// +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "lib/jxl/cms/transfer_functions_test.cc" +#include "lib/jxl/cms/transfer_functions.h" + +#include <cstdio> +#include <hwy/foreach_target.h> + +#include "lib/jxl/base/random.h" +#include "lib/jxl/cms/transfer_functions-inl.h" +#include "lib/jxl/testing.h" + +// Test utils +#include <hwy/highway.h> +#include <hwy/tests/hwy_gtest.h> +HWY_BEFORE_NAMESPACE(); +namespace jxl { +namespace HWY_NAMESPACE { +namespace { + +HWY_NOINLINE void TestPqEncodedFromDisplay() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + double intensity = 11000.0 + rng.UniformF(-150.0f, 150.0f); + TF_PQ tf_pq(intensity); + const float f = rng.UniformF(0.0f, 1.0f); + const float actual = GetLane(tf_pq.EncodedFromDisplay(d, Set(d, f))); + const float expected = TF_PQ_Base::EncodedFromDisplay(intensity, f); + const float abs_err = std::abs(expected - actual); + EXPECT_LT(abs_err, 5e-7) << "f = " << f; + max_abs_err = std::max(max_abs_err, abs_err); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +HWY_NOINLINE void TestHlgEncodedFromDisplay() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + const float f = rng.UniformF(0.0f, 1.0f); + const float actual = GetLane(TF_HLG().EncodedFromDisplay(d, Set(d, f))); + const float expected = TF_HLG_Base::EncodedFromDisplay(f); + const float abs_err = std::abs(expected - actual); + EXPECT_LT(abs_err, 4e-7) << "f = " << f; + max_abs_err = std::max(max_abs_err, abs_err); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +HWY_NOINLINE void TestPqDisplayFromEncoded() { + constexpr size_t kNumTrials = 1 << 23; + Rng rng(1); + float max_abs_err = 0; + HWY_FULL(float) d; + for (size_t i = 0; i < kNumTrials; i++) { + double intensity = 11000.0 + rng.UniformF(-150.0f, 150.0f); + TF_PQ tf_pq(intensity); + const float f = rng.UniformF(0.0f, 1.0f); + const float actual = GetLane(tf_pq.DisplayFromEncoded(d, Set(d, f))); + const float expected = TF_PQ_Base::DisplayFromEncoded(intensity, f); + const float abs_err = std::abs(expected - actual); + EXPECT_LT(abs_err, 3E-6) << "f = " << f; + max_abs_err = std::max(max_abs_err, abs_err); + } + printf("max abs err %e\n", static_cast<double>(max_abs_err)); +} + +} // namespace +// NOLINTNEXTLINE(google-readability-namespace-comments) +} // namespace HWY_NAMESPACE +} // namespace jxl +HWY_AFTER_NAMESPACE(); + +#if HWY_ONCE +namespace jxl { + +class TransferFunctionsTargetTest : public hwy::TestWithParamTarget {}; +HWY_TARGET_INSTANTIATE_TEST_SUITE_P(TransferFunctionsTargetTest); + +HWY_EXPORT_AND_TEST_P(TransferFunctionsTargetTest, TestPqEncodedFromDisplay); +HWY_EXPORT_AND_TEST_P(TransferFunctionsTargetTest, TestHlgEncodedFromDisplay); +HWY_EXPORT_AND_TEST_P(TransferFunctionsTargetTest, TestPqDisplayFromEncoded); + +} // namespace jxl +#endif // HWY_ONCE |