adding nearest and distance to point methods to lookup3d

src/mol-math/geometry/lookup3d/common.ts

@@ -41,8 +41,10 @@ export namespace Result {
 export interface Lookup3D<T = number> {
     // The result is mutated with each call to find.
     find(x: number, y: number, z: number, radius: number, result?: Result<T>): Result<T>,
+    nearest(x: number, y: number, z: number): { index: number, squaredDistance: number } | undefined,
+    distanceTo(x: number, y: number, z: number): number,
     check(x: number, y: number, z: number, radius: number): boolean,
     readonly boundary: { readonly box: Box3D, readonly sphere: Sphere3D }
     /** transient result */
     readonly result: Result<T>
-}
+}

src/mol-math/geometry/lookup3d/grid.ts

@@ -40,6 +40,40 @@ class GridLookup3DImpl<T extends number = number> implements GridLookup3D<T> {
         return ret;
     }
 
+    nearest(x: number, y: number, z: number): { index: number, squaredDistance: number } | undefined {
+        if (!OrderedSet.size(this.ctx.grid.data.indices)) return undefined;
+        this.ctx.x = x;
+        this.ctx.y = y;
+        this.ctx.z = z;
+        const radiusIncrement = 0.1; // how to choose a good increment?
+        const startingRadius = this.distanceTo(x, y, z);
+        this.ctx.radius = startingRadius < 0 ? radiusIncrement : startingRadius + radiusIncrement;
+        this.ctx.isCheck = false;
+        const result = this.result;
+        query(this.ctx, result);
+        while (!result.count) { // necessary to check if grid is empty to avoid infinite loop
+            this.ctx.radius += radiusIncrement;
+            query(this.ctx, result);
+        }
+        const { indices, squaredDistances } = result;
+        let index = indices[0], nearestDist = squaredDistances[0];
+
+        for (let i = 1, l = indices.length; i < l; i++) {
+            const dist = squaredDistances[i];
+            if (dist < nearestDist) {
+                nearestDist = dist;
+                index = indices[i];
+            }
+        }
+        return { index: index, squaredDistance: nearestDist };
+    }
+
+    distanceTo(x: number, y: number, z: number): number {
+        // distance between a sphere and a point
+        const { center, radius } = this.boundary.sphere;
+        return Vec3.distance(Vec3.create(x, y, z), center) - radius; // if negative, point is inside sphere
+    }
+
     check(x: number, y: number, z: number, radius: number): boolean {
         this.ctx.x = x;
         this.ctx.y = y;
@@ -277,4 +311,4 @@ function query<T extends number = number>(ctx: QueryContext, result: Result<T>):
         }
     }
     return result.count > 0;
-}
+}

src/mol-model/structure/structure/util/lookup3d.ts

@@ -86,6 +86,55 @@ export class StructureLookup3D {
         return ctx.result;
     }
 
+    nearest(x: number, y: number, z: number): { index: number, unit: Unit, squaredDistance: number } | undefined {
+        return this._nearest(x, y, z);
+    }
+
+    _nearest(x: number, y: number, z: number): { index: number, unit: Unit, squaredDistance: number } | undefined {
+        const ctx = this.findContext;
+        const { units, elementCount } = this.structure;
+
+        if (!elementCount) return undefined;
+
+        const radiusIncrement = 0.1; // how to choose a good increment?
+        const startingRadius = this.distanceTo(x, y, z);
+        let radius = startingRadius < 0 ? radiusIncrement : startingRadius + radiusIncrement;
+
+        while (true) {
+            const closeUnits = this.unitLookup.find(x, y, z, radius, ctx.closeUnitsResult);
+            radius += radiusIncrement;
+            console.log(radius)
+            if (closeUnits.count) {
+                const { indices, squaredDistances, count } = closeUnits;
+                let nearestIndex: number, nearestUnit: Unit, nearestDist = Number.MAX_SAFE_INTEGER;
+                for (let t = 1, _t = closeUnits.count; t < _t; t++) {
+                    const unit = units[closeUnits.indices[t]];
+                    Vec3.set(this.pivot, x, y, z);
+                    if (!unit.conformation.operator.isIdentity) {
+                        Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse);
+                    }
+                    const unitLookup = unit.lookup3d;
+                    const nearestResult = unitLookup.nearest(this.pivot[0], this.pivot[1], this.pivot[2]);
+                    if (nearestResult) {
+                        const { index: unitNearestIndex, squaredDistance: unitNearestDist } = nearestResult;
+                        if (unitNearestDist < nearestDist) {
+                            nearestDist = unitNearestDist;
+                            nearestIndex = unitNearestIndex;
+                            nearestUnit = unit;
+                        }
+                    }
+                }
+                return { index: nearestIndex!, unit: nearestUnit!, squaredDistance: nearestDist };
+            }
+        }
+    }
+
+    distanceTo(x: number, y: number, z: number): number {
+        // distance between a sphere and a point
+        const { center, radius } = this.boundary.sphere;
+        return Vec3.distance(Vec3.create(x, y, z), center) - radius; // if negative, point is inside sphere
+    }
+
     findIntoBuilder(x: number, y: number, z: number, radius: number, builder: StructureUniqueSubsetBuilder) {
         const { units } = this.structure;
         const closeUnits = this.unitLookup.find(x, y, z, radius);
@@ -217,4 +266,4 @@ export class StructureLookup3D {
         const position = { x: xs, y: ys, z: zs, radius, indices: OrderedSet.ofBounds(0, unitCount) };
         this.unitLookup = GridLookup3D(position, boundary);
     }
-}
+}