Merge branch 'idw_interpolation'

Author: alexlib

bump_version.py

@@ -53,16 +53,31 @@ def bump_version(version_type):
     new_version = f"{major}.{minor}.{patch}"
     print(f"New version: {new_version}")
 
-    # Replace the version in the file
+    # Replace the version in __init__.py
     new_content = re.sub(
         r"__version__\s*=\s*['\"]([^'\"]*)['\"]",
         f"__version__ = '{new_version}'",
         content
     )
-    
-    # Write the updated content back to the file
     init_file.write_text(new_content)
-    print(f"Version bumped to {new_version}")
+    print(f"Updated flowtracks/__init__.py to version {new_version}")
+
+    # Also update version in pyproject.toml
+    pyproject_file = Path('pyproject.toml')
+    if pyproject_file.exists():
+        pyproject_content = pyproject_file.read_text()
+        if '[project]' in pyproject_content:
+            pyproject_new = re.sub(
+                r'version\s*=\s*["\"][^"\"]*["\"]',
+                f'version = "{new_version}"',
+                pyproject_content
+            )
+            pyproject_file.write_text(pyproject_new)
+            print(f"Updated pyproject.toml to version {new_version}")
+        else:
+            print("Warning: No [project] table found in pyproject.toml, version not updated there.")
+    else:
+        print("Warning: pyproject.toml not found, version not updated there.")
 
 if __name__ == "__main__":
     if len(sys.argv) != 2:

flowtracks/__init__.py

@@ -1,3 +1,3 @@
 # -*- coding: utf-8 -*-
 
-__version__ = '1.1.0'
+__version__ = '1.1.1'

flowtracks/future_idea_xarray_dask_zarr.py

