Add specialized matrix transform for Vector3 and Matrix3x3 (#585)

* Add specialized matrix transform for Vector3 and Matrix3x3

Adds a specialized matrix multiplication function that
transforms an array of length 3 with a 3x3 matrix. The
specialized function can also transform the array in place
which can eliminate one array allocation.

The https://github.com/texel-org/color color library
has a benchmark comparing itself to Color.js.

Profiling the conversion benchmark showed that matrix
multiplication was one of the areas that was taking
the most amount of time.

This PR speeds up the conversion benchmark by more
than 3x and some of the other benchmarks by 2-3x

Benchmark numbers before this PR:
---------------------------------

conversion (Colorjs.io procedural API) --
Colorjs.io: 6213.74 ms
Ours: 733.75 ms
Speedup: 8.5x faster

conversion (Colorjs.io main API) --
Colorjs.io: 25444.47 ms
Ours: 703.50 ms
Speedup: 36.2x faster

gamut mapping OKLCH - sRGB (Colorjs.io procedural API) --
Colorjs.io: 13302.80 ms
Ours: 72.97 ms
Speedup: 182.3x faster

gamut mapping OKLCH - sRGB (Colorjs.io main API) --
Colorjs.io: 14900.11 ms
Ours: 67.78 ms
Speedup: 219.8x faster

gamut mapping all spaces to P3 (Colorjs.io procedural API) --
Colorjs.io: 10878.26 ms
Ours: 201.85 ms
Speedup: 53.9x faster

gamut mapping all spaces to P3 (Colorjs.io main API) --
Colorjs.io: 13731.21 ms
Ours: 201.76 ms
Speedup: 68.1x faster

After this PR:
--------------

conversion (Colorjs.io procedural API) --
Colorjs.io: 1793.97 ms
Ours: 713.03 ms
Speedup: 2.5x faster

conversion (Colorjs.io main API) --
Colorjs.io: 19884.18 ms
Ours: 704.59 ms
Speedup: 28.2x faster

gamut mapping OKLCH - sRGB (Colorjs.io procedural API) --
Colorjs.io: 4981.04 ms
Ours: 74.01 ms
Speedup: 67.3x faster

gamut mapping OKLCH - sRGB (Colorjs.io main API) --
Colorjs.io: 6477.12 ms
Ours: 68.02 ms
Speedup: 95.2x faster

gamut mapping all spaces to P3 (Colorjs.io procedural API) --
Colorjs.io: 4205.36 ms
Ours: 202.27 ms
Speedup: 20.8x faster

gamut mapping all spaces to P3 (Colorjs.io main API) --
Colorjs.io: 6137.99 ms
Ours: 199.91 ms
Speedup: 30.7x faster

* Fix formatting

* Fix type test

* Rename transform to multiply_v3_m3x3

* Fix tests

src/RGBColorSpace.js

@@ -1,5 +1,5 @@
 import ColorSpace from "./ColorSpace.js";
-import {multiplyMatrices} from "./util.js";
+import {multiply_v3_m3x3} from "./util.js";
 import adapt from "./adapt.js";
 import XYZ_D65 from "./spaces/xyz-d65.js";
 
@@ -39,7 +39,7 @@ export default class RGBColorSpace extends ColorSpace {
 
 		if (options.toXYZ_M && options.fromXYZ_M) {
 			options.toBase ??= rgb => {
-				let xyz = /** @type {[number, number, number]} */ (multiplyMatrices(options.toXYZ_M, rgb));
+				let xyz = multiply_v3_m3x3(rgb, options.toXYZ_M);
 
 				if (this.white !== this.base.white) {
 					// Perform chromatic adaptation
@@ -51,7 +51,7 @@ export default class RGBColorSpace extends ColorSpace {
 
 			options.fromBase ??= xyz => {
 				xyz = adapt(this.base.white, this.white, xyz);
-				return multiplyMatrices(options.fromXYZ_M, xyz);
+				return multiply_v3_m3x3(xyz, options.fromXYZ_M);
 			};
 		}
 

src/adapt.js

@@ -1,5 +1,5 @@
 import hooks from "./hooks.js";
-import {multiplyMatrices} from "./util.js";
+import {multiply_v3_m3x3} from "./util.js";
 
 // Type "imports"
 /** @typedef {import("./types.js").White} White */
@@ -70,7 +70,7 @@ export default function adapt (W1, W2, XYZ, options = {}) {
 	hooks.run("chromatic-adaptation-end", env);
 
 	if (env.M) {
-		return /** @type {[number, number, number]} */ (multiplyMatrices(/** @type {number[][]}*/ (env.M), env.XYZ));
+		return multiply_v3_m3x3(env.XYZ, env.M);
 	}
 	else {
 		throw new TypeError("Only Bradford CAT with white points D50 and D65 supported for now.");

src/multiply-matrices.js

@@ -1,3 +1,8 @@
+// Type "imports"
+/** @typedef {import("./types.js").Matrix3x3} Matrix3x3 */
+/** @typedef {import("./types.js").Vector3} Vector3 */
+
+
 /**
  * A is m x n. B is n x p. product is m x p.
  *
@@ -85,3 +90,58 @@ export default function multiplyMatrices (A, B) {
 
 	return product;
 }
+
+
+// dot3 and transform functions adapted from https://github.com/texel-org/color/blob/9793c7d4d02b51f068e0f3fd37131129a4270396/src/core.js
+//
+// The MIT License (MIT)
+// Copyright (c) 2024 Matt DesLauriers
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+// DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+// OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
+// OR OTHER DEALINGS IN THE SOFTWARE.
+
+
+/**
+ * Returns the dot product of two vectors each with a length of 3.
+ *
+ * @param {Vector3} a
+ * @param {Vector3} b
+ * @returns {number}
+ */
+function dot3 (a, b) {
+	return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+/**
+ * Transforms a vector of length 3 by a 3x3 matrix. Specify the same input and output
+ * vector to transform in place.
+ *
+ * @param {Vector3} input
+ * @param {Matrix3x3} matrix
+ * @param {Vector3} [out]
+ * @returns {Vector3}
+*/
+export function multiply_v3_m3x3 (input, matrix, out = [0, 0, 0]) {
+	const x = dot3(input, matrix[0]);
+	const y = dot3(input, matrix[1]);
+	const z = dot3(input, matrix[2]);
+	out[0] = x;
+	out[1] = y;
+	out[2] = z;
+	return out;
+}

src/spaces/a98rgb-linear.js

@@ -1,17 +1,25 @@
 import RGBColorSpace from "../RGBColorSpace.js";
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
 // convert an array of linear-light a98-rgb values to CIE XYZ
 // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
 // has greater numerical precision than section 4.3.5.3 of
 // https://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf
 // but the values below were calculated from first principles
 // from the chromaticity coordinates of R G B W
+
+
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[ 0.5766690429101305,   0.1855582379065463,   0.1882286462349947  ],
 	[ 0.29734497525053605,  0.6273635662554661,   0.07529145849399788 ],
 	[ 0.02703136138641234,  0.07068885253582723,  0.9913375368376388  ],
 ];
 
+
+/** @type {Matrix3x3} */
 const fromXYZ_M = [
 	[  2.0415879038107465,    -0.5650069742788596,   -0.34473135077832956 ],
 	[ -0.9692436362808795,     1.8759675015077202,    0.04155505740717557 ],

src/spaces/acescg.js

@@ -2,18 +2,23 @@ import RGBColorSpace from "../RGBColorSpace.js";
 import {WHITES} from "../adapt.js";
 import "../CATs.js"; // because of the funky whitepoint
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
 // The ACES whitepoint
 // see TB-2018-001 Derivation of the ACES White Point CIE Chromaticity Coordinates
 // also https://github.com/ampas/aces-dev/blob/master/documents/python/TB-2018-001/aces_wp.py
 // Similar to D60
 WHITES.ACES = [0.32168 / 0.33767, 1.00000, (1.00000 - 0.32168 - 0.33767) / 0.33767];
 
 // convert an array of linear-light ACEScc values to CIE XYZ
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[  0.6624541811085053,   0.13400420645643313,  0.1561876870049078  ],
 	[  0.27222871678091454,  0.6740817658111484,   0.05368951740793705 ],
 	[ -0.005574649490394108, 0.004060733528982826, 1.0103391003129971  ],
 ];
+/** @type {Matrix3x3} */
 const fromXYZ_M = [
 	[  1.6410233796943257,   -0.32480329418479,    -0.23642469523761225  ],
 	[ -0.6636628587229829,    1.6153315916573379,   0.016756347685530137 ],

src/spaces/oklab.js

@@ -1,26 +1,35 @@
 import ColorSpace from "../ColorSpace.js";
-import {multiplyMatrices} from "../util.js";
+import {multiply_v3_m3x3} from "../util.js";
 import XYZ_D65 from "./xyz-d65.js";
 
+
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
+
 // Recalculated for consistent reference white
 // see https://github.com/w3c/csswg-drafts/issues/6642#issuecomment-943521484
+/** @type {Matrix3x3} */
 const XYZtoLMS_M = [
 	[ 0.8190224379967030, 0.3619062600528904, -0.1288737815209879 ],
 	[ 0.0329836539323885, 0.9292868615863434,  0.0361446663506424 ],
 	[ 0.0481771893596242, 0.2642395317527308,  0.6335478284694309 ],
 ];
 // inverse of XYZtoLMS_M
+/** @type {Matrix3x3} */
 const LMStoXYZ_M = [
 	[  1.2268798758459243, -0.5578149944602171,  0.2813910456659647 ],
 	[ -0.0405757452148008,  1.1122868032803170, -0.0717110580655164 ],
 	[ -0.0763729366746601, -0.4214933324022432,  1.5869240198367816 ],
 ];
+/** @type {Matrix3x3} */
 export const LMStoLab_M = [
 	[ 0.2104542683093140,  0.7936177747023054, -0.0040720430116193 ],
 	[ 1.9779985324311684, -2.4285922420485799,  0.4505937096174110 ],
 	[ 0.0259040424655478,  0.7827717124575296, -0.8086757549230774 ],
 ];
 // LMStoIab_M inverted
+/** @type {Matrix3x3} */
 export const LabtoLMS_M = [
 	[ 1.0000000000000000,  0.3963377773761749,  0.2158037573099136 ],
 	[ 1.0000000000000000, -0.1055613458156586, -0.0638541728258133 ],
@@ -48,22 +57,25 @@ export default new ColorSpace({
 	base: XYZ_D65,
 	fromBase (XYZ) {
 		// move to LMS cone domain
-		let LMS = multiplyMatrices(XYZtoLMS_M, XYZ);
+		let LMS = multiply_v3_m3x3(XYZ, XYZtoLMS_M);
 
 		// non-linearity
-		let LMSg = LMS.map(val => Math.cbrt(val));
-
-		return multiplyMatrices(LMStoLab_M, LMSg);
+		LMS[0] = Math.cbrt(LMS[0]);
+		LMS[1] = Math.cbrt(LMS[1]);
+		LMS[2] = Math.cbrt(LMS[2]);
 
+		return multiply_v3_m3x3(LMS, LMStoLab_M, LMS);
 	},
 	toBase (OKLab) {
 		// move to LMS cone domain
-		let LMSg = multiplyMatrices(LabtoLMS_M, OKLab);
+		let LMSg = multiply_v3_m3x3(OKLab, LabtoLMS_M);
 
 		// restore linearity
-		let LMS = LMSg.map(val => val ** 3);
+		LMSg[0] = LMSg[0] ** 3;
+		LMSg[1] = LMSg[1] ** 3;
+		LMSg[2] = LMSg[2] ** 3;
 
-		return multiplyMatrices(LMStoXYZ_M, LMS);
+		return multiply_v3_m3x3(LMSg, LMStoXYZ_M, LMSg);
 	},
 
 	formats: {

src/spaces/p3-linear.js

@@ -1,11 +1,16 @@
 import RGBColorSpace from "../RGBColorSpace.js";
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[0.4865709486482162, 0.26566769316909306, 0.1982172852343625],
 	[0.2289745640697488, 0.6917385218365064,  0.079286914093745],
 	[0.0000000000000000, 0.04511338185890264, 1.043944368900976],
 ];
 
+/** @type {Matrix3x3} */
 const fromXYZ_M = [
 	[ 2.493496911941425,   -0.9313836179191239, -0.40271078445071684],
 	[-0.8294889695615747,   1.7626640603183463,  0.023624685841943577],

src/spaces/prophoto-linear.js

@@ -1,16 +1,21 @@
 import RGBColorSpace from "../RGBColorSpace.js";
 import XYZ_D50 from "./xyz-d50.js";
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
 // convert an array of  prophoto-rgb values to CIE XYZ
 // using  D50 (so no chromatic adaptation needed afterwards)
 // matrix cannot be expressed in rational form, but is calculated to 64 bit accuracy
 // see https://github.com/w3c/csswg-drafts/issues/7675
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[ 0.79776664490064230,  0.13518129740053308,  0.03134773412839220 ],
 	[ 0.28807482881940130,  0.71183523424187300,  0.00008993693872564 ],
 	[ 0.00000000000000000,  0.00000000000000000,  0.82510460251046020 ],
 ];
 
+/** @type {Matrix3x3} */
 const fromXYZ_M = [
 	[  1.34578688164715830, -0.25557208737979464, -0.05110186497554526 ],
 	[ -0.54463070512490190,  1.50824774284514680,  0.02052744743642139 ],

src/spaces/rec2020-linear.js

@@ -1,16 +1,21 @@
 import RGBColorSpace from "../RGBColorSpace.js";
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
 // convert an array of linear-light rec2020 values to CIE XYZ
 // using  D65 (no chromatic adaptation)
 // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
 // 0 is actually calculated as  4.994106574466076e-17
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[ 0.6369580483012914, 0.14461690358620832,  0.1688809751641721  ],
 	[ 0.2627002120112671, 0.6779980715188708,   0.05930171646986196 ],
 	[ 0.000000000000000,  0.028072693049087428, 1.060985057710791   ],
 ];
 
 // from ITU-R BT.2124-0 Annex 2 p.3
+/** @type {Matrix3x3} */
 const fromXYZ_M = [
 	[  1.716651187971268,  -0.355670783776392, -0.253366281373660  ],
 	[ -0.666684351832489,   1.616481236634939,  0.0157685458139111 ],

src/spaces/srgb-linear.js

@@ -1,12 +1,16 @@
 import RGBColorSpace from "../RGBColorSpace.js";
 
+// Type "imports"
+/** @typedef {import("../types.js").Matrix3x3} Matrix3x3 */
+
 // This is the linear-light version of sRGB
 // as used for example in SVG filters
 // or in Canvas
 
 // This matrix was calculated directly from the RGB and white chromaticities
 // when rounded to 8 decimal places, it agrees completely with the official matrix
 // see https://github.com/w3c/csswg-drafts/issues/5922
+/** @type {Matrix3x3} */
 const toXYZ_M = [
 	[ 0.41239079926595934, 0.357584339383878,   0.1804807884018343  ],
 	[ 0.21263900587151027, 0.715168678767756,   0.07219231536073371 ],
@@ -15,6 +19,7 @@ const toXYZ_M = [
 
 // This matrix is the inverse of the above;
 // again it agrees with the official definition when rounded to 8 decimal places
+/** @type {Matrix3x3} */
 export const fromXYZ_M = [
 	[  3.2409699419045226,  -1.537383177570094,   -0.4986107602930034  ],
 	[ -0.9692436362808796,   1.8759675015077202,   0.04155505740717559 ],

src/types.d.ts

@@ -2,6 +2,11 @@
  * @packageDocumentation
  * Defines and re-exports many types for use throughout the library.
  */
+
+// muliply-matricies.js
+export type Matrix3x3 = [[number, number, number], [number, number, number], [number, number, number]];
+export type Vector3 = [number, number, number];
+
 // contrast/
 export type * from "./contrast/index.js";
 
@@ -104,8 +109,8 @@ export interface ParseFunctionReturn {
 
 // rgbspace.js
 export interface RGBOptions extends SpaceOptions {
-	toXYZ_M?: number[][] | undefined;
-	fromXYZ_M?: number[][] | undefined;
+	toXYZ_M?: Matrix3x3 | undefined;
+	fromXYZ_M?: Matrix3x3 | undefined;
 }
 
 // serialize.js

src/util.js

@@ -2,7 +2,7 @@
  * Various utility functions
  */
 
-export {default as multiplyMatrices} from "./multiply-matrices.js";
+export {default as multiplyMatrices, multiply_v3_m3x3} from "./multiply-matrices.js";
 
 /**
  * Check if a value is a string (including a String object)

test/multiply_matrices.js

@@ -1,5 +1,5 @@
 import * as math from "mathjs"; // Used as test oracle
-import multiplyMatrices from "../src/multiply-matrices.js";
+import {multiplyMatrices, multiply_v3_m3x3} from "../src/util.js";
 import * as check from "../node_modules/htest.dev/src/check.js";
 
 // Used to collect expected results from oracle
@@ -93,5 +93,25 @@ export default {
 				args: [[], []],
 			},
 		].map(expectThrows),
+	},
+	{
+		name: "Transform",
+		run: multiply_v3_m3x3,
+		tests: [
+			{
+				name: "3x3 matrix with vector",
+				args: [[1, .5, 0], M_lin_sRGB_to_XYZ],
+				expect: math.multiply(math.matrix(M_lin_sRGB_to_XYZ), math.matrix([1, .5, 0])).valueOf(),
+			},
+			{
+				name: "3x3 matrix with vector in place",
+				run (A, B) {
+					multiply_v3_m3x3(A, B, A);
+					return A;
+				},
+				args: [[1, .5, 0], M_lin_sRGB_to_XYZ],
+				expect: math.multiply(math.matrix(M_lin_sRGB_to_XYZ), math.matrix([1, .5, 0])).valueOf(),
+			},
+		],
 	}],
 };

types/test/RGBColorSpace.ts

@@ -13,6 +13,6 @@ new RGBColorSpace({
 new RGBColorSpace({
 	name: "RGBSpace",
 	id: "rgbspace",
-	toXYZ_M: [[1, 2, 3]],
-	fromXYZ_M: [[3, 2, 1]],
+	toXYZ_M: [[1, 2, 3], [4, 5, 6], [7, 8, 9]],
+	fromXYZ_M: [[3, 2, 1], [6, 5, 4], [9, 8, 7]],
 });