Math

Introduction

vtkMath provides methods to perform common math operations. These include
providing constants such as Pi; conversion from degrees to radians; vector
operations such as dot and cross products and vector norm; matrix determinant
for 2x2 and 3x3 matrices; univariate polynomial solvers; and for random
number generation (for backward compatibility only).

Methods

LUFactor3x3

LU Factorization of a 3x3 matrix.

Argument Type Required Description
mat_3x3 Matrix3x3 Yes
index_3 Vector3 Yes

LUSolve3x3

LU back substitution for a 3x3 matrix.

Argument Type Required Description
mat_3x3 Matrix3x3 Yes
index_3 Vector3 Yes
x_3 Vector3 Yes

Pi

Get the number π.

add

Addition of two 3-vectors.

Argument Type Required Description
a Vector3 Yes The first 3D vector.
b Vector3 Yes The second 3D vector.
out Vector3 Yes The output 3D vector.

angleBetweenVectors

Angle between 3D vectors.

Argument Type Required Description
v1 Vector3 Yes The first 3D vector.
v2 Vector3 Yes The second 3D vector.

areBoundsInitialized

Argument Type Required Description
bounds Bounds Yes The bounds to check.

areEquals

Returns true if elements of both arrays are equals.

Argument Type Required Description
a Array. Yes An array of numbers (vector, point, matrix…)
b Array. Yes An array of numbers (vector, point, matrix…)
eps Number No The tolerance value.

arrayMax

Get the maximum of the array.

Argument Type Required Description
arr Array. Yes The array.
offset Number Yes The offset.
stride Number Yes The stride.

arrayMin

Get the minimum of the array.

Argument Type Required Description
arr Array. Yes The array.
offset Number Yes The offset.
stride Number Yes The stride.

arrayRange

Argument Type Required Description
arr Array. Yes The array.
offset Number Yes The offset.
stride Number Yes The stride.

beginCombination

Start iterating over “m choose n” objects.

Argument Type Required Description
m Number No The first numeric expression.
n Number No The second numeric expression.

binomial

The number of combinations of n objects from a pool of m objects (m>n).

Argument Type Required Description
m Number Yes The first numeric expression.
n Number Yes The second numeric expression.

boundsIsWithinOtherBounds

Check if first 3D bounds is within the second 3D bounds.

Argument Type Required Description
bounds1_6 Bounds Yes The first bounds.
bounds2_6 Bounds Yes The second bounds.
delta_3 Vector3 Yes The error margin along each axis.

ceil

Same as Math.ceil().

Argument Type Required Description
param1 Number Yes A numeric expression.

ceilLog2

Gives the exponent of the lowest power of two not less than x.

clampAndNormalizeValue

Argument Type Required Description
value Number Yes
range Range Yes

clampValue

Clamp some value against a range.

Argument Type Required Description
value Number Yes The value to clamp.
minValue Number Yes The minimum value.
maxValue Number Yes The maximum value.

clampVector

Clamp some vector against a range.

Argument Type Required Description
vector Vector3 Yes The vector to clamp.
minVector Vector3 Yes The minimum vector.
maxVector Vector3 Yes The maximum vector.
out Vector3 Yes The output vector.

columnsToMat3

Fill a 3x3 matrix with the given column vectors

Argument Type Required Description
column0 Vector3 Yes
column1 Vector3 Yes
column2 Vector3 Yes
mat Matrix Yes

columnsToMat4

Fill a 4x4 matrix with the given column vectors

Argument Type Required Description
column0 Vector4 Yes
column1 Vector4 Yes
column2 Vector4 Yes
column3 Vector4 Yes
mat Matrix Yes

createArray

Argument Type Required Description
size Number Yes The size of the array.

createUninitializedBounds

cross

Computes cross product of 3D vectors x and y.

Argument Type Required Description
x Vector3 Yes The first 3D vector.
y Vector3 Yes The second 3D vector.
out Vector3 Yes The output 3D vector.

degreesFromRadians

Convert radians to degrees.

Argument Type Required Description
rad Number Yes The value in radians.

determinant2x2

Argument Type Required Description
args Array. Yes

determinant3x3

Calculate the determinant of a 3x3 matrix.

Argument Type Required Description
mat_3x3 Matrix3x3 Yes The input 3x3 matrix.

diagonalize3x3

Argument Type Required Description
a_3x3 Matrix3x3 Yes
w_3 Vector3 Yes
v_3x3 Matrix3x3 Yes

distance2BetweenPoints

Compute distance squared between two points p1 and p2.

Argument Type Required Description
x Vector3 Yes The first 3D vector.
y Vector3 Yes The second 3D vector.

dot

Argument Type Required Description
x Vector3 Yes
y Vector3 Yes

dot2D

Argument Type Required Description
x Vector2 Yes
y Vector2 Yes

estimateMatrixCondition

Argument Type Required Description
A Matrix Yes
size Number Yes

extentIsWithinOtherExtent

Check if first 3D extent is within second 3D extent.

Argument Type Required Description
extent1 Extent Yes The first extent.
extent2 Extent Yes The second extent.

factorial

Compute N factorial, N! = N*(N-1) * (N-2)…32*1.

float2CssRGBA

Convert RGB or RGBA color array to CSS representation

Argument Type Required Description
rgbArray RGBColor or RGBAColor Yes The color array.

floatRGB2HexCode

Argument Type Required Description
rgbArray RGBColor Yes
prefix string No

floatToHex2

Argument Type Required Description
value Number Yes The value to convert.

floor

Same as Math.floor().

Argument Type Required Description
param1 Number Yes A numeric expression.

gaussian

Generate pseudo-random numbers distributed according to the standard normal
distribution.

gaussianAmplitude

Compute the amplitude of a Gaussian function with mean=0 and specified variance.

Argument Type Required Description
mean Number Yes The mean value.
variance Number Yes The variance value.
position Number Yes The position value.

gaussianWeight

Compute the amplitude of an unnormalized Gaussian function with mean=0 and
specified variance.

Argument Type Required Description
mean Number Yes The mean value.
variance Number Yes The variance value.
position Number Yes The position value.

getAdjustedScalarRange

getMajorAxisIndex

Argument Type Required Description
vector Array. Yes

getScalarTypeFittingRange

Get the scalar type that is most likely to have enough precision to store a
given range of data once it has been scaled and shifted

getSeed

Return the current seed used by the random number generator.

getSparseOrthogonalMatrix

Return the closest orthogonal matrix of 1, -1 and 0
It works for both column major and row major matrices
This function iteratively associate a column with a row by choosing
the greatest absolute value from the remaining row and columns
For each association, a -1 or a 1 is set in the output, depending on
the sign of the value in the original matrix

Argument Type Required Description
matrix Array. Yes The matrix of size nxn
n Array. Yes The size of the square matrix, defaults to 3

hex2float

Argument Type Required Description
hexStr String Yes
outFloatArray Array. No

hsv2rgb

Argument Type Required Description
hsv HSVColor Yes An Array of the HSV color.
rgb RGBColor Yes An Array of the RGB color.

identity

Set mat to the identity matrix.

Argument Type Required Description
n Number Yes The size of the matrix.
mat Array. Yes The output matrix.

identity3x3

Set mat_3x3 to the identity matrix.

Argument Type Required Description
mat_3x3 Matrix3x3 Yes The input 3x3 matrix.

invert3x3

Invert a 3x3 matrix.

Argument Type Required Description
in_3x3 Matrix3x3 Yes The input 3x3 matrix.
outI_3x3 Matrix3x3 Yes The output 3x3 matrix.

invertMatrix

Argument Type Required Description
A Matrix Yes The input matrix. It is modified during the inversion.
AI Matrix Yes The output inverse matrix. Can be the same as input matrix.
size Number No The square size of the matrix to invert : 4 for a 4x4
index Array. No
column Array. No

isFinite

Determines whether the passed value is a finite number.

Argument Type Required Description
value Yes The value to check.

isIdentity

Returns true if provided matrix is the identity matrix.

Argument Type Required Description
mat Array. Yes The 3x3 matrix to check
eps Number No The tolerance value.

isIdentity3x3

Returns true if provided 3x3 matrix is the identity matrix.

Argument Type Required Description
mat Matrix3x3 Yes The 3x3 matrix to check
eps Number No The tolerance value.

isInf

Determines whether the passed value is a infinite number.

Argument Type Required Description
value Number Yes The value to check.

isNaN

Determines whether the passed value is a NaN.

Argument Type Required Description
value Number Yes The value to check.

isNan

Determines whether the passed value is a NaN.

Argument Type Required Description
value Number Yes The value to check.

isPowerOfTwo

Returns true if integer is a power of two.

Argument Type Required Description
x Number Yes A numeric expression.

jacobi

Argument Type Required Description
a_3x3 Matrix3x3 Yes
w Array. Yes
v Array. Yes

jacobiN

Jacobi iteration for the solution of eigenvectors/eigenvalues. Input matrix a is modified (the upper triangle is filled with zeros)

Argument Type Required Description
a Matrix Yes real symetric nxn matrix
n Number Yes matrix size
w Array. Yes vector of size n to store eigenvalues (stored in decreasing order)
v Array. Yes matrix of size nxn to store eigenvectors (stored in decreasing order, normalized)

lab2rgb

Argument Type Required Description
lab Vector3 Yes
rgb RGBColor Yes An Array of the RGB color.

lab2xyz

Argument Type Required Description
lab Vector3 Yes
xyz Vector3 Yes

linearSolve3x3

Solve mat_3x3y_3 = x_3 for y and place the result in y.

Argument Type Required Description
mat_3x3 Matrix3x3 Yes
x_3 Vector3 Yes
y_3 Vector3 Yes

