Math

Pi() : Math.PI;

Return PI.

radiansFromDegrees(degree) : radian

Convert degrees to radians.

degreesFromRadians(radian) : degree

Convert radians to degrees.

round(float): int

Same as Math.round().

floor(float) : int

Same as Math.floor().

ceil(float) : int

Same as Math.ceil().

ceilLog2()

NOT IMPLEMENTED

min(a, b)

Same as Math.min().

max(a, b)

Same as Math.max().

arrayMin(array)

Minimum of the array.

arrayMax(a, b)

Maximum of the array.

factorial(n) : number

NOT IMPLEMENTED

binomial(m, n) : int

beginCombination(m, n) : Array or null

nextCombination(m, n, r) : Boolean

randomSeed()

NOT IMPLEMENTED

getSeed()

NOT IMPLEMENTED

random(minValue = 0, maxValue = 1) : Number

gaussian()

NOT IMPLEMENTED

add(a, b, out)

a[3] + b[3] => out[3]

subtract(a, b, out)

a[3] - b[3] => out[3]

multiplyScalar(vec, scalar) {

vec[3] * scalar => vec[3]

multiplyScalar2D(vec, scalar)

vec[2] * scalar => vec[2]

dot(x, y)

outer(x, y, out_3x3)

cross(x, y, out)

norm(x, n = 3)

normalize(x)

Normalize in place

perpendiculars(x, y, z, theta)

projectVector(a, b, projection)

dot2D(x, y)

projectVector2D(a, b, projection)

distance2BetweenPoints(x, y)

angleBetweenVectors(v1, v2)

gaussianAmplitude(mean, variance, position)

gaussianWeight(mean, variance, position)

outer2D(x, y, out_2x2)

norm2D(x2D)

normalize2D(x)

determinant2x2(c1, c2)

LUFactor3x3(mat_3x3, index_3)

LUSolve3x3(mat_3x3, index_3, x_3)

linearSolve3x3(mat_3x3, x_3, y_3)

multiply3x3_vect3(mat_3x3, in_3, out_3)

multiply3x3_mat3(a_3x3, b_3x3, out_3x3)

multiplyMatrix(a, b, rowA, colA, rowB, colB, out_rowXcol)

transpose3x3(in_3x3, outT_3x3)

invert3x3(in_3x3, outI_3x3)

identity3x3(mat_3x3)

determinant3x3(mat_3x3)

quaternionToMatrix3x3(quat_4, mat_3x3)

jacobiN(a, n, w, v)

matrix3x3ToQuaternion(mat_3x3, quat_4)

multiplyQuaternion(quat_1, quat_2, quat_out)

orthogonalize3x3(a_3x3, out_3x3)

diagonalize3x3(a_3x3, w_3, v_3x3)

singularValueDecomposition3x3(a_3x3, u_3x3, w_3, vT_3x3)

luFactorLinearSystem(A, index, size)

luSolveLinearSystem(A, index, x, size)

solveLinearSystem(A, x, size)

invertMatrix(A, AI, size, index = null, column = null)

estimateMatrixCondition(A, size)

jacobi(a_3x3, w, v)

jacobiN(a_3x3, 3, w, v);

solveHomogeneousLeastSquares(numberOfSamples, xt, xOrder, mt)

solveLeastSquares(numberOfSamples, xt, xOrder, yt, yOrder, mt, checkHomogeneous = true)

rgb2hsv(rgb, hsv)

hsv2rgb(hsv, rgb)

lab2xyz(lab, xyz)

xyz2lab(xyz, lab)

xyz2rgb(xyz, rgb)

rgb2xyz(rgb, xyz)

rgb2lab(rgb, lab)

LabToRGB(lab, rgb)

uninitializeBounds(bounds)

areBoundsInitialized(bounds)

clampValue(value, minValue, maxValue)

clampAndNormalizeValue(value, range)

getScalarTypeFittingRange

NOT IMPLEMENTED

getAdjustedScalarRange

NOT IMPLEMENTED

extentIsWithinOtherExtent(extent1, extent2)

boundsIsWithinOtherBounds(bounds1_6, bounds2_6, delta_3)

pointIsWithinBounds(point_3, bounds_6, delta_3)

solve3PointCircle(p1, p2, p3, center)

inf()

negInf()

isInf(number) : Boolean

isNan(?) : Boolean

isFinite(number) : Boolean

createUninitializedBouds() : new [6]

Source

index.js
import seedrandom from 'seedrandom';
import macro from 'vtk.js/Sources/macro';

const { vtkErrorMacro, vtkWarningMacro } = macro;

// ----------------------------------------------------------------------------
/* eslint-disable camelcase */
/* eslint-disable no-cond-assign */
/* eslint-disable no-bitwise */
/* eslint-disable no-multi-assign */
// ----------------------------------------------------------------------------
let randomSeedValue = 0;
const VTK_MAX_ROTATIONS = 20;
const VTK_SMALL_NUMBER = 1.0e-12;

function notImplemented(method) {
return () => vtkErrorMacro(`vtkMath::${method} - NOT IMPLEMENTED`);
}

function vtkSwapVectors3(v1, v2) {
for (let i = 0; i < 3; i++) {
const tmp = v1[i];
v1[i] = v2[i];
v2[i] = tmp;
}
}

function createArray(size = 3) {
const array = [];
while (array.length < size) {
array.push(0);
}
return array;
}

// ----------------------------------------------------------------------------
// Global methods
// ----------------------------------------------------------------------------

export const Pi = () => Math.PI;

export function radiansFromDegrees(deg) {
return (deg / 180) * Math.PI;
}

export function degreesFromRadians(rad) {
return (rad * 180) / Math.PI;
}

export const { round, floor, ceil, min, max } = Math;

export function arrayMin(arr, offset = 0, stride = 1) {
let minValue = Infinity;
for (let i = offset, len = arr.length; i < len; i += stride) {
if (arr[i] < minValue) {
minValue = arr[i];
}
}

return minValue;
}

export function arrayMax(arr, offset = 0, stride = 1) {
let maxValue = -Infinity;
for (let i = offset, len = arr.length; i < len; i += stride) {
if (maxValue < arr[i]) {
maxValue = arr[i];
}
}

return maxValue;
}

export function arrayRange(arr, offset = 0, stride = 1) {
let minValue = Infinity;
let maxValue = -Infinity;
for (let i = offset, len = arr.length; i < len; i += stride) {
if (arr[i] < minValue) {
minValue = arr[i];
}
if (maxValue < arr[i]) {
maxValue = arr[i];
}
}

return [minValue, maxValue];
}

export const ceilLog2 = notImplemented('ceilLog2');
export const factorial = notImplemented('factorial');

export function nearestPowerOfTwo(xi) {
let v = 1;
while (v < xi) {
v *= 2;
}
return v;
}

export function isPowerOfTwo(x) {
return x === nearestPowerOfTwo(x);
}

export function binomial(m, n) {
let r = 1;
for (let i = 1; i <= n; ++i) {
r *= (m - i + 1) / i;
}
return Math.floor(r);
}

export function beginCombination(m, n) {
if (m < n) {
return 0;
}

const r = createArray(n);
for (let i = 0; i < n; ++i) {
r[i] = i;
}
return r;
}

export function nextCombination(m, n, r) {
let status = 0;
for (let i = n - 1; i >= 0; --i) {
if (r[i] < m - n + i) {
let j = r[i] + 1;
while (i < n) {
r[i++] = j++;
}
status = 1;
break;
}
}
return status;
}

export function randomSeed(seed) {
seedrandom(`${seed}`, { global: true });
randomSeedValue = seed;
}

export function getSeed() {
return randomSeedValue;
}

export function random(minValue = 0, maxValue = 1) {
const delta = maxValue - minValue;
return minValue + delta * Math.random();
}

export const gaussian = notImplemented('gaussian');

// Vect3 operations
export function add(a, b, out) {
out[0] = a[0] + b[0];
out[1] = a[1] + b[1];
out[2] = a[2] + b[2];
}

export function subtract(a, b, out) {
out[0] = a[0] - b[0];
out[1] = a[1] - b[1];
out[2] = a[2] - b[2];
}

export function multiplyScalar(vec, scalar) {
vec[0] *= scalar;
vec[1] *= scalar;
vec[2] *= scalar;
}

export function multiplyScalar2D(vec, scalar) {
vec[0] *= scalar;
vec[1] *= scalar;
}

export function dot(x, y) {
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
}

export function outer(x, y, out_3x3) {
for (let i = 0; i < 3; i++) {
for (let j = 0; j < 3; j++) {
out_3x3[i][j] = x[i] * y[j];
}
}
}

export function cross(x, y, out) {
const Zx = x[1] * y[2] - x[2] * y[1];
const Zy = x[2] * y[0] - x[0] * y[2];
const Zz = x[0] * y[1] - x[1] * y[0];
out[0] = Zx;
out[1] = Zy;
out[2] = Zz;
}

export function norm(x, n = 3) {
switch (n) {
case 1:
return Math.abs(x);
case 2:
return Math.sqrt(x[0] * x[0] + x[1] * x[1]);
case 3:
return Math.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
default: {
let sum = 0;
for (let i = 0; i < n; i++) {
sum += x[i] * x[i];
}
return Math.sqrt(sum);
}
}
}

export function normalize(x) {
const den = norm(x);
if (den !== 0.0) {
x[0] /= den;
x[1] /= den;
x[2] /= den;
}
return den;
}

export function perpendiculars(x, y, z, theta) {
const x2 = x[0] * x[0];
const y2 = x[1] * x[1];
const z2 = x[2] * x[2];
const r = Math.sqrt(x2 + y2 + z2);

let dx;
let dy;
let dz;

// transpose the vector to avoid divide-by-zero error
if (x2 > y2 && x2 > z2) {
dx = 0;
dy = 1;
dz = 2;
} else if (y2 > z2) {
dx = 1;
dy = 2;
dz = 0;
} else {
dx = 2;
dy = 0;
dz = 1;
}

const a = x[dx] / r;
const b = x[dy] / r;
const c = x[dz] / r;
const tmp = Math.sqrt(a * a + c * c);

if (theta !== 0) {
const sintheta = Math.sin(theta);
const costheta = Math.cos(theta);

if (y) {
y[dx] = (c * costheta - a * b * sintheta) / tmp;
y[dy] = sintheta * tmp;
y[dz] = (-(a * costheta) - b * c * sintheta) / tmp;
}

if (z) {
z[dx] = (-(c * sintheta) - a * b * costheta) / tmp;
z[dy] = costheta * tmp;
z[dz] = (a * sintheta - b * c * costheta) / tmp;
}
} else {
if (y) {
y[dx] = c / tmp;
y[dy] = 0;
y[dz] = -a / tmp;
}

if (z) {
z[dx] = (-a * b) / tmp;
z[dy] = tmp;
z[dz] = (-b * c) / tmp;
}
}
}

export function projectVector(a, b, projection) {
const bSquared = dot(b, b);

if (bSquared === 0) {
projection[0] = 0;
projection[1] = 0;
projection[2] = 0;
return false;
}

const scale = dot(a, b) / bSquared;

for (let i = 0; i < 3; i++) {
projection[i] = b[i];
}
multiplyScalar(projection, scale);

return true;
}

export function dot2D(x, y) {
return x[0] * y[0] + x[1] * y[1];
}

export function projectVector2D(a, b, projection) {
const bSquared = dot2D(b, b);

if (bSquared === 0) {
projection[0] = 0;
projection[1] = 0;
return false;
}

const scale = dot2D(a, b) / bSquared;

for (let i = 0; i < 2; i++) {
projection[i] = b[i];
}
multiplyScalar2D(projection, scale);

return true;
}

export function distance2BetweenPoints(x, y) {
return (
(x[0] - y[0]) * (x[0] - y[0]) +
(x[1] - y[1]) * (x[1] - y[1]) +
(x[2] - y[2]) * (x[2] - y[2])
);
}

export function angleBetweenVectors(v1, v2) {
const crossVect = [0, 0, 0];
cross(v1, v2, crossVect);
return Math.atan2(norm(crossVect), dot(v1, v2));
}

export function gaussianAmplitude(mean, variance, position) {
const distanceFromMean = Math.abs(mean - position);
return (
(1 / Math.sqrt(2 * Math.PI * variance)) *
Math.exp(-(distanceFromMean ** 2) / (2 * variance))
);
}

export function gaussianWeight(mean, variance, position) {
const distanceFromMean = Math.abs(mean - position);
return Math.exp(-(distanceFromMean ** 2) / (2 * variance));
}

export function outer2D(x, y, out_2x2) {
for (let i = 0; i < 2; i++) {
for (let j = 0; j < 2; j++) {
out_2x2[i][j] = x[i] * y[j];
}
}
}

export function norm2D(x2D) {
return Math.sqrt(x2D[0] * x2D[0] + x2D[1] * x2D[1]);
}

export function normalize2D(x) {
const den = norm2D(x);
if (den !== 0.0) {
x[0] /= den;
x[1] /= den;
}
return den;
}

export function determinant2x2(...args) {
if (args.length === 2) {
return args[0][0] * args[1][1] - args[1][0] * args[0][1];
}
if (args.length === 4) {
return args[0] * args[3] - args[1] * args[2];
}
return Number.NaN;
}

export function LUFactor3x3(mat_3x3, index_3) {
let maxI;
let tmp;
let largest;
const scale = [0, 0, 0];

// Loop over rows to get implicit scaling information
for (let i = 0; i < 3; i++) {
largest = Math.abs(mat_3x3[i][0]);
if ((tmp = Math.abs(mat_3x3[i][1])) > largest) {
largest = tmp;
}
if ((tmp = Math.abs(mat_3x3[i][2])) > largest) {
largest = tmp;
}
scale[i] = 1 / largest;
}

// Loop over all columns using Crout's method

// first column
largest = scale[0] * Math.abs(mat_3x3[0][0]);
maxI = 0;
if ((tmp = scale[1] * Math.abs(mat_3x3[1][0])) >= largest) {
largest = tmp;
maxI = 1;
}
if ((tmp = scale[2] * Math.abs(mat_3x3[2][0])) >= largest) {
maxI = 2;
}
if (maxI !== 0) {
vtkSwapVectors3(mat_3x3[maxI], mat_3x3[0]);
scale[maxI] = scale[0];
}
index_3[0] = maxI;

mat_3x3[1][0] /= mat_3x3[0][0];
mat_3x3[2][0] /= mat_3x3[0][0];

// second column
mat_3x3[1][1] -= mat_3x3[1][0] * mat_3x3[0][1];
mat_3x3[2][1] -= mat_3x3[2][0] * mat_3x3[0][1];
largest = scale[1] * Math.abs(mat_3x3[1][1]);
maxI = 1;
if ((tmp = scale[2] * Math.abs(mat_3x3[2][1])) >= largest) {
maxI = 2;
vtkSwapVectors3(mat_3x3[2], mat_3x3[1]);
scale[2] = scale[1];
}
index_3[1] = maxI;
mat_3x3[2][1] /= mat_3x3[1][1];

// third column
mat_3x3[1][2] -= mat_3x3[1][0] * mat_3x3[0][2];
mat_3x3[2][2] -=
mat_3x3[2][0] * mat_3x3[0][2] + mat_3x3[2][1] * mat_3x3[1][2];
index_3[2] = 2;
}

export function LUSolve3x3(mat_3x3, index_3, x_3) {
// forward substitution
let sum = x_3[index_3[0]];
x_3[index_3[0]] = x_3[0];
x_3[0] = sum;

sum = x_3[index_3[1]];
x_3[index_3[1]] = x_3[1];
x_3[1] = sum - mat_3x3[1][0] * x_3[0];

sum = x_3[index_3[2]];
x_3[index_3[2]] = x_3[2];
x_3[2] = sum - mat_3x3[2][0] * x_3[0] - mat_3x3[2][1] * x_3[1];

// back substitution
x_3[2] /= mat_3x3[2][2];
x_3[1] = (x_3[1] - mat_3x3[1][2] * x_3[2]) / mat_3x3[1][1];
x_3[0] =
(x_3[0] - mat_3x3[0][1] * x_3[1] - mat_3x3[0][2] * x_3[2]) / mat_3x3[0][0];
}

export function linearSolve3x3(mat_3x3, x_3, y_3) {
const a1 = mat_3x3[0][0];
const b1 = mat_3x3[0][1];
const c1 = mat_3x3[0][2];
const a2 = mat_3x3[1][0];
const b2 = mat_3x3[1][1];
const c2 = mat_3x3[1][2];
const a3 = mat_3x3[2][0];
const b3 = mat_3x3[2][1];
const c3 = mat_3x3[2][2];

// Compute the adjoint
const d1 = +determinant2x2(b2, b3, c2, c3);
const d2 = -determinant2x2(a2, a3, c2, c3);
const d3 = +determinant2x2(a2, a3, b2, b3);

const e1 = -determinant2x2(b1, b3, c1, c3);
const e2 = +determinant2x2(a1, a3, c1, c3);
const e3 = -determinant2x2(a1, a3, b1, b3);

const f1 = +determinant2x2(b1, b2, c1, c2);
const f2 = -determinant2x2(a1, a2, c1, c2);
const f3 = +determinant2x2(a1, a2, b1, b2);

// Compute the determinant
const det = a1 * d1 + b1 * d2 + c1 * d3;

// Multiply by the adjoint
const v1 = d1 * x_3[0] + e1 * x_3[1] + f1 * x_3[2];
const v2 = d2 * x_3[0] + e2 * x_3[1] + f2 * x_3[2];
const v3 = d3 * x_3[0] + e3 * x_3[1] + f3 * x_3[2];

// Divide by the determinant
y_3[0] = v1 / det;
y_3[1] = v2 / det;
y_3[2] = v3 / det;
}

export function multiply3x3_vect3(mat_3x3, in_3, out_3) {
const x =
mat_3x3[0][0] * in_3[0] + mat_3x3[0][1] * in_3[1] + mat_3x3[0][2] * in_3[2];
const y =
mat_3x3[1][0] * in_3[0] + mat_3x3[1][1] * in_3[1] + mat_3x3[1][2] * in_3[2];
const z =
mat_3x3[2][0] * in_3[0] + mat_3x3[2][1] * in_3[1] + mat_3x3[2][2] * in_3[2];

out_3[0] = x;
out_3[1] = y;
out_3[2] = z;
}

export function multiply3x3_mat3(a_3x3, b_3x3, out_3x3) {
const tmp = [[0, 0, 0], [0, 0, 0], [0, 0, 0]];

for (let i = 0; i < 3; i++) {
tmp[0][i] =
a_3x3[0][0] * b_3x3[0][i] +
a_3x3[0][1] * b_3x3[1][i] +
a_3x3[0][2] * b_3x3[2][i];
tmp[1][i] =
a_3x3[1][0] * b_3x3[0][i] +
a_3x3[1][1] * b_3x3[1][i] +
a_3x3[1][2] * b_3x3[2][i];
tmp[2][i] =
a_3x3[2][0] * b_3x3[0][i] +
a_3x3[2][1] * b_3x3[1][i] +
a_3x3[2][2] * b_3x3[2][i];
}

for (let j = 0; j < 3; j++) {
out_3x3[j][0] = tmp[j][0];
out_3x3[j][1] = tmp[j][1];
out_3x3[j][2] = tmp[j][2];
}
}

export function multiplyMatrix(a, b, rowA, colA, rowB, colB, out_rowXcol) {
// we need colA == rowB
if (colA !== rowB) {
vtkErrorMacro('Number of columns of A must match number of rows of B.');
}

// output matrix is rowA*colB
// output row
for (let i = 0; i < rowA; i++) {
// output col
for (let j = 0; j < colB; j++) {
out_rowXcol[i][j] = 0;
// sum for this point
for (let k = 0; k < colA; k++) {
out_rowXcol[i][j] += a[i][k] * b[k][j];
}
}
}
}

export function transpose3x3(in_3x3, outT_3x3) {
let tmp;
tmp = in_3x3[1][0];
outT_3x3[1][0] = in_3x3[0][1];
outT_3x3[0][1] = tmp;
tmp = in_3x3[2][0];
outT_3x3[2][0] = in_3x3[0][2];
outT_3x3[0][2] = tmp;
tmp = in_3x3[2][1];
outT_3x3[2][1] = in_3x3[1][2];
outT_3x3[1][2] = tmp;

outT_3x3[0][0] = in_3x3[0][0];
outT_3x3[1][1] = in_3x3[1][1];
outT_3x3[2][2] = in_3x3[2][2];
}

export function invert3x3(in_3x3, outI_3x3) {
const a1 = in_3x3[0][0];
const b1 = in_3x3[0][1];
const c1 = in_3x3[0][2];
const a2 = in_3x3[1][0];
const b2 = in_3x3[1][1];
const c2 = in_3x3[1][2];
const a3 = in_3x3[2][0];
const b3 = in_3x3[2][1];
const c3 = in_3x3[2][2];

// Compute the adjoint
const d1 = +determinant2x2(b2, b3, c2, c3);
const d2 = -determinant2x2(a2, a3, c2, c3);
const d3 = +determinant2x2(a2, a3, b2, b3);

const e1 = -determinant2x2(b1, b3, c1, c3);
const e2 = +determinant2x2(a1, a3, c1, c3);
const e3 = -determinant2x2(a1, a3, b1, b3);

const f1 = +determinant2x2(b1, b2, c1, c2);
const f2 = -determinant2x2(a1, a2, c1, c2);
const f3 = +determinant2x2(a1, a2, b1, b2);

// Divide by the determinant
const det = a1 * d1 + b1 * d2 + c1 * d3;

outI_3x3[0][0] = d1 / det;
outI_3x3[1][0] = d2 / det;
outI_3x3[2][0] = d3 / det;

outI_3x3[0][1] = e1 / det;
outI_3x3[1][1] = e2 / det;
outI_3x3[2][1] = e3 / det;

outI_3x3[0][2] = f1 / det;
outI_3x3[1][2] = f2 / det;
outI_3x3[2][2] = f3 / det;
}

export function identity3x3(mat_3x3) {
for (let i = 0; i < 3; i++) {
mat_3x3[i][0] = mat_3x3[i][1] = mat_3x3[i][2] = 0;
mat_3x3[i][i] = 1;
}
}

export function determinant3x3(mat_3x3) {
return (
mat_3x3[0][0] * mat_3x3[1][1] * mat_3x3[2][2] +
mat_3x3[1][0] * mat_3x3[2][1] * mat_3x3[0][2] +
mat_3x3[2][0] * mat_3x3[0][1] * mat_3x3[1][2] -
mat_3x3[0][0] * mat_3x3[2][1] * mat_3x3[1][2] -
mat_3x3[1][0] * mat_3x3[0][1] * mat_3x3[2][2] -
mat_3x3[2][0] * mat_3x3[1][1] * mat_3x3[0][2]
);
}

export function quaternionToMatrix3x3(quat_4, mat_3x3) {
const ww = quat_4[0] * quat_4[0];
const wx = quat_4[0] * quat_4[1];
const wy = quat_4[0] * quat_4[2];
const wz = quat_4[0] * quat_4[3];

const xx = quat_4[1] * quat_4[1];
const yy = quat_4[2] * quat_4[2];
const zz = quat_4[3] * quat_4[3];

const xy = quat_4[1] * quat_4[2];
const xz = quat_4[1] * quat_4[3];
const yz = quat_4[2] * quat_4[3];

const rr = xx + yy + zz;
// normalization factor, just in case quaternion was not normalized
let f = 1 / (ww + rr);
const s = (ww - rr) * f;
f *= 2;

mat_3x3[0][0] = xx * f + s;
mat_3x3[1][0] = (xy + wz) * f;
mat_3x3[2][0] = (xz - wy) * f;

mat_3x3[0][1] = (xy - wz) * f;
mat_3x3[1][1] = yy * f + s;
mat_3x3[2][1] = (yz + wx) * f;

mat_3x3[0][2] = (xz + wy) * f;
mat_3x3[1][2] = (yz - wx) * f;
mat_3x3[2][2] = zz * f + s;
}

export function areMatricesEqual(matA, matB) {
if (!matA.length === matB.length) {
return false;
}

function isEqual(element, index) {
return element === matB[index];
}
return matA.every(isEqual);
}

export function jacobiN(a, n, w, v) {
let i;
let j;
let k;
let iq;
let ip;
let numPos;
let tresh;
let theta;
let t;
let tau;
let sm;
let s;
let h;
let g;
let c;
let tmp;
const b = createArray(n);
const z = createArray(n);

const vtkROTATE = (aa, ii, jj, kk, ll) => {
g = aa[ii][jj];
h = aa[kk][ll];
aa[ii][jj] = g - s * (h + g * tau);
aa[kk][ll] = h + s * (g - h * tau);
};

// initialize
for (ip = 0; ip < n; ip++) {
for (iq = 0; iq < n; iq++) {
v[ip][iq] = 0.0;
}
v[ip][ip] = 1.0;
}
for (ip = 0; ip < n; ip++) {
b[ip] = w[ip] = a[ip][ip];
z[ip] = 0.0;
}

// begin rotation sequence
for (i = 0; i < VTK_MAX_ROTATIONS; i++) {
sm = 0.0;
for (ip = 0; ip < n - 1; ip++) {
for (iq = ip + 1; iq < n; iq++) {
sm += Math.abs(a[ip][iq]);
}
}
if (sm === 0.0) {
break;
}

// first 3 sweeps
if (i < 3) {
tresh = (0.2 * sm) / (n * n);
} else {
tresh = 0.0;
}

for (ip = 0; ip < n - 1; ip++) {
for (iq = ip + 1; iq < n; iq++) {
g = 100.0 * Math.abs(a[ip][iq]);

// after 4 sweeps
if (
i > 3 &&
Math.abs(w[ip]) + g === Math.abs(w[ip]) &&
Math.abs(w[iq]) + g === Math.abs(w[iq])
) {
a[ip][iq] = 0.0;
} else if (Math.abs(a[ip][iq]) > tresh) {
h = w[iq] - w[ip];
if (Math.abs(h) + g === Math.abs(h)) {
t = a[ip][iq] / h;
} else {
theta = (0.5 * h) / a[ip][iq];
t = 1.0 / (Math.abs(theta) + Math.sqrt(1.0 + theta * theta));
if (theta < 0.0) {
t = -t;
}
}
c = 1.0 / Math.sqrt(1 + t * t);
s = t * c;
tau = s / (1.0 + c);
h = t * a[ip][iq];
z[ip] -= h;
z[iq] += h;
w[ip] -= h;
w[iq] += h;
a[ip][iq] = 0.0;

// ip already shifted left by 1 unit
for (j = 0; j <= ip - 1; j++) {
vtkROTATE(a, j, ip, j, iq);
}
// ip and iq already shifted left by 1 unit
for (j = ip + 1; j <= iq - 1; j++) {
vtkROTATE(a, ip, j, j, iq);
}
// iq already shifted left by 1 unit
for (j = iq + 1; j < n; j++) {
vtkROTATE(a, ip, j, iq, j);
}
for (j = 0; j < n; j++) {
vtkROTATE(v, j, ip, j, iq);
}
}
}
}

for (ip = 0; ip < n; ip++) {
b[ip] += z[ip];
w[ip] = b[ip];
z[ip] = 0.0;
}
}

// this is NEVER called
if (i >= VTK_MAX_ROTATIONS) {
vtkWarningMacro('vtkMath::Jacobi: Error extracting eigenfunctions');
return 0;
}

// sort eigenfunctions: these changes do not affect accuracy
for (j = 0; j < n - 1; j++) {
// boundary incorrect
k = j;
tmp = w[k];
for (i = j + 1; i < n; i++) {
// boundary incorrect, shifted already
if (w[i] >= tmp) {
// why exchange if same?
k = i;
tmp = w[k];
}
}
if (k !== j) {
w[k] = w[j];
w[j] = tmp;
for (i = 0; i < n; i++) {
tmp = v[i][j];
v[i][j] = v[i][k];
v[i][k] = tmp;
}
}
}
// ensure eigenvector consistency (i.e., Jacobi can compute vectors that
// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
// reek havoc in hyperstreamline/other stuff. We will select the most
// positive eigenvector.
const ceil_half_n = (n >> 1) + (n & 1);
for (j = 0; j < n; j++) {
for (numPos = 0, i = 0; i < n; i++) {
if (v[i][j] >= 0.0) {
numPos++;
}
}
// if ( numPos < ceil(double(n)/double(2.0)) )
if (numPos < ceil_half_n) {
for (i = 0; i < n; i++) {
v[i][j] *= -1.0;
}
}
}
return 1;
}

export function matrix3x3ToQuaternion(mat_3x3, quat_4) {
const tmp = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]];

