filters.gaussian does not work on 1D input

[stefanv]

import numpy as np
from skimage import filters

x = np.arange(50)
y = filters.gaussian(np.arange(50), sigma=(1,))

np.all(y < 1e-10)  # True

[sciunto]

There is actually two different aspects.

• 1D case was not handled. Fixed in a PR.
• Arrays are converted to float, and preserve_range must be used to manipulate this conversion in the original example.

[michalkrawczyk]

Hi,
I was wondering if problem don't lies in img_as_float from convert_to_float function:
I run test for np.arrange(9) with that and:

• result on filters.gaussian(np.arange(9), sigma=(1,))

[9.36517148e-20 1.69488592e-19 2.71560066e-19 3.79485270e-19
4.87890978e-19 5.96296685e-19 7.04221889e-19 8.06293363e-19
8.82130240e-19]

• result on filters.gaussian(np.arange(9), sigma=1, preserve_range=True)

[0.36378461 1.06325634 2.00469952 3.00013383 4. 4.99986617
5.99530048 6.93674366 7.63621539]

I might be wrong but it seems to be ok. Going deeper:

• convert_to_float(np.arange(9), True)

[0. 1. 2. 3. 4. 5. 6. 7. 8.] # OK

• convert_to_float(np.arange(9), False)

[5.42101086e-20 1.62630326e-19 2.71050543e-19 3.79470760e-19
4.87890978e-19 5.96311195e-19 7.04731412e-19 8.13151629e-19
9.21571847e-19] # <- What is That?

[jni]

@michalkrawczyk, see this page:

https://scikit-image.org/docs/dev/user_guide/data_types.html

convert to float with preserve range = False renormalises based on the assumed range of the data type, as shown in the table on the top of that page. So in your case, you start with int64.

In [3]: 1 / np.iinfo(np.int64).max                                         
Out[3]: 1.0842021724855044e-19

That's not exactly your number but it's about right. I'm not sure why part of the int64 range is truncated here, but that's the gist.

[michalkrawczyk]

@jni
I guess I found why is 'truncated'. See skimage.util.dtype.convert:

# signed/unsigned int -> float
if kind_out == 'f':
    # use float type that can exactly represent input integers
    computation_type = _dtype_itemsize(itemsize_in, dtype_out,
                                       np.float32, np.float64)

    if kind_in == 'u':
        # using np.divide or np.multiply doesn't copy the data
        # until the computation time
        image = np.multiply(image, 1. / imax_in,
                            dtype=computation_type)
        # DirectX uses this conversion also for signed ints
        # if imin_in:
        #     np.maximum(image, -1.0, out=image)
    else:
        image = np.add(image, 0.5, dtype=computation_type)               # here
        image *= 2 / (imax_in - imin_in)

Edit v2:

It seems according to Data Conversion Rules (and that comment from code) that also int conversions should be normalized with
image = np.multiply(image, 1. / imax_in, dtype=computation_type) not from else block