feat(ngff zarr): better multichannel support

src/Context/ImageActorContext.js

@@ -5,10 +5,10 @@ class ImageActorContext {
   // The rendered image / label image scale
   renderedScale = null
 
-  // MultiscaleChunked label image to be visualized
+  // MultiscaleSpatialImage label image to be visualized
   labelImage = null
 
-  // MultiscaleChunked label image to be visualized for use with
+  // MultiscaleSpatialImage label image to be visualized for use with
   // interactive, manual editing as opposed to stored or algorithmic results
   editorLabelImage = null
 

src/IO/ImageDataFromChunks.worker.js

@@ -1,9 +1,19 @@
 import registerWebworker from 'webworker-promise/lib/register'
 import componentTypeToTypedArray from './componentTypeToTypedArray'
-import { toDimensionArray } from './dimensionUtils'
 
 const haveSharedArrayBuffer = typeof self.SharedArrayBuffer === 'function'
 
+const validateIndices = ({ chunkStart, chunkEnd, roiStart, roiEnd }) => {
+  if (
+    ['x', 'y', 'z'].some(
+      dim => chunkStart[dim] > roiEnd[dim] || chunkEnd[dim] < roiStart[dim]
+    )
+  ) {
+    // We should never get here...
+    console.error('Requested a chunk outside the region of interest!')
+  }
+}
+
 registerWebworker().operation(
   'imageDataFromChunks',
   ({
@@ -14,29 +24,15 @@ registerWebworker().operation(
     indexStart,
     indexEnd,
   }) => {
-    const chunkSize = toDimensionArray(['c', 'x', 'y', 'z'], info.chunkSize)
-    const chunkStrides = [
-      chunkSize[0],
-      chunkSize[0] * chunkSize[1],
-      chunkSize[0] * chunkSize[1] * chunkSize[2],
-      chunkSize[0] * chunkSize[1] * chunkSize[2] * chunkSize[3],
-    ] // c, x, y, z,
-
-    const size = toDimensionArray(['x', 'y', 'z'], info.arrayShape)
-    const components = imageType.components
-
-    const pixelStrides = [
-      components,
-      components * size[0],
-      components * size[0] * size[1],
-      components * size[0] * size[1] * size[2],
-    ] // c, x, y, z
+    info.arrayShape.set('c', imageType.components)
 
     const pixelArrayType = componentTypeToTypedArray.get(
       imageType.componentType
     )
     let pixelArray = null
-    const pixelArrayElements = size.reduce((a, b) => a * b) * components
+    const pixelArrayElements = Array.from(info.arrayShape.values()).reduce(
+      (a, b) => a * b
+    )
     if (haveSharedArrayBuffer) {
       const pixelArrayBytes =
         pixelArrayElements * pixelArrayType.BYTES_PER_ELEMENT
@@ -46,68 +42,75 @@ registerWebworker().operation(
       pixelArray = new pixelArrayType(pixelArrayElements)
     }
 
+    const arrayShape = Object.fromEntries(info.arrayShape)
+    const pixelStrides = {
+      z: arrayShape.c * arrayShape.x * arrayShape.y,
+      y: arrayShape.c * arrayShape.x,
+      x: arrayShape.c,
+    }
+
+    const chunkSize = Object.fromEntries(info.chunkSize)
+    const chunkStrides = {
+      c: chunkSize.x * chunkSize.y * chunkSize.z,
+      z: chunkSize.x * chunkSize.y,
+      y: chunkSize.x,
+      x: 1,
+    }
+
     for (let index = 0; index < chunkIndices.length; index++) {
       const chunk = chunks[index]
-      const [h, i, j, k, l] = chunkIndices[index]
+      const [c, x, y, z /*t*/] = chunkIndices[index]
 
-      const chunkStart = [
-        i * chunkSize[1],
-        j * chunkSize[2],
-        k * chunkSize[3],
-        l * chunkSize[4],
-      ]
-      const chunkEnd = [
-        (i + 1) * chunkSize[1],
-        (j + 1) * chunkSize[2],
-        (k + 1) * chunkSize[3],
-        (l + 1) * chunkSize[4],
-      ]
-      // Skip if the chunk lives outside the region of interest
-      if (
-        chunkStart[0] > indexEnd[0] ||
-        chunkEnd[0] < indexStart[0] ||
-        chunkStart[1] > indexEnd[1] ||
-        chunkEnd[1] < indexStart[1] ||
-        chunkStart[2] > indexEnd[2] ||
-        chunkEnd[2] < indexStart[2] ||
-        chunkStart[3] > indexEnd[3] ||
-        chunkEnd[3] < indexStart[3]
-      ) {
-        // We should never get here...
-        console.error('Requested a chunk outside the region of interest!')
+      const chunkStart = {
+        c: c * chunkSize.c,
+        z: z * chunkSize.z,
+        y: y * chunkSize.y,
+        x: x * chunkSize.x,
       }
-      const itStart = [
-        Math.max(chunkStart[0], indexStart[0]),
-        Math.max(chunkStart[1], indexStart[1]) - i,
-        Math.max(chunkStart[2], indexStart[2]),
-        Math.max(chunkStart[3], indexStart[3]),
-      ]
-      const itEnd = [
-        Math.min(chunkEnd[0], indexEnd[0]),
-        Math.min(chunkEnd[1], indexEnd[1]),
-        Math.min(chunkEnd[2], indexEnd[2]),
-        Math.min(chunkEnd[3], indexEnd[3]),
-      ]
-      const itChunkOffsets = [0, 0, 0, 0]
-      itChunkOffsets[3] = chunkStrides[3] * l
-      const itPixelOffsets = [0, 0, 0]
-      for (let kk = itStart[2]; kk < itEnd[2]; kk++) {
-        itChunkOffsets[2] = chunkStrides[2] * (kk - k * chunkSize[3])
-        itPixelOffsets[2] = pixelStrides[2] * (kk - indexStart[2])
-        for (let jj = itStart[1]; jj < itEnd[1]; jj++) {
-          itChunkOffsets[1] = chunkStrides[1] * (jj - j * chunkSize[2])
-          itPixelOffsets[1] = pixelStrides[1] * (i + jj - indexStart[1])
-          for (let ii = itStart[0]; ii < itEnd[0]; ii++) {
-            const begin =
-              ii + itChunkOffsets[1] + itChunkOffsets[2] + itChunkOffsets[3]
-            const offset =
-              h + ii * pixelStrides[0] + itPixelOffsets[1] + itPixelOffsets[2]
+      const chunkEnd = {
+        c: (c + 1) * chunkSize.c,
+        z: (z + 1) * chunkSize.z,
+        y: (y + 1) * chunkSize.y,
+        x: (x + 1) * chunkSize.x,
+      }
+      const roiStart = Object.fromEntries(indexStart)
+      const roiEnd = Object.fromEntries(indexEnd)
+      validateIndices({ chunkStart, chunkEnd, roiStart, roiEnd })
 
-            pixelArray[offset] = chunk[begin]
-          } // for every column
-        } // for every row
-      } // for every slice
-    }
+      const itStart = {
+        c: Math.max(chunkStart.c, roiStart.c),
+        z: Math.max(chunkStart.z, roiStart.z),
+        y: Math.max(chunkStart.y, roiStart.y) - x,
+        x: Math.max(chunkStart.x, roiStart.x),
+      }
+      const itEnd = {
+        c: Math.min(chunkEnd.c, roiEnd.c),
+        z: Math.min(chunkEnd.z, roiEnd.z),
+        y: Math.min(chunkEnd.y, roiEnd.y),
+        x: Math.min(chunkEnd.x, roiEnd.x),
+      }
+
+      for (let cc = itStart.c; cc < itEnd.c; cc++) {
+        for (let zz = itStart.z; zz < itEnd.z; zz++) {
+          for (let yy = itStart.y; yy < itEnd.y; yy++) {
+            for (let xx = itStart.x; xx < itEnd.x; xx++) {
+              pixelArray[
+                cc +
+                  zz * pixelStrides.z +
+                  (yy + x) * pixelStrides.y +
+                  xx * pixelStrides.x
+              ] =
+                chunk[
+                  xx +
+                    (yy - y * chunkSize.y) * chunkStrides.y +
+                    (zz - z * chunkSize.z) * chunkStrides.z +
+                    (cc - c * chunkSize.c) * chunkStrides.c
+                ]
+            } // for every column
+          } // for every row
+        } // for every slice
+      } // for every channel
+    } // for every chunk
 
     let response = pixelArray
     if (!haveSharedArrayBuffer) {

src/IO/MultiscaleSpatialImage.js

@@ -133,14 +133,12 @@ class MultiscaleSpatialImage {
 
     const info = this.scaleInfo[scale]
 
-    const size = toDimensionArray(['x', 'y', 'z'], info.arrayShape)
-    const start = [0, 0, 0, 0] // x, y, z, t
-    const end = [
-      start[0] + size[0],
-      start[1] + size[1],
-      start[2] + size[2],
-      start[3] + 1,
-    ] // x, y, z, t
+    const start = new Map(Object.entries({ t: 0, c: 0, z: 0, y: 0, x: 0 }))
+    const end = Array.from(start).reduce(
+      (end, [dim, startIndex]) =>
+        end.set(dim, startIndex + info.arrayShape.get(dim)),
+      new Map()
+    )
 
     const numChunks = toDimensionArray(CXYZT, info.chunkCount)
     const l = 0
@@ -198,7 +196,9 @@ class MultiscaleSpatialImage {
       origin,
       spacing,
       direction: this.direction,
-      size: size.slice(0, this.imageType.dimension),
+      size: ['x', 'y', 'z']
+        .slice(0, this.imageType.dimension)
+        .map(dim => info.arrayShape.get(dim)),
       data: pixelArray,
     }
 

src/IO/ZarrMultiscaleSpatialImage.js

@@ -4,7 +4,7 @@ import MultiscaleSpatialImage from './MultiscaleSpatialImage'
 import bloscZarrDecompress from '../Compression/bloscZarrDecompress'
 import ZarrStore from './ZarrStore'
 import HttpStore from './HttpStore'
-import { CXYZT, toDimensionArray, toDimensionMap } from './dimensionUtils'
+import { CXYZT, toDimensionMap } from './dimensionUtils'
 
 // ends with zarr and optional nested image name like foo.zarr/image1
 export const isZarr = url => /zarr((\/)[\w-]+\/?)?$/.test(url)
@@ -26,7 +26,8 @@ const dtypeToComponentType = new Map([
   ['<f8', FloatTypes.Float64],
 ])
 
-const CONTIGUOUS_CHANNEL_INDEXING = Object.freeze(['t', 'c', 'z', 'y', 'x'])
+const TCZYX = Object.freeze(['t', 'c', 'z', 'y', 'x'])
+const toDimensionMapTCZYX = (dims, array) => toDimensionMap(dims, array, TCZYX)
 
 const composeTransforms = (transforms = [], dimCount) =>
   transforms.reduce(
@@ -105,12 +106,23 @@ const computeScaleSpacing = ({
   pixelArrayMetadata,
   dataset,
 }) => {
-  // "axis" metadata not defined in ngff V0.1 so fallback to CONTIGUOUS_CHANNEL_INDEXING
-  const dims =
-    multiscaleImage.axes?.map(({ name }) => name) ?? CONTIGUOUS_CHANNEL_INDEXING
+  // "axis" metadata not defined in ngff V0.1 so fallback to TCZYX
+  const dims = multiscaleImage.axes?.map(({ name }) => name) ?? TCZYX
 
   const { shape, chunks } = pixelArrayMetadata
 
+  const chunkSize = toDimensionMapTCZYX(dims, chunks)
+  const arrayShape = toDimensionMapTCZYX(dims, shape)
+
+  const componentsInData = arrayShape.get('c')
+  const components = Math.min(componentsInData, 3)
+  if (componentsInData !== components) {
+    console.warn(
+      `itk-vtk-viewer: ${componentsInData} components are not supported.`
+    )
+    arrayShape.set('c', components)
+  }
+
   return {
     dims,
     pixelArrayMetadata,
@@ -119,12 +131,12 @@ const computeScaleSpacing = ({
     coords: makeCoords(dims, shape, multiscaleImage, dataset),
     ranges: zattrs.ranges ?? undefined,
     direction: zattrs.direction ?? undefined,
-    chunkCount: toDimensionMap(
+    chunkCount: toDimensionMapTCZYX(
       dims,
-      dims.map((_, i) => Math.ceil(shape[i] / chunks[i]))
+      dims.map(dim => Math.ceil(arrayShape.get(dim) / chunkSize.get(dim)))
     ),
-    chunkSize: toDimensionMap(dims, chunks),
-    arrayShape: toDimensionMap(dims, shape),
+    chunkSize,
+    arrayShape,
   }
 }
 
@@ -150,19 +162,12 @@ const extractScaleSpacing = async store => {
 
   const info = scaleInfo[0]
 
-  const componentsInData = info.arrayShape.get('c') ?? 1
-  const components = Math.min(componentsInData, 3)
-  if (componentsInData !== components) {
-    console.warn(
-      `itk-vtk-viewer: ${componentsInData} components are not supported.  Maximum 3 components are supported.`
-    )
-  }
+  const components = info.arrayShape.get('c')
 
   const imageType = {
     // How many spatial dimensions?  Count greater than 1, X Y Z elements because "axis" metadata not defined in ngff V0.1
-    dimension: toDimensionArray(['x', 'y', 'z'], info.arrayShape).filter(
-      size => size > 1
-    ).length,
+    dimension: ['x', 'y', 'z'].filter(dim => info.arrayShape.get(dim) > 1)
+      .length,
     pixelType:
       components === 1 ? PixelTypes.Scalar : PixelTypes.VariableLengthVector,
     componentType: dtypeToComponentType.get(info.pixelArrayMetadata.dtype),

src/IO/dimensionUtils.js

@@ -1,7 +1,16 @@
 export const CXYZT = Object.freeze(['c', 'x', 'y', 'z', 't']) // viewer indexing
 
-export const toDimensionMap = (dims, array) =>
+const default1 = (dimMap, ensuredDims) =>
+  ensuredDims.reduce(
+    (map, dim) => (map.has(dim) ? map : map.set(dim, 1)),
+    dimMap
+  )
+
+const makeDimMap = (dims, array) =>
   new Map(dims.map((dim, i) => [dim, array[i]]))
 
+export const toDimensionMap = (dims, array, minimumDimensions = []) =>
+  default1(makeDimMap(dims, array), minimumDimensions)
+
 export const toDimensionArray = (dims, map) =>
   dims.map(dim => map.get(dim) ?? 1)