fixed DoG function for correct scaling

[cmarasinou]

Description

This is an fix of the function for generating the difference of Gaussians in blob.py

The correct approximation should divide (the existing implementation) by the change in the standard deviation between the two Gaussians. Simplifying this we can remove the multiplication by the standard deviation and divide by (sigma_ratio-1).

This is an important fix since it affects differently each scale and therefore it affects the local maxima in the scale-space representation. This can be seen in images with blobs at different scales.

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.
• Consider backporting the PR with @meeseeksdev backport to v0.14.x

[pep8speaks]

Hello @cmarasinou! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file skimage/feature/tests/test_blob.py:

Line 53:25: E225 missing whitespace around operator

Comment last updated at 2021-08-05 13:58:41 UTC

[cmarasinou]

Here are two examples:

Example 1

With blob_dog so far:

With fix:

Example 2

With blob_dog so far:

With fix:

[ericPrince]

Hi, just wanted to point out that LoG blob detection uses the same sigma scaling as DoG. If this change is approved for DoG, should LoG also be updated to be consistent? Or is that scaling correct for LoG?

https://github.com/scikit-image/scikit-image/blob/main/skimage/feature/blob.py#L501

    gl_images = [-gaussian_laplace(image, s) * np.mean(s) ** 2
                 for s in sigma_list]

[grlee77]

Thank you @cmarasinou . Apologies that this has gone neglected so long, but I think this does indeed look correct.

Here is a simple 1d example involving synthetic sinusoidal signals at two different frequencies. The peak amplitude across scales is only comparable with scaling as adjusted here:

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter, gaussian_laplace

n = 256
s0 = np.zeros(n)
s2 = np.cos(np.linspace(0, 4 * np.pi, 2 * n))
s4 = np.cos(np.linspace(0, 16 * np.pi, 2 * n))
sig = np.concatenate((s0, s0, s0, s0, s2, s2, s0, s0, s0, s4, s0, s0, s0))

impulse_sig = False
if impulse_sig:
    sig = np.zeros_like(sig)
    sig[sig.size // 2] = 1

fig, axes = plt.subplots(3, 1, figsize=[12, 14])
axes[0].plot(sig, 'k')
axes[0].set_title('signal')

min_sigma = 1
max_sigma = 250
sigma_ratio = 1.6
# k such that min_sigma*(sigma_ratio**k) > max_sigma
k = int(np.mean(np.log(max_sigma / min_sigma) / np.log(sigma_ratio) + 1))

# a geometric progression of standard deviations for gaussian kernels
sigma_list = np.array([min_sigma * (sigma_ratio ** i)
                       for i in range(k + 1)])


gaussian_signals = [gaussian_filter(sig, s) for s in sigma_list]


# compute difference of Gaussians using the current scaling
dog_signals_current_scaling = [(gaussian_signals[i] - gaussian_signals[i + 1])
                                * np.mean(sigma_list[i]) for i in range(k)]

dog_current = np.stack(dog_signals_current_scaling, axis=-1)
axes[1].plot(dog_current)
axes[1].set_title('DoG (current scaling)')


# compute difference of Gaussians using the proposed scaling
dog_signals_proposed = [(gaussian_signals[i] - gaussian_signals[i + 1])
                        for i in range(k)]
dog_proposed = np.stack(dog_signals_proposed, axis=-1)
axes[2].plot(dog_proposed)
axes[2].set_title('DoG (proposed scaling)')

plt.tight_layout()
plt.show()

Here the different colors in the two lower plots correspond to differences of Gaussians at different scales. In the lower plot where the scaling is as proposed in this PR, the peak amplitude across scales looks the same for both sinusoids.

[grlee77]

If this change is approved for DoG, should LoG also be updated to be consistent? Or is that scaling correct for LoG?

No, the LoG case seems to be correct as-is. Using the same signal from the example above, but with:

gl_images = [-gaussian_laplace(sig, s) * np.mean(s) ** 2
             for s in sigma_list]
gl_stack = np.stack(gl_images, axis=-1)
plt.figure()
plt.plot(gl_stack)

Gives a figure with the same peak amplitude across scales as expected.

[grlee77]

I fixed the reported PEP8 issue and updated tolerances in one test that was failing locally to see if we can pass the automated tests now

[mkcor]

Thank you so much, @cmarasinou! (Plus I'm not sure what it meant to take the mean of a single element: np.mean(sigma_list[i])?)

But where is the division by (sigma_ratio-1), please?
I'm referring to your PR description:

Simplifying this we can remove the multiplication by the standard deviation and divide by (sigma_ratio-1).

and to the bottom of page 6 in the paper you linked to (https://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf).

[grlee77]

I think whether or not we divide by (sigma_ratio - 1) or any global scale factor will not make any difference. This is because the following lines just use peak_local_max to return the coordinates of the peaks (the absolute magnitude of the peaks are not ever used).

[mkcor]

Oh, right, right! Thanks, @grlee77 😄
Waiting for CI to complete and merging

[grlee77]

Hi @mkcor, I will need to make a followup to this revisiting the absolute scaling. I had not noticed that there is also an absolute threshold being passed onto peak_local_max as well, so it will be important in terms of the appropriate setting of the user-parameter threshold. Also, the docstring example now longer matches after this PR, but this was not caught due to #5502!

For reference this is what image_cube looks like for the docstring example prior to this PR:

vs. after this PR where you can see that there is less bias toward larger scales

To avoid small scale edges being detected we will probably have to increase min_sigma substantially for the docstring example and fix the threshold.

[mkcor]

Oops! I'm diving deeper and it's a little scary, indeed... Reviewing your PR #5503 now!