Volseg extenison tweaks (#671)

* improve currentVolume behavior

* use different iso-value if 1sigma isnt present in the data

src/extensions/volumes-and-segmentations/entry-root.ts

@@ -15,7 +15,7 @@ import { PluginStateObject } from '../../mol-plugin-state/objects';
 import { PluginBehavior } from '../../mol-plugin/behavior';
 import { PluginCommands } from '../../mol-plugin/commands';
 import { PluginContext } from '../../mol-plugin/context';
-import { StateObjectCell, StateTransform } from '../../mol-state';
+import { StateObjectCell, StateSelection, StateTransform } from '../../mol-state';
 import { shallowEqualObjects } from '../../mol-util';
 import { ParamDefinition } from '../../mol-util/param-definition';
 import { MeshlistData } from '../meshes/mesh-extension';
@@ -151,6 +151,7 @@ export class VolsegEntryData extends PluginBehavior.WithSubscribers<VolsegEntryP
             // do nothing
         }
 
+        let volumeRef: string | undefined;
         this.subscribeObservable(this.plugin.state.data.events.cell.stateUpdated, e => {
             try { (this.getStateNode()); } catch { return; } // if state not does not exist yet
             if (e.cell.transform.ref === this.getStateNode().transform.ref) {
@@ -158,14 +159,20 @@ export class VolsegEntryData extends PluginBehavior.WithSubscribers<VolsegEntryP
                 if (newState && !shallowEqualObjects(newState, this.currentState.value)) { // avoid repeated update
                     this.currentState.next(newState);
                 }
+            } else if (e.cell.transform.tags?.includes(VOLUME_VISUAL_TAG)) {
+                if (e.ref === volumeRef) {
+                    this.currentVolume.next(e.cell.transform);
+                } else if (StateSelection.findAncestor(this.plugin.state.data.tree, this.plugin.state.data.cells, e.ref, a => a.transform.ref === ref)) {
+                    volumeRef = e.ref;
+                    this.currentVolume.next(e.cell.transform);
+                }
             }
         });
 
