From b2809ad631c7635383f64a62d76f51f193e27c64 Mon Sep 17 00:00:00 2001
From: Alexander Rose <alexander.rose@weirdbyte.de>
Date: Sat, 24 Aug 2019 16:24:54 -0700
Subject: [PATCH] handle curves/fibers in mesoscale models
---
src/apps/viewer/extensions/cellpack/curve.ts | 193 +++++++++++++++++++
src/apps/viewer/extensions/cellpack/data.ts | 7 +-
src/apps/viewer/extensions/cellpack/model.ts | 53 +++--
3 files changed, 241 insertions(+), 12 deletions(-)
create mode 100644 src/apps/viewer/extensions/cellpack/curve.ts
diff --git a/src/apps/viewer/extensions/cellpack/curve.ts b/src/apps/viewer/extensions/cellpack/curve.ts
new file mode 100644
index 000000000..aeda0b1ea
--- /dev/null
+++ b/src/apps/viewer/extensions/cellpack/curve.ts
@@ -0,0 +1,193 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Ludovic Autin <autin@scripps.edu>
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Vec3, Quat, Mat4 } from '../../../../mol-math/linear-algebra';
+import { NumberArray } from '../../../../mol-util/type-helpers';
+
+interface Frame {
+ t: Vec3,
+ r: Vec3,
+ s: Vec3,
+}
+
+function CubicInterpolate(y0: Vec3, y1: Vec3, y2: Vec3, y3: Vec3, mu: number): Vec3 {
+ const out = Vec3.zero()
+ const mu2 = mu * mu;
+ const a0 = Vec3()
+ const a1 = Vec3()
+ const a2 = Vec3()
+ const a3 = Vec3()
+ Vec3.sub(a0, y3, y2)
+ Vec3.sub(a0, a0, y0)
+ Vec3.add(a0, a0, y1)
+
+ Vec3.sub(a1, y0, y1)
+ Vec3.sub(a1, a1, a0)
+
+ Vec3.sub(a2, y2, y0)
+
+ Vec3.copy(a3, y1)
+
+ out[0] = a0[0] * mu * mu2 + a1[0] * mu2 + a2[0] * mu + a3[0]
+ out[1] = a0[1] * mu * mu2 + a1[1] * mu2 + a2[1] * mu + a3[1]
+ out[2] = a0[2] * mu * mu2 + a1[2] * mu2 + a2[2] * mu + a3[2]
+
+ return out
+}
+
+function ResampleControlPoints(points: NumberArray, segmentLength: number) {
+ const nP = points.length / 3
+ // insert a point at the end and at the begining
+ // controlPoints.Insert(0, controlPoints[0] + (controlPoints[0] - controlPoints[1]) / 2.0f);
+ // controlPoints.Add(controlPoints[nP - 1] + (controlPoints[nP - 1] - controlPoints[nP - 2]) / 2.0f);
+
+ let resampledControlPoints: Vec3[] = []
+ // resampledControlPoints.Add(controlPoints[0]);
+ // resampledControlPoints.Add(controlPoints[1]);
+
+ let idx = 1
+ let currentPosition = Vec3.create(points[idx * 3], points[idx * 3 + 1], points[idx * 3 + 2])
+
+ let lerpValue = 0.0
+
+ // Normalize the distance between control points
+ while (true) {
+ if (idx + 2 >= nP) break
+ const cp0 = Vec3.create(points[(idx-1)*3], points[(idx-1)*3+1], points[(idx-1)*3+2]) // controlPoints[currentPointId - 1];
+ const cp1 = Vec3.create(points[idx*3], points[idx*3+1], points[idx*3+2]) // controlPoints[currentPointId];
+ const cp2 = Vec3.create(points[(idx+1)*3], points[(idx+1)*3+1], points[(idx+1)*3+2]) // controlPoints[currentPointId + 1];
+ const cp3 = Vec3.create(points[(idx+2)*3], points[(idx+2)*3+1], points[(idx+2)*3+2]); // controlPoints[currentPointId + 2];
+ let found = false
+ for (; lerpValue <= 1; lerpValue += 0.01) {
+ // lerp?slerp
+ // let candidate:Vec3 = Vec3.lerp(Vec3.zero(), cp0, cp1, lerpValue);
+ // const candidate:Vec3 = Vec3.bezier(Vec3.zero(), cp0, cp1, cp2, cp3, lerpValue);
+ const candidate = CubicInterpolate(cp0, cp1, cp2, cp3, lerpValue)
+ const d = Vec3.distance(currentPosition, candidate);
+ if (d > segmentLength) {
+ resampledControlPoints.push(candidate)
+ Vec3.copy(currentPosition, candidate)
+ found = true
+ break
+ }
+ }
+ if (!found) {
+ lerpValue = 0
+ idx += 1
+ }
+ }
+ return resampledControlPoints
+}
+
+// easier to align to theses normals
+function GetSmoothNormals(points: Vec3[]) {
+ const nP: number = points.length;
+ const smoothNormals: Vec3[] = []
+ if (points.length < 3) {
+ for (let i = 0; i < points.length; ++i)
+ smoothNormals.push(Vec3.normalize(Vec3(), points[i]))
+ return smoothNormals;
+ }
+ let p0 = Vec3.copy(Vec3(), points[0]) // undefined?
+ let p1 = Vec3.copy(Vec3(), points[1])
+ let p2 = Vec3.copy(Vec3(), points[2])
+ const p21 = Vec3.sub(Vec3(), p2, p1)
+ const p01 = Vec3.sub(Vec3(), p0, p1)
+ const p0121 = Vec3.cross(Vec3(), p01, p21)
+ let last = Vec3.normalize(Vec3(), p0121)
+ smoothNormals.push(last)
+ for (let i = 1; i < points.length - 1; ++i) {
+ p0 = points[i - 1]
+ p1 = points[i]
+ p2 = points[i + 1]
+ const t = Vec3.normalize(Vec3(), Vec3.sub(Vec3(), p2 , p0))
+ const b = Vec3.normalize(Vec3(), Vec3.cross(Vec3(), t, last))
+ let n = Vec3.normalize(Vec3(), Vec3.cross(Vec3(), t, b))
+ n = Vec3.scale(n, n, -1.0)
+ last = Vec3.copy(last, n)
+ smoothNormals.push(n)
+ }
+ last = Vec3.normalize(Vec3(), Vec3.cross(Vec3(), Vec3.sub(Vec3(), points[nP - 3], points[nP-2]), Vec3.sub(Vec3(), points[nP-2] , points[nP-1])))
+ smoothNormals.push(last)
+ return smoothNormals;
+}
+
+function getFrame(reference: Vec3, tangent: Vec3) {
+ const t: Vec3 = Vec3.normalize(Vec3(), tangent);
+ // make reference vector orthogonal to tangent
+ const proj_r_to_t: Vec3 = Vec3.scale(
+ Vec3(), tangent, Vec3.dot(reference, tangent) / Vec3.dot(tangent, tangent)
+ )
+ const r: Vec3 = Vec3.normalize(Vec3(), Vec3.sub(Vec3(), reference , proj_r_to_t))
+ // make bitangent vector orthogonal to the others
+ const s: Vec3 = Vec3.normalize(Vec3(), Vec3.cross(Vec3(), t, r))
+ return { t, r, s }
+}
+
+// easier to align to theses normals
+// https://github.com/bzamecnik/gpg/blob/master/rotation-minimizing-frame/rmf.py
+function GetMiniFrame(points: Vec3[], normals: Vec3[]) {
+ const frames: Frame[] = [];
+ const t0: Vec3 = Vec3.normalize(Vec3(), Vec3.sub(Vec3(), points[1], points[0]))
+ frames.push(getFrame(normals[0], t0))
+
+ for (let i = 0; i< points.length-2; ++i) {
+ const t2 = Vec3.normalize(Vec3(), Vec3.sub(Vec3(), points[i+2], points[i+1]))
+ const v1: Vec3 = Vec3.sub(Vec3(), points[i + 1], points[i]) // this is tangeant
+ const c1 = Vec3.dot(v1, v1)
+ // compute r_i^L = R_1 * r_i
+ const v1r = Vec3.scale(Vec3(), v1, (2.0/c1)*Vec3.dot(v1, frames[i].r))
+ const ref_L_i: Vec3 = Vec3.sub(Vec3(), frames[i].r, v1r)
+ // compute t_i^L = R_1 * t_i
+ const v1t = Vec3.scale(Vec3(), v1, (2.0/c1) * Vec3.dot(v1, frames[i].t))
+ const tan_L_i: Vec3 = Vec3.sub(Vec3(), frames[i].t, v1t)
+ // # compute reflection vector of R_2
+ const v2: Vec3 = Vec3.sub(Vec3(), t2 , tan_L_i)
+ const c2 = Vec3.dot(v2, v2)
+ // compute r_(i+1) = R_2 * r_i^L
+ const v2l = Vec3.scale(Vec3(), v1, (2.0/c2) * Vec3.dot(v2, ref_L_i))
+ const ref_next = Vec3.sub(Vec3(), ref_L_i, v2l) // ref_L_i - (2 / c2) * v2.dot(ref_L_i) * v2
+ frames.push(getFrame(ref_next, t2)) // frames.append(Frame(ref_next, tangents[i+1]))
+ }
+ return frames;
+}
+
+export function getMatFromResamplePoints(points: NumberArray) {
+ let segmentLength = 3.4
+ let new_points: Vec3[] = ResampleControlPoints(points, 3.4)
+ const npoints = new_points.length
+ let new_normal: Vec3[] = GetSmoothNormals(new_points)
+ let frames: Frame[] = GetMiniFrame(new_points, new_normal)
+ const limit = npoints
+ let transforms: Mat4[] = []
+ let pti: Vec3 = Vec3.copy(Vec3(), new_points[0]);
+ // console.log(new_points.length)
+ // console.log(points.length/3)
+ // console.log(limit)
+ // console.log(segmentLength)
+ for (let i = 0; i<npoints-2; ++i) {
+ const pti1: Vec3 = new_points[i+1] // Vec3.create(points[(i+1)*3],points[(i+1)*3+1],points[(i+1)*3+2]);
+ const d = Vec3.distance(pti, pti1)
+ if (d >= segmentLength) {
+ // use twist or random?
+ const quat: Quat = Quat.rotationTo(Quat.zero(), Vec3.create(0, 0, 1), frames[i].t) // Quat.rotationTo(Quat.zero(), Vec3.create(0,0,1),new_normal[i]);//Quat.rotationTo(Quat.zero(), Vec3.create(0,0,1),direction);new_normal
+ const rq: Quat = Quat.setAxisAngle(Quat.zero(), frames[i].t, Math.random()*3.60 ) // Quat.setAxisAngle(Quat.zero(),direction, Math.random()*3.60 );//Quat.identity();//
+ let m: Mat4 = Mat4.fromQuat(Mat4.zero(), Quat.multiply(Quat.zero(), rq, quat)) // Mat4.fromQuat(Mat4.zero(),Quat.multiply(Quat.zero(),quat1,quat2));//Mat4.fromQuat(Mat4.zero(),quat);//Mat4.identity();//Mat4.fromQuat(Mat4.zero(),Quat.multiply(Quat.zero(),rq,quat));
+ // let pos:Vec3 = Vec3.add(Vec3.zero(),pti1,pti)
+ // pos = Vec3.scale(pos,pos,1.0/2.0);
+ // Vec3.makeRotation(Mat4.zero(),Vec3.create(0,0,1),frames[i].t);//
+ Mat4.setTranslation(m, pti1)
+ // let m2:Mat4 = GetTubePropertiesMatrix(pti,pti1);
+ // let q:Quat = Quat.rotationTo(Quat.zero(), Vec3.create(0,1,0),Vec3.create(0,0,1))
+ // m2=Mat4.mul(Mat4.identity(),Mat4.fromQuat(Mat4.zero(),q),m2);
+ transforms.push(m)
+ pti = Vec3.copy(pti, pti1)
+ }
+ if (transforms.length >= limit) break
+ }
+ return transforms
+}
diff --git a/src/apps/viewer/extensions/cellpack/data.ts b/src/apps/viewer/extensions/cellpack/data.ts
index 006b38976..183dfab4f 100644
--- a/src/apps/viewer/extensions/cellpack/data.ts
+++ b/src/apps/viewer/extensions/cellpack/data.ts
@@ -27,7 +27,8 @@ export interface Cell {
export interface Recipe {
setupfile: string
- paths: [string, string][] // [name: string, path: string][]
+ /** First entry is name, secound is path: [name: string, path: string][] */
+ paths: [string, string][]
version: string
name: string
}
@@ -47,7 +48,11 @@ export interface Ingredient {
name: string
positions?: [Vec3[]] // why wrapped in an extra array?
radii?: [number[]] // why wrapped in an extra array?
+
+ /** Number of `curveX` properties in the object where `X` is a 0-indexed number */
nbCurve?: number
+ /** Curve properties are Vec3[] but that is not expressable in TypeScript */
+ [curveX: string]: unknown
}
export interface IngredientSource {
diff --git a/src/apps/viewer/extensions/cellpack/model.ts b/src/apps/viewer/extensions/cellpack/model.ts
index 18e80099c..d7e144071 100644
--- a/src/apps/viewer/extensions/cellpack/model.ts
+++ b/src/apps/viewer/extensions/cellpack/model.ts
@@ -10,7 +10,7 @@ import { PluginStateObject as PSO } from '../../../../mol-plugin/state/objects';
import { ParamDefinition as PD } from '../../../../mol-util/param-definition';
import { Ingredient, CellPacking } from './data';
import { getFromPdb, getFromCellPackDB } from './util';
-import { Model, Structure, StructureSymmetry, StructureSelection, Queries, QueryContext } from '../../../../mol-model/structure';
+import { Model, Structure, StructureSymmetry, StructureSelection, QueryContext } from '../../../../mol-model/structure';
import { trajectoryFromMmCIF } from '../../../../mol-model-formats/structure/mmcif';
import { trajectoryFromPDB } from '../../../../mol-model-formats/structure/pdb';
import { Mat4, Vec3, Quat } from '../../../../mol-math/linear-algebra';
@@ -24,6 +24,8 @@ import { Hcl } from '../../../../mol-util/color/spaces/hcl';
import { ParseCellPack, StructureFromCellpack } from './state';
import { formatMolScript } from '../../../../mol-script/language/expression-formatter';
import { MolScriptBuilder as MS } from '../../../../mol-script/language/builder';
+import { getMatFromResamplePoints } from './curve';
+import { compile } from '../../../../mol-script/runtime/query/compiler';
function getCellPackModelUrl(fileName: string, baseUrl: string) {
return `${baseUrl}/cellPACK_database_1.1.0/results/${fileName}`
@@ -50,8 +52,13 @@ async function getStructure(model: Model, props: { assembly?: string }) {
structure = await StructureSymmetry.buildAssembly(structure, assembly).run()
}
- const query = Queries.internal.atomicSequence()
- const result = query(new QueryContext(structure))
+ const query = MS.struct.modifier.union([
+ MS.struct.generator.atomGroups({
+ 'entity-test': MS.core.rel.eq([MS.ammp('entityType'), 'polymer'])
+ })
+ ])
+ const compiled = compile<StructureSelection>(query)
+ const result = compiled(new QueryContext(structure))
structure = StructureSelection.unionStructure(result)
return structure
@@ -66,10 +73,34 @@ function getTransform(trans: Vec3, rot: Quat) {
return m
}
-function getTransforms(results: Ingredient['results']) {
+function getResultTransforms(results: Ingredient['results']) {
return results.map((r: Ingredient['results'][0]) => getTransform(r[0], r[1]))
}
+function getCurveTransforms(ingredient: Ingredient) {
+ const n = ingredient.nbCurve || 0
+ const instances: Mat4[] = []
+
+ for (let i = 0; i < n; ++i) {
+ const cname = `curve${i}`
+ if (!(cname in ingredient)) {
+ // console.warn(`Expected '${cname}' in ingredient`)
+ continue
+ }
+ const _points = ingredient[cname] as Vec3[]
+ if (_points.length <= 2) {
+ // TODO handle curve with 2 or less points
+ continue
+ }
+ const points = new Float32Array(_points.length * 3)
+ for (let i = 0, il = _points.length; i < il; ++i) Vec3.toArray(_points[i], points, i * 3)
+ const newInstances = getMatFromResamplePoints(points)
+ instances.push(...newInstances)
+ }
+
+ return instances
+}
+
function getAssembly(transforms: Mat4[], structure: Structure) {
const builder = Structure.Builder()
const { units } = structure;
@@ -92,16 +123,16 @@ async function getIngredientStructure(ingredient: Ingredient, baseUrl: string) {
if (name === 'HIV1_CAhex_0_1_0') return
if (name === 'HIV1_CAhexCyclophilA_0_1_0') return
if (name === 'iLDL') return
- if (source.pdb === 'None') return
+ if (name === 'peptides') return
+ if (name === 'lypoglycane') return
- // TODO handle fibers
- if (nbCurve) return
+ if (source.pdb === 'None') return
const model = await getModel(source.pdb || name, baseUrl)
if (!model) return
const structure = await getStructure(model, { assembly: source.biomt ? '1' : undefined })
- const transforms = getTransforms(results)
+ const transforms = nbCurve ? getCurveTransforms(ingredient) : getResultTransforms(results)
const assembly = getAssembly(transforms, structure)
return assembly
}
@@ -139,12 +170,12 @@ export const LoadCellPackModel = StateAction.build({
id: PD.Select('influenza_model1.json', [
['blood_hiv_immature_inside.json', 'blood_hiv_immature_inside'],
['BloodHIV1.0_mixed_fixed_nc1.cpr', 'BloodHIV1.0_mixed_fixed_nc1'],
- ['BloodPlasma1.2.apr.json', 'BloodPlasma1.2'],
- ['BloodSerumfillResult.apr', 'BloodSerumfillResult'],
+ // ['BloodPlasma1.2.apr.json', 'BloodPlasma1.2'],
+ // ['BloodSerumfillResult.apr', 'BloodSerumfillResult'],
['HIV-1_0.1.6-8_mixed_radii_pdb.cpr', 'HIV-1_0.1.6-8_mixed_radii_pdb'],
['influenza_model1.json', 'influenza_model1'],
['Mycoplasma1.5_mixed_pdb_fixed.cpr', 'Mycoplasma1.5_mixed_pdb_fixed'],
- ['NM_Analysis_FigureC1.4.cpr.json', 'NM_Analysis_FigureC1.4']
+ // ['NM_Analysis_FigureC1.4.cpr.json', 'NM_Analysis_FigureC1.4']
]),
baseUrl: PD.Text('https://cdn.jsdelivr.net/gh/mesoscope/cellPACK_data@master/'),
preset: PD.Group({
--
GitLab