shape_index has numerical errors with nearly equal eigenvalues

[scottstanie]

Description

When shape_index comes across an image that has a section where the hessian has identical eigenvalues, the output can flip from -1 to 1 at random pixels.

This produces a warning, but no traceback /home/scott/envs/data/bin/ipython:1: RuntimeWarning: divide by zero encountered in true_divide

Way to reproduce

• [ ] Code example

This can be seen when looking at the shape index of a quadratic bowl:

import numpy as np
import matplotlib.pyplot as plt
from skimage.feature import shape_index

def make_bowl(shape):
    rows, cols = shape
    x = np.linspace(-1, 1, cols)
    y = np.linspace(-1, 1, rows)
    xx, yy = np.meshgrid(x, y)
    z = xx**2 + yy**2
    return z / np.max(z)

bowl = make_bowl((50, 50))
plt.figure(); plt.imshow(shape_index(bowl)); plt.colorbar()  

The answer should be -1 in most of the image except edges, which corresponds to a perfect bowl upward, but there are speckles of +1 (a bowl downward)

• [ ] Relevant images (if any)

• [ ] Operating system and version: Ubuntu 18.04
• [ ] Python version: python 3.6
• [ ] scikit-image version (run skimage.__version__) skimage.__version__ # '0.14.1'

[scottstanie]

Digging in, this happens in feature.shape_index on the line (2.0 / np.pi) * np.arctan((l2 + l1) / (l2 - l1))

This will happen when l2 is exactly l1, as it flips the expression inside arctan from a large negative to np.inf, which is treated as positive (hence the speckles of +1 within the real answer of -1 in the image above).

A simple fix could be add an eps=1e-16 argument to the function and in the denominator do

l2_safe = l2 + eps
return (2.0 / np.pi) * np.arctan((l2_safe + l1) / (l2_safe - l1))

testing this works for me. Let me know if that solution makes sense or if something else would work better overall

[stefanv]

Thanks for the well-formatted issue, @scottstanie. I've checked, and it looks like using arctan2 instead of arctan fixes the instability. Would you like to submit a pull request?

[scottstanie]

Oh awesome, sounds even better. Sure i'll submit one using arctan2

[scottstanie]

Ah, it seems like arctan2 may not be what we want- it returns a value between -pi and pi, but the shape index formula from the paper is 2 / pi * arctan (...), since it is explicitly saying "k1 >= k2" (where those are the eigenvalues l1 and l2.

When I change the code to 2.0 / np.pi) * np.arctan2(l2 + l1, l2 - l1) on my bowl example, I actually get all values of > 1 or <-1, since having a negative denominator means that pi is either added to or subtracted from the the normal arctan result https://en.wikipedia.org/wiki/Atan2#Definition_and_computation

Also, the center of the box becomes all +1 instead of -1, which is what it should be.

[stefanv]

Did you flip the arguments to arctan2? That one often catches me, it is denomenator, then numerator:

https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.arctan2.html

[scottstanie]

Haha you are right I had them switched, that is kinda unexpected.

changing it to

return (2.0 / np.pi) * np.arctan2(l2 - l1, l2 + l1)

I get values ranging from -2 to 0...

Now I can't tell if it's me doing something wrong or if arctan2 doesn't work with the paper's way of implementing this. I'll check again tomorrow