diff --git a/src/mol-math/linear-algebra/matrix/matrix.ts b/src/mol-math/linear-algebra/matrix/matrix.ts
new file mode 100644
index 0000000000000000000000000000000000000000..549635053c5ce1ff03ecd97490d6cba7a76ea9ee
--- /dev/null
+++ b/src/mol-math/linear-algebra/matrix/matrix.ts
@@ -0,0 +1,75 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+interface Matrix { data: Helpers.NumberArray, size: number, cols: number, rows: number }
+
+namespace Matrix {
+ export function create(cols: number, rows: number, ctor: { new (size: number): Helpers.NumberArray } = Float32Array): Matrix {
+ const size = cols * rows
+ return { data: new ctor(size), size, cols, rows }
+ }
+
+ export function fromArray(data: Helpers.NumberArray, cols: number, rows: number): Matrix {
+ return { data, size: cols * rows, cols, rows }
+ }
+
+ export function transpose(out: Matrix, mat: Matrix): Matrix {
+ const nrows = mat.rows, ncols = mat.cols
+ const md = mat.data, mtd = out.data
+
+ for (let i = 0, mi = 0, mti = 0; i < nrows; mti += 1, mi += ncols, ++i) {
+ let ri = mti
+ for (let j = 0; j < ncols; ri += nrows, j++) mtd[ri] = md[mi + j]
+ }
+ return mat
+ }
+
+ /** out = matA * matB' */
+ export function multiplyABt (out: Matrix, matA: Matrix, matB: Matrix) {
+ const ncols = matA.cols, nrows = matA.rows, mrows = matB.rows
+ const ad = matA.data, bd = matB.data, cd = out.data
+
+ for (let i = 0, matAp = 0, outP = 0; i < nrows; matAp += ncols, i++) {
+ for (let pB = 0, j = 0; j < mrows; outP++, j++) {
+ let sum = 0.0
+ let pMatA = matAp
+ for (let k = 0; k < ncols; pMatA++, pB++, k++) {
+ sum += ad[pMatA] * bd[pB]
+ }
+ cd[outP] = sum
+ }
+ }
+ return out
+ }
+
+ /** Get the mean of rows in `mat` */
+ export function meanRows (mat: Matrix) {
+ const nrows = mat.rows, ncols = mat.cols
+ const md = mat.data
+ const mean = new Array(ncols)
+
+ for (let j = 0; j < ncols; ++j) mean[ j ] = 0.0
+ for (let i = 0, p = 0; i < nrows; ++i) {
+ for (let j = 0; j < ncols; ++j, ++p) mean[ j ] += md[ p ]
+ }
+ for (let j = 0; j < ncols; ++j) mean[ j ] /= nrows
+
+ return mean
+ }
+
+ /** Subtract `row` from all rows in `mat` */
+ export function subRows (mat: Matrix, row: Helpers.NumberArray) {
+ const nrows = mat.rows, ncols = mat.cols
+ const md = mat.data
+
+ for (let i = 0, p = 0; i < nrows; ++i) {
+ for (let j = 0; j < ncols; ++j, ++p) md[ p ] -= row[ j ]
+ }
+ return mat
+ }
+}
+
+export default Matrix
\ No newline at end of file
diff --git a/src/mol-math/linear-algebra/matrix/principal-axes.ts b/src/mol-math/linear-algebra/matrix/principal-axes.ts
new file mode 100644
index 0000000000000000000000000000000000000000..575668f62985a6fe95da8aa124729e5923bd7b1c
--- /dev/null
+++ b/src/mol-math/linear-algebra/matrix/principal-axes.ts
@@ -0,0 +1,191 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Matrix from './matrix.js';
+import { Vec3 } from '../3d.js';
+// import { Vec3, Mat4 } from '../3d.js';
+import { svd } from './svd.js';
+
+// const negateVector = Vec3.create(-1, -1, -1)
+// const tmpMatrix = Mat4.identity()
+
+/**
+ * Principal axes
+ */
+class PrincipalAxes {
+ begA: Vec3
+ endA: Vec3
+ begB: Vec3
+ endB: Vec3
+ begC: Vec3
+ endC: Vec3
+
+ center: Vec3
+
+ vecA: Vec3
+ vecB: Vec3
+ vecC: Vec3
+
+ normVecA: Vec3
+ normVecB: Vec3
+ normVecC: Vec3
+
+ /**
+ * points is a 3xN matrix
+ */
+ constructor(points: Matrix) {
+ const n = points.rows
+ const n3 = n / 3
+ const pointsT = Matrix.create(n, 3)
+ const A = Matrix.create(3, 3)
+ const W = Matrix.create(1, 3)
+ const U = Matrix.create(3, 3)
+ const V = Matrix.create(3, 3)
+
+ // calculate
+ const mean = Matrix.meanRows(points)
+ Matrix.subRows(points, mean)
+ Matrix.transpose(pointsT, points)
+ Matrix.multiplyABt(A, pointsT, pointsT)
+ svd(A, W, U, V)
+
+ // center
+ const vm = Vec3.create(mean[0], mean[1], mean[2])
+
+ // normalized
+ const van = Vec3.create(U.data[0], U.data[3], U.data[6])
+ const vbn = Vec3.create(U.data[1], U.data[4], U.data[7])
+ const vcn = Vec3.create(U.data[2], U.data[5], U.data[8])
+
+ // scaled
+ const va = Vec3.scale(Vec3.zero(), van, Math.sqrt(W.data[0] / n3))
+ const vb = Vec3.scale(Vec3.zero(), vbn, Math.sqrt(W.data[1] / n3))
+ const vc = Vec3.scale(Vec3.zero(), vcn, Math.sqrt(W.data[2] / n3))
+ // const va = van.clone().multiplyScalar(Math.sqrt(W.data[0] / n3))
+ // const vb = vbn.clone().multiplyScalar(Math.sqrt(W.data[1] / n3))
+ // const vc = vcn.clone().multiplyScalar(Math.sqrt(W.data[2] / n3))
+
+ // points
+ this.begA = Vec3.sub(Vec3.clone(vm), vm, va)
+ this.endA = Vec3.add(Vec3.clone(vm), vm, va)
+ this.begB = Vec3.sub(Vec3.clone(vm), vm, vb)
+ this.endB = Vec3.add(Vec3.clone(vm), vm, vb)
+ this.begC = Vec3.sub(Vec3.clone(vm), vm, vc)
+ this.endC = Vec3.add(Vec3.clone(vm), vm, vc)
+ // this.begA = vm.clone().sub(va)
+ // this.endA = vm.clone().add(va)
+ // this.begB = vm.clone().sub(vb)
+ // this.endB = vm.clone().add(vb)
+ // this.begC = vm.clone().sub(vc)
+ // this.endC = vm.clone().add(vc)
+
+ //
+
+ this.center = vm
+
+ this.vecA = va
+ this.vecB = vb
+ this.vecC = vc
+
+ this.normVecA = van
+ this.normVecB = vbn
+ this.normVecC = vcn
+ }
+
+ // TODO
+ // /**
+ // * Get the basis matrix descriping the axes
+ // * @param {Matrix4} [optionalTarget] - target object
+ // * @return {Matrix4} the basis
+ // */
+ // getBasisMatrix(optionalTarget = new Matrix4()) {
+ // const basis = optionalTarget
+
+ // basis.makeBasis(this.normVecB, this.normVecA, this.normVecC)
+ // if (basis.determinant() < 0) {
+ // basis.scale(negateVector)
+ // }
+
+ // return basis
+ // }
+
+ // TODO
+ // /**
+ // * Get a quaternion descriping the axes rotation
+ // * @param {Quaternion} [optionalTarget] - target object
+ // * @return {Quaternion} the rotation
+ // */
+ // getRotationQuaternion(optionalTarget = new Quaternion()) {
+ // const q = optionalTarget
+ // q.setFromRotationMatrix(this.getBasisMatrix(tmpMatrix))
+
+ // return q.inverse()
+ // }
+
+ // TODO
+ // /**
+ // * Get the scale/length for each dimension for a box around the axes
+ // * to enclose the atoms of a structure
+ // * @param {Structure|StructureView} structure - the structure
+ // * @return {{d1a: Number, d2a: Number, d3a: Number, d1b: Number, d2b: Number, d3b: Number}} scale
+ // */
+ // getProjectedScaleForAtoms(structure: Structure) {
+ // let d1a = -Infinity
+ // let d1b = -Infinity
+ // let d2a = -Infinity
+ // let d2b = -Infinity
+ // let d3a = -Infinity
+ // let d3b = -Infinity
+
+ // const p = new Vector3()
+ // const t = new Vector3()
+
+ // const center = this.center
+ // const ax1 = this.normVecA
+ // const ax2 = this.normVecB
+ // const ax3 = this.normVecC
+
+ // structure.eachAtom(function (ap: AtomProxy) {
+ // projectPointOnVector(p.copy(ap as any), ax1, center) // TODO
+ // const dp1 = t.subVectors(p, center).normalize().dot(ax1)
+ // const dt1 = p.distanceTo(center)
+ // if (dp1 > 0) {
+ // if (dt1 > d1a) d1a = dt1
+ // } else {
+ // if (dt1 > d1b) d1b = dt1
+ // }
+
+ // projectPointOnVector(p.copy(ap as any), ax2, center)
+ // const dp2 = t.subVectors(p, center).normalize().dot(ax2)
+ // const dt2 = p.distanceTo(center)
+ // if (dp2 > 0) {
+ // if (dt2 > d2a) d2a = dt2
+ // } else {
+ // if (dt2 > d2b) d2b = dt2
+ // }
+
+ // projectPointOnVector(p.copy(ap as any), ax3, center)
+ // const dp3 = t.subVectors(p, center).normalize().dot(ax3)
+ // const dt3 = p.distanceTo(center)
+ // if (dp3 > 0) {
+ // if (dt3 > d3a) d3a = dt3
+ // } else {
+ // if (dt3 > d3b) d3b = dt3
+ // }
+ // })
+
+ // return {
+ // d1a: d1a,
+ // d2a: d2a,
+ // d3a: d3a,
+ // d1b: -d1b,
+ // d2b: -d2b,
+ // d3b: -d3b
+ // }
+ // }
+}
+
+export default PrincipalAxes
\ No newline at end of file
diff --git a/src/mol-math/linear-algebra/matrix/svd.ts b/src/mol-math/linear-algebra/matrix/svd.ts
new file mode 100644
index 0000000000000000000000000000000000000000..e7ace3f959ae173d49c9869ce0015cce89c18b3f
--- /dev/null
+++ b/src/mol-math/linear-algebra/matrix/svd.ts
@@ -0,0 +1,292 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import Matrix from './matrix';
+
+// svd method adapted from http://inspirit.github.io/jsfeat/ MIT Eugene Zatepyakin
+
+export function swap(A: Helpers.NumberArray, i0: number, i1: number, t: number) {
+ t = A[i0]
+ A[i0] = A[i1]
+ A[i1] = t
+}
+
+export function hypot(a: number, b: number) {
+ a = Math.abs(a)
+ b = Math.abs(b)
+ if (a > b) {
+ b /= a
+ return a * Math.sqrt(1.0 + b * b)
+ }
+ if (b > 0) {
+ a /= b
+ return b * Math.sqrt(1.0 + a * a)
+ }
+ return 0.0
+}
+
+const EPSILON = 0.0000001192092896
+const FLT_MIN = 1E-37
+
+export function JacobiSVDImpl(At: Helpers.NumberArray, astep: number, _W: Helpers.NumberArray, Vt: Helpers.NumberArray, vstep: number, m: number, n: number, n1: number) {
+ const eps = EPSILON * 2.0
+ const minval = FLT_MIN
+ let i = 0
+ let j = 0
+ let k = 0
+ let iter = 0
+ const maxIter = Math.max(m, 30)
+ let Ai = 0
+ let Aj = 0
+ let Vi = 0
+ let Vj = 0
+ let changed = 0
+ let c = 0.0
+ let s = 0.0
+ let t = 0.0
+ let t0 = 0.0
+ let t1 = 0.0
+ let sd = 0.0
+ let beta = 0.0
+ let gamma = 0.0
+ let delta = 0.0
+ let a = 0.0
+ let p = 0.0
+ let b = 0.0
+ let seed = 0x1234
+ let val = 0.0
+ let val0 = 0.0
+ let asum = 0.0
+
+ const W = new Float64Array(n << 3)
+
+ for (; i < n; i++) {
+ for (k = 0, sd = 0; k < m; k++) {
+ t = At[i * astep + k]
+ sd += t * t
+ }
+ W[i] = sd
+
+ if (Vt) {
+ for (k = 0; k < n; k++) {
+ Vt[i * vstep + k] = 0
+ }
+ Vt[i * vstep + i] = 1
+ }
+ }
+
+ for (; iter < maxIter; iter++) {
+ changed = 0
+
+ for (i = 0; i < n - 1; i++) {
+ for (j = i + 1; j < n; j++) {
+ Ai = (i * astep) | 0
+ Aj = (j * astep) | 0
+ a = W[i]
+ p = 0
+ b = W[j]
+
+ k = 2
+ p += At[Ai] * At[Aj]
+ p += At[Ai + 1] * At[Aj + 1]
+
+ for (; k < m; k++) { p += At[Ai + k] * At[Aj + k] }
+
+ if (Math.abs(p) <= eps * Math.sqrt(a * b)) continue
+
+ p *= 2.0
+ beta = a - b
+ gamma = hypot(p, beta)
+ if (beta < 0) {
+ delta = (gamma - beta) * 0.5
+ s = Math.sqrt(delta / gamma)
+ c = (p / (gamma * s * 2.0))
+ } else {
+ c = Math.sqrt((gamma + beta) / (gamma * 2.0))
+ s = (p / (gamma * c * 2.0))
+ }
+
+ a = 0.0
+ b = 0.0
+
+ k = 2 // unroll
+ t0 = c * At[Ai] + s * At[Aj]
+ t1 = -s * At[Ai] + c * At[Aj]
+ At[Ai] = t0; At[Aj] = t1
+ a += t0 * t0; b += t1 * t1
+
+ t0 = c * At[Ai + 1] + s * At[Aj + 1]
+ t1 = -s * At[Ai + 1] + c * At[Aj + 1]
+ At[Ai + 1] = t0; At[Aj + 1] = t1
+ a += t0 * t0; b += t1 * t1
+
+ for (; k < m; k++) {
+ t0 = c * At[Ai + k] + s * At[Aj + k]
+ t1 = -s * At[Ai + k] + c * At[Aj + k]
+ At[Ai + k] = t0; At[Aj + k] = t1
+
+ a += t0 * t0; b += t1 * t1
+ }
+
+ W[i] = a
+ W[j] = b
+
+ changed = 1
+
+ if (Vt) {
+ Vi = (i * vstep) | 0
+ Vj = (j * vstep) | 0
+
+ k = 2
+ t0 = c * Vt[Vi] + s * Vt[Vj]
+ t1 = -s * Vt[Vi] + c * Vt[Vj]
+ Vt[Vi] = t0; Vt[Vj] = t1
+
+ t0 = c * Vt[Vi + 1] + s * Vt[Vj + 1]
+ t1 = -s * Vt[Vi + 1] + c * Vt[Vj + 1]
+ Vt[Vi + 1] = t0; Vt[Vj + 1] = t1
+
+ for (; k < n; k++) {
+ t0 = c * Vt[Vi + k] + s * Vt[Vj + k]
+ t1 = -s * Vt[Vi + k] + c * Vt[Vj + k]
+ Vt[Vi + k] = t0; Vt[Vj + k] = t1
+ }
+ }
+ }
+ }
+ if (changed === 0) break
+ }
+
+ for (i = 0; i < n; i++) {
+ for (k = 0, sd = 0; k < m; k++) {
+ t = At[i * astep + k]
+ sd += t * t
+ }
+ W[i] = Math.sqrt(sd)
+ }
+
+ for (i = 0; i < n - 1; i++) {
+ j = i
+ for (k = i + 1; k < n; k++) {
+ if (W[j] < W[k]) { j = k }
+ }
+ if (i !== j) {
+ swap(W, i, j, sd)
+ if (Vt) {
+ for (k = 0; k < m; k++) {
+ swap(At, i * astep + k, j * astep + k, t)
+ }
+
+ for (k = 0; k < n; k++) {
+ swap(Vt, i * vstep + k, j * vstep + k, t)
+ }
+ }
+ }
+ }
+
+ for (i = 0; i < n; i++) {
+ _W[i] = W[i]
+ }
+
+ if (!Vt) {
+ return
+ }
+
+ for (i = 0; i < n1; i++) {
+ sd = i < n ? W[i] : 0
+
+ while (sd <= minval) {
+ // if we got a zero singular value, then in order to get the corresponding left singular vector
+ // we generate a random vector, project it to the previously computed left singular vectors,
+ // subtract the projection and normalize the difference.
+ val0 = (1.0 / m)
+ for (k = 0; k < m; k++) {
+ seed = (seed * 214013 + 2531011)
+ val = (((seed >> 16) & 0x7fff) & 256) !== 0 ? val0 : -val0
+ At[i * astep + k] = val
+ }
+ for (iter = 0; iter < 2; iter++) {
+ for (j = 0; j < i; j++) {
+ sd = 0
+ for (k = 0; k < m; k++) {
+ sd += At[i * astep + k] * At[j * astep + k]
+ }
+ asum = 0.0
+ for (k = 0; k < m; k++) {
+ t = (At[i * astep + k] - sd * At[j * astep + k])
+ At[i * astep + k] = t
+ asum += Math.abs(t)
+ }
+ asum = asum ? 1.0 / asum : 0
+ for (k = 0; k < m; k++) {
+ At[i * astep + k] *= asum
+ }
+ }
+ }
+ sd = 0
+ for (k = 0; k < m; k++) {
+ t = At[i * astep + k]
+ sd += t * t
+ }
+ sd = Math.sqrt(sd)
+ }
+
+ s = (1.0 / sd)
+ for (k = 0; k < m; k++) {
+ At[i * astep + k] *= s
+ }
+ }
+}
+
+export function svd(A: Matrix, W: Matrix, U: Matrix, V: Matrix) {
+ let at = 0
+ let i = 0
+ const _m = A.rows
+ const _n = A.cols
+ let m = _m
+ let n = _n
+
+ if (m < n) {
+ at = 1
+ i = m
+ m = n
+ n = i
+ }
+
+ const amt = Matrix.create(m, m)
+ const wmt = Matrix.create(1, n)
+ const vmt = Matrix.create(n, n)
+
+ if (at === 0) {
+ Matrix.transpose(amt, A)
+ } else {
+ for (i = 0; i < _n * _m; i++) {
+ amt.data[i] = A.data[i]
+ }
+ for (; i < n * m; i++) {
+ amt.data[i] = 0
+ }
+ }
+
+ JacobiSVDImpl(amt.data, m, wmt.data, vmt.data, n, m, n, m)
+
+ if (W) {
+ for (i = 0; i < n; i++) {
+ W.data[i] = wmt.data[i]
+ }
+ for (; i < _n; i++) {
+ W.data[i] = 0
+ }
+ }
+
+ if (at === 0) {
+ if (U) Matrix.transpose(U, amt)
+ if (V) Matrix.transpose(V, vmt)
+ } else {
+ if (U) Matrix.transpose(U, vmt)
+ if (V) Matrix.transpose(V, amt)
+ }
+}
\ No newline at end of file