// on-diagonal elements
tmp[0][0] = mat_3x3[0][0] + mat_3x3[1][1] + mat_3x3[2][2];
tmp[1][1] = mat_3x3[0][0] - mat_3x3[1][1] - mat_3x3[2][2];
tmp[2][2] = -mat_3x3[0][0] + mat_3x3[1][1] - mat_3x3[2][2];
tmp[3][3] = -mat_3x3[0][0] - mat_3x3[1][1] + mat_3x3[2][2];

// off-diagonal elements
tmp[0][1] = tmp[1][0] = mat_3x3[2][1] - mat_3x3[1][2];
tmp[0][2] = tmp[2][0] = mat_3x3[0][2] - mat_3x3[2][0];
tmp[0][3] = tmp[3][0] = mat_3x3[1][0] - mat_3x3[0][1];

tmp[1][2] = tmp[2][1] = mat_3x3[1][0] + mat_3x3[0][1];
tmp[1][3] = tmp[3][1] = mat_3x3[0][2] + mat_3x3[2][0];
tmp[2][3] = tmp[3][2] = mat_3x3[2][1] + mat_3x3[1][2];

const eigenvectors = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]];
const eigenvalues = [0, 0, 0, 0];

// convert into format that JacobiN can use,
// then use Jacobi to find eigenvalues and eigenvectors
const NTemp = [0, 0, 0, 0];
const eigenvectorsTemp = [0, 0, 0, 0];
for (let i = 0; i < 4; i++) {
NTemp[i] = tmp[i];
eigenvectorsTemp[i] = eigenvectors[i];
}
jacobiN(NTemp, 4, eigenvalues, eigenvectorsTemp);

