Feature/segment histogram management

src/image/histogram.js

@@ -0,0 +1,165 @@
+/**
+ * Histogram utilities.
+ */
+
+/**
+ * Percentage of a histogram that falls in one interval.
+ *
+ * @typedef {object} HistogramIntervalPercentage
+ * @property {string} name Interval name.
+ * @property {number} percentage Percentage of total counts in [0, 100].
+ */
+
+/**
+ * Per-volume histogram (for example from
+ * {@link MaskImage#getHistogramBySegment}).
+ *
+ * @typedef {object} VolumeHistogram
+ * @property {number} volume Volume index.
+ * @property {number[]} histogram Sparse/dense counts indexed by intensity.
+ */
+
+/**
+ * Per-volume interval percentages.
+ *
+ * @typedef {object} VolumeHistogramIntervalPercentages
+ * @property {number} volume Volume index.
+ * @property {HistogramIntervalPercentage[]} intervals Percentages in
+ *   input order.
+ */
+
+/**
+ * Sum counts for intensities in the semi-open range [min, max).
+ *
+ * @param {number[]} histogram Intensity-indexed counts.
+ * @param {number} min Inclusive lower bound.
+ * @param {number} max Exclusive upper bound.
+ * @returns {number} Sum of counts in the range.
+ */
+function sumCountsInRange(histogram, min, max) {
+  let sum = 0;
+  for (let v = min; v < max; ++v) {
+    sum += histogram[v] || 0;
+  }
+  return sum;
+}
+
+/**
+ * Sum all counts in a histogram.
+ *
+ * @param {number[]} histogram Intensity-indexed counts.
+ * @returns {number} Total count.
+ */
+function sumHistogramCounts(histogram) {
+  let total = 0;
+  for (let i = 0; i < histogram.length; ++i) {
+    total += histogram[i] || 0;
+  }
+  return total;
+}
+
+/**
+ * Get the inclusive min intensity and exclusive max intensity present in a
+ * histogram (indices that are defined, including explicit zeros).
+ *
+ * @param {number[]} histogram Intensity-indexed counts.
+ * @returns {{min: number, maxExclusive: number}|undefined} Range, or
+ *   undefined if the histogram has no defined bins.
+ */
+function getHistogramIntensityRange(histogram) {
+  let min;
+  let max;
+  for (let i = 0; i < histogram.length; ++i) {
+    if (typeof histogram[i] !== 'undefined') {
+      if (typeof min === 'undefined') {
+        min = i;
+      }
+      max = i;
+    }
+  }
+  if (typeof min === 'undefined') {
+    return undefined;
+  }
+  return {min, maxExclusive: max + 1};
+}
+
+/**
+ * Build semi-open intervals [min, t0), [t0, t1), ..., [tLast, max) from
+ * interior thresholds and the histogram intensity range.
+ *
+ * @param {number} min Inclusive histogram minimum intensity.
+ * @param {number} maxExclusive Exclusive histogram maximum intensity.
+ * @param {number[]} thresholds Interior cut points (ascending).
+ * @param {string[]} names One name per resulting interval
+ *   (`thresholds.length + 1`).
+ * @returns {{name: string, min: number, max: number}[]} Named ranges.
+ */
+function buildIntervalsFromThresholds(min, maxExclusive, thresholds, names) {
+  const bounds = [min, ...thresholds, maxExclusive];
+  const intervals = [];
+  for (let i = 0; i < names.length; ++i) {
+    intervals.push({
+      name: names[i],
+      min: bounds[i],
+      max: bounds[i + 1]
+    });
+  }
+  return intervals;
+}
+
+/**
+ * Compute, per volume, the percentage of histogram counts that fall in each
+ * intensity interval.
+ *
+ * The caller only passes interior cut points (e.g. `[1, 5]`). For each volume
+ * the code builds semi-open intervals covering the histogram range:
+ * `[min, 1)`, `[1, 5)`, `[5, max)` where `min` / `max` come from that
+ * volume's histogram. Interval names are provided separately and must have
+ * length `thresholds.length + 1`.
+ *
+ * @param {VolumeHistogram[]} volumeHistograms Per-volume histograms (for
+ *   example the output of {@link MaskImage#getHistogramBySegment}).
+ * @param {number[]} thresholds Interior intensity cut points in ascending
+ *   order (e.g. `[1, 5]`).
+ * @param {string[]} names Names for the resulting intervals
+ *   (`thresholds.length + 1` entries).
+ * @returns {VolumeHistogramIntervalPercentages[]} One entry per volume, with
+ *   interval percentages in the same order as `names`.
+ */
+export function getHistogramIntervalPercentages(
+  volumeHistograms, thresholds, names) {
+  if (names.length !== thresholds.length + 1) {
+    throw new Error(
+      'getHistogramIntervalPercentages: names length must be ' +
+      'thresholds.length + 1.');
+  }
+
+  return volumeHistograms.map((volumeHistogram) => {
+    const {volume, histogram} = volumeHistogram;
+    const total = sumHistogramCounts(histogram);
+
+    if (total === 0) {
+      return {
+        volume,
+        intervals: names.map((name) => ({name, percentage: 0}))
+      };
+    }
+
+    const range = getHistogramIntensityRange(histogram);
+    const built = buildIntervalsFromThresholds(
+      range.min, range.maxExclusive, thresholds, names);
+
+    const intervalResults = built.map((interval) => {
+      const count = sumCountsInRange(histogram, interval.min, interval.max);
+      return {
+        name: interval.name,
+        percentage: (count * 100) / total
+      };
+    });
+
+    return {
+      volume,
+      intervals: intervalResults
+    };
+  });
+}

