Skip to content

Commit 12b54f7

Browse files
PaulHaxPaulHax
committed
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.
1 parent f482c6d commit 12b54f7

7 files changed

Lines changed: 888 additions & 59 deletions

File tree

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

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
import { describe, it, expect, vi } from 'vitest';
2+
import { Chunk } from '@/src/core/streaming/chunk';
3+
import { ChunkState } from '@/src/core/streaming/chunkStateMachine';
4+
import { Tags } from '@/src/core/dicomTags';
5+
import DicomChunkImage, {
6+
DicomChunkImageInit,
7+
} from '@/src/core/streaming/dicomChunkImage';
8+
import { ChunkStatus } from '@/src/core/streaming/chunkImage';
9+
10+
const ROWS = 2;
11+
const COLUMNS = 2;
12+
const PIXELS_PER_SLICE = ROWS * COLUMNS;
13+
const PUBLIC_DSC_SLOPE = 112067.85375182;
14+
15+
function metadataFor(z: number, overrides: Record<string, string> = {}) {
16+
const metadata = [
17+
[Tags.SOPInstanceUID, `1.2.3.${z}`],
18+
[Tags.ImagePositionPatient, `0\\0\\${z}`],
19+
[Tags.ImageOrientationPatient, '1\\0\\0\\0\\1\\0'],
20+
[Tags.Rows, String(ROWS)],
21+
[Tags.Columns, String(COLUMNS)],
22+
[Tags.PixelSpacing, '1\\1'],
23+
[Tags.BitsStored, '16'],
24+
[Tags.PixelRepresentation, '0'],
25+
[Tags.SamplesPerPixel, '1'],
26+
] as Array<[string, string]>;
27+
Object.entries(overrides).forEach(([tag, value]) => {
28+
const existing = metadata.find((entry) => entry[0] === tag);
29+
if (existing) existing[1] = value;
30+
else metadata.push([tag, value]);
31+
});
32+
return metadata;
33+
}
34+
35+
// The slice's z position is also its pixel value, so the decoded contents of a
36+
// slice identify which chunk it came from.
37+
async function makeLoadedChunk(
38+
z: number,
39+
overrides: Record<string, string> = {}
40+
) {
41+
const meta = metadataFor(z, overrides);
42+
const chunk = new Chunk({
43+
metaLoader: {
44+
meta,
45+
metaBlob: new Blob([`meta-${z}`]),
46+
load: () => {},
47+
stop: () => {},
48+
},
49+
dataLoader: {
50+
data: new Blob([String(z)]),
51+
load: () => {},
52+
stop: () => {},
53+
},
54+
});
55+
await chunk.loadMeta();
56+
await chunk.loadData();
57+
expect(chunk.state).toBe(ChunkState.Loaded);
58+
return chunk;
59+
}
60+
61+
function zOf(chunk: Chunk) {
62+
const meta = Object.fromEntries(chunk.metadata!);
63+
return Number(meta[Tags.ImagePositionPatient].split('\\')[2]);
64+
}
65+
66+
function splitAndSortByPosition(chunks: Chunk[]) {
67+
return Promise.resolve({
68+
volume: [...chunks].sort((a, b) => zOf(a) - zOf(b)),
69+
});
70+
}
71+
72+
// Decodes a chunk to a constant frame, letting the test choose the array type.
73+
function decodeTo(dataFor: (value: number) => ArrayLike<number>) {
74+
const read: DicomChunkImageInit['readDicomImage'] = async (file) => {
75+
const value = Number(await file.text());
76+
return {
77+
image: {
78+
size: [COLUMNS, ROWS, 1],
79+
data: dataFor(value) as Uint16Array,
80+
imageType: { components: 1 },
81+
},
82+
};
83+
};
84+
return read;
85+
}
86+
87+
function sliceOf(image: DicomChunkImage, index: number) {
88+
const scalars = image.getVtkImageData().getPointData().getScalars();
89+
const data = scalars.getData();
90+
return Array.from(
91+
data.slice(index * PIXELS_PER_SLICE, (index + 1) * PIXELS_PER_SLICE)
92+
);
93+
}
94+
95+
async function loadRejectingSeries(
96+
readDicomImage: DicomChunkImageInit['readDicomImage']
97+
) {
98+
const image = new DicomChunkImage({
99+
splitAndSort: splitAndSortByPosition,
100+
readDicomImage,
101+
});
102+
const errors: unknown[] = [];
103+
image.addEventListener('chunkError', ({ error }) => {
104+
errors.push(error);
105+
});
106+
const [valid, invalid] = await Promise.all([
107+
makeLoadedChunk(1, { [Tags.BitsStored]: '8' }),
108+
makeLoadedChunk(2, { [Tags.BitsStored]: '8' }),
109+
]);
110+
111+
await image.addChunks([valid, invalid]);
112+
await vi.waitFor(() =>
113+
expect(image.getChunkStatuses()).toEqual([
114+
ChunkStatus.Loaded,
115+
ChunkStatus.Errored,
116+
])
117+
);
118+
119+
expect(image.status.value).toBe('complete');
120+
expect(errors).toHaveLength(1);
121+
expect(sliceOf(image, 0)).toEqual(Array(PIXELS_PER_SLICE).fill(1));
122+
expect(sliceOf(image, 1)).toEqual(Array(PIXELS_PER_SLICE).fill(0));
123+
image.dispose();
124+
return String(errors[0]);
125+
}
126+
127+
describe('DicomChunkImage', () => {
128+
it('preserves exact modality-rescaled pixels from the public DSC series', async () => {
129+
// The public frames are 200x230; reduced geometry keeps the exact encoding,
130+
// rescale, and an observed stored-pixel maximum in a focused volume test.
131+
const decoded = Float64Array.from([
132+
0,
133+
PUBLIC_DSC_SLOPE,
134+
2 * PUBLIC_DSC_SLOPE,
135+
65131 * PUBLIC_DSC_SLOPE,
136+
]);
137+
const image = new DicomChunkImage({
138+
splitAndSort: splitAndSortByPosition,
139+
readDicomImage: decodeTo(() => decoded),
140+
});
141+
const frame = await makeLoadedChunk(1, {
142+
[Tags.SeriesInstanceUID]:
143+
'1.3.6.1.4.1.9590.100.1.2.284777661700890778225181143863199482857',
144+
[Tags.RescaleSlope]: String(PUBLIC_DSC_SLOPE),
145+
[Tags.RescaleIntercept]: '0',
146+
});
147+
148+
await image.addChunks([frame]);
149+
await vi.waitFor(() =>
150+
expect(image.getChunkStatuses()).toEqual([ChunkStatus.Loaded])
151+
);
152+
153+
const data = image.getVtkImageData().getPointData().getScalars().getData();
154+
expect(data).toBeInstanceOf(Float64Array);
155+
expect(Array.from(data)).toEqual(Array.from(decoded));
156+
157+
image.dispose();
158+
});
159+
160+
it('settles after rejecting decoded values its integer buffer cannot hold', async () => {
161+
const message = await loadRejectingSeries(
162+
decodeTo((value) =>
163+
value === 2
164+
? new Uint16Array(PIXELS_PER_SLICE).fill(5000)
165+
: new Uint8Array(PIXELS_PER_SLICE).fill(value)
166+
)
167+
);
168+
expect(message).toContain('5000');
169+
expect(message).toContain('Uint8Array');
170+
});
171+
172+
it('settles after rejecting fractional samples bound for an integer buffer', async () => {
173+
const message = await loadRejectingSeries(
174+
decodeTo((value) =>
175+
value === 2
176+
? new Float64Array(PIXELS_PER_SLICE).fill(2.5)
177+
: new Uint8Array(PIXELS_PER_SLICE).fill(value)
178+
)
179+
);
180+
expect(message).toContain('fractional');
181+
expect(message).toContain('Uint8Array');
182+
});
183+
});

