micarrsim/packages/acoustics-core/evaluation.js

import { scanIndexToPosition } from "./math.js";

export const RESULT_METRIC_KEYS = Object.freeze([
  "peakToSidelobeDb",
  "energyInsideTruthRatio",
  "precision",
  "recall",
  "mainLobeWidthM",
]);

const RESULT_METRIC_KEY_SET = new Set(RESULT_METRIC_KEYS);
const DEFAULT_MAXIMUM_DETECTIONS = 128;
const DEFAULT_DYNAMIC_RANGE_DB = 120;

function requireFinite(name, value) {
  if (!Number.isFinite(value)) throw new RangeError(`${name} must be finite.`);
  return value;
}

function requireNonNegative(name, value) {
  requireFinite(name, value);
  if (value < 0) throw new RangeError(`${name} must not be negative.`);
  return value;
}

function requireRatio(name, value, { positive = false } = {}) {
  requireFinite(name, value);
  if (value < (positive ? Number.EPSILON : 0) || value > 1) {
    throw new RangeError(`${name} must be ${positive ? "greater than zero and " : ""}at most one.`);
  }
  return value;
}

function validateScan(scan) {
  if (!scan || typeof scan !== "object") throw new TypeError("scan is required.");
  for (const key of ["xMin", "xMax", "yMin", "yMax", "z"]) requireFinite(`scan.${key}`, scan[key]);
  for (const key of ["columns", "rows"]) {
    if (!Number.isInteger(scan[key]) || scan[key] < 1) {
      throw new RangeError(`scan.${key} must be a positive integer.`);
    }
  }
  if (scan.xMin > scan.xMax || scan.yMin > scan.yMax) {
    throw new RangeError("scan minimum bounds must not exceed maximum bounds.");
  }
}

function normalizeRawMap(rawMap, scan) {
  const descriptor = rawMap && !Array.isArray(rawMap) && !ArrayBuffer.isView(rawMap) ? rawMap : null;
  const values = descriptor?.values ?? rawMap;
  if (!Array.isArray(values) && !ArrayBuffer.isView(values)) {
    throw new TypeError("rawMap must be an Array, a typed array, or a map descriptor with values.");
  }
  if (descriptor?.scale !== undefined && descriptor.scale !== "power") {
    throw new RangeError("rawMap scale must be power.");
  }
  if (descriptor?.rows !== undefined && descriptor.rows !== scan.rows) {
    throw new RangeError("rawMap rows do not match scan.rows.");
  }
  if (descriptor?.columns !== undefined && descriptor.columns !== scan.columns) {
    throw new RangeError("rawMap columns do not match scan.columns.");
  }
  const expectedLength = scan.rows * scan.columns;
  if (values.length !== expectedLength) {
    throw new RangeError(`rawMap length must be scan.rows × scan.columns (${expectedLength}).`);
  }
  for (let index = 0; index < values.length; index += 1) {
    requireNonNegative(`rawMap[${index}]`, values[index]);
  }
  return values;
}

function normalizeSources(sources) {
  if (!Array.isArray(sources)) throw new TypeError("sources must be an array.");
  const ids = new Set();
  return sources.map((source, index) => {
    if (!source || typeof source !== "object") throw new TypeError(`sources[${index}] must be an object.`);
    if (typeof source.id !== "string" || source.id.length === 0) {
      throw new TypeError(`sources[${index}].id must be a non-empty string.`);
    }
    if (ids.has(source.id)) throw new RangeError(`duplicate source id: ${source.id}`);
    ids.add(source.id);
    if (!Array.isArray(source.position) || source.position.length !== 3) {
      throw new TypeError(`sources[${index}].position must contain three coordinates.`);
    }
    source.position.forEach((coordinate, coordinateIndex) => {
      requireFinite(`sources[${index}].position[${coordinateIndex}]`, coordinate);
    });
    return source;
  });
}

function gridResolution(scan) {
  const x = scan.columns > 1 ? Math.abs(scan.xMax - scan.xMin) / (scan.columns - 1) : 0;
  const y = scan.rows > 1 ? Math.abs(scan.yMax - scan.yMin) / (scan.rows - 1) : 0;
  return { x, y, diagonal: Math.hypot(x, y) };
}