src/image/index.js

@@ -45,6 +45,7 @@ import {MaskSegmentViewHelper} from './maskSegmentViewHelper.js';
 import {SegmentCollection} from './segmentCollection.js';
 import {Diameter, Diameters, Label} from './label.js';
 import {WindowPreset} from './windowPreset.js';
+import {getHistogramIntervalPercentages} from './histogram.js';
 
 export {
   Annotation,
@@ -80,6 +81,7 @@ export {
   equalWl,
   getDefaultDicomSegJson,
   getDefaultDicomRTStructJson,
+  getHistogramIntervalPercentages,
   mergeMaskImages,
   bufferToPolygons,
   simplifyPolygon,

src/image/maskImage.js

@@ -1,15 +1,19 @@
 import {Index} from '../math/index.js';
+import {REAL_WORLD_EPSILON} from '../math/number.js';
 import {Size} from './size.js';
 import {valueRange} from './iterator.js';
 import {LabelingThread} from './labelingThread.js';
 import {SegmentCollection} from './segmentCollection.js';
 import {ImageContour} from './imageContour.js';
 import {Image} from './image.js';
+import {logger} from '../utils/logger.js';
+import {MaskSegmentHelper} from './maskSegmentHelper.js';
 
 /**
  * @import {Geometry} from './geometry.js';
  * @import {RGB} from '../utils/colour.js';
  * @import {Label} from './label.js';
+ * @import {MaskSegment} from '../dicom/dicomSegment.js';
  */
 
 const ML_PER_MM = 0.001; // ml/mm^3
@@ -529,4 +533,149 @@ export class MaskImage extends Image {
     return copy;
   }
 
+  /**
+   * Get per-volume histograms of reference-image intensities for one segment.
+   *
+   * For each volume of `image`, counts how often each rescaled intensity
+   * appears in voxels belonging to the given segment of this mask.
+   * Segmentations are 3D: the same mask is applied to every volume.
+   * Histogram arrays are dense from the min to max intensity found in
+   * that volume (`histogram[value] === 0` means that intensity is absent).
+   *
+   * @param {Image} image The reference image (intensities).
+   * @param {string} segmentId Segment tracking UID, or segment number as
+   *   string.
+   * @returns {{volume: number, histogram: number[]}[]} One entry per volume.
+   */
+  getHistogramBySegment(image, segmentId) {
+    const segment = this.#findSegmentById(segmentId);
+    if (typeof segment === 'undefined') {
+      logger.warn(
+        `getHistogramBySegment: unknown segment id: ${segmentId}`);
+      return [];
+    }
+
+    const imageGeometry = image.getGeometry();
+    const maskGeometry = this.getGeometry();
+    this.#assertCompatibleHistogramGeometry(maskGeometry, imageGeometry);
+    const imageSize = imageGeometry.getSize();
+
+    const segmentValue = typeof segment.displayValue !== 'undefined'
+      ? segment.displayValue
+      : segment.number;
+
+    const numberOfVolumes = imageSize.length() === 4
+      ? imageSize.get(3)
+      : 1;
+    const volumeSize = imageSize.length() === 4
+      ? imageSize.getDimSize(3)
+      : imageSize.getTotalSize();
+
+    const maskBuffer = this.getBuffer();
+    const result = [];
+
+    for (let v = 0; v < numberOfVolumes; ++v) {
+      const start = v * volumeSize;
+      const end = start + volumeSize;
+      /** @type {number[]} */
+      const counts = [];
+      let rmin;
+      let rmax;
+
+      for (let offset = start; offset < end; ++offset) {
+        // mask is always 3D spatial: reuse it for each volume
+        const maskOffset = offset - start;
+        if (maskBuffer[maskOffset] !== segmentValue) {
+          continue;
+        }
+        const value = image.getRescaledValueAtOffset(offset);
+        counts[value] = (counts[value] || 0) + 1;
+        if (typeof rmin === 'undefined' || value < rmin) {
+          rmin = value;
+        }
+        if (typeof rmax === 'undefined' || value > rmax) {
+          rmax = value;
+        }
+      }
+
+      /** @type {number[]} */
+      const histogram = [];
+      if (typeof rmin !== 'undefined') {
+        for (let b = rmin; b <= rmax; ++b) {
+          histogram[b] = counts[b] || 0;
+        }
+      }
+
+      result.push({volume: v, histogram});
+    }
+
+    return result;
+  }
+
+  /**
+   * Find a segment by tracking UID or by number as string.
+   *
+   * @param {string} segmentId The segment id.
+   * @returns {MaskSegment|undefined} The segment, or undefined if not found.
+   */
+  #findSegmentById(segmentId) {
+    const helper = new MaskSegmentHelper(this);
+    const segments = this.getMeta()?.custom?.segments ?? [];
+
+    for (const segment of segments) {
+      if (segment.trackingUid === segmentId) {
+        return segment;
+      }
+    }
+    for (const segment of segments) {
+      if (String(segment.number) === segmentId) {
+        return segment;
+      }
+    }
+    // also allow lookup when segments are only known via the helper
+    // (e.g. freshly added) using the same number-as-string rule
+    const asNumber = Number(segmentId);
+    if (!Number.isNaN(asNumber) && helper.hasSegment(asNumber)) {
+      return helper.getSegment(asNumber);
+    }
+    return undefined;
+  }
+
+  /**
+   * Ensure mask and image share spatial size (dims 0-2), origins and
+   * orientation. Segmentations are 3D; dim 3 of a 4D reference image
+   * is ignored here.
+   *
+   * @param {Geometry} maskGeometry The mask geometry.
+   * @param {Geometry} imageGeometry The reference image geometry.
+   */
+  #assertCompatibleHistogramGeometry(maskGeometry, imageGeometry) {
+    const maskSize = maskGeometry.getSize();
+    const imageSize = imageGeometry.getSize();
+    for (let i = 0; i < 3; ++i) {
+      if (maskSize.get(i) !== imageSize.get(i)) {
+        throw new Error(
+          'getHistogramBySegment: mask and image spatial sizes differ.');
+      }
+    }
+    if (!maskGeometry.getOrientation().isSimilar(
+      imageGeometry.getOrientation(), REAL_WORLD_EPSILON)) {
+      throw new Error(
+        'getHistogramBySegment: mask and image orientations differ.');
+    }
+    const maskOrigins = maskGeometry.getOrigins();
+    const imageOrigins = imageGeometry.getOrigins();
+    if (maskOrigins.length !== imageOrigins.length) {
+      throw new Error(
+        'getHistogramBySegment: mask and image origins differ.');
+    }
+    for (let i = 0; i < maskOrigins.length; ++i) {
+      if (!maskOrigins[i].isSimilar(
+        imageOrigins[i], REAL_WORLD_EPSILON)) {
+        throw new Error(
+          'getHistogramBySegment: mask and image origins differ.');
+      }
+    }
+  }
+
 } // MaskImage class