summaryrefslogtreecommitdiffstats
path: root/gfx/skia/skia/modules/skcms
diff options
context:
space:
mode:
authorDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 19:33:14 +0000
committerDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 19:33:14 +0000
commit36d22d82aa202bb199967e9512281e9a53db42c9 (patch)
tree105e8c98ddea1c1e4784a60a5a6410fa416be2de /gfx/skia/skia/modules/skcms
parentInitial commit. (diff)
downloadfirefox-esr-36d22d82aa202bb199967e9512281e9a53db42c9.tar.xz
firefox-esr-36d22d82aa202bb199967e9512281e9a53db42c9.zip
Adding upstream version 115.7.0esr.upstream/115.7.0esr
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'gfx/skia/skia/modules/skcms')
-rw-r--r--gfx/skia/skia/modules/skcms/README.chromium5
-rw-r--r--gfx/skia/skia/modules/skcms/skcms.cc3064
-rw-r--r--gfx/skia/skia/modules/skcms/skcms.gni20
-rw-r--r--gfx/skia/skia/modules/skcms/skcms.h418
-rw-r--r--gfx/skia/skia/modules/skcms/skcms_internal.h56
-rw-r--r--gfx/skia/skia/modules/skcms/src/Transform_inl.h1628
-rwxr-xr-xgfx/skia/skia/modules/skcms/version.sha11
7 files changed, 5192 insertions, 0 deletions
diff --git a/gfx/skia/skia/modules/skcms/README.chromium b/gfx/skia/skia/modules/skcms/README.chromium
new file mode 100644
index 0000000000..046f6b1d19
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/README.chromium
@@ -0,0 +1,5 @@
+Name: skcms
+URL: https://skia.org/
+Version: unknown
+Security Critical: yes
+License: BSD
diff --git a/gfx/skia/skia/modules/skcms/skcms.cc b/gfx/skia/skia/modules/skcms/skcms.cc
new file mode 100644
index 0000000000..246c08af94
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/skcms.cc
@@ -0,0 +1,3064 @@
+/*
+ * Copyright 2018 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "skcms.h"
+#include "skcms_internal.h"
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include <stdlib.h>
+#include <string.h>
+
+#if defined(__ARM_NEON)
+ #include <arm_neon.h>
+#elif defined(__SSE__)
+ #include <immintrin.h>
+
+ #if defined(__clang__)
+ // That #include <immintrin.h> is usually enough, but Clang's headers
+ // "helpfully" skip including the whole kitchen sink when _MSC_VER is
+ // defined, because lots of programs on Windows would include that and
+ // it'd be a lot slower. But we want all those headers included so we
+ // can use their features after runtime checks later.
+ #include <smmintrin.h>
+ #include <avxintrin.h>
+ #include <avx2intrin.h>
+ #include <avx512fintrin.h>
+ #include <avx512dqintrin.h>
+ #endif
+#endif
+
+static bool runtime_cpu_detection = true;
+void skcms_DisableRuntimeCPUDetection() {
+ runtime_cpu_detection = false;
+}
+
+// sizeof(x) will return size_t, which is 32-bit on some machines and 64-bit on others.
+// We have better testing on 64-bit machines, so force 32-bit machines to behave like 64-bit.
+//
+// Please do not use sizeof() directly, and size_t only when required.
+// (We have no way of enforcing these requests...)
+#define SAFE_SIZEOF(x) ((uint64_t)sizeof(x))
+
+// Same sort of thing for _Layout structs with a variable sized array at the end (named "variable").
+#define SAFE_FIXED_SIZE(type) ((uint64_t)offsetof(type, variable))
+
+static const union {
+ uint32_t bits;
+ float f;
+} inf_ = { 0x7f800000 };
+#define INFINITY_ inf_.f
+
+#if defined(__clang__) || defined(__GNUC__)
+ #define small_memcpy __builtin_memcpy
+#else
+ #define small_memcpy memcpy
+#endif
+
+static float log2f_(float x) {
+ // The first approximation of log2(x) is its exponent 'e', minus 127.
+ int32_t bits;
+ small_memcpy(&bits, &x, sizeof(bits));
+
+ float e = (float)bits * (1.0f / (1<<23));
+
+ // If we use the mantissa too we can refine the error signficantly.
+ int32_t m_bits = (bits & 0x007fffff) | 0x3f000000;
+ float m;
+ small_memcpy(&m, &m_bits, sizeof(m));
+
+ return (e - 124.225514990f
+ - 1.498030302f*m
+ - 1.725879990f/(0.3520887068f + m));
+}
+static float logf_(float x) {
+ const float ln2 = 0.69314718f;
+ return ln2*log2f_(x);
+}
+
+static float exp2f_(float x) {
+ float fract = x - floorf_(x);
+
+ float fbits = (1.0f * (1<<23)) * (x + 121.274057500f
+ - 1.490129070f*fract
+ + 27.728023300f/(4.84252568f - fract));
+
+ // Before we cast fbits to int32_t, check for out of range values to pacify UBSAN.
+ // INT_MAX is not exactly representable as a float, so exclude it as effectively infinite.
+ // Negative values are effectively underflow - we'll end up returning a (different) negative
+ // value, which makes no sense. So clamp to zero.
+ if (fbits >= (float)INT_MAX) {
+ return INFINITY_;
+ } else if (fbits < 0) {
+ return 0;
+ }
+
+ int32_t bits = (int32_t)fbits;
+ small_memcpy(&x, &bits, sizeof(x));
+ return x;
+}
+
+// Not static, as it's used by some test tools.
+float powf_(float x, float y) {
+ assert (x >= 0);
+ return (x == 0) || (x == 1) ? x
+ : exp2f_(log2f_(x) * y);
+}
+
+static float expf_(float x) {
+ const float log2_e = 1.4426950408889634074f;
+ return exp2f_(log2_e * x);
+}
+
+static float fmaxf_(float x, float y) { return x > y ? x : y; }
+static float fminf_(float x, float y) { return x < y ? x : y; }
+
+static bool isfinitef_(float x) { return 0 == x*0; }
+
+static float minus_1_ulp(float x) {
+ int32_t bits;
+ memcpy(&bits, &x, sizeof(bits));
+ bits = bits - 1;
+ memcpy(&x, &bits, sizeof(bits));
+ return x;
+}
+
+// Most transfer functions we work with are sRGBish.
+// For exotic HDR transfer functions, we encode them using a tf.g that makes no sense,
+// and repurpose the other fields to hold the parameters of the HDR functions.
+struct TF_PQish { float A,B,C,D,E,F; };
+struct TF_HLGish { float R,G,a,b,c,K_minus_1; };
+// We didn't originally support a scale factor K for HLG, and instead just stored 0 in
+// the unused `f` field of skcms_TransferFunction for HLGish and HLGInvish transfer functions.
+// By storing f=K-1, those old unusued f=0 values now mean K=1, a noop scale factor.
+
+static float TFKind_marker(skcms_TFType kind) {
+ // We'd use different NaNs, but those aren't guaranteed to be preserved by WASM.
+ return -(float)kind;
+}
+
+static skcms_TFType classify(const skcms_TransferFunction& tf, TF_PQish* pq = nullptr
+ , TF_HLGish* hlg = nullptr) {
+ if (tf.g < 0 && static_cast<float>(static_cast<int>(tf.g)) == tf.g) {
+ // TODO: soundness checks for PQ/HLG like we do for sRGBish?
+ switch ((int)tf.g) {
+ case -skcms_TFType_PQish:
+ if (pq) {
+ memcpy(pq , &tf.a, sizeof(*pq ));
+ }
+ return skcms_TFType_PQish;
+ case -skcms_TFType_HLGish:
+ if (hlg) {
+ memcpy(hlg, &tf.a, sizeof(*hlg));
+ }
+ return skcms_TFType_HLGish;
+ case -skcms_TFType_HLGinvish:
+ if (hlg) {
+ memcpy(hlg, &tf.a, sizeof(*hlg));
+ }
+ return skcms_TFType_HLGinvish;
+ }
+ return skcms_TFType_Invalid;
+ }
+
+ // Basic soundness checks for sRGBish transfer functions.
+ if (isfinitef_(tf.a + tf.b + tf.c + tf.d + tf.e + tf.f + tf.g)
+ // a,c,d,g should be non-negative to make any sense.
+ && tf.a >= 0
+ && tf.c >= 0
+ && tf.d >= 0
+ && tf.g >= 0
+ // Raising a negative value to a fractional tf->g produces complex numbers.
+ && tf.a * tf.d + tf.b >= 0) {
+ return skcms_TFType_sRGBish;
+ }
+
+ return skcms_TFType_Invalid;
+}
+
+skcms_TFType skcms_TransferFunction_getType(const skcms_TransferFunction* tf) {
+ return classify(*tf);
+}
+bool skcms_TransferFunction_isSRGBish(const skcms_TransferFunction* tf) {
+ return classify(*tf) == skcms_TFType_sRGBish;
+}
+bool skcms_TransferFunction_isPQish(const skcms_TransferFunction* tf) {
+ return classify(*tf) == skcms_TFType_PQish;
+}
+bool skcms_TransferFunction_isHLGish(const skcms_TransferFunction* tf) {
+ return classify(*tf) == skcms_TFType_HLGish;
+}
+
+bool skcms_TransferFunction_makePQish(skcms_TransferFunction* tf,
+ float A, float B, float C,
+ float D, float E, float F) {
+ *tf = { TFKind_marker(skcms_TFType_PQish), A,B,C,D,E,F };
+ assert(skcms_TransferFunction_isPQish(tf));
+ return true;
+}
+
+bool skcms_TransferFunction_makeScaledHLGish(skcms_TransferFunction* tf,
+ float K, float R, float G,
+ float a, float b, float c) {
+ *tf = { TFKind_marker(skcms_TFType_HLGish), R,G, a,b,c, K-1.0f };
+ assert(skcms_TransferFunction_isHLGish(tf));
+ return true;
+}
+
+float skcms_TransferFunction_eval(const skcms_TransferFunction* tf, float x) {
+ float sign = x < 0 ? -1.0f : 1.0f;
+ x *= sign;
+
+ TF_PQish pq;
+ TF_HLGish hlg;
+ switch (classify(*tf, &pq, &hlg)) {
+ case skcms_TFType_Invalid: break;
+
+ case skcms_TFType_HLGish: {
+ const float K = hlg.K_minus_1 + 1.0f;
+ return K * sign * (x*hlg.R <= 1 ? powf_(x*hlg.R, hlg.G)
+ : expf_((x-hlg.c)*hlg.a) + hlg.b);
+ }
+
+ // skcms_TransferFunction_invert() inverts R, G, and a for HLGinvish so this math is fast.
+ case skcms_TFType_HLGinvish: {
+ const float K = hlg.K_minus_1 + 1.0f;
+ x /= K;
+ return sign * (x <= 1 ? hlg.R * powf_(x, hlg.G)
+ : hlg.a * logf_(x - hlg.b) + hlg.c);
+ }
+
+ case skcms_TFType_sRGBish:
+ return sign * (x < tf->d ? tf->c * x + tf->f
+ : powf_(tf->a * x + tf->b, tf->g) + tf->e);
+
+ case skcms_TFType_PQish: return sign * powf_(fmaxf_(pq.A + pq.B * powf_(x, pq.C), 0)
+ / (pq.D + pq.E * powf_(x, pq.C)),
+ pq.F);
+ }
+ return 0;
+}
+
+
+static float eval_curve(const skcms_Curve* curve, float x) {
+ if (curve->table_entries == 0) {
+ return skcms_TransferFunction_eval(&curve->parametric, x);
+ }
+
+ float ix = fmaxf_(0, fminf_(x, 1)) * static_cast<float>(curve->table_entries - 1);
+ int lo = (int) ix ,
+ hi = (int)(float)minus_1_ulp(ix + 1.0f);
+ float t = ix - (float)lo;
+
+ float l, h;
+ if (curve->table_8) {
+ l = curve->table_8[lo] * (1/255.0f);
+ h = curve->table_8[hi] * (1/255.0f);
+ } else {
+ uint16_t be_l, be_h;
+ memcpy(&be_l, curve->table_16 + 2*lo, 2);
+ memcpy(&be_h, curve->table_16 + 2*hi, 2);
+ uint16_t le_l = ((be_l << 8) | (be_l >> 8)) & 0xffff;
+ uint16_t le_h = ((be_h << 8) | (be_h >> 8)) & 0xffff;
+ l = le_l * (1/65535.0f);
+ h = le_h * (1/65535.0f);
+ }
+ return l + (h-l)*t;
+}
+
+float skcms_MaxRoundtripError(const skcms_Curve* curve, const skcms_TransferFunction* inv_tf) {
+ uint32_t N = curve->table_entries > 256 ? curve->table_entries : 256;
+ const float dx = 1.0f / static_cast<float>(N - 1);
+ float err = 0;
+ for (uint32_t i = 0; i < N; i++) {
+ float x = static_cast<float>(i) * dx,
+ y = eval_curve(curve, x);
+ err = fmaxf_(err, fabsf_(x - skcms_TransferFunction_eval(inv_tf, y)));
+ }
+ return err;
+}
+
+bool skcms_AreApproximateInverses(const skcms_Curve* curve, const skcms_TransferFunction* inv_tf) {
+ return skcms_MaxRoundtripError(curve, inv_tf) < (1/512.0f);
+}
+
+// Additional ICC signature values that are only used internally
+enum {
+ // File signature
+ skcms_Signature_acsp = 0x61637370,
+
+ // Tag signatures
+ skcms_Signature_rTRC = 0x72545243,
+ skcms_Signature_gTRC = 0x67545243,
+ skcms_Signature_bTRC = 0x62545243,
+ skcms_Signature_kTRC = 0x6B545243,
+
+ skcms_Signature_rXYZ = 0x7258595A,
+ skcms_Signature_gXYZ = 0x6758595A,
+ skcms_Signature_bXYZ = 0x6258595A,
+
+ skcms_Signature_A2B0 = 0x41324230,
+ skcms_Signature_B2A0 = 0x42324130,
+
+ skcms_Signature_CHAD = 0x63686164,
+ skcms_Signature_WTPT = 0x77747074,
+
+ skcms_Signature_CICP = 0x63696370,
+
+ // Type signatures
+ skcms_Signature_curv = 0x63757276,
+ skcms_Signature_mft1 = 0x6D667431,
+ skcms_Signature_mft2 = 0x6D667432,
+ skcms_Signature_mAB = 0x6D414220,
+ skcms_Signature_mBA = 0x6D424120,
+ skcms_Signature_para = 0x70617261,
+ skcms_Signature_sf32 = 0x73663332,
+ // XYZ is also a PCS signature, so it's defined in skcms.h
+ // skcms_Signature_XYZ = 0x58595A20,
+};
+
+static uint16_t read_big_u16(const uint8_t* ptr) {
+ uint16_t be;
+ memcpy(&be, ptr, sizeof(be));
+#if defined(_MSC_VER)
+ return _byteswap_ushort(be);
+#else
+ return __builtin_bswap16(be);
+#endif
+}
+
+static uint32_t read_big_u32(const uint8_t* ptr) {
+ uint32_t be;
+ memcpy(&be, ptr, sizeof(be));
+#if defined(_MSC_VER)
+ return _byteswap_ulong(be);
+#else
+ return __builtin_bswap32(be);
+#endif
+}
+
+static int32_t read_big_i32(const uint8_t* ptr) {
+ return (int32_t)read_big_u32(ptr);
+}
+
+static float read_big_fixed(const uint8_t* ptr) {
+ return static_cast<float>(read_big_i32(ptr)) * (1.0f / 65536.0f);
+}
+
+// Maps to an in-memory profile so that fields line up to the locations specified
+// in ICC.1:2010, section 7.2
+typedef struct {
+ uint8_t size [ 4];
+ uint8_t cmm_type [ 4];
+ uint8_t version [ 4];
+ uint8_t profile_class [ 4];
+ uint8_t data_color_space [ 4];
+ uint8_t pcs [ 4];
+ uint8_t creation_date_time [12];
+ uint8_t signature [ 4];
+ uint8_t platform [ 4];
+ uint8_t flags [ 4];
+ uint8_t device_manufacturer [ 4];
+ uint8_t device_model [ 4];
+ uint8_t device_attributes [ 8];
+ uint8_t rendering_intent [ 4];
+ uint8_t illuminant_X [ 4];
+ uint8_t illuminant_Y [ 4];
+ uint8_t illuminant_Z [ 4];
+ uint8_t creator [ 4];
+ uint8_t profile_id [16];
+ uint8_t reserved [28];
+ uint8_t tag_count [ 4]; // Technically not part of header, but required
+} header_Layout;
+
+typedef struct {
+ uint8_t signature [4];
+ uint8_t offset [4];
+ uint8_t size [4];
+} tag_Layout;
+
+static const tag_Layout* get_tag_table(const skcms_ICCProfile* profile) {
+ return (const tag_Layout*)(profile->buffer + SAFE_SIZEOF(header_Layout));
+}
+
+// s15Fixed16ArrayType is technically variable sized, holding N values. However, the only valid
+// use of the type is for the CHAD tag that stores exactly nine values.
+typedef struct {
+ uint8_t type [ 4];
+ uint8_t reserved [ 4];
+ uint8_t values [36];
+} sf32_Layout;
+
+bool skcms_GetCHAD(const skcms_ICCProfile* profile, skcms_Matrix3x3* m) {
+ skcms_ICCTag tag;
+ if (!skcms_GetTagBySignature(profile, skcms_Signature_CHAD, &tag)) {
+ return false;
+ }
+
+ if (tag.type != skcms_Signature_sf32 || tag.size < SAFE_SIZEOF(sf32_Layout)) {
+ return false;
+ }
+
+ const sf32_Layout* sf32Tag = (const sf32_Layout*)tag.buf;
+ const uint8_t* values = sf32Tag->values;
+ for (int r = 0; r < 3; ++r)
+ for (int c = 0; c < 3; ++c, values += 4) {
+ m->vals[r][c] = read_big_fixed(values);
+ }
+ return true;
+}
+
+// XYZType is technically variable sized, holding N XYZ triples. However, the only valid uses of
+// the type are for tags/data that store exactly one triple.
+typedef struct {
+ uint8_t type [4];
+ uint8_t reserved [4];
+ uint8_t X [4];
+ uint8_t Y [4];
+ uint8_t Z [4];
+} XYZ_Layout;
+
+static bool read_tag_xyz(const skcms_ICCTag* tag, float* x, float* y, float* z) {
+ if (tag->type != skcms_Signature_XYZ || tag->size < SAFE_SIZEOF(XYZ_Layout)) {
+ return false;
+ }
+
+ const XYZ_Layout* xyzTag = (const XYZ_Layout*)tag->buf;
+
+ *x = read_big_fixed(xyzTag->X);
+ *y = read_big_fixed(xyzTag->Y);
+ *z = read_big_fixed(xyzTag->Z);
+ return true;
+}
+
+bool skcms_GetWTPT(const skcms_ICCProfile* profile, float xyz[3]) {
+ skcms_ICCTag tag;
+ return skcms_GetTagBySignature(profile, skcms_Signature_WTPT, &tag) &&
+ read_tag_xyz(&tag, &xyz[0], &xyz[1], &xyz[2]);
+}
+
+static bool read_to_XYZD50(const skcms_ICCTag* rXYZ, const skcms_ICCTag* gXYZ,
+ const skcms_ICCTag* bXYZ, skcms_Matrix3x3* toXYZ) {
+ return read_tag_xyz(rXYZ, &toXYZ->vals[0][0], &toXYZ->vals[1][0], &toXYZ->vals[2][0]) &&
+ read_tag_xyz(gXYZ, &toXYZ->vals[0][1], &toXYZ->vals[1][1], &toXYZ->vals[2][1]) &&
+ read_tag_xyz(bXYZ, &toXYZ->vals[0][2], &toXYZ->vals[1][2], &toXYZ->vals[2][2]);
+}
+
+typedef struct {
+ uint8_t type [4];
+ uint8_t reserved_a [4];
+ uint8_t function_type [2];
+ uint8_t reserved_b [2];
+ uint8_t variable [1/*variable*/]; // 1, 3, 4, 5, or 7 s15.16, depending on function_type
+} para_Layout;
+
+static bool read_curve_para(const uint8_t* buf, uint32_t size,
+ skcms_Curve* curve, uint32_t* curve_size) {
+ if (size < SAFE_FIXED_SIZE(para_Layout)) {
+ return false;
+ }
+
+ const para_Layout* paraTag = (const para_Layout*)buf;
+
+ enum { kG = 0, kGAB = 1, kGABC = 2, kGABCD = 3, kGABCDEF = 4 };
+ uint16_t function_type = read_big_u16(paraTag->function_type);
+ if (function_type > kGABCDEF) {
+ return false;
+ }
+
+ static const uint32_t curve_bytes[] = { 4, 12, 16, 20, 28 };
+ if (size < SAFE_FIXED_SIZE(para_Layout) + curve_bytes[function_type]) {
+ return false;
+ }
+
+ if (curve_size) {
+ *curve_size = SAFE_FIXED_SIZE(para_Layout) + curve_bytes[function_type];
+ }
+
+ curve->table_entries = 0;
+ curve->parametric.a = 1.0f;
+ curve->parametric.b = 0.0f;
+ curve->parametric.c = 0.0f;
+ curve->parametric.d = 0.0f;
+ curve->parametric.e = 0.0f;
+ curve->parametric.f = 0.0f;
+ curve->parametric.g = read_big_fixed(paraTag->variable);
+
+ switch (function_type) {
+ case kGAB:
+ curve->parametric.a = read_big_fixed(paraTag->variable + 4);
+ curve->parametric.b = read_big_fixed(paraTag->variable + 8);
+ if (curve->parametric.a == 0) {
+ return false;
+ }
+ curve->parametric.d = -curve->parametric.b / curve->parametric.a;
+ break;
+ case kGABC:
+ curve->parametric.a = read_big_fixed(paraTag->variable + 4);
+ curve->parametric.b = read_big_fixed(paraTag->variable + 8);
+ curve->parametric.e = read_big_fixed(paraTag->variable + 12);
+ if (curve->parametric.a == 0) {
+ return false;
+ }
+ curve->parametric.d = -curve->parametric.b / curve->parametric.a;
+ curve->parametric.f = curve->parametric.e;
+ break;
+ case kGABCD:
+ curve->parametric.a = read_big_fixed(paraTag->variable + 4);
+ curve->parametric.b = read_big_fixed(paraTag->variable + 8);
+ curve->parametric.c = read_big_fixed(paraTag->variable + 12);
+ curve->parametric.d = read_big_fixed(paraTag->variable + 16);
+ break;
+ case kGABCDEF:
+ curve->parametric.a = read_big_fixed(paraTag->variable + 4);
+ curve->parametric.b = read_big_fixed(paraTag->variable + 8);
+ curve->parametric.c = read_big_fixed(paraTag->variable + 12);
+ curve->parametric.d = read_big_fixed(paraTag->variable + 16);
+ curve->parametric.e = read_big_fixed(paraTag->variable + 20);
+ curve->parametric.f = read_big_fixed(paraTag->variable + 24);
+ break;
+ }
+ return skcms_TransferFunction_isSRGBish(&curve->parametric);
+}
+
+typedef struct {
+ uint8_t type [4];
+ uint8_t reserved [4];
+ uint8_t value_count [4];
+ uint8_t variable [1/*variable*/]; // value_count, 8.8 if 1, uint16 (n*65535) if > 1
+} curv_Layout;
+
+static bool read_curve_curv(const uint8_t* buf, uint32_t size,
+ skcms_Curve* curve, uint32_t* curve_size) {
+ if (size < SAFE_FIXED_SIZE(curv_Layout)) {
+ return false;
+ }
+
+ const curv_Layout* curvTag = (const curv_Layout*)buf;
+
+ uint32_t value_count = read_big_u32(curvTag->value_count);
+ if (size < SAFE_FIXED_SIZE(curv_Layout) + value_count * SAFE_SIZEOF(uint16_t)) {
+ return false;
+ }
+
+ if (curve_size) {
+ *curve_size = SAFE_FIXED_SIZE(curv_Layout) + value_count * SAFE_SIZEOF(uint16_t);
+ }
+
+ if (value_count < 2) {
+ curve->table_entries = 0;
+ curve->parametric.a = 1.0f;
+ curve->parametric.b = 0.0f;
+ curve->parametric.c = 0.0f;
+ curve->parametric.d = 0.0f;
+ curve->parametric.e = 0.0f;
+ curve->parametric.f = 0.0f;
+ if (value_count == 0) {
+ // Empty tables are a shorthand for an identity curve
+ curve->parametric.g = 1.0f;
+ } else {
+ // Single entry tables are a shorthand for simple gamma
+ curve->parametric.g = read_big_u16(curvTag->variable) * (1.0f / 256.0f);
+ }
+ } else {
+ curve->table_8 = nullptr;
+ curve->table_16 = curvTag->variable;
+ curve->table_entries = value_count;
+ }
+
+ return true;
+}
+
+// Parses both curveType and parametricCurveType data. Ensures that at most 'size' bytes are read.
+// If curve_size is not nullptr, writes the number of bytes used by the curve in (*curve_size).
+static bool read_curve(const uint8_t* buf, uint32_t size,
+ skcms_Curve* curve, uint32_t* curve_size) {
+ if (!buf || size < 4 || !curve) {
+ return false;
+ }
+
+ uint32_t type = read_big_u32(buf);
+ if (type == skcms_Signature_para) {
+ return read_curve_para(buf, size, curve, curve_size);
+ } else if (type == skcms_Signature_curv) {
+ return read_curve_curv(buf, size, curve, curve_size);
+ }
+
+ return false;
+}
+
+// mft1 and mft2 share a large chunk of data
+typedef struct {
+ uint8_t type [ 4];
+ uint8_t reserved_a [ 4];
+ uint8_t input_channels [ 1];
+ uint8_t output_channels [ 1];
+ uint8_t grid_points [ 1];
+ uint8_t reserved_b [ 1];
+ uint8_t matrix [36];
+} mft_CommonLayout;
+
+typedef struct {
+ mft_CommonLayout common [1];
+
+ uint8_t variable [1/*variable*/];
+} mft1_Layout;
+
+typedef struct {
+ mft_CommonLayout common [1];
+
+ uint8_t input_table_entries [2];
+ uint8_t output_table_entries [2];
+ uint8_t variable [1/*variable*/];
+} mft2_Layout;
+
+static bool read_mft_common(const mft_CommonLayout* mftTag, skcms_A2B* a2b) {
+ // MFT matrices are applied before the first set of curves, but must be identity unless the
+ // input is PCSXYZ. We don't support PCSXYZ profiles, so we ignore this matrix. Note that the
+ // matrix in skcms_A2B is applied later in the pipe, so supporting this would require another
+ // field/flag.
+ a2b->matrix_channels = 0;
+ a2b-> input_channels = mftTag-> input_channels[0];
+ a2b->output_channels = mftTag->output_channels[0];
+
+ // We require exactly three (ie XYZ/Lab/RGB) output channels
+ if (a2b->output_channels != ARRAY_COUNT(a2b->output_curves)) {
+ return false;
+ }
+ // We require at least one, and no more than four (ie CMYK) input channels
+ if (a2b->input_channels < 1 || a2b->input_channels > ARRAY_COUNT(a2b->input_curves)) {
+ return false;
+ }
+
+ for (uint32_t i = 0; i < a2b->input_channels; ++i) {
+ a2b->grid_points[i] = mftTag->grid_points[0];
+ }
+ // The grid only makes sense with at least two points along each axis
+ if (a2b->grid_points[0] < 2) {
+ return false;
+ }
+ return true;
+}
+
+// All as the A2B version above, except where noted.
+static bool read_mft_common(const mft_CommonLayout* mftTag, skcms_B2A* b2a) {
+ // Same as A2B.
+ b2a->matrix_channels = 0;
+ b2a-> input_channels = mftTag-> input_channels[0];
+ b2a->output_channels = mftTag->output_channels[0];
+
+
+ // For B2A, exactly 3 input channels (XYZ) and 3 (RGB) or 4 (CMYK) output channels.
+ if (b2a->input_channels != ARRAY_COUNT(b2a->input_curves)) {
+ return false;
+ }
+ if (b2a->output_channels < 3 || b2a->output_channels > ARRAY_COUNT(b2a->output_curves)) {
+ return false;
+ }
+
+ // Same as A2B.
+ for (uint32_t i = 0; i < b2a->input_channels; ++i) {
+ b2a->grid_points[i] = mftTag->grid_points[0];
+ }
+ if (b2a->grid_points[0] < 2) {
+ return false;
+ }
+ return true;
+}
+
+template <typename A2B_or_B2A>
+static bool init_tables(const uint8_t* table_base, uint64_t max_tables_len, uint32_t byte_width,
+ uint32_t input_table_entries, uint32_t output_table_entries,
+ A2B_or_B2A* out) {
+ // byte_width is 1 or 2, [input|output]_table_entries are in [2, 4096], so no overflow
+ uint32_t byte_len_per_input_table = input_table_entries * byte_width;
+ uint32_t byte_len_per_output_table = output_table_entries * byte_width;
+
+ // [input|output]_channels are <= 4, so still no overflow
+ uint32_t byte_len_all_input_tables = out->input_channels * byte_len_per_input_table;
+ uint32_t byte_len_all_output_tables = out->output_channels * byte_len_per_output_table;
+
+ uint64_t grid_size = out->output_channels * byte_width;
+ for (uint32_t axis = 0; axis < out->input_channels; ++axis) {
+ grid_size *= out->grid_points[axis];
+ }
+
+ if (max_tables_len < byte_len_all_input_tables + grid_size + byte_len_all_output_tables) {
+ return false;
+ }
+
+ for (uint32_t i = 0; i < out->input_channels; ++i) {
+ out->input_curves[i].table_entries = input_table_entries;
+ if (byte_width == 1) {
+ out->input_curves[i].table_8 = table_base + i * byte_len_per_input_table;
+ out->input_curves[i].table_16 = nullptr;
+ } else {
+ out->input_curves[i].table_8 = nullptr;
+ out->input_curves[i].table_16 = table_base + i * byte_len_per_input_table;
+ }
+ }
+
+ if (byte_width == 1) {
+ out->grid_8 = table_base + byte_len_all_input_tables;
+ out->grid_16 = nullptr;
+ } else {
+ out->grid_8 = nullptr;
+ out->grid_16 = table_base + byte_len_all_input_tables;
+ }
+
+ const uint8_t* output_table_base = table_base + byte_len_all_input_tables + grid_size;
+ for (uint32_t i = 0; i < out->output_channels; ++i) {
+ out->output_curves[i].table_entries = output_table_entries;
+ if (byte_width == 1) {
+ out->output_curves[i].table_8 = output_table_base + i * byte_len_per_output_table;
+ out->output_curves[i].table_16 = nullptr;
+ } else {
+ out->output_curves[i].table_8 = nullptr;
+ out->output_curves[i].table_16 = output_table_base + i * byte_len_per_output_table;
+ }
+ }
+
+ return true;
+}
+
+template <typename A2B_or_B2A>
+static bool read_tag_mft1(const skcms_ICCTag* tag, A2B_or_B2A* out) {
+ if (tag->size < SAFE_FIXED_SIZE(mft1_Layout)) {
+ return false;
+ }
+
+ const mft1_Layout* mftTag = (const mft1_Layout*)tag->buf;
+ if (!read_mft_common(mftTag->common, out)) {
+ return false;
+ }
+
+ uint32_t input_table_entries = 256;
+ uint32_t output_table_entries = 256;
+
+ return init_tables(mftTag->variable, tag->size - SAFE_FIXED_SIZE(mft1_Layout), 1,
+ input_table_entries, output_table_entries, out);
+}
+
+template <typename A2B_or_B2A>
+static bool read_tag_mft2(const skcms_ICCTag* tag, A2B_or_B2A* out) {
+ if (tag->size < SAFE_FIXED_SIZE(mft2_Layout)) {
+ return false;
+ }
+
+ const mft2_Layout* mftTag = (const mft2_Layout*)tag->buf;
+ if (!read_mft_common(mftTag->common, out)) {
+ return false;
+ }
+
+ uint32_t input_table_entries = read_big_u16(mftTag->input_table_entries);
+ uint32_t output_table_entries = read_big_u16(mftTag->output_table_entries);
+
+ // ICC spec mandates that 2 <= table_entries <= 4096
+ if (input_table_entries < 2 || input_table_entries > 4096 ||
+ output_table_entries < 2 || output_table_entries > 4096) {
+ return false;
+ }
+
+ return init_tables(mftTag->variable, tag->size - SAFE_FIXED_SIZE(mft2_Layout), 2,
+ input_table_entries, output_table_entries, out);
+}
+
+static bool read_curves(const uint8_t* buf, uint32_t size, uint32_t curve_offset,
+ uint32_t num_curves, skcms_Curve* curves) {
+ for (uint32_t i = 0; i < num_curves; ++i) {
+ if (curve_offset > size) {
+ return false;
+ }
+
+ uint32_t curve_bytes;
+ if (!read_curve(buf + curve_offset, size - curve_offset, &curves[i], &curve_bytes)) {
+ return false;
+ }
+
+ if (curve_bytes > UINT32_MAX - 3) {
+ return false;
+ }
+ curve_bytes = (curve_bytes + 3) & ~3U;
+
+ uint64_t new_offset_64 = (uint64_t)curve_offset + curve_bytes;
+ curve_offset = (uint32_t)new_offset_64;
+ if (new_offset_64 != curve_offset) {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+// mAB and mBA tags use the same encoding, including color lookup tables.
+typedef struct {
+ uint8_t type [ 4];
+ uint8_t reserved_a [ 4];
+ uint8_t input_channels [ 1];
+ uint8_t output_channels [ 1];
+ uint8_t reserved_b [ 2];
+ uint8_t b_curve_offset [ 4];
+ uint8_t matrix_offset [ 4];
+ uint8_t m_curve_offset [ 4];
+ uint8_t clut_offset [ 4];
+ uint8_t a_curve_offset [ 4];
+} mAB_or_mBA_Layout;
+
+typedef struct {
+ uint8_t grid_points [16];
+ uint8_t grid_byte_width [ 1];
+ uint8_t reserved [ 3];
+ uint8_t variable [1/*variable*/];
+} CLUT_Layout;
+
+static bool read_tag_mab(const skcms_ICCTag* tag, skcms_A2B* a2b, bool pcs_is_xyz) {
+ if (tag->size < SAFE_SIZEOF(mAB_or_mBA_Layout)) {
+ return false;
+ }
+
+ const mAB_or_mBA_Layout* mABTag = (const mAB_or_mBA_Layout*)tag->buf;
+
+ a2b->input_channels = mABTag->input_channels[0];
+ a2b->output_channels = mABTag->output_channels[0];
+
+ // We require exactly three (ie XYZ/Lab/RGB) output channels
+ if (a2b->output_channels != ARRAY_COUNT(a2b->output_curves)) {
+ return false;
+ }
+ // We require no more than four (ie CMYK) input channels
+ if (a2b->input_channels > ARRAY_COUNT(a2b->input_curves)) {
+ return false;
+ }
+
+ uint32_t b_curve_offset = read_big_u32(mABTag->b_curve_offset);
+ uint32_t matrix_offset = read_big_u32(mABTag->matrix_offset);
+ uint32_t m_curve_offset = read_big_u32(mABTag->m_curve_offset);
+ uint32_t clut_offset = read_big_u32(mABTag->clut_offset);
+ uint32_t a_curve_offset = read_big_u32(mABTag->a_curve_offset);
+
+ // "B" curves must be present
+ if (0 == b_curve_offset) {
+ return false;
+ }
+
+ if (!read_curves(tag->buf, tag->size, b_curve_offset, a2b->output_channels,
+ a2b->output_curves)) {
+ return false;
+ }
+
+ // "M" curves and Matrix must be used together
+ if (0 != m_curve_offset) {
+ if (0 == matrix_offset) {
+ return false;
+ }
+ a2b->matrix_channels = a2b->output_channels;
+ if (!read_curves(tag->buf, tag->size, m_curve_offset, a2b->matrix_channels,
+ a2b->matrix_curves)) {
+ return false;
+ }
+
+ // Read matrix, which is stored as a row-major 3x3, followed by the fourth column
+ if (tag->size < matrix_offset + 12 * SAFE_SIZEOF(uint32_t)) {
+ return false;
+ }
+ float encoding_factor = pcs_is_xyz ? (65535 / 32768.0f) : 1.0f;
+ const uint8_t* mtx_buf = tag->buf + matrix_offset;
+ a2b->matrix.vals[0][0] = encoding_factor * read_big_fixed(mtx_buf + 0);
+ a2b->matrix.vals[0][1] = encoding_factor * read_big_fixed(mtx_buf + 4);
+ a2b->matrix.vals[0][2] = encoding_factor * read_big_fixed(mtx_buf + 8);
+ a2b->matrix.vals[1][0] = encoding_factor * read_big_fixed(mtx_buf + 12);
+ a2b->matrix.vals[1][1] = encoding_factor * read_big_fixed(mtx_buf + 16);
+ a2b->matrix.vals[1][2] = encoding_factor * read_big_fixed(mtx_buf + 20);
+ a2b->matrix.vals[2][0] = encoding_factor * read_big_fixed(mtx_buf + 24);
+ a2b->matrix.vals[2][1] = encoding_factor * read_big_fixed(mtx_buf + 28);
+ a2b->matrix.vals[2][2] = encoding_factor * read_big_fixed(mtx_buf + 32);
+ a2b->matrix.vals[0][3] = encoding_factor * read_big_fixed(mtx_buf + 36);
+ a2b->matrix.vals[1][3] = encoding_factor * read_big_fixed(mtx_buf + 40);
+ a2b->matrix.vals[2][3] = encoding_factor * read_big_fixed(mtx_buf + 44);
+ } else {
+ if (0 != matrix_offset) {
+ return false;
+ }
+ a2b->matrix_channels = 0;
+ }
+
+ // "A" curves and CLUT must be used together
+ if (0 != a_curve_offset) {
+ if (0 == clut_offset) {
+ return false;
+ }
+ if (!read_curves(tag->buf, tag->size, a_curve_offset, a2b->input_channels,
+ a2b->input_curves)) {
+ return false;
+ }
+
+ if (tag->size < clut_offset + SAFE_FIXED_SIZE(CLUT_Layout)) {
+ return false;
+ }
+ const CLUT_Layout* clut = (const CLUT_Layout*)(tag->buf + clut_offset);
+
+ if (clut->grid_byte_width[0] == 1) {
+ a2b->grid_8 = clut->variable;
+ a2b->grid_16 = nullptr;
+ } else if (clut->grid_byte_width[0] == 2) {
+ a2b->grid_8 = nullptr;
+ a2b->grid_16 = clut->variable;
+ } else {
+ return false;
+ }
+
+ uint64_t grid_size = a2b->output_channels * clut->grid_byte_width[0]; // the payload
+ for (uint32_t i = 0; i < a2b->input_channels; ++i) {
+ a2b->grid_points[i] = clut->grid_points[i];
+ // The grid only makes sense with at least two points along each axis
+ if (a2b->grid_points[i] < 2) {
+ return false;
+ }
+ grid_size *= a2b->grid_points[i];
+ }
+ if (tag->size < clut_offset + SAFE_FIXED_SIZE(CLUT_Layout) + grid_size) {
+ return false;
+ }
+ } else {
+ if (0 != clut_offset) {
+ return false;
+ }
+
+ // If there is no CLUT, the number of input and output channels must match
+ if (a2b->input_channels != a2b->output_channels) {
+ return false;
+ }
+
+ // Zero out the number of input channels to signal that we're skipping this stage
+ a2b->input_channels = 0;
+ }
+
+ return true;
+}
+
+// Exactly the same as read_tag_mab(), except where there are comments.
+// TODO: refactor the two to eliminate common code?
+static bool read_tag_mba(const skcms_ICCTag* tag, skcms_B2A* b2a, bool pcs_is_xyz) {
+ if (tag->size < SAFE_SIZEOF(mAB_or_mBA_Layout)) {
+ return false;
+ }
+
+ const mAB_or_mBA_Layout* mBATag = (const mAB_or_mBA_Layout*)tag->buf;
+
+ b2a->input_channels = mBATag->input_channels[0];
+ b2a->output_channels = mBATag->output_channels[0];
+
+ // Require exactly 3 inputs (XYZ) and 3 (RGB) or 4 (CMYK) outputs.
+ if (b2a->input_channels != ARRAY_COUNT(b2a->input_curves)) {
+ return false;
+ }
+ if (b2a->output_channels < 3 || b2a->output_channels > ARRAY_COUNT(b2a->output_curves)) {
+ return false;
+ }
+
+ uint32_t b_curve_offset = read_big_u32(mBATag->b_curve_offset);
+ uint32_t matrix_offset = read_big_u32(mBATag->matrix_offset);
+ uint32_t m_curve_offset = read_big_u32(mBATag->m_curve_offset);
+ uint32_t clut_offset = read_big_u32(mBATag->clut_offset);
+ uint32_t a_curve_offset = read_big_u32(mBATag->a_curve_offset);
+
+ if (0 == b_curve_offset) {
+ return false;
+ }
+
+ // "B" curves are our inputs, not outputs.
+ if (!read_curves(tag->buf, tag->size, b_curve_offset, b2a->input_channels,
+ b2a->input_curves)) {
+ return false;
+ }
+
+ if (0 != m_curve_offset) {
+ if (0 == matrix_offset) {
+ return false;
+ }
+ // Matrix channels is tied to input_channels (3), not output_channels.
+ b2a->matrix_channels = b2a->input_channels;
+
+ if (!read_curves(tag->buf, tag->size, m_curve_offset, b2a->matrix_channels,
+ b2a->matrix_curves)) {
+ return false;
+ }
+
+ if (tag->size < matrix_offset + 12 * SAFE_SIZEOF(uint32_t)) {
+ return false;
+ }
+ float encoding_factor = pcs_is_xyz ? (32768 / 65535.0f) : 1.0f; // TODO: understand
+ const uint8_t* mtx_buf = tag->buf + matrix_offset;
+ b2a->matrix.vals[0][0] = encoding_factor * read_big_fixed(mtx_buf + 0);
+ b2a->matrix.vals[0][1] = encoding_factor * read_big_fixed(mtx_buf + 4);
+ b2a->matrix.vals[0][2] = encoding_factor * read_big_fixed(mtx_buf + 8);
+ b2a->matrix.vals[1][0] = encoding_factor * read_big_fixed(mtx_buf + 12);
+ b2a->matrix.vals[1][1] = encoding_factor * read_big_fixed(mtx_buf + 16);
+ b2a->matrix.vals[1][2] = encoding_factor * read_big_fixed(mtx_buf + 20);
+ b2a->matrix.vals[2][0] = encoding_factor * read_big_fixed(mtx_buf + 24);
+ b2a->matrix.vals[2][1] = encoding_factor * read_big_fixed(mtx_buf + 28);
+ b2a->matrix.vals[2][2] = encoding_factor * read_big_fixed(mtx_buf + 32);
+ b2a->matrix.vals[0][3] = encoding_factor * read_big_fixed(mtx_buf + 36);
+ b2a->matrix.vals[1][3] = encoding_factor * read_big_fixed(mtx_buf + 40);
+ b2a->matrix.vals[2][3] = encoding_factor * read_big_fixed(mtx_buf + 44);
+ } else {
+ if (0 != matrix_offset) {
+ return false;
+ }
+ b2a->matrix_channels = 0;
+ }
+
+ if (0 != a_curve_offset) {
+ if (0 == clut_offset) {
+ return false;
+ }
+
+ // "A" curves are our output, not input.
+ if (!read_curves(tag->buf, tag->size, a_curve_offset, b2a->output_channels,
+ b2a->output_curves)) {
+ return false;
+ }
+
+ if (tag->size < clut_offset + SAFE_FIXED_SIZE(CLUT_Layout)) {
+ return false;
+ }
+ const CLUT_Layout* clut = (const CLUT_Layout*)(tag->buf + clut_offset);
+
+ if (clut->grid_byte_width[0] == 1) {
+ b2a->grid_8 = clut->variable;
+ b2a->grid_16 = nullptr;
+ } else if (clut->grid_byte_width[0] == 2) {
+ b2a->grid_8 = nullptr;
+ b2a->grid_16 = clut->variable;
+ } else {
+ return false;
+ }
+
+ uint64_t grid_size = b2a->output_channels * clut->grid_byte_width[0];
+ for (uint32_t i = 0; i < b2a->input_channels; ++i) {
+ b2a->grid_points[i] = clut->grid_points[i];
+ if (b2a->grid_points[i] < 2) {
+ return false;
+ }
+ grid_size *= b2a->grid_points[i];
+ }
+ if (tag->size < clut_offset + SAFE_FIXED_SIZE(CLUT_Layout) + grid_size) {
+ return false;
+ }
+ } else {
+ if (0 != clut_offset) {
+ return false;
+ }
+
+ if (b2a->input_channels != b2a->output_channels) {
+ return false;
+ }
+
+ // Zero out *output* channels to skip this stage.
+ b2a->output_channels = 0;
+ }
+ return true;
+}
+
+// If you pass f, we'll fit a possibly-non-zero value for *f.
+// If you pass nullptr, we'll assume you want *f to be treated as zero.
+static int fit_linear(const skcms_Curve* curve, int N, float tol,
+ float* c, float* d, float* f = nullptr) {
+ assert(N > 1);
+ // We iteratively fit the first points to the TF's linear piece.
+ // We want the cx + f line to pass through the first and last points we fit exactly.
+ //
+ // As we walk along the points we find the minimum and maximum slope of the line before the
+ // error would exceed our tolerance. We stop when the range [slope_min, slope_max] becomes
+ // emtpy, when we definitely can't add any more points.
+ //
+ // Some points' error intervals may intersect the running interval but not lie fully
+ // within it. So we keep track of the last point we saw that is a valid end point candidate,
+ // and once the search is done, back up to build the line through *that* point.
+ const float dx = 1.0f / static_cast<float>(N - 1);
+
+ int lin_points = 1;
+
+ float f_zero = 0.0f;
+ if (f) {
+ *f = eval_curve(curve, 0);
+ } else {
+ f = &f_zero;
+ }
+
+
+ float slope_min = -INFINITY_;
+ float slope_max = +INFINITY_;
+ for (int i = 1; i < N; ++i) {
+ float x = static_cast<float>(i) * dx;
+ float y = eval_curve(curve, x);
+
+ float slope_max_i = (y + tol - *f) / x,
+ slope_min_i = (y - tol - *f) / x;
+ if (slope_max_i < slope_min || slope_max < slope_min_i) {
+ // Slope intervals would no longer overlap.
+ break;
+ }
+ slope_max = fminf_(slope_max, slope_max_i);
+ slope_min = fmaxf_(slope_min, slope_min_i);
+
+ float cur_slope = (y - *f) / x;
+ if (slope_min <= cur_slope && cur_slope <= slope_max) {
+ lin_points = i + 1;
+ *c = cur_slope;
+ }
+ }
+
+ // Set D to the last point that met our tolerance.
+ *d = static_cast<float>(lin_points - 1) * dx;
+ return lin_points;
+}
+
+// If this skcms_Curve holds an identity table, rewrite it as an identity skcms_TransferFunction.
+static void canonicalize_identity(skcms_Curve* curve) {
+ if (curve->table_entries && curve->table_entries <= (uint32_t)INT_MAX) {
+ int N = (int)curve->table_entries;
+
+ float c = 0.0f, d = 0.0f, f = 0.0f;
+ if (N == fit_linear(curve, N, 1.0f/static_cast<float>(2*N), &c,&d,&f)
+ && c == 1.0f
+ && f == 0.0f) {
+ curve->table_entries = 0;
+ curve->table_8 = nullptr;
+ curve->table_16 = nullptr;
+ curve->parametric = skcms_TransferFunction{1,1,0,0,0,0,0};
+ }
+ }
+}
+
+static bool read_a2b(const skcms_ICCTag* tag, skcms_A2B* a2b, bool pcs_is_xyz) {
+ bool ok = false;
+ if (tag->type == skcms_Signature_mft1) { ok = read_tag_mft1(tag, a2b); }
+ if (tag->type == skcms_Signature_mft2) { ok = read_tag_mft2(tag, a2b); }
+ if (tag->type == skcms_Signature_mAB ) { ok = read_tag_mab(tag, a2b, pcs_is_xyz); }
+ if (!ok) {
+ return false;
+ }
+
+ if (a2b->input_channels > 0) { canonicalize_identity(a2b->input_curves + 0); }
+ if (a2b->input_channels > 1) { canonicalize_identity(a2b->input_curves + 1); }
+ if (a2b->input_channels > 2) { canonicalize_identity(a2b->input_curves + 2); }
+ if (a2b->input_channels > 3) { canonicalize_identity(a2b->input_curves + 3); }
+
+ if (a2b->matrix_channels > 0) { canonicalize_identity(a2b->matrix_curves + 0); }
+ if (a2b->matrix_channels > 1) { canonicalize_identity(a2b->matrix_curves + 1); }
+ if (a2b->matrix_channels > 2) { canonicalize_identity(a2b->matrix_curves + 2); }
+
+ if (a2b->output_channels > 0) { canonicalize_identity(a2b->output_curves + 0); }
+ if (a2b->output_channels > 1) { canonicalize_identity(a2b->output_curves + 1); }
+ if (a2b->output_channels > 2) { canonicalize_identity(a2b->output_curves + 2); }
+
+ return true;
+}
+
+static bool read_b2a(const skcms_ICCTag* tag, skcms_B2A* b2a, bool pcs_is_xyz) {
+ bool ok = false;
+ if (tag->type == skcms_Signature_mft1) { ok = read_tag_mft1(tag, b2a); }
+ if (tag->type == skcms_Signature_mft2) { ok = read_tag_mft2(tag, b2a); }
+ if (tag->type == skcms_Signature_mBA ) { ok = read_tag_mba(tag, b2a, pcs_is_xyz); }
+ if (!ok) {
+ return false;
+ }
+
+ if (b2a->input_channels > 0) { canonicalize_identity(b2a->input_curves + 0); }
+ if (b2a->input_channels > 1) { canonicalize_identity(b2a->input_curves + 1); }
+ if (b2a->input_channels > 2) { canonicalize_identity(b2a->input_curves + 2); }
+
+ if (b2a->matrix_channels > 0) { canonicalize_identity(b2a->matrix_curves + 0); }
+ if (b2a->matrix_channels > 1) { canonicalize_identity(b2a->matrix_curves + 1); }
+ if (b2a->matrix_channels > 2) { canonicalize_identity(b2a->matrix_curves + 2); }
+
+ if (b2a->output_channels > 0) { canonicalize_identity(b2a->output_curves + 0); }
+ if (b2a->output_channels > 1) { canonicalize_identity(b2a->output_curves + 1); }
+ if (b2a->output_channels > 2) { canonicalize_identity(b2a->output_curves + 2); }
+ if (b2a->output_channels > 3) { canonicalize_identity(b2a->output_curves + 3); }
+
+ return true;
+}
+
+typedef struct {
+ uint8_t type [4];
+ uint8_t reserved [4];
+ uint8_t color_primaries [1];
+ uint8_t transfer_characteristics [1];
+ uint8_t matrix_coefficients [1];
+ uint8_t video_full_range_flag [1];
+} CICP_Layout;
+
+static bool read_cicp(const skcms_ICCTag* tag, skcms_CICP* cicp) {
+ if (tag->type != skcms_Signature_CICP || tag->size < SAFE_SIZEOF(CICP_Layout)) {
+ return false;
+ }
+
+ const CICP_Layout* cicpTag = (const CICP_Layout*)tag->buf;
+
+ cicp->color_primaries = cicpTag->color_primaries[0];
+ cicp->transfer_characteristics = cicpTag->transfer_characteristics[0];
+ cicp->matrix_coefficients = cicpTag->matrix_coefficients[0];
+ cicp->video_full_range_flag = cicpTag->video_full_range_flag[0];
+ return true;
+}
+
+void skcms_GetTagByIndex(const skcms_ICCProfile* profile, uint32_t idx, skcms_ICCTag* tag) {
+ if (!profile || !profile->buffer || !tag) { return; }
+ if (idx > profile->tag_count) { return; }
+ const tag_Layout* tags = get_tag_table(profile);
+ tag->signature = read_big_u32(tags[idx].signature);
+ tag->size = read_big_u32(tags[idx].size);
+ tag->buf = read_big_u32(tags[idx].offset) + profile->buffer;
+ tag->type = read_big_u32(tag->buf);
+}
+
+bool skcms_GetTagBySignature(const skcms_ICCProfile* profile, uint32_t sig, skcms_ICCTag* tag) {
+ if (!profile || !profile->buffer || !tag) { return false; }
+ const tag_Layout* tags = get_tag_table(profile);
+ for (uint32_t i = 0; i < profile->tag_count; ++i) {
+ if (read_big_u32(tags[i].signature) == sig) {
+ tag->signature = sig;
+ tag->size = read_big_u32(tags[i].size);
+ tag->buf = read_big_u32(tags[i].offset) + profile->buffer;
+ tag->type = read_big_u32(tag->buf);
+ return true;
+ }
+ }
+ return false;
+}
+
+static bool usable_as_src(const skcms_ICCProfile* profile) {
+ return profile->has_A2B
+ || (profile->has_trc && profile->has_toXYZD50);
+}
+
+bool skcms_ParseWithA2BPriority(const void* buf, size_t len,
+ const int priority[], const int priorities,
+ skcms_ICCProfile* profile) {
+ assert(SAFE_SIZEOF(header_Layout) == 132);
+
+ if (!profile) {
+ return false;
+ }
+ memset(profile, 0, SAFE_SIZEOF(*profile));
+
+ if (len < SAFE_SIZEOF(header_Layout)) {
+ return false;
+ }
+
+ // Byte-swap all header fields
+ const header_Layout* header = (const header_Layout*)buf;
+ profile->buffer = (const uint8_t*)buf;
+ profile->size = read_big_u32(header->size);
+ uint32_t version = read_big_u32(header->version);
+ profile->data_color_space = read_big_u32(header->data_color_space);
+ profile->pcs = read_big_u32(header->pcs);
+ uint32_t signature = read_big_u32(header->signature);
+ float illuminant_X = read_big_fixed(header->illuminant_X);
+ float illuminant_Y = read_big_fixed(header->illuminant_Y);
+ float illuminant_Z = read_big_fixed(header->illuminant_Z);
+ profile->tag_count = read_big_u32(header->tag_count);
+
+ // Validate signature, size (smaller than buffer, large enough to hold tag table),
+ // and major version
+ uint64_t tag_table_size = profile->tag_count * SAFE_SIZEOF(tag_Layout);
+ if (signature != skcms_Signature_acsp ||
+ profile->size > len ||
+ profile->size < SAFE_SIZEOF(header_Layout) + tag_table_size ||
+ (version >> 24) > 4) {
+ return false;
+ }
+
+ // Validate that illuminant is D50 white
+ if (fabsf_(illuminant_X - 0.9642f) > 0.0100f ||
+ fabsf_(illuminant_Y - 1.0000f) > 0.0100f ||
+ fabsf_(illuminant_Z - 0.8249f) > 0.0100f) {
+ return false;
+ }
+
+ // Validate that all tag entries have sane offset + size
+ const tag_Layout* tags = get_tag_table(profile);
+ for (uint32_t i = 0; i < profile->tag_count; ++i) {
+ uint32_t tag_offset = read_big_u32(tags[i].offset);
+ uint32_t tag_size = read_big_u32(tags[i].size);
+ uint64_t tag_end = (uint64_t)tag_offset + (uint64_t)tag_size;
+ if (tag_size < 4 || tag_end > profile->size) {
+ return false;
+ }
+ }
+
+ if (profile->pcs != skcms_Signature_XYZ && profile->pcs != skcms_Signature_Lab) {
+ return false;
+ }
+
+ bool pcs_is_xyz = profile->pcs == skcms_Signature_XYZ;
+
+ // Pre-parse commonly used tags.
+ skcms_ICCTag kTRC;
+ if (profile->data_color_space == skcms_Signature_Gray &&
+ skcms_GetTagBySignature(profile, skcms_Signature_kTRC, &kTRC)) {
+ if (!read_curve(kTRC.buf, kTRC.size, &profile->trc[0], nullptr)) {
+ // Malformed tag
+ return false;
+ }
+ profile->trc[1] = profile->trc[0];
+ profile->trc[2] = profile->trc[0];
+ profile->has_trc = true;
+
+ if (pcs_is_xyz) {
+ profile->toXYZD50.vals[0][0] = illuminant_X;
+ profile->toXYZD50.vals[1][1] = illuminant_Y;
+ profile->toXYZD50.vals[2][2] = illuminant_Z;
+ profile->has_toXYZD50 = true;
+ }
+ } else {
+ skcms_ICCTag rTRC, gTRC, bTRC;
+ if (skcms_GetTagBySignature(profile, skcms_Signature_rTRC, &rTRC) &&
+ skcms_GetTagBySignature(profile, skcms_Signature_gTRC, &gTRC) &&
+ skcms_GetTagBySignature(profile, skcms_Signature_bTRC, &bTRC)) {
+ if (!read_curve(rTRC.buf, rTRC.size, &profile->trc[0], nullptr) ||
+ !read_curve(gTRC.buf, gTRC.size, &profile->trc[1], nullptr) ||
+ !read_curve(bTRC.buf, bTRC.size, &profile->trc[2], nullptr)) {
+ // Malformed TRC tags
+ return false;
+ }
+ profile->has_trc = true;
+ }
+
+ skcms_ICCTag rXYZ, gXYZ, bXYZ;
+ if (skcms_GetTagBySignature(profile, skcms_Signature_rXYZ, &rXYZ) &&
+ skcms_GetTagBySignature(profile, skcms_Signature_gXYZ, &gXYZ) &&
+ skcms_GetTagBySignature(profile, skcms_Signature_bXYZ, &bXYZ)) {
+ if (!read_to_XYZD50(&rXYZ, &gXYZ, &bXYZ, &profile->toXYZD50)) {
+ // Malformed XYZ tags
+ return false;
+ }
+ profile->has_toXYZD50 = true;
+ }
+ }
+
+ for (int i = 0; i < priorities; i++) {
+ // enum { perceptual, relative_colormetric, saturation }
+ if (priority[i] < 0 || priority[i] > 2) {
+ return false;
+ }
+ uint32_t sig = skcms_Signature_A2B0 + static_cast<uint32_t>(priority[i]);
+ skcms_ICCTag tag;
+ if (skcms_GetTagBySignature(profile, sig, &tag)) {
+ if (!read_a2b(&tag, &profile->A2B, pcs_is_xyz)) {
+ // Malformed A2B tag
+ return false;
+ }
+ profile->has_A2B = true;
+ break;
+ }
+ }
+
+ for (int i = 0; i < priorities; i++) {
+ // enum { perceptual, relative_colormetric, saturation }
+ if (priority[i] < 0 || priority[i] > 2) {
+ return false;
+ }
+ uint32_t sig = skcms_Signature_B2A0 + static_cast<uint32_t>(priority[i]);
+ skcms_ICCTag tag;
+ if (skcms_GetTagBySignature(profile, sig, &tag)) {
+ if (!read_b2a(&tag, &profile->B2A, pcs_is_xyz)) {
+ // Malformed B2A tag
+ return false;
+ }
+ profile->has_B2A = true;
+ break;
+ }
+ }
+
+ skcms_ICCTag cicp_tag;
+ if (skcms_GetTagBySignature(profile, skcms_Signature_CICP, &cicp_tag)) {
+ if (!read_cicp(&cicp_tag, &profile->CICP)) {
+ // Malformed CICP tag
+ return false;
+ }
+ profile->has_CICP = true;
+ }
+
+ return usable_as_src(profile);
+}
+
+
+const skcms_ICCProfile* skcms_sRGB_profile() {
+ static const skcms_ICCProfile sRGB_profile = {
+ nullptr, // buffer, moot here
+
+ 0, // size, moot here
+ skcms_Signature_RGB, // data_color_space
+ skcms_Signature_XYZ, // pcs
+ 0, // tag count, moot here
+
+ // We choose to represent sRGB with its canonical transfer function,
+ // and with its canonical XYZD50 gamut matrix.
+ true, // has_trc, followed by the 3 trc curves
+ {
+ {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
+ {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
+ {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
+ },
+
+ true, // has_toXYZD50, followed by 3x3 toXYZD50 matrix
+ {{
+ { 0.436065674f, 0.385147095f, 0.143066406f },
+ { 0.222488403f, 0.716873169f, 0.060607910f },
+ { 0.013916016f, 0.097076416f, 0.714096069f },
+ }},
+
+ false, // has_A2B, followed by A2B itself, which we don't care about.
+ {
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ {0,0,0,0},
+ nullptr,
+ nullptr,
+
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ {{
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ }},
+
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ },
+
+ false, // has_B2A, followed by B2A itself, which we also don't care about.
+ {
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+
+ 0,
+ {{
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ }},
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+
+ 0,
+ {0,0,0,0},
+ nullptr,
+ nullptr,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ },
+
+ false, // has_CICP, followed by cicp itself which we don't care about.
+ { 0, 0, 0, 0 },
+ };
+ return &sRGB_profile;
+}
+
+const skcms_ICCProfile* skcms_XYZD50_profile() {
+ // Just like sRGB above, but with identity transfer functions and toXYZD50 matrix.
+ static const skcms_ICCProfile XYZD50_profile = {
+ nullptr, // buffer, moot here
+
+ 0, // size, moot here
+ skcms_Signature_RGB, // data_color_space
+ skcms_Signature_XYZ, // pcs
+ 0, // tag count, moot here
+
+ true, // has_trc, followed by the 3 trc curves
+ {
+ {{0, {1,1, 0,0,0,0,0}}},
+ {{0, {1,1, 0,0,0,0,0}}},
+ {{0, {1,1, 0,0,0,0,0}}},
+ },
+
+ true, // has_toXYZD50, followed by 3x3 toXYZD50 matrix
+ {{
+ { 1,0,0 },
+ { 0,1,0 },
+ { 0,0,1 },
+ }},
+
+ false, // has_A2B, followed by A2B itself, which we don't care about.
+ {
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ {0,0,0,0},
+ nullptr,
+ nullptr,
+
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ {{
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ }},
+
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ },
+
+ false, // has_B2A, followed by B2A itself, which we also don't care about.
+ {
+ 0,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+
+ 0,
+ {{
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ { 0,0,0,0 },
+ }},
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+
+ 0,
+ {0,0,0,0},
+ nullptr,
+ nullptr,
+ {
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ {{0, {0,0, 0,0,0,0,0}}},
+ },
+ },
+
+ false, // has_CICP, followed by cicp itself which we don't care about.
+ { 0, 0, 0, 0 },
+ };
+
+ return &XYZD50_profile;
+}
+
+const skcms_TransferFunction* skcms_sRGB_TransferFunction() {
+ return &skcms_sRGB_profile()->trc[0].parametric;
+}
+
+const skcms_TransferFunction* skcms_sRGB_Inverse_TransferFunction() {
+ static const skcms_TransferFunction sRGB_inv =
+ {0.416666657f, 1.137283325f, -0.0f, 12.920000076f, 0.003130805f, -0.054969788f, -0.0f};
+ return &sRGB_inv;
+}
+
+const skcms_TransferFunction* skcms_Identity_TransferFunction() {
+ static const skcms_TransferFunction identity = {1,1,0,0,0,0,0};
+ return &identity;
+}
+
+const uint8_t skcms_252_random_bytes[] = {
+ 8, 179, 128, 204, 253, 38, 134, 184, 68, 102, 32, 138, 99, 39, 169, 215,
+ 119, 26, 3, 223, 95, 239, 52, 132, 114, 74, 81, 234, 97, 116, 244, 205, 30,
+ 154, 173, 12, 51, 159, 122, 153, 61, 226, 236, 178, 229, 55, 181, 220, 191,
+ 194, 160, 126, 168, 82, 131, 18, 180, 245, 163, 22, 246, 69, 235, 252, 57,
+ 108, 14, 6, 152, 240, 255, 171, 242, 20, 227, 177, 238, 96, 85, 16, 211,
+ 70, 200, 149, 155, 146, 127, 145, 100, 151, 109, 19, 165, 208, 195, 164,
+ 137, 254, 182, 248, 64, 201, 45, 209, 5, 147, 207, 210, 113, 162, 83, 225,
+ 9, 31, 15, 231, 115, 37, 58, 53, 24, 49, 197, 56, 120, 172, 48, 21, 214,
+ 129, 111, 11, 50, 187, 196, 34, 60, 103, 71, 144, 47, 203, 77, 80, 232,
+ 140, 222, 250, 206, 166, 247, 139, 249, 221, 72, 106, 27, 199, 117, 54,
+ 219, 135, 118, 40, 79, 41, 251, 46, 93, 212, 92, 233, 148, 28, 121, 63,
+ 123, 158, 105, 59, 29, 42, 143, 23, 0, 107, 176, 87, 104, 183, 156, 193,
+ 189, 90, 188, 65, 190, 17, 198, 7, 186, 161, 1, 124, 78, 125, 170, 133,
+ 174, 218, 67, 157, 75, 101, 89, 217, 62, 33, 141, 228, 25, 35, 91, 230, 4,
+ 2, 13, 73, 86, 167, 237, 84, 243, 44, 185, 66, 130, 110, 150, 142, 216, 88,
+ 112, 36, 224, 136, 202, 76, 94, 98, 175, 213
+};
+
+bool skcms_ApproximatelyEqualProfiles(const skcms_ICCProfile* A, const skcms_ICCProfile* B) {
+ // Test for exactly equal profiles first.
+ if (A == B || 0 == memcmp(A,B, sizeof(skcms_ICCProfile))) {
+ return true;
+ }
+
+ // For now this is the essentially the same strategy we use in test_only.c
+ // for our skcms_Transform() smoke tests:
+ // 1) transform A to XYZD50
+ // 2) transform B to XYZD50
+ // 3) return true if they're similar enough
+ // Our current criterion in 3) is maximum 1 bit error per XYZD50 byte.
+
+ // skcms_252_random_bytes are 252 of a random shuffle of all possible bytes.
+ // 252 is evenly divisible by 3 and 4. Only 192, 10, 241, and 43 are missing.
+
+ // We want to allow otherwise equivalent profiles tagged as grayscale and RGB
+ // to be treated as equal. But CMYK profiles are a totally different ballgame.
+ const auto CMYK = skcms_Signature_CMYK;
+ if ((A->data_color_space == CMYK) != (B->data_color_space == CMYK)) {
+ return false;
+ }
+
+ // Interpret as RGB_888 if data color space is RGB or GRAY, RGBA_8888 if CMYK.
+ // TODO: working with RGBA_8888 either way is probably fastest.
+ skcms_PixelFormat fmt = skcms_PixelFormat_RGB_888;
+ size_t npixels = 84;
+ if (A->data_color_space == skcms_Signature_CMYK) {
+ fmt = skcms_PixelFormat_RGBA_8888;
+ npixels = 63;
+ }
+
+ // TODO: if A or B is a known profile (skcms_sRGB_profile, skcms_XYZD50_profile),
+ // use pre-canned results and skip that skcms_Transform() call?
+ uint8_t dstA[252],
+ dstB[252];
+ if (!skcms_Transform(
+ skcms_252_random_bytes, fmt, skcms_AlphaFormat_Unpremul, A,
+ dstA, skcms_PixelFormat_RGB_888, skcms_AlphaFormat_Unpremul, skcms_XYZD50_profile(),
+ npixels)) {
+ return false;
+ }
+ if (!skcms_Transform(
+ skcms_252_random_bytes, fmt, skcms_AlphaFormat_Unpremul, B,
+ dstB, skcms_PixelFormat_RGB_888, skcms_AlphaFormat_Unpremul, skcms_XYZD50_profile(),
+ npixels)) {
+ return false;
+ }
+
+ // TODO: make sure this final check has reasonable codegen.
+ for (size_t i = 0; i < 252; i++) {
+ if (abs((int)dstA[i] - (int)dstB[i]) > 1) {
+ return false;
+ }
+ }
+ return true;
+}
+
+bool skcms_TRCs_AreApproximateInverse(const skcms_ICCProfile* profile,
+ const skcms_TransferFunction* inv_tf) {
+ if (!profile || !profile->has_trc) {
+ return false;
+ }
+
+ return skcms_AreApproximateInverses(&profile->trc[0], inv_tf) &&
+ skcms_AreApproximateInverses(&profile->trc[1], inv_tf) &&
+ skcms_AreApproximateInverses(&profile->trc[2], inv_tf);
+}
+
+static bool is_zero_to_one(float x) {
+ return 0 <= x && x <= 1;
+}
+
+typedef struct { float vals[3]; } skcms_Vector3;
+
+static skcms_Vector3 mv_mul(const skcms_Matrix3x3* m, const skcms_Vector3* v) {
+ skcms_Vector3 dst = {{0,0,0}};
+ for (int row = 0; row < 3; ++row) {
+ dst.vals[row] = m->vals[row][0] * v->vals[0]
+ + m->vals[row][1] * v->vals[1]
+ + m->vals[row][2] * v->vals[2];
+ }
+ return dst;
+}
+
+bool skcms_AdaptToXYZD50(float wx, float wy,
+ skcms_Matrix3x3* toXYZD50) {
+ if (!is_zero_to_one(wx) || !is_zero_to_one(wy) ||
+ !toXYZD50) {
+ return false;
+ }
+
+ // Assumes that Y is 1.0f.
+ skcms_Vector3 wXYZ = { { wx / wy, 1, (1 - wx - wy) / wy } };
+
+ // Now convert toXYZ matrix to toXYZD50.
+ skcms_Vector3 wXYZD50 = { { 0.96422f, 1.0f, 0.82521f } };
+
+ // Calculate the chromatic adaptation matrix. We will use the Bradford method, thus
+ // the matrices below. The Bradford method is used by Adobe and is widely considered
+ // to be the best.
+ skcms_Matrix3x3 xyz_to_lms = {{
+ { 0.8951f, 0.2664f, -0.1614f },
+ { -0.7502f, 1.7135f, 0.0367f },
+ { 0.0389f, -0.0685f, 1.0296f },
+ }};
+ skcms_Matrix3x3 lms_to_xyz = {{
+ { 0.9869929f, -0.1470543f, 0.1599627f },
+ { 0.4323053f, 0.5183603f, 0.0492912f },
+ { -0.0085287f, 0.0400428f, 0.9684867f },
+ }};
+
+ skcms_Vector3 srcCone = mv_mul(&xyz_to_lms, &wXYZ);
+ skcms_Vector3 dstCone = mv_mul(&xyz_to_lms, &wXYZD50);
+
+ *toXYZD50 = {{
+ { dstCone.vals[0] / srcCone.vals[0], 0, 0 },
+ { 0, dstCone.vals[1] / srcCone.vals[1], 0 },
+ { 0, 0, dstCone.vals[2] / srcCone.vals[2] },
+ }};
+ *toXYZD50 = skcms_Matrix3x3_concat(toXYZD50, &xyz_to_lms);
+ *toXYZD50 = skcms_Matrix3x3_concat(&lms_to_xyz, toXYZD50);
+
+ return true;
+}
+
+bool skcms_PrimariesToXYZD50(float rx, float ry,
+ float gx, float gy,
+ float bx, float by,
+ float wx, float wy,
+ skcms_Matrix3x3* toXYZD50) {
+ if (!is_zero_to_one(rx) || !is_zero_to_one(ry) ||
+ !is_zero_to_one(gx) || !is_zero_to_one(gy) ||
+ !is_zero_to_one(bx) || !is_zero_to_one(by) ||
+ !is_zero_to_one(wx) || !is_zero_to_one(wy) ||
+ !toXYZD50) {
+ return false;
+ }
+
+ // First, we need to convert xy values (primaries) to XYZ.
+ skcms_Matrix3x3 primaries = {{
+ { rx, gx, bx },
+ { ry, gy, by },
+ { 1 - rx - ry, 1 - gx - gy, 1 - bx - by },
+ }};
+ skcms_Matrix3x3 primaries_inv;
+ if (!skcms_Matrix3x3_invert(&primaries, &primaries_inv)) {
+ return false;
+ }
+
+ // Assumes that Y is 1.0f.
+ skcms_Vector3 wXYZ = { { wx / wy, 1, (1 - wx - wy) / wy } };
+ skcms_Vector3 XYZ = mv_mul(&primaries_inv, &wXYZ);
+
+ skcms_Matrix3x3 toXYZ = {{
+ { XYZ.vals[0], 0, 0 },
+ { 0, XYZ.vals[1], 0 },
+ { 0, 0, XYZ.vals[2] },
+ }};
+ toXYZ = skcms_Matrix3x3_concat(&primaries, &toXYZ);
+
+ skcms_Matrix3x3 DXtoD50;
+ if (!skcms_AdaptToXYZD50(wx, wy, &DXtoD50)) {
+ return false;
+ }
+
+ *toXYZD50 = skcms_Matrix3x3_concat(&DXtoD50, &toXYZ);
+ return true;
+}
+
+
+bool skcms_Matrix3x3_invert(const skcms_Matrix3x3* src, skcms_Matrix3x3* dst) {
+ double a00 = src->vals[0][0],
+ a01 = src->vals[1][0],
+ a02 = src->vals[2][0],
+ a10 = src->vals[0][1],
+ a11 = src->vals[1][1],
+ a12 = src->vals[2][1],
+ a20 = src->vals[0][2],
+ a21 = src->vals[1][2],
+ a22 = src->vals[2][2];
+
+ double b0 = a00*a11 - a01*a10,
+ b1 = a00*a12 - a02*a10,
+ b2 = a01*a12 - a02*a11,
+ b3 = a20,
+ b4 = a21,
+ b5 = a22;
+
+ double determinant = b0*b5
+ - b1*b4
+ + b2*b3;
+
+ if (determinant == 0) {
+ return false;
+ }
+
+ double invdet = 1.0 / determinant;
+ if (invdet > +FLT_MAX || invdet < -FLT_MAX || !isfinitef_((float)invdet)) {
+ return false;
+ }
+
+ b0 *= invdet;
+ b1 *= invdet;
+ b2 *= invdet;
+ b3 *= invdet;
+ b4 *= invdet;
+ b5 *= invdet;
+
+ dst->vals[0][0] = (float)( a11*b5 - a12*b4 );
+ dst->vals[1][0] = (float)( a02*b4 - a01*b5 );
+ dst->vals[2][0] = (float)( + b2 );
+ dst->vals[0][1] = (float)( a12*b3 - a10*b5 );
+ dst->vals[1][1] = (float)( a00*b5 - a02*b3 );
+ dst->vals[2][1] = (float)( - b1 );
+ dst->vals[0][2] = (float)( a10*b4 - a11*b3 );
+ dst->vals[1][2] = (float)( a01*b3 - a00*b4 );
+ dst->vals[2][2] = (float)( + b0 );
+
+ for (int r = 0; r < 3; ++r)
+ for (int c = 0; c < 3; ++c) {
+ if (!isfinitef_(dst->vals[r][c])) {
+ return false;
+ }
+ }
+ return true;
+}
+
+skcms_Matrix3x3 skcms_Matrix3x3_concat(const skcms_Matrix3x3* A, const skcms_Matrix3x3* B) {
+ skcms_Matrix3x3 m = { { { 0,0,0 },{ 0,0,0 },{ 0,0,0 } } };
+ for (int r = 0; r < 3; r++)
+ for (int c = 0; c < 3; c++) {
+ m.vals[r][c] = A->vals[r][0] * B->vals[0][c]
+ + A->vals[r][1] * B->vals[1][c]
+ + A->vals[r][2] * B->vals[2][c];
+ }
+ return m;
+}
+
+#if defined(__clang__)
+ [[clang::no_sanitize("float-divide-by-zero")]] // Checked for by classify() on the way out.
+#endif
+bool skcms_TransferFunction_invert(const skcms_TransferFunction* src, skcms_TransferFunction* dst) {
+ TF_PQish pq;
+ TF_HLGish hlg;
+ switch (classify(*src, &pq, &hlg)) {
+ case skcms_TFType_Invalid: return false;
+ case skcms_TFType_sRGBish: break; // handled below
+
+ case skcms_TFType_PQish:
+ *dst = { TFKind_marker(skcms_TFType_PQish), -pq.A, pq.D, 1.0f/pq.F
+ , pq.B, -pq.E, 1.0f/pq.C};
+ return true;
+
+ case skcms_TFType_HLGish:
+ *dst = { TFKind_marker(skcms_TFType_HLGinvish), 1.0f/hlg.R, 1.0f/hlg.G
+ , 1.0f/hlg.a, hlg.b, hlg.c
+ , hlg.K_minus_1 };
+ return true;
+
+ case skcms_TFType_HLGinvish:
+ *dst = { TFKind_marker(skcms_TFType_HLGish), 1.0f/hlg.R, 1.0f/hlg.G
+ , 1.0f/hlg.a, hlg.b, hlg.c
+ , hlg.K_minus_1 };
+ return true;
+ }
+
+ assert (classify(*src) == skcms_TFType_sRGBish);
+
+ // We're inverting this function, solving for x in terms of y.
+ // y = (cx + f) x < d
+ // (ax + b)^g + e x ≥ d
+ // The inverse of this function can be expressed in the same piecewise form.
+ skcms_TransferFunction inv = {0,0,0,0,0,0,0};
+
+ // We'll start by finding the new threshold inv.d.
+ // In principle we should be able to find that by solving for y at x=d from either side.
+ // (If those two d values aren't the same, it's a discontinuous transfer function.)
+ float d_l = src->c * src->d + src->f,
+ d_r = powf_(src->a * src->d + src->b, src->g) + src->e;
+ if (fabsf_(d_l - d_r) > 1/512.0f) {
+ return false;
+ }
+ inv.d = d_l; // TODO(mtklein): better in practice to choose d_r?
+
+ // When d=0, the linear section collapses to a point. We leave c,d,f all zero in that case.
+ if (inv.d > 0) {
+ // Inverting the linear section is pretty straightfoward:
+ // y = cx + f
+ // y - f = cx
+ // (1/c)y - f/c = x
+ inv.c = 1.0f/src->c;
+ inv.f = -src->f/src->c;
+ }
+
+ // The interesting part is inverting the nonlinear section:
+ // y = (ax + b)^g + e.
+ // y - e = (ax + b)^g
+ // (y - e)^1/g = ax + b
+ // (y - e)^1/g - b = ax
+ // (1/a)(y - e)^1/g - b/a = x
+ //
+ // To make that fit our form, we need to move the (1/a) term inside the exponentiation:
+ // let k = (1/a)^g
+ // (1/a)( y - e)^1/g - b/a = x
+ // (ky - ke)^1/g - b/a = x
+
+ float k = powf_(src->a, -src->g); // (1/a)^g == a^-g
+ inv.g = 1.0f / src->g;
+ inv.a = k;
+ inv.b = -k * src->e;
+ inv.e = -src->b / src->a;
+
+ // We need to enforce the same constraints here that we do when fitting a curve,
+ // a >= 0 and ad+b >= 0. These constraints are checked by classify(), so they're true
+ // of the source function if we're here.
+
+ // Just like when fitting the curve, there's really no way to rescue a < 0.
+ if (inv.a < 0) {
+ return false;
+ }
+ // On the other hand we can rescue an ad+b that's gone slightly negative here.
+ if (inv.a * inv.d + inv.b < 0) {
+ inv.b = -inv.a * inv.d;
+ }
+
+ // That should usually make classify(inv) == sRGBish true, but there are a couple situations
+ // where we might still fail here, like non-finite parameter values.
+ if (classify(inv) != skcms_TFType_sRGBish) {
+ return false;
+ }
+
+ assert (inv.a >= 0);
+ assert (inv.a * inv.d + inv.b >= 0);
+
+ // Now in principle we're done.
+ // But to preserve the valuable invariant inv(src(1.0f)) == 1.0f, we'll tweak
+ // e or f of the inverse, depending on which segment contains src(1.0f).
+ float s = skcms_TransferFunction_eval(src, 1.0f);
+ if (!isfinitef_(s)) {
+ return false;
+ }
+
+ float sign = s < 0 ? -1.0f : 1.0f;
+ s *= sign;
+ if (s < inv.d) {
+ inv.f = 1.0f - sign * inv.c * s;
+ } else {
+ inv.e = 1.0f - sign * powf_(inv.a * s + inv.b, inv.g);
+ }
+
+ *dst = inv;
+ return classify(*dst) == skcms_TFType_sRGBish;
+}
+
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ //
+
+// From here below we're approximating an skcms_Curve with an skcms_TransferFunction{g,a,b,c,d,e,f}:
+//
+// tf(x) = cx + f x < d
+// tf(x) = (ax + b)^g + e x ≥ d
+//
+// When fitting, we add the additional constraint that both pieces meet at d:
+//
+// cd + f = (ad + b)^g + e
+//
+// Solving for e and folding it through gives an alternate formulation of the non-linear piece:
+//
+// tf(x) = cx + f x < d
+// tf(x) = (ax + b)^g - (ad + b)^g + cd + f x ≥ d
+//
+// Our overall strategy is then:
+// For a couple tolerances,
+// - fit_linear(): fit c,d,f iteratively to as many points as our tolerance allows
+// - invert c,d,f
+// - fit_nonlinear(): fit g,a,b using Gauss-Newton given those inverted c,d,f
+// (and by constraint, inverted e) to the inverse of the table.
+// Return the parameters with least maximum error.
+//
+// To run Gauss-Newton to find g,a,b, we'll also need the gradient of the residuals
+// of round-trip f_inv(x), the inverse of the non-linear piece of f(x).
+//
+// let y = Table(x)
+// r(x) = x - f_inv(y)
+//
+// ∂r/∂g = ln(ay + b)*(ay + b)^g
+// - ln(ad + b)*(ad + b)^g
+// ∂r/∂a = yg(ay + b)^(g-1)
+// - dg(ad + b)^(g-1)
+// ∂r/∂b = g(ay + b)^(g-1)
+// - g(ad + b)^(g-1)
+
+// Return the residual of roundtripping skcms_Curve(x) through f_inv(y) with parameters P,
+// and fill out the gradient of the residual into dfdP.
+static float rg_nonlinear(float x,
+ const skcms_Curve* curve,
+ const skcms_TransferFunction* tf,
+ float dfdP[3]) {
+ const float y = eval_curve(curve, x);
+
+ const float g = tf->g, a = tf->a, b = tf->b,
+ c = tf->c, d = tf->d, f = tf->f;
+
+ const float Y = fmaxf_(a*y + b, 0.0f),
+ D = a*d + b;
+ assert (D >= 0);
+
+ // The gradient.
+ dfdP[0] = logf_(Y)*powf_(Y, g)
+ - logf_(D)*powf_(D, g);
+ dfdP[1] = y*g*powf_(Y, g-1)
+ - d*g*powf_(D, g-1);
+ dfdP[2] = g*powf_(Y, g-1)
+ - g*powf_(D, g-1);
+
+ // The residual.
+ const float f_inv = powf_(Y, g)
+ - powf_(D, g)
+ + c*d + f;
+ return x - f_inv;
+}
+
+static bool gauss_newton_step(const skcms_Curve* curve,
+ skcms_TransferFunction* tf,
+ float x0, float dx, int N) {
+ // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing.
+ //
+ // Let P = [ tf->g, tf->a, tf->b ] (the three terms that we're adjusting).
+ //
+ // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P),
+ // where r(P) is the residual vector
+ // and Jf is the Jacobian matrix of f(), ∂r/∂P.
+ //
+ // Let's review the shape of each of these expressions:
+ // r(P) is [N x 1], a column vector with one entry per value of x tested
+ // Jf is [N x 3], a matrix with an entry for each (x,P) pair
+ // Jf^T is [3 x N], the transpose of Jf
+ //
+ // Jf^T Jf is [3 x N] * [N x 3] == [3 x 3], a 3x3 matrix,
+ // and so is its inverse (Jf^T Jf)^-1
+ // Jf^T r(P) is [3 x N] * [N x 1] == [3 x 1], a column vector with the same shape as P
+ //
+ // Our implementation strategy to get to the final ∆P is
+ // 1) evaluate Jf^T Jf, call that lhs
+ // 2) evaluate Jf^T r(P), call that rhs
+ // 3) invert lhs
+ // 4) multiply inverse lhs by rhs
+ //
+ // This is a friendly implementation strategy because we don't have to have any
+ // buffers that scale with N, and equally nice don't have to perform any matrix
+ // operations that are variable size.
+ //
+ // Other implementation strategies could trade this off, e.g. evaluating the
+ // pseudoinverse of Jf ( (Jf^T Jf)^-1 Jf^T ) directly, then multiplying that by
+ // the residuals. That would probably require implementing singular value
+ // decomposition, and would create a [3 x N] matrix to be multiplied by the
+ // [N x 1] residual vector, but on the upside I think that'd eliminate the
+ // possibility of this gauss_newton_step() function ever failing.
+
+ // 0) start off with lhs and rhs safely zeroed.
+ skcms_Matrix3x3 lhs = {{ {0,0,0}, {0,0,0}, {0,0,0} }};
+ skcms_Vector3 rhs = { {0,0,0} };
+
+ // 1,2) evaluate lhs and evaluate rhs
+ // We want to evaluate Jf only once, but both lhs and rhs involve Jf^T,
+ // so we'll have to update lhs and rhs at the same time.
+ for (int i = 0; i < N; i++) {
+ float x = x0 + static_cast<float>(i)*dx;
+
+ float dfdP[3] = {0,0,0};
+ float resid = rg_nonlinear(x,curve,tf, dfdP);
+
+ for (int r = 0; r < 3; r++) {
+ for (int c = 0; c < 3; c++) {
+ lhs.vals[r][c] += dfdP[r] * dfdP[c];
+ }
+ rhs.vals[r] += dfdP[r] * resid;
+ }
+ }
+
+ // If any of the 3 P parameters are unused, this matrix will be singular.
+ // Detect those cases and fix them up to indentity instead, so we can invert.
+ for (int k = 0; k < 3; k++) {
+ if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 &&
+ lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0) {
+ lhs.vals[k][k] = 1;
+ }
+ }
+
+ // 3) invert lhs
+ skcms_Matrix3x3 lhs_inv;
+ if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) {
+ return false;
+ }
+
+ // 4) multiply inverse lhs by rhs
+ skcms_Vector3 dP = mv_mul(&lhs_inv, &rhs);
+ tf->g += dP.vals[0];
+ tf->a += dP.vals[1];
+ tf->b += dP.vals[2];
+ return isfinitef_(tf->g) && isfinitef_(tf->a) && isfinitef_(tf->b);
+}
+
+static float max_roundtrip_error_checked(const skcms_Curve* curve,
+ const skcms_TransferFunction* tf_inv) {
+ skcms_TransferFunction tf;
+ if (!skcms_TransferFunction_invert(tf_inv, &tf) || skcms_TFType_sRGBish != classify(tf)) {
+ return INFINITY_;
+ }
+
+ skcms_TransferFunction tf_inv_again;
+ if (!skcms_TransferFunction_invert(&tf, &tf_inv_again)) {
+ return INFINITY_;
+ }
+
+ return skcms_MaxRoundtripError(curve, &tf_inv_again);
+}
+
+// Fit the points in [L,N) to the non-linear piece of tf, or return false if we can't.
+static bool fit_nonlinear(const skcms_Curve* curve, int L, int N, skcms_TransferFunction* tf) {
+ // This enforces a few constraints that are not modeled in gauss_newton_step()'s optimization.
+ auto fixup_tf = [tf]() {
+ // a must be non-negative. That ensures the function is monotonically increasing.
+ // We don't really know how to fix up a if it goes negative.
+ if (tf->a < 0) {
+ return false;
+ }
+ // ad+b must be non-negative. That ensures we don't end up with complex numbers in powf.
+ // We feel just barely not uneasy enough to tweak b so ad+b is zero in this case.
+ if (tf->a * tf->d + tf->b < 0) {
+ tf->b = -tf->a * tf->d;
+ }
+ assert (tf->a >= 0 &&
+ tf->a * tf->d + tf->b >= 0);
+
+ // cd+f must be ~= (ad+b)^g+e. That ensures the function is continuous. We keep e as a free
+ // parameter so we can guarantee this.
+ tf->e = tf->c*tf->d + tf->f
+ - powf_(tf->a*tf->d + tf->b, tf->g);
+
+ return true;
+ };
+
+ if (!fixup_tf()) {
+ return false;
+ }
+
+ // No matter where we start, dx should always represent N even steps from 0 to 1.
+ const float dx = 1.0f / static_cast<float>(N-1);
+
+ skcms_TransferFunction best_tf = *tf;
+ float best_max_error = INFINITY_;
+
+ // Need this or several curves get worse... *sigh*
+ float init_error = max_roundtrip_error_checked(curve, tf);
+ if (init_error < best_max_error) {
+ best_max_error = init_error;
+ best_tf = *tf;
+ }
+
+ // As far as we can tell, 1 Gauss-Newton step won't converge, and 3 steps is no better than 2.
+ for (int j = 0; j < 8; j++) {
+ if (!gauss_newton_step(curve, tf, static_cast<float>(L)*dx, dx, N-L) || !fixup_tf()) {
+ *tf = best_tf;
+ return isfinitef_(best_max_error);
+ }
+
+ float max_error = max_roundtrip_error_checked(curve, tf);
+ if (max_error < best_max_error) {
+ best_max_error = max_error;
+ best_tf = *tf;
+ }
+ }
+
+ *tf = best_tf;
+ return isfinitef_(best_max_error);
+}
+
+bool skcms_ApproximateCurve(const skcms_Curve* curve,
+ skcms_TransferFunction* approx,
+ float* max_error) {
+ if (!curve || !approx || !max_error) {
+ return false;
+ }
+
+ if (curve->table_entries == 0) {
+ // No point approximating an skcms_TransferFunction with an skcms_TransferFunction!
+ return false;
+ }
+
+ if (curve->table_entries == 1 || curve->table_entries > (uint32_t)INT_MAX) {
+ // We need at least two points, and must put some reasonable cap on the maximum number.
+ return false;
+ }
+
+ int N = (int)curve->table_entries;
+ const float dx = 1.0f / static_cast<float>(N - 1);
+
+ *max_error = INFINITY_;
+ const float kTolerances[] = { 1.5f / 65535.0f, 1.0f / 512.0f };
+ for (int t = 0; t < ARRAY_COUNT(kTolerances); t++) {
+ skcms_TransferFunction tf,
+ tf_inv;
+
+ // It's problematic to fit curves with non-zero f, so always force it to zero explicitly.
+ tf.f = 0.0f;
+ int L = fit_linear(curve, N, kTolerances[t], &tf.c, &tf.d);
+
+ if (L == N) {
+ // If the entire data set was linear, move the coefficients to the nonlinear portion
+ // with G == 1. This lets use a canonical representation with d == 0.
+ tf.g = 1;
+ tf.a = tf.c;
+ tf.b = tf.f;
+ tf.c = tf.d = tf.e = tf.f = 0;
+ } else if (L == N - 1) {
+ // Degenerate case with only two points in the nonlinear segment. Solve directly.
+ tf.g = 1;
+ tf.a = (eval_curve(curve, static_cast<float>(N-1)*dx) -
+ eval_curve(curve, static_cast<float>(N-2)*dx))
+ / dx;
+ tf.b = eval_curve(curve, static_cast<float>(N-2)*dx)
+ - tf.a * static_cast<float>(N-2)*dx;
+ tf.e = 0;
+ } else {
+ // Start by guessing a gamma-only curve through the midpoint.
+ int mid = (L + N) / 2;
+ float mid_x = static_cast<float>(mid) / static_cast<float>(N - 1);
+ float mid_y = eval_curve(curve, mid_x);
+ tf.g = log2f_(mid_y) / log2f_(mid_x);
+ tf.a = 1;
+ tf.b = 0;
+ tf.e = tf.c*tf.d + tf.f
+ - powf_(tf.a*tf.d + tf.b, tf.g);
+
+
+ if (!skcms_TransferFunction_invert(&tf, &tf_inv) ||
+ !fit_nonlinear(curve, L,N, &tf_inv)) {
+ continue;
+ }
+
+ // We fit tf_inv, so calculate tf to keep in sync.
+ // fit_nonlinear() should guarantee invertibility.
+ if (!skcms_TransferFunction_invert(&tf_inv, &tf)) {
+ assert(false);
+ continue;
+ }
+ }
+
+ // We'd better have a sane, sRGB-ish TF by now.
+ // Other non-Bad TFs would be fine, but we know we've only ever tried to fit sRGBish;
+ // anything else is just some accident of math and the way we pun tf.g as a type flag.
+ // fit_nonlinear() should guarantee this, but the special cases may fail this test.
+ if (skcms_TFType_sRGBish != classify(tf)) {
+ continue;
+ }
+
+ // We find our error by roundtripping the table through tf_inv.
+ //
+ // (The most likely use case for this approximation is to be inverted and
+ // used as the transfer function for a destination color space.)
+ //
+ // We've kept tf and tf_inv in sync above, but we can't guarantee that tf is
+ // invertible, so re-verify that here (and use the new inverse for testing).
+ // fit_nonlinear() should guarantee this, but the special cases that don't use
+ // it may fail this test.
+ if (!skcms_TransferFunction_invert(&tf, &tf_inv)) {
+ continue;
+ }
+
+ float err = skcms_MaxRoundtripError(curve, &tf_inv);
+ if (*max_error > err) {
+ *max_error = err;
+ *approx = tf;
+ }
+ }
+ return isfinitef_(*max_error);
+}
+
+// ~~~~ Impl. of skcms_Transform() ~~~~
+
+typedef enum {
+ Op_load_a8,
+ Op_load_g8,
+ Op_load_8888_palette8,
+ Op_load_4444,
+ Op_load_565,
+ Op_load_888,
+ Op_load_8888,
+ Op_load_1010102,
+ Op_load_101010x_XR,
+ Op_load_161616LE,
+ Op_load_16161616LE,
+ Op_load_161616BE,
+ Op_load_16161616BE,
+ Op_load_hhh,
+ Op_load_hhhh,
+ Op_load_fff,
+ Op_load_ffff,
+
+ Op_swap_rb,
+ Op_clamp,
+ Op_invert,
+ Op_force_opaque,
+ Op_premul,
+ Op_unpremul,
+ Op_matrix_3x3,
+ Op_matrix_3x4,
+
+ Op_lab_to_xyz,
+ Op_xyz_to_lab,
+
+ Op_tf_r,
+ Op_tf_g,
+ Op_tf_b,
+ Op_tf_a,
+
+ Op_pq_r,
+ Op_pq_g,
+ Op_pq_b,
+ Op_pq_a,
+
+ Op_hlg_r,
+ Op_hlg_g,
+ Op_hlg_b,
+ Op_hlg_a,
+
+ Op_hlginv_r,
+ Op_hlginv_g,
+ Op_hlginv_b,
+ Op_hlginv_a,
+
+ Op_table_r,
+ Op_table_g,
+ Op_table_b,
+ Op_table_a,
+
+ Op_clut_A2B,
+ Op_clut_B2A,
+
+ Op_store_a8,
+ Op_store_g8,
+ Op_store_4444,
+ Op_store_565,
+ Op_store_888,
+ Op_store_8888,
+ Op_store_1010102,
+ Op_store_161616LE,
+ Op_store_16161616LE,
+ Op_store_161616BE,
+ Op_store_16161616BE,
+ Op_store_101010x_XR,
+ Op_store_hhh,
+ Op_store_hhhh,
+ Op_store_fff,
+ Op_store_ffff,
+} Op;
+
+#if defined(__clang__)
+ template <int N, typename T> using Vec = T __attribute__((ext_vector_type(N)));
+#elif defined(__GNUC__)
+ // For some reason GCC accepts this nonsense, but not the more straightforward version,
+ // template <int N, typename T> using Vec = T __attribute__((vector_size(N*sizeof(T))));
+ template <int N, typename T>
+ struct VecHelper { typedef T __attribute__((vector_size(N*sizeof(T)))) V; };
+
+ template <int N, typename T> using Vec = typename VecHelper<N,T>::V;
+#endif
+
+// First, instantiate our default exec_ops() implementation using the default compiliation target.
+
+namespace baseline {
+#if defined(SKCMS_PORTABLE) || !(defined(__clang__) || defined(__GNUC__)) \
+ || (defined(__EMSCRIPTEN_major__) && !defined(__wasm_simd128__))
+ #define N 1
+ template <typename T> using V = T;
+ using Color = float;
+#elif defined(__AVX512F__) && defined(__AVX512DQ__)
+ #define N 16
+ template <typename T> using V = Vec<N,T>;
+ using Color = float;
+#elif defined(__AVX__)
+ #define N 8
+ template <typename T> using V = Vec<N,T>;
+ using Color = float;
+#elif defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(SKCMS_OPT_INTO_NEON_FP16)
+ #define N 8
+ template <typename T> using V = Vec<N,T>;
+ using Color = _Float16;
+#else
+ #define N 4
+ template <typename T> using V = Vec<N,T>;
+ using Color = float;
+#endif
+
+ #include "src/Transform_inl.h"
+ #undef N
+}
+
+// Now, instantiate any other versions of run_program() we may want for runtime detection.
+#if !defined(SKCMS_PORTABLE) && \
+ !defined(SKCMS_NO_RUNTIME_CPU_DETECTION) && \
+ (( defined(__clang__) && __clang_major__ >= 5) || \
+ (!defined(__clang__) && defined(__GNUC__))) \
+ && defined(__x86_64__)
+
+ #if !defined(__AVX2__)
+ #if defined(__clang__)
+ #pragma clang attribute push(__attribute__((target("avx2,f16c"))), apply_to=function)
+ #elif defined(__GNUC__)
+ #pragma GCC push_options
+ #pragma GCC target("avx2,f16c")
+ #endif
+
+ namespace hsw {
+ #define USING_AVX
+ #define USING_AVX_F16C
+ #define USING_AVX2
+ #define N 8
+ template <typename T> using V = Vec<N,T>;
+ using Color = float;
+
+ #include "src/Transform_inl.h"
+
+ // src/Transform_inl.h will undefine USING_* for us.
+ #undef N
+ }
+
+ #if defined(__clang__)
+ #pragma clang attribute pop
+ #elif defined(__GNUC__)
+ #pragma GCC pop_options
+ #endif
+
+ #define TEST_FOR_HSW
+ #endif
+
+ #if !defined(__AVX512F__) || !defined(__AVX512DQ__)
+ #if defined(__clang__)
+ #pragma clang attribute push(__attribute__((target("avx512f,avx512dq,avx512cd,avx512bw,avx512vl"))), apply_to=function)
+ #elif defined(__GNUC__)
+ #pragma GCC push_options
+ #pragma GCC target("avx512f,avx512dq,avx512cd,avx512bw,avx512vl")
+ #endif
+
+ namespace skx {
+ #define USING_AVX512F
+ #define N 16
+ template <typename T> using V = Vec<N,T>;
+ using Color = float;
+
+ #include "src/Transform_inl.h"
+
+ // src/Transform_inl.h will undefine USING_* for us.
+ #undef N
+ }
+
+ #if defined(__clang__)
+ #pragma clang attribute pop
+ #elif defined(__GNUC__)
+ #pragma GCC pop_options
+ #endif
+
+ #define TEST_FOR_SKX
+ #endif
+
+ #if defined(TEST_FOR_HSW) || defined(TEST_FOR_SKX)
+ enum class CpuType { None, HSW, SKX };
+ static CpuType cpu_type() {
+ static const CpuType type = []{
+ if (!runtime_cpu_detection) {
+ return CpuType::None;
+ }
+ // See http://www.sandpile.org/x86/cpuid.htm
+
+ // First, a basic cpuid(1) lets us check prerequisites for HSW, SKX.
+ uint32_t eax, ebx, ecx, edx;
+ __asm__ __volatile__("cpuid" : "=a"(eax), "=b"(ebx), "=c"(ecx), "=d"(edx)
+ : "0"(1), "2"(0));
+ if ((edx & (1u<<25)) && // SSE
+ (edx & (1u<<26)) && // SSE2
+ (ecx & (1u<< 0)) && // SSE3
+ (ecx & (1u<< 9)) && // SSSE3
+ (ecx & (1u<<12)) && // FMA (N.B. not used, avoided even)
+ (ecx & (1u<<19)) && // SSE4.1
+ (ecx & (1u<<20)) && // SSE4.2
+ (ecx & (1u<<26)) && // XSAVE
+ (ecx & (1u<<27)) && // OSXSAVE
+ (ecx & (1u<<28)) && // AVX
+ (ecx & (1u<<29))) { // F16C
+
+ // Call cpuid(7) to check for AVX2 and AVX-512 bits.
+ __asm__ __volatile__("cpuid" : "=a"(eax), "=b"(ebx), "=c"(ecx), "=d"(edx)
+ : "0"(7), "2"(0));
+ // eax from xgetbv(0) will tell us whether XMM, YMM, and ZMM state is saved.
+ uint32_t xcr0, dont_need_edx;
+ __asm__ __volatile__("xgetbv" : "=a"(xcr0), "=d"(dont_need_edx) : "c"(0));
+
+ if ((xcr0 & (1u<<1)) && // XMM register state saved?
+ (xcr0 & (1u<<2)) && // YMM register state saved?
+ (ebx & (1u<<5))) { // AVX2
+ // At this point we're at least HSW. Continue checking for SKX.
+ if ((xcr0 & (1u<< 5)) && // Opmasks state saved?
+ (xcr0 & (1u<< 6)) && // First 16 ZMM registers saved?
+ (xcr0 & (1u<< 7)) && // High 16 ZMM registers saved?
+ (ebx & (1u<<16)) && // AVX512F
+ (ebx & (1u<<17)) && // AVX512DQ
+ (ebx & (1u<<28)) && // AVX512CD
+ (ebx & (1u<<30)) && // AVX512BW
+ (ebx & (1u<<31))) { // AVX512VL
+ return CpuType::SKX;
+ }
+ return CpuType::HSW;
+ }
+ }
+ return CpuType::None;
+ }();
+ return type;
+ }
+ #endif
+
+#endif
+
+typedef struct {
+ Op op;
+ const void* arg;
+} OpAndArg;
+
+static OpAndArg select_curve_op(const skcms_Curve* curve, int channel) {
+ static const struct { Op sRGBish, PQish, HLGish, HLGinvish, table; } ops[] = {
+ { Op_tf_r, Op_pq_r, Op_hlg_r, Op_hlginv_r, Op_table_r },
+ { Op_tf_g, Op_pq_g, Op_hlg_g, Op_hlginv_g, Op_table_g },
+ { Op_tf_b, Op_pq_b, Op_hlg_b, Op_hlginv_b, Op_table_b },
+ { Op_tf_a, Op_pq_a, Op_hlg_a, Op_hlginv_a, Op_table_a },
+ };
+ const auto& op = ops[channel];
+
+ if (curve->table_entries == 0) {
+ const OpAndArg noop = { Op_load_a8/*doesn't matter*/, nullptr };
+
+ const skcms_TransferFunction& tf = curve->parametric;
+
+ if (tf.g == 1 && tf.a == 1 &&
+ tf.b == 0 && tf.c == 0 && tf.d == 0 && tf.e == 0 && tf.f == 0) {
+ return noop;
+ }
+
+ switch (classify(tf)) {
+ case skcms_TFType_Invalid: return noop;
+ case skcms_TFType_sRGBish: return OpAndArg{op.sRGBish, &tf};
+ case skcms_TFType_PQish: return OpAndArg{op.PQish, &tf};
+ case skcms_TFType_HLGish: return OpAndArg{op.HLGish, &tf};
+ case skcms_TFType_HLGinvish: return OpAndArg{op.HLGinvish, &tf};
+ }
+ }
+ return OpAndArg{op.table, curve};
+}
+
+static size_t bytes_per_pixel(skcms_PixelFormat fmt) {
+ switch (fmt >> 1) { // ignore rgb/bgr
+ case skcms_PixelFormat_A_8 >> 1: return 1;
+ case skcms_PixelFormat_G_8 >> 1: return 1;
+ case skcms_PixelFormat_RGBA_8888_Palette8 >> 1: return 1;
+ case skcms_PixelFormat_ABGR_4444 >> 1: return 2;
+ case skcms_PixelFormat_RGB_565 >> 1: return 2;
+ case skcms_PixelFormat_RGB_888 >> 1: return 3;
+ case skcms_PixelFormat_RGBA_8888 >> 1: return 4;
+ case skcms_PixelFormat_RGBA_8888_sRGB >> 1: return 4;
+ case skcms_PixelFormat_RGBA_1010102 >> 1: return 4;
+ case skcms_PixelFormat_RGB_101010x_XR >> 1: return 4;
+ case skcms_PixelFormat_RGB_161616LE >> 1: return 6;
+ case skcms_PixelFormat_RGBA_16161616LE >> 1: return 8;
+ case skcms_PixelFormat_RGB_161616BE >> 1: return 6;
+ case skcms_PixelFormat_RGBA_16161616BE >> 1: return 8;
+ case skcms_PixelFormat_RGB_hhh_Norm >> 1: return 6;
+ case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: return 8;
+ case skcms_PixelFormat_RGB_hhh >> 1: return 6;
+ case skcms_PixelFormat_RGBA_hhhh >> 1: return 8;
+ case skcms_PixelFormat_RGB_fff >> 1: return 12;
+ case skcms_PixelFormat_RGBA_ffff >> 1: return 16;
+ }
+ assert(false);
+ return 0;
+}
+
+static bool prep_for_destination(const skcms_ICCProfile* profile,
+ skcms_Matrix3x3* fromXYZD50,
+ skcms_TransferFunction* invR,
+ skcms_TransferFunction* invG,
+ skcms_TransferFunction* invB) {
+ // skcms_Transform() supports B2A destinations...
+ if (profile->has_B2A) { return true; }
+ // ...and destinations with parametric transfer functions and an XYZD50 gamut matrix.
+ return profile->has_trc
+ && profile->has_toXYZD50
+ && profile->trc[0].table_entries == 0
+ && profile->trc[1].table_entries == 0
+ && profile->trc[2].table_entries == 0
+ && skcms_TransferFunction_invert(&profile->trc[0].parametric, invR)
+ && skcms_TransferFunction_invert(&profile->trc[1].parametric, invG)
+ && skcms_TransferFunction_invert(&profile->trc[2].parametric, invB)
+ && skcms_Matrix3x3_invert(&profile->toXYZD50, fromXYZD50);
+}
+
+bool skcms_Transform(const void* src,
+ skcms_PixelFormat srcFmt,
+ skcms_AlphaFormat srcAlpha,
+ const skcms_ICCProfile* srcProfile,
+ void* dst,
+ skcms_PixelFormat dstFmt,
+ skcms_AlphaFormat dstAlpha,
+ const skcms_ICCProfile* dstProfile,
+ size_t npixels) {
+ return skcms_TransformWithPalette(src, srcFmt, srcAlpha, srcProfile,
+ dst, dstFmt, dstAlpha, dstProfile,
+ npixels, nullptr);
+}
+
+bool skcms_TransformWithPalette(const void* src,
+ skcms_PixelFormat srcFmt,
+ skcms_AlphaFormat srcAlpha,
+ const skcms_ICCProfile* srcProfile,
+ void* dst,
+ skcms_PixelFormat dstFmt,
+ skcms_AlphaFormat dstAlpha,
+ const skcms_ICCProfile* dstProfile,
+ size_t nz,
+ const void* palette) {
+ const size_t dst_bpp = bytes_per_pixel(dstFmt),
+ src_bpp = bytes_per_pixel(srcFmt);
+ // Let's just refuse if the request is absurdly big.
+ if (nz * dst_bpp > INT_MAX || nz * src_bpp > INT_MAX) {
+ return false;
+ }
+ int n = (int)nz;
+
+ // Null profiles default to sRGB. Passing null for both is handy when doing format conversion.
+ if (!srcProfile) {
+ srcProfile = skcms_sRGB_profile();
+ }
+ if (!dstProfile) {
+ dstProfile = skcms_sRGB_profile();
+ }
+
+ // We can't transform in place unless the PixelFormats are the same size.
+ if (dst == src && dst_bpp != src_bpp) {
+ return false;
+ }
+ // TODO: more careful alias rejection (like, dst == src + 1)?
+
+ if (needs_palette(srcFmt) && !palette) {
+ return false;
+ }
+
+ Op program [32];
+ const void* arguments[32];
+
+ Op* ops = program;
+ const void** args = arguments;
+
+ // These are always parametric curves of some sort.
+ skcms_Curve dst_curves[3];
+ dst_curves[0].table_entries =
+ dst_curves[1].table_entries =
+ dst_curves[2].table_entries = 0;
+
+ skcms_Matrix3x3 from_xyz;
+
+ switch (srcFmt >> 1) {
+ default: return false;
+ case skcms_PixelFormat_A_8 >> 1: *ops++ = Op_load_a8; break;
+ case skcms_PixelFormat_G_8 >> 1: *ops++ = Op_load_g8; break;
+ case skcms_PixelFormat_ABGR_4444 >> 1: *ops++ = Op_load_4444; break;
+ case skcms_PixelFormat_RGB_565 >> 1: *ops++ = Op_load_565; break;
+ case skcms_PixelFormat_RGB_888 >> 1: *ops++ = Op_load_888; break;
+ case skcms_PixelFormat_RGBA_8888 >> 1: *ops++ = Op_load_8888; break;
+ case skcms_PixelFormat_RGBA_1010102 >> 1: *ops++ = Op_load_1010102; break;
+ case skcms_PixelFormat_RGB_101010x_XR >> 1: *ops++ = Op_load_101010x_XR; break;
+ case skcms_PixelFormat_RGB_161616LE >> 1: *ops++ = Op_load_161616LE; break;
+ case skcms_PixelFormat_RGBA_16161616LE >> 1: *ops++ = Op_load_16161616LE; break;
+ case skcms_PixelFormat_RGB_161616BE >> 1: *ops++ = Op_load_161616BE; break;
+ case skcms_PixelFormat_RGBA_16161616BE >> 1: *ops++ = Op_load_16161616BE; break;
+ case skcms_PixelFormat_RGB_hhh_Norm >> 1: *ops++ = Op_load_hhh; break;
+ case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: *ops++ = Op_load_hhhh; break;
+ case skcms_PixelFormat_RGB_hhh >> 1: *ops++ = Op_load_hhh; break;
+ case skcms_PixelFormat_RGBA_hhhh >> 1: *ops++ = Op_load_hhhh; break;
+ case skcms_PixelFormat_RGB_fff >> 1: *ops++ = Op_load_fff; break;
+ case skcms_PixelFormat_RGBA_ffff >> 1: *ops++ = Op_load_ffff; break;
+
+ case skcms_PixelFormat_RGBA_8888_Palette8 >> 1: *ops++ = Op_load_8888_palette8;
+ *args++ = palette;
+ break;
+ case skcms_PixelFormat_RGBA_8888_sRGB >> 1:
+ *ops++ = Op_load_8888;
+ *ops++ = Op_tf_r; *args++ = skcms_sRGB_TransferFunction();
+ *ops++ = Op_tf_g; *args++ = skcms_sRGB_TransferFunction();
+ *ops++ = Op_tf_b; *args++ = skcms_sRGB_TransferFunction();
+ break;
+ }
+ if (srcFmt == skcms_PixelFormat_RGB_hhh_Norm ||
+ srcFmt == skcms_PixelFormat_RGBA_hhhh_Norm) {
+ *ops++ = Op_clamp;
+ }
+ if (srcFmt & 1) {
+ *ops++ = Op_swap_rb;
+ }
+ skcms_ICCProfile gray_dst_profile;
+ if ((dstFmt >> 1) == (skcms_PixelFormat_G_8 >> 1)) {
+ // When transforming to gray, stop at XYZ (by setting toXYZ to identity), then transform
+ // luminance (Y) by the destination transfer function.
+ gray_dst_profile = *dstProfile;
+ skcms_SetXYZD50(&gray_dst_profile, &skcms_XYZD50_profile()->toXYZD50);
+ dstProfile = &gray_dst_profile;
+ }
+
+ if (srcProfile->data_color_space == skcms_Signature_CMYK) {
+ // Photoshop creates CMYK images as inverse CMYK.
+ // These happen to be the only ones we've _ever_ seen.
+ *ops++ = Op_invert;
+ // With CMYK, ignore the alpha type, to avoid changing K or conflating CMY with K.
+ srcAlpha = skcms_AlphaFormat_Unpremul;
+ }
+
+ if (srcAlpha == skcms_AlphaFormat_Opaque) {
+ *ops++ = Op_force_opaque;
+ } else if (srcAlpha == skcms_AlphaFormat_PremulAsEncoded) {
+ *ops++ = Op_unpremul;
+ }
+
+ if (dstProfile != srcProfile) {
+
+ if (!prep_for_destination(dstProfile,
+ &from_xyz,
+ &dst_curves[0].parametric,
+ &dst_curves[1].parametric,
+ &dst_curves[2].parametric)) {
+ return false;
+ }
+
+ if (srcProfile->has_A2B) {
+ if (srcProfile->A2B.input_channels) {
+ for (int i = 0; i < (int)srcProfile->A2B.input_channels; i++) {
+ OpAndArg oa = select_curve_op(&srcProfile->A2B.input_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ *ops++ = Op_clamp;
+ *ops++ = Op_clut_A2B;
+ *args++ = &srcProfile->A2B;
+ }
+
+ if (srcProfile->A2B.matrix_channels == 3) {
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(&srcProfile->A2B.matrix_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+
+ static const skcms_Matrix3x4 I = {{
+ {1,0,0,0},
+ {0,1,0,0},
+ {0,0,1,0},
+ }};
+ if (0 != memcmp(&I, &srcProfile->A2B.matrix, sizeof(I))) {
+ *ops++ = Op_matrix_3x4;
+ *args++ = &srcProfile->A2B.matrix;
+ }
+ }
+
+ if (srcProfile->A2B.output_channels == 3) {
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(&srcProfile->A2B.output_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ }
+
+ if (srcProfile->pcs == skcms_Signature_Lab) {
+ *ops++ = Op_lab_to_xyz;
+ }
+
+ } else if (srcProfile->has_trc && srcProfile->has_toXYZD50) {
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(&srcProfile->trc[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ } else {
+ return false;
+ }
+
+ // A2B sources are in XYZD50 by now, but TRC sources are still in their original gamut.
+ assert (srcProfile->has_A2B || srcProfile->has_toXYZD50);
+
+ if (dstProfile->has_B2A) {
+ // B2A needs its input in XYZD50, so transform TRC sources now.
+ if (!srcProfile->has_A2B) {
+ *ops++ = Op_matrix_3x3;
+ *args++ = &srcProfile->toXYZD50;
+ }
+
+ if (dstProfile->pcs == skcms_Signature_Lab) {
+ *ops++ = Op_xyz_to_lab;
+ }
+
+ if (dstProfile->B2A.input_channels == 3) {
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(&dstProfile->B2A.input_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ }
+
+ if (dstProfile->B2A.matrix_channels == 3) {
+ static const skcms_Matrix3x4 I = {{
+ {1,0,0,0},
+ {0,1,0,0},
+ {0,0,1,0},
+ }};
+ if (0 != memcmp(&I, &dstProfile->B2A.matrix, sizeof(I))) {
+ *ops++ = Op_matrix_3x4;
+ *args++ = &dstProfile->B2A.matrix;
+ }
+
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(&dstProfile->B2A.matrix_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ }
+
+ if (dstProfile->B2A.output_channels) {
+ *ops++ = Op_clamp;
+ *ops++ = Op_clut_B2A;
+ *args++ = &dstProfile->B2A;
+ for (int i = 0; i < (int)dstProfile->B2A.output_channels; i++) {
+ OpAndArg oa = select_curve_op(&dstProfile->B2A.output_curves[i], i);
+ if (oa.arg) {
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ }
+ } else {
+ // This is a TRC destination.
+ // We'll concat any src->xyz matrix with our xyz->dst matrix into one src->dst matrix.
+ // (A2B sources are already in XYZD50, making that src->xyz matrix I.)
+ static const skcms_Matrix3x3 I = {{
+ { 1.0f, 0.0f, 0.0f },
+ { 0.0f, 1.0f, 0.0f },
+ { 0.0f, 0.0f, 1.0f },
+ }};
+ const skcms_Matrix3x3* to_xyz = srcProfile->has_A2B ? &I : &srcProfile->toXYZD50;
+
+ // There's a chance the source and destination gamuts are identical,
+ // in which case we can skip the gamut transform.
+ if (0 != memcmp(&dstProfile->toXYZD50, to_xyz, sizeof(skcms_Matrix3x3))) {
+ // Concat the entire gamut transform into from_xyz,
+ // now slightly misnamed but it's a handy spot to stash the result.
+ from_xyz = skcms_Matrix3x3_concat(&from_xyz, to_xyz);
+ *ops++ = Op_matrix_3x3;
+ *args++ = &from_xyz;
+ }
+
+ // Encode back to dst RGB using its parametric transfer functions.
+ for (int i = 0; i < 3; i++) {
+ OpAndArg oa = select_curve_op(dst_curves+i, i);
+ if (oa.arg) {
+ assert (oa.op != Op_table_r &&
+ oa.op != Op_table_g &&
+ oa.op != Op_table_b &&
+ oa.op != Op_table_a);
+ *ops++ = oa.op;
+ *args++ = oa.arg;
+ }
+ }
+ }
+ }
+
+ // Clamp here before premul to make sure we're clamping to normalized values _and_ gamut,
+ // not just to values that fit in [0,1].
+ //
+ // E.g. r = 1.1, a = 0.5 would fit fine in fixed point after premul (ra=0.55,a=0.5),
+ // but would be carrying r > 1, which is really unexpected for downstream consumers.
+ if (dstFmt < skcms_PixelFormat_RGB_hhh) {
+ *ops++ = Op_clamp;
+ }
+
+ if (dstProfile->data_color_space == skcms_Signature_CMYK) {
+ // Photoshop creates CMYK images as inverse CMYK.
+ // These happen to be the only ones we've _ever_ seen.
+ *ops++ = Op_invert;
+
+ // CMYK has no alpha channel, so make sure dstAlpha is a no-op.
+ dstAlpha = skcms_AlphaFormat_Unpremul;
+ }
+
+ if (dstAlpha == skcms_AlphaFormat_Opaque) {
+ *ops++ = Op_force_opaque;
+ } else if (dstAlpha == skcms_AlphaFormat_PremulAsEncoded) {
+ *ops++ = Op_premul;
+ }
+ if (dstFmt & 1) {
+ *ops++ = Op_swap_rb;
+ }
+ switch (dstFmt >> 1) {
+ default: return false;
+ case skcms_PixelFormat_A_8 >> 1: *ops++ = Op_store_a8; break;
+ case skcms_PixelFormat_G_8 >> 1: *ops++ = Op_store_g8; break;
+ case skcms_PixelFormat_ABGR_4444 >> 1: *ops++ = Op_store_4444; break;
+ case skcms_PixelFormat_RGB_565 >> 1: *ops++ = Op_store_565; break;
+ case skcms_PixelFormat_RGB_888 >> 1: *ops++ = Op_store_888; break;
+ case skcms_PixelFormat_RGBA_8888 >> 1: *ops++ = Op_store_8888; break;
+ case skcms_PixelFormat_RGBA_1010102 >> 1: *ops++ = Op_store_1010102; break;
+ case skcms_PixelFormat_RGB_161616LE >> 1: *ops++ = Op_store_161616LE; break;
+ case skcms_PixelFormat_RGBA_16161616LE >> 1: *ops++ = Op_store_16161616LE; break;
+ case skcms_PixelFormat_RGB_161616BE >> 1: *ops++ = Op_store_161616BE; break;
+ case skcms_PixelFormat_RGBA_16161616BE >> 1: *ops++ = Op_store_16161616BE; break;
+ case skcms_PixelFormat_RGB_hhh_Norm >> 1: *ops++ = Op_store_hhh; break;
+ case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: *ops++ = Op_store_hhhh; break;
+ case skcms_PixelFormat_RGB_101010x_XR >> 1: *ops++ = Op_store_101010x_XR; break;
+ case skcms_PixelFormat_RGB_hhh >> 1: *ops++ = Op_store_hhh; break;
+ case skcms_PixelFormat_RGBA_hhhh >> 1: *ops++ = Op_store_hhhh; break;
+ case skcms_PixelFormat_RGB_fff >> 1: *ops++ = Op_store_fff; break;
+ case skcms_PixelFormat_RGBA_ffff >> 1: *ops++ = Op_store_ffff; break;
+
+ case skcms_PixelFormat_RGBA_8888_sRGB >> 1:
+ *ops++ = Op_tf_r; *args++ = skcms_sRGB_Inverse_TransferFunction();
+ *ops++ = Op_tf_g; *args++ = skcms_sRGB_Inverse_TransferFunction();
+ *ops++ = Op_tf_b; *args++ = skcms_sRGB_Inverse_TransferFunction();
+ *ops++ = Op_store_8888;
+ break;
+ }
+
+ auto run = baseline::run_program;
+#if defined(TEST_FOR_HSW)
+ switch (cpu_type()) {
+ case CpuType::None: break;
+ case CpuType::HSW: run = hsw::run_program; break;
+ case CpuType::SKX: run = hsw::run_program; break;
+ }
+#endif
+#if defined(TEST_FOR_SKX)
+ switch (cpu_type()) {
+ case CpuType::None: break;
+ case CpuType::HSW: break;
+ case CpuType::SKX: run = skx::run_program; break;
+ }
+#endif
+ run(program, arguments, (const char*)src, (char*)dst, n, src_bpp,dst_bpp);
+ return true;
+}
+
+static void assert_usable_as_destination(const skcms_ICCProfile* profile) {
+#if defined(NDEBUG)
+ (void)profile;
+#else
+ skcms_Matrix3x3 fromXYZD50;
+ skcms_TransferFunction invR, invG, invB;
+ assert(prep_for_destination(profile, &fromXYZD50, &invR, &invG, &invB));
+#endif
+}
+
+bool skcms_MakeUsableAsDestination(skcms_ICCProfile* profile) {
+ if (!profile->has_B2A) {
+ skcms_Matrix3x3 fromXYZD50;
+ if (!profile->has_trc || !profile->has_toXYZD50
+ || !skcms_Matrix3x3_invert(&profile->toXYZD50, &fromXYZD50)) {
+ return false;
+ }
+
+ skcms_TransferFunction tf[3];
+ for (int i = 0; i < 3; i++) {
+ skcms_TransferFunction inv;
+ if (profile->trc[i].table_entries == 0
+ && skcms_TransferFunction_invert(&profile->trc[i].parametric, &inv)) {
+ tf[i] = profile->trc[i].parametric;
+ continue;
+ }
+
+ float max_error;
+ // Parametric curves from skcms_ApproximateCurve() are guaranteed to be invertible.
+ if (!skcms_ApproximateCurve(&profile->trc[i], &tf[i], &max_error)) {
+ return false;
+ }
+ }
+
+ for (int i = 0; i < 3; ++i) {
+ profile->trc[i].table_entries = 0;
+ profile->trc[i].parametric = tf[i];
+ }
+ }
+ assert_usable_as_destination(profile);
+ return true;
+}
+
+bool skcms_MakeUsableAsDestinationWithSingleCurve(skcms_ICCProfile* profile) {
+ // Call skcms_MakeUsableAsDestination() with B2A disabled;
+ // on success that'll return a TRC/XYZ profile with three skcms_TransferFunctions.
+ skcms_ICCProfile result = *profile;
+ result.has_B2A = false;
+ if (!skcms_MakeUsableAsDestination(&result)) {
+ return false;
+ }
+
+ // Of the three, pick the transfer function that best fits the other two.
+ int best_tf = 0;
+ float min_max_error = INFINITY_;
+ for (int i = 0; i < 3; i++) {
+ skcms_TransferFunction inv;
+ if (!skcms_TransferFunction_invert(&result.trc[i].parametric, &inv)) {
+ return false;
+ }
+
+ float err = 0;
+ for (int j = 0; j < 3; ++j) {
+ err = fmaxf_(err, skcms_MaxRoundtripError(&profile->trc[j], &inv));
+ }
+ if (min_max_error > err) {
+ min_max_error = err;
+ best_tf = i;
+ }
+ }
+
+ for (int i = 0; i < 3; i++) {
+ result.trc[i].parametric = result.trc[best_tf].parametric;
+ }
+
+ *profile = result;
+ assert_usable_as_destination(profile);
+ return true;
+}
diff --git a/gfx/skia/skia/modules/skcms/skcms.gni b/gfx/skia/skia/modules/skcms/skcms.gni
new file mode 100644
index 0000000000..aa7daa2cf4
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/skcms.gni
@@ -0,0 +1,20 @@
+# DO NOT EDIT: This is a generated file.
+# See //bazel/exporter_tool/README.md for more information.
+#
+# The source of truth is //modules/skcms/BUILD.bazel
+
+# To update this file, run make -C bazel generate_gni
+
+_modules = get_path_info("../../modules", "abspath")
+
+# Generated by Bazel rule //modules/skcms:public_hdrs
+skcms_public_headers = [ "$_modules/skcms/skcms.h" ]
+
+# List generated by Bazel rules:
+# //modules/skcms:srcs
+# //modules/skcms:textual_hdrs
+skcms_sources = [
+ "$_modules/skcms/skcms.cc",
+ "$_modules/skcms/skcms_internal.h",
+ "$_modules/skcms/src/Transform_inl.h",
+]
diff --git a/gfx/skia/skia/modules/skcms/skcms.h b/gfx/skia/skia/modules/skcms/skcms.h
new file mode 100644
index 0000000000..322549b38f
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/skcms.h
@@ -0,0 +1,418 @@
+/*
+ * Copyright 2018 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#pragma once
+
+// skcms.h contains the entire public API for skcms.
+
+#ifndef SKCMS_API
+ #define SKCMS_API
+#endif
+
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// A row-major 3x3 matrix (ie vals[row][col])
+typedef struct skcms_Matrix3x3 {
+ float vals[3][3];
+} skcms_Matrix3x3;
+
+// It is _not_ safe to alias the pointers to invert in-place.
+SKCMS_API bool skcms_Matrix3x3_invert(const skcms_Matrix3x3*, skcms_Matrix3x3*);
+SKCMS_API skcms_Matrix3x3 skcms_Matrix3x3_concat(const skcms_Matrix3x3*, const skcms_Matrix3x3*);
+
+// A row-major 3x4 matrix (ie vals[row][col])
+typedef struct skcms_Matrix3x4 {
+ float vals[3][4];
+} skcms_Matrix3x4;
+
+// A transfer function mapping encoded values to linear values,
+// represented by this 7-parameter piecewise function:
+//
+// linear = sign(encoded) * (c*|encoded| + f) , 0 <= |encoded| < d
+// = sign(encoded) * ((a*|encoded| + b)^g + e), d <= |encoded|
+//
+// (A simple gamma transfer function sets g to gamma and a to 1.)
+typedef struct skcms_TransferFunction {
+ float g, a,b,c,d,e,f;
+} skcms_TransferFunction;
+
+SKCMS_API float skcms_TransferFunction_eval (const skcms_TransferFunction*, float);
+SKCMS_API bool skcms_TransferFunction_invert(const skcms_TransferFunction*,
+ skcms_TransferFunction*);
+
+typedef enum skcms_TFType {
+ skcms_TFType_Invalid,
+ skcms_TFType_sRGBish,
+ skcms_TFType_PQish,
+ skcms_TFType_HLGish,
+ skcms_TFType_HLGinvish,
+} skcms_TFType;
+
+// Identify which kind of transfer function is encoded in an skcms_TransferFunction
+SKCMS_API skcms_TFType skcms_TransferFunction_getType(const skcms_TransferFunction*);
+
+// We can jam a couple alternate transfer function forms into skcms_TransferFunction,
+// including those matching the general forms of the SMPTE ST 2084 PQ function or HLG.
+//
+// PQish:
+// max(A + B|encoded|^C, 0)
+// linear = sign(encoded) * (------------------------) ^ F
+// D + E|encoded|^C
+SKCMS_API bool skcms_TransferFunction_makePQish(skcms_TransferFunction*,
+ float A, float B, float C,
+ float D, float E, float F);
+// HLGish:
+// { K * sign(encoded) * ( (R|encoded|)^G ) when 0 <= |encoded| <= 1/R
+// linear = { K * sign(encoded) * ( e^(a(|encoded|-c)) + b ) when 1/R < |encoded|
+SKCMS_API bool skcms_TransferFunction_makeScaledHLGish(skcms_TransferFunction*,
+ float K, float R, float G,
+ float a, float b, float c);
+
+// Compatibility shim with K=1 for old callers.
+static inline bool skcms_TransferFunction_makeHLGish(skcms_TransferFunction* fn,
+ float R, float G,
+ float a, float b, float c) {
+ return skcms_TransferFunction_makeScaledHLGish(fn, 1.0f, R,G, a,b,c);
+}
+
+// PQ mapping encoded [0,1] to linear [0,1].
+static inline bool skcms_TransferFunction_makePQ(skcms_TransferFunction* tf) {
+ return skcms_TransferFunction_makePQish(tf, -107/128.0f, 1.0f, 32/2523.0f
+ , 2413/128.0f, -2392/128.0f, 8192/1305.0f);
+}
+// HLG mapping encoded [0,1] to linear [0,12].
+static inline bool skcms_TransferFunction_makeHLG(skcms_TransferFunction* tf) {
+ return skcms_TransferFunction_makeHLGish(tf, 2.0f, 2.0f
+ , 1/0.17883277f, 0.28466892f, 0.55991073f);
+}
+
+// Is this an ordinary sRGB-ish transfer function, or one of the HDR forms we support?
+SKCMS_API bool skcms_TransferFunction_isSRGBish(const skcms_TransferFunction*);
+SKCMS_API bool skcms_TransferFunction_isPQish (const skcms_TransferFunction*);
+SKCMS_API bool skcms_TransferFunction_isHLGish (const skcms_TransferFunction*);
+
+// Unified representation of 'curv' or 'para' tag data, or a 1D table from 'mft1' or 'mft2'
+typedef union skcms_Curve {
+ struct {
+ uint32_t alias_of_table_entries;
+ skcms_TransferFunction parametric;
+ };
+ struct {
+ uint32_t table_entries;
+ const uint8_t* table_8;
+ const uint8_t* table_16;
+ };
+} skcms_Curve;
+
+// Complex transforms between device space (A) and profile connection space (B):
+// A2B: device -> [ "A" curves -> CLUT ] -> [ "M" curves -> matrix ] -> "B" curves -> PCS
+// B2A: device <- [ "A" curves <- CLUT ] <- [ "M" curves <- matrix ] <- "B" curves <- PCS
+
+typedef struct skcms_A2B {
+ // Optional: N 1D "A" curves, followed by an N-dimensional CLUT.
+ // If input_channels == 0, these curves and CLUT are skipped,
+ // Otherwise, input_channels must be in [1, 4].
+ uint32_t input_channels;
+ skcms_Curve input_curves[4];
+ uint8_t grid_points[4];
+ const uint8_t* grid_8;
+ const uint8_t* grid_16;
+
+ // Optional: 3 1D "M" curves, followed by a color matrix.
+ // If matrix_channels == 0, these curves and matrix are skipped,
+ // Otherwise, matrix_channels must be 3.
+ uint32_t matrix_channels;
+ skcms_Curve matrix_curves[3];
+ skcms_Matrix3x4 matrix;
+
+ // Required: 3 1D "B" curves. Always present, and output_channels must be 3.
+ uint32_t output_channels;
+ skcms_Curve output_curves[3];
+} skcms_A2B;
+
+typedef struct skcms_B2A {
+ // Required: 3 1D "B" curves. Always present, and input_channels must be 3.
+ uint32_t input_channels;
+ skcms_Curve input_curves[3];
+
+ // Optional: a color matrix, followed by 3 1D "M" curves.
+ // If matrix_channels == 0, this matrix and these curves are skipped,
+ // Otherwise, matrix_channels must be 3.
+ uint32_t matrix_channels;
+ skcms_Matrix3x4 matrix;
+ skcms_Curve matrix_curves[3];
+
+ // Optional: an N-dimensional CLUT, followed by N 1D "A" curves.
+ // If output_channels == 0, this CLUT and these curves are skipped,
+ // Otherwise, output_channels must be in [1, 4].
+ uint32_t output_channels;
+ uint8_t grid_points[4];
+ const uint8_t* grid_8;
+ const uint8_t* grid_16;
+ skcms_Curve output_curves[4];
+} skcms_B2A;
+
+typedef struct skcms_CICP {
+ uint8_t color_primaries;
+ uint8_t transfer_characteristics;
+ uint8_t matrix_coefficients;
+ uint8_t video_full_range_flag;
+} skcms_CICP;
+
+typedef struct skcms_ICCProfile {
+ const uint8_t* buffer;
+
+ uint32_t size;
+ uint32_t data_color_space;
+ uint32_t pcs;
+ uint32_t tag_count;
+
+ // skcms_Parse() will set commonly-used fields for you when possible:
+
+ // If we can parse red, green and blue transfer curves from the profile,
+ // trc will be set to those three curves, and has_trc will be true.
+ bool has_trc;
+ skcms_Curve trc[3];
+
+ // If this profile's gamut can be represented by a 3x3 transform to XYZD50,
+ // skcms_Parse() sets toXYZD50 to that transform and has_toXYZD50 to true.
+ bool has_toXYZD50;
+ skcms_Matrix3x3 toXYZD50;
+
+ // If the profile has a valid A2B0 or A2B1 tag, skcms_Parse() sets A2B to
+ // that data, and has_A2B to true. skcms_ParseWithA2BPriority() does the
+ // same following any user-provided prioritization of A2B0, A2B1, or A2B2.
+ bool has_A2B;
+ skcms_A2B A2B;
+
+ // If the profile has a valid B2A0 or B2A1 tag, skcms_Parse() sets B2A to
+ // that data, and has_B2A to true. skcms_ParseWithA2BPriority() does the
+ // same following any user-provided prioritization of B2A0, B2A1, or B2A2.
+ bool has_B2A;
+ skcms_B2A B2A;
+
+ // If the profile has a valid CICP tag, skcms_Parse() sets CICP to that data,
+ // and has_CICP to true.
+ bool has_CICP;
+ skcms_CICP CICP;
+} skcms_ICCProfile;
+
+// The sRGB color profile is so commonly used that we offer a canonical skcms_ICCProfile for it.
+SKCMS_API const skcms_ICCProfile* skcms_sRGB_profile(void);
+// Ditto for XYZD50, the most common profile connection space.
+SKCMS_API const skcms_ICCProfile* skcms_XYZD50_profile(void);
+
+SKCMS_API const skcms_TransferFunction* skcms_sRGB_TransferFunction(void);
+SKCMS_API const skcms_TransferFunction* skcms_sRGB_Inverse_TransferFunction(void);
+SKCMS_API const skcms_TransferFunction* skcms_Identity_TransferFunction(void);
+
+// Practical equality test for two skcms_ICCProfiles.
+// The implementation is subject to change, but it will always try to answer
+// "can I substitute A for B?" and "can I skip transforming from A to B?".
+SKCMS_API bool skcms_ApproximatelyEqualProfiles(const skcms_ICCProfile* A,
+ const skcms_ICCProfile* B);
+
+// Practical test that answers: Is curve roughly the inverse of inv_tf? Typically used by passing
+// the inverse of a known parametric transfer function (like sRGB), to determine if a particular
+// curve is very close to sRGB.
+SKCMS_API bool skcms_AreApproximateInverses(const skcms_Curve* curve,
+ const skcms_TransferFunction* inv_tf);
+
+// Similar to above, answering the question for all three TRC curves of the given profile. Again,
+// passing skcms_sRGB_InverseTransferFunction as inv_tf will answer the question:
+// "Does this profile have a transfer function that is very close to sRGB?"
+SKCMS_API bool skcms_TRCs_AreApproximateInverse(const skcms_ICCProfile* profile,
+ const skcms_TransferFunction* inv_tf);
+
+// Parse an ICC profile and return true if possible, otherwise return false.
+// Selects an A2B profile (if present) according to priority list (each entry 0-2).
+// The buffer is not copied; it must remain valid as long as the skcms_ICCProfile will be used.
+SKCMS_API bool skcms_ParseWithA2BPriority(const void*, size_t,
+ const int priority[], int priorities,
+ skcms_ICCProfile*);
+
+static inline bool skcms_Parse(const void* buf, size_t len, skcms_ICCProfile* profile) {
+ // For continuity of existing user expectations,
+ // prefer A2B0 (perceptual) over A2B1 (relative colormetric), and ignore A2B2 (saturation).
+ const int priority[] = {0,1};
+ return skcms_ParseWithA2BPriority(buf, len,
+ priority, sizeof(priority)/sizeof(*priority),
+ profile);
+}
+
+SKCMS_API bool skcms_ApproximateCurve(const skcms_Curve* curve,
+ skcms_TransferFunction* approx,
+ float* max_error);
+
+SKCMS_API bool skcms_GetCHAD(const skcms_ICCProfile*, skcms_Matrix3x3*);
+SKCMS_API bool skcms_GetWTPT(const skcms_ICCProfile*, float xyz[3]);
+
+// These are common ICC signature values
+enum {
+ // data_color_space
+ skcms_Signature_CMYK = 0x434D594B,
+ skcms_Signature_Gray = 0x47524159,
+ skcms_Signature_RGB = 0x52474220,
+
+ // pcs
+ skcms_Signature_Lab = 0x4C616220,
+ skcms_Signature_XYZ = 0x58595A20,
+};
+
+typedef enum skcms_PixelFormat {
+ skcms_PixelFormat_A_8,
+ skcms_PixelFormat_A_8_,
+ skcms_PixelFormat_G_8,
+ skcms_PixelFormat_G_8_,
+ skcms_PixelFormat_RGBA_8888_Palette8,
+ skcms_PixelFormat_BGRA_8888_Palette8,
+
+ skcms_PixelFormat_RGB_565,
+ skcms_PixelFormat_BGR_565,
+
+ skcms_PixelFormat_ABGR_4444,
+ skcms_PixelFormat_ARGB_4444,
+
+ skcms_PixelFormat_RGB_888,
+ skcms_PixelFormat_BGR_888,
+ skcms_PixelFormat_RGBA_8888,
+ skcms_PixelFormat_BGRA_8888,
+ skcms_PixelFormat_RGBA_8888_sRGB, // Automatic sRGB encoding / decoding.
+ skcms_PixelFormat_BGRA_8888_sRGB, // (Generally used with linear transfer functions.)
+
+ skcms_PixelFormat_RGBA_1010102,
+ skcms_PixelFormat_BGRA_1010102,
+
+ skcms_PixelFormat_RGB_161616LE, // Little-endian. Pointers must be 16-bit aligned.
+ skcms_PixelFormat_BGR_161616LE,
+ skcms_PixelFormat_RGBA_16161616LE,
+ skcms_PixelFormat_BGRA_16161616LE,
+
+ skcms_PixelFormat_RGB_161616BE, // Big-endian. Pointers must be 16-bit aligned.
+ skcms_PixelFormat_BGR_161616BE,
+ skcms_PixelFormat_RGBA_16161616BE,
+ skcms_PixelFormat_BGRA_16161616BE,
+
+ skcms_PixelFormat_RGB_hhh_Norm, // 1-5-10 half-precision float in [0,1]
+ skcms_PixelFormat_BGR_hhh_Norm, // Pointers must be 16-bit aligned.
+ skcms_PixelFormat_RGBA_hhhh_Norm,
+ skcms_PixelFormat_BGRA_hhhh_Norm,
+
+ skcms_PixelFormat_RGB_hhh, // 1-5-10 half-precision float.
+ skcms_PixelFormat_BGR_hhh, // Pointers must be 16-bit aligned.
+ skcms_PixelFormat_RGBA_hhhh,
+ skcms_PixelFormat_BGRA_hhhh,
+
+ skcms_PixelFormat_RGB_fff, // 1-8-23 single-precision float (the normal kind).
+ skcms_PixelFormat_BGR_fff, // Pointers must be 32-bit aligned.
+ skcms_PixelFormat_RGBA_ffff,
+ skcms_PixelFormat_BGRA_ffff,
+
+ skcms_PixelFormat_RGB_101010x_XR, // Note: This is located here to signal no clamping.
+ skcms_PixelFormat_BGR_101010x_XR, // Compatible with MTLPixelFormatBGR10_XR.
+} skcms_PixelFormat;
+
+// We always store any alpha channel linearly. In the chart below, tf-1() is the inverse
+// transfer function for the given color profile (applying the transfer function linearizes).
+
+// We treat opaque as a strong requirement, not just a performance hint: we will ignore
+// any source alpha and treat it as 1.0, and will make sure that any destination alpha
+// channel is filled with the equivalent of 1.0.
+
+// We used to offer multiple types of premultiplication, but now just one, PremulAsEncoded.
+// This is the premul you're probably used to working with.
+
+typedef enum skcms_AlphaFormat {
+ skcms_AlphaFormat_Opaque, // alpha is always opaque
+ // tf-1(r), tf-1(g), tf-1(b), 1.0
+ skcms_AlphaFormat_Unpremul, // alpha and color are unassociated
+ // tf-1(r), tf-1(g), tf-1(b), a
+ skcms_AlphaFormat_PremulAsEncoded, // premultiplied while encoded
+ // tf-1(r)*a, tf-1(g)*a, tf-1(b)*a, a
+} skcms_AlphaFormat;
+
+// Convert npixels pixels from src format and color profile to dst format and color profile
+// and return true, otherwise return false. It is safe to alias dst == src if dstFmt == srcFmt.
+SKCMS_API bool skcms_Transform(const void* src,
+ skcms_PixelFormat srcFmt,
+ skcms_AlphaFormat srcAlpha,
+ const skcms_ICCProfile* srcProfile,
+ void* dst,
+ skcms_PixelFormat dstFmt,
+ skcms_AlphaFormat dstAlpha,
+ const skcms_ICCProfile* dstProfile,
+ size_t npixels);
+
+// As skcms_Transform(), supporting srcFmts with a palette.
+SKCMS_API bool skcms_TransformWithPalette(const void* src,
+ skcms_PixelFormat srcFmt,
+ skcms_AlphaFormat srcAlpha,
+ const skcms_ICCProfile* srcProfile,
+ void* dst,
+ skcms_PixelFormat dstFmt,
+ skcms_AlphaFormat dstAlpha,
+ const skcms_ICCProfile* dstProfile,
+ size_t npixels,
+ const void* palette);
+
+// If profile can be used as a destination in skcms_Transform, return true. Otherwise, attempt to
+// rewrite it with approximations where reasonable. If successful, return true. If no reasonable
+// approximation exists, leave the profile unchanged and return false.
+SKCMS_API bool skcms_MakeUsableAsDestination(skcms_ICCProfile* profile);
+
+// If profile can be used as a destination with a single parametric transfer function (ie for
+// rasterization), return true. Otherwise, attempt to rewrite it with approximations where
+// reasonable. If successful, return true. If no reasonable approximation exists, leave the
+// profile unchanged and return false.
+SKCMS_API bool skcms_MakeUsableAsDestinationWithSingleCurve(skcms_ICCProfile* profile);
+
+// Returns a matrix to adapt XYZ color from given the whitepoint to D50.
+SKCMS_API bool skcms_AdaptToXYZD50(float wx, float wy,
+ skcms_Matrix3x3* toXYZD50);
+
+// Returns a matrix to convert RGB color into XYZ adapted to D50, given the
+// primaries and whitepoint of the RGB model.
+SKCMS_API bool skcms_PrimariesToXYZD50(float rx, float ry,
+ float gx, float gy,
+ float bx, float by,
+ float wx, float wy,
+ skcms_Matrix3x3* toXYZD50);
+
+// Call before your first call to skcms_Transform() to skip runtime CPU detection.
+SKCMS_API void skcms_DisableRuntimeCPUDetection(void);
+
+// Utilities for programmatically constructing profiles
+static inline void skcms_Init(skcms_ICCProfile* p) {
+ memset(p, 0, sizeof(*p));
+ p->data_color_space = skcms_Signature_RGB;
+ p->pcs = skcms_Signature_XYZ;
+}
+
+static inline void skcms_SetTransferFunction(skcms_ICCProfile* p,
+ const skcms_TransferFunction* tf) {
+ p->has_trc = true;
+ for (int i = 0; i < 3; ++i) {
+ p->trc[i].table_entries = 0;
+ p->trc[i].parametric = *tf;
+ }
+}
+
+static inline void skcms_SetXYZD50(skcms_ICCProfile* p, const skcms_Matrix3x3* m) {
+ p->has_toXYZD50 = true;
+ p->toXYZD50 = *m;
+}
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/gfx/skia/skia/modules/skcms/skcms_internal.h b/gfx/skia/skia/modules/skcms/skcms_internal.h
new file mode 100644
index 0000000000..cc6d578ba0
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/skcms_internal.h
@@ -0,0 +1,56 @@
+/*
+ * Copyright 2018 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#pragma once
+
+// skcms_internal.h contains APIs shared by skcms' internals and its test tools.
+// Please don't use this header from outside the skcms repo.
+
+#include "skcms.h"
+#include <stdbool.h>
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// ~~~~ General Helper Macros ~~~~
+ #define ARRAY_COUNT(arr) (int)(sizeof((arr)) / sizeof(*(arr)))
+
+ typedef struct skcms_ICCTag {
+ uint32_t signature;
+ uint32_t type;
+ uint32_t size;
+ const uint8_t* buf;
+ } skcms_ICCTag;
+
+ void skcms_GetTagByIndex (const skcms_ICCProfile*, uint32_t idx, skcms_ICCTag*);
+ bool skcms_GetTagBySignature(const skcms_ICCProfile*, uint32_t sig, skcms_ICCTag*);
+
+ float skcms_MaxRoundtripError(const skcms_Curve* curve, const skcms_TransferFunction* inv_tf);
+
+ // 252 of a random shuffle of all possible bytes.
+ // 252 is evenly divisible by 3 and 4. Only 192, 10, 241, and 43 are missing.
+ // Used for ICC profile equivalence testing.
+ extern const uint8_t skcms_252_random_bytes[252];
+
+// ~~~~ Portable Math ~~~~
+ static inline float floorf_(float x) {
+ float roundtrip = (float)((int)x);
+ return roundtrip > x ? roundtrip - 1 : roundtrip;
+ }
+ static inline float fabsf_(float x) { return x < 0 ? -x : x; }
+ float powf_(float, float);
+
+// ~~~~ Does this pixel format need a palette pointer to be usable? ~~~~
+ static inline bool needs_palette(skcms_PixelFormat fmt) {
+ return (fmt >> 1) == (skcms_PixelFormat_RGBA_8888_Palette8 >> 1);
+ }
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/gfx/skia/skia/modules/skcms/src/Transform_inl.h b/gfx/skia/skia/modules/skcms/src/Transform_inl.h
new file mode 100644
index 0000000000..350f6a20a6
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/src/Transform_inl.h
@@ -0,0 +1,1628 @@
+/*
+ * Copyright 2018 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+// Intentionally NO #pragma once... included multiple times.
+
+// This file is included from skcms.cc in a namespace with some pre-defines:
+// - N: depth of all vectors, 1,4,8, or 16 (preprocessor define)
+// - V<T>: a template to create a vector of N T's.
+
+using F = V<Color>; // Called F for historic reasons... maybe rename C?
+using I32 = V<int32_t>;
+using U64 = V<uint64_t>;
+using U32 = V<uint32_t>;
+using U16 = V<uint16_t>;
+using U8 = V<uint8_t>;
+
+
+#if defined(__GNUC__) && !defined(__clang__)
+ // Once again, GCC is kind of weird, not allowing vector = scalar directly.
+ static constexpr F F0 = F() + 0.0f,
+ F1 = F() + 1.0f,
+ FInfBits = F() + 0x7f800000; // equals 2139095040, the bit pattern of +Inf
+#else
+ static constexpr F F0 = 0.0f,
+ F1 = 1.0f,
+ FInfBits = 0x7f800000; // equals 2139095040, the bit pattern of +Inf
+#endif
+
+// Instead of checking __AVX__ below, we'll check USING_AVX.
+// This lets skcms.cc set USING_AVX to force us in even if the compiler's not set that way.
+// Same deal for __F16C__ and __AVX2__ ~~~> USING_AVX_F16C, USING_AVX2.
+
+#if !defined(USING_AVX) && N == 8 && defined(__AVX__)
+ #define USING_AVX
+#endif
+#if !defined(USING_AVX_F16C) && defined(USING_AVX) && defined(__F16C__)
+ #define USING AVX_F16C
+#endif
+#if !defined(USING_AVX2) && defined(USING_AVX) && defined(__AVX2__)
+ #define USING_AVX2
+#endif
+#if !defined(USING_AVX512F) && N == 16 && defined(__AVX512F__) && defined(__AVX512DQ__)
+ #define USING_AVX512F
+#endif
+
+// Similar to the AVX+ features, we define USING_NEON and USING_NEON_F16C.
+// This is more for organizational clarity... skcms.cc doesn't force these.
+#if N > 1 && defined(__ARM_NEON)
+ #define USING_NEON
+ #if __ARM_FP & 2
+ #define USING_NEON_F16C
+ #endif
+ #if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(SKCMS_OPT_INTO_NEON_FP16)
+ #define USING_NEON_FP16
+ #endif
+#endif
+
+// These -Wvector-conversion warnings seem to trigger in very bogus situations,
+// like vst3q_f32() expecting a 16x char rather than a 4x float vector. :/
+#if defined(USING_NEON) && defined(__clang__)
+ #pragma clang diagnostic ignored "-Wvector-conversion"
+#endif
+
+// GCC & Clang (but not clang-cl) warn returning U64 on x86 is larger than a register.
+// You'd see warnings like, "using AVX even though AVX is not enabled".
+// We stifle these warnings; our helpers that return U64 are always inlined.
+#if defined(__SSE__) && defined(__GNUC__)
+ #if !defined(__has_warning)
+ #pragma GCC diagnostic ignored "-Wpsabi"
+ #elif __has_warning("-Wpsabi")
+ #pragma GCC diagnostic ignored "-Wpsabi"
+ #endif
+#endif
+
+#if defined(__clang__)
+ #define FALLTHROUGH [[clang::fallthrough]]
+#else
+ #define FALLTHROUGH
+#endif
+
+// We tag most helper functions as SI, to enforce good code generation
+// but also work around what we think is a bug in GCC: when targeting 32-bit
+// x86, GCC tends to pass U16 (4x uint16_t vector) function arguments in the
+// MMX mm0 register, which seems to mess with unrelated code that later uses
+// x87 FP instructions (MMX's mm0 is an alias for x87's st0 register).
+//
+// It helps codegen to call __builtin_memcpy() when we know the byte count at compile time.
+#if defined(__clang__) || defined(__GNUC__)
+ #define SI static inline __attribute__((always_inline))
+#else
+ #define SI static inline
+#endif
+
+template <typename T, typename P>
+SI T load(const P* ptr) {
+ T val;
+ small_memcpy(&val, ptr, sizeof(val));
+ return val;
+}
+template <typename T, typename P>
+SI void store(P* ptr, const T& val) {
+ small_memcpy(ptr, &val, sizeof(val));
+}
+
+// (T)v is a cast when N == 1 and a bit-pun when N>1,
+// so we use cast<T>(v) to actually cast or bit_pun<T>(v) to bit-pun.
+template <typename D, typename S>
+SI D cast(const S& v) {
+#if N == 1
+ return (D)v;
+#elif defined(__clang__)
+ return __builtin_convertvector(v, D);
+#else
+ D d;
+ for (int i = 0; i < N; i++) {
+ d[i] = v[i];
+ }
+ return d;
+#endif
+}
+
+template <typename D, typename S>
+SI D bit_pun(const S& v) {
+ static_assert(sizeof(D) == sizeof(v), "");
+ return load<D>(&v);
+}
+
+// When we convert from float to fixed point, it's very common to want to round,
+// and for some reason compilers generate better code when converting to int32_t.
+// To serve both those ends, we use this function to_fixed() instead of direct cast().
+#if defined(USING_NEON_FP16)
+ // NEON's got a F16 -> U16 instruction, so this should be fine without going via I16.
+ SI U16 to_fixed(F f) { return cast<U16>(f + 0.5f); }
+#else
+ SI U32 to_fixed(F f) { return (U32)cast<I32>(f + 0.5f); }
+#endif
+
+
+// Sometimes we do something crazy on one branch of a conditonal,
+// like divide by zero or convert a huge float to an integer,
+// but then harmlessly select the other side. That trips up N==1
+// sanitizer builds, so we make if_then_else() a macro to avoid
+// evaluating the unused side.
+
+#if N == 1
+ #define if_then_else(cond, t, e) ((cond) ? (t) : (e))
+#else
+ template <typename C, typename T>
+ SI T if_then_else(C cond, T t, T e) {
+ return bit_pun<T>( ( cond & bit_pun<C>(t)) |
+ (~cond & bit_pun<C>(e)) );
+ }
+#endif
+
+
+SI F F_from_Half(U16 half) {
+#if defined(USING_NEON_FP16)
+ return bit_pun<F>(half);
+#elif defined(USING_NEON_F16C)
+ return vcvt_f32_f16((float16x4_t)half);
+#elif defined(USING_AVX512F)
+ return (F)_mm512_cvtph_ps((__m256i)half);
+#elif defined(USING_AVX_F16C)
+ typedef int16_t __attribute__((vector_size(16))) I16;
+ return __builtin_ia32_vcvtph2ps256((I16)half);
+#else
+ U32 wide = cast<U32>(half);
+ // A half is 1-5-10 sign-exponent-mantissa, with 15 exponent bias.
+ U32 s = wide & 0x8000,
+ em = wide ^ s;
+
+ // Constructing the float is easy if the half is not denormalized.
+ F norm = bit_pun<F>( (s<<16) + (em<<13) + ((127-15)<<23) );
+
+ // Simply flush all denorm half floats to zero.
+ return if_then_else(em < 0x0400, F0, norm);
+#endif
+}
+
+#if defined(__clang__)
+ // The -((127-15)<<10) underflows that side of the math when
+ // we pass a denorm half float. It's harmless... we'll take the 0 side anyway.
+ __attribute__((no_sanitize("unsigned-integer-overflow")))
+#endif
+SI U16 Half_from_F(F f) {
+#if defined(USING_NEON_FP16)
+ return bit_pun<U16>(f);
+#elif defined(USING_NEON_F16C)
+ return (U16)vcvt_f16_f32(f);
+#elif defined(USING_AVX512F)
+ return (U16)_mm512_cvtps_ph((__m512 )f, _MM_FROUND_CUR_DIRECTION );
+#elif defined(USING_AVX_F16C)
+ return (U16)__builtin_ia32_vcvtps2ph256(f, 0x04/*_MM_FROUND_CUR_DIRECTION*/);
+#else
+ // A float is 1-8-23 sign-exponent-mantissa, with 127 exponent bias.
+ U32 sem = bit_pun<U32>(f),
+ s = sem & 0x80000000,
+ em = sem ^ s;
+
+ // For simplicity we flush denorm half floats (including all denorm floats) to zero.
+ return cast<U16>(if_then_else(em < 0x38800000, (U32)F0
+ , (s>>16) + (em>>13) - ((127-15)<<10)));
+#endif
+}
+
+// Swap high and low bytes of 16-bit lanes, converting between big-endian and little-endian.
+#if defined(USING_NEON_FP16)
+ SI U16 swap_endian_16(U16 v) {
+ return (U16)vrev16q_u8((uint8x16_t) v);
+ }
+#elif defined(USING_NEON)
+ SI U16 swap_endian_16(U16 v) {
+ return (U16)vrev16_u8((uint8x8_t) v);
+ }
+#endif
+
+SI U64 swap_endian_16x4(const U64& rgba) {
+ return (rgba & 0x00ff00ff00ff00ff) << 8
+ | (rgba & 0xff00ff00ff00ff00) >> 8;
+}
+
+#if defined(USING_NEON_FP16)
+ SI F min_(F x, F y) { return (F)vminq_f16((float16x8_t)x, (float16x8_t)y); }
+ SI F max_(F x, F y) { return (F)vmaxq_f16((float16x8_t)x, (float16x8_t)y); }
+#elif defined(USING_NEON)
+ SI F min_(F x, F y) { return (F)vminq_f32((float32x4_t)x, (float32x4_t)y); }
+ SI F max_(F x, F y) { return (F)vmaxq_f32((float32x4_t)x, (float32x4_t)y); }
+#else
+ SI F min_(F x, F y) { return if_then_else(x > y, y, x); }
+ SI F max_(F x, F y) { return if_then_else(x < y, y, x); }
+#endif
+
+SI F floor_(F x) {
+#if N == 1
+ return floorf_(x);
+#elif defined(USING_NEON_FP16)
+ return vrndmq_f16(x);
+#elif defined(__aarch64__)
+ return vrndmq_f32(x);
+#elif defined(USING_AVX512F)
+ // Clang's _mm512_floor_ps() passes its mask as -1, not (__mmask16)-1,
+ // and integer santizer catches that this implicit cast changes the
+ // value from -1 to 65535. We'll cast manually to work around it.
+ // Read this as `return _mm512_floor_ps(x)`.
+ return _mm512_mask_floor_ps(x, (__mmask16)-1, x);
+#elif defined(USING_AVX)
+ return __builtin_ia32_roundps256(x, 0x01/*_MM_FROUND_FLOOR*/);
+#elif defined(__SSE4_1__)
+ return _mm_floor_ps(x);
+#else
+ // Round trip through integers with a truncating cast.
+ F roundtrip = cast<F>(cast<I32>(x));
+ // If x is negative, truncating gives the ceiling instead of the floor.
+ return roundtrip - if_then_else(roundtrip > x, F1, F0);
+
+ // This implementation fails for values of x that are outside
+ // the range an integer can represent. We expect most x to be small.
+#endif
+}
+
+SI F approx_log2(F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ return x;
+#else
+ // The first approximation of log2(x) is its exponent 'e', minus 127.
+ I32 bits = bit_pun<I32>(x);
+
+ F e = cast<F>(bits) * (1.0f / (1<<23));
+
+ // If we use the mantissa too we can refine the error signficantly.
+ F m = bit_pun<F>( (bits & 0x007fffff) | 0x3f000000 );
+
+ return e - 124.225514990f
+ - 1.498030302f*m
+ - 1.725879990f/(0.3520887068f + m);
+#endif
+}
+
+SI F approx_log(F x) {
+ const float ln2 = 0.69314718f;
+ return ln2 * approx_log2(x);
+}
+
+SI F approx_exp2(F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ return x;
+#else
+ F fract = x - floor_(x);
+
+ F fbits = (1.0f * (1<<23)) * (x + 121.274057500f
+ - 1.490129070f*fract
+ + 27.728023300f/(4.84252568f - fract));
+ I32 bits = cast<I32>(min_(max_(fbits, F0), FInfBits));
+
+ return bit_pun<F>(bits);
+#endif
+}
+
+SI F approx_pow(F x, float y) {
+ return if_then_else((x == F0) | (x == F1), x
+ , approx_exp2(approx_log2(x) * y));
+}
+
+SI F approx_exp(F x) {
+ const float log2_e = 1.4426950408889634074f;
+ return approx_exp2(log2_e * x);
+}
+
+// Return tf(x).
+SI F apply_tf(const skcms_TransferFunction* tf, F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ (void)tf;
+ return x;
+#else
+ // Peel off the sign bit and set x = |x|.
+ U32 bits = bit_pun<U32>(x),
+ sign = bits & 0x80000000;
+ x = bit_pun<F>(bits ^ sign);
+
+ // The transfer function has a linear part up to d, exponential at d and after.
+ F v = if_then_else(x < tf->d, tf->c*x + tf->f
+ , approx_pow(tf->a*x + tf->b, tf->g) + tf->e);
+
+ // Tack the sign bit back on.
+ return bit_pun<F>(sign | bit_pun<U32>(v));
+#endif
+}
+
+SI F apply_pq(const skcms_TransferFunction* tf, F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ (void)tf;
+ return x;
+#else
+ U32 bits = bit_pun<U32>(x),
+ sign = bits & 0x80000000;
+ x = bit_pun<F>(bits ^ sign);
+
+ F v = approx_pow(max_(tf->a + tf->b * approx_pow(x, tf->c), F0)
+ / (tf->d + tf->e * approx_pow(x, tf->c)),
+ tf->f);
+
+ return bit_pun<F>(sign | bit_pun<U32>(v));
+#endif
+}
+
+SI F apply_hlg(const skcms_TransferFunction* tf, F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ (void)tf;
+ return x;
+#else
+ const float R = tf->a, G = tf->b,
+ a = tf->c, b = tf->d, c = tf->e,
+ K = tf->f + 1;
+ U32 bits = bit_pun<U32>(x),
+ sign = bits & 0x80000000;
+ x = bit_pun<F>(bits ^ sign);
+
+ F v = if_then_else(x*R <= 1, approx_pow(x*R, G)
+ , approx_exp((x-c)*a) + b);
+
+ return K*bit_pun<F>(sign | bit_pun<U32>(v));
+#endif
+}
+
+SI F apply_hlginv(const skcms_TransferFunction* tf, F x) {
+#if defined(USING_NEON_FP16)
+ // TODO(mtklein)
+ (void)tf;
+ return x;
+#else
+ const float R = tf->a, G = tf->b,
+ a = tf->c, b = tf->d, c = tf->e,
+ K = tf->f + 1;
+ U32 bits = bit_pun<U32>(x),
+ sign = bits & 0x80000000;
+ x = bit_pun<F>(bits ^ sign);
+ x /= K;
+
+ F v = if_then_else(x <= 1, R * approx_pow(x, G)
+ , a * approx_log(x - b) + c);
+
+ return bit_pun<F>(sign | bit_pun<U32>(v));
+#endif
+}
+
+
+// Strided loads and stores of N values, starting from p.
+template <typename T, typename P>
+SI T load_3(const P* p) {
+#if N == 1
+ return (T)p[0];
+#elif N == 4
+ return T{p[ 0],p[ 3],p[ 6],p[ 9]};
+#elif N == 8
+ return T{p[ 0],p[ 3],p[ 6],p[ 9], p[12],p[15],p[18],p[21]};
+#elif N == 16
+ return T{p[ 0],p[ 3],p[ 6],p[ 9], p[12],p[15],p[18],p[21],
+ p[24],p[27],p[30],p[33], p[36],p[39],p[42],p[45]};
+#endif
+}
+
+template <typename T, typename P>
+SI T load_4(const P* p) {
+#if N == 1
+ return (T)p[0];
+#elif N == 4
+ return T{p[ 0],p[ 4],p[ 8],p[12]};
+#elif N == 8
+ return T{p[ 0],p[ 4],p[ 8],p[12], p[16],p[20],p[24],p[28]};
+#elif N == 16
+ return T{p[ 0],p[ 4],p[ 8],p[12], p[16],p[20],p[24],p[28],
+ p[32],p[36],p[40],p[44], p[48],p[52],p[56],p[60]};
+#endif
+}
+
+template <typename T, typename P>
+SI void store_3(P* p, const T& v) {
+#if N == 1
+ p[0] = v;
+#elif N == 4
+ p[ 0] = v[ 0]; p[ 3] = v[ 1]; p[ 6] = v[ 2]; p[ 9] = v[ 3];
+#elif N == 8
+ p[ 0] = v[ 0]; p[ 3] = v[ 1]; p[ 6] = v[ 2]; p[ 9] = v[ 3];
+ p[12] = v[ 4]; p[15] = v[ 5]; p[18] = v[ 6]; p[21] = v[ 7];
+#elif N == 16
+ p[ 0] = v[ 0]; p[ 3] = v[ 1]; p[ 6] = v[ 2]; p[ 9] = v[ 3];
+ p[12] = v[ 4]; p[15] = v[ 5]; p[18] = v[ 6]; p[21] = v[ 7];
+ p[24] = v[ 8]; p[27] = v[ 9]; p[30] = v[10]; p[33] = v[11];
+ p[36] = v[12]; p[39] = v[13]; p[42] = v[14]; p[45] = v[15];
+#endif
+}
+
+template <typename T, typename P>
+SI void store_4(P* p, const T& v) {
+#if N == 1
+ p[0] = v;
+#elif N == 4
+ p[ 0] = v[ 0]; p[ 4] = v[ 1]; p[ 8] = v[ 2]; p[12] = v[ 3];
+#elif N == 8
+ p[ 0] = v[ 0]; p[ 4] = v[ 1]; p[ 8] = v[ 2]; p[12] = v[ 3];
+ p[16] = v[ 4]; p[20] = v[ 5]; p[24] = v[ 6]; p[28] = v[ 7];
+#elif N == 16
+ p[ 0] = v[ 0]; p[ 4] = v[ 1]; p[ 8] = v[ 2]; p[12] = v[ 3];
+ p[16] = v[ 4]; p[20] = v[ 5]; p[24] = v[ 6]; p[28] = v[ 7];
+ p[32] = v[ 8]; p[36] = v[ 9]; p[40] = v[10]; p[44] = v[11];
+ p[48] = v[12]; p[52] = v[13]; p[56] = v[14]; p[60] = v[15];
+#endif
+}
+
+
+SI U8 gather_8(const uint8_t* p, I32 ix) {
+#if N == 1
+ U8 v = p[ix];
+#elif N == 4
+ U8 v = { p[ix[0]], p[ix[1]], p[ix[2]], p[ix[3]] };
+#elif N == 8
+ U8 v = { p[ix[0]], p[ix[1]], p[ix[2]], p[ix[3]],
+ p[ix[4]], p[ix[5]], p[ix[6]], p[ix[7]] };
+#elif N == 16
+ U8 v = { p[ix[ 0]], p[ix[ 1]], p[ix[ 2]], p[ix[ 3]],
+ p[ix[ 4]], p[ix[ 5]], p[ix[ 6]], p[ix[ 7]],
+ p[ix[ 8]], p[ix[ 9]], p[ix[10]], p[ix[11]],
+ p[ix[12]], p[ix[13]], p[ix[14]], p[ix[15]] };
+#endif
+ return v;
+}
+
+SI U16 gather_16(const uint8_t* p, I32 ix) {
+ // Load the i'th 16-bit value from p.
+ auto load_16 = [p](int i) {
+ return load<uint16_t>(p + 2*i);
+ };
+#if N == 1
+ U16 v = load_16(ix);
+#elif N == 4
+ U16 v = { load_16(ix[0]), load_16(ix[1]), load_16(ix[2]), load_16(ix[3]) };
+#elif N == 8
+ U16 v = { load_16(ix[0]), load_16(ix[1]), load_16(ix[2]), load_16(ix[3]),
+ load_16(ix[4]), load_16(ix[5]), load_16(ix[6]), load_16(ix[7]) };
+#elif N == 16
+ U16 v = { load_16(ix[ 0]), load_16(ix[ 1]), load_16(ix[ 2]), load_16(ix[ 3]),
+ load_16(ix[ 4]), load_16(ix[ 5]), load_16(ix[ 6]), load_16(ix[ 7]),
+ load_16(ix[ 8]), load_16(ix[ 9]), load_16(ix[10]), load_16(ix[11]),
+ load_16(ix[12]), load_16(ix[13]), load_16(ix[14]), load_16(ix[15]) };
+#endif
+ return v;
+}
+
+SI U32 gather_32(const uint8_t* p, I32 ix) {
+ // Load the i'th 32-bit value from p.
+ auto load_32 = [p](int i) {
+ return load<uint32_t>(p + 4*i);
+ };
+#if N == 1
+ U32 v = load_32(ix);
+#elif N == 4
+ U32 v = { load_32(ix[0]), load_32(ix[1]), load_32(ix[2]), load_32(ix[3]) };
+#elif N == 8
+ U32 v = { load_32(ix[0]), load_32(ix[1]), load_32(ix[2]), load_32(ix[3]),
+ load_32(ix[4]), load_32(ix[5]), load_32(ix[6]), load_32(ix[7]) };
+#elif N == 16
+ U32 v = { load_32(ix[ 0]), load_32(ix[ 1]), load_32(ix[ 2]), load_32(ix[ 3]),
+ load_32(ix[ 4]), load_32(ix[ 5]), load_32(ix[ 6]), load_32(ix[ 7]),
+ load_32(ix[ 8]), load_32(ix[ 9]), load_32(ix[10]), load_32(ix[11]),
+ load_32(ix[12]), load_32(ix[13]), load_32(ix[14]), load_32(ix[15]) };
+#endif
+ // TODO: AVX2 and AVX-512 gathers (c.f. gather_24).
+ return v;
+}
+
+SI U32 gather_24(const uint8_t* p, I32 ix) {
+ // First, back up a byte. Any place we're gathering from has a safe junk byte to read
+ // in front of it, either a previous table value, or some tag metadata.
+ p -= 1;
+
+ // Load the i'th 24-bit value from p, and 1 extra byte.
+ auto load_24_32 = [p](int i) {
+ return load<uint32_t>(p + 3*i);
+ };
+
+ // Now load multiples of 4 bytes (a junk byte, then r,g,b).
+#if N == 1
+ U32 v = load_24_32(ix);
+#elif N == 4
+ U32 v = { load_24_32(ix[0]), load_24_32(ix[1]), load_24_32(ix[2]), load_24_32(ix[3]) };
+#elif N == 8 && !defined(USING_AVX2)
+ U32 v = { load_24_32(ix[0]), load_24_32(ix[1]), load_24_32(ix[2]), load_24_32(ix[3]),
+ load_24_32(ix[4]), load_24_32(ix[5]), load_24_32(ix[6]), load_24_32(ix[7]) };
+#elif N == 8
+ (void)load_24_32;
+ // The gather instruction here doesn't need any particular alignment,
+ // but the intrinsic takes a const int*.
+ const int* p4 = bit_pun<const int*>(p);
+ I32 zero = { 0, 0, 0, 0, 0, 0, 0, 0},
+ mask = {-1,-1,-1,-1, -1,-1,-1,-1};
+ #if defined(__clang__)
+ U32 v = (U32)__builtin_ia32_gatherd_d256(zero, p4, 3*ix, mask, 1);
+ #elif defined(__GNUC__)
+ U32 v = (U32)__builtin_ia32_gathersiv8si(zero, p4, 3*ix, mask, 1);
+ #endif
+#elif N == 16
+ (void)load_24_32;
+ // The intrinsic is supposed to take const void* now, but it takes const int*, just like AVX2.
+ // And AVX-512 swapped the order of arguments. :/
+ const int* p4 = bit_pun<const int*>(p);
+ U32 v = (U32)_mm512_i32gather_epi32((__m512i)(3*ix), p4, 1);
+#endif
+
+ // Shift off the junk byte, leaving r,g,b in low 24 bits (and zero in the top 8).
+ return v >> 8;
+}
+
+#if !defined(__arm__)
+ SI void gather_48(const uint8_t* p, I32 ix, U64* v) {
+ // As in gather_24(), with everything doubled.
+ p -= 2;
+
+ // Load the i'th 48-bit value from p, and 2 extra bytes.
+ auto load_48_64 = [p](int i) {
+ return load<uint64_t>(p + 6*i);
+ };
+
+ #if N == 1
+ *v = load_48_64(ix);
+ #elif N == 4
+ *v = U64{
+ load_48_64(ix[0]), load_48_64(ix[1]), load_48_64(ix[2]), load_48_64(ix[3]),
+ };
+ #elif N == 8 && !defined(USING_AVX2)
+ *v = U64{
+ load_48_64(ix[0]), load_48_64(ix[1]), load_48_64(ix[2]), load_48_64(ix[3]),
+ load_48_64(ix[4]), load_48_64(ix[5]), load_48_64(ix[6]), load_48_64(ix[7]),
+ };
+ #elif N == 8
+ (void)load_48_64;
+ typedef int32_t __attribute__((vector_size(16))) Half_I32;
+ typedef long long __attribute__((vector_size(32))) Half_I64;
+
+ // The gather instruction here doesn't need any particular alignment,
+ // but the intrinsic takes a const long long*.
+ const long long int* p8 = bit_pun<const long long int*>(p);
+
+ Half_I64 zero = { 0, 0, 0, 0},
+ mask = {-1,-1,-1,-1};
+
+ ix *= 6;
+ Half_I32 ix_lo = { ix[0], ix[1], ix[2], ix[3] },
+ ix_hi = { ix[4], ix[5], ix[6], ix[7] };
+
+ #if defined(__clang__)
+ Half_I64 lo = (Half_I64)__builtin_ia32_gatherd_q256(zero, p8, ix_lo, mask, 1),
+ hi = (Half_I64)__builtin_ia32_gatherd_q256(zero, p8, ix_hi, mask, 1);
+ #elif defined(__GNUC__)
+ Half_I64 lo = (Half_I64)__builtin_ia32_gathersiv4di(zero, p8, ix_lo, mask, 1),
+ hi = (Half_I64)__builtin_ia32_gathersiv4di(zero, p8, ix_hi, mask, 1);
+ #endif
+ store((char*)v + 0, lo);
+ store((char*)v + 32, hi);
+ #elif N == 16
+ (void)load_48_64;
+ const long long int* p8 = bit_pun<const long long int*>(p);
+ __m512i lo = _mm512_i32gather_epi64(_mm512_extracti32x8_epi32((__m512i)(6*ix), 0), p8, 1),
+ hi = _mm512_i32gather_epi64(_mm512_extracti32x8_epi32((__m512i)(6*ix), 1), p8, 1);
+ store((char*)v + 0, lo);
+ store((char*)v + 64, hi);
+ #endif
+
+ *v >>= 16;
+ }
+#endif
+
+SI F F_from_U8(U8 v) {
+ return cast<F>(v) * (1/255.0f);
+}
+
+SI F F_from_U16_BE(U16 v) {
+ // All 16-bit ICC values are big-endian, so we byte swap before converting to float.
+ // MSVC catches the "loss" of data here in the portable path, so we also make sure to mask.
+ U16 lo = (v >> 8),
+ hi = (v << 8) & 0xffff;
+ return cast<F>(lo|hi) * (1/65535.0f);
+}
+
+SI U16 U16_from_F(F v) {
+ // 65535 == inf in FP16, so promote to FP32 before converting.
+ return cast<U16>(cast<V<float>>(v) * 65535 + 0.5f);
+}
+
+SI F minus_1_ulp(F v) {
+#if defined(USING_NEON_FP16)
+ return bit_pun<F>( bit_pun<U16>(v) - 1 );
+#else
+ return bit_pun<F>( bit_pun<U32>(v) - 1 );
+#endif
+}
+
+SI F table(const skcms_Curve* curve, F v) {
+ // Clamp the input to [0,1], then scale to a table index.
+ F ix = max_(F0, min_(v, F1)) * (float)(curve->table_entries - 1);
+
+ // We'll look up (equal or adjacent) entries at lo and hi, then lerp by t between the two.
+ I32 lo = cast<I32>( ix ),
+ hi = cast<I32>(minus_1_ulp(ix+1.0f));
+ F t = ix - cast<F>(lo); // i.e. the fractional part of ix.
+
+ // TODO: can we load l and h simultaneously? Each entry in 'h' is either
+ // the same as in 'l' or adjacent. We have a rough idea that's it'd always be safe
+ // to read adjacent entries and perhaps underflow the table by a byte or two
+ // (it'd be junk, but always safe to read). Not sure how to lerp yet.
+ F l,h;
+ if (curve->table_8) {
+ l = F_from_U8(gather_8(curve->table_8, lo));
+ h = F_from_U8(gather_8(curve->table_8, hi));
+ } else {
+ l = F_from_U16_BE(gather_16(curve->table_16, lo));
+ h = F_from_U16_BE(gather_16(curve->table_16, hi));
+ }
+ return l + (h-l)*t;
+}
+
+SI void sample_clut_8(const uint8_t* grid_8, I32 ix, F* r, F* g, F* b) {
+ U32 rgb = gather_24(grid_8, ix);
+
+ *r = cast<F>((rgb >> 0) & 0xff) * (1/255.0f);
+ *g = cast<F>((rgb >> 8) & 0xff) * (1/255.0f);
+ *b = cast<F>((rgb >> 16) & 0xff) * (1/255.0f);
+}
+
+SI void sample_clut_8(const uint8_t* grid_8, I32 ix, F* r, F* g, F* b, F* a) {
+ // TODO: don't forget to optimize gather_32().
+ U32 rgba = gather_32(grid_8, ix);
+
+ *r = cast<F>((rgba >> 0) & 0xff) * (1/255.0f);
+ *g = cast<F>((rgba >> 8) & 0xff) * (1/255.0f);
+ *b = cast<F>((rgba >> 16) & 0xff) * (1/255.0f);
+ *a = cast<F>((rgba >> 24) & 0xff) * (1/255.0f);
+}
+
+SI void sample_clut_16(const uint8_t* grid_16, I32 ix, F* r, F* g, F* b) {
+#if defined(__arm__)
+ // This is up to 2x faster on 32-bit ARM than the #else-case fast path.
+ *r = F_from_U16_BE(gather_16(grid_16, 3*ix+0));
+ *g = F_from_U16_BE(gather_16(grid_16, 3*ix+1));
+ *b = F_from_U16_BE(gather_16(grid_16, 3*ix+2));
+#else
+ // This strategy is much faster for 64-bit builds, and fine for 32-bit x86 too.
+ U64 rgb;
+ gather_48(grid_16, ix, &rgb);
+ rgb = swap_endian_16x4(rgb);
+
+ *r = cast<F>((rgb >> 0) & 0xffff) * (1/65535.0f);
+ *g = cast<F>((rgb >> 16) & 0xffff) * (1/65535.0f);
+ *b = cast<F>((rgb >> 32) & 0xffff) * (1/65535.0f);
+#endif
+}
+
+SI void sample_clut_16(const uint8_t* grid_16, I32 ix, F* r, F* g, F* b, F* a) {
+ // TODO: gather_64()-based fast path?
+ *r = F_from_U16_BE(gather_16(grid_16, 4*ix+0));
+ *g = F_from_U16_BE(gather_16(grid_16, 4*ix+1));
+ *b = F_from_U16_BE(gather_16(grid_16, 4*ix+2));
+ *a = F_from_U16_BE(gather_16(grid_16, 4*ix+3));
+}
+
+static void clut(uint32_t input_channels, uint32_t output_channels,
+ const uint8_t grid_points[4], const uint8_t* grid_8, const uint8_t* grid_16,
+ F* r, F* g, F* b, F* a) {
+
+ const int dim = (int)input_channels;
+ assert (0 < dim && dim <= 4);
+ assert (output_channels == 3 ||
+ output_channels == 4);
+
+ // For each of these arrays, think foo[2*dim], but we use foo[8] since we know dim <= 4.
+ I32 index [8]; // Index contribution by dimension, first low from 0, then high from 4.
+ F weight[8]; // Weight for each contribution, again first low, then high.
+
+ // O(dim) work first: calculate index,weight from r,g,b,a.
+ const F inputs[] = { *r,*g,*b,*a };
+ for (int i = dim-1, stride = 1; i >= 0; i--) {
+ // x is where we logically want to sample the grid in the i-th dimension.
+ F x = inputs[i] * (float)(grid_points[i] - 1);
+
+ // But we can't index at floats. lo and hi are the two integer grid points surrounding x.
+ I32 lo = cast<I32>( x ), // i.e. trunc(x) == floor(x) here.
+ hi = cast<I32>(minus_1_ulp(x+1.0f));
+ // Notice how we fold in the accumulated stride across previous dimensions here.
+ index[i+0] = lo * stride;
+ index[i+4] = hi * stride;
+ stride *= grid_points[i];
+
+ // We'll interpolate between those two integer grid points by t.
+ F t = x - cast<F>(lo); // i.e. fract(x)
+ weight[i+0] = 1-t;
+ weight[i+4] = t;
+ }
+
+ *r = *g = *b = F0;
+ if (output_channels == 4) {
+ *a = F0;
+ }
+
+ // We'll sample 2^dim == 1<<dim table entries per pixel,
+ // in all combinations of low and high in each dimension.
+ for (int combo = 0; combo < (1<<dim); combo++) { // This loop can be done in any order.
+
+ // Each of these upcoming (combo&N)*K expressions here evaluates to 0 or 4,
+ // where 0 selects the low index contribution and its weight 1-t,
+ // or 4 the high index contribution and its weight t.
+
+ // Since 0<dim≤4, we can always just start off with the 0-th channel,
+ // then handle the others conditionally.
+ I32 ix = index [0 + (combo&1)*4];
+ F w = weight[0 + (combo&1)*4];
+
+ switch ((dim-1)&3) { // This lets the compiler know there are no other cases to handle.
+ case 3: ix += index [3 + (combo&8)/2];
+ w *= weight[3 + (combo&8)/2];
+ FALLTHROUGH;
+ // fall through
+
+ case 2: ix += index [2 + (combo&4)*1];
+ w *= weight[2 + (combo&4)*1];
+ FALLTHROUGH;
+ // fall through
+
+ case 1: ix += index [1 + (combo&2)*2];
+ w *= weight[1 + (combo&2)*2];
+ }
+
+ F R,G,B,A=F0;
+ if (output_channels == 3) {
+ if (grid_8) { sample_clut_8 (grid_8 ,ix, &R,&G,&B); }
+ else { sample_clut_16(grid_16,ix, &R,&G,&B); }
+ } else {
+ if (grid_8) { sample_clut_8 (grid_8 ,ix, &R,&G,&B,&A); }
+ else { sample_clut_16(grid_16,ix, &R,&G,&B,&A); }
+ }
+ *r += w*R;
+ *g += w*G;
+ *b += w*B;
+ *a += w*A;
+ }
+}
+
+static void clut(const skcms_A2B* a2b, F* r, F* g, F* b, F a) {
+ clut(a2b->input_channels, a2b->output_channels,
+ a2b->grid_points, a2b->grid_8, a2b->grid_16,
+ r,g,b,&a);
+}
+static void clut(const skcms_B2A* b2a, F* r, F* g, F* b, F* a) {
+ clut(b2a->input_channels, b2a->output_channels,
+ b2a->grid_points, b2a->grid_8, b2a->grid_16,
+ r,g,b,a);
+}
+
+static void exec_ops(const Op* ops, const void** args,
+ const char* src, char* dst, int i) {
+ F r = F0, g = F0, b = F0, a = F1;
+ while (true) {
+ switch (*ops++) {
+ case Op_load_a8:{
+ a = F_from_U8(load<U8>(src + 1*i));
+ } break;
+
+ case Op_load_g8:{
+ r = g = b = F_from_U8(load<U8>(src + 1*i));
+ } break;
+
+ case Op_load_4444:{
+ U16 abgr = load<U16>(src + 2*i);
+
+ r = cast<F>((abgr >> 12) & 0xf) * (1/15.0f);
+ g = cast<F>((abgr >> 8) & 0xf) * (1/15.0f);
+ b = cast<F>((abgr >> 4) & 0xf) * (1/15.0f);
+ a = cast<F>((abgr >> 0) & 0xf) * (1/15.0f);
+ } break;
+
+ case Op_load_565:{
+ U16 rgb = load<U16>(src + 2*i);
+
+ r = cast<F>(rgb & (uint16_t)(31<< 0)) * (1.0f / (31<< 0));
+ g = cast<F>(rgb & (uint16_t)(63<< 5)) * (1.0f / (63<< 5));
+ b = cast<F>(rgb & (uint16_t)(31<<11)) * (1.0f / (31<<11));
+ } break;
+
+ case Op_load_888:{
+ const uint8_t* rgb = (const uint8_t*)(src + 3*i);
+ #if defined(USING_NEON_FP16)
+ // See the explanation under USING_NEON below. This is that doubled up.
+ uint8x16x3_t v = {{ vdupq_n_u8(0), vdupq_n_u8(0), vdupq_n_u8(0) }};
+ v = vld3q_lane_u8(rgb+ 0, v, 0);
+ v = vld3q_lane_u8(rgb+ 3, v, 2);
+ v = vld3q_lane_u8(rgb+ 6, v, 4);
+ v = vld3q_lane_u8(rgb+ 9, v, 6);
+
+ v = vld3q_lane_u8(rgb+12, v, 8);
+ v = vld3q_lane_u8(rgb+15, v, 10);
+ v = vld3q_lane_u8(rgb+18, v, 12);
+ v = vld3q_lane_u8(rgb+21, v, 14);
+
+ r = cast<F>((U16)v.val[0]) * (1/255.0f);
+ g = cast<F>((U16)v.val[1]) * (1/255.0f);
+ b = cast<F>((U16)v.val[2]) * (1/255.0f);
+ #elif defined(USING_NEON)
+ // There's no uint8x4x3_t or vld3 load for it, so we'll load each rgb pixel one at
+ // a time. Since we're doing that, we might as well load them into 16-bit lanes.
+ // (We'd even load into 32-bit lanes, but that's not possible on ARMv7.)
+ uint8x8x3_t v = {{ vdup_n_u8(0), vdup_n_u8(0), vdup_n_u8(0) }};
+ v = vld3_lane_u8(rgb+0, v, 0);
+ v = vld3_lane_u8(rgb+3, v, 2);
+ v = vld3_lane_u8(rgb+6, v, 4);
+ v = vld3_lane_u8(rgb+9, v, 6);
+
+ // Now if we squint, those 3 uint8x8_t we constructed are really U16s, easy to
+ // convert to F. (Again, U32 would be even better here if drop ARMv7 or split
+ // ARMv7 and ARMv8 impls.)
+ r = cast<F>((U16)v.val[0]) * (1/255.0f);
+ g = cast<F>((U16)v.val[1]) * (1/255.0f);
+ b = cast<F>((U16)v.val[2]) * (1/255.0f);
+ #else
+ r = cast<F>(load_3<U32>(rgb+0) ) * (1/255.0f);
+ g = cast<F>(load_3<U32>(rgb+1) ) * (1/255.0f);
+ b = cast<F>(load_3<U32>(rgb+2) ) * (1/255.0f);
+ #endif
+ } break;
+
+ case Op_load_8888:{
+ U32 rgba = load<U32>(src + 4*i);
+
+ r = cast<F>((rgba >> 0) & 0xff) * (1/255.0f);
+ g = cast<F>((rgba >> 8) & 0xff) * (1/255.0f);
+ b = cast<F>((rgba >> 16) & 0xff) * (1/255.0f);
+ a = cast<F>((rgba >> 24) & 0xff) * (1/255.0f);
+ } break;
+
+ case Op_load_8888_palette8:{
+ const uint8_t* palette = (const uint8_t*) *args++;
+ I32 ix = cast<I32>(load<U8>(src + 1*i));
+ U32 rgba = gather_32(palette, ix);
+
+ r = cast<F>((rgba >> 0) & 0xff) * (1/255.0f);
+ g = cast<F>((rgba >> 8) & 0xff) * (1/255.0f);
+ b = cast<F>((rgba >> 16) & 0xff) * (1/255.0f);
+ a = cast<F>((rgba >> 24) & 0xff) * (1/255.0f);
+ } break;
+
+ case Op_load_1010102:{
+ U32 rgba = load<U32>(src + 4*i);
+
+ r = cast<F>((rgba >> 0) & 0x3ff) * (1/1023.0f);
+ g = cast<F>((rgba >> 10) & 0x3ff) * (1/1023.0f);
+ b = cast<F>((rgba >> 20) & 0x3ff) * (1/1023.0f);
+ a = cast<F>((rgba >> 30) & 0x3 ) * (1/ 3.0f);
+ } break;
+
+ case Op_load_101010x_XR:{
+ static constexpr float min = -0.752941f;
+ static constexpr float max = 1.25098f;
+ static constexpr float range = max - min;
+ U32 rgba = load<U32>(src + 4*i);
+ r = cast<F>((rgba >> 0) & 0x3ff) * (1/1023.0f) * range + min;
+ g = cast<F>((rgba >> 10) & 0x3ff) * (1/1023.0f) * range + min;
+ b = cast<F>((rgba >> 20) & 0x3ff) * (1/1023.0f) * range + min;
+ } break;
+
+ case Op_load_161616LE:{
+ uintptr_t ptr = (uintptr_t)(src + 6*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgb = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = vld3q_u16(rgb);
+ r = cast<F>((U16)v.val[0]) * (1/65535.0f);
+ g = cast<F>((U16)v.val[1]) * (1/65535.0f);
+ b = cast<F>((U16)v.val[2]) * (1/65535.0f);
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = vld3_u16(rgb);
+ r = cast<F>((U16)v.val[0]) * (1/65535.0f);
+ g = cast<F>((U16)v.val[1]) * (1/65535.0f);
+ b = cast<F>((U16)v.val[2]) * (1/65535.0f);
+ #else
+ r = cast<F>(load_3<U32>(rgb+0)) * (1/65535.0f);
+ g = cast<F>(load_3<U32>(rgb+1)) * (1/65535.0f);
+ b = cast<F>(load_3<U32>(rgb+2)) * (1/65535.0f);
+ #endif
+ } break;
+
+ case Op_load_16161616LE:{
+ uintptr_t ptr = (uintptr_t)(src + 8*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgba = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = vld4q_u16(rgba);
+ r = cast<F>((U16)v.val[0]) * (1/65535.0f);
+ g = cast<F>((U16)v.val[1]) * (1/65535.0f);
+ b = cast<F>((U16)v.val[2]) * (1/65535.0f);
+ a = cast<F>((U16)v.val[3]) * (1/65535.0f);
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = vld4_u16(rgba);
+ r = cast<F>((U16)v.val[0]) * (1/65535.0f);
+ g = cast<F>((U16)v.val[1]) * (1/65535.0f);
+ b = cast<F>((U16)v.val[2]) * (1/65535.0f);
+ a = cast<F>((U16)v.val[3]) * (1/65535.0f);
+ #else
+ U64 px = load<U64>(rgba);
+
+ r = cast<F>((px >> 0) & 0xffff) * (1/65535.0f);
+ g = cast<F>((px >> 16) & 0xffff) * (1/65535.0f);
+ b = cast<F>((px >> 32) & 0xffff) * (1/65535.0f);
+ a = cast<F>((px >> 48) & 0xffff) * (1/65535.0f);
+ #endif
+ } break;
+
+ case Op_load_161616BE:{
+ uintptr_t ptr = (uintptr_t)(src + 6*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgb = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = vld3q_u16(rgb);
+ r = cast<F>(swap_endian_16((U16)v.val[0])) * (1/65535.0f);
+ g = cast<F>(swap_endian_16((U16)v.val[1])) * (1/65535.0f);
+ b = cast<F>(swap_endian_16((U16)v.val[2])) * (1/65535.0f);
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = vld3_u16(rgb);
+ r = cast<F>(swap_endian_16((U16)v.val[0])) * (1/65535.0f);
+ g = cast<F>(swap_endian_16((U16)v.val[1])) * (1/65535.0f);
+ b = cast<F>(swap_endian_16((U16)v.val[2])) * (1/65535.0f);
+ #else
+ U32 R = load_3<U32>(rgb+0),
+ G = load_3<U32>(rgb+1),
+ B = load_3<U32>(rgb+2);
+ // R,G,B are big-endian 16-bit, so byte swap them before converting to float.
+ r = cast<F>((R & 0x00ff)<<8 | (R & 0xff00)>>8) * (1/65535.0f);
+ g = cast<F>((G & 0x00ff)<<8 | (G & 0xff00)>>8) * (1/65535.0f);
+ b = cast<F>((B & 0x00ff)<<8 | (B & 0xff00)>>8) * (1/65535.0f);
+ #endif
+ } break;
+
+ case Op_load_16161616BE:{
+ uintptr_t ptr = (uintptr_t)(src + 8*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgba = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = vld4q_u16(rgba);
+ r = cast<F>(swap_endian_16((U16)v.val[0])) * (1/65535.0f);
+ g = cast<F>(swap_endian_16((U16)v.val[1])) * (1/65535.0f);
+ b = cast<F>(swap_endian_16((U16)v.val[2])) * (1/65535.0f);
+ a = cast<F>(swap_endian_16((U16)v.val[3])) * (1/65535.0f);
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = vld4_u16(rgba);
+ r = cast<F>(swap_endian_16((U16)v.val[0])) * (1/65535.0f);
+ g = cast<F>(swap_endian_16((U16)v.val[1])) * (1/65535.0f);
+ b = cast<F>(swap_endian_16((U16)v.val[2])) * (1/65535.0f);
+ a = cast<F>(swap_endian_16((U16)v.val[3])) * (1/65535.0f);
+ #else
+ U64 px = swap_endian_16x4(load<U64>(rgba));
+
+ r = cast<F>((px >> 0) & 0xffff) * (1/65535.0f);
+ g = cast<F>((px >> 16) & 0xffff) * (1/65535.0f);
+ b = cast<F>((px >> 32) & 0xffff) * (1/65535.0f);
+ a = cast<F>((px >> 48) & 0xffff) * (1/65535.0f);
+ #endif
+ } break;
+
+ case Op_load_hhh:{
+ uintptr_t ptr = (uintptr_t)(src + 6*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgb = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = vld3q_u16(rgb);
+ U16 R = (U16)v.val[0],
+ G = (U16)v.val[1],
+ B = (U16)v.val[2];
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = vld3_u16(rgb);
+ U16 R = (U16)v.val[0],
+ G = (U16)v.val[1],
+ B = (U16)v.val[2];
+ #else
+ U16 R = load_3<U16>(rgb+0),
+ G = load_3<U16>(rgb+1),
+ B = load_3<U16>(rgb+2);
+ #endif
+ r = F_from_Half(R);
+ g = F_from_Half(G);
+ b = F_from_Half(B);
+ } break;
+
+ case Op_load_hhhh:{
+ uintptr_t ptr = (uintptr_t)(src + 8*i);
+ assert( (ptr & 1) == 0 ); // src must be 2-byte aligned for this
+ const uint16_t* rgba = (const uint16_t*)ptr; // cast to const uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = vld4q_u16(rgba);
+ U16 R = (U16)v.val[0],
+ G = (U16)v.val[1],
+ B = (U16)v.val[2],
+ A = (U16)v.val[3];
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = vld4_u16(rgba);
+ U16 R = (U16)v.val[0],
+ G = (U16)v.val[1],
+ B = (U16)v.val[2],
+ A = (U16)v.val[3];
+ #else
+ U64 px = load<U64>(rgba);
+ U16 R = cast<U16>((px >> 0) & 0xffff),
+ G = cast<U16>((px >> 16) & 0xffff),
+ B = cast<U16>((px >> 32) & 0xffff),
+ A = cast<U16>((px >> 48) & 0xffff);
+ #endif
+ r = F_from_Half(R);
+ g = F_from_Half(G);
+ b = F_from_Half(B);
+ a = F_from_Half(A);
+ } break;
+
+ case Op_load_fff:{
+ uintptr_t ptr = (uintptr_t)(src + 12*i);
+ assert( (ptr & 3) == 0 ); // src must be 4-byte aligned for this
+ const float* rgb = (const float*)ptr; // cast to const float* to be safe.
+ #if defined(USING_NEON_FP16)
+ float32x4x3_t lo = vld3q_f32(rgb + 0),
+ hi = vld3q_f32(rgb + 12);
+ r = (F)vcombine_f16(vcvt_f16_f32(lo.val[0]), vcvt_f16_f32(hi.val[0]));
+ g = (F)vcombine_f16(vcvt_f16_f32(lo.val[1]), vcvt_f16_f32(hi.val[1]));
+ b = (F)vcombine_f16(vcvt_f16_f32(lo.val[2]), vcvt_f16_f32(hi.val[2]));
+ #elif defined(USING_NEON)
+ float32x4x3_t v = vld3q_f32(rgb);
+ r = (F)v.val[0];
+ g = (F)v.val[1];
+ b = (F)v.val[2];
+ #else
+ r = load_3<F>(rgb+0);
+ g = load_3<F>(rgb+1);
+ b = load_3<F>(rgb+2);
+ #endif
+ } break;
+
+ case Op_load_ffff:{
+ uintptr_t ptr = (uintptr_t)(src + 16*i);
+ assert( (ptr & 3) == 0 ); // src must be 4-byte aligned for this
+ const float* rgba = (const float*)ptr; // cast to const float* to be safe.
+ #if defined(USING_NEON_FP16)
+ float32x4x4_t lo = vld4q_f32(rgba + 0),
+ hi = vld4q_f32(rgba + 16);
+ r = (F)vcombine_f16(vcvt_f16_f32(lo.val[0]), vcvt_f16_f32(hi.val[0]));
+ g = (F)vcombine_f16(vcvt_f16_f32(lo.val[1]), vcvt_f16_f32(hi.val[1]));
+ b = (F)vcombine_f16(vcvt_f16_f32(lo.val[2]), vcvt_f16_f32(hi.val[2]));
+ a = (F)vcombine_f16(vcvt_f16_f32(lo.val[3]), vcvt_f16_f32(hi.val[3]));
+ #elif defined(USING_NEON)
+ float32x4x4_t v = vld4q_f32(rgba);
+ r = (F)v.val[0];
+ g = (F)v.val[1];
+ b = (F)v.val[2];
+ a = (F)v.val[3];
+ #else
+ r = load_4<F>(rgba+0);
+ g = load_4<F>(rgba+1);
+ b = load_4<F>(rgba+2);
+ a = load_4<F>(rgba+3);
+ #endif
+ } break;
+
+ case Op_swap_rb:{
+ F t = r;
+ r = b;
+ b = t;
+ } break;
+
+ case Op_clamp:{
+ r = max_(F0, min_(r, F1));
+ g = max_(F0, min_(g, F1));
+ b = max_(F0, min_(b, F1));
+ a = max_(F0, min_(a, F1));
+ } break;
+
+ case Op_invert:{
+ r = F1 - r;
+ g = F1 - g;
+ b = F1 - b;
+ a = F1 - a;
+ } break;
+
+ case Op_force_opaque:{
+ a = F1;
+ } break;
+
+ case Op_premul:{
+ r *= a;
+ g *= a;
+ b *= a;
+ } break;
+
+ case Op_unpremul:{
+ F scale = if_then_else(F1 / a < INFINITY_, F1 / a, F0);
+ r *= scale;
+ g *= scale;
+ b *= scale;
+ } break;
+
+ case Op_matrix_3x3:{
+ const skcms_Matrix3x3* matrix = (const skcms_Matrix3x3*) *args++;
+ const float* m = &matrix->vals[0][0];
+
+ F R = m[0]*r + m[1]*g + m[2]*b,
+ G = m[3]*r + m[4]*g + m[5]*b,
+ B = m[6]*r + m[7]*g + m[8]*b;
+
+ r = R;
+ g = G;
+ b = B;
+ } break;
+
+ case Op_matrix_3x4:{
+ const skcms_Matrix3x4* matrix = (const skcms_Matrix3x4*) *args++;
+ const float* m = &matrix->vals[0][0];
+
+ F R = m[0]*r + m[1]*g + m[ 2]*b + m[ 3],
+ G = m[4]*r + m[5]*g + m[ 6]*b + m[ 7],
+ B = m[8]*r + m[9]*g + m[10]*b + m[11];
+
+ r = R;
+ g = G;
+ b = B;
+ } break;
+
+ case Op_lab_to_xyz:{
+ // The L*a*b values are in r,g,b, but normalized to [0,1]. Reconstruct them:
+ F L = r * 100.0f,
+ A = g * 255.0f - 128.0f,
+ B = b * 255.0f - 128.0f;
+
+ // Convert to CIE XYZ.
+ F Y = (L + 16.0f) * (1/116.0f),
+ X = Y + A*(1/500.0f),
+ Z = Y - B*(1/200.0f);
+
+ X = if_then_else(X*X*X > 0.008856f, X*X*X, (X - (16/116.0f)) * (1/7.787f));
+ Y = if_then_else(Y*Y*Y > 0.008856f, Y*Y*Y, (Y - (16/116.0f)) * (1/7.787f));
+ Z = if_then_else(Z*Z*Z > 0.008856f, Z*Z*Z, (Z - (16/116.0f)) * (1/7.787f));
+
+ // Adjust to XYZD50 illuminant, and stuff back into r,g,b for the next op.
+ r = X * 0.9642f;
+ g = Y ;
+ b = Z * 0.8249f;
+ } break;
+
+ // As above, in reverse.
+ case Op_xyz_to_lab:{
+ F X = r * (1/0.9642f),
+ Y = g,
+ Z = b * (1/0.8249f);
+
+ X = if_then_else(X > 0.008856f, approx_pow(X, 1/3.0f), X*7.787f + (16/116.0f));
+ Y = if_then_else(Y > 0.008856f, approx_pow(Y, 1/3.0f), Y*7.787f + (16/116.0f));
+ Z = if_then_else(Z > 0.008856f, approx_pow(Z, 1/3.0f), Z*7.787f + (16/116.0f));
+
+ F L = Y*116.0f - 16.0f,
+ A = (X-Y)*500.0f,
+ B = (Y-Z)*200.0f;
+
+ r = L * (1/100.f);
+ g = (A + 128.0f) * (1/255.0f);
+ b = (B + 128.0f) * (1/255.0f);
+ } break;
+
+ case Op_tf_r:{ r = apply_tf((const skcms_TransferFunction*)*args++, r); } break;
+ case Op_tf_g:{ g = apply_tf((const skcms_TransferFunction*)*args++, g); } break;
+ case Op_tf_b:{ b = apply_tf((const skcms_TransferFunction*)*args++, b); } break;
+ case Op_tf_a:{ a = apply_tf((const skcms_TransferFunction*)*args++, a); } break;
+
+ case Op_pq_r:{ r = apply_pq((const skcms_TransferFunction*)*args++, r); } break;
+ case Op_pq_g:{ g = apply_pq((const skcms_TransferFunction*)*args++, g); } break;
+ case Op_pq_b:{ b = apply_pq((const skcms_TransferFunction*)*args++, b); } break;
+ case Op_pq_a:{ a = apply_pq((const skcms_TransferFunction*)*args++, a); } break;
+
+ case Op_hlg_r:{ r = apply_hlg((const skcms_TransferFunction*)*args++, r); } break;
+ case Op_hlg_g:{ g = apply_hlg((const skcms_TransferFunction*)*args++, g); } break;
+ case Op_hlg_b:{ b = apply_hlg((const skcms_TransferFunction*)*args++, b); } break;
+ case Op_hlg_a:{ a = apply_hlg((const skcms_TransferFunction*)*args++, a); } break;
+
+ case Op_hlginv_r:{ r = apply_hlginv((const skcms_TransferFunction*)*args++, r); } break;
+ case Op_hlginv_g:{ g = apply_hlginv((const skcms_TransferFunction*)*args++, g); } break;
+ case Op_hlginv_b:{ b = apply_hlginv((const skcms_TransferFunction*)*args++, b); } break;
+ case Op_hlginv_a:{ a = apply_hlginv((const skcms_TransferFunction*)*args++, a); } break;
+
+ case Op_table_r: { r = table((const skcms_Curve*)*args++, r); } break;
+ case Op_table_g: { g = table((const skcms_Curve*)*args++, g); } break;
+ case Op_table_b: { b = table((const skcms_Curve*)*args++, b); } break;
+ case Op_table_a: { a = table((const skcms_Curve*)*args++, a); } break;
+
+ case Op_clut_A2B: {
+ const skcms_A2B* a2b = (const skcms_A2B*) *args++;
+ clut(a2b, &r,&g,&b,a);
+
+ if (a2b->input_channels == 4) {
+ // CMYK is opaque.
+ a = F1;
+ }
+ } break;
+
+ case Op_clut_B2A: {
+ const skcms_B2A* b2a = (const skcms_B2A*) *args++;
+ clut(b2a, &r,&g,&b,&a);
+ } break;
+
+ // Notice, from here on down the store_ ops all return, ending the loop.
+
+ case Op_store_a8: {
+ store(dst + 1*i, cast<U8>(to_fixed(a * 255)));
+ } return;
+
+ case Op_store_g8: {
+ // g should be holding luminance (Y) (r,g,b ~~~> X,Y,Z)
+ store(dst + 1*i, cast<U8>(to_fixed(g * 255)));
+ } return;
+
+ case Op_store_4444: {
+ store<U16>(dst + 2*i, cast<U16>(to_fixed(r * 15) << 12)
+ | cast<U16>(to_fixed(g * 15) << 8)
+ | cast<U16>(to_fixed(b * 15) << 4)
+ | cast<U16>(to_fixed(a * 15) << 0));
+ } return;
+
+ case Op_store_565: {
+ store<U16>(dst + 2*i, cast<U16>(to_fixed(r * 31) << 0 )
+ | cast<U16>(to_fixed(g * 63) << 5 )
+ | cast<U16>(to_fixed(b * 31) << 11 ));
+ } return;
+
+ case Op_store_888: {
+ uint8_t* rgb = (uint8_t*)dst + 3*i;
+ #if defined(USING_NEON_FP16)
+ // See the explanation under USING_NEON below. This is that doubled up.
+ U16 R = to_fixed(r * 255),
+ G = to_fixed(g * 255),
+ B = to_fixed(b * 255);
+
+ uint8x16x3_t v = {{ (uint8x16_t)R, (uint8x16_t)G, (uint8x16_t)B }};
+ vst3q_lane_u8(rgb+ 0, v, 0);
+ vst3q_lane_u8(rgb+ 3, v, 2);
+ vst3q_lane_u8(rgb+ 6, v, 4);
+ vst3q_lane_u8(rgb+ 9, v, 6);
+
+ vst3q_lane_u8(rgb+12, v, 8);
+ vst3q_lane_u8(rgb+15, v, 10);
+ vst3q_lane_u8(rgb+18, v, 12);
+ vst3q_lane_u8(rgb+21, v, 14);
+ #elif defined(USING_NEON)
+ // Same deal as load_888 but in reverse... we'll store using uint8x8x3_t, but
+ // get there via U16 to save some instructions converting to float. And just
+ // like load_888, we'd prefer to go via U32 but for ARMv7 support.
+ U16 R = cast<U16>(to_fixed(r * 255)),
+ G = cast<U16>(to_fixed(g * 255)),
+ B = cast<U16>(to_fixed(b * 255));
+
+ uint8x8x3_t v = {{ (uint8x8_t)R, (uint8x8_t)G, (uint8x8_t)B }};
+ vst3_lane_u8(rgb+0, v, 0);
+ vst3_lane_u8(rgb+3, v, 2);
+ vst3_lane_u8(rgb+6, v, 4);
+ vst3_lane_u8(rgb+9, v, 6);
+ #else
+ store_3(rgb+0, cast<U8>(to_fixed(r * 255)) );
+ store_3(rgb+1, cast<U8>(to_fixed(g * 255)) );
+ store_3(rgb+2, cast<U8>(to_fixed(b * 255)) );
+ #endif
+ } return;
+
+ case Op_store_8888: {
+ store(dst + 4*i, cast<U32>(to_fixed(r * 255)) << 0
+ | cast<U32>(to_fixed(g * 255)) << 8
+ | cast<U32>(to_fixed(b * 255)) << 16
+ | cast<U32>(to_fixed(a * 255)) << 24);
+ } return;
+
+ case Op_store_101010x_XR: {
+ static constexpr float min = -0.752941f;
+ static constexpr float max = 1.25098f;
+ static constexpr float range = max - min;
+ store(dst + 4*i, cast<U32>(to_fixed(((r - min) / range) * 1023)) << 0
+ | cast<U32>(to_fixed(((g - min) / range) * 1023)) << 10
+ | cast<U32>(to_fixed(((b - min) / range) * 1023)) << 20);
+ return;
+ }
+ case Op_store_1010102: {
+ store(dst + 4*i, cast<U32>(to_fixed(r * 1023)) << 0
+ | cast<U32>(to_fixed(g * 1023)) << 10
+ | cast<U32>(to_fixed(b * 1023)) << 20
+ | cast<U32>(to_fixed(a * 3)) << 30);
+ } return;
+
+ case Op_store_161616LE: {
+ uintptr_t ptr = (uintptr_t)(dst + 6*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgb = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = {{
+ (uint16x8_t)U16_from_F(r),
+ (uint16x8_t)U16_from_F(g),
+ (uint16x8_t)U16_from_F(b),
+ }};
+ vst3q_u16(rgb, v);
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = {{
+ (uint16x4_t)U16_from_F(r),
+ (uint16x4_t)U16_from_F(g),
+ (uint16x4_t)U16_from_F(b),
+ }};
+ vst3_u16(rgb, v);
+ #else
+ store_3(rgb+0, U16_from_F(r));
+ store_3(rgb+1, U16_from_F(g));
+ store_3(rgb+2, U16_from_F(b));
+ #endif
+
+ } return;
+
+ case Op_store_16161616LE: {
+ uintptr_t ptr = (uintptr_t)(dst + 8*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgba = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = {{
+ (uint16x8_t)U16_from_F(r),
+ (uint16x8_t)U16_from_F(g),
+ (uint16x8_t)U16_from_F(b),
+ (uint16x8_t)U16_from_F(a),
+ }};
+ vst4q_u16(rgba, v);
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = {{
+ (uint16x4_t)U16_from_F(r),
+ (uint16x4_t)U16_from_F(g),
+ (uint16x4_t)U16_from_F(b),
+ (uint16x4_t)U16_from_F(a),
+ }};
+ vst4_u16(rgba, v);
+ #else
+ U64 px = cast<U64>(to_fixed(r * 65535)) << 0
+ | cast<U64>(to_fixed(g * 65535)) << 16
+ | cast<U64>(to_fixed(b * 65535)) << 32
+ | cast<U64>(to_fixed(a * 65535)) << 48;
+ store(rgba, px);
+ #endif
+ } return;
+
+ case Op_store_161616BE: {
+ uintptr_t ptr = (uintptr_t)(dst + 6*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgb = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = {{
+ (uint16x8_t)swap_endian_16(U16_from_F(r)),
+ (uint16x8_t)swap_endian_16(U16_from_F(g)),
+ (uint16x8_t)swap_endian_16(U16_from_F(b)),
+ }};
+ vst3q_u16(rgb, v);
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = {{
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(r))),
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(g))),
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(b))),
+ }};
+ vst3_u16(rgb, v);
+ #else
+ U32 R = to_fixed(r * 65535),
+ G = to_fixed(g * 65535),
+ B = to_fixed(b * 65535);
+ store_3(rgb+0, cast<U16>((R & 0x00ff) << 8 | (R & 0xff00) >> 8) );
+ store_3(rgb+1, cast<U16>((G & 0x00ff) << 8 | (G & 0xff00) >> 8) );
+ store_3(rgb+2, cast<U16>((B & 0x00ff) << 8 | (B & 0xff00) >> 8) );
+ #endif
+
+ } return;
+
+ case Op_store_16161616BE: {
+ uintptr_t ptr = (uintptr_t)(dst + 8*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgba = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = {{
+ (uint16x8_t)swap_endian_16(U16_from_F(r)),
+ (uint16x8_t)swap_endian_16(U16_from_F(g)),
+ (uint16x8_t)swap_endian_16(U16_from_F(b)),
+ (uint16x8_t)swap_endian_16(U16_from_F(a)),
+ }};
+ vst4q_u16(rgba, v);
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = {{
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(r))),
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(g))),
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(b))),
+ (uint16x4_t)swap_endian_16(cast<U16>(U16_from_F(a))),
+ }};
+ vst4_u16(rgba, v);
+ #else
+ U64 px = cast<U64>(to_fixed(r * 65535)) << 0
+ | cast<U64>(to_fixed(g * 65535)) << 16
+ | cast<U64>(to_fixed(b * 65535)) << 32
+ | cast<U64>(to_fixed(a * 65535)) << 48;
+ store(rgba, swap_endian_16x4(px));
+ #endif
+ } return;
+
+ case Op_store_hhh: {
+ uintptr_t ptr = (uintptr_t)(dst + 6*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgb = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+
+ U16 R = Half_from_F(r),
+ G = Half_from_F(g),
+ B = Half_from_F(b);
+ #if defined(USING_NEON_FP16)
+ uint16x8x3_t v = {{
+ (uint16x8_t)R,
+ (uint16x8_t)G,
+ (uint16x8_t)B,
+ }};
+ vst3q_u16(rgb, v);
+ #elif defined(USING_NEON)
+ uint16x4x3_t v = {{
+ (uint16x4_t)R,
+ (uint16x4_t)G,
+ (uint16x4_t)B,
+ }};
+ vst3_u16(rgb, v);
+ #else
+ store_3(rgb+0, R);
+ store_3(rgb+1, G);
+ store_3(rgb+2, B);
+ #endif
+ } return;
+
+ case Op_store_hhhh: {
+ uintptr_t ptr = (uintptr_t)(dst + 8*i);
+ assert( (ptr & 1) == 0 ); // The dst pointer must be 2-byte aligned
+ uint16_t* rgba = (uint16_t*)ptr; // for this cast to uint16_t* to be safe.
+
+ U16 R = Half_from_F(r),
+ G = Half_from_F(g),
+ B = Half_from_F(b),
+ A = Half_from_F(a);
+ #if defined(USING_NEON_FP16)
+ uint16x8x4_t v = {{
+ (uint16x8_t)R,
+ (uint16x8_t)G,
+ (uint16x8_t)B,
+ (uint16x8_t)A,
+ }};
+ vst4q_u16(rgba, v);
+ #elif defined(USING_NEON)
+ uint16x4x4_t v = {{
+ (uint16x4_t)R,
+ (uint16x4_t)G,
+ (uint16x4_t)B,
+ (uint16x4_t)A,
+ }};
+ vst4_u16(rgba, v);
+ #else
+ store(rgba, cast<U64>(R) << 0
+ | cast<U64>(G) << 16
+ | cast<U64>(B) << 32
+ | cast<U64>(A) << 48);
+ #endif
+
+ } return;
+
+ case Op_store_fff: {
+ uintptr_t ptr = (uintptr_t)(dst + 12*i);
+ assert( (ptr & 3) == 0 ); // The dst pointer must be 4-byte aligned
+ float* rgb = (float*)ptr; // for this cast to float* to be safe.
+ #if defined(USING_NEON_FP16)
+ float32x4x3_t lo = {{
+ vcvt_f32_f16(vget_low_f16(r)),
+ vcvt_f32_f16(vget_low_f16(g)),
+ vcvt_f32_f16(vget_low_f16(b)),
+ }}, hi = {{
+ vcvt_f32_f16(vget_high_f16(r)),
+ vcvt_f32_f16(vget_high_f16(g)),
+ vcvt_f32_f16(vget_high_f16(b)),
+ }};
+ vst3q_f32(rgb + 0, lo);
+ vst3q_f32(rgb + 12, hi);
+ #elif defined(USING_NEON)
+ float32x4x3_t v = {{
+ (float32x4_t)r,
+ (float32x4_t)g,
+ (float32x4_t)b,
+ }};
+ vst3q_f32(rgb, v);
+ #else
+ store_3(rgb+0, r);
+ store_3(rgb+1, g);
+ store_3(rgb+2, b);
+ #endif
+ } return;
+
+ case Op_store_ffff: {
+ uintptr_t ptr = (uintptr_t)(dst + 16*i);
+ assert( (ptr & 3) == 0 ); // The dst pointer must be 4-byte aligned
+ float* rgba = (float*)ptr; // for this cast to float* to be safe.
+ #if defined(USING_NEON_FP16)
+ float32x4x4_t lo = {{
+ vcvt_f32_f16(vget_low_f16(r)),
+ vcvt_f32_f16(vget_low_f16(g)),
+ vcvt_f32_f16(vget_low_f16(b)),
+ vcvt_f32_f16(vget_low_f16(a)),
+ }}, hi = {{
+ vcvt_f32_f16(vget_high_f16(r)),
+ vcvt_f32_f16(vget_high_f16(g)),
+ vcvt_f32_f16(vget_high_f16(b)),
+ vcvt_f32_f16(vget_high_f16(a)),
+ }};
+ vst4q_f32(rgba + 0, lo);
+ vst4q_f32(rgba + 16, hi);
+ #elif defined(USING_NEON)
+ float32x4x4_t v = {{
+ (float32x4_t)r,
+ (float32x4_t)g,
+ (float32x4_t)b,
+ (float32x4_t)a,
+ }};
+ vst4q_f32(rgba, v);
+ #else
+ store_4(rgba+0, r);
+ store_4(rgba+1, g);
+ store_4(rgba+2, b);
+ store_4(rgba+3, a);
+ #endif
+ } return;
+ }
+ }
+}
+
+
+static void run_program(const Op* program, const void** arguments,
+ const char* src, char* dst, int n,
+ const size_t src_bpp, const size_t dst_bpp) {
+ int i = 0;
+ while (n >= N) {
+ exec_ops(program, arguments, src, dst, i);
+ i += N;
+ n -= N;
+ }
+ if (n > 0) {
+ char tmp[4*4*N] = {0};
+
+ memcpy(tmp, (const char*)src + (size_t)i*src_bpp, (size_t)n*src_bpp);
+ exec_ops(program, arguments, tmp, tmp, 0);
+ memcpy((char*)dst + (size_t)i*dst_bpp, tmp, (size_t)n*dst_bpp);
+ }
+}
+
+// Clean up any #defines we may have set so that we can be #included again.
+#if defined(USING_AVX)
+ #undef USING_AVX
+#endif
+#if defined(USING_AVX_F16C)
+ #undef USING_AVX_F16C
+#endif
+#if defined(USING_AVX2)
+ #undef USING_AVX2
+#endif
+#if defined(USING_AVX512F)
+ #undef USING_AVX512F
+#endif
+
+#if defined(USING_NEON)
+ #undef USING_NEON
+#endif
+#if defined(USING_NEON_F16C)
+ #undef USING_NEON_F16C
+#endif
+#if defined(USING_NEON_FP16)
+ #undef USING_NEON_FP16
+#endif
+
+#undef FALLTHROUGH
diff --git a/gfx/skia/skia/modules/skcms/version.sha1 b/gfx/skia/skia/modules/skcms/version.sha1
new file mode 100755
index 0000000000..ae57df7d06
--- /dev/null
+++ b/gfx/skia/skia/modules/skcms/version.sha1
@@ -0,0 +1 @@
+ba39d81f9797aa973bdf01aa6b0363b280352fba
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-10 21:06:30 +0000