color_transformations.dart

// Copyright 2022 Google Inc. Use of this source code is governed by an
// MIT-style license that can be found at https://opensource.org/licenses/MIT.

import 'dart:typed_data';

import 'package:rational/rational.dart';

// Matrix values from https://www.w3.org/TR/css-color-4/#color-conversion-code.

final d65 = chromaToXyz(Rational.parse('0.3127'), Rational.parse('0.3290'));
final d50 = chromaToXyz(Rational.parse('0.3457'), Rational.parse('0.3585'));

final linearToXyzD65 = {
  "srgb": linearLightRgbToXyz(0.640, 0.330, 0.300, 0.600, 0.150, 0.060, d65),
  "display-p3":
      linearLightRgbToXyz(0.680, 0.320, 0.265, 0.690, 0.150, 0.060, d65),
  "a98-rgb":
      linearLightRgbToXyz(0.6400, 0.3300, 0.2100, 0.7100, 0.1500, 0.0600, d65),
  "rec2020": linearLightRgbToXyz(0.708, 0.292, 0.170, 0.797, 0.131, 0.046, d65)
};

final linearToXyzD50 = {
  "prophoto-rgb": linearLightRgbToXyz(
      0.734699, 0.265301, 0.159597, 0.840403, 0.036598, 0.000105, d50),
};

final bradford = RationalMatrix.fromFloat64List(Float64List.fromList([
  00.8951000, 00.2664000, -0.1614000, //
  -0.7502000, 01.7135000, 00.0367000,
  00.0389000, -0.0685000, 01.0296000
]));

/// The transformation matrix for converting D65 XYZ colors to D50 XYZ.
final RationalMatrix d65XyzToD50 = () {
  // Algorithm from http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html
  var source = bradford.timesVector(d65);
  var destination = bradford.timesVector(d50);
  return bradford.invert() *
      RationalMatrix([
        [destination[0] / source[0], Rational.zero, Rational.zero],
        [Rational.zero, destination[1] / source[1], Rational.zero],
        [Rational.zero, Rational.zero, destination[2] / source[2]]
      ]) *
      bradford;
}();

void main() {
  for (var src in linearToXyzD65.entries) {
    print("// ${src.key} to xyz-d65");
    print('var M = ${src.value.toJSString()};\n');

    print("// xyz-d65 to ${src.key}");
    print('var M = ${src.value.invert().toJSString()};\n');
  }

  for (var src in linearToXyzD50.entries) {
    print("// ${src.key} to xyz-d50");
    print('var M = ${src.value.toJSString()};\n');

    print("// xyz-d50 to ${src.key}");
    print('var M = ${src.value.invert().toJSString()};\n');
  }

  print("// xyz-d65 to xyz-d50");
  print('var M = ${d65XyzToD50.toJSString()};\n');

  print("// xyz-d50 to xyz-d65");
  print('var M = ${d65XyzToD50.invert().toJSString()};\n');
}

class RationalMatrix {
  final List<List<Rational>> contents;

  RationalMatrix(Iterable<Iterable<Rational>> contents)
      : contents = List.unmodifiable(
            contents.map((iter) => List<Rational>.unmodifiable(iter)));

  RationalMatrix.empty()
      : contents = List.generate(3, (_) => List.filled(3, Rational.zero));

  factory RationalMatrix.fromFloat64List(Float64List list) =>
      RationalMatrix(List.generate(
          3,
          (i) => List.generate(
              3, (j) => Rational.parse(list[i * 3 + j].toString()))));

  RationalMatrix operator *(RationalMatrix other) => RationalMatrix([
        for (var i = 0; i < 3; i++)
          [
            for (var j = 0; j < 3; j++)
              [for (var k = 0; k < 3; k++) get(i, k) * other.get(k, j)].sum
          ]
      ]);

  List<Rational> timesVector(List<Rational> vector) => List.generate(
      3, (i) => Iterable.generate(3, (j) => get(i, j) * vector[j]).sum);