// the first eigenvector is the one we want
quat_4[0] = eigenvectors[0][0];
quat_4[1] = eigenvectors[1][0];
quat_4[2] = eigenvectors[2][0];
quat_4[3] = eigenvectors[3][0];
}

export function multiplyQuaternion(quat_1, quat_2, quat_out) {
const ww = quat_1[0] * quat_2[0];
const wx = quat_1[0] * quat_2[1];
const wy = quat_1[0] * quat_2[2];
const wz = quat_1[0] * quat_2[3];

const xw = quat_1[1] * quat_2[0];
const xx = quat_1[1] * quat_2[1];
const xy = quat_1[1] * quat_2[2];
const xz = quat_1[1] * quat_2[3];

const yw = quat_1[2] * quat_2[0];
const yx = quat_1[2] * quat_2[1];
const yy = quat_1[2] * quat_2[2];
const yz = quat_1[2] * quat_2[3];

const zw = quat_1[3] * quat_2[0];
const zx = quat_1[3] * quat_2[1];
const zy = quat_1[3] * quat_2[2];
const zz = quat_1[3] * quat_2[3];

quat_out[0] = ww - xx - yy - zz;
quat_out[1] = wx + xw + yz - zy;
quat_out[2] = wy - xz + yw + zx;
quat_out[3] = wz + xy - yx + zw;
}