function planarDistance(a, b) {
  return Math.hypot(a[0] - b[0], a[1] - b[1]);
}

function comparePeaks(a, b) {
  return b.value - a.value || a.index - b.index;
}

/**
 * Return positive local maxima with each equal-valued 8-connected plateau
 * represented by its lowest row-major index. A positive flat map falls back
 * to the deterministic global peak at index zero.
 */
export function findLocalPeaks(rawMap, scan) {
  validateScan(scan);
  const values = normalizeRawMap(rawMap, scan);
  let globalIndex = -1;
  let globalValue = 0;
  for (let index = 0; index < values.length; index += 1) {
    if (values[index] > globalValue) {
      globalValue = values[index];
      globalIndex = index;
    }
  }
  if (!(globalValue > 0)) return [];

  const visited = new Uint8Array(values.length);
  const queue = new Int32Array(values.length);
  const peaks = [];
  for (let startIndex = 0; startIndex < values.length; startIndex += 1) {
    if (visited[startIndex] || !(values[startIndex] > 0)) continue;
    const plateauValue = values[startIndex];
    let queueStart = 0;
    let queueEnd = 1;
    let representativeIndex = startIndex;
    let hasHigherNeighbor = false;
    let hasLowerNeighbor = false;
    queue[0] = startIndex;
    visited[startIndex] = 1;

    while (queueStart < queueEnd) {
      const index = queue[queueStart];
      queueStart += 1;
      representativeIndex = Math.min(representativeIndex, index);
      const row = Math.floor(index / scan.columns);
      const column = index % scan.columns;
      const rowStart = Math.max(0, row - 1);
      const rowEnd = Math.min(scan.rows - 1, row + 1);
      const columnStart = Math.max(0, column - 1);
      const columnEnd = Math.min(scan.columns - 1, column + 1);
      for (let neighborRow = rowStart; neighborRow <= rowEnd; neighborRow += 1) {
        for (let neighborColumn = columnStart; neighborColumn <= columnEnd; neighborColumn += 1) {
          if (neighborRow === row && neighborColumn === column) continue;
          const neighborIndex = neighborRow * scan.columns + neighborColumn;
          const neighborValue = values[neighborIndex];
          if (neighborValue === plateauValue) {
            if (!visited[neighborIndex]) {
              visited[neighborIndex] = 1;
              queue[queueEnd] = neighborIndex;
              queueEnd += 1;
            }
          } else if (neighborValue > plateauValue) {
            hasHigherNeighbor = true;
          } else {
            hasLowerNeighbor = true;
          }
        }
      }
    }

    if (!hasHigherNeighbor && hasLowerNeighbor) {
      peaks.push({
        index: representativeIndex,
        position: scanIndexToPosition(representativeIndex, scan),
        value: plateauValue,
      });
    }
  }

  if (!peaks.some(({ index }) => index === globalIndex)) {
    peaks.push({ index: globalIndex, position: scanIndexToPosition(globalIndex, scan), value: globalValue });
  }
  peaks.sort(comparePeaks);
  return peaks;
}

function selectDetections(localPeaks, thresholdValue, minimumPeakDistanceM, maximumDetections) {
  const detections = [];
  for (const peak of localPeaks) {
    if (peak.value < thresholdValue) continue;
    if (
      minimumPeakDistanceM > 0 &&
      detections.some((selected) => planarDistance(selected.position, peak.position) <= minimumPeakDistanceM)
    ) {
      continue;
    }
    detections.push({ ...peak });
    if (detections.length === maximumDetections) break;
  }
  return detections;
}

