fix(dicom): match the volume buffer to the decoded pixel values

src/core/streaming/__tests__/dicomChunkImage.spec.ts

@@ -0,0 +1,183 @@
+import { describe, it, expect, vi } from 'vitest';
+import { Chunk } from '@/src/core/streaming/chunk';
+import { ChunkState } from '@/src/core/streaming/chunkStateMachine';
+import { Tags } from '@/src/core/dicomTags';
+import DicomChunkImage, {
+  DicomChunkImageInit,
+} from '@/src/core/streaming/dicomChunkImage';
+import { ChunkStatus } from '@/src/core/streaming/chunkImage';
+
+const ROWS = 2;
+const COLUMNS = 2;
+const PIXELS_PER_SLICE = ROWS * COLUMNS;
+const PUBLIC_DSC_SLOPE = 112067.85375182;
+
+function metadataFor(z: number, overrides: Record<string, string> = {}) {
+  const metadata = [
+    [Tags.SOPInstanceUID, `1.2.3.${z}`],
+    [Tags.ImagePositionPatient, `0\\0\\${z}`],
+    [Tags.ImageOrientationPatient, '1\\0\\0\\0\\1\\0'],
+    [Tags.Rows, String(ROWS)],
+    [Tags.Columns, String(COLUMNS)],
+    [Tags.PixelSpacing, '1\\1'],
+    [Tags.BitsStored, '16'],
+    [Tags.PixelRepresentation, '0'],
+    [Tags.SamplesPerPixel, '1'],
+  ] as Array<[string, string]>;
+  Object.entries(overrides).forEach(([tag, value]) => {
+    const existing = metadata.find((entry) => entry[0] === tag);
+    if (existing) existing[1] = value;
+    else metadata.push([tag, value]);
+  });
+  return metadata;
+}
+
+// The slice's z position is also its pixel value, so the decoded contents of a
+// slice identify which chunk it came from.
+async function makeLoadedChunk(
+  z: number,
+  overrides: Record<string, string> = {}
+) {
+  const meta = metadataFor(z, overrides);
+  const chunk = new Chunk({
+    metaLoader: {
+      meta,
+      metaBlob: new Blob([`meta-${z}`]),
+      load: () => {},
+      stop: () => {},
+    },
+    dataLoader: {
+      data: new Blob([String(z)]),
+      load: () => {},
+      stop: () => {},
+    },
+  });
+  await chunk.loadMeta();
+  await chunk.loadData();
+  expect(chunk.state).toBe(ChunkState.Loaded);
+  return chunk;
+}
+
+function zOf(chunk: Chunk) {
+  const meta = Object.fromEntries(chunk.metadata!);
+  return Number(meta[Tags.ImagePositionPatient].split('\\')[2]);
+}
+
+function splitAndSortByPosition(chunks: Chunk[]) {
+  return Promise.resolve({
+    volume: [...chunks].sort((a, b) => zOf(a) - zOf(b)),
+  });
+}
+
+// Decodes a chunk to a constant frame, letting the test choose the array type.
+function decodeTo(dataFor: (value: number) => ArrayLike<number>) {
+  const read: DicomChunkImageInit['readDicomImage'] = async (file) => {
+    const value = Number(await file.text());
+    return {
+      image: {
+        size: [COLUMNS, ROWS, 1],
+        data: dataFor(value) as Uint16Array,
+        imageType: { components: 1 },
+      },
+    };
+  };
+  return read;
+}
+
+function sliceOf(image: DicomChunkImage, index: number) {
+  const scalars = image.getVtkImageData().getPointData().getScalars();
+  const data = scalars.getData();
+  return Array.from(
+    data.slice(index * PIXELS_PER_SLICE, (index + 1) * PIXELS_PER_SLICE)
+  );
+}
+
+async function loadRejectingSeries(
+  readDicomImage: DicomChunkImageInit['readDicomImage']
+) {
+  const image = new DicomChunkImage({
+    splitAndSort: splitAndSortByPosition,
+    readDicomImage,
+  });
+  const errors: unknown[] = [];
+  image.addEventListener('chunkError', ({ error }) => {
+    errors.push(error);
+  });
+  const [valid, invalid] = await Promise.all([
+    makeLoadedChunk(1, { [Tags.BitsStored]: '8' }),
+    makeLoadedChunk(2, { [Tags.BitsStored]: '8' }),
+  ]);
+
+  await image.addChunks([valid, invalid]);
+  await vi.waitFor(() =>
+    expect(image.getChunkStatuses()).toEqual([
+      ChunkStatus.Loaded,
+      ChunkStatus.Errored,
+    ])
+  );
+
+  expect(image.status.value).toBe('complete');
+  expect(errors).toHaveLength(1);
+  expect(sliceOf(image, 0)).toEqual(Array(PIXELS_PER_SLICE).fill(1));
+  expect(sliceOf(image, 1)).toEqual(Array(PIXELS_PER_SLICE).fill(0));
+  image.dispose();
+  return String(errors[0]);
+}
+
+describe('DicomChunkImage', () => {
+  it('preserves exact modality-rescaled pixels from the public DSC series', async () => {
+    // The public frames are 200x230; reduced geometry keeps the exact encoding,
+    // rescale, and an observed stored-pixel maximum in a focused volume test.
+    const decoded = Float64Array.from([
+      0,
+      PUBLIC_DSC_SLOPE,
+      2 * PUBLIC_DSC_SLOPE,
+      65131 * PUBLIC_DSC_SLOPE,
+    ]);
+    const image = new DicomChunkImage({
+      splitAndSort: splitAndSortByPosition,
+      readDicomImage: decodeTo(() => decoded),
+    });
+    const frame = await makeLoadedChunk(1, {
+      [Tags.SeriesInstanceUID]:
+        '1.3.6.1.4.1.9590.100.1.2.284777661700890778225181143863199482857',
+      [Tags.RescaleSlope]: String(PUBLIC_DSC_SLOPE),
+      [Tags.RescaleIntercept]: '0',
+    });
+
+    await image.addChunks([frame]);
+    await vi.waitFor(() =>
+      expect(image.getChunkStatuses()).toEqual([ChunkStatus.Loaded])
+    );
+
+    const data = image.getVtkImageData().getPointData().getScalars().getData();
+    expect(data).toBeInstanceOf(Float64Array);
+    expect(Array.from(data)).toEqual(Array.from(decoded));
+
+    image.dispose();
+  });
+
+  it('settles after rejecting decoded values its integer buffer cannot hold', async () => {
+    const message = await loadRejectingSeries(
+      decodeTo((value) =>
+        value === 2
+          ? new Uint16Array(PIXELS_PER_SLICE).fill(5000)
+          : new Uint8Array(PIXELS_PER_SLICE).fill(value)
+      )
+    );
+    expect(message).toContain('5000');
+    expect(message).toContain('Uint8Array');
+  });
+
+  it('settles after rejecting fractional samples bound for an integer buffer', async () => {
+    const message = await loadRejectingSeries(
+      decodeTo((value) =>
+        value === 2
+          ? new Float64Array(PIXELS_PER_SLICE).fill(2.5)
+          : new Uint8Array(PIXELS_PER_SLICE).fill(value)
+      )
+    );
+    expect(message).toContain('fractional');
+    expect(message).toContain('Uint8Array');
+  });
+});