export function orthogonalize3x3(a_3x3, out_3x3) {
// copy the matrix
for (let i = 0; i < 3; i++) {
out_3x3[0][i] = a_3x3[0][i];
out_3x3[1][i] = a_3x3[1][i];
out_3x3[2][i] = a_3x3[2][i];
}

// Pivot the matrix to improve accuracy
const scale = createArray(3);
const index = createArray(3);
let largest;

// Loop over rows to get implicit scaling information
for (let i = 0; i < 3; i++) {
const x1 = Math.abs(out_3x3[i][0]);
const x2 = Math.abs(out_3x3[i][1]);
const x3 = Math.abs(out_3x3[i][2]);
largest = x2 > x1 ? x2 : x1;
largest = x3 > largest ? x3 : largest;
scale[i] = 1;
if (largest !== 0) {
scale[i] /= largest;
}
}

// first column
const x1 = Math.abs(out_3x3[0][0]) * scale[0];
const x2 = Math.abs(out_3x3[1][0]) * scale[1];
const x3 = Math.abs(out_3x3[2][0]) * scale[2];
index[0] = 0;
largest = x1;
if (x2 >= largest) {
largest = x2;
index[0] = 1;
}
if (x3 >= largest) {
index[0] = 2;
}
if (index[0] !== 0) {
vtkSwapVectors3(out_3x3[index[0]], out_3x3[0]);
scale[index[0]] = scale[0];
}

// second column
const y2 = Math.abs(out_3x3[1][1]) * scale[1];
const y3 = Math.abs(out_3x3[2][1]) * scale[2];
index[1] = 1;
largest = y2;
if (y3 >= largest) {
index[1] = 2;
vtkSwapVectors3(out_3x3[2], out_3x3[1]);
}

// third column
index[2] = 2;

// A quaternion can only describe a pure rotation, not
// a rotation with a flip, therefore the flip must be
// removed before the matrix is converted to a quaternion.
let flip = 0;
if (determinant3x3(out_3x3) < 0) {
flip = 1;
for (let i = 0; i < 3; i++) {
out_3x3[0][i] = -out_3x3[0][i];
out_3x3[1][i] = -out_3x3[1][i];
out_3x3[2][i] = -out_3x3[2][i];
}
}

// Do orthogonalization using a quaternion intermediate
// (this, essentially, does the orthogonalization via
// diagonalization of an appropriately constructed symmetric
// 4x4 matrix rather than by doing SVD of the 3x3 matrix)
const quat = createArray(4);
matrix3x3ToQuaternion(out_3x3, quat);
quaternionToMatrix3x3(quat, out_3x3);

// Put the flip back into the orthogonalized matrix.
if (flip) {
for (let i = 0; i < 3; i++) {
out_3x3[0][i] = -out_3x3[0][i];
out_3x3[1][i] = -out_3x3[1][i];
out_3x3[2][i] = -out_3x3[2][i];
}
}

// Undo the pivoting
if (index[1] !== 1) {
vtkSwapVectors3(out_3x3[index[1]], out_3x3[1]);
}
if (index[0] !== 0) {
vtkSwapVectors3(out_3x3[index[0]], out_3x3[0]);
}
}