function matchTruths(truths, detections, matchToleranceM) {
  const truthToDetection = new Int32Array(truths.length);
  truthToDetection.fill(-1);

  // First maximize valid tolerance-based matches. Stronger detections are
  // considered first; each one's candidate truths are ordered by distance.
  const adjacency = detections.map((detection) =>
    truths
      .map((truth, truthIndex) => ({
        truthIndex,
        distance: planarDistance(detection.position, truth.position),
      }))
      .filter(({ distance }) => distance <= matchToleranceM)
      .sort((a, b) => a.distance - b.distance || a.truthIndex - b.truthIndex),
  );

  function assign(detectionIndex, visitedTruths) {
    for (const { truthIndex } of adjacency[detectionIndex]) {
      if (visitedTruths[truthIndex]) continue;
      visitedTruths[truthIndex] = 1;
      const previousDetection = truthToDetection[truthIndex];
      if (previousDetection === -1 || assign(previousDetection, visitedTruths)) {
        truthToDetection[truthIndex] = detectionIndex;
        return true;
      }
    }
    return false;
  }

  detections.forEach((_, detectionIndex) => assign(detectionIndex, new Uint8Array(truths.length)));
  const usedDetections = new Uint8Array(detections.length);
  truthToDetection.forEach((detectionIndex) => {
    if (detectionIndex >= 0) usedDetections[detectionIndex] = 1;
  });

  // Preserve a nearest-peak localization measurement for still-unmatched
  // truths without counting a distant assignment as a true positive.
  const remainingPairs = [];
  truths.forEach((truth, truthIndex) => {
    if (truthToDetection[truthIndex] >= 0) return;
    detections.forEach((detection, detectionIndex) => {
      if (!usedDetections[detectionIndex]) {
        remainingPairs.push({
          truthIndex,
          detectionIndex,
          distance: planarDistance(detection.position, truth.position),
        });
      }
    });
  });
  remainingPairs.sort(
    (a, b) => a.distance - b.distance || a.truthIndex - b.truthIndex || a.detectionIndex - b.detectionIndex,
  );
  for (const pair of remainingPairs) {
    if (truthToDetection[pair.truthIndex] >= 0 || usedDetections[pair.detectionIndex]) continue;
    truthToDetection[pair.truthIndex] = pair.detectionIndex;
    usedDetections[pair.detectionIndex] = 1;
  }

  const sourceMatches = truths.map((truth, truthIndex) => {
    const detectionIndex = truthToDetection[truthIndex];
    if (detectionIndex < 0) {
      return {
        sourceId: truth.sourceId,
        truthPosition: [...truth.position],
        peak: null,
        localizationErrorM: null,
        withinTolerance: false,
      };
    }
    const detection = detections[detectionIndex];
    const localizationErrorM = planarDistance(detection.position, truth.position);
    return {
      sourceId: truth.sourceId,
      truthPosition: [...truth.position],
      peak: { ...detection, position: [...detection.position] },
      localizationErrorM,
      withinTolerance: localizationErrorM <= matchToleranceM,
    };
  });

  return {
    sourceMatches,
    truePositiveCount: sourceMatches.filter(({ withinTolerance }) => withinTolerance).length,
  };
}

function energyInsideRatio(values, scan, truths, radiusM, globalValue) {
  if (!(globalValue > 0) || truths.length === 0) return 0;
  let total = 0;
  let inside = 0;
  for (let index = 0; index < values.length; index += 1) {
    const scaledValue = values[index] / globalValue;
    total += scaledValue;
    if (truths.some((truth) => planarDistance(scanIndexToPosition(index, scan), truth.position) <= radiusM)) {
      inside += scaledValue;
    }
  }
  return total > 0 ? Math.min(1, Math.max(0, inside / total)) : 0;
}

function peakToSidelobeDb(values, scan, truths, globalPeak, exclusionRadiusM, dynamicRangeDb) {
  if (!globalPeak) return 0;
  const exclusionCenters = truths.length > 0 ? truths.map(({ position }) => position) : [globalPeak.position];
  let mainPeak = truths.length > 0 ? 0 : globalPeak.value;
  let sidelobePeak = 0;
  for (let index = 0; index < values.length; index += 1) {
    const position = scanIndexToPosition(index, scan);
    const isMainRegion = exclusionCenters.some((center) => planarDistance(position, center) <= exclusionRadiusM);
    if (isMainRegion) mainPeak = Math.max(mainPeak, values[index]);
    else sidelobePeak = Math.max(sidelobePeak, values[index]);
  }
  if (!(mainPeak > 0)) return sidelobePeak > 0 ? -dynamicRangeDb : 0;
  if (!(sidelobePeak > 0)) return dynamicRangeDb;
  return Math.max(
    -dynamicRangeDb,
    Math.min(dynamicRangeDb, 10 * (Math.log10(mainPeak) - Math.log10(sidelobePeak))),
  );
}

