Additions to interp_to_regular_grid!

Added assertion that the determined dim has no non-zero entries!
Added distance_upper_bound param to customize non-convex interpolations!

src/empyre/utils/misc.py

@@ -24,7 +24,7 @@ def levi_civita(i, j, k):
     return (j-i) * (k-i) * (k-j) / (np.abs(j-i) * np.abs(k-i) * np.abs(k-j))
 
 
-def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex=True):
+def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex=True, distance_upper_bound=None):
     """Interpolate values on points to regular grid.
 
     Parameters
@@ -45,6 +45,11 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     convex : bool, optional
         By default True. If this is set to False, additional measures are taken to find holes in the point cloud.
         WARNING: this is an experimental feature that should be used with caution!
+    distance_upper_bound: float, optional
+        Only used if `convex=False`. Set the upper bound, determining if a point of the new (interpolated) grid is too
+        far away from any original point. They are assumed to be in a "hole" and their values are set to zero. Set this
+        value in nm, it will be converted to the local unit of the original points internally.
+        If not set and `convex=True`, double of the the mean of `scale` is assumed.
 
     Returns
     -------
@@ -67,12 +72,12 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
     _log.info(f'z-range: {z_min:.2g} <-> {z_max:.2g} ({z_diff:.2g})')
     # Determine dimensions from given grid spacing a:
     dim = tuple(np.round(np.asarray((z_diff/scale[0], y_diff/scale[1], x_diff/scale[2]), dtype=int)))
-    x = x_min + scale[2] * (np.arange(dim[2]) + 0.5)  # +0.5: shift to pixel center!
-    y = y_min + scale[1] * (np.arange(dim[1]) + 0.5)  # +0.5: shift to pixel center!
+    assert all(dim) > 0, f'All dimensions of dim={dim} need to be > 0, please adjust the scale accordingly!'
     z = z_min + scale[0] * (np.arange(dim[0]) + 0.5)  # +0.5: shift to pixel center!
+    y = y_min + scale[1] * (np.arange(dim[1]) + 0.5)  # +0.5: shift to pixel center!
+    x = x_min + scale[2] * (np.arange(dim[2]) + 0.5)  # +0.5: shift to pixel center!
     # Create points for new Euclidian grid; fliplr for (x, y, z) order:
-    test = np.asarray(list(itertools.product(z, y, x)))
-    points_euc = np.fliplr(test)
+    points_euc = np.fliplr(np.asarray(list(itertools.product(z, y, x))))
     # Make values 2D (if not already); double .T so that a new axis is added at the END (n, 1):
     values = np.atleast_2d(values.T).T
     # Prepare interpolated grid:
@@ -96,7 +101,11 @@ def interp_to_regular_grid(points, values, scale, scale_factor=1, step=1, convex
         tree = cKDTree(points)
         # Query the tree for nearest neighbors, x: points to query, k: number of neighbors, p: norm
         # to use (here: 2 - Euclidean), distance_upper_bound: maximum distance that is searched!
-        data, leafsize = tree.query(x=points_euc, k=1, p=2, distance_upper_bound=2*np.mean(scale))
+        if distance_upper_bound is None:  # Take the mean of the scale for the upper bound:
+            distance_upper_bound = 2 * np.mean(scale)  # NOTE: could be problematic for wildly varying scale numbers.
+        else:  # Convert to local scale:
+            distance_upper_bound *= scale_factor
+        data, leafsize = tree.query(x=points_euc, k=1, p=2, distance_upper_bound=distance_upper_bound)
         # Create boolean mask that determines which interpolation points have no neighbor near enough:
         mask = np.isinf(data).reshape(dim)  # Points further away than upper bound were marked 'inf'!
         for i in tqdm(range(values.shape[-1])):  # Set these points to zero (NOTE: This can take a looooong time):