fix gaussian factor

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;
 }
 
 /**

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);
 });

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,15 +13,15 @@ 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;
 
   expect(data[center]).toBe(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];
 

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);
 

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) {

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;