# 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