feat: visualize beam spread around peaks

packages/acoustics-core/evaluation.js

@@ -386,6 +386,129 @@ export function measureMainLobe(rawMap, scan, peakIndex, halfPowerRatio = 0.5) {
   return { peakIndex, thresholdValue, pixelCount: queueEnd, widthXM, widthYM, equivalentDiameterM };
 }
 
+function collectHalfPowerIndices(values, scan, peakIndex, halfPowerRatio) {
+  const peakValue = values[peakIndex];
+  if (!(peakValue > 0)) return { indices: new Int32Array(0), thresholdValue: 0 };
+
+  const thresholdValue = peakValue * halfPowerRatio;
+  const visited = new Uint8Array(values.length);
+  const queue = new Int32Array(values.length);
+  let queueStart = 0;
+  let queueEnd = 1;
+  queue[0] = peakIndex;
+  visited[peakIndex] = 1;
+  while (queueStart < queueEnd) {
+    const index = queue[queueStart];
+    queueStart += 1;
+    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
+      ) {
+        const neighborIndex = neighborRow * scan.columns + neighborColumn;
+        if (!visited[neighborIndex] && values[neighborIndex] >= thresholdValue) {
+          visited[neighborIndex] = 1;
+          queue[queueEnd] = neighborIndex;
+          queueEnd += 1;
+        }
+      }
+    }
+  }
+  return { indices: queue.subarray(0, queueEnd), thresholdValue };
+}
+
+function measureBeamSpreadFromValues(values, scan, peakIndex, halfPowerRatio) {
+  const component = collectHalfPowerIndices(values, scan, peakIndex, halfPowerRatio);
+  const peakPosition = scanIndexToPosition(peakIndex, scan);
+  const base = {
+    calibrated: false,
+    method: "connected-half-power-power-weighted-covariance",
+    halfPowerRatio,
+    thresholdValue: component.thresholdValue,
+    pixelCount: component.indices.length,
+    centroid: [...peakPosition],
+    covarianceM2: { xx: 0, xy: 0, yy: 0 },
+    ellipse: {
+      axisScale: "rms",
+      angleReference: "counterclockwise-from-positive-x",
+      majorRmsRadiusM: 0,
+      minorRmsRadiusM: 0,
+      angleDeg: 0,
+    },
+  };
+  if (component.indices.length === 0) return base;
+
+  // Normalizing by the component maximum changes no weighted moments and
+  // prevents a sum of otherwise finite power values from overflowing.
+  let maximum = 0;
+  for (const index of component.indices) maximum = Math.max(maximum, values[index]);
+  let totalWeight = 0;
+  let weightedX = 0;
+  let weightedY = 0;
+  for (const index of component.indices) {
+    const weight = values[index] / maximum;
+    const position = scanIndexToPosition(index, scan);
+    totalWeight += weight;
+    weightedX += weight * position[0];
+    weightedY += weight * position[1];
+  }
+  const centroidX = scan.columns === 1 ? scan.xMin : weightedX / totalWeight;
+  const centroidY = scan.rows === 1 ? scan.yMax : weightedY / totalWeight;
+  let xx = 0;
+  let xy = 0;
+  let yy = 0;
+  for (const index of component.indices) {
+    const weight = values[index] / maximum;
+    const position = scanIndexToPosition(index, scan);
+    const x = position[0] - centroidX;
+    const y = position[1] - centroidY;
+    xx += weight * x * x;
+    xy += weight * x * y;
+    yy += weight * y * y;
+  }
+  xx = scan.columns === 1 ? 0 : Math.max(0, xx / totalWeight);
+  xy = scan.columns === 1 || scan.rows === 1 ? 0 : xy / totalWeight;
+  yy = scan.rows === 1 ? 0 : Math.max(0, yy / totalWeight);
+
+  const trace = xx + yy;
+  const discriminant = Math.hypot(xx - yy, 2 * xy);
+  const majorVariance = Math.max(0, (trace + discriminant) / 2);
+  const minorVariance = Math.max(0, (trace - discriminant) / 2);
+  const covarianceScale = Math.max(Math.abs(xx), Math.abs(xy), Math.abs(yy));
+  const directionDefined = covarianceScale > 0 && discriminant > Number.EPSILON * covarianceScale * 8;
+  let angleDeg = directionDefined ? (Math.atan2(2 * xy, xx - yy) * 90) / Math.PI : 0;
+  if (Object.is(angleDeg, -0)) angleDeg = 0;
+
+  return {
+    ...base,
+    centroid: [centroidX, centroidY, scan.z],
+    covarianceM2: { xx, xy, yy },
+    ellipse: {
+      ...base.ellipse,
+      majorRmsRadiusM: Math.sqrt(majorVariance),
+      minorRmsRadiusM: Math.sqrt(minorVariance),
+      angleDeg,
+    },
+  };
+}
+
+/**
+ * Fit a deterministic, power-weighted covariance ellipse to the connected
+ * half-power lobe. This describes beam spread, not probability or confidence.
+ */
+export function measureBeamSpread(rawMap, scan, peakIndex, halfPowerRatio = 0.5) {
+  validateScan(scan);
+  const values = normalizeRawMap(rawMap, scan);
+  if (!Number.isInteger(peakIndex) || peakIndex < 0 || peakIndex >= values.length) {
+    throw new RangeError("peakIndex is outside rawMap.");
+  }
+  requireRatio("halfPowerRatio", halfPowerRatio, { positive: true });
+  return measureBeamSpreadFromValues(values, scan, peakIndex, halfPowerRatio);
+}
+
 /**
  * Evaluate one linear-power DAS band. Source z coordinates are projected onto
  * scan.z, matching the image-plane ground-truth convention.
@@ -432,7 +555,12 @@ export function evaluateDasBand(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);
+  const detections = selectDetections(localPeaks, thresholdValue, minimumPeakDistanceM, maximumDetections).map(
+    (detection) => ({
+      ...detection,
+      beamSpread: measureBeamSpreadFromValues(values, scan, detection.index, halfPowerRatio),
+    }),
+  );
   const { sourceMatches, truePositiveCount } = matchTruths(truths, detections, matchToleranceM);
   const precision = detections.length > 0 ? truePositiveCount / detections.length : truths.length === 0 ? 1 : 0;
   const recall = truths.length > 0 ? truePositiveCount / truths.length : 1;
@@ -513,6 +641,18 @@ export function exportBandEvaluation(evaluation, options = {}) {
   }
   const peaks = (evaluation.detections ?? []).slice(0, maximumPeaks).map((peak) => {
     const resultPeak = { position: [...peak.position], value: peak.value };
+    if (peak.beamSpread) {
+      resultPeak.beamSpread = {
+        calibrated: peak.beamSpread.calibrated,
+        method: peak.beamSpread.method,
+        halfPowerRatio: peak.beamSpread.halfPowerRatio,
+        thresholdValue: peak.beamSpread.thresholdValue,
+        pixelCount: peak.beamSpread.pixelCount,
+        centroid: [...peak.beamSpread.centroid],
+        covarianceM2: { ...peak.beamSpread.covarianceM2 },
+        ellipse: { ...peak.beamSpread.ellipse },
+      };
+    }
     const match = matchByPeakIndex.get(peak.index);
     if (match) {
       resultPeak.sourceId = match.sourceId;