Add experimental YCgCo-R support

Using the latest version of the spec
https://jvet-experts.org/doc_end_user/documents/29_Teleconference/wg11/JVET-AC1008-v2.zip

Add support for MC 15 and 16 as described in the spec linked above.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f942aee..21a494f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,8 @@
 
 option(AVIF_ENABLE_WERROR "Treat all compiler warnings as errors" ON)
 
+option(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R "Enable experimental YCgCo-R matrix code" OFF)
+
 option(AVIF_CODEC_AOM "Use the AOM codec for encoding/decoding (see AVIF_CODEC_AOM_DECODE/AVIF_CODEC_AOM_ENCODE)" OFF)
 option(AVIF_CODEC_DAV1D "Use the dav1d codec for decoding" OFF)
 option(AVIF_CODEC_LIBGAV1 "Use the libgav1 codec for decoding" OFF)
@@ -205,6 +207,10 @@
     endif()
 endif()
 
+if(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    add_compile_definitions(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+endif()
+
 set(AVIF_SRCS
     src/alpha.c
     src/avif.c
diff --git a/apps/avifenc.c b/apps/avifenc.c
index 0b4ba99..9d06223 100644
--- a/apps/avifenc.c
+++ b/apps/avifenc.c
@@ -1479,11 +1479,23 @@
             returnCode = 1;
         }
         // Matrix coefficients.
-        if (cicpExplicitlySet && settings.matrixCoefficients != AVIF_MATRIX_COEFFICIENTS_IDENTITY) {
-            fprintf(stderr, "Matrix coefficients have to be identity in lossless mode.\n");
-            returnCode = 1;
+        if (cicpExplicitlySet) {
+            avifBool incompatibleMC = (settings.matrixCoefficients != AVIF_MATRIX_COEFFICIENTS_IDENTITY);
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+            incompatibleMC &= (settings.matrixCoefficients != AVIF_MATRIX_COEFFICIENTS_YCGCO_RE &&
+                               settings.matrixCoefficients != AVIF_MATRIX_COEFFICIENTS_YCGCO_RO);
+#endif
+            if (incompatibleMC) {
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+                fprintf(stderr, "Matrix coefficients have to be identity, YCgCo-Re, or YCgCo-Ro in lossless mode.\n");
+#else
+                fprintf(stderr, "Matrix coefficients have to be identity in lossless mode.\n");
+#endif
+                returnCode = 1;
+            }
+        } else {
+            settings.matrixCoefficients = AVIF_MATRIX_COEFFICIENTS_IDENTITY;
         }
-        settings.matrixCoefficients = AVIF_MATRIX_COEFFICIENTS_IDENTITY;
         if (returnCode == 1)
             goto cleanup;
     } else {
@@ -1696,7 +1708,6 @@
     avifBool hasAlpha = (image->alphaPlane && image->alphaRowBytes);
     avifBool usingLosslessColor = (settings.quality == AVIF_QUALITY_LOSSLESS);
     avifBool usingLosslessAlpha = (settings.qualityAlpha == AVIF_QUALITY_LOSSLESS);
-    avifBool depthMatches = (sourceDepth == image->depth);
     avifBool using400 = (image->yuvFormat == AVIF_PIXEL_FORMAT_YUV400);
     avifBool using444 = (image->yuvFormat == AVIF_PIXEL_FORMAT_YUV444);
     avifBool usingFullRange = (image->yuvRange == AVIF_RANGE_FULL);
@@ -1725,9 +1736,9 @@
             lossless = AVIF_FALSE;
         }
 
-        if (!depthMatches) {
+        if (usingIdentityMatrix && (sourceDepth != image->depth)) {
             fprintf(stderr,
-                    "WARNING: [--lossless] Input depth (%d) does not match output depth (%d). Output might not be lossless.\n",
+                    "WARNING: [--lossless] Identity matrix is used but input depth (%d) does not match output depth (%d). Output might not be lossless.\n",
                     sourceDepth,
                     image->depth);
             lossless = AVIF_FALSE;
@@ -1744,8 +1755,17 @@
                 lossless = AVIF_FALSE;
             }
 
-            if (!usingIdentityMatrix && !using400) {
+            avifBool matrixCoefficientsAreLosslessCompatible = usingIdentityMatrix;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+            matrixCoefficientsAreLosslessCompatible |= (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE ||
+                                                        image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO);
+#endif
+            if (!matrixCoefficientsAreLosslessCompatible && !using400) {
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+                fprintf(stderr, "WARNING: [--lossless] Input data was RGB and matrixCoefficients isn't set to identity (--cicp x/x/0) or YCgCo-Re/Ro (--cicp x/x/15 or x/x/16); Output might not be lossless.\n");
+#else
                 fprintf(stderr, "WARNING: [--lossless] Input data was RGB and matrixCoefficients isn't set to identity (--cicp x/x/0); Output might not be lossless.\n");
+#endif
                 lossless = AVIF_FALSE;
             }
         }
diff --git a/apps/shared/avifjpeg.c b/apps/shared/avifjpeg.c
index 87ff123..b599b3e 100644
--- a/apps/shared/avifjpeg.c
+++ b/apps/shared/avifjpeg.c
@@ -356,12 +356,32 @@
 
         avif->width = cinfo.output_width;
         avif->height = cinfo.output_height;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+        const avifBool useYCgCoR = (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE ||
+                                    avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO);
+#else
+        const avifBool useYCgCoR = AVIF_FALSE;
+#endif
         if (avif->yuvFormat == AVIF_PIXEL_FORMAT_NONE) {
-            // Identity is only valid with YUV444.
-            avif->yuvFormat = (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_IDENTITY) ? AVIF_PIXEL_FORMAT_YUV444
-                                                                                              : AVIF_APP_DEFAULT_PIXEL_FORMAT;
+            // Identity and YCgCo-R are only valid with YUV444.
+            avif->yuvFormat = (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_IDENTITY || useYCgCoR)
+                                  ? AVIF_PIXEL_FORMAT_YUV444
+                                  : AVIF_APP_DEFAULT_PIXEL_FORMAT;
         }
         avif->depth = requestedDepth ? requestedDepth : 8;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+        if (useYCgCoR) {
+            if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO) {
+                fprintf(stderr, "AVIF_MATRIX_COEFFICIENTS_YCGCO_RO cannot be used with JPEG because it has an even bit depth.\n");
+                goto cleanup;
+            }
+            if (requestedDepth && requestedDepth != 10) {
+                fprintf(stderr, "Cannot request %u bits for YCgCo-Re as it uses 2 extra bits.\n", requestedDepth);
+                goto cleanup;
+            }
+            avif->depth = 10;
+        }
+#endif
         avifRGBImageSetDefaults(&rgb, avif);
         rgb.format = AVIF_RGB_FORMAT_RGB;
         rgb.chromaDownsampling = chromaDownsampling;
diff --git a/apps/shared/avifpng.c b/apps/shared/avifpng.c
index bcd908a..bb18ebc 100644
--- a/apps/shared/avifpng.c
+++ b/apps/shared/avifpng.c
@@ -313,9 +313,18 @@
     avif->width = rawWidth;
     avif->height = rawHeight;
     avif->yuvFormat = requestedFormat;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO) {
+        fprintf(stderr, "AVIF_MATRIX_COEFFICIENTS_YCGCO_RO cannot be used with PNG because it has an even bit depth.\n");
+        goto cleanup;
+    }
+    const avifBool useYCgCoR = (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE);
+#else
+    const avifBool useYCgCoR = AVIF_FALSE;
+#endif
     if (avif->yuvFormat == AVIF_PIXEL_FORMAT_NONE) {
-        if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_IDENTITY) {
-            // Identity is only valid with YUV444.
+        if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_IDENTITY || useYCgCoR) {
+            // Identity and YCgCo-R are only valid with YUV444.
             avif->yuvFormat = AVIF_PIXEL_FORMAT_YUV444;
         } else if ((rawColorType == PNG_COLOR_TYPE_GRAY) || (rawColorType == PNG_COLOR_TYPE_GRAY_ALPHA)) {
             avif->yuvFormat = AVIF_PIXEL_FORMAT_YUV400;
@@ -331,6 +340,19 @@
             avif->depth = 12;
         }
     }
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    if (useYCgCoR) {
+        if (imgBitDepth != 8) {
+            fprintf(stderr, "AVIF_MATRIX_COEFFICIENTS_YCGCO_RE cannot be used on 16 bit input because it adds two bits.\n");
+            goto cleanup;
+        }
+        if (requestedDepth && requestedDepth != 10) {
+            fprintf(stderr, "Cannot request %u bits for YCgCo-Re as it uses 2 extra bits.\n", requestedDepth);
+            goto cleanup;
+        }
+        avif->depth = 10;
+    }
+#endif
 
     avifRGBImageSetDefaults(&rgb, avif);
     rgb.chromaDownsampling = chromaDownsampling;
