diff --git a/src/mol-math/geometry/symmetry-operator.ts b/src/mol-math/geometry/symmetry-operator.ts
new file mode 100644
index 0000000000000000000000000000000000000000..1e3a57de38f5a782220061e838ca5ee2dafcdc28
--- /dev/null
+++ b/src/mol-math/geometry/symmetry-operator.ts
@@ -0,0 +1,141 @@
+/**
+ * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Vec3, Mat4 } from '../linear-algebra/3d'
+
+interface SymmetryOperator {
+ readonly name: string,
+ readonly hkl: Vec3,
+
+ readonly matrix: Mat4,
+ // cache the inverse of the transform
+ readonly inverse: Mat4,
+ // optimize the identity case
+ readonly isIdentity: boolean
+}
+
+namespace SymmetryOperator {
+ export const DefaultName = '1_555'
+ export const Default: SymmetryOperator = create(DefaultName, Mat4.identity());
+
+ export function create(name: string, matrix: Mat4, hkl?: number[]): SymmetryOperator {
+ const _hkl = hkl ? Vec3.create(hkl[0], hkl[1], hkl[2]) : Vec3.zero();
+ if (Mat4.isIdentity(matrix)) return { name, matrix, inverse: Mat4.identity(), isIdentity: true, hkl: _hkl };
+ if (!Mat4.isRotationAndTranslation(matrix)) throw new Error('Symmetry operator must be a composition of rotation and translation.');
+ return { name, matrix, inverse: Mat4.invert(Mat4.zero(), matrix), isIdentity: false, hkl: _hkl };
+ }
+
+ export interface CoordinateMapper { (index: number, slot: Vec3): Vec3 }
+ export interface ArrayMapping {
+ readonly invariantPosition: CoordinateMapper,
+ readonly position: CoordinateMapper,
+ x(index: number): number,
+ y(index: number): number,
+ z(index: number): number
+ }
+
+ export interface Coordinates { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> }
+
+ export function createMapping(operator: SymmetryOperator, coords: Coordinates) {
+ const invariantPosition = SymmetryOperator.createCoordinateMapper(SymmetryOperator.Default, coords);
+ const position = operator.isIdentity ? invariantPosition : SymmetryOperator.createCoordinateMapper(operator, coords);
+ const { x, y, z } = createProjections(operator, coords);
+ return { invariantPosition, position, x, y, z };
+ }
+
+ export function createCoordinateMapper(t: SymmetryOperator, coords: Coordinates): CoordinateMapper {
+ if (t.isIdentity) return identityPosition(coords);
+ return generalPosition(t, coords);
+ }
+}
+
+export default SymmetryOperator
+
+interface Projections { x(index: number): number, y(index: number): number, z(index: number): number }
+
+function createProjections(t: SymmetryOperator, coords: SymmetryOperator.Coordinates): Projections {
+ if (t.isIdentity) return { x: projectCoord(coords.x), y: projectCoord(coords.y), z: projectCoord(coords.z) };
+ return { x: projectX(t, coords), y: projectY(t, coords), z: projectZ(t, coords) };
+}
+
+function projectCoord(xs: ArrayLike<number>) {
+ return (i: number) => xs[i];
+}
+
+function isW1(m: Mat4) {
+ return m[3] === 0 && m[7] === 0 && m[11] === 0 && m[15] === 1;
+}
+
+function projectX({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+ const xx = m[0], yy = m[4], zz = m[8], tx = m[12];
+
+ if (isW1(m)) {
+ // this should always be the case.
+ return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + tx;
+ }
+
+ return (i: number) => {
+ const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+ return (xx * x + yy * y + zz * z + tx) / w;
+ }
+}
+
+function projectY({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+ const xx = m[1], yy = m[5], zz = m[9], ty = m[13];
+
+ if (isW1(m)) {
+ // this should always be the case.
+ return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + ty;
+ }
+
+ return (i: number) => {
+ const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+ return (xx * x + yy * y + zz * z + ty) / w;
+ }
+}
+
+function projectZ({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs}: SymmetryOperator.Coordinates) {
+ const xx = m[2], yy = m[6], zz = m[10], tz = m[14];
+
+ if (isW1(m)) {
+ // this should always be the case.
+ return (i: number) => xx * xs[i] + yy * ys[i] + zz * zs[i] + tz;
+ }
+
+ return (i: number) => {
+ const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
+ return (xx * x + yy * y + zz * z + tz) / w;
+ }
+}
+
+function identityPosition({ x, y, z }: SymmetryOperator.Coordinates): SymmetryOperator.CoordinateMapper {
+ return (i, s) => {
+ s[0] = x[i];
+ s[1] = y[i];
+ s[2] = z[i];
+ return s;
+ }
+}
+
+function generalPosition({ matrix: m }: SymmetryOperator, { x: xs, y: ys, z: zs }: SymmetryOperator.Coordinates) {
+ if (isW1(m)) {
+ // this should always be the case.
+ return (i: number, r: Vec3): Vec3 => {
+ const x = xs[i], y = ys[i], z = zs[i];
+ r[0] = m[0] * x + m[4] * y + m[8] * z + m[12];
+ r[1] = m[1] * x + m[5] * y + m[9] * z + m[13];
+ r[2] = m[2] * x + m[6] * y + m[10] * z + m[14];
+ return r;
+ }
+ }
+ return (i: number, r: Vec3): Vec3 => {
+ r[0] = xs[i];
+ r[1] = ys[i];
+ r[2] = zs[i];
+ Vec3.transformMat4(r, r, m);
+ return r;
+ }
+}
\ No newline at end of file
diff --git a/src/mol-math/geometry/transform.ts b/src/mol-math/geometry/transform.ts
deleted file mode 100644
index cf9aaa6d937b378c12edb14b9ae98376a23f8e99..0000000000000000000000000000000000000000
--- a/src/mol-math/geometry/transform.ts
+++ /dev/null
@@ -1,27 +0,0 @@
-/**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import { Mat4 } from '../linear-algebra/3d'
-
-interface Transform extends Readonly<{
- transform: Mat4,
- // cache the inverse of the transform
- inverse: Mat4,
- // optimize the identity case
- isIdentity: boolean
-}> { }
-
-namespace Transform {
- export function isIdentity(m: Mat4) {
- throw 'nyi'
- }
-
- export function isRotationAndTranslation(m: Mat4) {
- throw 'nyi'
- }
-}
-
-export default Transform
\ No newline at end of file
diff --git a/src/mol-math/linear-algebra/3d.ts b/src/mol-math/linear-algebra/3d.ts
index fccdb52b784f181d6d856963f875940a96447646..f989a61c08071c119239843a33b756dea1cbd267 100644
--- a/src/mol-math/linear-algebra/3d.ts
+++ b/src/mol-math/linear-algebra/3d.ts
@@ -152,7 +152,8 @@ export namespace Mat4 {
let det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
if (!det) {
- return Number.NaN;
+ console.warn('non-invertible matrix.', a);
+ return out;
}
det = 1.0 / det;
@@ -430,6 +431,18 @@ export namespace Mat4 {
// Calculate the determinant
return b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
}
+
+ export function isRotationAndTranslation(a: Mat4) {
+ const a00 = a[0], a01 = a[1], a02 = a[2], a03 = a[3],
+ a10 = a[4], a11 = a[5], a12 = a[6], a13 = a[7],
+ a20 = a[8], a21 = a[9], a22 = a[10], a23 = a[11],
+ /* a30 = a[12], a31 = a[13], a32 = a[14],*/ a33 = a[15];
+
+ if (a33 !== 1 || a03 !== 0 || a13 !== 0 || a23 !== 0) return false;
+ const det3x3 = a00 * (a11 * a22 - a21 * a21) - a01 * (a10 * a22 - a12 * a20) + a02 * (a10 * a21 - a11 * a20);
+ if (det3x3 < 1 - EPSILON.Value || det3x3 > 1 + EPSILON.Value) return false;
+ return true;
+ }
}
export namespace Vec3 {
diff --git a/src/mol-model/structure/model/formats/mmcif.ts b/src/mol-model/structure/model/formats/mmcif.ts
index dee0fddc23278dbb9ec3b999c3c8ff2670eb9037..aca855e52a55d203f85d7038679c1925ed68385e 100644
--- a/src/mol-model/structure/model/formats/mmcif.ts
+++ b/src/mol-model/structure/model/formats/mmcif.ts
@@ -70,9 +70,9 @@ function getConformation({ data }: mmCIF_Format, bounds: Interval): Conformation
atomId: Column.window(atom_site.id, start, end),
occupancy: Column.window(atom_site.occupancy, start, end),
B_iso_or_equiv: Column.window(atom_site.B_iso_or_equiv, start, end),
- __x: atom_site.Cartn_x.toArray({ array: Float32Array, start, end }),
- __y: atom_site.Cartn_y.toArray({ array: Float32Array, start, end }),
- __z: atom_site.Cartn_z.toArray({ array: Float32Array, start, end }),
+ x: atom_site.Cartn_x.toArray({ array: Float32Array, start, end }),
+ y: atom_site.Cartn_y.toArray({ array: Float32Array, start, end }),
+ z: atom_site.Cartn_z.toArray({ array: Float32Array, start, end }),
}
}
diff --git a/src/mol-model/structure/model/properties/conformation.ts b/src/mol-model/structure/model/properties/conformation.ts
index 66a975d62bf73b3b17c7873544fea784a959a04f..ce94b8597d2d1920b469f86d3942898e403fa5de 100644
--- a/src/mol-model/structure/model/properties/conformation.ts
+++ b/src/mol-model/structure/model/properties/conformation.ts
@@ -20,9 +20,9 @@ interface Conformation {
// Coordinates. Generally, not to be accessed directly because the coordinate might be
// transformed by an operator. Use Unit.getPosition instead.
- __x: ArrayLike<number>,
- __y: ArrayLike<number>,
- __z: ArrayLike<number>
+ x: ArrayLike<number>,
+ y: ArrayLike<number>,
+ z: ArrayLike<number>
}
export default Conformation
\ No newline at end of file
diff --git a/src/mol-model/structure/structure.ts b/src/mol-model/structure/structure.ts
index 1ceba8e5a8e0486cd28a7186f41ad85df1ae9987..4fe835cfaffe28eeb60b503fc14b6e713f2074b5 100644
--- a/src/mol-model/structure/structure.ts
+++ b/src/mol-model/structure/structure.ts
@@ -7,7 +7,6 @@
import Atom from './structure/atom'
import AtomSet from './structure/atom/set'
import Structure from './structure/structure'
-import Operator from './structure/operator'
import Unit from './structure/unit'
-export { Atom, AtomSet, Structure, Operator, Unit }
\ No newline at end of file
+export { Atom, AtomSet, Structure, Unit }
\ No newline at end of file
diff --git a/src/mol-model/structure/structure/operator.ts b/src/mol-model/structure/structure/operator.ts
deleted file mode 100644
index 10ba47cb8c94e46f587f5e39cb7958e980d2bcbb..0000000000000000000000000000000000000000
--- a/src/mol-model/structure/structure/operator.ts
+++ /dev/null
@@ -1,23 +0,0 @@
-/**
- * Copyright (c) 2017 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import { Mat4 } from 'mol-math/linear-algebra'
-
-interface Operator extends Readonly<{
- name: string,
- hkl: number[], // defaults to [0, 0, 0] for non symmetry entries
- matrix: Mat4,
- // cache the inverse of the transform
- inverse: Mat4,
- // optimize the identity case
- isIdentity: boolean
-}> { }
-
-namespace Operator {
- export const Identity: Operator = { name: '1_555', hkl: [0, 0, 0], matrix: Mat4.identity(), inverse: Mat4.identity(), isIdentity: true };
-}
-
-export default Operator
\ No newline at end of file
diff --git a/src/mol-model/structure/structure/structure.ts b/src/mol-model/structure/structure/structure.ts
index 0d09fa4db4b252e52620feaac8050b1c1fb60aea..0d262041252f50b998b7465ce7625921028560c0 100644
--- a/src/mol-model/structure/structure/structure.ts
+++ b/src/mol-model/structure/structure/structure.ts
@@ -6,9 +6,9 @@
import { OrderedSet, Iterator } from 'mol-data/int'
import UniqueArray from 'mol-data/util/unique-array'
+import SymmetryOperator from 'mol-math/geometry/symmetry-operator'
import { Model, Format } from '../model'
import Unit from './unit'
-import Operator from './operator'
import AtomSet from './atom/set'
import Atom from './atom'
@@ -35,7 +35,7 @@ namespace Structure {
const builder = Builder();
for (let c = 0; c < chains.count; c++) {
- const unit = Unit.create(model, Operator.Identity);
+ const unit = Unit.create(model, SymmetryOperator.Default);
builder.addUnit(unit);
builder.addAtoms(unit.id, OrderedSet.ofBounds(chains.segments[c], chains.segments[c + 1]));
}
diff --git a/src/mol-model/structure/structure/unit.ts b/src/mol-model/structure/structure/unit.ts
index c819836d99d6fe8e471ded60f6509d6c6dcf04b5..be83164a38f62a4d8dc99d4c190aa27d04c58343 100644
--- a/src/mol-model/structure/structure/unit.ts
+++ b/src/mol-model/structure/structure/unit.ts
@@ -4,45 +4,34 @@
* @author David Sehnal <david.sehnal@gmail.com>
*/
+import SymmetryOperator from 'mol-math/geometry/symmetry-operator'
import { Model } from '../model'
-import Operator from './operator'
-import { Vec3, Mat4 } from 'mol-math/linear-algebra'
-interface Unit extends Readonly<{
+interface Unit extends SymmetryOperator.ArrayMapping {
// Structure-level unique identifier of the unit.
- id: number,
+ readonly id: number,
// Provides access to the underlying data.
- model: Model,
+ readonly model: Model,
// Determines the operation applied to this unit.
// The transform and and inverse are baked into the "getPosition" function
- operator: Operator,
+ readonly operator: SymmetryOperator,
// Reference some commonly accessed things for faster access.
- residueIndex: ArrayLike<number>,
- chainIndex: ArrayLike<number>,
- hierarchy: Model['hierarchy'],
- conformation: Model['conformation']
-}> {
- // returns the untransformed position. Used for spatial queries.
- getInvariantPosition(atom: number, slot: Vec3): Vec3,
+ readonly residueIndex: ArrayLike<number>,
+ readonly chainIndex: ArrayLike<number>,
+ readonly hierarchy: Model['hierarchy'],
+ readonly conformation: Model['conformation']
- // gets the transformed position of the specified atom.
- getPosition(atom: number, slot: Vec3): Vec3,
-
- // optimized x/y/z coordinate projections for your query needs.
- x(atom: number): number,
- y(atom: number): number,
- z(atom: number): number
+ // TODO: add velocity?
}
namespace Unit {
- export function create(model: Model, operator: Operator): Unit {
+ export function create(model: Model, operator: SymmetryOperator): Unit {
const h = model.hierarchy;
- const { __x: x, __y: y, __z: z } = model.conformation;
- const getInvariantPosition = makeGetInvariantPosition(x, y, z);
- const getPosition = operator.isIdentity ? getInvariantPosition : makeGetPosition(operator.matrix, x, y, z);
+ const { invariantPosition, position, x, y, z } = SymmetryOperator.createMapping(operator, model.conformation);
+
return {
id: nextUnitId(),
model,
@@ -51,11 +40,9 @@ namespace Unit {
chainIndex: h.chainSegments.segmentMap,
hierarchy: model.hierarchy,
conformation: model.conformation,
- getInvariantPosition,
- getPosition,
- x: operator.isIdentity ? makeGetCoord(x) : makeX(operator.matrix, x, y, z),
- y: operator.isIdentity ? makeGetCoord(y) : makeY(operator.matrix, x, y, z),
- z: operator.isIdentity ? makeGetCoord(z) : makeZ(operator.matrix, x, y, z)
+ invariantPosition,
+ position,
+ x, y, z
};
}
}
@@ -67,73 +54,4 @@ function nextUnitId() {
const ret = _id;
_id = (_id + 1) % 0x3fffffff;
return ret;
-}
-
-function makeGetInvariantPosition(x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>) {
- return (a: number, r: Vec3): Vec3 => {
- r[0] = x[a];
- r[1] = y[a];
- r[2] = z[a];
- return r;
- }
-}
-
-function makeGetCoord(xs: ArrayLike<number>) {
- return (i: number) => xs[i];
-}
-
-function isWConst(m: Mat4) {
- return m[3] === 0 && m[7] === 0 && m[11] === 0;
-}
-
-function makeX(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
- const xx = m[0], yy = m[4], zz = m[8], ww = m[12];
-
- if (isWConst(m)) {
- const w = m[15] || 1.0;
- return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
- }
-
- return (i: number) => {
- const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
- return (xx * x + yy * y + zz * z + ww) / w;
- }
-}
-
-function makeY(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
- const xx = m[1], yy = m[5], zz = m[9], ww = m[13];
-
- if (isWConst(m)) {
- const w = m[15] || 1.0;
- return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
- }
-
- return (i: number) => {
- const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
- return (xx * x + yy * y + zz * z + ww) / w;
- }
-}
-
-function makeZ(m: Mat4, xs: ArrayLike<number>, ys: ArrayLike<number>, zs: ArrayLike<number>) {
- const xx = m[2], yy = m[6], zz = m[10], ww = m[14];
-
- if (isWConst(m)) {
- const w = m[15] || 1.0;
- return (i: number) => (xx * xs[i] + yy * ys[i] + zz * zs[i] + ww) / w;
- }
-
- return (i: number) => {
- const x = xs[i], y = ys[i], z = zs[i], w = (m[3] * x + m[7] * y + m[11] * z + m[15]) || 1.0;
- return (xx * x + yy * y + zz * z + ww) / w;
- }
-}
-
-function makeGetPosition(m: Mat4, x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>) {
- return (a: number, r: Vec3): Vec3 => {
- r[0] = x[a];
- r[1] = y[a];
- r[2] = z[a];
- Vec3.transformMat4(r, r, m);
- return r;
- }
}
\ No newline at end of file
diff --git a/src/perf-tests/structure.ts b/src/perf-tests/structure.ts
index 3ae723fc755aa4eae7250d397a0bed8c21113a44..156d441f8393095f4d3236f42e561262d88b7719 100644
--- a/src/perf-tests/structure.ts
+++ b/src/perf-tests/structure.ts
@@ -242,8 +242,8 @@ export namespace PropertyAccess {
}
export async function run() {
- //const { structures, models } = await readCIF('./examples/1cbs_full.bcif');
- const { structures, models, mmcif } = await readCIF('e:/test/quick/3j3q_full.bcif');
+ const { structures, models, mmcif } = await readCIF('./examples/1cbs_full.bcif');
+ //const { structures, models, mmcif } = await readCIF('e:/test/quick/3j3q_full.bcif');
//const { structures, models, mmcif } = await readCIF('e:/test/quick/1cbs_updated.cif');
console.log(mmcif.pdbx_struct_oper_list.matrix.toArray());
@@ -272,7 +272,7 @@ export namespace PropertyAccess {
console.log('atom.x', sumProperty(structures[0], Q.props.atom.x));
console.timeEnd('atom.x');
console.time('__x')
- console.log('__x', sumProperty(structures[0], l => l.unit.conformation.__x[l.atom]));
+ console.log('__x', sumProperty(structures[0], l => l.unit.conformation.x[l.atom]));
console.timeEnd('__x')
//const authSeqId = Atom.property(l => l.unit.hierarchy.residues.auth_seq_id.value(l.unit.residueIndex[l.atom]));