Merge branch 'topic/default/piv-various' into 'branch/default'

Various PIV improvements: perf + better correl and fft

See merge request fluiddyn/fluidimage!95

doc/examples/piv_try_params.py

@@ -6,16 +6,16 @@
 
 params.series.path = "../../image_samples/Oseen/Images"
 
-params.mask.strcrop = "30:250, 100:"
+params.mask.strcrop = "30:250, 40:"
 
 params.piv0.shape_crop_im0 = 32
-params.piv0.displacement_max = 5
+params.piv0.displacement_max = 7
 
 params.fix.correl_min = 0.2
 params.fix.threshold_diff_neighbour = 8
 
 params.multipass.number = 2
-params.multipass.use_tps = True
+params.multipass.use_tps = "last"
 
 work = Work(params=params)
 

doc/examples/profile_piv_work.py

@@ -27,15 +27,15 @@
 params.piv0.shape_crop_im0 = 40
 params.piv0.displacement_max = 14
 
-params.piv0.nb_peaks_to_search = 2
+params.piv0.nb_peaks_to_search = 1
 params.piv0.particle_radius = 3
 
 params.mask.strcrop = ":, :1500"
 
 params.multipass.number = 2
 
-params.multipass.use_tps = "last"
-# params.multipass.use_tps = False
+# params.multipass.use_tps = "last"
+params.multipass.use_tps = False
 params.multipass.subdom_size = 200
 params.multipass.smoothing_coef = 10.0
 params.multipass.threshold_tps = 0.5
@@ -56,8 +56,10 @@
             unicode=True,
             color=True,
             show_all=False,
-            # time="percent_of_total",
+            time="percent_of_total",
             # flat=True,  # buggy with pyinstrument 4.6.2!
         )
     )
 )
+
+# piv.display(show_interp=False, scale=1, show_error=True)

src/fluidimage/calcul/__init__.py

@@ -7,10 +7,10 @@
 .. autosummary::
    :toctree:
 
-   fft
    correl
-   subpix
+   fft
    interpolate
-   smooth_clean
+   mean_neighbors
+   subpix
 
 """

src/fluidimage/calcul/correl.py

@@ -43,16 +43,21 @@
 """
 
 from abc import ABC, abstractmethod
+from math import sqrt
 
 import numpy as np
-from numpy.fft import fft2, ifft2
 from scipy.ndimage import correlate
 from scipy.signal import correlate2d
 from transonic import Array, Type, boost
 
 from .correl_pycuda import correl_pycuda
 from .errors import PIVError
-from .fft import CUFFT2DReal2Complex, FFTW2DReal2Complex, SKCUFFT2DReal2Complex
+from .fft import (
+    CUFFT2DReal2Complex,
+    FFTW2DReal2Complex,
+    NumpyFFT2DReal2Complex,
+    SKCUFFT2DReal2Complex,
+)
 from .subpix import SubPix
 
 
@@ -516,7 +521,7 @@ def _finalize_init(self):
 
             for indices, v in np.ndenumerate(where_large_displacement):
                 dx, dy = self.compute_displacement_from_indices(*indices[::-1])
-                displacement = np.sqrt(dx**2 + dy**2)
+                displacement = sqrt(dx**2 + dy**2)
                 if displacement > self.displacement_max:
                     where_large_displacement[indices] = True
 
@@ -569,7 +574,7 @@ def _norm_images(im0: A2df32, im1: A2df32):
 
     Should return something close to:
 
-    np.sqrt(np.sum(im1**2) * np.sum(im0**2))
+    sqrt(np.sum(im1**2) * np.sum(im0**2))
 
     """
     im0 = im0.ravel()
@@ -585,7 +590,7 @@ def _norm_images(im0: A2df32, im1: A2df32):
         for idx in range(1, im0.size):
             tmp0 += im0[idx] ** 2
             tmp1 += im1[idx] ** 2
-    return np.sqrt(tmp0 * tmp1)
+    return sqrt(tmp0 * tmp1)
 
 
 @boost
@@ -631,18 +636,6 @@ def _like_fftshift(arr: A2dC):
     return arr
 
 
-class CorrelFFTNumpy(CorrelFFTBase):
-    """Correlations using numpy.fft."""
-
-    _tag = "np.fft"
-
-    def __call__(self, im0, im1):
-        """Compute the correlation from images."""
-        norm = _norm_images(im0, im1)
-        correl = ifft2(fft2(im0).conj() * fft2(im1)).real
-        return _like_fftshift(np.ascontiguousarray(correl)), norm
-
-
 class CorrelFFTWithOperBase(CorrelFFTBase):
 
     FFTClass: object
@@ -658,12 +651,40 @@ def __call__(self, im0, im1):
         Warning: important for perf, so use Pythran
 
         """
-        norm = _norm_images(im0, im1) * im0.size
-        oper = self.oper
-        correl = oper.ifft(oper.fft(im0).conj() * oper.fft(im1))
+        fft_im0 = self.oper.fft(im0)
+        fft_im0[0, 0] = 0.0
+        fft_im1 = self.oper.fft(im1)
+        fft_im1[0, 0] = 0.0
+        energy0 = self.oper.compute_energy_from_fourier(fft_im0)
+        energy1 = self.oper.compute_energy_from_fourier(fft_im1)
+        norm = sqrt(4 * energy0 * energy1) * self.oper.coef_norm_correl
+        correl = self.oper.ifft(fft_im0.conj() * fft_im1)
         return _like_fftshift(correl), norm
 
 
+class CorrelFFTNumpy(CorrelFFTWithOperBase):
+    """Correlations using numpy.fft."""
+
+    FFTClass = NumpyFFT2DReal2Complex
+    _tag = "np.fft"
+
+    def __call__(self, im0, im1):
+        """Compute the correlation from images.
+
+        Warning: important for perf, so use Pythran
+
+        """
+        fft_im0 = self.oper.fft(im0)
+        fft_im0[0, 0] = 0.0
+        fft_im1 = self.oper.fft(im1)
+        fft_im1[0, 0] = 0.0
+        correl = self.oper.ifft(fft_im0.conj() * fft_im1).real
+        energy0 = self.oper.compute_energy_from_fourier(fft_im0)
+        energy1 = self.oper.compute_energy_from_fourier(fft_im1)
+        norm = sqrt(4 * energy0 * energy1) * self.oper.coef_norm_correl
+        return _like_fftshift(np.ascontiguousarray(correl)), norm
+
+
 class CorrelFFTW(CorrelFFTWithOperBase):
     """Correlations using fluidimage.fft.FFTW2DReal2Complex"""