@@ -396,12 +418,22 @@
 
     volatile int rgbDepth = requestedDepth;
     if (rgbDepth == 0) {
-        if (avif->depth > 8) {
-            rgbDepth = 16;
-        } else {
-            rgbDepth = 8;
-        }
+        rgbDepth = (avif->depth > 8) ? 16 : 8;
     }
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO) {
+        fprintf(stderr, "AVIF_MATRIX_COEFFICIENTS_YCGCO_RO cannot be used with PNG because it has an even bit depth.\n");
+        goto cleanup;
+    }
+    if (avif->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE) {
+        if (avif->depth != 10 || (requestedDepth && requestedDepth != 8)) {
+            fprintf(stderr, "Cannot request %u bits for YCgCo-Re as it only works for 8 bits.\n", requestedDepth);
+            goto cleanup;
+        }
+
+        rgbDepth = 8;
+    }
+#endif
 
     volatile avifBool monochrome8bit = (avif->yuvFormat == AVIF_PIXEL_FORMAT_YUV400) && !avif->alphaPlane && (avif->depth == 8) &&
                                        (rgbDepth == 8);
diff --git a/include/avif/avif.h b/include/avif/avif.h
index 5a59e22..12fbec4 100644
--- a/include/avif/avif.h
+++ b/include/avif/avif.h
@@ -331,7 +331,12 @@
     AVIF_MATRIX_COEFFICIENTS_SMPTE2085 = 11,
     AVIF_MATRIX_COEFFICIENTS_CHROMA_DERIVED_NCL = 12,
     AVIF_MATRIX_COEFFICIENTS_CHROMA_DERIVED_CL = 13,
