diff --git a/.vscode/extensions.json b/.vscode/extensions.json index 966b67347d36bc37436d3ba78525db8a5a4275fc..93c043ab7c334b348ea377d93146bf0506bebd6b 100644 --- a/.vscode/extensions.json +++ b/.vscode/extensions.json @@ -6,7 +6,6 @@ "recommendations": [ "dbaeumer.vscode-eslint", "firsttris.vscode-jest-runner", - "msjsdiag.debugger-for-chrome", "slevesque.shader", "stpn.vscode-graphql", "wayou.vscode-todo-highlight" diff --git a/CHANGELOG.md b/CHANGELOG.md index db5514c0ce0cadc427001db6bab2180ce220e745..dd6600dfda724293d007181cfd210cd9244d8e06 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,11 @@ Note that since we don't clearly distinguish between a public and private interf ## [Unreleased] -- [Breaking] Renamed the ``model-index`` color theme to ``trajectory-index`` -- Added a new ``model-index`` color theme that uses the `Model.Index` property -- Added the new ``model-index`` color theme as an option for the carbon color in the ``element-symbol`` and ``ilustrative`` color themes +- [Breaking] Rename the ``model-index`` color theme to ``trajectory-index`` +- Add a new ``model-index`` color theme that uniquely colors each loaded model +- Add the new ``model-index`` color theme as an option for the carbon color in the ``element-symbol`` and ``ilustrative`` color themes +- Add nearest method to lookup3d +- Add mipmap-based blur for skybox backgrounds ## [v3.19.0] - 2022-10-01 diff --git a/src/extensions/backgrounds/index.ts b/src/extensions/backgrounds/index.ts index c87c8899f367a38bf377f6d266586faf0fe32ae0..c5bbdc99e1f013f3658ffa660cc0df4889db5e44 100644 --- a/src/extensions/backgrounds/index.ts +++ b/src/extensions/backgrounds/index.ts @@ -72,6 +72,7 @@ export const Backgrounds = PluginBehavior.create<{ }>({ lightness: 0, saturation: 0, opacity: 1, + blur: 0.3, } } }, 'Purple Nebula Skybox'], diff --git a/src/mol-canvas3d/passes/background.ts b/src/mol-canvas3d/passes/background.ts index d4bfb3e59aa3fb02cfdc4e958b61618dbc500235..14d194b6baa418f1e381cf47b182735590df2be2 100644 --- a/src/mol-canvas3d/passes/background.ts +++ b/src/mol-canvas3d/passes/background.ts @@ -49,6 +49,7 @@ const SkyboxParams = { pz: PD.File({ label: 'Positive Z / Front', accept: 'image/*' }), }, { isExpanded: true, label: 'Files' }), }), + blur: PD.Numeric(0, { min: 0.0, max: 1.0, step: 0.01 }, { description: 'Note, this only works in WebGL2 or when "EXT_shader_texture_lod" is available.' }), ...SharedParams, }; type SkyboxProps = PD.Values<typeof SkyboxParams> @@ -170,6 +171,7 @@ export class BackgroundPass { Mat4.invert(m, m); ValueCell.update(this.renderable.values.uViewDirectionProjectionInverse, m); + ValueCell.updateIfChanged(this.renderable.values.uBlur, props.blur); ValueCell.updateIfChanged(this.renderable.values.uOpacity, props.opacity); ValueCell.updateIfChanged(this.renderable.values.uSaturation, props.saturation); ValueCell.updateIfChanged(this.renderable.values.uLightness, props.lightness); @@ -367,7 +369,7 @@ function getSkyboxTexture(ctx: WebGLContext, assetManager: AssetManager, faces: const cubeAssets = getCubeAssets(assetManager, faces); const cubeFaces = getCubeFaces(assetManager, cubeAssets); const assets = [cubeAssets.nx, cubeAssets.ny, cubeAssets.nz, cubeAssets.px, cubeAssets.py, cubeAssets.pz]; - const texture = ctx.resources.cubeTexture(cubeFaces, false, onload); + const texture = ctx.resources.cubeTexture(cubeFaces, true, onload); return { texture, assets }; } @@ -424,12 +426,15 @@ const BackgroundSchema = { uGradientColorA: UniformSpec('v3'), uGradientColorB: UniformSpec('v3'), uGradientRatio: UniformSpec('f'), + uBlur: UniformSpec('f'), uOpacity: UniformSpec('f'), uSaturation: UniformSpec('f'), uLightness: UniformSpec('f'), dVariant: DefineSpec('string', ['skybox', 'image', 'verticalGradient', 'horizontalGradient', 'radialGradient']), }; -const SkyboxShaderCode = ShaderCode('background', background_vert, background_frag); +const SkyboxShaderCode = ShaderCode('background', background_vert, background_frag, { + shaderTextureLod: 'optional' +}); type BackgroundRenderable = ComputeRenderable<Values<typeof BackgroundSchema>> function getBackgroundRenderable(ctx: WebGLContext, width: number, height: number): BackgroundRenderable { @@ -448,6 +453,7 @@ function getBackgroundRenderable(ctx: WebGLContext, width: number, height: numbe uGradientColorA: ValueCell.create(Vec3()), uGradientColorB: ValueCell.create(Vec3()), uGradientRatio: ValueCell.create(0.5), + uBlur: ValueCell.create(0), uOpacity: ValueCell.create(1), uSaturation: ValueCell.create(0), uLightness: ValueCell.create(0), diff --git a/src/mol-canvas3d/passes/smaa.ts b/src/mol-canvas3d/passes/smaa.ts index 4b1de5e41389eddec14bfbadc39d90ca95a3bc15..79752e4da69232a69081ab6795baf71c7488ba42 100644 --- a/src/mol-canvas3d/passes/smaa.ts +++ b/src/mol-canvas3d/passes/smaa.ts @@ -1,5 +1,5 @@ /** - * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info. + * Copyright (c) 2020-2022 mol* contributors, licensed under MIT, See LICENSE file for more info. * * @author Alexander Rose <alexander.rose@weirdbyte.de> */ @@ -11,7 +11,7 @@ import { ShaderCode } from '../../mol-gl/shader-code'; import { WebGLContext } from '../../mol-gl/webgl/context'; import { createComputeRenderItem } from '../../mol-gl/webgl/render-item'; import { RenderTarget } from '../../mol-gl/webgl/render-target'; -import { createTexture, loadImageTexture, Texture } from '../../mol-gl/webgl/texture'; +import { loadImageTexture, Texture } from '../../mol-gl/webgl/texture'; import { Vec2, Vec4 } from '../../mol-math/linear-algebra'; import { ValueCell } from '../../mol-util'; import { ParamDefinition as PD } from '../../mol-util/param-definition'; @@ -192,8 +192,8 @@ function getWeightsRenderable(ctx: WebGLContext, edgesTexture: Texture): Weights const width = edgesTexture.getWidth(); const height = edgesTexture.getHeight(); - const areaTexture = createTexture(ctx.gl, ctx.extensions, 'image-uint8', 'rgb', 'ubyte', 'linear'); - const searchTexture = createTexture(ctx.gl, ctx.extensions, 'image-uint8', 'rgba', 'ubyte', 'nearest'); + const areaTexture = ctx.resources.texture('image-uint8', 'rgb', 'ubyte', 'linear'); + const searchTexture = ctx.resources.texture('image-uint8', 'rgba', 'ubyte', 'nearest'); const values: Values<typeof WeightsSchema> = { ...QuadValues, diff --git a/src/mol-geo/geometry/text/font-atlas.ts b/src/mol-geo/geometry/text/font-atlas.ts index 609a6985a50ab450799a602ed1ae91caaa95beb3..834b8ecbd49a1f278644496d67946c9554a0501e 100644 --- a/src/mol-geo/geometry/text/font-atlas.ts +++ b/src/mol-geo/geometry/text/font-atlas.ts @@ -1,5 +1,5 @@ /** - * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info. + * Copyright (c) 2019-2022 mol* contributors, licensed under MIT, See LICENSE file for more info. * * @author Alexander Rose <alexander.rose@weirdbyte.de> */ @@ -88,7 +88,7 @@ export class FontAtlas { this.scratchCanvas.width = this.maxWidth; this.scratchCanvas.height = this.lineHeight; - this.scratchContext = this.scratchCanvas.getContext('2d')!; + this.scratchContext = this.scratchCanvas.getContext('2d', { willReadFrequently: true })!; this.scratchContext.font = `${p.fontStyle} ${p.fontVariant} ${p.fontWeight} ${fontSize}px ${p.fontFamily}`; this.scratchContext.fillStyle = 'black'; this.scratchContext.textBaseline = 'middle'; diff --git a/src/mol-gl/shader/background.frag.ts b/src/mol-gl/shader/background.frag.ts index a764a9ad8237ec635bb53ad45c7bb7e08f297c4f..835314ff237890841d3d5c57ffccd73eed72911b 100644 --- a/src/mol-gl/shader/background.frag.ts +++ b/src/mol-gl/shader/background.frag.ts @@ -6,6 +6,7 @@ precision mediump sampler2D; #if defined(dVariant_skybox) uniform samplerCube tSkybox; uniform mat4 uViewDirectionProjectionInverse; + uniform float uBlur; uniform float uOpacity; uniform float uSaturation; uniform float uLightness; @@ -49,7 +50,11 @@ vec3 lightenColor(vec3 c, float amount) { void main() { #if defined(dVariant_skybox) vec4 t = uViewDirectionProjectionInverse * vPosition; - gl_FragColor = textureCube(tSkybox, normalize(t.xyz / t.w)); + #ifdef enabledShaderTextureLod + gl_FragColor = textureCubeLodEXT(tSkybox, normalize(t.xyz / t.w), uBlur * 8.0); + #else + gl_FragColor = textureCube(tSkybox, normalize(t.xyz / t.w)); + #endif gl_FragColor.a = uOpacity; gl_FragColor.rgb = lightenColor(saturateColor(gl_FragColor.rgb, uSaturation), uLightness); #elif defined(dVariant_image) diff --git a/src/mol-math/geometry/_spec/lookup3d.spec.ts b/src/mol-math/geometry/_spec/lookup3d.spec.ts index 1641a2c6e85a431a265c4d8708baf453976b872a..9b0fe2d35f617048544926f20f6e4e61416aa308 100644 --- a/src/mol-math/geometry/_spec/lookup3d.spec.ts +++ b/src/mol-math/geometry/_spec/lookup3d.spec.ts @@ -2,6 +2,7 @@ * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info. * * @author David Sehnal <david.sehnal@gmail.com> + * @author Gianluca Tomasello <giagitom@gmail.com> */ import { GridLookup3D } from '../../geometry'; @@ -24,9 +25,17 @@ describe('GridLookup3d', () => { expect(r.count).toBe(1); expect(r.indices[0]).toBe(0); + r = grid.nearest(0, 0, 0, 1); + expect(r.count).toBe(1); + expect(r.indices[0]).toBe(0); + r = grid.find(0, 0, 0, 1); expect(r.count).toBe(3); expect(sortArray(r.indices)).toEqual([0, 1, 2]); + + r = grid.nearest(0, 0, 0, 3); + expect(r.count).toBe(3); + expect(sortArray(r.indices)).toEqual([0, 1, 2]); }); it('radius', () => { @@ -38,9 +47,17 @@ describe('GridLookup3d', () => { expect(r.count).toBe(1); expect(r.indices[0]).toBe(0); + r = grid.nearest(0, 0, 0, 1); + expect(r.count).toBe(1); + expect(r.indices[0]).toBe(0); + r = grid.find(0, 0, 0, 0.5); expect(r.count).toBe(2); expect(sortArray(r.indices)).toEqual([0, 1]); + + r = grid.nearest(0, 0, 0, 3); + expect(r.count).toBe(3); + expect(sortArray(r.indices)).toEqual([0, 1, 2]); }); it('indexed', () => { @@ -51,8 +68,15 @@ describe('GridLookup3d', () => { let r = grid.find(0, 0, 0, 0); expect(r.count).toBe(0); + r = grid.nearest(0, 0, 0, 1); + expect(r.count).toBe(1); + r = grid.find(0, 0, 0, 0.5); expect(r.count).toBe(1); expect(sortArray(r.indices)).toEqual([0]); + + r = grid.nearest(0, 0, 0, 3); + expect(r.count).toBe(1); + expect(sortArray(r.indices)).toEqual([0]); }); -}); \ No newline at end of file +}); diff --git a/src/mol-math/geometry/lookup3d/common.ts b/src/mol-math/geometry/lookup3d/common.ts index 24b457d77d5207956cf997c1bd346e5dc8af6ade..2c6c76cc2635963617d02ac10d6f5fa47e9ca0f3 100644 --- a/src/mol-math/geometry/lookup3d/common.ts +++ b/src/mol-math/geometry/lookup3d/common.ts @@ -41,8 +41,9 @@ 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, k: number, stopIf?: Function, result?: Result<T>): Result<T>, check(x: number, y: number, z: number, radius: number): boolean, readonly boundary: { readonly box: Box3D, readonly sphere: Sphere3D } /** transient result */ readonly result: Result<T> -} \ No newline at end of file +} diff --git a/src/mol-math/geometry/lookup3d/grid.ts b/src/mol-math/geometry/lookup3d/grid.ts index 970cce5881913c9d5447ca613dba7b26619b5378..b4f96dd0f9aee2c8bc1c4aaec740ed3e5fbbde03 100644 --- a/src/mol-math/geometry/lookup3d/grid.ts +++ b/src/mol-math/geometry/lookup3d/grid.ts @@ -3,6 +3,7 @@ * * @author David Sehnal <david.sehnal@gmail.com> * @author Alexander Rose <alexander.rose@weirdbyte.de> + * @author Gianluca Tomasello <giagitom@gmail.com> */ import { Result, Lookup3D } from './common'; @@ -12,6 +13,7 @@ import { PositionData } from '../common'; import { Vec3 } from '../../linear-algebra'; import { OrderedSet } from '../../../mol-data/int'; import { Boundary } from '../boundary'; +import { FibonacciHeap } from '../../../mol-util/fibonacci-heap'; interface GridLookup3D<T = number> extends Lookup3D<T> { readonly buckets: { readonly offset: ArrayLike<number>, readonly count: ArrayLike<number>, readonly array: ArrayLike<number> } @@ -40,6 +42,17 @@ class GridLookup3DImpl<T extends number = number> implements GridLookup3D<T> { return ret; } + nearest(x: number, y: number, z: number, k: number = 1, stopIf?: Function, result?: Result<T>): Result<T> { + this.ctx.x = x; + this.ctx.y = y; + this.ctx.z = z; + this.ctx.k = k; + this.ctx.stopIf = stopIf; + const ret = result ?? this.result; + queryNearest(this.ctx, ret); + return ret; + } + check(x: number, y: number, z: number, radius: number): boolean { this.ctx.x = x; this.ctx.y = y; @@ -221,12 +234,14 @@ interface QueryContext { x: number, y: number, z: number, + k: number, + stopIf?: Function, radius: number, isCheck: boolean } function createContext(grid: Grid3D): QueryContext { - return { grid, x: 0.1, y: 0.1, z: 0.1, radius: 0.1, isCheck: false }; + return { grid, x: 0.1, y: 0.1, z: 0.1, k: 1, stopIf: undefined, radius: 0.1, isCheck: false }; } function query<T extends number = number>(ctx: QueryContext, result: Result<T>): boolean { @@ -277,4 +292,152 @@ function query<T extends number = number>(ctx: QueryContext, result: Result<T>): } } return result.count > 0; -} \ No newline at end of file +} + +const tmpDirVec = Vec3(); +const tmpVec = Vec3(); +const tmpSetG = new Set<number>(); +const tmpSetG2 = new Set<number>(); +const tmpArrG1 = [0.1]; +const tmpArrG2 = [0.1]; +const tmpArrG3 = [0.1]; +const tmpHeapG = new FibonacciHeap(); +function queryNearest<T extends number = number>(ctx: QueryContext, result: Result<T>): boolean { + const { min, expandedBox: box, boundingSphere: { center }, size: [sX, sY, sZ], bucketOffset, bucketCounts, bucketArray, grid, data: { x: px, y: py, z: pz, indices, radius }, delta, maxRadius } = ctx.grid; + const { x, y, z, k, stopIf } = ctx; + const indicesCount = OrderedSet.size(indices); + Result.reset(result); + if (indicesCount === 0 || k <= 0) return false; + let gX, gY, gZ, stop = false, gCount = 1, expandGrid = true, nextGCount = 0, arrG = tmpArrG1, nextArrG = tmpArrG2, maxRange = 0, expandRange = true, gridId: number, gridPointsFinished = false; + const expandedArrG = tmpArrG3, sqMaxRadius = maxRadius * maxRadius; + arrG.length = 0; + expandedArrG.length = 0; + tmpSetG.clear(); + tmpHeapG.clear(); + Vec3.set(tmpVec, x, y, z); + if (!Box3D.containsVec3(box, tmpVec)) { + // intersect ray pointing to box center + Box3D.nearestIntersectionWithRay(tmpVec, box, tmpVec, Vec3.normalize(tmpDirVec, Vec3.sub(tmpDirVec, center, tmpVec))); + gX = Math.max(0, Math.min(sX - 1, Math.floor((tmpVec[0] - min[0]) / delta[0]))); + gY = Math.max(0, Math.min(sY - 1, Math.floor((tmpVec[1] - min[1]) / delta[1]))); + gZ = Math.max(0, Math.min(sZ - 1, Math.floor((tmpVec[2] - min[2]) / delta[2]))); + } else { + gX = Math.floor((x - min[0]) / delta[0]); + gY = Math.floor((y - min[1]) / delta[1]); + gZ = Math.floor((z - min[2]) / delta[2]); + } + const dX = maxRadius !== 0 ? Math.max(1, Math.min(sX - 1, Math.ceil(maxRadius / delta[0]))) : 1; + const dY = maxRadius !== 0 ? Math.max(1, Math.min(sY - 1, Math.ceil(maxRadius / delta[1]))) : 1; + const dZ = maxRadius !== 0 ? Math.max(1, Math.min(sZ - 1, Math.ceil(maxRadius / delta[2]))) : 1; + arrG.push(gX, gY, gZ, (((gX * sY) + gY) * sZ) + gZ); + while (result.count < indicesCount) { + const arrGLen = gCount * 4; + for (let ig = 0; ig < arrGLen; ig += 4) { + gridId = arrG[ig + 3]; + if (!tmpSetG.has(gridId)) { + tmpSetG.add(gridId); + gridPointsFinished = tmpSetG.size >= grid.length; + const bucketIdx = grid[gridId]; + if (bucketIdx !== 0) { + const _maxRange = maxRange; + const ki = bucketIdx - 1; + const offset = bucketOffset[ki]; + const count = bucketCounts[ki]; + const end = offset + count; + for (let i = offset; i < end; i++) { + const bIdx = bucketArray[i]; + const idx = OrderedSet.getAt(indices, bIdx); + const dx = px[idx] - x; + const dy = py[idx] - y; + const dz = pz[idx] - z; + let distSq = dx * dx + dy * dy + dz * dz; + if (maxRadius !== 0) { + const r = radius![idx]; + distSq -= r * r; + } + if (expandRange && distSq > maxRange) { + maxRange = distSq; + } + tmpHeapG.insert(distSq, bIdx); + } + if (_maxRange < maxRange) expandRange = false; + } + } + } + // find next grid points + nextArrG.length = 0; + nextGCount = 0; + tmpSetG2.clear(); + for (let ig = 0; ig < arrGLen; ig += 4) { + gX = arrG[ig]; + gY = arrG[ig + 1]; + gZ = arrG[ig + 2]; + // fill grid points array with valid adiacent positions + for (let ix = -dX; ix <= dX; ix++) { + const xPos = gX + ix; + if (xPos < 0 || xPos >= sX) continue; + for (let iy = -dY; iy <= dY; iy++) { + const yPos = gY + iy; + if (yPos < 0 || yPos >= sY) continue; + for (let iz = -dZ; iz <= dZ; iz++) { + const zPos = gZ + iz; + if (zPos < 0 || zPos >= sZ) continue; + gridId = (((xPos * sY) + yPos) * sZ) + zPos; + if (tmpSetG2.has(gridId)) continue; // already scanned + tmpSetG2.add(gridId); + if (tmpSetG.has(gridId)) continue; // already visited + if (!expandGrid) { + const xP = min[0] + xPos * delta[0] - x; + const yP = min[1] + yPos * delta[1] - y; + const zP = min[2] + zPos * delta[2] - z; + const distSqG = (xP * xP) + (yP * yP) + (zP * zP) - sqMaxRadius; // is sqMaxRadius necessary? + if (distSqG > maxRange) { + expandedArrG.push(xPos, yPos, zPos, gridId); + continue; + } + } + nextArrG.push(xPos, yPos, zPos, gridId); + nextGCount++; + } + } + } + } + expandGrid = false; + if (nextGCount === 0) { + if (k === 1) { + const node = tmpHeapG.findMinimum(); + if (node) { + const { key: squaredDistance, value: index } = node!; + // const squaredDistance = node!.key, index = node!.value; + Result.add(result, index as number, squaredDistance as number); + return true; + } + } else { + while (!tmpHeapG.isEmpty() && (gridPointsFinished || tmpHeapG.findMinimum()!.key as number <= maxRange) && result.count < k) { + const node = tmpHeapG.extractMinimum(); + const squaredDistance = node!.key, index = node!.value; + Result.add(result, index as number, squaredDistance as number); + if (stopIf && !stop) { + stop = stopIf(index, squaredDistance); + } + } + } + if (result.count >= k || stop || result.count >= indicesCount) return result.count > 0; + expandGrid = true; + expandRange = true; + if (expandedArrG.length > 0) { + for (let i = 0, l = expandedArrG.length; i < l; i++) { + arrG.push(expandedArrG[i]); + } + expandedArrG.length = 0; + gCount = arrG.length; + } + } else { + const tmp = arrG; + arrG = nextArrG; + nextArrG = tmp; + gCount = nextGCount; + } + } + return result.count > 0; +} diff --git a/src/mol-math/geometry/primitives/box3d.ts b/src/mol-math/geometry/primitives/box3d.ts index fff4923615107fddbbe9be77febc27b070182889..0176be2aed2a120beb44c398d851539a3d1c1f43 100644 --- a/src/mol-math/geometry/primitives/box3d.ts +++ b/src/mol-math/geometry/primitives/box3d.ts @@ -138,6 +138,48 @@ namespace Box3D { a.max[2] < b.min[2] || a.min[2] > b.max[2] ); } + + // const tmpTransformV = Vec3(); + export function nearestIntersectionWithRay(out: Vec3, box: Box3D, origin: Vec3, dir: Vec3): Vec3 { + const [minX, minY, minZ] = box.min; + const [maxX, maxY, maxZ] = box.max; + const [x, y, z] = origin; + const invDirX = 1.0 / dir[0]; + const invDirY = 1.0 / dir[1]; + const invDirZ = 1.0 / dir[2]; + let tmin, tmax, tymin, tymax, tzmin, tzmax; + if (invDirX >= 0) { + tmin = (minX - x) * invDirX; + tmax = (maxX - x) * invDirX; + } else { + tmin = (maxX - x) * invDirX; + tmax = (minX - x) * invDirX; + } + if (invDirY >= 0) { + tymin = (minY - y) * invDirY; + tymax = (maxY - y) * invDirY; + } else { + tymin = (maxY - y) * invDirY; + tymax = (minY - y) * invDirY; + } + if (invDirZ >= 0) { + tzmin = (minZ - z) * invDirZ; + tzmax = (maxZ - z) * invDirZ; + } else { + tzmin = (maxZ - z) * invDirZ; + tzmax = (minZ - z) * invDirZ; + } + if (tymin > tmin) + tmin = tymin; + if (tymax < tmax) + tmax = tymax; + if (tzmin > tmin) + tmin = tzmin; + if (tzmax < tmax) + tmax = tzmax; + Vec3.scale(out, dir, tmin); + return Vec3.set(out, out[0] + x, out[1] + y, out[2] + z); + } } -export { Box3D }; \ No newline at end of file +export { Box3D }; diff --git a/src/mol-math/geometry/primitives/sphere3d.ts b/src/mol-math/geometry/primitives/sphere3d.ts index 6f09cce9503363b2c5da245f0198bf8a478c63f7..69953f360f467ce8a3ae6200b4b2517c135a7fe4 100644 --- a/src/mol-math/geometry/primitives/sphere3d.ts +++ b/src/mol-math/geometry/primitives/sphere3d.ts @@ -277,6 +277,12 @@ namespace Sphere3D { export function distance(a: Sphere3D, b: Sphere3D) { return Vec3.distance(a.center, b.center) - a.radius + b.radius; } + + /** Get the distance of v from sphere. If negative, v is inside sphere */ + export function distanceToVec(sphere: Sphere3D, v: Vec3): number { + const { center, radius } = sphere; + return Vec3.distance(v, center) - radius; + } } -export { Sphere3D }; \ No newline at end of file +export { Sphere3D }; diff --git a/src/mol-model/structure/structure/util/lookup3d.ts b/src/mol-model/structure/structure/util/lookup3d.ts index 49de19cbd20e0ffa6790767068e6f519aabc48fe..693d17aca13aaa4d607ab3826a3a6b834deaca40 100644 --- a/src/mol-model/structure/structure/util/lookup3d.ts +++ b/src/mol-model/structure/structure/util/lookup3d.ts @@ -3,6 +3,7 @@ * * @author David Sehnal <david.sehnal@gmail.com> * @author Alexander Rose <alexander.rose@weirdbyte.de> + * @author Gianluca Tomasello <giagitom@gmail.com> */ import { Structure } from '../structure'; @@ -13,6 +14,7 @@ import { StructureUniqueSubsetBuilder } from './unique-subset-builder'; import { StructureElement } from '../element'; import { Unit } from '../unit'; import { UnitIndex } from '../element/util'; +import { FibonacciHeap } from '../../../../mol-util/fibonacci-heap'; export interface StructureResult extends Result<StructureElement.UnitIndex> { units: Unit[] @@ -54,6 +56,7 @@ export function StructureLookup3DResultContext(): StructureLookup3DResultContext export class StructureLookup3D { private unitLookup: Lookup3D; private pivot = Vec3(); + private heap = new FibonacciHeap(); findUnitIndices(x: number, y: number, z: number, radius: number): Result<number> { return this.unitLookup.find(x, y, z, radius); @@ -86,6 +89,54 @@ export class StructureLookup3D { return ctx.result; } + nearest(x: number, y: number, z: number, k: number = 1, ctx?: StructureLookup3DResultContext): StructureResult { + return this._nearest(x, y, z, k, ctx ?? this.findContext); + } + + _nearest(x: number, y: number, z: number, k: number, ctx: StructureLookup3DResultContext): StructureResult { + const result = ctx.result, heap = this.heap; + Result.reset(result); + heap.clear(); + const { units } = this.structure; + let elementsCount = 0; + const closeUnits = this.unitLookup.nearest(x, y, z, units.length, (uid: number) => (elementsCount += units[uid].elements.length) >= k, ctx.closeUnitsResult); // sort units based on distance to the point + if (closeUnits.count === 0) return result; + let totalCount = 0, maxDistResult = -Number.MAX_VALUE; + for (let t = 0, _t = closeUnits.count; t < _t; t++) { + const unitSqDist = closeUnits.squaredDistances[t]; + if (totalCount >= k && maxDistResult < unitSqDist) break; + Vec3.set(this.pivot, x, y, z); + const unit = units[closeUnits.indices[t]]; + if (!unit.conformation.operator.isIdentity) { + Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse); + } + const unitLookup = unit.lookup3d; + const groupResult = unitLookup.nearest(this.pivot[0], this.pivot[1], this.pivot[2], k, void 0, ctx.unitGroupResult); + if (groupResult.count === 0) continue; + totalCount += groupResult.count; + maxDistResult = Math.max(maxDistResult, groupResult.squaredDistances[groupResult.count - 1]); + for (let j = 0, _j = groupResult.count; j < _j; j++) { + heap.insert(groupResult.squaredDistances[j], { index: groupResult.indices[j], unit: unit }); + } + } + if (k === 1) { + const node = heap.findMinimum(); + if (node) { + const { key: squaredDistance } = node; + const { unit, index } = node.value as { index: UnitIndex, unit: Unit }; + StructureResult.add(result, unit as Unit, index as UnitIndex, squaredDistance as number); + } + } else { + while (!heap.isEmpty() && result.count < k) { + const node = heap.extractMinimum(); + const { key: squaredDistance } = node!; + const { unit, index } = node!.value as { index: UnitIndex, unit: Unit }; + StructureResult.add(result, unit as Unit, index as UnitIndex, squaredDistance as number); + } + } + return result; + } + 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 +268,4 @@ export class StructureLookup3D { const position = { x: xs, y: ys, z: zs, radius, indices: OrderedSet.ofBounds(0, unitCount) }; this.unitLookup = GridLookup3D(position, boundary); } -} \ No newline at end of file +} diff --git a/src/mol-util/_spec/fibonacci-heap.spec.ts b/src/mol-util/_spec/fibonacci-heap.spec.ts new file mode 100644 index 0000000000000000000000000000000000000000..4016e8faf9655ff72f6d95f01174fa0798cd55be --- /dev/null +++ b/src/mol-util/_spec/fibonacci-heap.spec.ts @@ -0,0 +1,22 @@ +/** + * Copyright (c) 2022 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author Gianluca Tomasello <giagitom@gmail.com> + */ + +import { FibonacciHeap } from '../fibonacci-heap'; + +describe('fibonacci-heap', () => { + it('basic', () => { + const heap = new FibonacciHeap(); + heap.insert(1, 2); + heap.insert(4); + heap.insert(2); + heap.insert(3); + expect(heap.size()).toBe(4); + const node = heap.extractMinimum(); + expect(node!.key).toBe(1); + expect(node!.value).toBe(2); + expect(heap.size()).toBe(3); + }); +}); diff --git a/src/mol-util/fibonacci-heap.ts b/src/mol-util/fibonacci-heap.ts new file mode 100644 index 0000000000000000000000000000000000000000..5602bd4b837117db57bed611d2ae282fa2e8fb33 --- /dev/null +++ b/src/mol-util/fibonacci-heap.ts @@ -0,0 +1,407 @@ +/** + * Copyright (c) 2022 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author Gianluca Tomasello <giagitom@gmail.com> + * + * Adapted from https://github.com/gwtw/ts-fibonacci-heap, Copyright (c) 2014 Daniel Imms, MIT + */ + +interface INode<K, V> { + key: K; + value?: V; +} + +type CompareFunction<K, V> = (a: INode<K, V>, b: INode<K, V>) => number; + +class Node<K, V> implements INode<K, V> { + public key: K; + public value: V | undefined; + public prev: Node<K, V>; + public next: Node<K, V>; + public parent: Node<K, V> | null = null; + public child: Node<K, V> | null = null; + + public degree: number = 0; + public isMarked: boolean = false; + + constructor(key: K, value?: V) { + this.key = key; + this.value = value; + this.prev = this; + this.next = this; + } +} + +class NodeListIterator<K, V> { + private _index: number; + private _items: Node<K, V>[]; + private _len: number; + /** + * Creates an Iterator used to simplify the consolidate() method. It works by + * making a shallow copy of the nodes in the root list and iterating over the + * shallow copy instead of the source as the source will be modified. + * @param start A node from the root list. + */ + constructor(start?: Node<K, V>) { + this._index = -1; + this._items = []; + this._len = 0; + if (start) { + let current = start, l = 0; + do { + this._items[l++] = current; + current = current.next; + } while (start !== current); + this._len = l; + } + } + + /** + * @return Whether there is a next node in the iterator. + */ + public hasNext(): boolean { + return this._index < this._len - 1; + } + + /** + * @return The next node. + */ + public next(): Node<K, V> { + return this._items[++this._index]; + } + + /** + * @return Resets iterator to reuse it. + */ + public reset(start: Node<K, V>) { + this._index = -1; + this._len = 0; + let current = start, l = 0; + do { + this._items[l++] = current; + current = current.next; + } while (start !== current); + this._len = l; + } +} + +const tmpIt = new NodeListIterator<any, any>(); +/** + * A Fibonacci heap data structure with a key and optional value. +*/ +export class FibonacciHeap<K, V> { + private _minNode: Node<K, V> | null = null; + private _nodeCount: number = 0; + private _compare: CompareFunction<K, V>; + + constructor( + compare?: CompareFunction<K, V> + ) { + this._compare = compare ? compare : this._defaultCompare; + } + + /** + * Clears the heap's data, making it an empty heap. + */ + public clear(): void { + this._minNode = null; + this._nodeCount = 0; + } + + /** + * Decreases a key of a node. + * @param node The node to decrease the key of. + * @param newKey The new key to assign to the node. + */ + public decreaseKey(node: Node<K, V>, newKey: K): void { + if (!node) { + throw new Error('Cannot decrease key of non-existent node'); + } + if (this._compare({ key: newKey }, { key: node.key }) > 0) { + throw new Error('New key is larger than old key'); + } + + node.key = newKey; + const parent = node.parent; + if (parent && this._compare(node, parent) < 0) { + this._cut(node, parent, <Node<K, V>> this._minNode); + this._cascadingCut(parent, <Node<K, V>> this._minNode); + } + if (this._compare(node, <Node<K, V>> this._minNode) < 0) { + this._minNode = node; + } + } + + /** + * Deletes a node. + * @param node The node to delete. + */ + public delete(node: Node<K, V>): void { + // This is a special implementation of decreaseKey that sets the argument to + // the minimum value. This is necessary to make generic keys work, since there + // is no MIN_VALUE constant for generic types. + const parent = node.parent; + if (parent) { + this._cut(node, parent, <Node<K, V>> this._minNode); + this._cascadingCut(parent, <Node<K, V>> this._minNode); + } + this._minNode = node; + + this.extractMinimum(); + } + + /** + * Extracts and returns the minimum node from the heap. + * @return The heap's minimum node or null if the heap is empty. + */ + public extractMinimum(): Node<K, V> | null { + const extractedMin = this._minNode; + if (extractedMin) { + // Set parent to null for the minimum's children + if (extractedMin.child) { + let child = extractedMin.child; + do { + child.parent = null; + child = child.next; + } while (child !== extractedMin.child); + } + + let nextInRootList = null; + if (extractedMin.next !== extractedMin) { + nextInRootList = extractedMin.next; + } + // Remove min from root list + this._removeNodeFromList(extractedMin); + this._nodeCount--; + + // Merge the children of the minimum node with the root list + this._minNode = this._mergeLists(nextInRootList, extractedMin.child); + if (this._minNode) { + this._minNode = this._consolidate(this._minNode); + } + } + return extractedMin; + } + + /** + * Returns the minimum node from the heap. + * @return The heap's minimum node or null if the heap is empty. + */ + public findMinimum(): Node<K, V> | null { + return this._minNode; + } + + /** + * Inserts a new key-value pair into the heap. + * @param key The key to insert. + * @param value The value to insert. + * @return node The inserted node. + */ + public insert(key: K, value?: V): Node<K, V> { + const node = new Node(key, value); + this._minNode = this._mergeLists(this._minNode, node); + this._nodeCount++; + return node; + } + + /** + * @return Whether the heap is empty. + */ + public isEmpty(): boolean { + return this._minNode === null; + } + + /** + * @return The size of the heap. + */ + public size(): number { + if (this._minNode === null) { + return 0; + } + return this._getNodeListSize(this._minNode); + } + + /** + * Joins another heap to this heap. + * @param other The other heap. + */ + public union(other: FibonacciHeap<K, V>): void { + this._minNode = this._mergeLists(this._minNode, other._minNode); + this._nodeCount += other._nodeCount; + } + + /** + * Compares two nodes with each other. + * @param a The first key to compare. + * @param b The second key to compare. + * @return -1, 0 or 1 if a < b, a == b or a > b respectively. + */ + private _defaultCompare(a: INode<K, V>, b: INode<K, V>): number { + if (a.key > b.key) { + return 1; + } + if (a.key < b.key) { + return -1; + } + return 0; + } + + /** + * Cut the link between a node and its parent, moving the node to the root list. + * @param node The node being cut. + * @param parent The parent of the node being cut. + * @param minNode The minimum node in the root list. + * @return The heap's new minimum node. + */ + private _cut(node: Node<K, V>, parent: Node<K, V>, minNode: Node<K, V>): Node<K, V> | null { + node.parent = null; + parent.degree--; + if (node.next === node) { + parent.child = null; + } else { + parent.child = node.next; + } + this._removeNodeFromList(node); + const newMinNode = this._mergeLists(minNode, node); + node.isMarked = false; + return newMinNode; + } + + /** + * Perform a cascading cut on a node; mark the node if it is not marked, + * otherwise cut the node and perform a cascading cut on its parent. + * @param node The node being considered to be cut. + * @param minNode The minimum node in the root list. + * @return The heap's new minimum node. + */ + private _cascadingCut(node: Node<K, V>, minNode: Node<K, V> | null): Node<K, V> | null { + const parent = node.parent; + if (parent) { + if (node.isMarked) { + minNode = this._cut(node, parent, <Node<K, V>>minNode); + minNode = this._cascadingCut(parent, minNode); + } else { + node.isMarked = true; + } + } + return minNode; + } + + /** + * Merge all trees of the same order together until there are no two trees of + * the same order. + * @param minNode The current minimum node. + * @return The new minimum node. + */ + private _consolidate(minNode: Node<K, V>): Node<K, V> | null { + + const aux = []; + tmpIt.reset(minNode); + while (tmpIt.hasNext()) { + let current = tmpIt.next(); + + // If there exists another node with the same degree, merge them + let auxCurrent = aux[current.degree]; + while (auxCurrent) { + if (this._compare(current, auxCurrent) > 0) { + const temp = current; + current = auxCurrent; + auxCurrent = temp; + } + this._linkHeaps(auxCurrent, current); + aux[current.degree] = null; + current.degree++; + auxCurrent = aux[current.degree]; + } + + aux[current.degree] = current; + } + + let newMinNode = null; + for (let i = 0; i < aux.length; i++) { + const node = aux[i]; + if (node) { + // Remove siblings before merging + node.next = node; + node.prev = node; + newMinNode = this._mergeLists(newMinNode, node); + } + } + return newMinNode; + } + + /** + * Removes a node from a node list. + * @param node The node to remove. + */ + private _removeNodeFromList(node: Node<K, V>): void { + const prev = node.prev; + const next = node.next; + prev.next = next; + next.prev = prev; + node.next = node; + node.prev = node; + } + + /** + * Links two heaps of the same order together. + * + * @private + * @param max The heap with the larger root. + * @param min The heap with the smaller root. + */ + private _linkHeaps(max: Node<K, V>, min: Node<K, V>): void { + this._removeNodeFromList(max); + min.child = this._mergeLists(max, min.child); + max.parent = min; + max.isMarked = false; + } + + /** + * Merge two lists of nodes together. + * + * @private + * @param a The first list to merge. + * @param b The second list to merge. + * @return The new minimum node from the two lists. + */ + private _mergeLists(a: Node<K, V> | null, b: Node<K, V> | null): Node<K, V> | null { + if (!a) { + if (!b) { + return null; + } + return b; + } + if (!b) { + return a; + } + + const temp = a.next; + a.next = b.next; + a.next.prev = a; + b.next = temp; + b.next.prev = b; + + return this._compare(a, b) < 0 ? a : b; + } + + /** + * Gets the size of a node list. + * @param node A node within the node list. + * @return The size of the node list. + */ + private _getNodeListSize(node: Node<K, V>): number { + let count = 0; + let current = node; + + do { + count++; + if (current.child) { + count += this._getNodeListSize(current.child); + } + current = current.next; + } while (current !== node); + + return count; + } +}