export function diagonalize3x3(a_3x3, w_3, v_3x3) {
let i;
let j;
let k;
let maxI;
let tmp;
let maxVal;

// do the matrix[3][3] to **matrix conversion for Jacobi
const C = [createArray(3), createArray(3), createArray(3)];
const ATemp = createArray(3);
const VTemp = createArray(3);
for (i = 0; i < 3; i++) {
C[i][0] = a_3x3[i][0];
C[i][1] = a_3x3[i][1];
C[i][2] = a_3x3[i][2];
ATemp[i] = C[i];
VTemp[i] = v_3x3[i];
}

// diagonalize using Jacobi
jacobiN(ATemp, 3, w_3, VTemp);

// if all the eigenvalues are the same, return identity matrix
if (w_3[0] === w_3[1] && w_3[0] === w_3[2]) {
identity3x3(v_3x3);
return;
}

// transpose temporarily, it makes it easier to sort the eigenvectors
transpose3x3(v_3x3, v_3x3);

// if two eigenvalues are the same, re-orthogonalize to optimally line
// up the eigenvectors with the x, y, and z axes
for (i = 0; i < 3; i++) {
// two eigenvalues are the same
if (w_3[(i + 1) % 3] === w_3[(i + 2) % 3]) {
// find maximum element of the independent eigenvector
maxVal = Math.abs(v_3x3[i][0]);
maxI = 0;
for (j = 1; j < 3; j++) {
if (maxVal < (tmp = Math.abs(v_3x3[i][j]))) {
maxVal = tmp;
maxI = j;
}
}
// swap the eigenvector into its proper position
if (maxI !== i) {
tmp = w_3[maxI];
w_3[maxI] = w_3[i];
w_3[i] = tmp;
vtkSwapVectors3(v_3x3[i], v_3x3[maxI]);
}
// maximum element of eigenvector should be positive
if (v_3x3[maxI][maxI] < 0) {
v_3x3[maxI][0] = -v_3x3[maxI][0];
v_3x3[maxI][1] = -v_3x3[maxI][1];
v_3x3[maxI][2] = -v_3x3[maxI][2];
}

// re-orthogonalize the other two eigenvectors
j = (maxI + 1) % 3;
k = (maxI + 2) % 3;

v_3x3[j][0] = 0.0;
v_3x3[j][1] = 0.0;
v_3x3[j][2] = 0.0;
v_3x3[j][j] = 1.0;
cross(v_3x3[maxI], v_3x3[j], v_3x3[k]);
normalize(v_3x3[k]);
cross(v_3x3[k], v_3x3[maxI], v_3x3[j]);

// transpose vectors back to columns
transpose3x3(v_3x3, v_3x3);
return;
}
}

// the three eigenvalues are different, just sort the eigenvectors
// to align them with the x, y, and z axes

// find the vector with the largest x element, make that vector
// the first vector
maxVal = Math.abs(v_3x3[0][0]);
maxI = 0;
for (i = 1; i < 3; i++) {
if (maxVal < (tmp = Math.abs(v_3x3[i][0]))) {
maxVal = tmp;
maxI = i;
}
}
// swap eigenvalue and eigenvector
if (maxI !== 0) {
tmp = w_3[maxI];
w_3[maxI] = w_3[0];
w_3[0] = tmp;
vtkSwapVectors3(v_3x3[maxI], v_3x3[0]);
}
// do the same for the y element
if (Math.abs(v_3x3[1][1]) < Math.abs(v_3x3[2][1])) {
tmp = w_3[2];
w_3[2] = w_3[1];
w_3[1] = tmp;
vtkSwapVectors3(v_3x3[2], v_3x3[1]);
}

// ensure that the sign of the eigenvectors is correct
for (i = 0; i < 2; i++) {
if (v_3x3[i][i] < 0) {
v_3x3[i][0] = -v_3x3[i][0];
v_3x3[i][1] = -v_3x3[i][1];
v_3x3[i][2] = -v_3x3[i][2];
}
}
// set sign of final eigenvector to ensure that determinant is positive
if (determinant3x3(v_3x3) < 0) {
v_3x3[2][0] = -v_3x3[2][0];
v_3x3[2][1] = -v_3x3[2][1];
v_3x3[2][2] = -v_3x3[2][2];
}

// transpose the eigenvectors back again
transpose3x3(v_3x3, v_3x3);
}

export function singularValueDecomposition3x3(a_3x3, u_3x3, w_3, vT_3x3) {
let i;
const B = [createArray(3), createArray(3), createArray(3)];

// copy so that A can be used for U or VT without risk
for (i = 0; i < 3; i++) {
B[0][i] = a_3x3[0][i];
B[1][i] = a_3x3[1][i];
B[2][i] = a_3x3[2][i];
}

// temporarily flip if determinant is negative
const d = determinant3x3(B);
if (d < 0) {
for (i = 0; i < 3; i++) {
B[0][i] = -B[0][i];
B[1][i] = -B[1][i];
B[2][i] = -B[2][i];
}
}

// orthogonalize, diagonalize, etc.
orthogonalize3x3(B, u_3x3);
transpose3x3(B, B);
multiply3x3_mat3(B, u_3x3, vT_3x3);
diagonalize3x3(vT_3x3, w_3, vT_3x3);
multiply3x3_mat3(u_3x3, vT_3x3, u_3x3);
transpose3x3(vT_3x3, vT_3x3);

// re-create the flip
if (d < 0) {
w_3[0] = -w_3[0];
w_3[1] = -w_3[1];
w_3[2] = -w_3[2];
}
}

export function luFactorLinearSystem(A, index, size) {
let i;
let j;
let k;
let largest;
let maxI = 0;
let sum;
let temp1;
let temp2;
const scale = createArray(size);

//
// Loop over rows to get implicit scaling information
//
for (i = 0; i < size; i++) {
for (largest = 0.0, j = 0; j < size; j++) {
if ((temp2 = Math.abs(A[i][j])) > largest) {
largest = temp2;
}
}

if (largest === 0.0) {
vtkWarningMacro('Unable to factor linear system');
return 0;
}
scale[i] = 1.0 / largest;
}
//
// Loop over all columns using Crout's method
//
for (j = 0; j < size; j++) {
for (i = 0; i < j; i++) {
sum = A[i][j];
for (k = 0; k < i; k++) {
sum -= A[i][k] * A[k][j];
}
A[i][j] = sum;
}
//
// Begin search for largest pivot element
//
for (largest = 0.0, i = j; i < size; i++) {
sum = A[i][j];
for (k = 0; k < j; k++) {
sum -= A[i][k] * A[k][j];
}
A[i][j] = sum;

if ((temp1 = scale[i] * Math.abs(sum)) >= largest) {
largest = temp1;
maxI = i;
}
}
//
// Check for row interchange
//
if (j !== maxI) {
for (k = 0; k < size; k++) {
temp1 = A[maxI][k];
A[maxI][k] = A[j][k];
A[j][k] = temp1;
}
scale[maxI] = scale[j];
}
//
// Divide by pivot element and perform elimination
//
index[j] = maxI;

if (Math.abs(A[j][j]) <= VTK_SMALL_NUMBER) {
vtkWarningMacro('Unable to factor linear system');
return 0;
}

if (j !== size - 1) {
temp1 = 1.0 / A[j][j];
for (i = j + 1; i < size; i++) {
A[i][j] *= temp1;
}
}
}
return 1;
}

export function luSolveLinearSystem(A, index, x, size) {
let i;
let j;
let ii;
let idx;
let sum;
//
// Proceed with forward and backsubstitution for L and U
// matrices. First, forward substitution.
//
for (ii = -1, i = 0; i < size; i++) {
idx = index[i];
sum = x[idx];
x[idx] = x[i];

if (ii >= 0) {
for (j = ii; j <= i - 1; j++) {
sum -= A[i][j] * x[j];
}
} else if (sum !== 0.0) {
ii = i;
}

x[i] = sum;
}
//
// Now, back substitution
//
for (i = size - 1; i >= 0; i--) {
sum = x[i];
for (j = i + 1; j < size; j++) {
sum -= A[i][j] * x[j];
}
x[i] = sum / A[i][i];
}
}