luFactorLinearSystem

Argument Type Required Description
A Matrix Yes
index Array. Yes
size Number Yes

luSolveLinearSystem

Argument Type Required Description
A Matrix Yes
index Array. Yes
x Array. Yes
size Number Yes

matrix3x3ToQuaternion

Argument Type Required Description
mat_3x3 Matrix3x3 Yes
quat_4 Vector4 Yes

max

Get the maximum of the two arguments provided. If either argument is NaN,
the first argument will always be returned.

Argument Type Required Description
param1 Number Yes The first numeric expression.
param2 Number Yes The second numeric expression.

min

Get the minimum of the two arguments provided. If either argument is NaN,
the first argument will always be returned.

Argument Type Required Description
param1 Number Yes The first numeric expression.
param2 Number Yes The second numeric expression.

multiply3x3_mat3

Argument Type Required Description
a_3x3 Matrix3x3 Yes
b_3x3 Matrix3x3 Yes
out_3x3 Matrix3x3 Yes

multiply3x3_vect3

Argument Type Required Description
mat_3x3 Matrix3x3 Yes
in_3 Vector3 Yes
out_3 Vector3 Yes

multiplyAccumulate

Argument Type Required Description
a Vector3 Yes
b Vector3 Yes
scalar Number Yes
out Vector3 Yes

multiplyAccumulate2D

Argument Type Required Description
a Vector2 Yes
b Vector2 Yes
scalar Number Yes
out Vector2 Yes

multiplyMatrix

Multiply two matrices.

Argument Type Required Description
a Matrix Yes
b Matrix Yes
rowA Number Yes
colA Number Yes
rowB Number Yes
colB Number Yes
outRowAColB Matrix Yes

multiplyQuaternion

Argument Type Required Description
quat_1 Vector4 Yes
quat_2 Vector4 Yes
quat_out Vector4 Yes

multiplyScalar

Argument Type Required Description
vec Vector3 Yes
scalar Number Yes

multiplyScalar2D

Argument Type Required Description
vec Vector2 Yes
scalar Number Yes

nearestPowerOfTwo

Compute the nearest power of two that is not less than x.

Argument Type Required Description
xi Number Yes A numeric expression.

nextCombination