  RationalMatrix invert() {
    // Using the same naming convention used in
    // https://en.wikipedia.org/wiki/Determinant and
    // https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices.
    var a = get(0, 0);
    var b = get(0, 1);
    var c = get(0, 2);
    var d = get(1, 0);
    var e = get(1, 1);
    var f = get(1, 2);
    var g = get(2, 0);
    var h = get(2, 1);
    var i = get(2, 2);

    var idet = Rational.one /
        (a * e * i + b * f * g + c * d * h - c * e * g - b * d * i - a * f * h);

    return RationalMatrix([
      [(e * i - f * h) * idet, -(b * i - c * h) * idet, (b * f - c * e) * idet],
      [
        -(d * i - f * g) * idet,
        (a * i - c * g) * idet,
        -(a * f - c * d) * idet
      ],
      [(d * h - e * g) * idet, -(a * h - b * g) * idet, (a * e - b * d) * idet],
    ]);
  }

  RationalMatrix transpose() => RationalMatrix(
      List.generate(3, (i) => List.generate(3, (j) => get(j, i))));

  Rational get(int i, int j) => contents[i][j];

  Rational set(int i, int j, Rational value) => contents[i][j] = value;

  String toString() =>
      '[ ' +
      contents
          .map((row) => row.map((number) => number.toDoubleString()).join(' '))
          .join('\n  ') +
      ' ]';

  String toExactString() =>
      '[ ' +
      contents
          .map((row) => row.map((number) => number.toExactString()).join(' '))
          .join('\n  ') +
      ' ]';

  String toJSString() =>
      '[\n\t' +
      contents
          .map((row) => row.map((number) => number.toDoubleString()).join(', '))
          .join(',\n\t') +
      '\n]';
}

const precision = 17;

extension on Rational {
  String toDoubleString() {
    var doubleString = double.parse(toExactString()).toString();
    if (!doubleString.startsWith('-')) doubleString = ' $doubleString';
    return doubleString.toString().padRight(precision + 3, '0');
  }

  String toExactString() {
    var newNum = (Rational(numerator) *
            Rational(BigInt.from(10).pow(precision), denominator))
        .truncate();

    var numString = newNum.abs().toString();
    if (numString.length == precision + 1) {
      numString = '${numString[0]}.${numString.substring(1)}';
    } else {
      numString = '0.${numString.padLeft(precision, '0')}';
    }

    return '${newNum.isNegative ? '-' : '0'}$numString';
  }
}

extension on Iterable<Rational> {
  Rational get sum => reduce((a, b) => a + b);
}

RationalMatrix linearLightRgbToXyz(
        double redChromaX,
        double redChromaY,
        double greenChromaX,
        double greenChromaY,
        double blueChromaX,
        double blueChromaY,
        List<Rational> white) =>
    _linearLightRgbToXyz(
        Rational.parse(redChromaX.toString()),
        Rational.parse(redChromaY.toString()),
        Rational.parse(greenChromaX.toString()),
        Rational.parse(greenChromaY.toString()),
        Rational.parse(blueChromaX.toString()),
        Rational.parse(blueChromaY.toString()),
        white);

RationalMatrix _linearLightRgbToXyz(
    Rational redChromaX,
    Rational redChromaY,
    Rational greenChromaX,
    Rational greenChromaY,
    Rational blueChromaX,
    Rational blueChromaY,
    List<Rational> white) {
  // Algorithm from http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
  var xyzRed = chromaToXyz(redChromaX, redChromaY);
  var xyzGreen = chromaToXyz(greenChromaX, greenChromaY);
  var xyzBlue = chromaToXyz(blueChromaX, blueChromaY);

  var s = RationalMatrix([xyzRed, xyzGreen, xyzBlue])
      .transpose()
      .invert()
      .timesVector(white);
  var sRed = s[0];
  var sGreen = s[1];
  var sBlue = s[2];

  return RationalMatrix([
    xyzRed.map((value) => sRed * value),
    xyzGreen.map((value) => sGreen * value),
    xyzBlue.map((value) => sBlue * value)
  ]).transpose();
}

/// Convert a two-dimensional chroma coordinates into a point in XYZ space.
List<Rational> chromaToXyz(Rational chromaX, Rational chromaY) => [
      chromaX / chromaY,
      Rational.one,
      (Rational.one - chromaX - chromaY) / chromaY
    ];