/** Measure the 8-connected half-power component containing a selected peak. */
export function measureMainLobe(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 });
  const peakValue = values[peakIndex];
  if (!(peakValue > 0)) {
    return {
      peakIndex,
      thresholdValue: 0,
      pixelCount: 0,
      widthXM: 0,
      widthYM: 0,
      equivalentDiameterM: 0,
    };
  }

  const thresholdValue = peakValue * halfPowerRatio;
  const visited = new Uint8Array(values.length);
  const queue = new Int32Array(values.length);
  let queueStart = 0;
  let queueEnd = 1;
  let minimumRow = Math.floor(peakIndex / scan.columns);
  let maximumRow = minimumRow;
  let minimumColumn = peakIndex % scan.columns;
  let maximumColumn = minimumColumn;
  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;
    minimumRow = Math.min(minimumRow, row);
    maximumRow = Math.max(maximumRow, row);
    minimumColumn = Math.min(minimumColumn, column);
    maximumColumn = Math.max(maximumColumn, column);
    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;
        }
      }
    }
  }

  const resolution = gridResolution(scan);
  const widthXM = scan.columns > 1 ? (maximumColumn - minimumColumn + 1) * resolution.x : 0;
  const widthYM = scan.rows > 1 ? (maximumRow - minimumRow + 1) * resolution.y : 0;
  let equivalentDiameterM = 0;
  if (resolution.x > 0 && resolution.y > 0) {
    equivalentDiameterM = 2 * Math.sqrt((queueEnd * resolution.x * resolution.y) / Math.PI);
  } else {
    equivalentDiameterM = Math.max(widthXM, widthYM);
  }
  return { peakIndex, thresholdValue, pixelCount: queueEnd, widthXM, widthYM, equivalentDiameterM };
}

/**
 * Evaluate one linear-power DAS band. Source z coordinates are projected onto
 * scan.z, matching the image-plane ground-truth convention.
 */
