// Copyright 2023 Yuan Tong. All rights reserved.
// SPDX-License-Identifier: BSD-2-Clause

#include "iccmaker.h"

#include <math.h>
#include <string.h>

// ICCv2 profile specification: https://www.color.org/icc32.pdf

/**
 * Color Profile Structure
 *
 * Header:
 *  size         = 376 bytes (*1)
 *  CMM          = 'lcms' (*2)
 *  Version      = 2.2.0
 *  Device Class = Display
 *  Color Space  = RGB
 *  Conn. Space  = XYZ
 *  Date, Time   = 1 Jan 2000, 0:00:00
 *  Platform     = Microsoft
 *  Flags        = Not Embedded Profile, Use anywhere
 *  Dev. Mnfctr. = 0x0
 *  Dev. Model   = 0x0
 *  Dev. Attrbts = Reflective, Glossy, Positive, Color
 *  Rndrng Intnt = Perceptual
 *  Illuminant   = 0.96420288, 1.00000000, 0.82490540    [Lab 100.000000, 0.000000, 0.000000]
 *  Creator      = 'avif'
 *
 * Profile Tags:
 *                    Tag    ID      Offset         Size                 Value
 *                   ----  ------    ------         ----                 -----
 *  profileDescriptionTag  'desc'       240           95                  avif
 *     mediaWhitePointTag  'wtpt'       268 (*3)      20        (to be filled)
 *         redColorantTag  'rXYZ'       288           20        (to be filled)
 *       greenColorantTag  'gXYZ'       308           20        (to be filled)
 *        blueColorantTag  'bXYZ'       328           20        (to be filled)
 *              redTRCTag  'rTRC'       348 (*4)      16        (to be filled)
 *            greenTRCTag  'gTRC'       348           16        (to be filled)
 *             blueTRCTag  'bTRC'       348           16        (to be filled)
 *           copyrightTag  'cprt'       364           12                   CC0
 *
 * (*1): The template data is padded to 448 bytes according to MD5 specification, so that computation can be applied
 *       directly on it. The actual ICC profile data is the first 376 bytes.
 * (*2): 6.1.2 CMM Type: The signatures must be registered in order to avoid conflicts.
 *       The registry can be found at https://www.color.org/signatures2.xalter (Private and ICC tag and CMM registry)
 *       Therefore we are using the signature of Little CMS.
 * (*3): The profileDescriptionTag requires 95 bytes of data, but with some trick, the content of the last 67 bytes
 *       can be anything. Therefore we are placing the following tags in this region to reduce profile size.
 * (*4): The transfer characteristic (gamma) of the 3 channels are the same, so the data can be shared.
 */

static const uint8_t iccColorTemplate[448] = {
    0x00, 0x00, 0x01, 0x78, 0x6c, 0x63, 0x6d, 0x73, 0x02, 0x20, 0x00, 0x00, 0x6d, 0x6e, 0x74, 0x72, 0x52, 0x47, 0x42, 0x20, 0x58,
    0x59, 0x5a, 0x20, 0x07, 0xd0, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x61, 0x63, 0x73, 0x70, 0x4d, 0x53,
    0x46, 0x54, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf6, 0xd6, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0xd3, 0x2d, 0x61, 0x76, 0x69, 0x66,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x09, 0x64, 0x65, 0x73, 0x63, 0x00, 0x00, 0x00, 0xf0, 0x00, 0x00, 0x00, 0x5f, 0x77, 0x74, 0x70,
    0x74, 0x00, 0x00, 0x01, 0x0c, 0x00, 0x00, 0x00, 0x14, 0x72, 0x58, 0x59, 0x5a, 0x00, 0x00, 0x01, 0x20, 0x00, 0x00, 0x00, 0x14,
    0x67, 0x58, 0x59, 0x5a, 0x00, 0x00, 0x01, 0x34, 0x00, 0x00, 0x00, 0x14, 0x62, 0x58, 0x59, 0x5a, 0x00, 0x00, 0x01, 0x48, 0x00,
    0x00, 0x00, 0x14, 0x72, 0x54, 0x52, 0x43, 0x00, 0x00, 0x01, 0x5c, 0x00, 0x00, 0x00, 0x10, 0x67, 0x54, 0x52, 0x43, 0x00, 0x00,
    0x01, 0x5c, 0x00, 0x00, 0x00, 0x10, 0x62, 0x54, 0x52, 0x43, 0x00, 0x00, 0x01, 0x5c, 0x00, 0x00, 0x00, 0x10, 0x63, 0x70, 0x72,
    0x74, 0x00, 0x00, 0x01, 0x6c, 0x00, 0x00, 0x00, 0x0c, 0x64, 0x65, 0x73, 0x63, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x05,
    0x61, 0x76, 0x69, 0x66, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x58, 0x59, 0x5a, 0x20, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0xf3, 0x54, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x16, 0xc9, 0x58, 0x59, 0x5a, 0x20, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x6f, 0xa0, 0x00, 0x00, 0x38, 0xf2, 0x00, 0x00, 0x03, 0x8f, 0x58, 0x59, 0x5a, 0x20, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x62, 0x96, 0x00, 0x00, 0xb7, 0x89, 0x00, 0x00, 0x18, 0xda, 0x58, 0x59, 0x5a, 0x20, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x24, 0xa0, 0x00, 0x00, 0x0f, 0x85, 0x00, 0x00, 0xb6, 0xc4, 0x63, 0x75, 0x72, 0x76, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x74, 0x65, 0x78, 0x74, 0x00, 0x00, 0x00, 0x00, 0x43, 0x43, 0x30, 0x00, 0x80, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xc0,
    0x0b, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};

