gaussian filter produces pixels outside valid range

[ghost]

Description

filters.gaussian produces pixel values outside of the valid range, -1 to 1.

This unwanted behavior occurs for certain sigma and truncate input values when the input image type is uint8 or float64. The behavior does not occur when the input image type is float32.
I bet that the cause involves floating-point errors, but I am uncertain of how to address it.

Way to reproduce

The following code block reproduces the behavior and demonstrates how it can lead to errors in subsequent functions.

from skimage import data, filters,util
import numpy as np 

image = data.checkerboard()
blurimg = filters.gaussian(image,sigma=2.0)
print(blurimg.max()) # prints 1.0000000000000004
util.img_as_ubyte(blurimg) # will raise a ValueError

The ValueError raised by img_as_ubyte reads.

ValueError: Images of type float must be between -1 and 1.

I tested gaussian with varying sigma values, truncate values, and image types. That testing script is pasted below.
The defined function conducts the tests. It returns a list of $N$ failed input combinations. The code that follows the function definition, applies the function to three image types.

def find_failed_gaussians(image):
    # list of sigmas to test in gaussian filter
    sigma_list = np.linspace(1.0,9.0,num=97) 
    # list of truncate values to test in gaussian filter
    truncate_list = np.arange(1,6) 
    # list of inputs which fail the 'try: assert' block
    failed_inputs = []
    for sigma in sigma_list:
        for trunc in truncate_list:
            imgblur = filters.gaussian(image,sigma=gsigma,truncate=trunc)
            try:
                assert(imgblur.max()<=1.0)
            except:
                failed_inputs.append((gsigma,trunc))
    return failed_inputs

image = data.checkerboard()
print(image.dtype)
print(len(find_failed_gaussians(image)))

image = util.img_as_float64(data.checkerboard())
print(image.dtype)
print(len(find_failed_gaussians(image)))

image = util.img_as_float32(data.checkerboard())
print(image.dtype)
print(len(find_failed_gaussians(image)))

Version information

3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 07:04:59) [GCC 10.3.0]
Linux-5.18.10-76051810-generic-x86_64-with-glibc2.33
scikit-image version: 0.19.3
numpy version: 1.23.0

[lagru]

I think this is caused by the upstream (SciPy 1.9.0, NumPy 1.23.2) function scipy.ndimage.gaussian_filter1d that is called by scipy.ndimage.gaussian_filter and in turn skimage.filter.gaussian:

import scipy as sp
import skimage as ski
image = ski.img_as_float(ski.data.checkerboard())
sp.ndimage.gaussian_filter(image, sigma=2.).max()  # returns 1.0000000000000002
sp.ndimage.gaussian_filter1d(image.ravel(), sigma=2).max() # returns 1.0000000000000002

I am not sure how to best approach this, too. For the current state of skimage expecting scaled input at a lot of places, it is indeed an annoying and surprising problem! I don't think this is really buggy behavior on SciPy's side; floating / precision errors are a fact of life sadly. Maybe one could start digging into SciPy's internals if there is anything that can be tweaked but I'm not sure how likely that will help. Clipping the output or raising a warning is a performance cost we likely don't want to pay in general for this function.

Summoning the experience of @scikit-image/core for sage advice. 🧙

[lagru]

Oh, and I don't think this (from the docstring)

The multi-dimensional filter is implemented as a sequence of one-dimensional convolution filters. The intermediate arrays are stored in the same data type as the output. Therefore, for output types with a limited precision, the results may be imprecise because intermediate results may be stored with insufficient precision.

is the cause as this behavior also happens with the one-dimensional filter.

[github-actions]

Hey, there hasn't been any activity on this issue for more than 180 days. For now, we have marked it as "dormant" until there is some new activity. You are welcome to reach out to people by mentioning them here or on our forum if you need more feedback! If you think that this issue is no longer relevant, you may close it by yourself; otherwise, we may do it at some point (either way, it will be done manually). In any case, thank you for your contributions so far!

[jni]

Hi all, apologies for the silence here.

Here are potential suggestions for course of action:

• clip the output of gaussian filter to the min and max of the input, as this is a theoretical guarantee of the filter in the absence of floating point errors.
• clip the input to img_as_ubyte and warn, rather than raise a ValueError
• same as 2, but only when the min/max are within some tolerance, say 1e-7, otherwise raise as we do now

It's tempting to ignore the problem because we are going to do away with automatic scaling in the future, but, I don't think we can, because this is actually an issue with our manual scaling functions — ie even in the non-auto-scaling future, we still want users to be able to rescale their images and for that to be robust to floating point errors.

In fact, on writing that last paragraph, I think that we might want to do both option 1 and option 3.

CC @stefanv