Given m, n, and a valid combination of n integers in the range [0,m[, this
function alters the integers into the next combination in a sequence of all
combinations of n items from a pool of m.

Argument Type Required Description
m Number Yes The first numeric expression.
n Number Yes The second numeric expression.
r Array. Yes

norm

Argument Type Required Description
x Array. Yes
n Number Yes

norm2D

Compute the norm of a 2-vector.

Argument Type Required Description
x2D Vector2 Yes x The 2D vector.

normalize

Normalize in place. Returns norm.

Argument Type Required Description
x Vector3 Yes The vector to normlize.

normalize2D

Normalize (in place) a 2-vector.

Argument Type Required Description
x Vector2 Yes The 2D vector.

orthogonalize3x3

Argument Type Required Description
a_3x3 Matrix3x3 Yes
out_3x3 Matrix3x3 Yes

outer

Outer product of two 3-vectors.

Argument Type Required Description
x Vector3 Yes The first 3D vector.
y Vector3 Yes The second 3D vector.
out_3x3 Matrix3x3 Yes The output 3x3 matrix.

outer2D

Outer product of two 2-vectors.

Argument Type Required Description
x Vector2 Yes The first 2D vector.
y Vector2 Yes The second 2D vector.
out_2x2 Matrix Yes The output 2x2 matrix.

perpendiculars

Given a unit vector v1, find two unit vectors v2 and v3 such that v1 cross v2 = v3

Argument Type Required Description
x Vector3 Yes The first vector.
y Vector3 Yes The second vector.
z Vector3 Yes The third vector.
theta Number Yes

pointIsWithinBounds

Check if point is within the given 3D bounds.

Argument Type Required Description
point_3 Vector3 Yes The coordinate of the point.
bounds_6 Bounds Yes The bounds.
delta_3 Vector3 Yes The error margin along each axis.

projectVector

Argument Type Required Description
a Vector3 Yes
b Vector3 Yes
projection Vector3 Yes

projectVector2D

Compute the projection of 2D vector a on 2D vector b and returns the
result in projection vector.

Argument Type Required Description
a Vector2 Yes The first 2D vector.
b Vector2 Yes The second 2D vector.
projection Vector2 Yes The projection 2D vector.

quaternionToMatrix3x3

Argument Type Required Description
quat_4 Vector4 Yes
mat_3x3 Matrix3x3 Yes

radiansFromDegrees

Convert degrees to radians.

Argument Type Required Description
deg Number Yes The value in degrees.

random

Generate pseudo-random numbers distributed according to the uniform
distribution between min and max.

Argument Type Required Description
minValue Number Yes
maxValue Number Yes

randomSeed

Initialize seed value.

Argument Type Required Description
seed Number Yes

rgb2hsv

Argument Type Required Description
rgb RGBColor Yes An Array of the RGB color.
hsv HSVColor Yes An Array of the HSV color.

rgb2lab

Argument Type Required Description
rgb RGBColor Yes
lab Vector3 Yes

rgb2xyz

Argument Type Required Description
rgb RGBColor Yes An Array of the RGB color.
xyz Vector3 Yes

round

Same as Math.round().

Argument Type Required Description
param1 Number Yes The value to be rounded to the nearest integer.

roundNumber

Argument Type Required Description
num Number Yes
digits Number No

roundVector

Argument Type Required Description
vector Vector3 Yes
out Vector3 Yes

roundVector

Argument Type Required Description
vector Vector3 Yes
out Vector3 No
digits Number No

rowsToMat3

Fill a 3x3 matrix with the given row vectors

Argument Type Required Description
row0 Vector3 Yes
row1 Vector3 Yes
row2 Vector3 Yes
mat Matrix Yes

rowsToMat4

Fill a 4x4 matrix with the given row vectors

Argument Type Required Description
row0 Vector4 Yes
row1 Vector4 Yes
row2 Vector4 Yes
row3 Vector4 Yes
mat Matrix Yes

signedAngleBetweenVectors

Signed angle between v1 and v2 with regards to plane defined by normal vN.
angle between v1 and v2 with regards to plane defined by normal vN.Signed
angle between v1 and v2 with regards to plane defined by normal
vN.t3(mat_3x3, in_3, out_3)

Argument Type Required Description
v1 Vector3 Yes The first 3D vector.
v2 Vector3 Yes The second 3D vector.
vN Vector3 Yes

singularValueDecomposition3x3

Argument Type Required Description
a_3x3 Matrix3x3 Yes
u_3x3 Matrix3x3 Yes
w_3 Vector3 Yes
vT_3x3 Matrix3x3 Yes

solve3PointCircle

In Euclidean space, there is a unique circle passing through any given three
non-collinear points P1, P2, and P3.

Using Cartesian coordinates to represent these points as spatial vectors, it
is possible to use the dot product and cross product to calculate the radius
and center of the circle. See:
http://en.wikipedia.org/wiki/Circumscribed_circle and more specifically the
section Barycentric coordinates from cross- and dot-products

Argument Type Required Description
p1 Vector3 Yes The coordinate of the first point.
p2 Vector3 Yes The coordinate of the second point.
p3 Vector3 Yes The coordinate of the third point.
center Vector3 Yes The coordinate of the center point.

solveHomogeneousLeastSquares

Solves for the least squares best fit matrix for the homogeneous equation
X’M’ = 0’. Uses the method described on pages 40-41 of Computer Vision by
Forsyth and Ponce, which is that the solution is the eigenvector associated
with the minimum eigenvalue of T(X)X, where T(X) is the transpose of X. The
inputs and output are transposed matrices. Dimensions: X’ is numberOfSamples
by xOrder, M’ dimension is xOrder by yOrder. M’ should be pre-allocated. All
matrices are row major. The resultant matrix M’ should be pre-multiplied to
X’ to get 0’, or transposed and then post multiplied to X to get 0

Argument Type Required Description
numberOfSamples Number Yes
xt Matrix Yes
xOrder Number Yes
mt Matrix Yes

solveLeastSquares

Solves for the least squares best fit matrix for the equation X’M’ = Y’. Uses
pseudoinverse to get the ordinary least squares. The inputs and output are
transposed matrices. Dimensions: X’ is numberOfSamples by xOrder, Y’ is
numberOfSamples by yOrder, M’ dimension is xOrder by yOrder. M’ should be
pre-allocated. All matrices are row major. The resultant matrix M’ should be
pre-multiplied to X’ to get Y’, or transposed and then post multiplied to X
to get Y By default, this method checks for the homogeneous condition where
Y==0, and if so, invokes SolveHomogeneousLeastSquares. For better performance
when the system is known not to be homogeneous, invoke with
checkHomogeneous=0.

Argument Type Required Description
numberOfSamples Number Yes
xt Matrix Yes
xOrder Number Yes
yt Matrix Yes
yOrder Number Yes
mt Matrix Yes
checkHomogeneous Boolean No

solveLinearSystem

Argument Type Required Description
A Matrix Yes
x Array. Yes
size Number Yes

subtract

Subtraction of two 3-vectors.

Argument Type Required Description
a Vector3 Yes The first 3D vector.
b Vector3 Yes The second 3D vector.
out Vector3 Yes The output 3D vector.

swapColumnsMatrix_nxn

Given two columns indices, swap the two columns of a nxn matrix

Argument Type Required Description
matrix Array. Yes The n by n matrix in wich we want to swap the vectors.
n Number Yes size of the matrix.
column1 Number Yes index of first col to swap with the other.
column2 Number Yes index of second col to swap with the other.

swapRowsMatrix_nxn

Given two rows indices, swap the two rows of a nxn matrix

Argument Type Required Description
matrix Array. Yes The n by n matrix in wich we want to swap the vectors.
n Number Yes size of the matrix.
row1 Number Yes index of first row to swap with the other.
row2 Number Yes index of second row to swap with the other.

transpose3x3

Transpose a 3x3 matrix.

Argument Type Required Description
in_3x3 Matrix3x3 Yes The input 3x3 matrix.
outT_3x3 Matrix3x3 Yes The output 3x3 matrix.

uninitializeBounds

Returns bounds.

Argument Type Required Description
bounds Bounds Yes Output array that hold bounds, optionally empty.

xyz2lab

Argument Type Required Description
xyz Vector3 Yes
lab Vector3 Yes

xyz2rgb

Argument Type Required Description
xyz Vector3 Yes
rgb RGBColor Yes An Array of the RGB color.

Source

Constants.js
export const IDENTITY = [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1];
export const IDENTITY_3X3 = [1, 0, 0, 0, 1, 0, 0, 0, 1];

export const EPSILON = 1e-6;
export const VTK_SMALL_NUMBER = 1.0e-12;

export default {
IDENTITY,
IDENTITY_3X3,
EPSILON,
VTK_SMALL_NUMBER,
};
index.d.ts
import { Bounds, Extent, HSVColor, RGBAColor, RGBColor, Matrix, Matrix3x3, Range, Vector2, Vector3, Vector4 } from "../../../types";

/**
*
* @param {Number} size The size of the array.
*/
export function createArray(size?: number): number[];

/**
* Given two rows indices, swap the two rows of a nxn matrix
* @param {Number[]} matrix The n by n matrix in wich we want to swap the vectors.
* @param {Number} n size of the matrix.
* @param {Number} row1 index of first row to swap with the other.
* @param {Number} row2 index of second row to swap with the other.
*/
export function swapRowsMatrix_nxn(matrix: number[], n: number, row1: number, row2: number): void;

/**
* Given two columns indices, swap the two columns of a nxn matrix
* @param {Number[]} matrix The n by n matrix in wich we want to swap the vectors.
* @param {Number} n size of the matrix.
* @param {Number} column1 index of first col to swap with the other.
* @param {Number} column2 index of second col to swap with the other.
*/
export function swapColumnsMatrix_nxn(matrix: number[], n: number, column1: number, column2: number): void;

/**
* Get the number π.
*/
export function Pi(): number;

/**
* Convert degrees to radians.
* @param {Number} deg The value in degrees.
*/
export function radiansFromDegrees(deg: number): number;

/**
* Convert radians to degrees.
* @param {Number} rad The value in radians.
*/
export function degreesFromRadians(rad: number): number;

/**
* Same as Math.round().
* @param {Number} param1 The value to be rounded to the nearest integer.
*/
export function round(param1: number): number;

/**
* Same as Math.floor().
* @param {Number} param1 A numeric expression.
*/
export function floor(param1: number): number;

/**
* Same as Math.ceil().
* @param {Number} param1 A numeric expression.
*/
export function ceil(param1: number): number;

/**
* Get the minimum of the two arguments provided. If either argument is NaN,
* the first argument will always be returned.
* @param {Number} param1 The first numeric expression.
* @param {Number} param2 The second numeric expression.
*/
export function min(param1: number, param2: number): number;

/**
* Get the maximum of the two arguments provided. If either argument is NaN,
* the first argument will always be returned.
* @param {Number} param1 The first numeric expression.
* @param {Number} param2 The second numeric expression.
*/
export function max(param1: number, param2: number): number;

/**
* Get the minimum of the array.
* @param {Number[]} arr The array.
* @param {Number} offset The offset.
* @param {Number} stride The stride.
*/
export function arrayMin(arr: number[], offset: number, stride: number): number;

/**
* Get the maximum of the array.
* @param {Number[]} arr The array.
* @param {Number} offset The offset.
* @param {Number} stride The stride.
*/
export function arrayMax(arr: number[], offset: number, stride: number): number;

/**
*
* @param {Number[]} arr The array.
* @param {Number} offset The offset.
* @param {Number} stride The stride.
*/
export function arrayRange(arr: number[], offset: number, stride: number): number[];

/**
* Gives the exponent of the lowest power of two not less than x.
*
*/
export function ceilLog2(): void;

/**
* Compute N factorial, N! = N*(N-1) * (N-2)...*3*2*1.
*/
export function factorial(): void;

/**
* Generate pseudo-random numbers distributed according to the standard normal
* distribution.
*/
export function gaussian(): void;

/**
* Compute the nearest power of two that is not less than x.
* @param {Number} xi A numeric expression.
*/
export function nearestPowerOfTwo(xi: number): number;

/**
* Returns true if integer is a power of two.
* @param {Number} x A numeric expression.
*/
export function isPowerOfTwo(x: number): boolean;

/**
* The number of combinations of n objects from a pool of m objects (m>n).
* @param {Number} m The first numeric expression.
* @param {Number} n The second numeric expression.
*/
export function binomial(m: number, n: number): number;

/**
* Start iterating over "m choose n" objects.
* @param {Number} [m] The first numeric expression.
* @param {Number} [n] The second numeric expression.
*/
export function beginCombination(m?: number, n?: number): number;

/**
* Given m, n, and a valid combination of n integers in the range [0,m[, this
* function alters the integers into the next combination in a sequence of all
* combinations of n items from a pool of m.
* @param {Number} m The first numeric expression.
* @param {Number} n The second numeric expression.
* @param {Number[]} r
*/
export function nextCombination(m: number, n: number, r: number[]): number;

/**
* Initialize seed value.
* @param {Number} seed
*/
export function randomSeed(seed: number): void;

/**
* Return the current seed used by the random number generator.
*/
export function getSeed(): number;

/**
* Generate pseudo-random numbers distributed according to the uniform
* distribution between min and max.
* @param {Number} minValue
* @param {Number} maxValue
*/
export function random(minValue: number, maxValue: number): number;

/**
* Addition of two 3-vectors.
* @param {Vector3} a The first 3D vector.
* @param {Vector3} b The second 3D vector.
* @param {Vector3} out The output 3D vector.
* @example
* ```js
* a[3] + b[3] => out[3]
* ```
*/
export function add(a: Vector3, b: Vector3, out: Vector3): Vector3;

/**
* Subtraction of two 3-vectors.
* @param {Vector3} a The first 3D vector.
* @param {Vector3} b The second 3D vector.
* @param {Vector3} out The output 3D vector.
* @example
* ```js
* a[3] - b[3] => out[3]
* ```
*/
export function subtract(a: Vector3, b: Vector3, out: Vector3): Vector3;

/**
*
* @param {Vector3} vec
* @param {Number} scalar
* @example
* ```js
* vec[3] * scalar => vec[3]
* ```
*/
export function multiplyScalar(vec: Vector3, scalar: number): Vector3;

/**
*
* @param {Vector2} vec
* @param {Number} scalar
* @example
* ```js
* vec[3] * scalar => vec[3]
* ```
*/
export function multiplyScalar2D(vec: Vector2, scalar: number): Vector2;

/**
*
* @param {Vector3} a
* @param {Vector3} b
* @param {Number} scalar
* @param {Vector3} out
* @example
* ```js
* a[3] + b[3] * scalar => out[3]
* ```
*/
export function multiplyAccumulate(a: Vector3, b: Vector3, scalar: number, out: Vector3): Vector3;

/**
*
* @param {Vector2} a
* @param {Vector2} b
* @param {Number} scalar
* @param {Vector2} out
* @example
* ```js
* a[2] + b[2] * scalar => out[2]
* ```
*/
export function multiplyAccumulate2D(a: Vector2, b: Vector2, scalar: number, out: Vector2): Vector2;

/**
*
* @param {Vector3} x
* @param {Vector3} y
* @example
* ```js
* a[2] + b[2] * scalar => out[2]
* ```
*/
export function dot(x: Vector3, y: Vector3): number;

/**
* Outer product of two 3-vectors.
* @param {Vector3} x The first 3D vector.
* @param {Vector3} y The second 3D vector.
* @param {Matrix3x3} out_3x3 The output 3x3 matrix.
*/
export function outer(x: Vector3, y: Vector3, out_3x3: Matrix3x3): void;

/**
* Computes cross product of 3D vectors x and y.
* @param {Vector3} x The first 3D vector.
* @param {Vector3} y The second 3D vector.
* @param {Vector3} out The output 3D vector.
*/
export function cross(x: Vector3, y: Vector3, out: Vector3): Vector3;

/**
*
* @param {Number[]} x
* @param {Number} n
*/
export function norm(x: number[], n: number): number;

/**
* Normalize in place. Returns norm.
* @param {Vector3} x The vector to normlize.
*/
export function normalize(x: Vector3): number;

/**
* Given a unit vector v1, find two unit vectors v2 and v3 such that v1 cross v2 = v3
* @param {Vector3} x The first vector.
* @param {Vector3} y The second vector.
* @param {Vector3} z The third vector.
* @param {Number} theta
*/
export function perpendiculars(x: Vector3, y: Vector3, z: Vector3, theta: number): void;

/**
*
* @param {Vector3} a
* @param {Vector3} b
* @param {Vector3} projection
*/
export function projectVector(a: Vector3, b: Vector3, projection: Vector3): boolean;

/**
*
* @param {Vector2} x
* @param {Vector2} y
*/
export function dot2D(x: Vector2, y: Vector2): number;

/**
* Compute the projection of 2D vector `a` on 2D vector `b` and returns the
* result in projection vector.
* @param {Vector2} a The first 2D vector.
* @param {Vector2} b The second 2D vector.
* @param {Vector2} projection The projection 2D vector.
*/
export function projectVector2D(a: Vector2, b: Vector2, projection: Vector2): boolean;

/**
* Compute distance squared between two points p1 and p2.
* @param {Vector3} x The first 3D vector.
* @param {Vector3} y The second 3D vector.
*/
export function distance2BetweenPoints(x: Vector3, y: Vector3): number;

/**
* Angle between 3D vectors.
* @param {Vector3} v1 The first 3D vector.
* @param {Vector3} v2 The second 3D vector.
*/
export function angleBetweenVectors(v1: Vector3, v2: Vector3): number;


/**
* Signed angle between v1 and v2 with regards to plane defined by normal vN.
* angle between v1 and v2 with regards to plane defined by normal vN.Signed
* angle between v1 and v2 with regards to plane defined by normal
* vN.t3(mat_3x3, in_3, out_3)
* @param {Vector3} v1 The first 3D vector.
* @param {Vector3} v2 The second 3D vector.
* @param {Vector3} vN
*/
export function signedAngleBetweenVectors(v1: Vector3, v2: Vector3, vN: Vector3): number;


/**
* Compute the amplitude of a Gaussian function with mean=0 and specified variance.
* @param {Number} mean The mean value.
* @param {Number} variance The variance value.
* @param {Number} position The position value.
*/
export function gaussianAmplitude(mean: number, variance: number, position: number): number;

/**
* Compute the amplitude of an unnormalized Gaussian function with mean=0 and
* specified variance.
* @param {Number} mean The mean value.
* @param {Number} variance The variance value.
* @param {Number} position The position value.
*/
export function gaussianWeight(mean: number, variance: number, position: number): number;

/**
* Outer product of two 2-vectors.
* @param {Vector2} x The first 2D vector.
* @param {Vector2} y The second 2D vector.
* @param {Matrix} out_2x2 The output 2x2 matrix.
*/
export function outer2D(x: Vector2, y: Vector2, out_2x2: Matrix): void;

/**
* Compute the norm of a 2-vector.
* @param {Vector2} x2D x The 2D vector.
*/
export function norm2D(x2D: Vector2): number;

/**
* Normalize (in place) a 2-vector.
* @param {Vector2} x The 2D vector.
*/
export function normalize2D(x: Vector2): number;

/**
*
* @param {Number[]} args
*/
export function determinant2x2(args: number[]): number;

/**
* Fill a 4x4 matrix with the given row vectors
* @param {Vector4} row0
* @param {Vector4} row1
* @param {Vector4} row2
* @param {Vector4} row3
* @param {Matrix} mat
*/
export function rowsToMat4(row0: Vector4, row1: Vector4, row2: Vector4, row3: Vector4, mat: Matrix): Matrix;

/**
* Fill a 4x4 matrix with the given column vectors
* @param {Vector4} column0
* @param {Vector4} column1
* @param {Vector4} column2
* @param {Vector4} column3
* @param {Matrix} mat
*/
export function columnsToMat4(column0: Vector4, column1: Vector4, column2: Vector4, column3: Vector4, mat: Matrix): Matrix;

/**
* Fill a 3x3 matrix with the given row vectors
* @param {Vector3} row0
* @param {Vector3} row1
* @param {Vector3} row2
* @param {Matrix} mat
*/
export function rowsToMat3(row0: Vector3, row1: Vector3, row2: Vector3, mat: Matrix): Matrix;

/**
* Fill a 3x3 matrix with the given column vectors
* @param {Vector3} column0
* @param {Vector3} column1
* @param {Vector3} column2
* @param {Matrix} mat
*/
export function columnsToMat3(column0: Vector3, column1: Vector3, column2: Vector3, mat: Matrix): Matrix;

/**
* LU Factorization of a 3x3 matrix.
* @param {Matrix3x3} mat_3x3
* @param {Vector3} index_3
*/
export function LUFactor3x3(mat_3x3: Matrix3x3, index_3: Vector3): void;

/**
* LU back substitution for a 3x3 matrix.
* @param {Matrix3x3} mat_3x3
* @param {Vector3} index_3
* @param {Vector3} x_3
*/
export function LUSolve3x3(mat_3x3: Matrix3x3, index_3: Vector3, x_3: Vector3): void;

/**
* Solve mat_3x3y_3 = x_3 for y and place the result in y.
* @param {Matrix3x3} mat_3x3
* @param {Vector3} x_3
* @param {Vector3} y_3
*/
export function linearSolve3x3(mat_3x3: Matrix3x3, x_3: Vector3, y_3: Vector3): void;

/**
*
* @param {Matrix3x3} mat_3x3
* @param {Vector3} in_3
* @param {Vector3} out_3
*/
export function multiply3x3_vect3(mat_3x3: Matrix3x3, in_3: Vector3, out_3: Vector3): void;

/**
*
* @param {Matrix3x3} a_3x3
* @param {Matrix3x3} b_3x3
* @param {Matrix3x3} out_3x3
*/
export function multiply3x3_mat3(a_3x3: Matrix3x3, b_3x3: Matrix3x3, out_3x3: Matrix3x3): void;

/**
* Multiply two matrices.
* @param {Matrix} a
* @param {Matrix} b
* @param {Number} rowA
* @param {Number} colA
* @param {Number} rowB
* @param {Number} colB
* @param {Matrix} outRowAColB
*/
export function multiplyMatrix(a: Matrix, b: Matrix, rowA: number, colA: number, rowB: number, colB: number, outRowAColB: Matrix): void;

/**
* Transpose a 3x3 matrix.
* @param {Matrix3x3} in_3x3 The input 3x3 matrix.
* @param {Matrix3x3} outT_3x3 The output 3x3 matrix.
*/
export function transpose3x3(in_3x3: Matrix3x3, outT_3x3: Matrix3x3): void;

/**
* Invert a 3x3 matrix.
* @param {Matrix3x3} in_3x3 The input 3x3 matrix.
* @param {Matrix3x3} outI_3x3 The output 3x3 matrix.
*/
export function invert3x3(in_3x3: Matrix3x3, outI_3x3: Matrix3x3): void;

/**
* Set mat to the identity matrix.
* @param {Number} n The size of the matrix.
* @param {Number[]} mat The output matrix.
* @see isIdentity()
* @see identity3x3()
*/
export function identity(n: number, mat: number[]): void;

/**
* Set mat_3x3 to the identity matrix.
* @param {Matrix3x3} mat_3x3 The input 3x3 matrix.
* @see isIdentity3x3()
* @see identity()
*/
export function identity3x3(mat_3x3: Matrix3x3): void;

/**
* Returns true if provided matrix is the identity matrix.
* @param {Number[]} mat The 3x3 matrix to check
* @param {Number} [eps] The tolerance value.
* @see isIdentity()
* @see identity()
*/
export function isIdentity(mat: Matrix3x3, eps?: number): boolean;

/**
* Returns true if provided 3x3 matrix is the identity matrix.
* @param {Matrix3x3} mat The 3x3 matrix to check
* @param {Number} [eps] The tolerance value.
* @see isIdentity()
* @see identity3x3()
*/
export function isIdentity3x3(mat: Matrix3x3, eps?: number): boolean;

/**
* Calculate the determinant of a 3x3 matrix.
* @param {Matrix3x3} mat_3x3 The input 3x3 matrix.
*/
export function determinant3x3(mat_3x3: Matrix3x3): number;

/**
*
* @param {Vector4} quat_4
* @param {Matrix3x3} mat_3x3
*/
export function quaternionToMatrix3x3(quat_4: Vector4, mat_3x3: Matrix3x3): void;

/**
* Returns true if elements of both arrays are equals.
* @param {Number[]} a An array of numbers (vector, point, matrix...)
* @param {Number[]} b An array of numbers (vector, point, matrix...)
* @param {Number} [eps] The tolerance value.
*/
export function areEquals(a: number[], b: number[], eps?: number): boolean;

/**
*
* @param {Number} num
* @param {Number} [digits]
*/
export function roundNumber(num: number, digits?: number): number;

/**
*
* @param {Vector3} vector
* @param {Vector3} [out]
* @param {Number} [digits]
*/
export function roundVector(vector: Vector3, out?: Vector3, digits?: number): Vector3;

/**
* Jacobi iteration for the solution of eigenvectors/eigenvalues. Input matrix a is modified (the upper triangle is filled with zeros)
* @param {Matrix} a real symetric nxn matrix
* @param {Number} n matrix size
* @param {Number[]} w vector of size n to store eigenvalues (stored in decreasing order)
* @param {Number[]} v matrix of size nxn to store eigenvectors (stored in decreasing order, normalized)
*/
export function jacobiN(a: Matrix, n: number, w: number[], v: number[]): number;

/**
*
* @param {Matrix3x3} mat_3x3
* @param {Vector4} quat_4
*/
export function matrix3x3ToQuaternion(mat_3x3: Matrix3x3, quat_4: Vector4): void;

/**
*
* @param {Vector4} quat_1
* @param {Vector4} quat_2
* @param {Vector4} quat_out
*/
export function multiplyQuaternion(quat_1: Vector4, quat_2: Vector4, quat_out: Vector4): void;

/**
*
* @param {Matrix3x3} a_3x3
* @param {Matrix3x3} out_3x3
*/
export function orthogonalize3x3(a_3x3: Matrix3x3, out_3x3: Matrix3x3): void;

/**
*
* @param {Matrix3x3} a_3x3
* @param {Vector3} w_3
* @param {Matrix3x3} v_3x3
*/
export function diagonalize3x3(a_3x3: Matrix3x3, w_3: Vector3, v_3x3: Matrix3x3): void;

/**
*
* @param {Matrix3x3} a_3x3
* @param {Matrix3x3} u_3x3
* @param {Vector3} w_3
* @param {Matrix3x3} vT_3x3
*/
export function singularValueDecomposition3x3(a_3x3: Matrix3x3, u_3x3: Matrix3x3, w_3: Vector3, vT_3x3: Matrix3x3): void;

/**
*
* @param {Matrix} A
* @param {Number[]} index
* @param {Number} size
*/
export function luFactorLinearSystem(A: Matrix, index: number[], size: number): number;

/**
*
* @param {Matrix} A
* @param {Number[]} index
* @param {Number[]} x
* @param {Number} size
*/
export function luSolveLinearSystem(A: Matrix, index: number[], x: number[], size: number): void;

/**
*
* @param {Matrix} A
* @param {Number[]} x
* @param {Number} size
*/
export function solveLinearSystem(A: Matrix, x: number[], size: number): number;

/**
*
* @param {Matrix} A The input matrix. It is modified during the inversion.
* @param {Matrix} AI The output inverse matrix. Can be the same as input matrix.
* @param {Number} [size] The square size of the matrix to invert : 4 for a 4x4
* @param {Number[]} [index]
* @param {Number[]} [column]
* @return AI on success, null otherwise
*/
export function invertMatrix(A: Matrix, AI: Matrix, size?: number, index?: number[], column?: number[]): Matrix|null;

/**
*
* @param {Matrix} A
* @param {Number} size
*/
export function estimateMatrixCondition(A: Matrix, size: number): number;

/**
*
* @param {Matrix3x3} a_3x3
* @param {Number[]} w
* @param {Number[]} v
*/
export function jacobi(a_3x3: Matrix3x3, w: number[], v: number[]): number;

/**
* Solves for the least squares best fit matrix for the homogeneous equation
* X'M' = 0'. Uses the method described on pages 40-41 of Computer Vision by
* Forsyth and Ponce, which is that the solution is the eigenvector associated
* with the minimum eigenvalue of T(X)X, where T(X) is the transpose of X. The
* inputs and output are transposed matrices. Dimensions: X' is numberOfSamples
* by xOrder, M' dimension is xOrder by yOrder. M' should be pre-allocated. All
* matrices are row major. The resultant matrix M' should be pre-multiplied to
* X' to get 0', or transposed and then post multiplied to X to get 0
* @param {Number} numberOfSamples
* @param {Matrix} xt
* @param {Number} xOrder
* @param {Matrix} mt
*/
export function solveHomogeneousLeastSquares(numberOfSamples: number, xt: Matrix, xOrder: number, mt: Matrix): number;

/**
* Solves for the least squares best fit matrix for the equation X'M' = Y'. Uses
* pseudoinverse to get the ordinary least squares. The inputs and output are
* transposed matrices. Dimensions: X' is numberOfSamples by xOrder, Y' is
* numberOfSamples by yOrder, M' dimension is xOrder by yOrder. M' should be
* pre-allocated. All matrices are row major. The resultant matrix M' should be
* pre-multiplied to X' to get Y', or transposed and then post multiplied to X
* to get Y By default, this method checks for the homogeneous condition where
* Y==0, and if so, invokes SolveHomogeneousLeastSquares. For better performance
* when the system is known not to be homogeneous, invoke with
* checkHomogeneous=0.
* @param {Number} numberOfSamples
* @param {Matrix} xt
* @param {Number} xOrder
* @param {Matrix} yt
* @param {Number} yOrder
* @param {Matrix} mt
* @param {Boolean} [checkHomogeneous]
*/
export function solveLeastSquares(numberOfSamples: number, xt: Matrix, xOrder: number, yt: Matrix, yOrder: number, mt: Matrix, checkHomogeneous?: boolean): number;

/**
*
* @param {String} hexStr
* @param {Number[]} [outFloatArray]
*/
export function hex2float(hexStr: string, outFloatArray?: number[]): number[];

/**
*
* @param {RGBColor} rgb An Array of the RGB color.
* @param {HSVColor} hsv An Array of the HSV color.
*/
export function rgb2hsv(rgb: RGBColor, hsv: HSVColor): void;

/**
*
* @param {HSVColor} hsv An Array of the HSV color.
* @param {RGBColor} rgb An Array of the RGB color.
*/
export function hsv2rgb(hsv: HSVColor, rgb: RGBColor): void;

/**
*
* @param {Vector3} lab
* @param {Vector3} xyz
*/
export function lab2xyz(lab: Vector3, xyz: Vector3): void;

/**
*
* @param {Vector3} xyz
* @param {Vector3} lab
*/
export function xyz2lab(xyz: Vector3, lab: Vector3): void;

/**
*
* @param {Vector3} xyz
* @param {RGBColor} rgb An Array of the RGB color.
*/
export function xyz2rgb(xyz: Vector3, rgb: RGBColor): void;

/**
*
* @param {RGBColor} rgb An Array of the RGB color.
* @param {Vector3} xyz
*/
export function rgb2xyz(rgb: RGBColor, xyz: Vector3): void;

/**
*
* @param {RGBColor} rgb
* @param {Vector3} lab
*/
export function rgb2lab(rgb: RGBColor, lab: Vector3): void;

/**
*
* @param {Vector3} lab
* @param {RGBColor} rgb An Array of the RGB color.
*/
export function lab2rgb(lab: Vector3, rgb: RGBColor): void;

/**
* Returns bounds.
* @param {Bounds} bounds Output array that hold bounds, optionally empty.
*/
export function uninitializeBounds(bounds: Bounds): Bounds;

/**
*
* @param {Bounds} bounds The bounds to check.
*/
export function areBoundsInitialized(bounds: Bounds): boolean;

/**
* Compute the bounds from points.
* @param {Vector3} point1 The coordinate of the first point.
* @param {Vector3} point2 The coordinate of the second point.
* @param {Bounds} bounds Output array that hold bounds, optionally empty.
* @deprecated please use vtkBoundingBox.addPoints(vtkBoundingBox.reset([]), points)
*/
export function computeBoundsFromPoints(point1: Vector3, point2: Vector3, bounds: Bounds): Bounds;

/**
* Clamp some value against a range.
* @param {Number} value The value to clamp.
* @param {Number} minValue The minimum value.
* @param {Number} maxValue The maximum value.
*/
export function clampValue(value: number, minValue: number, maxValue: number): number;

/**
* Clamp some vector against a range.
* @param {Vector3} vector The vector to clamp.
* @param {Vector3} minVector The minimum vector.
* @param {Vector3} maxVector The maximum vector.
* @param {Vector3} out The output vector.
*/
export function clampVector(vector: Vector3, minVector: Vector3, maxVector: Vector3, out: Vector3): Vector3;

/**
*
* @param {Vector3} vector
* @param {Vector3} out
*/
export function roundVector(vector: Vector3, out: Vector3): Vector3;

/**
*
* @param {Number} value
* @param {Range} range
*/
export function clampAndNormalizeValue(value: number, range: Range): number;

/**
* Get the scalar type that is most likely to have enough precision to store a
* given range of data once it has been scaled and shifted
*/
export function getScalarTypeFittingRange(): void;

/**
*
*/
export function getAdjustedScalarRange(): void;

/**
* Check if first 3D extent is within second 3D extent.
* @param {Extent} extent1 The first extent.
* @param {Extent} extent2 The second extent.
*/
export function extentIsWithinOtherExtent(extent1: Extent, extent2: Extent): number;

/**
* Check if first 3D bounds is within the second 3D bounds.
* @param {Bounds} bounds1_6 The first bounds.
* @param {Bounds} bounds2_6 The second bounds.
* @param {Vector3} delta_3 The error margin along each axis.
*/
export function boundsIsWithinOtherBounds(bounds1_6: Bounds, bounds2_6: Bounds, delta_3: Vector3): number;

/**
* Check if point is within the given 3D bounds.
* @param {Vector3} point_3 The coordinate of the point.
* @param {Bounds} bounds_6 The bounds.
* @param {Vector3} delta_3 The error margin along each axis.
*/
export function pointIsWithinBounds(point_3: Vector3, bounds_6: Bounds, delta_3: Vector3): number;

/**
* In Euclidean space, there is a unique circle passing through any given three
* non-collinear points P1, P2, and P3.
*
* Using Cartesian coordinates to represent these points as spatial vectors, it
* is possible to use the dot product and cross product to calculate the radius
* and center of the circle. See:
* http://en.wikipedia.org/wiki/Circumscribed_circle and more specifically the
* section Barycentric coordinates from cross- and dot-products
* @param {Vector3} p1 The coordinate of the first point.
* @param {Vector3} p2 The coordinate of the second point.
* @param {Vector3} p3 The coordinate of the third point.
* @param {Vector3} center The coordinate of the center point.
*/
export function solve3PointCircle(p1: Vector3, p2: Vector3, p3: Vector3, center: Vector3): number;

/**
* Determines whether the passed value is a infinite number.
* @param {Number} value The value to check.
*/
export function isInf(value: number): boolean;

/**
*
*/
export function createUninitializedBounds(): Bounds;

/**
*
* @param {Number[]} vector
*/
export function getMajorAxisIndex(vector: number[]): number;

/**
* Return the closest orthogonal matrix of 1, -1 and 0
* It works for both column major and row major matrices
* This function iteratively associate a column with a row by choosing
* the greatest absolute value from the remaining row and columns
* For each association, a -1 or a 1 is set in the output, depending on
* the sign of the value in the original matrix
*
* @param {Number[]} matrix The matrix of size nxn
* @param {Number[]} n The size of the square matrix, defaults to 3
*/
export function getSparseOrthogonalMatrix(matrix: number[], n: number): number[];

/**
*
* @param {Number} value The value to convert.
*/
export function floatToHex2(value: number): string;

/**
*
* @param {RGBColor} rgbArray
* @param {string} [prefix]
*/
export function floatRGB2HexCode(rgbArray: RGBColor | RGBAColor, prefix?: string): string;

/**
* Convert RGB or RGBA color array to CSS representation
* @param {RGBColor|RGBAColor} rgbArray The color array.
*/
export function float2CssRGBA(rgbArray: RGBColor | RGBAColor): string;

/**
* Determines whether the passed value is a NaN.
* @param {Number} value The value to check.
*/
export function isNan(value: number): boolean;

/**
* Determines whether the passed value is a NaN.
* @param {Number} value The value to check.
*/
export function isNaN(value: number): boolean;

/**
* Determines whether the passed value is a finite number.
* @param value The value to check.
*/
export function isFinite(value: any): boolean;

/**
* vtkMath provides methods to perform common math operations. These include
* providing constants such as Pi; conversion from degrees to radians; vector
* operations such as dot and cross products and vector norm; matrix determinant
* for 2x2 and 3x3 matrices; univariate polynomial solvers; and for random
* number generation (for backward compatibility only).
*/
export declare const vtkMath: {
createArray: typeof createArray;
swapRowsMatrix_nxn: typeof swapRowsMatrix_nxn;
swapColumnsMatrix_nxn: typeof swapColumnsMatrix_nxn;
Pi: typeof Pi;
radiansFromDegrees: typeof radiansFromDegrees;
degreesFromRadians: typeof degreesFromRadians;
round: typeof round;
floor: typeof floor;
ceil: typeof ceil;
min: typeof min;
max: typeof max;
arrayMin: typeof arrayMin;
arrayMax: typeof arrayMax;
arrayRange: typeof arrayRange;
ceilLog2: typeof ceilLog2;
factorial: typeof factorial;
gaussian: typeof gaussian;
nearestPowerOfTwo: typeof nearestPowerOfTwo;
isPowerOfTwo: typeof isPowerOfTwo;
binomial: typeof binomial;
beginCombination: typeof beginCombination;
nextCombination: typeof nextCombination;
randomSeed: typeof randomSeed;
getSeed: typeof getSeed;
random: typeof random;
add: typeof add;
subtract: typeof subtract;
multiplyScalar: typeof multiplyScalar;
multiplyScalar2D: typeof multiplyScalar2D;
multiplyAccumulate: typeof multiplyAccumulate;
multiplyAccumulate2D: typeof multiplyAccumulate2D;
dot: typeof dot;
outer: typeof outer;
cross: typeof cross;
norm: typeof norm;
normalize: typeof normalize;
perpendiculars: typeof perpendiculars;
projectVector: typeof projectVector;
dot2D: typeof dot2D;
projectVector2D: typeof projectVector2D;
distance2BetweenPoints: typeof distance2BetweenPoints;
angleBetweenVectors: typeof angleBetweenVectors;
gaussianAmplitude: typeof gaussianAmplitude;
gaussianWeight: typeof gaussianWeight;
outer2D: typeof outer2D;
norm2D: typeof norm2D;
normalize2D: typeof normalize2D;
determinant2x2: typeof determinant2x2;
rowsToMat4: typeof rowsToMat4;
columnsToMat4: typeof columnsToMat4;
rowsToMat3: typeof rowsToMat3;
columnsToMat3: typeof columnsToMat3;
LUFactor3x3: typeof LUFactor3x3;
LUSolve3x3: typeof LUSolve3x3;
linearSolve3x3: typeof linearSolve3x3;
multiply3x3_vect3: typeof multiply3x3_vect3;
multiply3x3_mat3: typeof multiply3x3_mat3;
multiplyMatrix: typeof multiplyMatrix;
transpose3x3: typeof transpose3x3;
invert3x3: typeof invert3x3;
identity3x3: typeof identity3x3;
determinant3x3: typeof determinant3x3;
quaternionToMatrix3x3: typeof quaternionToMatrix3x3;
areEquals: typeof areEquals;
areMatricesEqual: typeof areEquals;
roundNumber: typeof roundNumber;
roundVector: typeof roundVector;
jacobiN: typeof jacobiN;
matrix3x3ToQuaternion: typeof matrix3x3ToQuaternion;
multiplyQuaternion: typeof multiplyQuaternion;
orthogonalize3x3: typeof orthogonalize3x3;
diagonalize3x3: typeof diagonalize3x3;
singularValueDecomposition3x3: typeof singularValueDecomposition3x3;
luFactorLinearSystem: typeof luFactorLinearSystem;
luSolveLinearSystem: typeof luSolveLinearSystem;
solveLinearSystem: typeof solveLinearSystem;
invertMatrix: typeof invertMatrix;
estimateMatrixCondition: typeof estimateMatrixCondition;
jacobi: typeof jacobi;
solveHomogeneousLeastSquares: typeof solveHomogeneousLeastSquares;
solveLeastSquares: typeof solveLeastSquares;
hex2float: typeof hex2float;
rgb2hsv: typeof rgb2hsv;
hsv2rgb: typeof hsv2rgb;
lab2xyz: typeof lab2xyz;
xyz2lab: typeof xyz2lab;
xyz2rgb: typeof xyz2rgb;
rgb2xyz: typeof rgb2xyz;
rgb2lab: typeof rgb2lab;
lab2rgb: typeof lab2rgb;
uninitializeBounds: typeof uninitializeBounds;
areBoundsInitialized: typeof areBoundsInitialized;
computeBoundsFromPoints: typeof computeBoundsFromPoints;
clampValue: typeof clampValue;
clampVector: typeof clampVector;
clampAndNormalizeValue: typeof clampAndNormalizeValue;
getScalarTypeFittingRange: typeof getScalarTypeFittingRange;
getAdjustedScalarRange: typeof getAdjustedScalarRange;
extentIsWithinOtherExtent: typeof extentIsWithinOtherExtent;
boundsIsWithinOtherBounds: typeof boundsIsWithinOtherBounds;
pointIsWithinBounds: typeof pointIsWithinBounds;
solve3PointCircle: typeof solve3PointCircle;
isInf: typeof isInf;
createUninitializedBounds: typeof createUninitializedBounds;
getMajorAxisIndex: typeof getMajorAxisIndex;
getSparseOrthogonalMatrix: typeof getSparseOrthogonalMatrix;
floatToHex2: typeof floatToHex2;
floatRGB2HexCode: typeof floatRGB2HexCode;
float2CssRGBA: typeof float2CssRGBA;
inf: number;
negInf: number;
isNan: typeof isNaN,
isNaN: typeof isNaN;
isFinite: typeof isFinite
}
export default vtkMath;
index.js
import seedrandom from 'seedrandom';
import macro from 'vtk.js/Sources/macros';
import {
IDENTITY,
IDENTITY_3X3,
EPSILON,
VTK_SMALL_NUMBER,
} from 'vtk.js/Sources/Common/Core/Math/Constants';

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;

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

// Swap rows for n by n matrix
function swapRowsMatrix_nxn(matrix, n, row1, row2) {
let tmp;
for (let i = 0; i < n; i++) {
tmp = matrix[row1 * n + i];
matrix[row1 * n + i] = matrix[row2 * n + i];
matrix[row2 * n + i] = tmp;
}
}

// Swap columns for n by n matrix
function swapColumnsMatrix_nxn(matrix, n, column1, column2) {
let tmp;
for (let i = 0; i < n; i++) {
tmp = matrix[i * n + column1];
matrix[i * n + column1] = matrix[i * n + column2];
matrix[i * n + column2] = tmp;
}
}

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

export function createArray(size = 3) {
// faster than Array.from and/or while loop
const res = Array(size);
for (let i = 0; i < size; ++i) {
res[i] = 0;
}
return res;
}

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];
return out;
}

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

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

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