static const size_t iccColorLength = 376;

static const ptrdiff_t colorWhiteOffset = 0x114;
static const ptrdiff_t colorRedOffset = 0x128;
static const ptrdiff_t colorGreenOffset = 0x13c;
static const ptrdiff_t colorBlueOffset = 0x150;
static const ptrdiff_t colorGammaOffset = 0x168;

/**
 * Gray Profile Structure
 *
 * Header:
 *  size         = 275 bytes
 *  CMM          = 'lcms'
 *  Version      = 2.2.0
 *  Device Class = Display
 *  Color Space  = Gray
 *  Conn. Space  = XYZ
 *  Date, Time   = 1 Jan 2000, 0:00:00
 *  Platform     = Microsoft
 *  Flags        = Not Embedded Profile, Use anywhere
 *  Dev. Mnfctr. = 0x0
 *  Dev. Model   = 0x0
 *  Dev. Attrbts = Reflective, Glossy, Positive, Color
 *  Rndrng Intnt = Perceptual
 *  Illuminant   = 0.96420288, 1.00000000, 0.82490540    [Lab 100.000000, 0.000000, 0.000000]
 *  Creator      = 'avif'
 *
 * Profile Tags:
 *                    Tag    ID      Offset         Size                 Value
 *                   ----  ------    ------         ----                 -----
 *  profileDescriptionTag  'desc'       180           95                  avif
 *     mediaWhitePointTag  'wtpt'       208           20        (to be filled)
 *             grayTRCTag  'kTRC'       228           16        (to be filled)
 *           copyrightTag  'cprt'       244           12                   CC0
 */

static const uint8_t iccGrayTemplate[320] = {
    0x00, 0x00, 0x01, 0x13, 0x6c, 0x63, 0x6d, 0x73, 0x02, 0x20, 0x00, 0x00, 0x6d, 0x6e, 0x74, 0x72, 0x47, 0x52, 0x41, 0x59,
    0x58, 0x59, 0x5a, 0x20, 0x07, 0xd0, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x61, 0x63, 0x73, 0x70,
    0x4d, 0x53, 0x46, 0x54, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf6, 0xd6, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0xd3, 0x2d,
    0x61, 0x76, 0x69, 0x66, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x04, 0x64, 0x65, 0x73, 0x63, 0x00, 0x00, 0x00, 0xb4,
    0x00, 0x00, 0x00, 0x5f, 0x77, 0x74, 0x70, 0x74, 0x00, 0x00, 0x00, 0xd0, 0x00, 0x00, 0x00, 0x14, 0x6b, 0x54, 0x52, 0x43,
    0x00, 0x00, 0x00, 0xe4, 0x00, 0x00, 0x00, 0x10, 0x63, 0x70, 0x72, 0x74, 0x00, 0x00, 0x00, 0xf4, 0x00, 0x00, 0x00, 0x0c,
    0x64, 0x65, 0x73, 0x63, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x05, 0x61, 0x76, 0x69, 0x66, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x58, 0x59, 0x5a, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf3, 0x54,
    0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x16, 0xc9, 0x63, 0x75, 0x72, 0x76, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01,
    0x01, 0x00, 0x00, 0x00, 0x74, 0x65, 0x78, 0x74, 0x00, 0x00, 0x00, 0x00, 0x43, 0x43, 0x30, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x80, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x98, 0x08, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};

