cgcostume · November 18, 2020 08:23 · cgcostume · Nov 18, 2020
diff --git a/erf-kernel.js b/erf-kernel.js

 factorize = (x) => x == 0 ? 1 : x * factorize(x - 1);

 // Computation of the gauss-error (which gives the accumulated area below the normal 
 // distribution) using its Maclaurin series (https://en.wikipedia.org/wiki/Error_function)

 erf = (z, n) => {

    const lower = (n) => factorize(n) * (2.0 * n + 1.0);
    const upper = (n, z) => Math.pow(-1, n) * Math.pow(z, 2.0 * n + 1.0);

    let sum = 0.0;
    for(let i = 0; i <= n; ++i) {
        sum += upper(i, z) / lower(i);
    }

    // area of normal distribution is sqrt(PI) and erf is single-sided
    return sum * 2.0 / Math.sqrt(Math.PI); 
 }

 weight = (i, k) => {
    // some sigma to scale kernel size into
    const sigma = 2.75106392; // empirically done, trying to get a rest area of < 0.0001
    
    // calculate the number of buckets
    const kreal = 2.0 * k + 1;
  
    // start and end position of ith bucket
    const x0 = (i < 1 ? 0.0 : (-1.0 + 2.0 * i) / kreal * sigma);  
    const x1 = (+1.0 + 2.0 * i) / kreal * sigma;

    // areas for both positions
    const p0 = erf(x0, 32);
    const p1 = erf(x1, 32);
  
    // difference is area within bucket, note that all except the first buckets are applied twice
    return (p1 - p0) * (i < 1 ? 1.0 : 0.5);
 }

 erfKernel = (n) => {

    const size = Math.max(0, n);
    const weights = new Float32Array(size + 1);

    let sum = 0.0; // check sum for remaining area computation
   
    for(let i = 0; i <= size; ++i) {
        weights[i] = weight(i, size);
     	sum += weights[i] * (i > 0 ? 2.0 : 1.0);
    }
   
    // account for the remaining area (should be approximatelly the 
    // single-sided error of 0.0001 tweaked towards above ...)
    weights[size] += (1.0 - sum) * (size > 0 ? 0.5 : 1.0);

    return weights;
 }

 // TEST: area of normal distribution is sqrt(PI) and erf is single-sided
 for(let i = 0; i < 8; ++i) {

    const k = i * 2 + 1;
    const kernel = erfKernel(i);
    const checksum = kernel.reduce((s, x, i) => s + x * (i > 0 ? 2.0 : 1.0));

    console.log(`erf${String(k).padStart(2, '0')} (checksum = ${checksum}):`)
    console.log(kernel);
 }

	factorize = (x) => x == 0 ? 1 : x * factorize(x - 1);

	// Computation of the gauss-error (which gives the accumulated area below the normal
	// distribution) using its Maclaurin series (https://en.wikipedia.org/wiki/Error_function)

	erf = (z, n) => {

	const lower = (n) => factorize(n) * (2.0 * n + 1.0);
	const upper = (n, z) => Math.pow(-1, n) * Math.pow(z, 2.0 * n + 1.0);

	let sum = 0.0;
	for(let i = 0; i <= n; ++i) {
	sum += upper(i, z) / lower(i);
	}

	// area of normal distribution is sqrt(PI) and erf is single-sided
	return sum * 2.0 / Math.sqrt(Math.PI);
	}

	weight = (i, k) => {
	// some sigma to scale kernel size into
	const sigma = 2.75106392; // empirically done, trying to get a rest area of < 0.0001

	// calculate the number of buckets
	const kreal = 2.0 * k + 1;

	// start and end position of ith bucket
	const x0 = (i < 1 ? 0.0 : (-1.0 + 2.0 * i) / kreal * sigma);
	const x1 = (+1.0 + 2.0 * i) / kreal * sigma;

	// areas for both positions
	const p0 = erf(x0, 32);
	const p1 = erf(x1, 32);

	// difference is area within bucket, note that all except the first buckets are applied twice
	return (p1 - p0) * (i < 1 ? 1.0 : 0.5);
	}

	erfKernel = (n) => {

	const size = Math.max(0, n);
	const weights = new Float32Array(size + 1);

	let sum = 0.0; // check sum for remaining area computation

	for(let i = 0; i <= size; ++i) {
	weights[i] = weight(i, size);
	sum += weights[i] * (i > 0 ? 2.0 : 1.0);
	}

	// account for the remaining area (should be approximatelly the
	// single-sided error of 0.0001 tweaked towards above ...)
	weights[size] += (1.0 - sum) * (size > 0 ? 0.5 : 1.0);

	return weights;
	}

	// TEST: area of normal distribution is sqrt(PI) and erf is single-sided
	for(let i = 0; i < 8; ++i) {

	const k = i * 2 + 1;
	const kernel = erfKernel(i);
	const checksum = kernel.reduce((s, x, i) => s + x * (i > 0 ? 2.0 : 1.0));

	console.log(`erf${String(k).padStart(2, '0')} (checksum = ${checksum}):`)
	console.log(kernel);
	}