export function multiplyAccumulate(a, b, scalar, out) {
out[0] = a[0] + b[0] * scalar;
out[1] = a[1] + b[1] * scalar;
out[2] = a[2] + b[2] * scalar;
return out;
}

export function multiplyAccumulate2D(a, b, scalar, out) {
out[0] = a[0] + b[0] * scalar;
out[1] = a[1] + b[1] * scalar;
return out;
}

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) {
out_3x3[0] = x[0] * y[0];
out_3x3[1] = x[0] * y[1];
out_3x3[2] = x[0] * y[2];
out_3x3[3] = x[1] * y[0];
out_3x3[4] = x[1] * y[1];
out_3x3[5] = x[1] * y[2];
out_3x3[6] = x[2] * y[0];
out_3x3[7] = x[2] * y[1];
out_3x3[8] = x[2] * y[2];
}

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;
return out;
}

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 signedAngleBetweenVectors(v1, v2, vN) {
const crossVect = [0, 0, 0];
cross(v1, v2, crossVect);
const angle = Math.atan2(norm(crossVect), dot(v1, v2));
return dot(crossVect, vN) >= 0 ? angle : -angle;
}

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) {
out_2x2[0] = x[0] * y[0];
out_2x2[1] = x[0] * y[1];
out_2x2[2] = x[1] * y[0];
out_2x2[3] = x[1] * y[1];
}

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 rowsToMat4(row0, row1, row2, row3, mat) {
for (let i = 0; i < 4; i++) {
mat[i] = row0[i];
mat[4 + i] = row1[i];
mat[8 + i] = row2[i];
mat[12 + i] = row3[i];
}
return mat;
}

