# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1713800256 -7200
#      Mon Apr 22 17:37:36 2024 +0200
# Node ID a17adc60f680bfbaa9cc500b0cb01c7cb3bf7600
# Parent  d1d2300c6ef2d3eaa0b0bfff2a85d6bb5f0694ce
TPS: self.norm_coefs

diff --git a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
--- a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
+++ b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
@@ -21,6 +21,7 @@
 """
 
 from logging import debug
+from typing import List
 
 import numpy as np
 from transonic import Array
@@ -31,6 +32,10 @@
 class ThinPlateSplineSubdom:
     """Helper class for thin plate interpolation."""
 
+    new_positions: "np.int64[:,:]"
+    ind_new_positions_subdom: List["np.int64[:]"]
+    norm_coefs: "float[:]"
+
     def __init__(
         self,
         centers,
@@ -137,6 +142,15 @@
         return U_smooth_tmp, U_tps, summary
 
     def init_with_new_positions(self, new_positions):
+        """Initialize with the new positions
+
+        Parameters
+        ----------
+
+        new_positions: 2d array of int64
+          new_positions[0] and new_positions[1] correspond to the x and y values, respectively.
+
+        """
         npos = self.new_positions = new_positions
 
         ind_new_positions_subdom = []
@@ -161,6 +175,15 @@
                 i_subdom += 1
 
         self.ind_new_positions_subdom = ind_new_positions_subdom
+
+        self.norm_coefs = np.zeros(self.new_positions.shape[1])
+        for i_subdom in range(self.nb_subdom):
+            # TODO: replace 1 by an appropriate function of
+            # ind_new_positions_subdom[i_subdom]
+            # + save another list of 1d array like ind_new_positions_subdom
+            # containing the norm coefficients for each subdomain
+            self.norm_coefs[ind_new_positions_subdom[i_subdom]] += 1
+
         self._init_EM_subdom()
 
     def _init_EM_subdom(self):
@@ -176,16 +199,14 @@
         self.EM = EM
 
     def compute_eval(self, U_tps):
-        U_eval = np.zeros(self.new_positions[1].shape)
-        nb_tps = np.zeros(self.new_positions[1].shape, dtype=int)
+        U_eval = np.zeros(self.new_positions.shape[1])
 
         for i in range(self.nb_subdom):
             U_eval[self.ind_new_positions_subdom[i]] += np.dot(
                 U_tps[i], self.EM[i]
             )
-            nb_tps[self.ind_new_positions_subdom[i]] += 1
 
-        U_eval /= nb_tps
+        U_eval /= self.norm_coefs
 
         return U_eval
 
# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1714121049 -7200
#      Fri Apr 26 10:44:09 2024 +0200
# Node ID 04051f3e2f7dda60770bcbbc0db812c6c684d994
# Parent  a17adc60f680bfbaa9cc500b0cb01c7cb3bf7600
Start to clean up calcul/interpolate/thin_plate_spline_subdom.py

diff --git a/doc/examples/piv_try_params_Karman.py b/doc/examples/piv_try_params_Karman.py
--- a/doc/examples/piv_try_params_Karman.py
+++ b/doc/examples/piv_try_params_Karman.py
@@ -6,7 +6,7 @@
 
 params.series.path = "../../image_samples/Karman/Images"
 
-params.mask.strcrop = "50:350, 0:380"
+params.mask.strcrop = "20:380, 0:450"
 
 params.fix.correl_min = 0.4
 params.fix.threshold_diff_neighbour = 2.0
@@ -17,12 +17,14 @@
 
 params.multipass.number = 2
 params.multipass.use_tps = "last"
+params.multipass.subdom_size = 400
+params.multipass.smoothing_coef = 2.0
 
 work = Work(params=params)
 
 piv = work.process_1_serie()
 
-# piv.display(show_interp=True, scale=0.3, show_error=False)
-piv.display(show_interp=False, scale=1, show_error=True)
+piv.display(show_interp=True, scale=0.3, show_error=False, show_correl=False)
+# piv.display(show_interp=False, scale=1, show_error=True)
 
 # piv.save()
diff --git a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
--- a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
+++ b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
@@ -42,7 +42,7 @@
         subdom_size,
         smoothing_coef,
         threshold=None,
-        percent_buffer_area=0.2,
+        percent_buffer_area=20,
     ):
         self.centers = centers
         self.subdom_size = subdom_size
@@ -50,13 +50,17 @@
         self.threshold = threshold
         self.compute_indices(percent_buffer_area)
 
-    def compute_indices(self, percent_buffer_area=0.25):
+    def compute_indices(self, percent_buffer_area):
+        # note: `centers = np.vstack([ys, xs])`
         xs = self.centers[1]
         ys = self.centers[0]
-        max_coord = np.max(self.centers, 1)
-        min_coord = np.min(self.centers, 1)
-        range_coord = max_coord - min_coord
-        aspect_ratio = range_coord[0] / range_coord[1]
+        x_min = xs.min()
+        y_min = ys.min()
+        x_max = xs.max()
+        y_max = ys.max()
+        range_x = x_max - x_min
+        range_y = y_max - y_min
+        aspect_ratio = range_y / range_x
 
         nb_subdom = xs.size / self.subdom_size
         nb_subdomx = int(np.floor(np.sqrt(nb_subdom / aspect_ratio)))
@@ -71,45 +75,33 @@
         self.nb_subdomx = nb_subdomx
         self.nb_subdomy = nb_subdomy
 
-        x_dom = np.linspace(min_coord[1], max_coord[1], nb_subdomx + 1)
-        y_dom = np.linspace(min_coord[0], max_coord[0], nb_subdomy + 1)
+        x_dom = np.linspace(x_min, x_max, nb_subdomx + 1)
+        y_dom = np.linspace(y_min, y_max, nb_subdomy + 1)
+
+        coef_buffer = percent_buffer_area / 100
+        buffer_length_x = coef_buffer * range_x / nb_subdomx
+        buffer_length_y = coef_buffer * range_y / nb_subdomy
 
-        buffer_area_x = (
-            range_coord[1]
-            / (nb_subdomx)
-            * percent_buffer_area
-            * np.ones_like(x_dom)
-        )
-        buffer_area_y = (
-            range_coord[0]
-            / (nb_subdomy)
-            * percent_buffer_area
-            * np.ones_like(y_dom)
-        )
+        self.xmin_limits = x_dom[:-1] - buffer_length_x
+        self.xmax_limits = x_dom[1:] + buffer_length_x
 
-        self.x_dom = x_dom
-        self.y_dom = y_dom
-        self.buffer_area_x = buffer_area_x
-        self.buffer_area_y = buffer_area_y
+        self.ymin_limits = y_dom[:-1] - buffer_length_y
+        self.ymax_limits = y_dom[1:] + buffer_length_y
 
-        ind_subdom = np.zeros([nb_subdom, 2])
         ind_v_subdom = []
 
-        i_subdom = 0
-        for i in range(nb_subdomx):
-            for j in range(nb_subdomy):
-                ind_subdom[i_subdom, :] = [i, j]
+        for iy in range(nb_subdomy):
+            for ix in range(nb_subdomx):
 
                 ind_v_subdom.append(
                     np.where(
-                        (xs >= x_dom[i] - buffer_area_x[i])
-                        & (xs < x_dom[i + 1] + buffer_area_x[i + 1])
-                        & (ys >= y_dom[j] - buffer_area_y[j])
-                        & (ys < y_dom[j + 1] + buffer_area_y[j + 1])
+                        (xs >= self.xmin_limits[ix])
+                        & (xs < self.xmax_limits[ix])
+                        & (ys >= self.ymin_limits[iy])
+                        & (ys < self.ymin_limits[iy])
                     )[0]
                 )
 
-                i_subdom += 1
         self.ind_v_subdom = ind_v_subdom
         self.nb_subdom = nb_subdom
 
@@ -120,7 +112,6 @@
 
         for i in range(self.nb_subdom):
             centers_tmp = self.centers[:, self.ind_v_subdom[i]]
-
             U_tmp = U[self.ind_v_subdom[i]]
             U_smooth[i], U_tps[i], summaries[i] = self.compute_tps_coeff_iter(
                 centers_tmp, U_tmp
@@ -148,32 +139,26 @@
         ----------
 
         new_positions: 2d array of int64
-          new_positions[0] and new_positions[1] correspond to the x and y values, respectively.
+          new_positions[1] and new_positions[0] correspond to the x and y values, respectively.
 
         """
