From 2b12e1631e76f273bbe6450171885e03864b30c7 Mon Sep 17 00:00:00 2001 From: jobo322 Date: Tue, 4 Aug 2026 10:13:48 -0500 Subject: [PATCH 1/4] feat(baselineRecognition): implement automatic baseline recognition function and corresponding tests --- .../xreimAutomaticBaselineRecognition.test.ts | 27 +++ .../xreimAutomaticBaselineRecognition.ts | 224 ++++++++++++++++++ 2 files changed, 251 insertions(+) create mode 100644 src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts create mode 100644 src/xreim/xreimAutomaticBaselineRecognition.ts diff --git a/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts b/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts new file mode 100644 index 000000000..465a0f516 --- /dev/null +++ b/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts @@ -0,0 +1,27 @@ +import { expect, test } from 'vitest'; + +import { xreimAutomaticBaselineRecognition } from '../xreimAutomaticBaselineRecognition.ts'; + +test('detects baseline regions around a Gaussian peak', () => { + const length = 201; + const x = new Float64Array(length); + const re = new Float64Array(length); + const im = new Float64Array(length); + + for (let i = 0; i < length; i++) { + x[i] = i; + re[i] = 0.05 * i + 25 * Math.exp(-((i - 100) ** 2) / (2 * 15 ** 2)); + } + + const mask = xreimAutomaticBaselineRecognition( + { x, re, im }, + { scale: 4, thresholdFactor: 0.5, erosionRadius: 0, component: 're' }, + ); + + const maskArray = Array.from(mask); + expect(maskArray[0]).toBe(1); + expect(maskArray[200]).toBe(1); + expect(maskArray[100]).toBe(0); + expect(maskArray.some((value) => value === 1)).toBe(true); + expect(maskArray.some((value) => value === 0)).toBe(true); +}); diff --git a/src/xreim/xreimAutomaticBaselineRecognition.ts b/src/xreim/xreimAutomaticBaselineRecognition.ts new file mode 100644 index 000000000..03992fa67 --- /dev/null +++ b/src/xreim/xreimAutomaticBaselineRecognition.ts @@ -0,0 +1,224 @@ +import type { DoubleArray } from 'cheminfo-types'; + +import type { DataXReIm } from '../types/index.ts'; + +export interface AutomaticBaselineRecognitionOptions { + /** + * Scale parameter for the discrete CWT-Haar derivative approximation. + * Larger values reduce noise but widen the detected peaks. + * A value of 'auto' computes a heuristic scale based on the spectrum length. + * @default 'auto' + */ + scale?: number | 'auto'; + + /** + * Multiplicative factor used in the iterative thresholding rule. + * The threshold is the mean plus factor × standard deviation. + * A value close to 0.5 is often more robust for the unnormalized CWT derivative + * used here than the 3 used in the original Dietrich procedure. + * @default 0.5 + */ + thresholdFactor?: number; + + /** + * Radius of the 1D erosion filter used to remove isolated spikes in the mask. + * @default 1 + */ + erosionRadius?: number; + + /** + * Signal component used for the analysis. + * - 're': use the real component. + * - 'im': use the imaginary component. + * - 'magnitude': use the magnitude of the complex signal. + * @default 're' + */ + component?: 're' | 'im' | 'magnitude'; + + /** + * Deprecated alias for component. + */ + mode?: 're' | 'im' | 'magnitude'; +} + +/** + * Automatically detects signal-free regions in a 1D spectrum using an + * approximate CWT-Haar derivative followed by iterative thresholding. + * + * The output is a binary mask where 1 marks points belonging to the baseline + * (signal-free regions) and 0 marks points associated with peaks or signal. + * + * @param data - object of kind {x:[], re:[], im:[]} + * @param options - recognition options + * @returns a binary mask as a Uint8Array + */ +export function xreimAutomaticBaselineRecognition< + ArrayType extends DoubleArray = DoubleArray, +>( + data: DataXReIm, + options: AutomaticBaselineRecognitionOptions = {}, +): Uint8Array { + const { + scale = 'auto', + thresholdFactor = 0.5, + erosionRadius = 1, + component, + mode, + } = options; + + const length = data.x.length; + if (data.re.length !== length || data.im.length !== length) { + throw new TypeError('length of x, re and im must be identical'); + } + + const signal = getSignal(data, component ?? mode ?? 're'); + const actualScale = resolveScale(length, scale); + const derivative = computeCwtHaarDerivative(signal, actualScale); + const power = new Float64Array(length); + for (let i = 0; i < length; i++) { + const value = derivative[i]; + power[i] = value * value; + } + + const threshold = iterativeThreshold(power, thresholdFactor); + const mask = new Uint8Array(length); + for (let i = 0; i < length; i++) { + mask[i] = power[i] <= threshold ? 1 : 0; + } + + if (erosionRadius > 0) { + return erodeMask(mask, erosionRadius); + } + return mask; +} + +function getSignal( + data: DataXReIm, + component: 're' | 'im' | 'magnitude', +): DoubleArray { + const { re, im } = data; + switch (component) { + case 'im': + return im; + case 'magnitude': { + const magnitude = new Float64Array(re.length); + for (let i = 0; i < re.length; i++) { + magnitude[i] = Math.hypot(re[i], im[i]); + } + return magnitude; + } + case 're': + return re; + default: + return re; + } +} + +function resolveScale( + length: number, + scale: AutomaticBaselineRecognitionOptions['scale'], +): number { + if (typeof scale === 'number') { + return Math.max( + 1, + Math.min(Math.floor(scale), Math.floor((length - 1) / 2)), + ); + } + return Math.max(1, Math.floor(length / 512)); +} + +function computeCwtHaarDerivative( + signal: DoubleArray, + scale: number, +): Float64Array { + const length = signal.length; + const derivative = new Float64Array(length); + + for (let i = scale; i < length - scale; i++) { + let sum = 0; + for (let j = 1; j <= scale; j++) { + sum -= signal[i - j]; + sum += signal[i + j]; + } + derivative[i] = sum / scale; + } + + return derivative; +} + +function iterativeThreshold(values: Float64Array, factor: number): number { + let threshold = getThreshold(values, factor); + let previousThreshold = Number.POSITIVE_INFINITY; + + while (Math.abs(previousThreshold - threshold) > 1e-12) { + const valuesBelow = new Float64Array(values.length); + let count = 0; + for (let i = 0; i < values.length; i++) { + if (values[i] <= threshold) { + valuesBelow[count++] = values[i]; + } + } + + if (count === 0) { + return threshold; + } + + previousThreshold = threshold; + threshold = getThreshold(valuesBelow.subarray(0, count), factor); + } + + return threshold; +} + +function getThreshold(values: DoubleArray, factor: number): number { + const mean = getMean(values); + const std = getStandardDeviation(values, mean); + return mean + factor * std; +} + +function getMean(values: DoubleArray): number { + let sum = 0; + for (let i = 0; i < values.length; i++) { + sum += values[i]; + } + return values.length === 0 ? 0 : sum / values.length; +} + +function getStandardDeviation(values: DoubleArray, mean: number): number { + if (values.length < 2) { + return 0; + } + let sumSquared = 0; + for (let i = 0; i < values.length; i++) { + const diff = values[i] - mean; + sumSquared += diff * diff; + } + return Math.sqrt(sumSquared / values.length); +} + +function erodeMask(mask: Uint8Array, radius: number): Uint8Array { + const result = new Uint8Array(mask.length); + for (let i = 0; i < mask.length; i++) { + if (mask[i] === 0) { + continue; + } + + let trueCount = 1; + let falseCount = 0; + for (let offset = -radius; offset <= radius; offset++) { + if (offset === 0) continue; + const index = i + offset; + if (index < 0 || index >= mask.length) { + falseCount++; + } else if (mask[index] === 1) { + trueCount++; + } else { + falseCount++; + } + } + + result[i] = falseCount >= trueCount ? 0 : 1; + } + + return result; +} From 8bd80d67912a280767e70575bdf3505cb5660dae Mon Sep 17 00:00:00 2001 From: jobo322 Date: Fri, 7 Aug 2026 10:50:57 -0500 Subject: [PATCH 2/4] feat(xreim): export automatic baseline recognition function and update dependencies --- src/xreim/index.ts | 1 + .../xreimAutomaticBaselineRecognition.ts | 40 ++++--------------- 2 files changed, 8 insertions(+), 33 deletions(-) diff --git a/src/xreim/index.ts b/src/xreim/index.ts index f987f0711..a7c71ea2a 100644 --- a/src/xreim/index.ts +++ b/src/xreim/index.ts @@ -1,2 +1,3 @@ export * from './xreimSortX.ts'; export * from './xreimZeroFilling.ts'; +export * from './xreimAutomaticBaselineRecognition.ts'; diff --git a/src/xreim/xreimAutomaticBaselineRecognition.ts b/src/xreim/xreimAutomaticBaselineRecognition.ts index 03992fa67..80d428c14 100644 --- a/src/xreim/xreimAutomaticBaselineRecognition.ts +++ b/src/xreim/xreimAutomaticBaselineRecognition.ts @@ -1,5 +1,6 @@ import type { DoubleArray } from 'cheminfo-types'; +import { xMean, xStandardDeviation } from '../index.ts'; import type { DataXReIm } from '../types/index.ts'; export interface AutomaticBaselineRecognitionOptions { @@ -34,11 +35,6 @@ export interface AutomaticBaselineRecognitionOptions { * @default 're' */ component?: 're' | 'im' | 'magnitude'; - - /** - * Deprecated alias for component. - */ - mode?: 're' | 'im' | 'magnitude'; } /** @@ -47,7 +43,6 @@ export interface AutomaticBaselineRecognitionOptions { * * The output is a binary mask where 1 marks points belonging to the baseline * (signal-free regions) and 0 marks points associated with peaks or signal. - * * @param data - object of kind {x:[], re:[], im:[]} * @param options - recognition options * @returns a binary mask as a Uint8Array @@ -63,7 +58,6 @@ export function xreimAutomaticBaselineRecognition< thresholdFactor = 0.5, erosionRadius = 1, component, - mode, } = options; const length = data.x.length; @@ -71,7 +65,7 @@ export function xreimAutomaticBaselineRecognition< throw new TypeError('length of x, re and im must be identical'); } - const signal = getSignal(data, component ?? mode ?? 're'); + const signal = getSignal(data, component ?? 're'); const actualScale = resolveScale(length, scale); const derivative = computeCwtHaarDerivative(signal, actualScale); const power = new Float64Array(length); @@ -153,9 +147,9 @@ function iterativeThreshold(values: Float64Array, factor: number): number { while (Math.abs(previousThreshold - threshold) > 1e-12) { const valuesBelow = new Float64Array(values.length); let count = 0; - for (let i = 0; i < values.length; i++) { - if (values[i] <= threshold) { - valuesBelow[count++] = values[i]; + for (const value of values) { + if (value <= threshold) { + valuesBelow[count++] = value; } } @@ -171,31 +165,11 @@ function iterativeThreshold(values: Float64Array, factor: number): number { } function getThreshold(values: DoubleArray, factor: number): number { - const mean = getMean(values); - const std = getStandardDeviation(values, mean); + const mean = xMean(values); + const std = xStandardDeviation(values, { mean }); return mean + factor * std; } -function getMean(values: DoubleArray): number { - let sum = 0; - for (let i = 0; i < values.length; i++) { - sum += values[i]; - } - return values.length === 0 ? 0 : sum / values.length; -} - -function getStandardDeviation(values: DoubleArray, mean: number): number { - if (values.length < 2) { - return 0; - } - let sumSquared = 0; - for (let i = 0; i < values.length; i++) { - const diff = values[i] - mean; - sumSquared += diff * diff; - } - return Math.sqrt(sumSquared / values.length); -} - function erodeMask(mask: Uint8Array, radius: number): Uint8Array { const result = new Uint8Array(mask.length); for (let i = 0; i < mask.length; i++) { From b3fde1ef0e0f22677042cc4204fc17b66b29f70e Mon Sep 17 00:00:00 2001 From: jobo322 Date: Tue, 25 Aug 2026 13:55:10 -0500 Subject: [PATCH 3/4] chore: update snapshot --- src/__tests__/__snapshots__/index.test.ts.snap | 1 + 1 file changed, 1 insertion(+) diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index b7eb3cb7d..73653ede1 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -136,6 +136,7 @@ exports[`existence of exported functions 1`] = ` "xy2ToXY", "xreimSortX", "xreimZeroFilling", + "xreimAutomaticBaselineRecognition", "xyArrayAlign", "xyArrayAlignToFirst", "xyArrayMerge", From f5a291e9ad042953f6479521e3317f0b509e03f1 Mon Sep 17 00:00:00 2001 From: jobo322 Date: Wed, 26 Aug 2026 09:54:25 -0500 Subject: [PATCH 4/4] feat: improve speed from 2 second to 96 ms for 1D spectra --- src/utils/__tests__/createRandomArray.test.ts | 1 - .../xreimAutomaticBaselineRecognition.test.ts | 6 +- .../xreimAutomaticBaselineRecognition.ts | 210 +++++++++++------- 3 files changed, 131 insertions(+), 86 deletions(-) diff --git a/src/utils/__tests__/createRandomArray.test.ts b/src/utils/__tests__/createRandomArray.test.ts index 907fdd37c..0505160e2 100644 --- a/src/utils/__tests__/createRandomArray.test.ts +++ b/src/utils/__tests__/createRandomArray.test.ts @@ -63,7 +63,6 @@ test('Testing in conjunction with spectra-fitting', () => { }); expect(fittedPeaks.peaks[0].x).toBeDeepCloseTo(10, 2); - //@ts-expect-error it is a gaussian shape expect(fittedPeaks.peaks[0].shape.fwhm).toBeDeepCloseTo( 2 * Math.sqrt(2 * Math.log(2)), 1, diff --git a/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts b/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts index 465a0f516..c53799c2a 100644 --- a/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts +++ b/src/xreim/__tests__/xreimAutomaticBaselineRecognition.test.ts @@ -19,9 +19,11 @@ test('detects baseline regions around a Gaussian peak', () => { ); const maskArray = Array.from(mask); + expect(maskArray[0]).toBe(1); expect(maskArray[200]).toBe(1); expect(maskArray[100]).toBe(0); - expect(maskArray.some((value) => value === 1)).toBe(true); - expect(maskArray.some((value) => value === 0)).toBe(true); + + expect(maskArray).toContain(1); + expect(maskArray).toContain(0); }); diff --git a/src/xreim/xreimAutomaticBaselineRecognition.ts b/src/xreim/xreimAutomaticBaselineRecognition.ts index 80d428c14..a416b133e 100644 --- a/src/xreim/xreimAutomaticBaselineRecognition.ts +++ b/src/xreim/xreimAutomaticBaselineRecognition.ts @@ -1,6 +1,5 @@ import type { DoubleArray } from 'cheminfo-types'; -import { xMean, xStandardDeviation } from '../index.ts'; import type { DataXReIm } from '../types/index.ts'; export interface AutomaticBaselineRecognitionOptions { @@ -47,34 +46,43 @@ export interface AutomaticBaselineRecognitionOptions { * @param options - recognition options * @returns a binary mask as a Uint8Array */ -export function xreimAutomaticBaselineRecognition< - ArrayType extends DoubleArray = DoubleArray, ->( - data: DataXReIm, +export function xreimAutomaticBaselineRecognition( + data: Omit & { im?: DoubleArray }, options: AutomaticBaselineRecognitionOptions = {}, ): Uint8Array { const { scale = 'auto', thresholdFactor = 0.5, erosionRadius = 1, - component, + component = 're', } = options; const length = data.x.length; - if (data.re.length !== length || data.im.length !== length) { + if (data.re.length !== length || (data.im && data.im.length !== length)) { throw new TypeError('length of x, re and im must be identical'); } - const signal = getSignal(data, component ?? 're'); + if (!data.im && component !== 're') { + throw new TypeError( + `component '${component}' requires im array to be defined`, + ); + } + + const signal = getSignal({ im: data.re, ...data }, component); const actualScale = resolveScale(length, scale); + + // OPTIMIZATION 1: O(N) derivative computation const derivative = computeCwtHaarDerivative(signal, actualScale); + const power = new Float64Array(length); for (let i = 0; i < length; i++) { const value = derivative[i]; power[i] = value * value; } + // OPTIMIZATION 2: Iterative threshold without array allocations const threshold = iterativeThreshold(power, thresholdFactor); + const mask = new Uint8Array(length); for (let i = 0; i < length; i++) { mask[i] = power[i] <= threshold ? 1 : 0; @@ -86,113 +94,149 @@ export function xreimAutomaticBaselineRecognition< return mask; } -function getSignal( - data: DataXReIm, - component: 're' | 'im' | 'magnitude', -): DoubleArray { - const { re, im } = data; - switch (component) { - case 'im': - return im; - case 'magnitude': { - const magnitude = new Float64Array(re.length); - for (let i = 0; i < re.length; i++) { - magnitude[i] = Math.hypot(re[i], im[i]); - } - return magnitude; - } - case 're': - return re; - default: - return re; - } -} - -function resolveScale( - length: number, - scale: AutomaticBaselineRecognitionOptions['scale'], -): number { - if (typeof scale === 'number') { - return Math.max( - 1, - Math.min(Math.floor(scale), Math.floor((length - 1) / 2)), - ); - } - return Math.max(1, Math.floor(length / 512)); -} - function computeCwtHaarDerivative( signal: DoubleArray, scale: number, ): Float64Array { const length = signal.length; const derivative = new Float64Array(length); + if (scale <= 0) return derivative; - for (let i = scale; i < length - scale; i++) { - let sum = 0; - for (let j = 1; j <= scale; j++) { - sum -= signal[i - j]; - sum += signal[i + j]; - } - derivative[i] = sum / scale; + // We want: sum_{j=1 to scale} (signal[i+j] - signal[i-j]) + // This is: (sum of right window) - (sum of left window) + let leftSum = 0; + let rightSum = 0; + + // Initialize windows for the first valid i (i = scale) + for (let j = 1; j <= scale; j++) { + leftSum += signal[scale - j]; + rightSum += signal[scale + j]; + } + derivative[scale] = (rightSum - leftSum) / scale; + + for (let i = scale + 1; i < length - scale; i++) { + // Slide windows: subtract the element leaving and add the element entering + leftSum = leftSum - signal[i - scale - 1] + signal[i - 1]; + rightSum = rightSum - signal[i + 1] + signal[i + scale + 1]; + derivative[i] = (rightSum - leftSum) / scale; } return derivative; } function iterativeThreshold(values: Float64Array, factor: number): number { - let threshold = getThreshold(values, factor); + let threshold = calculateThresholdFiltered( + values, + factor, + Number.POSITIVE_INFINITY, + ); let previousThreshold = Number.POSITIVE_INFINITY; - while (Math.abs(previousThreshold - threshold) > 1e-12) { - const valuesBelow = new Float64Array(values.length); - let count = 0; - for (const value of values) { - if (value <= threshold) { - valuesBelow[count++] = value; - } - } - - if (count === 0) { - return threshold; - } - + // Max iterations safety cap to prevent infinite loops in edge cases + let iterations = 0; + while (Math.abs(previousThreshold - threshold) > 1e-12 && iterations < 100) { previousThreshold = threshold; - threshold = getThreshold(valuesBelow.subarray(0, count), factor); + threshold = calculateThresholdFiltered(values, factor, threshold); + iterations++; } return threshold; } -function getThreshold(values: DoubleArray, factor: number): number { - const mean = xMean(values); - const std = xStandardDeviation(values, { mean }); +/** + * Calculates mean + factor * std, but only considers values <= currentThreshold + * This removes the need to create a new filtered array every iteration. + * @param values + * @param factor + * @param currentThreshold + */ +function calculateThresholdFiltered( + values: Float64Array, + factor: number, + currentThreshold: number, +): number { + let sum = 0; + let count = 0; + + for (const val of values) { + if (val <= currentThreshold) { + sum += val; + count++; + } + } + + if (count === 0) return currentThreshold; + + const mean = sum / count; + let varianceSum = 0; + + for (const val of values) { + if (val <= currentThreshold) { + const diff = val - mean; + varianceSum += diff * diff; + } + } + + const std = Math.sqrt(varianceSum / count); return mean + factor * std; } function erodeMask(mask: Uint8Array, radius: number): Uint8Array { - const result = new Uint8Array(mask.length); - for (let i = 0; i < mask.length; i++) { - if (mask[i] === 0) { - continue; - } + const length = mask.length; + const result = new Uint8Array(length); + const windowSize = 2 * radius + 1; + const threshold = windowSize / 2; - let trueCount = 1; - let falseCount = 0; + for (let i = 0; i < length; i++) { + if (mask[i] === 0) continue; + + let trueCount = 0; for (let offset = -radius; offset <= radius; offset++) { - if (offset === 0) continue; const index = i + offset; - if (index < 0 || index >= mask.length) { - falseCount++; - } else if (mask[index] === 1) { + if (index >= 0 && index < length && mask[index] === 1) { trueCount++; - } else { - falseCount++; } } - - result[i] = falseCount >= trueCount ? 0 : 1; + // Result is 1 only if majority of window is 1 + result[i] = trueCount > threshold ? 1 : 0; } return result; } + +function getSignal( + data: DataXReIm, + component: 're' | 'im' | 'magnitude', +): DoubleArray { + // Validate component when im is not provided + + const { re, im } = data; + switch (component) { + case 'im': + return im; + case 'magnitude': { + const magnitude = new Float64Array(re.length); + for (let i = 0; i < re.length; i++) { + magnitude[i] = Math.hypot(re[i], im[i]); + } + return magnitude; + } + case 're': + return re; + default: + return re; + } +} + +function resolveScale( + length: number, + scale: AutomaticBaselineRecognitionOptions['scale'], +): number { + if (typeof scale === 'number') { + return Math.max( + 1, + Math.min(Math.floor(scale), Math.floor((length - 1) / 2)), + ); + } + return Math.max(1, Math.floor(length / 512)); +}