export function columnsToMat4(column0, column1, column2, column3, mat) {
for (let i = 0; i < 4; i++) {
mat[4 * i] = column0[i];
mat[4 * i + 1] = column1[i];
mat[4 * i + 2] = column2[i];
mat[4 * i + 3] = column3[i];
}
return mat;
}

export function rowsToMat3(row0, row1, row2, mat) {
for (let i = 0; i < 3; i++) {
mat[i] = row0[i];
mat[3 + i] = row1[i];
mat[6 + i] = row2[i];
}
return mat;
}

export function columnsToMat3(column0, column1, column2, mat) {
for (let i = 0; i < 3; i++) {
mat[3 * i] = column0[i];
mat[3 * i + 1] = column1[i];
mat[3 * i + 2] = column2[i];
}
return mat;
}

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 * 3]);
if ((tmp = Math.abs(mat_3x3[i * 3 + 1])) > largest) {
largest = tmp;
}
if ((tmp = Math.abs(mat_3x3[i * 3 + 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]);
maxI = 0;
if ((tmp = scale[1] * Math.abs(mat_3x3[3])) >= largest) {
largest = tmp;
maxI = 1;
}
if ((tmp = scale[2] * Math.abs(mat_3x3[6])) >= largest) {
maxI = 2;
}
if (maxI !== 0) {
swapRowsMatrix_nxn(mat_3x3, 3, maxI, 0);
scale[maxI] = scale[0];
}
index_3[0] = maxI;

mat_3x3[3] /= mat_3x3[0];
mat_3x3[6] /= mat_3x3[0];

// second column
mat_3x3[4] -= mat_3x3[3] * mat_3x3[1];
mat_3x3[7] -= mat_3x3[6] * mat_3x3[1];
largest = scale[1] * Math.abs(mat_3x3[4]);
maxI = 1;
if ((tmp = scale[2] * Math.abs(mat_3x3[7])) >= largest) {
maxI = 2;
swapRowsMatrix_nxn(mat_3x3, 3, 1, 2);
scale[2] = scale[1];
}
index_3[1] = maxI;
mat_3x3[7] /= mat_3x3[4];

// third column
mat_3x3[5] -= mat_3x3[3] * mat_3x3[2];
mat_3x3[8] -= mat_3x3[6] * mat_3x3[2] + mat_3x3[7] * mat_3x3[5];
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[3] * x_3[0];

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

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

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

// 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] * in_3[0] + mat_3x3[1] * in_3[1] + mat_3x3[2] * in_3[2];
const y = mat_3x3[3] * in_3[0] + mat_3x3[4] * in_3[1] + mat_3x3[5] * in_3[2];
const z = mat_3x3[6] * in_3[0] + mat_3x3[7] * in_3[1] + mat_3x3[8] * 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 copyA = [...a_3x3];
const copyB = [...b_3x3];
for (let i = 0; i < 3; i++) {
out_3x3[i] =
copyA[0] * copyB[i] + copyA[1] * copyB[i + 3] + copyA[2] * copyB[i + 6];
out_3x3[i + 3] =
copyA[3] * copyB[i] + copyA[4] * copyB[i + 3] + copyA[5] * copyB[i + 6];
out_3x3[i + 6] =
copyA[6] * copyB[i] + copyA[7] * copyB[i + 3] + copyA[8] * copyB[i + 6];
}
}

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

// If a or b is used to store the result, copying them is required
const copyA = [...a];
const copyB = [...b];
// output matrix is rowA*colB
// output row
for (let i = 0; i < rowA; i++) {
// output col
for (let j = 0; j < colB; j++) {
outRowAColB[i * colB + j] = 0;
// sum for this point
for (let k = 0; k < colA; k++) {
outRowAColB[i * colB + j] += copyA[i * colA + k] * copyB[j + colB * k];
}
}
}
}