export function solveLinearSystem(A, x, size) {
// if we solving something simple, just solve it
if (size === 2) {
const y = createArray(2);
const det = determinant2x2(A[0][0], A[0][1], A[1][0], A[1][1]);

if (det === 0.0) {
// Unable to solve linear system
return 0;
}

y[0] = (A[1][1] * x[0] - A[0][1] * x[1]) / det;
y[1] = (-(A[1][0] * x[0]) + A[0][0] * x[1]) / det;

x[0] = y[0];
x[1] = y[1];
return 1;
}

if (size === 1) {
if (A[0][0] === 0.0) {
// Unable to solve linear system
return 0;
}

x[0] /= A[0][0];
return 1;
}

//
// System of equations is not trivial, use Crout's method
//

// Check on allocation of working vectors
const index = createArray(size);

// Factor and solve matrix
if (luFactorLinearSystem(A, index, size) === 0) {
return 0;
}
luSolveLinearSystem(A, index, x, size);

return 1;
}

export function invertMatrix(A, AI, size, index = null, column = null) {
const tmp1Size = index || createArray(size);
const tmp2Size = column || createArray(size);

// Factor matrix; then begin solving for inverse one column at a time.
// Note: tmp1Size returned value is used later, tmp2Size is just working
// memory whose values are not used in LUSolveLinearSystem
if (luFactorLinearSystem(A, tmp1Size, size, tmp2Size) === 0) {
return 0;
}

for (let j = 0; j < size; j++) {
for (let i = 0; i < size; i++) {
tmp2Size[i] = 0.0;
}
tmp2Size[j] = 1.0;

luSolveLinearSystem(A, tmp1Size, tmp2Size, size);

for (let i = 0; i < size; i++) {
AI[i][j] = tmp2Size[i];
}
}

return 1;
}

export function estimateMatrixCondition(A, size) {
let minValue = +Number.MAX_VALUE;
let maxValue = -Number.MAX_VALUE;

// find the maximum value
for (let i = 0; i < size; i++) {
for (let j = i; j < size; j++) {
if (Math.abs(A[i][j]) > max) {
maxValue = Math.abs(A[i][j]);
}
}
}

// find the minimum diagonal value
for (let i = 0; i < size; i++) {
if (Math.abs(A[i][i]) < min) {
minValue = Math.abs(A[i][i]);
}
}

if (minValue === 0.0) {
return Number.MAX_VALUE;
}
return maxValue / minValue;
}

export function jacobi(a_3x3, w, v) {
return jacobiN(a_3x3, 3, w, v);
}

export function solveHomogeneousLeastSquares(numberOfSamples, xt, xOrder, mt) {
// check dimensional consistency
if (numberOfSamples < xOrder) {
vtkWarningMacro('Insufficient number of samples. Underdetermined.');
return 0;
}

let i;
let j;
let k;

// set up intermediate variables
// Allocate matrix to hold X times transpose of X
const XXt = createArray(xOrder); // size x by x
// Allocate the array of eigenvalues and eigenvectors
const eigenvals = createArray(xOrder);
const eigenvecs = createArray(xOrder);

// Clear the upper triangular region (and btw, allocate the eigenvecs as well)
for (i = 0; i < xOrder; i++) {
eigenvecs[i] = createArray(xOrder);
XXt[i] = createArray(xOrder);
for (j = 0; j < xOrder; j++) {
XXt[i][j] = 0.0;
}
}

// Calculate XXt upper half only, due to symmetry
for (k = 0; k < numberOfSamples; k++) {
for (i = 0; i < xOrder; i++) {
for (j = i; j < xOrder; j++) {
XXt[i][j] += xt[k][i] * xt[k][j];
}
}
}

// now fill in the lower half of the XXt matrix
for (i = 0; i < xOrder; i++) {
for (j = 0; j < i; j++) {
XXt[i][j] = XXt[j][i];
}
}

// Compute the eigenvectors and eigenvalues
jacobiN(XXt, xOrder, eigenvals, eigenvecs);

// Smallest eigenval is at the end of the list (xOrder-1), and solution is
// corresponding eigenvec.
for (i = 0; i < xOrder; i++) {
mt[i][0] = eigenvecs[i][xOrder - 1];
}

return 1;
}

export function solveLeastSquares(
numberOfSamples,
xt,
xOrder,
yt,
yOrder,
mt,
checkHomogeneous = true
) {
// check dimensional consistency
if (numberOfSamples < xOrder || numberOfSamples < yOrder) {
vtkWarningMacro('Insufficient number of samples. Underdetermined.');
return 0;
}

const homogenFlags = createArray(yOrder);
let allHomogeneous = 1;
let hmt;
let homogRC = 0;
let i;
let j;
let k;
let someHomogeneous = 0;

// Ok, first init some flags check and see if all the systems are homogeneous
if (checkHomogeneous) {
// If Y' is zero, it's a homogeneous system and can't be solved via
// the pseudoinverse method. Detect this case, warn the user, and
// invoke SolveHomogeneousLeastSquares instead. Note that it doesn't
// really make much sense for yOrder to be greater than one in this case,
// since that's just yOrder occurrences of a 0 vector on the RHS, but
// we allow it anyway. N

// Initialize homogeneous flags on a per-right-hand-side basis
for (j = 0; j < yOrder; j++) {
homogenFlags[j] = 1;
}
for (i = 0; i < numberOfSamples; i++) {
for (j = 0; j < yOrder; j++) {
if (Math.abs(yt[i][j]) > VTK_SMALL_NUMBER) {
allHomogeneous = 0;
homogenFlags[j] = 0;
}
}
}

// If we've got one system, and it's homogeneous, do it and bail out quickly.
if (allHomogeneous && yOrder === 1) {
vtkWarningMacro(
'Detected homogeneous system (Y=0), calling SolveHomogeneousLeastSquares()'
);
return solveHomogeneousLeastSquares(numberOfSamples, xt, xOrder, mt);
}

// Ok, we've got more than one system of equations.
// Figure out if we need to calculate the homogeneous equation solution for
// any of them.
if (allHomogeneous) {
someHomogeneous = 1;
} else {
for (j = 0; j < yOrder; j++) {
if (homogenFlags[j]) {
someHomogeneous = 1;
}
}
}
}

// If necessary, solve the homogeneous problem
if (someHomogeneous) {
// hmt is the homogeneous equation version of mt, the general solution.
hmt = createArray(xOrder);
for (j = 0; j < xOrder; j++) {
// Only allocate 1 here, not yOrder, because here we're going to solve
// just the one homogeneous equation subset of the entire problem
hmt[j] = [0];
}

// Ok, solve the homogeneous problem
homogRC = solveHomogeneousLeastSquares(numberOfSamples, xt, xOrder, hmt);
}

// set up intermediate variables
const XXt = createArray(xOrder); // size x by x
const XXtI = createArray(xOrder); // size x by x
const XYt = createArray(xOrder); // size x by y
for (i = 0; i < xOrder; i++) {
XXt[i] = createArray(xOrder);
XXtI[i] = createArray(xOrder);

for (j = 0; j < xOrder; j++) {
XXt[i][j] = 0.0;
XXtI[i][j] = 0.0;
}

XYt[i] = createArray(yOrder);
for (j = 0; j < yOrder; j++) {
XYt[i][j] = 0.0;
}
}

// first find the pseudoinverse matrix
for (k = 0; k < numberOfSamples; k++) {
for (i = 0; i < xOrder; i++) {
// first calculate the XXt matrix, only do the upper half (symmetrical)
for (j = i; j < xOrder; j++) {
XXt[i][j] += xt[k][i] * xt[k][j];
}

// now calculate the XYt matrix
for (j = 0; j < yOrder; j++) {
XYt[i][j] += xt[k][i] * yt[k][j];
}
}
}

// now fill in the lower half of the XXt matrix
for (i = 0; i < xOrder; i++) {
for (j = 0; j < i; j++) {
XXt[i][j] = XXt[j][i];
}
}

const successFlag = invertMatrix(XXt, XXtI, xOrder);

// next get the inverse of XXt
if (successFlag) {
for (i = 0; i < xOrder; i++) {
for (j = 0; j < yOrder; j++) {
mt[i][j] = 0.0;
for (k = 0; k < xOrder; k++) {
mt[i][j] += XXtI[i][k] * XYt[k][j];
}
}
}
}

// Fix up any of the solutions that correspond to the homogeneous equation
// problem.
if (someHomogeneous) {
for (j = 0; j < yOrder; j++) {
if (homogenFlags[j]) {
// Fix this one
for (i = 0; i < xOrder; i++) {
mt[i][j] = hmt[i][0];
}
}
}
}

if (someHomogeneous) {
return homogRC && successFlag;
}

return successFlag;
}

