feat: add standard peak grouping methods

2026-08-29 09:00:40 +09:00
parent 447f77dfa7
commit e048ccbadc
7 changed files with 299 additions and 4 deletions
@@ -7,6 +7,11 @@ export const RESULT_METRIC_KEYS = Object.freeze([
  "recall",
  "mainLobeWidthM",
 ]);
+export const PEAK_DETECTION_METHODS = Object.freeze([
+  "local-max",
+  "connected-components",
+  "watershed",
+]);

 const RESULT_METRIC_KEY_SET = new Set(RESULT_METRIC_KEYS);
 const DEFAULT_MAXIMUM_DETECTIONS = 128;
@@ -114,6 +119,214 @@ function comparePeaks(a, b) {
  return b.value - a.value || a.index - b.index;
 }

+function requireDetectionMethod(value = "local-max") {
+  if (!PEAK_DETECTION_METHODS.includes(value)) {
+    throw new RangeError(
+      "detectionMethod must be local-max, connected-components, or watershed.",
+    );
+  }
+  return value;
+}
+
+function forEachNeighbor(index, scan, callback) {
+  const row = Math.floor(index / scan.columns);
+  const column = index % scan.columns;
+  for (
+    let neighborRow = Math.max(0, row - 1);
+    neighborRow <= Math.min(scan.rows - 1, row + 1);
+    neighborRow += 1
+  ) {
+    for (
+      let neighborColumn = Math.max(0, column - 1);
+      neighborColumn <= Math.min(scan.columns - 1, column + 1);
+      neighborColumn += 1
+    ) {
+      if (neighborRow === row && neighborColumn === column) continue;
+      callback(neighborRow * scan.columns + neighborColumn);
+    }
+  }
+}
+
+function centroidCandidate(values, scan, indices, peakIndex, extra = {}) {
+  let totalWeight = 0;
+  let weightedX = 0;
+  let weightedY = 0;
+  for (const index of indices) {
+    const weight = values[index];
+    const position = scanIndexToPosition(index, scan);
+    totalWeight += weight;
+    weightedX += position[0] * weight;
+    weightedY += position[1] * weight;
+  }
+  const fallback = scanIndexToPosition(peakIndex, scan);
+  return {
+    index: peakIndex,
+    position: totalWeight > 0
+      ? [weightedX / totalWeight, weightedY / totalWeight, scan.z]
+      : fallback,
+    value: values[peakIndex],
+    clusterPixelCount: indices.length,
+    ...extra,
+  };
+}
+
+function findConnectedComponentCandidates(values, scan, thresholdValue) {
+  const visited = new Uint8Array(values.length);
+  const queue = new Int32Array(values.length);
+  const candidates = [];
+  const active = (index) => values[index] > 0 && values[index] >= thresholdValue;
+  for (let startIndex = 0; startIndex < values.length; startIndex += 1) {
+    if (visited[startIndex] || !active(startIndex)) continue;
+    let queueStart = 0;
+    let queueEnd = 1;
+    let peakIndex = startIndex;
+    queue[0] = startIndex;
+    visited[startIndex] = 1;
+    while (queueStart < queueEnd) {
+      const index = queue[queueStart];
+      queueStart += 1;
+      if (
+        values[index] > values[peakIndex] ||
+        (values[index] === values[peakIndex] && index < peakIndex)
+      ) {
+        peakIndex = index;
+      }
+      forEachNeighbor(index, scan, (neighborIndex) => {
+        if (visited[neighborIndex] || !active(neighborIndex)) return;
+        visited[neighborIndex] = 1;
+        queue[queueEnd] = neighborIndex;
+        queueEnd += 1;
+      });
+    }
+    candidates.push(
+      centroidCandidate(
+        values,
+        scan,
+        queue.subarray(0, queueEnd),
+        peakIndex,
+        { clusterMethod: "connected-components" },
+      ),
+    );
+  }
+  candidates.sort(comparePeaks);
+  return candidates;
+}
+
+function watershedEntryHigher(left, right) {
+  return (
+    left.priority > right.priority ||
+    (
+      left.priority === right.priority &&
+      (
+        left.seedOrder < right.seedOrder ||
+        (left.seedOrder === right.seedOrder && left.index < right.index)
+      )
+    )
+  );
+}
+
+function pushWatershedEntry(heap, entry) {
+  heap.push(entry);
+  let index = heap.length - 1;
+  while (index > 0) {
+    const parent = Math.floor((index - 1) / 2);
+    if (!watershedEntryHigher(heap[index], heap[parent])) break;
+    [heap[index], heap[parent]] = [heap[parent], heap[index]];
+    index = parent;
+  }
+}
+
+function popWatershedEntry(heap) {
+  const first = heap[0];
+  const last = heap.pop();
+  if (heap.length === 0) return first;
+  heap[0] = last;
+  let index = 0;
+  while (true) {
+    const left = index * 2 + 1;
+    const right = left + 1;
+    let best = index;
+    if (left < heap.length && watershedEntryHigher(heap[left], heap[best])) best = left;
+    if (right < heap.length && watershedEntryHigher(heap[right], heap[best])) best = right;
+    if (best === index) break;
+    [heap[index], heap[best]] = [heap[best], heap[index]];
+    index = best;
+  }
+  return first;
+}
+
+function findWatershedCandidates(values, scan, thresholdValue, localPeaks) {
+  const seeds = localPeaks.filter((peak) => peak.value >= thresholdValue);
+  if (seeds.length === 0) return [];
+  const active = (index) => values[index] > 0 && values[index] >= thresholdValue;
+  const assignments = new Int32Array(values.length);
+  assignments.fill(-1);
+  const priorities = new Float64Array(values.length);
+  priorities.fill(Number.NEGATIVE_INFINITY);
+  const heap = [];
+  seeds.forEach((seed, seedOrder) => {
+    assignments[seed.index] = seedOrder;
+    priorities[seed.index] = seed.value;
+    pushWatershedEntry(heap, {
+      index: seed.index,
+      priority: seed.value,
+      seedOrder,
+    });
+  });
+
+  while (heap.length > 0) {
+    const current = popWatershedEntry(heap);
+    if (
+      priorities[current.index] !== current.priority ||
+      assignments[current.index] !== current.seedOrder
+    ) {
+      continue;
+    }
+    forEachNeighbor(current.index, scan, (neighborIndex) => {
+      if (!active(neighborIndex)) return;
+      const nextPriority = Math.min(current.priority, values[neighborIndex]);
+      const previousPriority = priorities[neighborIndex];
+      const previousSeed = assignments[neighborIndex];
+      if (
+        nextPriority < previousPriority ||
+        (nextPriority === previousPriority && previousSeed >= 0 && previousSeed <= current.seedOrder)
+      ) {
+        return;
+      }
+      priorities[neighborIndex] = nextPriority;
+      assignments[neighborIndex] = current.seedOrder;
+      pushWatershedEntry(heap, {
+        index: neighborIndex,
+        priority: nextPriority,
+        seedOrder: current.seedOrder,
+      });
+    });
+  }
+
+  const basinIndices = seeds.map(() => []);
+  assignments.forEach((seedOrder, index) => {
+    if (seedOrder >= 0) basinIndices[seedOrder].push(index);
+  });
+  const candidates = seeds.map((seed, seedOrder) =>
+    centroidCandidate(
+      values,
+      scan,
+      basinIndices[seedOrder],
+      seed.index,
+      { clusterMethod: "watershed", basin: seedOrder + 1 },
+    ));
+  candidates.sort(comparePeaks);
+  return candidates;
+}
+
+function detectionCandidates(values, scan, localPeaks, thresholdValue, detectionMethod) {
+  if (detectionMethod === "local-max") return localPeaks;
+  if (detectionMethod === "connected-components") {
+    return findConnectedComponentCandidates(values, scan, thresholdValue);
+  }
+  return findWatershedCandidates(values, scan, thresholdValue, localPeaks);
+}
+
 /**
 * Return positive local maxima with each equal-valued 8-connected plateau
 * represented by its lowest row-major index. A positive flat map falls back
@@ -557,6 +770,7 @@ export function evaluateLocalizationBand(options) {
    "detectionThresholdRatio",
    options.detectionThresholdRatio ?? 0.5,
  );
+  const detectionMethod = requireDetectionMethod(options.detectionMethod ?? "local-max");
  const minimumPeakDistanceM = requireNonNegative(
    "minimumPeakDistanceM",
    options.minimumPeakDistanceM ?? resolution.diagonal,
@@ -578,7 +792,14 @@ export function evaluateLocalizationBand(options) {
  const localPeaks = findLocalPeaks(values, scan);
  const globalPeak = localPeaks.length > 0 ? { ...localPeaks[0], position: [...localPeaks[0].position] } : null;
  const thresholdValue = (globalPeak?.value ?? 0) * detectionThresholdRatio;
-  const detections = selectDetections(localPeaks, thresholdValue, minimumPeakDistanceM, maximumDetections).map(
+  const candidates = detectionCandidates(
+    values,
+    scan,
+    localPeaks,
+    thresholdValue,
+    detectionMethod,
+  );
+  const detections = selectDetections(candidates, thresholdValue, minimumPeakDistanceM, maximumDetections).map(
    (detection) => ({
      ...detection,
      beamSpread: measureBeamSpreadFromValues(
@@ -612,6 +833,7 @@ export function evaluateLocalizationBand(options) {
      matchToleranceM,
      detectionThresholdRatio,
      detectionThresholdValue: thresholdValue,
+      detectionMethod,
      minimumPeakDistanceM,
      sidelobeExclusionRadiusM,
      halfPowerRatio,
@@ -634,6 +856,7 @@ export function evaluateLocalizationBand(options) {
    },
    diagnostics: {
      mapScale,
+      detectionCandidateCount: candidates.length,
      mainLobe,
      mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0,
      gridResolutionM: resolution,
@@ -711,6 +934,7 @@ export function exportBandEvaluation(evaluation, options = {}) {
      matchToleranceM: evaluation.parameters.matchToleranceM,
      detectionThresholdRatio: evaluation.parameters.detectionThresholdRatio,
      detectionThresholdValue: evaluation.parameters.detectionThresholdValue,
+      detectionMethod: evaluation.parameters.detectionMethod ?? "local-max",
      minimumPeakDistanceM: evaluation.parameters.minimumPeakDistanceM,
      sidelobeExclusionRadiusM: evaluation.parameters.sidelobeExclusionRadiusM,
      halfPowerRatio: evaluation.parameters.halfPowerRatio,