export function evaluateDasBand(options) {
  if (!options || typeof options !== "object") throw new TypeError("evaluation options are required.");
  const { rawMap, scan } = options;
  validateScan(scan);
  const values = normalizeRawMap(rawMap, scan);
  const sources = normalizeSources(options.sources ?? options.scene?.sources ?? []);
  const truths = sources.map((source) => ({
    sourceId: source.id,
    position: [source.position[0], source.position[1], scan.z],
  }));
  const resolution = gridResolution(scan);
  const defaultRadius = resolution.diagonal > 0 ? resolution.diagonal * 1.5 : 0;
  const truthRadiusM = requireNonNegative("truthRadiusM", options.truthRadiusM ?? defaultRadius);
  const matchToleranceM = requireNonNegative(
    "matchToleranceM",
    options.matchToleranceM ?? Math.max(defaultRadius, truthRadiusM),
  );
  const detectionThresholdRatio = requireRatio(
    "detectionThresholdRatio",
    options.detectionThresholdRatio ?? 0.5,
  );
  const minimumPeakDistanceM = requireNonNegative(
    "minimumPeakDistanceM",
    options.minimumPeakDistanceM ?? resolution.diagonal,
  );
  const sidelobeExclusionRadiusM = requireNonNegative(
    "sidelobeExclusionRadiusM",
    options.sidelobeExclusionRadiusM ?? truthRadiusM,
  );
  const halfPowerRatio = requireRatio("halfPowerRatio", options.halfPowerRatio ?? 0.5, { positive: true });
  const dynamicRangeDb = requireNonNegative(
    "dynamicRangeDb",
    options.dynamicRangeDb ?? DEFAULT_DYNAMIC_RANGE_DB,
  );
  const maximumDetections = options.maximumDetections ?? DEFAULT_MAXIMUM_DETECTIONS;
  if (!Number.isInteger(maximumDetections) || maximumDetections < 1 || maximumDetections > 128) {
    throw new RangeError("maximumDetections must be an integer from 1 through 128.");
  }

  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 { 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;
  const mainLobe = globalPeak ? measureMainLobe(values, scan, globalPeak.index, halfPowerRatio) : null;

  return {
    globalPeak,
    localPeaks,
    detections,
    sourceMatches,
    counts: {
      truth: truths.length,
      detected: detections.length,
      truePositive: truePositiveCount,
      falsePositive: detections.length - truePositiveCount,
      falseNegative: truths.length - truePositiveCount,
    },
    parameters: {
      truthRadiusM,
      matchToleranceM,
      detectionThresholdRatio,
      detectionThresholdValue: thresholdValue,
      minimumPeakDistanceM,
      sidelobeExclusionRadiusM,
      halfPowerRatio,
      dynamicRangeDb,
      maximumDetections,
    },
    metrics: {
      peakToSidelobeDb: peakToSidelobeDb(
        values,
        scan,
        truths,
        globalPeak,
        sidelobeExclusionRadiusM,
        dynamicRangeDb,
      ),
      energyInsideTruthRatio: energyInsideRatio(values, scan, truths, truthRadiusM, globalPeak?.value ?? 0),
      precision,
      recall,
      mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0,
    },
    diagnostics: {
      mainLobe,
      mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0,
      gridResolutionM: resolution,
    },
  };
}

/** Select only metric names accepted by result.schema.json. */
export function selectResultMetrics(evaluationOrMetrics, metricNames = RESULT_METRIC_KEYS) {
  const metrics = evaluationOrMetrics?.metrics ?? evaluationOrMetrics;
  if (!metrics || typeof metrics !== "object") throw new TypeError("evaluation metrics are required.");
  if (!Array.isArray(metricNames)) throw new TypeError("metricNames must be an array.");
  const selected = {};
  for (const metricName of metricNames) {
    if (!RESULT_METRIC_KEY_SET.has(metricName)) {
      throw new RangeError(`result.schema.json does not allow metric: ${metricName}`);
    }
    const value = metrics[metricName];
    requireFinite(`metrics.${metricName}`, value);
    selected[metricName] = value;
  }
  return selected;
}

/** Convert an evaluation to the schema-safe `peaks` and `metrics` band fields. */
export function exportBandEvaluation(evaluation, options = {}) {
  if (!evaluation || typeof evaluation !== "object") throw new TypeError("evaluation is required.");
  const maximumPeaks = options.maximumPeaks ?? 128;
  if (!Number.isInteger(maximumPeaks) || maximumPeaks < 0 || maximumPeaks > 128) {
    throw new RangeError("maximumPeaks must be an integer from zero through 128.");
  }
  const matchByPeakIndex = new Map();
  for (const match of evaluation.sourceMatches ?? []) {
    if (match.withinTolerance && match.peak) matchByPeakIndex.set(match.peak.index, match);
  }
  const peaks = (evaluation.detections ?? []).slice(0, maximumPeaks).map((peak) => {
    const resultPeak = { position: [...peak.position], value: peak.value };
    const match = matchByPeakIndex.get(peak.index);
    if (match) {
      resultPeak.sourceId = match.sourceId;
      resultPeak.localizationErrorM = match.localizationErrorM;
      const resolution = evaluation.diagnostics?.gridResolutionM;
      if (resolution?.x > 0 && resolution?.y > 0) {
        resultPeak.localizationErrorPixels = Math.hypot(
          (match.peak.position[0] - match.truthPosition[0]) / resolution.x,
          (match.peak.position[1] - match.truthPosition[1]) / resolution.y,
        );
      }
    }
    return resultPeak;
  });
  return {
    peaks,
    metrics: selectResultMetrics(evaluation, options.metricNames ?? RESULT_METRIC_KEYS),
    evaluation: {
      version: "1.0.0",
      truthRadiusM: evaluation.parameters.truthRadiusM,
      matchToleranceM: evaluation.parameters.matchToleranceM,
      detectionThresholdRatio: evaluation.parameters.detectionThresholdRatio,
      detectionThresholdValue: evaluation.parameters.detectionThresholdValue,
      minimumPeakDistanceM: evaluation.parameters.minimumPeakDistanceM,
      sidelobeExclusionRadiusM: evaluation.parameters.sidelobeExclusionRadiusM,
      halfPowerRatio: evaluation.parameters.halfPowerRatio,
      dynamicRangeDb: evaluation.parameters.dynamicRangeDb,
      maximumDetections: evaluation.parameters.maximumDetections,
    },
  };
}