-        npos = self.new_positions = new_positions
+        self.new_positions = new_positions
+        xs = new_positions[1]
+        ys = new_positions[0]
 
         ind_new_positions_subdom = []
 
-        x_dom = self.x_dom
-        y_dom = self.y_dom
-        buffer_area_x = self.buffer_area_x
-        buffer_area_y = self.buffer_area_y
-
-        i_subdom = 0
-        for i in range(self.nb_subdomx):
-            for j in range(self.nb_subdomy):
+        for iy in range(self.nb_subdomy):
+            for ix in range(self.nb_subdomx):
                 ind_new_positions_subdom.append(
                     np.where(
-                        (npos[1] >= x_dom[i] - buffer_area_x[i])
-                        & (npos[1] < x_dom[i + 1] + buffer_area_x[i + 1])
-                        & (npos[0] >= y_dom[j] - buffer_area_y[j])
-                        & (npos[0] < y_dom[j + 1] + buffer_area_y[j + 1])
+                        (xs >= self.xmin_limits[ix])
+                        & (xs < self.xmax_limits[ix])
+                        & (ys >= self.ymin_limits[iy])
+                        & (ys < self.ymin_limits[iy])
                     )[0]
                 )
 
-                i_subdom += 1
-
         self.ind_new_positions_subdom = ind_new_positions_subdom
 
         self.norm_coefs = np.zeros(self.new_positions.shape[1])
@@ -184,21 +169,16 @@
             # containing the norm coefficients for each subdomain
             self.norm_coefs[ind_new_positions_subdom[i_subdom]] += 1
 
-        self._init_EM_subdom()
-
-    def _init_EM_subdom(self):
         EM = [None] * self.nb_subdom
 
         for i in range(self.nb_subdom):
             centers_tmp = self.centers[:, self.ind_v_subdom[i]]