‎src/core/streaming/dicomChunkImage.ts‎

Lines changed: 80 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -2,12 +2,21 @@ import {
22
buildSegmentGroups,
33
ReadOverlappingSegmentationMeta,
44
readVolumeSlice,
5-
splitAndSort,
5+
splitAndSort as splitAndSortChunks,
66
} from '@/src/io/dicom';
77
import { Chunk, waitForChunkState } from '@/src/core/streaming/chunk';
8-
import { Image, JsonCompatible, readImage } from '@itk-wasm/image-io';
8+
import {
9+
Image,
10+
JsonCompatible,
11+
readImage as readItkImage,
12+
} from '@itk-wasm/image-io';
913
import { getWorker } from '@/src/io/itk/worker';
10-
import { allocateImageFromChunks } from '@/src/utils/allocateImageFromChunks';
14+
import {
15+
allocateImageFromChunks,
16+
getBufferValueRange,
17+
samplesAreIntegral,
18+
valuesFitBuffer,
19+
} from '@/src/utils/allocateImageFromChunks';
1120
import { TypedArray } from '@kitware/vtk.js/types';
1221
import { Tags } from '@/src/core/dicomTags';
1322
import vtkDataArray from '@kitware/vtk.js/Common/Core/DataArray';
@@ -71,10 +80,28 @@ async function dicomSliceToImageUri(blob: Blob) {
7180
return itkImageToURI(itkImage);
7281
}
7382