export function transpose3x3(in_3x3, outT_3x3) {
let tmp;

// off-diagonal elements
tmp = in_3x3[3];
outT_3x3[3] = in_3x3[1];
outT_3x3[1] = tmp;
tmp = in_3x3[6];
outT_3x3[6] = in_3x3[2];
outT_3x3[2] = tmp;
tmp = in_3x3[7];
outT_3x3[7] = in_3x3[5];
outT_3x3[5] = tmp;

// on-diagonal elements
outT_3x3[0] = in_3x3[0];
outT_3x3[4] = in_3x3[4];
outT_3x3[8] = in_3x3[8];
}

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

// 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;
if (det === 0) {
vtkWarningMacro('Matrix has 0 determinant');
}

outI_3x3[0] = d1 / det;
outI_3x3[3] = d2 / det;
outI_3x3[6] = d3 / det;

outI_3x3[1] = e1 / det;
outI_3x3[4] = e2 / det;
outI_3x3[7] = e3 / det;

outI_3x3[2] = f1 / det;
outI_3x3[5] = f2 / det;
outI_3x3[8] = f3 / det;
}

export function determinant3x3(mat_3x3) {
return (
mat_3x3[0] * mat_3x3[4] * mat_3x3[8] +
mat_3x3[3] * mat_3x3[7] * mat_3x3[2] +
mat_3x3[6] * mat_3x3[1] * mat_3x3[5] -
mat_3x3[0] * mat_3x3[7] * mat_3x3[5] -
mat_3x3[3] * mat_3x3[1] * mat_3x3[8] -
mat_3x3[6] * mat_3x3[4] * mat_3x3[2]
);
}