-    AVIF_MATRIX_COEFFICIENTS_ICTCP = 14
+    AVIF_MATRIX_COEFFICIENTS_ICTCP = 14,
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    AVIF_MATRIX_COEFFICIENTS_YCGCO_RE = 15,
+    AVIF_MATRIX_COEFFICIENTS_YCGCO_RO = 16,
+#endif
+    AVIF_MATRIX_COEFFICIENTS_LAST
 };
 typedef uint16_t avifMatrixCoefficients; // AVIF_MATRIX_COEFFICIENTS_*
 
diff --git a/include/avif/internal.h b/include/avif/internal.h
index ae0763f..050c64f 100644
--- a/include/avif/internal.h
+++ b/include/avif/internal.h
@@ -120,7 +120,11 @@
 {
     AVIF_REFORMAT_MODE_YUV_COEFFICIENTS = 0, // Normal YUV conversion using coefficients
     AVIF_REFORMAT_MODE_IDENTITY,             // Pack GBR directly into YUV planes (AVIF_MATRIX_COEFFICIENTS_IDENTITY)
-    AVIF_REFORMAT_MODE_YCGCO                 // YUV conversion using AVIF_MATRIX_COEFFICIENTS_YCGCO
+    AVIF_REFORMAT_MODE_YCGCO,                // YUV conversion using AVIF_MATRIX_COEFFICIENTS_YCGCO
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    AVIF_REFORMAT_MODE_YCGCO_RE, // YUV conversion using AVIF_MATRIX_COEFFICIENTS_YCGCO_RE
+    AVIF_REFORMAT_MODE_YCGCO_RO, // YUV conversion using AVIF_MATRIX_COEFFICIENTS_YCGCO_RO
+#endif
 } avifReformatMode;
 
 typedef enum avifAlphaMultiplyMode
