From 735e497a1d1a8d021635f6f20e204d1869c6904d Mon Sep 17 00:00:00 2001
From: Alexander Rose <alexander.rose@weirdbyte.de>
Date: Wed, 10 Apr 2019 00:59:58 -0700
Subject: [PATCH] wip, molsurf
---
src/mol-math/geometry/molecular-surface.ts | 41 ++++++++++++-------
.../visual/util/molecular-surface.ts | 9 ++--
2 files changed, 31 insertions(+), 19 deletions(-)
diff --git a/src/mol-math/geometry/molecular-surface.ts b/src/mol-math/geometry/molecular-surface.ts
index 254c34815..96303db45 100644
--- a/src/mol-math/geometry/molecular-surface.ts
+++ b/src/mol-math/geometry/molecular-surface.ts
@@ -49,6 +49,8 @@ function getAngleTables (probePositions: number): AnglesTables {
* Ignore indices a and b (these are the relevant atoms in projectPoints/Torii)
*
* Cache the last clipped atom (as very often the same one in subsequent calls)
+ *
+ * `a` and `b` must be resolved indices
*/
function obscured (state: MolSurfCalcState, x: number, y: number, z: number, a: number, b: number) {
if (state.lastClip !== -1) {
@@ -61,7 +63,7 @@ function obscured (state: MolSurfCalcState, x: number, y: number, z: number, a:
}
for (let j = 0, jl = state.neighbours.count; j < jl; ++j) {
- const ai = state.neighbours.indices[j]
+ const ai = OrderedSet.getAt(state.position.indices, state.neighbours.indices[j])
if (ai !== a && ai !== b && singleAtomObscures(state, ai, x, y, z)) {
state.lastClip = ai
return ai
@@ -71,12 +73,14 @@ function obscured (state: MolSurfCalcState, x: number, y: number, z: number, a:
return -1
}
+/**
+ * `ai` must be a resolved index
+ */
function singleAtomObscures (state: MolSurfCalcState, ai: number, x: number, y: number, z: number) {
- const j = OrderedSet.getAt(state.position.indices, ai)
- const r = state.position.radius[j]
- const dx = state.position.x[j] - x
- const dy = state.position.y[j] - y
- const dz = state.position.z[j] - z
+ const r = state.position.radius[ai]
+ const dx = state.position.x[ai] - x
+ const dy = state.position.y[ai] - y
+ const dz = state.position.z[ai] - z
const dSq = dx * dx + dy * dy + dz * dz
return dSq < (r * r)
}
@@ -92,7 +96,7 @@ function singleAtomObscures (state: MolSurfCalcState, ai: number, x: number, y:
* itself and delta
*/
async function projectPoints (ctx: RuntimeContext, state: MolSurfCalcState) {
- const { position, lookup3d, min, space, data, idData, scaleFactor } = state
+ const { position, lookup3d, min, space, data, idData, scaleFactor, updateChunk } = state
const { gridx, gridy, gridz } = state
const { indices, x, y, z, radius } = position
@@ -154,7 +158,7 @@ async function projectPoints (ctx: RuntimeContext, state: MolSurfCalcState) {
const spy = dy * ap + vy
const spz = dz * ap + vz
- if (obscured(state, spx, spy, spz, i, -1) === -1) {
+ if (obscured(state, spx, spy, spz, j, -1) === -1) {
const dd = rad - d
if (dd < data[idx]) {
data[idx] = dd
@@ -166,7 +170,7 @@ async function projectPoints (ctx: RuntimeContext, state: MolSurfCalcState) {
}
}
- if (i % 10000 === 0 && ctx.shouldUpdate) {
+ if (i % updateChunk === 0 && ctx.shouldUpdate) {
await ctx.update({ message: 'projecting points', current: i, max: n })
}
}
@@ -177,6 +181,9 @@ const atob = Vec3()
const mid = Vec3()
const n1 = Vec3()
const n2 = Vec3()
+/**
+ * `a` and `b` must be resolved indices
+ */
function projectTorus (state: MolSurfCalcState, a: number, b: number) {
const { position, min, space, data, idData } = state
const { cosTable, sinTable, probePositions, probeRadius, scaleFactor } = state
@@ -195,7 +202,7 @@ function projectTorus (state: MolSurfCalcState, a: number, b: number) {
const dSq = dx * dx + dy * dy + dz * dz
// This check now redundant as already done in AVHash.withinRadii
- if (dSq > ((rA + rB) * (rA + rB))) { return }
+ // if (dSq > ((rA + rB) * (rA + rB))) { return }
const d = Math.sqrt(dSq)
@@ -271,7 +278,7 @@ function projectTorus (state: MolSurfCalcState, a: number, b: number) {
// Is this grid point closer to a or b?
// Take dot product of atob and gridpoint->p (dx, dy, dz)
const dp = dx * atob[0] + dy * atob[1] + dz * atob[2]
- idData[idx] = dp < 0.0 ? b : a
+ idData[idx] = OrderedSet.indexOf(position.indices, dp < 0.0 ? b : a)
}
}
}
@@ -281,17 +288,17 @@ function projectTorus (state: MolSurfCalcState, a: number, b: number) {
}
async function projectTorii (ctx: RuntimeContext, state: MolSurfCalcState) {
- const { n, lookup3d, position } = state
+ const { n, lookup3d, position, updateChunk } = state
const { x, y, z, indices, radius } = position
for (let i = 0; i < n; ++i) {
const k = OrderedSet.getAt(indices, i)
state.neighbours = lookup3d.find(x[k], y[k], z[k], radius[k])
for (let j = 0, jl = state.neighbours.count; j < jl; ++j) {
- const l = state.neighbours.indices[j]
+ const l = OrderedSet.getAt(indices, state.neighbours.indices[j])
if (k < l) projectTorus(state, k, l)
}
- if (i % 10000 === 0 && ctx.shouldUpdate) {
+ if (i % updateChunk === 0 && ctx.shouldUpdate) {
await ctx.update({ message: 'projecting torii', current: i, max: n })
}
}
@@ -308,6 +315,7 @@ interface MolSurfCalcState {
lookup3d: Lookup3D
position: Required<PositionData>
min: Vec3
+ updateChunk: number
maxRadius: number
@@ -336,7 +344,9 @@ async function createState(ctx: RuntimeContext, position: Required<PositionData>
const { resolution, probeRadius, probePositions } = props
const scaleFactor = 1 / resolution
- const lookup3d = GridLookup3D(position)
+ const cellSize = Vec3.create(maxRadius, maxRadius, maxRadius)
+ Vec3.scale(cellSize, cellSize, 2)
+ const lookup3d = GridLookup3D(position, cellSize)
const box = lookup3d.boundary.box
const { indices } = position
const n = OrderedSet.size(indices)
@@ -370,6 +380,7 @@ async function createState(ctx: RuntimeContext, position: Required<PositionData>
lookup3d,
position,
min,
+ updateChunk,
maxRadius,
diff --git a/src/mol-repr/structure/visual/util/molecular-surface.ts b/src/mol-repr/structure/visual/util/molecular-surface.ts
index 35f7c3697..9a045e1d1 100644
--- a/src/mol-repr/structure/visual/util/molecular-surface.ts
+++ b/src/mol-repr/structure/visual/util/molecular-surface.ts
@@ -17,15 +17,16 @@ export function computeUnitMolecularSurface(unit: Unit, props: MolecularSurfaceC
return Task.create('Molecular Surface', async ctx => {
const { indices } = position
const n = OrderedSet.size(indices)
- const radii = new Float32Array(n)
+ const radii = new Float32Array(OrderedSet.max(indices))
let maxRadius = 0
for (let i = 0; i < n; ++i) {
- const r = radius(OrderedSet.getAt(indices, i))
+ const j = OrderedSet.getAt(indices, i)
+ const r = radius(j)
if (maxRadius < r) maxRadius = r
- radii[i] = r + props.probeRadius
+ radii[j] = r + props.probeRadius
- if (i % 10000 === 0 && ctx.shouldUpdate) {
+ if (i % 100000 === 0 && ctx.shouldUpdate) {
await ctx.update({ message: 'calculating max radius', current: i, max: n })
}
}
--
GitLab