@@ -0,0 +1,176 @@
+#!/usr/bin/env python3
+"""
+Boilerplate for 3D-PTV trajectory post-processing
+using xarray + dask + zarr
+
+Features:
+- Ragged array encoding for trajectories of different lengths
+- Vector-style storage (obs, component)
+- Spline smoothing and derivative calculation (velocity, acceleration)
+- Resampling onto a uniform time base
+- Streaming mode: append processed trajectories to a Zarr store
+"""
+
+import numpy as np
+import xarray as xr
+from scipy.interpolate import UnivariateSpline
+import zarr
+
+# ---------------------------------------------------------------------
+# 1. Example: Ragged array encoding for variable-length trajectories
+# ---------------------------------------------------------------------
+
+def build_ragged_example():
+    traj1_t = np.array([0, 1, 2])
+    traj1_x = np.array([0, 1, 2])
+    traj1_y = np.array([0, 1, 4])
+    traj1_z = np.array([0, 0, 0])
+
+    traj2_t = np.array([0, 2, 4, 6, 8])
+    traj2_x = np.array([0, 2, 4, 6, 8])
+    traj2_y = np.array([0, -1, -2, -3, -4])
+    traj2_z = np.array([0, 1, 0, -1, 0])
+
+    # Concatenate into ragged structure
+    times = np.concatenate([traj1_t, traj2_t])
+    positions = np.vstack([
+        np.stack([traj1_x, traj1_y, traj1_z], axis=-1),
+        np.stack([traj2_x, traj2_y, traj2_z], axis=-1)
+    ])
+    trajectory_id = np.concatenate([
+        np.full(traj1_t.shape, 0),
+        np.full(traj2_t.shape, 1)
+    ])
+
+    ds = xr.Dataset(
+        {
+            "t": ("obs", times),
+            "pos": (("obs", "component"), positions),
+            "trajectory": ("obs", trajectory_id),
+        },
+        coords={"component": ["x", "y", "z"], "obs": np.arange(len(times))}
+    )
+
+    return ds
+
+
+# ---------------------------------------------------------------------
+# 2. Compute derivatives (velocity, acceleration) in ragged array
+# ---------------------------------------------------------------------
+
+def compute_derivatives_ragged(ds):
+    def _derivs(sub):
+        dt = np.gradient(sub.t.values)
+        dpos = np.gradient(sub.pos.values, axis=0)
+        vel = dpos / dt[:, None]
+        acc = np.gradient(vel, axis=0) / dt[:, None]
+        return xr.Dataset({
+            "vel": (("obs", "component"), vel),
+            "acc": (("obs", "component"), acc)
+        })
+
+    derivs = ds.groupby("trajectory").map(_derivs)
+    return xr.merge([ds, derivs])
+
+
+# ---------------------------------------------------------------------
+# 3. Spline smoothing + resampling on uniform time base
+# ---------------------------------------------------------------------
+
+def smooth_and_resample(t, pos, t_uniform, s=0.0):
+    """Smooth trajectory with spline, resample to uniform time base.
+    Returns position, velocity, acceleration arrays of shape (len(t_uniform), 3).
+    """
+    comps = []
+    vels = []
+    accs = []
+    for d in range(pos.shape[1]):  # loop over x,y,z
+        spline = UnivariateSpline(t, pos[:, d], s=s)
+        p = spline(t_uniform)
+        v = spline.derivative(1)(t_uniform)
+        a = spline.derivative(2)(t_uniform)
+        comps.append(p)
+        vels.append(v)
+        accs.append(a)
+
+    pos_u = np.stack(comps, axis=-1)
+    vel_u = np.stack(vels, axis=-1)
+    acc_u = np.stack(accs, axis=-1)
+    return pos_u, vel_u, acc_u
+
+
+# ---------------------------------------------------------------------
+# 4. Streaming mode: append processed trajectories to Zarr
+# ---------------------------------------------------------------------
+
+def init_zarr_store(store_path, t_uniform):
+    components = ["x", "y", "z"]
+
+    ds = xr.Dataset(
+        data_vars={
+            "position": (("trajectory", "time", "component"),
+                         np.empty((0, len(t_uniform), len(components)))),
+            "velocity": (("trajectory", "time", "component"),
+                         np.empty((0, len(t_uniform), len(components)))),
+            "acceleration": (("trajectory", "time", "component"),
+                             np.empty((0, len(t_uniform), len(components))))
+        },
+        coords={
+            "time": t_uniform,
+            "component": components,
+            "trajectory": []
+        }
+    )
+    ds.to_zarr(store_path, mode="w")
+
+
+def append_to_zarr(store_path, traj_id, t, pos, t_uniform, s=0.0):
+    pos_u, vel_u, acc_u = smooth_and_resample(t, pos, t_uniform, s=s)
+
+    new = xr.Dataset(
+        {
+            "position": (("trajectory", "time", "component"), pos_u[np.newaxis, ...]),
+            "velocity": (("trajectory", "time", "component"), vel_u[np.newaxis, ...]),
+            "acceleration": (("trajectory", "time", "component"), acc_u[np.newaxis, ...]),
+        },
+        coords={
+            "trajectory": [traj_id],
+            "time": t_uniform,
+            "component": ["x", "y", "z"]
+        }
+    )
+    new.to_zarr(store_path, mode="a", append_dim="trajectory")
+
+
+# ---------------------------------------------------------------------
+# 5. Example usage
+# ---------------------------------------------------------------------
+
+if __name__ == "__main__":
+    # Step 1: Build ragged example
+    ds_ragged = build_ragged_example()
+    print("Ragged dataset:")
+    print(ds_ragged)
+
+    # Step 2: Compute derivatives in ragged array
+    ds_with_derivs = compute_derivatives_ragged(ds_ragged)
+    print("\nRagged dataset with velocity and acceleration:")
+    print(ds_with_derivs)
+
+    # Step 3+4: Streaming to Zarr
+    t_uniform = np.linspace(0, 8, 81)  # uniform time base (0.1s step)
+    store = "trajectories.zarr"
+    init_zarr_store(store, t_uniform)
+
+    # Add first trajectory
+    obs0 = ds_ragged.where(ds_ragged.trajectory == 0, drop=True)
+    append_to_zarr(store, traj_id=0, t=obs0.t.values, pos=obs0.pos.values, t_uniform=t_uniform, s=0.1)
+
+    # Add second trajectory
+    obs1 = ds_ragged.where(ds_ragged.trajectory == 1, drop=True)
+    append_to_zarr(store, traj_id=1, t=obs1.t.values, pos=obs1.pos.values, t_uniform=t_uniform, s=0.1)
+
+    # Open the final store lazily with xarray+dask
+    ds_zarr = xr.open_zarr(store)
+    print("\nZarr dataset (streamed trajectories):")
+    print(ds_zarr)

flowtracks/interpolation.py

@@ -640,13 +640,17 @@ def weights(self, dists, use_parts, unused_marker=None):
             interpolation point 1..m.
 
         Returns:
-        weights - an (m,n) array.
+        weights - an (m,k) array.
         """
-        if unused_marker is None:
-            unused_marker = self.field_positions().shape[0]
-
         weights = dists**-self._par
-        weights[use_parts == unused_marker] = 0.
+        # If use_parts is boolean, just mask out non-neighbors
+        if use_parts.dtype == bool:
+            weights[~use_parts] = 0.
+        else:
+            # If use_parts is indices, mask out invalid indices
+            if unused_marker is None:
+                unused_marker = dists.shape[1]
+            weights[use_parts == unused_marker] = 0.
         return weights
 
     def set_scene(self, tracer_pos, interp_points,
@@ -712,7 +716,42 @@ def __call__(self, tracer_pos, interp_points, data, companionship=None):
         vel_interp - an (m,3) array with the interpolated value at the position
             of each particle, [m/s].
         """
