From 0b038559eb69625b740f36b9adeb53da9cd5d0aa Mon Sep 17 00:00:00 2001 From: David Sehnal <david.sehnal@gmail.com> Date: Thu, 19 Apr 2018 18:09:30 +0200 Subject: [PATCH] structure boundary --- src/mol-math/geometry/lookup3d/common.ts | 2 +- src/mol-math/linear-algebra/3d/vec3.ts | 8 +++ .../structure/structure/structure.ts | 4 ++ .../structure/structure/util/boundary.ts | 70 +++++++++++++++++++ .../structure/structure/util/lookup3d.ts | 31 +++++--- src/perf-tests/lookup3d.ts | 3 + 6 files changed, 109 insertions(+), 9 deletions(-) create mode 100644 src/mol-model/structure/structure/util/boundary.ts diff --git a/src/mol-math/geometry/lookup3d/common.ts b/src/mol-math/geometry/lookup3d/common.ts index 2cf669f88..0aa58601d 100644 --- a/src/mol-math/geometry/lookup3d/common.ts +++ b/src/mol-math/geometry/lookup3d/common.ts @@ -32,5 +32,5 @@ export interface Lookup3D<T = number> { // The result is mutated with each call to find. find(x: number, y: number, z: number, radius: number): Result<T>, check(x: number, y: number, z: number, radius: number): boolean, - boundary: { box: Box3D, sphere: Sphere3D } + readonly boundary: { readonly box: Box3D, readonly sphere: Sphere3D } } \ No newline at end of file diff --git a/src/mol-math/linear-algebra/3d/vec3.ts b/src/mol-math/linear-algebra/3d/vec3.ts index bbf28f3ac..5a48c82e4 100644 --- a/src/mol-math/linear-algebra/3d/vec3.ts +++ b/src/mol-math/linear-algebra/3d/vec3.ts @@ -65,6 +65,14 @@ namespace Vec3 { return out; } + export function ofArray(array: ArrayLike<number>) { + const out = zero(); + out[0] = array[0]; + out[1] = array[1]; + out[2] = array[2]; + return out; + } + export function set(out: Vec3, x: number, y: number, z: number): Vec3 { out[0] = x; out[1] = y; diff --git a/src/mol-model/structure/structure/structure.ts b/src/mol-model/structure/structure/structure.ts index 60ddb8f65..5dbcd62ec 100644 --- a/src/mol-model/structure/structure/structure.ts +++ b/src/mol-model/structure/structure/structure.ts @@ -91,6 +91,10 @@ namespace Structure { return s.__lookup3d__; } + export function getBoundary(s: Structure) { + return getLookup3d(s).boundary; + } + // TODO: "lift" atom set operators? // TODO: "diff" } diff --git a/src/mol-model/structure/structure/util/boundary.ts b/src/mol-model/structure/structure/util/boundary.ts new file mode 100644 index 000000000..83655a247 --- /dev/null +++ b/src/mol-model/structure/structure/util/boundary.ts @@ -0,0 +1,70 @@ +/** + * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author David Sehnal <david.sehnal@gmail.com> + */ + +import Structure from '../structure' +import ElementSet from '../element/set' +import { ElementGroup } from '../../structure' +import { Box3D, Sphere3D } from 'mol-math/geometry'; +import { Vec3 } from 'mol-math/linear-algebra'; + +function computeStructureBoundary(s: Structure): { box: Box3D, sphere: Sphere3D } { + const min = [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE]; + const max = [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE]; + + const { units, elements } = s; + + let cx = 0, cy = 0, cz = 0; + let radiusSq = 0; + let size = 0; + + for (let i = 0, _i = ElementSet.unitCount(elements); i < _i; i++) { + const group = ElementSet.unitGetByIndex(elements, i); + const { x, y, z } = units[ElementSet.unitGetId(elements, i)]; + + size += ElementGroup.size(group); + for (let j = 0, _j = ElementGroup.size(group); j < _j; j++) { + const e = ElementGroup.getAt(group, j); + const xx = x(e), yy = y(e), zz = z(e); + + min[0] = Math.min(xx, min[0]); + min[1] = Math.min(yy, min[1]); + min[2] = Math.min(zz, min[2]); + max[0] = Math.max(xx, max[0]); + max[1] = Math.max(yy, max[1]); + max[2] = Math.max(zz, max[2]); + + cx += xx; + cy += yy; + cz += zz; + } + } + + if (size > 0) { + cx /= size; + cy /= size; + cz /= size; + } + + for (let i = 0, _i = ElementSet.unitCount(elements); i < _i; i++) { + const group = ElementSet.unitGetByIndex(elements, i); + const { x, y, z } = units[ElementSet.unitGetId(elements, i)]; + + size += ElementGroup.size(group); + for (let j = 0, _j = ElementGroup.size(group); j < _j; j++) { + const e = ElementGroup.getAt(group, j); + const dx = x(e) - cx, dy = y(e) - cy, dz = z(e) - cz; + const d = dx * dx + dy * dy + dz * dz; + if (d > radiusSq) radiusSq = d; + } + } + + return { + box: { min: Vec3.ofArray(min), max: Vec3.ofArray(max) }, + sphere: { center: Vec3.create(cx, cy, cz), radius: Math.sqrt(radiusSq) } + }; +} + +export { computeStructureBoundary } \ No newline at end of file diff --git a/src/mol-model/structure/structure/util/lookup3d.ts b/src/mol-model/structure/structure/util/lookup3d.ts index b34f5e947..7434852bc 100644 --- a/src/mol-model/structure/structure/util/lookup3d.ts +++ b/src/mol-model/structure/structure/util/lookup3d.ts @@ -10,6 +10,7 @@ import { Lookup3D, GridLookup3D, Result, Box3D, Sphere3D } from 'mol-math/geomet import { ElementGroup, ElementSet } from '../../structure'; import { Vec3 } from 'mol-math/linear-algebra'; import { OrderedSet } from 'mol-data/int'; +import { computeStructureBoundary } from './boundary'; interface StructureLookup3D extends Lookup3D<Element> {} @@ -25,8 +26,6 @@ namespace StructureLookup3D { const closeUnits = this.unitLookup.find(x, y, z, radius); if (closeUnits.count === 0) return this.result; - console.log({ closeUnits }); - for (let t = 0, _t = closeUnits.count; t < _t; t++) { const i = closeUnits.indices[t]; const unitId = ElementSet.unitGetId(elements, i); @@ -38,8 +37,6 @@ namespace StructureLookup3D { } const groupLookup = ElementGroup.getLookup3d(unit, group); const groupResult = groupLookup.find(this.pivot[0], this.pivot[1], this.pivot[2], radius); - //console.log(groupLookup); - //console.log({ groupCount: groupResult.count }); for (let j = 0, _j = groupResult.count; j < _j; j++) { Result.add(this.result, Element.create(unitId, groupResult.indices[j]), groupResult.squaredDistances[j]); } @@ -47,10 +44,29 @@ namespace StructureLookup3D { return this.result; } + check(x: number, y: number, z: number, radius: number): boolean { - throw new Error("Method not implemented."); + const { units, elements } = this.structure; + const closeUnits = this.unitLookup.find(x, y, z, radius); + if (closeUnits.count === 0) return false; + + for (let t = 0, _t = closeUnits.count; t < _t; t++) { + const i = closeUnits.indices[t]; + const unitId = ElementSet.unitGetId(elements, i); + const group = ElementSet.unitGetByIndex(elements, i); + const unit = units[unitId]; + Vec3.set(this.pivot, x, y, z); + if (!unit.operator.isIdentity) { + Vec3.transformMat4(this.pivot, this.pivot, unit.operator.inverse); + } + const groupLookup = ElementGroup.getLookup3d(unit, group); + if (groupLookup.check(this.pivot[0], this.pivot[1], this.pivot[2], radius)) return true; + } + + return false; } - boundary: { box: Box3D; sphere: Sphere3D; } = 0 as any; + + boundary: { box: Box3D; sphere: Sphere3D; }; constructor(private structure: Structure) { const { units, elements } = structure; @@ -75,9 +91,8 @@ namespace StructureLookup3D { radius[i] = s.radius; } - console.log({ xs, ys, zs, radius, unitCount }) - this.unitLookup = GridLookup3D({ x: xs, y: ys, z: zs, radius, indices: OrderedSet.ofBounds(0, unitCount) }); + this.boundary = computeStructureBoundary(structure); } } diff --git a/src/perf-tests/lookup3d.ts b/src/perf-tests/lookup3d.ts index 574eebaa6..cadc871cc 100644 --- a/src/perf-tests/lookup3d.ts +++ b/src/perf-tests/lookup3d.ts @@ -56,6 +56,9 @@ export async function test() { const sl = Structure.getLookup3d(structures[0]); const result1 = sl.find(-30.07, 8.178, -13.897, 10); console.log(result1.count);//, result1.indices); + + console.log(Structure.getBoundary(structures[0])); + console.log(lookup.boundary); } test(); \ No newline at end of file -- GitLab