Calcul de la corrélation 2D enveloppée avec fftconvolve

Question

J'ai un ensemble de tableaux 2D que j'ai pour calculer la corrélation 2D de. J'ai essayé beaucoup de choses différentes (même en le programmant dans Fortran), mais je pense que le moyen le plus rapide sera de le calculer en utilisant FFT.

Basé sur mes tests et this answer je peux utiliser scipy.signal.fftconvolve et il fonctionne très bien si je suis en train de reproduire la sortie de scipy.signal.correlate2d avec boundary='fill'. Donc, fondamentalement, ce

scipy.signal.fftconvolve(a, a[::-1, ::-1], mode='same') 

est égal à ce (à l'exception d'un léger décalage)

scipy.signal.correlate2d(a, a, boundary='fill', mode='same') 

La chose est que les tableaux doivent être calculés en mode enveloppé, car ils sont des tableaux périodiques 2D (c'est-à-dire boundary='wrap'). Donc, si je suis en train de reproduire la sortie de

scipy.signal.correlate2d(a, a, boundary='wrap', mode='same') 

Je ne peux pas, ou du moins je ne vois pas comment le faire. (Et je veux utiliser la méthode FFT, car il est beaucoup plus rapide.)

Apparemment Scipy avait l'habitude d'avoir something like that qui aurait pu faire l'affaire, mais apparemment il a été laissé de côté et je ne peux pas le trouver, donc je pense Scipy pourrait avoir abandonné le support pour cela.

Quoi qu'il en soit, existe-t-il un moyen d'utiliser les routines FFT scipy ou numpy pour calculer cette corrélation des tableaux de périodes?

Answer 1

La corrélation enveloppée peut être implémentée en utilisant la FFT. Voici un code pour montrer comment:

In [276]: import numpy as np 

In [277]: from scipy.signal import correlate2d 

Créez un tableau aléatoire a à travailler avec:

In [278]: a = np.random.randn(200, 200) 

Calculer la corrélation 2D à l'aide scipy.signal.correlate2d:

In [279]: c = correlate2d(a, a, boundary='wrap', mode='same') 

calcule maintenant le même résultat, en utilisant les fonctions 2D FFT de numpy.fft. (Ce code suppose a est carré.)

In [280]: from numpy.fft import fft2, ifft2 

In [281]: fc = np.roll(ifft2(fft2(a).conj()*fft2(a)).real, (a.shape[0] - 1)//2, axis=(0,1)) 

Vérifiez que les deux méthodes donnent le même résultat:

In [282]: np.allclose(c, fc) 
Out[282]: True 

Et comme vous le soulignez, en utilisant la FFT est beaucoup plus rapide. Pour cet exemple, il est d'environ 1000 fois plus vite:

In [283]: %timeit c = correlate2d(a, a, boundary='wrap', mode='same') 
1 loop, best of 3: 3.2 s per loop 

In [284]: %timeit fc = np.roll(ifft2(fft2(a).conj()*fft2(a)).real, (a.shape[0] - 1)//2, axis=(0,1)) 
100 loops, best of 3: 3.19 ms per loop 

Et qui comprend le calcul de dupliquée fft2(a). Bien sûr, fft2(a) ne doit être calculé qu'une fois:

In [285]: fta = fft2(a) 

In [286]: fc = np.roll(ifft2(fta.conj()*fta).real, (a.shape[0] - 1)//2, axis=(0,1))