1078 lines
38 KiB
JavaScript
1078 lines
38 KiB
JavaScript
import { scanIndexToPosition } from "./math.js";
|
||
|
||
export const RESULT_METRIC_KEYS = Object.freeze([
|
||
"peakToSidelobeDb",
|
||
"energyInsideTruthRatio",
|
||
"precision",
|
||
"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;
|
||
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" &&
|
||
descriptor.scale !== "score" &&
|
||
descriptor.scale !== "noise-projection"
|
||
) {
|
||
throw new RangeError("rawMap scale must be power, score, or noise-projection.");
|
||
}
|
||
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 rawMapScale(rawMap) {
|
||
if (rawMap && !Array.isArray(rawMap) && !ArrayBuffer.isView(rawMap)) {
|
||
return rawMap.scale ?? "power";
|
||
}
|
||
return "power";
|
||
}
|
||
|
||
function rawMapObjectiveDirection(rawMap, mapScale) {
|
||
const descriptor =
|
||
rawMap && !Array.isArray(rawMap) && !ArrayBuffer.isView(rawMap)
|
||
? rawMap
|
||
: null;
|
||
const direction =
|
||
descriptor?.objectiveDirection ??
|
||
(mapScale === "noise-projection" ? "minimum" : "maximum");
|
||
if (direction !== "maximum" && direction !== "minimum") {
|
||
throw new RangeError('rawMap objectiveDirection must be "maximum" or "minimum".');
|
||
}
|
||
if (mapScale === "noise-projection" && direction !== "minimum") {
|
||
throw new RangeError("noise-projection maps must use objectiveDirection minimum.");
|
||
}
|
||
if (mapScale !== "noise-projection" && direction !== "maximum") {
|
||
throw new RangeError("power and score maps must use objectiveDirection maximum.");
|
||
}
|
||
return direction;
|
||
}
|
||
|
||
export function normalizedMusicPseudospectrum(values) {
|
||
if (!Array.isArray(values) && !ArrayBuffer.isView(values)) {
|
||
throw new TypeError("noise projection values must be an array or typed array.");
|
||
}
|
||
if (values.length === 0) {
|
||
throw new RangeError("noise projection values must not be empty.");
|
||
}
|
||
let rawMinimum = Number.POSITIVE_INFINITY;
|
||
let rawMaximum = Number.NEGATIVE_INFINITY;
|
||
for (let index = 0; index < values.length; index += 1) {
|
||
const value = requireNonNegative(
|
||
`noiseProjection[${index}]`,
|
||
values[index],
|
||
);
|
||
rawMinimum = Math.min(rawMinimum, value);
|
||
rawMaximum = Math.max(rawMaximum, value);
|
||
}
|
||
const pseudospectrum = new Float64Array(values.length);
|
||
if (rawMinimum === rawMaximum) {
|
||
pseudospectrum.fill(1);
|
||
return {
|
||
values: pseudospectrum,
|
||
transform: "normalized-reciprocal-noise-projection",
|
||
rawMinimum,
|
||
flat: true,
|
||
};
|
||
}
|
||
if (rawMinimum === 0) {
|
||
for (let index = 0; index < values.length; index += 1) {
|
||
pseudospectrum[index] = values[index] === 0 ? 1 : 0;
|
||
}
|
||
} else {
|
||
for (let index = 0; index < values.length; index += 1) {
|
||
pseudospectrum[index] = rawMinimum / values[index];
|
||
}
|
||
}
|
||
return {
|
||
values: pseudospectrum,
|
||
transform: "normalized-reciprocal-noise-projection",
|
||
rawMinimum,
|
||
flat: false,
|
||
};
|
||
}
|
||
|
||
function objectiveMerit(values, direction) {
|
||
if (direction === "maximum") {
|
||
return { values, transform: "identity", rawMinimum: null, flat: false };
|
||
}
|
||
return normalizedMusicPseudospectrum(values);
|
||
}
|
||
|
||
function restoreObjectiveValue(candidate, rawValues) {
|
||
if (!candidate) return null;
|
||
return {
|
||
...candidate,
|
||
value: rawValues[candidate.index],
|
||
meritValue: candidate.value,
|
||
};
|
||
}
|
||
|
||
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;
|
||
}
|
||
|
||
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
|
||
* 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 };
|
||
}
|
||
|
||
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, mapScale = "power") {
|
||
const component = collectHalfPowerIndices(values, scan, peakIndex, halfPowerRatio);
|
||
const peakPosition = scanIndexToPosition(peakIndex, scan);
|
||
const base = {
|
||
calibrated: false,
|
||
method: mapScale === "score"
|
||
? "connected-half-maximum-score-weighted-covariance"
|
||
: "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,
|
||
rawMapScale(rawMap),
|
||
);
|
||
}
|
||
|
||
/**
|
||
* Evaluate one nonnegative spatial map. Noise-projection minima are detected
|
||
* through the standard MUSIC reciprocal pseudospectrum normalized as
|
||
* P/Pmax=qmin/q, while exported peak values remain the unmodified q values.
|
||
* Source z coordinates are projected onto scan.z, matching the image-plane
|
||
* ground-truth convention.
|
||
*/
|
||
export function evaluateLocalizationBand(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 mapScale = options.mapScale ?? rawMapScale(rawMap);
|
||
if (mapScale !== "power" && mapScale !== "score" && mapScale !== "noise-projection") {
|
||
throw new RangeError("mapScale must be power, score, or noise-projection.");
|
||
}
|
||
const objectiveDirection =
|
||
options.objectiveDirection ?? rawMapObjectiveDirection(rawMap, mapScale);
|
||
if (
|
||
(mapScale === "noise-projection" && objectiveDirection !== "minimum") ||
|
||
(mapScale !== "noise-projection" && objectiveDirection !== "maximum")
|
||
) {
|
||
throw new RangeError(
|
||
"noise-projection uses a minimum objective; power and score use a maximum objective.",
|
||
);
|
||
}
|
||
const merit = objectiveMerit(values, objectiveDirection);
|
||
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 detectionMethod = requireDetectionMethod(options.detectionMethod ?? "local-max");
|
||
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 localMeritPeaks = merit.flat
|
||
? []
|
||
: findLocalPeaks(merit.values, scan);
|
||
const globalMeritPeak = localMeritPeaks[0] ?? null;
|
||
const thresholdValue = (globalMeritPeak?.value ?? 0) * detectionThresholdRatio;
|
||
const meritCandidates = detectionCandidates(
|
||
merit.values,
|
||
scan,
|
||
localMeritPeaks,
|
||
thresholdValue,
|
||
detectionMethod,
|
||
);
|
||
const selectedMeritDetections = selectDetections(
|
||
meritCandidates,
|
||
thresholdValue,
|
||
minimumPeakDistanceM,
|
||
maximumDetections,
|
||
);
|
||
const localPeaks = localMeritPeaks.map((peak) =>
|
||
restoreObjectiveValue(peak, values));
|
||
const globalPeak = restoreObjectiveValue(globalMeritPeak, values);
|
||
const detections = selectedMeritDetections.map((detection) => {
|
||
const restored = restoreObjectiveValue(detection, values);
|
||
if (objectiveDirection === "minimum") return restored;
|
||
return {
|
||
...restored,
|
||
beamSpread: measureBeamSpreadFromValues(
|
||
values,
|
||
scan,
|
||
detection.index,
|
||
halfPowerRatio,
|
||
mapScale,
|
||
),
|
||
};
|
||
});
|
||
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 =
|
||
objectiveDirection === "maximum" && globalPeak
|
||
? measureMainLobe(values, scan, globalPeak.index, halfPowerRatio)
|
||
: null;
|
||
const metrics = { precision, recall };
|
||
if (objectiveDirection === "maximum") {
|
||
Object.assign(metrics, {
|
||
peakToSidelobeDb: peakToSidelobeDb(
|
||
values,
|
||
scan,
|
||
truths,
|
||
globalPeak,
|
||
sidelobeExclusionRadiusM,
|
||
dynamicRangeDb,
|
||
),
|
||
energyInsideTruthRatio: energyInsideRatio(
|
||
values,
|
||
scan,
|
||
truths,
|
||
truthRadiusM,
|
||
globalPeak?.value ?? 0,
|
||
),
|
||
mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0,
|
||
});
|
||
}
|
||
|
||
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,
|
||
detectionThresholdDomain:
|
||
objectiveDirection === "minimum"
|
||
? "normalized-music-pseudospectrum"
|
||
: "raw-map-value",
|
||
detectionMethod,
|
||
minimumPeakDistanceM,
|
||
sidelobeExclusionRadiusM,
|
||
halfPowerRatio,
|
||
dynamicRangeDb,
|
||
maximumDetections,
|
||
},
|
||
metrics,
|
||
diagnostics: {
|
||
mapScale,
|
||
objectiveDirection,
|
||
detectionTransform: merit.transform,
|
||
...(merit.rawMinimum === null ? {} : { rawMinimum: merit.rawMinimum }),
|
||
flatObjectiveMap: merit.flat,
|
||
detectionCandidateCount: meritCandidates.length,
|
||
mainLobe,
|
||
mainLobeWidthM: mainLobe?.equivalentDiameterM ?? 0,
|
||
gridResolutionM: resolution,
|
||
},
|
||
};
|
||
}
|
||
|
||
/** Backwards-compatible name for callers that evaluate conventional DAS maps. */
|
||
export function evaluateDasBand(options) {
|
||
return evaluateLocalizationBand(options);
|
||
}
|
||
|
||
/** Select only metric names accepted by result.schema.json. */
|
||
export function selectResultMetrics(evaluationOrMetrics, metricNames) {
|
||
const metrics = evaluationOrMetrics?.metrics ?? evaluationOrMetrics;
|
||
if (!metrics || typeof metrics !== "object") throw new TypeError("evaluation metrics are required.");
|
||
const selectedMetricNames =
|
||
metricNames ?? RESULT_METRIC_KEYS.filter((metricName) => metricName in metrics);
|
||
if (!Array.isArray(selectedMetricNames)) throw new TypeError("metricNames must be an array.");
|
||
const selected = {};
|
||
for (const metricName of selectedMetricNames) {
|
||
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 };
|
||
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;
|
||
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),
|
||
evaluation: {
|
||
version: "1.0.0",
|
||
truthRadiusM: evaluation.parameters.truthRadiusM,
|
||
matchToleranceM: evaluation.parameters.matchToleranceM,
|
||
detectionThresholdRatio: evaluation.parameters.detectionThresholdRatio,
|
||
detectionThresholdValue: evaluation.parameters.detectionThresholdValue,
|
||
objectiveDirection: evaluation.diagnostics?.objectiveDirection ?? "maximum",
|
||
detectionThresholdDomain:
|
||
evaluation.parameters.detectionThresholdDomain ?? "raw-map-value",
|
||
detectionMethod: evaluation.parameters.detectionMethod ?? "local-max",
|
||
minimumPeakDistanceM: evaluation.parameters.minimumPeakDistanceM,
|
||
sidelobeExclusionRadiusM: evaluation.parameters.sidelobeExclusionRadiusM,
|
||
halfPowerRatio: evaluation.parameters.halfPowerRatio,
|
||
dynamicRangeDb: evaluation.parameters.dynamicRangeDb,
|
||
maximumDetections: evaluation.parameters.maximumDetections,
|
||
},
|
||
};
|
||
}
|