feat(loaders): Support mutifile & multiscale OME-TIFF (#748)

.changeset/popular-steaks-end.md

@@ -0,0 +1,13 @@
+---
+'@vivjs/loaders': minor
+---
+
+feat: Support multiscale multifile-OME-TIFFs
+
+This release extends Viv's multifile OME-TIFF data-loading capabilities to multiscale TIFFs as well. The `loadOmeTiff` utility now recognizes and loads multiresolution images described in a `companion.ome` metadata file.
+
+```js
+import { loadOmeTiff } from '@vivjs/loaders';
+
+let loader = await loadOmeTiff("http://localhost:8080/data.companion.ome");
+```

packages/loaders/src/omexml.ts

@@ -165,18 +165,18 @@ const ImageSchema = z
   })
   .transform(flattenAttributes);
 
-const OmeSchema = z
-  .object({
-    Image: z.preprocess(ensureArray, ImageSchema.array())
-  })
-  .extend({
-    attr: z.object({
-      xmlns: z.string(),
-      'xmlns:xsi': z.string(),
-      'xsi:schemaLocation': z.string()
-    })
-  })
-  .transform(flattenAttributes);
+const OmeSchema = z.object({
+  Image: z.preprocess(ensureArray, ImageSchema.array())
+});
+// TODO: Verify that these attributes are always present
+// .extend({
+//   attr: z.object({
+//     'xmlns': z.string(),
+//     'xmlns:xsi': z.string(),
+//     'xsi:schemaLocation': z.string()
+//   })
+// })
+// .transform(flattenAttributes);
 
 export function fromString(str: string) {
   const raw = parseXML(str);

packages/loaders/src/tiff/lib/indexers.ts

@@ -1,97 +1,69 @@
 /* eslint-disable no-use-before-define */
 import { GeoTIFFImage, GeoTIFF } from 'geotiff';
 import type { OmeTiffSelection } from './utils';
-import type { OmeXml } from '../../omexml';
 import type { MultiTiffImage } from '../multi-tiff';
 
+type ImageFileDirectory = Awaited<ReturnType<GeoTIFF['parseFileDirectoryAt']>>;
+
 export type OmeTiffIndexer = (
   sel: OmeTiffSelection,
-  z: number
+  resolutionLevel: number
 ) => Promise<GeoTIFFImage>;
 
+/**
+ * A function that resolves a given selection to a GeoTIFF and IFD index.
+ *
+ * Ideally we could just return the GeoTIFFImage, but for the legacy
+ * bioformats case we need to return the GeoTIFF object and the IFD index.
+ */
+export type OmeTiffResolver = {
+  (sel: OmeTiffSelection):
+    | { tiff: GeoTIFF; ifdIndex: number }
+    | Promise<{ tiff: GeoTIFF; ifdIndex: number }>;
+};
+
 /*
  * An "indexer" for a GeoTIFF-based source is a function that takes a
  * "selection" (e.g. { z, t, c }) and returns a Promise for the GeoTIFFImage
  * object corresponding to that selection.
  *
  * For OME-TIFF images, the "selection" object is the same regardless of
  * the format version. However, modern version of Bioformats have a different
- * memory layout for pyramidal resolutions. Thus, we have two different "indexers"
+ * memory layout for pyramidal resolutions. Thus, we have different "indexers"
  * depending on which format version is detected.
- *
- * TODO: We currently only support indexing the first image in the OME-TIFF with
- * our indexers. There can be multiple images in an OME-TIFF, so supporting these
- * images will require extending these indexers or creating new methods.
- */
-
-/*
- * Returns an indexer for legacy Bioformats images. This assumes that
- * downsampled resolutions are stored sequentially in the OME-TIFF.
- */
-export function getOmeLegacyIndexer(
-  tiff: GeoTIFF,
-  rootMeta: OmeXml
-): OmeTiffIndexer {
-  const { SizeT, SizeC, SizeZ } = rootMeta[0].Pixels;
-  const ifdIndexer = getOmeIFDIndexer(rootMeta, 0);
-
-  return (sel: OmeTiffSelection, pyramidLevel: number) => {
-    // Get IFD index at base pyramid level
-    const index = ifdIndexer(sel);
-    // Get index of first image at pyramidal level
-    const pyramidIndex = pyramidLevel * SizeZ * SizeT * SizeC;
-    // Return image at IFD index for pyramidal level
-    return tiff.getImage(index + pyramidIndex);
-  };
-}
-
-/*
- * Returns an indexer for modern Bioforamts images that store multiscale
- * resolutions using SubIFDs.
- *
- * The ifdIndexer returns the 'index' to the base resolution for a
- * particular 'selection'. The SubIFDs to the downsampled resolutions
- * of the 'selection' are stored within the `baseImage.fileDirectory`.
- * We use the SubIFDs to get the IFD for the corresponding sub-resolution.
- *
- * NOTE: This function create a custom IFD cache rather than mutating
- * `GeoTIFF.ifdRequests` with a random offset. The IFDs are cached in
- * an ES6 Map that maps a string key that identifies the selection uniquely
- * to the corresponding IFD.
  */
-export function getOmeSubIFDIndexer(
-  tiff: GeoTIFF,
-  rootMeta: OmeXml,
-  image = 0
-): OmeTiffIndexer {
-  const ifdIndexer = getOmeIFDIndexer(rootMeta, image);
-  const ifdCache: Map<
-    string,
-    ReturnType<GeoTIFF['parseFileDirectoryAt']>
-  > = new Map();
-
+export function createOmeImageIndexerFromResolver(
+  resolveBaseResolutionImageLocation: OmeTiffResolver,
+  image: {
+    size: { z: number; t: number; c: number };
+  }
+) {
+  const ifdCache: ImageFileDirectory[] = [];
   return async (sel: OmeTiffSelection, pyramidLevel: number) => {
-    const index = ifdIndexer(sel);
-    const baseImage = await tiff.getImage(index);
+    const { tiff, ifdIndex } = await resolveBaseResolutionImageLocation(sel);
+    const baseImage = await tiff.getImage(ifdIndex);
 
     // It's the highest resolution, no need to look up SubIFDs.
     if (pyramidLevel === 0) {
       return baseImage;
     }
 
-    const { SubIFDs } = baseImage.fileDirectory;
-    if (!SubIFDs) {
-      throw Error('Indexing Error: OME-TIFF is missing SubIFDs.');
+    let index;
+    if (baseImage.fileDirectory.SubIFDs) {
+      index = baseImage.fileDirectory.SubIFDs[pyramidLevel - 1];
+    } else {
+      // Legacy Bioformats OME-TIFFs don't have SubIFDs, and instead
+      // store each resolution level in a separate IFD at the end
+      // of the base image.
+      const resolutionOffset =
+        pyramidLevel * image.size.z * image.size.t * image.size.c;
+      index = ifdIndex + resolutionOffset;
     }
 
-    // Get IFD for the selection at the pyramidal level
-    const key = `${sel.t}-${sel.c}-${sel.z}-${pyramidLevel}`;
-    if (!ifdCache.has(key)) {
-      // Only create a new request if we don't have the key.
-      const subIfdOffset = SubIFDs[pyramidLevel - 1];
-      ifdCache.set(key, tiff.parseFileDirectoryAt(subIfdOffset));
+    if (!ifdCache[index]) {
+      ifdCache[index] = await tiff.parseFileDirectoryAt(index);
     }
-    const ifd = await ifdCache.get(key)!;
+    const ifd = ifdCache[index];
 
     // Create a new image object manually from IFD
     return new GeoTIFFImage(
@@ -105,53 +77,6 @@ export function getOmeSubIFDIndexer(
   };
 }
 
-/*
- * Returns a function that computes the image index based on the dimension
- * order and dimension sizes.
- */
-function getOmeIFDIndexer(
-  rootMeta: OmeXml,
-  image = 0
-): (sel: OmeTiffSelection) => number {
-  const { SizeC, SizeZ, SizeT, DimensionOrder } = rootMeta[image].Pixels;
-  // For multi-image OME-TIFF files, we need to offset by the full dimensions
-  // of the previous images dimensions i.e Z * C * T of image - 1 + that of image - 2 etc.
-  let imageOffset = 0;
-  if (image > 0) {
-    for (let i = 0; i < image; i += 1) {
-      const {
-        SizeC: prevSizeC,
-        SizeZ: prevSizeZ,
-        SizeT: prevSizeT
-      } = rootMeta[i].Pixels;
-      imageOffset += prevSizeC * prevSizeZ * prevSizeT;
-    }
-  }
-  switch (DimensionOrder) {
-    case 'XYZCT': {
-      return ({ t, c, z }) => imageOffset + t * SizeZ * SizeC + c * SizeZ + z;
-    }
-    case 'XYZTC': {
-      return ({ t, c, z }) => imageOffset + c * SizeZ * SizeT + t * SizeZ + z;
-    }
-    case 'XYCTZ': {
-      return ({ t, c, z }) => imageOffset + z * SizeC * SizeT + t * SizeC + c;
-    }
-    case 'XYCZT': {
-      return ({ t, c, z }) => imageOffset + t * SizeC * SizeZ + z * SizeC + c;
-    }
-    case 'XYTCZ': {
-      return ({ t, c, z }) => imageOffset + z * SizeT * SizeC + c * SizeT + t;
-    }
-    case 'XYTZC': {
-      return ({ t, c, z }) => imageOffset + c * SizeT * SizeZ + z * SizeT + t;
-    }
-    default: {
-      throw new Error(`Invalid OME-XML DimensionOrder, got ${DimensionOrder}.`);
-    }
-  }
-}
-
 export function getMultiTiffIndexer(tiffs: MultiTiffImage[]) {
   function selectionToKey({ c = 0, t = 0, z = 0 }: OmeTiffSelection): string {
     return `${c}-${t}-${z}`;

packages/loaders/src/tiff/lib/utils.ts

@@ -1,14 +1,9 @@
-import {
-  type GeoTIFF,
-  type GeoTIFFImage,
-  fromFile,
-  fromUrl,
-  fromBlob
-} from 'geotiff';
-import { getLabels, DTYPE_LOOKUP } from '../../utils';
+import { fromFile, fromUrl, fromBlob } from 'geotiff';
+import { getLabels, DTYPE_LOOKUP, prevPowerOf2, assert } from '../../utils';
+import { createOffsetsProxy } from './proxies';
+import type { GeoTIFF, GeoTIFFImage } from 'geotiff';
 import type { OmeXml, PhysicalUnit, DimensionOrder } from '../../omexml';
 import type { MultiTiffImage } from '../multi-tiff';
-import { createOffsetsProxy } from './proxies';
 
 // TODO: Remove the fancy label stuff
 export type OmeTiffDims =
@@ -36,7 +31,7 @@ type PhysicalSizes = {
   z?: PhysicalSize;
 };
 
-function extractPhysicalSizesfromOmeXml(
+export function extractPhysicalSizesfromPixels(
   d: OmeXml[number]['Pixels']
 ): undefined | PhysicalSizes {
   if (
@@ -60,51 +55,67 @@ function extractPhysicalSizesfromOmeXml(
   return physicalSizes;
 }
 
-export function getOmePixelSourceMeta({ Pixels }: OmeXml[0]) {
+export function parsePixelDataType(dtype: string) {
+  assert(dtype in DTYPE_LOOKUP, `Pixel type ${dtype} not supported.`);
+  return DTYPE_LOOKUP[dtype as keyof typeof DTYPE_LOOKUP];
+}
+
+export function extractAxesFromPixels(d: OmeXml[number]['Pixels']) {
   // e.g. 'XYZCT' -> ['t', 'c', 'z', 'y', 'x']
-  const labels = getLabels(Pixels.DimensionOrder);
+  const labels = getLabels(d['DimensionOrder']);
 
   // Compute "shape" of image
   const shape: number[] = Array(labels.length).fill(0);
-  shape[labels.indexOf('t')] = Pixels.SizeT;
-  shape[labels.indexOf('c')] = Pixels.SizeC;
-  shape[labels.indexOf('z')] = Pixels.SizeZ;
+  shape[labels.indexOf('t')] = d['SizeT'];
+  shape[labels.indexOf('c')] = d['SizeC'];
+  shape[labels.indexOf('z')] = d['SizeZ'];
+  shape[labels.indexOf('y')] = d['SizeY'];
+  shape[labels.indexOf('x')] = d['SizeX'];
 
   // Push extra dimension if data are interleaved.
-  if (Pixels.Interleaved) {
+  if (d['Interleaved']) {
     // @ts-expect-error private, unused dim name for selection
     labels.push('_c');
     shape.push(3);
   }
 
-  // Creates a new shape for different level of pyramid.
-  // Assumes factor-of-two downsampling.
-  const getShape = (level: number = 0) => {
-    const s = [...shape];
-    s[labels.indexOf('x')] = Pixels.SizeX >> level;
-    s[labels.indexOf('y')] = Pixels.SizeY >> level;
-    return s;
-  };
+  return { labels, shape };
+}
 
-  if (!(Pixels.Type in DTYPE_LOOKUP)) {
-    throw Error(`Pixel type ${Pixels.Type} not supported.`);
-  }
+/**
+ * Compute the shape of the image at a given resolution level.
+ *
+ * Assumes that the image is downsampled by a factor of 2 for each
+ * pyramid level.
+ */
+export function getShapeForBinaryDownsampleLevel(options: {
+  axes: { shape: number[]; labels: string[] };
+  level: number;
+}) {
+  const { axes, level } = options;
+  const xIndex = axes.labels.indexOf('x');
+  assert(xIndex !== -1, 'x dimension not found');
+  const yIndex = axes.labels.indexOf('y');
+  assert(yIndex !== -1, 'y dimension not found');
+  const resolutionShape = axes.shape.slice();
+  resolutionShape[xIndex] = axes.shape[xIndex] >> level;
+  resolutionShape[yIndex] = axes.shape[yIndex] >> level;
+  return resolutionShape;
+}
 
-  const dtype = DTYPE_LOOKUP[Pixels.Type as keyof typeof DTYPE_LOOKUP];
-  const maybePhysicalSizes = extractPhysicalSizesfromOmeXml(Pixels);
-  if (maybePhysicalSizes) {
-    return { labels, getShape, dtype, physicalSizes: maybePhysicalSizes };
-  }
-  return { labels, getShape, dtype };
+export function getTiffTileSize(image: GeoTIFFImage) {
+  const tileWidth = image.getTileWidth();
+  const tileHeight = image.getTileHeight();
+  const size = Math.min(tileWidth, tileHeight);
+  // deck.gl requirement for power-of-two tile size.
+  return prevPowerOf2(size);
 }
 
 // Inspired by/borrowed from https://geotiffjs.github.io/geotiff.js/geotiffimage.js.html#line297
 function guessImageDataType(image: GeoTIFFImage) {
   // Assuming these are flat TIFFs, just grab the info for the first image/sample.
   const sampleIndex = 0;
-  const format = image.fileDirectory.SampleFormat
-    ? image.fileDirectory.SampleFormat[sampleIndex]
-    : 1;
+  const format = image.fileDirectory?.SampleFormat?.[sampleIndex] ?? 1;
   const bitsPerSample = image.fileDirectory.BitsPerSample[sampleIndex];
   switch (format) {
     case 1: // unsigned integer data

packages/loaders/src/tiff/multi-tiff.ts

@@ -1,10 +1,10 @@
 import type { GeoTIFFImage } from 'geotiff';
 
 import TiffPixelSource from './pixel-source';
-import { guessTiffTileSize } from '../utils';
 import {
   getMultiTiffMetadata,
   getMultiTiffMeta,
+  getTiffTileSize,
   type OmeTiffSelection
 } from './lib/utils';
 import type Pool from './lib/Pool';
@@ -56,7 +56,7 @@ export async function load(
     firstImage.fileDirectory;
   // Not sure if we need this or if the order matters for this use case.
   const dimensionOrder = 'XYZCT';
-  const tileSize = guessTiffTileSize(firstImage);
+  const tileSize = getTiffTileSize(firstImage);
   const meta = { photometricInterpretation };
   const indexer = getMultiTiffIndexer(images);
   const { shape, labels, dtype } = getMultiTiffMeta(dimensionOrder, images);