diff --git a/src/reformat.c b/src/reformat.c
index c66f240..a8e8441 100644
--- a/src/reformat.c
+++ b/src/reformat.c
@@ -47,12 +47,20 @@
     //
     // YCgCo performs limited-full range adjustment on R,G,B but the current implementation performs range adjustment
     // on Y,U,V. So YCgCo with limited range is unsupported.
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    const avifBool useYCgCoRe = (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE);
+    const avifBool useYCgCoRo = (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO);
+#else
+    const avifBool useYCgCoRe = AVIF_FALSE;
+    const avifBool useYCgCoRo = AVIF_FALSE;
+#endif
     if ((image->matrixCoefficients == 3 /* CICP reserved */) ||
-        ((image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO) && (image->yuvRange == AVIF_RANGE_LIMITED)) ||
+        ((image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO || useYCgCoRe || useYCgCoRo) &&
+         (image->yuvRange == AVIF_RANGE_LIMITED)) ||
         (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_BT2020_CL) ||
         (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_SMPTE2085) ||
         (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_CHROMA_DERIVED_CL) ||
-        (image->matrixCoefficients >= AVIF_MATRIX_COEFFICIENTS_ICTCP)) { // Note the >= catching "future" CICP values here too
+        (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_ICTCP) || (image->matrixCoefficients >= AVIF_MATRIX_COEFFICIENTS_LAST)) {
         return AVIF_FALSE;
     }
 