-            new_positions_tmp = self.new_positions[
-                :, self.ind_new_positions_subdom[i]
-            ]
+            new_positions_tmp = new_positions[:, ind_new_positions_subdom[i]]
             EM[i] = compute_tps_matrix(new_positions_tmp, centers_tmp)
 
         self.EM = EM
 
-    def compute_eval(self, U_tps):
+    def interpolate(self, U_tps):
         U_eval = np.zeros(self.new_positions.shape[1])
 
         for i in range(self.nb_subdom):
diff --git a/src/fluidimage/works/piv/multipass.py b/src/fluidimage/works/piv/multipass.py
--- a/src/fluidimage/works/piv/multipass.py
+++ b/src/fluidimage/works/piv/multipass.py
@@ -62,7 +62,7 @@
                 "coeff_zoom": 2,
                 "use_tps": "last",
                 "subdom_size": 200,
-                "smoothing_coef": 0.5,
+                "smoothing_coef": 2.0,
                 "threshold_tps": 1.0,
             },
         )
diff --git a/src/fluidimage/works/piv/singlepass.py b/src/fluidimage/works/piv/singlepass.py
--- a/src/fluidimage/works/piv/singlepass.py
+++ b/src/fluidimage/works/piv/singlepass.py
@@ -483,7 +483,7 @@
                 subdom_size,
                 smoothing_coef,
                 threshold=threshold,
-                percent_buffer_area=0.25,
+                percent_buffer_area=20,
             )
             try:
                 (
@@ -514,8 +514,8 @@
 
                 tps.init_with_new_positions(new_positions)
 
-                deltaxs_approx = tps.compute_eval(deltaxs_tps)
-                deltays_approx = tps.compute_eval(deltays_tps)
+                deltaxs_approx = tps.interpolate(deltaxs_tps)
+                deltays_approx = tps.interpolate(deltays_tps)
 
                 print("TPS summary: ", end="")
                 nb_fixed_vectors_tot = summary.pop("nb_fixed_vectors_tot")
# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1714134511 -7200
#      Fri Apr 26 14:28:31 2024 +0200
# Node ID 8dca5014bda76cd23bd656815219ba9eaef83977
# Parent  04051f3e2f7dda60770bcbbc0db812c6c684d994
Refactoring and renaming TPS

diff --git a/src/fluidimage/calcul/interpolate/test_thin_plate_spline.py b/src/fluidimage/calcul/interpolate/test_thin_plate_spline.py
--- a/src/fluidimage/calcul/interpolate/test_thin_plate_spline.py
+++ b/src/fluidimage/calcul/interpolate/test_thin_plate_spline.py
@@ -3,7 +3,7 @@
 from .thin_plate_spline import (
     ThinPlateSpline,
     ThinPlateSplineNumpy,
-    compute_tps_coeff,
+    compute_tps_weights,
 )
 
 pi = np.pi
@@ -29,7 +29,7 @@
     new_positions = np.vstack([XI, YI])
 
     # calculate tps coeff
-    U_smooth, U_tps = compute_tps_coeff(centers, U, 0)
+    U_smooth, U_tps = compute_tps_weights(centers, U, 0)
     # evaluate interpolation on the new grid
     tps = ThinPlateSpline(new_positions, centers)
     tps.compute_field(U_tps)
diff --git a/src/fluidimage/calcul/interpolate/thin_plate_spline.py b/src/fluidimage/calcul/interpolate/thin_plate_spline.py
--- a/src/fluidimage/calcul/interpolate/thin_plate_spline.py
+++ b/src/fluidimage/calcul/interpolate/thin_plate_spline.py
@@ -8,14 +8,14 @@
 squared difference from the initial data.
 
 We first need to compute tps coefficients ``U_tps`` (function
-``compute_tps_coeff``). Interpolated data can then be obtained as the
+``compute_tps_weights``). Interpolated data can then be obtained as the
 matrix product ``dot(U_tps, EM)`` where the matrix ``EM`` is obtained
 by the function ``compute_tps_matrix``.  The spatial derivatives are
 obtained as ``dot(U_tps, EMDX)`` and ``dot(U_tps, EMDY)``, where
 ``EMDX`` and ``EMDY`` are obtained from the function
 ``compute_tps_matrix_dxy``. A helper class is also provided.
 
-.. autofunction:: compute_tps_coeff
+.. autofunction:: compute_tps_weights
 
 .. autoclass:: ThinPlateSpline
    :members:
@@ -165,7 +165,7 @@
     compute_tps_matrix = compute_tps_matrix_numpy
 
 
-def compute_tps_coeff(centers, U, smoothing_coef):
+def compute_tps_weights(centers, values, smoothing_coef):
     """Calculate the thin plate spline (tps) coefficients
 
     Parameters
@@ -175,7 +175,7 @@
         ``[nb_dim,  N]`` array representing the positions of the N centers,
         sources of the TPS (nb_dim = space dimension).
 
-    U : np.array
+    values : np.array
         ``[N]`` array representing the values of the considered
         scalar measured at the centres ``centers``.
 
@@ -185,25 +185,25 @@
     Returns
     -------
 
-    U_smooth : np.array
+    smoothed_values : np.array
          Values of the quantity U at the N centres after smoothing.
 
-    U_tps : np.array
+    tps_weights : np.array
          TPS weights of the centres and columns of the linear.
 
     """