export function hex2float(hexStr, outFloatArray = [0, 0.5, 1]) {
switch (hexStr.length) {
case 3: // abc => #aabbcc
outFloatArray[0] = (parseInt(hexStr[0], 16) * 17) / 255;
outFloatArray[1] = (parseInt(hexStr[1], 16) * 17) / 255;
outFloatArray[2] = (parseInt(hexStr[2], 16) * 17) / 255;
return outFloatArray;
case 4: // #abc => #aabbcc
outFloatArray[0] = (parseInt(hexStr[1], 16) * 17) / 255;
outFloatArray[1] = (parseInt(hexStr[2], 16) * 17) / 255;
outFloatArray[2] = (parseInt(hexStr[3], 16) * 17) / 255;
return outFloatArray;
case 6: // ab01df => #ab01df
outFloatArray[0] = parseInt(hexStr.substr(0, 2), 16) / 255;
outFloatArray[1] = parseInt(hexStr.substr(2, 2), 16) / 255;
outFloatArray[2] = parseInt(hexStr.substr(4, 2), 16) / 255;
return outFloatArray;
case 7: // #ab01df
outFloatArray[0] = parseInt(hexStr.substr(1, 2), 16) / 255;
outFloatArray[1] = parseInt(hexStr.substr(3, 2), 16) / 255;
outFloatArray[2] = parseInt(hexStr.substr(5, 2), 16) / 255;
return outFloatArray;
case 9: // #ab01df00
outFloatArray[0] = parseInt(hexStr.substr(1, 2), 16) / 255;
outFloatArray[1] = parseInt(hexStr.substr(3, 2), 16) / 255;
outFloatArray[2] = parseInt(hexStr.substr(5, 2), 16) / 255;
outFloatArray[3] = parseInt(hexStr.substr(7, 2), 16) / 255;
return outFloatArray;
default:
return outFloatArray;
}
}

export function rgb2hsv(rgb, hsv) {
let h;
let s;
const [r, g, b] = rgb;
const onethird = 1.0 / 3.0;
const onesixth = 1.0 / 6.0;
const twothird = 2.0 / 3.0;

let cmax = r;
let cmin = r;

if (g > cmax) {
cmax = g;
} else if (g < cmin) {
cmin = g;
}
if (b > cmax) {
cmax = b;
} else if (b < cmin) {
cmin = b;
}
const v = cmax;

if (v > 0.0) {
s = (cmax - cmin) / cmax;
} else {
s = 0.0;
}
if (s > 0) {
if (r === cmax) {
h = (onesixth * (g - b)) / (cmax - cmin);
} else if (g === cmax) {
h = onethird + (onesixth * (b - r)) / (cmax - cmin);
} else {
h = twothird + (onesixth * (r - g)) / (cmax - cmin);
}
if (h < 0.0) {
h += 1.0;
}
} else {
h = 0.0;
}

// Set the values back to the array
hsv[0] = h;
hsv[1] = s;
hsv[2] = v;
}

export function hsv2rgb(hsv, rgb) {
const [h, s, v] = hsv;
const onethird = 1.0 / 3.0;
const onesixth = 1.0 / 6.0;
const twothird = 2.0 / 3.0;
const fivesixth = 5.0 / 6.0;
let r;
let g;
let b;

// compute RGB from HSV
if (h > onesixth && h <= onethird) {
// green/red
g = 1.0;
r = (onethird - h) / onesixth;
b = 0.0;
} else if (h > onethird && h <= 0.5) {
// green/blue
g = 1.0;
b = (h - onethird) / onesixth;
r = 0.0;
} else if (h > 0.5 && h <= twothird) {
// blue/green
b = 1.0;
g = (twothird - h) / onesixth;
r = 0.0;
} else if (h > twothird && h <= fivesixth) {
// blue/red
b = 1.0;
r = (h - twothird) / onesixth;
g = 0.0;
} else if (h > fivesixth && h <= 1.0) {
// red/blue
r = 1.0;
b = (1.0 - h) / onesixth;
g = 0.0;
} else {
// red/green
r = 1.0;
g = h / onesixth;
b = 0.0;
}

// add Saturation to the equation.
r = s * r + (1.0 - s);
g = s * g + (1.0 - s);
b = s * b + (1.0 - s);

r *= v;
g *= v;
b *= v;

// Assign back to the array
rgb[0] = r;
rgb[1] = g;
rgb[2] = b;
}

export function lab2xyz(lab, xyz) {
// LAB to XYZ
const [L, a, b] = lab;
let var_Y = (L + 16) / 116;
let var_X = a / 500 + var_Y;
let var_Z = var_Y - b / 200;

if (var_Y ** 3 > 0.008856) {
var_Y **= 3;
} else {
var_Y = (var_Y - 16.0 / 116.0) / 7.787;
}

if (var_X ** 3 > 0.008856) {
var_X **= 3;
} else {
var_X = (var_X - 16.0 / 116.0) / 7.787;
}

if (var_Z ** 3 > 0.008856) {
var_Z **= 3;
} else {
var_Z = (var_Z - 16.0 / 116.0) / 7.787;
}
const ref_X = 0.9505;
const ref_Y = 1.0;
const ref_Z = 1.089;
xyz[0] = ref_X * var_X; // ref_X = 0.9505 Observer= 2 deg Illuminant= D65
xyz[1] = ref_Y * var_Y; // ref_Y = 1.000
xyz[2] = ref_Z * var_Z; // ref_Z = 1.089
}

export function xyz2lab(xyz, lab) {
const [x, y, z] = xyz;
const ref_X = 0.9505;
const ref_Y = 1.0;
const ref_Z = 1.089;
let var_X = x / ref_X; // ref_X = 0.9505 Observer= 2 deg, Illuminant= D65
let var_Y = y / ref_Y; // ref_Y = 1.000
let var_Z = z / ref_Z; // ref_Z = 1.089

if (var_X > 0.008856) var_X **= 1.0 / 3.0;
else var_X = 7.787 * var_X + 16.0 / 116.0;
if (var_Y > 0.008856) var_Y **= 1.0 / 3.0;
else var_Y = 7.787 * var_Y + 16.0 / 116.0;
if (var_Z > 0.008856) var_Z **= 1.0 / 3.0;
else var_Z = 7.787 * var_Z + 16.0 / 116.0;

lab[0] = 116 * var_Y - 16;
lab[1] = 500 * (var_X - var_Y);
lab[2] = 200 * (var_Y - var_Z);
}

export function xyz2rgb(xyz, rgb) {
const [x, y, z] = xyz;
let r = x * 3.2406 + y * -1.5372 + z * -0.4986;
let g = x * -0.9689 + y * 1.8758 + z * 0.0415;
let b = x * 0.0557 + y * -0.204 + z * 1.057;

// The following performs a "gamma correction" specified by the sRGB color
// space. sRGB is defined by a canonical definition of a display monitor and
// has been standardized by the International Electrotechnical Commission (IEC
// 61966-2-1). The nonlinearity of the correction is designed to make the
// colors more perceptually uniform. This color space has been adopted by
// several applications including Adobe Photoshop and Microsoft Windows color
// management. OpenGL is agnostic on its RGB color space, but it is reasonable
// to assume it is close to this one.
if (r > 0.0031308) r = 1.055 * r ** (1 / 2.4) - 0.055;
else r *= 12.92;
if (g > 0.0031308) g = 1.055 * g ** (1 / 2.4) - 0.055;
else g *= 12.92;
if (b > 0.0031308) b = 1.055 * b ** (1 / 2.4) - 0.055;
else b *= 12.92;

// Clip colors. ideally we would do something that is perceptually closest
// (since we can see colors outside of the display gamut), but this seems to
// work well enough.
let maxVal = r;
if (maxVal < g) maxVal = g;
if (maxVal < b) maxVal = b;
if (maxVal > 1.0) {
r /= maxVal;
g /= maxVal;
b /= maxVal;
}
if (r < 0) r = 0;
if (g < 0) g = 0;
if (b < 0) b = 0;

// Push values back to array
rgb[0] = r;
rgb[1] = g;
rgb[2] = b;
}

export function rgb2xyz(rgb, xyz) {
let [r, g, b] = rgb;
// The following performs a "gamma correction" specified by the sRGB color
// space. sRGB is defined by a canonical definition of a display monitor and
// has been standardized by the International Electrotechnical Commission (IEC
// 61966-2-1). The nonlinearity of the correction is designed to make the
// colors more perceptually uniform. This color space has been adopted by
// several applications including Adobe Photoshop and Microsoft Windows color
// management. OpenGL is agnostic on its RGB color space, but it is reasonable
// to assume it is close to this one.
if (r > 0.04045) r = ((r + 0.055) / 1.055) ** 2.4;
else r /= 12.92;
if (g > 0.04045) g = ((g + 0.055) / 1.055) ** 2.4;
else g /= 12.92;
if (b > 0.04045) b = ((b + 0.055) / 1.055) ** 2.4;
else b /= 12.92;

// Observer. = 2 deg, Illuminant = D65
xyz[0] = r * 0.4124 + g * 0.3576 + b * 0.1805;
xyz[1] = r * 0.2126 + g * 0.7152 + b * 0.0722;
xyz[2] = r * 0.0193 + g * 0.1192 + b * 0.9505;
}

