Add float32 support to moments_hu

[rfezzani]

Description

Related to #5199.
Fused types are used to allow processing single precision values.

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[rfezzani]

I wonder if maintaining the cython implementation of moments_hu worth it 😕.

I made some tests by implementing a pure python equivalent function:

def py_moments_hu(x):
    offset = x.shape[1]
    _x = x.reshape(-1)

    x03, x12, x21, x30 = _x[3: 3 * offset + 1: offset - 1]
    t0 = x30 + x12
    t1 = x21 + x03
    q0 = t0 * t0
    q1 = t1 * t1

    x02, x11, x20 = _x[2: 2 * offset + 1: offset - 1]
    n4 = 4 * x11
    s = x20 + x02
    d = x20 - x02

    hu = np.empty((7, ), x.dtype)
    hu[0] = s
    hu[1] = d * d + n4 * x11

    hu[3] = q0 + q1
    hu[5] = d * (q0 - q1) + n4 * t0 * t1

    t0 *= q0 - 3 * q1
    t1 *= 3 * q0 - q1
    q0 = x30 - 3 * x12
    q1 = 3 * x21 - x03

    hu[2] = q0 * q0 + q1 * q1
    hu[4] = q0 * t0 + q1 * t1
    hu[6] = q1 * t0 - q0 * t1

    return hu

I then performed some timing:

In [2]: x = np.random.rand(4, 4)                                                                                   

In [3]: %timeit py_moments_hu(x)                                                                    
7.3 µs ± 37.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [4]: from skimage.measure import moments_hu                                                           

In [5]: %timeit moments_hu(x)                                                                            
4.33 µs ± 20.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Pure python version is as expected slower then cython implementation (1.69 times), but we are talking about µseconds here...

[rfezzani]

This PR replicate #4307 (I forgot about it 😅).

[grlee77]

Thanks @rfezzani, is there any benefit to this PR over the prior one? (i.e. which would you like to have reviewed?).

I am also interested in improving preservation of float32 in general. In working on cupyimg I found quite a few places on the python side where scikit-image converts to float64 when not strictly necessary. I didn't see a currently open issue on this, so I will open one where we can list functions that should have single precision support and keep track of open PRs needing review.

[rfezzani]

is there any benefit to this PR over the prior one?

As you already noticed, #4307 use the dtype parameter but this one not. Choose the option you prefer, I will close the other 😉

[grlee77]

This looks good to me. I prefer this approach to the one in #4307 (I think it is cleaner without adding a new dtype argument to the API)

[grlee77]

This seems fine, but I wonder if we should explicitly raise on error on complex-valued inputs? I think there are probably a ton of other places where this is just assumed and not explicitly checked, thgouh. not sure it is worth it

if nu.dtype.kind = 'c':
    raise TValueError("complex-valued inputs are not suppoted.")

[rfezzani]

I have no strict opinion on this, but most of skimage functions do not explicitly check for complex input values, I am not sure it worth it 😕

[jni]

I'm ok with the PR as is. Thanks @rfezzani!

Yeah maybe we should use pure Python but this is working and is faster so 🤷

[rfezzani]

Thank you @grlee77 🎉