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(