View 2D: IJK slicing

packages/element/examples/image-io-read-image.ts

@@ -12,7 +12,7 @@ document.addEventListener('DOMContentLoaded', async function () {
   const viewerElement = document.querySelector('#viewer')! as ItkViewer2d;
   const viewer = viewerElement.getActor();
 
-  const path = 'HeadMRVolume.nrrd';
+  const path = 'prostate-slight-rotation.nrrd';
   const url = new URL(path, document.location.origin);
   const response = await fetch(url.href);
   const data = new Uint8Array(await response.arrayBuffer());

packages/io/src/MultiscaleSpatialImage.ts

@@ -23,6 +23,7 @@ import {
 } from './dimensionUtils.js';
 import { transformBounds } from './transformBounds.js';
 import {
+  ChunkParameter,
   Dimension,
   Extent,
   ReadOnlyDimensionBounds,
@@ -42,6 +43,11 @@ function setMatrixElement(
   matrixData[column + row * columns] = value;
 }
 
+type ReadonlyFloatArray =
+  | Readonly<Float32Array>
+  | Readonly<Float64Array>
+  | ReadonlyMat3;
+
 type ImageDataFromChunksWorkerArgs = {
   scaleInfo: {
     chunkSize: Map<Dimension, number>;
@@ -80,7 +86,8 @@ function inflate(bounds: Bounds, delta: number) {
 }
 
 // code modified from vtk.js/ImageData
-const extentToBounds = (ex: Extent, indexToWorld: mat4) => {
+const extentToBounds = (extent: Extent, indexToWorld: ReadonlyMat4) => {
+  const ex = extent;
   // prettier-ignore
   const corners = [
     ex[0], ex[2], ex[4],
@@ -129,7 +136,11 @@ const extentToBounds = (ex: Extent, indexToWorld: mat4) => {
   return bounds;
 };
 
-export const ensure3dDirection = (d: Float64Array): ReadonlyMat3 => {
+// returns a copy
+export const ensure3dDirection = (
+  maybe2dDirection: ReadonlyFloatArray,
+): mat3 => {
+  const d = maybe2dDirection;
   if (d.length >= 9) {
     return mat3.fromValues(
       d[0],
@@ -172,21 +183,28 @@ const makeMat4 = ({
   return mat4.scale(mat, mat, spacing);
 };
 
-const makeIndexToWorld = ({
+type Maybe2dVec =
+  | [number, number]
+  | [number, number, number]
+  | vec3
+  | Array<number>;
+
+const maybe2dSpatialToMat4 = ({
   direction: inDirection,
   origin,
   spacing,
 }: {
-  direction: ReadonlyMat3;
-  origin: Array<number>;
-  spacing: Array<number>;
+  direction: ReadonlyFloatArray;
+  origin: Maybe2dVec;
+  spacing: Maybe2dVec;
 }) => {
+  const inDirection3d = ensure3dDirection(inDirection);
   // ITK (and VTKMath) uses row-major index axis, but gl-matrix uses column-major. Transpose.
   const DIMENSIONS = 3;
-  const direction = Array(inDirection.length) as mat3;
+  const direction = Array(inDirection3d.length) as mat3;
   for (let idx = 0; idx < DIMENSIONS; ++idx) {
     for (let col = 0; col < DIMENSIONS; ++col) {
-      direction[col + idx * 3] = inDirection[idx + col * DIMENSIONS];
+      direction[col + idx * 3] = inDirection3d[idx + col * DIMENSIONS];
     }
   }
 
@@ -239,6 +257,53 @@ export const worldBoundsToIndexBounds = ({
   return new Map([...spaceBounds, ...ctBounds]);
 };
 
+// Ensures CXYZT dimensions are present
+const ensureBoundsCXYZT = ({
+  indexBounds,
+  fullIndexBounds,
+}: {
+  indexBounds: ReadonlyBounds;
+  fullIndexBounds: ReadOnlyDimensionBounds;
+}) => {
+  const fullIndexBoundsWithZCT = ensuredDims(
+    [0, 1] as [number, number],
+    CXYZT,
+    fullIndexBounds,
+  );
+  // clamp to existing integer indexes
+  const imageBoundsByDim = chunk(2, indexBounds);
+  const spaceBounds = (['x', 'y', 'z'] as const).map((dim, idx) => {
+    // eslint-disable-next-line @typescript-eslint/no-non-null-assertion
+    const [min, max] = fullIndexBoundsWithZCT.get(dim)!;
+    const [bmin, bmax] = imageBoundsByDim[idx];
+    return [
+      dim,
+      [
+        Math.floor(Math.min(max, Math.max(min, bmin))),
+        Math.ceil(Math.min(max, Math.max(min, bmax))),
+      ],
+    ] as const;
+  });
+  const ctBounds = (['c', 't'] as const).map(
+    (dim) => [dim, fullIndexBoundsWithZCT.get(dim)!] as const,
+  );
+  return new Map([...spaceBounds, ...ctBounds]);
+};
+
+const normalizedImageBoundsToIndexBounds = (
+  arrayShape: ChunkParameter,
+  normalizedImageBounds: ReadonlyBounds,
+) => {
+  const toIndexScale = XYZ.map((axis) => arrayShape.get(axis) ?? 1);
+  const normalizedToIndex = maybe2dSpatialToMat4({
+    direction: mat3.identity(mat3.create()),
+    origin: vec3.create(),
+    spacing: toIndexScale,
+  }) as ReadonlyMat4;
+  const indexBounds = transformBounds(normalizedToIndex, normalizedImageBounds);
+  return indexBounds;
+};
+
 function isContained(
   benchmarkBounds: ReadOnlyDimensionBounds,
   testedBounds: ReadOnlyDimensionBounds,
@@ -327,7 +392,7 @@ export class MultiscaleSpatialImage {
     const info = this.scaleInfos[scale];
     if (info.origin) return info.origin;
 
-    const origin = new Array(this.spatialDims.length);
+    const origin = new Array<number>(this.spatialDims.length);
     for (let index = 0; index < this.spatialDims.length; index++) {
       const dim = this.spatialDims[index];
       if (info.coords.has(dim)) {
@@ -346,7 +411,7 @@ export class MultiscaleSpatialImage {
     const info = this.scaleInfos[scale];
     if (info.spacing) return info.spacing;
 
-    const spacing = new Array(this.spatialDims.length);
+    const spacing = new Array<number>(this.spatialDims.length);
     for (let index = 0; index < this.spatialDims.length; index++) {
       const dim = this.spatialDims[index];
       const dimCoords = await info.coords.get(dim);
@@ -389,6 +454,7 @@ export class MultiscaleSpatialImage {
         setMatrixElement(direction, dimension, d, d, 1.0);
       }
     }
+
     return direction;
   }
 
@@ -497,25 +563,47 @@ export class MultiscaleSpatialImage {
   async scaleIndexToWorld(requestedScale: number) {
     const scale = Math.min(requestedScale, this.scaleInfos.length - 1);
     if (this.scaleInfos[scale].indexToWorld)
-      return this.scaleInfos[scale].indexToWorld as mat4;
+      return this.scaleInfos[scale].indexToWorld as ReadonlyMat4;
 
     // compute and cache origin/scale on info
     const [origin, spacing] = await Promise.all([
       this.scaleOrigin(scale),
       this.scaleSpacing(scale),
     ]);
 
-    const direction = ensure3dDirection(this.direction);
-
-    const indexToWorld = makeIndexToWorld({
-      direction,
+    const indexToWorld = maybe2dSpatialToMat4({
+      direction: this.direction,
       origin,
       spacing,
-    });
+    }) as ReadonlyMat4;
     this.scaleInfos[scale].indexToWorld = indexToWorld;
     return indexToWorld;
   }
 
+  async buildAndCacheImage(scale: number, indexBounds: Bounds) {
+    const indexBoundsCXYZT = ensureBoundsCXYZT({
+      indexBounds,
+      fullIndexBounds: this.getIndexBounds(scale),
+    });
+    const { dims } = this.scaleInfos[scale];
+    const indexBoundsByDimension = orderBy(dims)(indexBoundsCXYZT);
+    const cachedImage = findImageInBounds({
+      cache: this.cachedImages,
+      scale,
+      bounds: indexBoundsByDimension,
+    });
+    if (cachedImage) return cachedImage;
+
+    const image = await this.buildImage(scale, indexBoundsByDimension);
+    storeImage({
+      cache: this.cachedImages,
+      scale,
+      bounds: indexBoundsByDimension,
+      image,
+    });
+    return image;
+  }
+
   /* Retrieve bounded image at scale. */
   async getImage(
     requestedScale: number,
@@ -524,25 +612,29 @@ export class MultiscaleSpatialImage {
     const scale = Math.min(requestedScale, this.scaleInfos.length - 1);
     const indexToWorld = await this.scaleIndexToWorld(scale);
 
-    const { dims } = this.scaleInfos[scale];
-    const indexBounds = orderBy(dims)(
-      worldBoundsToIndexBounds({
-        bounds: worldBounds,
-        fullIndexBounds: this.getIndexBounds(scale),
-        worldToIndex: mat4.invert(mat4.create(), indexToWorld),
-      }),
+    const fullIndexBounds = ensuredDims(
+      [0, 0],
+      XYZ,
+      this.getIndexBounds(scale),
     );
+    let indexBounds = XYZ.flatMap((dim) => fullIndexBounds.get(dim)) as Bounds;
+    if (worldBounds) {
+      const worldToIndex = mat4.invert(mat4.create(), indexToWorld);
+      indexBounds = transformBounds(worldToIndex, worldBounds);
+    }
 
-    const cachedImage = findImageInBounds({
-      cache: this.cachedImages,
-      scale,
-      bounds: indexBounds,
-    });
-    if (cachedImage) return cachedImage;
+    return this.buildAndCacheImage(scale, indexBounds);
+  }
 
-    const image = await this.buildImage(scale, indexBounds);
-    storeImage({ cache: this.cachedImages, scale, bounds: indexBounds, image });
-    return image;
+  async getImageInImageSpace(
+    scale: number,
+    normalizedImageBounds: ReadonlyBounds = [0, 1, 0, 1, 0, 1],
+  ) {
+    const indexBounds = normalizedImageBoundsToIndexBounds(
+      this.scaleInfos[scale].arrayShape,
+      normalizedImageBounds,
+    );
+    return this.buildAndCacheImage(scale, indexBounds);
   }
 
   getIndexBounds(scale: number) {
@@ -555,8 +647,7 @@ export class MultiscaleSpatialImage {
     );
   }
 
-  async getWorldBounds(scale: number) {
-    const indexToWorld = await this.scaleIndexToWorld(scale);
+  getIndexExtent(scale: number) {
     const imageBounds = ensuredDims(
       [0, 1],
       ['x', 'y', 'z'],
@@ -566,6 +657,12 @@ export class MultiscaleSpatialImage {
       imageBounds.get(dim),
     ) as Bounds;
     inflate(bounds, 0.5);
+    return bounds;
+  }
+
+  async getWorldBounds(scale: number) {
+    const indexToWorld = await this.scaleIndexToWorld(scale);
+    const bounds = this.getIndexExtent(scale);
     return extentToBounds(bounds, indexToWorld);
   }
 }

packages/io/src/dimensionUtils.ts

@@ -1,6 +1,7 @@
 import { Dimension, SpatialDimensions } from './types.js';
 
 export const CXYZT = Object.freeze(['c', 'x', 'y', 'z', 't'] as const); // viewer indexing
+export const XYZ = Object.freeze(['x', 'y', 'z']) as SpatialDimensions;
 
 export const ensuredDims = <T>(
   defaultValue: T,
@@ -30,7 +31,7 @@ export const orderBy =
       }),
     );
 
-export const chunk = <T>(chunkSize: number, array: Array<T>) => {
+export const chunk = <T>(chunkSize: number, array: ReadonlyArray<T>) => {
   const chunks = [];
   for (let i = 0; i < array.length; i += chunkSize) {
     chunks.push(array.slice(i, i + chunkSize));
@@ -41,5 +42,3 @@ export const chunk = <T>(chunkSize: number, array: Array<T>) => {
 export const nonNullable = <T>(value: T): value is NonNullable<T> => {
   return value != null;
 };
-
-export const XYZ = Object.freeze(['x', 'y', 'z']) as SpatialDimensions;

packages/io/src/types.ts

@@ -1,4 +1,4 @@
-import { mat4 } from 'gl-matrix';
+import { ReadonlyMat4 } from 'gl-matrix';
 import { TypedArray } from 'itk-wasm';
 
 export type ValueOf<T> = T[keyof T];
@@ -21,7 +21,7 @@ export type TypedArrayConstructor = Constructor<TypedArray>;
 
 export type SpatialDimensions = ['x'] | ['x', 'y'] | ['x', 'y', 'z'];
 
-type ChunkParameter = Map<Dimension, number>;
+export type ChunkParameter = Map<Dimension, number>;
 type LazyCoords = {
   get: (dim: Dimension) => Awaitable<Float32Array | undefined>;
   has: (dim: Dimension) => boolean;
@@ -40,7 +40,7 @@ export type ScaleInfo = {
   ranges?: Array<Array<number>>; // Map('1': [0, 140], '2': [3, 130]) or null if unknown. Range of values for each component
   name?: string; // dataset name
 
-  indexToWorld?: mat4;
+  indexToWorld?: ReadonlyMat4;
 
   // Zarr specific
   pixelArrayMetadata?: {

packages/viewer/src/camera.ts

@@ -2,7 +2,6 @@ import {
   Bounds,
   addPoint,
   createBounds,
-  getCorners,
   getLength,
 } from '@itk-viewer/utils/bounding-box.js';
 import { ReadonlyVec3, mat4, vec3, quat, ReadonlyQuat } from 'gl-matrix';
@@ -209,24 +208,29 @@ export const reset3d = (
 export const reset2d = (
   pose: Pose,
   verticalFieldOfView: number,
-  bounds: Bounds,
+  pointsToFit: Array<ReadonlyVec3>,
   aspect: number,
 ) => {
-  const center = vec3.fromValues(
-    (bounds[0] + bounds[1]) / 2.0,
-    (bounds[2] + bounds[3]) / 2.0,
-    (bounds[4] + bounds[5]) / 2.0,
-  );
+  const center = vec3.create();
+  for (let i = 0; i < pointsToFit.length; ++i) {
+    vec3.add(center, center, pointsToFit[i]);
+  }
+  vec3.scale(center, center, 1.0 / pointsToFit.length);
 
   // Get the bounds in view coordinates
   const viewBounds = createBounds();
-  const visiblePoints = getCorners(bounds);
   const viewMat = mat4.create();
   toMat4(viewMat, pose);
-  for (let i = 0; i < visiblePoints.length; ++i) {
-    const point = visiblePoints[i];
-    vec3.transformMat4(point, point, viewMat);
-    addPoint(viewBounds, ...point);
+  const viewSpacePoint = vec3.create();
+  for (let i = 0; i < pointsToFit.length; ++i) {
+    const point = pointsToFit[i];
+    vec3.transformMat4(viewSpacePoint, point, viewMat);
+    addPoint(
+      viewBounds,
+      viewSpacePoint[0],
+      viewSpacePoint[1],
+      viewSpacePoint[2],
+    );
   }
 
   const xLength = getLength(viewBounds, 0);