diff --git a/src/mol-math/linear-algebra/3d/minimize-rmsd.ts b/src/mol-math/linear-algebra/3d/minimize-rmsd.ts index 84444b6672419ee985876795e108dc5a320cec25..1296e10f3aa316031ac5a048b26c702b9269ebe8 100644 --- a/src/mol-math/linear-algebra/3d/minimize-rmsd.ts +++ b/src/mol-math/linear-algebra/3d/minimize-rmsd.ts @@ -19,6 +19,11 @@ namespace MinimizeRmsd { } export interface Positions { x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number> } + export namespace Positions { + export function empty(n: number) { + return { x: new Float64Array(n), y: new Float64Array(n), z: new Float64Array(n) }; + } + } export interface Input { a: Positions, @@ -54,10 +59,10 @@ class RmsdTransformState { this.b = data.b; if (data.centerA) this.centerA = data.centerA; - else this.centerA = CentroidHelper.compute(data.a, Vec3.zero()); + else this.centerA = data.centerA = CentroidHelper.compute(data.a, Vec3.zero()); if (data.centerB) this.centerB = data.centerB; - else this.centerB = CentroidHelper.compute(data.b, Vec3.zero()); + else this.centerB = data.centerB = CentroidHelper.compute(data.b, Vec3.zero()); this.result = into; } diff --git a/src/mol-model/structure/structure/element.ts b/src/mol-model/structure/structure/element.ts index 3e0ec4ec82004ed7b7d30d7d4a56716a61a35acd..da1e5a307f76c3a70415a4d2c4822d81cac747ae 100644 --- a/src/mol-model/structure/structure/element.ts +++ b/src/mol-model/structure/structure/element.ts @@ -118,6 +118,12 @@ namespace StructureElement { } export namespace Loci { + export function size(loci: Loci) { + let s = 0; + for (const u of loci.elements) s += OrderedSet.size(u.indices); + return s; + } + export function all(structure: Structure): Loci { return Loci(structure, structure.units.map(unit => ({ unit, diff --git a/src/mol-model/structure/structure/util/superposition.ts b/src/mol-model/structure/structure/util/superposition.ts new file mode 100644 index 0000000000000000000000000000000000000000..8c41ba996510e03eadbbf815c56715787fd41718 --- /dev/null +++ b/src/mol-model/structure/structure/util/superposition.ts @@ -0,0 +1,55 @@ +/** + * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author David Sehnal <david.sehnal@gmail.com> + */ + +import { MinimizeRmsd } from 'mol-math/linear-algebra'; +import StructureElement from '../element'; +import { OrderedSet } from 'mol-data/int'; + +export function superposeStructures(xs: StructureElement.Loci[]): MinimizeRmsd.Result[] { + const ret: MinimizeRmsd.Result[] = []; + if (xs.length <= 0) return ret; + + const n = getMinSize(xs); + const input: MinimizeRmsd.Input = { a: getPositionTable(xs[0], n), b: getPositionTable(xs[1], n) }; + ret[0] = MinimizeRmsd.compute(input); + for (let i = 2; i < xs.length; i++) { + input.b = getPositionTable(xs[i], n); + input.centerB = void 0; + ret.push(MinimizeRmsd.compute(input)); + } + + return ret; +} + +function getPositionTable(xs: StructureElement.Loci, n: number): MinimizeRmsd.Positions { + const ret = MinimizeRmsd.Positions.empty(n); + let o = 0; + for (const u of xs.elements) { + const { unit, indices } = u; + const { elements } = unit; + const { x, y, z } = unit.conformation; + for (let i = 0, _i = OrderedSet.size(indices); i < _i; i++) { + const e = elements[OrderedSet.getAt(indices, i)]; + ret.x[o] = x(e); + ret.y[o] = y(e); + ret.z[o] = z(e); + o++; + if (o >= n) break; + } + if (o >= n) break; + } + return ret; +} + +function getMinSize(xs: StructureElement.Loci[]) { + if (xs.length === 0) return 0; + let s = StructureElement.Loci.size(xs[0]); + for (let i = 1; i < xs.length; i++) { + const t = StructureElement.Loci.size(xs[i]); + if (t < s) s = t; + } + return s; +} \ No newline at end of file