-        raise NotImplementedError()
+        if len(tracer_pos) == 0:
+            warnings.warn("No tracers in frame, interpolation returned zeros.")
+            ret_shape = data.shape[-1] if data.ndim > 1 else 1
+            return np.zeros((interp_points.shape[0], ret_shape))
+
+        dists, use_parts = select_neighbs(tracer_pos, interp_points,
+                                          self._radius, self._neighbs,
+                                          companionship)
+
+        m, n = dists.shape
+        if data.ndim == 1:
+            data = data[:, None]
+        matched_data = np.zeros((m, n, data.shape[1]))
+        for i in range(m):
+            matched_data[i] = data
+
+        # Handle exact matches: if any distance is zero, set result to that data value
+        exact_match = (dists == 0)
+        has_exact = exact_match.any(axis=1)
+        vel_interp = np.empty((m, data.shape[1]), dtype=data.dtype)
+        weights = self.weights(dists, use_parts)
+        for i in range(m):
+            if has_exact[i]:
+                # Use the first exact match (should only be one)
+                idx = np.where(exact_match[i])[0][0]
+                vel_interp[i] = matched_data[i, idx]
+            else:
+                sum_weights = weights[i].sum()
+                if sum_weights == 0:
+                    vel_interp[i] = 0
+                else:
+                    vel_interp[i] = (weights[i][:, None] * matched_data[i]).sum(axis=0) / sum_weights
+        # Always return a 2D array (m, d) even for single-point or 1D data
+        if vel_interp.ndim == 1:
+            vel_interp = vel_interp[:, None]
+        return vel_interp
 
     def _meth_interp(self, act_neighbs=None):
         """

pyproject.toml

@@ -1,3 +1,33 @@
 [build-system]
 requires = ["setuptools>=42", "wheel"]
 build-backend = "setuptools.build_meta"
+
+[project]
+name = "flowtracks"
+version = "1.1.1"
+description = "Library for handling of PTV trajectory database."
+readme = "README.md"
+requires-python = ">=3.9"
+authors = [{ name = "Yosef Meller", email = "yosefm@gmail.com" }]
+license = { file = "LICENSE.txt" }
+dependencies = [
+    "numpy",
+    "scipy",
+    "tables"
+]
+classifiers = [
+    "Development Status :: 5 - Production/Stable",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: BSD License",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.8",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
+    "Topic :: Scientific/Engineering"
+]
+
+[dependency-groups]
+dev = [
+    "pytest>=8.4.2",
+]

tests/test_idw_call.py

@@ -0,0 +1,63 @@
+import unittest
+import numpy as np
+from flowtracks import interpolation
+
+class TestIDWCall(unittest.TestCase):
+    def setUp(self):
+        # Simple 2D case for easy checking
+        self.tracer_pos = np.array([
+            [0.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [1.0, 1.0, 0.0],
+        ])
+        self.data = np.array([
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0],
+            [1.0, 1.0, 1.0],
+        ])
+        self.interp_points = np.array([
+            [0.5, 0.5, 0.0],
+        ])
+
+    def test_basic_idw(self):
+        idw = interpolation.InverseDistanceWeighter(num_neighbs=4, param=1)
+        result = idw(self.tracer_pos, self.interp_points, self.data)
+        self.assertEqual(result.shape, (1, 3))
+        # Should be a weighted average, check sum is reasonable
+        self.assertTrue(np.all(result >= 0))
+        self.assertTrue(np.all(result <= 1))
+
+    def test_empty_tracers(self):
+        import warnings
+        idw = interpolation.InverseDistanceWeighter(num_neighbs=2, param=1)
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter("always")
+            result = idw(np.empty((0, 3)), self.interp_points, np.empty((0, 3)))
+            self.assertTrue(any(issubclass(warn.category, UserWarning) for warn in w))
+        self.assertEqual(result.shape, (1, 3))
+        self.assertTrue(np.all(result == 0))
+
+    def test_1d_data(self):
+        data_1d = np.array([1.0, 2.0, 3.0, 4.0])
+        idw = interpolation.InverseDistanceWeighter(num_neighbs=4, param=1)
+        result = idw(self.tracer_pos, self.interp_points, data_1d)
+        self.assertEqual(result.shape, (1, 1))
+        self.assertTrue(np.all(result >= 1))
+        self.assertTrue(np.all(result <= 4))
+
+    def test_companionship(self):
+        # Exclude the first tracer from the interpolation
+        import warnings
+        companions = np.array([0])
+        idw = interpolation.InverseDistanceWeighter(num_neighbs=4, param=1)
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter("always")
+            result = idw(self.tracer_pos, self.interp_points, self.data, companionship=companions)
+            self.assertTrue(any(issubclass(warn.category, RuntimeWarning) for warn in w))
+        self.assertEqual(result.shape, (1, 3))
+        # The result should not be equal to the average including the first tracer
+        idw_all = interpolation.InverseDistanceWeighter(num_neighbs=4, param=1)
+        result_all = idw_all(self.tracer_pos, self.interp_points, self.data)
+        self.assertFalse(np.allclose(result, result_all))