83+
function readDicomImage(file: File) {
84+
return readItkImage(file, { webWorker: getWorker() });
85+
}
86+
87+
export interface DicomChunkImageInit {
88+
splitAndSort: (
89+
chunks: Chunk[],
90+
mapToBlob: (chunk: Chunk, index: number) => Blob
91+
) => Promise<Record<string, Chunk[]>>;
92+
readDicomImage: (file: File) => Promise<{
93+
image: Pick<Image, 'size' | 'data'> & {
94+
imageType: Pick<Image['imageType'], 'components'>;
95+
};
96+
}>;
97+
}
98+
7499
export default class DicomChunkImage
75100
extends BaseProgressiveImage
76101
implements ChunkImage
77102
{
103+
private splitAndSort: DicomChunkImageInit['splitAndSort'];
104+
private readDicomImage: DicomChunkImageInit['readDicomImage'];
78105
protected chunks: Chunk[];
79106
private chunkListeners: Array<() => void>;
80107
private thumbnailCache: WeakMap<Chunk, Promise<string>>;
@@ -85,9 +112,12 @@ export default class DicomChunkImage
85112
| (JsonCompatible & ReadOverlappingSegmentationMeta)
86113
| null;
87114

88-
constructor() {
115+
constructor(init: Partial<DicomChunkImageInit> = {}) {
89116
super();
90117

118+
this.splitAndSort = init.splitAndSort ?? splitAndSortChunks;
119+
this.readDicomImage = init.readDicomImage ?? readDicomImage;
120+
91121
this.status.value = 'incomplete';
92122
this.loaded = computed(() => {
93123
return !this.loading.value && this.status.value === 'complete';
@@ -179,7 +209,7 @@ export default class DicomChunkImage
179209
});
180210

181211
await Promise.all(chunks.map((chunk) => chunk.loadMeta()));
182-
const chunksByVolume = await splitAndSort(
212+
const chunksByVolume = await this.splitAndSort(
183213
this.chunks,
184214
(chunk) => chunk.metaBlob!
185215
);
@@ -385,11 +415,8 @@ export default class DicomChunkImage
385415
throw new Error(`Chunk ${chunkIndex} does not have data`);
386416

387417
const chunkId = chunk.metadata ? getChunkId(chunk) : `index-${chunkIndex}`;
388-
const result = await readImage(
389-
new File([chunk.dataBlob], `file-${chunkIndex}.dcm`),
390-
{
391-
webWorker: getWorker(),
392-
}
418+
const result = await this.readDicomImage(
419+
new File([chunk.dataBlob], `file-${chunkIndex}.dcm`)
393420
);
394421

395422
if (!result.image.data)
@@ -405,13 +432,12 @@ export default class DicomChunkImage
405432

406433
const scalars = this.vtkImageData.value.getPointData().getScalars();
407434
const pixelData = scalars.getData() as TypedArray;
435+
const componentCount = scalars.getNumberOfComponents();
408436

409437
const dims = this.vtkImageData.value.getDimensions();
410-
const components = scalars.getNumberOfComponents();
411438

412-
// The volume buffer is sized from the first chunk's metadata, so each
413-
// chunk gets a fixed slot: one frame per chunk in a multi-file volume,
414-
// the whole volume when a single multi-frame chunk fills it.
439+
// Each chunk gets a fixed slot: one frame per chunk in a multi-file
440+
// volume, or the whole volume when a single multi-frame chunk fills it.
415441
const framesPerChunk = this.chunks.length > 1 ? 1 : dims[2];
416442
const [chunkWidth, chunkHeight] = result.image.size;
417443
const chunkFrames = result.image.size[2] ?? 1;
@@ -420,7 +446,7 @@ export default class DicomChunkImage
420446
chunkWidth !== dims[0] ||
421447
chunkHeight !== dims[1] ||
422448
chunkFrames !== framesPerChunk ||
423-
chunkComponents !== components
449+
chunkComponents !== componentCount
424450
) {
425451
// A lone chunk defines the volume it fails to fit, so advice about
426452
// agreeing with the other files only makes sense for a multi-file volume.
@@ -431,34 +457,60 @@ export default class DicomChunkImage
431457
throw new Error(
432458
`File ${chunkId} (chunk ${chunkIndex}) does not fit the volume it belongs to. ` +
433459
`It decoded to ${chunkWidth}x${chunkHeight}x${chunkFrames} with ${chunkComponents} component(s), ` +
434-
`but the volume has room for ${dims[0]}x${dims[1]}x${framesPerChunk} with ${components} component(s).` +
460+
`but the volume has room for ${dims[0]}x${dims[1]}x${framesPerChunk} with ${componentCount} component(s).` +
435461
advice
436462
);
437463
}
438464

439-
const offset = dims[0] * dims[1] * components * chunkIndex;
440-
pixelData.set(result.image.data as TypedArray, offset);
441-
442-
const rangeAlreadyInitialized = this.chunkStatus.some(
443-
(status) => status === ChunkStatus.Loaded
444-
);
445-
446-
// update the data range
447465
const chunkDataRange: Array<[number, number]> = [];
448-
for (let comp = 0; comp < scalars.getNumberOfComponents(); comp++) {
466+
for (let comp = 0; comp < componentCount; comp++) {
449467
const { min, max } = fastComputeRange(
450468
result.image.data as unknown as number[],
451469
comp,
452-
scalars.getNumberOfComponents()
470+
componentCount
453471
);
454472
chunkDataRange.push([min, max]);
473+
}
455474

456-
const curRange = scalars.getRange(comp);
475+
// The buffer is allocated for the range every chunk's tags declare, so a
476+
// chunk only fails here when its decoded values disagree with its tags.
477+
// TypedArray.set raises nothing for such values: integers wrap and
478+
// fractions truncate.
479+
const chunkMin = Math.min(...chunkDataRange.map(([min]) => min));
480+
const chunkMax = Math.max(...chunkDataRange.map(([, max]) => max));
481+
const decoded = result.image.data as unknown as ArrayLike<number>;
482+
if (!valuesFitBuffer({ min: chunkMin, max: chunkMax }, pixelData)) {
483+
const bufferRange = getBufferValueRange(pixelData)!;
484+
throw new Error(
485+
`File ${chunkId} (chunk ${chunkIndex}) has pixel values the volume it belongs to cannot represent. ` +
486+
`Its pixel values run from ${chunkMin} to ${chunkMax}, but the volume's buffer is ` +
487+
`${pixelData.constructor.name}, holding values from ${bufferRange.min} to ${bufferRange.max}. ` +
488+
`Every file in a volume must decode to values its buffer can hold without conversion.`
489+
);
490+
}
491+
if (!samplesAreIntegral(decoded, pixelData)) {
492+
throw new Error(
493+
`File ${chunkId} (chunk ${chunkIndex}) has fractional pixel values the volume it belongs to cannot represent. ` +
494+
`Its pixel values run from ${chunkMin} to ${chunkMax}, but the volume's buffer is ` +
495+
`${pixelData.constructor.name}, which holds only whole numbers. ` +
496+
`Every file in a volume must decode to values its buffer can hold without conversion.`
497+
);
498+
}
499+
500+
const offset = dims[0] * dims[1] * componentCount * chunkIndex;
501+
pixelData.set(result.image.data as TypedArray, offset);
457502

503+
const rangeAlreadyInitialized = this.chunkStatus.some(
504+
(status) => status === ChunkStatus.Loaded
505+
);
506+
507+
// update the data range
508+
chunkDataRange.forEach(([min, max], comp) => {
509+
const curRange = scalars.getRange(comp);
458510
const newMin = rangeAlreadyInitialized ? Math.min(min, curRange[0]) : min;
459511
const newMax = rangeAlreadyInitialized ? Math.max(max, curRange[1]) : max;
460512
scalars.setRange({ min: newMin, max: newMax }, comp);
461-
}
513+
});
462514
scalars.modified(); // so image-stats will trigger update of range
463515

464516
chunk.setUserData(DATA_RANGE_KEY, chunkDataRange);

0 commit comments

Comments
 (0)