From 05a15030c95c4237313897c3ecf8c6de6e466059 Mon Sep 17 00:00:00 2001 From: Pierre-antoine Comby Date: Fri, 2 Feb 2024 10:17:56 +0100 Subject: [PATCH] fix: robust voronoi estimation. --- src/mrinufft/density/geometry_based.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/mrinufft/density/geometry_based.py b/src/mrinufft/density/geometry_based.py index 238bad95..4dc0ecc5 100644 --- a/src/mrinufft/density/geometry_based.py +++ b/src/mrinufft/density/geometry_based.py @@ -77,14 +77,22 @@ def voronoi_unique(traj, *args, **kwargs): idx_vertices = v.regions[v.point_region[mm]] if np.all([i != -1 for i in idx_vertices]): wi[mm] = vol(v.vertices[idx_vertices]) + else: + wi[mm] = np.inf + # some voronoi cell are considered closed, but have a too big area. + # (They are closing near infinity). + # we classify them as open cells as well. + outlier_thresh = np.percentile(wi, 95) + wi[wi > outlier_thresh] = np.inf + # For edge point (infinite voronoi cells) we extrapolate from neighbours # Initial implementation in Jeff Fessler's MIRT rho = np.sum(traj**2, axis=1) - igood = rho > 0.6 * np.max(rho) + igood = (rho > 0.6 * np.max(rho)) & ~np.isinf(wi) if len(igood) < 10: print("dubious extrapolation with", len(igood), "points") poly = np.polynomial.Polynomial.fit(rho[igood], wi[igood], 3) - wi[wi == 0] = poly(rho[wi == 0]) + wi[np.isinf(wi)] = poly(rho[np.isinf(wi)]) return wi