static const size_t iccGrayLength = 275;

static const ptrdiff_t grayWhiteOffset = 0xd8;
static const ptrdiff_t grayGammaOffset = 0xf0;

static const ptrdiff_t checksumOffset = 0x54;

static const double small = 1e-12;

static uint32_t readLittleEndianU32(const uint8_t * data)
{
    return ((uint32_t)data[0] << 0) | ((uint32_t)data[1] << 8) | ((uint32_t)data[2] << 16) | ((uint32_t)data[3] << 24);
}

static void writeLittleEndianU32(uint8_t * data, uint32_t value)
{
    data[0] = (value >> 0) & 0xff;
    data[1] = (value >> 8) & 0xff;
    data[2] = (value >> 16) & 0xff;
    data[3] = (value >> 24) & 0xff;
}

static void writeBigEndianU16(uint8_t * data, uint16_t value)
{
    data[0] = (value >> 8) & 0xff;
    data[1] = (value >> 0) & 0xff;
}

static void writeBigEndianU32(uint8_t * data, uint32_t value)
{
    data[0] = (value >> 24) & 0xff;
    data[1] = (value >> 16) & 0xff;
    data[2] = (value >> 8) & 0xff;
    data[3] = (value >> 0) & 0xff;
}

static avifBool putS15Fixed16(uint8_t * data, double value)
{
    value = round(value * 65536);
    if (value > INT32_MAX || value < INT32_MIN) {
        return AVIF_FALSE;
    }

    int32_t fixed = (int32_t)value;
    // reinterpret into uint32_t to ensure the exact bits are written.
    writeBigEndianU32(data, *(uint32_t *)&fixed);

    return AVIF_TRUE;
}

static avifBool putU8Fixed8(uint8_t * data, float value)
{
    value = roundf(value * 256);
    if (value > UINT16_MAX || value < 1) {
        return AVIF_FALSE;
    }

    uint16_t fixed = (uint16_t)value;
    writeBigEndianU16(data, fixed);
    return AVIF_TRUE;
}

static avifBool putColorant(uint8_t * data, const double XYZ[3])
{
    if (!putS15Fixed16(data, XYZ[0])) {
        return AVIF_FALSE;
    }

    if (!putS15Fixed16(data + 4, XYZ[1])) {
        return AVIF_FALSE;
    }

    if (!putS15Fixed16(data + 8, XYZ[2])) {
        return AVIF_FALSE;
    }

    return AVIF_TRUE;
}

static avifBool xyToXYZ(const float xy[2], double XYZ[3])
{
    if (fabsf(xy[1]) < small) {
        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 matInv(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) < small) {
        return AVIF_FALSE;
    }
    det = 1 / 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 matMul(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 matDiag(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];
}

static void swap(double * a, double * b)
{
    double tmp = *a;
    *a = *b;
    *b = tmp;
}

// Transpose M
static void matTrans(double M[3][3])
{
    swap(&M[0][1], &M[1][0]);
    swap(&M[0][2], &M[2][0]);
    swap(&M[1][2], &M[2][1]);
}

// Computes y = M.x
static void vecMul(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];
}

// MD5 algorithm. See https://www.ietf.org/rfc/rfc1321.html#appendix-A.3
// This function writes the MD5 checksum in place at offset `checksumOffset` of `data`.
// This function shall only be called with a copy of iccColorTemplate or iccGrayTemplate, and sizeof(icc*Template).
static void computeMD5(uint8_t * data, size_t length)
{
    static const uint32_t sineparts[64] = {
        0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501, 0x698098d8, 0x8b44f7af,
        0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa,
        0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8,
        0x676f02d9, 0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70,
        0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665, 0xf4292244, 0x432aff97,
        0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1,
        0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391,
    };
    static const uint8_t shift[64] = {
        7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 5, 9,  14, 20, 5, 9,  14, 20, 5, 9,  14, 20, 5, 9,  14, 20,
        4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21,
    };

    uint32_t a0 = 0x67452301, b0 = 0xefcdab89, c0 = 0x98badcfe, d0 = 0x10325476;

    for (uint32_t i = 0; i < length; i += 64) {
        uint32_t a = a0, b = b0, c = c0, d = d0, f, g;
        for (uint32_t j = 0; j < 64; j++) {
            if (j < 16) {
                f = (b & c) | ((~b) & d);
                g = j;
            } else if (j < 32) {
                f = (d & b) | ((~d) & c);
                g = (5 * j + 1) & 0xf;
            } else if (j < 48) {
                f = b ^ c ^ d;
                g = (3 * j + 5) & 0xf;
            } else {
                f = c ^ (b | (~d));
                g = (7 * j) & 0xf;
            }
            uint32_t u = readLittleEndianU32(data + i + g * 4);
            f += a + sineparts[j] + u;
            a = d;
            d = c;
            c = b;
            b += (f << shift[j]) | (f >> (32u - shift[j]));
        }
        a0 += a;
        b0 += b;
        c0 += c;
        d0 += d;
    }

    uint8_t * output = data + checksumOffset;
    writeLittleEndianU32(output, a0);
    writeLittleEndianU32(output + 4, b0);
    writeLittleEndianU32(output + 8, c0);
    writeLittleEndianU32(output + 12, d0);
}

