From 9dfd670a5b3d0ed1d55e52a8bd8bd81f3b2bd7cb Mon Sep 17 00:00:00 2001 From: Paul Elliott Date: Tue, 1 Sep 2026 16:56:41 -0400 Subject: [PATCH] fix(dicom): match the volume buffer to the decoded pixel values The volume buffer's element type was derived from the first instance's stored range put through its rescale, floored and ceiled to integers. ITK applies the rescale while decoding and picks the output type from it: an integral slope and intercept keep an integer type, anything fractional produces float64. A fractional rescale whose endpoints landed on integers got an integer buffer, and every fraction was dropped on write. An integral range past 32 bits threw at allocation. A later instance decoding outside the first instance's type wrote through TypedArray.set with no error, wrapping every value the type cannot hold. Allocate for the union of every instance's declared range, choosing Float64Array whenever the slope or intercept is fractional or the range is not a safe 32-bit integer range, and the narrowest integer type otherwise. Ordinary CT and MR allocation is unchanged. An instance whose tags give no finite range is left out of the union and judged on its decoded values. Check each chunk after it decodes and before it is written. Its actual range must fit the buffer's element type, and float output bound for an integer buffer must be whole numbers. Both rejections name the file, the values it decoded to, and the buffer it has to fit. Only decoder output that disagrees with the tags can fail these checks. Reproduced with a public DSC perfusion series whose slope is 112067.85375182 over 16-bit unsigned storage. --- .../__tests__/dicomChunkImage.spec.ts | 183 ++++++++++ src/core/streaming/dicomChunkImage.ts | 108 ++++-- .../__tests__/allocateImageFromChunks.spec.ts | 316 ++++++++++++++++++ src/utils/allocateImageFromChunks.ts | 168 ++++++++-- tests/specs/dicom-modality-rescale.e2e.ts | 71 ++++ tests/specs/imageCacheUtils.ts | 45 ++- tests/specs/syntheticDicom.ts | 56 +++- 7 files changed, 888 insertions(+), 59 deletions(-) create mode 100644 src/core/streaming/__tests__/dicomChunkImage.spec.ts create mode 100644 tests/specs/dicom-modality-rescale.e2e.ts diff --git a/src/core/streaming/__tests__/dicomChunkImage.spec.ts b/src/core/streaming/__tests__/dicomChunkImage.spec.ts new file mode 100644 index 000000000..9ddc38851 --- /dev/null +++ b/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 = {}) { + 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 = {} +) { + 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) { + 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'); + }); +}); diff --git a/src/core/streaming/dicomChunkImage.ts b/src/core/streaming/dicomChunkImage.ts index 2d2df56f7..3e71a2167 100644 --- a/src/core/streaming/dicomChunkImage.ts +++ b/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>; + readDicomImage: (file: File) => Promise<{ + image: Pick & { + imageType: Pick; + }; + }>; +} + 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>; @@ -85,9 +112,12 @@ export default class DicomChunkImage | (JsonCompatible & ReadOverlappingSegmentationMeta) | null; - constructor() { + constructor(init: Partial = {}) { 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; + 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); diff --git a/src/utils/__tests__/allocateImageFromChunks.spec.ts b/src/utils/__tests__/allocateImageFromChunks.spec.ts index f6c675246..a47bb2be8 100644 --- a/src/utils/__tests__/allocateImageFromChunks.spec.ts +++ b/src/utils/__tests__/allocateImageFromChunks.spec.ts @@ -2,7 +2,14 @@ import type { Chunk } from '@/src/core/streaming/chunk'; import { Tags } from '@/src/core/dicomTags'; import { allocateImageFromChunks, + getPixelFormat, + getRescaledValueRange, + getBufferValueRange, getTypedArrayForDataRange, + getTypedArrayValueRange, + getVolumeBufferType, + samplesAreIntegral, + valuesFitBuffer, } from '@/src/utils/allocateImageFromChunks'; import { describe, it, expect } from 'vitest'; @@ -41,7 +48,316 @@ describe('getTypedArrayForDataRange', () => { }); }); +// Builds a format from exactly the tags given, with no defaults filled in, so +// each test states the whole instance it is describing. +function format(tags: Record) { + return getPixelFormat(Object.entries(tags)); +} + +const UNSIGNED_16 = format({ + [Tags.BitsStored]: '16', + [Tags.PixelRepresentation]: '0', +}); +const SIGNED_16 = format({ + [Tags.BitsStored]: '16', + [Tags.PixelRepresentation]: '1', +}); +const UNSIGNED_8 = format({ + [Tags.BitsStored]: '8', + [Tags.PixelRepresentation]: '0', +}); +// A CT stored as 12-bit unsigned and rescaled to Hounsfield units. +const RESCALED_CT = format({ + [Tags.BitsStored]: '12', + [Tags.PixelRepresentation]: '0', + [Tags.RescaleSlope]: '1', + [Tags.RescaleIntercept]: '-1024', +}); +const NEGATIVE_SLOPE_8 = format({ + [Tags.BitsStored]: '8', + [Tags.PixelRepresentation]: '0', + [Tags.RescaleSlope]: '-1', +}); +const PUBLIC_DSC_SLOPE = '112067.85375182'; +const PUBLIC_DSC = format({ + [Tags.BitsStored]: '16', + [Tags.PixelRepresentation]: '0', + [Tags.RescaleSlope]: PUBLIC_DSC_SLOPE, + [Tags.RescaleIntercept]: '0', +}); + +describe('getPixelFormat', () => { + it('uses the DICOM defaults for absent rescale tags', () => { + expect(UNSIGNED_16).toEqual({ + bitsStored: 16, + pixelRepresentation: 0, + rescaleSlope: 1, + rescaleIntercept: 0, + }); + }); + + it('reports bitsStored 0 when the tag is absent', () => { + expect(format({ [Tags.PixelRepresentation]: '0' }).bitsStored).toBe(0); + }); + + it('treats a zero-length element as absent rather than as zero', () => { + expect( + format({ [Tags.BitsStored]: '16', [Tags.RescaleSlope]: '' }).rescaleSlope + ).toBe(1); + expect( + format({ [Tags.BitsStored]: '16', [Tags.RescaleSlope]: ' ' }) + .rescaleSlope + ).toBe(1); + expect( + format({ [Tags.BitsStored]: '', [Tags.PixelRepresentation]: '' }) + ).toEqual({ + bitsStored: 0, + pixelRepresentation: 0, + rescaleSlope: 1, + rescaleIntercept: 0, + }); + }); + + it('parses padded numeric strings', () => { + expect( + format({ [Tags.BitsStored]: '16', [Tags.RescaleIntercept]: ' -1024 ' }) + .rescaleIntercept + ).toBe(-1024); + }); +}); + +describe('getRescaledValueRange', () => { + it('spans the stored range put through slope and intercept', () => { + expect(getRescaledValueRange(UNSIGNED_16)).toEqual({ min: 0, max: 65535 }); + expect(getRescaledValueRange(SIGNED_16)).toEqual({ + min: -32768, + max: 32767, + }); + expect(getRescaledValueRange(RESCALED_CT)).toEqual({ + min: -1024, + max: 3071, + }); + }); + + it('orders the endpoints when the slope is negative', () => { + expect(getRescaledValueRange(NEGATIVE_SLOPE_8)).toEqual({ + min: -255, + max: 0, + }); + }); + + it('is null when BitsStored is missing', () => { + expect(getRescaledValueRange(format({}))).toBeNull(); + }); +}); + +describe('getVolumeBufferType', () => { + it('widens to hold the rescaled range', () => { + expect(getVolumeBufferType([UNSIGNED_8])).toBe(Uint8Array); + expect(getVolumeBufferType([UNSIGNED_16])).toBe(Uint16Array); + expect(getVolumeBufferType([SIGNED_16])).toBe(Int16Array); + expect(getVolumeBufferType([RESCALED_CT])).toBe(Int16Array); + expect(getVolumeBufferType([NEGATIVE_SLOPE_8])).toBe(Int16Array); + expect( + getVolumeBufferType([ + format({ [Tags.BitsStored]: '16', [Tags.RescaleIntercept]: '-1024' }), + ]) + ).toBe(Int32Array); + }); + + it('throws when BitsStored is missing', () => { + expect(() => getVolumeBufferType([format({})])).toThrow(); + }); + + it('matches the Float64 output observed in the public DSC perfusion series', () => { + expect(getVolumeBufferType([PUBLIC_DSC])).toBe(Float64Array); + }); + + it('uses Float64 when fractional rescale has integral endpoints', () => { + expect( + getVolumeBufferType([ + format({ + [Tags.BitsStored]: '12', + [Tags.RescaleSlope]: '0.2', + [Tags.RescaleIntercept]: '0', + }), + ]) + ).toBe(Float64Array); + }); + + it('uses Float64 when an integral modality range exceeds 32 bits', () => { + expect( + getVolumeBufferType([ + format({ + [Tags.BitsStored]: '16', + [Tags.RescaleSlope]: '112068', + [Tags.RescaleIntercept]: '0', + }), + ]) + ).toBe(Float64Array); + }); + + it('rejects a non-finite modality range', () => { + expect(() => + getVolumeBufferType([ + format({ + [Tags.BitsStored]: '16', + [Tags.PixelRepresentation]: '0', + [Tags.RescaleSlope]: '1e308', + }), + ]) + ).toThrow('No instance declares a finite modality rescale range'); + }); + + it('leaves an instance with no usable range out of the union', () => { + // The containment guard rejects such an instance when it decodes; its + // faulty tags must not stop the rest of the series from allocating. + expect(getVolumeBufferType([format({}), UNSIGNED_16])).toBe(Uint16Array); + expect( + getVolumeBufferType([ + UNSIGNED_16, + format({ [Tags.BitsStored]: '16', [Tags.RescaleSlope]: '1e308' }), + ]) + ).toBe(Uint16Array); + }); + + it('chooses one type that represents every instance in the volume', () => { + expect(getVolumeBufferType([UNSIGNED_8, UNSIGNED_16])).toBe(Uint16Array); + expect(getVolumeBufferType([UNSIGNED_16, PUBLIC_DSC])).toBe(Float64Array); + }); +}); + +describe('getTypedArrayValueRange', () => { + it('reports what each element type holds', () => { + expect(getTypedArrayValueRange(Uint8Array)).toEqual({ min: 0, max: 255 }); + expect(getTypedArrayValueRange(Int16Array)).toEqual({ + min: -32768, + max: 32767, + }); + }); + + it('has no range to report for element types the allocator never makes', () => { + expect(getTypedArrayValueRange(Float32Array)).toBeUndefined(); + }); +}); + +describe('getBufferValueRange', () => { + it('reads the range off the buffer itself', () => { + expect(getBufferValueRange(new Uint16Array(4))).toEqual({ + min: 0, + max: 65535, + }); + expect(getBufferValueRange(new Int32Array(4))).toEqual({ + min: -(2 ** 31), + max: 2 ** 31 - 1, + }); + }); +}); + +describe('valuesFitBuffer', () => { + it('accepts a range inside the buffer type', () => { + expect(valuesFitBuffer({ min: 0, max: 255 }, new Uint16Array(1))).toBe( + true + ); + expect(valuesFitBuffer({ min: -1024, max: 3071 }, new Int16Array(1))).toBe( + true + ); + }); + + it('accepts the exact bounds of the buffer type', () => { + expect(valuesFitBuffer({ min: 0, max: 65535 }, new Uint16Array(1))).toBe( + true + ); + expect( + valuesFitBuffer({ min: -32768, max: 32767 }, new Int16Array(1)) + ).toBe(true); + }); + + it('rejects negative values in an unsigned buffer', () => { + expect( + valuesFitBuffer({ min: -2048, max: -2048 }, new Uint16Array(1)) + ).toBe(false); + }); + + it('rejects values wider than the buffer type', () => { + expect(valuesFitBuffer({ min: 0, max: 5000 }, new Uint8Array(1))).toBe( + false + ); + expect(valuesFitBuffer({ min: 0, max: 65536 }, new Uint16Array(1))).toBe( + false + ); + }); + + it('accepts a slice whose declared range would overflow but whose values do not', () => { + // A slice declaring BitsStored 16 could reach 65535, which an Int16Array + // volume cannot hold, but what it actually decoded to fits. + const declaredRange = getRescaledValueRange(UNSIGNED_16)!; + expect(valuesFitBuffer(declaredRange, new Int16Array(1))).toBe(false); + expect(valuesFitBuffer({ min: 2000, max: 2000 }, new Int16Array(1))).toBe( + true + ); + }); + + it('judges the range only, leaving fractions to samplesAreIntegral', () => { + expect(valuesFitBuffer({ min: 0.5, max: 100.5 }, new Uint16Array(1))).toBe( + true + ); + }); + + it('constrains nothing when the buffer has no integer range', () => { + expect( + valuesFitBuffer({ min: -1e30, max: 1e30 }, new Float64Array(1)) + ).toBe(true); + }); +}); + +describe('samplesAreIntegral', () => { + it('trusts integer typed arrays without scanning them', () => { + expect( + samplesAreIntegral(new Uint16Array([1, 5000]), new Uint8Array(1)) + ).toBe(true); + }); + + it('rejects fractional float samples bound for an integer buffer', () => { + expect( + samplesAreIntegral(new Float64Array([0, 0.5, 1]), new Uint16Array(1)) + ).toBe(false); + expect(samplesAreIntegral([1, 1.5], new Int16Array(1))).toBe(false); + }); + + it('accepts float samples that are whole numbers', () => { + expect( + samplesAreIntegral(new Float64Array([0, 1, 2]), new Uint16Array(1)) + ).toBe(true); + expect( + samplesAreIntegral(new Float32Array([-3, 4]), new Int16Array(1)) + ).toBe(true); + }); + + it('accepts anything bound for a float buffer', () => { + expect( + samplesAreIntegral(new Float64Array([0.5, 1e30]), new Float64Array(1)) + ).toBe(true); + }); +}); + describe('allocateImageFromChunks', () => { + it('allocates for the modality range of every chunk', () => { + const image = allocateImageFromChunks([ + positionedChunk(0), + positionedChunk(1, { + [Tags.BitsStored]: '16', + [Tags.PixelRepresentation]: '0', + [Tags.RescaleSlope]: PUBLIC_DSC_SLOPE, + [Tags.RescaleIntercept]: '0', + }), + ]); + const data = image.getPointData().getScalars().getData(); + + expect(data).toBeInstanceOf(Float64Array); + expect(data).toHaveLength(3 * 4 * 2); + }); + it('matches ITK spacing order for single-slice images with SpacingBetweenSlices', () => { const image = allocateImageFromChunks([ chunk({ diff --git a/src/utils/allocateImageFromChunks.ts b/src/utils/allocateImageFromChunks.ts index d46f59ed7..77ffd0f70 100644 --- a/src/utils/allocateImageFromChunks.ts +++ b/src/utils/allocateImageFromChunks.ts @@ -36,6 +36,77 @@ function getBitStorageSize(num: number, signed: boolean) { return 2 ** Math.ceil(Math.log2(nbits)); } +type TypedArrayConstructor = + | typeof Int8Array + | typeof Uint8Array + | typeof Int16Array + | typeof Uint16Array + | typeof Int32Array + | typeof Uint32Array; + +const TYPED_ARRAY_VALUE_RANGES = new Map< + TypedArrayConstructor, + { min: number; max: number } +>([ + [Int8Array, { min: -(2 ** 7), max: 2 ** 7 - 1 }], + [Uint8Array, { min: 0, max: 2 ** 8 - 1 }], + [Int16Array, { min: -(2 ** 15), max: 2 ** 15 - 1 }], + [Uint16Array, { min: 0, max: 2 ** 16 - 1 }], + [Int32Array, { min: -(2 ** 31), max: 2 ** 31 - 1 }], + [Uint32Array, { min: 0, max: 2 ** 32 - 1 }], +]); + +/** + * The values a buffer of the given element type can hold without wrapping. + * Undefined for element types the allocator never produces, floats included. + */ +export function getTypedArrayValueRange(ctor: unknown) { + return TYPED_ARRAY_VALUE_RANGES.get(ctor as TypedArrayConstructor); +} + +/** + * The values `buffer` can hold without wrapping, or undefined if its element + * type has no fixed integer range to enforce. + */ +export function getBufferValueRange(buffer: ArrayBufferView) { + return getTypedArrayValueRange(buffer.constructor); +} + +/** + * Whether a decoded value range is representable in `buffer`'s element type. + * A buffer with no enforceable range accepts everything. + */ +export function valuesFitBuffer( + range: { min: number; max: number }, + buffer: ArrayBufferView +) { + const bufferRange = getBufferValueRange(buffer); + if (!bufferRange) return true; + return range.min >= bufferRange.min && range.max <= bufferRange.max; +} + +/** + * Whether every sample is a whole number, which an integer-element buffer + * needs. Integer typed arrays hold nothing else, so only float or plain + * arrays are scanned. A buffer with no integer range accepts everything. + */ +export function samplesAreIntegral( + values: ArrayLike, + buffer: ArrayBufferView +) { + if (!getBufferValueRange(buffer)) return true; + if ( + ArrayBuffer.isView(values) && + !(values instanceof Float32Array) && + !(values instanceof Float64Array) + ) { + return true; + } + return Array.prototype.every.call(values, (value: number) => + Number.isInteger(value) + ); +} + export function getTypedArrayForDataRange(min: number, max: number) { if (!Number.isSafeInteger(min) || !Number.isSafeInteger(max)) throw new Error('Input must be integers'); @@ -53,24 +124,85 @@ export function getTypedArrayForDataRange(min: number, max: number) { throw new Error(`Cannot handle ${nbits}-bit sized ranges`); } -function getTypedArrayConstructor( - bitsStored: number, - pixelRepresentation: number, - rescaleIntercept: number, - rescaleSlope: number -) { - if (bitsStored === 0) throw new Error('bits stored is zero!'); +function numberOr(value: Maybe, fallback: number) { + const text = value?.trim(); + if (!text) return fallback; + const num = Number(text); + return Number.isFinite(num) ? num : fallback; +} + +/** + * The tags that decide what values an instance decodes to. + * + * A bitsStored of 0 means the tag was absent or unparseable. It is not a DICOM + * default: BitsStored is Type 1, so a conforming instance always carries one. + */ +export function getPixelFormat(metadata: Maybe>) { + const meta = new Map(metadata ?? []); + return { + bitsStored: numberOr(meta.get(BitsStoredTag), 0), + pixelRepresentation: numberOr(meta.get(PixelRepresentationTag), 0), + rescaleSlope: numberOr(meta.get(RescaleSlope), 1), + rescaleIntercept: numberOr(meta.get(RescaleIntercept), 0), + }; +} + +type PixelFormat = ReturnType; + +/** + * The values an instance can decode to. ITK/GDCM applies RescaleSlope and + * RescaleIntercept while decoding, so this is the range after rescaling. + * + * Null when BitsStored is missing, since nothing can be derived without it. + */ +export function getRescaledValueRange(format: PixelFormat) { + const { bitsStored, pixelRepresentation, rescaleSlope, rescaleIntercept } = + format; + if (!Number.isInteger(bitsStored) || bitsStored <= 0) return null; - // Maybe constrain bitsAllocated to allowed values of 8, 16, 32? const isSigned = pixelRepresentation === 1; const storedMin = isSigned ? -(2 ** (bitsStored - 1)) : 0; const storedMax = 2 ** (bitsStored - (isSigned ? 1 : 0)) - 1; - const outputMin = Math.floor(storedMin * rescaleSlope + rescaleIntercept); - const outputMax = Math.ceil(storedMax * rescaleSlope + rescaleIntercept); + const a = storedMin * rescaleSlope + rescaleIntercept; + const b = storedMax * rescaleSlope + rescaleIntercept; + + return { min: Math.min(a, b), max: Math.max(a, b) }; +} + +/** + * The element type a volume allocated from instances of these formats holds. + * + * An instance whose tags give no finite range is left out of the union and + * judged on its decoded values by the containment guard instead. + */ +export function getVolumeBufferType(formats: PixelFormat[]) { + const usable = formats.flatMap((format) => { + const range = getRescaledValueRange(format); + return range && Number.isFinite(range.min) && Number.isFinite(range.max) + ? [{ format, range }] + : []; + }); + if (usable.length === 0) + throw new Error('No instance declares a finite modality rescale range'); + + const needsFloat64 = usable.some( + ({ format, range }) => + !Number.isInteger(format.rescaleSlope) || + !Number.isInteger(format.rescaleIntercept) || + !Number.isSafeInteger(range.min) || + !Number.isSafeInteger(range.max) + ); + if (needsFloat64) return Float64Array; + + const min = Math.min(...usable.map(({ range }) => range.min)); + const max = Math.max(...usable.map(({ range }) => range.max)); + const exceedsSigned32 = min < 0 && (min < -(2 ** 31) || max > 2 ** 31 - 1); + const exceedsUnsigned32 = min >= 0 && max > 2 ** 32 - 1; + if (exceedsSigned32 || exceedsUnsigned32) return Float64Array; // NOTE(fli): might be better to assume (u)int16 and re-allocate to (u)int32 // if needed, since the data range might actually fit in a smaller datatype. - return getTypedArrayForDataRange(outputMin, outputMax); + return getTypedArrayForDataRange(min, max); } export function allocateImageFromChunks(sortedChunks: Chunk[]) { @@ -86,11 +218,10 @@ export function allocateImageFromChunks(sortedChunks: Chunk[]) { const spacingBetweenSlices = Number(meta.get(SpacingBetweenSlicesTag)); const rows = Number(meta.get(RowsTag) ?? 0); const columns = Number(meta.get(ColumnsTag) ?? 0); - const bitsStored = Number(meta.get(BitsStoredTag) ?? 0); - const pixelRepresentation = Number(meta.get(PixelRepresentationTag)); + const volumeFormats = sortedChunks.map((chunk) => + getPixelFormat(chunk.metadata) + ); const samplesPerPixel = Number(meta.get(SamplesPerPixelTag) ?? 1); - const rescaleIntercept = Number(meta.get(RescaleIntercept) ?? 0); - const rescaleSlope = Number(meta.get(RescaleSlope) ?? 1); const numberOfFrames = meta.has(NumberOfFrames) ? Number(meta.get(NumberOfFrames)) : null; @@ -109,12 +240,7 @@ export function allocateImageFromChunks(sortedChunks: Chunk[]) { // Some CT modality series have NumberOfFrames === 1, so use the number of chunks if more than 1 chunk. const slices = sortedChunks.length > 1 ? sortedChunks.length : (numberOfFrames ?? 1); - const TypedArrayCtor = getTypedArrayConstructor( - bitsStored, - pixelRepresentation, - rescaleIntercept, - rescaleSlope - ); + const TypedArrayCtor = getVolumeBufferType(volumeFormats); const pixelData = new TypedArrayCtor( rows * columns * slices * samplesPerPixel ); diff --git a/tests/specs/dicom-modality-rescale.e2e.ts b/tests/specs/dicom-modality-rescale.e2e.ts new file mode 100644 index 000000000..9b7225ee2 --- /dev/null +++ b/tests/specs/dicom-modality-rescale.e2e.ts @@ -0,0 +1,71 @@ +import * as fs from 'fs'; +import * as path from 'path'; +import { cleanuptotal } from 'wdio-cleanuptotal-service'; + +import { TEMP_DIR } from '../../wdio.shared.conf'; +import { volViewPage } from '../pageobjects/volview.page'; +import { buildSyntheticDicom, newUid } from './syntheticDicom'; +import { waitForFirstCompleteCachedImageScalars } from './imageCacheUtils'; +import { writeManifestToFile } from './utils'; + +const PUBLIC_DSC_SERIES_UID = + '1.3.6.1.4.1.9590.100.1.2.284777661700890778225181143863199482857'; +const PUBLIC_DSC_SLOPE = 112067.85375182; +const STORED_VALUES = [0, 2, 65131]; +const ROWS = 4; +const COLUMNS = 4; + +async function writeRescaledSeries() { + const dirName = `modality-rescale-${Date.now()}`; + const dir = path.join(TEMP_DIR, dirName); + fs.mkdirSync(dir, { recursive: true }); + cleanuptotal.addCleanup(async () => { + fs.rmSync(dir, { recursive: true, force: true }); + }); + + const studyUid = newUid(); + const resources = STORED_VALUES.map((pixelValue, index) => { + const filename = `slice-${index}.dcm`; + fs.writeFileSync( + path.join(dir, filename), + buildSyntheticDicom({ + studyUid, + seriesUid: PUBLIC_DSC_SERIES_UID, + sopUid: newUid(), + instanceNumber: index + 1, + imageOrientationPatient: [1, 0, 0, 0, 1, 0], + imagePositionPatient: [0, 0, index], + rows: ROWS, + cols: COLUMNS, + bitsAllocated: 16, + bitsStored: 16, + highBit: 15, + pixelRepresentation: 0, + rescaleSlope: PUBLIC_DSC_SLOPE, + rescaleIntercept: 0, + pixelValue, + }) + ); + return { url: `tmp/${dirName}/${filename}`, name: filename }; + }); + + const manifestName = `modality-rescale-${Date.now()}.json`; + await writeManifestToFile({ resources }, manifestName); + return manifestName; +} + +describe('DICOM modality rescale', () => { + it('preserves the public DSC series Float64 output through volume loading', async () => { + const manifestName = await writeRescaledSeries(); + await volViewPage.open(`?urls=[tmp/${manifestName}]`); + await volViewPage.waitForViews(); + const scalars = await waitForFirstCompleteCachedImageScalars(); + + expect(await volViewPage.getNotificationsCount()).toBe(0); + expect(scalars.type).toBe('Float64Array'); + const expected = STORED_VALUES.flatMap((stored) => + Array(ROWS * COLUMNS).fill(stored * PUBLIC_DSC_SLOPE) + ); + expect(scalars.values).toEqual(expected); + }); +}); diff --git a/tests/specs/imageCacheUtils.ts b/tests/specs/imageCacheUtils.ts index e006ee030..78f680e11 100644 --- a/tests/specs/imageCacheUtils.ts +++ b/tests/specs/imageCacheUtils.ts @@ -1,5 +1,7 @@ -export async function getFirstCachedImageSpacing() { - return browser.execute(() => { +type CachedImageScalars = { type: string; values: number[] }; + +async function readFirstCachedImage(property: 'spacing' | 'complete-scalars') { + return browser.execute((requestedProperty) => { const app = (document.querySelector('#app') as any)?.__vue_app__; const pinia = app?.config?.globalProperties?.$pinia ?? @@ -15,8 +17,22 @@ export async function getFirstCachedImageSpacing() { const id = imageCache?.imageIds?.[0]; const imageData = imageCache?.getVtkImageData(id); if (!imageData) return null; - return Array.from(imageData.getSpacing()).map(Number); - }); + if (requestedProperty === 'spacing') { + return Array.from(imageData.getSpacing()).map(Number); + } + if (imageCache.imageStatus[id] !== 'complete') return null; + + const data = imageData.getPointData().getScalars()?.getData(); + if (!data) return null; + return { + type: data.constructor.name, + values: Array.from(data as ArrayLike), + }; + }, property); +} + +export function getFirstCachedImageSpacing() { + return readFirstCachedImage('spacing') as Promise; } export async function waitForFirstCachedImageSpacing() { @@ -33,3 +49,24 @@ export async function waitForFirstCachedImageSpacing() { ); return spacing!; } + +export function getFirstCompleteCachedImageScalars() { + return readFirstCachedImage( + 'complete-scalars' + ) as Promise; +} + +export async function waitForFirstCompleteCachedImageScalars() { + let scalars: CachedImageScalars | null = null; + await browser.waitUntil( + async () => { + scalars = await getFirstCompleteCachedImageScalars(); + return scalars !== null; + }, + { + timeout: 30_000, + timeoutMsg: 'Expected first cached image scalars to become available', + } + ); + return scalars!; +} diff --git a/tests/specs/syntheticDicom.ts b/tests/specs/syntheticDicom.ts index 70744bcee..f1750e94a 100644 --- a/tests/specs/syntheticDicom.ts +++ b/tests/specs/syntheticDicom.ts @@ -2,7 +2,7 @@ // Emits just enough tags for ITK/GDCM to categorize and load a series: // SOP Class/Instance UIDs, Study/SeriesInstanceUID, SeriesNumber, Modality, // Patient identifiers, ImageOrientationPatient, ImagePositionPatient, -// PixelSpacing, SliceThickness, image geometry, and zeroed PixelData. +// PixelSpacing, SliceThickness, image geometry, and constant-valued PixelData. const SOP_CLASS_MR = '1.2.840.10008.5.1.4.1.1.4'; const SOP_CLASS_ULTRASOUND_MULTIFRAME = '1.2.840.10008.5.1.4.1.1.3.1'; @@ -92,6 +92,22 @@ const ds = (g: number, e: number, v: string) => const us = (g: number, e: number, v: number) => elemShort(g, e, 'US', writeShort(v)); +// DICOM element values must have an even length, so an odd-sized 8-bit frame +// gets a trailing pad byte. +const frameBytes8 = (sampleCount: number, value: number) => { + const bytes = new Uint8Array(sampleCount + (sampleCount % 2)); + bytes.fill(value & 0xff, 0, sampleCount); + return bytes; +}; + +const frameBytes16 = (sampleCount: number, value: number) => { + const bytes = new Uint8Array(sampleCount * 2); + const view = new DataView(bytes.buffer); + for (let i = 0; i < sampleCount; i++) + view.setUint16(i * 2, value & 0xffff, true); + return bytes; +}; + export type SyntheticSliceOptions = { studyUid: string; seriesUid: string; @@ -108,6 +124,15 @@ export type SyntheticSliceOptions = { imagePositionPatient: readonly [number, number, number]; rows?: number; cols?: number; + bitsAllocated?: number; + bitsStored?: number; + highBit?: number; + pixelRepresentation?: number; + rescaleSlope?: number; + rescaleIntercept?: number; + // Stored value written to every sample. Signed values are written as two's + // complement, so they pair with pixelRepresentation 1. + pixelValue?: number; pixelSpacing?: readonly [number, number]; spacingBetweenSlices?: number; sliceThickness?: number; @@ -128,6 +153,13 @@ export function buildSyntheticDicom(opts: SyntheticSliceOptions): Uint8Array { imagePositionPatient, rows = 4, cols = 4, + bitsAllocated = 16, + bitsStored = bitsAllocated, + highBit = bitsStored - 1, + pixelRepresentation = 0, + rescaleSlope, + rescaleIntercept, + pixelValue = 0, pixelSpacing = [1, 1] as const, spacingBetweenSlices, sliceThickness = 1, @@ -138,6 +170,12 @@ export function buildSyntheticDicom(opts: SyntheticSliceOptions): Uint8Array { studyDate = '20260101', } = opts; + if (bitsAllocated !== 8 && bitsAllocated !== 16) { + throw new Error( + `bitsAllocated must be 8 or 16, got ${bitsAllocated}. Other widths would need a pixel data VR this helper does not emit.` + ); + } + const dataset = combine( ui(0x0008, 0x0016, SOP_CLASS_MR), ui(0x0008, 0x0018, sopUid), @@ -172,11 +210,17 @@ export function buildSyntheticDicom(opts: SyntheticSliceOptions): Uint8Array { us(0x0028, 0x0010, rows), us(0x0028, 0x0011, cols), ds(0x0028, 0x0030, pixelSpacing.map((n) => n.toString()).join('\\')), - us(0x0028, 0x0100, 16), - us(0x0028, 0x0101, 16), - us(0x0028, 0x0102, 15), - us(0x0028, 0x0103, 0), - elemLong(0x7fe0, 0x0010, 'OW', new Uint8Array(rows * cols * 2)) + us(0x0028, 0x0100, bitsAllocated), + us(0x0028, 0x0101, bitsStored), + us(0x0028, 0x0102, highBit), + us(0x0028, 0x0103, pixelRepresentation), + ...(rescaleIntercept == null + ? [] + : [ds(0x0028, 0x1052, String(rescaleIntercept))]), + ...(rescaleSlope == null ? [] : [ds(0x0028, 0x1053, String(rescaleSlope))]), + bitsAllocated === 8 + ? elemLong(0x7fe0, 0x0010, 'OB', frameBytes8(rows * cols, pixelValue)) + : elemLong(0x7fe0, 0x0010, 'OW', frameBytes16(rows * cols, pixelValue)) ); const fileMetaBody = combine(