-    nb_dim, N = centers.shape
-    U = np.hstack([U, np.zeros(nb_dim + 1)])
-    U = U.reshape([U.size, 1])
+    nb_dim, num_values = centers.shape
+    values = np.hstack([values, np.zeros(nb_dim + 1)])
+    values = values.reshape([values.size, 1])
     try:
         EM = compute_tps_matrix(centers, centers).T
     except TypeError as e:
         print(centers.dtype, centers.shape)
         raise e
 
-    smoothing_mat = smoothing_coef * np.eye(N, N)
-    smoothing_mat = np.hstack([smoothing_mat, np.zeros([N, nb_dim + 1])])
-    PM = np.hstack([np.ones([N, 1]), centers.T])
+    smoothing_mat = smoothing_coef * np.eye(num_values, num_values)
+    smoothing_mat = np.hstack([smoothing_mat, np.zeros([num_values, nb_dim + 1])])
+    PM = np.hstack([np.ones([num_values, 1]), centers.T])
     IM = np.vstack(
         [
             EM + smoothing_mat,
@@ -211,9 +211,9 @@
         ]
     )
 
-    U_tps = np.linalg.solve(IM, U)
-    U_smooth = np.dot(EM, U_tps)
-    return U_smooth.ravel(), U_tps.ravel()
+    tps_weights = np.linalg.solve(IM, values)
+    smoothed_values = np.dot(EM, tps_weights)
+    return smoothed_values.ravel(), tps_weights.ravel()
 
 
 def compute_tps_matrices_dxy(dsites, centers):
diff --git a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
--- a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
+++ b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
@@ -7,14 +7,6 @@
 1988) minimises a linear combination of the squared curvature and
 squared difference from the initial data.
 
-We first need to compute tps coefficients ``U_tps`` (function
-``compute_tps_coeff``). Interpolated data can then be obtained as the
-matrix product ``dot(U_tps, EM)`` where the matrix ``EM`` is obtained
-by the function ``compute_tps_matrix``.  The spatial derivatives are
-obtained as ``dot(U_tps, EMDX)`` and ``dot(U_tps, EMDY)``, where
-``EMDX`` and ``EMDY`` are obtained from the function
-``compute_tps_matrix_dxy``. A helper class is also provided.
-
 .. autoclass:: ThinPlateSplineSubdom
    :members:
 
@@ -26,15 +18,17 @@
 import numpy as np
 from transonic import Array
 
-from .thin_plate_spline import compute_tps_coeff, compute_tps_matrix
+from .thin_plate_spline import compute_tps_matrix, compute_tps_weights
 
 
 class ThinPlateSplineSubdom:
     """Helper class for thin plate interpolation."""
 
-    new_positions: "np.int64[:,:]"
-    ind_new_positions_subdom: List["np.int64[:]"]
+    ind_new_positions_domains: List["np.int64[:]"]
     norm_coefs: "float[:]"
+    num_new_positions: int
+    num_centers: int
+    tps_matrices: List["float[:,:]"]
 
     def __init__(
         self,
@@ -46,18 +40,14 @@
     ):
         self.centers = centers
         self.subdom_size = subdom_size
-        self.smoothing_coef = smoothing_coef
         self.threshold = threshold
-        self.compute_indices(percent_buffer_area)
 
-    def compute_indices(self, percent_buffer_area):
         # note: `centers = np.vstack([ys, xs])`
         xs = self.centers[1]
         ys = self.centers[0]
-        x_min = xs.min()
-        y_min = ys.min()
-        x_max = xs.max()
-        y_max = ys.max()
+        self.num_centers = xs.size
+        x_min, x_max = xs.min(), xs.max()
+        y_min, y_max = ys.min(), ys.max()
         range_x = x_max - x_min
         range_y = y_max - y_min
         aspect_ratio = range_y / range_x
@@ -70,14 +60,19 @@
 
         debug(f"nb_subdomx: {nb_subdomx} ; nb_subdomy: {nb_subdomy}")
 
-        nb_subdom = nb_subdomx * nb_subdomy
-
         self.nb_subdomx = nb_subdomx
         self.nb_subdomy = nb_subdomy
+        self.nb_subdom = nb_subdomx * nb_subdomy
 
         x_dom = np.linspace(x_min, x_max, nb_subdomx + 1)
         y_dom = np.linspace(y_min, y_max, nb_subdomy + 1)
 
+        # normalization as UVmat so that the effect of the filter do not depends
+        # too much on the size of the domains
+        self.smoothing_coef = (
+            smoothing_coef * (x_dom[1] - x_dom[0]) * (y_dom[1] - y_dom[0]) / 1000
+        )
+
         coef_buffer = percent_buffer_area / 100
         buffer_length_x = coef_buffer * range_x / nb_subdomx
         buffer_length_y = coef_buffer * range_y / nb_subdomy
@@ -88,12 +83,10 @@
         self.ymin_limits = y_dom[:-1] - buffer_length_y
         self.ymax_limits = y_dom[1:] + buffer_length_y
 
