diff --git a/doc/examples/piv_try_params.py b/doc/examples/piv_try_params.py
index 1fa40400..64fdf961 100644
--- a/doc/examples/piv_try_params.py
+++ b/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)
 
diff --git a/doc/examples/profile_piv_work.py b/doc/examples/profile_piv_work.py
index 78b606d1..1c8ebea5 100644
--- a/doc/examples/profile_piv_work.py
+++ b/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)
diff --git a/src/fluidimage/calcul/__init__.py b/src/fluidimage/calcul/__init__.py
index 22327d8b..3221680a 100644
--- a/src/fluidimage/calcul/__init__.py
+++ b/src/fluidimage/calcul/__init__.py
@@ -7,10 +7,10 @@
 .. autosummary::
    :toctree:
 
-   fft
    correl
-   subpix
+   fft
    interpolate
-   smooth_clean
+   mean_neighbors
+   subpix
 
 """
diff --git a/src/fluidimage/calcul/correl.py b/src/fluidimage/calcul/correl.py
index 762683e1..a006ba0f 100644
--- a/src/fluidimage/calcul/correl.py
+++ b/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"""
 
diff --git a/src/fluidimage/calcul/fft.py b/src/fluidimage/calcul/fft.py
index 7a6d4724..6261676e 100644
--- a/src/fluidimage/calcul/fft.py
+++ b/src/fluidimage/calcul/fft.py
@@ -26,7 +26,11 @@
 
 """
 
+from abc import ABC, abstractmethod
+
 import numpy as np
+import numpy.fft as np_fft
+from transonic import Array, Type, boost
 
 try:
     import pyfftw
@@ -62,7 +66,70 @@
 nthreads = 1
 
 
-class CUFFT2DReal2Complex:
+A2d_complex = Array[Type(np.complex64, np.complex128), "2d"]
+
+
+@boost
+def _compute_energy_from_fourier(field_fft: A2d_complex, coef_norm: int):
+    """Simple Pythran implementation of
+
+    (
+        0.5 / coef_norm
+        * (
+            np.sum(abs(field_fft[:, 0]) ** 2 + abs(field_fft[:, -1]) ** 2)
+            + 2 * np.sum(abs(field_fft[:, 1:-1]) ** 2)
+        )
+    )
+    """
+    n0, n1 = field_fft.shape
+    result = 0.0
+    for i0 in range(n0):
+        result += abs(field_fft[i0, 0]) ** 2 + abs(field_fft[i0, n1 - 1]) ** 2
+        for i1 in range(1, n1 - 1):
+            result += 2 * abs(field_fft[i0, i1]) ** 2
+    return 0.5 / coef_norm * result
+
+
+class OperatorFFTBase(ABC):
+    """Abstract class for FFT operators"""
+
+    coef_norm: int
+    type_real: str = "float32"
+    type_complex: str = "complex64"
+
+    def __init__(self, nx, ny):
+        shapeX = [ny, nx]
+        shapeK = [ny, nx // 2 + 1]
+
+        self.shapeX = shapeX
+        self.shapeK = shapeK
+
+        self.coef_norm = nx * ny
+        self.coef_norm_correl = self.coef_norm**2
+        self.coef_norm_energy = self.coef_norm**2
+
+    @abstractmethod
+    def fft(self, field):
+        """Forwards Fast Fourier Transform"""
+
+    @abstractmethod
+    def ifft(self, field_fft):
+        """Inverse Fast Fourier Transform"""
+
+    def compute_energy_from_fourier(self, field_fft):
+        """Compute the energy from a field in Fourier space"""
+        return _compute_energy_from_fourier(field_fft, self.coef_norm_energy)
+
+    def compute_energy_from_spatial(self, field):
+        """Compute energy from a field in real space"""
+        return np.mean(abs(field) ** 2) / 2
+
+    def project_fft_on_real(self, field_fft):
+        """Project a field in Fourier space onto the real manifold"""
+        return self.fft(self.ifft(field_fft))
+
+
+class CUFFT2DReal2Complex(OperatorFFTBase):
     """A class to use cufft with float32."""
 
     type_real = "float32"
@@ -90,24 +157,18 @@ def __init__(self, nx, ny):
 
         self.coef_norm = nx * ny
 
-    def fft(self, ff):
-        arr_dev = self.thr.to_device(ff.astype(self.type_complex))
+    def fft(self, field):
+        arr_dev = self.thr.to_device(field.astype(self.type_complex))
         self.fftplan(arr_dev, arr_dev, 1.0 / self.coef_norm)
         return arr_dev.get()
 
-    def ifft(self, ff_fft):
-        arr_dev = self.thr.to_device(ff_fft)
+    def ifft(self, field_fft):
+        arr_dev = self.thr.to_device(field_fft)
         self.fftplan(arr_dev, arr_dev, self.coef_norm, inverse=True)
         return arr_dev.get()
 
-    def compute_energy_from_Fourier(self, ff_fft):
-        return np.sum(abs(ff_fft) ** 2) / 2
-
-    def compute_energy_from_spatial(self, ff):
-        return np.mean(abs(ff) ** 2) / 2
-
-    def project_fft_on_realX(self, ff_fft):
-        return self.fft(self.ifft(ff_fft))
+    def compute_energy_from_fourier(self, field_fft):
+        return np.sum(abs(field_fft) ** 2) / 2 / self.coef_norm_energy
 
 
 class CUFFT2DReal2ComplexFloat64(CUFFT2DReal2Complex):
@@ -117,50 +178,31 @@ class CUFFT2DReal2ComplexFloat64(CUFFT2DReal2Complex):
     type_complex = "complex128"
 
 
-class SKCUFFT2DReal2Complex:
+class SKCUFFT2DReal2Complex(OperatorFFTBase):
     """A class to use skcuda-cufft with float32."""
 
     type_real = "float32"
     type_complex = "complex64"
 
     def __init__(self, nx, ny):
-        shapeX = [ny, nx]
-        shapeK = [ny, nx // 2 + 1]
-
-        self.shapeX = shapeX
-        self.shapeK = shapeK
-
+        super().__init__(nx, ny)
         self.fftplan = skfft.Plan(self.shapeX, np.float32, np.complex64)
         self.ifftplan = skfft.Plan(self.shapeX, np.complex64, np.float32)
 
-        self.coef_norm = nx * ny
-
-    def fft(self, ff):
-        x_gpu = gpuarray.to_gpu(ff)
+    def fft(self, field):
+        x_gpu = gpuarray.to_gpu(field)
         xf_gpu = gpuarray.empty(self.shapeK, np.complex64)
         skfft.fft(x_gpu, xf_gpu, self.fftplan, False)
         return xf_gpu.get()
 
-    def ifft(self, ff_fft):
-        xf_gpu = gpuarray.to_gpu(ff_fft)
+    def ifft(self, field_fft):
+        xf_gpu = gpuarray.to_gpu(field_fft)
         x_gpu = gpuarray.empty(self.shapeX, np.float32)
         skfft.ifft(xf_gpu, x_gpu, self.ifftplan, False)
         return x_gpu.get()
 
-    def compute_energy_from_Fourier(self, ff_fft):
-        return (
-            np.sum(abs(ff_fft[:, 0]) ** 2 + abs(ff_fft[:, -1]) ** 2)
-            + 2 * np.sum(abs(ff_fft[:, 1:-1]) ** 2)
-        ) / 2
-
-    def compute_energy_from_spatial(self, ff):
-        return np.mean(abs(ff) ** 2) / 2
-
-    def project_fft_on_realX(self, ff_fft):
-        return self.fft(self.ifft(ff_fft))
-
 
-class FFTW2DReal2Complex:
+class FFTW2DReal2Complex(OperatorFFTBase):
     """A class to use fftw with float32.
 
     These ffts are NOT normalized (faster)!
