From 337a2ca072289f656add1dda5736f8592bba2e13 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Tue, 28 Jul 2026 11:38:28 +0200 Subject: [PATCH] fix: get gaussian factor was overestimating the value --- src/shapes/1d/gaussian/Gaussian.ts | 8 +- .../gaussian/__tests__/Gaussian.fct.test.ts | 2 +- .../gaussian/__tests__/Gaussian.shape.test.ts | 74 +++++++++++++++---- src/shapes/2d/gaussian2D/Gaussian2D.ts | 21 +++++- .../__tests__/Gaussian2D.shapes.test.ts | 47 ++++++++---- src/util/constants.ts | 1 + 6 files changed, 120 insertions(+), 33 deletions(-) diff --git a/src/shapes/1d/gaussian/Gaussian.ts b/src/shapes/1d/gaussian/Gaussian.ts index 3c6a91d..593c983 100644 --- a/src/shapes/1d/gaussian/Gaussian.ts +++ b/src/shapes/1d/gaussian/Gaussian.ts @@ -1,6 +1,7 @@ import { GAUSSIAN_EXP_FACTOR, ROOT_2LN2, + ROOT_LN2, ROOT_PI_OVER_LN2, } from '../../../util/constants.ts'; import erfinv from '../../../util/erfinv.ts'; @@ -178,12 +179,15 @@ export function getGaussianArea(options: GetGaussianAreaOptions) { } /** - * Calculate the half-width factor corresponding to a given area coverage fraction. + * Calculate the width factor corresponding to a given area coverage fraction. * @param area - target area fraction (0–1). Defaults to `0.9999`. * @returns the factor by which to multiply fwhm to cover the given area. */ export function getGaussianFactor(area = 0.9999) { - return Math.sqrt(2) * erfinv(area); + if (area >= 1) { + throw new Error('area should be (0 - 1)'); + } + return erfinv(area) / ROOT_LN2; } /** diff --git a/src/shapes/1d/gaussian/__tests__/Gaussian.fct.test.ts b/src/shapes/1d/gaussian/__tests__/Gaussian.fct.test.ts index 0160126..b6c9aa1 100644 --- a/src/shapes/1d/gaussian/__tests__/Gaussian.fct.test.ts +++ b/src/shapes/1d/gaussian/__tests__/Gaussian.fct.test.ts @@ -17,5 +17,5 @@ test('new Gaussian() gaussian.fct', () => { }); test('Guassian functions', () => { - expect(getGaussianFactor()).toBeCloseTo(3.8833175701198104); + expect(getGaussianFactor()).toBeCloseTo(3.298186269744253); }); diff --git a/src/shapes/1d/gaussian/__tests__/Gaussian.shape.test.ts b/src/shapes/1d/gaussian/__tests__/Gaussian.shape.test.ts index 415c455..ffb3a6c 100644 --- a/src/shapes/1d/gaussian/__tests__/Gaussian.shape.test.ts +++ b/src/shapes/1d/gaussian/__tests__/Gaussian.shape.test.ts @@ -1,7 +1,7 @@ -import erfinv from 'compute-erfinv'; import { expect, test } from 'vitest'; import { ROOT_2LN2, ROOT_PI_OVER_LN2 } from '../../../../util/constants.ts'; +import { pseudoVoigtFindFactor } from '../../pseudoVoigt/computeFactor.ts'; import { Gaussian, calculateGaussianHeight, @@ -13,7 +13,7 @@ test('height 1', () => { const gaussian = new Gaussian({ fwhm: 10 }); const data = gaussian.getData({ height: 1 }); - expect(data).toHaveLength(39); + expect(data).toHaveLength(33); const center = (data.length - 1) / 2; @@ -21,7 +21,7 @@ test('height 1', () => { const area = data.reduce((a, b) => a + b, 0); - expect(area).toBeCloseTo((ROOT_PI_OVER_LN2 * 10) / 2, 3); + expect(area).toBeCloseTo(((ROOT_PI_OVER_LN2 * 10) / 2) * 0.9999, 3); expect(gaussian.getParameters()).toStrictEqual(['fwhm']); }); @@ -37,11 +37,11 @@ test('fwhm fixed and normalized', () => { const gaussian = new Gaussian({ fwhm: 50 }); const data = gaussian.getData(); - expect(data).toHaveLength(195); + expect(data).toHaveLength(165); const area = data.reduce((a, b) => a + b, 0); - expect(area).toBeCloseTo(0.9999, 2); + expect(area).toBeCloseTo(0.9999, 4); const computedArea = gaussian.getArea(); @@ -59,7 +59,7 @@ test('sd fixed', () => { const area = data.reduce((a, b) => a + b, 0); - expect(area).toBeCloseTo(height * Math.sqrt(2 * Math.PI) * sd, 2); + expect(area).toBeCloseTo(height * Math.sqrt(2 * Math.PI) * sd * 0.9999, 2); }); test('odd fwhm', () => { @@ -125,18 +125,64 @@ test('height calculations', () => { expect(height).toBeCloseTo(expectedHeight, 4); }); -test('factor should be close', () => { - const gaussian = new Gaussian({ fwhm: 100 }); - for (let i = 1; i < 11; i++) { - const area = i * 0.1; +test('factor covers the requested area', () => { + const gaussian = new Gaussian({ fwhm: 1000 }); + const areas = [0.98, 0.96, 0.7, 0.4, 0.2]; + for (const area of areas) { + const data = gaussian.getData({ factor: gaussian.getFactor(area) }); + const sum = data.reduce((a, b) => a + b, 0); + + expect(sum).toBeCloseTo(area, 2); + } +}); - expect(gaussian.getFactor(area)).toBeCloseTo( - Math.sqrt(2) * erfinv(area), - 1, - ); +test('factor integrates to the requested area', () => { + const gaussian = new Gaussian({ fwhm: 1 }); + const total = gaussian.getArea(1); + const areas = [0.2, 0.4, 0.7, 0.9, 0.96, 0.98, 0.995, 0.9999]; + for (const area of areas) { + const halfWidth = gaussian.getFactor(area) / 2; + + expect(integrate(gaussian, halfWidth) / total).toBeCloseTo(area, 3); } }); +test('factor is consistent with the pseudo-Voigt gaussian limit', () => { + const gaussian = new Gaussian({ fwhm: 1 }); + const areas = [0.2, 0.7, 0.98]; + for (const area of areas) { + const factor = gaussian.getFactor(area); + const bisected = pseudoVoigtFindFactor(area, 0.999999); + + // the residual is the accuracy of the erfinv approximation + expect(Math.abs(factor - bisected) / factor).toBeLessThan(1e-3); + } +}); + +test('factor throws when the area cannot be reached', () => { + const gaussian = new Gaussian({ fwhm: 1 }); + + expect(() => gaussian.getFactor(1)).toThrow('area should be (0 - 1)'); + expect(() => gaussian.getFactor(1.2)).toThrow('area should be (0 - 1)'); +}); + +/** + * Numerically integrate the shape over [-halfWidth, halfWidth] with the + * trapezoidal rule, independently of any closed-form area formula. + * @param gaussian - the shape to integrate. + * @param halfWidth - half of the integration window. + * @returns the integral over the window. + */ +function integrate(gaussian: Gaussian, halfWidth: number) { + const nbSteps = 100000; + const step = (2 * halfWidth) / nbSteps; + let sum = (gaussian.fct(-halfWidth) + gaussian.fct(halfWidth)) / 2; + for (let index = 1; index < nbSteps; index++) { + sum += gaussian.fct(-halfWidth + index * step); + } + return sum * step; +} + function getNbChanges(y: Float64Array) { const yPrime = [0]; diff --git a/src/shapes/2d/gaussian2D/Gaussian2D.ts b/src/shapes/2d/gaussian2D/Gaussian2D.ts index 2bebd01..9871672 100644 --- a/src/shapes/2d/gaussian2D/Gaussian2D.ts +++ b/src/shapes/2d/gaussian2D/Gaussian2D.ts @@ -88,8 +88,8 @@ export class Gaussian2D implements Shape2DClass { ); } - public getFactor(volume = 1) { - return getGaussianFactor(volume); + public getFactor(volume?: number) { + return getGaussian2DFactor(volume); } public getVolume( @@ -126,6 +126,21 @@ export class Gaussian2D implements Shape2DClass { } } +/** + * Calculate the width factor corresponding to a given volume coverage fraction. + * The factor delimits a rectangular window, and the shape is separable, so the + * enclosed volume is the product of the coverage of both axes. Each axis + * therefore has to cover the square root of the target volume. + * @param volume - target volume fraction (0–1). Defaults to `0.9999`. + * @returns the factor by which to multiply fwhm, on each axis, to cover the given volume. + */ +export const getGaussian2DFactor = (volume = 0.9999) => { + if (volume >= 1) { + throw new Error('volume should be (0 - 1)'); + } + return getGaussianFactor(Math.sqrt(volume)); +}; + export const gaussian2DFct = ( x: number, y: number, @@ -145,7 +160,7 @@ export const getGaussian2DData = ( fwhm = ensureFWHM2D(fwhm, sd); const { height = calculateGaussian2DHeight({ fwhm, volume: 1 }) } = options; - let { factor = getGaussianFactor(), length = { x: 0, y: 0 } } = options; + let { factor = getGaussian2DFactor(), length = { x: 0, y: 0 } } = options; factor = ensureXYNumber(factor); diff --git a/src/shapes/2d/gaussian2D/__tests__/Gaussian2D.shapes.test.ts b/src/shapes/2d/gaussian2D/__tests__/Gaussian2D.shapes.test.ts index 3ef9957..a69c0c5 100644 --- a/src/shapes/2d/gaussian2D/__tests__/Gaussian2D.shapes.test.ts +++ b/src/shapes/2d/gaussian2D/__tests__/Gaussian2D.shapes.test.ts @@ -1,10 +1,10 @@ -import erfinv from 'compute-erfinv'; import { expect, test } from 'vitest'; import { ROOT_2LN2 } from '../../../../util/constants.ts'; import { gaussianFwhmToWidth, gaussianWidthToFWHM, + getGaussianFactor, } from '../../../1d/gaussian/Gaussian.ts'; import { Gaussian2D } from '../Gaussian2D.ts'; @@ -14,7 +14,7 @@ test('height 1', () => { }); const data = gaussian2D.getData({ height: 1 }); - expect(data).toHaveLength(39); + expect(data).toHaveLength(35); const xCenter = (data.length - 1) / 2; const yCenter = (data[0].length - 1) / 2; @@ -23,7 +23,7 @@ test('height 1', () => { const volume = getVolume(data); - expect(volume).toBeCloseTo((100 * Math.PI) / Math.LN2 / 4, 2); + expect(volume).toBeCloseTo(((100 * Math.PI) / Math.LN2 / 4) * 0.9999, 2); }); test('check gaussian2D continuous', () => { @@ -38,11 +38,11 @@ test('fwhm fixed and normalized', () => { const gaussian2D = new Gaussian2D({ fwhm: 50 }); const data = gaussian2D.getData(); - expect(data).toHaveLength(195); + expect(data).toHaveLength(173); const volume = getVolume(data); - expect(volume).toBeCloseTo(0.9999, 2); + expect(volume).toBeCloseTo(0.9999, 4); const computeSurface = gaussian2D.getVolume(); @@ -62,7 +62,7 @@ test('sd fixed', () => { const volume = getVolume(data); - expect(volume).toBeCloseTo(height * 2 * Math.PI * sd * sd, 0); + expect(volume).toBeCloseTo(height * 2 * Math.PI * sd * sd * 0.9999, 0); }); test('odd fwhm', () => { @@ -118,18 +118,39 @@ test('change height should change area', () => { expect(gaussian2D.getVolume(2)).toBeCloseTo(2 * volume, 4); }); -test('factor should be close', () => { +test('factor covers the requested volume', () => { const gaussian2D = new Gaussian2D({ fwhm: 100 }); - for (let i = 1; i < 11; i++) { - const volume: number = i * 0.1; + const volumes = [0.98, 0.96, 0.7, 0.4, 0.2]; + for (const volume of volumes) { + const data = gaussian2D.getData({ + factor: gaussian2D.getFactor(volume), + }); - expect(gaussian2D.getFactor(volume)).toBeCloseTo( - Math.sqrt(2) * erfinv(volume), - 1, - ); + expect(getVolume(data)).toBeCloseTo(volume, 2); } }); +test('factor covers the square root of the volume on each axis', () => { + const gaussian2D = new Gaussian2D({ fwhm: 100 }); + + expect(gaussian2D.getFactor(0.9)).toBeCloseTo( + getGaussianFactor(Math.sqrt(0.9)), + 10, + ); +}); + +test('default factor is finite', () => { + const gaussian2D = new Gaussian2D({ fwhm: 100 }); + + expect(gaussian2D.getFactor()).toBeCloseTo(3.4383781513674276); +}); + +test('factor throws when the volume cannot be reached', () => { + const gaussian2D = new Gaussian2D({ fwhm: 100 }); + + expect(() => gaussian2D.getFactor(1)).toThrow('volume should be (0 - 1)'); +}); + function getVolume(data: Float64Array[]) { let volume = 0; for (const row of data) { diff --git a/src/util/constants.ts b/src/util/constants.ts index 5ee458a..f6982c4 100644 --- a/src/util/constants.ts +++ b/src/util/constants.ts @@ -1,5 +1,6 @@ export const GAUSSIAN_EXP_FACTOR = -4 * Math.LN2; export const ROOT_PI_OVER_LN2 = Math.sqrt(Math.PI / Math.LN2); +export const ROOT_LN2 = Math.sqrt(Math.LN2); export const ROOT_THREE = Math.sqrt(3); export const ROOT_2LN2 = Math.sqrt(2 * Math.LN2); export const ROOT_2LN2_MINUS_ONE = Math.sqrt(2 * Math.LN2) - 1;