handle principal axes of points in a plane

CHANGELOG.md

@@ -7,6 +7,7 @@ Note that since we don't clearly distinguish between a public and private interf
 ## [Unreleased]
 
 - Fix: only update camera state if manualReset is off (#494)
+- Improve handling principal axes of points in a plane
 
 ## [v3.12.1] - 2022-07-20
 

src/mol-math/linear-algebra/matrix/principal-axes.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2022 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -9,6 +9,7 @@ import { Vec3 } from '../3d/vec3';
 import { svd } from './svd';
 import { NumberArray } from '../../../mol-util/type-helpers';
 import { Axes3D } from '../../geometry';
+import { EPSILON } from '../3d/common';
 
 export { PrincipalAxes };
 
@@ -58,10 +59,15 @@ namespace PrincipalAxes {
         return Axes3D.create(origin, dirA, dirB, dirC);
     }
 
+    export function calculateNormalizedAxes(momentsAxes: Axes3D): Axes3D {
+        const a = Axes3D.clone(momentsAxes);
+        if (Vec3.magnitude(a.dirC) < EPSILON) {
+            Vec3.cross(a.dirC, a.dirA, a.dirB);
+        }
+        return Axes3D.normalize(a, a);
+    }
+
     const tmpBoxVec = Vec3();
-    const tmpBoxVecA = Vec3();
-    const tmpBoxVecB = Vec3();
-    const tmpBoxVecC = Vec3();
     /**
      * Get the scale/length for each dimension for a box around the axes
      * to enclose the given positions
@@ -82,13 +88,11 @@ namespace PrincipalAxes {
         const t = Vec3();
 
         const center = momentsAxes.origin;
-        const normVecA = Vec3.normalize(tmpBoxVecA, momentsAxes.dirA);
-        const normVecB = Vec3.normalize(tmpBoxVecB, momentsAxes.dirB);
-        const normVecC = Vec3.normalize(tmpBoxVecC, momentsAxes.dirC);
+        const a = calculateNormalizedAxes(momentsAxes);
 
         for (let i = 0, il = positions.length; i < il; i += 3) {
-            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecA, center);
-            const dp1 = Vec3.dot(normVecA, Vec3.normalize(t, Vec3.sub(t, p, center)));
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), a.dirA, center);
+            const dp1 = Vec3.dot(a.dirA, Vec3.normalize(t, Vec3.sub(t, p, center)));
             const dt1 = Vec3.distance(p, center);
             if (dp1 > 0) {
                 if (dt1 > d1a) d1a = dt1;
@@ -96,8 +100,8 @@ namespace PrincipalAxes {
                 if (dt1 > d1b) d1b = dt1;
             }
 
-            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecB, center);
-            const dp2 = Vec3.dot(normVecB, Vec3.normalize(t, Vec3.sub(t, p, center)));
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), a.dirB, center);
+            const dp2 = Vec3.dot(a.dirB, Vec3.normalize(t, Vec3.sub(t, p, center)));
             const dt2 = Vec3.distance(p, center);
             if (dp2 > 0) {
                 if (dt2 > d2a) d2a = dt2;
@@ -105,8 +109,8 @@ namespace PrincipalAxes {
                 if (dt2 > d2b) d2b = dt2;
             }
 
-            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecC, center);
-            const dp3 = Vec3.dot(normVecC, Vec3.normalize(t, Vec3.sub(t, p, center)));
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), a.dirC, center);
+            const dp3 = Vec3.dot(a.dirC, Vec3.normalize(t, Vec3.sub(t, p, center)));
             const dt3 = Vec3.distance(p, center);
             if (dp3 > 0) {
                 if (dt3 > d3a) d3a = dt3;
@@ -115,16 +119,16 @@ namespace PrincipalAxes {
             }
         }
 
-        const dirA = Vec3.setMagnitude(Vec3(), normVecA, (d1a + d1b) / 2);
-        const dirB = Vec3.setMagnitude(Vec3(), normVecB, (d2a + d2b) / 2);
-        const dirC = Vec3.setMagnitude(Vec3(), normVecC, (d3a + d3b) / 2);
+        const dirA = Vec3.setMagnitude(Vec3(), a.dirA, (d1a + d1b) / 2);
+        const dirB = Vec3.setMagnitude(Vec3(), a.dirB, (d2a + d2b) / 2);
+        const dirC = Vec3.setMagnitude(Vec3(), a.dirC, (d3a + d3b) / 2);
 
         const origin = Vec3();
         const addCornerHelper = function (d1: number, d2: number, d3: number) {
             Vec3.copy(tmpBoxVec, center);
-            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecA, d1);
-            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecB, d2);
-            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecC, d3);
+            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, a.dirA, d1);
+            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, a.dirB, d2);
+            Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, a.dirC, d3);
             Vec3.add(origin, origin, tmpBoxVec);
         };
         addCornerHelper(d1a, d2a, d3a);

src/mol-model/structure/structure/carbohydrates/compute.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2022 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  * @author David Sehnal <david.sehnal@gmail.com>
@@ -141,7 +141,7 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
         Vec3.normalize(elements[iA].geometry.direction, elements[iA].geometry.direction);
     }
 
-    const tmpV = Vec3.zero();
+    const tmpV = Vec3();
     function fixTerminalLinkDirection(iA: number, indexB: number, unitB: Unit.Atomic) {
         const pos = unitB.conformation.position, geo = elements[iA].geometry;
         Vec3.sub(geo.direction, pos(unitB.elements[indexB], tmpV), geo.center);
@@ -189,9 +189,10 @@ export function computeCarbohydrates(structure: Structure): Carbohydrates {
                     const anomericCarbon = getAnomericCarbon(unit, ringAtoms);
 
                     const ma = PrincipalAxes.calculateMomentsAxes(getPositions(unit, ringAtoms));
-                    const center = Vec3.copy(Vec3.zero(), ma.origin);
-                    const normal = Vec3.copy(Vec3.zero(), ma.dirC);
-                    const direction = getDirection(Vec3.zero(), unit, anomericCarbon, center);
+                    const a = PrincipalAxes.calculateNormalizedAxes(ma);
+                    const center = Vec3.copy(Vec3(), a.origin);
+                    const normal = Vec3.copy(Vec3(), a.dirC);
+                    const direction = getDirection(Vec3(), unit, anomericCarbon, center);
                     Vec3.orthogonalize(direction, normal, direction);
 
                     const ringAltId = UnitRing.getAltId(unit, ringAtoms);