-        this.subscribeObservable(this.plugin.state.data.events.cell.stateUpdated, e => {
-            // TODO: subscribe cell.removed event as well to set the current volume to undefined
-            if (e.cell.transform.tags?.includes(VOLUME_VISUAL_TAG)) {
-                // TODO: make sure this belongs to the "current entry subtree"
-                this.currentVolume.next(e.cell.transform);
+        this.subscribeObservable(this.plugin.state.data.events.cell.removed, e => {
+            if (e.ref === volumeRef) {
+                volumeRef = undefined;
+                this.currentVolume.next(undefined);
             }
         });
 

src/extensions/volumes-and-segmentations/entry-volume.ts

@@ -52,7 +52,7 @@ export class VolsegVolumeData {
     async loadVolume() {
         const hasVolumes = this.entryData.metadata.raw.grid.volumes.volume_downsamplings.length > 0;
         if (hasVolumes) {
-            const isoLevelPromise = ExternalAPIs.getIsovalue(this.entryData.metadata.raw.grid.general.source_db_id ?? this.entryData.entryId);
+            const isoLevelPromise = ExternalAPIs.tryGetIsovalue(this.entryData.metadata.raw.grid.general.source_db_id ?? this.entryData.entryId);
             let group = this.entryData.findNodesByTags(GROUP_TAG)[0]?.transform.ref;
             if (!group) {
                 const newGroupNode = await this.entryData.newUpdate().apply(CreateGroup, { label: 'Volume' }, { tags: [GROUP_TAG], state: { isCollapsed: true } }).commit();
@@ -65,7 +65,17 @@ export class VolsegVolumeData {
             const volumeData = volumeNode.cell!.obj!.data;
 
             const volumeType = VolsegStateParams.volumeType.defaultValue;
-            const isovalue = await isoLevelPromise;
+            let isovalue = await isoLevelPromise;
+            if (!isovalue) {
+                const stats = volumeData.grid.stats;
+                const maxRelative = (stats.max - stats.mean) / stats.sigma;
+                if (maxRelative > 1) {
+                    isovalue = { kind: 'relative', value: 1.0 };
+                } else {
+                    isovalue = { kind: 'relative', value: maxRelative * 0.5 };
+                }
+            }
+
             const adjustedIsovalue = Volume.adjustedIsoValue(volumeData, isovalue.value, isovalue.kind);
             const visualParams = this.createVolumeVisualParams(volumeData, volumeType);
             this.changeIsovalueInVolumeVisualParams(visualParams, adjustedIsovalue, volumeData.grid.stats);

src/extensions/volumes-and-segmentations/external-api.ts

@@ -8,7 +8,7 @@ import { splitEntryId } from './helpers';
 
 
 /** Try to get author-defined contour value for isosurface from EMDB API. Return relative value 1.0, if not applicable or fails.  */
-export async function getIsovalue(entryId: string): Promise<{ kind: 'absolute' | 'relative', value: number }> {
+export async function tryGetIsovalue(entryId: string): Promise<{ kind: 'absolute' | 'relative', value: number } | undefined> {
     const split = splitEntryId(entryId);
     if (split.source === 'emdb') {
         try {
@@ -24,7 +24,7 @@ export async function getIsovalue(entryId: string): Promise<{ kind: 'absolute' |
             // do nothing
         }
     }
-    return { kind: 'relative', value: 1.0 };
+    return undefined;
 }
 
 export async function getPdbIdsForEmdbEntry(entryId: string): Promise<string[]> {

src/extensions/volumes-and-segmentations/lattice-segmentation.ts

@@ -1,274 +0,0 @@
-/**
- * Copyright (c) 2018-2022 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Adam Midlik <midlik@gmail.com>
- */
-
-import { CIF, CifBlock } from '../../mol-io/reader/cif';
-import { Box3D } from '../../mol-math/geometry';
-import { Tensor, Vec3 } from '../../mol-math/linear-algebra';
-import { volumeFromDensityServerData } from '../../mol-model-formats/volume/density-server';
-import { CustomProperties } from '../../mol-model/custom-property';
-import { Grid, Volume } from '../../mol-model/volume';
-import { Segment } from './volseg-api/data';
-import { lazyGetter } from './helpers';
-
-
-export class LatticeSegmentation {
-    private segments: number[];
-    private sets: number[];
-    /** Maps setId to a set of segmentIds*/
-    private segmentMap: Map<number, Set<number>>; // computations with objects might be actually faster than with Maps and Sets?
-    /** Maps segmentId to a set of setIds*/
-    private inverseSegmentMap: Map<number, Set<number>>;
-    private grid: Grid;
-
-    private constructor(segmentationDataBlock: CifBlock, grid: Grid) {
-        const segmentationValues = segmentationDataBlock!.categories['segmentation_data_3d'].getField('values')?.toIntArray()!;
-        this.segmentMap = LatticeSegmentation.makeSegmentMap(segmentationDataBlock);
-        this.inverseSegmentMap = LatticeSegmentation.invertMultimap(this.segmentMap);
-        this.segments = Array.from(this.inverseSegmentMap.keys());
-        this.sets = Array.from(this.segmentMap.keys());
-        this.grid = grid;
-        this.grid.cells.data = Tensor.Data1(segmentationValues);
-    }
-
-    public static async fromCifBlock(segmentationDataBlock: CifBlock) {
-        const densityServerCif = CIF.schema.densityServer(segmentationDataBlock);
-        const volume = await volumeFromDensityServerData(densityServerCif).run();
-        const grid = volume.grid;
-        return new LatticeSegmentation(segmentationDataBlock, grid);
-    }
-
-    public createSegment_old(segId: number): Volume {
-        // console.time('createSegment_old');
-        const n = this.grid.cells.data.length;
-        const newData = new Float32Array(n);
-
-        for (let i = 0; i < n; i++) {
-            newData[i] = this.segmentMap.get(this.grid.cells.data[i])?.has(segId) ? 1 : 0;
-        }
-
-        const result: Volume = {
-            sourceData: { kind: 'custom', name: 'test', data: newData as any },
-            customProperties: new CustomProperties(),
-            _propertyData: {},
-            grid: {
-                ...this.grid,
-                // stats: { min: 0, max: 1, mean: newMean, sigma: arrayRms(newData) },
-                stats: { min: 0, max: 1, mean: 0, sigma: 1 },
-                cells: {
-                    ...this.grid.cells,
-                    data: newData as any,
-                }
-            }
-        };
-        // console.timeEnd('createSegment_old');
-        return result;
-    }
-
-    public hasSegment(segId: number): boolean {
-        return this.inverseSegmentMap.has(segId);
-    }
-    public createSegment(seg: Segment, propertyData?: {[key: string]: any}): Volume {
-        const { space, data }: Tensor = this.grid.cells;
-        const [nx, ny, nz] = space.dimensions;
-        const axisOrder = [...space.axisOrderSlowToFast];
-        const get = space.get;
-        const cell = Box.create(0, nx, 0, ny, 0, nz);
-
-        const EXPAND_START = 2; // We need to add 2 layers of zeros, probably because of a bug in GPU marching cubes implementation
-        const EXPAND_END = 1;
-        let bbox = this.getSegmentBoundingBoxes()[seg.id];
-        bbox = Box.expand(bbox, EXPAND_START, EXPAND_END);
-        bbox = Box.confine(bbox, cell);
-        const [ox, oy, oz] = Box.origin(bbox);
-        const [mx, my, mz] = Box.size(bbox);
-        // n, i refer to original box; m, j to the new box
-
-        const newSpace = Tensor.Space([mx, my, mz], axisOrder, Uint8Array);
-        const newTensor = Tensor.create(newSpace, newSpace.create());
-        const newData = newTensor.data;
-        const newSet = newSpace.set;
-
-        const sets = this.inverseSegmentMap.get(seg.id);
-        if (!sets) throw new Error(`This LatticeSegmentation does not contain segment ${seg.id}`);
-
-        for (let jz = 0; jz < mz; jz++) {
-            for (let jy = 0; jy < my; jy++) {
-                for (let jx = 0; jx < mx; jx++) {
-                    // Iterating in ZYX order is faster (probably fewer cache misses)
-                    const setId = get(data, ox + jx, oy + jy, oz + jz);
-                    const value = sets.has(setId) ? 1 : 0;
-                    newSet(newData, jx, jy, jz, value);
-                }
-            }
-        }
-
-        const transform = this.grid.transform;
-        let newTransform: Grid.Transform;
-        if (transform.kind === 'matrix') {
-            throw new Error('Not implemented for transform of kind "matrix"'); // TODO ask if this is really needed
-        } else if (transform.kind === 'spacegroup') {
-            const newFractionalBox = Box.toFractional(bbox, cell);
-            const origFractSize = Vec3.sub(Vec3.zero(), transform.fractionalBox.max, transform.fractionalBox.min);
-            Vec3.mul(newFractionalBox.min, newFractionalBox.min, origFractSize);
-            Vec3.mul(newFractionalBox.max, newFractionalBox.max, origFractSize);
-            Vec3.add(newFractionalBox.min, newFractionalBox.min, transform.fractionalBox.min);
-            Vec3.add(newFractionalBox.max, newFractionalBox.max, transform.fractionalBox.min);
-            newTransform = { ...transform, fractionalBox: newFractionalBox };
-        } else {
-            throw new Error(`Unknown transform kind: ${transform}`);
-        }
-        const result = {
-            sourceData: { kind: 'custom', name: 'test', data: newTensor.data as any },
-            label: seg.biological_annotation.name ?? `Segment ${seg.id}`,
-            customProperties: new CustomProperties(),
-            _propertyData: propertyData ?? {},
-            grid: {
-                stats: { min: 0, max: 1, mean: 0, sigma: 1 },
-                cells: newTensor,
-                transform: newTransform,
-            }
-        } as Volume;
-        return result;
-    }
-
-    private static _getSegmentBoundingBoxes(self: LatticeSegmentation) {
-        const { space, data }: Tensor = self.grid.cells;
-        const [nx, ny, nz] = space.dimensions;
-        const get = space.get;
-
-        const setBoxes: { [setId: number]: Box } = {}; // with object this is faster than with Map
-        self.sets.forEach(setId => setBoxes[setId] = Box.create(nx, -1, ny, -1, nz, -1));
-
-        for (let iz = 0; iz < nz; iz++) {
-            for (let iy = 0; iy < ny; iy++) {
-                for (let ix = 0; ix < nx; ix++) {
-                    // Iterating in ZYX order is faster (probably fewer cache misses)
-                    const setId = get(data, ix, iy, iz);
-                    Box.addPoint_InclusiveEnd(setBoxes[setId], ix, iy, iz);
-                }
-            }
-        }
-
-        const segmentBoxes: { [segmentId: number]: Box } = {};
-        self.segments.forEach(segmentId => segmentBoxes[segmentId] = Box.create(nx, -1, ny, -1, nz, -1));
-        self.inverseSegmentMap.forEach((setIds, segmentId) => {
-            setIds.forEach(setId => {
-                segmentBoxes[segmentId] = Box.cover(segmentBoxes[segmentId], setBoxes[setId]);
-            });
-        });
-
-        for (const segmentId in segmentBoxes) {
-            if (segmentBoxes[segmentId][5] === -1) { // segment's box left unchanged -> contains no voxels
-                segmentBoxes[segmentId] = Box.create(0, 1, 0, 1, 0, 1);
-            } else {
-                segmentBoxes[segmentId] = Box.expand(segmentBoxes[segmentId], 0, 1); // inclusive end -> exclusive end
-            }
-        }
-        return segmentBoxes;
-    }
-    private getSegmentBoundingBoxes = lazyGetter(() => LatticeSegmentation._getSegmentBoundingBoxes(this));
-
-    private static invertMultimap<K, V>(map: Map<K, Set<V>>): Map<V, Set<K>> {
-        const inverted = new Map<V, Set<K>>();
-        map.forEach((values, key) => {
-            values.forEach(value => {
-                if (!inverted.has(value)) inverted.set(value, new Set<K>());
-                inverted.get(value)?.add(key);
-            });
-        });
-        return inverted;
-    }
-
-    private static makeSegmentMap(segmentationDataBlock: CifBlock): Map<number, Set<number>> {
-        const setId = segmentationDataBlock.categories['segmentation_data_table'].getField('set_id')?.toIntArray()!;
-        const segmentId = segmentationDataBlock.categories['segmentation_data_table'].getField('segment_id')?.toIntArray()!;
-        const map = new Map<number, Set<number>>();
-        for (let i = 0; i < segmentId.length; i++) {
-            if (!map.has(setId[i])) {
-                map.set(setId[i], new Set());
-            }
-            map.get(setId[i])!.add(segmentId[i]);
-        }
-        return map;
-    }
-
-    public benchmark(segment: Segment) {
-        const N = 100;
-
-        console.time(`createSegment ${segment.id} ${N}x`);
-        for (let i = 0; i < N; i++) {
-            this.getSegmentBoundingBoxes = lazyGetter(() => LatticeSegmentation._getSegmentBoundingBoxes(this));
-            this.createSegment(segment);
-        }
-        console.timeEnd(`createSegment ${segment.id} ${N}x`);
-    }
-}
-
-
-type Box = [number, number, number, number, number, number];
-
-/** Represents a 3D box in integer coordinates. xFrom... is inclusive, xTo... is exclusive. */
-namespace Box {
-    export function create(xFrom: number, xTo: number, yFrom: number, yTo: number, zFrom: number, zTo: number): Box {
-        return [xFrom, xTo, yFrom, yTo, zFrom, zTo];
-    }
-    export function expand(box: Box, expandFrom: number, expandTo: number): Box {
-        const [xFrom, xTo, yFrom, yTo, zFrom, zTo] = box;
-        return [xFrom - expandFrom, xTo + expandTo, yFrom - expandFrom, yTo + expandTo, zFrom - expandFrom, zTo + expandTo];
-    }
-    export function confine(box1: Box, box2: Box): Box {
-        const [xFrom1, xTo1, yFrom1, yTo1, zFrom1, zTo1] = box1;
-        const [xFrom2, xTo2, yFrom2, yTo2, zFrom2, zTo2] = box2;
-        return [
-            Math.max(xFrom1, xFrom2), Math.min(xTo1, xTo2),
-            Math.max(yFrom1, yFrom2), Math.min(yTo1, yTo2),
-            Math.max(zFrom1, zFrom2), Math.min(zTo1, zTo2)
-        ];
-    }
-    export function cover(box1: Box, box2: Box): Box {
-        const [xFrom1, xTo1, yFrom1, yTo1, zFrom1, zTo1] = box1;
-        const [xFrom2, xTo2, yFrom2, yTo2, zFrom2, zTo2] = box2;
-        return [
-            Math.min(xFrom1, xFrom2), Math.max(xTo1, xTo2),
-            Math.min(yFrom1, yFrom2), Math.max(yTo1, yTo2),
-            Math.min(zFrom1, zFrom2), Math.max(zTo1, zTo2)
-        ];
-    }
-    export function size(box: Box): [number, number, number] {
-        const [xFrom, xTo, yFrom, yTo, zFrom, zTo] = box;
-        return [xTo - xFrom, yTo - yFrom, zTo - zFrom];
-    }
-    export function origin(box: Box): [number, number, number] {
-        const xFrom = box[0];
-        const yFrom = box[2];
-        const zFrom = box[4];
-        return [xFrom, yFrom, zFrom];
-    }
-    export function log(name: string, box: Box): void {
-        const [xFrom, xTo, yFrom, yTo, zFrom, zTo] = box;
-        console.log(`Box ${name}: [${xFrom}:${xTo}, ${yFrom}:${yTo}, ${zFrom}:${zTo}], size: ${size(box)}`);
-    }
-    export function toFractional(box: Box, relativeTo: Box): Box3D {
-        const [xFrom, xTo, yFrom, yTo, zFrom, zTo] = box;
-        const [x0, y0, z0] = origin(relativeTo);
-        const [sizeX, sizeY, sizeZ] = size(relativeTo);
-        const min = Vec3.create((xFrom - x0) / sizeX, (yFrom - y0) / sizeY, (zFrom - z0) / sizeZ);
-        const max = Vec3.create((xTo - x0) / sizeX, (yTo - y0) / sizeY, (zTo - z0) / sizeZ);
-        return Box3D.create(min, max);
-    }
-    export function addPoint_InclusiveEnd(box: Box, x: number, y: number, z: number): void {
-        if (x < box[0]) box[0] = x;
-        if (x > box[1]) box[1] = x;
-        if (y < box[2]) box[2] = y;
-        if (y > box[3]) box[3] = y;
-        if (z < box[4]) box[4] = z;
-        if (z > box[5]) box[5] = z;
-    }
-    export function equal(box1: Box, box2: Box): boolean {
-        return box1.every((value, i) => value === box2[i]);
-    }
-}
-