@@ -171,14 +213,10 @@ class FFTW2DReal2Complex:
     type_complex = "complex64"
 
     def __init__(self, nx, ny):
-        shapeX = [ny, nx]
-        shapeK = [ny, nx // 2 + 1]
+        super().__init__(nx, ny)
 
-        self.shapeX = shapeX
-        self.shapeK = shapeK
-
-        self.arrayX = pyfftw.empty_aligned(shapeX, self.type_real)
-        self.arrayK = pyfftw.empty_aligned(shapeK, self.type_complex)
+        self.arrayX = pyfftw.empty_aligned(self.shapeX, self.type_real)
+        self.arrayK = pyfftw.empty_aligned(self.shapeK, self.type_complex)
 
         self.fftplan = pyfftw.FFTW(
             input_array=self.arrayX,
@@ -195,33 +233,52 @@ def __init__(self, nx, ny):
             threads=nthreads,
         )
 
-        self.coef_norm = nx * ny
-
-    def fft(self, ff):
-        self.arrayX[:] = ff
+    def fft(self, field):
+        self.arrayX[:] = field
         self.fftplan()
         return self.arrayK.copy()
 
-    def ifft(self, ff_fft):
-        self.arrayK[:] = ff_fft
+    def ifft(self, field_fft):
+        self.arrayK[:] = field_fft
         self.ifftplan(normalise_idft=False)
         return self.arrayX.copy()
 
-    def compute_energy_from_Fourier(self, ff_fft):
-        return (
-            np.sum(abs(ff_fft[:, 0]) ** 2 + abs(ff_fft[:, -1]) ** 2)
-            + 2 * np.sum(abs(ff_fft[:, 1:-1]) ** 2)
-        ) / 2
 
-    def compute_energy_from_spatial(self, ff):
-        return np.mean(abs(ff) ** 2) / 2
+class FFTW2DReal2ComplexFloat64(FFTW2DReal2Complex):
+    """A class to use fftw with float64."""
+
+    type_real = "float64"
+    type_complex = "complex128"
 
-    def project_fft_on_realX(self, ff_fft):
-        return self.fft(self.ifft(ff_fft))
 
+class NumpyFFT2DReal2Complex(OperatorFFTBase):
+    """FFT operator using numpy.fft"""
 
-class FFTW2DReal2ComplexFloat64(FFTW2DReal2Complex):
-    """A class to use fftw with float64."""
+    def __init__(self, nx, ny):
+        super().__init__(nx, ny)
+        self.coef_norm_correl = self.coef_norm
+        self.coef_norm = 1
+
+    def fft(self, field):
+        return np_fft.fft2(field)
+
+    def ifft(self, field_fft):
+        return np_fft.ifft2(field_fft)
+
+    def compute_energy_from_fourier(self, field_fft):
+        return np.sum(abs(field_fft) ** 2) / 2 / self.coef_norm_energy
+
+
+class NumpyFFT2DReal2ComplexFloat64(NumpyFFT2DReal2Complex):
+    """FFT operator using numpy.fft"""
 
     type_real = "float64"
     type_complex = "complex128"
+
+
+classes = [
+    FFTW2DReal2Complex,
+    FFTW2DReal2ComplexFloat64,
+    NumpyFFT2DReal2Complex,
+    NumpyFFT2DReal2ComplexFloat64,
+]
diff --git a/src/fluidimage/calcul/interpolate/griddata.py b/src/fluidimage/calcul/interpolate/griddata.py
index aa587518..169031e0 100644
--- a/src/fluidimage/calcul/interpolate/griddata.py
+++ b/src/fluidimage/calcul/interpolate/griddata.py
@@ -5,27 +5,18 @@
 """
 
 import numpy as np
-from matplotlib import mlab
 from scipy import interpolate
 
 
-def griddata(centers, values, new_positions, using="scipy"):
+def griddata(centers, values, new_positions):
     xnew = new_positions[0]
     ynew = new_positions[1]
     grid_x, grid_y = np.meshgrid(xnew, ynew)
 
-    if using == "scipy":
-        # 'linear' or 'cubic'
-        values_new = interpolate.griddata(
-            centers.T, values, (grid_x, grid_y), "cubic"
-        )
-    elif using == "matplotlib":
-        x = centers[0]
-        y = centers[1]
-        values_new = mlab.griddata(x, y, values, xnew, ynew, "linear")
-        values_new[values_new.mask] = np.nan
-    else:
-        raise ValueError
+    # 'linear' or 'cubic'
+    values_new = interpolate.griddata(
+        centers.T, values, (grid_x, grid_y), "linear"
+    )
 
     inds = np.where(np.isnan(values_new))
     values_nearest = interpolate.griddata(
diff --git a/src/fluidimage/calcul/mean_neighbors.py b/src/fluidimage/calcul/mean_neighbors.py
new file mode 100644
index 00000000..bc1bc066
--- /dev/null
+++ b/src/fluidimage/calcul/mean_neighbors.py
@@ -0,0 +1,155 @@
+import numpy as np
+from transonic import boost
+
+
+def _mean_neighbors(arr, arr_is_not_nan):
+    """Nan has to be zeros in `arr`!"""
+
+    arr_a = arr
+    arr_b = arr_is_not_nan.astype(np.uint8)
+
+    n0, n1 = arr.shape
+    arr_output = np.zeros_like(arr)
+    arr_nb_pts = np.zeros(arr.shape, dtype=np.uint8)
+
+    # corners
+    i0, i1 = 0, 0
+    if arr_is_not_nan[i0, i1]:
+        arr_output[i0, i1] = (
+            arr_a[i0 + 1, i1] + arr_a[i0 + 1, i1 + 1] + arr_a[i0, i1 + 1]
+        )
+        arr_nb_pts[i0, i1] = (
+            arr_b[i0 + 1, i1] + arr_b[i0 + 1, i1 + 1] + arr_b[i0, i1 + 1]
+        )
+
+    i0, i1 = 0, n1 - 1
+    if arr_is_not_nan[i0, i1]:
+        arr_output[i0, i1] = (
+            arr_a[i0 + 1, i1] + arr_a[i0 + 1, i1 - 1] + arr_a[i0, i1 - 1]
+        )
+        arr_nb_pts[i0, i1] = (
+            arr_b[i0 + 1, i1] + arr_b[i0 + 1, i1 - 1] + arr_b[i0, i1 - 1]
+        )
+
+    i0, i1 = n0 - 1, 0
+    if arr_is_not_nan[i0, i1]:
+        arr_output[i0, i1] = (
+            arr_a[i0 - 1, i1] + arr_a[i0 - 1, i1 + 1] + arr_a[i0, i1 + 1]
+        )
+        arr_nb_pts[i0, i1] = (
+            arr_b[i0 - 1, i1] + arr_b[i0 - 1, i1 + 1] + arr_b[i0, i1 + 1]
+        )
+
+    i0, i1 = n0 - 1, n1 - 1
+    if arr_is_not_nan[i0, i1]:
+        arr_output[i0, i1] = (
+            arr_a[i0 - 1, i1] + arr_a[i0 - 1, i1 - 1] + arr_a[i0, i1 - 1]
+        )
+        arr_nb_pts[i0, i1] = (
+            arr_b[i0 - 1, i1] + arr_b[i0 - 1, i1 - 1] + arr_b[i0, i1 - 1]
+        )
+
+    # upper raw
+    i0 = 0
+    for i1 in range(1, n1 - 1):
+        if not arr_is_not_nan[i0, i1]:
+            continue
+        tmp_a = arr_a[i0, i1 - 1] + arr_a[i0, i1 + 1]
+        tmp_b = arr_b[i0, i1 - 1] + arr_b[i0, i1 + 1]
+        for idx in range(-1, 2):
+            tmp_a += arr_a[i0 + 1, i1 - idx]
+            tmp_b += arr_b[i0 + 1, i1 - idx]
+        arr_output[i0, i1] = tmp_a
+        arr_nb_pts[i0, i1] = tmp_b
+
+    # bottom raw
+    i0 = n0 - 1
+    for i1 in range(1, n1 - 1):
+        if not arr_is_not_nan[i0, i1]:
+            continue
+        tmp_a = arr_a[i0, i1 - 1] + arr_a[i0, i1 + 1]
+        tmp_b = arr_b[i0, i1 - 1] + arr_b[i0, i1 + 1]
+        for idx in range(-1, 2):
+            tmp_a += arr_a[i0 - 1, i1 - idx]
+            tmp_b += arr_b[i0 - 1, i1 - idx]
+        arr_output[i0, i1] = tmp_a
+        arr_nb_pts[i0, i1] = tmp_b
+
+    # left column
+    i1 = 0
+    for i0 in range(1, n0 - 1):
+        if not arr_is_not_nan[i0, i1]:
+            continue
+        tmp_a = arr_a[i0 - 1, i1] + arr_a[i0 + 1, i1]
+        tmp_b = arr_b[i0 - 1, i1] + arr_b[i0 + 1, i1]
+        for idx in range(-1, 2):
+            tmp_a += arr_a[i0 + idx, i1 + 1]
+            tmp_b += arr_b[i0 + idx, i1 + 1]
+        arr_output[i0, i1] = tmp_a
+        arr_nb_pts[i0, i1] = tmp_b
+
+    # right column
+    i1 = n1 - 1
+    for i0 in range(1, n0 - 1):
+        if not arr_is_not_nan[i0, i1]:
+            continue
+        tmp_a = arr_a[i0 - 1, i1] + arr_a[i0 + 1, i1]
+        tmp_b = arr_b[i0 - 1, i1] + arr_b[i0 + 1, i1]
+        for idx in range(-1, 2):
+            tmp_a += arr_a[i0 + idx, i1 - 1]
+            tmp_b += arr_b[i0 + idx, i1 - 1]
+        arr_output[i0, i1] = tmp_a
+        arr_nb_pts[i0, i1] = tmp_b
+
+    # core
+    for i0 in range(1, n0 - 1):
+        for i1 in range(1, n1 - 1):
+            if not arr_is_not_nan[i0, i1]:
+                continue
+            tmp_a = arr_a[i0, i1 - 1] + arr_a[i0, i1 + 1]
+            tmp_b = arr_b[i0, i1 - 1] + arr_b[i0, i1 + 1]
+            for idx in range(-1, 2):
+                tmp_a += arr_a[i0 - 1, i1 - idx] + arr_a[i0 + 1, i1 - idx]
+                tmp_b += arr_b[i0 - 1, i1 - idx] + arr_b[i0 + 1, i1 - idx]
+            arr_output[i0, i1] = tmp_a
+            arr_nb_pts[i0, i1] = tmp_b
+
+    # normalize
+    for i0 in range(n0):
+        for i1 in range(n1):
+            if not arr_is_not_nan[i0, i1] or arr_nb_pts[i0, i1] == 0:
+                arr_output[i0, i1] = np.nan
+                continue
+            arr_output[i0, i1] /= arr_nb_pts[i0, i1]
+
+    return arr_output
+
+
+@boost
+def mean_neighbors_xy(
+    deltaxs: "float32[]",
+    deltays: "float32[]",
+    iyvecs: "int64[]",
+    ixvecs: "int64[]",
+):
+    """Compute the mean over the nearest neighbors"""
+
+    n1 = nx = len(ixvecs)
+    n0 = ny = len(iyvecs)
+
+    shape = (ny, nx)
+    arr_is_not_nan = np.ones(shape, dtype=bool)
+
+    deltaxs2d = deltaxs.copy().reshape(shape)
+    deltays2d = deltays.copy().reshape(shape)
+
+    for i0 in range(n0):
+        for i1 in range(n1):
+            if np.isnan(deltaxs2d[i0, i1]) or np.isnan(deltays2d[i0, i1]):
+                arr_is_not_nan[i0, i1] = False
+                deltaxs2d[i0, i1] = 0.0
+                deltays2d[i0, i1] = 0.0
+
+    mean_x = _mean_neighbors(deltaxs2d, arr_is_not_nan)
+    mean_y = _mean_neighbors(deltays2d, arr_is_not_nan)
+    return mean_x.reshape(n0 * n1), mean_y.reshape(n0 * n1)
diff --git a/src/fluidimage/calcul/meson.build b/src/fluidimage/calcul/meson.build
index 2b1cf574..cf2f07f4 100644
--- a/src/fluidimage/calcul/meson.build
+++ b/src/fluidimage/calcul/meson.build
@@ -6,7 +6,7 @@ python_sources = [
   'correl_pycuda.py',
   'errors.py',
   'fft.py',
-  'smooth_clean.py',
+  'mean_neighbors.py',
   'subpix.py',
   'test_correl.py',
   'test_fft.py',
@@ -19,7 +19,7 @@ py.install_sources(
 
 subdir('interpolate')
 
-run_command(['transonic', '--meson', '--backend', backend, 'correl.py', 'subpix.py'], check: true)
+run_command(['transonic', '--meson', '--backend', backend, 'correl.py', 'subpix.py', 'mean_neighbors.py', 'fft.py'], check: true)
 
 foreach be : backends
   subdir('__' + be + '__')
diff --git a/src/fluidimage/calcul/smooth_clean.py b/src/fluidimage/calcul/smooth_clean.py
deleted file mode 100644
index 9a1ef0e8..00000000
--- a/src/fluidimage/calcul/smooth_clean.py
+++ /dev/null
@@ -1,75 +0,0 @@
-"""Smooth 2d fields
-===================
-
-
-"""
-
-import itertools
-
-import numpy as np
-from scipy.interpolate import LinearNDInterpolator
-from scipy.ndimage import convolve
-
-from fluidimage.calcul.interpolate.griddata import griddata
-
-weights = np.ones([3, 3])
-
-
-def _smooth(a, for_norm):
-    norm = convolve(for_norm, weights, mode="nearest")
-    ind = np.where(norm == 0)
-    norm[ind] = 1
-    return convolve(a, weights, mode="nearest") / norm
-
-
-def smooth_clean(xs, ys, deltaxs, deltays, iyvecs, ixvecs, threshold):
-    """Smooth and clean the displacements
-
-    Warning: important for perf (~25% for PIV)
-
-    """
-    nx = len(ixvecs)
-    ny = len(iyvecs)
-
-    shape = [ny, nx]
-    for_norm = np.ones(shape)
-
-    selection = ~(np.isnan(deltaxs) | np.isnan(deltays))
-    if not selection.any():
-        return deltaxs, deltays
-
-    centers = np.vstack([xs[selection], ys[selection]])
-    dxs_select = deltaxs[selection]
-    dys_select = deltays[selection]
-    dxs = griddata(centers, dxs_select, (ixvecs, iyvecs)).reshape([ny, nx])
-    dys = griddata(centers, dys_select, (ixvecs, iyvecs)).reshape([ny, nx])
-
-    dxs2 = _smooth(dxs, for_norm)
-    dys2 = _smooth(dys, for_norm)
-
-    inds = (abs(dxs2 - dxs) + abs(dys2 - dys) > threshold).nonzero()
-
-    dxs[inds] = 0
-    dys[inds] = 0
-    for_norm[inds] = 0
-
-    dxs_smooth = _smooth(dxs, for_norm)
-    dys_smooth = _smooth(dys, for_norm)
-
-    # come back to the unstructured grid
-    xy = list(itertools.product(ixvecs, iyvecs))
-    interpolator_x = LinearNDInterpolator(xy, dxs_smooth.T.flat)
-    interpolator_y = LinearNDInterpolator(xy, dys_smooth.T.flat)
-    out_dxs = interpolator_x(xs, ys)
-    out_dys = interpolator_y(xs, ys)
-
-    # previous implementation (interp2d is depreciated)
-    # fxs = interp2d(ixvecs, iyvecs, dxs_smooth, kind="linear")
-    # fys = interp2d(ixvecs, iyvecs, dys_smooth, kind="linear")
-    # out_dxs = np.empty_like(deltaxs)
-    # out_dys = np.empty_like(deltays)
-    # for i, (x, y) in enumerate(zip(xs, ys)):
-    #     out_dxs[i] = fxs(x, y)[0]
-    #     out_dys[i] = fys(x, y)[0]
-
-    return out_dxs, out_dys
diff --git a/src/fluidimage/calcul/test_correl.py b/src/fluidimage/calcul/test_correl.py
index 3d368ada..1823ba40 100644
--- a/src/fluidimage/calcul/test_correl.py
+++ b/src/fluidimage/calcul/test_correl.py
@@ -17,7 +17,7 @@
 logger = logging.getLogger("fluidimage")
 
 classes = {k.replace(".", "_"): v for k, v in correlation_classes.items()}
-classes2 = {
+classes_real_space = {
     "signal": CorrelScipySignal,
     "pycuda": CorrelPyCuda,
     "pythran": CorrelPythran,
@@ -34,7 +34,7 @@
     import pycuda
 except ImportError:
     classes.pop("pycuda")
-    classes2.pop("pycuda")
+    classes_real_space.pop("pycuda")
 
 try:
     import skcuda
@@ -86,14 +86,10 @@ def _test(self, cls=cls, k=k):
         ) = correl.compute_displacements_from_correl(c, norm=norm)
         displacement_computed = np.array([dx, dy])
 
-        logger.debug(
-            k
-            + ", displacement = [0, 0]\t error= {}\n".format(
-                abs(displacement_computed)
-            )
-        )
+        print(f"{k}, displacement = [0, 0]\t error= {abs(displacement_computed)}")
 
-        self.assertTrue(np.allclose([0, 0], displacement_computed, atol=1e-03))
+        assert np.allclose([0, 0], displacement_computed, atol=1e-03)
+        assert np.allclose(correl_max, 1.0), correl_max
 
         # then, with the 2 figures with displacements
         c, norm = correl(self.im0, self.im1)
@@ -159,7 +155,8 @@ def _test1(self, cls=cls, k=k):
             )
         )
 
-        self.assertTrue(np.allclose([0, 0], displacement_computed, atol=1e-03))
+        assert np.allclose(correl_max, 1.0), correl_max
+        assert np.allclose([0, 0], displacement_computed, atol=1e-03)
 
         # then, with the 2 figures with displacements
         c, norm = correl(self.im0, self.im1)
@@ -214,7 +211,7 @@ def setUpClass(cls):
         cls.im1 = cls.im1.astype("float32")
 
 
-for k, cls in classes2.items():
+for k, cls in classes_real_space.items():
 
     def _test2(self, cls=cls, k=k):
         correl = cls(self.im0.shape, self.im1.shape, mode="valid")
diff --git a/src/fluidimage/calcul/test_fft.py b/src/fluidimage/calcul/test_fft.py
index 0529b6a1..14124790 100644
--- a/src/fluidimage/calcul/test_fft.py
+++ b/src/fluidimage/calcul/test_fft.py
@@ -1,134 +1,89 @@
-import unittest
-from time import time
-
 import numpy as np
+import pytest
+from numpy.random import PCG64, Generator
 
-from fluidimage.calcul.correl import CUFFT2DReal2Complex, FFTW2DReal2Complex
-
-# from scipy.misc import lena
-# from correl import calcul_correl_norm_scipy, CorrelWithFFT
-
-
-class TestFFTW2DReal2Complex(unittest.TestCase):
-    # def test_correl(self):
-
-    #     rtime, ntime = 0., 0.
-
-    #     im1 = lena()
-    #     im2 = im1 + np.random.randn(*im1.shape) * 50  # add noise
-    #     im1 = im1.astype('float32')
-    #     im2 = im2.astype('float32')
-    #     nx, ny = im1.shape
-
-    #     t0 = time()
-    #     correl = CorrelWithFFT(nx, ny)
-    #     correlfft = correl.calcul_correl_norm(im1, im2)
-    #     ntime += time() - t0
+from fluidimage.calcul.fft import _compute_energy_from_fourier, classes
 
-    #     t0 = time()
-    #     correl_dir = calcul_correl_norm_scipy(im1, im2)
-    #     rtime += time() - t0
 
-    #     print 'correl fft speedup = %g' % (rtime / ntime)
-    #     # ax1 = plt.subplot(121)
-    #     # ax2 = plt.subplot(122)
-    #     # ax1.imshow(correlfft)
-    #     # ax2.imshow(correl_dir)
-    #     # plt.show(block=True)
+@pytest.mark.parametrize("cls", classes)
+def test_fft_random(cls):
 
-    #     rtol = 8e-0
-    #     atol = 3e-04
-    #     # tmp = np.absolute(correlfft - correl_dir) #- atol - rtol*np.abs(back)
+    nx, ny = 20, 24
+    oper = cls(nx, ny)
 
-    #     # print(tmp.max())
-    #     self.assertTrue(np.allclose(correlfft,
-    #                                 correl_dir, rtol=rtol, atol=atol))
+    generator = Generator(PCG64())
 
-    def test_fft(self):
-        """simple"""
-        nx = 4
-        ny = 2
-        op = FFTW2DReal2Complex(nx, ny)
+    arr = generator.random(nx * ny, dtype=oper.type_real).reshape(oper.shapeX)
+    arr_fft = oper.fft(arr)
 
-        func_fft = np.zeros(op.shapeK, dtype=op.type_complex)
-        func_fft[0, 1] = 1
+    energyX = oper.compute_energy_from_spatial(arr)
+    energyK = oper.compute_energy_from_fourier(arr_fft)
+    back = oper.ifft(arr_fft) / oper.coef_norm
+    energyX = oper.compute_energy_from_spatial(back)
+    arr_fft_2 = oper.fft(back)
+    energyKback = oper.compute_energy_from_fourier(arr_fft_2)
+    rtol = 8e-05
+    atol = 1e-04
+    assert np.allclose(arr_fft, arr_fft_2, rtol=rtol, atol=atol)
+    assert np.allclose(arr, back, rtol=rtol, atol=atol)
+    assert energyK == pytest.approx(energyKback)
+    assert energyX == pytest.approx(energyK)
 
-        self.compute_and_check(func_fft, op)
+    correl = oper.ifft(arr_fft.conj() * arr_fft)
+    assert np.allclose(correl.max() / (2 * energyK * oper.coef_norm_correl), 1)
 
-    def compute_and_check(self, func_fft, op):
-        energyK = op.compute_energy_from_Fourier(func_fft)
 
-        func = op.ifft(func_fft)
-        energyX = op.compute_energy_from_spatial(func)
+@pytest.mark.parametrize("cls", classes)
+def test_fft_simple(cls):
+    """simple"""
+    nx = 4
+    ny = 2
+    oper = cls(nx, ny)
 
-        back_fft = op.fft(func) / op.coef_norm
-        back = op.ifft(back_fft)
+    arr_fft = np.zeros(oper.shapeK, dtype=cls.type_complex)
+    arr_fft[0, 1] = 1
 
-        rtol = 8e-05
-        atol = 1e-04
-        self.assertTrue(np.allclose(func_fft, back_fft, rtol=rtol, atol=atol))
-        self.assertTrue(np.allclose(func, back, rtol=rtol, atol=atol))
+    energyK = oper.compute_energy_from_fourier(arr_fft)
 
-        self.assertAlmostEqual(energyX / energyK, 1.0, places=3)
+    func = oper.ifft(arr_fft)
+    energyX = oper.compute_energy_from_spatial(func)
 
-        energyKback = op.compute_energy_from_Fourier(back_fft)
-        self.assertAlmostEqual(energyK / energyKback, 1.0, places=3)
+    back_fft = oper.fft(func) / oper.coef_norm
+    back = oper.ifft(back_fft)
 
-    def compute_and_check2(self, func, op):
-        energyX = op.compute_energy_from_spatial(func)
-        func_fft = op.fft(func)
-        energyK = op.compute_energy_from_Fourier(func_fft)
-        back = op.ifft(func_fft)
-        energyX = op.compute_energy_from_spatial(back)
-        func_fft_2 = op.fft(back)
-        energyKback = op.compute_energy_from_Fourier(func_fft_2)
-        rtol = 8e-05
-        atol = 1e-04
-        self.assertTrue(np.allclose(func_fft, func_fft_2, rtol=rtol, atol=atol))
+    rtol = 8e-05
+    atol = 1e-04
+    assert np.allclose(arr_fft, back_fft, rtol=rtol, atol=atol)
+    assert np.allclose(func, back, rtol=rtol, atol=atol)
 
-        # tmp = np.absolute(func - back) - atol - rtol*np.abs(back)
+    assert energyX, pytest.approx(energyK)
 
-        # print(tmp.max())
+    energyKback = oper.compute_energy_from_fourier(back_fft)
 
-        self.assertTrue(np.allclose(func, back, rtol=rtol, atol=atol))
+    assert energyK, pytest.approx(energyKback)
 
-        self.assertAlmostEqual(energyX / energyK, 1.0, places=3)
-        self.assertAlmostEqual(energyK / energyKback, 1.0, places=3)
 
-    def bench_fft_random(self):
-        """random"""
-        nx = 128 * 16
-        ny = 64 * 4 * 8
+@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])
+def test_compute_energy_from_fourier(dtype):
 
-        rtime, ntime = 0.0, 0.0
-        Nloops = 2
-        for nloop in range(Nloops):
-            op = FFTW2DReal2Complex(nx, ny)
+    coef_norm = 10
+    n0 = 12
+    n1 = 8
 
-            func = np.random.random(op.shapeX)
+    generator = Generator(PCG64())
+    field_fft = generator.random(n0 * n1) + 1j * generator.random(n0 * n1)
+    field_fft = field_fft.reshape((n0, n1)).astype(np.complex64)
 
-            func = np.array(func, dtype=op.type_real)
+    assert field_fft.shape == (n0, n1)
 
-            t0 = time()
-            func_fft = op.fft(func)
-            func = op.ifft(func_fft)
-
-            self.compute_and_check2(func, op)
-            ntime += time() - t0
-
-            op = CUFFT2DReal2Complex(nx, ny)
-
-            t0 = time()
-            func_fft = op.fft(func)
-            func1 = op.ifft(func_fft)
-
-            self.compute_and_check2(func1, op)
-            rtime += time() - t0
-
-        print(
-            "array size = %5d x %5d : gpu speedup = %g" % (nx, ny, ntime / rtime)
+    expected = (
+        0.5
+        / coef_norm
+        * (
+            np.sum(abs(field_fft[:, 0]) ** 2 + abs(field_fft[:, -1]) ** 2)
+            + 2 * np.sum(abs(field_fft[:, 1:-1]) ** 2)
         )
+    )
 
-
-if __name__ == "__main__":
-    unittest.main()
+    energy = _compute_energy_from_fourier(field_fft, coef_norm)
+    assert energy == pytest.approx(expected)
diff --git a/src/fluidimage/gui/imviewer.py b/src/fluidimage/gui/imviewer.py
index dd93e679..4b9b27c8 100644
--- a/src/fluidimage/gui/imviewer.py
+++ b/src/fluidimage/gui/imviewer.py
@@ -15,6 +15,7 @@
 
 import matplotlib.pyplot as plt
 
+import fluidimage
 from fluiddyn.io.image import imread
 from fluiddyn.util import time_as_str
 from fluiddyn.util.serieofarrays import SerieOfArraysFromFiles
@@ -62,7 +63,12 @@ def parse_args():
         "-cm", "--colormap", help="colormap", type=str, default=""
     )
     parser.add_argument("-v", "--verbose", help="verbose mode", action="count")
-
+    parser.add_argument(
+        "-V",
+        "--version",
+        help="Print fluidimage version and exit",
+        action="count",
+    )
     return parser.parse_args()
 
 
@@ -309,4 +315,9 @@ def onclick(self, event):
 
 def main():
     args = parse_args()
+
+    if args.version:
+        print(f"fluidimage {fluidimage.__version__}")
+        return
+
     ImageViewer(args)
diff --git a/src/fluidimage/gui/monitor.py b/src/fluidimage/gui/monitor.py
index cc983a6e..d6085d8f 100644
--- a/src/fluidimage/gui/monitor.py
+++ b/src/fluidimage/gui/monitor.py
@@ -30,6 +30,7 @@
     Tree,
 )
 
+import fluidimage
 from fluidimage import ParamContainer
 
 
@@ -98,7 +99,12 @@ def parse_args(cls):
         parser.add_argument(
             "-v", "--verbose", help="verbose mode", action="count"
         )
-
+        parser.add_argument(
+            "-V",
+            "--version",
+            help="Print fluidimage version and exit",
+            action="count",
+        )
         return parser.parse_args()
 
     def __init__(self, args):
@@ -270,5 +276,10 @@ def on_tree_node_selected(self, event):
 def main():
     """Main function for fluidimage-monitor"""
     args = MonitorApp.parse_args()
+
+    if args.version:
+        print(f"fluidimage {fluidimage.__version__}")
+        return
+
     app = MonitorApp(args)
     app.run()
diff --git a/src/fluidimage/gui/piv_viewer.py b/src/fluidimage/gui/piv_viewer.py
index bcc13dce..5ba0ef0b 100644
--- a/src/fluidimage/gui/piv_viewer.py
+++ b/src/fluidimage/gui/piv_viewer.py
@@ -16,6 +16,7 @@
 import h5py
 import matplotlib.pyplot as plt
 
+import fluidimage
 from fluiddyn.util.serieofarrays import SerieOfArraysFromFiles
 from fluidimage.gui.base_matplotlib import AppMatplotlibWidgets
 
@@ -84,7 +85,12 @@ def parse_args(cls):
         parser.add_argument(
             "-v", "--verbose", help="verbose mode", action="count"
         )
-
+        parser.add_argument(
+            "-V",
+            "--version",
+            help="Print fluidimage version and exit",
+            action="count",
+        )
         return parser.parse_args()
 
     def __init__(self, args):
@@ -208,4 +214,9 @@ def _update_fig(self):
 
 def main():
     args = VectorFieldsViewer.parse_args()
+
+    if args.version:
+        print(f"fluidimage {fluidimage.__version__}")
+        return
+
     return VectorFieldsViewer(args)
diff --git a/src/fluidimage/gui/test_imviewer.py b/src/fluidimage/gui/test_imviewer.py
index d1b24f4d..cd30fe01 100644
--- a/src/fluidimage/gui/test_imviewer.py
+++ b/src/fluidimage/gui/test_imviewer.py
@@ -1,42 +1,46 @@
 import sys
-import unittest
 
-from fluiddyn.io import stdout_redirected
 from fluidimage import get_path_image_samples
-from fluidimage.gui.imviewer import ImageViewer, parse_args
+from fluidimage.gui.imviewer import ImageViewer, main, parse_args
 
 path_image_samples = get_path_image_samples()
 
 
-class TestImageViewer(unittest.TestCase):
-    def test_main(self):
-        with stdout_redirected():
-            command = "fluidimviewer"
-            args = command.split()
-            args.append(str(path_image_samples / "Karman/Images"))
-            sys.argv = args
-            args = parse_args()
-            ImageViewer(args)
+def test_fluidimviewer_version(monkeypatch):
+    command = "fluidimviewer --version"
+    with monkeypatch.context() as ctx:
+        ctx.setattr(sys, "argv", command.split())
+        main()
 
-            args = command.split()
-            args.append(str(path_image_samples / "Karman/Images/*"))
-            sys.argv = args
 
-            args = parse_args()
-            self = ImageViewer(args)
+def test_main(monkeypatch):
 
-            self.set_autoclim(None)
+    words = ["fluidimviewer", str(path_image_samples / "Karman/Images")]
+    with monkeypatch.context() as ctx:
+        ctx.setattr(sys, "argv", words)
+        args = parse_args()
 
-            self._switch()
-            self._switch()
+    ImageViewer(args)
 
-            self._increase_ifile()
-            self._decrease_ifile()
+    words = ["fluidimviewer", str(path_image_samples / "Karman/Images/*")]
+    with monkeypatch.context() as ctx:
+        ctx.setattr(sys, "argv", words)
+        args = parse_args()
 
-            self._submit_n("2")
+    viewer = ImageViewer(args)
 
-            self._increase_ifile_n()
-            self._decrease_ifile_n()
+    viewer.set_autoclim(None)
 
-            self._change_cmin("1")
-            self._change_cmax("2")
+    viewer._switch()
+    viewer._switch()
+
+    viewer._increase_ifile()
+    viewer._decrease_ifile()
+
+    viewer._submit_n("2")
+
+    viewer._increase_ifile_n()
+    viewer._decrease_ifile_n()
+
+    viewer._change_cmin("1")
+    viewer._change_cmax("2")
diff --git a/src/fluidimage/gui/test_monitor.py b/src/fluidimage/gui/test_monitor.py
index 7b3cacf8..c159ede6 100644
--- a/src/fluidimage/gui/test_monitor.py
+++ b/src/fluidimage/gui/test_monitor.py
@@ -24,6 +24,14 @@ def test_monitor_help(monkeypatch):
             raise RuntimeError
 
 
+def test_monitor_version(monkeypatch):
+    command = "fluidimage-monitor --version"
+
+    with monkeypatch.context() as ctx:
+        ctx.setattr(sys, "argv", command.split())
+        main()
+
+
 def test_monitor_bad_paths(monkeypatch, tmp_path):
 
     command = "fluidimage-monitor __dir_does_not_exist__"
diff --git a/src/fluidimage/gui/test_piv_viewer.py b/src/fluidimage/gui/test_piv_viewer.py
index e38d92a6..df04b613 100644
--- a/src/fluidimage/gui/test_piv_viewer.py
+++ b/src/fluidimage/gui/test_piv_viewer.py
@@ -4,6 +4,13 @@
 from fluidimage.gui.piv_viewer import main
 
 
+def test_fluidpivviewer_version(monkeypatch):
+    command = "fluidpivviewer --version"
+    with monkeypatch.context() as ctx:
+        ctx.setattr(sys, "argv", command.split())
+        main()
+
+
 def test_main(monkeypatch):
 
     path_image_samples = get_path_image_samples()
diff --git a/src/fluidimage/topologies/test_bos.py b/src/fluidimage/topologies/test_bos.py
index 32a198e4..5452072f 100644
--- a/src/fluidimage/topologies/test_bos.py
+++ b/src/fluidimage/topologies/test_bos.py
@@ -25,7 +25,7 @@ def test_bos(tmp_path_karman, executor):
     params.piv0.shape_crop_im0 = 16
 
     # compute only few vectors
-    params.piv0.grid.overlap = -8
+    params.piv0.grid.overlap = -2
 
     params.saving.how = "recompute"
     params.saving.postfix = postfix
diff --git a/src/fluidimage/works/piv/fix.py b/src/fluidimage/works/piv/fix.py
index 155e55a9..6457dc10 100644
--- a/src/fluidimage/works/piv/fix.py
+++ b/src/fluidimage/works/piv/fix.py
@@ -7,12 +7,14 @@
 
 """
 
+from math import sqrt
+
 import numpy as np
 
 from fluiddyn.util.paramcontainer import ParamContainer
 from fluidimage.calcul.errors import PIVError
+from fluidimage.calcul.mean_neighbors import mean_neighbors_xy
 
-from ...calcul.smooth_clean import smooth_clean
 from .. import BaseWork
 
 
@@ -104,32 +106,28 @@ def put_to_nan(inds, explanation):
             threshold = self.params_fix.threshold_diff_neighbour
             ixvecs = self.piv_work.ixvecs
             iyvecs = self.piv_work.iyvecs
-            xs = piv_results.xs
-            ys = piv_results.ys
 
             ixvecs, iyvecs = self.piv_work._xyoriginalimage_from_xymasked(
                 ixvecs, iyvecs
             )
 
-            dxs_smooth, dys_smooth = smooth_clean(
-                xs, ys, deltaxs, deltays, iyvecs, ixvecs, threshold
+            dxs_neighbors, dys_neighbors = mean_neighbors_xy(
+                deltaxs, deltays, iyvecs, ixvecs
             )
-            piv_results.dxs_smooth_clean = dxs_smooth
-            piv_results.dys_smooth_clean = dys_smooth
 
-            differences = np.sqrt(
-                (dxs_smooth - deltaxs) ** 2 + (dys_smooth - deltays) ** 2
-            )
+            diff_x = dxs_neighbors - deltaxs
+            diff_y = dys_neighbors - deltays
+            differences2 = diff_x**2 + diff_y**2
 
             with np.errstate(invalid="ignore"):
-                inds = (differences > threshold).nonzero()[0]
+                inds = (differences2 > threshold**2).nonzero()[0]
 
             put_to_nan(inds, "diff neighbour too large")
 
             for ivec in inds:
-                piv_results.errors[ivec] += " (diff = {:.2f})".format(
-                    differences[ivec]
-                )
+                piv_results.errors[
+                    ivec
+                ] += f" (diff = {sqrt(differences2[ivec]):.2f})"
 
             piv_results.deltaxs_wrong = deltaxs_wrong
             piv_results.deltays_wrong = deltays_wrong
@@ -169,11 +167,12 @@ def put_to_nan(inds, explanation):
                 continue
 
             diff_neighbours = np.empty(len(other_peaks_good) + 1)
-            diff_neighbours[0] = differences[ivec]
+            diff_neighbours[0] = sqrt(differences_2[ivec])
 
             for i, (dx, dy, corr) in enumerate(other_peaks_good):
                 diff_neighbours[i + 1] = np.sqrt(
-                    (dxs_smooth[ivec] - dx) ** 2 + (dys_smooth[ivec] - dy) ** 2
+                    (dxs_neighbors[ivec] - dx) ** 2
+                    + (dys_neighbors[ivec] - dy) ** 2
                 )
 
             argmin = diff_neighbours.argmin()
diff --git a/src/fluidimage/works/piv/singlepass.py b/src/fluidimage/works/piv/singlepass.py
index 45a02884..a8dd54f0 100644
--- a/src/fluidimage/works/piv/singlepass.py
+++ b/src/fluidimage/works/piv/singlepass.py
@@ -155,8 +155,10 @@ def calcul(self, couple):
         couple.apply_mask(self.params.mask)
 
         im0, im1 = couple.get_arrays()
+
         if not hasattr(self, "ixvecs_grid"):
             self._prepare_with_image(im0)
+
         (
             deltaxs,
             deltays,
@@ -201,8 +203,8 @@ def _pad_images(self, im0, im1):
         """
         npad = self.npad = max(self._start_for_crop0 + self._stop_for_crop0)
         tmp = [(npad, npad), (npad, npad)]
-        im0pad = np.pad(im0 - im0.min(), tmp, "constant")
-        im1pad = np.pad(im1 - im1.min(), tmp, "constant")
+        im0pad = np.pad(im0 - im0.min(), tmp, "constant", constant_values=0)
+        im1pad = np.pad(im1 - im1.min(), tmp, "constant", constant_values=0)
         return im0pad, im1pad
 
     def _calcul_positions_vectors_subimages(
@@ -308,22 +310,8 @@ def _loop_vectors(self, im0, im1, deltaxs_input=None, deltays_input=None):
                 errors[ivec] = "Bad im_crop shape."
                 continue
 
-            if np.isnan(im0crop).any() or np.isnan(im1crop).any():
-                deltaxs[ivec] = np.nan
-                deltays[ivec] = np.nan
-                correls_max[ivec] = np.nan
-                errors[ivec] = "Nan(s) in cropped images."
-                continue
-
             # compute and store correlation map
             correl, norm = self.correl(im0crop, im1crop)
-            if (
-                self.index_pass == 0
-                and self.params.piv0.coef_correl_no_displ is not None
-            ):
-                correl[
-                    self.correl.get_indices_no_displacement()
-                ] *= self.params.piv0.coef_correl_no_displ
 
             correls[ivec] = correl
 
@@ -401,16 +389,13 @@ def _init_crop(self):
         self._stop_for_crop1 = tuple(_stop_for_crop1)
 
     def _crop_im0(self, ixvec, iyvec, im):
-        """Crop image 0.
-
-        Warning: important for perf (~12% for PIV with _crop_im1)
-        """
+        """Crop image 0."""
         subim = im[
             iyvec - self._start_for_crop0[0] : iyvec + self._stop_for_crop0[0],
             ixvec - self._start_for_crop0[1] : ixvec + self._stop_for_crop0[1],
         ]
-        subim = np.array(subim, dtype=np.float32)
-        return subim - subim.mean()
+        subim = np.ascontiguousarray(subim, dtype=np.float32)
+        return subim
 
     def _crop_im1(self, ixvec, iyvec, im):
         """Crop image 1."""
@@ -418,8 +403,8 @@ def _crop_im1(self, ixvec, iyvec, im):
             iyvec - self._start_for_crop1[0] : iyvec + self._stop_for_crop1[0],
             ixvec - self._start_for_crop1[1] : ixvec + self._stop_for_crop1[1],
         ]
-        subim = np.array(subim, dtype=np.float32)
-        return subim - subim.mean()
+        subim = np.ascontiguousarray(subim, dtype=np.float32)
+        return subim
 
     def apply_interp(self, piv_results, last=False):
         """Interpolate a PIV result object on the grid of the PIV work.
@@ -592,7 +577,6 @@ def _complete_params_with_default(cls, params):
                 "method_correl": "fftw",
                 "method_subpix": "2d_gaussian2",
                 "nsubpix": None,
-                "coef_correl_no_displ": None,
                 "nb_peaks_to_search": 1,
                 "particle_radius": 3,
             },
@@ -639,11 +623,6 @@ def _complete_params_with_default(cls, params):
   particles. It has to be increased in case of peak locking (plot the
   histograms of the displacements).
 
-- coef_correl_no_displ : None, number
-
-  If this coefficient is not None, the correlation of the point corresponding
-  to no displacement is multiplied by this coefficient (for the first pass).
-
 - nb_peaks_to_search : 1, int
 
   Number of peaks to search. Secondary peaks can be used during the fix step.
diff --git a/src/fluidimage/works/piv/test_piv.py b/src/fluidimage/works/piv/test_piv.py
index 382dbfae..7e3db013 100644
--- a/src/fluidimage/works/piv/test_piv.py
+++ b/src/fluidimage/works/piv/test_piv.py
@@ -25,7 +25,7 @@ def test_minimal_piv(tmp_path):
 
     params.multipass.number = 2
 
-    params.fix.displacement_max = 2
+    params.fix.displacement_max = 3
     params.fix.threshold_diff_neighbour = 2
 
     params.series.path = str(path_images / "Oseen*")