@@ -68,6 +76,12 @@
         state->mode = AVIF_REFORMAT_MODE_IDENTITY;
     } else if (image->matrixCoefficients == AVIF_MATRIX_COEFFICIENTS_YCGCO) {
         state->mode = AVIF_REFORMAT_MODE_YCGCO;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    } else if (useYCgCoRe) {
+        state->mode = AVIF_REFORMAT_MODE_YCGCO_RE;
+    } else if (useYCgCoRo) {
+        state->mode = AVIF_REFORMAT_MODE_YCGCO_RO;
+#endif
     }
 
     if (state->mode != AVIF_REFORMAT_MODE_YUV_COEFFICIENTS) {
@@ -158,7 +172,13 @@
 
     // YCgCo performs limited-full range adjustment on R,G,B but the current implementation performs range adjustment
     // on Y,U,V. So YCgCo with limited range is unsupported.
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+    assert((state->mode != AVIF_REFORMAT_MODE_YCGCO && state->mode != AVIF_REFORMAT_MODE_YCGCO_RE &&
+            state->mode != AVIF_REFORMAT_MODE_YCGCO_RO) ||
+           (state->yuvRange == AVIF_RANGE_FULL));
+#else
     assert((state->mode != AVIF_REFORMAT_MODE_YCGCO) || (state->yuvRange == AVIF_RANGE_FULL));
+#endif
 
     if (state->mode == AVIF_REFORMAT_MODE_IDENTITY) {
         unorm = (int)avifRoundf(v * state->rangeY + state->biasY);
@@ -313,6 +333,19 @@
                             yuvBlock[bI][bJ].y = 0.5f * rgbPixel[1] + 0.25f * (rgbPixel[0] + rgbPixel[2]);
                             yuvBlock[bI][bJ].u = 0.5f * rgbPixel[1] - 0.25f * (rgbPixel[0] + rgbPixel[2]);
                             yuvBlock[bI][bJ].v = 0.5f * (rgbPixel[0] - rgbPixel[2]);
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+                        } else if (state.mode == AVIF_REFORMAT_MODE_YCGCO_RE || state.mode == AVIF_REFORMAT_MODE_YCGCO_RO) {
+                            // Formulas from JVET-U0093.
+                            const int R = (int)avifRoundf(AVIF_CLAMP(rgbPixel[0] * rgbMaxChannelF, 0.f, rgbMaxChannelF));
+                            const int G = (int)avifRoundf(AVIF_CLAMP(rgbPixel[1] * rgbMaxChannelF, 0.f, rgbMaxChannelF));
+                            const int B = (int)avifRoundf(AVIF_CLAMP(rgbPixel[2] * rgbMaxChannelF, 0.f, rgbMaxChannelF));
+                            const int Co = R - B;
+                            const int t = B + (Co >> 1);
+                            const int Cg = G - t;
+                            yuvBlock[bI][bJ].y = (t + (Cg >> 1)) / state.rangeY;
+                            yuvBlock[bI][bJ].u = Cg / state.rangeUV;
+                            yuvBlock[bI][bJ].v = Co / state.rangeUV;
+#endif
                         } else {
                             float Y = (kr * rgbPixel[0]) + (kg * rgbPixel[1]) + (kb * rgbPixel[2]);
                             yuvBlock[bI][bJ].y = Y;
@@ -701,6 +734,22 @@
                     G = Y + Cb;
                     B = t - Cr;
                     R = t + Cr;
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+                } else if (state->mode == AVIF_REFORMAT_MODE_YCGCO_RE || state->mode == AVIF_REFORMAT_MODE_YCGCO_RO) {
+                    const int bitOffset = (state->mode == AVIF_REFORMAT_MODE_YCGCO_RE) ? 2 : 1;
+                    const int maxValue = (1 << (state->yuvDepth - bitOffset)) - 1;
+                    assert((float)maxValue == rgbMaxChannelF);
+                    const int YY = unormY;
+                    const int Cg = (int)avifRoundf(Cb * yuvMaxChannel);
+                    const int Co = (int)avifRoundf(Cr * yuvMaxChannel);
+                    const int t = YY - (Cg >> 1);
+                    G = AVIF_CLAMP(t + Cg, 0, maxValue);
+                    B = AVIF_CLAMP(t - (Co >> 1), 0, maxValue);
+                    R = AVIF_CLAMP(B + Co, 0, maxValue);
+                    G /= rgbMaxChannelF;
+                    B /= rgbMaxChannelF;
+                    R /= rgbMaxChannelF;
+#endif
                 } else {
                     // Normal YUV
                     R = Y + (2 * (1 - kr)) * Cr;
diff --git a/tests/gtest/aviflosslesstest.cc b/tests/gtest/aviflosslesstest.cc
index 82c446c..a4e31c5 100644
--- a/tests/gtest/aviflosslesstest.cc
+++ b/tests/gtest/aviflosslesstest.cc
@@ -21,8 +21,12 @@
     const testutil::AvifImagePtr ground_truth_image =
         testutil::ReadImage(data_path, file_name);
 
-    for (auto matrix_coefficient :
-         {AVIF_MATRIX_COEFFICIENTS_IDENTITY, AVIF_MATRIX_COEFFICIENTS_YCGCO}) {
+    for (auto matrix_coefficient : {
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+           AVIF_MATRIX_COEFFICIENTS_YCGCO_RE, AVIF_MATRIX_COEFFICIENTS_YCGCO_RO,
+#endif
+               AVIF_MATRIX_COEFFICIENTS_IDENTITY, AVIF_MATRIX_COEFFICIENTS_YCGCO
+         }) {
       // Read a ground truth image but ask for certain matrix coefficients.
       testutil::AvifImagePtr image(avifImageCreateEmpty(), avifImageDestroy);
       ASSERT_NE(image, nullptr);
@@ -36,6 +40,14 @@
           /*ignoreXMP=*/false, image.get(),
           /*outDepth=*/nullptr, /*sourceTiming=*/nullptr,
           /*frameIter=*/nullptr);
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+      if (matrix_coefficient == AVIF_MATRIX_COEFFICIENTS_YCGCO_RO) {
+        // AVIF_MATRIX_COEFFICIENTS_YCGCO_RO does not work because the input
+        // depth is not odd.
+        ASSERT_EQ(file_format, AVIF_APP_FILE_FORMAT_UNKNOWN);
+        continue;
+      }
+#endif
       ASSERT_NE(file_format, AVIF_APP_FILE_FORMAT_UNKNOWN);
 
       // Encode.
@@ -70,7 +82,12 @@
       // Verify that the ground truth and decoded images are the same.
       const bool are_images_equal =
           testutil::AreImagesEqual(*ground_truth_image, *decoded_default);
-      if (matrix_coefficient == AVIF_MATRIX_COEFFICIENTS_IDENTITY) {
+
+      if (matrix_coefficient == AVIF_MATRIX_COEFFICIENTS_IDENTITY
+#if defined(AVIF_ENABLE_EXPERIMENTAL_YCGCO_R)
+          || matrix_coefficient == AVIF_MATRIX_COEFFICIENTS_YCGCO_RE
+#endif
+      ) {
         ASSERT_TRUE(are_images_equal);
       } else {
         // AVIF_MATRIX_COEFFICIENTS_YCGCO is not lossless because it does not