// Bradford chromatic adaptation matrix
// from https://www.researchgate.net/publication/253799640_A_uniform_colour_space_based_upon_CIECAM97s
static const double bradford[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 lmsD50[3] = { 0.996284, 1.02043, 0.818644 };

avifBool avifGenerateRGBICC(avifRWData * icc, float gamma, const float primaries[8])
{
    uint8_t buffer[sizeof(iccColorTemplate)];
    memcpy(buffer, iccColorTemplate, sizeof(iccColorTemplate));

    double whitePointXYZ[3];
    if (!xyToXYZ(&primaries[6], whitePointXYZ)) {
        return AVIF_FALSE;
    }

    if (!putColorant(buffer + colorWhiteOffset, whitePointXYZ)) {
        return AVIF_FALSE;
    }

    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];
    if (!matInv(rgbPrimaries, rgbPrimariesInv)) {
        return AVIF_FALSE;
    }

    double rgbCoefficients[3];
    vecMul(rgbPrimariesInv, whitePointXYZ, rgbCoefficients);

    double rgbCoefficientsMat[3][3];
    matDiag(rgbCoefficients, rgbCoefficientsMat);

    double rgbXYZ[3][3];
    matMul(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];
    vecMul(bradford, whitePointXYZ, lms);
    for (int i = 0; i < 3; ++i) {
        if (fabs(lms[i]) < small) {
            return AVIF_FALSE;
        }
        lms[i] = lmsD50[i] / lms[i];
    }

    double adaptation[3][3];
    matDiag(lms, adaptation);

    double tmp[3][3];
    matMul(adaptation, bradford, tmp);

    double bradfordInv[3][3];
    if (!matInv(bradford, bradfordInv)) {
        return AVIF_FALSE;
    }
    matMul(bradfordInv, tmp, adaptation);

    double rgbXYZD50[3][3];
    matMul(adaptation, rgbXYZ, rgbXYZD50);
    matTrans(rgbXYZD50);

    if (!putColorant(buffer + colorRedOffset, rgbXYZD50[0])) {
        return AVIF_FALSE;
    }

    if (!putColorant(buffer + colorGreenOffset, rgbXYZD50[1])) {
        return AVIF_FALSE;
    }

    if (!putColorant(buffer + colorBlueOffset, rgbXYZD50[2])) {
        return AVIF_FALSE;
    }

    if (!putU8Fixed8(buffer + colorGammaOffset, gamma)) {
        return AVIF_FALSE;
    }

    computeMD5(buffer, sizeof(iccColorTemplate));
    if (avifRWDataSet(icc, buffer, iccColorLength) != AVIF_RESULT_OK) {
        return AVIF_FALSE;
    }

    return AVIF_TRUE;
}

avifBool avifGenerateGrayICC(avifRWData * icc, float gamma, const float white[2])
{
    uint8_t buffer[sizeof(iccGrayTemplate)];
    memcpy(buffer, iccGrayTemplate, sizeof(iccGrayTemplate));

    double whitePointXYZ[3];
    if (!xyToXYZ(white, whitePointXYZ)) {
        return AVIF_FALSE;
    }

    if (!putColorant(buffer + grayWhiteOffset, whitePointXYZ)) {
        return AVIF_FALSE;
    }

    if (!putU8Fixed8(buffer + grayGammaOffset, gamma)) {
        return AVIF_FALSE;
    }

    computeMD5(buffer, sizeof(iccGrayTemplate));
    if (avifRWDataSet(icc, buffer, iccGrayLength) != AVIF_RESULT_OK) {
        return AVIF_FALSE;
    }

    return AVIF_TRUE;
}