handle curves/fibers in mesoscale models

src/apps/viewer/extensions/cellpack/curve.ts

+/**
+ * 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
+}

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 {

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({