From d00331a5c364412b9aa48714009c3284e55ff8c1 Mon Sep 17 00:00:00 2001 From: ezwelty Date: Fri, 16 Feb 2024 16:39:18 +0100 Subject: [PATCH] Add surface-tangent motion model (Brinkerhoff 2017) Includes both cartesian and cylindrical versions. --- src/glimpse/track/__init__.py | 17 ++- src/glimpse/track/motion.py | 211 ++++++++++++++++++++++++++++++++++ 2 files changed, 226 insertions(+), 2 deletions(-) diff --git a/src/glimpse/track/__init__.py b/src/glimpse/track/__init__.py index abd23b40..9f5326c1 100644 --- a/src/glimpse/track/__init__.py +++ b/src/glimpse/track/__init__.py @@ -1,8 +1,21 @@ """Track image features through time using particle filtering.""" -from .motion import CartesianMotion, CylindricalMotion +from .motion import ( + CartesianMotion, + CylindricalMotion, + TangentCartesianMotion, + TangentCylindricalMotion, +) from .observer import Observer from .tracker import Tracker from .tracks import Tracks -__all__ = ["CartesianMotion", "CylindricalMotion", "Observer", "Tracker", "Tracks"] +__all__ = [ + "CartesianMotion", + "CylindricalMotion", + "TangentCartesianMotion", + "TangentCylindricalMotion", + "Observer", + "Tracker", + "Tracks", +] diff --git a/src/glimpse/track/motion.py b/src/glimpse/track/motion.py index 7b27970a..20a7c962 100644 --- a/src/glimpse/track/motion.py +++ b/src/glimpse/track/motion.py @@ -309,3 +309,214 @@ def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> Non time_units * particles[:, 3:6] + 0.5 * axyz * time_units ** 2 ) particles[:, 3:6] += time_units * axyz + + +class TangentCartesianMotion(Motion): + """ + Evolves particles tangent to a mean surface. + + Initial particle positions and velocities, and random accelerations, are + specified by independent and normally distributed x and y components. + Temporal arguments (e.g. :attr:`vxy`, :attr:`axy`) are assumed to be + in :attr:`time_unit` time units. + + Particle heights (z) are initialized based on a mean surface (:attr:`dem`) + and its uncertainty (:attr:`dem_sigma`), then maintain the same relative + height adjusted by a random walk proportional to the particle's horizontal + distance and a characteristic slope of small-scale features + (:attr:`slope_sigma`). + + This is the model used in Brinkerhoff (2017): Bayesian methods in glaciology + (http://hdl.handle.net/11122/8113), chapter 4. + + Attributes: + xy (iterable): Mean initial position (x, y). + time_unit (timedelta): Length of time unit for temporal arguments. + dem: Elevation of the surface on which to track points. + dem_sigma: Elevation standard deviations. + n (int): Number of particles. + xy_sigma (iterable): Standard deviation of initial position (x, y). + vxy (iterable): Mean initial velocity (dx/dt, dy/dt). + vxy_sigma (iterable): Standard deviation of initial velocity + (dx/dt, dy/dt). + axy (iterable): Mean acceleration (d^2x/dt^2, d^2y/dt^2). + axy_sigma (iterable): Standard deviation of acceleration + (d^2x/dt^2, d^2y/dt^2). + slope_sigma (float): Standard deviation of the characteristic slope + of small-scale features. + """ + + def __init__( + self, + xy: Iterable[Number], + time_unit: datetime.timedelta, + dem: Union[Number, Raster], + dem_sigma: Union[Number, Raster] = 0, + n: int = 1000, + xy_sigma: Iterable[Number] = (0, 0), + vxy: Iterable[Number] = (0, 0), + vxy_sigma: Iterable[Number] = (0, 0), + axy: Iterable[Number] = (0, 0), + axy_sigma: Iterable[Number] = (0, 0), + slope_sigma: Number = 0, + ) -> None: + self.xy = xy + self.time_unit = time_unit + if not isinstance(dem, Raster): + dem = Raster(dem) + self.dem = dem + if not isinstance(dem_sigma, Raster): + dem_sigma = Raster(dem_sigma) + self.dem_sigma = dem_sigma + self.n = n + self.xy_sigma = xy_sigma + self.vxy = vxy + self.vxy_sigma = vxy_sigma + self.axy = axy + self.axy_sigma = axy_sigma + self.slope_sigma = slope_sigma + + def initialize_particles(self) -> np.ndarray: + """ + Initialize particles around an initial mean position. + + Returns: + particles: Particle positions and velocities (x, y, z, vx, vy, vz). + """ + particles = np.zeros((self.n, 6), dtype=float) + particles[:, 0:2] = self.xy + self.xy_sigma * np.random.randn(self.n, 2) + z_offsets = self.dem_sigma.sample(particles[:, 0:2]) * np.random.randn(self.n) + particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets + particles[:, 3:5] = self.vxy + self.vxy_sigma * np.random.randn(self.n, 2) + return particles + + def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> None: + """ + Evolve particles through time by stochastic differentiation. + + Arguments: + particles: Particle positions and velocities (x, y, z, vx, vy, vz). + dt: Time step to evolve particles forward or backward. + """ + n = len(particles) + time_units = dt.total_seconds() / self.time_unit.total_seconds() + axy = self.axy + self.axy_sigma * np.random.randn(n, 2) + dxy = time_units * particles[:, 3:5] + 0.5 * axy * time_units ** 2 + # HACK: Recover z_offsets (since particles may have been resampled) + z_offsets = particles[:, 2] - self.dem.sample(particles[:, 0:2]) + z_offsets += ( + self.slope_sigma * np.random.randn(n) * (dxy ** 2).sum(axis=1) ** 0.5 + ) + particles[:, 0:2] += dxy + particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets + particles[:, 3:5] += time_units * axy + + +class TangentCylindricalMotion(Motion): + """ + Evolves particles tangent to a mean surface using cylindrical motion. + + Identical to :class:`TangentCartesianMotion`, except that particle motion is + specified by independent and normally distributed components of magnitude + (radius) and direction (theta). Angular arguments are + assumed to be in radians counterclockwise from the +x axis. + + Attributes: + xy (iterable): Mean initial position (x, y). + time_unit (timedelta): Length of time unit for temporal arguments. + dem: Elevation of the surface on which to track points. + dem_sigma: Elevation standard deviations. + n (int): Number of particles. + xy_sigma (iterable): Standard deviation of initial position (x, y). + vrth (iterable): Mean initial velocity (d radius/dt, theta). + vrth_sigma (iterable): Standard deviation of initial velocity + (d radius/dt, theta). + arth (iterable): Mean acceleration + (d^2 radius/dt^2, d theta/dt). + arth_sigma (iterable): Standard deviation of acceleration + (d^2 radius/dt^2, d theta/dt). + """ + + def __init__( + self, + xy: Iterable[Number], + time_unit: datetime.timedelta, + dem: Union[Number, Raster], + dem_sigma: Union[Number, Raster] = None, + n: int = 1000, + xy_sigma: Iterable[Number] = (0, 0), + vrth: Iterable[Number] = (0, 0), + vrth_sigma: Iterable[Number] = (0, 0), + arth: Iterable[Number] = (0, 0), + arth_sigma: Iterable[Number] = (0, 0), + slope_sigma: Number = 0, + ) -> None: + self.xy = xy + self.time_unit = time_unit + if not isinstance(dem, Raster): + dem = Raster(dem) + self.dem = dem + if not isinstance(dem_sigma, Raster): + dem_sigma = Raster(dem_sigma) + self.dem_sigma = dem_sigma + self.n = n + self.xy_sigma = xy_sigma + self.vrth = vrth + self.vrth_sigma = vrth_sigma + self.arth = arth + self.arth_sigma = arth_sigma + self.slope_sigma = slope_sigma + + def initialize_particles(self) -> np.ndarray: + """ + Initialize particles around an initial mean position. + + Returns: + particles: Particle positions and velocities (x, y, z, vx, vy, vz). + """ + particles = np.zeros((self.n, 6), dtype=float) + particles[:, 0:2] = self.xy + self.xy_sigma * np.random.randn(self.n, 2) + z_offsets = self.dem_sigma.sample(particles[:, 0:2]) * np.random.randn(self.n) + particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets + vrth = self.vrth + self.vrth_sigma * np.random.randn(self.n, 2) + particles[:, 3:5] = np.column_stack( + ( + # r' * cos(th) + vrth[:, 0] * np.cos(vrth[:, 1]), + # r' * sin(th) + vrth[:, 0] * np.sin(vrth[:, 1]), + ) + ) + return particles + + def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> None: + """ + Evolve particles through time by stochastic differentiation. + + Arguments: + particles: Particle positions and velocities (x, y, z, vx, vy, vz). + dt: Time step to evolve particles forward or backward. + """ + n = len(particles) + time_units = dt.total_seconds() / self.time_unit.total_seconds() + vx = particles[:, 3] + vy = particles[:, 4] + vr = np.sqrt(vx ** 2 + vy ** 2) + arth = self.arth + self.arth_sigma * np.random.randn(n, 2) + axy = np.column_stack( + ( + # r'' * cos(th) - r' * sin(th) * th' + arth[:, 0] * (vx / vr) - vy * arth[:, 1], + # r'' * sin(th) - r' * cos(th) * th' + arth[:, 0] * (vy / vr) + vx * arth[:, 1], + ) + ) + dxy = time_units * particles[:, 3:5] + 0.5 * axy * time_units ** 2 + # HACK: Recover z_offsets (since particles may have been resampled) + z_offsets = particles[:, 2] - self.dem.sample(particles[:, 0:2]) + z_offsets += ( + self.slope_sigma * np.random.randn(n) * (dxy ** 2).sum(axis=1) ** 0.5 + ) + particles[:, 0:2] += dxy + particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets + particles[:, 3:5] += time_units * axy