src/core/streaming/dicomChunkImage.ts

@@ -2,12 +2,21 @@ import {
   buildSegmentGroups,
   ReadOverlappingSegmentationMeta,
   readVolumeSlice,
-  splitAndSort,
+  splitAndSort as splitAndSortChunks,
 } from '@/src/io/dicom';
 import { Chunk, waitForChunkState } from '@/src/core/streaming/chunk';
-import { Image, JsonCompatible, readImage } from '@itk-wasm/image-io';
+import {
+  Image,
+  JsonCompatible,
+  readImage as readItkImage,
+} from '@itk-wasm/image-io';
 import { getWorker } from '@/src/io/itk/worker';
-import { allocateImageFromChunks } from '@/src/utils/allocateImageFromChunks';
+import {
+  allocateImageFromChunks,
+  getBufferValueRange,
+  samplesAreIntegral,
+  valuesFitBuffer,
+} from '@/src/utils/allocateImageFromChunks';
 import { TypedArray } from '@kitware/vtk.js/types';
 import { Tags } from '@/src/core/dicomTags';
 import vtkDataArray from '@kitware/vtk.js/Common/Core/DataArray';
@@ -71,10 +80,28 @@ async function dicomSliceToImageUri(blob: Blob) {
   return itkImageToURI(itkImage);
 }
 
