feat: improve min level determination contours

src/data/data2d/Spectrum2D/contours.ts

@@ -2,10 +2,12 @@ import type { Spectrum2D, Spectrum } from '@zakodium/nmrium-core';
 import { isSpectrum2DFt } from '@zakodium/nmrium-core';
 import type { NmrData2DFt } from 'cheminfo-types';
 import { Conrec } from 'ml-conrec';
-import { xMaxAbsoluteValue } from 'ml-spectra-processing';
+import { matrixToArray } from 'ml-spectra-processing';
 
 import type { SpectrumFTData } from '../../../component/hooks/use2DReducer.tsx';
-import { calculateSanPlot } from '../../utilities/calculateSanPlot.js';
+
+import { estimateNoiseLevel } from './findBestMinContour/estimateNoiseLevel.js';
+import { findAutomaticContourLevels } from './findBestMinContour/findAutomaticContourLevels.ts';
 
 interface Level {
   positive: ContourItem;
@@ -26,6 +28,7 @@ interface BaseWheelOptions {
   altKey: boolean;
   invertScroll?: boolean;
 }
+
 interface WheelOptions extends BaseWheelOptions {
   contourOptions: ContourOptions;
 }
@@ -40,6 +43,7 @@ const DEFAULT_CONTOURS_OPTIONS: ContourOptions = {
     numberOfLayers: 10,
   },
 };
+
 type LevelSign = keyof Level;
 
 const LEVEL_SIGNS: Readonly<[LevelSign, LevelSign]> = ['positive', 'negative'];
@@ -56,21 +60,43 @@ function getDefaultContoursLevel(spectrum: Spectrum2D, quadrant = 'rr') {
   // @ts-expect-error type of NmrData2D should have a discriminator field to separate fid and ft
   const quadrantData = data[quadrant];
 
-  //@ts-expect-error will be included in nexts versions
-  const { noise = calculateSanPlot('2D', quadrantData) } = info;
-
-  const { positive = 0, negative = 0 } = noise;
-  const max = Math.max(
-    Math.abs(quadrantData.minZ),
-    Math.abs(quadrantData.maxZ),
+  const {
+    experiment = '',
+    // @ts-expect-error will be included in nexts versions
+    noise,
+  } = info;
+  const matrix = matrixToArray(quadrantData.z);
+
+  const noiseLevel = (noise as { positive: number; negative: number })
+    ? Math.max(noise.positive, noise.negative)
+    : estimateNoiseLevel(matrix);
+
+  const bestMinLevel = findAutomaticContourLevels(
+    matrix,
+    quadrantData.z.length,
+    quadrantData.z[0].length,
+    noiseLevel,
+    {
+      maxFdr: 0.005,
+      diagnostics: false,
+      maxOccupancy: 0.005,
+      persistenceLevels: 5,
+      contourRatio: 1.8,
+      minPersistence: 0.5,
+      maxVerticalRidgeScore: 0.05,
+      maxHorizontalRidgeScore: 0.05,
+      ridgeCoverageThreshold: experiment.includes('jres') ? 0.9 : 0.1,
+    },
   );
 
-  const minAbsPeakBase = 0.005 * max;
-  const minAllowed = 3 * xMaxAbsoluteValue([positive, negative]);
-
-  const minLevel = Math.max(minAbsPeakBase, minAllowed);
+  const { hasT1Noise, minLevelWithoutT1Noise, minLevel, maxAbsoluteValue } =
+    bestMinLevel;
   const minContourLevel = Math.min(
-    calculateValueOfLevel(minLevel, max, true),
+    calculateValueOfLevel(
+      hasT1Noise ? minLevelWithoutT1Noise : minLevel,
+      maxAbsoluteValue,
+      true,
+    ),
     DEFAULT_CONTOURS_OPTIONS.positive.contourLevels[1] -
       DEFAULT_CONTOURS_OPTIONS.positive.numberOfLayers,
   );

src/data/data2d/Spectrum2D/findBestMinContour/estimateNoiseLevel.ts

@@ -0,0 +1,67 @@
+import { xMedian } from 'ml-spectra-processing';
+
+const GAUSSIAN_MAD_SCALE = 1 / 0.6744897501960817;
+const DEFAULT_MAX_SAMPLES = 65_536;
+
+export interface EstimateNoiseLevelOptions {
+  /**
+   * Maximum number of evenly spaced finite values used for estimation.
+   * @default 65536
+   */
+  maxSamples?: number;
+}
+
+/**
+ * Estimates the standard deviation of signed Gaussian noise using a
+ * deterministically downsampled median absolute deviation.
+ */
+export function estimateNoiseLevel(
+  matrix: ArrayLike<number>,
+  options: EstimateNoiseLevelOptions = {},
+): number {
+  const { maxSamples = DEFAULT_MAX_SAMPLES } = options;
+
+  if (!Number.isInteger(maxSamples) || maxSamples <= 0) {
+    throw new Error('maxSamples must be a positive integer');
+  }
+
+  const sample = collectFiniteSample(matrix, maxSamples);
+
+  if (sample.length === 0) {
+    throw new Error('matrix must contain at least one finite value');
+  }
+
+  const median = xMedian(sample);
+
+  if (median === 0) {
+    return 1; //assume NUS reconstruction.
+  }
+
+  const averageDeviations = new Float64Array(sample.length);
+  for (let i = 0; i < sample.length; i++) {
+    averageDeviations[i] = Math.abs(sample[i] - median);
+  }
+
+  const noiseLevel = xMedian(averageDeviations) * GAUSSIAN_MAD_SCALE;
+
+  if (!Number.isFinite(noiseLevel) || noiseLevel < 0) {
+    throw new Error('unable to estimate a positive noise level');
+  }
+
+  return noiseLevel;
+}
+
+function collectFiniteSample(
+  matrix: ArrayLike<number>,
+  maxSamples: number,
+): Float64Array {
+  const step = Math.max(1, Math.ceil(matrix.length / maxSamples));
+  const sample = new Float64Array(Math.ceil(matrix.length / step));
+  let sampleIndex = 0;
+
+  for (let index = 0; index < matrix.length; index += step) {
+    sample[sampleIndex++] = matrix[index];
+  }
+
+  return sample;
+}