blur.Gaussian: radius is non-standard and not commensurate with sigma

[rcoreilly]

The Gaussian kernel used in blur.Gaussian, in blur/blur.go, is computed using a radius value which is non-standard relative to e.g., the scipy image blurring code, and doesn't map onto the sigma (standard deviation) parameter that defines a Gaussian.

The code uses:

math.Exp(-(x * x / 4 / radius))

But a standard gaussian is defined in terms of sigma, where sigma^2 goes in the denominator, and there is a 1/2 factor, not 1/4: https://en.wikipedia.org/wiki/Gaussian_blur

sigma2 := sigma * sigma
...
math.Exp(-(x * x) / 2 / sigma2))

And in the scipy implementation, you use a radius that is some multiple, default 4, of the sigma value, not sigma^2. Thus, there is no way to "square" the single radius parameter used in the current implementation with other standard implementations. It would be better to have an explicit radius multiplier, in addition to the sigma parameter.

Here's an alternative implementation (in a separate codebase, not a fork) that reproduces the scipy test results from here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html (except for apparent edge / rounding issues), and splits out the kernel for independent inspection during tests. If you'd accept a PR with this, I could do that. However, this would be a breaking change, so it is not clear how best to handle that?

// scipy impl:
// https://github.com/scipy/scipy/blob/4bfc152f6ee1ca48c73c06e27f7ef021d729f496/scipy/ndimage/filters.py#L136
// #L214 has the invocation: radius = Ceil(sigma)

// bild uses:
// math.Exp(-0.5 * (x * x / (2 * radius))
// so   sigma = sqrt(radius) / 2
// and radius = sigma * sigma * 2

// GaussianBlurKernel1D returns a 1D Gaussian kernel.
// Sigma is the standard deviation,
// and the radius of the kernel is 4 * sigma.
func GaussianBlurKernel1D(sigma float64) *convolution.Kernel {
	sigma2 := sigma * sigma
	sfactor := -0.5 / sigma2
	radius := math.Ceil(4 * sigma) // truncate = 4 in scipy
	length := 2*int(radius) + 1

	// Create the 1-d gaussian kernel
	k := convolution.NewKernel(length, 1)
	for i, x := 0, -radius; i < length; i, x = i+1, x+1 {
		k.Matrix[i] = math.Exp(sfactor * (x * x))
	}
	return k
}

// GaussianBlur returns a smoothly blurred version of the image using
// a Gaussian function. Sigma is the standard deviation of the Gaussian
// function, and a kernel of radius = 4 * Sigma is used.
func GaussianBlur(src image.Image, sigma float64) *image.RGBA {
	if sigma <= 0 {
		return clone.AsRGBA(src)
	}

	k := GaussianBlurKernel1D(sigma).Normalized()

	// Perform separable convolution
	options := convolution.Options{Bias: 0, Wrap: false, KeepAlpha: false}
	result := convolution.Convolve(src, k, &options)
	result = convolution.Convolve(result, k.Transposed(), &options)

	return result
}