Highlight occupied depth in en face views

eslint.config.js

@@ -313,12 +313,14 @@ export default tseslint.config(
     {
       dir: 'enface',
       publicSurface: false,
-      pure: ['src/enface/sliceGeometry.ts'],
+      pure: ['src/enface/thickness.ts', 'src/enface/sliceGeometry.ts'],
     },
     {
       dir: 'segmentation',
       pure: [
         'src/segmentation/geometry.ts',
+        'src/segmentation/measurements/occupiedDepth.ts',
+        'src/segmentation/measurements/occupiedDepth.worker.ts',
         'src/segmentation/color.ts',
         'src/segmentation/model.ts',
         'src/segmentation/segment.ts',

scripts/samples/jhu-oct-provenance.json

@@ -93,7 +93,7 @@
     "axial_spacing_microns": 3.867251332849264,
     "lateral_sampling": "Retain native lateral indices 0,2,...,1022; native index 1023 is omitted. ScaleX is doubled.",
     "slow_sampling": "Linearly resample transformed intensities and unrounded boundaries from per-B-scan physical positions to 49 evenly spaced positions across 6.149238953366875 mm. Native adjacent pitches range from 0.12010203066263614 to 0.13211284885730468 mm.",
-    "axial_sampling": "Native depth samples and ScaleZ are preserved at 3.867251332849264 microns per voxel.",
+    "axial_sampling": "Native depth samples and ScaleZ are preserved; thickness is occupied depth voxels multiplied by 3.867251332849264 microns.",
     "coordinates": "Shared local zero origin, identity direction, and millimetre spacing for image and segmentation. Patient anatomical orientation is not inferred.",
     "intensities": "Nested float32 square roots of source reflectivity before slow-axis interpolation, then nearest-even rounding of 255 times the regularized values to uint8. Stored values range from 0 to 255; divide by 255 for normalized display intensity. This is a display transform, not unmodified linear reflectivity.",
     "missing_reflectivity": "Non-finite or out-of-range native values are missing. A resampled voxel is valid only when every source with nonzero interpolation weight is valid. Missing output intensity is zero.",

src/core/slab/geometry.ts

@@ -21,6 +21,11 @@ export function deriveSlabGeometry(
   spec: SlabSpecification
 ) {
   const axes = [0, 1, 2].filter((axis) => axis !== spec.axis);
+  const dimensions = image.getDimensions();
+  const sourceSpacing = image.getSpacing();
+  const width = dimensions[axes[0]];
+  const height = dimensions[axes[1]];
+  const spacing = axes.map((axis) => Math.abs(sourceSpacing[axis]));
   const planeIndex: Vector3 = [0, 0, 0];
   planeIndex[spec.axis] = (spec.first + spec.last) / 2;
   const planeOrigin = image.indexToWorld(planeIndex) as Vector3;
@@ -48,6 +53,9 @@ export function deriveSlabGeometry(
     ...metadataFromVtkImageData(image),
   };
   return {
+    width,
+    height,
+    spacing,
     metadata,
     pixelToWorld,
     planeOrigin,

src/core/streaming/dicom/frameGeometry.ts

@@ -0,0 +1,70 @@
+import { vec3 } from 'gl-matrix';
+import { areEquals } from '@kitware/vtk.js/Common/Core/Math';
+
+export const near = (a: number, b: number, tolerance: number) =>
+  Math.abs(a - b) <= tolerance;
+
+export function imageDirection(basis: number[] | null | undefined) {
+  if (!basis || basis.length !== 6 || !basis.every(Number.isFinite))
+    return null;
+  const row = basis.slice(0, 3) as vec3;
+  const column = basis.slice(3) as vec3;
+  if (
+    !near(vec3.length(row), 1, 1e-4) ||
+    !near(vec3.length(column), 1, 1e-4) ||
+    // Rounding two unit vectors to four decimals can shift their dot by 1.8e-4.
+    !near(vec3.dot(row, column), 0, 2e-4)
+  )
+    return null;
+  return [...basis, ...vec3.cross([0, 0, 0], row, column)];
+}
+
+/** Validate equally spaced slices along the positive image normal. */
+export function frameGridGeometry(
+  positions: number[][],
+  direction: number[] | null
+) {
+  if (positions.some((p) => p.length !== 3 || !p.every(Number.isFinite)))
+    return {
+      reason:
+        'DICOM slice positions are invalid; spatial reconstruction is required.',
+    };
+  if (positions.length < 2) return { reason: null };
+  const step = positions
+    .at(-1)!
+    .map((n, axis) => (n - positions[0][axis]) / (positions.length - 1));
+  const distance = Math.hypot(...step);
+  // Patient positions may be printed to four decimal places in millimetres.
+  const tolerance = Math.max(1e-4, distance * 1e-3);
+  if (distance <= 1e-5)
+    return {
+      reason: 'Repeated DICOM slice locations need spatial reconstruction.',
+    };
+  if (
+    direction &&
+    !areEquals(
+      step,
+      direction.slice(6).map((n) => n * distance),
+      tolerance
+    )
+  )
+    return {
+      reason:
+        'DICOM slice positions do not follow the image grid; spatial reconstruction is required.',
+    };
+  if (
+    positions.some(
+      (position, i) =>
+        !areEquals(
+          position,
+          positions[0].map((n, axis) => n + i * step[axis]),
+          tolerance
+        )
+    )
+  )
+    return {
+      reason:
+        'DICOM slices have nonuniform positions; spatial reconstruction is required.',
+    };
+  return { reason: null, distance };
+}

src/core/streaming/dicom/octMetadata.ts

@@ -1,74 +1,25 @@
 import dicomParser, { type DataSet } from 'dicom-parser';
 import { inflateSync } from 'fflate';
-import { vec3 } from 'gl-matrix';
 import { areEquals } from '@kitware/vtk.js/Common/Core/Math';
 import { Tags, toDicomParserKey } from '@/src/core/dicomTags';
+import {
+  frameGridGeometry,
+  imageDirection,
+} from '@/src/core/streaming/dicom/frameGeometry';
 import {
   isOCTMetadata,
   OCT_GEOMETRY_ERROR,
   OCT_SCAN_PATTERN,
 } from '@/src/core/streaming/dicom/octMetadataTags';
 import type { Maybe } from '@/src/types';
 
-const close = (a: number, b: number, tolerance: number) =>
-  Number.isFinite(a) && Number.isFinite(b) && Math.abs(a - b) <= tolerance;
 const finiteVector = (
   value: number[] | null | undefined,
   length: number
 ): value is number[] =>
   value?.length === length && value.every(Number.isFinite);
 const positiveNumber = (value: number) => Number.isFinite(value) && value > 0;
 
-function imageDirection(basis: number[] | null | undefined) {
-  if (!finiteVector(basis, 6)) return null;
-  const row = basis.slice(0, 3) as vec3;
-  const column = basis.slice(3) as vec3;
-  if (
-    !close(vec3.length(row), 1, 1e-4) ||
-    !close(vec3.length(column), 1, 1e-4) ||
-    // Rounding two unit vectors to four decimals can shift their dot by 1.8e-4.
-    !close(vec3.dot(row, column), 0, 2e-4)
-  )
-    return null;
-  return [...basis, ...vec3.cross([0, 0, 0], row, column)];
-}
-
-/** Validate equally spaced slices along the positive image normal. */
-function frameGridError(positions: number[][], direction: number[] | null) {
-  if (positions.some((p) => p.length !== 3 || !p.every(Number.isFinite)))
-    return 'DICOM slice positions are invalid; spatial reconstruction is required.';
-  if (positions.length < 2) return null;
-  const step = positions
-    .at(-1)!
-    .map((n, axis) => (n - positions[0][axis]) / (positions.length - 1));
-  const distance = Math.hypot(...step);
-  // Patient positions may be printed to four decimal places in millimetres.
-  const tolerance = Math.max(1e-4, distance * 1e-3);
-  if (distance <= 1e-5)
-    return 'Repeated DICOM slice locations need spatial reconstruction.';
-  if (
-    direction &&
-    !areEquals(
-      step,
-      direction.slice(6).map((n) => n * distance),
-      tolerance
-    )
-  )
-    return 'DICOM slice positions do not follow the image grid; spatial reconstruction is required.';
-  if (
-    positions.some(
-      (position, i) =>
-        !areEquals(
-          position,
-          positions[0].map((n, axis) => n + i * step[axis]),
-          tolerance
-        )
-    )
-  )
-    return 'DICOM slices have nonuniform positions; spatial reconstruction is required.';
-  return null;
-}
-
 const SHARED = 'x52009229';
 const PER_FRAME = 'x52009230';
 const MEASURES = 'x00289110';
@@ -161,7 +112,7 @@ function positionError(geometry: ReturnType<typeof frameGeometry>) {
     !positiveNumber(Number(geometry.sliceSpacing))
   )
     return 'OCT slice spacing is unavailable; spatial reconstruction is required.';
-  return frameGridError(positions, direction);
+  return frameGridGeometry(positions, direction).reason;
 }
 
 function geometryError(geometry: ReturnType<typeof frameGeometry>) {

src/core/streaming/dicom/seriesGeometry.ts

@@ -0,0 +1,142 @@
+import { areEquals } from '@kitware/vtk.js/Common/Core/Math';
+import { Tags } from '@/src/core/dicomTags';
+import {
+  frameGridGeometry,
+  imageDirection,
+  near,
+} from '@/src/core/streaming/dicom/frameGeometry';
+import { OCT_GEOMETRY_ERROR } from '@/src/core/streaming/dicom/octMetadataTags';
+
+const TOLERANCE = 1e-5;
+
+type Grid = {
+  spacing: ArrayLike<number>;
+  direction: ArrayLike<number>;
+  origin: ArrayLike<number>;
+};
+
+function vector(tags: Map<string, string>, tag: string, length: number) {
+  const result = tags.get(tag)?.split('\\').map(Number);
+  return result?.length === length && result.every(Number.isFinite)
+    ? result
+    : null;
+}
+
+function declaredSeriesReason(tags: Map<string, string>[]) {
+  const reason = tags
+    .map((entry) => entry.get(OCT_GEOMETRY_ERROR))
+    .find(Boolean);
+  if (reason) return reason;
+  const sop = tags[0].get(Tags.SOPClassUID)?.trim();
+  if (tags.some((entry) => entry.get(Tags.SOPClassUID)?.trim() !== sop))
+    return 'The DICOM files have incompatible image types; spatial reconstruction is required.';
+  const spacing = vector(tags[0], Tags.PixelSpacing, 2);
+  if (
+    tags.some((entry) => {
+      const other = vector(entry, Tags.PixelSpacing, 2);
+      return spacing && other
+        ? !areEquals(spacing, other, TOLERANCE)
+        : spacing !== other;
+    })
+  )
+    return 'The DICOM files have inconsistent pixel spacing; spatial reconstruction is required.';
+  return null;
+}
+
+function seriesPositions(tags: Map<string, string>[]) {
+  const positions = tags.map((entry) =>
+    vector(entry, Tags.ImagePositionPatient, 3)
+  );
+  return positions.every((position) => position !== null) ? positions : null;
+}
+
+function loadedGridReason(
+  grid: Grid | undefined,
+  direction: number[],
+  origin: number[],
+  distance: number
+) {
+  if (
+    grid &&
+    (!areEquals(Array.from(grid.direction), direction, TOLERANCE) ||
+      !areEquals(Array.from(grid.origin), origin, TOLERANCE) ||
+      !near(grid.spacing[2], distance, TOLERANCE))
+  )
+    return 'The loaded DICOM grid does not match its file geometry; spatial reconstruction is required.';
+  return null;
+}
+
+function seriesGeometry(tags: Map<string, string>[], grid?: Grid) {
+  const declaredReason = declaredSeriesReason(tags);
+  if (declaredReason) return { reason: declaredReason };
+  const directions = tags.map((entry) =>
+    imageDirection(vector(entry, Tags.ImageOrientationPatient, 6))
+  );
+  const [direction] = directions;
+  if (
+    !direction ||
+    directions.some(
+      (other) => !other || !areEquals(direction, other, TOLERANCE)
+    )
+  )
+    return {
+      reason:
+        'The DICOM files have missing or inconsistent image orientation; spatial reconstruction is required.',
+    };
+  const positions = seriesPositions(tags);
+  if (!positions)
+    return {
+      reason:
+        'The DICOM files have incomplete or invalid slice positions; spatial reconstruction is required.',
+    };
+  const geometry = frameGridGeometry(positions, direction);
+  if (geometry.reason) return geometry;
+  const distance = geometry.distance!;
+  const reason = loadedGridReason(grid, direction, positions[0], distance);
+  return { reason, distance };
+}
+
+function declaredPixelSpacing(tags: Map<string, string>, component: number) {
+  const value = vector(tags, Tags.PixelSpacing, 2)?.[component] ?? NaN;
+  return value > 0 ? value : NaN;
+}
+
+function calibratedAxes(
+  tags: Map<string, string>[],
+  grid: Grid | undefined,
+  distance: number | undefined
+) {
+  const declared = [
+    declaredPixelSpacing(tags[0], 1),
+    declaredPixelSpacing(tags[0], 0),
+    tags.length > 1 ? distance : Number(tags[0].get(Tags.SpacingBetweenSlices)),
+  ];
+  return declared.map((spacing, axis) => {
+    if (!spacing || !Number.isFinite(spacing) || spacing <= 0) return false;
+    if (grid && !near(Math.abs(grid.spacing[axis]), spacing, TOLERANCE))
+      return false;
+    return (
+      axis === 2 ||
+      tags.every((entry) =>
+        near(declaredPixelSpacing(entry, 1 - axis), spacing, TOLERANCE)
+      )
+    );
+  });
+}
+
+/** Validate source slice arrangement against the loaded grid. */
+export function dicomSeriesGeometry(tags: Map<string, string>[], grid?: Grid) {
+  if (!tags.length)
+    return { reason: null, calibratedAxes: [false, false, false] };
+  const geometry =
+    tags.length > 1
+      ? seriesGeometry(tags, grid)
+      : {
+          reason: tags[0].get(OCT_GEOMETRY_ERROR) ?? null,
+          distance: undefined,
+        };
+  return {
+    ...geometry,
+    calibratedAxes: calibratedAxes(tags, grid, geometry.distance),
+  };
+}

src/enface/EnFaceControlPanel.vue

@@ -31,6 +31,15 @@ const controlState = computed(() => {
     axes: [],
     depthMaximum: 0,
     projectionReason: reason,
+    thickness: {
+      thresholdMaximum: 1,
+      segmentationReason: reason,
+      highlightReason: reason,
+      thicknessReason: reason,
+      segments: [],
+      segmentId: null,
+      highlightThin: false,
+    },
   } satisfies EnFaceViewControlState;
 });
 function patch(patchSettings: Partial<EnFaceViewConfig>) {

src/enface/EnFaceControls.vue

@@ -3,6 +3,7 @@ import { ref, watch } from 'vue';
 import ReasonedAction from '@/src/components/ReasonedAction.vue';
 import { rangeThumbLabels } from '@/src/utils/rangeThumbLabels';
 import { SLAB_METHODS } from '@/src/core/slab/types';
+import EnFaceThicknessControls from './EnFaceThicknessControls.vue';
 import type { EnFaceViewConfig, EnFaceViewControlState } from './store';
 
 const props = defineProps<
@@ -11,7 +12,7 @@ const props = defineProps<
 const emit = defineEmits<{
   patch: [settings: Partial<EnFaceViewConfig>];
 }>();
-const expandedSections = ref([0]);
+const expandedSections = ref([0, 1]);
 const slabDraft = ref<number[] | null>(null);
 watch(
   () => [
@@ -152,6 +153,12 @@ function updateAxis(axis: number) {
             </ReasonedAction>
           </v-expansion-panel-text>
         </v-expansion-panel>
+        <EnFaceThicknessControls
+          v-bind="thickness"
+          :settings="settings"
+          :embedded="embedded"
+          @patch="emit('patch', $event)"
+        />
       </v-expansion-panels>
     </v-card-text>
   </v-card>