+function readDicomImage(file: File) {
+  return readItkImage(file, { webWorker: getWorker() });
+}
+
+export interface DicomChunkImageInit {
+  splitAndSort: (
+    chunks: Chunk[],
+    mapToBlob: (chunk: Chunk, index: number) => Blob
+  ) => Promise<Record<string, Chunk[]>>;
+  readDicomImage: (file: File) => Promise<{
+    image: Pick<Image, 'size' | 'data'> & {
+      imageType: Pick<Image['imageType'], 'components'>;
+    };
+  }>;
+}
+
 export default class DicomChunkImage
   extends BaseProgressiveImage
   implements ChunkImage
 {
+  private splitAndSort: DicomChunkImageInit['splitAndSort'];
+  private readDicomImage: DicomChunkImageInit['readDicomImage'];
   protected chunks: Chunk[];
   private chunkListeners: Array<() => void>;
   private thumbnailCache: WeakMap<Chunk, Promise<string>>;
@@ -85,9 +112,12 @@ export default class DicomChunkImage
     | (JsonCompatible & ReadOverlappingSegmentationMeta)
     | null;
 
-  constructor() {
+  constructor(init: Partial<DicomChunkImageInit> = {}) {
     super();
 
+    this.splitAndSort = init.splitAndSort ?? splitAndSortChunks;
+    this.readDicomImage = init.readDicomImage ?? readDicomImage;
+
     this.status.value = 'incomplete';
     this.loaded = computed(() => {
       return !this.loading.value && this.status.value === 'complete';
@@ -179,7 +209,7 @@ export default class DicomChunkImage
     });
 
     await Promise.all(chunks.map((chunk) => chunk.loadMeta()));
-    const chunksByVolume = await splitAndSort(
+    const chunksByVolume = await this.splitAndSort(
       this.chunks,
       (chunk) => chunk.metaBlob!
     );
@@ -385,11 +415,8 @@ export default class DicomChunkImage
       throw new Error(`Chunk ${chunkIndex} does not have data`);
 
     const chunkId = chunk.metadata ? getChunkId(chunk) : `index-${chunkIndex}`;
-    const result = await readImage(
-      new File([chunk.dataBlob], `file-${chunkIndex}.dcm`),
-      {
-        webWorker: getWorker(),
-      }
+    const result = await this.readDicomImage(
+      new File([chunk.dataBlob], `file-${chunkIndex}.dcm`)
     );
 
     if (!result.image.data)
@@ -405,13 +432,12 @@ export default class DicomChunkImage
 
     const scalars = this.vtkImageData.value.getPointData().getScalars();
     const pixelData = scalars.getData() as TypedArray;
+    const componentCount = scalars.getNumberOfComponents();
 
     const dims = this.vtkImageData.value.getDimensions();
-    const components = scalars.getNumberOfComponents();
 
-    // The volume buffer is sized from the first chunk's metadata, so each
-    // chunk gets a fixed slot: one frame per chunk in a multi-file volume,
-    // the whole volume when a single multi-frame chunk fills it.
+    // Each chunk gets a fixed slot: one frame per chunk in a multi-file
+    // volume, or the whole volume when a single multi-frame chunk fills it.
     const framesPerChunk = this.chunks.length > 1 ? 1 : dims[2];
     const [chunkWidth, chunkHeight] = result.image.size;
     const chunkFrames = result.image.size[2] ?? 1;
@@ -420,7 +446,7 @@ export default class DicomChunkImage
       chunkWidth !== dims[0] ||
       chunkHeight !== dims[1] ||
       chunkFrames !== framesPerChunk ||
-      chunkComponents !== components
+      chunkComponents !== componentCount
     ) {
       // A lone chunk defines the volume it fails to fit, so advice about
       // agreeing with the other files only makes sense for a multi-file volume.
@@ -431,34 +457,60 @@ export default class DicomChunkImage
       throw new Error(
         `File ${chunkId} (chunk ${chunkIndex}) does not fit the volume it belongs to. ` +
           `It decoded to ${chunkWidth}x${chunkHeight}x${chunkFrames} with ${chunkComponents} component(s), ` +
-          `but the volume has room for ${dims[0]}x${dims[1]}x${framesPerChunk} with ${components} component(s).` +
+          `but the volume has room for ${dims[0]}x${dims[1]}x${framesPerChunk} with ${componentCount} component(s).` +
           advice
       );
     }
 
-    const offset = dims[0] * dims[1] * components * chunkIndex;
-    pixelData.set(result.image.data as TypedArray, offset);
-
-    const rangeAlreadyInitialized = this.chunkStatus.some(
-      (status) => status === ChunkStatus.Loaded
-    );
-
-    // update the data range
     const chunkDataRange: Array<[number, number]> = [];
-    for (let comp = 0; comp < scalars.getNumberOfComponents(); comp++) {
+    for (let comp = 0; comp < componentCount; comp++) {
       const { min, max } = fastComputeRange(
         result.image.data as unknown as number[],
         comp,
-        scalars.getNumberOfComponents()
+        componentCount
       );
       chunkDataRange.push([min, max]);
+    }
 
-      const curRange = scalars.getRange(comp);
+    // The buffer is allocated for the range every chunk's tags declare, so a
+    // chunk only fails here when its decoded values disagree with its tags.
+    // TypedArray.set raises nothing for such values: integers wrap and
+    // fractions truncate.
+    const chunkMin = Math.min(...chunkDataRange.map(([min]) => min));
+    const chunkMax = Math.max(...chunkDataRange.map(([, max]) => max));
+    const decoded = result.image.data as unknown as ArrayLike<number>;
+    if (!valuesFitBuffer({ min: chunkMin, max: chunkMax }, pixelData)) {
+      const bufferRange = getBufferValueRange(pixelData)!;
+      throw new Error(
+        `File ${chunkId} (chunk ${chunkIndex}) has pixel values the volume it belongs to cannot represent. ` +
+          `Its pixel values run from ${chunkMin} to ${chunkMax}, but the volume's buffer is ` +
+          `${pixelData.constructor.name}, holding values from ${bufferRange.min} to ${bufferRange.max}. ` +
+          `Every file in a volume must decode to values its buffer can hold without conversion.`
+      );
+    }
+    if (!samplesAreIntegral(decoded, pixelData)) {
+      throw new Error(
+        `File ${chunkId} (chunk ${chunkIndex}) has fractional pixel values the volume it belongs to cannot represent. ` +
+          `Its pixel values run from ${chunkMin} to ${chunkMax}, but the volume's buffer is ` +
+          `${pixelData.constructor.name}, which holds only whole numbers. ` +
+          `Every file in a volume must decode to values its buffer can hold without conversion.`
+      );
+    }
+
+    const offset = dims[0] * dims[1] * componentCount * chunkIndex;
+    pixelData.set(result.image.data as TypedArray, offset);
 
+    const rangeAlreadyInitialized = this.chunkStatus.some(
+      (status) => status === ChunkStatus.Loaded
+    );
+
+    // update the data range
+    chunkDataRange.forEach(([min, max], comp) => {
+      const curRange = scalars.getRange(comp);
       const newMin = rangeAlreadyInitialized ? Math.min(min, curRange[0]) : min;
       const newMax = rangeAlreadyInitialized ? Math.max(max, curRange[1]) : max;
       scalars.setRange({ min: newMin, max: newMax }, comp);
-    }
+    });
     scalars.modified(); // so image-stats will trigger update of range
 
     chunk.setUserData(DATA_RANGE_KEY, chunkDataRange);