gauss-error-based separated blur kernel

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

const erf01 = new Float32Array [ 
	1.0000000000000000000 ];

const erf03 = new Float32Array [ 
	0.8053219914436340000, 
	0.0973390042781829800 ];

const erf05 = new Float32Array(3) [ 
	0.5635007619857788000, 
	0.2084610462188720700, 
	0.0097885727882385250 ];

const erf07 = new Float32Array(4) [ 
	0.4216516315937042000, 
	0.2414563894271850600, 
	0.0449914149940013900, 
	0.0027263797819614410 ];

const erf09 = new Float32Array(5) [ 
	0.3344678878784179700, 
	0.2354270666837692300, 
	0.0820083096623420700, 
	0.0140916891396045680, 
	0.0012389905750751495 ];

const erf11 = new Float32Array(6) [ 
	0.2764289081096649000, 
	0.2174566239118576000, 
	0.1058361455798149100, 
	0.0318464599549770360, 
	0.0059180646203458310, 
	0.0007282518781721592 ];

const erf13 = new Float32Array(7) [ 
	0.2352707684040069600, 
	0.1977262347936630200, 
	0.1173612549901008600, 
	0.0491887070238590240, 
	0.0145529769361019130, 
	0.0030381297692656517, 
	0.0004973122850060463 ];

const erf15 = new Float32Array(8) [ 
	0.2046523839235305800, 
	0.1794241964817047000, 
	0.1209106147289276100, 
	0.0626238510012626600, 
	0.0249265767633914950, 
	0.0076238219626247880, 
	0.0017914142226800323, 
	0.0003733328776434064 ];