Add color conversion utilities. (#1861)

The color conversion matrix code is copied from apps/shared/iccmaker.c

Add tests.

CMakeLists.txt

@@ -324,6 +324,7 @@ set(AVIF_SRCS
     src/alpha.c
     src/avif.c
     src/colr.c
+    src/colrconvert.c
     src/diag.c
     src/exif.c
     src/io.c

apps/shared/iccmaker.c

@@ -217,7 +217,7 @@ static avifBool xyToXYZ(const float xy[2], double XYZ[3])
         return AVIF_FALSE;
     }
 
-    const double factor = 1 / xy[1];
+    const double factor = 1.0 / xy[1];
     XYZ[0] = xy[0] * factor;
     XYZ[1] = 1;
     XYZ[2] = (1 - xy[0] - xy[1]) * factor;
@@ -381,7 +381,7 @@ avifBool avifGenerateRGBICC(avifRWData * icc, float gamma, const float primaries
         return AVIF_FALSE;
     }
 
-    double rgbPrimaries[3][3] = {
+    const double rgbPrimaries[3][3] = {
         { primaries[0], primaries[2], primaries[4] },
         { primaries[1], primaries[3], primaries[5] },
         { 1.0 - primaries[0] - primaries[1], 1.0 - primaries[2] - primaries[3], 1.0 - primaries[4] - primaries[5] }

include/avif/internal.h

@@ -96,6 +96,19 @@ avifTransferFunction avifTransferCharacteristicsGetLinearToGammaFunction(avifTra
 // Computes the RGB->YUV conversion coefficients kr, kg, kb, such that Y=kr*R+kg*G+kb*B.
 void avifColorPrimariesComputeYCoeffs(avifColorPrimaries colorPrimaries, float coeffs[3]);
 
+// Computes a conversion matrix from RGB to XYZ with a D50 white point.
+AVIF_NODISCARD avifBool avifColorPrimariesComputeRGBToXYZD50Matrix(avifColorPrimaries colorPrimaries, double coeffs[3][3]);
+// Computes a conversion matrix from XYZ with a D50 white point to RGB.
+AVIF_NODISCARD avifBool avifColorPrimariesComputeXYZD50ToRGBMatrix(avifColorPrimaries colorPrimaries, double coeffs[3][3]);
+// Computes the RGB->RGB conversion matrix to convert from one set of RGB primaries to another.
+AVIF_NODISCARD avifBool avifColorPrimariesComputeRGBToRGBMatrix(avifColorPrimaries srcColorPrimaries,
+                                                                avifColorPrimaries dstColorPrimaries,
+                                                                double coeffs[3][3]);
+// Converts the given linear RGB pixel from one color space to another using the provided coefficients.
+// The coefficients can be obtained with avifColorPrimariesComputeRGBToRGBMatrix().
+// The output values are not clamped and may be < 0 or > 1.
+void avifLinearRGBConvertColorSpace(float rgb[4], const double coeffs[3][3]);
+
 #define AVIF_ARRAY_DECLARE(TYPENAME, ITEMSTYPE, ITEMSNAME) \
     typedef struct TYPENAME                                \
     {                                                      \

src/colrconvert.c

@@ -0,0 +1,186 @@
+// Copyright 2023 Google LLC
+// SPDX-License-Identifier: BSD-2-Clause
+
+#include "avif/internal.h"
+
+#include <math.h>
+
+static const double epsilon = 1e-12;
+
+static avifBool avifXyToXYZ(const float xy[2], double XYZ[3])
+{
+    if (fabsf(xy[1]) < epsilon) {
+        return AVIF_FALSE;
+    }
+
+    const double factor = 1.0 / xy[1];
+    XYZ[0] = xy[0] * factor;
+    XYZ[1] = 1;
+    XYZ[2] = (1 - xy[0] - xy[1]) * factor;
+
+    return AVIF_TRUE;
+}
+
+// Computes I = M^-1. Returns false if M seems to be singular.
+static avifBool avifMatInv(const double M[3][3], double I[3][3])
+{
+    double det = M[0][0] * (M[1][1] * M[2][2] - M[2][1] * M[1][2]) - M[0][1] * (M[1][0] * M[2][2] - M[1][2] * M[2][0]) +
+                 M[0][2] * (M[1][0] * M[2][1] - M[1][1] * M[2][0]);
+    if (fabs(det) < epsilon) {
+        return AVIF_FALSE;
+    }
+    det = 1.0 / det;
+
+    I[0][0] = (M[1][1] * M[2][2] - M[2][1] * M[1][2]) * det;
+    I[0][1] = (M[0][2] * M[2][1] - M[0][1] * M[2][2]) * det;
+    I[0][2] = (M[0][1] * M[1][2] - M[0][2] * M[1][1]) * det;
+    I[1][0] = (M[1][2] * M[2][0] - M[1][0] * M[2][2]) * det;
+    I[1][1] = (M[0][0] * M[2][2] - M[0][2] * M[2][0]) * det;
+    I[1][2] = (M[1][0] * M[0][2] - M[0][0] * M[1][2]) * det;
+    I[2][0] = (M[1][0] * M[2][1] - M[2][0] * M[1][1]) * det;
+    I[2][1] = (M[2][0] * M[0][1] - M[0][0] * M[2][1]) * det;
+    I[2][2] = (M[0][0] * M[1][1] - M[1][0] * M[0][1]) * det;
+
+    return AVIF_TRUE;
+}
+
+// Computes C = A*B
+static void avifMatMul(const double A[3][3], const double B[3][3], double C[3][3])
+{
+    C[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0] + A[0][2] * B[2][0];
+    C[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1] + A[0][2] * B[2][1];
+    C[0][2] = A[0][0] * B[0][2] + A[0][1] * B[1][2] + A[0][2] * B[2][2];
+    C[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0] + A[1][2] * B[2][0];
+    C[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1] + A[1][2] * B[2][1];
+    C[1][2] = A[1][0] * B[0][2] + A[1][1] * B[1][2] + A[1][2] * B[2][2];
+    C[2][0] = A[2][0] * B[0][0] + A[2][1] * B[1][0] + A[2][2] * B[2][0];
+    C[2][1] = A[2][0] * B[0][1] + A[2][1] * B[1][1] + A[2][2] * B[2][1];
+    C[2][2] = A[2][0] * B[0][2] + A[2][1] * B[1][2] + A[2][2] * B[2][2];
+}
+
+// Set M to have values of d on the leading diagonal, and zero elsewhere.
+static void avifMatDiag(const double d[3], double M[3][3])
+{
+    M[0][0] = d[0];
+    M[0][1] = 0;
+    M[0][2] = 0;
+    M[1][0] = 0;
+    M[1][1] = d[1];
+    M[1][2] = 0;
+    M[2][0] = 0;
+    M[2][1] = 0;
+    M[2][2] = d[2];
+}
+
+// Computes y = M.x
+static void avifVecMul(const double M[3][3], const double x[3], double y[3])
+{
+    y[0] = M[0][0] * x[0] + M[0][1] * x[1] + M[0][2] * x[2];
+    y[1] = M[1][0] * x[0] + M[1][1] * x[1] + M[1][2] * x[2];
+    y[2] = M[2][0] * x[0] + M[2][1] * x[1] + M[2][2] * x[2];
+}
+
+// Bradford chromatic adaptation matrix
+// from https://www.researchgate.net/publication/253799640_A_uniform_colour_space_based_upon_CIECAM97s
+static const double avifBradford[3][3] = {
+    { 0.8951, 0.2664, -0.1614 },
+    { -0.7502, 1.7135, 0.0367 },
+    { 0.0389, -0.0685, 1.0296 },
+};
+
+// LMS values for D50 whitepoint
+static const double avifLmsD50[3] = { 0.996284, 1.02043, 0.818644 };
+
+avifBool avifColorPrimariesComputeRGBToXYZD50Matrix(avifColorPrimaries colorPrimaries, double coeffs[3][3])
+{
+    float primaries[8];
+    avifColorPrimariesGetValues(colorPrimaries, primaries);
+
+    double whitePointXYZ[3];
+    AVIF_CHECK(avifXyToXYZ(&primaries[6], whitePointXYZ));
+
+    const double rgbPrimaries[3][3] = {
+        { primaries[0], primaries[2], primaries[4] },
+        { primaries[1], primaries[3], primaries[5] },
+        { 1.0 - primaries[0] - primaries[1], 1.0 - primaries[2] - primaries[3], 1.0 - primaries[4] - primaries[5] }
+    };
+
+    double rgbPrimariesInv[3][3];
+    AVIF_CHECK(avifMatInv(rgbPrimaries, rgbPrimariesInv));
+
+    double rgbCoefficients[3];
+    avifVecMul(rgbPrimariesInv, whitePointXYZ, rgbCoefficients);
+
+    double rgbCoefficientsMat[3][3];
+    avifMatDiag(rgbCoefficients, rgbCoefficientsMat);
+
+    double rgbXYZ[3][3];
+    avifMatMul(rgbPrimaries, rgbCoefficientsMat, rgbXYZ);
+
+    // ICC stores primaries XYZ under PCS.
+    // Adapt using linear bradford transform
+    // from https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119021780.app3
+    double lms[3];
+    avifVecMul(avifBradford, whitePointXYZ, lms);
+    for (int i = 0; i < 3; ++i) {
+        if (fabs(lms[i]) < epsilon) {
+            return AVIF_FALSE;
+        }
+        lms[i] = avifLmsD50[i] / lms[i];
+    }
+
+    double adaptation[3][3];
+    avifMatDiag(lms, adaptation);
+
+    double tmp[3][3];
+    avifMatMul(adaptation, avifBradford, tmp);
+
+    double bradfordInv[3][3];
+    if (!avifMatInv(avifBradford, bradfordInv)) {
+        return AVIF_FALSE;
+    }
+    avifMatMul(bradfordInv, tmp, adaptation);
+
+    avifMatMul(adaptation, rgbXYZ, coeffs);
+
+    return AVIF_TRUE;
+}
+
+avifBool avifColorPrimariesComputeXYZD50ToRGBMatrix(avifColorPrimaries colorPrimaries, double coeffs[3][3])
+{
+    double rgbToXyz[3][3];
+    AVIF_CHECK(avifColorPrimariesComputeRGBToXYZD50Matrix(colorPrimaries, rgbToXyz));
+    AVIF_CHECK(avifMatInv(rgbToXyz, coeffs));
+    return AVIF_TRUE;
+}
+
+avifBool avifColorPrimariesComputeRGBToRGBMatrix(avifColorPrimaries srcColorPrimaries,
+                                                 avifColorPrimaries dstColorPrimaries,
+                                                 double coeffs[3][3])
+{
+    // Note: no special casing for srcColorPrimaries == dstColorPrimaries to allow
+    // testing that the computation actually produces the identity matrix.
+    double srcRGBToXYZ[3][3];
+    AVIF_CHECK(avifColorPrimariesComputeRGBToXYZD50Matrix(srcColorPrimaries, srcRGBToXYZ));
+    double xyzToDstRGB[3][3];
+    AVIF_CHECK(avifColorPrimariesComputeXYZD50ToRGBMatrix(dstColorPrimaries, xyzToDstRGB));
+    // coeffs = xyzToDstRGB * srcRGBToXYZ
+    // i.e. srcRGB -> XYZ -> dstRGB
+    avifMatMul(xyzToDstRGB, srcRGBToXYZ, coeffs);
+    return AVIF_TRUE;
+}
+
+// Converts a linear RGBA pixel to a different color space. This function actually works for gamma encoded
+// RGB as well but linear gives better results. Also, for gamma encoded values, it would be
+// better to clamp the output to [0, 1]. Linear values don't need clamping because values
+// > 1.0 are valid for HDR transfer curves, and the gamma compression function will do the
+// clamping as necessary.
+void avifLinearRGBConvertColorSpace(float rgb[4], const double coeffs[3][3])
+{
+    const double rgbDouble[3] = { rgb[0], rgb[1], rgb[2] };
+    double converted[3];
+    avifVecMul(coeffs, rgbDouble, converted);
+    rgb[0] = (float)converted[0];
+    rgb[1] = (float)converted[1];
+    rgb[2] = (float)converted[2];
+}

tests/CMakeLists.txt

@@ -118,6 +118,7 @@ if(AVIF_ENABLE_GTEST)
     add_test(NAME avifincrtest COMMAND avifincrtest ${CMAKE_CURRENT_SOURCE_DIR}/data/)
 
     add_avif_gtest(avifcolrtest)
+    add_avif_gtest_with_data(avifcolrconverttest)
     add_avif_gtest_with_data(avifiostatstest)
     add_avif_gtest_with_data(aviflosslesstest)
     add_avif_gtest_with_data(avifmetadatatest)
@@ -316,6 +317,7 @@ if(AVIF_CODEC_AVM_ENABLED)
                 avifbasictest
                 avifchangesettingtest
                 avifcllitest
+                avifcolrconverttest
                 avifgridapitest
                 avifincrtest
                 avifiostatstest

tests/data/README.md

@@ -470,6 +470,8 @@ Source : created from a personal photo, converted to HDR using Photoshop, see
 https://helpx.adobe.com/camera-raw/using/hdr-output.html and
 https://gregbenzphotography.com/hdr-images/jpg-hdr-gain-maps-in-adobe-camera-raw/
 
+See also sources/seine.psd
+
 HDR image using the PQ transfer curve. Contains a gain map in
 [Adobe's format](https://helpx.adobe.com/camera-raw/using/gain-map.html) that is not recognized by
 libavif and ignored by the tests.
@@ -484,6 +486,14 @@ Source : created from a personal photo, converted to HDR using Photoshop, then s
 see https://helpx.adobe.com/camera-raw/using/hdr-output.html and
 https://gregbenzphotography.com/hdr-images/jpg-hdr-gain-maps-in-adobe-camera-raw/
 
+### File [seine_hdr_rec2020.avif](seine_hdr_rec2020.avif)
+
+![](seine_hdr_rec2020.avif)
+
+License: [same as libavif](https://github.com/AOMediaCodec/libavif/blob/main/LICENSE)
+
+Source : same as seine_hdr_srgb but saved with the Rec. 2020 color space in Photoshop (Camera Raw 16.0.1.1683).
+
 ### File [seine_sdr_gainmap_srgb.avif](seine_sdr_gainmap_srgb.avif)
 
 ![](seine_sdr_gainmap_srgb.avif)
@@ -532,6 +542,25 @@ See `CreateTestImages` in `avifgainmaptest.cc` (set kUpdateTestImages to update
 
 SDR image with a gain map to allow tone mapping to HDR. The gain map's width and height are halved compared to the base image.
 
+## Files colors*_hdr_*.avif and colors*_sdr_srgb.avif
+
+SDR and HDR (PQ) AVIF images in various colorspaces.
+The HDR versions all show the same colors: all colors fit in sRGB but are encoded in various colorspaces.
+
+The colors_text* files have text on them. They are not currently used in tests but can be used for manual
+testing of gain maps, as they make it easy to see which version the browser is displaying.
+
+Source : created with Photoshop 25.1.0 (Camera Raw 16.0.1.1683), see sources/colors.psd and
+https://helpx.adobe.com/camera-raw/using/hdr-output.html,
+https://gregbenzphotography.com/hdr-images/jpg-hdr-gain-maps-in-adobe-camera-raw/
+
+Basic process: create a 32bit image, export it as png for the SDR version.
+Then open the Camera Raw filter (Filter > Camera Raw Filter...), click HDR at the top right, and drag
+the histogram towards the right to create brighter pixels.
+Click the save icon on the top right. Select AVIF as output format and check "HDR output" then save.
+
+To export more images from sources/colors.psd, flatten desired layers before opening the Camera Raw dialog.
+
 ## Animated Images
 
 ### File [colors-animated-8bpc.avif](colors-animated-8bpc.avif)