/**
* Returns true if elements of both arrays are equals.
* @param {Array} a an array of numbers (vector, point, matrix...)
* @param {Array} b an array of numbers (vector, point, matrix...)
* @param {Number} eps tolerance
*/
export function areEquals(a, b, eps = EPSILON) {
if (a.length !== b.length) {
return false;
}

function isEqual(element, index) {
return Math.abs(element - b[index]) <= eps;
}
return a.every(isEqual);
}

export const areMatricesEqual = areEquals;

export function identity3x3(mat_3x3) {
for (let i = 0; i < 3; i++) {
/* eslint-disable-next-line no-multi-assign */
mat_3x3[i * 3] = mat_3x3[i * 3 + 1] = mat_3x3[i * 3 + 2] = 0;
mat_3x3[i * 3 + i] = 1;
}
}

export function identity(n, mat) {
for (let i = 0; i < n; i++) {
for (let j = 0; j < n; j++) {
mat[i * n + j] = 0;
}
mat[i * n + i] = 1;
}
return mat;
}

export function isIdentity(mat, eps = EPSILON) {
return areMatricesEqual(mat, IDENTITY, eps);
}

export function isIdentity3x3(mat, eps = EPSILON) {
return areMatricesEqual(mat, IDENTITY_3X3, eps);
}

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] = xx * f + s;
mat_3x3[3] = (xy + wz) * f;
mat_3x3[6] = (xz - wy) * f;

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

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

export function roundNumber(num, digits = 0) {
if (!`${num}`.includes('e')) {
return +`${Math.round(`${num}e+${digits}`)}e-${digits}`;
}
const arr = `${num}`.split('e');
let sig = '';
if (+arr[1] + digits > 0) {
sig = '+';
}
return +`${Math.round(`${+arr[0]}e${sig}${+arr[1] + digits}`)}e-${digits}`;
}

export function roundVector(vector, out = [0, 0, 0], digits = 0) {
out[0] = roundNumber(vector[0], digits);
out[1] = roundNumber(vector[1], digits);
out[2] = roundNumber(vector[2], digits);

return out;
}

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) => {
g = aa[ii];
h = aa[jj];
aa[ii] = g - s * (h + g * tau);
aa[jj] = h + s * (g - h * tau);
};

// initialize
identity(n, v);
for (ip = 0; ip < n; ip++) {
b[ip] = w[ip] = a[ip + ip * n];
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 * n + 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 * n + 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 * n + iq] = 0.0;
} else if (Math.abs(a[ip * n + iq]) > tresh) {
h = w[iq] - w[ip];
if (Math.abs(h) + g === Math.abs(h)) {
t = a[ip * n + iq] / h;
} else {
theta = (0.5 * h) / a[ip * n + 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 * n + iq];
z[ip] -= h;
z[iq] += h;
w[ip] -= h;
w[iq] += h;
a[ip * n + iq] = 0.0;

// ip already shifted left by 1 unit
for (j = 0; j <= ip - 1; j++) {
vtkROTATE(a, j * n + ip, j * n + iq);
}
// ip and iq already shifted left by 1 unit
for (j = ip + 1; j <= iq - 1; j++) {
vtkROTATE(a, ip * n + j, j * n + iq);
}
// iq already shifted left by 1 unit
for (j = iq + 1; j < n; j++) {
vtkROTATE(a, ip * n + j, iq * n + j);
}
for (j = 0; j < n; j++) {
vtkROTATE(v, j * n + ip, j * n + 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 || Math.abs(w[i] - tmp) < VTK_SMALL_NUMBER) {
// why exchange if same?
k = i;
tmp = w[k];
}
}
if (k !== j) {
w[k] = w[j];
w[j] = tmp;
swapColumnsMatrix_nxn(v, n, j, k);
}
}
// 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 (numPos = 0, i = 0; i < n * n; i++) {
if (v[i] >= 0.0) {
numPos++;
}
}
// if ( numPos < ceil(double(n)/double(2.0)) )
if (numPos < ceil_half_n) {
for (i = 0; i < n; i++) {
v[i * n + 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] = mat_3x3[0] + mat_3x3[4] + mat_3x3[8];
tmp[5] = mat_3x3[0] - mat_3x3[4] - mat_3x3[8];
tmp[10] = -mat_3x3[0] + mat_3x3[4] - mat_3x3[8];
tmp[15] = -mat_3x3[0] - mat_3x3[4] + mat_3x3[8];

// off-diagonal elements
tmp[1] = tmp[4] = mat_3x3[7] - mat_3x3[5];
tmp[2] = tmp[8] = mat_3x3[2] - mat_3x3[6];
tmp[3] = tmp[12] = mat_3x3[3] - mat_3x3[1];

tmp[6] = tmp[9] = mat_3x3[3] + mat_3x3[1];
tmp[7] = tmp[13] = mat_3x3[2] + mat_3x3[6];
tmp[11] = tmp[14] = mat_3x3[7] + mat_3x3[5];

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
// tmp is copied because jacobiN may modify it
const NTemp = [...tmp];
jacobiN(NTemp, 4, eigenvalues, eigenvectors);

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

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 < 9; i++) {
out_3x3[i] = a_3x3[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 * 3]);
const x2 = Math.abs(out_3x3[i * 3 + 1]);
const x3 = Math.abs(out_3x3[i * 3 + 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]) * scale[0];
const x2 = Math.abs(out_3x3[3]) * scale[1];
const x3 = Math.abs(out_3x3[6]) * 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) {
// swap vectors
swapColumnsMatrix_nxn(out_3x3, 3, index[0], 0);
scale[index[0]] = scale[0];
}