export function rgb2lab(rgb, lab) {
const xyz = [0, 0, 0];
rgb2xyz(rgb, xyz);
xyz2lab(xyz, lab);
}

export function lab2rgb(lab, rgb) {
const xyz = [0, 0, 0];
lab2xyz(lab, xyz);
xyz2rgb(xyz, rgb);
}

export function uninitializeBounds(bounds) {
bounds[0] = 1.0;
bounds[1] = -1.0;
bounds[2] = 1.0;
bounds[3] = -1.0;
bounds[4] = 1.0;
bounds[5] = -1.0;
}

export function areBoundsInitialized(bounds) {
return !(bounds[1] - bounds[0] < 0.0);
}

export function computeBoundsFromPoints(point1, point2, bounds) {
bounds[0] = Math.min(point1[0], point2[0]);
bounds[1] = Math.max(point1[0], point2[0]);
bounds[2] = Math.min(point1[1], point2[1]);
bounds[3] = Math.max(point1[1], point2[1]);
bounds[4] = Math.min(point1[2], point2[2]);
bounds[5] = Math.max(point1[2], point2[2]);
}

export function clampValue(value, minValue, maxValue) {
if (value < minValue) {
return minValue;
}
if (value > maxValue) {
return maxValue;
}
return value;
}

export function clampVector(vector, minVector, maxVector, out = []) {
out[0] = clampValue(vector[0], minVector[0], maxVector[0]);
out[1] = clampValue(vector[1], minVector[1], maxVector[1]);
out[2] = clampValue(vector[2], minVector[2], maxVector[2]);

return out;
}

export function roundVector(vector, out = []) {
out[0] = Math.round(vector[0]);
out[1] = Math.round(vector[1]);
out[2] = Math.round(vector[2]);

return out;
}

export function clampAndNormalizeValue(value, range) {
let result = 0;
if (range[0] !== range[1]) {
// clamp
if (value < range[0]) {
result = range[0];
} else if (value > range[1]) {
result = range[1];
} else {
result = value;
}
// normalize
result = (result - range[0]) / (range[1] - range[0]);
}

return result;
}

export const getScalarTypeFittingRange = notImplemented(
'GetScalarTypeFittingRange'
);
export const getAdjustedScalarRange = notImplemented('GetAdjustedScalarRange');

export function extentIsWithinOtherExtent(extent1, extent2) {
if (!extent1 || !extent2) {
return 0;
}

for (let i = 0; i < 6; i += 2) {
if (
extent1[i] < extent2[i] ||
extent1[i] > extent2[i + 1] ||
extent1[i + 1] < extent2[i] ||
extent1[i + 1] > extent2[i + 1]
) {
return 0;
}
}

return 1;
}

export function boundsIsWithinOtherBounds(bounds1_6, bounds2_6, delta_3) {
if (!bounds1_6 || !bounds2_6) {
return 0;
}
for (let i = 0; i < 6; i += 2) {
if (
bounds1_6[i] + delta_3[i / 2] < bounds2_6[i] ||
bounds1_6[i] - delta_3[i / 2] > bounds2_6[i + 1] ||
bounds1_6[i + 1] + delta_3[i / 2] < bounds2_6[i] ||
bounds1_6[i + 1] - delta_3[i / 2] > bounds2_6[i + 1]
) {
return 0;
}
}
return 1;
}

export function pointIsWithinBounds(point_3, bounds_6, delta_3) {
if (!point_3 || !bounds_6 || !delta_3) {
return 0;
}
for (let i = 0; i < 3; i++) {
if (
point_3[i] + delta_3[i] < bounds_6[2 * i] ||
point_3[i] - delta_3[i] > bounds_6[2 * i + 1]
) {
return 0;
}
}
return 1;
}

export function solve3PointCircle(p1, p2, p3, center) {
const v21 = createArray(3);
const v32 = createArray(3);
const v13 = createArray(3);
const v12 = createArray(3);
const v23 = createArray(3);
const v31 = createArray(3);

for (let i = 0; i < 3; ++i) {
v21[i] = p1[i] - p2[i];
v32[i] = p2[i] - p3[i];
v13[i] = p3[i] - p1[i];
v12[i] = -v21[i];
v23[i] = -v32[i];
v31[i] = -v13[i];
}

const norm12 = norm(v12);
const norm23 = norm(v23);
const norm13 = norm(v13);

const crossv21v32 = createArray(3);
cross(v21, v32, crossv21v32);
const normCross = norm(crossv21v32);

const radius = (norm12 * norm23 * norm13) / (2 * normCross);

const normCross22 = 2 * normCross * normCross;
const alpha = (norm23 * norm23 * dot(v21, v31)) / normCross22;
const beta = (norm13 * norm13 * dot(v12, v32)) / normCross22;
const gamma = (norm12 * norm12 * dot(v13, v23)) / normCross22;

for (let i = 0; i < 3; ++i) {
center[i] = alpha * p1[i] + beta * p2[i] + gamma * p3[i];
}
return radius;
}

export const inf = Infinity;
export const negInf = -Infinity;

export const isInf = (value) => !Number.isFinite(value);
export const { isFinite, isNaN } = Number;
export const isNan = isNaN;

// JavaScript - add-on ----------------------

export function createUninitializedBounds() {
return [].concat([
Number.MAX_VALUE,
-Number.MAX_VALUE, // X
Number.MAX_VALUE,
-Number.MAX_VALUE, // Y
Number.MAX_VALUE,
-Number.MAX_VALUE, // Z
]);
}

export function getMajorAxisIndex(vector) {
let maxValue = -1;
let axisIndex = -1;
for (let i = 0; i < vector.length; i++) {
const value = Math.abs(vector[i]);
if (value > maxValue) {
axisIndex = i;
maxValue = value;
}
}

return axisIndex;
}

export function floatToHex2(value) {
const integer = Math.floor(value * 255);
if (integer > 15) {
return integer.toString(16);
}
return `0${integer.toString(16)}`;
}

export function floatRGB2HexCode(rgbArray, prefix = '#') {
return `${prefix}${rgbArray.map(floatToHex2).join('')}`;
}

function floatToChar(f) {
return Math.round(f * 255);
}

export function float2CssRGBA(rgbArray) {
if (rgbArray.length === 3) {
return `rgb(${rgbArray.map(floatToChar).join(', ')})`;
}
return `rgba(${floatToChar(rgbArray[0] || 0)}, ${floatToChar(
rgbArray[1] || 0
)}, ${floatToChar(rgbArray[2] || 0)}, ${rgbArray[3] || 0})`;
}

// ----------------------------------------------------------------------------
// Only Static API
// ----------------------------------------------------------------------------

export default {
Pi,
radiansFromDegrees,
degreesFromRadians,
round,
floor,
ceil,
ceilLog2,
min,
max,
arrayMin,
arrayMax,
arrayRange,
isPowerOfTwo,
nearestPowerOfTwo,
factorial,
binomial,
beginCombination,
nextCombination,
randomSeed,
getSeed,
random,
gaussian,
add,
subtract,
multiplyScalar,
multiplyScalar2D,
dot,
outer,
cross,
norm,
normalize,
perpendiculars,
projectVector,
projectVector2D,
distance2BetweenPoints,
angleBetweenVectors,
gaussianAmplitude,
gaussianWeight,
dot2D,
outer2D,
norm2D,
normalize2D,
determinant2x2,
LUFactor3x3,
LUSolve3x3,
linearSolve3x3,
multiply3x3_vect3,
multiply3x3_mat3,
multiplyMatrix,
transpose3x3,
invert3x3,
identity3x3,
determinant3x3,
quaternionToMatrix3x3,
areMatricesEqual,
matrix3x3ToQuaternion,
multiplyQuaternion,
orthogonalize3x3,
diagonalize3x3,
singularValueDecomposition3x3,
solveLinearSystem,
invertMatrix,
luFactorLinearSystem,
luSolveLinearSystem,
estimateMatrixCondition,
jacobi,
jacobiN,
solveHomogeneousLeastSquares,
solveLeastSquares,
hex2float,
rgb2hsv,
hsv2rgb,
lab2xyz,
xyz2lab,
xyz2rgb,
rgb2xyz,
rgb2lab,
lab2rgb,
uninitializeBounds,
areBoundsInitialized,
computeBoundsFromPoints,
clampValue,
clampVector,
clampAndNormalizeValue,
getScalarTypeFittingRange,
getAdjustedScalarRange,
extentIsWithinOtherExtent,
boundsIsWithinOtherBounds,
pointIsWithinBounds,
solve3PointCircle,
inf,
negInf,
isInf,
isNan: isNaN,
isNaN,
isFinite,

// JS add-on
createUninitializedBounds,
getMajorAxisIndex,
floatToHex2,
floatRGB2HexCode,
float2CssRGBA,
};