scikit-image icon indicating copy to clipboard operation
scikit-image copied to clipboard

phase_cross_correlation is no more correct in 0.19.3 (instead of 0.18.1)

Open patquem opened this issue 2 years ago • 6 comments

Description

Shift returned by phase_cross_correlative is no more correct when calculating it with 0.19.3 (instead of 0.18.1) on 'smoothed' image

Way to reproduce

import skimage
from skimage.registration import phase_cross_correlation
from skimage.data import binary_blobs
from scipy.ndimage import gaussian_filter


blobs = binary_blobs(256, blob_size_fraction=0.05, volume_fraction=0.7)
ref0 = skimage.img_as_float(blobs)
ref1 = gaussian_filter(ref0.copy(), 3)

shift = [5, 3]

def shift_calculation(ref):
    mov = ref.copy()
    mov = mov[10 - shift[0]:-10 - shift[0], 10 - shift[1]: -10 - shift[1]]
    ref = ref[10:-10, 10: -10]
    return phase_cross_correlation(ref, mov)[0]

shift_calc0 = shift_calculation(ref0)
shift_calc1 = shift_calculation(ref1)

print(f"Shift (binary image) : {shift_calc0}")
print(f"Shift (smoothed image) : {shift_calc1}")

Version information

3.7.7 (tags/v3.7.7:d7c567b08f, Mar 10 2020, 10:41:24) [MSC v.1900 64 bit (AMD64)] Windows-10-10.0.19041-SP0 scikit-image version: 0.19.3 numpy version: 1.21.6

Shift (binary image) : [-5. -3.]
Shift (smoothed image) : [0. 0.]

Version information

3.7.7 (tags/v3.7.7:d7c567b08f, Mar 10 2020, 10:41:24) [MSC v.1900 64 bit (AMD64)] Windows-10-10.0.14393-SP0 scikit-image version: 0.18.1 numpy version: 1.19.5

Shift (binary image) : [-5. -3.]
Shift (smoothed image) : [-5. -3.]

patquem avatar Aug 05 '22 13:08 patquem

I confirm, I can reproduce:

dev branch

3.8.2 (default, Mar 26 2020, 15:53:00) 
[GCC 7.3.0]
Linux-5.13.0-7614-generic-x86_64-with-glibc2.10
scikit-image version: 0.20.0.dev0+git20220725.38b595d60
numpy version: 1.19.2

Shift (binary image) : [-5. -3.]
Shift (smoothed image) : [-1.  0.]

Version 0.18

3.10.4 (main, Mar 31 2022, 08:41:55) [GCC 7.5.0]
Linux-5.13.0-7614-generic-x86_64-with-glibc2.32
scikit-image version: 0.18.3
numpy version: 1.23.1

Shift (binary image) : [-5. -3.]
Shift (smoothed image) : [-5. -3.]

Thank you @patquem for the report :wink:

rfezzani avatar Aug 05 '22 13:08 rfezzani

Well after some investigation, setting normalization parameter to "None" fixes the issue. Reading the note in the doc

In a high noise scenario, the unnormalized method may be preferable.

rfezzani avatar Aug 05 '22 14:08 rfezzani

Hello @rfezzani, Thanks for your anwser. My test case has no noise, so the recommandation in the documentation is not really suitable The problem occured by normalization has already been raised in #6376. The @vincefn 's proposal to go back to the unormalized method by default is maybe the good solution ? Patrick

patquem avatar Aug 05 '22 14:08 patquem

I agree that the normalization by default broke the backward compatibility :grimacing:

rfezzani avatar Aug 05 '22 14:08 rfezzani

Please see #5459 and #5460 for cross ref

rfezzani avatar Aug 05 '22 14:08 rfezzani

My image registration code broke and after some hours and heartache I tracked it down to #5459.

Registering two 2D gaussians now produces the wrong answer:

image

(This is actually the second time I've had to update my code for this one little function, the other one being https://github.com/scikit-image/scikit-image/pull/4502. Maybe there should be a register_translation function that isn't a phase cross-correlation? Or lacks the normalization that seems to be causing regressions?)

joseph-long avatar Aug 18 '22 21:08 joseph-long

This line in phase_cross_correlation does not seem right: https://github.com/scikit-image/scikit-image/blob/74e0652fbbd2bf5600aacc3937e2b566dddb7374/skimage/registration/_phase_cross_correlation.py#L243

This takes the maximum of all |image_product| elements, which is a single scalar, and then the maximum of this scalar with 100*eps, resulting in a single scalar. But I think normalization should divide each product element by their corresponding magnitude, element-wise. $R_{jk}=(G_{a,jk}\circ G_{b,jk}^{*}) / |G_{a,jk}\circ G_{b,jk}|$

A possible code fix:

        product_amp = np.abs(src_freq * target_freq)
        product_amp[product_amp < 100*eps] = 100*eps
        image_product /= product_amp

dirkpitt2050 avatar Nov 30 '22 13:11 dirkpitt2050