// second column
const y2 = Math.abs(out_3x3[4]) * scale[1];
const y3 = Math.abs(out_3x3[7]) * scale[2];
index[1] = 1;
largest = y2;
if (y3 >= largest) {
index[1] = 2;
// swap vectors
swapColumnsMatrix_nxn(out_3x3, 3, 1, 2);
}

// 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 < 9; i++) {
out_3x3[i] = -out_3x3[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 < 9; i++) {
out_3x3[i] = -out_3x3[i];
}
}

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

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

// a is copied because jacobiN may modify it
const copyA = [...a_3x3];

// diagonalize using Jacobi
jacobiN(copyA, 3, w_3, v_3x3);

// 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 * 3]);
maxI = 0;
for (j = 1; j < 3; j++) {
if (maxVal < (tmp = Math.abs(v_3x3[i * 3 + 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;
swapRowsMatrix_nxn(v_3x3, 3, i, maxI);
}
// maximum element of eigenvector should be positive
if (v_3x3[maxI * 3 + maxI] < 0) {
v_3x3[maxI * 3] = -v_3x3[maxI * 3];
v_3x3[maxI * 3 + 1] = -v_3x3[maxI * 3 + 1];
v_3x3[maxI * 3 + 2] = -v_3x3[maxI * 3 + 2];
}

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

v_3x3[j * 3] = 0.0;
v_3x3[j * 3 + 1] = 0.0;
v_3x3[j * 3 + 2] = 0.0;
v_3x3[j * 3 + j] = 1.0;
const vectTmp1 = cross(
[v_3x3[maxI * 3], v_3x3[maxI * 3 + 1], v_3x3[maxI * 3 + 2]],
[v_3x3[j * 3], v_3x3[j * 3 + 1], v_3x3[j * 3 + 2]],
[]
);
normalize(vectTmp1);
const vectTmp2 = cross(
vectTmp1,
[v_3x3[maxI * 3], v_3x3[maxI * 3 + 1], v_3x3[maxI * 3 + 2]],
[]
);
for (let t = 0; t < 3; t++) {
v_3x3[k * 3 + t] = vectTmp1[t];
v_3x3[j * 3 + t] = vectTmp2[t];
}

// 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]);
maxI = 0;
for (i = 1; i < 3; i++) {
if (maxVal < (tmp = Math.abs(v_3x3[i * 3]))) {
maxVal = tmp;
maxI = i;
}
}
// swap eigenvalue and eigenvector
if (maxI !== 0) {
const eigenValTmp = w_3[maxI];
w_3[maxI] = w_3[0];
w_3[0] = eigenValTmp;
swapRowsMatrix_nxn(v_3x3, 3, maxI, 0);
}
// do the same for the y element
if (Math.abs(v_3x3[4]) < Math.abs(v_3x3[7])) {
const eigenValTmp = w_3[2];
w_3[2] = w_3[1];
w_3[1] = eigenValTmp;
swapRowsMatrix_nxn(v_3x3, 3, 1, 2);
}

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

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

export function singularValueDecomposition3x3(a_3x3, u_3x3, w_3, vT_3x3) {
let i;
// copy so that A can be used for U or VT without risk
const B = [...a_3x3];

// temporarily flip if determinant is negative
const d = determinant3x3(B);
if (d < 0) {
for (i = 0; i < 9; i++) {
B[i] = -B[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];
}
}

/**
* Factor linear equations Ax = b using LU decomposition A = LU. Output factorization LU is in matrix A.
* @param {Matrix} A square matrix
* @param {Number} index integer array of pivot indices index[0->n-1]
* @param {Number} size matrix size
*/
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 * size + 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 * size + j];
for (k = 0; k < i; k++) {
sum -= A[i * size + k] * A[k * size + j];
}
A[i * size + j] = sum;
}
//
// Begin search for largest pivot element
//
for (largest = 0.0, i = j; i < size; i++) {
sum = A[i * size + j];
for (k = 0; k < j; k++) {
sum -= A[i * size + k] * A[k * size + j];
}
A[i * size + 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 * size + k];
A[maxI * size + k] = A[j * size + k];
A[j * size + k] = temp1;
}
scale[maxI] = scale[j];
}
//
// Divide by pivot element and perform elimination
//
index[j] = maxI;

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

if (j !== size - 1) {
temp1 = 1.0 / A[j * size + j];
for (i = j + 1; i < size; i++) {
A[i * size + 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 * size + 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 * size + j] * x[j];
}
x[i] = sum / A[i * size + 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], A[1], A[2], A[3]);

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

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

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

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

x[0] /= A[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;
}

// Note that A is modified during the inversion !
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 null;
}

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 * size + j] = tmp2Size[i];
}
}

return AI;
}

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 * size + j]) > maxValue) {
maxValue = Math.abs(A[i * size + j]);
}
}
}

// find the minimum diagonal value
for (let i = 0; i < size; i++) {
if (Math.abs(A[i * size + i]) < minValue) {
minValue = Math.abs(A[i * size + 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 * xOrder); // size x by x
// Allocate the array of eigenvalues and eigenvectors
const eigenvals = createArray(xOrder);
const eigenvecs = createArray(xOrder * xOrder);

// 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 * xOrder + j] += xt[k * xOrder + i] * xt[k * xOrder + 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 * xOrder + j] = XXt[j * xOrder + 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] = eigenvecs[i * xOrder + 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 * yOrder + 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 should be xOrder x yOrder, but since we are solving only the homogeneous part, here it is xOrder x 1
hmt = createArray(xOrder);

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

// set up intermediate variables
const XXt = createArray(xOrder * xOrder); // size x by x
const XXtI = createArray(xOrder * xOrder); // size x by x
const XYt = createArray(xOrder * yOrder); // size x by y

// 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 * xOrder + j] += xt[k * xOrder + i] * xt[k * xOrder + j];
}

// now calculate the XYt matrix
for (j = 0; j < yOrder; j++) {
XYt[i * yOrder + j] += xt[k * xOrder + i] * yt[k * yOrder + 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 * xOrder + j] = XXt[j * xOrder + 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 * yOrder + j] = 0.0;
for (k = 0; k < xOrder; k++) {
mt[i * yOrder + j] += XXtI[i * xOrder + k] * XYt[k * yOrder + 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 * yOrder + j] = hmt[i * yOrder];
}
}
}
}

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;
return bounds;
}

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

/**
* @deprecated please use vtkBoundingBox.addPoints(vtkBoundingBox.reset([]), points)
*/
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]);
return bounds;
}

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 = [0, 0, 0]) {
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 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;
}

// Return the closest orthogonal matrix of 1, -1 and 0
// It works for both column major and row major matrices
// This function iteratively associate a column with a row by choosing
// the greatest absolute value from the remaining row and columns
// For each association, a -1 or a 1 is set in the output, depending on
// the sign of the value in the original matrix
export function getSparseOrthogonalMatrix(matrix, n = 3) {
// Initialize rows and columns to available indices
const rows = new Array(n);
const cols = new Array(n);
for (let i = 0; i < n; ++i) {
rows[i] = i;
cols[i] = i;
}
// No need for the last iteration: i = 0
for (let i = n - 1; i > 0; i--) {
// Loop invariant:
// rows[0:i] and cols[0:i] contain the remaining rows and columns
// rows]i:n[ and cols]i:n[ contain the associations found (rows[k] is associated with cols[k])
let bestValue = -Infinity;
let bestRowI = 0;
let bestColI = 0;
for (let rowI = 0; rowI <= i; ++rowI) {
const row = rows[rowI];
for (let colI = 0; colI <= i; ++colI) {
const col = cols[colI];
const absVal = Math.abs(matrix[row + n * col]);
if (absVal > bestValue) {
bestValue = absVal;
bestRowI = rowI;
bestColI = colI;
}
}
}
// Found an association between rows[bestRowI] and cols[bestColI]
// Put both at the end of their array by swapping with i
[rows[i], rows[bestRowI]] = [rows[bestRowI], rows[i]];
[cols[i], cols[bestColI]] = [cols[bestColI], cols[i]];
}

// Convert row/column association to a matrix
const output = new Array(n * n).fill(0);
for (let i = 0; i < n; ++i) {
const matIdx = rows[i] + n * cols[i];
output[matIdx] = matrix[matIdx] < 0 ? -1 : 1;
}

return output;
}

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,
multiplyAccumulate,
multiplyAccumulate2D,
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,
identity,
isIdentity,
isIdentity3x3,
determinant3x3,
quaternionToMatrix3x3,
areEquals,
areMatricesEqual,
roundNumber,
roundVector,
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,
getSparseOrthogonalMatrix,
floatToHex2,
floatRGB2HexCode,
float2CssRGBA,
};