-        ind_v_subdom = []
-
+        self.indices_domains = []
         for iy in range(nb_subdomy):
             for ix in range(nb_subdomx):
-
-                ind_v_subdom.append(
+                self.indices_domains.append(
                     np.where(
                         (xs >= self.xmin_limits[ix])
                         & (xs < self.xmax_limits[ix])
@@ -102,35 +95,39 @@
                     )[0]
                 )
 
-        self.ind_v_subdom = ind_v_subdom
-        self.nb_subdom = nb_subdom
-
-    def compute_tps_coeff_subdom(self, U):
-        U_smooth = [None] * self.nb_subdom
-        U_tps = [None] * self.nb_subdom
+    def compute_tps_weights_subdom(self, values):
+        """Compute the TPS weights for all subdomains"""
+        smoothed_field_domains = [None] * self.nb_subdom
+        weights_domains = [None] * self.nb_subdom
         summaries = [None] * self.nb_subdom
 
-        for i in range(self.nb_subdom):
-            centers_tmp = self.centers[:, self.ind_v_subdom[i]]
-            U_tmp = U[self.ind_v_subdom[i]]
-            U_smooth[i], U_tps[i], summaries[i] = self.compute_tps_coeff_iter(
-                centers_tmp, U_tmp
+        for idx in range(self.nb_subdom):
+            centers_domain = self.centers[:, self.indices_domains[idx]]
+            values_domain = values[self.indices_domains[idx]]
+            (
+                smoothed_field_domains[idx],
+                weights_domains[idx],
+                summaries[idx],
+            ) = self.compute_tps_weights_iter(
+                centers_domain, values_domain, self.smoothing_coef
             )
 
-        U_smooth_tmp = np.zeros(self.centers[1].shape)
-        nb_tps = np.zeros(self.centers[1].shape, dtype=int)
+        smoothed_field = np.zeros(self.num_centers)
+        nb_tps = np.zeros(self.num_centers, dtype=int)
         summary = {"nb_fixed_vectors": [], "max(Udiff)": [], "nb_iterations": []}
 
-        for i in range(self.nb_subdom):
-            U_smooth_tmp[self.ind_v_subdom[i]] += U_smooth[i]
-            nb_tps[self.ind_v_subdom[i]] += 1
+        for idx in range(self.nb_subdom):
+            smoothed_field[self.indices_domains[idx]] += smoothed_field_domains[
+                idx
+            ]
+            nb_tps[self.indices_domains[idx]] += 1
             for key in ("nb_fixed_vectors", "max(Udiff)", "nb_iterations"):
-                summary[key].append(summaries[i][key])
+                summary[key].append(summaries[idx][key])
 
         summary["nb_fixed_vectors_tot"] = sum(summary["nb_fixed_vectors"])
-        U_smooth_tmp /= nb_tps
+        smoothed_field /= nb_tps
 
-        return U_smooth_tmp, U_tps, summary
+        return smoothed_field, weights_domains, summary
 
     def init_with_new_positions(self, new_positions):
         """Initialize with the new positions
@@ -142,15 +139,14 @@
           new_positions[1] and new_positions[0] correspond to the x and y values, respectively.
 
         """
-        self.new_positions = new_positions
         xs = new_positions[1]
         ys = new_positions[0]
+        self.num_new_positions = xs.size
 
-        ind_new_positions_subdom = []
-
+        self.ind_new_positions_domains = ind_new_positions_domains = []
         for iy in range(self.nb_subdomy):
             for ix in range(self.nb_subdomx):
-                ind_new_positions_subdom.append(
+                ind_new_positions_domains.append(
                     np.where(
                         (xs >= self.xmin_limits[ix])
                         & (xs < self.xmax_limits[ix])
@@ -159,56 +155,52 @@
                     )[0]
                 )
 
-        self.ind_new_positions_subdom = ind_new_positions_subdom
-
-        self.norm_coefs = np.zeros(self.new_positions.shape[1])
-        for i_subdom in range(self.nb_subdom):
+        self.norm_coefs = np.zeros(self.num_new_positions)
+        for i_domain in range(self.nb_subdom):
             # TODO: replace 1 by an appropriate function of
-            # ind_new_positions_subdom[i_subdom]
-            # + save another list of 1d array like ind_new_positions_subdom
+            # ind_new_positions_domains[i_subdom]
+            # + save another list of 1d array like ind_new_positions_domains
             # containing the norm coefficients for each subdomain
-            self.norm_coefs[ind_new_positions_subdom[i_subdom]] += 1
-
-        EM = [None] * self.nb_subdom
+            self.norm_coefs[ind_new_positions_domains[i_domain]] += 1
 
-        for i in range(self.nb_subdom):
-            centers_tmp = self.centers[:, self.ind_v_subdom[i]]
-            new_positions_tmp = new_positions[:, ind_new_positions_subdom[i]]
-            EM[i] = compute_tps_matrix(new_positions_tmp, centers_tmp)
-
-        self.EM = EM
-
-    def interpolate(self, U_tps):
-        U_eval = np.zeros(self.new_positions.shape[1])
-
-        for i in range(self.nb_subdom):
-            U_eval[self.ind_new_positions_subdom[i]] += np.dot(
-                U_tps[i], self.EM[i]
+        self.tps_matrices = [None] * self.nb_subdom
+        for i_domain in range(self.nb_subdom):
+            centers_tmp = self.centers[:, self.indices_domains[i_domain]]
+            new_positions_tmp = new_positions[
+                :, ind_new_positions_domains[i_domain]
+            ]
+            self.tps_matrices[i_domain] = compute_tps_matrix(
+                new_positions_tmp, centers_tmp
             )
 
-        U_eval /= self.norm_coefs
+    def interpolate(self, tps_weights_domains):
+        """Interpolate on new positions"""
+        values = np.zeros(self.num_new_positions)
+        for i_domain in range(self.nb_subdom):
+            values[self.ind_new_positions_domains[i_domain]] += np.dot(
+                tps_weights_domains[i_domain], self.tps_matrices[i_domain]
+            )
+        values /= self.norm_coefs
+        return values
 
-        return U_eval
-
-    def compute_tps_coeff_iter(self, centers, values: Array[np.float64, "1d"]):
+    def compute_tps_weights_iter(
+        self, centers, values: Array[np.float64, "1d"], smoothing_coef
+    ):
         """Compute the thin plate spline (tps) coefficients removing erratic
         vectors
 
-        It computes iteratively "compute_tps_coeff", compares the tps
+        It computes iteratively "compute_tps_weights", compares the tps
         result to the initial data and remove it if difference is
         larger than the given threshold
 
         """
         summary = {"nb_fixed_vectors": 0}
-
-        # normalization as UVmat so that the effect of the filter do not depends
-        # too much on the size of the domains
-        smoothing_coef = self.smoothing_coef * values.size / 1000
-
-        U_smooth, U_tps = compute_tps_coeff(centers, values, smoothing_coef)
+        smoothed_values, tps_weights = compute_tps_weights(
+            centers, values, smoothing_coef
+        )
         count = 1
         if self.threshold is not None:
-            differences = np.sqrt((U_smooth - values) ** 2)
+            differences = np.sqrt((smoothed_values - values) ** 2)
             ind_erratic_vector = np.argwhere(differences > self.threshold)
 
             summary["max(Udiff)"] = max(differences)
@@ -216,12 +208,12 @@
             nb_fixed_vectors = 0
             while ind_erratic_vector.size != 0:
                 nb_fixed_vectors += ind_erratic_vector.size
-                values[ind_erratic_vector] = U_smooth[ind_erratic_vector]
-                U_smooth, U_tps = compute_tps_coeff(
+                values[ind_erratic_vector] = smoothed_values[ind_erratic_vector]
+                smoothed_values, tps_weights = compute_tps_weights(
                     centers, values, smoothing_coef
                 )
 
-                differences = np.sqrt((U_smooth - values) ** 2)
+                differences = np.sqrt((smoothed_values - values) ** 2)
                 ind_erratic_vector = np.argwhere(differences > self.threshold)
                 count += 1
 
@@ -235,4 +227,4 @@
 
             summary["nb_fixed_vectors"] = nb_fixed_vectors
         summary["nb_iterations"] = count
-        return U_smooth, U_tps, summary
+        return smoothed_values, tps_weights, summary
diff --git a/src/fluidimage/works/piv/singlepass.py b/src/fluidimage/works/piv/singlepass.py
--- a/src/fluidimage/works/piv/singlepass.py
+++ b/src/fluidimage/works/piv/singlepass.py
@@ -490,12 +490,12 @@
                     deltaxs_smooth,
                     deltaxs_tps,
                     summary,
-                ) = tps.compute_tps_coeff_subdom(deltaxs)
+                ) = tps.compute_tps_weights_subdom(deltaxs)
                 (
                     deltays_smooth,
                     deltays_tps,
                     summary,
-                ) = tps.compute_tps_coeff_subdom(deltays)
+                ) = tps.compute_tps_weights_subdom(deltays)
             except np.linalg.LinAlgError:
                 print("LinAlgError in TPS: compute delta_approx with griddata")
                 deltaxs_approx = griddata(
# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1714141700 -7200
#      Fri Apr 26 16:28:20 2024 +0200
# Node ID e365c409ebdbf7607596470bf8b82fac8dfde28a
# Parent  8dca5014bda76cd23bd656815219ba9eaef83977
TPS with norm coefs

diff --git a/doc/examples/profile_piv_work.py b/doc/examples/profile_piv_work.py
--- a/doc/examples/profile_piv_work.py
+++ b/doc/examples/profile_piv_work.py
@@ -38,7 +38,7 @@
 params.multipass.use_tps = False
 params.multipass.subdom_size = 200
 params.multipass.smoothing_coef = 10.0
-params.multipass.threshold_tps = 0.5
+params.multipass.threshold_tps = 1.5
 
 params.fix.correl_min = 0.15
 params.fix.threshold_diff_neighbour = 3
diff --git a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
--- a/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
+++ b/src/fluidimage/calcul/interpolate/thin_plate_spline_subdom.py
@@ -24,11 +24,15 @@
 class ThinPlateSplineSubdom:
     """Helper class for thin plate interpolation."""
 
-    ind_new_positions_domains: List["np.int64[:]"]
-    norm_coefs: "float[:]"
-    num_new_positions: int
     num_centers: int
     tps_matrices: List["float[:,:]"]
+    norm_coefs: "float[:]"
+    norm_coefs_domains: List["float[:]"]
+
+    num_new_positions: int
+    ind_new_positions_domains: List["np.int64[:]"]
+    norm_coefs_new_pos: "float[:]"
+    norm_coefs_new_pos_domains: List["float[:]"]
 
     def __init__(
         self,
@@ -77,24 +81,48 @@
         buffer_length_x = coef_buffer * range_x / nb_subdomx
         buffer_length_y = coef_buffer * range_y / nb_subdomy
 
-        self.xmin_limits = x_dom[:-1] - buffer_length_x
-        self.xmax_limits = x_dom[1:] + buffer_length_x
+        self.limits_min_x = x_dom[:-1] - buffer_length_x
+        self.limits_max_x = x_dom[1:] + buffer_length_x
+
+        self.limits_min_y = y_dom[:-1] - buffer_length_y
+        self.limits_max_y = y_dom[1:] + buffer_length_y
 
-        self.ymin_limits = y_dom[:-1] - buffer_length_y
-        self.ymax_limits = y_dom[1:] + buffer_length_y
+        self.xmax_domains = np.empty(self.nb_subdom)
+        self.ymax_domains = np.empty(self.nb_subdom)
+        self.xc_domains = np.empty(self.nb_subdom)
+        self.yc_domains = np.empty(self.nb_subdom)
+        self.indices_domains = []
 
-        self.indices_domains = []
+        i_dom = 0
         for iy in range(nb_subdomy):
             for ix in range(nb_subdomx):
+                xmin = self.limits_min_x[ix]
+                xmax = self.limits_max_x[ix]
+                ymin = self.limits_min_y[iy]
+                ymax = self.limits_max_y[iy]
+
                 self.indices_domains.append(
                     np.where(
-                        (xs >= self.xmin_limits[ix])
-                        & (xs < self.xmax_limits[ix])
-                        & (ys >= self.ymin_limits[iy])
-                        & (ys < self.ymin_limits[iy])
+                        (xs >= xmin) & (xs < xmax) & (ys >= ymin) & (ys < ymax)
                     )[0]
                 )
 
+                self.xmax_domains[i_dom] = xmax
+                self.ymax_domains[i_dom] = ymax
+                self.xc_domains[i_dom] = 0.5 * (xmin + xmax)
+                self.yc_domains[i_dom] = 0.5 * (ymin + ymax)
+                i_dom += 1
+
+        self.norm_coefs = np.zeros(self.num_centers)
+        self.norm_coefs_domains = []
+        for i_dom in range(self.nb_subdom):
+            indices_domain = self.indices_domains[i_dom]
+            xs_domain = xs[indices_domain]
+            ys_domain = ys[indices_domain]
+            coefs = self._compute_coef_norm(xs_domain, ys_domain, i_dom)
+            self.norm_coefs_domains.append(coefs)
+            self.norm_coefs[indices_domain] += coefs
+
     def compute_tps_weights_subdom(self, values):
         """Compute the TPS weights for all subdomains"""
         smoothed_field_domains = [None] * self.nb_subdom
@@ -113,19 +141,18 @@
             )
 
         smoothed_field = np.zeros(self.num_centers)
-        nb_tps = np.zeros(self.num_centers, dtype=int)
         summary = {"nb_fixed_vectors": [], "max(Udiff)": [], "nb_iterations": []}
 
         for idx in range(self.nb_subdom):
-            smoothed_field[self.indices_domains[idx]] += smoothed_field_domains[
-                idx
-            ]
-            nb_tps[self.indices_domains[idx]] += 1
+            indices_domain = self.indices_domains[idx]
+            smoothed_field[indices_domain] += (
+                self.norm_coefs_domains[idx] * smoothed_field_domains[idx]
+            )
             for key in ("nb_fixed_vectors", "max(Udiff)", "nb_iterations"):
                 summary[key].append(summaries[idx][key])
 
         summary["nb_fixed_vectors_tot"] = sum(summary["nb_fixed_vectors"])
-        smoothed_field /= nb_tps
+        smoothed_field /= self.norm_coefs
 
         return smoothed_field, weights_domains, summary
 
@@ -148,20 +175,23 @@
             for ix in range(self.nb_subdomx):
                 ind_new_positions_domains.append(
                     np.where(
-                        (xs >= self.xmin_limits[ix])
-                        & (xs < self.xmax_limits[ix])
-                        & (ys >= self.ymin_limits[iy])
-                        & (ys < self.ymin_limits[iy])
+                        (xs >= self.limits_min_x[ix])
+                        & (xs < self.limits_max_x[ix])
+                        & (ys >= self.limits_min_y[iy])
+                        & (ys < self.limits_max_y[iy])
                     )[0]
                 )
 
-        self.norm_coefs = np.zeros(self.num_new_positions)
+        self.norm_coefs_new_pos = np.zeros(self.num_new_positions)
+        self.norm_coefs_new_pos_domains = []
+
         for i_domain in range(self.nb_subdom):
-            # TODO: replace 1 by an appropriate function of
-            # ind_new_positions_domains[i_subdom]
-            # + save another list of 1d array like ind_new_positions_domains
-            # containing the norm coefficients for each subdomain
-            self.norm_coefs[ind_new_positions_domains[i_domain]] += 1
+            indices_domain = ind_new_positions_domains[i_domain]
+            xs_domain = xs[indices_domain]
+            ys_domain = ys[indices_domain]
+            coefs = self._compute_coef_norm(xs_domain, ys_domain, i_domain)
+            self.norm_coefs_new_pos_domains.append(coefs)
+            self.norm_coefs_new_pos[indices_domain] += coefs
 
         self.tps_matrices = [None] * self.nb_subdom
         for i_domain in range(self.nb_subdom):
@@ -173,14 +203,32 @@
                 new_positions_tmp, centers_tmp
             )
 
-    def interpolate(self, tps_weights_domains):
+    def _compute_coef_norm(self, xs_domain, ys_domain, i_domain):
+
+        x_center_domain = self.xc_domains[i_domain]
+        y_center_domain = self.yc_domains[i_domain]
+
+        x_max_domain = self.xmax_domains[i_domain]
+        y_max_domain = self.ymax_domains[i_domain]
+
+        dx_max = x_max_domain - x_center_domain
+        dy_max = y_max_domain - y_center_domain
+
+        dx2_normed = (xs_domain - x_center_domain) ** 2 / dx_max**2
+        dy2_normed = (ys_domain - y_center_domain) ** 2 / dy_max**2
+
+        return (1 - dx2_normed) * (1 - dy2_normed) + 0.001
+
+    def interpolate(self, weights_domains):
         """Interpolate on new positions"""
         values = np.zeros(self.num_new_positions)
         for i_domain in range(self.nb_subdom):
-            values[self.ind_new_positions_domains[i_domain]] += np.dot(
-                tps_weights_domains[i_domain], self.tps_matrices[i_domain]
+            values[
+                self.ind_new_positions_domains[i_domain]
+            ] += self.norm_coefs_new_pos_domains[i_domain] * np.dot(
+                weights_domains[i_domain], self.tps_matrices[i_domain]
             )
-        values /= self.norm_coefs
+        values /= self.norm_coefs_new_pos
         return values
 
     def compute_tps_weights_iter(
diff --git a/src/fluidimage/works/piv/test_piv.py b/src/fluidimage/works/piv/test_piv.py
--- a/src/fluidimage/works/piv/test_piv.py
+++ b/src/fluidimage/works/piv/test_piv.py
@@ -23,6 +23,7 @@
     params.piv0.nb_peaks_to_search = 2
 
     params.multipass.number = 2
+    params.multipass.use_tps = True
 
     params.fix.displacement_max = 3
     params.fix.threshold_diff_neighbour = 2
# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1714143087 -7200
#      Fri Apr 26 16:51:27 2024 +0200
# Node ID 967222aa08b0b059a9ef9fefbaf6cc12bf2ce163
# Parent  e365c409ebdbf7607596470bf8b82fac8dfde28a
Save deltaxs_smooth, deltays_smooth when available

diff --git a/src/fluidimage/data_objects/piv.py b/src/fluidimage/data_objects/piv.py
--- a/src/fluidimage/data_objects/piv.py
+++ b/src/fluidimage/data_objects/piv.py
@@ -246,6 +246,10 @@
       Equivalent to the `_approx` variables but for the last pass
       and of size ``num_vectors``.
 
+    deltaxs_smooth, deltays_smooth:
+
+      Smoothed displacements at the positions ``xs``, ``ys`` (present only for
+      TPS interpolation).
     """
 
     _keys_to_be_saved = [
@@ -260,6 +264,8 @@
         "iyvecs_approx",
         "deltaxs_final",
         "deltays_final",
+        "deltaxs_smooth",
+        "deltays_smooth",
         "ixvecs_final",
         "iyvecs_final",
     ]
# HG changeset patch
# User paugier <pierre.augier@univ-grenoble-alpes.fr>
# Date 1714144918 -7200
#      Fri Apr 26 17:21:58 2024 +0200
# Node ID 65002c3eafc3ed29572e0b20cd60c3e5b10681cf
# Parent  967222aa08b0b059a9ef9fefbaf6cc12bf2ce163
threshold_tps, change default to 1.5 (like UVmat)

diff --git a/src/fluidimage/works/piv/multipass.py b/src/fluidimage/works/piv/multipass.py
--- a/src/fluidimage/works/piv/multipass.py
+++ b/src/fluidimage/works/piv/multipass.py
@@ -63,7 +63,7 @@
                 "use_tps": "last",
                 "subdom_size": 200,
                 "smoothing_coef": 2.0,
-                "threshold_tps": 1.0,
+                "threshold_tps": 1.5,
             },
         )
 
@@ -93,7 +93,7 @@
 - smoothing_coef : float
 
   Coefficient used for the TPS method. The result is smoother for larger
-  smoothing_coef.
+  smoothing_coef. 2 is often reasonable. Can typically be between 0 to 40.
 
 - threshold_tps :  float