From bcc3c7fc7dde1d1e6369c356666076b786540bdf Mon Sep 17 00:00:00 2001 From: Paul Elliott Date: Sun, 20 Sep 2026 12:32:06 +0000 Subject: [PATCH 1/3] feat(segmentation): add the segment model, masks, algorithms and render helpers Adds the framework-free layer of the segmentation module: the segment and mask model, geometry helpers, colour parsing, bounded mask storage with overlap planning and voxel access, the fill-holes and gaussian smooth algorithms with their workers, the shared identity registry, the reference coordinator, and the display, projection and render-mask helpers. Nothing consumes them yet; the store, persistence, tools and UI that do follow in the next change. Also adds reorientLabelImage, which ensureSameSpace uses to reorder an equivalent voxel grid onto its image without interpolating label boundaries. --- package-lock.json | 9 + package.json | 2 + .../__tests__/ensureSameSpace.spec.ts} | 19 +- .../__tests__/reorientLabelImage.spec.ts | 103 +++++ src/io/resample/reorientLabelImage.ts | 82 ++++ src/io/resample/resample.ts | 8 +- src/segmentation/__tests__/model.spec.ts | 297 ++++++++++++++ .../__tests__/reframeMaskScalars.spec.ts | 44 ++ src/segmentation/color.ts | 41 ++ .../__tests__/gaussianSmooth.spec.ts | 122 ++++++ .../__tests__/gaussianSmoothGolden.spec.ts | 237 +++++++++++ .../editing/algorithms/fillHoles.ts | 165 ++++++++ .../editing/algorithms/fillHoles.worker.ts | 17 + .../algorithms/gaussianSmooth.worker.ts | 384 ++++++++++++++++++ src/segmentation/editing/coordinator.ts | 29 ++ src/segmentation/geometry.ts | 174 ++++++++ src/segmentation/io/__tests__/export.spec.ts | 58 +++ src/segmentation/io/export.ts | 34 ++ src/segmentation/io/maskFileNaming.ts | 27 ++ src/segmentation/masks/labelValue.ts | 33 ++ src/segmentation/masks/overlap.ts | 197 +++++++++ src/segmentation/masks/storage.ts | 101 +++++ src/segmentation/masks/voxelAccess.ts | 216 ++++++++++ src/segmentation/model.ts | 127 ++++++ .../__tests__/segmentDisplay.spec.ts | 129 ++++++ .../__tests__/segmentRenderMask.spec.ts | 114 ++++++ .../__tests__/segmentSliceVisibility.spec.ts | 87 ++++ src/segmentation/rendering/display.ts | 87 ++++ src/segmentation/rendering/projection.ts | 60 +++ src/segmentation/rendering/renderMask.ts | 83 ++++ src/segmentation/segment.ts | 97 +++++ src/segmentation/segmentReferences.ts | 40 ++ src/segmentation/segmentRegistry.ts | 276 +++++++++++++ src/utils/color.ts | 10 + 34 files changed, 3492 insertions(+), 17 deletions(-) rename src/{store/__tests__/segmentImportPath.spec.ts => io/resample/__tests__/ensureSameSpace.spec.ts} (75%) create mode 100644 src/io/resample/__tests__/reorientLabelImage.spec.ts create mode 100644 src/io/resample/reorientLabelImage.ts create mode 100644 src/segmentation/__tests__/model.spec.ts create mode 100644 src/segmentation/__tests__/reframeMaskScalars.spec.ts create mode 100644 src/segmentation/color.ts create mode 100644 src/segmentation/editing/algorithms/__tests__/gaussianSmooth.spec.ts create mode 100644 src/segmentation/editing/algorithms/__tests__/gaussianSmoothGolden.spec.ts create mode 100644 src/segmentation/editing/algorithms/fillHoles.ts create mode 100644 src/segmentation/editing/algorithms/fillHoles.worker.ts create mode 100644 src/segmentation/editing/algorithms/gaussianSmooth.worker.ts create mode 100644 src/segmentation/editing/coordinator.ts create mode 100644 src/segmentation/geometry.ts create mode 100644 src/segmentation/io/__tests__/export.spec.ts create mode 100644 src/segmentation/io/export.ts create mode 100644 src/segmentation/io/maskFileNaming.ts create mode 100644 src/segmentation/masks/labelValue.ts create mode 100644 src/segmentation/masks/overlap.ts create mode 100644 src/segmentation/masks/storage.ts create mode 100644 src/segmentation/masks/voxelAccess.ts create mode 100644 src/segmentation/model.ts create mode 100644 src/segmentation/rendering/__tests__/segmentDisplay.spec.ts create mode 100644 src/segmentation/rendering/__tests__/segmentRenderMask.spec.ts create mode 100644 src/segmentation/rendering/__tests__/segmentSliceVisibility.spec.ts create mode 100644 src/segmentation/rendering/display.ts create mode 100644 src/segmentation/rendering/projection.ts create mode 100644 src/segmentation/rendering/renderMask.ts create mode 100644 src/segmentation/segment.ts create mode 100644 src/segmentation/segmentReferences.ts create mode 100644 src/segmentation/segmentRegistry.ts diff --git a/package-lock.json b/package-lock.json index c50b7d12a..8c1c0594f 100644 --- a/package-lock.json +++ b/package-lock.json @@ -23,6 +23,7 @@ "@sentry/vue": "^10.27.0", "@thi.ng/api": "^8.12.9", "@thi.ng/rasterize": "^1.0.171", + "@types/color-name": "^2.0.0", "@types/cors": "^2.8.19", "@types/deep-equal": "^1.0.4", "@types/express": "^5.0.5", @@ -39,6 +40,7 @@ "@wdio/spec-reporter": "^9.20.0", "@wdio/static-server-service": "^9.20.0", "@wdio/visual-service": "^10.1.0", + "color-name": "^1.1.4", "comlink": "^4.4.2", "concurrently": "^10.0.4", "core-js": "3.47.0", @@ -5142,6 +5144,13 @@ "assertion-error": "^2.0.1" } }, + "node_modules/@types/color-name": { + "version": "2.0.0", + "resolved": "https://registry.npmjs.org/@types/color-name/-/color-name-2.0.0.tgz", + "integrity": "sha512-63mTjolMJv75upGaUbT6J3lRDWl6pETPQsaWni9w3dMArhNBpgtHkX8ISb9zLV3YYLPA/SMk8ZGALa3k9WY/aQ==", + "dev": true, + "license": "MIT" + }, "node_modules/@types/connect": { "version": "3.4.38", "resolved": "https://registry.npmjs.org/@types/connect/-/connect-3.4.38.tgz", diff --git a/package.json b/package.json index 71cfe5912..93fe9da79 100644 --- a/package.json +++ b/package.json @@ -49,6 +49,7 @@ "@sentry/vue": "^10.27.0", "@thi.ng/api": "^8.12.9", "@thi.ng/rasterize": "^1.0.171", + "@types/color-name": "^2.0.0", "@types/cors": "^2.8.19", "@types/deep-equal": "^1.0.4", "@types/express": "^5.0.5", @@ -65,6 +66,7 @@ "@wdio/spec-reporter": "^9.20.0", "@wdio/static-server-service": "^9.20.0", "@wdio/visual-service": "^10.1.0", + "color-name": "^1.1.4", "comlink": "^4.4.2", "concurrently": "^10.0.4", "core-js": "3.47.0", diff --git a/src/store/__tests__/segmentImportPath.spec.ts b/src/io/resample/__tests__/ensureSameSpace.spec.ts similarity index 75% rename from src/store/__tests__/segmentImportPath.spec.ts rename to src/io/resample/__tests__/ensureSameSpace.spec.ts index 043e2e18e..910b1d77d 100644 --- a/src/store/__tests__/segmentImportPath.spec.ts +++ b/src/io/resample/__tests__/ensureSameSpace.spec.ts @@ -1,11 +1,8 @@ import { afterEach, describe, expect, it, vi } from 'vitest'; -import { createPinia, setActivePinia } from 'pinia'; import { resolve } from 'node:path'; import { InterfaceTypes, runPipelineNode, type Image } from 'itk-wasm'; import vtkDataArray from '@kitware/vtk.js/Common/Core/DataArray'; import vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; -import { useSegmentGroupStore } from '@/src/store/segmentGroups'; -import { useImageCacheStore } from '@/src/store/image-cache'; import { ensureSameSpace } from '@/src/io/resample/resample'; import * as wasm from '@/src/io/resample/itkWasmUtils'; @@ -53,26 +50,16 @@ describe('labelmap import index alignment', () => { expect(runWasm).not.toHaveBeenCalled(); }); - it('resamples a reflected child onto the parent grid with label interpolation', async () => { - setActivePinia(createPinia()); + it('reorients a reflected label child onto the parent grid without resampling', async () => { const runWasm = useNodeResampling(); const parent = makeImage(); const child = makeImage(); child.setOrigin([0, 2, 0]); child.setDirection(1, 0, 0, 0, -1, 0, 0, 0, 1); child.getPointData().getScalars().setComponent(0, 0, 7); - const cache = useImageCacheStore(); - cache.addVTKImageData(parent, 'CTA', { id: 'parent' }); - cache.addVTKImageData(child, 'cta-head-neck-total.seg.nii.gz', { - id: 'child', - }); - - const store = useSegmentGroupStore(); - const [id] = await store.convertImageToLabelmap('child', 'parent'); - const imported = store.dataIndex[id]; + const imported = await ensureSameSpace(parent, child, true); - expect(runWasm).toHaveBeenCalledOnce(); - expect(runWasm.mock.calls[0][1]).toContain('--label'); + expect(runWasm).not.toHaveBeenCalled(); expect(imported.getDirection()).toEqual(parent.getDirection()); expect(imported.getOrigin()).toEqual(parent.getOrigin()); expect(Array.from(imported.getPointData().getScalars().getData())).toEqual([ diff --git a/src/io/resample/__tests__/reorientLabelImage.spec.ts b/src/io/resample/__tests__/reorientLabelImage.spec.ts new file mode 100644 index 000000000..16b90b43a --- /dev/null +++ b/src/io/resample/__tests__/reorientLabelImage.spec.ts @@ -0,0 +1,103 @@ +import { describe, expect, it } from 'vitest'; +import vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; +import vtkDataArray from '@kitware/vtk.js/Common/Core/DataArray'; +import { reorientLabelImage } from '../reorientLabelImage'; + +const image = (dimensions: [number, number, number]) => { + const result = vtkImageData.newInstance(); + result.setDimensions(dimensions); + result.getPointData().setScalars( + vtkDataArray.newInstance({ + values: Uint16Array.from( + { length: dimensions.reduce((a, b) => a * b, 1) }, + (_, i) => i + 1 + ), + }) + ); + return result; +}; + +describe('label-grid reorientation', () => { + it.each([0, 1, 2, 3, 4, 5, 6, 7])( + 'preserves every label under axis flips %i', + (flips) => { + const source = image([3, 4, 5]); + const target = image([3, 4, 5]); + const direction: [ + number, + number, + number, + number, + number, + number, + number, + number, + number, + ] = [1, 0, 0, 0, 1, 0, 0, 0, 1]; + const origin: [number, number, number] = [0, 0, 0]; + [3, 4, 5].forEach((size, axis) => { + if (flips & (1 << axis)) { + direction[axis * 4] = -1; + origin[axis] = size - 1; + } + }); + source.setDirection(direction); + source.setOrigin(origin); + const output = reorientLabelImage(target, source)!; + const values = output.getPointData().getScalars().getData(); + for (let k = 0; k < 5; k++) + for (let j = 0; j < 4; j++) + for (let i = 0; i < 3; i++) { + const x = flips & 1 ? 2 - i : i, + y = flips & 2 ? 3 - j : j, + z = flips & 4 ? 4 - k : k; + expect(values[i + 3 * (j + 4 * k)]).toBe(1 + x + 3 * (y + 4 * z)); + } + expect(output.getDirection()).toEqual(target.getDirection()); + expect(source.getPointData().getScalars().getData()[0]).toBe(1); + } + ); + + it('permutes unequal axes', () => { + const source = image([3, 4, 5]); + const target = image([4, 3, 5]); + target.setDirection([0, 1, 0, 1, 0, 0, 0, 0, 1]); + const output = reorientLabelImage(target, source)!; + expect(output.getExtent()).toEqual(target.getExtent()); + const values = output.getPointData().getScalars().getData(); + for (let k = 0; k < 5; k++) + for (let j = 0; j < 3; j++) + for (let i = 0; i < 4; i++) + expect(values[i + 4 * (j + 3 * k)]).toBe(1 + j + 3 * (i + 4 * k)); + }); + + it('defers nonzero extents to the general path', () => { + const source = image([3, 4, 5]); + source.setExtent(1, 3, 2, 5, 3, 7); + expect(reorientLabelImage(source, source)).toBeNull(); + }); + + it('accepts DICOM precision differences without changing labels', () => { + const source = image([3, 4, 5]); + const target = image([3, 4, 5]); + source.setOrigin([0.000004, 0.000004, 0.000012]); + expect([ + ...reorientLabelImage(target, source)! + .getPointData() + .getScalars() + .getData(), + ]).toEqual([...source.getPointData().getScalars().getData()]); + }); + + it('defers fractional shifts, different sampling, and cropping to interpolation', () => { + const source = image([3, 4, 5]); + const target = image([3, 4, 5]); + source.setOrigin([0.25, 0, 0]); + expect(reorientLabelImage(target, source)).toBeNull(); + source.setOrigin([0, 0, 0]); + source.setSpacing([0.5, 1, 1]); + expect(reorientLabelImage(target, source)).toBeNull(); + source.setSpacing([1, 1, 1]); + expect(reorientLabelImage(image([1, 4, 5]), source)).toBeNull(); + }); +}); diff --git a/src/io/resample/reorientLabelImage.ts b/src/io/resample/reorientLabelImage.ts new file mode 100644 index 000000000..bb9d85922 --- /dev/null +++ b/src/io/resample/reorientLabelImage.ts @@ -0,0 +1,82 @@ +import vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; +import vtkImageReslice from '@kitware/vtk.js/Imaging/Core/ImageReslice'; +import { InterpolationMode } from '@kitware/vtk.js/Imaging/Core/AbstractImageInterpolator/Constants'; +import { mat4, vec3 } from 'gl-matrix'; + +/** Reorder an equivalent voxel grid without interpolating label boundaries. */ +export function reorientLabelImage(target: vtkImageData, source: vtkImageData) { + // ImageReslice produces zero-based output extents. + if ( + [target, source].some((image) => + [0, 2, 4].some((axis) => image.getExtent()[axis] !== 0) + ) + ) + return null; + const matrix = mat4.multiply( + mat4.create(), + source.getWorldToIndex(), + target.getIndexToWorld() + ); + const tolerance = 1e-3; + const axes = [0, 1, 2].map((column) => { + const values = [0, 1, 2].map((row) => matrix[4 * column + row]); + const axis = values.findIndex((value) => Math.abs(value) > 0.5); + if ( + axis < 0 || + values.some((value, row) => + row === axis + ? Math.abs(Math.abs(value) - 1) > tolerance + : Math.abs(value) > tolerance + ) + ) + return -1; + return axis; + }); + if (axes.includes(-1) || new Set(axes).size !== 3) return null; + const sourceSize = source.getDimensions(); + if ( + target.getDimensions().some((size, axis) => size !== sourceSize[axes[axis]]) + ) + return null; + + // Check the entire extent so rounding error cannot accumulate into a shift + // at the far edge. Matching physical bounds alone does not imply equal grids. + const from = target.getExtent(); + const to = source.getExtent(); + for (let corner = 0; corner < 8; corner++) { + const point = vec3.fromValues( + from[corner & 1], + from[2 + ((corner >> 1) & 1)], + from[4 + ((corner >> 2) & 1)] + ); + vec3.transformMat4(point, point, matrix); + if ( + [0, 1, 2].some( + (axis) => + Math.min( + Math.abs(point[axis] - to[axis * 2]), + Math.abs(point[axis] - to[axis * 2 + 1]) + ) > tolerance + ) + ) + return null; + } + + if (matrix.every((value, index) => value === (index % 5 === 0 ? 1 : 0))) + return source; + + const filter = vtkImageReslice.newInstance(); + filter.setOutputOrigin(target.getOrigin()); + filter.setOutputSpacing(target.getSpacing()); + filter.setOutputDirection(target.getDirection()); + filter.setOutputExtent(target.getExtent()); + filter.setOutputDimensionality(3); + filter.setTransformInputSampling(false); + filter.setInterpolationMode(InterpolationMode.NEAREST); + try { + filter.setInputData(source); + return filter.getOutputData() as vtkImageData; + } finally { + filter.delete(); + } +} diff --git a/src/io/resample/resample.ts b/src/io/resample/resample.ts index 41b302ec2..2d96fddaf 100644 --- a/src/io/resample/resample.ts +++ b/src/io/resample/resample.ts @@ -4,6 +4,7 @@ import vtkITKHelper from '@kitware/vtk.js/Common/DataModel/ITKHelper'; import { compareImageIndexGrids } from '@/src/utils/imageSpace'; import { shallowCopyImageData } from '@/src/utils/vtk-helpers'; import { runWasm } from './itkWasmUtils'; +import { reorientLabelImage } from './reorientLabelImage'; export async function resample(fixed: Image, moving: Image, label = false) { const labelFlag = label ? ['--label'] : []; @@ -30,7 +31,12 @@ export async function ensureSameSpace( ) { // Callers own what they get back and may hand it to something that disposes // it, so never return the candidate itself. - if (compareImageIndexGrids(target, resampleCandidate)) { + if (label) { + const reoriented = reorientLabelImage(target, resampleCandidate); + if (reoriented === resampleCandidate) + return shallowCopyImageData(resampleCandidate); + if (reoriented) return reoriented; + } else if (compareImageIndexGrids(target, resampleCandidate)) { return shallowCopyImageData(resampleCandidate); } const itkImage = await resample( diff --git a/src/segmentation/__tests__/model.spec.ts b/src/segmentation/__tests__/model.spec.ts new file mode 100644 index 000000000..f2fe69850 --- /dev/null +++ b/src/segmentation/__tests__/model.spec.ts @@ -0,0 +1,297 @@ +import { describe, expect, it } from 'vitest'; +import type { RGBAColor } from '@kitware/vtk.js/types'; + +import { TOOL_COLORS } from '@/src/config'; +import { + emptyExtent, + isEmptyExtent, + markedExtent, +} from '@/src/segmentation/geometry'; +import { + cssColorToRGBA, + tryCssColorToRGBA, + rgbaToCssColor, +} from '@/src/segmentation/color'; +import type { + LabelmapBinding, + SegmentMask, + Segmentation, +} from '@/src/segmentation/model'; +import type { Extent3D } from '@/src/segmentation/geometry'; +import { resolveSegmentAppearance } from '@/src/segmentation/segment'; +import vtkLabelMap from '@/src/vtk/LabelMap'; + +describe('emptyExtent', () => { + it('is the pinned empty sentinel', () => { + expect(emptyExtent()).toEqual([0, -1, 0, -1, 0, -1]); + }); + + it('returns a fresh extent per call', () => { + const first = emptyExtent(); + first[1] = 10; + + expect(emptyExtent()).toEqual([0, -1, 0, -1, 0, -1]); + }); +}); + +describe('isEmptyExtent', () => { + it('accepts the empty sentinel', () => { + expect(isEmptyExtent(emptyExtent())).toBe(true); + }); + + it('rejects an extent covering a whole image', () => { + expect(isEmptyExtent([0, 9, 0, 19, 0, 29])).toBe(false); + }); + + it('rejects a single voxel extent', () => { + // vtk.js extents are inclusive, so min === max is one voxel, not empty. + expect(isEmptyExtent([4, 4, 5, 5, 6, 6])).toBe(false); + }); + + it.each([ + ['i', [5, 4, 0, 9, 0, 9] as Extent3D], + ['j', [0, 9, 5, 4, 0, 9] as Extent3D], + ['k', [0, 9, 0, 9, 5, 4] as Extent3D], + ])('is empty when the %s axis is inverted', (_axis, extent) => { + expect(isEmptyExtent(extent)).toBe(true); + }); +}); + +describe('cssColorToRGBA', () => { + it('parses a tool color hex string', () => { + expect(cssColorToRGBA('#58f24c')).toEqual([88, 242, 76, 255]); + }); + + it('parses an alpha channel when present', () => { + expect(cssColorToRGBA('#58f24c80')).toEqual([88, 242, 76, 128]); + }); + + it('parses the named color used by the vector tool label defaults', () => { + expect(cssColorToRGBA('red')).toEqual([255, 0, 0, 255]); + }); + + it.each([ + ['orange', [255, 165, 0, 255]], + ['rebeccapurple', [102, 51, 153, 255]], + ])('parses the CSS color keyword %s', (name, expected) => { + expect(cssColorToRGBA(name)).toEqual(expected); + }); +}); + +describe('rgbaToCssColor', () => { + it.each([ + [[88, 242, 76, 255] as RGBAColor], + [[0, 0, 0, 255] as RGBAColor], + [[214, 0, 0, 128] as RGBAColor], + ])('round trips %j', (rgba) => { + expect(cssColorToRGBA(rgbaToCssColor(rgba))).toEqual(rgba); + }); + + it('emits a css color literal', () => { + expect(rgbaToCssColor([88, 242, 76, 255])).toMatch( + /^(#[0-9a-fA-F]{6,8}|rgba?\(.+\))$/ + ); + }); +}); + +describe('color conversion of existing label colors', () => { + it.each(TOOL_COLORS)('round trips %s', (css) => { + const rgba = cssColorToRGBA(css); + + expect(cssColorToRGBA(rgbaToCssColor(rgba))).toEqual(rgba); + }); +}); + +describe('mask model', () => { + it('holds no labelmap binding until voxels are allocated', () => { + const segment: SegmentMask = { + id: 'segment-1', + segmentId: 'segment-1', + representations: {}, + }; + + expect(segment.representations.labelmap).toBeUndefined(); + }); + + it('binds a segment to its own voxels and the box they cover', () => { + const image = vtkLabelMap.newInstance(); + const binding: LabelmapBinding = { + image, + extent: [0, 9, 0, 19, 0, 29], + name: 'Tumor', + }; + const segment: SegmentMask = { + id: 'segment-1', + segmentId: 'segment-1', + representations: { labelmap: binding }, + }; + + expect(segment.representations.labelmap?.image).toBe(image); + expect(isEmptyExtent(binding.extent)).toBe(false); + }); + + it('keeps record order separate from the records themselves', () => { + const makeMask = (id: string, segmentId: string): SegmentMask => ({ + id, + segmentId, + representations: {}, + }); + const segmentation: Segmentation = { + id: 'segmentation-1', + name: 'Segmentation', + parentImageId: 'image-1', + masks: { + 'segment-1': makeMask('segment-1', 'segment-1'), + 'segment-2': makeMask('segment-2', 'segment-2'), + }, + order: ['segment-2', 'segment-1'], + fillOpacity: 1, + outlineOpacity: 1, + outlineThickness: 2, + }; + + expect(segmentation.order).toEqual(['segment-2', 'segment-1']); + expect(Object.keys(segmentation.masks).sort()).toEqual([ + 'segment-1', + 'segment-2', + ]); + }); +}); + +describe('the appearance resolver', () => { + it('fills the app defaults for what a type leaves unset', () => { + const resolved = resolveSegmentAppearance({ + id: 'segment-1', + name: 'Tumor', + color: [255, 0, 0, 255], + visible: true, + locked: false, + }); + + expect(resolved).toMatchObject({ + name: 'Tumor', + cssColor: '#ff0000', + fillOpacity: 1, + outlineOpacity: 1, + }); + }); + + it('keeps what a type does state', () => { + const resolved = resolveSegmentAppearance({ + id: 'segment-1', + name: 'Tumor', + color: [255, 0, 0, 255], + visible: false, + locked: true, + fillOpacity: 0.25, + strokeWidth: 3, + }); + + expect(resolved.fillOpacity).toBe(0.25); + expect(resolved.strokeWidth).toBe(3); + expect(resolved.visible).toBe(false); + expect(resolved.locked).toBe(true); + }); + + it('answers for a type that is gone', () => { + const resolved = resolveSegmentAppearance(undefined); + + expect(resolved.name).toBe(''); + expect(resolved.fillOpacity).toBe(1); + expect(resolved.visible).toBe(true); + expect(resolved.locked).toBe(false); + }); + + describe('CSS colors', () => { + // Functional notation is not supported. A config using it is told so at the + // boundary that reads the file, rather than resolving to a plausible black. + it.each([ + 'rgb(0, 255, 0)', + 'rgb(0 255 0)', + 'rgba(255, 0, 0, 0.5)', + 'hsl(120, 100%, 50%)', + ])('does not parse %s', (css) => { + expect(tryCssColorToRGBA(css)).toBeUndefined(); + }); + + it('treats transparent as fully transparent, not black', () => { + expect(cssColorToRGBA('transparent')).toEqual([0, 0, 0, 0]); + }); + + it('reports unparseable input rather than silently blackening it', () => { + expect(tryCssColorToRGBA('not-a-color')).toBeUndefined(); + expect(cssColorToRGBA('not-a-color')).toEqual([0, 0, 0, 255]); + }); + + // A config.json label color and a 6.x state file both reach here unvalidated, + // and the migration that calls this is not wrapped in a try. + it.each(['constructor', '__proto__', 'toString', 'valueOf'])( + 'treats the inherited property name %s as unparseable', + (name) => { + expect(tryCssColorToRGBA(name)).toBeUndefined(); + expect(cssColorToRGBA(name)).toEqual([0, 0, 0, 255]); + } + ); + }); +}); + +// A binding's extent is the allocation: paint pads it on growth and an erase +// never shrinks it, so the box a segment occupies has to be read off the voxels. +describe('markedExtent', () => { + const VALUE = 3; + + // A mask bounded to `extent`, holding VALUE at each parent index in `marks`. + const scalarsOf = ( + extent: Extent3D, + marks: Array<[number, number, number]> + ) => { + const si = extent[1] - extent[0] + 1; + const sj = extent[3] - extent[2] + 1; + const sk = extent[5] - extent[4] + 1; + const scalars = new Uint8Array(si * sj * sk); + marks.forEach(([i, j, k]) => { + scalars[ + i - extent[0] + (j - extent[2]) * si + (k - extent[4]) * si * sj + ] = VALUE; + }); + return scalars; + }; + + it('bounds the marked voxels, not the allocation they sit in', () => { + const extent: Extent3D = [0, 5, 0, 5, 0, 5]; + const scalars = scalarsOf(extent, [ + [1, 2, 3], + [4, 2, 3], + [2, 5, 1], + ]); + + expect(markedExtent(scalars, extent, VALUE)).toEqual([1, 4, 2, 5, 1, 3]); + }); + + it('reads a single voxel as its own box', () => { + const extent: Extent3D = [2, 4, 2, 4, 2, 4]; + const scalars = scalarsOf(extent, [[3, 3, 3]]); + + expect(markedExtent(scalars, extent, VALUE)).toEqual([3, 3, 3, 3, 3, 3]); + }); + + it('reports an empty box for a mask holding nothing of that value', () => { + const extent: Extent3D = [0, 3, 0, 3, 0, 3]; + const scalars = scalarsOf(extent, [[1, 1, 1]]); + + expect(isEmptyExtent(markedExtent(scalars, extent, VALUE + 1))).toBe(true); + expect(isEmptyExtent(markedExtent(new Uint8Array(64), extent, VALUE))).toBe( + true + ); + }); + + it('reads a mask whose own origin is not the image origin', () => { + const extent: Extent3D = [4, 6, 7, 9, 1, 2]; + const scalars = scalarsOf(extent, [ + [5, 8, 1], + [6, 9, 2], + ]); + + expect(markedExtent(scalars, extent, VALUE)).toEqual([5, 6, 8, 9, 1, 2]); + }); +}); diff --git a/src/segmentation/__tests__/reframeMaskScalars.spec.ts b/src/segmentation/__tests__/reframeMaskScalars.spec.ts new file mode 100644 index 000000000..ba4554f60 --- /dev/null +++ b/src/segmentation/__tests__/reframeMaskScalars.spec.ts @@ -0,0 +1,44 @@ +import { describe, expect, it } from 'vitest'; +import { reframeMaskScalars } from '@/src/segmentation/masks/storage'; +import type { Extent3D } from '@/src/segmentation/geometry'; + +describe('reframing into a reusable mask buffer', () => { + it('overwrites old labels and padding when the source moves or is erased', () => { + const to: Extent3D = [4, 6, 5, 7, 6, 8]; + const output = new Uint8Array(27).fill(9); + const source = new Uint8Array([1]); + expect(reframeMaskScalars(source, [5, 5, 6, 6, 7, 7], to, output)).toBe( + output + ); + expect([...output]).toEqual( + Array.from({ length: 27 }, (_, i) => (i === 13 ? 1 : 0)) + ); + reframeMaskScalars(source, [6, 6, 7, 7, 8, 8], to, output); + expect([...output]).toEqual( + Array.from({ length: 27 }, (_, i) => (i === 26 ? 1 : 0)) + ); + source[0] = 0; + reframeMaskScalars(source, [6, 6, 7, 7, 8, 8], to, output); + expect([...output]).toEqual(Array(27).fill(0)); + }); + + it('clips source rows at the destination and clears disjoint copies', () => { + const from: Extent3D = [2, 5, 3, 4, 4, 4]; + const to: Extent3D = [3, 4, 4, 5, 4, 4]; + const source = new Uint8Array([1, 2, 3, 4, 5, 6, 7, 8]); + const output = new Uint8Array(4).fill(9); + reframeMaskScalars(source, from, to, output); + expect([...output]).toEqual([6, 7, 0, 0]); + expect([...source]).toEqual([1, 2, 3, 4, 5, 6, 7, 8]); + reframeMaskScalars(source, from, [3, 4, 9, 10, 4, 4], output); + expect([...output]).toEqual([0, 0, 0, 0]); + }); + + it('rejects a destination buffer of the wrong size before changing it', () => { + const output = new Uint8Array([9]); + expect(() => + reframeMaskScalars([1], [0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0], output) + ).toThrow('Mask output size mismatch'); + expect([...output]).toEqual([9]); + }); +}); diff --git a/src/segmentation/color.ts b/src/segmentation/color.ts new file mode 100644 index 000000000..832a59fb5 --- /dev/null +++ b/src/segmentation/color.ts @@ -0,0 +1,41 @@ +import type { RGBAColor } from '@kitware/vtk.js/types'; +import colorNames from 'color-name'; +import { hexaToRGBA, rgbaToHexa } from '@/src/utils/color'; + +// A Map, not the package's plain object: prototype keys such as 'constructor' +// must not be mistaken for colors. +const NAMED_COLORS = new Map(Object.entries(colorNames)); + +const HEX_COLOR = /^#?([0-9a-f]{3,4}|[0-9a-f]{6}|[0-9a-f]{8})$/; + +const expandShorthandHex = (hex: string) => + hex.length <= 4 + ? hex + .split('') + .map((digit) => digit.repeat(2)) + .join('') + : hex; + +/** Parses `transparent`, the CSS colour keywords, and hex. */ +export function tryCssColorToRGBA(css: string): RGBAColor | undefined { + const value = css.trim().toLowerCase(); + if (value === 'transparent') return [0, 0, 0, 0]; + + const named = NAMED_COLORS.get(value); + if (named) return [...named, 255] as RGBAColor; + + const hex = HEX_COLOR.exec(value)?.[1]; + return hex ? hexaToRGBA(expandShorthandHex(hex)) : undefined; +} + +/** Falls back to opaque black for unparseable label colors. */ +export function cssColorToRGBA(css: string): RGBAColor { + return tryCssColorToRGBA(css) ?? [0, 0, 0, 255]; +} + +// Opaque colors keep the 6-digit form label colors are written in, so a color +// that round trips through a segment comes back byte-identical. +export function rgbaToCssColor(rgba: RGBAColor) { + const hexa = rgbaToHexa(rgba); + return rgba[3] === 255 ? hexa.slice(0, 7) : hexa; +} diff --git a/src/segmentation/editing/algorithms/__tests__/gaussianSmooth.spec.ts b/src/segmentation/editing/algorithms/__tests__/gaussianSmooth.spec.ts new file mode 100644 index 000000000..4a809e805 --- /dev/null +++ b/src/segmentation/editing/algorithms/__tests__/gaussianSmooth.spec.ts @@ -0,0 +1,122 @@ +import { describe, it, expect } from 'vitest'; +import { fullExtent, type Extent3D } from '@/src/segmentation/geometry'; +import { gaussianSmoothLabelMapWorker } from '@/src/segmentation/editing/algorithms/gaussianSmooth.worker'; + +const LABEL = 3; +type Dims = [number, number, number]; + +/** Reconstructs the whole parent, including growth outside the input mask. */ +function smooth( + data: Uint8Array, + dimensions: Dims, + maskExtent: Extent3D = fullExtent(dimensions), + parentDimensions: Dims = dimensions +) { + const { scalars, extent } = gaussianSmoothLabelMapWorker({ + data, + dimensions, + spacing: [1, 1, 1], + maskExtent, + parentDimensions, + params: { sigma: 1, label: LABEL }, + }); + const [pi, pj, pk] = parentDimensions; + const parent = new Uint8Array(pi * pj * pk); + let offset = 0; + for (let k = extent[4]; k <= extent[5]; k += 1) { + for (let j = extent[2]; j <= extent[3]; j += 1) { + for (let i = extent[0]; i <= extent[1]; i += 1) { + parent[i + pi * (j + pj * k)] = scalars[offset++]; + } + } + } + return parent; +} + +function boxInParent(extent: Extent3D, dimensions: Dims) { + const [pi, pj, pk] = dimensions; + const parent = new Uint8Array(pi * pj * pk); + for (let k = extent[4]; k <= extent[5]; k += 1) { + for (let j = extent[2]; j <= extent[3]; j += 1) { + for (let i = extent[0]; i <= extent[1]; i += 1) { + parent[i + pi * (j + pj * k)] = LABEL; + } + } + } + return parent; +} + +describe('gaussianSmoothLabelMapWorker', () => { + it('smooths an isolated voxel away even when the mask is that one voxel', () => { + const smoothed = smooth( + new Uint8Array([LABEL]), + [1, 1, 1], + [1, 1, 1, 1, 1, 1], + [3, 3, 3] + ); + expect(smoothed.every((value) => value === 0)).toBe(true); + }); + + it.each([ + [2, 4, 2, 4, 2, 4], + [1, 3, 1, 3, 1, 3], + [0, 2, 2, 4, 2, 4], + [1, 4, 1, 4, 1, 4], + [4, 7, 4, 7, 4, 7], + ] as Extent3D[])( + 'preserves the full-parent output for a tight mask at [%i, %i, %i, %i, %i, %i]', + (...extent) => { + const parentDimensions: Dims = [9, 9, 9]; + const dimensions: Dims = [ + extent[1] - extent[0] + 1, + extent[3] - extent[2] + 1, + extent[5] - extent[4] + 1, + ]; + const data = new Uint8Array( + dimensions[0] * dimensions[1] * dimensions[2] + ).fill(LABEL); + const croppedResult = smooth(data, dimensions, extent, parentDimensions); + const parentResult = smooth( + boxInParent(extent, parentDimensions), + parentDimensions + ); + expect(croppedResult).toEqual(parentResult); + } + ); + + it('retains all 62 voxels when mirroring grows a box toward parent faces', () => { + const smoothed = smooth( + new Uint8Array(64).fill(LABEL), + [4, 4, 4], + [1, 4, 1, 4, 1, 4], + [9, 9, 9] + ); + expect(smoothed.filter((value) => value === LABEL)).toHaveLength(62); + expect(smoothed[0 + 2 * 9 + 2 * 81]).toBe(LABEL); + }); + + it('erodes a cropped cube at its corners and keeps its centre', () => { + const smoothed = smooth( + new Uint8Array(27).fill(LABEL), + [3, 3, 3], + [1, 3, 1, 3, 1, 3], + [5, 5, 5] + ); + expect(smoothed[1 + 1 * 5 + 1 * 25]).toBe(0); + expect(smoothed[2 + 2 * 5 + 2 * 25]).toBe(LABEL); + }); + + it('mirrors only at parent faces, not at the mask allocation', () => { + const data = new Uint8Array(27).fill(LABEL); + const againstFace = smooth(data, [3, 3, 3], [0, 2, 2, 4, 2, 4], [9, 9, 9]); + const awayFromFace = smooth(data, [3, 3, 3], [1, 3, 2, 4, 2, 4], [9, 9, 9]); + expect(againstFace[0 + 3 * 9 + 3 * 81]).toBe(LABEL); + expect(awayFromFace[0 + 3 * 9 + 3 * 81]).toBe(0); + }); + + it('leaves a buffer with none of the label alone', () => { + expect(Array.from(smooth(new Uint8Array([0, 1, 0, 1]), [4, 1, 1]))).toEqual( + [0, 1, 0, 1] + ); + }); +}); diff --git a/src/segmentation/editing/algorithms/__tests__/gaussianSmoothGolden.spec.ts b/src/segmentation/editing/algorithms/__tests__/gaussianSmoothGolden.spec.ts new file mode 100644 index 000000000..be65c6a3d --- /dev/null +++ b/src/segmentation/editing/algorithms/__tests__/gaussianSmoothGolden.spec.ts @@ -0,0 +1,237 @@ +import { describe, it, expect } from 'vitest'; +import { gaussianSmoothLabelMapWorker } from '@/src/segmentation/editing/algorithms/gaussianSmooth.worker'; + +// A byte-exact record of the filter's output. The rest of the smoothing suite +// asserts on shape properties, which a change in the arithmetic can satisfy +// while every voxel moves; this catches that. +// +// Each entry is one row of the volume, '1' where the label survives, so a +// failure names the rows that moved. To refresh one after a deliberate change +// in behaviour, print `asRows(smooth(...))`. + +const LABEL = 3; +const DIMENSIONS: [number, number, number] = [9, 8, 7]; + +const SIGMA_0_6 = [ + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '011100100', + '011101100', + '011011100', + '000111100', + '001111100', + '000000000', + '000000000', + '000000000', + '011101100', + '011111100', + '011111100', + '001111100', + '011111100', + '000000000', + '000000000', + '000000000', + '011011100', + '011111100', + '001111100', + '011111100', + '011111000', + '000000000', + '000000000', + '000000000', + '000111100', + '001111100', + '011111100', + '011111000', + '011110000', + '000000000', + '000000000', + '000000000', + '001111100', + '011111100', + '011111000', + '011110000', + '011100100', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', +]; + +const SIGMA_1_0 = [ + '000000000', + '000000000', + '000000000', + '000001000', + '000011000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '001111000', + '001111100', + '000111000', + '000010000', + '000000000', + '000000000', + '000000000', + '001111000', + '011111100', + '011111100', + '001111100', + '001111000', + '000000000', + '000000000', + '000001000', + '001111100', + '011111100', + '011111100', + '011111100', + '001111000', + '000000000', + '000000000', + '000011000', + '000111000', + '001111100', + '011111100', + '011111000', + '001110000', + '000000000', + '000000000', + '000000000', + '000111000', + '001111000', + '011111000', + '011110000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000110000', + '001110000', + '001100000', + '000000000', + '000000000', + '000000000', +]; + +const ANISOTROPIC = [ + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '001100000', + '011110000', + '111111100', + '111111100', + '001111100', + '000111000', + '000000000', + '000000000', + '001000000', + '011111100', + '111111100', + '011111100', + '001111100', + '001111000', + '000000000', + '000000000', + '000001000', + '000111100', + '001111100', + '011111100', + '011111100', + '011111000', + '000000000', + '000000000', + '000011000', + '000111100', + '001111100', + '011111100', + '111111000', + '011110000', + '000000000', + '000000000', + '000111000', + '001111100', + '011111100', + '111111000', + '111111000', + '011100000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', + '000000000', +]; +/** An asymmetric blob with a pitted interior and two lone corner voxels. */ +const blob = () => { + const [di, dj, dk] = DIMENSIONS; + const data = new Uint8Array(di * dj * dk); + const mark = (i: number, j: number, k: number) => { + data[i + j * di + k * di * dj] = LABEL; + }; + for (let k = 1; k <= 5; k += 1) { + for (let j = 1; j <= 5; j += 1) { + for (let i = 1; i <= 6; i += 1) { + if ((i + j + k) % 7 !== 0) mark(i, j, k); + } + } + } + mark(8, 0, 0); + mark(0, 7, 6); + return data; +}; + +const smooth = (sigma: number, spacing: [number, number, number]) => + gaussianSmoothLabelMapWorker({ + data: blob(), + dimensions: DIMENSIONS, + spacing, + maskExtent: [0, 8, 0, 7, 0, 6], + parentDimensions: DIMENSIONS, + params: { sigma, label: LABEL }, + }).scalars; + +const asRows = (output: ArrayLike) => { + const bits = Array.from(output, (value) => (value ? '1' : '0')).join(''); + return bits.match(new RegExp(`.{${DIMENSIONS[0]}}`, 'g')); +}; + +describe('gaussian smooth golden output', () => { + it('matches byte for byte at a sigma below one voxel', () => { + expect(asRows(smooth(0.6, [1, 1, 1]))).toEqual(SIGMA_0_6); + }); + + it('matches byte for byte at a one-voxel sigma', () => { + expect(asRows(smooth(1.0, [1, 1, 1]))).toEqual(SIGMA_1_0); + }); + + it('matches byte for byte through anisotropic spacing', () => { + expect(asRows(smooth(1.0, [1, 1, 3]))).toEqual(ANISOTROPIC); + }); +}); diff --git a/src/segmentation/editing/algorithms/fillHoles.ts b/src/segmentation/editing/algorithms/fillHoles.ts new file mode 100644 index 000000000..f8f24deba --- /dev/null +++ b/src/segmentation/editing/algorithms/fillHoles.ts @@ -0,0 +1,165 @@ +import { TypedArray } from '@kitware/vtk.js/types'; + +// 4-connected neighbor offsets, shared so the flood-fill loops never allocate +// a neighbor array per visited voxel. +const NEIGHBOR_DU = [-1, 1, 0, 0]; +const NEIGHBOR_DV = [0, 0, -1, 1]; + +type MaskData = TypedArray | number[]; + +export type FillHolesOptions = { + // Flat label-map scalar array, indexed as i + j*dimI + k*dimI*dimJ. + data: MaskData; + // Label-map IJK dimensions [dimI, dimJ, dimK]. + dimensions: [number, number, number]; + // IJK axis perpendicular to the fill plane (the slice axis). + axis: 0 | 1 | 2; + // When set, only this slice index along `axis` is processed. + // When omitted, every slice along `axis` is processed. + sliceIndex?: number; + // The one label treated as foreground, and the one enclosed background is + // filled with. Storage is a mask per segment, so a fill is always one + // segment's own; running over several segments is one call each. + label: number; +}; + +type Plane = { + uDim: number; + vDim: number; + offset: (u: number, v: number) => number; +}; + +// `visited` and `stack` are reused across slices, so the whole run allocates +// them once. +// `visited`: 0 = unvisited, 1 = outside (border-connected), 2 = hole. +type Flood = { + plane: Plane; + out: MaskData; + label: number; + visited: Uint8Array; + stack: number[]; +}; + +const inPlane = (plane: Plane, u: number, v: number) => + u >= 0 && u < plane.uDim && v >= 0 && v < plane.vDim; + +// Drain `stack`, expanding the region into unvisited non-foreground neighbors +// (each marked with `mark`). `collect`, when given, receives the flat offset of +// every region cell. +function drain(flood: Flood, mark: number, collect?: number[]) { + const { plane, out, label, visited, stack } = flood; + while (stack.length) { + const p = stack.pop()!; + const u = p % plane.uDim; + const v = (p - u) / plane.uDim; + if (collect) collect.push(plane.offset(u, v)); + for (let n = 0; n < 4; n++) { + const nu = u + NEIGHBOR_DU[n]; + const nv = v + NEIGHBOR_DV[n]; + if (!inPlane(plane, nu, nv)) continue; + const np = nu + nv * plane.uDim; + if (out[plane.offset(nu, nv)] !== label && visited[np] === 0) { + visited[np] = mark; + stack.push(np); + } + } + } +} + +// Flood the non-foreground cells reachable from the slice border, marking them +// "outside". Whatever it does not reach is enclosed. +function markOutside(flood: Flood) { + const { plane, out, label, visited, stack } = flood; + const seed = (u: number, v: number) => { + const p = u + v * plane.uDim; + if (visited[p] === 0 && out[plane.offset(u, v)] !== label) { + visited[p] = 1; + stack.push(p); + } + }; + for (let u = 0; u < plane.uDim; u++) { + seed(u, 0); + seed(u, plane.vDim - 1); + } + for (let v = 0; v < plane.vDim; v++) { + seed(0, v); + seed(plane.uDim - 1, v); + } + drain(flood, 1); +} + +// Only fill background. A voxel another segment holds is not this segment's to +// take here; the write path decides that on confirm. +function fillBackground(out: MaskData, cells: number[], label: number) { + for (let c = 0; c < cells.length; c++) { + if (out[cells[c]] === 0) { + out[cells[c]] = label; + } + } +} + +// Any non-foreground cell not marked "outside" is part of a hole. Group each +// hole into a connected component and fill it. +function fillEnclosed(flood: Flood) { + const { plane, out, label, visited, stack } = flood; + for (let v = 0; v < plane.vDim; v++) { + for (let u = 0; u < plane.uDim; u++) { + const p = u + v * plane.uDim; + if (visited[p] !== 0 || out[plane.offset(u, v)] === label) continue; + + const holeCells: number[] = []; + visited[p] = 2; + stack.push(p); + drain(flood, 2, holeCells); + fillBackground(out, holeCells, label); + } + } +} + +// Fills enclosed background regions ("holes") on 2D slices of a label map. +// A hole is background that does not connect to the slice border. Only +// background (0) voxels are filled. Returns a copy of `data`; the input is left +// untouched. +export function fillHoles(opts: FillHolesOptions) { + const { data, dimensions, axis, sliceIndex, label } = opts; + const out = data.slice(); + + const strides = [1, dimensions[0], dimensions[0] * dimensions[1]]; + const sliceStride = strides[axis]; + const sliceCount = dimensions[axis]; + + // The two in-plane axes (everything that isn't the slice axis). + const [uAxis, vAxis] = [0, 1, 2].filter((a) => a !== axis); + const uDim = dimensions[uAxis]; + const vDim = dimensions[vAxis]; + const uStride = strides[uAxis]; + const vStride = strides[vAxis]; + + const visited = new Uint8Array(uDim * vDim); + const stack: number[] = []; + + const firstSlice = sliceIndex ?? 0; + const lastSlice = sliceIndex ?? sliceCount - 1; + + for (let slice = firstSlice; slice <= lastSlice; slice++) { + const base = slice * sliceStride; + visited.fill(0); + + const flood: Flood = { + plane: { + uDim, + vDim, + offset: (u, v) => base + u * uStride + v * vStride, + }, + out, + label, + visited, + stack, + }; + + markOutside(flood); + fillEnclosed(flood); + } + + return out; +} diff --git a/src/segmentation/editing/algorithms/fillHoles.worker.ts b/src/segmentation/editing/algorithms/fillHoles.worker.ts new file mode 100644 index 000000000..bf548bb4b --- /dev/null +++ b/src/segmentation/editing/algorithms/fillHoles.worker.ts @@ -0,0 +1,17 @@ +import * as Comlink from 'comlink'; +import { + fillHoles, + FillHolesOptions, +} from '@/src/segmentation/editing/algorithms/fillHoles'; + +// Runs the pure flood-fill off the main thread so whole-volume fills do not +// freeze the UI, mirroring gaussianSmooth.worker.ts. +export function fillHolesWorker(input: FillHolesOptions) { + return fillHoles(input); +} + +const workerApi = { + fillHolesWorker, +}; + +Comlink.expose(workerApi); diff --git a/src/segmentation/editing/algorithms/gaussianSmooth.worker.ts b/src/segmentation/editing/algorithms/gaussianSmooth.worker.ts new file mode 100644 index 000000000..54a130925 --- /dev/null +++ b/src/segmentation/editing/algorithms/gaussianSmooth.worker.ts @@ -0,0 +1,384 @@ +import * as Comlink from 'comlink'; +import { TypedArray } from '@kitware/vtk.js/types'; +import { createTypedArrayLike } from '@/src/utils'; +import { + extentSize, + extentUnion, + type Extent3D, +} from '@/src/segmentation/geometry'; + +export interface GaussianSmoothParams { + sigma: number; + label: number; +} + +export interface GaussianSmoothInput { + data: TypedArray | number[]; + dimensions: number[]; + spacing: [number, number, number]; + maskExtent: [number, number, number, number, number, number]; + parentDimensions: [number, number, number]; + params: GaussianSmoothParams; +} + +function generateGaussianKernel(sigma: number, radiusFactor = 1.5) { + const radius = Math.ceil(sigma * radiusFactor); + const size = 2 * radius + 1; + const kernel = new Float32Array(size); + const center = radius; + let sum = 0; + + for (let i = 0; i < size; i++) { + const x = i - center; + // VTK formula: exp(-(x * x) / (std * std * 2.0)) + const value = Math.exp(-(x * x) / (sigma * sigma * 2.0)); + kernel[i] = value; + sum += value; + } + + // Normalize kernel + for (let i = 0; i < size; i++) { + kernel[i] /= sum; + } + + return kernel; +} + +// Helper for robust boundary handling (mirroring) +function mirrorCoord(sampleCoord: number, axisDim: number) { + let finalCoord = sampleCoord; + if (sampleCoord < 0) { + finalCoord = -sampleCoord; // Reflect + } else if (sampleCoord >= axisDim) { + finalCoord = 2 * axisDim - sampleCoord - 2; // Reflect + } + // Clamp to ensure it's within bounds, useful if kernel is very large + return Math.max(0, Math.min(axisDim - 1, finalCoord)); +} + +/** + * What stays fixed for one axis pass: the two buffers, the kernel, and how a + * line along the convolved axis is addressed. Built once per pass, so walking + * a line allocates nothing. + */ +interface AxisPass { + inputData: TypedArray | number[]; + outputData: TypedArray | number[]; + kernel: Float32Array; + kernelCenter: number; + // Voxels along the convolved axis, and the index step between them. + count: number; + stride: number; +} + +// Convolves the one line of voxels that starts at `lineStart` and runs along +// the pass's axis. Every axis reduces to this, because a line differs only in +// where it starts, how long it is, and how far apart its voxels sit. +function convolveLine(pass: AxisPass, lineStart: number) { + const { inputData, outputData, kernel, kernelCenter, count, stride } = pass; + const kernelSize = kernel.length; + + for (let i = 0; i < count; i++) { + let sum = 0; + for (let k = 0; k < kernelSize; k++) { + const sample = mirrorCoord(i + k - kernelCenter, count); + sum += inputData[sample * stride + lineStart] * kernel[k]; + } + + outputData[i * stride + lineStart] = sum; + } +} + +function convolve1D( + inputData: TypedArray | number[], + outputData: TypedArray | number[], + kernel: Float32Array, + volume: { dimensions: number[]; axis: 0 | 1 | 2 } +) { + const { dimensions, axis } = volume; + const [dimX, dimY] = dimensions; + const strides = [1, dimX, dimX * dimY]; + const pass: AxisPass = { + inputData, + outputData, + kernel, + kernelCenter: Math.floor(kernel.length / 2), + count: dimensions[axis], + stride: strides[axis], + }; + + // The two axes the pass does not walk, the widest-striding one outermost: + // z, then y, then x. That is the loop order each axis wants for cache + // efficiency, so convolving along X visits z, y, x, along Y visits z, x, y, + // and along Z visits y, x, z. + const outer = axis === 2 ? 1 : 2; + const inner = axis === 0 ? 1 : 0; + + for (let o = 0; o < dimensions[outer]; o++) { + const outerOffset = o * strides[outer]; + for (let i = 0; i < dimensions[inner]; i++) { + convolveLine(pass, outerOffset + i * strides[inner]); + } + } +} + +function gaussianFilter3D( + inputData: TypedArray | number[], + dimensions: number[], + sigmaPixels: [number, number, number], + radiusFactor = 1.5 +) { + const totalSize = dimensions[0] * dimensions[1] * dimensions[2]; + const kernelX = generateGaussianKernel(sigmaPixels[0], radiusFactor); + const kernelY = generateGaussianKernel(sigmaPixels[1], radiusFactor); + const kernelZ = generateGaussianKernel(sigmaPixels[2], radiusFactor); + const temp = new Float32Array(totalSize); + const output = new Float32Array(totalSize); + + convolve1D(inputData, output, kernelX, { dimensions, axis: 0 }); + convolve1D(output, temp, kernelY, { dimensions, axis: 1 }); + convolve1D(temp, output, kernelZ, { dimensions, axis: 2 }); + + return output; +} + +// What a bounding-box scan holds still: the voxels being read, the bounds +// being widened, and the row addressing. Built once, so scanning a row +// allocates nothing. +interface RowScan { + data: TypedArray | number[]; + bounds: number[]; + dimX: number; + sliceSize: number; + label: number; +} + +// Its own function so that the per-voxel test sits two blocks deep rather than +// four. +function growBoundsOverRow(scan: RowScan, y: number, z: number) { + const { data, bounds, dimX, sliceSize, label } = scan; + const rowStart = y * dimX + z * sliceSize; + + for (let x = 0; x < dimX; x++) { + if (data[rowStart + x] !== label) continue; + bounds[0] = Math.min(bounds[0], x); + bounds[1] = Math.max(bounds[1], x); + bounds[2] = Math.min(bounds[2], y); + bounds[3] = Math.max(bounds[3], y); + bounds[4] = Math.min(bounds[4], z); + bounds[5] = Math.max(bounds[5], z); + } +} + +function calculateBoundingBox( + data: TypedArray | number[], + dimensions: number[], + label: number +) { + const [dimX, dimY, dimZ] = dimensions; + const bounds = [dimX, -1, dimY, -1, dimZ, -1]; + const scan: RowScan = { + data, + bounds, + dimX, + sliceSize: dimX * dimY, + label, + }; + + for (let z = 0; z < dimZ; z++) { + for (let y = 0; y < dimY; y++) { + growBoundsOverRow(scan, y, z); + } + } + + if (bounds[1] === -1) return null; + + return bounds; +} + +function expandBoundingBox({ + bounds, + maskExtent, + parentDimensions, + sigmaPixels, + radiusFactor = 1.5, +}: { + bounds: number[]; + maskExtent: GaussianSmoothInput['maskExtent']; + parentDimensions: GaussianSmoothInput['parentDimensions']; + sigmaPixels: [number, number, number]; + radiusFactor?: number; +}) { + return sigmaPixels.flatMap((sigma, axis) => { + const padding = Math.ceil(sigma * radiusFactor); + // The parent-image faces, stated in mask coordinates. Ending the + // convolution volume there keeps the established mirrored boundary + // wherever the mask sits, so the result does not depend on how much of + // the parent the mask happens to be allocated over. Crop faces are not + // clamped: outside the buffer reads as background. + const parentLow = -maskExtent[axis * 2]; + const parentHigh = parentDimensions[axis] - 1 - maskExtent[axis * 2]; + return [ + Math.max(parentLow, bounds[axis * 2] - padding), + Math.min(parentHigh, bounds[axis * 2 + 1] + padding), + ]; + }); +} + +/** + * Visits every voxel of `bounds` that the volume actually holds, giving each + * its offset in the volume and its offset in the padded sub-volume. The + * padding ring outside the volume is skipped by clipping the loops, not tested + * per voxel. + */ +function forEachClippedVoxel( + dimensions: number[], + bounds: number[], + visit: (origIndex: number, subIndex: number) => void +) { + const [dimX, dimY, dimZ] = dimensions; + const [minX, maxX, minY, maxY, minZ, maxZ] = bounds; + const subDimX = maxX - minX + 1; + const subDimY = maxY - minY + 1; + const lastX = Math.min(maxX, dimX - 1); + + for (let z = Math.max(minZ, 0); z <= Math.min(maxZ, dimZ - 1); z += 1) { + for (let y = Math.max(minY, 0); y <= Math.min(maxY, dimY - 1); y += 1) { + const rowOrig = y * dimX + z * dimX * dimY; + const rowSub = (y - minY) * subDimX + (z - minZ) * subDimX * subDimY; + for (let x = Math.max(minX, 0); x <= lastX; x += 1) { + visit(x + rowOrig, x - minX + rowSub); + } + } + } +} + +/** + * The label's own binary mask over `bounds`, which is the only thing the + * filter reads. Built in one pass rather than copying the labels out and + * thresholding them afterwards: the copy is a second volume-sized Float32 + * array, live at the same time as this one. + */ +function extractSubMask( + data: TypedArray | number[], + dimensions: number[], + bounds: number[], + label: number +) { + const [minX, maxX, minY, maxY, minZ, maxZ] = bounds; + const subDims = [maxX - minX + 1, maxY - minY + 1, maxZ - minZ + 1]; + // Zero filled, so everything outside the buffer stays background. + const subMask = new Float32Array(subDims[0] * subDims[1] * subDims[2]); + + forEachClippedVoxel(dimensions, bounds, (origIndex, subIndex) => { + subMask[subIndex] = data[origIndex] === label ? 255.0 : 0.0; + }); + + return { subMask, subDims }; +} + +// Output storage includes the padding ring: mirroring at a parent face can +// turn on voxels beyond the input mask's allocation. +function copySubVolumeBack( + subData: Float32Array, + originalData: TypedArray | number[], + region: { dimensions: number[]; bounds: number[] }, + label: number +) { + forEachClippedVoxel( + region.dimensions, + region.bounds, + (origIndex, subIndex) => { + const origLabel = originalData[origIndex]; + if (origLabel === label || origLabel === 0) { + originalData[origIndex] = subData[subIndex] > 127.5 ? label : 0; + } + } + ); +} + +export function gaussianSmoothLabelMapWorker(input: GaussianSmoothInput) { + const { + data: originalData, + dimensions, + spacing, + maskExtent, + parentDimensions, + params, + } = input; + const { sigma, label } = params; + + if (sigma <= 0) { + throw new Error('Sigma must be positive'); + } + + const sigmaPixels: [number, number, number] = [ + sigma / spacing[0], + sigma / spacing[1], + sigma / spacing[2], + ]; + + // Absent when the label is nowhere in the mask, which is also the whole + // answer for a mask with nothing to smooth: it comes back as it went in. + const bounds = calculateBoundingBox(originalData, dimensions, label); + if (!bounds) { + const outputData = createTypedArrayLike(originalData, originalData.length); + for (let i = 0; i < originalData.length; i++) { + outputData[i] = originalData[i]; + } + return { scalars: outputData, extent: maskExtent }; + } + + const expandedBounds = expandBoundingBox({ + bounds, + maskExtent, + parentDimensions, + sigmaPixels, + }); + const { subMask, subDims } = extractSubMask( + originalData, + dimensions, + expandedBounds, + label + ); + + const smoothedSubMask = gaussianFilter3D(subMask, subDims, sigmaPixels, 1.5); + + const expandedExtent = expandedBounds.map( + (value, axis) => value + maskExtent[axis - (axis % 2)] + ) as Extent3D; + const extent = extentUnion(maskExtent, expandedExtent); + const outputDimensions = extentSize(extent); + const outputData = createTypedArrayLike( + originalData, + outputDimensions[0] * outputDimensions[1] * outputDimensions[2] + ); + const outputBoundsInInput = extent.map( + (value, axis) => value - maskExtent[axis - (axis % 2)] + ); + forEachClippedVoxel( + dimensions, + outputBoundsInInput, + (origIndex, outIndex) => { + outputData[outIndex] = originalData[origIndex]; + } + ); + const smoothedBoundsInOutput = expandedExtent.map( + (value, axis) => value - extent[axis - (axis % 2)] + ); + + copySubVolumeBack( + smoothedSubMask, + outputData, + { dimensions: outputDimensions, bounds: smoothedBoundsInOutput }, + label + ); + + return { scalars: outputData, extent }; +} + +const workerApi = { + gaussianSmoothLabelMapWorker, +}; + +Comlink.expose(workerApi); diff --git a/src/segmentation/editing/coordinator.ts b/src/segmentation/editing/coordinator.ts new file mode 100644 index 000000000..0bb466539 --- /dev/null +++ b/src/segmentation/editing/coordinator.ts @@ -0,0 +1,29 @@ +import { defineStore } from 'pinia'; + +/** Coordinates temporary mask previews with competing edits and durable reads. */ +export const useSegmentationEditsStore = defineStore( + 'segmentationEdits', + () => { + let cancelPreview: (() => void) | undefined; + + function beforeEdit() { + const cancel = cancelPreview; + cancelPreview = undefined; + cancel?.(); + } + + function hold(cancel: () => void) { + beforeEdit(); + cancelPreview = cancel; + } + + function release(cancel: () => void) { + if (cancelPreview === cancel) cancelPreview = undefined; + } + + // Save and export read committed voxels, resolving an unconfirmed preview. + const beforeRead = beforeEdit; + + return { beforeEdit, beforeRead, hold, release }; + } +); diff --git a/src/segmentation/geometry.ts b/src/segmentation/geometry.ts new file mode 100644 index 000000000..7a2f4fe65 --- /dev/null +++ b/src/segmentation/geometry.ts @@ -0,0 +1,174 @@ +/** vtk.js index-space extent order: [iMin, iMax, jMin, jMax, kMin, kMax]. */ +export type Extent3D = [number, number, number, number, number, number]; + +/** Widens `box` in place to take in one more index. */ +export const growExtent = (box: Extent3D, i: number, j: number, k: number) => { + box[0] = Math.min(box[0], i); + box[1] = Math.max(box[1], i); + box[2] = Math.min(box[2], j); + box[3] = Math.max(box[3], j); + box[4] = Math.min(box[4], k); + box[5] = Math.max(box[5], k); +}; + +export function emptyExtent(): Extent3D { + return [0, -1, 0, -1, 0, -1]; +} + +/** vtk.js extents are inclusive, so an axis is empty only when max < min. */ +export function isEmptyExtent(extent: Extent3D) { + return ( + extent[1] < extent[0] || extent[3] < extent[2] || extent[5] < extent[4] + ); +} + +/** Voxel counts along i, j, k. Meaningless for an empty extent. */ +export function extentSize(extent: Extent3D) { + return [ + extent[1] - extent[0] + 1, + extent[3] - extent[2] + 1, + extent[5] - extent[4] + 1, + ] as [number, number, number]; +} + +export function extentContains(outer: Extent3D, inner: Extent3D) { + return ( + inner[0] >= outer[0] && + inner[1] <= outer[1] && + inner[2] >= outer[2] && + inner[3] <= outer[3] && + inner[4] >= outer[4] && + inner[5] <= outer[5] + ); +} + +export function extentContainsIndex( + extent: Extent3D, + i: number, + j: number, + k: number +) { + return ( + i >= extent[0] && + i <= extent[1] && + j >= extent[2] && + j <= extent[3] && + k >= extent[4] && + k <= extent[5] + ); +} + +/** A mask's extent with the row and plane strides that extent implies. */ +export type MaskBounds = { + extent: Extent3D; + mi: number; + mj: number; +}; + +/** Where a parent-index voxel sits in the buffer of a mask bounded that way. */ +export const maskOffset = ( + bounds: MaskBounds, + i: number, + j: number, + k: number +) => + i - + bounds.extent[0] + + (j - bounds.extent[2]) * bounds.mi + + (k - bounds.extent[4]) * bounds.mi * bounds.mj; + +/** + * The box `labelValue` actually occupies inside a mask bounded by `extent`, + * empty when it occupies nothing. A binding's extent is the allocation, padded + * and never shrunk by an erase, so it is not the segment's bounds. + */ +export function markedExtent( + scalars: ArrayLike, + extent: Extent3D, + labelValue: number +): Extent3D { + const ni = extent[1] - extent[0] + 1; + const nj = extent[3] - extent[2] + 1; + const nk = extent[5] - extent[4] + 1; + let bounds: Extent3D | undefined; + + const scanRow = (rowStart: number, j: number, k: number) => { + for (let index = 0; index < ni; index += 1) { + if (scalars[rowStart + index] !== labelValue) continue; + const i = extent[0] + index; + if (bounds) growExtent(bounds, i, j, k); + else bounds = [i, i, j, j, k, k]; + } + }; + + for (let row = 0; row < nj * nk; row += 1) { + scanRow(row * ni, extent[2] + (row % nj), extent[4] + Math.floor(row / nj)); + } + + return bounds ?? emptyExtent(); +} + +/** The parent-image slice indices containing `labelValue`, for i, j and k. */ +export function markedSlices( + scalars: ArrayLike, + extent: Extent3D, + labelValue: number +): [number[], number[], number[]] { + const ni = extent[1] - extent[0] + 1; + const nj = extent[3] - extent[2] + 1; + const occupied = [new Set(), new Set(), new Set()]; + + for (let offset = 0; offset < scalars.length; offset += 1) { + if (scalars[offset] !== labelValue) continue; + const i = extent[0] + (offset % ni); + const row = Math.floor(offset / ni); + const j = extent[2] + (row % nj); + const k = extent[4] + Math.floor(row / nj); + occupied[0].add(i); + occupied[1].add(j); + occupied[2].add(k); + } + + return occupied.map((slices) => [...slices]) as [ + number[], + number[], + number[], + ]; +} + +export function extentUnion(a: Extent3D, b: Extent3D): Extent3D { + return [ + Math.min(a[0], b[0]), + Math.max(a[1], b[1]), + Math.min(a[2], b[2]), + Math.max(a[3], b[3]), + Math.min(a[4], b[4]), + Math.max(a[5], b[5]), + ]; +} + +export function padExtent(extent: Extent3D, padding: number): Extent3D { + return [ + extent[0] - padding, + extent[1] + padding, + extent[2] - padding, + extent[3] + padding, + extent[4] - padding, + extent[5] + padding, + ]; +} + +export function clipExtent(extent: Extent3D, bounds: Extent3D): Extent3D { + return [ + Math.max(extent[0], bounds[0]), + Math.min(extent[1], bounds[1]), + Math.max(extent[2], bounds[2]), + Math.min(extent[3], bounds[3]), + Math.max(extent[4], bounds[4]), + Math.min(extent[5], bounds[5]), + ]; +} + +export function fullExtent(dimensions: number[] | Int32Array): Extent3D { + return [0, dimensions[0] - 1, 0, dimensions[1] - 1, 0, dimensions[2] - 1]; +} diff --git a/src/segmentation/io/__tests__/export.spec.ts b/src/segmentation/io/__tests__/export.spec.ts new file mode 100644 index 000000000..fbbb7a658 --- /dev/null +++ b/src/segmentation/io/__tests__/export.spec.ts @@ -0,0 +1,58 @@ +import { describe, expect, it } from 'vitest'; +import JSZip from 'jszip'; + +import { bundleExportFiles, layerFileName } from '@/src/segmentation/io/export'; + +// A labelmap file carries one label per voxel, so a segmentation with overlap +// leaves as several files. One file is the common case and stays the download +// it has always been: same name, same bytes, no archive around it. + +const bytes = (...values: number[]) => new Uint8Array(values); + +const readBlob = async (blob: Blob) => + Array.from(new Uint8Array(await blob.arrayBuffer())); + +describe('naming the file each group of segments writes', () => { + it('gives the first group the plain name', () => { + expect(layerFileName('Prostate', 'seg.nrrd', 0)).toBe('Prostate.seg.nrrd'); + }); + + it('numbers every later group after it', () => { + expect(layerFileName('Prostate', 'seg.nrrd', 1)).toBe( + 'Prostate_layer1.seg.nrrd' + ); + expect(layerFileName('Prostate', 'nii.gz', 2)).toBe( + 'Prostate_layer2.nii.gz' + ); + }); +}); + +describe('handing the written files to the browser', () => { + it('downloads a single file as itself', async () => { + const bundle = await bundleExportFiles('Prostate', [ + { name: 'Prostate.seg.nrrd', data: bytes(1, 2, 3) }, + ]); + + expect(bundle.name).toBe('Prostate.seg.nrrd'); + expect(await readBlob(bundle.blob)).toEqual([1, 2, 3]); + }); + + it('downloads several files as one archive of them', async () => { + const bundle = await bundleExportFiles('Prostate', [ + { name: 'Prostate.seg.nrrd', data: bytes(1, 2, 3) }, + { name: 'Prostate_layer1.seg.nrrd', data: bytes(4, 5) }, + ]); + + expect(bundle.name).toBe('Prostate.zip'); + const zip = await JSZip.loadAsync(bundle.blob); + expect(Object.keys(zip.files)).toEqual([ + 'Prostate.seg.nrrd', + 'Prostate_layer1.seg.nrrd', + ]); + expect( + Array.from( + await zip.files['Prostate_layer1.seg.nrrd'].async('uint8array') + ) + ).toEqual([4, 5]); + }); +}); diff --git a/src/segmentation/io/export.ts b/src/segmentation/io/export.ts new file mode 100644 index 000000000..08c1a7ed9 --- /dev/null +++ b/src/segmentation/io/export.ts @@ -0,0 +1,34 @@ +import JSZip from 'jszip'; + +export type ExportFile = { + name: string; + data: string | Uint8Array; +}; + +/** + * The name a group's file carries. The first keeps the plain stem, so a + * segmentation with no overlap saves as the one file it always did. + */ +export const layerFileName = (stem: string, format: string, layer: number) => + layer === 0 ? `${stem}.${format}` : `${stem}_layer${layer}.${format}`; + +export const archiveNameFor = (stem: string) => `${stem}.zip`; + +/** + * What a save hands to the browser: the single file itself, or every file in + * one archive. A labelmap file carries one label per voxel, so segments that + * overlap cannot share one and the save turns into several. + */ +export async function bundleExportFiles(stem: string, files: ExportFile[]) { + const [first] = files; + if (files.length === 1) { + return { name: first.name, blob: new Blob([first.data]) }; + } + + const zip = new JSZip(); + files.forEach((file) => zip.file(file.name, file.data)); + return { + name: archiveNameFor(stem), + blob: await zip.generateAsync({ type: 'blob' }), + }; +} diff --git a/src/segmentation/io/maskFileNaming.ts b/src/segmentation/io/maskFileNaming.ts new file mode 100644 index 000000000..873b4b76a --- /dev/null +++ b/src/segmentation/io/maskFileNaming.ts @@ -0,0 +1,27 @@ +const defaultName = (baseName: string, index: number) => + `Segment Group ${index} for ${baseName}`; + +/** + * Default names for mask files, counted per parent image. The + * count keeps rising so a deleted mask's name is not immediately handed to + * the next one, and `taken` skips a name something already holds. + */ +export function createMaskFileNamer(taken: () => Set) { + const nextIndex: Record = Object.create(null); + return { + pick(parentImageId: string, baseName: string) { + const held = taken(); + let name = ''; + do { + const index = nextIndex[parentImageId] ?? 1; + nextIndex[parentImageId] = index + 1; + name = defaultName(baseName, index); + } while (held.has(name)); + return name; + }, + /** Called by the deletion cascade, so a removed image stops counting. */ + forget(parentImageId: string) { + delete nextIndex[parentImageId]; + }, + }; +} diff --git a/src/segmentation/masks/labelValue.ts b/src/segmentation/masks/labelValue.ts new file mode 100644 index 000000000..5a72c3557 --- /dev/null +++ b/src/segmentation/masks/labelValue.ts @@ -0,0 +1,33 @@ +/** A composed labelmap is written as Uint8, so a label value has to fit in one byte. */ +export const LABELMAP_MAX_VALUE = 255; + +/** + * The value every mask marks its own voxels with. A mask holds one segment and + * nothing else, so the byte says only claimed or not; which segment it belongs + * to is the mask's identity, not its content. Export assigns its own values at + * write time, where one file does have to tell segments apart per voxel. + */ +export const SEGMENT_VALUE = 1; + +import { LABELMAP_BACKGROUND_VALUE } from '@/src/segmentation/model'; + +export function nextUnusedLabelValue( + used: Set, + maximum: number, + preferred?: number +) { + if ( + preferred !== undefined && + preferred > LABELMAP_BACKGROUND_VALUE && + preferred <= maximum && + !used.has(preferred) + ) { + return preferred; + } + let labelValue = LABELMAP_BACKGROUND_VALUE + 1; + while (used.has(labelValue)) labelValue += 1; + if (labelValue > maximum) { + throw new Error(`An image holds at most ${maximum} segments in a labelmap`); + } + return labelValue; +} diff --git a/src/segmentation/masks/overlap.ts b/src/segmentation/masks/overlap.ts new file mode 100644 index 000000000..5be07d464 --- /dev/null +++ b/src/segmentation/masks/overlap.ts @@ -0,0 +1,197 @@ +import type vtkLabelMap from '@/src/vtk/LabelMap'; +import { + LABELMAP_BACKGROUND_VALUE, + maskScalars, +} from '@/src/segmentation/model'; +import { + clipExtent, + extentContainsIndex, + extentSize, + isEmptyExtent, + maskOffset, + type Extent3D, + type MaskBounds, +} from '@/src/segmentation/geometry'; + +// Bounded masks read as one parent-shaped picture: how a mask is written into +// that picture, and which masks can share one without losing a voxel. + +export type BoundedScalars = MaskBounds & { + mask: vtkLabelMap; + scalars: Uint8Array; +}; + +/** + * The extent is copied because the callers read it per voxel and a segment's + * own copy lives in the reactive tree. + */ +export function boundScalars( + mask: vtkLabelMap | undefined, + bounds: Extent3D +): BoundedScalars | undefined { + if (!mask || isEmptyExtent(bounds)) return undefined; + const extent = [...bounds] as Extent3D; + const [mi, mj] = extentSize(extent); + return { mask, scalars: maskScalars(mask), extent, mi, mj }; +} + +/** + * Clipping once here keeps a mask that misses the box out of the per-voxel + * containment test, and lets the caller skip the walk when none is left. + */ +const masksReaching = (masks: BoundedScalars[], within: Extent3D) => + masks.filter((bounded) => !isEmptyExtent(clipExtent(bounded.extent, within))); + +/** + * Whether any of these masks holds the voxel at PARENT indices i, j, k, over + * the box the caller is about to walk. Absent when no mask reaches that box. + */ +export function masksHolding(masks: BoundedScalars[], within: Extent3D) { + const reaching = masksReaching(masks, within); + if (reaching.length === 0) return undefined; + return (i: number, j: number, k: number) => + reaching.some( + (bounded) => + extentContainsIndex(bounded.extent, i, j, k) && + bounded.scalars[maskOffset(bounded, i, j, k)] !== + LABELMAP_BACKGROUND_VALUE + ); +} + +/** + * Clears the voxel at PARENT indices i, j, k from every one of these masks and + * answers true: nothing left here can refuse the write, which is the same + * per-voxel answer the occupancy test gives. Absent when no mask reaches + * `within`, the box the caller is about to walk. A mask that does not reach the + * voxel has nothing there to clear, so nothing grows. Finish the operation in + * a finally block to publish each changed mask once, including partial writes. + */ +export function masksClearing(masks: BoundedScalars[], within: Extent3D) { + const reaching = masksReaching(masks, within); + if (reaching.length === 0) return undefined; + const changed = new Set(); + const claim = (i: number, j: number, k: number) => { + reaching.forEach((bounded) => { + if (!extentContainsIndex(bounded.extent, i, j, k)) return; + const offset = maskOffset(bounded, i, j, k); + if (bounded.scalars[offset] === LABELMAP_BACKGROUND_VALUE) return; + bounded.scalars[offset] = LABELMAP_BACKGROUND_VALUE; + changed.add(bounded.mask); + }); + return true; + }; + return { + claim, + finish: () => { + changed.forEach((mask) => mask.modified()); + changed.clear(); + }, + }; +} + +/** + * A mask's buffer, positioned where one row of the shared box starts in it. + */ +type MaskRow = { scalars: Uint8Array; from: number }; + +const rowAt = ( + bounded: BoundedScalars, + i: number, + j: number, + k: number +): MaskRow => ({ + scalars: bounded.scalars, + from: maskOffset(bounded, i, j, k), +}); + +/** + * Whether both rows hold a voxel at the same step along i. Background is 0, so + * a claimed voxel is a truthy one. + */ +function rowsIntersect(a: MaskRow, b: MaskRow, count: number) { + const { scalars: av, from: ai } = a; + const { scalars: bv, from: bi } = b; + for (let n = 0; n < count; n += 1) { + if (av[ai + n] && bv[bi + n]) return true; + } + return false; +} + +/** + * Whether two masks claim one voxel in common. A mask is bounded to the voxels + * it holds, so boxes that miss cannot share a voxel and neither buffer is read. + * Boxes that meet are swept a row at a time: two segments that touch nowhere + * usually still share a box, so the sweep is the common case, and each row + * costs one offset per mask with a plain step along i from there. + */ +export function masksIntersect(a: BoundedScalars, b: BoundedScalars) { + const shared = clipExtent(a.extent, b.extent); + if (isEmptyExtent(shared)) return false; + + const [ni] = extentSize(shared); + const i = shared[0]; + for (let k = shared[4]; k <= shared[5]; k += 1) { + for (let j = shared[2]; j <= shared[3]; j += 1) { + if (rowsIntersect(rowAt(a, i, j, k), rowAt(b, i, j, k), ni)) return true; + } + } + return false; +} + +/** + * Items grouped so no group holds two masks that claim a voxel in common. + * Greedy first fit: an item takes the lowest group it does not intersect, so a + * segmentation with no overlap stays one group, in order. An item with no mask + * claims nothing and joins the first group. + */ +export function groupByLayer( + items: T[], + maskOf: (item: T) => BoundedScalars | undefined +) { + const layers: Array<{ items: T[]; masks: BoundedScalars[] }> = []; + const fits = ( + layer: { masks: BoundedScalars[] }, + mask: BoundedScalars | undefined + ) => !mask || layer.masks.every((other) => !masksIntersect(mask, other)); + + items.forEach((item) => { + const mask = maskOf(item); + const found = layers.find((layer) => fits(layer, mask)); + const layer = found ?? { items: [], masks: [] }; + if (!found) layers.push(layer); + layer.items.push(item); + if (mask) layer.masks.push(mask); + }); + + return layers.map((layer) => layer.items); +} + +/** + * Marks a bounded mask's voxels in a parent-shaped buffer as `labelValue`. The + * mask's own bytes only say claimed or not, so the value is the caller's. + * Masks are written in `order`, so a later one takes a voxel an earlier one + * also claims. + */ +export function writeMaskInto( + values: Uint8Array, + dimensions: readonly number[], + bounded: BoundedScalars, + labelValue: number +) { + const { extent, scalars } = bounded; + const [dx, dy] = dimensions; + const [ni, nj, nk] = extentSize(extent); + // Rows flat in one loop: a j loop inside a k loop would nest deeper than the + // style allows once the row test is in it. + for (let row = 0; row < nj * nk; row += 1) { + const j = extent[2] + (row % nj); + const k = extent[4] + Math.floor(row / nj); + const to = extent[0] + j * dx + k * dx * dy; + const from = maskOffset(bounded, extent[0], j, k); + for (let n = 0; n < ni; n += 1) { + // Background is 0, so a voxel this mask leaves unclaimed keeps whatever + // the buffer already holds there. + if (scalars[from + n]) values[to + n] = labelValue; + } + } +} diff --git a/src/segmentation/masks/storage.ts b/src/segmentation/masks/storage.ts new file mode 100644 index 000000000..1e7846132 --- /dev/null +++ b/src/segmentation/masks/storage.ts @@ -0,0 +1,101 @@ +import vtkDataArray from '@kitware/vtk.js/Common/Core/DataArray'; +import type vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; +import type { TypedArray, Vector3 } from '@kitware/vtk.js/types'; + +import { maskScalars } from '@/src/segmentation/model'; +import { + clipExtent, + extentSize, + isEmptyExtent, + maskOffset, + type Extent3D, +} from '@/src/segmentation/geometry'; +import vtkLabelMap from '@/src/vtk/LabelMap'; + +export const setMaskScalars = (mask: vtkLabelMap, values: Uint8Array) => + mask + .getPointData() + .setScalars(vtkDataArray.newInstance({ numberOfComponents: 1, values })); + +/** Places a bounded mask on the same index grid as its parent image. */ +export function placeMask( + mask: vtkLabelMap, + parent: vtkImageData, + extent: Extent3D +) { + const dimensions = isEmptyExtent(extent) ? [0, 0, 0] : extentSize(extent); + const origin = isEmptyExtent(extent) + ? Array.from(parent.getOrigin()) + : Array.from( + parent.indexToWorld([extent[0], extent[2], extent[4]] as Vector3) + ); + mask.setOrigin(origin as Vector3); + mask.setDimensions(dimensions as Vector3); + mask.computeTransforms(); + return dimensions; +} + +export function allocateMask(parent: vtkImageData, extent: Extent3D) { + const mask = vtkLabelMap.newInstance( + parent.get('spacing', 'origin', 'direction') + ); + const dimensions = placeMask(mask, parent, extent); + setMaskScalars( + mask, + new Uint8Array(dimensions[0] * dimensions[1] * dimensions[2]) + ); + return mask; +} + +/** + * Copies a mask onto another extent of its parent grid, padding with zero. + * An output buffer must match the destination size and not alias the source. + */ +export function reframeMaskScalars( + scalars: TypedArray | number[], + from: Extent3D, + to: Extent3D, + output?: Uint8Array +) { + const [mi, mj, mk] = extentSize(to); + const size = isEmptyExtent(to) ? 0 : mi * mj * mk; + const values = output ?? new Uint8Array(size); + if (values.length !== size) throw new Error('Mask output size mismatch'); + if (output) values.fill(0); + const shared = clipExtent(from, to); + const [si, sj] = extentSize(from); + // Bounds can be reactive; read them once before the per-row copy loop. + const source = { extent: [...from] as Extent3D, mi: si, mj: sj }; + const destination = { extent: [...to] as Extent3D, mi, mj }; + if (isEmptyExtent(shared)) return values; + const count = shared[1] - shared[0] + 1; + for (let k = shared[4]; k <= shared[5]; k += 1) { + for (let j = shared[2]; j <= shared[3]; j += 1) { + const start = maskOffset(source, shared[0], j, k); + const end = maskOffset(destination, shared[0], j, k); + if (count === 1) { + values[end] = scalars[start]; + } else { + const row = Array.isArray(scalars) + ? scalars.slice(start, start + count) + : scalars.subarray(start, start + count); + values.set(row, end); + } + } + } + return values; +} + +/** Preserves the vtk image instance while replacing its scalar storage. */ +export function regrowMask( + mask: vtkLabelMap, + parent: vtkImageData, + from: Extent3D, + to: Extent3D +) { + const previous = maskScalars(mask); + const values = reframeMaskScalars(previous, from, to); + placeMask(mask, parent, to); + setMaskScalars(mask, values); + mask.modified(); +} diff --git a/src/segmentation/masks/voxelAccess.ts b/src/segmentation/masks/voxelAccess.ts new file mode 100644 index 000000000..f29e642cf --- /dev/null +++ b/src/segmentation/masks/voxelAccess.ts @@ -0,0 +1,216 @@ +import type { TypedArray } from '@kitware/vtk.js/types'; + +import type { Maybe } from '@/src/types'; +import type { VoxelGesture } from '@/src/segmentation/model'; +import type { useImageCacheStore } from '@/src/store/image-cache'; +import { regrowMask } from '@/src/segmentation/masks/storage'; +import { + boundScalars, + masksClearing, + masksHolding, +} from '@/src/segmentation/masks/overlap'; +import { + listMasks, + maskScalars, + type LabelmapBinding, + type MaskVoxelAccessor, + type SegmentMask, + type Segmentation, + type VoxelStorage, +} from '@/src/segmentation/model'; +import { + clipExtent, + extentContains, + extentUnion, + fullExtent, + isEmptyExtent, + padExtent, + type Extent3D, +} from '@/src/segmentation/geometry'; + +export type VoxelAccessDeps = { + imageCacheStore: ReturnType; + findMask: (maskId: string) => SegmentMask | undefined; + getMask: (maskId: string) => SegmentMask; + segmentationOfMask: (maskId: string) => Segmentation | undefined; + ensureLabelmapBinding: (maskId: string) => LabelmapBinding; + maskLocked: (mask: SegmentMask) => boolean; +}; + +/** + * Reading and growing the voxels behind a mask. Split out so the store holds + * the records; every accessor re-resolves its binding rather than capturing a + * buffer, so none of them outlive a mask they were made for. + */ +export function createVoxelAccess(deps: VoxelAccessDeps) { + const { + imageCacheStore, + findMask, + getMask, + segmentationOfMask, + ensureLabelmapBinding, + maskLocked, + } = deps; + + function requireParentImage(maskId: string) { + const segmentation = segmentationOfMask(maskId); + if (!segmentation) throw new Error('No such segment'); + const parent = imageCacheStore.getVtkImageData(segmentation.parentImageId); + if (!parent) throw new Error('No such parent image'); + return parent; + } + + /** + * Grows one mask, in place, to cover `extent` in parent index space, with + * `padding` voxels of room beyond it when it has to grow at all. + */ + function ensureMaskContains(maskId: string, extent: Extent3D, padding = 0) { + if (isEmptyExtent(extent)) return false; + + const binding = getMask(maskId).representations.labelmap; + if (!binding) throw new Error('No storage: call materialize() first'); + const parent = requireParentImage(maskId); + // Refused before anything is touched, so a rejected growth leaves the mask + // exactly as it was. + const parentExtent = fullExtent(parent.getDimensions()); + if (!extentContains(parentExtent, extent)) + throw new Error('Extent leaves the parent image'); + + const current = binding.extent; + if (!isEmptyExtent(current) && extentContains(current, extent)) + return false; + + const requested = clipExtent(padExtent(extent, padding), parentExtent); + const grown = isEmptyExtent(current) + ? requested + : extentUnion(current, requested); + regrowMask(binding.image, parent, current, grown); + binding.extent = grown; + return true; + } + + /** + * The voxel half of the accessor seam, over whichever mask `findBinding` + * resolves. Resolution is deferred to every call so a stale accessor sees + * deletion or growth done through another one. `onMissing` names why storage + * is unreachable, so `exists()` can answer without throwing. + */ + function voxelStorage( + maskId: string, + findBinding: () => Maybe, + onMissing: () => never + ): VoxelStorage { + const findImage = () => findBinding()?.image; + const requireImage = () => findImage() ?? onMissing(); + const requireScalars = () => maskScalars(requireImage()); + + return { + exists: () => !!findImage(), + image: requireImage, + scalars: requireScalars, + snapshot: () => requireScalars().slice(), + apply: (scalars: TypedArray | number[]) => { + const image = requireImage(); + const data = maskScalars(image); + if (scalars.length !== data.length) { + throw new Error('Scalar length does not match storage'); + } + data.set(scalars); + image.modified(); + }, + ensureContains: (extent: Extent3D, padding = 0) => { + requireImage(); + return ensureMaskContains(maskId, extent, padding); + }, + }; + } + + function maskVoxels(maskId: string): MaskVoxelAccessor { + // Validates eagerly: an accessor for a nonexistent segment is refused up + // front, not just on first use. + getMask(maskId); + + const binding = () => getMask(maskId).representations.labelmap; + + // Deliberately tolerant where binding() is not: the segment itself can be + // deleted out from under an accessor, and that is an absent storage, not a + // lookup error. + const findBinding = () => findMask(maskId)?.representations.labelmap; + + const onMissing = (): never => { + throw new Error('No storage: call materialize() first'); + }; + + return { + binding, + materialize: () => ensureLabelmapBinding(maskId), + ...voxelStorage(maskId, findBinding, onMissing), + }; + } + + /** + * The accessor for consumers holding an id a segment may already have left: + * the renderer and the paint widget are computeds keyed on one that can + * vanish a tick before they do, so this stays constructible either way. + */ + const findMaskVoxels = (maskId: string) => + voxelStorage( + maskId, + () => findMask(maskId)?.representations.labelmap, + () => { + throw new Error('No such segment'); + } + ); + + /** A bound segment's buffer, absent when it has none or holds nothing. */ + const boundedMask = (binding?: LabelmapBinding) => + binding && boundScalars(binding.image, binding.extent); + + /** + * The masks of an image's other segments that `gesture` may take a voxel + * from, resolved once per run because the caller below runs per voxel. A + * locked segment is not editable, so an aimed gesture is not offered its mask + * at all. + */ + function siblingMasks(maskId: string, gesture: VoxelGesture) { + const segmentation = segmentationOfMask(maskId); + if (!segmentation) return []; + return listMasks(segmentation).flatMap((segment) => { + if (segment.id === maskId) return []; + if (gesture === 'aimed' && maskLocked(segment)) return []; + const bounded = boundedMask(segment.representations.labelmap); + return bounded ? [bounded] : []; + }); + } + + /** + * Whether the voxel at PARENT indices i, j, k is this segment's to write, + * taking it from the neighbours that have to yield it. Absent when no other + * segment reaches `within`, the box the caller is about to walk: every voxel + * in it is then uncontested and the question need not be asked per voxel. + * + * `gesture` is the whole of the policy, so see {@link VoxelGesture}. An aimed + * operation must call finish in a finally block after its last voxel write. + */ + function voxelClaim(maskId: string, gesture: VoxelGesture, within: Extent3D) { + const masks = siblingMasks(maskId, gesture); + if (gesture === 'aimed') return masksClearing(masks, within); + const held = masksHolding(masks, within); + return ( + held && { + claim: (i: number, j: number, k: number) => !held(i, j, k), + finish: () => undefined, + } + ); + } + + return { + requireParentImage, + ensureMaskContains, + maskVoxels, + findMaskVoxels, + boundedMask, + siblingMasks, + voxelClaim, + }; +} diff --git a/src/segmentation/model.ts b/src/segmentation/model.ts new file mode 100644 index 000000000..76f62d905 --- /dev/null +++ b/src/segmentation/model.ts @@ -0,0 +1,127 @@ +import type { Extent3D } from '@/src/segmentation/geometry'; +import type { ProcessingResultSource } from '@/src/types'; +import type { RGBAColor, TypedArray } from '@kitware/vtk.js/types'; + +import type vtkLabelMap from '@/src/vtk/LabelMap'; + +/** A fresh segmentation tints the anatomy under it rather than hiding it. */ +export const DEFAULT_SEGMENTATION_FILL_OPACITY = 0.3; + +export type LabelmapBinding = { + /** + * This mask's voxels, and no other segment's. Held raw: a vtk object must + * not be proxied, so every writer of a binding marks it. + */ + image: vtkLabelMap; + extent: Extent3D; // the mask's own bounds, in parent image index space + /** Reaches the saved archive's entry path, so a round trip keeps it. */ + name: string; + source?: ProcessingResultSource; +}; + +/** + * One image's mask for one segment type. Its id is its own, distinct from the + * type id: everything the user sees or sets, visibility and lock included, + * lives on the type, so this record is storage and nothing else. + */ +export type SegmentMask = { + id: string; + segmentId: string; + representations: { + // absent until voxels are allocated + labelmap?: LabelmapBinding; + }; +}; + +export const LABELMAP_BACKGROUND_VALUE = 0; + +export const makeDefaultSegmentName = (value: number) => `Segment ${value}`; + +/** + * One mask's label descriptor, derived from the segment it delineates. + * Identity lives on `Segment`; this is the value-keyed view the labelmap + * renderer and the .seg.nrrd writer consume. + */ +export type LabelmapSegment = { + value: number; + name: string; + color: RGBAColor; + visible: boolean; + locked?: boolean; + // Absent on descriptors that come off a file rather than off a segment. + fillOpacity?: number; + outlineOpacity?: number; +}; + +/** vtk declares getData() as number[] | TypedArray; mask storage is typed. */ +export const maskScalars = (mask: vtkLabelMap) => + mask.getPointData().getScalars().getData() as Uint8Array; + +export type Segmentation = { + id: string; + name: string; + parentImageId: string; + masks: Record; + order: string[]; + fillOpacity: number; + outlineOpacity: number; + outlineThickness: number; +}; + +/** The display multipliers every segment of a segmentation is scaled by. */ +export type SegmentationDisplayPatch = Partial< + Pick +>; + +/** + * The voxel operations every labelmap consumer routes through. Storage is one + * bounded mask per segment, sized to the region that segment covers. + * + * `ensureContains` may replace the scalar array, dimensions and strides: + * anything that cached those from `image()` or `scalars()` must re-fetch after + * calling it. + */ +export type VoxelStorage = { + /** + * Whether the storage is still reachable. An accessor outlives what it + * points at, so callers holding one across a deletion check this before a + * read or a write; every other method throws while it is false. + */ + exists(): boolean; + image(): vtkLabelMap; + /** Live mask buffer. Writers publish changes through apply() or image().modified(). */ + scalars(): TypedArray; + snapshot(): TypedArray; + /** Bulk copy-in; keeps image() and scalars() identity, marks it modified. */ + apply(scalars: TypedArray | number[]): void; + /** + * Ensures storage covers `extent`, growing to the union of what it has and + * what it was asked for. Returns whether storage was invalidated + * (scalars/dimensions/strides changed). An empty extent is already covered. + * Throws when the extent leaves the parent image, so callers clip. When the + * extent is not already covered, the mask grows by `padding` voxels beyond + * it on every face (clipped to the parent), so nearby requests that follow + * grow nothing. + */ + ensureContains(extent: Extent3D, padding?: number): boolean; +}; + +/** + * Voxel access for one segment. Re-resolves the binding on every call rather + * than capturing it, so a caller that holds an accessor across a segment + * deletion or a growth sees the current state, not a stale one. `exists()` is + * false, and every storage method throws, before `materialize()`. + */ +export type MaskVoxelAccessor = VoxelStorage & { + binding(): LabelmapBinding | undefined; + /** Allocates storage if needed and returns the binding. Idempotent. */ + materialize(): LabelmapBinding; +}; + +/** Segments in display order. `order` is the authority, `segments` the store. */ +export function listMasks(segmentation: Segmentation) { + return segmentation.order.map((id) => segmentation.masks[id]); +} + +/** Aimed writes clear unlocked neighbors; sweeps only grow into unclaimed voxels. */ +export type VoxelGesture = 'aimed' | 'sweep'; diff --git a/src/segmentation/rendering/__tests__/segmentDisplay.spec.ts b/src/segmentation/rendering/__tests__/segmentDisplay.spec.ts new file mode 100644 index 000000000..bfb018c2a --- /dev/null +++ b/src/segmentation/rendering/__tests__/segmentDisplay.spec.ts @@ -0,0 +1,129 @@ +import { describe, expect, it } from 'vitest'; + +import { + SEGMENT_ACTOR_OPACITY, + segmentFillAlpha, + segmentOutlineTables, +} from '@/src/segmentation/rendering/display'; +import type { LabelmapSegment } from '@/src/segmentation/model'; + +const makeMask = ( + value: number, + overrides: Partial = {} +): LabelmapSegment => ({ + value, + name: `Segment ${value}`, + color: [255, 0, 0, 255], + visible: true, + ...overrides, +}); + +describe('segmentFillAlpha', () => { + it('is the segment alpha when the fill is fully opaque', () => { + expect(segmentFillAlpha(makeMask(1, { fillOpacity: 1 }))).toBe(1); + }); + + it('scales the segment alpha by the fill opacity', () => { + expect(segmentFillAlpha(makeMask(1, { fillOpacity: 0.5 }))).toBe(0.5); + }); + + it('hides a fill the user set to zero', () => { + expect(segmentFillAlpha(makeMask(1, { fillOpacity: 0 }))).toBe(0); + }); + + it('hides an invisible segment whatever its fill opacity', () => { + expect( + segmentFillAlpha(makeMask(1, { visible: false, fillOpacity: 1 })) + ).toBe(0); + }); + + it('treats a descriptor without a fill opacity as opaque', () => { + expect(segmentFillAlpha(makeMask(1))).toBe(1); + }); + + it('scales the segment alpha by the segmentation\u2019s fill opacity', () => { + expect(segmentFillAlpha(makeMask(1, { fillOpacity: 0.5 }), 0.5)).toBe(0.25); + }); + + it('hides every fill when the segmentation\u2019s fill opacity is zero', () => { + expect(segmentFillAlpha(makeMask(1, { fillOpacity: 1 }), 0)).toBe(0); + }); +}); + +describe('segmentOutlineTables', () => { + it('indexes both tables by label value minus one', () => { + const tables = segmentOutlineTables( + [ + makeMask(1, { outlineOpacity: 0.25 }), + makeMask(2, { outlineOpacity: 0.5 }), + ], + 2, + 1 + ); + + expect(tables.opacities).toEqual([0.25, 0.5]); + expect(tables.thicknesses).toEqual([2, 2]); + }); + + it('hides an outline the user set to zero', () => { + const tables = segmentOutlineTables( + [makeMask(1, { outlineOpacity: 0 })], + 2, + 1 + ); + + expect(tables.opacities).toEqual([0]); + }); + + it('scales every segment by the group outline opacity', () => { + const tables = segmentOutlineTables( + [makeMask(1, { outlineOpacity: 0.5 })], + 2, + 0.5 + ); + + expect(tables.opacities).toEqual([0.25]); + }); + + it('leaves values no segment claims at the group defaults', () => { + const tables = segmentOutlineTables( + [makeMask(3, { outlineOpacity: 0.5 })], + 2, + 1 + ); + + expect(tables.opacities).toEqual([1, 1, 0.5]); + expect(tables.thicknesses).toEqual([2, 2, 2]); + }); + + it('drops the thickness of an invisible segment', () => { + const tables = segmentOutlineTables( + [makeMask(1, { visible: false }), makeMask(2)], + 2, + 1 + ); + + expect(tables.thicknesses).toEqual([0, 2]); + }); + + it('has no entries for an artifact with no bound segments', () => { + expect(segmentOutlineTables([], 2, 1)).toEqual({ + thicknesses: [], + opacities: [], + }); + }); +}); + +describe('SEGMENT_ACTOR_OPACITY', () => { + it('leaves the fill to the transfer functions', () => { + // A fully opaque segment reaches the screen at its own alpha, so the actor + // must not scale it down. + expect(SEGMENT_ACTOR_OPACITY).toBeGreaterThan(0.999); + }); + + it('stays out of the opaque render pass', () => { + // vtk.js treats an image slice at opacity 1 as opaque and restacks it + // against the base image and the sibling segment actors. + expect(SEGMENT_ACTOR_OPACITY).toBeLessThan(1); + }); +}); diff --git a/src/segmentation/rendering/__tests__/segmentRenderMask.spec.ts b/src/segmentation/rendering/__tests__/segmentRenderMask.spec.ts new file mode 100644 index 000000000..efb2930fb --- /dev/null +++ b/src/segmentation/rendering/__tests__/segmentRenderMask.spec.ts @@ -0,0 +1,114 @@ +import { describe, expect, it } from 'vitest'; +import vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; +import { allocateMask, regrowMask } from '@/src/segmentation/masks/storage'; +import { maskScalars } from '@/src/segmentation/model'; +import { type Extent3D } from '@/src/segmentation/geometry'; +import { segmentRenderMask } from '@/src/segmentation/rendering/renderMask'; + +function scene(extent: Extent3D) { + const parent = vtkImageData.newInstance({ + spacing: [2, 3, 4], + origin: [10, 20, 30], + direction: [0, 1, 0, -1, 0, 0, 0, 0, 1], + }); + parent.setExtent(4, 9, 5, 10, 6, 11); + const source = allocateMask(parent, extent); + maskScalars(source).fill(1); + source.modified(); + return { parent, source }; +} + +describe('render-only segment slices', () => { + it.each([0, 1, 2])('pads only the displayed plane on axis %i', (axis) => { + const extent: Extent3D = [6, 6, 7, 7, 8, 8]; + const { parent, source } = scene(extent); + const original = maskScalars(source); + const origin = [...source.getOrigin()]; + const rendered = segmentRenderMask(source, parent, extent, { + axis: axis, + index: extent[axis * 2], + })!; + const dimensions = [3, 3, 3]; + dimensions[axis] = 1; + expect(rendered.getDimensions()).toEqual(dimensions); + const center: [number, number, number] = [1, 1, 1]; + center[axis] = 0; + expect(Array.from(rendered.indexToWorld(center))).toEqual(origin); + expect([...maskScalars(rendered)]).toEqual([0, 0, 0, 0, 1, 0, 0, 0, 0]); + expect(source.getDimensions()).toEqual([1, 1, 1]); + expect(source.getOrigin()).toEqual(origin); + expect(maskScalars(source)).toBe(original); + expect([...original]).toEqual([1]); + expect( + segmentRenderMask(source, parent, extent, { + axis: axis, + index: extent[axis * 2] + 1, + }) + ).toBeNull(); + expect( + segmentRenderMask(source, parent, [0, -1, 0, -1, 0, -1], { + axis: axis, + index: 0, + }) + ).toBeNull(); + }); + + it.each([0, 1, 2, 3, 4, 5])('clips padding at parent face %i', (face) => { + const extent: Extent3D = [6, 6, 7, 7, 8, 8]; + const bounds = [4, 9, 5, 10, 6, 11]; + const boundaryAxis = Math.floor(face / 2); + const axis = (boundaryAxis + 1) % 3; + extent[boundaryAxis * 2] = bounds[face]; + extent[boundaryAxis * 2 + 1] = bounds[face]; + const { parent, source } = scene(extent); + const rendered = segmentRenderMask(source, parent, extent, { + axis: axis, + index: extent[axis * 2], + })!; + expect(rendered.getDimensions()[boundaryAxis]).toBe(2); + const corner: [number, number, number] = [0, 0, 0]; + if (face % 2) corner[boundaryAxis] = 1; + expect( + parent.worldToIndex(rendered.indexToWorld(corner))[boundaryAxis] + ).toBeCloseTo(bounds[face]); + expect([...maskScalars(rendered)].reduce((a, b) => a + b, 0)).toBe(1); + }); + + it('refreshes after edits, slice changes and growth', () => { + const extent: Extent3D = [6, 6, 7, 7, 8, 9]; + const { parent, source } = scene(extent); + const first = segmentRenderMask(source, parent, extent, { + axis: 2, + index: 8, + })!; + expect( + segmentRenderMask(source, parent, extent, { axis: 2, index: 8 }) + ).toBe(first); + maskScalars(source)[0] = 0; + source.modified(); + expect( + segmentRenderMask(source, parent, extent, { axis: 2, index: 8 }) + ).toBe(first); + expect([...maskScalars(first)].every((v) => v === 0)).toBe(true); + expect(first.getPointData().getScalars().getRange()).toEqual([0, 0]); + const nextSlice = segmentRenderMask(source, parent, extent, { + axis: 2, + index: 9, + })!; + expect([...maskScalars(nextSlice)]).toEqual([0, 0, 0, 0, 1, 0, 0, 0, 0]); + expect( + Array.from(parent.worldToIndex(nextSlice.indexToWorld([1, 1, 0]))) + ).toEqual([6, 7, 9]); + const grown: Extent3D = [5, 6, 7, 7, 8, 9]; + regrowMask(source, parent, extent, grown); + maskScalars(source)[0] = 1; + source.modified(); + const next = segmentRenderMask(source, parent, grown, { + axis: 2, + index: 8, + })!; + expect(next.getDimensions()).toEqual([4, 3, 1]); + expect(next.getPointData().getScalars().getRange()).toEqual([0, 1]); + expect([...maskScalars(next)].reduce((a, b) => a + b, 0)).toBe(1); + }); +}); diff --git a/src/segmentation/rendering/__tests__/segmentSliceVisibility.spec.ts b/src/segmentation/rendering/__tests__/segmentSliceVisibility.spec.ts new file mode 100644 index 000000000..d9d3d39fe --- /dev/null +++ b/src/segmentation/rendering/__tests__/segmentSliceVisibility.spec.ts @@ -0,0 +1,87 @@ +import { describe, expect, it } from 'vitest'; + +import { + segmentCoincidentOffset, + sliceWithinExtent, +} from '@/src/segmentation/rendering/display'; +import { emptyExtent, type Extent3D } from '@/src/segmentation/geometry'; + +// --------------------------------------------------------------------------- +// The two view-layer rules a mask per segment needs. +// +// `sliceWithinExtent` answers whether a segment's actor has anything to draw on +// the slice being viewed. A bounded mask covers only part of the volume, and +// vtkImageMapper clamps a slice outside its input to the nearest one, so an +// actor left visible off its own extent paints a stale slice over the image. +// The slice and the extent are both in the PARENT image's index space, on the +// index axis the view's LPS axis maps to. +// +// `segmentCoincidentOffset` gives each segment its own coincident-topology +// polygon offset, by its back-to-front stack index. Overlap is +// representable, so segments sharing one offset would z-fight. +// Greater stack indices draw in front. Registry order is mapped in reverse. +// --------------------------------------------------------------------------- + +const EXTENT: Extent3D = [1, 2, 0, 3, 2, 5]; + +describe('sliceWithinExtent', () => { + it('is true for a slice inside the extent', () => { + expect(sliceWithinExtent(EXTENT, 2, 3)).toBe(true); + }); + + it('includes both ends of the extent', () => { + expect(sliceWithinExtent(EXTENT, 2, 2)).toBe(true); + expect(sliceWithinExtent(EXTENT, 2, 5)).toBe(true); + }); + + it('is false below the extent', () => { + expect(sliceWithinExtent(EXTENT, 2, 1)).toBe(false); + }); + + it('is false above the extent', () => { + expect(sliceWithinExtent(EXTENT, 2, 6)).toBe(false); + }); + + it('reads the axis it is given', () => { + expect(sliceWithinExtent(EXTENT, 0, 3)).toBe(false); + expect(sliceWithinExtent(EXTENT, 1, 3)).toBe(true); + }); + + it('is false for a mask that covers nothing, on every axis', () => { + expect(sliceWithinExtent(emptyExtent(), 0, 0)).toBe(false); + expect(sliceWithinExtent(emptyExtent(), 1, 0)).toBe(false); + expect(sliceWithinExtent(emptyExtent(), 2, 0)).toBe(false); + }); +}); + +describe('segmentCoincidentOffset', () => { + it('puts the first segment in front of the base image', () => { + const [factor, units] = segmentCoincidentOffset(0); + + expect(factor).toBeLessThan(0); + expect(units).toBeLessThan(0); + }); + + it('puts a greater stack index in front of a smaller one', () => { + const [earlierFactor, earlierUnits] = segmentCoincidentOffset(0); + const [laterFactor, laterUnits] = segmentCoincidentOffset(1); + + expect(laterUnits).toBeLessThan(earlierUnits); + expect(laterFactor).toBeLessThanOrEqual(earlierFactor); + }); + + it('keeps that order all the way down a long list', () => { + const offsets = Array.from({ length: 64 }, (_, index) => + segmentCoincidentOffset(index) + ); + + expect( + offsets.every( + ([factor, units]) => Number.isFinite(factor) && Number.isFinite(units) + ) + ).toBe(true); + offsets.slice(1).forEach(([, units], index) => { + expect(units).toBeLessThan(offsets[index][1]); + }); + }); +}); diff --git a/src/segmentation/rendering/display.ts b/src/segmentation/rendering/display.ts new file mode 100644 index 000000000..484cae0e3 --- /dev/null +++ b/src/segmentation/rendering/display.ts @@ -0,0 +1,87 @@ +import type { LabelmapSegment } from '@/src/segmentation/model'; +import type { Extent3D } from '@/src/segmentation/geometry'; +import { isEmptyExtent } from '@/src/segmentation/geometry'; + +/** + * Whether a segment's actor has anything to draw on the slice being viewed. + * Extent and slice are both in the parent image's index space, on the index + * axis the view's LPS axis maps to. vtkImageMapper clamps a slice outside its + * input to the nearest one, so an actor left visible off its own extent paints + * a stale slice over the image. + */ +export function sliceWithinExtent( + extent: Extent3D, + axisIndex: number, + slice: number +) { + if (isEmptyExtent(extent)) return false; + return slice >= extent[axisIndex * 2] && slice <= extent[axisIndex * 2 + 1]; +} + +const SEGMENT_OFFSET_FACTOR = -4; + +/** + * Actor opacity for a segment's slice representation. Per-segment and + * per-segmentation opacity live in the transfer functions, so the actor itself + * carries none. It must stay below 1: vtk.js puts an image slice in the opaque + * render pass at an opacity of 1, which restacks it against the base image and + * the sibling segment actors. + */ +export const SEGMENT_ACTOR_OPACITY = 0.9999; + +/** + * A mask's coincident-topology polygon offset, by back-to-front stack index. + * Overlapping segments need distinct offsets to avoid z-fighting. Greater + * stack indices sit closer to the viewer; the registry maps its first entry + * to the greatest index. + */ +export function segmentCoincidentOffset(stackIndex: number) { + return [SEGMENT_OFFSET_FACTOR, SEGMENT_OFFSET_FACTOR - stackIndex] as [ + number, + number, + ]; +} + +/** + * Fill alpha in 0..1 for the slice representation's piecewise function: the + * segment's own alpha scaled by its fill opacity and by the segmentation's, + * the same way the outline tables compose theirs. + */ +export const segmentFillAlpha = ( + segment: LabelmapSegment, + segmentationOpacity = 1 +) => + segment.visible + ? ((segment.color[3] || 0) / 255) * + (segment.fillOpacity ?? 1) * + segmentationOpacity + : 0; + +/** + * The label outline tables vtk.js indexes by label value minus one, so both run + * from value 1 to the largest value in use. A value no segment claims keeps the + * segmentation defaults. + */ +export const segmentOutlineTables = ( + segments: LabelmapSegment[], + segmentationThickness: number, + segmentationOpacity: number +) => { + const byValue = new Map(segments.map((segment) => [segment.value, segment])); + const largestValue = segments.reduce( + (largest, segment) => Math.max(largest, segment.value), + 0 + ); + const at = (index: number) => byValue.get(index + 1); + + return { + thicknesses: Array.from({ length: largestValue }, (_, index) => { + const segment = at(index); + return !segment || segment.visible ? segmentationThickness : 0; + }), + opacities: Array.from( + { length: largestValue }, + (_, index) => segmentationOpacity * (at(index)?.outlineOpacity ?? 1) + ), + }; +}; diff --git a/src/segmentation/rendering/projection.ts b/src/segmentation/rendering/projection.ts new file mode 100644 index 000000000..70c2d57ef --- /dev/null +++ b/src/segmentation/rendering/projection.ts @@ -0,0 +1,60 @@ +import { computed } from 'vue'; + +import type { SegmentRegistry } from '@/src/segmentation/segmentRegistry'; +import { + sameLabelmapSegment, + toLabelmapSegment, +} from '@/src/segmentation/segment'; +import { SEGMENT_VALUE } from '@/src/segmentation/masks/labelValue'; +import { + listMasks, + type LabelmapSegment, + type Segmentation, +} from '@/src/segmentation/model'; + +export type SegmentProjectionDeps = { + segmentations: Record; + segmentRegistry: SegmentRegistry; +}; + +/** One descriptor per bound mask; unchanged appearances retain object identity. */ +export function createSegmentProjection({ + segmentations, + segmentRegistry, +}: SegmentProjectionDeps) { + function project() { + const byMask: Record = {}; + Object.values(segmentations).forEach((segmentation) => { + listMasks(segmentation).forEach((segment) => { + if (!segment.representations.labelmap) return; + byMask[segment.id] = toLabelmapSegment( + segmentRegistry.getSegment(segment.segmentId), + SEGMENT_VALUE + ); + }); + }); + return byMask; + } + + let projected: Record = {}; + + return computed(() => { + const fresh = project(); + const stable = Object.fromEntries( + Object.entries(fresh).map(([maskId, list]) => { + const previous = projected[maskId]; + return [ + maskId, + previous && sameLabelmapSegment(previous, list) ? previous : list, + ]; + }) + ); + // The record's own identity is what a consumer of the whole projection + // watches, so it survives a change that left every mask alone. + const unchanged = + Object.keys(stable).length === Object.keys(projected).length && + Object.entries(stable).every(([id, list]) => projected[id] === list); + if (!unchanged) projected = stable; + return projected; + }); +} diff --git a/src/segmentation/rendering/renderMask.ts b/src/segmentation/rendering/renderMask.ts new file mode 100644 index 000000000..78a2343c4 --- /dev/null +++ b/src/segmentation/rendering/renderMask.ts @@ -0,0 +1,83 @@ +import type vtkImageData from '@kitware/vtk.js/Common/DataModel/ImageData'; +import type vtkLabelMap from '@/src/vtk/LabelMap'; +import { + allocateMask, + reframeMaskScalars, +} from '@/src/segmentation/masks/storage'; +import { SEGMENT_VALUE } from '@/src/segmentation/masks/labelValue'; +import { maskScalars } from '@/src/segmentation/model'; +import { + clipExtent, + isEmptyExtent, + padExtent, + type Extent3D, +} from '@/src/segmentation/geometry'; + +// Keep only the current slice per axis, outside segmentation storage and export. +const renderMasks = new WeakMap< + vtkLabelMap, + Map< + number, + { + image: vtkLabelMap; + key: string; + mtime: number; + } + > +>(); + +/** Add known background within the scan, without inventing data beyond it. */ +export function segmentRenderMask( + source: vtkLabelMap, + parent: vtkImageData, + extent: Extent3D, + { axis, index: slice }: { axis: number; index: number } +) { + if (isEmptyExtent(extent)) return null; + const padded = clipExtent( + padExtent(extent, 1), + parent.getExtent() as Extent3D + ); + const index = Math.round(slice); + if (index < extent[axis * 2] || index > extent[axis * 2 + 1]) return null; + padded[axis * 2] = index; + padded[axis * 2 + 1] = index; + const key = [ + ...extent, + ...padded, + ...parent.getOrigin(), + ...parent.getSpacing(), + ...parent.getDirection(), + ].join(','); + let slices = renderMasks.get(source); + if (!slices) { + slices = new Map(); + renderMasks.set(source, slices); + } + let cached = slices.get(axis); + if (!cached || cached.key !== key) { + cached = { image: allocateMask(parent, padded), key, mtime: -1 }; + slices.set(axis, cached); + } + if (cached.mtime !== source.getMTime()) { + const values = reframeMaskScalars( + maskScalars(source), + [...extent], + padded, + maskScalars(cached.image) + ); + const scalars = cached.image.getPointData().getScalars(); + scalars.dataChange(); + // Each stored mask is binary. Supply the range to avoid another scan. + scalars.setRange( + { + min: values.includes(0) ? 0 : SEGMENT_VALUE, + max: values.includes(SEGMENT_VALUE) ? SEGMENT_VALUE : 0, + }, + 0 + ); + cached.image.modified(); + cached.mtime = source.getMTime(); + } + return cached.image; +} diff --git a/src/segmentation/segment.ts b/src/segmentation/segment.ts new file mode 100644 index 000000000..68c089fc3 --- /dev/null +++ b/src/segmentation/segment.ts @@ -0,0 +1,97 @@ +import type { RGBAColor } from '@kitware/vtk.js/types'; + +import { + STROKE_WIDTH_ANNOTATION_TOOL_DEFAULT, + TOOL_COLORS, +} from '@/src/config'; +import type { Maybe } from '@/src/types'; +import { cleanUndefined } from '@/src/utils'; +import type { LabelmapSegment } from '@/src/segmentation/model'; +import { cssColorToRGBA, rgbaToCssColor } from '@/src/segmentation/color'; + +/** + * Identity and shared appearance for everything drawn as one thing: a paint + * mask on any image, a rectangle, a polygon. Appearance fields are absent + * until set and mean "app default" while they are, so a configured or imported + * type that states nothing follows the default and a ruler type carries no + * meaningless opacity. + */ +export type Segment = { + id: string; + name: string; + color: RGBAColor; + // State the user sets on the thing itself, so it holds on every image. + visible: boolean; + locked: boolean; + fillOpacity?: number; + outlineOpacity?: number; + strokeWidth?: number; +}; + +export type SegmentInit = Partial>; + +export const DEFAULT_SEGMENT_COLOR = cssColorToRGBA(TOOL_COLORS[0]); + +const APPEARANCE_DEFAULTS = { + name: '', + color: DEFAULT_SEGMENT_COLOR, + visible: true, + locked: false, + fillOpacity: 1, + outlineOpacity: 1, + strokeWidth: STROKE_WIDTH_ANNOTATION_TOOL_DEFAULT, +}; + +/** + * The one resolver. Every renderer, editor and encoder reads a type through + * it; nothing reads the optional fields directly, so an absent field means the + * app default in exactly one place. + */ +export const resolveSegmentAppearance = (type: Maybe) => { + const stated: Partial = type ?? {}; + const resolved = { + ...APPEARANCE_DEFAULTS, + ...cleanUndefined({ + name: stated.name, + color: stated.color, + visible: stated.visible, + locked: stated.locked, + fillOpacity: stated.fillOpacity, + outlineOpacity: stated.outlineOpacity, + strokeWidth: stated.strokeWidth, + }), + }; + return { ...resolved, cssColor: rgbaToCssColor(resolved.color) }; +}; + +export type SegmentAppearance = ReturnType; + +/** + * The descriptor a record projects onto the label value its mask holds, all of + * it resolved from the segment. The labelmap renderer and the .seg.nrrd writer + * consume it. + */ +export const toLabelmapSegment = ( + type: Maybe, + labelValue: number +): LabelmapSegment => { + const resolved = resolveSegmentAppearance(type); + return { + value: labelValue, + name: resolved.name, + color: [...resolved.color] as RGBAColor, + visible: resolved.visible, + locked: resolved.locked, + fillOpacity: resolved.fillOpacity, + outlineOpacity: resolved.outlineOpacity, + }; +}; + +export const sameLabelmapSegment = (a: LabelmapSegment, b: LabelmapSegment) => + a.value === b.value && + a.name === b.name && + a.visible === b.visible && + a.locked === b.locked && + a.fillOpacity === b.fillOpacity && + a.outlineOpacity === b.outlineOpacity && + a.color.every((channel, index) => channel === b.color[index]); diff --git a/src/segmentation/segmentReferences.ts b/src/segmentation/segmentReferences.ts new file mode 100644 index 000000000..62541be05 --- /dev/null +++ b/src/segmentation/segmentReferences.ts @@ -0,0 +1,40 @@ +import { getActivePinia, type Pinia } from 'pinia'; + +// Who points at a segment, declared by the stores that hold references. +// +// The registry knows nothing about masks or shapes: it asks these declarations +// whether a segment is still referenced and hands them the removal. Declarations +// are registered from store setup and scoped to the application instance that +// ran it, so two applications, or two tests, never answer for each other. +// Dependency-free on purpose: a store import here would close a cycle back +// through the registry. + +export type SegmentReferenceHolder = { + has: (segmentId: string) => boolean; + remove: (segmentId: string) => void; +}; + +const holdersByApp = new WeakMap>(); + +const holdersOf = () => { + const pinia = getActivePinia(); + if (!pinia) return undefined; + const existing = holdersByApp.get(pinia); + if (existing) return existing; + const holders = new Map(); + holdersByApp.set(pinia, holders); + return holders; +}; + +export function declareSegmentReferences( + name: string, + holder: SegmentReferenceHolder +) { + holdersOf()?.set(name, holder); +} + +export const segmentIsReferenced = (segmentId: string) => + [...(holdersOf()?.values() ?? [])].some((holder) => holder.has(segmentId)); + +export const removeSegmentReferences = (segmentId: string) => + holdersOf()?.forEach((holder) => holder.remove(segmentId)); diff --git a/src/segmentation/segmentRegistry.ts b/src/segmentation/segmentRegistry.ts new file mode 100644 index 000000000..14ef5b2c6 --- /dev/null +++ b/src/segmentation/segmentRegistry.ts @@ -0,0 +1,276 @@ +import { computed, ref, type Ref } from 'vue'; + +import { TOOL_COLORS } from '@/src/config'; +import { useIdStore } from '@/src/store/id'; +import type { Maybe } from '@/src/types'; +import { cssColorToRGBA } from '@/src/segmentation/color'; +import { + resolveSegmentAppearance, + type Segment, + type SegmentInit, +} from '@/src/segmentation/segment'; +import { omit } from '@/src/utils'; +import { cleanUndefined } from '@/src/utils'; + +export type ConfiguredSegment = { + color?: string; + fillOpacity?: number; + outlineOpacity?: number; + strokeWidth?: number; +}; + +export type ConfiguredSegments = Record; + +export type SegmentRegistryOptions = { + hasReferences?: (segmentId: string) => boolean; + removeReferences?: (segmentId: string) => void; +}; + +const fromConfigured = ( + name: string, + configured: ConfiguredSegment +): SegmentInit => + cleanUndefined({ + name, + color: configured.color ? cssColorToRGBA(configured.color) : undefined, + fillOpacity: configured.fillOpacity, + outlineOpacity: configured.outlineOpacity, + strokeWidth: configured.strokeWidth, + }); + +const configuredAppearance = ({ + color, + fillOpacity, + outlineOpacity, + strokeWidth, +}: Segment) => ({ color, fillOpacity, outlineOpacity, strokeWidth }); + +/** + * Identity and shared appearance for a family of segments: one instance backs + * paint, rectangles, polygons and rulers together. Explicit order drives the + * picker, shortcuts, serialization and labelmap stacking. + */ +export const createSegmentRegistry = ({ + hasReferences = () => false, + removeReferences = () => {}, +}: SegmentRegistryOptions = {}) => { + const segmentById = ref>({}) as Ref< + Record + >; + + const segmentOrder = ref([]); + const segmentList = computed(() => + segmentOrder.value.map((id) => segmentById.value[id]) + ); + + const selectedSegmentId = ref>(); + const selectionRevision = ref(0); + + // A type that is gone is not selected. + const selectedSegment = computed(() => + selectedSegmentId.value + ? segmentById.value[selectedSegmentId.value] + : undefined + ); + + const selectSegment = (id: Maybe) => { + selectedSegmentId.value = id && segmentById.value[id] ? id : undefined; + // Reselecting the same segment can still request that its row be revealed. + if (selectedSegmentId.value) selectionRevision.value += 1; + }; + + const getSegment = (id: Maybe) => + id ? segmentById.value[id] : undefined; + + const appearanceOf = (id: Maybe) => + resolveSegmentAppearance(getSegment(id)); + + // Cached: the renderer asks for one index per mask per re-render, and the + // export sort asks twice per comparison. + const orderIndex = computed( + () => new Map(segmentOrder.value.map((id, index) => [id, index])) + ); + + const orderIndexOf = (id: Maybe) => + (id === undefined || id === null ? undefined : orderIndex.value.get(id)) ?? + -1; + + const findSegmentByName = (name: Maybe) => + segmentList.value.find((type) => type.name === name); + + const uniqueName = (stem: string) => { + const taken = new Set(segmentList.value.map((type) => type.name.trim())); + if (!taken.has(stem)) return stem; + let index = 2; + while (taken.has(`${stem} (${index})`)) index += 1; + return `${stem} (${index})`; + }; + + const defaultName = () => { + const taken = new Set(segmentList.value.map((type) => type.name.trim())); + let index = 1; + while (taken.has(`Segment ${index}`)) index += 1; + return `Segment ${index}`; + }; + + let nextColorIndex = 0; + const nextColor = () => { + const color = cssColorToRGBA(TOOL_COLORS[nextColorIndex]); + nextColorIndex = (nextColorIndex + 1) % TOOL_COLORS.length; + return color; + }; + + /** Mints a segment without touching the selection. Allocates no voxels. */ + const mintSegment = (init: SegmentInit = {}) => { + const id = useIdStore().nextId(); + segmentById.value = { + ...segmentById.value, + [id]: { + name: defaultName(), + color: nextColor(), + visible: true, + locked: false, + ...cleanUndefined(init), + id, + }, + }; + segmentOrder.value = [...segmentOrder.value, id]; + return id; + }; + + const addSegment = (init: SegmentInit = {}) => { + const id = mintSegment(init); + selectSegment(id); + return id; + }; + + const updateSegment = (id: string, patch: SegmentInit) => { + const type = segmentById.value[id]; + if (!type) return; + segmentById.value = { + ...segmentById.value, + [id]: { ...type, ...patch, id }, + }; + }; + + // Deleting a referenced segment takes its masks and shapes with it; the + // caller owns the confirmation. + const deleteSegment = (id: string) => { + if (!segmentById.value[id]) return; + removeReferences(id); + segmentOrder.value = segmentOrder.value.filter((key) => key !== id); + segmentById.value = omit(segmentById.value, id); + if (selectedSegmentId.value === id) { + selectSegment(segmentOrder.value[0]); + } + }; + + const moveSegment = (id: string, target: string, after = false) => { + if (id === target || !getSegment(id) || !getSegment(target)) return; + const order = segmentOrder.value.filter((key) => key !== id); + order.splice(order.indexOf(target) + Number(after), 0, id); + segmentOrder.value = order; + }; + + const ensureSelectedSegment = () => { + if (selectedSegment.value) return selectedSegment.value.id; + const first = segmentList.value[0]; + if (!first) return addSegment(); + selectSegment(first.id); + return first.id; + }; + + /** + * Exact-name lookup, minting on a miss. Import binds descriptors this way, + * so a file's segment lands on the one already carrying that name and the + * registry's own color wins. + */ + const segmentNamed = (name: string, init: SegmentInit = {}) => { + const existing = findSegmentByName(name); + if (existing) return existing.id; + return mintSegment({ ...init, name }); + }; + + // --- config overlay --- // + + // Keep each key's identity and appearance beneath its config contribution. + // New segments begin with automatic color and default optional appearance; + // a restored segment begins with its session appearance. + const configEntries = new Map< + string, + { id: string; appearance: ReturnType } + >(); + + const replaceConfigSegments = (configured: Maybe) => { + const next = configured ?? {}; + + Object.entries(next).forEach(([name, props]) => { + let entry = configEntries.get(name); + if (!entry || !getSegment(entry.id)) { + const id = findSegmentByName(name)?.id ?? mintSegment({ name }); + entry = { id, appearance: configuredAppearance(getSegment(id)!) }; + } + updateSegment(entry.id, { + ...entry.appearance, + ...fromConfigured(name, props), + }); + configEntries.set(name, entry); + }); + + [...configEntries.entries()] + .filter(([name]) => !(name in next)) + .forEach(([name, { id }]) => { + configEntries.delete(name); + // Content keeps the last configured appearance as session state. + if (segmentById.value[id] && !hasReferences(id)) deleteSegment(id); + }); + + // A configured registry offers a selection from the start; selecting + // creates nothing, so the first edit lands in a configured type rather + // than minting one beside it. + if (!selectedSegment.value) selectSegment(segmentList.value[0]?.id); + }; + + // --- wire --- // + + const serialize = () => segmentList.value.map((type) => ({ ...type })); + + /** + * Seats restored segments beside the ones already here. Ids are minted fresh + * and every incoming reference is remapped through the returned map, so an + * import into a populated scene overwrites nothing. + */ + const adopt = (incoming: Maybe) => + Object.fromEntries( + (incoming ?? []).map(({ id, ...init }) => [ + id, + mintSegment(init as SegmentInit), // side effect in Array.map + ]) + ); + + return { + segmentById, + segmentList, + selectedSegmentId, + selectionRevision, + selectedSegment, + selectSegment, + getSegment, + appearanceOf, + orderIndexOf, + findSegmentByName, + uniqueName, + mintSegment, + addSegment, + updateSegment, + moveSegment, + deleteSegment, + ensureSelectedSegment, + segmentNamed, + replaceConfigSegments, + serialize, + adopt, + }; +}; + +export type SegmentRegistry = ReturnType; diff --git a/src/utils/color.ts b/src/utils/color.ts index e71f644ea..80ca63bba 100644 --- a/src/utils/color.ts +++ b/src/utils/color.ts @@ -1,5 +1,15 @@ import type { RGBAColor } from '@kitware/vtk.js/types'; +/** A cursor over a palette: each call hands out the next opaque colour. */ +export function cycleColors(palette: readonly (readonly number[])[]) { + let index = 0; + return () => { + const color = palette[index]; + index = (index + 1) % palette.length; + return [...color, 255] as RGBAColor; + }; +} + /** * Converts an RGBA tuple to a hex string with alpha. * From e8a21494b086c215c00d26a6449246fc9c31fe24 Mon Sep 17 00:00:00 2001 From: Paul Elliott Date: Sat, 19 Sep 2026 15:06:41 +0000 Subject: [PATCH 2/3] fix(resample): return a same-grid labelmap without reslicing The label path of ensureSameSpace only skipped vtkImageReslice when the product of the source worldToIndex and the target indexToWorld was a bit-exact identity matrix. Real geometry never satisfies that: with an origin and spacing such as [-120.1, -98.4, 33.7] and [0.7, 0.7, 3] the product carries float error of about 1e-5, so a labelmap that already sat on its parent's grid was copied through a full nearest-neighbour reslice on the main thread, once per imported or restored mask. The axis and eight-corner checks just above already prove the two grids coincide within the function's 1e-3 tolerance, so the identity test now only asks that the mapping neither permutes nor flips an axis, and the source image is returned as-is in that case. Permuted, flipped and genuinely different grids keep going through the reslice. --- .../__tests__/reorientLabelImage.spec.ts | 26 +++++++++++++++++++ src/io/resample/reorientLabelImage.ts | 6 ++++- 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/io/resample/__tests__/reorientLabelImage.spec.ts b/src/io/resample/__tests__/reorientLabelImage.spec.ts index 16b90b43a..df79b6a03 100644 --- a/src/io/resample/__tests__/reorientLabelImage.spec.ts +++ b/src/io/resample/__tests__/reorientLabelImage.spec.ts @@ -89,6 +89,32 @@ describe('label-grid reorientation', () => { ]).toEqual([...source.getPointData().getScalars().getData()]); }); + it('returns the source itself when it already sits on the target grid', () => { + const spacing: [number, number, number] = [0.7, 0.7, 3]; + const origin: [number, number, number] = [-120.1, -98.4, 33.7]; + const source = image([3, 4, 5]); + const target = image([3, 4, 5]); + [source, target].forEach((im) => { + im.setSpacing(spacing); + im.setOrigin(origin); + }); + expect(reorientLabelImage(target, source)).toBe(source); + + // The same geometry laid out along a flipped axis is a different grid and + // still has to go through the reslice. + const flipped = image([3, 4, 5]); + flipped.setSpacing(spacing); + flipped.setOrigin([origin[0] + spacing[0] * 2, origin[1], origin[2]]); + flipped.setDirection([-1, 0, 0, 0, 1, 0, 0, 0, 1]); + const output = reorientLabelImage(target, flipped)!; + expect(output).not.toBe(flipped); + const values = output.getPointData().getScalars().getData(); + for (let k = 0; k < 5; k++) + for (let j = 0; j < 4; j++) + for (let i = 0; i < 3; i++) + expect(values[i + 3 * (j + 4 * k)]).toBe(1 + (2 - i) + 3 * (j + 4 * k)); + }); + it('defers fractional shifts, different sampling, and cropping to interpolation', () => { const source = image([3, 4, 5]); const target = image([3, 4, 5]); diff --git a/src/io/resample/reorientLabelImage.ts b/src/io/resample/reorientLabelImage.ts index bb9d85922..c7ef31375 100644 --- a/src/io/resample/reorientLabelImage.ts +++ b/src/io/resample/reorientLabelImage.ts @@ -62,7 +62,11 @@ export function reorientLabelImage(target: vtkImageData, source: vtkImageData) { return null; } - if (matrix.every((value, index) => value === (index % 5 === 0 ? 1 : 0))) + // The corner check just proved both grids coincide within `tolerance`, so an + // unpermuted, unflipped mapping means the source already sits on the target + // grid. Demanding a bit-exact identity here instead would reslice every real + // image, whose transforms never multiply back to exactly one. + if (axes.every((axis, column) => axis === column && matrix[5 * column] > 0)) return source; const filter = vtkImageReslice.newInstance(); From a4ba32f044c35e4ad8a9ed23d60d7197c5d64d03 Mon Sep 17 00:00:00 2001 From: Paul Elliott Date: Sun, 20 Sep 2026 18:46:20 +0000 Subject: [PATCH 3/3] perf(segmentation): stop allocating a closure per claimed voxel masksHolding and masksClearing.claim swept the reaching masks with some/forEach, so every voxel a brush or polygon fill walked allocated a fresh callback capturing i, j and k. Both now sweep the reaching masks with a plain indexed loop, so the per-voxel path allocates nothing. The answers are unchanged. --- src/segmentation/masks/overlap.ts | 36 +++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/src/segmentation/masks/overlap.ts b/src/segmentation/masks/overlap.ts index 5be07d464..4a87514c5 100644 --- a/src/segmentation/masks/overlap.ts +++ b/src/segmentation/masks/overlap.ts @@ -45,17 +45,26 @@ const masksReaching = (masks: BoundedScalars[], within: Extent3D) => /** * Whether any of these masks holds the voxel at PARENT indices i, j, k, over * the box the caller is about to walk. Absent when no mask reaches that box. + * + * The answer is asked once per voxel the caller walks, so the sweep over the + * reaching masks is a plain indexed loop: a callback taking i, j, k would be a + * fresh closure per voxel. */ export function masksHolding(masks: BoundedScalars[], within: Extent3D) { const reaching = masksReaching(masks, within); if (reaching.length === 0) return undefined; - return (i: number, j: number, k: number) => - reaching.some( - (bounded) => + return (i: number, j: number, k: number) => { + for (let index = 0; index < reaching.length; index += 1) { + const bounded = reaching[index]; + if ( extentContainsIndex(bounded.extent, i, j, k) && bounded.scalars[maskOffset(bounded, i, j, k)] !== LABELMAP_BACKGROUND_VALUE - ); + ) + return true; + } + return false; + }; } /** @@ -70,14 +79,19 @@ export function masksClearing(masks: BoundedScalars[], within: Extent3D) { const reaching = masksReaching(masks, within); if (reaching.length === 0) return undefined; const changed = new Set(); + // Indexed loop, as in masksHolding: claim runs once per voxel the caller + // walks, and a callback over the reaching masks would allocate per voxel. const claim = (i: number, j: number, k: number) => { - reaching.forEach((bounded) => { - if (!extentContainsIndex(bounded.extent, i, j, k)) return; - const offset = maskOffset(bounded, i, j, k); - if (bounded.scalars[offset] === LABELMAP_BACKGROUND_VALUE) return; - bounded.scalars[offset] = LABELMAP_BACKGROUND_VALUE; - changed.add(bounded.mask); - }); + for (let index = 0; index < reaching.length; index += 1) { + const bounded = reaching[index]; + if (extentContainsIndex(bounded.extent, i, j, k)) { + const offset = maskOffset(bounded, i, j, k); + if (bounded.scalars[offset] !== LABELMAP_BACKGROUND_VALUE) { + bounded.scalars[offset] = LABELMAP_BACKGROUND_VALUE; + changed.add(bounded.mask); + } + } + } return true; }; return {